数值计算方法(Numerical Methods)MATLAB实现(2)--choleski分解

1. 乔列斯基分解
      乔列斯基分解当且仅当矩阵为正定对称阵时才可进行，运算量大约是LU分解的一半(n³/6)，另有n次取平方根的计算(n为方阵维度)。
      对于3*3矩阵A = LL^T。根据对称性，表示为

      将右侧乘入，得到

      现需根据以上矩阵求出每个L。由于对称，只需利用下三角产生的六个方程即可。
      先从上至下计算第一列。得出:
      $$A_{11}=L_{11}^2{L}_{11}=\sqrt{A_{11}}$$      $$A_{21}=L_{11}{L}_{21}{L}_{21}=A_{21}/L_{11}$$      $$A_{31}=L_{11}{L}_{31}{L}_{31}=A_{31}/L_{11}$$根据对称性，第二列只需从第二行开始计算。$$A_{22}=L_{21}^2+L_{22}^2{L}_{22}=\sqrt{A_{22}-L_{21}^2}$$$$A_{32}=L_{21}{L}_{31}+L_{22}{L}_{32}{L}_{32}=(A_{32}-L_{21}{L}_{31})/L_{22}$$最后从第三行求出L₃₃即可。$$A_{33}=L_{31}^2+L_{32}^2+L_{33}^2{L}_{33}=\sqrt{A_{33}-L_{31}^2-L_{32}^2}$$从三阶的情况可以看出，$$(L{L}^T{)}_{ij}=L_{i1}{L}_{j1}+L_{i2}{L}_{j2}+...+L_{ij}{L}_{jj}=\sum _{k=1}^j{L}_{ik}{L}_{jk},i\ge j$$将此式与对应A中的元素对应，得到$$A_{ij}=\sum _{k=1}^j{L}_{ik}{L}_{jk},i=j,j+1,...,n,j=1,2,...,n$$下标表示计算发生在下三角形矩阵。对于j=1(第一列)，有$$L_{11}=\sqrt{A_{11}}{L}_{i1}=A_{i1}/L_{i1},i=2,3,......,n$$推广可得，进行到A_ij得到的代求值应该是L_ij，(此时该行方程的其他需要求的L值都已求)。将含有L_ij的一项从和式中取出：$$A_{ij}=\sum _{k=1}^{j-1}{L}_{ik}{L}_{jk}+L_{ij}{L}_{jj}$$如果待求L在对角线上(i=j),L_jj的解为$$L_{jj}=\sqrt{A_{jj}-\sum _{k=1}^{j-1}{L}_{jk}^2},j=2,3,...,n$$若不是对角线元素，那么$$L_{ij}=(A_{ij}-\sum _{k=1}^{j-1}{L}_{ik}{L}_{jk})/L_{jj},j=2,3,....,n-1,i=j+1,j+2,.....,n$$
        分解阶段代码
        从以上分析中看出，A_ij仅在计算L_ij处使用。因此，若L_ij已被计算，A_ij就没有用了。这样就可以在程序中用L覆盖下三角矩阵中的A值。如果遇到Ljj²取负的情况，程序报错，终止。

function L = choleski(A)
%computes L in Coleski's decomposition A = LL'
n =size(A,1);
for j = 1:n
    temp = A(j,j) -dot(A(j,1:j-1),A(j,1:j-1));
    if temp < 0.0
        error("Matrix is not positive definite");
    end
    A(j,j)=sqrt(temp);
    for i =j+1:n
        A(i,j)=(A(i,j) - dot(A(i,1:j-1),A(j,1:j-1)))/A(j,j);
    end
end
L = tril(A);
end

      求解阶段代码
      与上篇Doolittle分解相似的代换方式即可解出:

function x = choleskiSol(L,b)
%solves [L]{L'} = {b}

n = length(b);
if size(b,2) >1;b=b';end %{b} must be column vector
for k = 1:n %Solution of [L]{y} = {b}
    b(k) = (b(k)-dot(L(k,1:k-1),b(1:k-1)'))/L(k,k);
end
for k = n:-1:1 %Solution of {L}'{x} = {y}
    b(k) = (b(k) - dot(L(k+1:n,k),b(k+1:n)))/L(k,k);
end
x=b;

      例如，使用乔列斯基分解法求解Ax=b并验证,其中，
$$\text{A}=\left(\begin{array}{cccc}1.44& -0.36& 5.52& 0.00\\ -0.36& 10.33& -7.78& 0.00\\ 5.52& -7.78& 28.40& 9.00\\ 0.00& 0.00& 9.00& 61.00\end{array}\right),\text{b}=\left(\begin{array}{c}0.04\\ -2.15\\ 0.00\\ 0.88\end{array}\right)$$    完整程序如下：

function L = choleski(A)
%computes L in Coleski's decomposition A = LL'
n =size(A,1);
for j = 1:n
    temp = A(j,j) -dot(A(j,1:j-1),A(j,1:j-1));
    if temp < 0.0
        error("Matrix is not positive definite");
    end
    A(j,j)=sqrt(temp);
    for i =j+1:n
        A(i,j)=(A(i,j) - dot(A(i,1:j-1),A(j,1:j-1)))/A(j,j);
    end
end
L = tril(A);
end

function x = choleskiSol(L,b)
%solves [L]{L'} = {b}

n = length(b);
if size(b,2) >1;b=b';end %{b} must be column vector
for k = 1:n %Solution of [L]{y} = {b}
    b(k) = (b(k)-dot(L(k,1:k-1),b(1:k-1)'))/L(k,k);
end
for k = n:-1:1 %Solution of {L}'{x} = {y}
    b(k) = (b(k) - dot(L(k+1:n,k),b(k+1:n)))/L(k,k);
end
x=b;

      在命令行中输入

A = [1.44 -0.36 5.52 0.00;
-0.36 10.33 -7.78 0.00;
5.52 -7.78 28.40 9.00;
0.00 0.00 9.00 61.00];
b = [0.04 -2.15 0 0.88];
L = choleski(A);
x = choleskiSol(L,b)
Check = A*x

验证结果如下：

      结果正确。