2021年10月21日 星期四

Python 學習筆記 : SciPy 的線性代數子套件 linalg 測試

SciPy 的線性代數子套件 linalg 是以 Numpy 的 ndarry 陣列為基礎建構的, 它提供了比 Numpy 更多更完整的線性代數函式庫, 事實上 SciPy 可視為 Numpy 的擴充集 (superset), 但 Numpy 基於歷史緣故仍保有核心的線性代數函式. 


參考書籍 :

基礎線性代數 (黃學亮, 五南, 2013)

使用 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]])
>>> linalg.det(A)    
0.0

矩陣的行列式有如下特性 :
  • 行列式值為 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 分解. 

參考 :


沒有留言:

張貼留言