注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

信息 灵感 创新

I? =Information,Inspiration,Innovation

 
 
 

日志

 
 
关于我

we are 5. Mathematics, Computation, Programming, Engineering, and Making fun of life.

网易考拉推荐

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

2013-07-10 14:40:18|  分类: M&M |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

前面的博客中我提到了如何实现正定矩阵的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函数并没有检验矩阵是否是正定的,所以计算结果有问题的时候,最好看看矩阵是否是正定的。
 

  评论这张
 
阅读(5741)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2016