😴 🚵 🕧 SciPy, optimasi 🛋️ 🐈 🤘🏿

SciPy (diucapkan sai pie) adalah paket aplikasi matematika berdasarkan ekstensi Numpy Python. Dengan SciPy, sesi Python interaktif berubah menjadi lingkungan pemrosesan data dan prototipe lengkap yang sama untuk sistem yang kompleks seperti MATLAB, IDL, Octave, R-Lab dan SciLab. Hari ini saya ingin berbicara singkat tentang cara menggunakan beberapa algoritma optimasi yang terkenal dalam paket scipy.optimize. Bantuan yang lebih terperinci dan terkini tentang penggunaan fungsi selalu dapat diperoleh dengan menggunakan perintah help () atau menggunakan Shift + Tab.

Pendahuluan

Untuk menyelamatkan diri Anda dan pembaca dari mencari dan membaca sumber, tautan ke deskripsi metode akan terutama ke Wikipedia. Sebagai aturan, informasi ini cukup untuk memahami metode secara umum dan kondisi untuk aplikasi mereka. Untuk memahami esensi metode matematika, kami mengikuti 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:

Minimalisasi bersyarat dan tanpa syarat dari fungsi skalar dari beberapa variabel (minimal) menggunakan berbagai algoritma (Nelder-Mead simplex, BFGS, gradien konjugat Newton, COBYLA dan SLSQP )
Optimalisasi global (mis: basinhopping , diff_evolution )
Meminimalkan residual kuadrat terkecil (least_squares) dan algoritma untuk memasang kurva ke kuadrat terkecil non-linear (curve_fit)
Meminimalkan fungsi skalar dari satu variabel (minim_scalar) dan menemukan root (root_scalar)
Pemecah multidimensi dari sistem persamaan (root) menggunakan berbagai algoritma (hybrid Powell, Levenberg-Marquardt atau metode skala besar, seperti Newton-Krylov ).

Dalam artikel ini kami hanya akan mempertimbangkan item pertama dari seluruh daftar ini.

Minimalisasi tanpa syarat dari fungsi skalar dari beberapa variabel

Fungsi minimal dari paket scipy.optimize menyediakan antarmuka umum untuk memecahkan masalah minimalisasi kondisional dan tanpa syarat dari fungsi skalar dari beberapa variabel. Untuk menunjukkan kerjanya, kita akan membutuhkan fungsi yang sesuai dari beberapa variabel, yang akan kita minimalkan dengan cara yang berbeda.

Untuk keperluan ini, fungsi Rosenbrock dari variabel N sempurna, yang memiliki bentuk:

$f \ kiri (\ mathbf {x} \ kanan) = \ sum_ {i = 1} ^ {N-1} [100 \ kiri (x_ {i + 1} -x_ {i} ^ {2} \ kanan) ^ {2} + \ kiri (1-x_ {i} \ kanan) ^ {2}]$

Terlepas dari kenyataan bahwa fungsi Rosenbrock dan matriks Jacobi dan Hessiannya (turunan pertama dan kedua, masing-masing) sudah didefinisikan dalam paket scipy.optimize, kami 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 kejelasan, kami menggambar dalam 3D nilai-nilai fungsi Rosenbrock dari dua variabel.

Kode untuk rendering

 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 minimum adalah 0 untuk $x_i = 1$ , perhatikan contoh cara menentukan nilai minimum fungsi Rosenbrock menggunakan berbagai prosedur scipy.optimize.

Metode Simpleks Nelder-Mead (Nelder-Mead)

Biarkan ada titik awal x0 dalam ruang 5 dimensi. Temukan titik minimum dari fungsi Rosenbrock yang paling dekat dengannya menggunakan algoritma simpleks Nelder-Mead (algoritma 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 termudah untuk meminimalkan fungsi yang jelas dan halus. Itu tidak memerlukan perhitungan turunan dari suatu fungsi, itu cukup untuk hanya menentukan nilainya. Metode Nelder-Mead adalah pilihan yang baik untuk masalah minimisasi sederhana. Namun, karena tidak menggunakan estimasi gradien, mungkin diperlukan waktu lebih lama untuk menemukan minimum.

Metode Powell

Algoritma optimasi lain di mana hanya nilai-nilai fungsi yang dihitung adalah metode Powell . Untuk menggunakannya, Anda perlu mengatur metode = '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-Channo (BFGS)

Untuk mendapatkan konvergensi yang lebih cepat ke solusi, prosedur BFGS menggunakan gradien dari fungsi tujuan. Gradien dapat ditentukan sebagai fungsi atau dihitung menggunakan perbedaan orde pertama. Bagaimanapun, biasanya metode BFGS membutuhkan lebih sedikit panggilan fungsi daripada metode simpleks.

Kami menemukan turunan dari fungsi Rosenbrock dalam bentuk analitis:

$\ frac {\ partial f} {\ partial x_j} = \ jumlah \ limit_ {i = 1} ^ N 200 (x_i - x_ {i-1} ^ 2) (\ delta_ {i, j} - 2x_ {i -1, j}) - 2 (1 - x_ {i-1}) \ delta_ {i-1, j} =$

$= 200 (x_j - x_ {j-1} ^ 2) - 400x_j (x_ {j + 1} - x_j ^ 2) - 2 (1-x_j)$

Ungkapan ini berlaku untuk turunan dari semua variabel kecuali yang pertama dan terakhir, yang didefinisikan sebagai:

$\ frac {\ partial f} {\ partial x_0} = -400 x_0 (x_1 - x_0 ^ 2) - 2 (1 - x_0),$

$\ frac {\ partial f} {\ partial x_ {N-1}} = 200 (x_ {N-1} - x_ {N-2} ^ 2).$

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 der

Fungsi 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]

Conjugate Gradient Algorithm (Newton)

Algoritma gradien konjugat Newton adalah metode Newton yang dimodifikasi.
Metode Newton didasarkan pada perkiraan suatu fungsi di wilayah lokal oleh polinomial tingkat kedua:

$f \ kiri (\ mathbf {x} \ kanan) \ approx f \ kiri (\ mathbf {x} _ {0} \ kanan) + \ nabla f \ kiri (\ mathbf {x} _ {0} \ kanan) \ cdot \ kiri (\ mathbf {x} - \ mathbf {x} _ {0} \ kanan) + \ frac {1} {2} \ kiri (\ mathbf {x} - \ mathbf {x} _ {0} \ kanan) ^ {T} \ mathbf {H} \ kiri (\ mathbf {x} _ {0} \ kanan) \ kiri (\ mathbf {x} - \ mathbf {x} _ {0} \ kanan)$

dimana $\ mathbf {H} \ kiri (\ mathbf {x} _ {0} \ kanan)$ adalah matriks turunan kedua (Hessian matrix, Hessian).
Jika Goni pasti positif, maka minimum lokal dari fungsi ini dapat ditemukan dengan menyamakan gradien nol dari bentuk kuadrat menjadi nol. Hasilnya adalah ekspresi:

$\ mathbf {x} _ {\ textrm {opt}} = \ mathbf {x} _ {0} - \ mathbf {H} ^ {- 1} \ nabla f$

Invers Hessian dihitung menggunakan metode gradien konjugat. Contoh menggunakan metode ini untuk meminimalkan fungsi Rosenbrock diberikan di bawah ini. Untuk menggunakan metode Newton-CG, Anda harus mendefinisikan fungsi yang mengevaluasi Hessian.
Fungsi Goni Rosenbrock dalam bentuk analitis sama dengan:

$H_ {i, j} = \ frac {\ partial ^ 2 f} {\ partial x_i x_j} = 200 (\ delta_ {i, j} - 2x_ {i-1} \ delta {i-1, j} - 400x_i (\ delta_ {i + 1, j} - 2x_i \ delta {i, j}) - 400 \ delta_ {i, j} (x_ {i + 1} - x_i ^ 2) + 2 \ delta_ {i, j } =$

$= (202 + 1200x_i ^ 2 - 400x_ {i + 1}) \ delta_ {i, j} - 400x_i \ delta_ {i + 1, j} - 400x_ {i-1} \ delta_ {i-1, j}$

dimana $i, j \ di \ kiri [1, N-2 \ kanan]$ dan $i, j \ di \ kiri [0, N-1 \ kanan]$ tentukan matriksnya $N \ kali N$ .

Elemen-elemen bukan nol yang tersisa dari matriks sama dengan:

$\ frac {\ partial ^ 2 f} {\ partial x_0 ^ 2} = 1200x_0 ^ 2 - 400x_1 +2$

$\ frac {\ partial ^ 2 f} {\ partial x_0 x_1} = \ frac {\ partial ^ 2 f} {\ partial x_1 x_0} = -400x_0$

$\ frac {\ partial ^ 2 f} {\ partial x_ {N-1} x_ {N-2}} = \ frac {\ partial ^ 2 f} {\ partial x_ {N-2} x_ {N-1 }} = -400x_ {N-2}$

$\ frac {\ partial ^ 2 f} {\ partial x_ {N-1} ^ 2} = 200x$

Misalnya, dalam ruang lima dimensi N = 5, matriks Hessian untuk fungsi Rosenbrock memiliki bentuk pita:

$\ tiny \ mathbf {H} = \ begin {bmatrix} 1200x_ {0} ^ {2} -400x_ {1} +2 & -400x_ {0} & 0 & 0 & 0 \\ -400x_ {0} & 202 + 1200x_ {1} ^ {2} -400x_ {2} & -400x_ {1} & 0 & 0 \\ 0 & -400x_ {1} & 202 + 1200x_ {2} ^ {2} -400x_ {3} & -400x_ {2} & 0 \\ 0 & & -400x_ {2} & 202 + 1200x_ {3} ^ {2} -400x_ {4} & -400x_ {3} \\ 0 & 0 & 0 & 0 & -400x_ { 3} & 200 \ end {bmatrix}$

Kode yang menghitung Goni ini bersama dengan kode untuk meminimalkan fungsi Rosenbrock menggunakan metode gradien konjugat (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 dengan definisi fungsi produk dari Hessian dan vektor sewenang-wenang

Dalam masalah-masalah dunia nyata, komputasi dan penyimpanan seluruh matriks Hessian mungkin membutuhkan sumber daya waktu dan memori yang signifikan. Selain itu, pada kenyataannya, tidak perlu menentukan matriks Hessian itu sendiri, karena prosedur minimalisasi hanya membutuhkan vektor yang sama dengan produk Hessian dengan vektor sewenang-wenang lainnya. Jadi, dari sudut pandang komputasi, jauh lebih disukai untuk segera menentukan fungsi yang mengembalikan hasil produk Hessian dengan vektor sewenang-wenang.

Pertimbangkan fungsi hess, yang menggunakan vektor minimalisasi sebagai argumen pertama, dan vektor arbitrer sebagai argumen kedua (bersama dengan argumen lain dari fungsi yang diperkecil). Dalam kasus kami, tidak terlalu sulit untuk menghitung produk Hessian dari fungsi Rosenbrock dengan vektor sewenang-wenang. Jika p adalah vektor sembarang, maka produk $H (x) \ cdot p$ memiliki bentuk:

$\ mathbf {H} \ kiri (\ mathbf {x} \ kanan) \ mathbf {p} = \ begin {bmatrix} \ kiri (1200x_ {0} ^ {2} -400x_ {1} +2 \ kanan) p_ {0} -400x_ {0} p_ {1} \\ \ vdots \\ -400x_ {i-1} p_ {i-1} + \ kiri (202 + 1200x_ {i} ^ {2} -400x_ {i + 1} \ kanan) p_ {i} -400x_ {i} p_ {i + 1} \\ \ vdots \\ -400x_ {N-2} p_ {N-2} + 200p_ {N-1} \ end {bmatrix }.$

Fungsi yang menghitung produk Hessian dan vektor sewenang-wenang dilewatkan sebagai nilai argumen hessp ke fungsi meminimalkan:

 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

Algoritme wilayah kepercayaan konjugasi gradien (Newton)

Kondisionalitas yang buruk dari matriks Hessian dan arah pencarian yang salah dapat menyebabkan fakta bahwa algoritma gradien konjugat Newton dapat menjadi tidak efisien. Dalam kasus seperti itu, preferensi diberikan pada metode trust-region dari gradien konjugat 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 dari Hessian dan vektor sewenang-wenang:

 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 jenis Krylovsky

Seperti metode trust-ncg, metode tipe Krylovsky sangat cocok untuk memecahkan masalah skala besar, karena mereka hanya menggunakan produk matriks-vektor. Esensi mereka adalah dalam menyelesaikan masalah di bidang rahasia yang dibatasi oleh subruang Krylov yang terpotong. Untuk tugas-tugas yang tidak jelas, lebih baik menggunakan metode ini, karena menggunakan lebih sedikit iterasi non-linier karena lebih sedikit produk vektor-matriks per subtugas, dibandingkan dengan metode trust-ncg. Selain itu, solusi untuk subtugas kuadrat lebih akurat daripada 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 dari Hessian dan vektor sewenang-wenang:

 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 Perkiraan Berbasis Keyakinan

Semua metode (Newton-CG, trust-ncg dan trust-krylov) sangat cocok untuk menyelesaikan tugas skala besar (dengan ribuan variabel). Hal ini disebabkan oleh fakta bahwa algoritma yang mendasari gradien konjugat menyiratkan penentuan perkiraan matriks Hessian terbalik. Solusinya adalah iteratif, tanpa dekomposisi eksplisit dari Hessian. Karena hanya perlu menentukan fungsi untuk produk Hessian dan vektor arbitrer, algoritma ini sangat baik untuk bekerja dengan matriks jarang (tape diagonal). Ini memberikan biaya memori yang rendah dan penghematan waktu yang signifikan.

Dalam masalah ukuran sedang, biaya penyimpanan dan faktorisasi Goni tidak kritis. Ini berarti bahwa solusi dapat diperoleh dalam iterasi yang lebih sedikit, menyelesaikan subtugas area trust hampir secara tepat. Untuk ini, beberapa persamaan nonlinier diselesaikan secara iteratif untuk setiap subtugas kuadratik. Solusi semacam itu biasanya membutuhkan 3 atau 4 dekomposisi dari matriks Holet Hessian. Akibatnya, metode ini konvergen dalam iterasi yang lebih sedikit dan membutuhkan lebih sedikit perhitungan fungsi obyektif daripada metode lain yang diimplementasikan domain kepercayaan. Algoritma ini hanya menyiratkan penentuan matriks Goni lengkap dan tidak mendukung kemampuan untuk menggunakan fungsi produk Goni dan vektor sewenang-wenang.

Contoh meminimalkan fungsi 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.])

Ini, mungkin, tinggal. Pada artikel selanjutnya saya akan mencoba menceritakan yang paling menarik tentang minimalisasi bersyarat, penerapan minimalisasi dalam menyelesaikan masalah perkiraan, meminimalkan fungsi satu variabel, minimizers sewenang-wenang dan menemukan akar sistem persamaan menggunakan paket scipy.optimize.

Sumber: https://docs.scipy.org/doc/scipy/reference/

SciPy, optimasi