[SciPy-Dev] optimize for nonsquare or degenerate systems

Dan Stahlke dstahlke@gmail....
Sun Jun 5 12:03:19 CDT 2011


I need to solve nonlinear equations in which the number of constraints 
is not equal to the number of variables, and in which some of the 
constraints may be redundant.  For example, see the code below which 
finds an isometry (I realize that this is not the easiest way to find an 
isometry, it is just an example to illustrate my point).  Matlab is able 
to handle such systems.  Is there support for this in scipy?  So far I 
have tried newton_krylov and leastsq with no success.  newton_krylov 
requires constraints=variables and leastsq requires constraints>variables.

If I pad my constraints with zeros to make the number of constraints 
match the number of variables, then leastsq seems to work.  Is there a 
reason that this is not done automatically?

I patched the KrylovJacobian function nonlin.py to give the option of 
using pinv for solving the Jacobian equation, and this seems to work for 
me.  I created a feature request ticket for this 
(http://projects.scipy.org/scipy/ticket/1454) but was asked to turn to 
the mailing list first for advice.

Another question somewhat related: is newton_krylov supposed to handle 
functions on complex numbers?  The documentation doesn't say either way 
as far as I can tell.  The solver seems to attempt such equations but 
doesn't converge upon the answer.  Wrapping the function such that two 
real variables map to one complex variable works though.

Thanks,
Dan

<pre>
# Example of a non-square optimization problem: find an isometry

import numpy as np
import scipy.optimize

(M, N) = (10, 4)

def F(x):
     arr = np.array(x).reshape(M, N)
     closure = np.dot(arr.T, arr)
     err = closure - np.eye(N)
     return err

start = np.random.normal(0, 1, N*M)

# Gives error:
# ValueError: expected square matrix, but got shape=(16, 40)
result = scipy.optimize.newton_krylov(F, start, verbose=True)

# Gives error:
# TypeError: Improper input: N=40 must not exceed M=4
result = scipy.optimize.leastsq(F, start)
</pre>


More information about the SciPy-Dev mailing list