LU分解指把一个矩阵分解成一个上三角和一个下三角相乘的形式。其中\(L\)指的是下三角矩阵，\(U\)指的是上三角矩阵。
\[L \cdot U = A\]
矩阵的形式如下：
\[
\begin{bmatrix} {\alpha}_{00} & 0 & 0 & 0 \\
{\alpha}_{10} & {\alpha}_{11} & 0 & 0 \\
{\alpha}_{20} & {\alpha}_{21} & {\alpha}_{22} & 0 \\
{\alpha}_{30} & {\alpha}_{31} & {\alpha}_{32} & {\alpha}_{33}\end{bmatrix} \cdot
\begin{bmatrix} {\beta}_{00} & {\beta}_{01} & {\beta}_{02} & {\beta}_{03} \\
0 & {\beta}_{11} & {\beta}_{12} & {\beta}_{13} \\
0 & 0 & {\beta}_{22} & {\beta}_{23} \\
0 & 0 & 0 & {\beta}_{33}\end{bmatrix} =
\begin{bmatrix} a_{00} & a_{01} & a_{02} & a_{03} \\
a_{10} & a_{11} & a_{12} & a_{13} \\
a_{20} & a_{21} & a_{22} & a_{23} \\
a_{30} & a_{31} & a_{32} & a_{33}\end{bmatrix}
\]
这样，求解线性方程组就可以分成两部完成：
\[
A \cdot x = \left( L \cdot U \right) \cdot x = L \cdot \left(U \cdot x \right) = L \cdot y = b
\]
通过\(b\)求\(y\)很容易。很明显有\(y_0 = \frac{b_0}{{\alpha}_{00}}\)，其余的\(y_i\)可以在知道前项值的情况下，逐次带入等式中求出：
\[
y_i = \frac{1}{{\alpha}_{ii}}\left[b_i – \sum_{j = 0}^{i-1} \alpha_{ij} y_j \right] \quad i = 1,2,\ldots,n-1
\]
然后,类似地，可以通过\(y\)求出\(x\)：
\[
x_{x-1} = \frac{y_{n-1}}{\beta_{n-1,n-1}}
\]
\[
x_i = \frac{1}{\beta_{ii}}\left[ y_i -\sum_{j=i+1}^{n-1} \beta_{ij} x_j \right] \quad i = n-2,N-3,\ldots,0
\]
所以，用\(LU\)分解求线性方程组的关键还是如何进行\(LU\)分解,一旦得到\(L\)，\(U\)矩阵，剩下的都很好解决。
用\(LU\)矩阵表示原矩阵，展开\(LU\)矩阵，我们可以得到\(n^2\)个方程，但是\(LU\)矩阵共有\(n^2+n\)个未知数。所以\(LU\)矩阵中未知数的个数实际上冗余。一个一般的做法是令$L$矩阵的对角线元素全为1；当然反过来令$U$的对角线元素全为1也可以。这样，很明显对\(L\)矩阵的第0列和\(U\)矩阵的第0行分别有\(\alpha_{i0} = \frac{a_{i0}}{\beta_{00}}\)，\(\beta_{0i} = a_{0i}\)。由矩阵乘法很明显可以得到以下等式：
\[
\beta_{ij}=a_{ij} – \sum_{k=0}^{i-1} \alpha_{ik} \beta_{kj} \tag{1.1}
\]
\[
\alpha_{ij} = \frac{1}{\beta_{jj}}\left(a_{ij} – \sum_{k=0}^{j-1} \alpha_{ik} \beta_{kj} \right) \tag{1.2}
\]
通过以上的公式可以看出，当我们知道较小的\(i\)，\(j\)下标对应的矩阵元素后可以求出下一个下标更大的矩阵元素。进而可以求出整个\(LU\)矩阵。这被称作求解\(LU\)矩阵的Crout算法。
因为\(LU\)分解中生成的两个矩阵中，未知数在矩阵中没有重合，所以两个矩阵可以存储在一个矩阵的存储空间里，如下所示：
\[
\begin{bmatrix} \beta_{00} & \beta_{01} & \beta_{02} & \beta_{03} \\
\alpha_{10} & \beta_{11} & \beta_{12} & \beta_{13} \\
\alpha_{20} & \alpha_{21} & \beta_{22} & \beta_{23} \\
\alpha_{30} & \alpha_{31} & \alpha_{32} & \beta_{33} \end{bmatrix}
\]
对\(LU\)分解中，\(\alpha\)元素的求取需要用到除法，所以可能会因为\(\beta_{ii}\)的值极端而造成误差过大甚至算法失效，所以选取主元是很有必要的。这同高斯消元法的具体做法很类似，可以通过行交换来选择合适的主元。通过保存行交换的次序，可以对线性方程组右端施以同样的交换，从而对等式的成立没有影响。
在作者写的程序中还有一处非常精妙的地方，作者对系数矩阵中的每一行都进行了比例变换。对比等式$1.1$和$1.2$，这两个等式非常相似，只是在右端相差了一个系数。如果我们对矩阵中的每一行进行比例变换，使得\(\beta_{jj}\)等于1，则这两个等式就可以用一个等式表示。程序写起来就会非常简单优美。因为\(\beta_{jj}\)已经被变换为1，是已知的确定值，所以每一行的变换系数可以存储在\(\beta_{jj}\)对应的存储位置（实际上就是\(\beta_{jj}\)原来的值。）
作者给出的代码如下。\(LU\)分解的原理很简单，这本书关于\(LU\)分解的代码，我觉得才是非常有意思的地方。当然代码写的感觉不是很好懂，作者把代码写成这样，真的是给学习算法的人看的么？但确实感觉很精妙就是了，也许对学习如何写代码有一定帮助。
[code lang=”cpp” title=”ludcmp.h”]struct LUdcmp
{
Int n;
MatDoub lu;
VecInt indx;
Doub d;
LUdcmp(MatDoub_I &a);
void solve(VecDoub_I &b, VecDoub_O &x);
void solve(MatDoub_I &b, MatDoub_O &x);
void inverse(MatDoub_O &ainv);
Doub det();
void mprove(VecDoub_I &b, VecDoub_IO &x);
MatDoub_I &aref;
};
[/code]
首先定义了一个\(LU\)分解类，这个类中n用来存储矩阵的大小，lu存储生成的lu矩阵，index用来存储矩阵中行的交换次序，d用来存储进行了奇数次还是偶数次交换（因为进行行变换时矩阵的行列式会取反，原矩阵被分解成两个三角矩阵的乘，而三角矩阵的行列式是对角线元素的乘积，所以进行LU变换后，原系数矩阵的行列式非常好求，）aref存储对原系数矩阵的引用。其余的函数包括构造函数，解方程函数，求逆矩阵函数，求行列式函数。
[code lang=”cpp” title=”ludcmp.h”]LUdcmp::LUdcmp(MatDoub_I &a) : n(a.nrows()), lu(a), aref(a), indx(n) {
const Doub TINY=1.0e-40;
Int i,imax,j,k;
Doub big,temp;
VecDoub vv(n);
d=1.0;
for (i=0;i<n;i++) {
big=0.0;
for (j=0;j<n;j++)
if ((temp=abs(lu[i][j])) > big) big=temp;
if (big == 0.0) throw("Singular matrix in LUdcmp");
vv[i]=1.0/big;
}
for (k=0;k<n;k++) {
big=0.0;
for (i=k;i<n;i++) {
temp=vv[i]*abs(lu[i][k]);
if (temp > big) {
big=temp;
imax=i;
}
}
if (k != imax) {
for (j=0;j<n;j++) {
temp=lu[imax][j];
lu[imax][j]=lu[k][j];
lu[k][j]=temp;
}
d = -d;
vv[imax]=vv[k];
}
indx[k]=imax;
if (lu[k][k] == 0.0) lu[k][k]=TINY;
for (i=k+1;i<n;i++) {
temp=lu[i][k] /= lu[k][k];
for (j=k+1;j<n;j++)
lu[i][j] -= temp*lu[k][j];
}
}
}[/code]
LUdcmp的构造函数，对原始矩阵进行LU分解有这个函数完成。首先初始化对象的变量。lu变量用来存储分解后的LU矩阵，但是因为对原始矩阵进行LU分解时，可以直接把分解的结果存储在原始矩阵里，所以可以先用原始矩阵初始化lu。a是原始矩阵的引用。index存储矩阵行交换的信息，长度要和矩阵的长度相同。
初始化之后，函数中首先定义了所有后面要用到的变量。这么写和Fortran类似，可能能方便编译器进行优化。这些变量的用处后面会说到。函数中首先有一个循环结构，用来找到每一行绝对值最大的元素，将这个最大绝对值存储在vv中相应的位置。用于后面挑选主元的时候，乘以相同行中可能被当做主元的元素。来判断应该和哪一行进行行交换。如果某一行所有的元素都是0，则原始矩阵很明显是奇异矩阵，函数抛出异常。
之后是另一个循环，这个循环就是用来进行LU分解。首先k层是控制最外围的循环。当k取某固定值时，相应的对角线上元素很明显只能在第k列中找。作者这里首先先把第k列中所有的元素除本行的最大值，然后查找除过之后的最大值。确定这个最大值所在的行，如果改行不是第k行，就用该行和第k行交换。相应的主元很明显就是相应行的第k个元素。至于为什么要除同一行的最大值，目的很明显是要使选取的主元和主元所在行最大元素比值尽量大，很明显这和直接查找最大元素做主元不太一样（尽管有些类似），这样是不是引起的误差会更小一些就不得而知了，也许将来能找到这么做的原因。
交换行之后，d取反，原因是行列式在交换行之后，行列式的值会取反（行列式的性质。）d用来保存行列式的正负性质。
之后index保存行交换的信息，用来解方程的时候，对等式右侧项做同样的交换。
之后是语句
[code lang=”cpp” gutter=”false”]if (lu[k][k] == 0.0) lu[k][k]=TINY;[/code]
如果主元是0的话，就令主元为最小值TINY。这里，和明显如果选定的主元是〇的话，原始矩阵应该是奇异矩阵。应该抛出异常才对，但作者似乎想让程序尽可能的能运行下去，在这种情况下令为〇的主元取一个非常小的值。实际上，语句的判断条件很难达到，对浮点数来说，即使主元真的是〇，但是因为误差的原因，也很难出现一个绝对等于〇的浮点数。
之后是这段代码真正关键的部分。
[code lang=”cpp” gutter=”false”]for (i=k+1;i<n;i++) {
temp=lu[i][k] /= lu[k][k];
for (j=k+1;j<n;j++)
lu[i][j] -= temp*lu[k][j];
}[/code]
当k循环到任一确定值时，第k行和第k列之前的行和列已经进行了彻底的LU分解，生成相应的LU矩阵的元素。
Visits: 1053