SciPy, 최적화

SciPy(sai Pie로 발음)는 Numpy Python 확장을 기반으로 하는 수학 애플리케이션 패키지입니다. SciPy를 사용하면 대화형 Python 세션이 MATLAB, IDL, Octave, R-Lab 및 SciLab과 동일한 완전한 데이터 과학 및 복잡한 시스템 프로토타입 제작 환경이 됩니다. 오늘 저는 scipy.optimize 패키지에서 잘 알려진 최적화 알고리즘을 사용하는 방법에 대해 간략하게 이야기하고 싶습니다. 기능 사용에 대한 더 자세한 최신 도움말은 언제든지 help() 명령을 사용하거나 Shift+Tab을 사용하여 얻을 수 있습니다.

소개

자신과 독자가 주요 소스를 검색하고 읽는 것을 방지하기 위해 방법 설명에 대한 링크는 주로 Wikipedia에 있습니다. 일반적으로 이 정보는 일반적인 용어 및 적용 조건을 이해하는 데 충분합니다. 수학적 방법의 본질을 이해하려면 각 기사의 끝이나 즐겨 사용하는 검색 엔진에서 찾을 수 있는 보다 권위 있는 출판물의 링크를 따르십시오.

따라서 scipy.optimize 모듈에는 다음 절차의 구현이 포함됩니다.

1. 다양한 알고리즘(Nelder-Mead Simplex, BFGS, Newton Conjugate Gradients, 코빌라 и SLSQP)
2. 전역 최적화(예: 유역 호핑, diff_evolution)
3. 잔차 최소화 MNC (least_squares) 및 비선형 최소 제곱을 사용한 곡선 피팅 알고리즘(curve_fit)
4. 한 변수(minim_scalar)의 스칼라 함수 최소화 및 루트 검색(root_scalar)
5. 다양한 알고리즘(하이브리드 Powell, Levenberg-Marquardt 또는 다음과 같은 대규모 방법 뉴턴-크릴로프).

이 기사에서는 전체 목록 중 첫 번째 항목만 고려할 것입니다.

여러 변수의 스칼라 함수를 무조건적으로 최소화

scipy.optimize 패키지의 minim 함수는 여러 변수의 스칼라 함수의 조건부 및 무조건적 최소화 문제를 해결하기 위한 일반 인터페이스를 제공합니다. 그것이 어떻게 작동하는지 보여주기 위해 우리는 여러 변수의 적절한 함수가 필요하며 이를 다양한 방법으로 최소화할 것입니다.

이러한 목적을 위해 N 변수의 Rosenbrock 함수는 완벽하며 다음과 같은 형식을 갖습니다.

Rosenbrock 함수와 해당 Jacobi 및 Hessian 행렬(각각 XNUMX차 및 XNUMX차 도함수)이 scipy.optimize 패키지에 이미 정의되어 있다는 사실에도 불구하고 우리는 이를 직접 정의하겠습니다.

import numpy as np

def rosen(x):
    """The Rosenbrock function"""
    return np.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0, axis=0)

명확성을 위해 두 변수의 Rosenbrock 함수 값을 3D로 그려 보겠습니다.

도면 코드

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter

# Настраиваем 3D график
fig = plt.figure(figsize=[15, 10])
ax = fig.gca(projection='3d')

# Задаем угол обзора
ax.view_init(45, 30)

# Создаем данные для графика
X = np.arange(-2, 2, 0.1)
Y = np.arange(-1, 3, 0.1)
X, Y = np.meshgrid(X, Y)
Z = rosen(np.array([X,Y]))

# Рисуем поверхность
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
plt.show()

최소값이 0이라는 것을 미리 알고 , 다양한 scipy.optimize 절차를 사용하여 Rosenbrock 함수의 최소값을 결정하는 방법에 대한 예를 살펴보겠습니다.

Nelder-Mead 단순법

0차원 공간에 초기점 x5이 있다고 가정합니다. 알고리즘을 이용하여 가장 가까운 Rosenbrock 함수의 최소점을 찾아보겠습니다. Nelder-Mead 단순형 (알고리즘은 메소드 매개변수의 값으로 지정됩니다):

from scipy.optimize import minimize
x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = minimize(rosen, x0, method='nelder-mead',
    options={'xtol': 1e-8, 'disp': True})
print(res.x)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 339
         Function evaluations: 571
[1. 1. 1. 1. 1.]

단순 방법은 명시적으로 정의되고 상당히 매끄러운 함수를 최소화하는 가장 간단한 방법입니다. 함수의 도함수를 계산할 필요가 없으며 해당 값만 지정하면 충분합니다. Nelder-Mead 방법은 간단한 최소화 문제에 적합한 선택입니다. 그러나 기울기 추정을 사용하지 않기 때문에 최소값을 찾는 데 시간이 더 오래 걸릴 수 있습니다.

파월 방법

함수값만 계산하는 또 다른 최적화 알고리즘은 파웰의 방법. 이를 사용하려면 minim 함수에서 method = 'powell'을 설정해야 합니다.

x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = minimize(rosen, x0, method='powell',
    options={'xtol': 1e-8, 'disp': True})
print(res.x)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 19
         Function evaluations: 1622
[1. 1. 1. 1. 1.]

Broyden-Fletcher-Goldfarb-Shanno(BFGS) 알고리즘

솔루션으로 더 빠르게 수렴하려면 다음 절차를 따르세요. BFGS 목적 함수의 기울기를 사용합니다. 기울기는 함수로 지정되거나 XNUMX차 차이를 사용하여 계산될 수 있습니다. 어떤 경우든 BFGS 메서드는 일반적으로 단순 메서드보다 더 적은 수의 함수 호출을 필요로 합니다.

분석적 형태로 Rosenbrock 함수의 미분을 찾아보겠습니다.

이 표현식은 다음과 같이 정의된 첫 번째 변수와 마지막 변수를 제외한 모든 변수의 파생물에 유효합니다.

이 기울기를 계산하는 Python 함수를 살펴보겠습니다.

def rosen_der (x):
    xm = x [1: -1]
    xm_m1 = x [: - 2]
    xm_p1 = x [2:]
    der = np.zeros_like (x)
    der [1: -1] = 200 * (xm-xm_m1 ** 2) - 400 * (xm_p1 - xm ** 2) * xm - 2 * (1-xm)
    der [0] = -400 * x [0] * (x [1] -x [0] ** 2) - 2 * (1-x [0])
    der [-1] = 200 * (x [-1] -x [-2] ** 2)
    return der

기울기 계산 함수는 아래와 같이 minim 함수의 jac 매개변수 값으로 지정됩니다.

res = minimize(rosen, x0, method='BFGS', jac=rosen_der, options={'disp': True})
print(res.x)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 25
         Function evaluations: 30
         Gradient evaluations: 30
[1.00000004 1.0000001  1.00000021 1.00000044 1.00000092]

켤레 기울기 알고리즘(뉴턴)

연산 뉴턴의 공액 기울기 수정된 뉴턴의 방법이다.
뉴턴의 방법은 XNUMX차 다항식으로 로컬 영역의 함수를 근사하는 데 기반을 두고 있습니다.

어디에서 는 XNUMX차 도함수 행렬(Hessian 행렬, Hessian)입니다.
헤세 행렬이 양의 정부호인 경우 이 함수의 국소 최소값은 XNUMX차 형식의 XNUMX 기울기를 XNUMX과 동일시하여 찾을 수 있습니다. 결과는 다음과 같습니다.

역 헤세 행렬은 켤레 기울기 방법을 사용하여 계산됩니다. Rosenbrock 함수를 최소화하기 위해 이 방법을 사용하는 예는 다음과 같습니다. Newton-CG 방법을 사용하려면 헤세 행렬을 계산하는 함수를 지정해야 합니다.
분석 형식의 Rosenbrock 함수의 헤세 행렬은 다음과 같습니다.

어디에서 и , 행렬을 정의 .

행렬의 XNUMX이 아닌 나머지 요소는 다음과 같습니다.

예를 들어, 5차원 공간 N = XNUMX에서 Rosenbrock 함수에 대한 헤세 행렬은 밴드 형태를 갖습니다.

켤레 기울기(뉴턴) 방법을 사용하여 Rosenbrock 함수를 최소화하는 코드와 함께 이 헤세 행렬을 계산하는 코드:

def rosen_hess(x):
    x = np.asarray(x)
    H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1)
    diagonal = np.zeros_like(x)
    diagonal[0] = 1200*x[0]**2-400*x[1]+2
    diagonal[-1] = 200
    diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]
    H = H + np.diag(diagonal)
    return H

res = minimize(rosen, x0, method='Newton-CG', 
               jac=rosen_der, hess=rosen_hess,
               options={'xtol': 1e-8, 'disp': True})
print(res.x)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 24
         Function evaluations: 33
         Gradient evaluations: 56
         Hessian evaluations: 24
[1.         1.         1.         0.99999999 0.99999999]

헤세 행렬의 곱 함수 정의와 임의 벡터의 예

실제 문제에서 전체 헤세 행렬을 계산하고 저장하려면 상당한 시간과 메모리 리소스가 필요할 수 있습니다. 이 경우 실제로 헤세 행렬 자체를 지정할 필요가 없습니다. 최소화 절차에는 헤세 행렬과 다른 임의 벡터의 곱과 동일한 벡터만 필요합니다. 따라서 계산적인 관점에서는 헤세 행렬과 임의의 벡터를 곱한 결과를 반환하는 함수를 즉시 정의하는 것이 훨씬 바람직합니다.

최소화 벡터를 첫 번째 인수로, 임의의 벡터를 두 번째 인수로(최소화할 함수의 다른 인수와 함께) 사용하는 hess 함수를 생각해 보세요. 우리의 경우 임의의 벡터를 사용하여 Rosenbrock 함수의 헤세 행렬의 곱을 계산하는 것은 그리 어렵지 않습니다. 만약에 p 는 임의의 벡터이면 곱은 형식은 다음과 같습니다.

헤세 행렬과 임의 벡터의 곱을 계산하는 함수는 hessp 인수의 값으로 최소화 함수에 전달됩니다.

def rosen_hess_p(x, p):
    x = np.asarray(x)
    Hp = np.zeros_like(x)
    Hp[0] = (1200*x[0]**2 - 400*x[1] + 2)*p[0] - 400*x[0]*p[1]
    Hp[1:-1] = -400*x[:-2]*p[:-2]+(202+1200*x[1:-1]**2-400*x[2:])*p[1:-1] 
    -400*x[1:-1]*p[2:]
    Hp[-1] = -400*x[-2]*p[-2] + 200*p[-1]
    return Hp

res = minimize(rosen, x0, method='Newton-CG',
               jac=rosen_der, hessp=rosen_hess_p,
               options={'xtol': 1e-8, 'disp': True})

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 24
         Function evaluations: 33
         Gradient evaluations: 56
         Hessian evaluations: 66

켤레 경사 신뢰 영역 알고리즘(뉴턴)

헤세 행렬의 조건이 좋지 않고 검색 방향이 올바르지 않으면 뉴턴의 켤레 기울기 알고리즘이 효과적이지 않을 수 있습니다. 그러한 경우에는 다음이 우선적으로 제공됩니다. 신뢰 영역 방법 (trust-region) 뉴턴 기울기를 공액화합니다.

헤세 행렬 정의의 예:

res = minimize(rosen, x0, method='trust-ncg',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})
print(res.x)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 20
         Function evaluations: 21
         Gradient evaluations: 20
         Hessian evaluations: 19
[1. 1. 1. 1. 1.]

헤세 행렬과 임의 벡터의 곱 함수의 예:

res = minimize(rosen, x0, method='trust-ncg', 
                jac=rosen_der, hessp=rosen_hess_p, 
                options={'gtol': 1e-8, 'disp': True})
print(res.x)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 20
         Function evaluations: 21
         Gradient evaluations: 20
         Hessian evaluations: 0
[1. 1. 1. 1. 1.]

Krylov 유형 방법

trust-ncg 방법과 마찬가지로 Krylov 유형 방법은 행렬-벡터 곱만 사용하므로 대규모 문제를 해결하는 데 매우 적합합니다. 그들의 본질은 잘린 Krylov 부분 공간으로 제한되는 신뢰 영역의 문제를 해결하는 것입니다. 불확실한 문제의 경우에는 trust-ncg 방법에 비해 하위 문제당 행렬-벡터 곱의 수가 적기 때문에 더 적은 수의 비선형 반복을 사용하므로 이 방법을 사용하는 것이 더 좋습니다. 또한, trust-ncg 방법을 사용하는 것보다 XNUMX차 하위 문제에 대한 해를 더 정확하게 구합니다.
헤세 행렬 정의의 예:

res = minimize(rosen, x0, method='trust-krylov',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 19
         Function evaluations: 20
         Gradient evaluations: 20
         Hessian evaluations: 18

print(res.x)

    [1. 1. 1. 1. 1.]

헤세 행렬과 임의 벡터의 곱 함수의 예:

res = minimize(rosen, x0, method='trust-krylov',
               jac=rosen_der, hessp=rosen_hess_p,
               options={'gtol': 1e-8, 'disp': True})

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 19
         Function evaluations: 20
         Gradient evaluations: 20
         Hessian evaluations: 0

print(res.x)

    [1. 1. 1. 1. 1.]

신뢰 영역의 근사해를 위한 알고리즘

모든 방법(Newton-CG, trust-ncg 및 trust-krylov)은 수천 개의 변수가 있는 대규모 문제를 해결하는 데 매우 적합합니다. 이는 기본 켤레 기울기 알고리즘이 역 헤시안 행렬의 대략적인 결정을 의미하기 때문입니다. 헤세 행렬을 명시적으로 확장하지 않고 반복적으로 해를 구합니다. 헤세 행렬과 임의 벡터의 곱에 대한 함수만 정의하면 되므로 이 알고리즘은 희소(대각선) 행렬 작업에 특히 좋습니다. 이는 낮은 메모리 비용과 상당한 시간 절약을 제공합니다.

중간 크기 문제의 경우 헤세 행렬을 저장하고 인수분해하는 비용은 중요하지 않습니다. 이는 더 적은 반복으로 솔루션을 얻을 수 있으며 신뢰 영역의 하위 문제를 거의 정확하게 해결할 수 있음을 의미합니다. 이를 위해 일부 비선형 방정식은 각 3차 하위 문제에 대해 반복적으로 해결됩니다. 이러한 솔루션에는 일반적으로 헤세 행렬의 촐레스키 분해가 4~XNUMX회 필요합니다. 결과적으로 이 방법은 구현된 다른 신뢰 영역 방법보다 더 적은 반복 횟수로 수렴하고 더 적은 목적 함수 계산이 필요합니다. 이 알고리즘은 완전한 헤세 행렬의 결정만을 의미하며 헤세 행렬과 임의 벡터의 곱 함수를 사용하는 기능은 지원하지 않습니다.

Rosenbrock 함수를 최소화한 예:

res = minimize(rosen, x0, method='trust-exact',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})
res.x

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 13
         Function evaluations: 14
         Gradient evaluations: 13
         Hessian evaluations: 14

array([1., 1., 1., 1., 1.])

우리는 아마 거기서 멈출 것입니다. 다음 기사에서는 조건부 최소화, 근사 문제 해결 시 최소화 적용, 한 변수의 함수 최소화, 임의 최소화, scipy.optimize를 사용하여 방정식 시스템의 근 찾기에 대해 가장 흥미로운 내용을 설명하려고 합니다. 패키지.

출처 : https://docs.scipy.org/doc/scipy/reference/

출처 : habr.com