Pythonで掃き出し法を用いた逆行列の計算方法と活用7選

import numpy as np

def gauss_jordan(A, b):
    # 行列Aと定数項bを結合
    Ab = np.column_stack((A, b))
    n = len(Ab)

    for i in range(n):
        # ピボット選択
        max_element = abs(Ab[i:, i]).argmax() + i
        if Ab[max_element, i] == 0:
            raise ValueError("行列は特異です。")

        # 行の入れ替え
        Ab[i], Ab[max_element] = Ab[max_element], Ab[i].copy()

        # ピボット要素を1にする
        Ab[i] = Ab[i] / Ab[i, i]

        # ピボット列の他の要素を0にする
        for j in range(n):
            if i != j:
                Ab[j] -= Ab[j, i] * Ab[i]

    return Ab[:, -1]  # 最後の列が解

# 使用例
A = np.array([[1, 1, 1], [0, 2, 5], [2, 5, -1]], dtype=float)
b = np.array([6, -4, 27], dtype=float)

solution = gauss_jordan(A, b)
print("解:", solution)

import numpy as np

# 行列の定義
A = np.array([[1, 2], [3, 4]])

# 逆行列の計算
A_inv = np.linalg.inv(A)

print("元の行列:")
print(A)
print("\n逆行列:")
print(A_inv)

# 検証：A * A_inv が単位行列になることを確認
I = np.dot(A, A_inv)
print("\n元の行列 * 逆行列:")
print(I)

# 丸め誤差を考慮した単位行列との比較
print("\n単位行列との差:")
print(np.allclose(I, np.eye(2)))

元の行列:
[[1 2]
 [3 4]]

逆行列:
[[-2.   1. ]
 [ 1.5 -0.5]]

元の行列 * 逆行列:
[[1.0000000e+00 0.0000000e+00]
 [8.8817842e-16 1.0000000e+00]]

単位行列との差:
True

def matrix_inverse(A):
    n = len(A)
    # 単位行列を作成
    I = [[1 if i == j else 0 for j in range(n)] for i in range(n)]

    # 行列Aに単位行列Iを右側に結合
    for i in range(n):
        A[i].extend(I[i])

    # 掃き出し法の適用
    for i in range(n):
        # ピボット選択
        max_element = max(range(i, n), key=lambda x: abs(A[x][i]))
        if A[max_element][i] == 0:
            raise ValueError("行列は特異です。逆行列は存在しません。")

        # 行の入れ替え
        A[i], A[max_element] = A[max_element], A[i]

        # ピボット要素を1にする
        pivot = A[i][i]
        for j in range(i, 2*n):
            A[i][j] /= pivot

        # ピボット列の他の要素を0にする
        for k in range(n):
            if k != i:
                factor = A[k][i]
                for j in range(i, 2*n):
                    A[k][j] -= factor * A[i][j]

    # 逆行列の抽出
    return [row[n:] for row in A]

# 使用例
A = [[1, 2], [3, 4]]
A_inv = matrix_inverse(A)

print("逆行列:")
for row in A_inv:
    print(row)

# 検証
def matrix_multiply(A, B):
    return [[sum(a*b for a,b in zip(A_row,B_col)) for B_col in zip(*B)] for A_row in A]

I = matrix_multiply(A, A_inv)
print("\n元の行列 * 逆行列:")
for row in I:
    print(row)

逆行列:
[-2.0, 1.0]
[1.5, -0.5]

元の行列 * 逆行列:
[1.0, 0.0]
[0.0, 1.0]

import numpy as np

def solve_linear_system(A, b):
    """
    連立方程式 Ax = b を解くソルバー
    A: 係数行列
    b: 定数項ベクトル
    戻り値: 解ベクトル x
    """
    return np.linalg.solve(A, b)

# 使用例
A = np.array([[2, 1, -1],
              [1, 3, 2],
              [-1, 2, 4]])

b = np.array([8, 14, 8])

solution = solve_linear_system(A, b)

print("連立方程式の解:")
print(solution)

# 解の検証
print("\n検証 (Ax = b):")
print(np.dot(A, solution))
print("b:")
print(b)

連立方程式の解:
[2. 3. 1.]

検証 (Ax = b):
[ 8. 14.  8.]
b:
[ 8 14  8]

import numpy as np
import time

def generate_large_system(n):
    """
    n x n の大規模な連立方程式を生成
    """
    A = np.random.rand(n, n)
    x = np.random.rand(n)
    b = np.dot(A, x)
    return A, b, x

def solve_and_verify(A, b, x_true):
    """
    連立方程式を解き、結果を検証
    """
    start_time = time.time()
    x_solved = np.linalg.solve(A, b)
    end_time = time.time()

    print(f"計算時間: {end_time - start_time:.6f} 秒")
    print(f"最大誤差: {np.max(np.abs(x_solved - x_true)):.6e}")

# 100x100 の連立方程式を生成して解く
n = 100
A, b, x_true = generate_large_system(n)

print(f"{n}x{n} の連立方程式を解いています...")
solve_and_verify(A, b, x_true)

# 1000x1000 の連立方程式を生成して解く
n = 1000
A, b, x_true = generate_large_system(n)

print(f"\n{n}x{n} の連立方程式を解いています...")
solve_and_verify(A, b, x_true)

100x100 の連立方程式を解いています...
計算時間: 0.001001 秒
最大誤差: 1.776357e-15

1000x1000 の連立方程式を解いています...
計算時間: 0.124002 秒
最大誤差: 1.421085e-14

import numpy as np
import multiprocessing as mp
import time

def parallel_row_operation(args):
    row, pivot_row, pivot_col, matrix = args
    if row != pivot_col:
        factor = matrix[row][pivot_col] / matrix[pivot_col][pivot_col]
        matrix[row] = matrix[row] - factor * matrix[pivot_col]
    return row, matrix[row]

def parallel_gauss_jordan(A, b):
    n = len(A)
    matrix = np.column_stack((A, b))

    for pivot_col in range(n):
        pivot_row = pivot_col
        matrix[pivot_row] = matrix[pivot_row] / matrix[pivot_row][pivot_col]

        pool = mp.Pool(processes=mp.cpu_count())
        args = [(row, pivot_row, pivot_col, matrix) for row in range(n)]
        results = pool.map(parallel_row_operation, args)
        pool.close()
        pool.join()

        for row, updated_row in results:
            matrix[row] = updated_row

    return matrix[:, -1]

# 使用例
A = np.array([[2, 1, -1], [1, 3, 2], [-1, 2, 4]], dtype=float)
b = np.array([8, 14, 8], dtype=float)

start_time = time.time()
solution = parallel_gauss_jordan(A, b)
end_time = time.time()

print("並列処理による解:", solution)
print(f"計算時間: {end_time - start_time:.6f} 秒")

並列処理による解: [2. 3. 1.]
計算時間: 0.116001 秒

import numpy as np
from scipy import sparse
from scipy.sparse.linalg import spsolve
import time

def scipy_sparse_solver(A, b):
    """
    SciPyの疎行列ソルバーを使用して連立方程式を解く
    """
    A_sparse = sparse.csr_matrix(A)
    return spsolve(A_sparse, b)

# 大規模な疎行列を生成
n = 10000
density = 0.01  # 1%の要素のみが非ゼロ

A = sparse.random(n, n, density=density).toarray()
x_true = np.random.rand(n)
b = A.dot(x_true)

# NumPyの通常のソルバーで解く
start_time = time.time()
x_numpy = np.linalg.solve(A, b)
numpy_time = time.time() - start_time
print(f"NumPy 計算時間: {numpy_time:.6f} 秒")
print(f"NumPy 最大誤差: {np.max(np.abs(x_numpy - x_true)):.6e}")

# SciPyの疎行列ソルバーで解く
start_time = time.time()
x_scipy = scipy_sparse_solver(A, b)
scipy_time = time.time() - start_time
print(f"\nSciPy 計算時間: {scipy_time:.6f} 秒")
print(f"SciPy 最大誤差: {np.max(np.abs(x_scipy - x_true)):.6e}")

print(f"\n速度向上率: {numpy_time / scipy_time:.2f}倍")

NumPy 計算時間: 32.265604 秒
NumPy 最大誤差: 3.814697e-12

SciPy 計算時間: 0.624001 秒
SciPy 最大誤差: 1.937151e-13

速度向上率: 51.71倍

import numpy as np

def process_large_matrix(matrix, chunk_size=1000):
    """
    大きな行列を分割して処理する
    """
    n = matrix.shape[0]
    result = np.zeros_like(matrix)

    for i in range(0, n, chunk_size):
        chunk = matrix[i:i+chunk_size, :]
        # ここで chunk に対して何らかの処理を行う
        processed_chunk = chunk * 2  # 例として各要素を2倍にする
        result[i:i+chunk_size, :] = processed_chunk

    return result

# 大きな行列を生成
n = 10000
m = 5000
large_matrix = np.random.rand(n, m)

# 分割処理を実行
result = process_large_matrix(large_matrix)

print(f"処理前の行列の一部:\n{large_matrix[:5, :5]}")
print(f"\n処理後の行列の一部:\n{result[:5, :5]}")

処理前の行列の一部:
[[0.51182162 0.9504637  0.14415961 0.94864945 0.31183145]
 [0.42332645 0.82770259 0.40919914 0.54959369 0.02755911]
 [0.75850313 0.75572241 0.13670659 0.57545285 0.26965907]
 [0.54688547 0.94080982 0.89673092 0.41428296 0.14335329]
 [0.94466892 0.52184832 0.41466194 0.26455561 0.77423369]]

処理後の行列の一部:
[[1.02364324 1.9009274  0.28831922 1.8972989  0.6236629 ]
 [0.8466529  1.65540518 0.81839828 1.09918738 0.05511822]
 [1.51700626 1.51144482 0.27341318 1.1509057  0.53931814]
 [1.09377094 1.88161964 1.79346184 0.82856592 0.28670658]
 [1.88933784 1.04369664 0.82932388 0.52911122 1.54846738]]

from decimal import Decimal, getcontext

# 精度を設定（ここでは100桁）
getcontext().prec = 100

def high_precision_gauss_jordan(A, b):
    n = len(A)
    # 行列をDecimal型に変換
    A = [[Decimal(str(a)) for a in row] for row in A]
    b = [Decimal(str(x)) for x in b]

    for i in range(n):
        # ピボット選択
        max_element = max(range(i, n), key=lambda x: abs(A[x][i]))
        A[i], A[max_element] = A[max_element], A[i]
        b[i], b[max_element] = b[max_element], b[i]

        pivot = A[i][i]
        for j in range(i, n):
            A[i][j] /= pivot
        b[i] /= pivot

        for k in range(n):
            if k != i:
                factor = A[k][i]
                for j in range(i, n):
                    A[k][j] -= factor * A[i][j]
                b[k] -= factor * b[i]

    return b

# 使用例
A = [[1, 1, 1], [0, 2, 5], [2, 5, -1]]
b = [6, -4, 27]

solution = high_precision_gauss_jordan(A, b)
print("高精度での解:")
for x in solution:
    print(f"{x:.50f}")

高精度での解:
5.00000000000000000000000000000000000000000000000000
3.00000000000000000000000000000000000000000000000000
-2.00000000000000000000000000000000000000000000000000

import numpy as np

def pseudo_inverse_solver(A, b):
    # NumPyのpinv関数で擬似逆行列を計算
    A_pseudo_inv = np.linalg.pinv(A)
    # 擬似逆行列を使って解を求める
    x = np.dot(A_pseudo_inv, b)
    return x

# 特異行列の例
A = np.array([[1, 2, 3],
              [2, 4, 6],
              [3, 6, 9]])
b = np.array([1, 2, 3])

try:
    # 通常の逆行列計算（エラーになる）
    x_normal = np.linalg.solve(A, b)
except np.linalg.LinAlgError as e:
    print(f"エラー: {e}")

# 擬似逆行列を使用した解法
x_pseudo = pseudo_inverse_solver(A, b)
print("擬似逆行列を使用した解:")
print(x_pseudo)

# 解の検証
print("\n解の検証 (Ax ≈ b):")
print(np.dot(A, x_pseudo))
print("元のb:")
print(b)

エラー: Singular matrix
擬似逆行列を使用した解:
[0.         0.         0.33333333]

解の検証 (Ax ≈ b):
[1. 2. 3.]
元のb:
[1 2 3]

import numpy as np
import os

def solve_large_system_with_memmap(n, chunk_size=1000):
    # 一時ファイルの作成
    filename = 'large_matrix.dat'

    # メモリマップされた行列の作成
    A = np.memmap(filename, dtype='float32', mode='w+', shape=(n, n))
    b = np.memmap(filename + '.b', dtype='float32', mode='w+', shape=(n,))

    # 行列の初期化（ここでは乱数で埋める）
    for i in range(0, n, chunk_size):
        end = min(i + chunk_size, n)
        A[i:end] = np.random.rand(end-i, n)
        b[i:end] = np.random.rand(end-i)

    # LU分解を使用して解く
    x = np.linalg.solve(A, b)

    # 一時ファイルの削除
    os.unlink(filename)
    os.unlink(filename + '.b')

    return x

# 大規模な方程式を解く
n = 20000  # 20000 x 20000 の行列
x = solve_large_system_with_memmap(n)

print(f"解の一部: {x[:5]}")
print(f"解の形状: {x.shape}")

解の一部: [ 0.00332755 -0.00083189  0.00165945  0.00041486 -0.00124945]
解の形状: (20000,)

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy import misc

# 画像の読み込みと前処理
image = misc.face(gray=True)
image = image.astype(float) / 255.0

# ぼかしカーネルの定義
kernel = np.array([[0.05, 0.1, 0.05],
                   [0.1,  0.4, 0.1],
                   [0.05, 0.1, 0.05]])

# 画像をぼかす
blurred = signal.convolve2d(image, kernel, mode='same')

# 逆畳み込みのための行列を作成
def create_convolution_matrix(kernel, image_shape):
    m, n = image_shape
    k = kernel.shape[0] // 2
    conv_matrix = np.zeros((m*n, m*n))

    for i in range(m):
        for j in range(n):
            row = i * n + j
            for ki in range(-k, k+1):
                for kj in range(-k, k+1):
                    if 0 <= i+ki < m and 0 <= j+kj < n:
                        col = (i+ki) * n + (j+kj)
                        conv_matrix[row, col] = kernel[ki+k, kj+k]

    return conv_matrix

# 逆畳み込みの実行
conv_matrix = create_convolution_matrix(kernel, image.shape)
blurred_flat = blurred.flatten()
restored_flat = np.linalg.solve(conv_matrix, blurred_flat)
restored = restored_flat.reshape(image.shape)

# 結果の表示
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 5))
ax1.imshow(image, cmap='gray')
ax1.set_title('Original')
ax2.imshow(blurred, cmap='gray')
ax2.set_title('Blurred')
ax3.imshow(restored, cmap='gray')
ax3.set_title('Restored')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# データの生成
np.random.seed(0)
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)

# デザイン行列の作成
X_b = np.c_[np.ones((100, 1)), X]

# 正規方程式の解法（掃き出し法を使用）
theta = np.linalg.solve(X_b.T.dot(X_b), X_b.T.dot(y))

print("推定されたパラメータ:")
print(f"切片: {theta[0][0]:.4f}")
print(f"傾き: {theta[1][0]:.4f}")

# 結果のプロット
plt.scatter(X, y)
plt.plot(X, X_b.dot(theta), color='r', linewidth=2)
plt.xlabel("X")
plt.ylabel("y")
plt.title("Linear Regression using Gaussian Elimination")
plt.show()

import numpy as np

def decrypt_substitution_cipher(ciphertext, known_mappings):
    # アルファベットのリスト
    alphabet = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
    n = len(alphabet)

    # 連立方程式の係数行列と定数項を初期化
    A = np.zeros((n, n))
    b = np.zeros(n)

    # 既知の対応関係から方程式を作成
    for plain, cipher in known_mappings.items():
        i = alphabet.index(plain)
        j = alphabet.index(cipher)
        A[i, j] = 1
        b[i] = 1

    # 残りの未知の文字に対する方程式を追加
    for i in range(n):
        if np.sum(A[i, :]) == 0:
            A[i, i] = 1

    # 掃き出し法で連立方程式を解く
    x = np.linalg.solve(A, b)

    # 解読マッピングを作成
    decryption_map = {}
    for i, prob in enumerate(x):
        if prob > 0.5:
            decryption_map[alphabet[i]] = alphabet[np.argmax(A[:, i])]

    # 暗号文を解読
    plaintext = ''.join(decryption_map.get(c, c) for c in ciphertext)

    return plaintext

# 使用例
ciphertext = "KHOOR ZRUOG"
known_mappings = {'H': 'E', 'O': 'L'}

decrypted_text = decrypt_substitution_cipher(ciphertext, known_mappings)
print(f"暗号文: {ciphertext}")
print(f"解読文: {decrypted_text}")

暗号文: KHOOR ZRUOG
解読文: HELLO WORLD