矩阵的LU分解初步：一个对角线上元素非零的方阵

上一篇我们对下三角矩阵的求解给出了一个方便的求解，利用消元代入可以在$\Theta (N^2)$ 的时间内完成，对于上三角矩阵，我们仍然可以利用类似的方法在相同的时间内求解。
对于一个非三角矩阵系统何对其进行求解，我们将在接下来几篇博文中进行讨论，而在这篇博里我们会对求解做一个浅显的分析和简易尝试。

---

【1】下三角矩阵的线性方程求解
【2】矩阵的LU分解初步：一个对角线上元素非零的方阵◀

---

对于如下的线性方程（其中矩阵对角线上的元素不为零，矩阵非奇异）
$$\left[\begin{array}{cccccc}{a}_{11}& a_{12}& .& .& .& a_{1n}\\ a_{21}& a_{22}& .& .& .& a_{2n}\\ .& .& .& & & .\\ .& .& & .& & .\\ .& .& & & .& .\\ a_{m1}& a_{m2}& a_{m3}& .& .& a_{mn}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ .\\ .\\ .\\ x_n\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\\ .\\ .\\ .\\ b_n\end{array}\right]$$
我们记$A=(a_{ij})$，$x=(x_i)$，$b=(b_i)$，于是我们可以简单地将方程记作$$Ax=b$$
为了求解这个方程，我们定义两个三角矩阵$L$、$U$其中，$L$为一个单位下三角矩阵，$U$为一个上三角矩阵。如果有$$LU=A$$成立，那么我们将这一对$L$、$U$称作$A$的一个$LU$分解，这样线性方程又可以记作$$LUx=b$$在这里，我们令$Ux=y$，这样上式变成$$Ly=b$$如此一来，好处显而易见：我们将一个非三角线性系统转化为两个三角线性系统，这样我们就可以采用上一篇博文中介绍的方法求解。
以上我们简要介绍了$LU$分解的思路。为了对$A$完成$LU$分解，我们可以尝试高斯消元(Gaussian elimination)的方式，实际上$LU$分解确实是一种特殊的高斯消元。

在《算法导论》(第三版，480页）我们可以看到其利用“主元”的方式（俺们小学时候解二元一次方程组的那种方法，同乘上一个系数，一减就完事了，这里写的好高大上啊…）参考其推导如下：

KaTeX parse error: Undefined control sequence: \ at position 281: …1 \end{bmatrix}\̲ ̲
而我们逆用矩阵乘法公式可以继续分解为两个矩阵相乘（可以试着验真）
$$\left[\begin{array}{cc}{a}_{11}& {\omega }^T\\ v& A^1\end{array}\right]=\left[\begin{array}{cc}1& 0\\ \frac{v}{a_{11}}& 1\end{array}\right]\left[\begin{array}{cc}{a}_{11}& {\omega }^T\\ 0& A^1-\frac{v{\omega }^T}{a_{11}}\end{array}\right]$$
这里我们采用$a_{11}$为主元，而我们又令$A^1-\frac{v{\omega }^T}{a_{11}}=L^1{U}^1$，从而又可以得到
$$A=\left[\begin{array}{cc}1& 0\\ \frac{v}{a_{11}}& L^1\end{array}\right]\left[\begin{array}{cc}{a}_{11}& {\omega }^T\\ 0& U^1\end{array}\right]$$
这个过程我们可以使用递归完成看出其递归性，我们能够使用递归算法来完成分解。递归是个老生常谈的问题，形式简洁但是使用起来可能会存在一些效率问题，我们可以将其重写为循环的形式，以下给出核心代码

	for (i = 0; i < Row - 1; i++)
	{
		for (j = i; j < Row - 1; j++)
			MatA[j + 1][i] /= MatA[i][i];
		for (j = i; j < Row - 1; j++)
		for (k = i; k < Row - 1; k++)
			MatA[j + 1][k + 1] -= MatA[i][k+1] * MatA[j+1][i];
	}

这段代码中我们得到了一个包含$LU$的新矩阵，实际上在这个矩阵中$L$和$U$以对角线分隔

    for (i = 0; i<Row; i++)
	{
		L[i][i] = 1;//L是一个单位下三角矩阵，因此对角线时1
		for (j = 0; j <=i; j++)
			L[i][j] = MatA[i][j];
	}
	
	for (i = 0; i<Row; i++)
	{
		for (j = 0; j <= i; j++)
			U[j][i] = MatA[j][i];
	}