scipy.optimize.least_squares에서 잘 최적화 할 수없는 예제와 그 해결 방법

scipy에서는 optimize.least_squares를 사용하여 비선형 함수의 매개 변수를 데이터에 맞출 수 있습니다. 그러나 비선형 함수의 형태에 따라 최적의 매개 변수를 찾지 못할 수 있습니다. 왜냐하면 optimize.least_squares 에서는 국소적인 최적해 밖에 구할 수 없기 때문입니다.

이번에는 optimize.least_squares 가 국소 최적해에 빠져 버리는 예를 들어, optimize.basinhopping 를 사용하여 전역적 최적해를 구해 보겠습니다.

버전:

Python 3.5.1

numpy (1.11.1)

scipy (0.18.0)

잘 최적화할 수 없는 예

$a$를 매개 변수로 사용하는 다음과 같은 함수를 생각해 봅시다.

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

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

반대로 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()

매개 변수의 초기 값을 $a_0 = -3$로 설정했는데 데이터에 잘 맞지 않았습니다.

잘 최적화할 수 없는 이유

매개변수의 초기값에 따라 결과가 어떻게 달라지는지 살펴보면,

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

초기값이 음수이면 국소적으로 최적의 파라미터에 빠져 버립니다. 이 이유는 매개 변수의 값과 잔차의 관계를 살펴보면 알 수 있습니다. 아래 그림과 같이, 파라미터에 대해서 극소치가 두 개 있는 탓으로, 초기치에 의존해 결과가 바뀌어 버리는 것입니다.

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

잘 최적화하는 방법

전역적으로 최적의 파라미터를 얻으려면 다양한 초기값에서 계산해 보면 된다는 것입니다. 이것을 좋은 느낌으로하는 방법으로 scipy에는 optimize.basinhopping가 있습니다. 그럼 해보자.

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

잘 매개 변수를 구했습니다. 요령은 인수의 stepsize 입니다. 이 인수가 얼마나 초기값을 크게 바꾸는지를 결정합니다.