Ottimizzazione di Scipy per il sistema non vincolato

Question

Spesso devo risolvere problemi non lineari in cui il numero di variabili supera il numero di vincoli (o talvolta il contrario). Di solito alcuni dei vincoli o delle variabili sono ridondanti in modo complicato. C'è un modo per risolvere tali problemi?

La maggior parte dei solutori di Scipy sembra supporre che il numero di vincoli sia uguale al numero di variabili e che il giacobino sia non singolare. leastsq Funziona a volte ma non prova nemmeno quando i vincoli sono meno del numero di variabili. Mi rendo conto che potrei semplicemente correre fmin Su linalg.norm(F), ma questo è molto meno efficiente di qualsiasi metodo che utilizza il giacobino.

Ecco un esempio di un problema che dimostra ciò di cui sto parlando. Ovviamente ha una soluzione, ma leastsq dà un errore. Naturalmente, questo esempio è facile da risolvere a mano, l'ho appena messo qui per dimostrare il problema.

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)

Il messaggio di errore che ricevo è:

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

Ho perlustrato la rete per una risposta e ho persino chiesto nella mailing list di Scipy e non ho ricevuto risposta. Per ora ho hackerato la fonte Scipy in modo che il newton_krylov usi del risolutore pinv(), ma non penso che questa sia una soluzione ottimale.

Answer 1

Che ne dici di ridimensionare l'array di restituzione da f () al numero di variabili:

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