Fix the n != 5 cases in the steepest descent tests.

Nerf the steepest descent test cases.
Take into account that pcg() uses the 2-norm in the steepest descent tests.

diff --git a/tests/steepest_descent_tests.m b/tests/steepest_descent_tests.m
index 7b9982d7d72510d47864ead036234d5abeb3de06..035c82d20bfafbc56021f1b2a11b7ffc1839dee2 100644 (file)
--- a/tests/steepest_descent_tests.m
+++ b/tests/steepest_descent_tests.m
@@ -6,7 +6,7 @@
 ## stopping condition, so we should too.
 ##
 max_iterations = 100000;
 ## stopping condition, so we should too.
 ##
 max_iterations = 100000;

-tolerance = 1e-11;
+tolerance = 1e-10;

 
 ## First a simple example.
 Q = [5,1,2; ...
 
 ## First a simple example.
 Q = [5,1,2; ...


@@ -23,7 +23,7 @@ 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.
 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));
+step_size = @(x) step_length_positive_definite(g(x), Q);

 sd = steepest_descent(g, x0, step_size, tolerance, max_iterations);
 
 diff = norm(cgm - sd, 'inf');
 sd = steepest_descent(g, x0, step_size, tolerance, max_iterations);
 
 diff = norm(cgm - sd, 'inf');


@@ -34,29 +34,33 @@ unit_test_equals("Steepest descent agrees with CGM", ...
 
 ## Test again Octave's pcg() function.
 for n = [ 5, 10, 25, 50, 100 ]
 
 ## 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);
+  Q = random_positive_definite_matrix(n, 100);
+  C = random_positive_definite_matrix(n, 100);

 
   ## Assumed by Octave's implementation when you don't supply a
   ## preconditioner.
 
   ## Assumed by Octave's implementation when you don't supply a
   ## preconditioner.

-  x0 = zeros(5, 1);
-  b  = unifrnd(-1000, 1000, 5, 1);
+  x0 = zeros(n, 1);
+  b  = unifrnd(-100, 100, n, 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.
 
   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));
+  step_size = @(x) step_length_positive_definite(g(x), Q);

 
 

+  ## pcg() stops when the /relative/ norm falls below tolerance. To
+  ## eliminate the relativity, we divide the tolerance by the
+  ## quantity that pcg() will divide by.

   [x_pcg, o_flag, o_relres, o_iter] = pcg(Q, ...
                                          b, ...
   [x_pcg, o_flag, o_relres, o_iter] = pcg(Q, ...
                                          b, ...

-                                         tolerance, ...
+                                         tolerance/norm(g(x0)), ...

                                          max_iterations, ...
                                          C, ...
                                          C');
   x_sd = steepest_descent(g, x0, step_size, tolerance, max_iterations);
 
                                          max_iterations, ...
                                          C, ...
                                          C');
   x_sd = steepest_descent(g, x0, step_size, tolerance, max_iterations);
 

-  diff = norm(x_pcg - x_sd, 'inf');
+  ## Note: pcg() uses the 2-norm.
+  diff = abs(norm(g(x_pcg)) - norm(g(x_sd), 'inf'));

   msg = sprintf("Our steepest descent agrees with Octave's pcg, n=%d.", n);
   msg = sprintf("Our steepest descent agrees with Octave's pcg, n=%d.", n);

-  unit_test_equals(msg, true, diff < 2*tolerance);
+  unit_test_equals(msg, true, diff <= tolerance);

 end
 end