Scipy won't find the optimum (simple cosine function)

Question

I am trying to estimate the argument of a cosine function using the scipy optimizer (yes I am aware arc cos could be used, but I don't want to do that).

The code + a demonstration:

import numpy
import scipy

def solver(data):
    Z=numpy.zeros(len(data))
    a=0.003
    for i in range(len(data)):
        def minimizer(b):
            return numpy.abs(data[i]-numpy.cos(b))
        Z[i]=scipy.optimize.minimize(minimizer,a,bounds=[(0,numpy.pi)],method="L-BFGS-B").x[0]
    return Z

Y=numpy.zeros(100)
for i in range(100):
   Y[i]=numpy.cos(i/25)

solver(Y)

The result is not good, when the argument of the cos function reaches values above 2, the estimation "skips over" the values and returns the maximum argument value instead.

array([0.        , 0.04      , 0.08      , 0.12      , 0.16      ,
       0.2       , 0.24      , 0.28      , 0.32      , 0.36      ,
       0.4       , 0.44      , 0.48      , 0.52      , 0.56      ,
       0.6       , 0.64      , 0.67999999, 0.72      , 0.75999999,
       0.8       , 0.83999999, 0.88      , 0.92      , 0.95999999,
       1.        , 1.04      , 1.08      , 1.12      , 1.16      ,
       1.2       , 1.24      , 1.28      , 1.32      , 1.36      ,
       1.4       , 1.44      , 1.48      , 1.52      , 1.56      ,
       1.6       , 1.64      , 1.68      , 1.72      , 1.76      ,
       1.8       , 1.84      , 1.88      , 1.91999999, 1.95999999,
       2.        , 2.04      , 3.14159265, 3.14159265, 3.14159265,
       3.14159265, 3.14159265, 3.14159265, 3.14159265, 3.14159265,...

What causes this phenomenon? Are there some other optimizers/ settings that could help with the issue?

Probably not the answer you are looking for, but preprocessing your inputs with `val %= 2 * numpy.pi` would probably do the trick. — Dillon Davis, Jul 06 '18 at 07:13
@Dole Apart from that the code is not runnable - I hope that other people will understand better than me what you are trying to achieve and what your strategy is to achieve this. — Mr. T, Jul 06 '18 at 09:27
@Mr. T Sorry, the array was not initialized. I would like to have an optimizer that can find the cosine argument. So if it is fed cos(2.4) it should returns 2.4. Instead as can be seen, it returns 3.14... — Dole, Jul 06 '18 at 11:49
@Dole Not sure why it doesn't work but this is definitely a job for [`scipy.optimize.root`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.root.html). What you want to find is the root of the function `data[i] - np.cos(b)`. — a_guest, Jul 06 '18 at 11:53

score 3 · Accepted Answer · answered Jul 06 '18 at 12:57

The reason is that for the function (for example) f = abs(cos(0.75*pi) - cos(z)) the gradient f' happens to vanish at z = pi, as can be seen from the following plot:

If you check the result the optimization procedure then you'll see that:

      fun: array([0.29289322])
 hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>
      jac: array([0.])
  message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     nfev: 16
      nit: 2
   status: 0
  success: True
        x: array([3.14159265])

So the optimization procedure reached one of its convergence criteria. More detailed information about the criterion can be found at the L-BFGS-B documentation. It says that

gtol : float

The iteration will stop when max{|proj g_i | i = 1, ..., n} <= gtol where pg_i is the i-th component of the projected gradient.

So it eventually reaches a point z >= pi which is then projected back to z = pi due to the constraint and at this point the gradient of the function is zero and hence it stops. You can observe that by registering a callback which prints the current parameter vector:

def new_callback():
    step = 1

    def callback(xk):
        nonlocal step
        print('Step #{}: xk = {}'.format(step, xk))
        step += 1

    return callback

scipy.optimize.minimize(..., callback=new_callback())

Which outputs:

Step #1: xk = [0.006]
Step #2: xk = [3.14159265]

So at the second step it hit z >= pi which is projected back to z = pi.

You can circumvent this problem by reducing the bounds to for example bounds=[(0, 0.99*np.pi)]. This will give you the expected result, however the method won't converge; you will see something like:

      fun: array([1.32930966e-09])
 hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>
      jac: array([0.44124484])
  message: b'ABNORMAL_TERMINATION_IN_LNSRCH'
     nfev: 160
      nit: 6
   status: 2
  success: False
        x: array([2.35619449])

Note the message ABNORMAL_TERMINATION_IN_LNSRCH. This is due to the nature of abs(x) and the fact that its derivative has a discontinuity at x = 0 (you can read more about that here).

Alternative approach (finding the root)

For all the lines above we were trying to find a value z for which cos(z) == cos(0.75*pi) (or abs(cos(z) - cos(0.75*pi)) < eps). This problem is actually finding the root of the function f = cos(z) - cos(0.75*pi) where we can make use of the fact that cos is a continuous function. We need to set the boundaries a, b such that f(a)*f(b) < 0 (i.e. they have opposite sign). For example using bisect method:

res = scipy.optimize.bisect(f, 0, np.pi)

I tried the bisect method. However, it complains that the signs are not different even though they are, is there a hack? The least squares worked well, but all the tolerances need to be maxed or the error grows rapidly. — Dole, Jul 08 '18 at 03:29
@Dole If it complains about similar signs, then that's probably because they are indeed the same. Could you provide the definition of `f` as well as the case for which it doesn't work? Better triple check the call to `bisect` as well. What is the exact error message? — a_guest, Jul 08 '18 at 08:37

score 3 · Answer 2 · answered Jul 06 '18 at 13:39

Besides the general minimize method, SciPy has minimize_scalar specifically for 1-dimensional problems like here, and least_squares for minimizing a particular kind of functions that measure the difference between two quantities (such as the difference between cos(b) and diff[i] here). The latter performs well here, even without fine-tuning.

for i in range(len(data)):
    Z[i] = scipy.optimize.least_squares(lambda b: data[i] - numpy.cos(b), a, bounds=(0, numpy.pi)).x[0]

The function passed to least_squares is the thing we'd like to be close to 0, without an absolute value on it. I'll add that a=0.003 seems a suboptimal choice for a starting point, being so close to the boundary; nonetheless it works.

Also, as a_guest already posted, a scalar root finding method should do the same thing while throwing fewer surprises here, given that we already have a nice bracketing interval [0, pi]. Bisection is reliable but slow; Brent's method is what I'd probably use.

for i in range(len(data)):
    Z[i] = scipy.optimize.brentq(lambda b: data[i] - numpy.cos(b), 0, numpy.pi)

Scipy won't find the optimum (simple cosine function)

2 Answers2

Alternative approach (finding the root)