今天來談談一個大家耳熟能詳的名詞：逆矩陣。今天來談談為什麼逆矩陣沒有想像中的好，以及為什麼在實際計算當中我們應該避免使用逆矩陣。
線性方程式求解
前一篇有提到，我們可以直接求 A 的逆矩陣，然後再求算 ，即可求得 x。
聽起來輕鬆寫意，然而在實際應用當中直接去求算逆矩陣解線性方程式是比較少見的。實際應用上，矩陣裡的元素不會全部都是整數，通常會充滿各種浮點數，維數也大很多。不會像大學的題目那樣，最多就是 3x3 矩陣，而且數字都盡可能設計得很好，不太會出現小數點的計算。
逆矩陣
求逆矩陣主要有兩個經典手法：
- 透過伴隨矩陣（Adj Matrix）與行列式來求解逆矩陣
- 透過高斯喬登（Gauss Jordan）消去法求解逆矩陣
直觀上來看，用伴隨矩陣跟行列式來求解逆矩陣比用高斯喬登消去法麻煩多了。實際上，伴隨矩陣與行列式也的確不利於計算機計算。
由於公式當中出現
這在行列式趨近於零的時候會讓值變得很大，進而導致誤差或甚至 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 分解
在計算機中求逆矩陣來解線性方程式其實並不是很常見的事，只是在大學時期手算為主的解題思維下，反而會去懷疑明明直接用逆矩陣就好，為什麼還要大費周章分解矩陣。
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 倍。
這邊比較不直觀的地方在於，明明 LU 分解多了 L * y 跟 U * x 的步驟，看起來比較麻煩，但還是優於逆矩陣？由於 LU 矩陣的特性，求算 Ly 與 Ux 的時候，計算量約為 O(n^2) 左右。
結論
從計算量來說，跟高斯喬登消去法求逆矩陣並沒有到非常巨大的差異，但在實際應用上，實際要求解的方程往往會是稀疏矩陣，也就是 0 佔據大部分的矩陣元素。
對於這類型的求解，如果直接求算逆矩陣，那麼這個逆矩陣不但非常醜，而且大部分的元素反而還不會是 0。然而如果改用 LU 分解，矩陣 A 的 0 元素大部分可以保留到 LU 上，讓計算變得比較簡單。
後記
有很多看起來很直觀的解法，不一定是最佳解。而找到最佳解的方式，往往可以從「換個角度看世界」來得到。例如之後會講到的快速傅立葉轉換、餘弦離散轉換，這些看起來有點繞的手法，卻是解決問題的妙計。
