0

I have a set of data; each column corresponds to a spectrum at a certain time. I want to fit the spectrum at a generic time (t_i) as a linear combination of the spectrum at time 0 (in the first column), at time 5 (in column 30) and time 35 (in column 210). So the equation I want to fit is:

S(t_i) = a * S(t_0) + b * S(t_5) + c * S(t_35)

where:

  1. 0 <= a, b, c <= 1

  2. a + b + c = 1

I found the solution at this question (Minimizing Least Squares with Algebraic Constraints and Bounds) super useful. But when I try it with my set of data the results are obviously wrong. I tried modifying the method to 'Nelder-Mead' but it doesn't respect my bound so I get negative values.

This is my script:

t0= df.iloc[:,0]    #Spectrum at time 0
t5 = df.iloc[:,30]  # Spectrum at time 5
t35 = df.iloc[:,120] # Spectrum at time 35
ti= df.iloc[:,20]
# Bounds that make every coefficient be between 0 and 1
bnds = [(0, 1), (0, 1), (0, 1)]
# Constrain the sum of the coefficient to 1
cons = [{"type": "eq", "fun": lambda x: x[0] + x[1] + x[2] - 1}]
xinit = np.array([1, 0, 0])
fun = lambda x: np.sum((ti -(x[0] * t0 + x[1] * t5 + x[2] * t35))**2)
res = minimize(fun, xinit,method='Nelder-Mead', bounds=bnds, constraints=cons)
print(res.x)

If I use the Nelder-Mead method I get: Out: [ 0.02732053 1.01961422 -0.04504698] , if I don't specify the method I get: [1. 0. 0.] (I believe that in this case the SLSQP method is being used).

The data I'm referring to is similar to the following:

0            3.333       5           35.001
0.001045089 0.001109701 0.001169798 0.000725486
0.001083051 0.001138815 0.001176665 0.000713021
0.001090994 0.001142676 0.001186642 0.000716149
0.001096258 0.001156476 0.001190218 0.00071286

Can you identify the problem? Can you suggest other ways to solve this problem? I have also tried using least_squares, but it failed.

Community
  • 1
  • 1
Francesca Mazzotta
  • 13
  • 2
  • Nelder-Mead does not take bounds or constraints. See the keywords at https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html – Joe May 21 '20 at 13:44
  • If you need both you have to use `SLSQP` or `trust-constr`. – Joe May 21 '20 at 13:46
  • Also see https://stackoverflow.com/a/52438791/7919597 – Joe May 21 '20 at 13:48
  • @Joe, I have tried SLSQP, but it gives wrong outputs. For example, by setting t_i = t5 which should give a=0, b=1 and c=0 but instead it gives 1, 0, 0. Using trust-constr I get better values but they are still wrong. I also read about this problem of Nelder-Mead that's why I'm having a hard time solving the problem. – Francesca Mazzotta May 21 '20 at 15:01
  • These are local optimizers, see my answer https://stackoverflow.com/a/52438791/7919597 The suggestions are the same. Try different initial values for `x0`. – Joe May 21 '20 at 15:47

1 Answers1

0

The result of a local optimization strongly depends on the initial values.

It might return [1, 0, 0] for the case you stated above because there simply was no possibility for the optimizer to find a "downhill-only" way to [0. 1. 0.]. In fact, you might have started in a local minima and all ways out of the dip went uphill. So the optimizer chose to stay. That's how these optimizers work.

Try

xinit = np.array([0.0, 1.0, 0.0])

for t_i = t5 and I am quite sure the optimizer will return the initial value.

For your case do what I stated here: Run the optimizer several times, each time pick random initial values inside your boundaries. You can pick the code posted there and just add your constraints, use SLSQP or trust-constr.

Joe
  • 6,758
  • 2
  • 26
  • 47