Opencv 三对角线矩阵(Tridiagonal Matrix)解法之(Thomas Algorithm)

1. 简介

三对角线矩阵(Tridiagonal Matrix)，结构如公式(1)所示：

aixi−1+bixi+cixx+1=di(1)

其中a1=0，cn=0。写成矩阵形式如(2):

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢b1a20c1b2a3c2b3⋱⋱⋱cn−1an0bn⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢x1x2x3⋮xn⎤⎦⎥⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢d1d2d3⋮dn⎤⎦⎥⎥⎥⎥⎥⎥⎥(2)

常用的解法为Thomas algorithm，又称为The Tridiagonal matrix algorithm(TDMA). 它是一种高斯消元法的解法。分为两个阶段：向前消元（Forward Elimination）和回代（Back Substitution）。

• 向前消元（Forward Elimination）：
  

  c′i=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪cibicibi−aic′i−1;i=1;i=2,3,…,n−1(3)
  
  
  d′i=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪dibidi−aid′i−1bi−aic′i−1;i=1;i=2,3,…,n.(4)

• 回代（Back Substitution）：
  

  xn=d′nxi=d′i−c′ixi+1;i=n−1,n−2,…,1.(5)

2.代码

• 维基百科提供的C语言版本：

void solve_tridiagonal_in_place_destructive(float * restrict const x, const size_t X, const float * restrict const a, const float * restrict const b, float * restrict const c) {    /*     solves Ax = v where A is a tridiagonal matrix consisting of vectors a, b, c     x - initially contains the input vector v, and returns the solution x. indexed from 0 to X - 1 inclusive     X - number of equations (length of vector x)     a - subdiagonal (means it is the diagonal below the main diagonal), indexed from 1 to X - 1 inclusive     b - the main diagonal, indexed from 0 to X - 1 inclusive     c - superdiagonal (means it is the diagonal above the main diagonal), indexed from 0 to X - 2 inclusive     Note: contents of input vector c will be modified, making this a one-time-use function (scratch space can be allocated instead for this purpose to make it reusable)     Note 2: We don't check for diagonal dominance, etc.; this is not guaranteed stable     */    /* index variable is an unsigned integer of same size as pointer */    size_t ix;    c[0] = c[0] / b[0];    x[0] = x[0] / b[0];    /* loop from 1 to X - 1 inclusive, performing the forward sweep */    for (ix = 1; ix < X; ix++) {        const float m = 1.0f / (b[ix] - a[ix] * c[ix - 1]);        c[ix] = c[ix] * m;        x[ix] = (x[ix] - a[ix] * x[ix - 1]) * m;    }    /* loop from X - 2 to 0 inclusive (safely testing loop condition for an unsigned integer), to perform the back substitution */    for (ix = X - 1; ix-- > 0; )        x[ix] = x[ix] - c[ix] * x[ix + 1];}

• 本人基于Opencv的版本：

bool caltridiagonalMatrices(     cv::Mat_<double> &input_a,     cv::Mat_<double> &input_b,     cv::Mat_<double> &input_c,    cv::Mat_<double> &input_d,    cv::Mat_<double> &output_x ){    /*     solves Ax = v where A is a tridiagonal matrix consisting of vectors input_a, input_b, input_c, and v is a vector consisting of input_d.     input_a - subdiagonal (means it is the diagonal below the main diagonal), indexed from 1 to X - 1 inclusive     input_b - the main diagonal, indexed from 0 to X - 1 inclusive     input_c - superdiagonal (means it is the diagonal above the main diagonal), indexed from 0 to X - 2 inclusive     input_d - the input vector v, indexed from 0 to X - 1 inclusive     output_x - returns the solution x. indexed from 0 to X - 1 inclusive     */    /* the size of input_a is 1*n or n*1 */    int rows = input_a.rows;    int cols = input_a.cols;    if ( ( rows == 1 && cols > rows ) ||         (cols == 1 && rows > cols ) )    {        const int count = ( rows > cols ? rows : cols ) - 1;        output_x = cv::Mat_<double>::zeros(rows, cols);        cv::Mat_<double> cCopy, dCopy;        input_c.copyTo(cCopy);        input_d.copyTo(dCopy);        if ( input_b(0) != 0 )        {            cCopy(0) /= input_b(0);            dCopy(0) /= input_b(0);        }        else        {            return false;        }        for ( int i=1; i < count; i++ )        {            double temp = input_b(i) - input_a(i) * cCopy(i-1);            if ( temp == 0.0 )            {                return false;            }            cCopy(i) /= temp;            dCopy(i) = ( dCopy(i) - dCopy(i-1)*input_a(i) ) / temp;        }        output_x(count) = dCopy(count);        for ( int i=count-2; i > 0; i-- )        {            output_x(i) = dCopy(i) - cCopy(i)*output_x(i+1);        }        return true;    }    else    {        return false;    }}

参考文献：https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm