-
Notifications
You must be signed in to change notification settings - Fork 326
/
perform_conjugate_gradient.m
129 lines (112 loc) · 3.12 KB
/
perform_conjugate_gradient.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
function [x,err,it] = perform_conjugate_gradient(A,y,options)
% perform_conjugate_gradient - perform (bi)-conjugate gradient
%
% [x,err,k] = perform_conjugate_gradient(A,y,options);
%
% Solves for A*x=y.
% Works both for vector x,y and matrix x,y (parralel soving of many
% linear equations).
%
% Important: the algorithm assumes that the matrix is symmetric definite
% positive. If it is not the case, then you must set options.is_sdp=0,
% and the algorithm will use a bi-conjugate gradient descent (x2 slower).
%
% A can be a matrix or a callback function y=A(x,options).
% In this case, you have to set options.ncols as the number of columns of
% A (if it is a callback).
% In the non-symmetric case, then options.transpose=1 for the computation
% of A*x and options.transpose=-1 for the computation of A'*x.
%
% err monitors the decreasing of |A*x-y|.
% k is the total number of iterations.
%
% You can set:
% options.x is an initial guess
% options.epsilon is maximum error
% options.niter_max is maximum number of iterations
%
% Copyright (c) 2007 Gabriel Peyre
options.null = 0;
niter = getoptions(options, 'niter_max', 100);
epsilon = getoptions(options, 'epsilon', 1e-5);
is_sdp = getoptions(options, 'is_sdp', 1);
x = getoptions(options, 'x', []);
use_callback = 0;
if not(isnumeric(A))
use_callback = 1;
end
if not(isempty(x))
x = options.x;
else
if use_callback==0
x = zeros(size(A,2),1);
else
if isfield(options, 'ncols')
x = zeros(options.ncols, 1);
else
error('You have to specify options.ncols');
end
end
end
if norm(y(:), 'fro') == 0
normb = epsilon;
else
normb = epsilon; % * sum(y(:).^2);
end
if use_callback==0
r = y - A*x;
else
options.transpose = 1;
r = y - feval(A,x,options);
end
p = r;
if is_sdp==0
rr = r;
pp = p;
end
r0 = sum(r(:).^2);
err = [sum(r0)];
for it=1:niter
% auxiliary vector
if use_callback==0
w = A*p;
else
options.transpose = 1;
% options.transpose = -1;
w = feval(A,p,options);
end
if is_sdp==1
d = sum(p(:) .* w(:));
else
d = sum(pp(:) .* w(:));
end
I = find(abs(d)<eps); d(I) = 1;
alpha = repmat( r0 ./ d, [size(x,1) 1] ); % found optimal alpha in line search
x = x + alpha.*p; % new guess
r = r - alpha.*w; % the residual is in fact r=b-A*x
if is_sdp==0
if use_callback==0
w = A'*pp;
else
options.transpose = -1;
w = feval(A,pp,options);
end
rr = rr-alpha.*w;
end
rprev = r0; % save norm of the old residual
if is_sdp==0
r0 = sum(r.*rr);
else
r0 = sum(r.^2); % compute norm of new residual
end
err(end+1) = sqrt( sum(r.^2) );
if err(end)<normb
return;
end
% search direction
beta = r0./rprev;
p = r + repmat(beta, [size(x,1) 1]).*p;
if is_sdp==0
pp = rr + repmat(beta, [size(x,1) 1]).*pp;
end
end