From: Michael Orlitzky <michael@orlitzky.com>
Date: Wed, 13 Mar 2013 18:05:21 +0000 (-0400)
Subject: Add first implementation of the conjugate gradient method and its tests.
X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=336f1181f7f065a07ee58bd772a875f7f4b39247;p=octave.git

Add first implementation of the conjugate gradient method and its tests.
---

diff --git a/optimization/conjugate_gradient_method.m b/optimization/conjugate_gradient_method.m
new file mode 100644
index 0000000..5b07e1c
--- /dev/null
+++ b/optimization/conjugate_gradient_method.m
@@ -0,0 +1,54 @@
+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);
+
+  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;
+end
diff --git a/tests/conjugate_gradient_method_tests.m b/tests/conjugate_gradient_method_tests.m
new file mode 100644
index 0000000..56987a5
--- /dev/null
+++ b/tests/conjugate_gradient_method_tests.m
@@ -0,0 +1,16 @@
+A = [5,1,2; ...
+     1,6,3;
+     2,3,7];
+
+b = [1;2;3];
+
+x0 = [1;1;1];
+
+## Solved over the rationals.
+expected = [2/73; 11/73; 26/73];
+actual = conjugate_gradient_method(A, b, x0, 1e-6);
+diff = norm(actual - expected);
+
+unit_test_equals("CGM works on an example", ...
+		 true, ...
+		 norm(diff) < 1e-6);