基于Cholesky分解的正定矩阵求逆矩阵

转载地址   http://379910987.blog.163.com/blog/static/3352379720136102212297/


前面的博客中我提到了如何实现正定矩阵的Cholesky分解,并提供了源代码,通过该代码可以将一个正定矩阵分解为一个上三角矩阵和其转置的乘积,在此基础上,对上三角矩阵进行求逆是十分简单的运算,在得到其逆矩阵之后,也就能求出原正定矩阵的逆矩阵了。

数学原理如下:

基于Cholesky的正定矩阵求逆矩阵 - Lemniscate - 信息,灵感,创新

 对于u的逆矩阵,可以使用下列函数进行计算:

function ut=invutri(u)
   [m,n]=size(u);
    if m~=n
       error( 'Error  using in invutri.Matrix must be square. ');
   end
    for i= 1:n
        for j= 1:i- 1
            if u(i,j)~= 0
               error( 'Error  using in invutri.Matrix must be up-triangle. ');
           end
       end
   end
    for i= 1:n
        if u(i,i)== 0
           error( 'Error  using in invutri.Matrix  is singular. ')
       end
   end
   mu=[u,eye(n)];
    for i=n:- 1: 1
       mu(i,:)=mu(i,:)/mu(i,i);
        for j=i- 1:- 1: 1
           mu(j,:)=mu(j,:)-mu(i,:)*mu(j,i);
       end
   end
   ut=mu(:,n+ 1: 2*n);
end

这样,结合前面的工作,我们就得到了两个m函数文件:cholesky.m和invutri.m文件,分别完成正定矩阵的Cholesky分解和上三角矩阵的求逆。接下来,我们实现正定矩阵的求逆过程,只需要三行代码即可完成工作:

u=cholesky(A);

tu=invutri(u);

B=tu*tu';

则B就是正定矩阵A的逆矩阵。

测试代码如下:

>> clear
>> A=[1 2 3 4;2  5 9 10;3 9 22 20;4 10 20 37]

A =

     1     2     3     4
     2     5     9    10
     3     9    22    20
     4    10    20    37

>> u=cholesky(A)

u =

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

>> tu=invutri(u)

tu =

    1.0000   -2.0000    1.5000   -0.3750
         0    1.0000   -1.5000   -0.1250
         0         0    0.5000   -0.1250
         0         0         0    0.2500

>> B=tu*tu'

B =

    7.3906   -4.2031    0.7969   -0.0938
   -4.2031    3.2656   -0.7344   -0.0313
    0.7969   -0.7344    0.2656   -0.0313
   -0.0938   -0.0313   -0.0313    0.0625

>> B*A

ans =

     1     0     0     0
     0     1     0     0
     0     0     1     0
     0     0     0     1

>> A*B

ans =

     1     0     0     0
     0     1     0     0
     0     0     1     0
     0     0     0     1

表明B确实是A的逆矩阵

这几段代码都是使用的MATLAB中的基本语句,而没有使用inv、chol等MATLAB中现成的函数,从数学原理上分析操作过程。但是也要小心的是,在使用时,cholesky函数并没有检验矩阵是否是正定的,所以计算结果有问题的时候,最好看看矩阵是否是正定的。




  • 0
    点赞
  • 0
    评论
  • 6
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

©️2021 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值