pythonscipysplinecubic-splinefminsearch

fmin_slsqp returns initial guess finding the minimum of cubic spline


I am trying to find the minimum of a natural cubic spline. I have written the following code to find the natural cubic spline. (I have been given test data and have confirmed this method is correct.) Now I can not figure out how to find the minimum of this function.

This is the data

xdata = np.linspace(0.25, 2, 8)
ydata = 10**(-12) * np.array([1,2,1,2,3,1,1,2])

This is the function

import scipy as sp
import numpy as np
import math
from numpy.linalg import inv
from scipy.optimize import fmin_slsqp
from scipy.optimize import minimize, rosen, rosen_der

def phi(x, xd,yd):
    n = len(xd)
    h = np.array(xd[1:n] - xd[0:n-1])
    f = np.divide(yd[1:n] - yd[0:(n-1)],h)

    q = [0]*(n-2)

    for i in range(n-2):
        q[i] = 3*(f[i+1] - f[i])

    A = np.zeros(((n-2),(n-2)))

    #define A for j=0
    A[0,0] = 2*(h[0] + h[1])
    A[0,1] = h[1]

    #define A for j = n-2
    A[-1,-2] = h[-2]
    A[-1,-1] = 2*(h[-2] + h[-1])

    #define A for in the middle
    for j in range(1,(n-3)):
        A[j,j-1] = h[j]
        A[j,j] = 2*(h[j] + h[j+1])
        A[j,j+1] = h[j+1]

    Ainv = inv(A)

    B = Ainv.dot(q)

    b = (n)*[0]
    b[1:(n-1)] = B

    # now we find a, b, c and d

    a = [0]*(n-1)
    c = [0]*(n-1)
    d = [0]*(n-1)

    s = [0]*(n-1)

    for r in range(n-1):
        a[r] = 1/(3*h[r]) * (b[r + 1] - b[r])
        c[r] = f[r] - h[r]*((2*b[r] + b[r+1])/3)
        d[r] = yd[r]

    #solution 1 start
    for m in range(n-1):
        if xd[m] <= x <= xd[m+1]:
            s = a[m]*(x - xd[m])**3 + b[m]*(x-xd[m])**2 + c[m]*(x-xd[m]) + d[m]

    return(s)
    #solution 1 end

I want to find the minimum on the domain of my xdata, so a fmin didn't work as you can not define bounds there. I tried both fmin_slsqp and minimize. They are not compatible with the phi function I wrote so I rewrote phi(x, xd,yd) and added an extra variable such that phi is phi(x, xd,yd, m). M indicates in which subfunction of the spline we are calculating a solution (from x_m to x_m+1). In the code we replaced #solution 1 by the following

# solution 2 start

return(a[m]*(x - xd[m])**3 + b[m]*(x-xd[m])**2 + c[m]*(x-xd[m]) + d[m])
# solution 2 end

To find the minimum in a domain x_m to x_(m+1) we use the following code: (we use an instance where m=0, so x from 0.25 to 0.5. The initial guess is 0.3)

fmin_slsqp(phi, x0 = 0.3, bounds=([(0.25,0.5)]), args=(xdata, ydata, 0))

What I would then do (I know it's crude), is iterate this with a for loop to find the minimum on all subdomains and then take the overall minimum. However, the function fmin_slsqp constantly returns the initial guess as the minimum. So there is something wrong, which I do not know how to fix. If you could help me this would be greatly appreciated. Thanks for reading this far.


Solution

  • When I plot your function phi and the data you feed in, I see that its range is of the order of 1e-12. However, fmin_slsqp is unable to handle that level of precision and fails to find any change in your objective.

    The solution I propose is scaling the return of your objective by the same order of precision like so: return(s*1e12)

    Then you get good results.

    >>> sol = fmin_slsqp(phi, x0=0.3, bounds=([(0.25, 0.5)]), args=(xdata, ydata))
    >>> print(sol)
    Optimization terminated successfully.    (Exit mode 0)
            Current function value: 1.0
            Iterations: 2
            Function evaluations: 6
            Gradient evaluations: 2
    [ 0.25]
    

    Plot of Phi