From 6c0b2e3b90fad51bbbb14f4755127a7920fd59c4 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Tue, 19 Mar 2013 09:28:17 -0400 Subject: [PATCH] Allow CGM to take a maximum number of iterations, return the amount needed. --- optimization/conjugate_gradient_method.m | 29 ++++++++++++++++-------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/optimization/conjugate_gradient_method.m b/optimization/conjugate_gradient_method.m index 9b2ddac..2c94487 100644 --- a/optimization/conjugate_gradient_method.m +++ b/optimization/conjugate_gradient_method.m @@ -1,4 +1,4 @@ -function x_star = conjugate_gradient_method(A, b, x0, tolerance) +function [x, k] = conjugate_gradient_method(A, b, x0, tolerance, max_iterations) % % Solve, % @@ -8,7 +8,8 @@ function x_star = conjugate_gradient_method(A, b, x0, tolerance) % % min [phi(x) = (1/2)* + ] % - % using Algorithm 5.2 in Nocedal and Wright. + % using the conjugate_gradient_method (Algorithm 5.2 in Nocedal and + % Wright). % % INPUT: % @@ -22,9 +23,14 @@ function x_star = conjugate_gradient_method(A, b, x0, tolerance) % - ``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_star`` - The solution to Ax=b. + % - ``x`` - The solution to Ax=b. + % + % - ``k`` - The ending value of k; that is, the number of iterations that + % were performed. % % NOTES: % @@ -33,22 +39,25 @@ function x_star = conjugate_gradient_method(A, b, x0, tolerance) 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 -- 2.43.2