5

I often have to solve nonlinear problems in which the number of variables exceeds the number of constraints (or sometimes the other way around). Usually some of the constraints or variables are redundant in a complicated way. Is there any way to solve such problems?

Most of the scipy solvers seem to assume that the number of constraints equals the number of variables, and that the Jacobian is nonsingular. leastsq works sometimes but it doesn't even try when the constraints are fewer than the number of variables. I realize that I could just run fmin on linalg.norm(F), but this is much less efficient than any method which makes use of the Jacobian.

Here is an example of a problem which demonstrates what I am talking about. It obviously has a solution, but leastsq gives an error. Of course, this example is easy to solve by hand, I just put it here to demonstrate the issue.

import numpy as np
import scipy.optimize

mat = np.random.randn(5, 7)

def F(x):
    y = np.dot(mat, x)
    return np.array([ y[0]**2 + y[1]**3 + 12, y[2] + 17 ])

x0 = np.random.randn(7)
scipy.optimize.leastsq(F, x0)

The error message I get is:

Traceback (most recent call last):
  File "question.py", line 13, in <module>
    scipy.optimize.leastsq(F, x0)
  File "/home/dstahlke/apps/scipy/lib64/python2.7/site-packages/scipy/optimize/minpack.py", line 278, in leastsq
    raise TypeError('Improper input: N=%s must not exceed M=%s' % (n,m))
TypeError: Improper input: N=7 must not exceed M=2

I have scoured the net for an answer and have even asked on the SciPy mailing list, and got no response. For now I hacked the SciPy source so that the newton_krylov solver uses pinv(), but I don't think this is an optimal solution.

Cairnarvon
  • 25,981
  • 9
  • 51
  • 65
Dan Stahlke
  • 1,417
  • 1
  • 15
  • 20
  • Is this really a scipy question or is it actually a mathematical one wearing a fake beard? – talonmies Feb 18 '12 at 16:48
  • I believe it is a scipy question. I can solve these types of problems using custom solvers I have written but would prefer to be able to use the existing scipy solvers. Additionally, matlab's fsolve seems to be able to solve these. This seems to be a common situation and it is kind of hard to believe that scipy can't handle it (seemingly). – Dan Stahlke Feb 18 '12 at 17:08
  • `fsolve` uses a trust region method IIRC. So are you really wanting to know whether there is an analogous scipy function to `fsolve`? – talonmies Feb 18 '12 at 17:27
  • I don't care so much about the details of the implementation (maybe I should), I just want to know if there are any solvers at all in scipy that can handle these types of problems, other than using the one dimensional minimizers on norm(F(x)). Or failing that, whether there are any other open source solvers in python that can do it. At this point, I'd even be willing to write a wrapper for a C library if there is a good one. I looked at a few, but none of them seemed to handle the cases of degenerate (i.e. redundant) variables. – Dan Stahlke Feb 18 '12 at 17:56

1 Answers1

3

How about resize the return array from F() to the number of variables:

import numpy as np
import scipy.optimize

mat = np.random.randn(5, 7)

def F(x):
    y = np.dot(mat, x)
    return np.resize(np.array([ y[0]**2 + y[1]**3 + 12, y[2] + 17]), 7)

while True:    
    x0 = np.random.randn(7)
    r = scipy.optimize.leastsq(F, x0)
    err = F(r[0])
    norm =  np.dot(err, err)
    if norm < 1e-6:
        break

print err
HYRY
  • 94,853
  • 25
  • 187
  • 187
  • Thanks, this seems to work. Is there a reason that this is not done automatically? If there are no side effects maybe I will suggest automatic padding as a patch to scipy. – Dan Stahlke Feb 19 '12 at 15:18