-function [x, k] = steepest_descent(g, x0, step_size, tolerance, max_iterations)
+function [x, k] = steepest_descent(g, ...
+ x0, ...
+ step_size, ...
+ tolerance, ...
+ max_iterations)
%
% An implementation of the steepest-descent algorithm, with the
% search direction d_{k} = -\nabla f(x_{k}).
%
% We should terminate when either,
%
- % a) The 2-norm of the gradient at x_{k} is greater than
- % ``tolerance``.
+ % a) The infinity-norm of the gradient at x_{k} is greater than
+ % ``tolerance``.
%
% b) We have performed ``max_iterations`` iterations.
%
% receive.
%
% * ``tolerance`` - the stopping tolerance. We stop when the norm
- % of the gradient falls below this value.
+ % of the gradient falls below this value.
%
% * ``max_iterations`` - a safety net; we return x_{k}
- % unconditionally if we exceed this number of iterations.
+ % unconditionally if we exceed this number of iterations.
%
% OUTPUT:
%
% * ``x`` - the solution at the final value of x_{k}.
%
% * ``k`` - the value of k when we stop; i.e. the number of
- % iterations.
+ % iterations.
+ %
+ % NOTES:
+ %
+ % A specialized implementation for solving e.g. Qx=b can avoid one
+ % matrix-vector multiplication below.
+ %
% The initial gradient at x_{0} is not supplied, so we compute it
- % here and begin the loop at k=1.
- x = x0;
- g_k = g(x);
+ % here and begin the loop at k=0.
+ k = 0;
+ xk = x0;
+ gk = g(xk);
- if (norm(g_k) < tolerance)
- % If x_0 is close enough to a solution, there's nothing for us to
- % do! We use g_k (the gradient of f at x_k) instead of d_k because
- % their 2-norms will be the same, and g_k is already stored.
- return;
- end
+ while (k <= max_iterations && norm(gk, 'inf') > tolerance)
+ % Loop until either of our stopping conditions are met. This
+ % should also work when x0 is a satisfactory point.
- for k = [1 : max_iterations]
- % Loop until either of our stopping conditions are met. If the
- % loop finishes, we have implicitly met the second stopping
- % condition (number of iterations).
- d_k = -g_k;
- alpha_k = step_size(x);
- x = x + (alpha_k * d_k);
- g_k = g(x);
+ dk = -gk;
+ alpha_k = step_size(xk);
+ x_next = xk + (alpha_k * dk);
- if (norm(g_k) < tolerance)
- return;
- end
+ % We potentially just performed one more iteration than necessary
+ % in order to simplify the loop. Note that due to the structure of
+ % our loop, we will have k > max_iterations when we fail to
+ % converge.
+ k = k + 1;
+ xk = x_next;
+ gk = g(x_next);
end
- % If we make it to the end of the loop, that means we've executed the
- % maximum allowed iterations. The caller should be able to examine the
- % return value ``k`` to determine what happened.
+ x = xk;
end