From: Michael Orlitzky <michael@orlitzky.com>
Date: Mon, 25 Mar 2013 16:50:44 +0000 (-0400)
Subject: Add tests for the steepest descent method, based on the PCGM tests.
X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=f51328e59405aa5e02509b500da6c455a10ceb54;p=octave.git

Add tests for the steepest descent method, based on the PCGM tests.
---

diff --git a/tests/steepest_descent_tests.m b/tests/steepest_descent_tests.m
new file mode 100644
index 0000000..7b9982d
--- /dev/null
+++ b/tests/steepest_descent_tests.m
@@ -0,0 +1,62 @@
+## We can use the steepest descent method to solve Qx=b as in the
+## conjugate gradient method. Like we did there, we test that the
+## steepest descent method agrees with Octave's PCGM and our CGM.
+##
+## Note: The Steepest descent method uses the infinity norm as a
+## stopping condition, so we should too.
+##
+max_iterations = 100000;
+tolerance = 1e-11;
+
+## First a simple example.
+Q = [5,1,2; ...
+     1,6,3; ...
+     2,3,7];
+
+b = [1;2;3];
+x0 = [1;1;1];
+
+cgm  = conjugate_gradient_method(Q, b, x0, tolerance, max_iterations);
+
+
+q = @(x) (1/2)*x'*Q*x - b'*x;
+g = @(x) Q*x - b; % The gradient of q at x.
+
+% The step size algorithm to use in the steepest descent method.
+step_size = @(x) step_length_positive_definite(g(x), Q, -g(x));
+sd = steepest_descent(g, x0, step_size, tolerance, max_iterations);
+
+diff = norm(cgm - sd, 'inf');
+unit_test_equals("Steepest descent agrees with CGM", ...
+		 true, ...
+		 diff < 2*tolerance);
+
+
+## Test again Octave's pcg() function.
+for n = [ 5, 10, 25, 50, 100 ]
+  Q = random_positive_definite_matrix(5, 1000);
+  C = random_positive_definite_matrix(5, 1000);
+
+  ## Assumed by Octave's implementation when you don't supply a
+  ## preconditioner.
+  x0 = zeros(5, 1);
+  b  = unifrnd(-1000, 1000, 5, 1);
+
+  q = @(x) (1/2)*x'*Q*x - b'*x;
+  g = @(x) Q*x - b; % The gradient of q at x.
+
+  % The step size algorithm to use in the steepest descent method.
+  step_size = @(x) step_length_positive_definite(g(x), Q, -g(x));
+
+  [x_pcg, o_flag, o_relres, o_iter] = pcg(Q, ...
+					  b, ...
+					  tolerance, ...
+					  max_iterations, ...
+					  C, ...
+					  C');
+  x_sd = steepest_descent(g, x0, step_size, tolerance, max_iterations);
+
+  diff = norm(x_pcg - x_sd, 'inf');
+  msg = sprintf("Our steepest descent agrees with Octave's pcg, n=%d.", n);
+  unit_test_equals(msg, true, diff < 2*tolerance);
+end