SciPy (diucapkan sai pie) adalah paket aplikasi matematika berdasarkan ekstensi Numpy Python. Dengan SciPy, sesi Python interaktif Anda menjadi ilmu data lengkap dan lingkungan pembuatan prototipe sistem yang kompleks seperti MATLAB, IDL, Octave, R-Lab, dan SciLab. Hari ini saya ingin berbicara singkat tentang cara menggunakan beberapa algoritma optimasi terkenal dalam paket scipy.optimize. Bantuan yang lebih detail dan terkini dalam menggunakan fungsi selalu dapat diperoleh dengan menggunakan perintah help() atau menggunakan Shift+Tab.
pengenalan
Untuk menyelamatkan diri Anda dan pembaca dari mencari dan membaca sumber utama, tautan ke deskripsi metode terutama ada di Wikipedia. Biasanya, informasi ini cukup untuk memahami metode secara umum dan ketentuan penerapannya. Untuk memahami inti dari metode matematika, ikuti tautan ke publikasi yang lebih otoritatif, yang dapat ditemukan di akhir setiap artikel atau di mesin pencari favorit Anda.
Jadi, modul scipy.optimize mencakup implementasi prosedur berikut:
- Minimisasi fungsi skalar beberapa variabel (minim) bersyarat dan tidak bersyarat menggunakan berbagai algoritma (Nelder-Mead simplex, BFGS, Newton conjugate gradien, ΠΈ )
- Pengoptimalan global (misalnya: , )
- Meminimalkan residu (least_squares) dan algoritma pemasangan kurva menggunakan kuadrat terkecil nonlinier (curve_fit)
- Meminimalkan fungsi skalar satu variabel (minim_scalar) dan mencari akar (root_scalar)
- Pemecah sistem persamaan multidimensi (root) menggunakan berbagai algoritma (hybrid Powell, atau metode skala besar seperti ).
Pada artikel ini kami hanya akan mempertimbangkan item pertama dari keseluruhan daftar ini.
Minimisasi tanpa syarat dari fungsi skalar beberapa variabel
Fungsi minim dari paket scipy.optimize menyediakan antarmuka umum untuk menyelesaikan masalah minimalisasi bersyarat dan tidak bersyarat dari fungsi skalar beberapa variabel. Untuk mendemonstrasikan cara kerjanya, kita memerlukan fungsi yang sesuai dari beberapa variabel, yang akan kita minimalkan dengan cara berbeda.
Untuk tujuan ini, fungsi Rosenbrock dari N variabel sempurna, yang berbentuk:
Meskipun fungsi Rosenbrock dan matriks Jacobi dan Hessiannya (turunan pertama dan kedua) sudah ditentukan dalam paket scipy.optimize, kami akan mendefinisikannya sendiri.
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)Untuk lebih jelasnya, mari kita gambarkan dalam 3D nilai fungsi Rosenbrock dari dua variabel.
Kode gambar
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()
Mengetahui sebelumnya bahwa minimumnya adalah 0 at , mari kita lihat contoh cara menentukan nilai minimum fungsi Rosenbrock menggunakan berbagai prosedur scipy.optimize.
Metode simpleks Nelder-Mead
Misalkan ada titik awal x0 dalam ruang 5 dimensi. Mari kita cari titik minimum dari fungsi Rosenbrock yang paling dekat dengannya menggunakan algoritma (algoritme ditentukan sebagai nilai parameter metode):
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.]Metode simpleks adalah cara paling sederhana untuk meminimalkan fungsi yang terdefinisi secara eksplisit dan cukup lancar. Tidak perlu menghitung turunan suatu fungsi; cukup menentukan nilainya saja. Metode Nelder-Mead adalah pilihan yang baik untuk permasalahan minimalisasi sederhana. Namun, karena tidak menggunakan perkiraan gradien, mungkin diperlukan waktu lebih lama untuk menemukan nilai minimumnya.
metode Powell
Algoritma optimasi lain yang hanya menghitung nilai fungsi adalah . Untuk menggunakannya, Anda perlu menyetel method = 'powell' di fungsi minim.
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.]Algoritma Broyden-Fletcher-Goldfarb-Shanno (BFGS).
Untuk mendapatkan konvergensi yang lebih cepat terhadap suatu solusi, prosedurnya menggunakan gradien fungsi tujuan. Gradien dapat ditentukan sebagai fungsi atau dihitung menggunakan perbedaan orde pertama. Bagaimanapun, metode BFGS biasanya memerlukan lebih sedikit pemanggilan fungsi dibandingkan metode simpleks.
Mari kita cari turunan fungsi Rosenbrock dalam bentuk analitik:
Ekspresi ini berlaku untuk turunan semua variabel kecuali variabel pertama dan terakhir, yang didefinisikan sebagai:
Mari kita lihat fungsi Python yang menghitung gradien ini:
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 derFungsi perhitungan gradien ditentukan sebagai nilai parameter jac dari fungsi minimal, seperti yang ditunjukkan di bawah ini.
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]Algoritma gradien konjugasi (Newton)
Algoritma adalah metode Newton yang dimodifikasi.
Metode Newton didasarkan pada perkiraan suatu fungsi di suatu daerah dengan polinomial derajat kedua:
dimana adalah matriks turunan kedua (matriks Hessian, Hessian).
Jika Hessian adalah definit positif, maka minimum lokal dari fungsi ini dapat dicari dengan menyamakan gradien nol bentuk kuadrat dengan nol. Hasilnya adalah ekspresi:
Invers Hessian dihitung menggunakan metode gradien konjugasi. Contoh penggunaan metode ini untuk meminimalkan fungsi Rosenbrock diberikan di bawah. Untuk menggunakan metode Newton-CG, Anda harus menentukan fungsi yang menghitung Hessian.
Fungsi Hessian Rosenbrock dalam bentuk analitik sama dengan:
dimana ΠΈ , tentukan matriksnya .
Elemen matriks bukan nol yang tersisa sama dengan:
Misalnya, dalam ruang lima dimensi N = 5, matriks Hessian untuk fungsi Rosenbrock berbentuk pita:
Kode yang menghitung Hessian ini beserta kode untuk meminimalkan fungsi Rosenbrock menggunakan metode gradien konjugasi (Newton):
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]Contoh definisi fungsi produk Hessian dan vektor sembarang
Dalam masalah dunia nyata, menghitung dan menyimpan seluruh matriks Hessian memerlukan waktu dan sumber daya memori yang signifikan. Dalam hal ini sebenarnya tidak perlu menentukan matriks Hessian itu sendiri, karena prosedur minimalisasi hanya memerlukan vektor yang sama dengan hasil kali Hessian dengan vektor sembarang lainnya. Jadi, dari sudut pandang komputasi, lebih baik jika segera mendefinisikan fungsi yang mengembalikan hasil perkalian Hessian dengan vektor sembarang.
Pertimbangkan fungsi hess, yang menggunakan vektor minimalisasi sebagai argumen pertama, dan vektor sembarang sebagai argumen kedua (bersama dengan argumen lain dari fungsi yang akan diminimalkan). Dalam kasus kami, menghitung produk fungsi Hessian dari Rosenbrock dengan vektor sembarang tidaklah terlalu sulit. Jika p adalah vektor sembarang, maka hasil kali memiliki bentuk:
Fungsi yang menghitung hasil kali Hessian dan vektor sembarang diteruskan sebagai nilai argumen hessp ke fungsi minimal:
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: 66Algoritma wilayah kepercayaan gradien konjugasi (Newton)
Pengondisian matriks Hessian yang buruk dan arah pencarian yang salah dapat menyebabkan algoritma gradien konjugasi Newton menjadi tidak efektif. Dalam kasus seperti itu, preferensi diberikan kepada (wilayah kepercayaan) konjugasi gradien Newton.
Contoh definisi matriks Hessian:
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.]Contoh dengan fungsi produk Hessian dan vektor sembarang:
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.]Metode tipe Krylov
Seperti metode trust-ncg, metode tipe Krylov sangat cocok untuk menyelesaikan masalah berskala besar karena hanya menggunakan produk matriks-vektor. Esensinya adalah untuk memecahkan masalah di wilayah kepercayaan yang dibatasi oleh subruang Krylov yang terpotong. Untuk masalah yang tidak pasti, lebih baik menggunakan metode ini, karena metode ini menggunakan jumlah iterasi nonlinier yang lebih sedikit karena jumlah produk vektor matriks per submasalah lebih sedikit dibandingkan dengan metode trust-ncg. Selain itu, penyelesaian submasalah kuadratik ditemukan lebih akurat dibandingkan menggunakan metode trust-ncg.
Contoh definisi matriks Hessian:
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.]
Contoh dengan fungsi produk Hessian dan vektor sembarang:
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.]
Algoritma untuk solusi perkiraan di wilayah kepercayaan
Semua metode (Newton-CG, trust-ncg dan trust-krylov) sangat cocok untuk memecahkan masalah skala besar (dengan ribuan variabel). Hal ini disebabkan oleh fakta bahwa algoritma gradien konjugasi yang mendasarinya menyiratkan perkiraan penentuan matriks Hessian terbalik. Solusinya ditemukan secara iteratif, tanpa perluasan Hessian secara eksplisit. Karena Anda hanya perlu mendefinisikan fungsi untuk hasil kali vektor Hessian dan vektor sembarang, algoritma ini sangat baik untuk bekerja dengan matriks renggang (diagonal pita). Hal ini memberikan biaya memori yang rendah dan penghematan waktu yang signifikan.
Untuk permasalahan berukuran sedang, biaya penyimpanan dan pemfaktoran Hessian tidak terlalu penting. Artinya, solusi dapat diperoleh dengan iterasi yang lebih sedikit, sehingga menyelesaikan sub-masalah di wilayah kepercayaan dengan hampir tepat. Untuk melakukan hal ini, beberapa persamaan nonlinier diselesaikan secara iteratif untuk setiap submasalah kuadrat. Solusi seperti itu biasanya memerlukan 3 atau 4 dekomposisi Cholesky dari matriks Hessian. Hasilnya, metode ini konvergen dalam iterasi yang lebih sedikit dan memerlukan penghitungan fungsi tujuan yang lebih sedikit dibandingkan metode wilayah kepercayaan lainnya yang diterapkan. Algoritma ini hanya menyiratkan penentuan matriks Hessian lengkap dan tidak mendukung kemampuan untuk menggunakan fungsi produk Hessian dan vektor sembarang.
Contoh dengan minimalisasi fungsi Rosenbrock:
res = minimize(rosen, x0, method='trust-exact',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
res.xOptimization terminated successfully.
Current function value: 0.000000
Iterations: 13
Function evaluations: 14
Gradient evaluations: 13
Hessian evaluations: 14
array([1., 1., 1., 1., 1.])Kami mungkin akan berhenti di situ. Pada artikel selanjutnya saya akan mencoba menceritakan hal-hal paling menarik tentang minimalisasi bersyarat, penerapan minimalisasi dalam menyelesaikan masalah aproksimasi, minimalisasi fungsi satu variabel, minimalisasi arbitrer, dan mencari akar sistem persamaan menggunakan scipy.optimize kemasan.
Sumber:
Sumber: www.habr.com