Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
NCTM or The Math Forum.



Nonlinear Conjugate Gradient
Posted:
Mar 16, 2011 8:18 PM


I have been trying to implement the nonlinear conjugate gradient algorithm and have been having some difficulty and am not sure what is going wrong.
On wikipedia there is a discussion of the conjugate gradient algorithm and a numerical example. http://en.wikipedia.org/wiki/Conjugate_gradient_method
This algorithm works in MATLAB as I modified the code provided by XinShe Yang while using the the Polak–Ribière formula
A=[4,1;1,3]; b=[1;2];
% Solve the system by the conjugate gradient method %u=rand(n,1); % initial guess u0 u=[2;1]; r=bA*u; % intial r0 d=r; % intial d0=r0 delta=10^(5); % accuracy
while max(abs(r))>delta alpha=r'*r/(d'*A*d); % alpha_n K=r'*r; % temporary value u=u+alpha*d; % u_{n+1} last_r=r; r=ralpha*A*d; % r_{n+1} %beta=r'*r/K; % beta_n Fletcher Reeves beta=(r'*(rlast_r))/K; % beta_n Polak Ribiere d=r+beta*d; % d_{n+1} u' % disp u on the screen end
The nonlinear conjugate gradient method is what I am interested in since it allows a generic function to be minimized instead of Ax=b. Discussion on Wikipedia is here http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient
There is also discussion here on the modifications from going from linear to nonlinear conjugate gradient. http://www.ee.ucla.edu/~vandenbe/236C/lectures/cg.pdf
My attempt at the nonlinear conjugate gradient method is
A=[4,1;1,3]; b=[1;2];
x=[2;1];
F=A*xb;
g=gradient(F,x);
d=g;
for k=1:10 a=linesrch(d,x,b,A); x = x + a*d; last_g=g; K=g'*g; F=A*xb;
g=gradient(F,x); beta=norm((g'*(glast_g)))^2/norm(K)^2; % beta_n Polak Ribiere d=g+beta*d; x' % disp x on the screen
end

function a = linesrch(d,x,b,A)
a = 1e4; while myfun(a,d,x,b,A) < 0 a = 1.6*a; end a = fminunc(@(a) myfun(a,d,x,b,A),1);

function [f] = myfun(a,p,x,b,A)
f=gradient(A*(x+a*p)b).'*p;
Right now the results blow up. I am rying to replicate the results of the Ax=b example using nonlinear conjugate gradient so I can move on to other functions.
I see some recent related discussion here as well http://www.mathworks.com/matlabcentral/newsreader/view_thread/304151 which I also tried. Any suggestions would be welcome. Thanks.



