## stopping condition, so we should too.
##
max_iterations = 100000;
-tolerance = 1e-11;
+tolerance = 1e-8;
## First a simple example.
Q = [5,1,2; ...
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));
+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');
## 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);
+ 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(5, 1);
- b = unifrnd(-1000, 1000, 5, 1);
+ 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, -g(x));
+ 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, ...
+ tolerance/norm(g(x0)), ...
max_iterations, ...
C, ...
C');
- x_sd = steepest_descent(g, x0, step_size, tolerance, max_iterations);
+ [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);
- unit_test_equals(msg, true, diff < 2*tolerance);
+
+ ## 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