In an effort to create an interactive gradient descent visualization I am using the following gradient descent function,

stepGradient[func_, start_, l_, steps_, vars_] := 
  stepGradient[func, start, l, steps, vars] = Module[{},
    vals = {start};
    grd[f_, v_][n_] := Grad[f, v] /. Thread[v -> n];
    grad = grd[func, vars];
    For[i = 0, i < steps, i++,
     vals = {vals, (Last@vals) - l*grad[Last@vals]}
    Cases[vals, {__}, {-2}]

It is used within a Manipulate[] so I'm trying to speed it up a bit. It returns a list of points, one from each step the gradient takes (i.e. {{1,2}, {3,2}, ...}). I understand using For[] is not optimal, but I'm not sure how to replace it.

Question: Is there a better, cleaner, and/or faster way to do this in Mathematica?

Also the function seems to work for most functions like $x^2+y^2$ and $\sin x\cos y$, but fails for a linear error function like below. The points grow rapidly up to the order of $10^{80}$ in only 10 steps. Why does it not work for this error function?

list = RandomReal[listMax, {listSize, 2}];
error = 1/listSize * Sum[(Last@list[[i]] - (m*First@list[[i]] + b))^2, {i, 1, listSize}];

The full source is here if you need something to test it.

  • $\begingroup$ In general, you need line search (e.g., Armijo backtracking) in order to make gradient methods convergent. $\endgroup$ Nov 6, 2017 at 21:58
  • $\begingroup$ @HenrikSchumacher why is my implementation any different than what's shown here? $\endgroup$
    – Dando18
    Nov 6, 2017 at 22:26
  • $\begingroup$ The code behind that link uses numerical differentiation. Apart from that, seemingly, no difference. But that does not make it working unless using tiny(!) step sizes (aka learning rates). $\endgroup$ Nov 7, 2017 at 7:03

1 Answer 1


For neural networks, one often prescribes a "learning rate", i.e. a constant step size. In is quite well known in optimization circles that this is a very, very bad idea as the gradient alone does not tell you how far you should travel without ascending the objective function (we want to descent!).

In the following, I show you an implementation of gradient descent with "Armijo step size rule with quadratic interpolation", applied to a linear regression problem. (Actually, with regression problems, it is often better to use the Gauss-Newton method.)

This is the code for the steepest descent. One has to supply a objective function f and a function generating its differential:

stepGradient[f_, Df_, start_, initialstepsize_, tolerance_, steps_] := 
  Module[{\[Sigma], \[Gamma], x, \[Phi]0, \[Phi]t, D\[Phi]0, DF, u, y, t, pts, iter, residual},
   \[Sigma] = 0.5; (*Armijo constant*)
   \[Gamma] = 0.5; (*shrinking factor for step sizes*)
   iter = 0;
   pts = {start};
   x = start;
   DF = Df[x];
   residual = Sqrt[DF.DF];
   While[residual > tolerance && iter < steps, (* gradient descent until Norm[Df[x]] < tolerance *)
    u = -DF;(*search direction*)

    (* Start Armijo backtracking *)
    (* merit function is given by \[Phi][t_] := f[x + t u ] *)
    \[Phi]0 = f[x]; (* \[Phi]0 = \[Phi][0] *)
    D\[Phi]0 = DF.u; (* D\[Phi]0 = \[Phi]'[0] *)
    y = x + t u;
    \[Phi]t = f[y]; (* \[Phi]t = \[Phi]'[t] *)
    t = initialstepsize;
    While[\[Phi]t > \[Phi]0 + t \[Sigma] D\[Phi]0,(* backtracking until f[x+t u]\[LessEqual]f[x]+\[Sigma] t Df[x].u *)
    tnew = - \[Sigma] t^2 D\[Phi]0/(\[Phi]t - \[Phi]0 - t D\[Phi]0);         (* compute minimum tnew of the quadratic interpolation through \[Phi][0], \[Phi][t] and \[Phi]'[0] *)
    t = Max[\[Gamma] t, tnew]; (* shrink step size *)
    y = x + t u;
    \[Phi]t = f[y];
   (* Armijo backtracking is done. Prepare for next iteration. *)

    x = y;
    DF = Df[x];
    residual = Sqrt[DF.DF];
    pts = {pts, x};
   Partition[Flatten[pts], Length[start]]

Simulating data sampled from a linear function perturbated by some Gaussian noise.

m0 = 1.5;
b0 = 4.;
listMax = {-1, 1};
listSize = 300;
list = Transpose[{Range[listSize]/listSize, m0 Range[listSize]/listSize + b0}] + 
   0.1 RandomVariate[NormalDistribution[], {listSize, 2}];

enter image description here

Setting up the fidelity function and its derivative. It is a good idea to do that outside of any loop since this involves symbolic computations which tend to be rather slow. So we perform them only once. (With big data, symbolic computations should be completely avoided. In this case, one should implement the derivative in a cleverer way.)

error = Total[(m list[[All, 1]] + b - list[[All, 2]])^2]/Length[list];
vars = {m, b};
Quiet[Block[{x, xx},
   xx = Table[x[[i]], {i, 1, Length[vars]}];
   f = x \[Function] Evaluate[Simplify[error] /. Thread[vars -> xx]];
   Df = x \[Function] Evaluate[D[f[xx], {xx, 1}]];
   (* DDf = x \[Function] Evaluate[D[f[xx], {xx, 2}]]; *)

Finally, we run steepest descent and plot the result:

pts = stepGradient[f, Df, {0., 0.}, 100., 0.0001, 100];
{m1, b1} = pts[[-1]]
 Plot[m1 t + b1, {t, 0, 1}],

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.