Allow CGM to take a maximum number of iterations, return the amount needed.
[octave.git] / optimization / conjugate_gradient_method.m
1 function [x, k] = conjugate_gradient_method(A, b, x0, tolerance, max_iterations)
2 %
3 % Solve,
4 %
5 % Ax = b
6 %
7 % or equivalently,
8 %
9 % min [phi(x) = (1/2)*<Ax,x> + <b,x>]
10 %
11 % using the conjugate_gradient_method (Algorithm 5.2 in Nocedal and
12 % Wright).
13 %
14 % INPUT:
15 %
16 % - ``A`` -- The coefficient matrix of the system to solve. Must
17 % be positive definite.
18 %
19 % - ``b`` -- The right-hand-side of the system to solve.
20 %
21 % - ``x0`` -- The starting point for the search.
22 %
23 % - ``tolerance`` -- How close ``Ax`` has to be to ``b`` (in
24 % magnitude) before we stop.
25 %
26 % - ``max_iterations`` -- The maximum number of iterations to perform.
27 %
28 % OUTPUT:
29 %
30 % - ``x`` - The solution to Ax=b.
31 %
32 % - ``k`` - The ending value of k; that is, the number of iterations that
33 % were performed.
34 %
35 % NOTES:
36 %
37 % All vectors are assumed to be *column* vectors.
38 %
39 zero_vector = zeros(length(x0), 1);
40
41 k = 0;
42 x = x0; % Eschew the 'k' suffix on 'x' for simplicity.
43 rk = A*x - b; % The first residual must be computed the hard way.
44 pk = -rk;
45
46 for k = [ 0 : max_iterations ]
47 if (norm(rk) < tolerance)
48 % Success.
49 return;
50 end
51
52 alpha_k = step_length_cgm(rk, A, pk);
53 x_next = x + alpha_k*pk;
54 r_next = rk + alpha_k*A*pk;
55 beta_next = (r_next' * r_next)/(rk' * rk);
56 p_next = -r_next + beta_next*pk;
57
58 k = k + 1;
59 x = x_next;
60 rk = r_next;
61 pk = p_next;
62 end
63 end