[SciPy-Dev] optimize for nonsquare or degenerate systems
Dan Stahlke
dstahlke at gmail.com
Sun Jun 5 13:03:19 EDT 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