逆矩陣沒有想像中的好

💡

質問やフィードバックがありましたら、フォームからお願いします

線形方程式の解法

Ax=b

前回の記事で述べたように、私たちは直接行列Aの逆行列を求め、その後 $A^{-1}b$ を計算することでxを求めることができます。

聞こえは簡単そうですが、実際の応用では逆行列を用いて線形方程式を解くことはあまり一般的ではありません。実際のデータでは、行列の要素がすべて整数であることは稀で、通常は様々な浮動小数点数が含まれており、次元も大きくなることが多いです。大学の問題のように、最大で3x3の行列で、数字もできるだけ整えられていることはほとんどありません。

逆行列

逆行列を求める主な二つの古典的手法があります：

伴随行列（Adj Matrix）と行列式を用いて逆行列を求める方法

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

ガウス・ジョルダン（Gauss Jordan）消去法を用いて逆行列を求める方法

直感的には、伴随行列と行列式を使って逆行列を求めるのは、ガウス・ジョルダン消去法を使うよりも手間がかかります。実際、伴随行列と行列式は計算機による計算には不向きです。

式に現れる

\frac{1}{det(A)}

は、行列式がゼロに近づくと値が非常に大きくなり、誤差やオーバーフローを引き起こす可能性があります。また、伴随行列は行列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):
        # 対角要素がゼロかどうかを確認
        if augmented[i][i] == 0:
            # 下の行と交換
            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 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で、対角線以下の要素がゼロでない行列、Uは対角線以上の要素がゼロでない行列です。

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  # Lの対角要素は1

        for j in range(i, n):
            # 上三角行列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):
            # 下三角行列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)

    # 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]

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

となり、ガウス・ジョルダン消去法のO(n^3)と比べて約3倍の差があります。

ここで直感的でない点は、LU分解にL * yとU * xのステップが追加され、手間がかかるように見えるのに、逆行列よりも優れているということです。LU行列の特性により、LyとUxを計算する際の計算量は約O(n^2)です。

結論

計算量の観点から見ると、ガウス・ジョルダン消去法で逆行列を求めることとの間に大きな違いはありませんが、実際の応用では、解くべき方程式はしばしば疎行列であり、つまり0が大部分を占める行列要素です。

このような解法の場合、逆行列を直接求めてしまうと、その逆行列は非常に複雑で、要素のほとんどがゼロでないことが多いです。しかし、LU分解を用いれば、行列Aのゼロ要素をLUに保持でき、計算がより簡単になります。

後記

直感的に見える多くの解法が、必ずしも最適解であるとは限りません。そして、最適解を見つける方法は、「別の視点から世界を見る」ことによって得られることが多いです。例えば、次回取り上げる高速フーリエ変換やコサイン離散変換など、少し複雑に見える手法が、問題解決の妙策であることがよくあります。

← 與程式有關的遊戲三選 (1) – A=B JPG 與離散餘弦轉換 →

この記事が役に立ったと思ったら、下のリンクからコーヒーを奢ってくれると嬉しいです ☕ 私の普通の一日が輝かしいものになります ✨

☕Buy me a coffee

逆矩陣沒有想像中的好

目次

線形方程式の解法

逆行列

LU分解

結論

後記

目次