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

Mr.Right

不顾一切的去想,于是我们有了梦想。脚踏实地的去做,于是梦想成了现实。

 
 
 

日志

 
 
关于我

人生一年又一年,只要每年都有所积累,有所成长,都有那么一次自己认为满意的花开时刻就好。即使一时不顺,也要敞开胸怀。生命的荣枯并不是简单的重复,一时的得失不是成败的尺度。花开不是荣耀,而是一个美丽的结束,花谢也不是耻辱,而是一个低调的开始。

网易考拉推荐

阿英讲如何编写共轭梯度法, matlab实现  

2013-11-02 23:22:18|  分类: 编程 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
即以此功德,庄严佛净土。上报四重恩,下济三途苦。惟愿见闻者,悉发菩提心。在世富贵全,往生极乐国。

缘起:共轭梯度法是解大规模正定线性方程组的有效方法之一,如何用程序实现是初学者遇到的困扰之一。周围愿意分享的人太少了,我编程过程中是多么希望有个高手给我指点下,节省我的摸索时间,可惜一直没有这样的人,这另我感到有点孤单和无助。今晚花了一晚上专门研究了下,将完整的代码贴出如下,希望对大家有帮助,节省初学者的摸索时间。

特别感谢Stanford University的Emmanuel Candes教授的指点!
阿英讲如何编写共轭梯度法, matlab实现 - 阿英 - Mr.Right
 

% Solve the positive definite system Ax = b by the CG method.

n = 400; % Problem size
[Q, R] = qr(randn(n)); % Make random orthonormal matrix
lmin = 1; % minimum eigenvalue
lmax = 10; % maximum eigenvalue
lambda = 1 + (lmax-lmin).*rand(n, 1); % Choose eigenvalues clustered in [1, 10]
A = Q*diag(lambda)*Q'; % Make A with given eigenvalues
b = randn(n,1); % Make b
x1 = A\b;

%% the standard form of the conjugate gradient method, Numerical Optimizatin 2e, P112
ErrTol = 1.e-9; % ErrTol is the desired accuracy
MaxIts = 100; % MaxIts is the maximum number of CG steps
x0 = zeros(length(b), 1); % x0 is the initial guess, the default value is x0 = 0

% Initialization
k = 0; x = x0; % Initialize x
% r is the gradient of objective function J = 0.5*x'*Ax - b'*x
r = A*x - b; % Initialize residuals 
p = -r;  % Initialize conjugate direction
deltaNew = r'*r; % deltaNew is the square error
delta0 = deltaNew; % initial square error 
% Main iteration
while (k < MaxIts) && (deltaNew > ErrTol^2*delta0)
q = A*p;
alpha = deltaNew/(p'*q); % update step size Eq.(5.24a)
x = x + alpha*p; % Update x Eq.(5.24b)
r = r + alpha*q; % update residuals Eq.(5.24c)
deltaOld = deltaNew; % save previous square error
deltaNew = r'*r; % Update square error Eq.(5.24e)
p = -r + (deltaNew/deltaOld)*p; % Update conjugate direction Eq.(5.24f)
k = k + 1; % Increment iteration count
end

% Check accuracy
norm(x1 - x)
norm(x1 - x)/norm(x1) 



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

历史上的今天

在LOFTER的更多文章

评论

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

页脚

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