在 Python 中使用 Levenberg-Marquardt 算法优化方程组

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] 给定结果的误差：

1*a+2*b+3*c - 2 = 0
4*a+5*b+6*c - 2 = 0

因此我写了一些代码。

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

到

def residual(x0,measdata,t):
    # we fix the value of `c` here
    return measdata - function(x0[0],x0[1],5,t)

# only two parameters for `x0`
erg = optimize.leastsq(func=residual,x0=(a0,b0),args=(measdata,t))

要回答如何做自己想做的事情，我不确定 scipy 是否可以做到。我发现这个 issue 说 scipy 无法处理不确定的系统：

有趣的是，我认为MINPACK 例程也可以处理m

因此，对于未定的最小二乘问题也添加一个小规模求解器，毕竟会有一些兴趣。

即使那篇文章是 3 年前的，我仍然无法在文档中找到任何证据表明 scipy 可以做你想做的事。但是，我找到了一个声称 you can solve for an underdetermined matrix 的 SO 答案，但我还没有完全掌握数学知识，以确定它是否适用于您的案例。因为我觉得很难总结这篇文章，所以我只引用似乎最重要的部分。

在A·x = b不确定的情况下，

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

将选择一个解决方案 x' 最小化 ||x||_L2时间> 服从 ||A·x - b||_L2 = 0。有时候是这样的 不是我们正在寻找的特定解决方案，但我们可以 线性变换它以获得我们想要的。为了做到这一点，我们将 首先计算 A 的right null space，即 表征 A·x 的所有可能解的空间 = b。我们可以使用rank-revealing QR decomposition