版权声明:本文为博主原创文章,未经博主允许不得转载。https://blog.csdn.net/weixin_38451800/article/details/88933683
1.Cholesky分解在MATLAB里面有专门的函数chol(其调用格式为:R=chol(X))这一函数功能就不介绍了,但有些场合只能用最基础的函数,所以,下面给出Cholesky分解的最基础形式。
%Cholesky分解法
function [X]=m_chol(A,b)
[N,N]=size(A);
X=zeros(N,1);
Y=zeros(N,1);
for i=1:N
A(i,i)=sqrt(A(i,i)-A(i,1:i-1)*A(i,1:i-1)');
if A(i,i)==0
'A is singular. no unique solution'
break
end
for j=i+1:N
A(j,i)=(A(j,i)-A(j,1:i-1)*A(i,1:i-1)')/A(i,i);
end
end
A
b
%前代法
for j=1:N
Y(j)=(b(j)-A(j,1:j-1)*Y(1:j-1))/A(j,j);
end
Y
%
A=A'
for k=N:-1:1
X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);
end
如果矩阵X是对称正定的,则Cholesky分解将矩阵X分解成一个下三角矩阵和上三角矩阵的乘积。设上三角矩阵为R,则下三角矩阵为其转置,即X=R’R。
另外,cholesky分解可能出现非正定问题, [R,p]=chol(X):这个命令格式将不输出出错信息。当X为对称正定的,则p=0,R与上述格式得到的结果相同;否则p为一个正整数。如果X为满秩矩阵,则R为一个阶数为q=p-1的上三角阵,且满足R’R=X(1:q,1:q)。 实现Cholesky分解后,线性方程组Ax=b变成R‘Rx=b,所以x=R(R’\b)。
2.下面给出两个博客讲解
(1)cholesky原理:https://blog.csdn.net/wangxiaojun911/article/details/7030964/
(2)C语言实现方式:https://blog.csdn.net/acdreamers/article/details/44656847