逆行列は想像ほど良くない

今日は、誰もが耳にしたことのある用語「逆行列」について話そうと思う。なぜ逆行列が想像ほど良くないのか、そして実際の計算においてなぜ逆行列の使用を避けるべきなのかについて解説する。

連立一次方程式の求解

前回の記事で触れたように、A の逆行列を直接求めてから $A^{-1}b$ を計算すれば、x を求めることができる。

一見すると簡単そうに思えるが、実際のアプリケーションにおいて、連立一次方程式を解くために逆行列を直接計算することは極めて稀だ。実践的な場面では、行列の要素がすべて整数であることはまずなく、通常はさまざまな浮動小数点数で溢れており、次元もはるかに大きい。大学の試験問題のように、せいぜい 3x3 の行列で、しかも小数の計算がほとんど出ないように数字が綺麗に設計されている、などということはないのだ。

逆行列

逆行列を求める古典的な手法は、主に 2 つある：

• 余因子行列（Adj Matrix）と行列式を使って逆行列を解く

• ガウス・ジョルダン（Gauss-Jordan）の消去法を用いて逆行列を解く

直感的に見ても、余因子行列と行列式を使って逆行列を解くのは、ガウス・ジョルダンの消去法を使うよりもはるかに面倒だ。実際、余因子行列と行列式はコンピュータによる計算には全く適していない。

公式の中に

が現れるため、行列式が 0 に近づくと値が非常に大きくなり、誤差やオーバーフロー（overflow）を引き起こす原因になる。さらに、余因子行列を求めるには行列 A をより小さな小行列に分解してそれぞれの行列式を計算する必要があり、その計算量は O(n!) に達する。そのため、公式自体はシンプルに見えても、実際の計算には向かない。

より一般的な逆行列の求め方は、ガウス・ジョルダンの消去法を用いるものだ。拡大係数行列

を行基本変形によって上三角行列に変形した後、対角以外の成分をすべて 0 に変形していく。このときの拡大行列は次のようになる。

このときの 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 を 2 つの行列 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 である行列になる。

このような問題を解く際、直接逆行列を求めてしまうと、その逆行列は数値的に非常に扱いづらいものになるだけでなく、要素の大部分が非ゼロになってしまう。しかし LU 分解を用いれば、行列 A の 0 の要素の多くを L と U にそのまま保持できるため、計算を大幅に単純化できるのだ。

あとがき

一見すると非常に直感的な解法が、必ずしも最適な解であるとは限らない。そして最適な解を見つける方法は、往々にして「視点を変えて世界を見る」ことから得られる。例えば今後取り上げる予定の高速フーリエ変換（FFT）や離散コサイン変換（DCT）なども、一見すると回りくどい手法に見えるが、実は問題を解決するための極めて優れた妙手なのである。