## We can use the steepest descent method to solve Qx=b as in the ## conjugate gradient method. Like we did there, we test that the ## steepest descent method agrees with Octave's PCGM and our CGM. ## ## Note: The Steepest descent method uses the infinity norm as a ## stopping condition, so we should too. ## max_iterations = 100000; tolerance = 1e-8; ## First a simple example. Q = [5,1,2; ... 1,6,3; ... 2,3,7]; b = [1;2;3]; x0 = [1;1;1]; cgm = conjugate_gradient_method(Q, b, x0, tolerance, max_iterations); 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); sd = steepest_descent(g, x0, step_size, tolerance, max_iterations); diff = norm(cgm - sd, 'inf'); unit_test_equals("Steepest descent agrees with CGM", ... true, ... diff < 2*tolerance); ## Test again Octave's pcg() function. for n = [ 5, 10, 25, 50, 100 ] 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. 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. 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, ... tolerance/norm(g(x0)), ... max_iterations, ... C, ... C'); [x_sd, k] = steepest_descent(g, x0, step_size, tolerance, max_iterations); diff = norm(x_pcg - x_sd, 'inf'); msg = sprintf("Our steepest descent agrees with Octave's pcg, n=%d.", n); ## There's no good way to choose the tolerance here, since each ## individual algorithm terminates based on the (2,infinity)-norm of ## the gradient. unit_test_equals(msg, true, diff <= sqrt(tolerance)); end