我正在使用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下观察到的。