Το SciPy (προφέρεται πίτα sai) είναι ένα πακέτο μαθηματικών εφαρμογών που βασίζεται στην επέκταση Numpy Python. Με το SciPy, η διαδραστική συνεδρία Python γίνεται το ίδιο πλήρες περιβάλλον επιστήμης δεδομένων και σύνθετων πρωτοτύπων συστημάτων όπως τα MATLAB, IDL, Octave, R-Lab και SciLab. Σήμερα θέλω να μιλήσω εν συντομία για τον τρόπο χρήσης μερικών γνωστών αλγορίθμων βελτιστοποίησης στο πακέτο scipy.optimize. Μπορείτε πάντα να λάβετε πιο λεπτομερή και ενημερωμένη βοήθεια σχετικά με τη χρήση συναρτήσεων χρησιμοποιώντας την εντολή help() ή χρησιμοποιώντας το Shift+Tab.
Εισαγωγή
Για να γλιτώσετε τον εαυτό σας και τους αναγνώστες από την αναζήτηση και την ανάγνωση πρωτογενών πηγών, οι σύνδεσμοι προς τις περιγραφές των μεθόδων θα βρίσκονται κυρίως στη Wikipedia. Κατά κανόνα, αυτές οι πληροφορίες είναι επαρκείς για την κατανόηση των μεθόδων σε γενικούς όρους και των προϋποθέσεων εφαρμογής τους. Για να κατανοήσετε την ουσία των μαθηματικών μεθόδων, ακολουθήστε τους συνδέσμους προς πιο έγκυρες δημοσιεύσεις, που μπορείτε να βρείτε στο τέλος κάθε άρθρου ή στην αγαπημένη σας μηχανή αναζήτησης.
Έτσι, η ενότητα scipy.optimize περιλαμβάνει την υλοποίηση των παρακάτω διαδικασιών:
- Υπό όρους και άνευ όρων ελαχιστοποίηση βαθμωτών συναρτήσεων πολλών μεταβλητών (ελάχιστη) με χρήση διαφόρων αλγορίθμων (Nelder-Mead simplex, BFGS, συζυγείς διαβαθμίσεις Newton, и )
- Παγκόσμια βελτιστοποίηση (για παράδειγμα: , )
- Ελαχιστοποίηση υπολειμμάτων (ελάχιστα τετράγωνα) και αλγόριθμοι προσαρμογής καμπύλης που χρησιμοποιούν μη γραμμικά ελάχιστα τετράγωνα (καμπύλη_προσαρμογή)
- Ελαχιστοποίηση βαθμωτών συναρτήσεων μιας μεταβλητής (minim_scalar) και αναζήτηση ριζών (root_scalar)
- Πολυδιάστατοι λύτες συστήματος εξισώσεων (ρίζα) χρησιμοποιώντας διάφορους αλγόριθμους (υβριδικό Powell, ή μεθόδους μεγάλης κλίμακας όπως ).
Σε αυτό το άρθρο θα εξετάσουμε μόνο το πρώτο στοιχείο από ολόκληρη τη λίστα.
Ελαχιστοποίηση άνευ όρων μιας κλιμακωτής συνάρτησης πολλών μεταβλητών
Η συνάρτηση minim από το πακέτο scipy.optimize παρέχει μια γενική διεπαφή για την επίλυση προβλημάτων ελαχιστοποίησης υπό όρους και άνευ όρων βαθμωτών συναρτήσεων πολλών μεταβλητών. Για να δείξουμε πώς λειτουργεί, θα χρειαστούμε μια κατάλληλη συνάρτηση πολλών μεταβλητών, την οποία θα ελαχιστοποιήσουμε με διαφορετικούς τρόπους.
Για αυτούς τους σκοπούς, η συνάρτηση Rosenbrock των N μεταβλητών είναι τέλεια, η οποία έχει τη μορφή:
Παρά το γεγονός ότι η συνάρτηση Rosenbrock και οι πίνακες Jacobi και Hessian (η πρώτη και η δεύτερη παράγωγος, αντίστοιχα) έχουν ήδη οριστεί στο πακέτο 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 δύο μεταβλητών.
Κώδικας σχεδίασης
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 στο , ας δούμε παραδείγματα για τον προσδιορισμό της ελάχιστης τιμής της συνάρτησης Rosenbrock χρησιμοποιώντας διάφορες διαδικασίες scipy.optimize.
Μέθοδος Simplex Nelder-Mead
Έστω ένα αρχικό σημείο x0 στον 5-διάστατο χώρο. Ας βρούμε το ελάχιστο σημείο της συνάρτησης Rosenbrock που βρίσκεται πιο κοντά σε αυτό χρησιμοποιώντας τον αλγόριθμο (ο αλγόριθμος καθορίζεται ως η τιμή της παραμέτρου της μεθόδου):
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.]Η μέθοδος simplex είναι ο απλούστερος τρόπος για την ελαχιστοποίηση μιας ρητά καθορισμένης και αρκετά ομαλής συνάρτησης. Δεν απαιτεί τον υπολογισμό των παραγώγων μιας συνάρτησης, αρκεί να καθορίσετε μόνο τις τιμές της. Η μέθοδος Nelder-Mead είναι μια καλή επιλογή για απλά προβλήματα ελαχιστοποίησης. Ωστόσο, δεδομένου ότι δεν χρησιμοποιεί εκτιμήσεις κλίσης, μπορεί να χρειαστεί περισσότερος χρόνος για να βρεθεί το ελάχιστο.
Μέθοδος Πάουελ
Ένας άλλος αλγόριθμος βελτιστοποίησης στον οποίο υπολογίζονται μόνο οι τιμές της συνάρτησης είναι . Για να το χρησιμοποιήσετε, πρέπει να ορίσετε τη μέθοδο = '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 απαιτεί συνήθως λιγότερες κλήσεις συναρτήσεων από τη μέθοδο simplex.
Ας βρούμε την παράγωγο της συνάρτησης 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Η συνάρτηση υπολογισμού κλίσης καθορίζεται ως η τιμή της παραμέτρου 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]Συζυγής αλγόριθμος κλίσης (Newton)
Αλγόριθμος είναι μια τροποποιημένη μέθοδος του Νεύτωνα.
Η μέθοδος του Νεύτωνα βασίζεται στην προσέγγιση μιας συνάρτησης σε μια τοπική περιοχή με ένα πολυώνυμο δεύτερου βαθμού:
όπου είναι ο πίνακας των δεύτερων παραγώγων (Hessian matrix, Hessian).
Εάν η Έσσια είναι θετική οριστική, τότε το τοπικό ελάχιστο αυτής της συνάρτησης μπορεί να βρεθεί εξισώνοντας τη μηδενική κλίση της τετραγωνικής μορφής με μηδέν. Το αποτέλεσμα θα είναι η έκφραση:
Ο αντίστροφος Έσσιος υπολογίζεται με τη μέθοδο της συζυγούς κλίσης. Ένα παράδειγμα χρήσης αυτής της μεθόδου για την ελαχιστοποίηση της συνάρτησης Rosenbrock δίνεται παρακάτω. Για να χρησιμοποιήσετε τη μέθοδο Newton-CG, πρέπει να καθορίσετε μια συνάρτηση που υπολογίζει την Hessian.
Το Hessian της συνάρτησης Rosenbrock σε αναλυτική μορφή ισούται με:
όπου и , ορίστε τον πίνακα .
Τα υπόλοιπα μη μηδενικά στοιχεία του πίνακα είναι ίσα με:
Για παράδειγμα, σε έναν πενταδιάστατο χώρο N = 5, ο πίνακας Hessian για τη συνάρτηση Rosenbrock έχει τη μορφή ζώνης:
Κώδικας που υπολογίζει αυτό το Hessian μαζί με κώδικα για την ελαχιστοποίηση της συνάρτησης Rosenbrock χρησιμοποιώντας τη μέθοδο συζυγούς κλίσης (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]Ένα παράδειγμα με τον ορισμό της συνάρτησης γινομένου της Έσσης και ένα αυθαίρετο διάνυσμα
Σε προβλήματα του πραγματικού κόσμου, ο υπολογισμός και η αποθήκευση ολόκληρου του πίνακα Hessian μπορεί να απαιτήσει σημαντικό χρόνο και πόρους μνήμης. Σε αυτήν την περίπτωση, στην πραγματικότητα δεν χρειάζεται να προσδιορίσετε τον ίδιο τον Hessian matrix, γιατί η διαδικασία ελαχιστοποίησης απαιτεί μόνο ένα διάνυσμα ίσο με το γινόμενο του Hessian με ένα άλλο αυθαίρετο διάνυσμα. Έτσι, από υπολογιστική σκοπιά, είναι πολύ προτιμότερο να ορίσουμε αμέσως μια συνάρτηση που επιστρέφει το αποτέλεσμα του γινόμενου της Έσσιας με ένα αυθαίρετο διάνυσμα.
Θεωρήστε τη συνάρτηση hess, η οποία παίρνει το διάνυσμα ελαχιστοποίησης ως πρώτο όρισμα και ένα αυθαίρετο διάνυσμα ως δεύτερο όρισμα (μαζί με άλλα ορίσματα της συνάρτησης που πρόκειται να ελαχιστοποιηθεί). Στην περίπτωσή μας, ο υπολογισμός του γινόμενου του Hessian της συνάρτησης Rosenbrock με ένα αυθαίρετο διάνυσμα δεν είναι πολύ δύσκολος. Αν p είναι ένα αυθαίρετο διάνυσμα, τότε το γινόμενο έχει τη μορφή:
Η συνάρτηση που υπολογίζει το γινόμενο του Hessian και ενός αυθαίρετου διανύσματος μεταβιβάζεται ως τιμή του ορίσματος 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Συζυγής αλγόριθμος αξιοπιστίας περιοχής κλίσης (Newton)
Η κακή ρύθμιση του πίνακα Hessian και οι λανθασμένες οδηγίες αναζήτησης μπορούν να προκαλέσουν αναποτελεσματικότητα του συζυγούς αλγόριθμου κλίσης του Newton. Σε τέτοιες περιπτώσεις, προτιμάται (εμπιστοσύνη-περιοχή) συζευγμένες διαβαθμίσεις Newton.
Παράδειγμα με τον ορισμό του Hessian matrix:
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.]Παράδειγμα με τη συνάρτηση γινομένου του Hessian και ένα αυθαίρετο διάνυσμα:
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. Επιπλέον, η λύση στο δευτεροβάθμιο υποπρόβλημα βρίσκεται με μεγαλύτερη ακρίβεια από τη χρήση της μεθόδου εμπιστοσύνης-ncg.
Παράδειγμα με τον ορισμό του Hessian matrix:
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.]
Παράδειγμα με τη συνάρτηση γινομένου του Hessian και ένα αυθαίρετο διάνυσμα:
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) είναι κατάλληλες για την επίλυση προβλημάτων μεγάλης κλίμακας (με χιλιάδες μεταβλητές). Αυτό οφείλεται στο γεγονός ότι ο υποκείμενος αλγόριθμος βαθμίδας συζυγούς υπονοεί έναν κατά προσέγγιση προσδιορισμό του αντίστροφου Hessian matrix. Η λύση βρίσκεται επαναληπτικά, χωρίς ρητή επέκταση της Έσσιας. Δεδομένου ότι χρειάζεται μόνο να ορίσετε μια συνάρτηση για το γινόμενο ενός Hessian και ενός αυθαίρετου διανύσματος, αυτός ο αλγόριθμος είναι ιδιαίτερα καλός για εργασία με αραιούς (διαγώνιο ζώνης) πίνακες. Αυτό παρέχει χαμηλό κόστος μνήμης και σημαντική εξοικονόμηση χρόνου.
Για προβλήματα μεσαίου μεγέθους, το κόστος αποθήκευσης και παραγοντοποίησης του Hessian δεν είναι κρίσιμο. Αυτό σημαίνει ότι είναι δυνατό να επιτευχθεί μια λύση σε λιγότερες επαναλήψεις, επιλύοντας σχεδόν ακριβώς τα υποπροβλήματα της περιοχής εμπιστοσύνης. Για να γίνει αυτό, μερικές μη γραμμικές εξισώσεις λύνονται επαναληπτικά για κάθε δευτεροβάθμιο υποπρόβλημα. Μια τέτοια λύση απαιτεί συνήθως 3 ή 4 αποσυνθέσεις Cholesky της μήτρας της Έσσης. Ως αποτέλεσμα, η μέθοδος συγκλίνει σε λιγότερες επαναλήψεις και απαιτεί λιγότερους υπολογισμούς αντικειμενικής συνάρτησης από άλλες εφαρμοσμένες μεθόδους περιοχής εμπιστοσύνης. Αυτός ο αλγόριθμος περιλαμβάνει μόνο τον προσδιορισμό του πλήρους πίνακα Hessian και δεν υποστηρίζει τη δυνατότητα χρήσης της συνάρτησης προϊόντος του Hessian και ενός αυθαίρετου διανύσματος.
Παράδειγμα με ελαχιστοποίηση της συνάρτησης 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.])Μάλλον θα σταματήσουμε εκεί. Στο επόμενο άρθρο θα προσπαθήσω να πω τα πιο ενδιαφέροντα πράγματα για την ελαχιστοποίηση υπό όρους, την εφαρμογή της ελαχιστοποίησης στην επίλυση προβλημάτων προσέγγισης, την ελαχιστοποίηση μιας συνάρτησης μιας μεταβλητής, τους αυθαίρετους ελαχιστοποιητές και την εύρεση των ριζών ενός συστήματος εξισώσεων χρησιμοποιώντας το scipy.optimize πακέτο.
Πηγή:
Πηγή: www.habr.com