🥇SciPy, אופטימיזציה

SciPy (מבוטא sai pie) היא חבילת אפליקציות מתמטית המבוססת על הרחבה Numpy Python. עם SciPy, הפעלת Python האינטראקטיבית שלך הופכת לאותה סביבת מדעי נתונים מלאה וסביבת אב-טיפוס מורכבת כמו MATLAB, IDL, Octave, R-Lab ו-SciLab. היום אני רוצה לדבר בקצרה על איך להשתמש בכמה אלגוריתמי אופטימיזציה ידועים בחבילת scipy.optimize. עזרה מפורטת ועדכנית יותר על שימוש בפונקציות תמיד אפשר לקבל באמצעות הפקודה help() או באמצעות Shift+Tab.

מבוא

על מנת להציל את עצמך ואת הקוראים מחיפוש וקריאת מקורות ראשוניים, קישורים לתיאורי שיטות יהיו בעיקר בוויקיפדיה. ככלל, מידע זה מספיק כדי להבין את השיטות במונחים כלליים ואת התנאים ליישומם. כדי להבין את מהות השיטות המתמטיות, עקוב אחר הקישורים לפרסומים מוסמכים יותר, אותם ניתן למצוא בסוף כל מאמר או במנוע החיפוש המועדף עליך.

אז, מודול scipy.optimize כולל יישום של ההליכים הבאים:

מזעור מותנה ובלתי מותנה של פונקציות סקלריות של מספר משתנים (מינימום) באמצעות אלגוריתמים שונים (Nelder-Mead simplex, BFGS, ניוטון מצומד גרדיאנטים, קובילה и SLSQP)
אופטימיזציה גלובלית (לדוגמה: אגן שופינג, diff_evolution)
מזעור שאריות MNC (לפחות_ריבועים) ואלגוריתמים של התאמת עקומה באמצעות ריבועים קטנים לא ליניאריים (Curve_fit)
מזעור פונקציות סקלריות של משתנה אחד (minim_scalar) וחיפוש שורשים (root_scalar)
פותרים רב מימדיים של מערכת משוואות (שורש) באמצעות אלגוריתמים שונים (פאוול היברידי, לבנברג-מרקורדט או שיטות בקנה מידה גדול כגון ניוטון קרילוב).

במאמר זה נשקול רק את הפריט הראשון מכל רשימה זו.

מזעור בלתי מותנה של פונקציה סקלרית של מספר משתנים

פונקציית המינימום מחבילת scipy.optimize מספקת ממשק כללי לפתרון בעיות מזעור מותנות ובלתי מותנות של פונקציות סקלריות של מספר משתנים. כדי להדגים איך זה עובד, נצטרך פונקציה מתאימה של כמה משתנים, אותם נמזער בדרכים שונות.

למטרות אלה, פונקציית Rosenbrock של N משתנים מושלמת, שיש לה את הצורה:

למרות העובדה שפונקציית Rosenbrock והמטריצות Jacobi וה-Hssian שלה (הנגזרת הראשונה והשנייה, בהתאמה) כבר מוגדרות בחבילת 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)

למען הבהירות, בואו נצייר בתלת מימד את הערכים של פונקציית רוזנברוק של שני משתנים.

קוד ציור

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.

שיטת Nelder-Mead Simplex

תהיה נקודה התחלתית x0 במרחב 5 מימדי. בואו נמצא את נקודת המינימום של פונקציית Rosenbrock הקרובה אליה באמצעות האלגוריתם Nelder-Mead simplex (האלגוריתם מצוין כערך של פרמטר השיטה):

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 היא בחירה טובה לבעיות מזעור פשוטות. עם זאת, מכיוון שהוא אינו משתמש בהערכות שיפוע, ייתכן שייקח זמן רב יותר למצוא את המינימום.

שיטת פאוול

אלגוריתם אופטימיזציה נוסף שבו רק ערכי הפונקציות מחושבים הוא שיטת פאוול. כדי להשתמש בו, עליך להגדיר שיטה = 'פאוול' בפונקציית המינימום.

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 משתמש בשיפוע של פונקציית המטרה. ניתן לציין את הגרדיאנט כפונקציה או לחשב באמצעות הפרשי מסדר ראשון. בכל מקרה, שיטת BFGS דורשת בדרך כלל פחות קריאות לפונקציה מאשר שיטת הסימפלקס.

הבה נמצא את הנגזרת של פונקציית רוזנברוק בצורה אנליטית:

ביטוי זה תקף עבור הנגזרות של כל המשתנים למעט הראשון והאחרון, המוגדרים כ:

בואו נסתכל על פונקציית 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]

אלגוריתם גרדיאנט מצומד (ניוטון)

אלגוריתם שיפועים מצומדים של ניוטון היא שיטת ניוטון שונה.
השיטה של ניוטון מבוססת על קירוב פונקציה באזור מקומי על ידי פולינום מהמעלה השנייה:

איפה היא המטריצה של נגזרות שניות (הסיאנית, הסיאנית).
אם ההסיאן מוגדר חיובי, אזי ניתן למצוא את המינימום המקומי של פונקציה זו על ידי השוואת שיפוע האפס של הצורה הריבועית לאפס. התוצאה תהיה הביטוי:

ההסיאן ההיפוך מחושב בשיטת שיפוע מצומד. דוגמה לשימוש בשיטה זו כדי למזער את פונקציית Rosenbrock ניתנת להלן. כדי להשתמש בשיטת ניוטון-CG, עליך לציין פונקציה שמחשבת את ההסיאן.
הפונקציה ההסינית של רוזנברוק בצורה אנליטית שווה ל:

איפה и , הגדר את המטריצה .

הרכיבים הנותרים שאינם אפס של המטריצה שווים ל:

לדוגמה, במרחב חמישה ממדי N = 5, המטריצה ההסיאנית עבור פונקציית רוזנברוק היא בצורת פס:

קוד שמחשב את ההסיאן הזה יחד עם קוד למזעור פונקציית רוזנברוק בשיטת שיפוע מצומד (ניוטון):

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, שלוקחת את וקטור המזעור כארגומנט הראשון, ווקטור שרירותי כארגומנט השני (יחד עם ארגומנטים אחרים של הפונקציה שיש למזער). במקרה שלנו, חישוב מכפלת ההסיאן של פונקציית רוזנברוק עם וקטור שרירותי אינו קשה במיוחד. אם p הוא וקטור שרירותי, ואז המוצר יש את הצורה:

הפונקציה שמחשבת את המכפלה של ה-Hssian ווקטור שרירותי מועברת כערך של הארגומנט 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

אלגוריתם אזור אמון גרדיאנט מצומד (ניוטון)

התניה לקויה של המטריצה ההסיאנית וכיווני חיפוש שגויים עלולים לגרום לאלגוריתם הגרדיאנט המצומד של ניוטון להיות לא יעיל. במקרים כאלה ניתנת עדיפות ל שיטת אמון באזור (אזור אמון) מצמידים גרדיאנטים של ניוטון.

דוגמה להגדרה של המטריצה ההסית:

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

שיטות מסוג קרילוב

בדומה לשיטת ה-trust-ncg, שיטות מסוג Krylov מתאימות היטב לפתרון בעיות בקנה מידה גדול מכיוון שהן משתמשות רק במוצרי וקטור מטריקס. המהות שלהם היא לפתור בעיה באזור ביטחון מוגבל על ידי תת מרחב קרילוב קטום. לבעיות לא ודאות, עדיף להשתמש בשיטה זו, מכיוון שהיא משתמשת במספר קטן יותר של איטרציות לא ליניאריות בגלל המספר הקטן יותר של מוצרי מטריקס-וקטור לכל תת-בעיה, בהשוואה לשיטת trust-ncg. בנוסף, הפתרון לבעיית המשנה הריבועית נמצא בצורה מדויקת יותר מאשר שימוש בשיטת trust-ncg.
דוגמה להגדרה של המטריצה ההסית:

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 פירוקים Cholesky של המטריצה ההסיאנית. כתוצאה מכך, השיטה מתכנסת בפחות איטרציות ודורשת פחות חישובי פונקציות אובייקטיביות מאשר שיטות אזורי ביטחון מיושמות אחרות. אלגוריתם זה מרמז רק על קביעת המטריצה ההסיאנית השלמה ואינו תומך ביכולת להשתמש בפונקציית המכפלה של ההסיאן ובוקטור שרירותי.

דוגמה עם מזעור של פונקציית רוזנברוק:

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/

מקור: www.habr.com

SciPy, אופטימיזציה