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