2016-05-13 12 views
0

私はHartleyとZissermanの多視点ジオメトリテキストを使って作業しており、Fundamental行列を計算するためのゴールド標準アルゴリズムを実装しています。これには、Levenberg-Marquardtを使用して非線形最小化問題を解く必要があります。Matlab lsqnonlin in Python

私はこれをscipy.optimize.least_squaresで実装しましたが、性能はlsqnonlinを使用する同様の(たとえば、同じ機能性の)MATLABコードの次数よりも遅いです。いずれにしても、私はヤコビ行列、またはヤコビ行列の希薄さのマスクを供給していません。

計算時間に関して、これは利用可能なscipyソルバの範囲に当てはまります。 MATLABに同様の性能(数値&の速度)を持つ代替案が存在するのか、ラップされコンパイルされたソルバーに移動する必要があるのだろうか?

コード要求コメントの編集。私は挿入されたコードの総量を制限しようとしています。

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));   

三角測量は、すべての点にわたって反復方程式Ax = 0を設定し、SVDを使用して解決する、標準的な方法です。

の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 

。三角測量は上記のとおりです。私はそれを加えることができますが、疑問の大きさを管理しやすいように要点を結びつけるでしょう。

+2

比較のために、両方のコードセットを追加する必要があります。 Pythonコードを書くよりも効率的なMATLABコードを書くことは可能でしょうか? – Dan

+1

関連:[MATLABはなぜ高速な行列乗算をするのですか?](0120-18753) – Adriaan

+0

MATLAB速度のポストは次のとおりです。 scipyコードはBLAS/LaPackを使用する必要があります。上記はパフォーマンスのためにベクトル化する必要があります。 – Jzl5325

答えて

0

least_squaresは非常に新しいです。 2015年秋、SciPyの土地には選択肢がありませんでした。それ以外の場合は、セレス。

確かに多くの機会があります。least_squares ---プルリクエストは喜んで受け入れられています:-)。最初にチェックするのは、SciPyがまともなLAPACKの実装にリンクされていることです。

+0

まず、比較が公正であることを確認します---デフォルトの許容誤差などはおそらく異なります。 –