tests/steepest_descent_tests.m

   1 ## We can use the steepest descent method to solve Qx=b as in the
   2 ## conjugate gradient method. Like we did there, we test that the
   3 ## steepest descent method agrees with Octave's PCGM and our CGM.
   4 ##
   5 ## Note: The Steepest descent method uses the infinity norm as a
   6 ## stopping condition, so we should too.
   7 ##
   8 max_iterations = 100000;
   9 tolerance = 1e-8;
  10
  11 ## First a simple example.
  12 Q = [5,1,2; ...
  13      1,6,3; ...
  14      2,3,7];
  15
  16 b = [1;2;3];
  17 x0 = [1;1;1];
  18
  19 cgm  = conjugate_gradient_method(Q, b, x0, tolerance, max_iterations);
  20
  21
  22 q = @(x) (1/2)*x'*Q*x - b'*x;
  23 g = @(x) Q*x - b; % The gradient of q at x.
  24
  25 % The step size algorithm to use in the steepest descent method.
  26 step_size = @(x) step_length_positive_definite(g(x), Q);
  27 sd = steepest_descent(g, x0, step_size, tolerance, max_iterations);
  28
  29 diff = norm(cgm - sd, 'inf');
  30 unit_test_equals("Steepest descent agrees with CGM", ...
  31                  true, ...
  32                  diff < 2*tolerance);
  33
  34
  35 ## Test again Octave's pcg() function.
  36 for n = [ 5, 10, 25, 50, 100 ]
  37   Q = random_positive_definite_matrix(n, 100);
  38   C = random_positive_definite_matrix(n, 100);
  39
  40   ## Assumed by Octave's implementation when you don't supply a
  41   ## preconditioner.
  42   x0 = zeros(n, 1);
  43   b  = unifrnd(-100, 100, n, 1);
  44
  45   q = @(x) (1/2)*x'*Q*x - b'*x;
  46   g = @(x) Q*x - b; % The gradient of q at x.
  47
  48   % The step size algorithm to use in the steepest descent method.
  49   step_size = @(x) step_length_positive_definite(g(x), Q);
  50
  51   ## pcg() stops when the /relative/ norm falls below tolerance. To
  52   ## eliminate the relativity, we divide the tolerance by the
  53   ## quantity that pcg() will divide by.
  54   [x_pcg, o_flag, o_relres, o_iter] = pcg(Q, ...
  55                                           b, ...
  56                                           tolerance/norm(g(x0)), ...
  57                                           max_iterations, ...
  58                                           C, ...
  59                                           C');
  60   [x_sd, k] = steepest_descent(g, x0, step_size, tolerance, max_iterations);
  61
  62   diff = norm(x_pcg - x_sd, 'inf');
  63   msg = sprintf("Our steepest descent agrees with Octave's pcg, n=%d.", n);
  64
  65   ## There's no good way to choose the tolerance here, since each
  66   ## individual algorithm terminates based on the (2,infinity)-norm of
  67   ## the gradient.
  68   unit_test_equals(msg, true, diff <= sqrt(tolerance));
  69 end