tests/preconditioned_conjugate_gradient_method_tests.m

   1 ## Used throughout. The PCGM uses the infinity norm as the stopping
   2 ## condition, so we had better also.
   3 max_iterations = 10000;
   4 tolerance = 1e-10;
   5
   6 ## First a simple example.
   7 A = [5,1,2; ...
   8      1,6,3; ...
   9      2,3,7];
  10
  11 M = eye(3);
  12 b = [1;2;3];
  13 x0 = [1;1;1];
  14
  15 cgm  = conjugate_gradient_method(A, b, x0, tolerance, max_iterations);
  16 pcgm = preconditioned_conjugate_gradient_method(A, ...
  17                                                 M, ...
  18                                                 b, ...
  19                                                 x0, ...
  20                                                 tolerance, ...
  21                                                 max_iterations);
  22 diff = norm(cgm - pcgm, 'inf');
  23
  24 unit_test_equals("PCGM agrees with CGM when M == I", ...
  25                  true, ...
  26                  diff < 2*tolerance);
  27
  28 pcgm_simple = simple_preconditioned_cgm(A, M, b, x0, tolerance, max_iterations);
  29 diff = norm(pcgm_simple - pcgm, 'inf');
  30
  31 unit_test_equals("PCGM agrees with SimplePCGM when M == I", ...
  32                  true, ...
  33                  diff < 2*tolerance);
  34
  35 ## Needs to be symmetric!
  36 M = [0.97466, 0.24345, 0.54850; ...
  37      0.24345, 0.73251, 0.76639; ...
  38      0.54850, 0.76639, 1.47581];
  39
  40 pcgm = preconditioned_conjugate_gradient_method(A, ...
  41                                                 M, ...
  42                                                 b, ...
  43                                                 x0, ...
  44                                                 tolerance, ...
  45                                                 max_iterations);
  46 diff = norm(cgm - pcgm, 'inf');
  47
  48 unit_test_equals("PCGM agrees with CGM when M != I", ...
  49                  true, ...
  50                  diff < 2*tolerance);
  51
  52
  53 pcgm_simple = simple_preconditioned_cgm(A, M, b, x0, tolerance, max_iterations);
  54 diff = norm(pcgm_simple - pcgm, 'inf');
  55
  56 unit_test_equals("PCGM agrees with Simple PCGM when M != I", ...
  57                  true, ...
  58                  diff < 2*tolerance);
  59
  60
  61 # Test again Octave's pcg() function.
  62 for n = [ 5, 10, 25, 50, 100 ]
  63   A = random_positive_definite_matrix(n, 100);
  64
  65   # Use the cholesky factorization as a preconditioner.
  66   Ct = perturb(chol(A));
  67   C = Ct';
  68   M = Ct*C;
  69
  70   # Assumed by Octave's implementation when you don't supply a
  71   # preconditioner.
  72   x0 = zeros(n, 1);
  73   b  = unifrnd(-100, 100, n, 1);
  74   g = @(x) A*x - b;
  75
  76   ## pcg() stops when the /relative/ norm falls below tolerance. To
  77   ## eliminate the relativity, we divide the tolerance by the
  78   ## quantity that pcg() will divide by.
  79   [o_x, o_flag, o_relres, o_iter] = pcg(A, b, tolerance/norm(g(x0)), ...
  80                                         max_iterations, C, C');
  81   [x, k] = preconditioned_conjugate_gradient_method(A,
  82                                                     M,
  83                                                     b,
  84                                                     x0,
  85                                                     tolerance,
  86                                                     max_iterations);
  87
  88   diff = norm(o_x - x, 'inf');
  89   msg = sprintf("Our PCGM agrees with Octave's, n=%d.", n);
  90   ## There's no good way to choose the tolerance here, since each
  91   ## individual algorithm terminates based on the (2,infinity)-norm of
  92   ## the gradient. So we use two orders of magnitude.
  93   unit_test_equals(msg, true, diff <= sqrt(tolerance));
  94 end