-function x_star = conjugate_gradient_method(A, b, x0, tolerance)
- ##
- ## Solve,
- ##
- ## Ax = b
- ##
- ## or equivalently,
- ##
- ## min [phi(x) = (1/2)*<Ax,x> + <b,x>]
- ##
- ## using Algorithm 5.2 in Nocedal and Wright.
- ##
- ## INPUT:
- ##
- ## - ``A`` -- The coefficient matrix of the system to solve. Must
- ## be positive definite.
- ##
- ## - ``b`` -- The right-hand-side of the system to solve.
- ##
- ## - ``x0`` -- The starting point for the search.
- ##
- ## - ``tolerance`` -- How close ``Ax`` has to be to ``b`` (in
- ## magnitude) before we stop.
- ##
- ## OUTPUT:
- ##
- ## - ``x_star`` - The solution to Ax=b.
- ##
- ## NOTES:
- ##
- ## All vectors are assumed to be *column* vectors.
- ##
+function [x, k] = conjugate_gradient_method(A, b, x0, tolerance, max_iterations)
+ %
+ % Solve,
+ %
+ % Ax = b
+ %
+ % or equivalently,
+ %
+ % min [phi(x) = (1/2)*<Ax,x> + <b,x>]
+ %
+ % using the conjugate_gradient_method (Algorithm 5.2 in Nocedal and
+ % Wright).
+ %
+ % INPUT:
+ %
+ % - ``A`` -- The coefficient matrix of the system to solve. Must
+ % be positive definite.
+ %
+ % - ``b`` -- The right-hand-side of the system to solve.
+ %
+ % - ``x0`` -- The starting point for the search.
+ %
+ % - ``tolerance`` -- How close ``Ax`` has to be to ``b`` (in
+ % magnitude) before we stop.
+ %
+ % - ``max_iterations`` -- The maximum number of iterations to perform.
+ %
+ % OUTPUT:
+ %
+ % - ``x`` - The solution to Ax=b.
+ %
+ % - ``k`` - The ending value of k; that is, the number of iterations that
+ % were performed.
+ %
+ % NOTES:
+ %
+ % All vectors are assumed to be *column* vectors.
+ %
zero_vector = zeros(length(x0), 1);
k = 0;
- xk = x0;
- rk = A*xk - b; # The first residual must be computed the hard way.
+ x = x0; % Eschew the 'k' suffix on 'x' for simplicity.
+ rk = A*x - b; % The first residual must be computed the hard way.
pk = -rk;
- while (norm(rk) > tolerance)
+ for k = [ 0 : max_iterations ]
+ if (norm(rk) < tolerance)
+ % Success.
+ return;
+ end
+
alpha_k = step_length_cgm(rk, A, pk);
- x_next = xk + alpha_k*pk;
+ x_next = x + alpha_k*pk;
r_next = rk + alpha_k*A*pk;
beta_next = (r_next' * r_next)/(rk' * rk);
p_next = -r_next + beta_next*pk;
k = k + 1;
- xk = x_next;
+ x = x_next;
rk = r_next;
pk = p_next;
end
-
- x_star = xk;
end