iT邦幫忙

2023 iThome 鐵人賽

DAY 4
0
自我挑戰組

從電子元件到傅立葉轉換 - 那些我有興趣的主題系列 第 4

[Day4] 逆矩陣沒有想像中的好

  • 分享至 

  • xImage
  •  

逆矩陣沒有想像中的好

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

線性方程式求解

picture

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

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

逆矩陣

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

  • 透過伴隨矩陣(Adj Matrix)與行列式來求解逆矩陣

picture2

  • 透過高斯喬登(Gauss Jordan)消去法求解逆矩陣

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

由於公式當中出現

picture3

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

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

picture4

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

picture5

此時的 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

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

picture7

用 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 * yU * x 的步驟,看起來比較麻煩,但還是優於逆矩陣?由於 LU 矩陣的特性,求算 LyUx 的時候,計算量約為 O(n^2) 左右。

結論

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

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

後記

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


上一篇
[Day3] 奇異值分解與應用
下一篇
[Day5] Seam Carving – 讓圖片不成比例地縮小
系列文
從電子元件到傅立葉轉換 - 那些我有興趣的主題30
圖片
  直播研討會
圖片
{{ item.channelVendor }} {{ item.webinarstarted }} |
{{ formatDate(item.duration) }}
直播中

尚未有邦友留言

立即登入留言