MATLAB之楚列斯基分解法（九）

楚列斯基分解法

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

$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]$

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

命令	说明
R=chol(A)	返回A的楚列斯基分解
[R,p]=chol(A)	同上，p是用来判断A是否正定。

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

例：利用楚列斯基分解法求解$\{\begin{array}{c}3x_1+3x_2-3x_3=1\\ 3x_1+5x_2-2x_3=2\\ -3x_1-2x_2+5x_3=3\end{array}$.

其基本思路是先对系数矩阵进行楚列斯基分解，即$A=R^{\prime }R$,然后解$R^{\prime }y=b$,最后解$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

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