## stopping condition, so we should too.
##
max_iterations = 100000;
-tolerance = 1e-11;
+tolerance = 1e-10;
## 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);
- diff = norm(x_pcg - x_sd, 'inf');
+ ## Note: pcg() uses the 2-norm.
+ diff = abs(norm(g(x_pcg)) - norm(g(x_sd), 'inf'));
msg = sprintf("Our steepest descent agrees with Octave's pcg, n=%d.", n);
- unit_test_equals(msg, true, diff < 2*tolerance);
+ unit_test_equals(msg, true, diff <= tolerance);
end