1

I was trying to minimize the objective function while using a for loop to set the constraints such that x1 = x2 = ... xn. However, the optimization doesn't seem to work. I.e. the end x still equals to the initial x. And I am getting an error message of 'Singular matrix C in LSQ subproblem'.

covariance_matrix = np.matrix([[0.159775519, 0.022286316, 0.00137635, -0.001861736],
                     [0.022286316, 0.180593862, -5.5578e-05, 0.00451056], 
                     [0.00137635, -5.5578e-05, 0.053093075, 0.02240866], 
                     [-0.001861736, 0.00451056, 0.02240866, 0.053778594]]) 

x0 = np.matrix([0.2,0.2,0.3,0.4])


fun = lambda x: x.dot(covariance_matrix).dot(x.transpose())
cons = np.array([])
for i in range(0,x0.size-1):
    con = {'type': 'eq', 'fun': lambda x:  x[i] - x[i+1]}    
    cons = np.append(cons, con)

con = {'type': 'eq', 'fun': lambda x:  sum(x)-1}   
cons = np.append(cons, con) 

solution = minimize(fun,x0,method='SLSQP',constraints = cons)


solution message:   Singular matrix C in LSQ subproblem
solution status:   6
solution success:   False

But if I append the constraints one by one, then it works perfectly, meaning the result gives me x1 = x2 = x3 = x4

con1 = {'type': 'eq', 'fun': lambda x:  sum(x)-1}   
con2 = {'type': 'eq', 'fun': lambda x:  x[1]-x[0]}   
con3 = {'type': 'eq', 'fun': lambda x:  x[2]-x[1]}   
con4 = {'type': 'eq', 'fun': lambda x:  x[3]-x[2]}   
cons = np.append(cons, con1) 
cons = np.append(cons, con2) 
cons = np.append(cons, con3) 
cons = np.append(cons, con4) 

solution message:   Optimization terminated successfully.
solution status:   0
solution success:   True
dickli2119
  • 145
  • 1
  • 14
  • "*... the optimization doesn't seem to work. I.e. the end x still equals to the initial x."* There is no need to guess that it "doesn't seem to work". The `solution` object *tells* you it didn't work by setting `solution.success` to False. A short description of the reason is in `solution.message`. Always check `solution.success` after calling `minimize()`. Could you edit the question to show `solution.success`, `solution.status` and `solution.message` after calling `minimize()`? – Warren Weckesser Aug 03 '17 at 18:01
  • Thank you Warren! I have edited the initial problem. And it appears that I am getting a 'Singular matrix C in LSQ subproblem' problem. Do you know where this is coming from? – dickli2119 Aug 03 '17 at 18:12
  • It looks like that message was copied from the source code of the underlying Fortran solver. Unfortunately, I don't see a description of the "matrix C" in the docstrings for `minimize` or `fmin_slsqp`, so someone would have to dig into the Fortran code (at least the comments in the Fortran code) to figure out what C is. The code in scipy is [slsqp_optmz.f](https://github.com/scipy/scipy/blob/master/scipy/optimize/slsqp/slsqp_optmz.f). – Warren Weckesser Aug 03 '17 at 18:44
  • I'll add, though, that your equality constraints uniquely determine the solution--there is nothing that can be optimized in your problem. – Warren Weckesser Aug 03 '17 at 18:45
  • Warren, could you elaborate a little bit more? I thought the set up is to minimize the objective functions. And my main problem is that, if I append those constraints one by one then it works. but if I use a loop then it doesn't. – dickli2119 Aug 03 '17 at 19:39
  • Your constraints say `x[0] = x[1]`, `x[1] = x[2]`, `x[2] = x[3]` and `sum(x) = 1`. The only solution to those equations is `x[0] = x[1] = x[2] = x[3] = 0.25`. So there is nothing to optimize--there is no freedom left in the variables. – Warren Weckesser Aug 03 '17 at 19:41
  • Yea, i am setting it up in this way because it's easier to test the code. The initial x's are not 0.25. And if the optimization is run properly then it should give me x's = 0.25. – dickli2119 Aug 03 '17 at 19:44
  • Show the code where you append the constraints one by one and it works. – Warren Weckesser Aug 03 '17 at 19:52
  • Just added the code at the bottom. Really appreciate the help Warren! – dickli2119 Aug 03 '17 at 20:03
  • Argh, I should have seen the problem sooner, considering that I answered a similar question over a year ago: https://stackoverflow.com/questions/37791680/scipy-optimize-minimize-slsqp-with-linear-constraints-fails – Warren Weckesser Aug 03 '17 at 20:12
  • In the loop, change the creation of `con` to ` con = {'type': 'eq', 'fun': lambda x, i=i: x[i] - x[i+1]}` – Warren Weckesser Aug 03 '17 at 20:13
  • Thank you very much!! It's working perfectly now. Also really appreciate for sharing the article!! – dickli2119 Aug 03 '17 at 20:25
  • @WarrenWeckesser Hi Warren. Thank you for your previous help. I have a follow-up question regarding optimization. I posted on stack overflow, but have yet solved the problem. My question is that in the case where I need to set one extra constraint for example x1 + x2 >0.8, then mathematically the original constraint x1 = x2 = ... = xn does not work anymore because x1 must be bigger than xn. I would love to hear your suggestion. You could refer to my original post: https://stackoverflow.com/questions/45517107/what-constraints-should-i-set-when-i-try-to-run-optimization-to-build-an-equal-w – dickli2119 Aug 05 '17 at 16:39
  • It sounds like you have to figure out the mathematical model of the problem that you are trying to solve. That's not a programming question, so stackoverflow isn't the best place for that question. – Warren Weckesser Aug 05 '17 at 23:53

1 Answers1

6

(Note: while the details are different, this question is about the same problem as Scipy.optimize.minimize SLSQP with linear constraints fails)

Your loop is

for i in range(0,x0.size-1):
    con = {'type': 'eq', 'fun': lambda x:  x[i] - x[i+1]}    
    cons = np.append(cons, con)

The problem is the use of i in the lambda expression. Python closures are "late binding". That means the value of i that is used when the function is called is not the same as the value of i when the function was created. After the loop, the value of i is 2, so the expression evaluated by all the functions created in the loop is x[2] - x[3]. (That also explains the "singular matrix C" referred to in the error message.)

One way to fix this is to make i an argument of the lambda expression whose default value is the current i:

    con = {'type': 'eq', 'fun': lambda x, i=i:  x[i] - x[i+1]}
Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214