numpy - scipy.optimize.leastsq 无法拟合简单模型

Question

我一直在尝试使用Python的scipy.optimize库来估计模型的参数，但到目前为止没有成功。我尝试使用 scipy.optimize.leastsq ，它使用 levenberg-marquardt 算法。不幸的是，即使我将初始参数猜测设置得非常接近最佳拟合，它总是无法找到模型函数的最小值。实际上，它总是返回初始猜测的参数。所以，我认为我做错了什么。我的模型是一个简单的圆，但为了让事情变得更简单，只有半径是实际参数，数据中圆的中心是已知的并且是硬编码的。该数据是一个 10x10 像素的 float 图像，其圆心为 5,5，半径为 4。实际上，数据是使用我试图拟合的模型生成的。因此，存在完美契合。这是我的代码:

import math
import numpy
import scipy.optimize

# ========================================================================

g_data_width = 10
g_data_height = 10
g_xo = 5.0
g_yo = 5.0

def evaluate_model01(x,y,r):
    x2 = x*x
    y2 = y*y
    r2 = r*r
    v = 0.0
    if(x2 + y2 <= r2):
        v = 20.0
    return v

def model01(params,data_o):
    data_r = numpy.zeros(g_data_height*g_data_width)
    r = params[0]
    for y in range(g_data_height):
        for x in range(g_data_width):
            xnew = x - g_xo
            ynew = y - g_yo
            data_r[y*g_data_width+x] = math.fabs(data_o[y,x]-evaluate_model01(xnew,ynew,r))
    return data_r

# ========================================================================

g_data_o = numpy.array([[ 0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
                        [ 0,  0,  0,  0,  0, 20,  0,  0,  0,  0],
                        [ 0,  0,  0, 20, 20, 20, 20, 20,  0,  0],
                        [ 0,  0, 20, 20, 20, 20, 20, 20, 20,  0],
                        [ 0,  0, 20, 20, 20, 20, 20, 20, 20,  0],
                        [ 0, 20, 20, 20, 20, 20, 20, 20, 20, 20],
                        [ 0,  0, 20, 20, 20, 20, 20, 20, 20,  0],
                        [ 0,  0, 20, 20, 20, 20, 20, 20, 20,  0],
                        [ 0,  0,  0, 20, 20, 20, 20, 20,  0,  0],
                        [ 0,  0,  0,  0,  0, 20,  0,  0,  0,  0]],dtype=numpy.float32)
g_params = numpy.array([8.0])
print(scipy.optimize.leastsq(model01,g_params,args=(g_data_o),full_output=1))

# ========================================================================

我对数据进行了硬编码，以消除任何数据依赖性，并允许代码在任何安装了 scipy 的计算机上开箱即用地运行。我不太明白的是 model01 函数应该返回什么。根据文档，它应该返回一个数组。什么的数组？在我的代码中，我假设我必须返回每个数据点的残差数组。那是对的吗？我的数据是一个二维数组，因为它是一个图像，但我的残差是一个扁平的二维残差数组。这可以吗？有人能告诉我我到底做错了什么吗？有人可以修改和修复代码吗？正如我上面提到的，代码应该在任何安装了 scipy 和 numpy 的机器上开箱即用。如果我想要实现的目标无法通过 scipy.optimize.leastsq 实现，您能否建议一些其他适合使用 levenberg-marquardt 算法的模型的库？

Answer 1

恐怕您无法使用leastsq解决您的特定问题。 leastsq 是 FORTRAN 库 minipack 的包装，它调用 MINPACK 的 lmdif 和 lmder算法。重要的是，它基于最小二乘目标函数的雅可比行列式和海森行列式。您的目标函数没有平滑导数，因为:

if(x2 + y2 <= r2):
    v = 20.0

和

math.fabs(......)

，因此leastsq将始终返回初始启动参数。

您应该尝试使用一些不需要梯度/导数的方法，例如 Powell 方法 fmin_powell 或 Nelder-Mead fmin。

关于您的 model101() 中发生的情况。它首先创建一个 g_data 大小为 width*height 的一维数组，称为 data_r。然后它迭代 g_data，计算每个元素的 math.fabs(data_o[y,x]-evaluate_model01(xnew,ynew,r)) 并将该值放入一维数组data_r。最后，它返回date_r。

您的解释是正确的，model101() 返回 2D 数据的展平一维残差数组。