title: SciPy optimisation: Newton-CG vs BFGS vs L-BFGS
body: <p>I am doing an optimisation problem using Scipy, where I am taking a flat network of vertices and bonds of size <code>NNxNN</code>, connecting two sides of it (i.e., making it periodic), and minimising an energy function, so that it curls up to form a cylinder. (See the links below.)</p>
  
  <p>Since I have the function <code>energy(xyz-position)</code> and it's gradient, I decided to use the three methods recommended in the Scipy manual -- <code>Newton-CG</code>, <code>BFGS</code>, <code>L-BFGS-B</code> -- and compare how they performed.</p>
  
  <p>I call the optimisation function as follows, and I merely replace <code>'Newton-CG'</code> with <code>'BFGS'</code> and <code>'L-BFGS-B'</code> according to case:</p>
  
  <pre><code>from scipy.optimize import minimize
  res = minimize(energy, xyzInit, method='Newton-CG', jac = energy_der,  options={'disp': True}) 
  </code></pre>
  
  <p>I found the following general behaviour (I am giving the output data for the case of <code>NN=9</code>, corresponding to a <code>3*9^2=243</code>-dimensional parameter space) -</p>
  
  <ol>
  <li><p>BFGS systematically failed to find the correct minimum (for low <code>NN</code>), and failed to converge at all for large <code>NN</code>. See <a href="https://plot.ly/~apal90/162/" rel="noreferrer">https://plot.ly/~apal90/162/</a> for end result.</p>
  
  <pre><code> NN=9
   Method: BFGS
   Warning: Desired error not necessarily achieved due to precision loss.
   Current function value: 204.465912
   Iterations: 1239
   Function evaluations: 1520
   Gradient evaluations: 1508
   Time taken for minimisation: 340.728140116
  </code></pre></li>
  <li><p>Newton-CG found the correct minimum for small <code>NN</code> (&lt;=8), but starting from NN=9, returned an incorrect minimum (viz., a cylinder squashed at one end), and for higher values stopped even converging. Note: This behaviour was for some reason aggravated for odd <code>NN</code>'s. See <a href="https://plot.ly/~apal90/164/" rel="noreferrer">https://plot.ly/~apal90/164/</a></p>
  
  <pre><code> NN=9
   Method: Newton-CG
   Optimization terminated successfully.
   Current function value: 7.954412
   Iterations: 49
   Function evaluations: 58
   Gradient evaluations: 1654
   Hessian evaluations: 0
   Time taken for minimisation: 294.203114033
  </code></pre></li>
  <li><p>L-BFGS-B found the correct minimum, and that too blazingly fast, for all <code>NN</code>'s that I tested (up to <code>NN=14</code>). See <a href="https://plot.ly/~apal90/160/" rel="noreferrer">https://plot.ly/~apal90/160/</a></p>
  
  <pre><code> NN=9
   Method: L-BFGS-B
   Time taken for minimisation: 36.3749790192
  </code></pre></li>
  </ol>
  
  <p><strong>Question</strong>: Why is <code>L-BFGS-B</code> superior in this case to the other two methods? In particular, why is it so much superior to <code>BFGS</code>, when both are supposed to be quasi-Newton methods that work (to my understanding), in exactly the same manner.</p>
  
  <p><strong>My thoughts on the situation</strong>: All three methods do quadratic approximations at every point x. For this, it needs a gradient and a Hessian. If the Hessian is not given, it must be calculated by the algorithm. In our case, where only the gradient is explicitly given, this is calculated at every step numerically by the algorithm. More specifically, what we require is the inverse of the Hessian, and this is a very expensive step, especially in higher dimensions.  Now, Newton-CG calculates this inverse Hessian explicitly, hence it's longer time requirements. The quasi-Newton methods like BFGS and L-BFGS calculate an approximation to the Hessian (i.e., the curvature) based on the gradient, which is cheaper on time, and which is also supposedly a better estimate of the curvature about a point. Thus, for quadratic functions, Newton-CG converges faster, whereas for non-quadratic functions, the quasi-Newton functions converge better. L-BFGS is a lower memory version of BFGS that stores far less memory at every step than the full NxN matrix, hence it is faster than BFGS.</p>
  
  <p>This explanation shows a divergence between Newton-CG and the quasi-Newton methods. What it does not explain is the inability of the algorithms to find the true minimum, and especially the disparity between BFGS and L-BFGS, which are both supposed to function in the same manner. </p>
  
  <p>My general hunch on the long convergence times is that the system is non-quadratic (i.e. flat) about the minimum, and thus the algorithm oscillates about with converging. Beyond that, if BFGS and L-BFGS truly work in the same manner, I believe there must be some difference between the convergence tolerance levels of the Scipy algorithms. Otherwise, BFGS and L-BFGS don't really work in the same manner, and the latter probably calculates the Hessian far more accurately.</p>
  
  <p>References --</p>
  
  <p><a href="http://www.scipy-lectures.org/advanced/mathematical_optimization/#newton-and-quasi-newton-methods" rel="noreferrer">http://www.scipy-lectures.org/advanced/mathematical_optimization/#newton-and-quasi-newton-methods</a></p>
  
  <p><a href="https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization" rel="noreferrer">https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization</a></p>
  
  <p><a href="https://en.wikipedia.org/wiki/Quasi-Newton_method" rel="noreferrer">https://en.wikipedia.org/wiki/Quasi-Newton_method</a></p>
  
  <p><a href="https://docs.scipy.org/doc/scipy-0.18.1/reference/optimize.minimize-bfgs.html#optimize-minimize-bfgs" rel="noreferrer">https://docs.scipy.org/doc/scipy-0.18.1/reference/optimize.minimize-bfgs.html#optimize-minimize-bfgs</a></p>
  
  <p><a href="https://docs.scipy.org/doc/scipy-0.18.1/reference/optimize.minimize-lbfgsb.html#optimize-minimize-lbfgsb" rel="noreferrer">https://docs.scipy.org/doc/scipy-0.18.1/reference/optimize.minimize-lbfgsb.html#optimize-minimize-lbfgsb</a></p>