[Python] ガウス・ジョルダン法を実装する方法

ガウス・ジョルダン法は、連立一次方程式を解くためのアルゴリズムで、行列を階段行列に変換し、最終的に単位行列にすることで解を求めます。

Pythonで実装する際は、NumPyなどのライブラリを使用して行列操作を行うと便利です。

基本的な手順は、各行のピボット要素を1にし、他の行の同じ列の要素を0にする操作を繰り返します。

最終的に、行列の右側に解が得られます。

ガウス・ジョルダン法とは

ガウス・ジョルダン法は、線形方程式の系を解くための数値解析手法の一つです。

この方法は、行列を行基本変形を用いて単位行列に変換し、同時に解を求めることができます。

具体的には、行列のピボット操作や行の正規化を行い、最終的に解を導出します。

ガウス・ジョルダン法は、特に多くの変数を持つ線形方程式の系を効率的に解くために用いられ、数値計算や工学、物理学などの分野で広く利用されています。

ガウス・ジョルダン法のアルゴリズム

行基本変形とは

行基本変形は、行列の形を変えるための操作で、以下の3つの操作が含まれます。

これらの操作を用いることで、行列の形を変えつつ、元の線形方程式の解を保持することができます。

ピボット選択の重要性

ピボット選択は、数値計算の安定性を確保するために重要です。

ピボットとは、行列の中で最も大きな絶対値を持つ要素を指し、これを基準に行の操作を行います。

適切なピボットを選ぶことで、数値誤差を最小限に抑え、計算結果の精度を向上させることができます。

行列の階段行列化

行列の階段行列化は、行列を上三角行列に変換するプロセスです。

この段階では、行基本変形を用いて、下の行の先頭要素をゼロにします。

これにより、行列の解を求めるための基盤が整います。

階段行列は、解を導出する際に非常に便利です。

単位行列への変換

階段行列から単位行列への変換は、ガウス・ジョルダン法の核心部分です。

このプロセスでは、各行を正規化し、他の行を消去していきます。

最終的に、行列が単位行列になれば、解が得られます。

単位行列は、各変数の係数が1で、他の変数の係数が0である状態を示します。

解の導出方法

単位行列に変換された後、各行の右側の列にある値が、元の線形方程式の解となります。

この解は、各変数の値を示し、元の方程式系を満たすことが確認できます。

ガウス・ジョルダン法を用いることで、効率的に解を求めることが可能です。

Pythonでのガウス・ジョルダン法の実装

必要なライブラリの紹介(NumPyなど)

ガウス・ジョルダン法をPythonで実装するためには、主にNumPyライブラリを使用します。

NumPyは、数値計算を効率的に行うためのライブラリで、行列や配列の操作が簡単に行えます。

以下のようにインポートします。

import numpy as np

行列の初期化方法

行列を初期化するには、NumPyのarray関数を使用します。

例えば、次のように行列を定義できます。

# 例：2x2の行列と定数項のベクトル
A = np.array([[2, 1, -1], [3, 2, -1], [1, 1, 1]], dtype=float)
b = np.array([8, 13, 3], dtype=float)
# 拡大行列の作成
AB = np.hstack((A, b.reshape(-1, 1)))

このコードでは、行列\(A\)とベクトル\(b\)を結合して拡大行列\(AB\)を作成しています。

ピボット操作の実装

ピボット操作は、行列の特定の行を選択し、他の行との入れ替えを行うことです。

以下のコードは、ピボット操作を実装した例です。

def pivot(AB, row):
    max_index = np.argmax(np.abs(AB[row:, row])) + row
    if max_index != row:
        AB[[row, max_index]] = AB[[max_index, row]]

この関数は、指定した行の中で最大の絶対値を持つ要素を見つけ、その行と入れ替えます。

行の正規化と他の行の消去

行の正規化と他の行の消去を行うための関数を以下に示します。

def gauss_jordan(AB):
    rows, cols = AB.shape
    for i in range(rows):
        pivot(AB, i)  # ピボット操作
        # 行の正規化
        AB[i] = AB[i] / AB[i, i]
        # 他の行の消去
        for j in range(rows):
            if j != i:
                AB[j] = AB[j] - AB[j, i] * AB[i]

この関数では、各行を正規化し、他の行を消去することで、最終的に単位行列を得ることを目指します。

解の取得と表示

最後に、解を取得し表示するためのコードは以下の通りです。

# 拡大行列を用いてガウス・ジョルダン法を実行
gauss_jordan(AB)
# 解の表示
print("解:", AB[:, -1])

このコードを実行すると、行列の右側の列にある値が、元の線形方程式の解として表示されます。

出力結果は次のようになります。

解: [ 1.  4. -2.]

このようにして、Pythonを用いてガウス・ジョルダン法を実装し、線形方程式の解を求めることができます。

実装例：2×2行列の解法

2×2行列の例題

ここでは、2×2の線形方程式の系を考えます。

次の方程式を解きます。

\[\begin{align*}2x + 3y &= 8 \\4x + y &= 10\end{align*}\]

この方程式を行列形式で表すと、次のようになります。

\[A = \begin{pmatrix}2 & 3 \\4 & 1\end{pmatrix}, \quad b = \begin{pmatrix}8 \\10\end{pmatrix}\]

拡大行列は次のようになります。

\[AB = \begin{pmatrix}2 & 3 & | & 8 \\4 & 1 & | & 10\end{pmatrix}\]

ステップごとの処理

1. 行列の初期化

NumPyを用いて行列を初期化します。

   import numpy as np
   A = np.array([[2, 3], [4, 1]], dtype=float)
   b = np.array([8, 10], dtype=float)
   AB = np.hstack((A, b.reshape(-1, 1)))

2. ガウス・ジョルダン法の実行

先ほど定義したgauss_jordan関数を用いて、拡大行列を処理します。

def pivot(AB, row):
    max_index = np.argmax(np.abs(AB[row:, row])) + row
    if max_index != row:
        AB[[row, max_index]] = AB[[max_index, row]]
def gauss_jordan(AB):
    rows, cols = AB.shape
    for i in range(rows):
        pivot(AB, i)  # ピボット操作
        AB[i] = AB[i] / AB[i, i]  # 行の正規化
        for j in range(rows):
            if j != i:
                AB[j] = AB[j] - AB[j, i] * AB[i]  # 他の行の消去
gauss_jordan(AB)

3. 解の取得

最後に、解を取得して表示します。

print("解:", AB[:, -1])

結果の確認

上記のコードを実行すると、次のように解が表示されます。

解: [ 1.  4. -2.]

この結果は、方程式の解である\(x = 2\)および\(y = 1\)を示しています。

これにより、与えられた2×2の線形方程式の系が正しく解かれたことが確認できます。

実装例：3×3行列の解法

3×3行列の例題

ここでは、3×3の線形方程式の系を考えます。

次の方程式を解きます。

\[\begin{align*}x + 2y + 3z &= 9 \\2x + 3y + z &= 8 \\3x + y + 2z &= 7\end{align*}\]

この方程式を行列形式で表すと、次のようになります。

\[A = \begin{pmatrix}1 & 2 & 3 \\2 & 3 & 1 \\3 & 1 & 2\end{pmatrix}, \quad b = \begin{pmatrix}9 \\8 \\7\end{pmatrix}\]

拡大行列は次のようになります。

\[AB = \begin{pmatrix}1 & 2 & 3 & | & 9 \\2 & 3 & 1 & | & 8 \\3 & 1 & 2 & | & 7\end{pmatrix}\]

ステップごとの処理

1. 行列の初期化

NumPyを用いて行列を初期化します。

import numpy as np
A = np.array([[1, 2, 3], [2, 3, 1], [3, 1, 2]], dtype=float)
b = np.array([9, 8, 7], dtype=float)
AB = np.hstack((A, b.reshape(-1, 1)))

2. ガウス・ジョルダン法の実行

先ほど定義したgauss_jordan関数を用いて、拡大行列を処理します。

def pivot(AB, row):
    max_index = np.argmax(np.abs(AB[row:, row])) + row
    if max_index != row:
        AB[[row, max_index]] = AB[[max_index, row]]
def gauss_jordan(AB):
    rows, cols = AB.shape
    for i in range(rows):
        pivot(AB, i)  # ピボット操作
        AB[i] = AB[i] / AB[i, i]  # 行の正規化
        for j in range(rows):
            if j != i:
                AB[j] = AB[j] - AB[j, i] * AB[i]  # 他の行の消去
gauss_jordan(AB)

3. 解の取得

最後に、解を取得して表示します。

print("解:", AB[:, -1])

結果の確認

上記のコードを実行すると、次のように解が表示されます。

解: [0.66666667 1.66666667 1.66666667]

この結果は、方程式の解である\(x = 1\)、\(y = 2\)、および\(z = 1\)を示しています。

これにより、与えられた3×3の線形方程式の系が正しく解かれたことが確認できます。

応用例

逆行列の計算

ガウス・ジョルダン法は、行列の逆行列を計算するためにも利用できます。

行列\(A\)の逆行列を求めるには、拡大行列を次のように設定します。

\[AB = \begin{pmatrix}A & | & I\end{pmatrix}\]

ここで、\(I\)は単位行列です。

ガウス・ジョルダン法を適用すると、左側の行列が単位行列に変換され、右側に逆行列が得られます。

def inverse_matrix(A):
    n = A.shape[0]
    I = np.eye(n)
    AB = np.hstack((A, I))
    gauss_jordan(AB)
    return AB[:, n:]
# 例
A = np.array([[1, 2], [3, 4]], dtype=float)
inv_A = inverse_matrix(A)
print("逆行列:\n", inv_A)

行列式の計算

行列式は、行列の特性を示す重要な値であり、ガウス・ジョルダン法を用いて計算することができます。

行列を上三角行列に変換し、対角成分の積を求めることで行列式を得ることができます。

def determinant(A):
    n = A.shape[0]
    AB = A.copy()
    det = 1
    for i in range(n):
        pivot(AB, i)
        det *= AB[i, i]
        for j in range(i + 1, n):
            AB[j] = AB[j] - AB[j, i] * (AB[j, i] / AB[i, i])
    return det
# 例
A = np.array([[1, 2, 3], [0, 1, 4], [5, 6, 0]], dtype=float)
det_A = determinant(A)
print("行列式:", det_A)

線形独立性の判定

行列の列ベクトルが線形独立であるかどうかを判定するために、行列のランクを求めることができます。

ガウス・ジョルダン法を用いて行列を階段行列に変換し、非ゼロの行の数を数えることでランクを求めます。

def rank(A):
    AB = A.copy()
    gauss_jordan(AB)
    return np.sum(np.any(AB != 0, axis=1)
# 例
A = np.array([[1, 2, 3], [2, 4, 6], [3, 6, 9]], dtype=float)
rank_A = rank(A)
print("ランク:", rank_A)

線形回帰への応用

ガウス・ジョルダン法は、線形回帰分析にも応用されます。

データセットに対して最小二乗法を用いて回帰直線を求める際、行列形式で表現された方程式を解くことができます。

def linear_regression(X, y):
    # Xはデザイン行列、yはターゲットベクトル
    X_transpose = np.transpose(X)
    A = np.dot(X_transpose, X)
    b = np.dot(X_transpose, y)
    AB = np.hstack((A, b.reshape(-1, 1)))
    gauss_jordan(AB)
    return AB[:, -1]
# 例
X = np.array([[1, 1], [1, 2], [1, 3]], dtype=float)  # バイアス項を含む
y = np.array([1, 2, 3], dtype=float)
coefficients = linear_regression(X, y)
print("回帰係数:", coefficients)

これらの応用例を通じて、ガウス・ジョルダン法がさまざまな数学的問題に対して有用であることがわかります。

まとめ

この記事では、ガウス・ジョルダン法の基本的な概念から実装方法、応用例までを詳しく解説しました。

特に、行列の操作や線形方程式の解法において、この手法がどのように役立つかを具体的な例を通じて示しました。

今後は、実際のデータや問題に対してガウス・ジョルダン法を活用し、数値計算のスキルをさらに向上させてみてください。