Scipy.optimize.leastsq называет целевую функцию с NaN

Я использую scipy.optimize.leastsq, чтобы попытаться сопоставить ряд параметров с реальными данными при наличии шума. Целевая функция иногда вызывается с NaN изнутри minpack. Это ожидаемое поведение 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 в первую очередь?

Это поведение наблюдалось в Python 2.7.9 x32, работающем в Windows 7.

Ответ 1

Так как в определении проблемы используется линейный решатель по нелинейной задаче (шуму), решатель будет волноваться, если уровень шума не будет ниже порога заметности.

Чтобы решить эту проблему, вы можете попробовать использовать нелинейный решатель. Например, используя алгоритм решения broyden1 вместо lesssq:

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]