LU 分解 （LU Decomposition）

LU分解是指将一个 NxN 矩阵 A 分解为一个上三角矩阵 U 和下三角矩阵 L 的过程， 即： ＬU＝Ａ。

比如我们可以将一个 ３ｘ３ 矩阵分解为：

 

如果我们需要求解方程 Ax = b，即求解 ＬU x = b。 那么令 Ux = y， 即求解 Ly=b， 得到y。接着求解Ux=y，得到x。由于L和U都是三角矩阵，极易使用追赶法得到解。在实际使用中，通常为了防止在分解过程中产生主元为零的情况，我们会带一个排列矩阵 P， 即： PA = LU。 

 

我们还可以使用LU分解来求矩阵的逆：设 A的逆为B， 那么  AB = I 即 A [b1, b2, …,bN] = [e1, e2, …,eN]
也就是说  A bj = ej； 其中 j = 1, 2, …, N。那么我们只有使用上面的解方程法解N次就可以求出逆矩阵的N列。

 

另外就是可以用来求行列式的值即det(A) = det(LU) = det(L)det(U).

 

下面我给出一个C++版本的实现：

 
 /*
　*　Decompose matrix: PA=LU
  *   
　*/
template< int M>
bool LUDecomposition::decompose(const Matrix<M, M> &matrix, 
	                            Matrix<M, M> &ll,
	                            Matrix<M, M> &uu,
	                            int pivot[M])
{
    uu = matrix;
    ll.identity();

    for (int iRow=0; iRow<M; ++iRow)
    {
        // find pivot rows
        pivot[iRow] = iRow;
        float maxValue = fabs(uu(iRow, iRow));
        for (int iNextRow=iRow+1; iNextRow<M; ++iNextRow)
        {
            if (maxValue + FLOAT_EPS < fabs(uu(iNextRow, iRow)) )
            {
                pivot[iRow] = iNextRow;
                maxValue = fabs(uu(iNextRow, iRow));
            }
        }

        // if the matrix is singular, return error
        if (almostZero(maxValue))
            return false;
        
        //  if the pivot row differs from the current row, then
        //  interchange the two rows.
        if (pivot[iRow] != iRow)
        {         
            uu.swapRows(iRow, pivot[iRow]);
        }
        
        // update lower matrix and upper matrix
        for (int iNextRow = iRow+1; iNextRow<M; ++iNextRow)
        {
            double factor = uu(iNextRow, iRow) / uu(iRow, iRow);
            // lower matrix
            ll(iNextRow, iRow) = (float)factor;
            
            // upper matrix
            uu(iNextRow, iRow) = 0.f;
            for (int iCol= iRow+1; iCol<M;++iCol)
            {
                uu(iNextRow, iCol) -= uu(iRow, iCol) * factor;
            }
        }
    }
    return true;
}

 /*
　*　Solve LUx=b
  *   
　*/
template<int M>
bool LUDecomposition::solve(const Matrix<M, M> &ll, 
                            const Matrix<M, M> &uu,
                            float bb[M],
                            float xx[M])
{
    // first, let y=Ux, solve Ly=b
    float yy[M];
    for (int ii=0; ii< M; ++ii)
    {
        yy[ii] = bb[ii];
        for(int jj=0; jj<ii; ++jj)
            yy[ii] -= yy[jj] * ll(ii, jj);
    }
    
    // then solve Ux = y
    for (int ii=M-1; ii>=0; --ii)
    {
        if (almostZero(uu(ii, ii)))
            return false;

        xx[ii] = yy[ii];

        for (int jj=ii+1; jj<M; ++jj)
        {
            xx[ii] -= xx[jj] * uu(ii, jj);
        }

        xx[ii] /=  uu(ii, ii);
    }
    return true;
}