LU分解是指将一个 NxN 矩阵 A 分解为一个上三角矩阵 U 和下三角矩阵 L 的过程, 即: LU=A。
比如我们可以将一个 3x3 矩阵分解为:
如果我们需要求解方程 Ax = b,即求解 LU 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;
}