Matlab lsqnonlin in Python

Question

Ich habe durch den Geometrietext der Mehrfachansicht von Hartley und Zisserman gearbeitet und den Goldstandardalgorithmus zur Berechnung der Fundamentalmatrix implementiert. Dies erfordert die Lösung eines nichtlinearen Minimierungsproblems mit Levenberg-Marquardt.

Ich implementiert dies mit scipy.optimize.least_squares, aber die Leistung ist Größenordnungen langsamer als ähnliche (z. B. gleiche Funktionalität) Matlab-Code, der lsqnonlin verwendet. In keinem Fall stelle ich die Jacobi oder eine Maske der Spärlichkeit der Jacobi zur Verfügung.

In Bezug auf Rechenzeiten gilt dies für die Palette der verfügbaren scipy Löser. Ich frage mich, ob es eine Alternative gibt, die eine ähnliche Leistung (numerisch & Geschwindigkeit) zu Matlab hat, oder ob der Wechsel zu einem eingepackten, kompilierten Solver notwendig sein wird?

Bearbeiten Sie für den Code Request Kommentar. Ich versuche die Gesamtmenge des eingefügten Codes zu begrenzen.

Matlab:

P2GS = lsqnonlin(@(h)ReprojErrGS(corres1,PF1,corres2,h),PF2); 

function REGS = ReprojErrGS(corres1,PF1,corres2,PF2) 
    %Find estimated 3D point by Triangulation method 
    XwEst = TriangulationGS(corres1,PF1,corres2,PF2); 
    %Reprojection Back to the image 
    x1hat = PF1*XwEst; 
    x1hat = x1hat ./ repmat(x1hat(3,:),3,1); 
    x2hat = PF2*XwEst; 
    x2hat = x2hat ./ repmat(x2hat(3,:),3,1); 
    %Find root mean squared distance error 
    dist = ((corres1 - x1hat).*(corres1 - x1hat)) + ((corres2 - x2hat).* (corres2 - x2hat)); 
    REGS = sqrt(sum(sum(dist))/size(corres1,2));   

Triangulation ist die Standardmethode, um alle Punkte zu iterieren, die Einrichtung von Ax = 0 und die Lösung von SVD verwendet wird.

Python:

# Using 'trf' for performance, swap to 'lm' for levenberg-marquardt 
result = optimize.least_squares(projection_error, p1.ravel(), args=(p, pt.values, pt1.values), method='trf') 
# Inputs are pandas dataframe, hence the .values 

# Triangulate the correspondences 
xw_est = triangulate(pt, pt1, p, p1) 
# SciPy does not like 2d multi-dimensional variables, so reshape 

if p1.shape != (3,4): 
    p1 = p1.reshape(3,4) 

xhat = p.dot(xw_est).T 
xhat /= xhat[:,-1][:,np.newaxis] 
x2hat = p1.dot(xw_est).T 
x2hat /= x2hat[:,-1][:,np.newaxis] 
# Compute error 
dist = (pt - xhat)**2 + (pt1 - x2hat)**2 
reproj_error = np.sqrt(np.sum(dist, axis=1)/len(pt)) 
# print(reproj_error) 
return reproj_error 

Dies sollte vollständig vektorisiert werden. Triangulation ist wie oben. Ich könnte das hinzufügen, aber würde wahrscheinlich einen Kern verknüpfen, um die Frage Größe handhabbar zu halten.

Answer 1

least_squares ist sehr neu. Ab Herbst 2015 gab es im SciPy-Land keine Alternativen mehr. Ansonsten gibt es z.B. Ceres.

Da gibt es sicherlich viele Möglichkeiten zur Beschleunigung least_squares --- Anfragen werden gerne angenommen :-). Die erste Sache zu überprüfen ist, dass SciPy mit einer anständigen LAPACK-Implementierung verbunden ist.