Exemples et solutions qui ne peuvent pas être optimisés avec scipy.optimize.least_squares

Dans scipy, vous pouvez utiliser ʻoptimize.least_squares pour adapter les paramètres d'une fonction non linéaire à vos données. Cependant, selon la forme de la fonction non linéaire, il peut ne pas être possible de trouver les paramètres optimaux. C'est parce que ʻoptimize.least_squares ne peut trouver que des solutions optimales locales.

Cette fois, je vais donner un exemple où ʻoptimize.least_squares tombe dans une solution localement optimale, et utiliser ʻoptimize.basinhopping pour trouver une solution optimale globale.

version:

• Python 3.5.1
• numpy (1.11.1)
• scipy (0.18.0)

Exemples qui ne peuvent pas être bien optimisés

Considérez la fonction suivante avec $ a $ comme paramètre.

y(x)=\frac{1}{100}(x-3a)(2x-a)(3x+a)(x+2a)

Supposons que vous obteniez des données bruyantes lorsque $ a = 2 $.

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('ggplot')

seed = 0
np.random.seed(seed)

def y(x, a):
    return (x-3.*a) * (2.*x-a) * (3.*x+a) * (x+2.*a) / 100.

a_orig = 2.
xs = np.linspace(-5, 7, 1000)
ys = y(xs,a_orig)

num_data = 30
data_x = np.random.uniform(-5, 5, num_data)
data_y = y(data_x, a_orig) + np.random.normal(0, 0.5, num_data)

plt.plot(xs, ys, label='true a = %.2f'%(a_orig))
plt.plot(data_x, data_y, 'o', label='data')
plt.legend()

D'autre part, essayez de trouver les paramètres avec ʻoptimize.least_squares`.

from scipy.optimize import least_squares

def calc_residuals(params, data_x, data_y):
    model_y = y(data_x, params[0])
    return model_y - data_y

a_init = -3
res = least_squares(calc_residuals, np.array([a_init]), args=(data_x, data_y))

a_fit = res.x[0]
ys_fit = y(xs,a_fit)

plt.plot(xs, ys, label='true a = %.2f'%(a_orig))
plt.plot(xs, ys_fit, label='fit a = %.2f'%(a_fit))
plt.plot(data_x, data_y, 'o')
plt.legend()

J'ai défini la valeur initiale du paramètre sur $ a_0 = -3 $, mais cela ne correspondait pas bien aux données.

Pourquoi vous ne pouvez pas bien optimiser

En regardant comment le résultat change en fonction de la valeur initiale du paramètre,

a_inits = np.linspace(-4, 4, 1000)
a_fits = np.zeros(1000)
for i, a_init in enumerate(a_inits):
    res = least_squares(calc_residuals, np.array([a_init]), args=(data_x, data_y))
    a_fits[i] = res.x[0]

plt.plot(a_inits, a_fits)
plt.xlabel("initial value")
plt.ylabel("optimized value")

Si la valeur initiale est négative, vous êtes tombé dans les paramètres localement optimaux. La raison à cela peut être vue en examinant la relation entre les valeurs des paramètres et les résidus. Comme le montre la figure ci-dessous, il existe deux valeurs minimales pour le paramètre, le résultat changera donc en fonction de la valeur initiale.

def calc_cost(params, data_x, data_y):
    residuals = calc_residuals(params, data_x, data_y)
    return (residuals * residuals).sum()

costs = np.zeros(1000)
for i, a in enumerate(a_inits):
    costs[i] = calc_cost(np.array([a]), data_x, data_y)
plt.plot(a_inits, costs)
plt.xlabel("parameter")
plt.ylabel("sum of squares")

Comment bien optimiser

Pour trouver globalement les paramètres optimaux, il suffit de calculer à partir de différentes valeurs initiales. Il y a ʻoptimize.basinhopping` dans scipy pour faire ça bien. Faisons le.

from scipy.optimize import basinhopping
a_init = -3.0
minimizer_kwargs = {"args":(data_x, data_y)}
res = basinhopping(calc_cost, np.array([a_init]),stepsize=2.,minimizer_kwargs=minimizer_kwargs)
print(res.x)

a_fit = res.x[0]
ys_fit = y(xs,a_fit)

plt.plot(xs, ys, label='true a = %.2f'%(a_orig))
plt.plot(xs, ys_fit, label='fit by basin-hopping a = %.2f'%(a_fit))
plt.plot(data_x, data_y, 'o')
plt.legend()

Les paramètres ont été obtenus avec succès. L'astuce est l'argument «stepize». Détermine dans quelle mesure cet argument modifie la valeur initiale.