SciPy optimization for under-constrained system

Question

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.

Answer 1

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