X-Git-Url: http://gitweb.michael.orlitzky.com/?p=octave.git;a=blobdiff_plain;f=optimization%2Fsteepest_descent.m;h=0a6475e4614c37822be1c36ce891c003aae6ca55;hp=32aad711457ec2d04aa5d0985cd4286769bb79cc;hb=8ff5cc31794c0f0d2e7c4c5e8c9cb9552c4719d2;hpb=b1aa6fb5c819eb0ba94dd9dd65e7711202cddeac

diff --git a/optimization/steepest_descent.m b/optimization/steepest_descent.m
index 32aad71..0a6475e 100644
--- a/optimization/steepest_descent.m
+++ b/optimization/steepest_descent.m
@@ -1,12 +1,16 @@
-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.
   %
@@ -24,45 +28,45 @@ function [x, k] = steepest_descent(g, x0, step_size, tolerance, max_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);
+    xk = xk + (alpha_k * dk);
+    gk = g(xk);
 
-    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;
   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