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;

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);

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(n, 100);

38 C = random_positive_definite_matrix(n, 100);

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);

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);

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);

62 diff = norm(x_pcg - x_sd, 'inf');

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

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