逆矩陣沒有想像中的好

今天來談談一個大家耳熟能詳的名詞：逆矩陣。今天來談談為什麼逆矩陣沒有想像中的好，以及為什麼在實際計算當中我們應該避免使用逆矩陣。

線性方程式求解

Ax=b

前一篇有提到，我們可以直接求 A 的逆矩陣，然後再求算 $A^{-1}b$ ，即可求得 x。

聽起來輕鬆寫意，然而在實際應用當中直接去求算逆矩陣解線性方程式是比較少見的。實際應用上，矩陣裡的元素不會全部都是整數，通常會充滿各種浮點數，維數也大很多。不會像大學的題目那樣，最多就是 3x3 矩陣，而且數字都盡可能設計得很好，不太會出現小數點的計算。

逆矩陣

求逆矩陣主要有兩個經典手法：

透過伴隨矩陣（Adj Matrix）與行列式來求解逆矩陣

A^{-1}=adj(A) / det(A)

透過高斯喬登（Gauss Jordan）消去法求解逆矩陣

直觀上來看，用伴隨矩陣跟行列式來求解逆矩陣比用高斯喬登消去法麻煩多了。實際上，伴隨矩陣與行列式也的確不利於計算機計算。

由於公式當中出現

\frac{1}{det(A)}

這在行列式趨近於零的時候會讓值變得很大，進而導致誤差或甚至 overflow。另外就是伴隨矩陣需要把矩陣 A 拆分成更小的矩陣並計算其行列式，複雜度為 O(n!)。因此儘管公式看起來很簡單，但實際上不好計算。

比較常見的求逆矩陣則是用高斯喬登消去法來算，將增廣矩陣

\begin{bmatrix} A|I \end{bmatrix}

透過列運算化簡為上三角矩陣後再將非對角元素全部化簡為零，此時的擴展矩陣會是

\begin{bmatrix} I|B \end{bmatrix}

此時的 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，對角線以下的元素全為零的矩陣。

A = LU

計算出 LU 之後，我們可以再透過：

Ax = b \\ LUx = b \\ 令y= Ux，求解 Ly=b \\ y = Ux 解出 x

用 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(\frac{n^3}{3})

跟高斯喬登消去法的 O(n^3) 比起來差了 3 倍。

這邊比較不直觀的地方在於，明明 LU 分解多了 L * y 跟 U * x 的步驟，看起來比較麻煩，但還是優於逆矩陣？由於 LU 矩陣的特性，求算 Ly 與 Ux 的時候，計算量約為 O(n^2) 左右。

結論

從計算量來說，跟高斯喬登消去法求逆矩陣並沒有到非常巨大的差異，但在實際應用上，實際要求解的方程往往會是稀疏矩陣，也就是 0 佔據大部分的矩陣元素。

對於這類型的求解，如果直接求算逆矩陣，那麼這個逆矩陣不但非常醜，而且大部分的元素反而還不會是 0。然而如果改用 LU 分解，矩陣 A 的 0 元素大部分可以保留到 LU 上，讓計算變得比較簡單。

後記

有很多看起來很直觀的解法，不一定是最佳解。而找到最佳解的方式，往往可以從「換個角度看世界」來得到。例如之後會講到的快速傅立葉轉換、餘弦離散轉換，這些看起來有點繞的手法，卻是解決問題的妙計。

Kalan's Blog