楚列斯基分解法

楚列斯基分解是专门针对对称正定矩阵的分解。设 A = a i j ∈ R n × n A=a_{ij}\in R^{n\times n} A=aij​∈Rn×n是对称正定矩阵， A = R T R A = R^TR A=RTR称为矩阵A的楚列斯基分解、其中 R ∈ R n × n R\in R^{n\times n} R∈Rn×n是一个具有正的上三角矩阵。

R = [ r 11 r 12 r 13 r 14 r 22 r 23 r 24 r 33 r 34 r 44 ] R=\left[\begin{array}{cccc}r_{11} & r_{12} & r_{13} & r_{14} \\ & r_{22} & r_{23} & r_{24} \\ & & r_{33} & r_{34} \\ & & & r_{44}\end{array}\right] R=⎣⎢⎢⎡​r11​​r12​r22​​r13​r23​r33​​r14​r24​r34​r44​​⎦⎥⎥⎤​

在MATLAB中，实现这种分解的命令是chol.

A = [1 1 1 1;1 2 3 4;1 3 6 10;1 4 10 20];
[R,p] = chol(A)

R =

     1     1     1     1
     0     1     2     3
     0     0     1     3
     0     0     0     1


p =

     0

>> B = R'*R

B =

     1     1     1     1
     1     2     3     4
     1     3     6    10
     1     4    10    20

>> eig(A)

ans =

    0.0380
    0.4538
    2.2034

   26.3047

例：利用楚列斯基分解法求解 { 3 x 1 + 3 x 2 − 3 x 3 = 1 3 x 1 + 5 x 2 − 2 x 3 = 2 − 3 x 1 − 2 x 2 + 5 x 3 = 3 \left\{\begin{array}{c}3 x_{1}+3 x_{2}-3 x_{3}=1 \\ 3 x_{1}+5 x_{2}-2 x_{3}=2 \\ -3 x_{1}-2 x_{2}+5 x_{3}=3\end{array}\right. ⎩⎨⎧​3x1​+3x2​−3x3​=13x1​+5x2​−2x3​=2−3x1​−2x2​+5x3​=3​.

其基本思路是先对系数矩阵进行楚列斯基分解，即 A = R ′ R A=R'R A=R′R,然后解 R ′ y = b R'y=b R′y=b,最后解 R x = y Rx = y Rx=y.流程图如下：

首先利用普通的方法进行计算：

A = [3,3,-3;3,5,-2;-3,-2,5];
b = [1;2;3];
lamba = eig(A)

lamba =

    0.3096
    3.0000
    9.6904

R = chol(A)

R =

    1.7321    1.7321   -1.7321
         0    1.4142    0.7071
         0         0    1.2247

r = R'

r =

    1.7321         0         0
    1.7321    1.4142         0

   -1.7321    0.7071    1.2247

y = inv(r)*b

y =

    0.5774
    0.7071
    2.8577

x = inv(R)*y

x =

    3.3333

   -0.6667
    2.3333

考虑到其通用性，可以将其编写为一个函数：

function x= solvebyCHOL(A,b)
%此函数用于利用楚列斯基分解法解方程组
lambda = eig(A);   %求A的特征值
if lambda > eps&isequal(A,A')     %判断A是否为正定对称矩阵
    [n,n] = size(A);
    R = chol(A);
    %解R'y = b
    y(1) = b(1)/R(1,1);
    if n>1
        for i=2:n
            y(i) = (b(i)-R(1:i-1,i)'*y(1:i-1)')/R(i,i);
        end
    end
    %解Rx = y
    x(n) = y(n)/R(n,n);
    if n>1
        for i=n-1:-1:1
            x(i) = (y(i)-R(i,i+1:n)*x(i+1:n)')/R(i,i);
        end
    end
    x=x';
else
    x=[];
    disp('该方法是适用于对称正定矩阵');
end

验证程序的正确性：

x = solvebyCHOL(A,b)

x =

    3.3333

   -0.6667
    2.3333

可以看出该函数的运行结果和刚开始的运行结果一样。