：科学和工程计算中的线性方程组系数矩阵(图)

在科学和工程计算中平方根公式，经常需要求解形如 Ax=b 的线性方程组，其中 A 是 ntimes m 矩阵，称为系数矩阵，b 是 n 维列向量，称为 a右手向量，x为要求解的m维列向量平方根公式，称为解向量。

在科学和工程的实际计算中，经常会遇到系数矩阵A是一个对称正定矩阵。如果

A=begin{} a_{11} \ a_{21} & a_{22} && 对称\ a_{31} & a_{32} & a_{33} \ vdots & vdots & vdots & ddots \ a_{n1} & a_{n2} & dots & a_{nn} end{} \

对于正定矩阵，有如下三角矩阵

L = begin{} l_{11} \ l_{21} & l_{22} & & bm 0\ l_{31} & l_{32} & l_{33} \ vdots & vdots & vdots & ddots \ l_{n1} & l_{n2} & l_{n3} & dots & l_{nn} end{} \

令 A = L cdot L^{T} 成立。如果 L 的主要对角元素取正值，则这种分解是唯一的。

直接展开矩阵关系A = LL^{T}，有

begin{align} a_{11} &= l_{11}^{2} \ a_{21} &= l_{21}l_{11},quad a_{22} = l_{21} ^{2}+l_{22}^{2}\a_{31} &= l_{31}l_{11},quad a_{32} = l_{31}l_{21}+l_{32} l_{22},quad a_{33}=l_{31}^{2}+l_{32}^{2}+l_{33}^{2}\dots end{align}\

据此，矩阵 L 的元素 l_{11} to l_{21} to l_{22} to l_{31} to l_{32} to dots 可以通过行，计算公式为

begin{cases} l_{ij} &= (a_{ij} - sum{k=1}^{j-1}l_{ik}l_{jk}) / l_{jj}, 四边形 & j = 1, 2, dots, i-1 \ l_{ii} &= (a_{ii} - sum{k=1}^{i-1}l_{ik}^2)^frac{1}{2}, quad & i = 1, 2, dots, n \ end{cases} \

基于矩阵分解 A = LL^{T} ，对称正定方程组 Ax = b 可以简化为两个三角方程 Ly = b 和 L^Tx = y 来求解。由 Ly = b 即

begin{cases} l_{11}y_{1} &= b_1 \ l_{21}y_{2} + l_{22}y_{2} &= b_2 \ dots dots dots l_{n1}y_{1} + l_{n2}y_{2} + dots + l_{nn}y_{n} &= b_n end{cases} \

y_1 to y_2 to dots to y_n 可以按顺序计算：

y_i = (b_i - sum{k=1}^{i-1}l_{ik}y_{k})/l_{ii}, quad i = 1, 2, dots,n

而 L^Tx = y 就是

begin{cases} begin{*}{2} l_{11}x_1 + l_{21}x_2 + dots + l_{n1}x_n &= y_1 \ l_{22}x_2 + dots + l_{n2}x_n &= y_2 \ dots dots dots \ l_{nn}x_n &= y_n \ end{*} end{cases} \

找到 x_n to x_{n-1} to dots to x_1 的可逆顺序：

x_i = (y_i - sum{k=i+1}^nl_{ki}x_{k})/l_{ii}, quad i = n, n-1, dots, 1 \

由于公式中包含矩阵分解时的平方根运算，所以该算法称为平方根法，也称为分解法。

根据上面的公式，可以通过编写程序求解方程：

subroutine cholesky_full(n, a, y)
    implicit none
    
    integer, intent(in) :: n
    real, intent(inout) :: a(n, n), y(n)
    
    integer :: i, j, k
    real :: temp
    
    ! 分解矩阵，生成下三角阵L
    ! 工程问题中的很多矩阵非常庞大，所以，计算过程中的数据应该直接存放在原始数组a中，
    ! 而不是新创建一个数组
    do i = 1, n
        ! 公式中，j的取值范围为1到j-1，此处换成1到j，可以将分解式统一起来，省去一次判断。
        ! 因为j=i时，j循环虽然会执行错误的操作、生成错误的a(i,j)结果，
        ! 但a(i,j)马上就会被最外层的i循环生成的正确数据替换
        do j = 1, i
            temp = a(i, j)
            do k = 1, j-1

                temp = temp - a(i, k) * a(j, k)
            end do
            a(i, j) = temp / a(j, j)
            ! a(j, i) = 0.  ! 对角线上方d的元素赋0，可有可无
        end do
        a(i, i) = sqrt(temp)
    end do
    
    ! 根据Ly=b求解出y
    do i = 1, n
        temp = y(i)
        do j = 1, i-1
            temp = temp - a(i, j) * y(j)
        end do
        y(i) = temp / a(i, i)
    end do
    
    ! 求解出x
    do i = n, 1, -1
        temp = y(i) / a(i, i)
        y(i) = temp
        ! 公式中k的范围为i+1到n，此处为1到i-1，因为下方a(i,k)的下标和公式中交换了顺序
        do k = 1, i-1
            y(k) = y(k) - temp * a(i, k)
        end do
    end do
end subroutine

上面代码的分解部分和前面的公式基本一样，容易理解，只是引入了一个临时变量temp来存储数据。但是，如果我们交换j和k两层循环的位置，然后稍微调整循环计数器的取值范围，就可以直接完成分解操作，而无需使用临时变量。代码如下：

do i = 1, n
    do k = 1, i - 1
        a(i, k) = a(i, k) / a(k, k)
        do j = k + 1, i
            a(i, j) = a(i, j) - a(i, k) * a(j, k)
        end do
    end do
    a(i, i) = sqrt(a(i, i))
end do

参考资料：

王能超。高等学校教材，数值分析简明教程，（第2版）[M]． 2003.

吴建平、王正华、李晓梅。稀疏线性方程组的高效求解与并行计算[M]．湖南科学技术出版社，2004.