6.13. אלגברה לינארית¶

אלגברה לינארית במצלמה היא עבודה עם מטריצות קטנות: סיבוב 3x3 שמשלב דגימת IMU למסגרת העולם, מטריצת כיול שמיישרת עדשה, עדכון שונות-מצב של מסנן קלמן, התאמת פולינום שהמשוואות הנורמליות שלה יוצאות כפתרון לינארי זעיר. תת-המודולים numpy.linalg ו-scipy.linalg מכסים בדיוק את הסקאלה הזו.

הפונקציות נמצאות בשני מודולים. פעולות מכפלת מטריצות נמצאות ברמה העליונה של numpy; פירוקים והופכי המטריצה נמצאים תחת numpy.linalg; פותרי המערכות הלינאריות הייעודיים נמצאים תחת scipy.linalg. כפל מטריצות 2-על-2, למשל:

from ulab import numpy as np

A = np.array([[1, 2], [3, 4]], dtype=np.float)
B = np.array([[5, 6], [7, 8]], dtype=np.float)
np.dot(A, B)
# array([[19.0, 22.0],
#        [43.0, 50.0]])

6.13.1. מה זמין¶

dot() – מכפלת מטריצות או וקטורים.
cross() – מכפלה וקטורית תלת-ממדית.
trace() – סכום האלכסון.
inv() – הופכי מטריצה.
det() – דטרמיננטה.
cholesky() – פירוק צ’ולסקי (קלט סימטרי חיובי-מובהק).
eig() – ערכים עצמיים ווקטורים עצמיים של מטריצה ממשית סימטרית.
norm() – נורמת-2 של וקטור או מטריצה.
qr() – פירוק QR עם mode='reduced' (ברירת מחדל) או mode='complete'.
solve_triangular() – פתרון A @ x = b כאשר A משולשת.
cho_solve() – פתרון A @ x = b בהינתן גורם צ’ולסקי של A.

6.13.2. dot, cross, trace¶

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

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

np.dot(a, b)             # 32.0 (scalar product)
np.cross(a, b)           # array([-3.0, 6.0, -3.0])

m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
np.dot(m, a)             # matrix-vector product
np.dot(m, m)             # matrix-matrix product
np.trace(m)              # 1 + 5 + 9 = 15.0

תוצאת dot() היא תמיד מסוג dtype float. cross() היא המכפלה הווקטורית של שני וקטורים תלת-ממדיים, ו-trace() היא סכום האלכסון הראשי של מטריצה ריבועית.

6.13.3. inv ו-det¶

m = np.array([[1, 2, 3, 4],
              [4, 5, 6, 4],
              [7, 9, 9, 4],
              [3, 4, 5, 6]])
print(np.linalg.inv(m))
print(np.linalg.det(m))

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

הדטרמיננטה משתמשת באותו חילוץ – זמן הריצה זהה במהותו להופכי.

כאשר מטרת היישום היא לפתור מערכת לינארית, אין להפוך ולכפול – עדיף להשתמש בפותרים הייעודיים שלהלן. שניהם מהירים יותר ומתנהגים טוב יותר מספרית.

6.13.4. cholesky¶

עבור מטריצה A סימטרית חיובית-מובהקת, cholesky() מחזירה מטריצה משולשת-תחתית L כך ש-A = L @ L.T

a = np.array([[25, 15, -5],
              [15, 18,  0],
              [-5,  0, 11]])
L = np.linalg.cholesky(a)

אם הקלט אינו חיובי-מובהק או אינו סימטרי, מועלית ValueError.

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

6.13.5. eig¶

eig() פועלת רק על מטריצות ממשיות סימטריות. מטריצות לא-סימטריות מעלות ValueError. היא מחזירה זוג (eigenvalues, eigenvectors)

a = np.array([[1, 2, 1, 4],
              [2, 5, 3, 5],
              [1, 3, 6, 1],
              [4, 5, 1, 7]], dtype=np.uint8)
x, y = np.linalg.eig(a)

הערות:

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

6.13.6. norm¶

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

v = np.array([1, 2, 3, 4, 5])
m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

np.linalg.norm(v)        # 7.416...
np.linalg.norm(m)        # 16.881...

מילת המפתח האופציונלית axis= לוקחת את הנורמה לאורך ציר יחיד במקום על פני המערך כולו.

6.13.7. qr¶

qr() מפרקת מטריצה מלבנית A (צורה (M, N)) ל-Q אורתונורמלי ול-R משולש-עליון כך ש-A == Q @ R

A = np.arange(6).reshape((3, 2))

q, r = np.linalg.qr(A)
# mode='reduced' (default): q is (3, 2), r is (2, 2)

q, r = np.linalg.qr(A, mode='complete')
# q is (3, 3), r is (3, 2)

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

6.13.8. פתרון מערכות¶

שני הפותרים הייעודיים תחת ulab.scipy.linalg מהירים ומדויקים יותר מ-np.dot(np.linalg.inv(A), b):

solve_triangular(a, b, lower=False)() – פותר a @ x = b בהנחה ש-a משולשת:

A = np.array([[3, 0, 0, 0],
              [2, 1, 0, 0],
              [1, 0, 1, 0],
              [1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])
x = sp.linalg.solve_triangular(A, b, lower=True)

cho_solve(L, b)() – בהינתן גורם צ’ולסקי L, פותר A @ x = b כאשר A = L @ L.T
```
L = np.linalg.cholesky(A)
x = sp.linalg.cho_solve(L, b)
```

השתמשו באלה במקום בהיפוך בכל פעם שמבנה A מאפשר זאת – הם חוסכים עבודת חילוץ וגם את ההופכי המפורש.

6.13.9. פתרון מערכת לינארית קטנה¶

A = np.array([[3, 0, 1, 1],
              [0, 1, 0, 2],
              [1, 0, 1, 1],
              [1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])

x = np.dot(np.linalg.inv(A), b)
print(x)
print(np.dot(A, x))         # should equal b

ניתן לבטא את אותה בעיה על ידי פירוק A וקריאה לפותר המתאים – מהיר ומדויק יותר כאשר ל-A יש את המבנה הנכון.

לעיון מלא ברמת הארגומנטים, ראו numpy.linalg — שגרות אלגברה לינארית ו-scipy.linalg — שגרות אלגברה לינארית.