はじめに
NTTドコモ サービスイノベーション部の中村圭佑です。普段の業務では画像認識や生成系AIに関する研究開発を行っています。ネット上の日本語記事で数値的な逆行列の計算方法を網羅した記事が少ないなぁと感じたので備忘録的に書いておこうと思います。今回はプログラムと簡単な説明を添えて紹介していきたいと思います。一般的に使われるロバストな手法とそうでないものを分けて示します。逆行列一つとっても様々な計算方法があるため、とても面白いですヨ。
何故、逆行列を計算する必要があるのか？？？
逆行列の計算は数学、工学、科学、経済学など多岐にわたる分野で重要な役割を果たしています。逆行列を理解し計算することは、連立方程式の解法、システムの安定性分析、最適化問題の解決など、多くの実用的応用に直結しています。以下に、逆行列の計算の使用例と重要性をわかりやすく説明します。
連立方程式の解法
逆行列は、連立一次方程式を解く際に重要な役割を果たします。例えば、行列形式で表された方程式 があるとき、
の逆行列が存在すれば、この方程式の解は
によって直接求めることができます。これにより、計算の手間を大幅に削減できます。特に複雑なシステムや大規模な問題における効率的な解法によく利用されます。
最適化問題
経済学や工学における最適化問題では、目的関数の最小化（または最大化）を求める際に逆行列が用いられます。特に、制約付き最適化問題を解く際にラグランジュ乗数法が用いられることがあり、この過程でヘッセ行列（二階偏微分の行列）の逆行列が必要となる場合があります。
数値的安定性と精度
数値的安定性と計算精度の観点から逆行列を計算する方法も重要です。実際の計算では、逆行列が存在するにもかかわらず、数値誤差によって正確な逆行列を得ることができない場合があります。そのため、LU分解・QR分解・特異値分解などの数値的に安定した方法を用いて逆行列を計算することが、実用的な問題解決において非常に重要になります。
紹介する手法はこちら
今回は大きく分けて6つの手法を紹介します。
数値的にロバストかつ一般的によく使われる手法
- LU分解
- QR分解
- 特異値分解 (SVD)
特殊な条件下や高度な数値解析で用いられる手法
- ガウスジョルダン法
- コレスキー分解
- シャーマン・モリソン公式
一般的によく使われる手法
LU分解
LU分解は行列を下三角行列と上三角行列
の積に分解する方法です。逆行列の計算にLU分解を利用するプロセスは以下の通りです。
1. 行列を
に分解します。
2. 単位行列の各列に対して、前進代入と後退代入を使って
と
を利用し、逆行列
の列を一つずつ求めます。
3. 得られた各列を組み合わせて、逆行列を構築します。
この方法は、分解が一度行われると、異なる右辺に対する線形方程式を効率的に解くことができる点で非常に有利です。
import numpy as np
def lu_decomposition_with_partial_pivoting(A):
n = A.shape[0]
L = np.zeros((n, n), dtype=float)
U = np.zeros((n, n), dtype=float)
P = np.eye(n, dtype=float) # 初期置換行列は単位行列
for i in range(n):
pivot = np.argmax(np.abs(A[i:, i])) + i
if pivot != i:
A[[i, pivot]] = A[[pivot, i]]
P[[i, pivot]] = P[[pivot, i]]
if i > 0:
L[[i, pivot], :i] = L[[pivot, i], :i]
L[i, i] = 1
for j in range(i, n):
U[i, j] = A[i, j] - np.dot(L[i, :i], U[:i, j])
for j in range(i+1, n):
L[j, i] = (A[j, i] - np.dot(L[j, :i], U[:i, i])) / U[i, i]
return P, L, U
def lu_inverse(P, L, U):
n = L.shape[0]
I = np.eye(n, dtype=float)
Y = np.zeros_like(L)
for i in range(n):
Y[i, :] = I[i, :] - np.dot(L[i, :i], Y[:i, :])
X = np.zeros_like(U)
for i in range(n-1, -1, -1):
X[i, :] = (Y[i, :] - np.dot(U[i, i+1:], X[i+1:, :])) / U[i, i]
return X
# 使用例
P, L, U = lu_decomposition_with_partial_pivoting(matrix)
matrix_inv = lu_inverse(P, L, U)
QR分解
QR分解は行列を直交行列と上三角行列
に分解する方法です。逆行列の計算には以下の手順が含まれます。
1. 行列を
に分解します。
2. 次に の関係を使って、
の逆行列を求め、それに
を乗じて逆行列
を計算します。
QR分解は数値的に非常に安定しており、特に条件数が大きな行列に対して有効です。
def qr_inverse(A):
"""
QR分解を用いて逆行列を計算する。
A: 逆行列を計算したい正方行列
"""
# QR分解
Q, R = np.linalg.qr(A)
# 単位行列（Aの逆行列を乗算すると得られる結果）
I = np.eye(A.shape[0])
# QR分解の結果を用いて逆行列を計算
# Rx = Q^T * b (ここでbは単位行列の列ベクトル)
# 逆行列の各列を求める
A_inv = np.zeros_like(A)
for i in range(A.shape[0]):
b = I[:, i]
# 後退代入を用いて上三角行列Rの逆行列を計算
x = np.linalg.solve(R, np.dot(Q.T, b))
A_inv[:, i] = x
return A_inv
# QR分解を用いて逆行列を計算
matrix_inv = qr_inverse(matrix)
特異値分解 (SVD)
特異値分解は行列をの形に分解します。逆行列の計算には以下の手順が含まれます。
1. 行列のSVDを計算します。
2. の非ゼロ要素の逆数を取り、逆行列
を構築します。
3. の式を用いて逆行列を計算します。
SVDは逆行列が存在しない場合や悪条件な行列に対しても擬似逆行列を求めるのに有効です。
def svd_inverse(A):
"""
特異値分解(SVD)を用いて逆行列を計算する。
A: 逆行列を計算したい正方行列
"""
# SVDを実行
U, s, VT = np.linalg.svd(A)
# sは特異値の配列。逆行列の対角成分を計算するために逆数を取る
S_inv = np.diag(1 / s)
# 完全な逆行列を計算
A_inv = VT.T @ S_inv @ U.T
return A_inv
# SVDを用いて逆行列を計算
matrix_inv = svd_inverse(matrix)
特殊な条件下や高度な数値解析で用いられることが多い手法
ガウス・ジョルダン法
ガウス-ジョルダン法は、拡大行列を使って行基本変形を行い、元の行列を単位行列に変換することで逆行列を直接求める方法です。
1. 行列と単位行列
を横に並べた拡大行列を作成します。
2. ガウス-ジョルダン消去法を用いて左側のを単位行列に変換します。
3. 同じ行基本変形を右側のにも適用すると、変換後の右側が逆行列
になります。
この手法は直接的で理解しやすいですが、大規模な行列に対しては計算コストが高くなる可能性があります。
def gauss_jordan_inverse(matrix):
n = len(matrix)
# 入力行列が正方行列であることを確認
if matrix.shape[0] != matrix.shape[1]:
raise ValueError("Matrix must be square")
# 拡大行列を作成（元の行列と同じサイズの単位行列を横に結合）
A = np.hstack((matrix, np.eye(n)))
for i in range(n):
# ピボットが0でないことを確認（必要であれば行を交換）
if np.fabs(A[i, i]) < 1e-10:
for j in range(i+1, n):
if np.fabs(A[j, i]) > np.fabs(A[i, i]):
A[[i, j]] = A[[j, i]] # 行を交換
break
else:
raise ValueError("Matrix is singular")
# ピボット行をピボット値で割って正規化
A[i] = A[i] / A[i, i]
# 他のすべての行からピボット行の定数倍を引いて、i列目を0にする
for j in range(n):
if i != j:
A[j] -= A[i] * A[j, i]
# 拡大行列の右半分が逆行列
return A[:, n:]
# 逆行列を計算
inverse_matrix = gauss_jordan_inverse(matrix)
コレスキー分解
コレスキー分解は、対称かつ正定値な行列にのみ適用可能な手法です。この条件を満たさない行列に対しては、他の分解方法や計算方法を検討する必要があります。
この手法は対称正定値行列に対して非常に効果的なアプローチです。数値的安定性と計算効率の良さから、工学や科学の分野で頻繁に利用されています。
def cholesky_decomposition(A):
n = A.shape[0]
L = np.zeros_like(A)
for i in range(n):
for j in range(i+1):
sum_k = sum(L[i][k] * L[j][k] for k in range(j))
if i == j: # 対角要素
L[i][j] = np.sqrt(A[i][i] - sum_k)
else:
L[i][j] = (1.0 / L[j][j] * (A[i][j] - sum_k))
return L
def forward_substitution(L, b):
n = L.shape[0]
y = np.zeros_like(b, dtype=np.float)
for i in range(n):
y[i] = (b[i] - np.dot(L[i, :i], y[:i])) / L[i, i]
return y
def backward_substitution(U, y):
n = U.shape[0]
x = np.zeros_like(y, dtype=np.float)
for i in reversed(range(n)):
x[i] = (y[i] - np.dot(U[i, i+1:], x[i+1:])) / U[i, i]
return x
def cholesky_inverse(A):
n = A.shape[0]
L = cholesky_decomposition(A)
L_inv = np.zeros_like(A)
I = np.eye(n)
for i in range(n):
e = I[:, i]
y = forward_substitution(L, e)
L_inv[:, i] = backward_substitution(L.T, y)
return L_inv
matrix = np.array([[4, 1, 1], [1, 3, 1], [1, 1, 2]], dtype=float)
matrix_inv = cholesky_inverse(matrix)
シャーマン・モリソン公式
シャーマン・モリソン公式は、ある行列の逆行列を既に知っている状況で、その行列にランク1の更新を行った新しい行列の逆行列を効率的に計算するための公式です。ランク1の更新とは、元の行列にベクトルの外積（ベクトルとその転置ベクトルの積）の形で加算することを指します。この公式は、大規模な計算において逆行列の再計算を避けることができるため、計算コストを削減するのに役立ちます。
def sherman_morrison_formula(A_inv, u, v):
"""
シャーマン・モリソン公式を用いて逆行列を更新する。
A_inv: 更新前の逆行列
u, v: ランク1更新に用いるベクトル
"""
numerator = np.dot(A_inv, np.dot(u, np.dot(v, A_inv)))
denominator = 1.0 + np.dot(v, np.dot(A_inv, u))
A_inv_updated = A_inv - numerator / denominator
return A_inv_updated
# テスト用の行列とランク1更新ベクトル
A = np.array([[2, 1], [1, 2]], dtype=float)
u = np.array([1, 1], dtype=float)
v = np.array([1, -1], dtype=float)
# Aの逆行列を計算
A_inv = np.linalg.inv(A)
# シャーマン・モリソン公式を用いて逆行列を更新
A_inv_updated = sherman_morrison_formula(A_inv, u, v)
まとめ
今回は逆行列の数値的計算手法をまとめました。ふとした時、逆行列の計算方法を思い出したい方/気になってしまった方にこの記事が届けば嬉しいです。実際に使う場合は計算量やロバスト性を考慮しましょう。
またNTTドコモでは画像生成・物体検出・姿勢推定・一般物体認識・特定物体認識・文字認識・類似画像検索等、様々なAI技術の研究開発をしています。インターンや新卒・キャリア採用を積極的に実施しているため、気になる方は以下のリンクをご参考にしてください。
- NTTドコモ新卒採用
- インターンシップ情報 | NTTドコモ新卒採用
- NTTドコモキャリア採用
