scipy.optimize.leastsq 使用 NaN 调用 objective 函数

scipy.optimize.leastsq calls objective function with NaN

我正在使用 scipy.optimize.leastsq 尝试在存在噪声的情况下将一些参数拟合到真实世界的数据中。 objective 函数偶尔会从 minpack 中用 NaN 调用。这是 scipy.optimize.leastsq 的预期行为吗?在这种情况下,有没有比 returning NaN 残差更好的选择?

以下代码演示了该行为:

import scipy.optimize
import numpy as np

xF = np.array([1.0, 2.0, 3.0, 4.0]) # Target value for fit
NOISE_LEVEL = 1e-6 # The random noise level
RETURN_LEN = 1000  # The objective function return vector length

def func(x):
    if np.isnan(np.sum(x)):
        raise ValueError('Invalid x: %s' % x)
    v = np.random.rand(RETURN_LEN) * NOISE_LEVEL
    v[:len(x)] += xF - x
    return v

iteration = 0
while (1):
    iteration += 1
    x = np.zeros(len(xF))
    y, cov = scipy.optimize.leastsq(func, x)
    print('%04d %s' % (iteration, y))

正在对雅可比行列式进行数值计算。在生产代码中,优化通常有效,除非初始估计太好,表面非常平坦,并且噪声压倒了用于数值计算雅可比行列式的增量。在这种情况下,objective 函数的残差就像上面的代码示例一样显示为随机噪声,我不希望优化收敛。

在此代码示例中,小 NOISE_LEVEL 值 (<1e-10) 始终会收敛。在 1e-6,ValueError 通常会在几百次尝试中抛出。

一种可能的解决方法是 return 高度惩罚的残差(NaN 或 INF),例如:

v = np.empty(RETURN_LEN)
v.fill(np.nan)
return v

如果变脏,此解决方法似乎有效。首先有没有更好的替代方法或方法来防止 NaN?

此行为是在 Windows 7 上的 Python 2.7.9 x32 运行 下观察到的。

由于您的问题定义在非线性问题(噪声)上使用线性求解器,求解器会崩溃,除非噪声水平低于明显阈值。

为了解决这个问题,您可以尝试使用非线性求解器。 例如使用 broyden1 求解算法而不是 leastsq:

import scipy.optimize
import numpy as np

xF = np.array([1.0, 2.0, 3.0, 4.0]) # Target value for fit
NOISE_LEVEL = 1.e-6 # The random noise level
RETURN_LEN = 1000  # The objective function return vector length

def func(x):
    if np.isnan(np.sum(x)):
        raise ValueError('Invalid x: %s' % x)
    v = np.random.rand(RETURN_LEN) * NOISE_LEVEL

    v[:len(x)] += xF - x

    return v[:len(x)]

iteration = 0
while iteration < 10:
    iteration += 1
    x = np.random.rand(len(xF))
    y = scipy.optimize.broyden1(func, x)
    print('%04d %s' % (iteration, y))

returns 作为输出:

0001 [ 1.00000092  2.00000068  3.00000051  4.00000097]
0002 [ 1.0000012   2.00000214  3.00000272  4.00000369]
0003 [ 0.99999991  1.99999931  2.99999815  3.9999978 ]
0004 [ 1.00000097  2.00000198  3.00000345  4.00000425]
0005 [ 1.00000047  1.99999983  2.99999938  3.99999922]
0006 [ 1.00000024  2.00000021  3.00000071  4.00000136]
0007 [ 1.00000116  2.00000102  3.00000225  4.00000357]
0008 [ 1.00000006  2.00000002  3.00000017  4.00000039]
0009 [ 1.0000002   2.00000034  3.00000062  4.00000051]
0010 [ 1.00000137  2.0000015   3.00000193  4.00000344]