From f51328e59405aa5e02509b500da6c455a10ceb54 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Mon, 25 Mar 2013 12:50:44 -0400 Subject: [PATCH] Add tests for the steepest descent method, based on the PCGM tests. --- tests/steepest_descent_tests.m | 62 ++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) create mode 100644 tests/steepest_descent_tests.m diff --git a/tests/steepest_descent_tests.m b/tests/steepest_descent_tests.m new file mode 100644 index 0000000..7b9982d --- /dev/null +++ b/tests/steepest_descent_tests.m @@ -0,0 +1,62 @@ +## 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-11; + +## 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, -g(x)); +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(5, 1000); + C = random_positive_definite_matrix(5, 1000); + + ## Assumed by Octave's implementation when you don't supply a + ## preconditioner. + x0 = zeros(5, 1); + b = unifrnd(-1000, 1000, 5, 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, -g(x)); + + [x_pcg, o_flag, o_relres, o_iter] = pcg(Q, ... + b, ... + tolerance, ... + max_iterations, ... + C, ... + C'); + x_sd = 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); + unit_test_equals(msg, true, diff < 2*tolerance); +end -- 2.33.1