1% Solve a linear system with the preconditioned conjugated-gradient method.
2%
3% [x] = conjugategradient(fun,b,tol,maxit,pc,pc2,x0);
4% [x,Q,T] = conjugategradient(fun,b,tol,maxit,pc,pc2,x0);
5%
6% solves for x
7% A x = b
8% using the preconditioned conjugated-gradient method.
9% It provides also an approximation of A:
10% A \sim Q*T*Q'
11%
12% J(x) = 1/2 (x' b - x' A x)
13% ∇ J = b - A x
14% A x = b - ∇ J
15% b = ∇ J(0)
16%
17% the columns of Q are the Lanczos vectors
18
19% Alexander Barth 2010,2012-2013
20
21function [x,Q,T,diag] = conjugategradient(fun,b,varargin)
22
23n = length(b);
24
25% default parameters
26maxit = min(n,20);
27minit = 0;
28tol = 1e-6;
29pc = [];
30x0 = [];
31renorm = 0;
32
33prop = varargin;
34
35for i=1:2:length(prop)
36    if strcmp(prop{i},'minit')
37        minit = prop{i+1};
38    elseif strcmp(prop{i},'maxit')
39        maxit = prop{i+1};
40    elseif strcmp(prop{i},'tol')
41        tol = prop{i+1};
42    elseif strcmp(prop{i},'x0')
43        x0 = prop{i+1};
44    elseif strcmp(prop{i},'pc')
45        pc = prop{i+1};
46    elseif strcmp(prop{i},'renorm')
47        renorm = prop{i+1};
48    else
49
50    end
51end
52
53
54if isempty(x0)
55    % random initial vector
56    x0 = randn(n,1);
57end
58
59
60if isempty(pc)
61    % no pre-conditioning
62    pc = @(x) x;
63elseif isnumeric(pc)
64    invM = pc;
65    pc = @(x)( invM *x);
66end
67
68
69tol2 = tol^2;
70
71delta = [];
72gamma = [];
73q = NaN*zeros(n,1);
74
75%M = inv(invM);
76%E = chol(M);
77
78
79% initial guess
80x = x0;
81
82% gradient at initial guess
83r = b - fun(x);
84
85% quick exit
86if r'*r < tol2
87    return
88end
89
90% apply preconditioner
91z = pc(r);
92
93% first search direction == gradient
94p = z;
95
96% compute: r' * inv(M) * z (we will need this product at several
97% occasions)
98
99zr_old = r'*z;
100
101% r_old: residual at previous iteration
102r_old = r;
103
104for k=1:maxit
105    if k <= n && nargout > 1
106        % keep at most n vectors
107        Q(:,k) = r/sqrt(zr_old);
108    end
109
110    % compute A*p
111    Ap = fun(p);
112    %maxdiff(A*p,Ap)
113
114    % how far do we need to go in direction p?
115    % alpha is determined by linesearch
116
117    % alpha z'*r / (p' * A * p)
118    alpha(k) = zr_old / ( p' * Ap);
119
120    % get new estimate of x
121    x = x + alpha(k)*p;
122
123    % recompute gradient at new x. Could be done by
124    % r = b-fun(x);
125    % but this does require an new call to fun
126    r = r - alpha(k)*Ap;
127
128    if renorm
129        r = r - Q(:,1:k) * Q(:,1:k)' * r ;
130    end
131    %maxdiff(r,b-fun(x))
132
133    % apply pre-conditionner
134    z = pc(r);
135
136
137    zr_new = r'*z;
138
139    if r'*r < tol2 && k >= minit
140        break
141    end
142
143    %Fletcher-Reeves
144    beta(k+1) = zr_new / zr_old;
145    %Polak-Ribiere
146    %beta(k+1) = r'*(r-r_old) / zr_old;
147    %Hestenes-Stiefel
148    %beta(k+1) = r'*(r-r_old) / (p'*(r-r_old));
149    %beta(k+1) = r'*(r-r_old) / (r_old'*r_old);
150
151
152    % norm(p)
153    p = z + beta(k+1)*p;
154    zr_old = zr_new;
155    r_old = r;
156end
157
158
159%disp('alpha and beta')
160%figure,plot(beta(2:end))
161%rg(alpha)
162%rg(beta(2:end))
163
164if nargout > 1
165    kmax = size(Q,2);
166
167    delta(1) = 1/alpha(1);
168
169    %delta(1) - Q(:,1)'*invM*A*invM*Q(:,1)
170
171
172    for k=1:kmax-1
173        delta(k+1) = 1/alpha(k+1) + beta(k+1)/alpha(k);
174        gamma(k) = -sqrt(beta(k+1))/alpha(k);
175    end
176
177    T = sparse([1:kmax   1:kmax-1 2:kmax  ],...
178        [1:kmax   2:kmax   1:kmax-1],...
179        [delta    gamma    gamma]);
180
181    diag.iter = k;
182    diag.relres = sqrt(r'*r);
183end
184
185% Copyright (C) 2004 Alexander Barth <a.barth@ulg.ac.be>
186%
187% This program is free software; you can redistribute it and/or modify it under
188% the terms of the GNU General Public License as published by the Free Software
189% Foundation; either version 2 of the License, or (at your option) any later
190% version.
191%
192% This program is distributed in the hope that it will be useful, but WITHOUT
193% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
194% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
195% details.
196%
197% You should have received a copy of the GNU General Public License along with
198% this program; if not, see <http://www.gnu.org/licenses/>.
199