על חישוב עקבה במטריצה

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

הבה נבין את המוטיבציה האלגברית מאחורי כל אחת מהשיטות, שלב אחר שלב.

מהי עקבה ומה המשמעות שלה?

העקבה של מטריצה ריבועית A מוגדרת כסכום איברי האלכסון שלה:

Tr(A) = Σ Aᵢᵢ (i=1 to n)

המשפט היסודי של עקבות קובע שהעקבה שווה לסכום הערכים העצמיים של המטריצה:

Tr(A) = Σ λᵢ (i=1 to n)

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

שיטת האצ'ינסון

בשיטת האצ'ינסון, עליה כתבתי בבלוג פוסט מוקדם יותר, מקבלים הערכה סטטיסטית של העקבה של מטריצה. השיטה לא דורשת גישה לאיברים בודדים בתוך המטריצה A. היא דורשת רק לדעת לבצע מכפלת מטריצה בוקטור. זה הופך אותה לאידיאלית עבור מטריצות שמוגדרות על ידי פונקציה או אופרטור. הבעיה היא שלשיטה זו יש שונות גבוהה. אם למטריצה יש ערכים עצמיים גדולים מאוד, הממוצע יקפוץ מצד לצד ותדרשנה דגימות רבות כדי להתכנס לתוצאה מדויקת, התכנסות של O(1/ε²).

פירוק לערכים סינגולריים, SVD, לצורך עקבה

אם האצ'ינסון סובל משונות בגלל ערכים עצמיים גדולים, למה שלא נחשב אותם ישירות?

המוטיבציה האלגברית היא שאם A היא מטריצה סימטרית ומוגדרת חיובית, פירוק SVD שלה חושף את הערכים העצמיים λᵢ. אנו יכולים לחשב את העקבה המדויקת על ידי סיכום שלהם:

Tr(A) = Σ λᵢ (i=1 to n)

נו אז מה הבעיה? ביצוע SVD מלא למטריצה מסדר n דורש סיבוכיות זמן ריצה של (O(n³. זה בלתי אפשרי למטריצות ענק. אולי נוכל למצוא רק את הערכים העצמיים הגדולים ביותר ולהשתמש בהם כדי לקרב את העקבה? כאן נכנסות לתמונה שיטות מתקדמות יותר.

שיטות מרחב קרילוב

שיטות קרילוב מנסות להבין את מבנה המטריצה בעזרת מספר קטן של הכפלות, ללא SVD מלא. המוטיבציה האלגברי היא שבמקום להכפיל וקטור אקראי z במטריצה A ולזרוק את התוצאה, אנו שומרים את ההיסטוריה ויוצרים מרחב נפרש:

Kₖ(A, v) = span{v, Av, A²v, ..., Aᵏ⁻¹v}

בעזרת אלגוריתם כמו לנצוש, Lanczos, אנו מטילים את המטריצה הגדולה A על מרחב קרילוב הקטן, בגודל k × k. המטריצה הקטנה שנוצרת משמרת בצורה מצוינת את הערכים העצמיים הקיצוניים הגדולים והקטנים של A. שילוב של מרחב קרילוב עם וקטורים אקראיים מאפשר להעריך עקבות של פונקציות של מטריצות, כמו Tr(eᴬ), ביעילות עצומה. שני המושגים האלו, לנצוש ו-SLQ, גרסה סטוכסטית שלו, הם בליבה של חישובים נומריים מתקדמים למטריצות ענק, והם משלבים אלגברה לינארית תיאורטית עם סטטיסטיקה בצורה יפהפייה. נצלול למוטיבציה האלגברית ולמנגנון של כל אחד מהם.

אלגוריתם לנצוש

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

כאן נכנס מרחב קרילוב. אם ניקח וקטור אקראי v, נכפיל אותו ב-A, ואז נכפיל את התוצאה שוב ב-A, נייצר סדרה של וקטורים: v, Av, A²v, A³v, ...

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

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

ניסוח אלגוריתם לנצוש

מתחילים עם וקטור אקראי מנורמל v₁ וקובעים v₀ = 0 ו- β₀ = 0.

בכל איטרציה j (מ-1 עד k):

1. מחשבים w = Avⱼ.

2. מוצאים את הרכיב על האלכסון: αⱼ = vⱼᵀ * w.

3. מחסרים את הרכיבים הידועים כדי לקבל וקטור מאונך: w = w - αⱼvⱼ - βⱼ₋₁vⱼ₋₁.

4. מחשבים את הנורמה החדשה: βⱼ = ||w||.

5. הווקטור הבא בבסיס הוא vⱼ₊₁ = w / βⱼ.

ונקבל מטריצה קטנה Tₖ, בגודל k × k, כאשר k << n, שהיא תלת-אלכסונית: האיברים αⱼ יושבים על האלכסון הראשי, והאיברים βⱼ יושבים מעליו ומתחתיו. הערכים העצמיים של Tₖ הקטנה, Ritz values, הם קירוב מצוין לערכים העצמיים הקיצוניים, הגדולים והקטנים ביותר, של A הענקית.

Stochastic Lanczos Quadrature (SLQ)

מה זה SLQ? מדובר באלגוריתם שמשלב את שיטת האצ'ינסון עם אלגוריתם לנצוש ושיטות אינטגרציה נומריות, כדי לחשב ביעילות עקבה של פונקציה של מטריצה, כלומר Tr(f(A)). למשל, פונקציית אקספוננט Tr(eᴬ) או לוגריתם Tr(log A) שמשמשות רבות בלמידת מכונה ופיזיקה סטטיסטית.

נזכור ששיטת האצ'ינסון מאפשרת לנו להעריך עקבה בעזרת וקטורים אקראיים v:

Tr(f(A)) ≈ E[vᵀ f(A) v]

הבעיה היא שלחשב את המכפלה vᵀ f(A) v ישירות זה קשה מאוד, כי אנחנו לא יכולים לחשב את f(A) המלאה.

כאן באה פריצת דרך אלגברית ואנליטית: ניתן לראות את הביטוי

vᵀ f(A) v

כאינטגרל של הפונקציה f(x) לפי פיזור הערכים העצמיים של A.

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

המטריצה התלת-אלכסונית Tₖ שמפיק אלגוריתם לנצוש מציעה את מה שאנחנו צריכים:

• הערכים העצמיים של Tₖ הם בדיוק ה-נקודות הדגימה, ה-Nodes, של הקוודרטורה.

• הרכיבים הראשונים של הווקטורים העצמיים של Tₖ מספקים בדיוק את המשקלים.

איך זה עובד בפועל? האלגוריתם מבצע את השלבים הבאים:

1. מגרילים וקטור אקראי מנורמל v, למשל מתפלג רדמאכר או גאוסי.

2. מריצים את אלגוריתם לנצוש עם המטריצה A ווקטור ההתחלה v, למשך k צעדים. התוצאה: המטריצה התלת-אלכסונית הקטנה Tₖ.

3. מחשבים את הפונקציה על המטריצה הקטנה בלבד: בגלל ש-Tₖ קטנה מאוד, קל לחשב עליה ישירות את f(Tₖ) בעזרת פירוק מלא לערכים עצמיים.

4. הקירוב שלנו לאיבר vᵀ f(A) v מתקבל פשוט על ידי הסתכלות על האיבר הראשון, השמאלי העליון, של f(Tₖ):

   
   

   vᵀ f(A) v ≈ ||v||² * [f(Tₖ)]₁,₁

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

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

שיטת ++Hutch

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

• החלק הכבד, מכיל את הערכים העצמיים הגדולים ומיוצג על ידי מטריצה מקרבת מדרגה נמוכה Aₖ.

• החלק הקל / השארית (A - Aₖ).

אלגוריתם עובד כך:

1. משתמשים באלגוריתם אקראי Randomized SVD כדי למצוא קירוב של התת-מרחב של הערכים העצמיים הגדולים.

2. מחשבים את העקבה של החלק הזה באופן מדויק.

3. מפעילים את שיטת האצ'ינסון הסטנדרטית רק על השארית.

מכיוון שהשארית מכילה כעת רק ערכים עצמיים קטנים, הספקטרום שוטח, השונות של שיטת האצ'ינסון צונחת לאפס כמעט. והתוצאה? קצב ההתכנסות משתפר דרמטית מ- (O(1/ε² ל- (O(1/ε.

איך hutch++ משתמש במרחבי קרילוב אם בכלל?

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

1. בניית תת-מרחב הדרגה הנמוכה

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

• הוא לוקח מטריצה אקראית S ומחשב את המכפלה: A כפול S.

• לפעמים הוא עושה זאת שוב: A כפול (A כפול S).

• הקבוצה הזו של וקטורים {S, AS, A²S} היא למעשה בסיס למרחב קרילוב.

2. קירוב המטריצה

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

3. פיצול העבודה

הנה התכסיס שמבדיל את Hutch++ משיטת Hutchinson הרגילה:

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

• הוא משתמש בשיטה אקראית, בהאצ'ינסון, רק כדי להעריך את השארית, כלומר, מה ש-Q לא הצליחה לתפוס מהמטריצה המקורית.

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

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

על חישוב עקבה במטריצה

אנחנו במתמטיקאי מחקר ופיתוח בע"מ מספקים שירותים של מחקר אלגוריתמי יישומי, שירותים של Fractional CTO / Chief Scientist ומייעצים בנושאים של מתמטיקה שימושית, תכן ניתוח שיפור ופיתוח של אלגוריתמים ומבני נתונים. דברו איתי:

שלמה יונה

מייסד ומדען ראשי, מתמטיקאי מחקר ופיתוח בע"מ

053-7326360

shlomo.yona@mathematic.ai 

https://mathematic.ai

פודקאסט על החברה ועליי, שלמה יונה, ואופן העבודה שלנו ואיתנו: A technical deep dive about Mathematic.ai