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-11;
11 ## First a simple example.
12 Q = [5,1,2; ...
13 1,6,3; ...
14 2,3,7];
16 b = [1;2;3];
17 x0 = [1;1;1];
19 cgm = conjugate_gradient_method(Q, b, x0, tolerance, max_iterations);
22 q = @(x) (1/2)*x'*Q*x - b'*x;
23 g = @(x) Q*x - b; % The gradient of q at x.
25 % The step size algorithm to use in the steepest descent method.
26 step_size = @(x) step_length_positive_definite(g(x), Q, -g(x));
27 sd = steepest_descent(g, x0, step_size, tolerance, max_iterations);
29 diff = norm(cgm - sd, 'inf');
30 unit_test_equals("Steepest descent agrees with CGM", ...
31 true, ...
32 diff < 2*tolerance);
35 ## Test again Octave's pcg() function.
36 for n = [ 5, 10, 25, 50, 100 ]
37 Q = random_positive_definite_matrix(5, 1000);
38 C = random_positive_definite_matrix(5, 1000);
40 ## Assumed by Octave's implementation when you don't supply a
41 ## preconditioner.
42 x0 = zeros(5, 1);
43 b = unifrnd(-1000, 1000, 5, 1);
45 q = @(x) (1/2)*x'*Q*x - b'*x;
46 g = @(x) Q*x - b; % The gradient of q at x.
48 % The step size algorithm to use in the steepest descent method.
49 step_size = @(x) step_length_positive_definite(g(x), Q, -g(x));
51 [x_pcg, o_flag, o_relres, o_iter] = pcg(Q, ...
52 b, ...
53 tolerance, ...
54 max_iterations, ...
55 C, ...
56 C');
57 x_sd = steepest_descent(g, x0, step_size, tolerance, max_iterations);
59 diff = norm(x_pcg - x_sd, 'inf');
60 msg = sprintf("Our steepest descent agrees with Octave's pcg, n=%d.", n);
61 unit_test_equals(msg, true, diff < 2*tolerance);
62 end