איך מחשבים דטרמיננטה עם אלגוריתם בריס?
חלק ראשון, שבו אנחנו לומדים שמינוס 32 הוא כמעט 256
בפוסט הקודם שלי על דטרמיננטות הסברתי איך מחשבים אותן. הצגתי שתי דרכים: אחת שהולכת דרך ההגדרה הרקורסיבית והיא מאוד לא יעילה לחישוב, ואחת שמבצעת פישוט למטריצה שרוצים לחשב את הדטרמיננטה שלה והיא משמעותית יותר יעילה. יש שיטות שהן אפילו יעילות יותר, אבל בפוסט הזה אני רוצה לקחת צעד אחד אחורה דווקא ולהציג שיטה יעילה פחות, אבל עם יתרון נחמד אחד: אם האיברים של המטריצה הם כולם מספרים שלמים, אנחנו לא נזדקק לשברים במהלך כל החישוב. למה זה טוב? תכף אתן דוגמא פשוטה.
ראשית, בואו נבין מה באלגוריתם ה”יעיל” מכריח את השברים להיכנס למשחק. בואו נסתכל במטריצה הבאה:
\( \left[\begin{array}{ccc} 2 & 1 & 3\\ 3 & 1 & 6\\ 4 & 1 & 8 \end{array}\right] \)
האלגוריתם מבוסס על כך שאנחנו עוברים עמודה-עמודה, ובכל פעם אנחנו משתמשים באיבר כלשהו מהעמודה כדי לאפס את כל האיברים שאחריו בעמודה - זה הופך את חישוב הדטרמיננטה לפשוט כי כאשר מפתחים את הדטרמיננטה לפי עמודות, יוצא שבכל פעם יש רק איבר אחד בכל עמודה שצריך לפתח לפיו, ולכן מה שהוא בדרך כלל חישוב רקורסיבי מסובך שמתפצל להרבה מקרים לא מתפצל בכלל. כדי לאפס איברים בעמודה משתמשים בתוצאה המרהיבה לפיה אם לוקחים שורה במטריצה, ומחברים אותה עם שורה אחרת כשהיא מוכפלת בסקלר כלשהו, התוצאה היא בעלת אותה דטרמיננטה כמו המטריצה המקורית. למשל, אם במטריצה לעיל אני אקח את השורה הראשונה, אכפול אותה ב-\( -2 \) ואחבר לשורה האחרונה, אני אקבל את המטריצה
\( \left[\begin{array}{ccc} 2 & 1 & 3\\ 3 & 1 & 6\\ 0 & -1 & 2 \end{array}\right] \)
והמטריצה הזו היא בעלת אותה דטרמיננטה בדיוק כמו המטריצה שהתחלתי ממנה, והנה - קיבלתי אפס בעמודה הראשונה. הבעיה היא שכדי להיפטר מה-3 שבשורה האמצעית אני צריך לכפול את השורה הראשונה במשהו שאינו מספר שלם - ב-\( -\frac{3}{2} \), מה שמוביל אותי למטריצה הבאה:
\( \left[\begin{array}{ccc} 2 & 1 & 3\\ 0 & -\frac{1}{2} & \frac{3}{2}\\ 0 & -1 & 2 \end{array}\right] \)
וזהו, עכשיו יש לי שברים בתוך המטריצה, למרות שהתחלתי עם מטריצה שכולה מספרים שלמים ולמרות שגם התוצאה תהיה מספר שלם. כדי לראות שהתוצאה תהיה מספר שלם מספיק להיזכר באופן שבו דטרמיננטה מוגדרת רקורסיבית: אנחנו מקבלים סכום של איברים שכל אחד מהם שייך למטריצה (כלומר הוא מספר שלם) שמוכפל ב-\( \pm1 \) (מספר שלם) וזה מוכפל בדטרמיננטה של תת-מטריצות שמתקבלות מהמטריצה המקורית על ידי מחיקת שורות ועמודות (ולכן גם בהן יש מספרים שלמים). לי הסיטואציה הזו טיפה מזכירה את המקרה של ה-Casus irreducibilis בפתרון משוואות ממעלה שלישית (יש לי פוסט על זה): שם העניין הוא שיש לנו משוואה פולינומית עם מקדמים שהם מספרים ממשיים, והפתרונות שלה כולם מספרים ממשיים, אבל בלי לעבור “דרך” המספרים המרוכבים (כלומר, להוציא שורש ריבועי למספר שלילי) לא נוכל לכתוב את הפתרונות הללו באמצעות נוסחה. אלא שיש הבדל מהותי - במקרה ההוא חייבים לעבור דרך המרוכבים (יש לזה הוכחה די מרהיבה באמצעות תורת גלואה שאני מראה בפוסט ההוא) אבל במקרה שלנו זה לא הכרחי. יש דרכים אחרות לחשב את הדטרמיננטה גם בלי להכניס שברים למשחק. כבר ראינו אחת: לחשב את הדטרמיננטה באמצעות האלגוריתם הרקורסיבי, אלא שזה כאמור פשוט לא יעיל ולכן בפוסט הזה אני אציג את מה שנקרא “האלגוריתם של ברייס” (Bareiss) ששומר על הכל שלם והנזק שהוא גורם לזמן הריצה הוא לא כזה משמעותי (אבל יש נזק; זה בהחלט אלגוריתם שבעייתי בשלל סיטואציות ולכן טוב שמכירים אותו וגם שיטות אחרות).
בואו נדבר שניה על המוטיבציה האישית שלי לכתוב את הפוסט כדי להבין למה לא כדאי לסמוך בעיניים עצומות על השיטות הקיימות ולמה לפעמים חשוב לדבוק במספרים שלמים במקום בשברים. נסיבות אלו ואחרות הובילו אותי לחשב דטרמיננטה של מטריצה \( 5\times5 \), שנותנת לנו את הנפח של טטרהדרון (פירמידה עם בסיס משולש) שאורכי הצלעות שלו נתונים על ידי \( d_{ij} \) עבור \( 1\le i,j\le4 \). הנפח \( V \) נתון באמצעות
\( 288V^{2}=\left|\begin{array}{ccccc} 0 & 1 & 1 & 1 & 1\\ 1 & 0 & d_{12}^{2} & d_{13}^{2} & d_{14}^{2}\\ 1 & d_{21}^{2} & 0 & d_{23}^{2} & d_{24}^{2}\\ 1 & d_{31}^{2} & d_{32}^{2} & 0 & d_{34}^{2}\\ 1 & d_{41}^{2} & d_{42}^{2} & d_{43}^{2} & 0 \end{array}\right| \)
הדטרמיננטה הזו נקראת דטרמיננטת קיילי-מנגר והיא ראויה לפוסט בפני עצמה (היא עושה יותר מאשר לחשב נפח של טטרהדרון) והאמת העצובה היא שהדרך הכי פשוטה לחשב אותה היא פשוט להשתמש בנוסחה הרקורסיבית - זו בסך הכל מטריצת \( 5\times5 \), הרקורסיה תסתיים כמעט מייד, ואין צורך באלגוריתם שאני הולך להציג - אבל את זה נזכרתי לנסות רק אחרי שכבר מימשתי את האלגוריתם (כי על מי אני עובד, כל כך התלהבתי מהאלגוריתם שהייתי חייב לממש אותו בכל מקרה). אבל למה היא גרמה לי להגיע אל אלגוריתם בריס מלכתחילה? ובכן, כי עשיתי את הטעות של לחשב אותה באמצעות numpy. מבלי להיכנס לפרטים, המטרה שלי הייתה למצוא אורכי צלעות שעבורם הדטרמיננטה תהיה שווה בדיוק ל-256. עכשיו, בואו נניח שיש לנו מספר שלם \( n \) כלשהו ואנחנו מסתכלים על הדטרמיננטה
\( \left|\begin{array}{ccccc} 0 & 1 & 1 & 1 & 1\\ 1 & 0 & 1 & n^{2} & n^{2}\\ 1 & 1 & 0 & \left(n-1\right)^{2} & \left(n-1\right)^{2}\\ 1 & n^{2} & \left(n-1\right)^{2} & 0 & 4\\ 1 & n^{2} & \left(n-1\right)^{2} & 4 & 0 \end{array}\right| \)
שהיא מקרה פרטי של הדטרמיננטה לעיל עבור סדרת האורכים \( 1,n,n,n-1,n-1,2 \). אם מחשבים את הדטרמיננטה בצורה סימבולית (כלומר, פשוט פותחים את הביטוי, מקבלים פולינום ב-\( n \) ומפשטים אותו) מקבלים \( -32 \), כלומר הערך של הדטרמיננטה הזו לא תלוי ב-\( n \). למי שתוהים איך עושים את זה במחשב בלי לעבור את התהליך המהנה של לחשב ידנית דטרמיננטה \( 5\times5 \), אפשר לעשות את זה עם ספריית הפייתון sympy:
import sympy as sp
n = sp.symbols('n')
matrix = sp.Matrix([
[0, 1, 1, 1, 1],
[1, 0, 1, n**2, n**2],
[1, 1, 0, (n-1)**2, (n-1)**2],
[1, n**2, (n-1)**2, 0, 4],
[1, n**2, (n-1)**2, 4, 0]
])
determinant = matrix.det()
print(sp.simplify(determinant))טוב ויפה, אלא שאני לא השתמשתי ב-sympy כי בדקתי שלל דטרמיננטות משלל סוגים ולא רק מהמבנה הספציפי הזה שאותו מצאתי אחר כך. ספציפית, אני בדקתי את הדטרמיננטה שמתקבלת עבור הערך \( n=524,283 \) וחישבתי אותו בעזרת הספריה הסטנדרטית לחישובים נומריים בפייתון - numpy. ומה ש-numpy נתנה הוא את התוצאה \( 255.99999991524982 \).
במבט ראשון, התוצאה הזו היא “זה 256, פשוט numpy משתמשת בייצוג שברים עם נקודה צפה ולכן יש אי דיוקים קטנים”. אבל לא! זו דרך טובה לעבוד על עצמנו! כי אם אני מכניס \( n=524,283 \) ואז עוד מעלה דברים בריבוע, אני עובד עם מספרי ענק ולכן יש לי שגיאות דיוק גדולות, וחישובים שקשורים במטריצות יכולים להיות מאוד רגישים לשגיאות דיוק, וזה מתנפח ומתנפח. התוצאה, כאמור, הייתה אמורה להיות \( -32 \); זה שהגענו אל כמעט 256 זה סתם מקרה (על ידי ערכים שונים של \( n \) אפשר להגיע קרוב לשלל מספרים שלמים לא קשורים, פשוט 256 היה מה שחיפשנו). הנה קוד פייתון שמשתמש ב-numpy ומאפשר לראות מה הולך פה. עבור \( n=1 \) החישוב מדויק כמעט לחלוטין, וכך גם עבור \( n=10,000 \); אבל כשהוא נשבר, הוא נשבר חזק.
import numpy as np
for n in [1, 10000, 524283]:
A = np.array([[0, 1, 1, 1, 1],
[1, 0, 1, n**2, n**2],
[1, 1, 0, (n-1)**2, (n-1)**2],
[1, n**2, (n-1)**2, 0, 4],
[1, n**2, (n-1)**2, 4, 0]]
)
print(np.linalg.det(A))הבעיה פה היא כאמור ש-numpy מייצג שברים בשיטת הנקודה הצפה, ולכן יש לו רמת דיוק מוגבלת, בעוד שבפייתון יש רמת דיוק בלתי מוגבלת לעבודה עם מספרים שלמים. אז מכיוון שחשדתי ש-numpy מרמה אותי החלטתי שאני לא צריך להתעצל ולהסתמך עליו או אפילו לחפש ספריה אחרת שמבצעת את החישוב (כי אולי גם היא תרמה אותי?) אלא פשוט לממש את זה בעצמי. עכשיו, הייתי יכול להשתמש באלגוריתם הרקורסיבי הנאיבי; וגם הייתי יכול להשתמש בשיטה מבוססת השברים שראינו, ופשוט להשתמש בספריה frac בפייתון שמאפשרת ייצוג של שברים עם דיוק לא מוגבל. אבל באותו הרגע הדבר הראשון שעלה לי לראש היה לחפש “אלגוריתם לחישוב דטרמיננטה שלא משתמש בשברים” ואני שמח שזה קרה כי אלגוריתם בריס הוא די מגניב גם אם הייתי יכול להסתדר בלעדיו.
אז בואו נדבר עליו סוף סוף.
חלק שני שבו אנו רואים את מה שלשמו התכנסנו פה
קודם כל, מה אלגוריתם בריס לא עושה: הוא לא מעביר את המטריצה \( A \) שלנו למטריצה אחרת \( B \) כך ש-\( \left|A\right|=\left|B\right| \), כמו האלגוריתם היעיל שהצגתי. זה לא הולך לקרות. אני כן הולך להפוך את \( A \) למטריצה אחרת \( B \), אבל יהיו להן דטרמיננטות שונות לגמרי והרעיון הוא שבסיום האלגוריתם, \( \left[B\right]_{nn} \) (הכניסה הימנית-תחתונה של \( B \)) תהיה שווה ל-\( \left|A\right| \). יותר מכך - אין בעצם סיבה להציג את האלגוריתם בתור אוסף של שינויים של המטריצה \( A \) (ובמאמר שלו עליו אני מתבסס כאן, Sylvester’s Identity and Multistep Integer Preserving Gaussian Elimination, בריס באמת לא מציג אותו ככה). אפשר לחשוב על האלגוריתם גם ככה: אם בהתחלה האיבר הכללי של \( A \) הוא \( a_{ij} \) (עבור \( 1\le i,j\le n \)) אז האלגוריתם מייצר סדרה של איברים \( a_{ij}^{\left(0\right)},a_{ij}^{\left(1\right)},\ldots a_{ij}^{\left(n-1\right)} \) כך שבסופו של דבר יוצא ש-\( a_{nn}^{\left(n-1\right)}=\left|A\right| \). ה”חזקה” של האיברים היא לא חזקה אלא אינדקס של השלב באלגוריתם שבו אנחנו נמצאים כרגע, והאתחול הוא \( a_{ij}^{\left(0\right)}=a_{ij} \). בפועל, כשמממשים את האלגוריתם, הכי נוח באמת לבצע שינויים מקומיים ב-\( A \) כדי לשמור את המספרים של השלב הבא - זה חוסך זיכרון.
אז איך עובד האלגוריתם? אין שום טעם להציג אותו כרגע כי לא נבין כלום ממה שהולך שם. זה לא הולך למנוע ממני להציג אותו בכל מקרה כי זה לפחות יוצר עניין - אני קודם כל מימשתי את האלגוריתם ורק אז תהיתי למה בעצם זה עובד. הרעיון, כאמור, הוא לעבוד בשלבים, כשבשלב מספר \( k \) אנחנו מעדכנים את המספרים \( a_{ij}^{\left(k\right)} \) שעדיין רלוונטיים לנו (כמו בחישוב דטרמיננטה רגיל, ככל שמתקדמים בשלבים כך פחות ופחות איברים הופכים לרלוונטיים לנו.
הנה האלגוריתם:
- אתחל את המשתנים \( a_{ij}^{\left(0\right)}=a_{ij} \) לכל \( 1\le i,j\le n \) ואת המשתנה המיוחד \( a_{00}^{\left(-1\right)}=1 \).
- לכל \( k=1,2,\ldots,n-1 \):
- לכל \( k+1\le i,j\le n \), קבעו \( a_{ij}^{\left(k\right)}\leftarrow\frac{a_{ij}^{\left(k-1\right)}a_{kk}^{\left(k-1\right)}-a_{ik}^{\left(k-1\right)}a_{kj}^{\left(k-1\right)}}{a_{k-1,k-1}^{\left(k-2\right)}} \)
בקיצור, האלגוריתם די דומה למה שקורה באלגוריתם הרגיל - גם פה פשוט מעדכנים באופן סדרתי כניסות עם כל מני מכפלות וחיסורים ו… רגע… האם זה רק אני או שיש סימן חילוק ענקי באמצע האלגוריתם? איך בדיוק פתרנו את הבעיה? ובכן, הרעיון הוא שבסימן החילוק שמופיע שם, המכנה מחלק את המונה בלי שארית. במילים אחרות - אנחנו צריכים לבצע פעולת חילוק, אבל פעולת חילוק שלמים שיכולה להיות מדויקת, ובשום שלב אנחנו לא צריכים לעבוד עם מספר שהוא שבר (אם תקראו את הפוסט עד כה תראו שנמנעתי בזהירות מלטעון דברים מופרכים כמו “האלגוריתם לא מבצע חילוק”).
אוקיי, אבל למה זה עובד? במאמר שלו שעליו אני מתבסס כאן, בריס מביא טיעון לא טריוויאלי, שאולי היה אפשר להחליף בטיעונים אלמנטריים יותר (יש לבריס מאמר מוקדם יותר שבו יש טיעונים כאלו) אבל אני דווקא מאוד אוהב את הטכניקה שלו, אז בואו נראה אותה במלואה.
חלק שלישי שבו קורים קסמים עם מטריצות בלוקים
באלגוריתם שהצגתי, \( a_{ij}^{\left(k\right)} \) חושב באמצעות תרגיל מפוקפק כלשהו (שבמסגרתו הבטחתי שתתבצע חלוקה ללא שארית). זה לא ממש עוזר לנו להבין מה הולך פה, אז מה שנרצה לעשות הוא למצוא דרך טובה יותר להגדיר את הערכים של אותם \( a_{ij}^{\left(k\right)} \) מלכתחילה. ההגדרה תיראה די משונה אם פשוט אציג אותה, אז בואו נוכיח תוצאה כללית כלשהי קודם.
משהו שלא דיברתי עליו בכלל עד כה בפוסטים הללו הוא שדרך נוחה להתמודד עם מטריצות היא לפעמים לחלק אותן לבלוקים. כל מטריצה ריבועית \( A \) מסדר \( n\times n \) אפשר להציג בתור \( A=\left(\begin{array}{cc} A_{11} & A_{12}\\ A_{21} & A_{22} \end{array}\right) \) כך ש-\( A_{11} \) היא מטריצה מסדר \( k\times k \) שכוללת את מה שמקבלים אם מעתיקים מ-\( A \) את כל הכניסות מהצורה \( a_{ij} \) עבור \( 1\le i,j\le k \). גם שאר המטריצות נקבעות בצורה דומה: \( A_{12} \) תהיה מסדר \( k\times\left(n-k\right) \), \( A_{21} \) תהיה מסדר \( \left(n-k\right)\times k \) ו-\( A_{22} \) תהיה מסדר \( \left(n-k\right)\times\left(n-k\right) \), אבל בואו נתמקד לרגע ב-\( A_{11} \). זו המטריצה
\( A_{11}=\left(\begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1k}\\ a_{21} & a_{22} & & a_{2k}\\ \vdots & & \ddots & \vdots\\ a_{k1} & a_{k2} & \cdots & a_{kk} \end{array}\right) \)
מטריצה כזו, שמתקבלת מ-\( A \) על ידי לקיחת \( k \) השורות והעמודות הראשונות ומחיקת כל היתר, נקראת לפעמים המינור העיקרי (principle) מסדר \( k \) של המטריצה. אני הולך להניח שהוא הפיך, כלומר \( \left|A_{11}\right|\ne0 \); בהמשך נראה מה קורה כשזה לא המצב (רמז: האלגוריתם שהצגתי קודם עדיין לא שלם). אם הוא לא הפיך, זה נותן לנו פירוק פשוט יחסית של \( A \) למכפלה של שתי מטריצות בלוקים אלכסוניות:
\( A=\left(\begin{array}{cc} A_{11} & A_{12}\\ A_{21} & A_{22} \end{array}\right)=\left(\begin{array}{cc} A_{11} & 0\\ A_{21} & I \end{array}\right)\left(\begin{array}{cc} I & A_{11}^{-1}A_{12}\\ 0 & A_{22}-A_{21}A_{11}^{-1}A_{12} \end{array}\right) \)
מה הולך כאן? במבט ראשון זה נראה מבעית, אבל בעצם יש כאן שתי שאלות פשוטות יחסית:
- כשכופלים מטריצות בלוקים, האם חוקי כפל המטריצות הרגילים עדיין תקפים כמו קודם?
- אם מניחים שכן, האם החישוב של המכפלה למעלה באמת מחזיר את \( A \)?
בואו נענה קודם לשאלה השניה - כמובן. כדי לוודא שאין בעיה עם זה, אני אחשב את כל ארבע הכניסות
- הכניסה של \( A_{11} \): מתקבלת מכפל השורה הראשונה בעמודה הראשונה, כלומר \( A_{11}\cdot I+0\cdot0=A_{11} \).
- הכניסה של \( A_{12} \): מתקבלת מכפל השורה הראשונה בעמודה השניה, כלומר \( A_{11}\cdot\left(A_{11}^{-1}A_{12}\right)+0\cdot\left(A_{22}-A_{21}A_{11}^{-1}A_{12}\right)=A_{12} \)
- הכניסה של \( A_{21} \): מתקבלת מכפל השורה השניה בעמודה הראשונה, כלומר \( A_{21}\cdot I+I\cdot0=A_{21} \).
- הכניסה של \( A_{22} \): מתקבלת מכפלה השורה השניה בעמודה השניה, כלומר \( A_{21}\cdot A_{11}^{-1}A_{12}+I\cdot A_{22}-A_{21}A_{11}^{-1}A_{12}=A_{22} \)
כל זה אולי נראה קצת מהונדס מדי כדי שהכל יעבוד, אבל זה כמובן עובד. השאלה היותר מהותית היא למה בכלל מותר לכפול מטריצות בלוקים - וכמו הרבה דברים באלגברה לינארית בסיסית זו תוצאה שאני גם מניח שכולנו מכירים והיא גם די טכנית להוכחה, אז אני לא אוכיח אותה כאן אבל אין כאן חוכמה מיוחדת - אם אתם סקפטיים, נסו להוכיח אותה על שתי מטריצות קטנות קונקרטיות ותראו מה קורה.
עכשיו, מה שנחמד בפירוק המוזר הזה של \( A \) למכפלה של שתי מטריצות בלוקים הוא שזה ממשיך לעבוד כשלוקחים דטרמיננטה. כזכור, דטרמיננטה של מכפלה היא מכפלת הדטרמיננטות; וכשיש לנו מטריצת בלוקים \( \left(\begin{array}{cc} A & 0\\ C & D \end{array}\right) \) הדטרמיננטה שלה היא \( \left|A\right|\left|D\right| \) (ובאופן דומה כשה-0 הוא במקום \( C \)), כך שנקבל:
\( \left|A\right|=\left|\left(\begin{array}{cc} A_{11} & 0\\ A_{21} & I \end{array}\right)\right|\left|\left(\begin{array}{cc} I & A_{11}^{-1}A_{12}\\ 0 & A_{22}-A_{21}A_{11}^{-1}A_{12} \end{array}\right)\right|= \)
\( \left|A_{11}\right|\left|A_{22}-A_{21}A_{11}^{-1}A_{12}\right| \)
עכשיו משתמשים בטריק: כופלים את שני האגפים בסקלר \( \left|A_{11}\right|^{n-k-1} \) ומקבלים
\( \left|A_{11}\right|^{n-k-1}\left|A\right|=\left|\left|A_{11}\right|\left(A_{22}-A_{21}A_{11}^{-1}A_{12}\right)\right| \)
זה יכול להיות די מבלבל, הנה מה שקרה פה: באופן כללי, אם יש לי מטריצה \( B \) ואני כופל שורה כלשהי שלה בסקלר, זה מכפיל את כל הדטרמיננטה באותו סקלר. כדי להוכיח את זה אפשר להשתמש באותה טכניקה שראינו בפוסט הקודם: כפל שורה בסקלר זה כמו כפל במטריצת יחידה \( I \) שאחד מה-1-ים על האלכסון שלה הוחלף בסקלר הבודד \( \lambda \), והדטרמיננטה של מטריצה כזו היא \( \lambda \) כי היא הרי מכפלת הערכים על האלכסון.
עכשיו, מה קורה אם כופלים את כל השורות של \( B \) באותו סקלר, כלומר אם פשוט כפלנו את כל המטריצה \( B \) בסקלר הזה? אם \( B \) היא מטריצה מסדר \( n\times n \) זה אומר שכפלנו את הדטרמיננטה \( n \) פעמים בסקלר \( \lambda \), אז יש לנו את המשוואה \( \left|\lambda B\right|=\lambda^{n}\left|B\right| \).
אם נחזור עכשיו למה שעשינו למעלה - יש לנו את המטריצה \( B=A_{22}-A_{21}A_{11}^{-1}A_{12} \). הסדר של המטריצה הזו הוא \( \left(n-k\right)\times\left(n-k\right) \) ולכן \( \lambda^{n-k}\left|B\right|=\left|\lambda B\right| \). אצלנו \( \lambda=\left|A_{11}\right| \), ואנחנו מקבלים אותו בחזקת \( n-k \) כי בהתחלה הוא הופיע פעם אחת באגף ימין ואז כפלנו ב-\( \lambda^{n-k-1} \). זה מסביר את המעבר הזה.
חלק רביעי, שבו אנו חוזים בשובם של המספרים מהאלגוריתם
עכשיו הגענו לקאץ’: אני אראה שהאיברים של המטריצה \( \left|A_{11}\right|\left(A_{22}-A_{21}A_{11}^{-1}A_{12}\right) \) שאליה הגענו הם בעצם ה-\( a_{ij}^{\left(k\right)} \) שמופיעים בשלבי הביניים של האלגוריתם. ראשית, בואו נגדיר את ה-\( a_{ij}^{\left(k\right)} \) הללו במפורש, מה שעד כה נמנעתי מלעשות:
\( a_{ij}^{\left(k\right)}=\left|\begin{array}{ccccc} a_{11} & a_{12} & \cdots & a_{1k} & a_{1j}\\ a_{21} & a_{22} & \cdots & a_{2k} & a_{2j}\\ \cdots & \cdots & \cdots & \cdots & \cdots\\ a_{k1} & a_{k2} & \cdots & a_{kk} & a_{kj}\\ a_{i1} & a_{i1} & \cdots & a_{ik} & a_{ij} \end{array}\right| \) עבור \( k<i,j\le n \)
מה קורה פה? ראשית לוקחים את המטריצה \( A_{11} \) שמתקבלת מ-\( A \) המקורית על ידי לקיחת \( k \) השורות והעמודות הראשונות. עכשיו, מוסיפים שורה חדשה - את שורה מס’ \( i \), כאשר \( k<i \) כלומר זו אחת מהשורות שלא הופיעו במטריצה קודם. אנחנו לא מוסיפים את כל השורה אלא רק את \( k \) העמודות הראשונות שלה. אחר כך אנחנו מוסיפים עמודה \( j \) כך ש-\( k<j \) וגם כאן - לא את כל העמודה, רק את \( k \) השורות הראשונות שלה. לסיום, אחרי שהוספנו את השורה והעמודה עוד נותר מקום ריק אחד בקצה הימני-תחתון של המטריצה, ושם אנחנו שמים את \( ij \).
למה שנגדיר ככה? במאמר המקורי שלו בריס מגיע אל המטריצה הזו בדרך שונה וטכנית יותר - אני דווקא מעדיף את הדרך שבה נוקטים כאן, למרות שהיא נראית כמו קסם (במתמטיקה אין קסמים; כל הגדרה קסומה מגיעה אחרי שמישהו שבר את הראש על הבעיה הרבה זמן וניסה הרבה דברים).
עכשיו, הדטרמיננטה הזו נראית קצת מפחידה אבל כבר ראינו דרך להתמודד איתה - אם שמים לב שיש לנו כאן דטרמיננטה של מטריצת בלוקים. הבלוקים הם בדיוק ארבעת הדברים שתיארתי: המטריצה \( A_{11} \) שממנה מתחילים, השורה שהוספתי למטה, העמודה שהוספתי משמאל, וה-\( a_{ij} \) שהוספתי למטה. כלומר יש לנו כאן את מטריצת הבלוקים
\( B=\left(\begin{array}{cc} B_{11} & B_{12}\\ B_{21} & B_{22} \end{array}\right) \)
כך ש-\( B_{11}=A_{11} \), \( B_{12}=\left(a_{1j},\ldots,a_{kj}\right)^{t},B_{21}=\left(a_{i1},\ldots,a_{ik}\right) \) ו-\( B_{22}=\left(a_{ij}\right) \).
ראינו את הנוסחה
\( \left|B\right|=\left|B_{11}\right|\left|B_{22}-B_{21}B_{11}^{-1}B_{12}\right| \)
כאן \( \left|B\right|=a_{ij}^{\left(k\right)} \) כי \( a_{ij}^{\left(k\right)} \) הוגדר בתור הדטרמיננטה של המטריצה הזו. קצת יותר מעניין מה זה \( B_{21}B_{11}^{-1}B_{12} \) - זה כפל של מטריצה (\( A_{11}^{-1} \)) במטריצת עמודה מימין ומטריצת שורה משמאל - התוצאה היא סקלר, וחישוב ישיר על פי ההגדרה של כפל מטריצות נותן ש-
\( B_{21}B_{11}^{-1}B_{12}=\sum_{r=1}^{k}\sum_{s=1}^{k}a_{ir}\left[A_{11}^{-1}\right]_{rs}a_{sj} \)
ולכן משילוב כל הדברים הללו אנחנו מקבלים
\( a_{ij}^{\left(k\right)}=\left|A_{11}\right|\left(a_{ij}-\sum_{r=1}^{k}\sum_{s=1}^{k}a_{ir}\left[A_{11}^{-1}\right]_{rs}a_{sj}\right) \)
מה שיש לנו באגף שמאל הוא בדיוק את האיבר ה-\( ij \) של המטריצה \( \left|A_{11}\right|\left(A_{22}-A_{21}A_{11}^{-1}A_{12}\right) \), אם כי צריך קצת להיזהר עם האינדקסים כדי לראות את זה כי המטריצה הזו מתחילה לא מהאינדקס \( 1,1 \) אלא מהאינדקס \( k+1,k+1 \).
איך זה עוזר לנו? ובכן, קודם כל שימו לב שכמו שאמרתי קודם - המינורים העיקריים של \( A \) מתקבלים ככה. ממש על פי ההגדרה, מתקיים \( \left|A_{11}\right|=a_{kk}^{\left(k-1\right)} \), כי \( A_{11} \) מתקבלת מלקיחת \( k \) השורות והעמודות הראשונות ואילו המטריצה שבה משתמשים לחישוב \( a_{kk}^{\left(k-1\right)} \) מתקבלת מלקיחת \( k-1 \) השורות והעמודות הראשונות ואז הוספה אליהן של השורה והעמודה הבאות בתור.
לכן אני יכול ללכת אל הנוסחה שהוכחתי קודם:
\( \left|A_{11}\right|^{n-k-1}\left|A\right|=\left|\left|A_{11}\right|\left(A_{22}-A_{21}A_{11}^{-1}A_{12}\right)\right| \)
ולהציב בתוכה את \( \left|A_{11}\right|=a_{kk}^{\left(k-1\right)} \) מצד אחד, ואילו בצד השני יש לנו כאמור את המטריצה \( \left|A_{11}\right|\left(A_{22}-A_{21}A_{11}^{-1}A_{12}\right) \) שהאיברים שלה הם בדיוק ה-\( a_{ij}^{\left(k\right)} \) עבור \( k<i,j\le n \), ולכן קיבלנו את הנוסחה
\( \left|A\right|\left[a_{kk}^{\left(k-1\right)}\right]^{n-k-1}=\left|\begin{array}{ccc} a_{k+1,k+1}^{\left(k\right)} & \cdots & a_{k+1,n}^{\left(k\right)}\\ \cdots & \cdots & \cdots\\ a_{n,k+1}^{\left(k\right)} & \cdots & a_{n,n}^{\left(k\right)} \end{array}\right| \)
את הנוסחה הזו הוכחנו לכל מטריצה \( A \) (עם ההנחה שהמינור העיקרי שלה, מה שסימנתי ב-\( A_{11} \), הוא מדטרמיננטה שונה מאפס), אז אפשר להשתמש בה גם עבור המטריצה שאיתה הגדרנו את \( a_{ij}^{\left(k\right)} \), כלומר המטריצה שמורכבת מהמינור העיקרי \( A_{11} \) ואז עוד שורה ועמודה. צריך קצת להיזהר עם האינדקסים: עד עכשיו השתמשתי ב-\( k \) כדי לציין אינדקס כלשהו, \( 1\le k\le n \) שמתאים למינור העיקרי שבו מטפלים באותו הרגע, אבל עכשיו כשאני רוצה להשתמש בנוסחה הזו כדי לטפל ב-\( a_{ij}^{\left(k\right)} \) ולכן יש לי חופש בחירה של המינור העיקרי שלו שאנחנו בוחרים, צריך להשתמש באינדקס אחר - במאמר שלו ברייס משתמש ב-\( l \), ואז מקבלים
\( a_{ij}^{\left(k\right)}\left[a_{ll}^{\left(l-1\right)}\right]^{\left(k+1\right)-l-1}=\left|\begin{array}{cccc} a_{l+1,l+1}^{\left(l\right)} & \cdots & a_{l+1,k}^{\left(l\right)} & a_{l+1,j}^{\left(l\right)}\\ \cdots & \cdots & \cdots & \cdots\\ a_{k,l+1}^{\left(l\right)} & \cdots & a_{k,k}^{\left(l\right)} & a_{k,j}^{\left(l\right)}\\ a_{i,l+1}^{\left(l\right)} & \cdots & a_{i,k}^{\left(l\right)} & a_{i,j}^{\left(l\right)} \end{array}\right| \)
הסיבה לכך שיש \( k+1 \) בחזקה באגף שמאל הוא שזה הסדר של המטריצה שמגדירה את \( a_{ij}^{\left(k\right)} \) (כי לקחנו מטריצה \( k\times k \) והוספנו לה שורה ועמודה). השינוי באיך שאני מציג את האיברים בתוך הדטרמיננטה באגף ימין (עוד שורה ועמודה) הוא רק כדי שיהיה יותר קל להבין מה קורה, בגלל שהמבנה של השורה והעמודה האחרונות במטריצה הוא שונה מאשר עבור יתר הכניסות.
עכשיו, תחת ההנחה שלנו על כך שהדטרמיננטה של המינורים העיקריים היא לא אפס אנחנו מקבלים ש-\( a_{ll}^{\left(l-1\right)} \) הוא לא אפס (כי הוא שווה למינור עיקרי שכזה) ולכן אפשר לחלק בו ולקבל
\( a_{ij}^{\left(k\right)}=\frac{1}{\left[a_{ll}^{\left(l-1\right)}\right]^{k-l}}\left|\begin{array}{cccc} a_{l+1,l+1}^{\left(l\right)} & \cdots & a_{l+1,k}^{\left(l\right)} & a_{l+1,j}^{\left(l\right)}\\ \cdots & \cdots & \cdots & \cdots\\ a_{k,l+1}^{\left(l\right)} & \cdots & a_{k,k}^{\left(l\right)} & a_{k,j}^{\left(l\right)}\\ a_{i,l+1}^{\left(l\right)} & \cdots & a_{i,k}^{\left(l\right)} & a_{i,j}^{\left(l\right)} \end{array}\right| \)
וזה נותן לנו דרך לחשב רקורסיבית את ה-\( a_{ij}^{\left(k\right)} \)-ים! נזכיר את תנאי ההתחלה שלנו:
\( a_{00}^{\left(-1\right)}=1,a_{ij}^{\left(0\right)}=a_{ij} \)
ובנוסחה למעלה, כדי לצמצם כמה שרק ניתן את הגודל של הדטרמיננטה אפשר לבחור \( l=k-1 \) ולקבל
\( a_{ij}^{\left(k\right)}=\frac{1}{a_{k-1,k-1}^{\left(k-2\right)}}\left|\begin{array}{cc} a_{kk}^{\left(k-1\right)} & a_{kj}^{\left(k-1\right)}\\ a_{ik}^{\left(k-1\right)} & a_{ij}^{\left(k-1\right)} \end{array}\right|=\frac{a_{ij}^{\left(k-1\right)}a_{kk}^{\left(k-1\right)}-a_{ik}^{\left(k-1\right)}a_{kj}^{\left(k-1\right)}}{a_{k-1,k-1}^{\left(k-2\right)}} \)
כשהשוויון האחרון נובע פשוט מחישוב ישיר של הדטרמיננטה. זה משלים את ההוכחה, כי קיבלנו בדיוק את הנוסחה שהופיעה באלגוריתם, ועכשיו אנחנו גם מבינים למה החלוקה ב-\( a_{k-1,k-1}^{\left(k-2\right)} \) “מצליחה” ואנחנו מקבלים מספר שלם - כי \( a_{ij}^{\left(k\right)} \) מלכתחילה הוגדר בתור דטרמיננטה של מטריצה שמכילה רק מספרים שלמים (זאת בהנחה שהמטריצה המקורית הכילה רק שלמים) ולכן אם יש שבר ששווה לו, המכנה חייב לחלק את המונה.
חלק חמישי ואחרון, שבו אנו מגיעים לגרסה המלאה של האלגוריתם
יפה, אז ראינו את ההוכחה לכך שהאלגוריתם הבא עובד:
- אתחל את המשתנים \( a_{ij}^{\left(0\right)}=a_{ij} \) לכל \( 1\le i,j\le n \) ואת המשתנה המיוחד \( a_{00}^{\left(-1\right)}=1 \).
- לכל \( k=1,2,\ldots,n-1 \):
- לכל \( k+1\le i,j\le n \), קבעו \( a_{ij}^{\left(k\right)}\leftarrow\frac{a_{ij}^{\left(k-1\right)}a_{kk}^{\left(k-1\right)}-a_{ik}^{\left(k-1\right)}a_{kj}^{\left(k-1\right)}}{a_{k-1,k-1}^{\left(k-2\right)}} \)
עכשיו, מהנוסחה הזו ברור שאם אנחנו בשלב \( k \), אז כל הערכים במטריצה שהשורה או העמודה שלהם היא לכל היותר \( k \) כבר לא מעניינים אותנו, למעט האיבר במקום \( k-1,k-1 \) - אבל גם עבורו, אנחנו נצטרך רק לקרוא את הערך שלו ולא נשנה אותו יותר. זה אומר שאם אנחנו מקבלים את המטריצה כקלט, אז בשביל לחסוך מקום אפשר לשמור את ערכי הביניים בתוכה (אם לא מפריע לנו שהיא תשתנה; אם לא נרצה שהפונקציה תשפיע על המטריצה המקורית יהיה צורך ליצור עותק שלה, אבל גם אז עדיף עותק אחד על סביבות \( n \) עותקים כאלו).
אם נוקטים בגישה הזו, זה מה שהאלגוריתם הופך להיות:
- אתחלו את המשתנה \( A_{00}=1 \)
- לכל \( k=1,2,\ldots,n-1 \):
- לכל \( k+1\le i,j\le n \), קבעו \( A_{ij}\leftarrow\frac{A_{ij}A_{kk}-A_{ik}A_{kj}}{A_{k-1,k-1}} \)
בפועל כשכותבים קוד לא כזה כיף להוסיף משתנה \( A_{00} \) כי צריך קוד מיוחד לטפל בו ולא במטריצה \( A \), ולכן הרבה יותר פשוט לבדוק אם \( k=1 \) ואם כן - פשוט לוותר על פעולת החילוק באיטרציה הזו של האלגוריתם. בקוד שאצרף עוד רגע אפשר יהיה לראות את זה.
יש רק דבר אחד שטרם התייחסתי אליו - מה קורה אם אחד מהמינורים העיקריים הוא עם דטרמיננטה 0? במקרה הזה, החישוב יניב בשלב כלשהו \( A_{kk}=0 \) ואי אפשר יהיה לחלק בו. בסיטואציה הזו עושים דבר דומה למה שהיה באלגוריתם ה”רגיל” לדטרמיננטה - מחפשים בעמודה ה-\( k \)-ית החל מהשורה ה-\( k+1 \) כניסה ששונה מאפס, ואם נמצאה כזו - מחליפים את השורות ומכפילים את הסימן של הדטרמיננטה ב-\( -1 \). זה האלגוריתם המלא (שימו לב שבפייתון // מייצג חילוק בשלמים - כלומר, חילוק שמתעלם מהשארית, אם יש כזו, ומחזיר רק את המנה השלמה):
def bareiss_det(A):
sign = 1
n = len(A)
for k in range(n-1):
if A[k][k] == 0:
for i in range(k+1,n):
if A[i][k] != 0:
A[k], A[i] = A[i], A[k]
sign *= -1
break
if A[k][k] == 0:
return 0
for i in range(k+1, n):
for j in range(k+1, n):
A[i][j] = (A[i][j] * A[k][k] - A[i][k] * A[k][j])
if k != 0:
A[i][j] = A[i][j] // A[k-1][k-1]
A[i][k]=0
return sign * A[-1][-1]נהניתם? התעניינתם? אם תרצו, אתם מוזמנים לתת טיפ:
