Given a function of multiple variables, and some initial conditions, I would like an efficient way to track the gradient to the local minimum of that function. Two options spring to mind — to either use NDSolve, or FindMinimum. I don't know an efficient way to make NDSolve identify when it has found a minimum, so I'm currently playing with FindMinimum. My problem is that FindMinimum is too efficient in the sense that it takes very large steps initially.

FindMinimum[x^2 + y^4, {x, 10}, {y, 10}, EvaluationMonitor :> Sow[{x, y}]
]]], PlotRange -> All, PlotMarkers -> {Automatic, 10}]

The initial steps are much larger than I would like

What I would like is a dense string of data points along the trajectory to the minimum as determined by gradient flow. Maybe FindMinimum is not the best option, but it seems like it may be, provided I can place a constraint on the maximum step size. Is this possible? Is there a better way?

  • $\begingroup$ Actually, does this method even follow the gradient?! Maybe I made a typo somewhere, or I am misunderstanding something, but I would expect the evolution to be primarily in the y-direction initially. $\endgroup$ – user41147 Apr 20 '17 at 19:55

From the tutorial Line Search Methods, there is an example similar to this:

Newton's method effectively uses a quadratic model and solves the equation $$H {\bar s} = - \nabla f(x,y)\,,$$ where $H$ is the Hessian $H = \nabla^2f(x,y)$, for the step $\bar s=(\Delta x, \Delta y)$. For an objective function that is a quadratic function like x^2 + 10 y^2, this will solve exactly for the minimum. With a step size reduction effected by "MaxRelativeStepSize" -> .1, this results in a straight-line trajectory toward the origin. To get a true gradient step, we need to set the Hessian to the identity matrix, which can be done via the Hessian option.

This shows the difference on the OP's original problem between the default Newton's method (purple) and the naive gradient descent method (red).

 ContourPlot[x^2 + y^4, {x, -0.5, 10.5}, {y, -0.5, 10.5}, 
  Contours -> 2^Range[15]],
    FindMinimum[x^2 + y^4, {x, 10}, {y, 10}, 
     Method -> {"Newton", 
       "StepControl" -> {"LineSearch", "MaxRelativeStepSize" -> .1}},
     EvaluationMonitor :> Sow[{x, y}]]]], PlotRange -> All, 
  PlotMarkers -> {Automatic, 10}, PlotStyle -> Purple],
    FindMinimum[x^2 + y^4, {x, 10}, {y, 10}, 
     Method -> {"Newton", Hessian -> {{1, 0}, {0, 1}},    (* <-- override Hessian *)
       "StepControl" -> {"LineSearch", "MaxRelativeStepSize" -> .1}},
     EvaluationMonitor :> Sow[{x, y}], MaxIterations -> 200]]], 
  PlotRange -> All, PlotMarkers -> {Automatic, 10}, PlotStyle -> Red]

Mathematica graphics

Note the gradient descent does not "converge to the requested accuracy or precision within 200 iterations" (FindMinimum::cvmit).

  • $\begingroup$ This is exactly what I was looking for. Thank you! $\endgroup$ – user41147 Apr 21 '17 at 11:32

You can use FindRoot to find zeros of the derivative. FindRoot as an option DampingFactor that should serve your purpose:

   FindRoot[{D[x^2 + y^4, x] == 0, 
     D[x^2 + y^4, y] == 0}, {{x, 10}, {y, 10}}, DampingFactor -> .1, 
    EvaluationMonitor :> Sow[{x, y}]]]], PlotRange -> All, 
 PlotMarkers -> {Automatic, 10}]

enter image description here


The algorithm for FindMinimum and FindRoot is Newton by default, not a gradient method. If one wants that, it's best to implement it yourself:

f[x_, y_] := x^2 + y^4
list = Module[{d = 10^-3., steps = 10000, x = 10, y = 10, dx, dy, xn, 
    yn, xp, yp},
   xn = x;
   yn = y;
    x = xn;
    y = yn;
    dx = D[f[xp, yp], xp] /. {xp -> x, yp -> y};
    dy = D[f[xp, yp], yp] /. {xp -> x, yp -> y};
    xn = x - d*dx;
    yn = y - d*dy;
    {x, y, dx, dy},
    {i, 0, steps}]];

ListLinePlot[list[[All, {1, 2}]], PlotRange -> All, PlotMarkers -> {Automatic, 10}]

enter image description here

Of course, there is plenty of room for improving the algorithm and making it more robust / choosing a smarter step size.

  • $\begingroup$ Thanks! I'm a bit confused though, this method does not seem to follow the gradient. If I change the equations to 'D[x^2 + 10 y^2, x] == 0, D[x^2 + 10 y^2, y] == 0' then the trajectory is a straight line. $\endgroup$ – user41147 Apr 20 '17 at 19:26
  • $\begingroup$ I thought your question is how to increase the number of points on the way to the minimum. Of course, you need to calculate the gradient at each of these points. Right now, it just plots the trajectory (x,y). $\endgroup$ – Felix Apr 20 '17 at 19:56
  • $\begingroup$ Given some initial conditions, I would like to plot points along a trajectory that follows the steepest descent, eventually ending up at the minimum. I thought this was how FindMinimum worked. $\endgroup$ – user41147 Apr 20 '17 at 20:00

Just for fun, here is a version using NDSolve:

sol = NDSolve[
        x'[t]==-2x[t] Exp[t],
        y'[t]==-4y[t]^3 Exp[t],
        WhenEvent[x'[t]^2 + y'[t]^2 < 10^-10, end=t; "StopIntegration"]

Note that I added a scaling factor (Exp[t]) so that the t range is not too large. Here's a plot of the result:

    ContourPlot[x^2 + y^4, {x, -1, 11}, {y, -1, 11}, Contours->2^Range[15]],
    ParametricPlot[{x[t], y[t]} /. sol, {t, 0, end}]

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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