From b1aa6fb5c819eb0ba94dd9dd65e7711202cddeac Mon Sep 17 00:00:00 2001
From: Michael Orlitzky
Date: Sun, 17 Mar 2013 22:21:31 -0400
Subject: [PATCH] Add first implementation of the steepest descent method.
---
optimization/steepest_descent.m | 68 +++++++++++++++++++++++++++++++++
1 file changed, 68 insertions(+)
create mode 100644 optimization/steepest_descent.m
diff --git a/optimization/steepest_descent.m b/optimization/steepest_descent.m
new file mode 100644
index 0000000..32aad71
--- /dev/null
+++ b/optimization/steepest_descent.m
@@ -0,0 +1,68 @@
+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``.
+ %
+ % b) We have performed ``max_iterations`` iterations.
+ %
+ % INPUT:
+ %
+ % * ``g`` - the gradient of the function ``f`` we're minimizing.
+ %
+ % * ``x0`` - the starting point for the search.
+ %
+ % * ``step_size`` - a function of x which returns the optimal step
+ % size alpha_{k}. This will generally require more information
+ % than just x; for example it might need the function ``f`` or the
+ % gradient ``g``. However, our caller has this information so it
+ % should be incorporated into the step size algorithm that we
+ % receive.
+ %
+ % * ``tolerance`` - the stopping tolerance. We stop when the norm
+ % of the gradient falls below this value.
+ %
+ % * ``max_iterations`` - a safety net; we return x_{k}
+ % 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.
+
+ % 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);
+
+ 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
+
+ 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);
+
+ if (norm(g_k) < tolerance)
+ return;
+ end
+ 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.
+end
--
2.33.1