אופטימיזציית הגחלילית

אופטימיזציית הגחלילית (באנגלית: Firefly optimization) היא שיטה מטה-היוריסטית מתחום הבינה המלאכותית לאופטימיזציה של פונקציות המקבלת את השראתה מהתנהגותן של גחליליות בטבע. השיטה הוצעה על ידי בשנת 2008, והיא שימושית בתחומים רבים, לרבות בעיבוד תמונה, תזמון תהליכים, ניתוח מבני של חלבונים, קינימטיקה הפוכה ועוד. סיבוכיות האלגוריתם היא <math>O(t)</math>, כאשר <math>t</math> הוא מספר האיטרציות שמבצע האלגוריתם.

==הקדמה==
===בעיית אופטימיזציה===
גחלילית מזן [[Photuris lucicrescens
]]

בעיית אופטימיזציה היא בעיה מתמטית העוסקת במציאת ערך אופטימלי (מינימלי או מקסימלי) עבור פונקציה נתונה, תחת אילוצים נתונים. כל בעיית אופטימיזציה ניתנת לתיאור ע"י:

<math> \begin{align}
& \text{optimize} && f(x) \\
& \text{subject to} && \mathbf{x\in S} \\
\end{align} </math>

כאשר <math>\text{S}</math> היא קבוצת הפתרונות האפשריים לבעיה. <math>s</math> הוא פתרון לבעיה אם ורק אם מתקיים <math>s\in S</math> וגם <math>f(x)\leq f(s)</math> לכל <math>x\in S</math>.

=== הגחליליות בטבע ===

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

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

==האלגוריתם ==
=== תיאור ===
בהינתן מספר איטרציות מרבי <math>t_{max}</math>, ניתן לתאר את האלגוריתם בצורה הכללית הבאה:
# אתחל אוכלוסייה אקראית <math>x_1,...,x_n</math> של <math>n</math> גחליליות על פני פונקציית מטרה <math>f(x)</math> ב- <math>d</math> ממדים, כאשר <math>f(x_i)</math>היא עוצמת ההארה של גחלילית <math>x_i</math>.
# לכל זוג גחליליות <math>x_i, x_j</math>, אם גחלילית <math>x_i</math> בהירה יותר מגחלילית <math>x_j</math>, הזז את גחלילית <math>x_j</math> לכיוון גחלילית <math>x_i</math> על פי משוואת התנועה: <math>\mathbf{x}_i^{t+1}=\mathbf{x}_i^t + \beta \exp[-\gamma r_{ij}^2] (\mathbf{x}_j^t - \mathbf{x}_i^t) +\alpha_t \boldsymbol{\epsilon}_t </math> כאשר <math>t</math> מייצג נקודה בזמן, <math>\ x_i, x_j</math> הן מיקום הגחליליות <math>i,j</math> בהתאמה, <math>\alpha_t </math> הוא מקדם אקראיות ו- <math>\boldsymbol{\epsilon}_t </math> הוא וקטור מספרים אקראיים המתפלגים נורמלית בזמן <math>t</math>. אחרת, הזז את הגחלילית <math>x_j</math> באופן אקראי.
# עדכן את תאורת הגחליליות ואת מידת המשיכה בין הגחליליות ע"פ המיקומים והמרחקים החדשים שהתקבלו.
# חזור על שלבים (2) ו- (3) למשך <math>t_{max}</math> פעמים. בסיום, הגחלילית המוארת ביותר היא קירוב לפתרון הבעיה.
בבעיות למציאת מקסימום, יתקיים <math>I \propto f(\mathbf{x})</math> ובבעיות למציאת מינימום יתקיים <math>I \propto f^{-1}(x)</math>.

=== פסאודו - קוד ===
<div class="mw-content-ltr">

 Begin
 Generate a random population of n fireflies x_1, x_2, ..., x_n
 While (t < MaxGeneration)
 update distances
 update lights intensities
 update attractiveness
 for i = 1 : n
 for j = 1 : n
 if (I_j > I_i),
 move firefly i to j
 end if
 end for j
 end for i
 set s as the firefly with the optimal light intensity.
 end while
 return s as solution
 end

</div>

== רקע תאורטי ==
=== מודל מתמטי ===
האלגוריתם מבוסס על אידיאליזצית תכונות ההארה של מאור הגחליליות. הגחלילית, בקרב מרבית זני הגחליליות בטבע, משתמשת בהארה על מנת למשוך אליה פרטי להקה אחרים לצורכי זיווג. התכונות האידיאליות הן כדלקמן:
# כל גחלילית נמשכת לגחליליות אחרות, ללא קשר למינן.
# מידת המשיכה של גחלילית <math>a</math> כלפי גחלילית <math>b</math> היא ביחס ישר למידת ההארה של גחלילית <math>b</math> (אנלוגיה להתנהגותן הטבעית של הגחליליות).
# עוצמת האור שגחלילית <math>a</math> קולטת מגחלילית <math>b</math> היא ביחס הפוך למרחק בין גחלילית <math>a</math> לבין גחלילית <math>b</math> (אנלוגיה לחוק ריבוע המרחק ההופכי, לפיו עוצמת האור נחלשת ביחס ריבועי הפוך למרחק ממקור האור).
# עוצמת ההארה של הגחלילית נקבעת ע"פ ערכה של פונקציית המטרה <math>f(x)</math>.
על פי תכונות (2) ו- (3), ניתן לגזור את משוואת המשיכה בין שתי גחליליות, כפונקציה מעריכית של המרחק ביניהן:
<math>\beta = \beta_0e^{-\gamma r^2}</math>
כאשר <math>\beta_0</math> הוא קבוע, <math>\gamma</math> הוא מקדם התפשטות האור באוויר, <math>\beta</math> היא עוצמת המשיכה בין שתי גחליליות ו- <math>r</math> הוא המרחק הגאומטרי בין שתי הגחליליות. תנועתה של גחלילית <math>i</math> מושפעת מגחלילית בהירה יותר <math>j</math> באופן הבא:

<math>\mathbf{x}_i^{t+1}=\mathbf{x}_i^t + \beta \exp[-\gamma r_{ij}^2] (\mathbf{x}_j^t - \mathbf{x}_i^t) +\alpha_t \boldsymbol{\epsilon}_t </math>

כאשר <math>t</math> מייצג נקודה בזמן, <math>\ x_i, x_j</math> הן מיקום הגחליליות <math>i,j</math> בהתאמה, <math>\alpha_t </math> הוא מקדם אקראיות ו- <math>\boldsymbol{\epsilon}_t </math> הוא וקטור מספרים אקראיים המתפלגים נורמלית בזמן <math>t</math>. עבור <math>\beta_0=0</math> מתקבלת תנועה אקראית פשוטה, ועבור <math>\gamma=0</math> מתקבלת משוואת תנועת חלקיקים של אופטימיזציית הנחיל.

=== חלוקה אוטומטית ===
במקרים רבים, אלגוריתם הגחליליות יעיל יותר מאלגוריתמים רבים אחרים מתחום בינת הנחיל. סיבה עיקרית לכך היא תופעת החלוקה האוטומטית. האלגוריתם מבוסס על כך שאטרקטיביות של פתרון פוחתת עם עליית המרחק, ולכן, נוצר מצב שבו האוכלוסייה מופרדת באופן אוטומטי לתתי-קבוצות, כאשר כל תת-קבוצה מרוכזת סביב אופטימום מקומי שונה. עבור אוכלוסייה הגדולה מאוד ממספר מוקדי האופטימום המקומיים, חלוקה זו מאפשרת למצוא את כל המוקדים המקומיים בו-זמנית. מבחינה מתמטית, <math>{\sqrt{\gamma}}^{-1}
</math> מאפיין את המרחק הממוצע בין שתי קבוצות גחליליות סמוכות. מכאן, כל אוכלוסייה נחלקת לתתי קבוצות סמוכות בעלות מרחק ממוצע קבוע. עובדה זו מאפשרת התמודדות עם בעיות שאינן ליניאריות.

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

<math>f(x,y)=e^{-x^{2}-y^{2}}({\left | x \right |}+{\left | y \right |})</math>

לפונקציה <math>f(x,y)</math> ישנן ארבע נקודות מקסימום <math>(\pm 0.5, \pm 0.5)</math>. האלגוריתם מאפשר את מציאת הקטבים על ידי <math>25</math> גחליליות ו- <math>20</math> איטרציות. הסכימה הבאה מתארת את מיקומן האקראי של הגחליליות (כנקודות) בתחילתה של ריצת האלגוריתם (משמאל) ואת מיקומן בסיום ריצת האלגוריתם:
מרכז

===התמקדות כנגד התרחבות===
==== רעיון כללי ====
שני המרכיבים העיקריים בכל שיטה מטה-היוריסטית הם התמקדות והתרחבות.
* התמקדות - חיפוש מקומי ממוקד בתחום שבו קיים פתרון אופטימלי מקומי. מידת התמקדות גבוהה נוטה להגדיל את קצב ההתכנסות של האלגוריתם, אך עשויה לגרום להתכנסות מקומית בלתי רצויה. 
* התרחבות - חיפוש גלובלי הבוחן את כל הפתרונות האפשריים, תוך הימנעות מהתכנסות לפתרון אופטימלי מקומי. מידת התרחבות גבוהה נוטה להקטין את קצב ההתכנסות של האלגוריתם, אך מעלה משמעותית את ההסתברות למציאת אופטימום גלובלי. 
אלגוריתם המחלק את פעולתו להתמקדות והתרחבות נקרא אלגוריתם חיפוש-לסירוגין. רמות ההתמקדות וההתרחבות של האלגוריתם, וכמו כן מידת האיזון בין שני המרכיבים, נקבעות על ידי מידת האקראיות וערכי הפרמטרים. יעילותו ודיוקו של האלגוריתם תלוי ביחס ההתמקדות-התרחבות <math>Q</math>:

<math>Q=\frac{\tau_a}{\tau_b}</math>

כאשר <math>\tau_a</math> ו- <math>\tau_b</math> הם פרקי הזמן המנוצלים לתהליך ההתמקדות ולתהליך ההתרחבות, בהתאמה.

==== מציאת יחס התמקדות-התרחבות אופטימלי ====
חשיבותו של יחס התמקדות-התרחבות מוביל לפיתוח המודל התאורטי הבא, עבור המקרים הדו-ממדיים והתלת-ממדיים: אלגוריתם חיפוש המתמקד במרחב פתרונות קטן בעל רדיוס <math>a</math> הנמצא במרחב פתרונות שרדיוסו <math>b</math>, כאשר <math>a\ll b</math>, ניתן למידול על ידי תהליך דיפוזיבי איטי המתואר על ידי מערכת המשוואות הסטוכסטיות הדיפרנציאליות הבאות:

<math>\left\{\begin{matrix}
D\nabla_r^2t_1+\frac{1}{2 \pi \tau_a } \int_{0}^{2 \pi}\left [ t_2(r)-t_1(r) \right ]d \theta +1 =0 \\
 \\
u\nabla_rt_2(r)-\frac{1}{\tau_b }\left [ r_2(r)-t_1(r) \right ] + 1 =0
\end{matrix}\right.</math>

כאשר <math>D</math> הוא קבוע הדיפוזיה, <math>\tau_a</math> ו- <math>\tau_b</math> הם תוחלות הזמן שנוצלו לתהליך ההתמקדות ולתהליך ההתרחבות, בהתאמה, <math>t_1</math> ו- <math>t_2</math> הם תוחלות זמני ההגעה למטרה בתהליך ההתמקדות ובתהליך ההתרחבות, בהתאמה, ו- <math>u</math> הוא קצב החיפוש הממוצע. תחת ההנחה שקצב החיפוש <math>u
</math> זהה בכל הצעדים, פרק הזמן המינימלי הדרוש לכל אחד מהשלבים הוא:

<math>\tau_{a}^{\mathrm{min}}\approx \frac{D}{2u^2}\cdot \frac{\ln^2\frac{b}{a}}{2\ln\frac{b}{a}-1} \quad \quad

\tau_{b}^{\mathrm{min}}\approx \frac{a}{u}\sqrt{\ln\frac{b}{a}-\frac{1}{2}}</math>

כאשר <math>u\rightarrow \infty </math>, מתקבל יחס ההתמקדות-התרחבות האופטימלי:

<math>(1) \quad Q_{\mathrm{optimal}}=\frac{\tau_a}{\tau_b^2} \approx \frac{D}{a^2}\frac{1}{\left [ 2- \frac{1}{\ln{\frac{b}{a}}} \right ]^2} </math>

התוצאות לעיל קבילות עבור המקרה הדו-מימדי, ושלא קיימות תוצאות כלליות עבור יותר מ-3 ממדים.

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

<math>D\approx \frac {s^2} {2}</math>

כאשר <math>D</math>הוא קבוע הדיפוזיה ו-<math>s</math> הוא אורך הצעד של המהלך בכל איטרציה. במקרה שבו <math>b / a\rightarrow \infty </math> נקבל מתוצאה (1):

<math>Q_{\mathrm{optimal}}=\frac{\tau_a}{\tau_b^2}\approx \frac{1}{8}</math>

מכאן ניתן להסיק שבמקרה זה, יש להשקיע זמן רב יותר בשלב ההתרחבות מאשר בשלב ההתמקדות.

==== דוגמה - מציאת יחס התמקדות-התרחבות אופטימלי של פונקציית גל עומד ====
גרף של גל הרמוני העומד על דיסק
נדגים את חישוב יחס ההתמקדות-התרחבות במקרה של פונקציית גל-עומד:

<math>f(x)=1+\left\{
\exp \left [
-\sum_{i=1}^{d}
\left ( \frac{x_i}{\beta } \right )^{10}
 \right ] - 2 \exp
\left [
-\sum_{i=1}^{d} \left ( x_i-\pi \right )^2
 \right ]
 \right\}
\cdot \prod_{i=1}^{d}\cos ^2x_i</math>

פונקציה זו מאופיינת במספר גדול של נקודות קיצון מקומיות, ובעלת מינימום גלובלית יחידה <math>f_{\mathrm {min}}=0</math> בנקודה <math>(\pi,\pi,\dots,\pi)</math> בתחום <math>-20\leq x_i \leq 20</math> כאשר <math>i=1,2, \dots ,d</math> ו- <math>\beta = 15</math>. נקבע ש- <math>\alpha = \pi/2</math>, ומכיוון ש- <math>b=20</math>, נובע ש- <math>b/a=12.7</math>, ועבור <math>d=2</math>נקבל:

<math>Q_{\mathrm {optimal}} \approx
\frac{1}{2\left [ 2-\ln\left ( \frac{b}{a} \right ) \right ]^2} \approx 0.19</math>

מכאן ניתן להסיק, שבהפעלת האלגוריתם על פונקציית גל-עומד, על האלגוריתם להשקיע 80% מזמן החיפוש בהתרחבות ו- 20% מזמן החיפוש בהתמקדות, על מנת לבצע את החישוב בזמן אופטימלי.

הטבלה הבאה מתארת רצף של חמישה ניסויים, בהם מופעל אלגוריתם הגחליליות למציאת מינימום עבור פונקציית גל-עומד, כאשר <math>n=15</math>, <math>t=1000</math>, <math>x \in [-20,20]^n</math>, <math>\beta = 15</math> ו- <math>Q</math> משתנה.
השוואה של ביצועי אלגוריתם הגחליליות מול האלגוריתם של בניצ'ו אטאל, עבור מספר ממדים משתנה.

{| class="wikitable" width="400px"
|<math>Q</math>
|<math>0.4</math>
|<math>0.3</math>
|<math>0.2</math>
|<math>0.1</math>
|<math>0.05</math>
|-
|תוצאה
|<math>9.4e-11</math>
|<math>9.4e-11</math>
|<math>9.4e-11</math>
|<math>9.4e-11</math>
|<math>9.4e-11</math>
|-
|}

קל להבחין שהתוצאה המדויקת ביותר מתקבלת כאשר <math>Q=0.2</math> - הבחנה אמפירית המתלכדת עם הטענה התאורטית.

==== הכללה לממדים גבוהים ====
אלגוריתם החיפוש-לסירוגין אשר הוצע על ידי בניצ'ו אטאל עבור ממדים גבוהים יותר, מתבסס על היחסים האסימפטוטיים הבאים, המתקבלים באמצעות אקסטרפולציה:

<math>\frac{\tau_a}{\tau_b^2} \sim O\left ( \frac{D}{a^2} \right ) \quad \quad \tau_m \sim O\left ( \frac{b}{u}\left ( \frac{b}{a} \right )^{d-1} \right )</math>

כאשר <math>\tau_m</math> הוא זמן החיפוש הכולל הממוצע. הודגם אמפירית כי בממדים גבוהים, אלגוריתם הגחליליות יעיל בהרבה מן האלגוריתם הנ"ל, אך לא נמצא מודל תאורטי המסביר את תוצאות הניסוי.

== קביעת הפרמטרים ==
מקדם האקראיות <math>\alpha_t </math> מווסת את מידת האקראיות של מהלכי הגחליליות. כיוון שהפחתה הדרגתית של מידת האקראיות תאפשר מניעה של התכנסות לאופטימום מקומי, התכנסות מהירה תתקבל עבור <math>\alpha_t=a_0 \delta^t</math>

כאשר <math> \delta\in[0,1]</math> ו- <math>a_0</math> הוא מקדם אקראיות התחלתי. סימולציות מראות שהתכנסות יעילה במיוחד מתקבלת עבור <math>\gamma=1/\sqrt{L}</math>, <math>n\in[25,40] </math>, <math> a_0=0.01L </math>, כאשר <math>L</math> הוא סדר הגודל הממוצע של פתרונות הבעיה.

==סיבוכיות זמן ==
סיבוכיות האלגוריתם היא <math>O(n^{2}t)</math>, כאשר <math>t</math> הוא מספר האיטרציות ו- <math>n</math> הוא מספר הגחליליות. עבור <math>n</math> גדול יחסית, ניתן להחליף את הלולאה הכפולה בלולאה יחידה, ולמיין את מידת המשיכה של הגחליליות באמצעות אלגוריתם מיון השוואתי. במקרה זה, סיבוכיות האלגוריתם היא <math>O(nt\log{n})</math>. בשני המקרים, הסיבוכיות ליניארית ב- <math>t</math>.

== הערכת ביצועים ==
להלן מספר האיטרציות שנדרשו לפתרון בעיות אופטימיזציה עבור שלוש פונקציות מרכזיות, בטווח נתון, וכתלות במספר הממדים של הפונקציה (<math>n</math>).

=== פונקציית ספירה (Sphere function) ===
גרף של פונקציית ספירה עבור n=2.
<math>\textrm{Sphere}(\boldsymbol{x}) = \sum_{i=1}^{n} x_{i}^{2} \quad\quad [-100,100]^n</math>
{|class="wikitable" width="400px"
! <math>n</math> !! מינימום !! מקסימום
|-
|15
|6239
|6440
|-
|30
|6202
|6452
|-
|60
|6149
|6398
|}

=== פונקציית אקליי (Ackley function) ===
גרף של פונקציית אקליי עבור n=2.
<math>\textrm{Ackley}(\boldsymbol{x}) =-20 \exp \left [ -\frac{1}{5}\sqrt{\frac{1}{n}\sum_{i=1}^{n}{x_{i}^{2}}} \right ]-
\exp \left [ \frac{1}{n} \sum_{i=1}^{n}{\cos (2\pi x_i)} \right]+20+e\quad\quad [-30,30]^n</math>
{| class="wikitable" width="400px"
!
<math>n</math>
!
מינימום
!
מקסימום

|-
|15
|1720
|2176

|-
|30
|1568
|1805

|-
|60
|1402
|1862
|}
=== פונקציית רוזנברוק ===
גרף של פונקציית רוזנברוק עבור n=2.
<math>\textrm{Rosenbrock}(\boldsymbol{x}) = \sum_{i=1}^{n-1} \left[ 100 \left(x_{i+1} - x_{i}^{2}\right)^{2} + \left(x_{i} - 1\right)^{2}\right]\quad\quad [-30,30]^n</math>
{| class="wikitable" width="400px"
!
<math>n</math>
!
מינימום
!
מקסימום
|-
|15
|2497
|2938
|-
|30
|2531
|2653
|-
|60
|2471
|2574
|}

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

== ראו גם==
* אופטימיזציית הנחיל
* בינה מלאכותית
* אלגוריתם קן הנמלים
* תנועה בראונית

== לקריאה נוספת ==

== קישורים חיצוניים ==
* 
* [http://code.google.com/p/csc6810project/ מימוש לאלגוריתם הגחליליות ב- Python]
* [https://github.com/firefly-cpp/Firefly-algorithm--FFA- מימוש לאלגוריתם הגחליליות ב- C]
* [http://www.mathworks.com/matlabcentral/fileexchange/29693-firefly-algorithm מימוש לאלגוריתם הגחליליות ב- Mathlab]

קטגוריה:אופטימיזציה רציפה
קטגוריה:אלגוריתמים