Optimal value of parameter/s from a set by scipy.optimize.minimize() method

Question

I have this function $y = \exp(-x)$. I have a list of $x$ values and corresponding $y$ values.

X = [-3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.,  6.]

Y = [20.085536923187668, 7.38905609893065, 2.718281828459045, 1.0, 0.36787944117144233, 0.1353352832366127, 0.049787068367863944, 0.01831563888873418, 0.006737946999085467, 0.0024787521766663585]

For each list element, I am getting $y-x = 0$, which is good.

for xi,yi in zip (X,Y):
    print(yi-np.exp(-xi))

What I am trying to do is to change my equation slightly and find the optimum value of the parameter $B$. So the modified equation is, $y = \exp(-B*x)$ Ideally, I should get a $B=1$ for each of the list elements of X and Y.

The $B$ value is 1 as expected.

def test_cost5(param):
    B=param
    
    #print(b)
    return Y[3]-np.exp(-beta*X[3])

test_initial_guess=1
results=optimize.minimize(test_cost5, test_initial_guess, method='BFGS')
results

Gives output,

fun: array([0.])
 hess_inv: array([[1]])
      jac: array([0.])
  message: 'Optimization terminated successfully.'
     nfev: 3
      nit: 0
     njev: 1
   status: 0
  success: True
        x: array([1.])

Can I have some help on the following problems?

1. Can I use one more parameter apart from $B$?
2. I want to optimize these multiple parameters for a set of data and find the optimal values of the parameters. Can it be done?

Answer 1

For your first question, you can optimize as many parameters. You simply need to pass in a list array as an initial guess to your function. For the second question, if you want to optimize across the full set rather than just single pair of $x/y$, you should pass these whole arrays in as arguments to the optimize routine.

So if you wanted to fit say $y=\exp(\beta x+\alpha)$, you could do something like

import numpy as np    
from scipy import optimize

X = np.array([-3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.,  6.])

Y = np.array([20.085536923187668, 7.38905609893065, 2.718281828459045,
              1.0, 0.36787944117144233, 0.1353352832366127,0.049787068367863944,
               0.01831563888873418, 0.006737946999085467, 0.0024787521766663585])

def shifted_exp_error(params,*args):
    beta,alpha=params
    x,y=args
    
    error=y-np.exp(beta*x+alpha)
    mse=np.mean(error**2)
    print(beta,alpha)
    print(mse)
    return mse

results=optimize.minimize(shifted_exp_error,[-1.0,0.0], args=(X,Y), method='BFGS',options={'eps':1e-5,'gtol':1e-2})


print(results)

Notice, since I'm trying to optimize $\alpha$ and $\beta$ for all the $y$ values, I need some single value to measure the error across the whole set, since optimize.minimize can only minimize scalar functions. So, my function computes the mean squared error (MSE) for all the values in $y$ and minimizes that error. In general, you will also need to play around with eps (the step size for numerical derivatives) and gtol (how small the gradient needs to be for the minimization to converge).

For the example given, there are better ways to solve it (e.g. take the $\ln(y)$ and solve for the parameters using linear least squares regression), but this example should show the basics for a more general nonlinear optimization. Scipy also has a function for nonlinear least-squares that works well for this problem. All you would need to change to use this routine is to have shifted_exp_error return error instead of mse. There is also a curve_fit function.

I left print statements in my error function so you can what the parameters/error are throughout the optimization. Note that since you don't pass in a function to compute the derivative, it does numerical derivatives. So if you look at the output for this, the first three steps are:

1. Take a small step to compute the derivative with respect to $\beta$
2. Take a small step to compute the derivative with respect to $\alpha$
3. Take a larger step along the gradient to get your next candidate solution.