今天來談談一個大家耳熟能詳的名詞:逆矩陣。今天來談談為什麼逆矩陣沒有想像中的好,以及為什麼在實際計算當中我們應該避免使用逆矩陣。
前一篇有提到,我們可以直接求 A 的逆矩陣,然後再求算 $A^{-1}b$,即可求得 x
。
聽起來輕鬆寫意,然而在實際應用當中直接去求算逆矩陣解線性方程式是比較少見的。實際應用上,矩陣裡的元素不會全部都是整數,通常會充滿各種浮點數,維數也大很多。不會像大學的題目那樣,最多就是 3x3 矩陣,而且數字都盡可能設計得很好,不太會出現小數點的計算。
求逆矩陣主要有兩個經典手法:
直觀上來看,用伴隨矩陣跟行列式來求解逆矩陣比用高斯喬登消去法麻煩多了。實際上,伴隨矩陣與行列式也的確不利於計算機計算。
由於公式當中出現
這在行列式趨近於零的時候會讓值變得很大,進而導致誤差或甚至 overflow。另外就是伴隨矩陣需要把矩陣 A 拆分成更小的矩陣並計算其行列式,複雜度為 O(n!)。因此儘管公式看起來很簡單,但實際上不好計算。
比較常見的求逆矩陣則是用高斯喬登消去法來算,將增廣矩陣
透過列運算化簡為上三角矩陣後再將非對角元素全部化簡為零,此時的擴展矩陣會是
此時的 B 就是 A 的逆矩陣。如果用 Python 來寫會像這樣(假設矩陣有逆矩陣的情況下):
import numpy as np
def gauss_jordan(m):
n = len(m)
identity = np.eye(n)
# 增廣矩陣
augmented = np.hstack((m, identity))
for i in range(n):
# checking if diagonal element is zero
if augmented[i][i] == 0:
# swapping with row below
for j in range(i+1, n):
if augmented[j][i] != 0:
augmented[[i, j]] = augmented[[j, i]]
break
augmented[i] = augmented[i] / augmented[i, i]
for j in range(0, n):
if i != j:
augmented[j] -= augmented[j, i] * augmented[i]
# return the right side of the augmented matrix
return augmented[:, n:]
m = np.array([[1., 2.],
[3., 4.]])
inv_m = gauss_jordan(m)
numpy 裡面有提供 np.linalg.inv
函數,這邊為了示範所以自己寫了簡單版本。如果從時間複雜度來看,可以發現為 O(n^3)
。
在計算機中求逆矩陣來解線性方程式其實並不是很常見的事,只是在大學時期手算為主的解題思維下,反而會去懷疑明明直接用逆矩陣就好,為什麼還要大費周章分解矩陣。
LU 分解是將矩陣 A 拆分為兩個矩陣 L 與 U。分別為下三角矩陣與上三角矩陣。其中,L 為對角元素為 1;對角線以下的元素不為零,且對角線以上的元素全為 0 的矩陣;U 為對角線以上的元素不為 0,對角線以下的元素全為零的矩陣。
計算出 LU 之後,我們可以再透過:
用 numpy 實作的話大概長這樣:
import numpy as np
def lu_decomposition(A):
n = len(A)
L = np.zeros((n, n))
U = np.zeros((n, n))
for i in range(n):
L[i][i] = 1 # Diagonal elements of L are 1
for j in range(i, n):
# Compute the upper triangular matrix U
summation = sum(L[i][k] * U[k][j] for k in range(i))
U[i][j] = A[i][j] - summation
for j in range(i + 1, n):
# Compute the lower triangular matrix L
summation = sum(L[j][k] * U[k][i] for k in range(i))
L[j][i] = (A[j][i] - summation) / U[i][i]
return L, U
def lu_solve(L, U, b):
n = len(L)
y = np.zeros(n)
x = np.zeros(n)
# Solving L * y = b
for i in range(n):
summation = sum(L[i][j] * y[j] for j in range(i))
y[i] = (b[i] - summation) / L[i][i]
# Solving U * x = y
for i in range(n - 1, -1, -1):
summation = sum(U[i][j] * x[j] for j in range(i + 1, n))
x[i] = (y[i] - summation) / U[i][i]
return x
A = np.array([[1, 1, 1], [2, 3, 1], [7,4,2]])
b = np.array([3, 6, 13])
L, U = lu_decomposition(A)
x = lu_solve(L, U, b)
LU 分解對比高斯喬登消去法,計算量約為:O(n^3 / 3)
,跟高斯喬登消去法的 O(n^3)
比起來差了 3 倍。
這邊比較不直觀的地方在於,明明 LU 分解多了 L * y
跟 U * x
的步驟,看起來比較麻煩,但還是優於逆矩陣?由於 LU 矩陣的特性,求算 Ly
與 Ux
的時候,計算量約為 O(n^2)
左右。
從計算量來說,跟高斯喬登消去法求逆矩陣並沒有到非常巨大的差異,但在實際應用上,實際要求解的方程往往會是稀疏矩陣,也就是 0
佔據大部分的矩陣元素。
對於這類型的求解,如果直接求算逆矩陣,那麼這個逆矩陣不但非常醜,而且大部分的元素反而還不會是 0。然而如果改用 LU 分解,矩陣 A 的 0 元素大部分可以保留到 LU 上,讓計算變得更簡單。
有很多看起來很直觀的解法,不一定是最佳解。而找到最佳解的方式,往往可以從「換個角度看世界」來得到。例如之後會講到的快速傅立葉轉換、餘弦離散轉換,這些看起來有點繞的手法,卻是解決問題的妙計。