Python中的Matlab lsqnonlin答案

【问题标题】：Matlab lsqnonlin in PythonPython中的Matlab lsqnonlin
【发布时间】：2016-09-09 17:30:11
【问题描述】：

我一直在研究 Hartley 和 Zisserman 多视图几何文本，并实现了计算基本矩阵的黄金标准算法。这需要使用 Levenberg-Marquardt 解决非线性最小化问题。

我使用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

这应该是完全矢量化的。三角剖分如上。我可以补充一点，但可能会链接一个要点以保持问题大小可控。

【问题讨论】：

您应该添加两组代码进行比较。有没有可能你写的 MATLAB 代码比你写的 Python 代码更高效？
相关：Why is MATLAB so fast in matrix multiplication?
MATLAB 速度帖子很有趣，但 scipy 代码也应该使用 BLAS/LaPack。上述内容应进行矢量化以提高性能。
函数求值次数相同还是不同？ lstnonlin 和 least_squares 的收敛容差默认值是否相同——您似乎没有在此处指定它们？性能瓶颈甚至是求解器本身，还是目标函数的实现？您应该检查这些以确定问题出在哪里。
我认为默认容差至少不同，least_squares 使用 sqrt(eps)~1e-8 而 lsqnonlin 的默认值是 1e-6 --- 因此，您会期望 least_squares 必须做更多的工作才能达到更严格的容忍度。所以上面的比较肯定是不公平的。

标签： python matlab scipy mathematical-optimization

【解决方案1】：

least_squares 非常新。截至 2015 年秋季，SciPy 领域没有其他选择。否则，例如谷神星。

肯定有很多机会可以加速least_squares --- 很高兴接受拉取请求 :-)。首先要检查的是 SciPy 是否链接到了一个不错的 LAPACK 实现。

【讨论】：

首先，我会确保比较是公平的——默认公差等可能不同。

【解决方案2】：

速度差异巨大的一个原因可能是 matlab 的 lsqnonlin 能够检测雅可比矩阵的稀疏结构，因此计算速度更快。另一方面，scipy 的 minimum_squares 不处理稀疏的 Jacobian 矩阵，而是像在标准情况下一样计算 Jacobian 的每个元素（也是不必要的条目）。

在这个特定的优化问题（确定基本矩阵的黄金标准算法）中，Jacobian 的稀疏计算对于获得 Hartley 和 Zisserman 中提到和描述的良好性能至关重要。（几分钟到>10，取决于用于计算F的点对应的数量）

这是稀疏 LM 实现的 Paython Wrapper 的链接。但是让它运行起来有点复杂：

http://www.bmp.ds.mpg.de/pysparselm.html

我希望 scipy 的 minimum_squares 能够很快升级以处理稀疏的 Jacobians。 :) （scipy 最小二乘的可信区域反射能够处理稀疏雅可比行列式，请查看tutorial。但是，由于 LM 实际上是 MINPACK 实现，您必须在其他地方寻找稀疏 LM 的实现。）

【讨论】：