解线性方程组——Cholesky分解 | 北太天元

一、问题描述

对于线性方程组
$$Ax=b,A=\left(\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right),b=\left(\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right)$$

A为对称正定矩阵，求向量$x$

二、对称正定矩阵的Cholesky分解

$$\left(\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right)=\left[\begin{array}{cccc}{l}_{11}& & & \\ l_{21}& l_{22}& & \\ ⋮& ⋮& \ddots & \\ l_{n1}& l_{n2}& \cdots & l_{nn}\end{array}\right]\left[\begin{array}{cccc}{l}_{11}& l_{21}& \cdots & l_{n1}\\ & l_{22}& \cdots & l_{n2}\\ & & \ddots & ⋮\\ & & & l_{nn}\end{array}\right].$$

上式可以看作，做第一步
$$[\begin{array}{cc}{a}_{11}& A_{21}^{\mathrm{T}}\\ A_{21}& A_{22}\end{array}]=\left[\begin{array}{cc}{l}_{11}& \\ L_{21}& L_{22}\end{array}\right]\left[\begin{array}{cc}{l}_{11}& L_{21}^{\mathrm{T}}\\ & L_{22}^{\mathrm{T}}\end{array}\right]=\left[\begin{array}{cc}{l}_{11}^2& l_{11}{L}_{21}^{\mathrm{T}}\\ l_{11}{L}_{21}& L_{21}{L}_{21}^{\mathrm{T}}+L_{22}{L}_{22}^{\mathrm{T}}\end{array}\right]$$

按顺序，求出 $l_{11},L_{21}$后，再得$L_{22}{L}_{22}^T$，而$L_{22}{L}_{22}^T$依然是对称正定的，继续按照第一步的思路进行计算。

举一个例子
$$\left[\begin{array}{ccc}4& 12& -16\\ 12& 37& -43\\ -16& -43& 98\end{array}\right]$$
对这个三维矩阵做Cholesky分解，自己试着算一下

三、算法

💖 Cholesky分解：[L]=Cholesky_fac(A)

输入：系数矩阵 A

输出：下三角矩阵L

实现过程：在原矩阵上做改变，

• for $i=1$ to $n-1$
• $a_{ii}=\sqrt{A_{ii}}$
• $A_{i+1:n,i}=A_{i+1:n,i}\ast \frac{1}{a_{ii}}$
• $A_{i+1:n,i+1:n}=A_{i+1:n,i+1:n}-A_{i+1:n,i}\ast A_{i+1:n,i}^T$

i = n 时，$A_{n,n}=\sqrt{A_{n,n}}$

四、北太天元源程序

Cholesky分解

function [L] = Cholesky_fac(A)
    % 对称正定矩阵 的 Cholesky分解 
    % A = LL'
    % 输入:
    %    A,对称正定
    % 输出: 
    %    分解后的 下三角 L 
    % 创建时间: 1/18/2024
    % 版本: 1.0
    n = length(A);
    for i = 1:1:n-1
        A(i,i) = sqrt(A(i,i));
        A(i+1:n,i) = A(i+1:n,i)/A(i,i);
        A(i+1:n,i+1:n) = A(i+1:n,i+1:n)-A(i+1:n,i)*A(i+1:n,i)';
    end
    A(n,n)= sqrt(A(n,n));
    L = tril(A);
end

将上述代码保存为 Cholesky_fac.m 文件

五、数值算例

例子1：

$$\left[\begin{array}{ccc}4& 12& -16\\ 12& 37& -43\\ -16& -43& 98\end{array}\right]$$
对其做Cholesky分解

%%
clc,clear all;
A1 = [4 12 -16;12 37 -43;-16 -43 98];
L1 = Cholesky_fac(A1)

得

L1 = 
   2   0   0
   6   1   0
  -8   5   3

例子2：

$$\left[\begin{array}{ccccc}10& 1& 2& 3& 4\\ 1& 9& -1& 2& -3\\ 2& -1& 7& 3& -5\\ 3& 2& 3& 12& -1\\ 4& -3& -5& -1& 15\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\end{array}\right]=\left[\begin{array}{c}12\\ -27\\ 14\\ -17\\ 12\end{array}\right]$$

用Cholesky分解求解向量 x

先用Cholesky分解求出 L2，再用两次回代得到 $x$

%%  对称正定矩阵下
A2 = [10 1 2 3 4; 1 9 -1 2 -3; 2 -1 7 3 -5; 3 2 3 12 -1; 4 -3 -5 -1 15];
b = [12; -27; 14; -17; 12];
L2 = Cholesky_fac(A2)
x = back_substitution_two(L2,L2',b)
x2 = gsem_column(A2,b)

得

L2 = 
    3.1623    0.0000    0.0000    0.0000    0.0000
    0.3162    2.9833    0.0000    0.0000    0.0000
    0.6325   -0.4022    2.5374    0.0000    0.0000
    0.9487    0.5698    1.0362    3.1147    0.0000
    1.2649   -1.1397   -2.4665    0.3227    2.4317

x = 
    1.0000
   -2.0000
    3.0000
   -2.0000
    1.0000

x2 = 
    1.0000
   -2.0000
    3.0000
   -2.0000
    1.0000

假如对非对称正定矩阵使用

%% 对称非正定矩阵下 结果不正确
A3 = [-10 1 2 3 4; 1 -5 -1 2 -3; 2 -1 7 3 -5; 3 2 3 12 -1; 4 -3 -5 -1 15];
b = [12; -27; 14; -17; 12];
L3 = Cholesky_fac(A3) %出现虚数
x31 = back_substitution_two(L3,L3',b)

x32 = gsem_column(A3,b)

运行后得

L3 = 
   0.0000 + 3.1623i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 - 0.3162i   0.0000 + 2.2583i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 - 0.6325i   0.0000 + 0.5314i   2.5135 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 - 0.9487i   0.0000 - 0.7528i   1.1140 + 0.0000i   3.0483 + 0.0000i   0.0000 + 0.0000i
   0.0000 - 1.2649i   0.0000 + 1.5055i  -2.6258 + 0.0000i   0.6097 + 0.0000i   1.9664 + 0.0000i


x31 = 
   8.1524 + 0.0000i
 -23.9180 + 0.0000i
  23.8006 + 0.0000i
  -6.7427 + 0.0000i
  16.5172 + 0.0000i


x32  是正确结果
x32 = 
    0.2250
    1.6622
    5.2829
   -2.8504
    2.6434

可以看到 L3 中出现虚数，显然不正确。