上一篇我们对下三角矩阵的求解给出了一个方便的求解,利用消元代入可以在
的时间内完成,对于上三角矩阵,我们仍然可以利用类似的方法在相同的时间内求解。
对于一个非三角矩阵系统何对其进行求解,我们将在接下来几篇博文中进行讨论,而在这篇博里我们会对求解做一个浅显的分析和简易尝试。
【1】下三角矩阵的线性方程求解
【2】矩阵的LU分解初步:一个对角线上元素非零的方阵◀
对于如下的线性方程(其中矩阵对角线上的元素不为零,矩阵非奇异)
我们记 , , ,于是我们可以简单地将方程记作
为了求解这个方程,我们定义两个三角矩阵 、 其中, 为一个单位下三角矩阵, 为一个上三角矩阵。如果有
以上我们简要介绍了 分解的思路。为了对 完成 分解,我们可以尝试高斯消元(Gaussian elimination)的方式,实际上 分解确实是一种特殊的高斯消元。
在《算法导论》(第三版,480页)我们可以看到其利用“主元”的方式(俺们小学时候解二元一次方程组的那种方法,同乘上一个系数,一减就完事了,这里写的好高大上啊…)参考其推导如下:
而我们逆用矩阵乘法公式可以继续分解为两个矩阵相乘(可以试着验真)
这里我们采用 为主元,而我们又令 ,从而又可以得到
这个过程我们可以使用递归完成看出其递归性,我们能够使用递归算法来完成分解。递归是个老生常谈的问题,形式简洁但是使用起来可能会存在一些效率问题,我们可以将其重写为循环的形式,以下给出核心代码
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];
}
这段代码中我们得到了一个包含 的新矩阵,实际上在这个矩阵中 和 以对角线分隔
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];
}