scipy.optimize.leastsq wywołuje funkcję celu za pomocą NaN

Question

Używam scipy.optimize.leastsq do próby dopasowania wielu parametrów do rzeczywistych danych w obecności szumu. Funkcja celu jest czasami wywoływana za pomocą NaNs z poziomu minpack. Czy jest to oczekiwane zachowanie scipy.optimize.leastsq? Czy jest jakaś lepsza opcja niż tylko zwracanie reszt NaN w tym stanie?

Poniższy kod demonstruje zachowanie:

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)) 

Jacobiego jest obliczone numerycznie. W kodzie produkcji optymalizacja działa normalnie, z wyjątkiem sytuacji, gdy początkowy szacunek jest zbyt dobry, powierzchnia jest wyjątkowo płaska, a hałas przytłacza delty używane do obliczenia liczbowego jakobianów. W takim przypadku reszty funkcji celu pojawiają się jako losowy szum, jak w powyższym przykładzie kodu, i nie oczekiwałbym, że optymalizacja się zbiegnie.

W tym przykładzie kodu małe wartości NOISE_LEVEL (< 1e-10) zawsze się zbiegają. W 1e-6, ValueError jest generowany zwykle w ciągu kilkuset prób.

Możliwym rozwiązaniem jest powrót bardzo karane szczątkowego (albo NaN lub INF), takie jak:

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

To rozwiązanie wydaje się być skuteczne, jeżeli brudne. Czy w ogóle istnieją lepsze alternatywy lub sposoby zapobiegania NaN?

Takie zachowanie obserwowano pod Python 2.7.9 x32 działa na Windows 7.

Answer 1

Ponieważ definicja problemu wykorzystuje liniowe solver na problem nieliniowej (hałas), solver będzie wariować, chyba że hałas poziom jest poniżej progu zauważalności.

Aby rozwiązać ten problem, można spróbować użyć nieliniowego solver. na przykład za pomocą algorytmu rozwiązywania broyden1 zamiast 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)) 

powraca wyjściowej:

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]