Move conjugate_gradient_method.m to vanilla_cgm.m.

[octave.git] / optimization / conjugate_gradient_method.m

diff --git a/optimization/conjugate_gradient_method.m b/optimization/conjugate_gradient_method.m
index 5b07e1c183f88f5f0764975a23edc1c27dad446e..a6401e5ec21ff963883368576deb59ff0e1b36d0 100644 (file)
--- a/optimization/conjugate_gradient_method.m
+++ b/optimization/conjugate_gradient_method.m
@@ -1,54 +1,49 @@
-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.
-  ##
-  zero_vector = zeros(length(x0), 1);
+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.
+  %
+  % 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.
+  %
+  n = length(x0);
+  M = eye(n);
 
-  k = 0;
-  xk = x0;
-  rk = A*xk - b; # The first residual must be computed the hard way.
-  pk = -rk;
-
-  while (norm(rk) > tolerance)
-    alpha_k = step_length_cgm(rk, A, pk);
-    x_next = xk + 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;
-    rk = r_next;
-    pk = p_next;
-  end
-
-  x_star = xk;
+  % The standard CGM is equivalent to the preconditioned CGM is you
+  % use the identity matrix as your preconditioner.
+  [x, k] = preconditioned_conjugate_gradient_method(A,
+                                                   M,
+                                                   b,
+                                                   x0,
+                                                   tolerance,
+                                                   max_iterations);
 end