scipy.optimize.leastsq会将NaN传递给目标函数

Question

我正在使用scipy.optimize.leastsq来拟合一组真实数据的多个参数，以适应噪声。但是在minpack中，目标函数偶尔会被NaN调用。这是scipy.optimize.leastsq的预期行为吗？ 如果出现这种情况，除了返回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))

雅可比矩阵正在通过数值计算得出。在实际生产代码中，通常情况下优化都能正常运行，但当起始估计值过于准确、曲面过于平坦且噪音掩盖了用于数值计算雅可比矩阵的增量时，优化将无法收敛，这时目标函数的残差会像上面的代码示例一样出现随机噪声。

在此代码示例中，小于NOISE_LEVEL值（＜1e-10）始终会收敛。而在1e-6时，ValueError通常会在几百次尝试后抛出。

一种可能的解决方法是返回高度受罚的残差（如NaN或INF），例如：

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

如果有脏数据，这种解决方法似乎是有效的。是否有更好的替代方案或方法可以在一开始就防止NaN？

这种行为是在运行于Windows 7的Python 2.7.9 x32下观察到的。

- Matt Liberty

---

似乎噪声和数值雅可比矩阵的组合会使 x 跑到远离区域，导致 NaN。这是优化器无法提供可靠结果的一个强烈迹象。通常需要通过施加额外的约束条件来重新制定优化问题，使其不那么病态。也许使用有约束的优化器（https://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html#constrained-minimization-of-multivariate-scalar-functions-minimize）可能是一个可行的解决方案。此外，放宽终止条件可能会有所帮助。 - Dietrich

Answer 1

由于您所提出的问题定义在非线性问题（噪声）上使用了线性求解器，除非噪声水平低于可感知性阈值，否则求解器会失控。

为了解决这个问题，您可以尝试使用非线性求解器。例如，使用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))

输出结果如下：

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]

- Uvar