Why does scipy.optimize.minimize (default) report success without moving with Skyfield?

Question

scipy.optimize.minimize using default method is returning the initial value as the result, without any error or warning messages. While using the Nelder-Mead method as suggested by this answer solves the problem, I would like to understand:

Why does the default method returns the wrong answer without warning the starting point as the answer - and is there a way I can protect against "wrong answer without warning" avoid this behavior in this case?

Note, the function separation uses the python package Skyfield to generate the values to be minimized which is not guaranteed smooth, which may be why Simplex is better here.

RESULTS:

test result: [ 2.14159739] 'correct': 2.14159265359 initial: 0.0

default result: [ 10000.] 'correct': 13054 initial: 10000

Nelder-Mead result: [ 13053.81011963] 'correct': 13054 initial: 10000

FULL OUTPUT using DEFAULT METHOD:
   status: 0
  success: True
     njev: 1
     nfev: 3
 hess_inv: array([[1]])
      fun: 1694.98753895812
        x: array([ 10000.])
  message: 'Optimization terminated successfully.'
      jac: array([ 0.])
      nit: 0

FULL OUTPUT using Nelder-Mead METHOD:
  status: 0
    nfev: 63
 success: True
     fun: 3.2179306044608054
       x: array([ 13053.81011963])
 message: 'Optimization terminated successfully.'
     nit: 28

Here is the full script:

def g(x, a, b):
    return np.cos(a*x + b)

def separation(seconds, lat, lon):
    lat, lon, seconds = float(lat), float(lon), float(seconds) # necessary it seems
    place = earth.topos(lat, lon)
    jd = JulianDate(utc=(2016, 3, 9, 0, 0, seconds))
    mpos = place.at(jd).observe(moon).apparent().position.km
    spos = place.at(jd).observe(sun).apparent().position.km
    mlen = np.sqrt((mpos**2).sum())
    slen = np.sqrt((spos**2).sum())
    sepa = ((3600.*180./np.pi) *
            np.arccos(np.dot(mpos, spos)/(mlen*slen)))
    return sepa

from skyfield.api import load, now, JulianDate
import numpy as np
from scipy.optimize import minimize

data = load('de421.bsp')

sun   = data['sun']
earth = data['earth']
moon  = data['moon']

x_init = 0.0
out_g = minimize(g, x_init, args=(1, 1))
print "test result: ", out_g.x, "'correct': ", np.pi-1, "initial: ", x_init    # gives right answer

sec_init = 10000
out_s_def = minimize(separation, sec_init, args=(32.5, 215.1))
print "default result: ", out_s_def.x, "'correct': ", 13054, "initial: ", sec_init

sec_init = 10000
out_s_NM = minimize(separation, sec_init, args=(32.5, 215.1),
                 method = "Nelder-Mead")
print "Nelder-Mead result: ", out_s_NM.x, "'correct': ", 13054, "initial: ", sec_init

print ""
print "FULL OUTPUT using DEFAULT METHOD:"
print out_s_def
print ""
print "FULL OUTPUT using Nelder-Mead METHOD:"
print out_s_NM

Answer 1

1)

Your function is piecewise constant (has small-scale "staircase" pattern). It is not everywhere differentiable.

The gradient of the function at the initial guess is zero.

The default BFGS optimizer sees the zero gradient and decides it is a local minimum by its criteria (which are based on assumptions about the input function that are not true in this case, such as differentiability).

Basically, the exactly flat regions bomb the optimizer. The optimizer probes the function in the small exactly flat region around the initial point, where everything looks like the function is just a constant, so it thinks you gave it a constant function. Because your function is not differentiable everywhere, it is possible that almost all initial points are inside such flat regions, so that this can happen without bad luck in the choice of the initial point.

Note also that Nelder-Mead is not immune to this --- it just happens its initial simplex is larger than the size of the staircase, so it probes the function around a larger spot. If the initial simplex would be smaller than the staircase size, the optimizer would behave similarly as BFGS.

2)

General answer: local optimizers return local optima. Whether these coincide with the true optimum depends on the properties of the function.

In general, to see if you're stuck in a local optimum, try different initial guesses.

Also, using a derivative-based optimizer on a non-differentiable function is not a good idea. If the function is differentiable on a "large" scale, you can adjust the step size of the numerical differentiation.

Because there is no cheap/general way to check numerically if a function is everywhere differentiable, no such check is done --- instead it is an assumption in the optimization method that must be ensured by whoever inputs the objective function and chooses the optimization method.

Answer 2

The accepted answer by @pv. explains that Skyfield has a "staircase" response, meaning that some values it returns are locally flat except for discrete jumps.

I did a little experiment on the first step - converting times to JulianDate objects, indeed it looks like roughly 40 microseconds per increment, or about 5E-10 days. That's reasonable, considering the JPL databases span thousands of years. While this is probably fine for almost any general astronomical-scale application, it's not actually smooth. As the answer points out - the local flatness will trigger "success" in some (probably many) minimizers. This is expected and reasonable and is not in any way a failure of the method.

from skyfield.api import load, now, JulianDate
import numpy as np
import matplotlib.pyplot as plt

t  = 10000 + np.logspace(-10, 2, 25)        # logarithmic spacing
jd = JulianDate(utc=(2016, 3, 9, 0, 0, t))

dt  = t[1:] - t[:-1]
djd = jd.tt[1:] - jd.tt[:-1]

t  = 10000 + np.linspace(0, 0.001, 1001)        # linear spacing
jd = JulianDate(utc=(2016, 3, 9, 0, 0, t))

plt.figure()

plt.subplot(1,2,1)

plt.plot(dt, djd)
plt.xscale('log')
plt.yscale('log')

plt.subplot(1,2,2)

plt.plot(t, jd.tt-jd.tt[0])

plt.show()