用Python中的LevenbergMarquardt算法优化方程组

Question

我希望使用scipy的scipy.optimize.leastsq()方法来优化三个参数a,b,c。我有这两个方程。在

1*a+2*b+3*c = x1
4*a+5*b+6*c = x2

从分析上讲，这组方程是不确定的，但在数值上我试图找到a,b,c，以将测量结果的误差最小化[2,2]：

^{pr2}$

所以我写了一些代码。在

def function(a,b,c,t):
    return np.array([1*a+2*b+3*c+t[1],4*a+5*b+6*c+t[1]])

a0 = 1
b0 = 1
c0 = 1
measdata = np.array([2,2])
t = [1,2]

def residual(x0,measdata,t):
    return measdata - function(x0[0],x0[1],x0[2],t)

erg = optimize.leastsq(func=residual,x0=(a0,b0,c0),args=(measdata,t))

它总是导致：

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-296-ab0fc90a2253> in <module>()
     14     return result - function(x0[0],x0[1],x0[2],t)
     15 
---> 16 erg = optimize.leastsq(func = residual, x0 = (a0,b0,c0) , args=(result,t), maxfev=10000)
     17 
     18 function(erg[0][0],erg[0][1])

    //anaconda/lib/python3.5/site-packages/scipy/optimize/minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
        378     m = shape[0]
        379     if n > m:
    --> 380         raise TypeError('Improper input: N=%s must not exceed M=%s' % (n, m))
        381     if epsfcn is None:
        382         epsfcn = finfo(dtype).eps
    TypeError: Improper input: N=3 must not exceed M=2

怎样才能找到最小值？我知道这只是当地的最低要求，但我会很高兴的。在

Answer 1

这个错误告诉你你已经知道的，也就是说，系统是欠定的，其中n是参数的数目，m是约束的数目。在

如果修复其中一个参数，使n > m为False，代码将停止抱怨。例如，改变

def residual(x0,measdata,t):
    return measdata - function(x0[0],x0[1],x0[2],t)

erg = optimize.leastsq(func=residual,x0=(a0,b0,c0),args=(measdata,t))

到

^{pr2}$

为了回答这个问题，你怎么能做你想做的事，我不确定能不能用scipy。我发现这个issue说明scipy不能处理不确定的系统：

Interesting, I assumed the MINPACK routines would handle also m < n, but apparently not. The reason why they don't is probably that for m < n the minimum is some manifold of points, which causes problems in the termination conditions.

So would be after all some interest in adding also a small-scale solver for underdetermined least-square problems.

尽管那篇文章是3年前的，但我仍然在文档中找不到任何证据证明scipy可以做你想做的事。然而，我找到了一个如此的答案，声称you can solve for an underdetermined matrix，但我还没有完全掌握数学，以确定它是否适用于你的情况。因为我发现很难总结这篇文章，所以我只引用看起来最重要的部分。在

In the case where A·x = b is underdetermined,

x1, res, rnk, s = np.linalg.lstsq(A, b)

will pick a solution x' that minimizes ||x||_L2 subject to ||A·x - b||_L2 = 0. This happens not to be the particular solution we are looking for, but we can linearly transform it to get what we want. In order to do that, we'll first compute the right null space of A, which characterizes the space of all possible solutions to A·x = b. We can get this using a rank-revealing QR decomposition