Kalan's Blog

本部落主要是關於前端、軟體開發以及我在日本的生活,也會寫一些對時事的觀察和雜感

目前主題 亮色

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

線性方程式求解

Ax=bAx=b

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

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

逆矩陣

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

  • 透過伴隨矩陣(Adj Matrix)與行列式來求解逆矩陣
A1=adj(A)/det(A)A^{-1}=adj(A) / det(A)
  • 透過高斯喬登(Gauss Jordan)消去法求解逆矩陣

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

由於公式當中出現

1det(A)\frac{1}{det(A)}

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

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

[AI]\begin{bmatrix} A|I \end{bmatrix}

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

[IB]\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,對角線以下的元素全為零的矩陣。

picture6

A=LUA = LU

計算出 LU 之後,我們可以再透過:

Ax=bLUx=by=Ux,求解Ly=by=Ux解出xAx = 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(n33)O(\frac{n^3}{3})

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

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

結論

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

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

後記

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

上一篇

與程式有關的遊戲三選 (1) – A=B

下一篇

JPG 與離散餘弦轉換

如果覺得這篇文章對你有幫助的話,可以考慮到下面的連結請我喝一杯 ☕️ 可以讓我平凡的一天變得閃閃發光 ✨

Buy me a coffee

作者

Kalan 頭像照片,在淡水拍攝,淺藍背景

愷開 | Kalan

愷開。台灣人,在 2019 年到日本工作,目前定居在福岡。除了熟悉前端之外對 IoT、App 開發、後端、電子電路領域都有涉略。最近開始玩電吉他。 歡迎 Email 諮詢或合作,聊聊音樂也可以,希望能透過這個部落格和更多的人交流。