SciPy 的線性代數子套件 linalg 是以 Numpy 的 ndarry 陣列為基礎建構的, 它提供了比 Numpy 更多更完整的線性代數函式庫, 事實上 SciPy 可視為 Numpy 的擴充集 (superset), 但 Numpy 基於歷史緣故仍保有核心的線性代數函式.
參考書籍 :
使用 SciPy 線性代數子套件 linalg 的方式之一是匯入所有 linalg 的函數 :
from scipy.linalg import *
此方式很方便, 因為可以直接呼叫此子套件的函數, 但要注意命名空間汙染問題.
比較安全的方式是 :
from scipy import linalg
然後用 linalg 去呼叫旗下的函數, 例如 linalg.inv(A).
常用函數如下表 :
| scipy.linalg 常用函數 | 說明 |
| det(A) | 傳回方陣 A 的行列式值 |
| inv(A) | 傳回方陣 A 的反矩陣 |
| solve(A, b) | 傳回線性方程組 Ax=b 的變數向量 x 之解 (串列) |
| eig(A) | 傳回方陣 A 的特徵值 (eigen value) 與特徵向量 (eigen vector) |
| qr(A) | 將方陣 A 進行 QR 分解, 傳回正交矩陣 Q 與上三角矩陣 (元組) |
| lu(A) | 將方陣分解為置換矩陣 P, 上三角矩陣 U, 下三角矩陣 L (元組) |
參考 scipy.linalg 子套件的說明文件 :
1. 用 det() 函數求方陣的行列式 (determinant) 值 :
一個二階方陣 A 的行列式的算法為 :
即左上右下對角線元素相乘後相加, 再減所有右上左下對角線元素之乘積. 同樣的算法在三階方陣 A 的行列式為 :
下面用 linalg 的 det() 函數來計算上面兩個方陣的行列式 :
>>> A=np.array([[1, -1],[1, 2]]) # 二階方陣
>>> linalg.det(A) # 計算方陣 A 的行列式值
3.0
>>> A=np.array([[3, 2, 3], [2, 5, -7], [1, 2, -2]]) # 三階方陣
>>> linalg.det(A) # 計算方陣 A 的行列式值
3.0000000000000013
下面這個三階方陣的行列式值為 0, 例如 :
>>> A=np.array([[2, 2, 2], [2, 2, 2], [2, 2, 2]]) # 三階方陣
>>> A
array([[2, 2, 2],
[2, 2, 2],
[2, 2, 2]])
矩陣的行列式有如下特性 :
- 行列式值為 0 的矩陣為不可逆, 亦即其反矩陣不存在.
- 若方陣 A 任何一列或一行均為 0, 則 det(A)=0.
- 若方陣 A 任何兩列值相同, 則 det(A)=0.
- 若方陣 A 為對角, 上三角或下三角矩陣, 則 det(A)=對角線元素的乘積.
- 任何階的單位矩陣其行列式值必為 0.
- 兩個方陣內積的行列式值等於個別方陣行列式之乘積 det(AB)=det(A)*det(B).
- 方陣 A 轉置後的行列式值等於方陣 A 的行列式值.
- 方陣 A 的 n 次方 (n 次內積) 的行列式值為 A 行列式值的 n 次方.
2. 用 inv() 求方陣的反矩陣 (inverse matrix) :
在線性代數中, 一個 n 階方陣 A 的反矩陣 B 的定義是, A 與 B 的內積為 n 階單元方陣 :
A 的反矩陣 B 記做 :
一個方陣的反矩陣若存在, 此矩陣稱為可逆 (invertible), 可逆的矩陣又稱為非奇異矩陣 (non-singular matrix); 反之, 若方陣的反矩陣不存在即為不可逆, 不可逆的方陣被稱為奇異矩陣 (singular matrix). 反矩陣若存在, 則它是唯一的, 即一個方陣不可能有兩個反矩陣.
線性代數中反矩陣是透過
伴隨矩陣 (adjugate matrix) 計算出來的, 反矩陣與伴隨矩陣存在一個比例關係, 將方陣 A 的伴隨矩陣 A* 除以 A 的行列式 |A| 即可得到反矩陣, 因此行列式值為 0 之矩陣無法求反矩陣 :
SciPy 子套件 linalg 的 inv() 函數可用來求方陣的反矩陣, 例如 :
>>> import numpy as np
>>> from scipy import linalg
>>> A=np.array([[1, -1],[1, 2]]) # 定義一個方陣 A
>>> A
array([[ 1, -1],
[ 1, 2]])
>>> B=linalg.inv(A) # 求 A 的反矩陣 B
>>> B
array([[ 0.66666667, 0.33333333],
[-0.33333333, 0.33333333]])
>>> np.dot(A, B) # 計算 A, B 的內積
array([[1.00000000e+00, 0.00000000e+00],
[1.11022302e-16, 1.00000000e+00]])
>>> np.dot(B, A) # 計算 B, A 的內積
array([[ 1.00000000e+00, -1.11022302e-16],
[ 0.00000000e+00, 1.00000000e+00]])
A, B 交換做內積結果都是單元矩陣, 可見 A 的反矩陣為 B.
反矩陣的用途之一是用來解線性聯立方程組 AX=b, 其中 A 為係數矩陣 (coefficient matrix), X 是變數向量, b 則是右手係數 (right-hand coefficients) 向量, 則方程組的解如下 :
例如下面這個三元線性聯立方程組 :
寫成矩陣方程式為 :
此方程組的解為係數矩陣之反矩陣與右手係數向量之內積 :
下面用 linalg 的 inv() 函數來計算其解 :
>>> A=np.array([[3, 2, 3], [2, 5, -7], [1, 2, -2]]) # 係數矩陣
>>> A
array([[ 3, 2, 3],
[ 2, 5, -7],
[ 1, 2, -2]])
>>> b=np.array([9, -12, -3]) # 右手係數向量
>>> b
array([ 9, -12, -3])
>>> B=linalg.inv(A) # 係數矩陣的反矩陣
>>> B
array([[ 1.33333333, 3.33333333, -9.66666667],
[-1. , -3. , 9. ],
[-0.33333333, -1.33333333, 3.66666667]])
>>> np.dot(B, b) # 聯立方程組的解
array([ 1.00000000e+00, -7.10542736e-15, 2.00000000e+00])
第二個變數之解近乎 0, 故其解為 [x, y, z]=[1, 0, 2].
注意, 行列式值為 0 的矩陣不可逆, 沒有反矩陣, 例如 :
>>> A=np.array([[2, 2, 2], [2, 2, 2], [2, 2, 2]])
>>> A
array([[2, 2, 2],
[2, 2, 2],
[2, 2, 2]])
>>> linalg.det(A) # 行列式值為 0 的矩陣不可逆
0.0
>>> linalg.inv(A)
Traceback (most recent call last):
File "<pyshell>", line 1, in <module>
File "C:\Python37\lib\site-packages\scipy\linalg\basic.py", line 979, in inv
raise LinAlgError("singular matrix")
numpy.linalg.LinAlgError: singular matrix
行列式值為 0 的矩陣為奇異矩陣, 沒有反矩陣.
3. 用 solve() 函數求解線性聯立方程組 :
上面我們使用反矩陣求解線性聯立方程組之解, 事實上呼叫 solve(A, b) 函數並將係數矩陣 A 與右手係數向量 b 傳入分別做為第一與第二參數亦可求解, 它會傳回以串列表示的變數向量, 例如 :
>>> A=np.array([[3, 2, 3], [2, 5, -7], [1, 2, -2]]) # 係數矩陣
>>> A
array([[ 3, 2, 3],
[ 2, 5, -7],
[ 1, 2, -2]])
>>> b=np.array([9, -12, -3]) # 右手係數向量
>>> b
array([ 9, -12, -3])
>>> x=linalg.solve(A, b) # 求解線性聯立方程組
>>> x
array([1., 0., 2.]) # 變數向量 [x, y, z] 之解
可見答案與上面用反矩陣計算的相同.
4. 用 eig() 函數求特徵值與特徵向量 :
一個 n 階方陣 A 的特徵值與特徵向量定義如下 : 存在一個非 0 向量 X 與純量 λ 使得 AX=λX, 則 λ 稱為方陣 A 的特徵值, 而向量 X 為其特徵向量. AX=λX 可看成向量 X 經過 A 的線性轉換後得到一個新向量 Y=AX, 而此新向量與原向量 X 呈平行關係, 特徵值即兩平行向量之比例係數.
AX=λX 經移項後為 (A-λI)X=0, 故 λ 為 A 的特徵值之充要條件為 |A-λI|=0, 將此式展開即可得到 λ 的特徵方程式, 解此方程式之根即可得到特徵值 λ. 呼叫 linalg 子套件的 eig(A) 函數會傳回兩個元素的 tuple, 第一個傳回值為特徵值 λ, 第二個傳回值為特徵向量 :
l, v=linalg.eig(A)
例如 :
>>> A=np.array([[1, 2],[3, 2]]) # 二階方陣
>>> A
array([[1, 2],
[3, 2]])
>>> l, v=linalg.eig(A) # 求特徵值與特徵向量
>>> l # 有兩個特徵值 : -1 與 4
array([-1.+0.j, 4.+0.j])
>>> v # 有兩個特徵向量 : 前者對應特徵值 -1, 後者對應 4
array([[-0.70710678, -0.5547002 ],
[ 0.70710678, -0.83205029]])
此例有兩組特徵向量, 分別對應兩個特徵值.
特徵方程式可能會有複數根, 例如 :
>>> A=np.array([[3, 2, 3], [2, 5, -7], [1, 2, -2]])
>>> A
array([[ 3, 2, 3],
[ 2, 5, -7],
[ 1, 2, -2]])
>>> l, v=linalg.eig(A)
>>> l # 有三個特徵值
array([4.90057187+0.j , 0.54971406+0.55676557j,
0.54971406-0.55676557j])
>>> v # 有三組特徵向量
array([[-0.86169137+0.j , -0.74295072+0.j ,
-0.74295072-0.j ],
[-0.44017844+0.j , 0.63299867-0.06720996j,
0.63299867+0.06720996j],
[-0.25244984+0.j , 0.18481478-0.09307649j,
0.18481478+0.09307649j]])
此例有三個特徵向量, 分別對應三個特徵值.
5. 用 qr() 函數對方陣進行 QR 分解 :
一個實數方陣 A 的 QR 分解是指它可拆解為一個正交矩陣 Q 與一個上三角矩陣 R 的內積, 所謂正交矩陣意指轉置矩陣等於反矩陣, 亦即正交矩陣與其轉置矩陣之內積為單元矩陣 :
參考 :
linalg 子套件的 qr(A) 函數會傳回方陣 A 的正交矩陣 Q 與一個上三角矩陣 R 組成之 tuple, 使得 A=QR, 例如 :
>>> A=np.array([[1, 2],[3, 2]]) # 二階方陣
>>> A
array([[1, 2],
[3, 2]])
>>> Q, R=linalg.qr(A) # 對方陣 A 進行 QR 分解
>>> Q # 正交矩陣 Q
array([[-0.31622777, -0.9486833 ],
[-0.9486833 , 0.31622777]])
>>> R # 上三角矩陣 R
array([[-3.16227766, -2.52982213],
[ 0. , -1.26491106]])
>>> np.dot(Q, Q.T) # 正交矩陣與其轉置矩陣之內積必為單元矩陣 I
array([[ 1.0000000e+00, -5.8339229e-18],
[-5.8339229e-18, 1.0000000e+00]])
>>> B=linalg.inv(Q) # 求正交矩陣 Q 的反矩陣 B
>>> B # 正交矩陣的反矩陣即為本身之轉置
array([[-0.31622777, -0.9486833 ],
[-0.9486833 , 0.31622777]])
>>> np.allclose(A, np.dot(Q, R)) # 檢查 A 與 QR 是否相等
True
由上面驗算可知, 經 QR 分解後得到的正交矩陣與其轉置矩陣之內積為單元矩陣 I. 此外也使用了 Numpy 的 allclose() 函數檢驗 A 確實等於 QR. 函數 np.allclose() 常用來檢驗兩個陣列 A, B 的每個對應元素是否都相等 (element-wise equal), 若全部相等傳回 True, 否則傳回 False :
np.allclose(A, B)
其預設判別條件如下 :
其中 a, b 分別是矩陣 A, B 中的相對應元素, atol 為絕對容忍參數, rtol 為相對容忍參數 (相對於 b 的絕對值). atol 預設值為 1e-8 (即 0.00000001), rtol 預設值為 1e-5 (即 0.00001). 當對應元素 a, b 差之絕對值小於等於右方條件值時傳回 True, 否則傳回 False, 參考 :
例如 :
>>> np.allclose([1e10,1e-7], [1.00001e10,1e-8])
False
>>> np.allclose([1e10,1e-8], [1.00001e10,1e-9])
True
下面是三階方陣的範例 :
>>> A=np.array([[3, 2, 3], [2, 5, -7], [1, 2, -2]]) # 三階方陣
>>> A
array([[ 3, 2, 3],
[ 2, 5, -7],
[ 1, 2, -2]])
>>> Q, R=linalg.qr(A) # 對方陣 A 進行 QR 分解
>>> Q # 正交矩陣 Q
array([[-0.80178373, 0.59152048, -0.08512565],
[-0.53452248, -0.77352678, -0.34050261],
[-0.26726124, -0.22750788, 0.93638218]])
>>> R # 上三角矩陣 R
array([[-3.74165739, -4.81070235, 1.87082869],
[ 0. , -3.13960871, 7.64426469],
[ 0. , 0. , 0.25537696]])
>>> np.dot(Q, Q.T) # 正交矩陣與其轉置矩陣之內積必為單元矩陣 I
array([[ 1.00000000e+00, -4.66827300e-17, 8.68933848e-17],
[-4.66827300e-17, 1.00000000e+00, 7.00177831e-17],
[ 8.68933848e-17, 7.00177831e-17, 1.00000000e+00]])
>>> B=linalg.inv(Q) # 求正交矩陣 Q 的反矩陣 B
>>> B # 正交矩陣的反矩陣即為本身之轉置
array([[-0.80178373, -0.53452248, -0.26726124],
[ 0.59152048, -0.77352678, -0.22750788],
[-0.08512565, -0.34050261, 0.93638218]])
>>> np.allclose(A, np.dot(Q, R)) # 檢查 A 與 QR 是否相等
True
從三階方陣更能清楚看出正交矩陣的特性.
6. 用 lu() 函數對方陣進行 LU 分解 :
LU 分解是將一個方陣 A 分解為一個下三角矩陣 L 與上三角矩陣 U 的內積 :
在線性代數中 LU 分解是高斯消去法的矩陣形式, 也是求反矩陣與計算行列式值的關鍵步驟, 主要用來解線性聯立方程組. 不過, 並非每個方陣都存在 LU 分解, 但若透過一個置換矩陣 (permutation matrix) P 先做行列順序調換, 則任何方陣都可以進行 LU 分解了 :
所謂置換矩陣是指每列與每行都只有一個 1, 其餘元素均為 0 的方陣, 參考 :
SciPy 的 linalg 所提供 lu() 函數的就是做 PLU 分解, 它會傳回 P, L, U 組成之元組 :
P, L, U=linalg.lu(A)
例如 :
>>> A=np.array([[3, 2, 3], [2, 5, -7], [1, 2, -2]]) # 三階方陣
>>> A
array([[ 3, 2, 3],
[ 2, 5, -7],
[ 1, 2, -2]])
>>> P, L, U=linalg.lu(A) # 對方陣 A 做 PLU 分解
>>> P # 置換矩陣 P
array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
>>> L # 下三角矩陣 L
array([[1. , 0. , 0. ],
[0.66666667, 1. , 0. ],
[0.33333333, 0.36363636, 1. ]])
>>> U # 上三角矩陣 U
array([[ 3. , 2. , 3. ],
[ 0. , 3.66666667, -9. ],
[ 0. , 0. , 0.27272727]])
>>> np.allclose(A, P.dot(L.dot(U))) # 檢查 A 是否等於 PLU 內積
True
下面範例使用 random 子套件來產生一個隨機三階方陣, 其 LU 分解如下 :
>>> from scipy import random, linalg # 匯入隨機子套件 random
>>> A=random.randn(3, 3) # 產生三階隨機方陣
>>> A
array([[ 1.05600878, 0.35537349, -0.04829833],
[ 1.33457398, 0.192862 , -0.03187792],
[ 0.57623532, -0.36047869, -0.31709815]])
>>> P, L, U=linalg.lu(A) # 對方陣 A 做 PLU 分解
>>> P # 置換矩陣 P
array([[0., 0., 1.],
[1., 0., 0.],
[0., 1., 0.]])
>>> L # 下三角矩陣 L
array([[ 1. , 0. , 0. ],
[ 0.43177473, 1. , 0. ],
[ 0.79127032, -0.4569392 , 1. ]])
>>> U # 上三角矩陣 U
array([[ 1.33457398, 0.192862 , -0.03187792],
[ 0. , -0.44375163, -0.30333407],
[ 0. , 0. , -0.16167951]])
>>> np.allclose(A, P.dot(L.dot(U)))
True
此例的 P 不是單元矩陣, 可見方陣 A 是經過行列順序調換才能做 LU 分解.
參考 :
