【数値解析】定常反復法の殴り書き実装

やかんです。

明後日、最適化手法と数値解析の試験があります。今更感はすごいんですけど、やっぱり勉強してると実装した方がいいなという局面があり、今がそれです。

というわけで、線型方程式（連立一次方程式）の反復法的求解を試みたいと思います。

いつ反復法を使うのか。

授業では、

• 小さくて密 → 直接法
• 大きくて疎 → 反復法

という説明がありました。小さくて密な行列についてAx = bを解く局面がいつ訪れるのかわからないので直接法をいつ使うのかピンときていませんが、とりあえず反復法を実装してみたいです。

Jacobi法

結論、うまく動きませんでした。不慣れなpythonをザザザっと勢いで実装したので、そりゃうまく動かないよねといったところではあります。

あと、「大きくて疎な行列に適しているよ！」という反復法ですが、プログラム上、小さくて密な行列で実験しちゃってます。色々欠陥だらけの実装ですね。。

import numpy as np

SIZE = 5
LOOP = 1000

matrix = np.random.randint(1, 10, (SIZE, SIZE)).astype(float)  
b = np.random.randint(1, 10, SIZE).astype(float).reshape(SIZE, 1)
print(matrix)
print(b)

まず、Ax = bのAとbを用意します。変数名もテキトーですごめんなさい。

U = np.zeros((SIZE,SIZE), dtype=float)
D = np.zeros((SIZE,SIZE), dtype=float)
L = np.zeros((SIZE,SIZE), dtype=float)

for i in range(len(matrix)):
  for j in range(len(matrix[i])):
    if(i < j):
      # 右上三角(U)
      U[i][j] = matrix[i][j]
    elif(i == j):
      # 対角成分(D)
      D[i][j] = matrix[i][j]
    else:
      # 左下三角(L)
      L[i][j] = matrix[i][j]

print("---- U ----")
print(U)

print("---- D ----")
print(D)

print("---- L ----")
print(L)

print("---- U + D+ L ----")
print(U + D + L)

Aから対角成分を取り出しました。

# こっちJacobi

D_inv = np.linalg.inv(D)
D_inv_b = D_inv @ b

base_matrix = -D_inv @ (L + U)

def Jacobi():
  x = np.zeros((SIZE, 1), dtype=float)
  
  for i in range(LOOP):
    x_new = np.zeros((SIZE, 1), dtype=float)
    for j in range(len(x)):
      tmp = 0
      for k in range(len(base_matrix)):
        tmp += base_matrix[j][k] * x[k]
      tmp += D_inv_b[j]
      x_new[j] = tmp

    x = x_new.copy()

  print(x)
  
  return x

はい、このJacobi関数が問題の関数ですが、もうすごい雑でごめんなさい。動きませんこのコードは。「逆行列使う時はLU分解を経由しましょう」と習いましたが、numpyが内部的にやってくれていたりしないんでしょうか。

Gauss-Siedel法

ほぼJacobi（もどき）と同じです。

# Gauss-Seidel

def GaussSeidel():
  x = np.zeros((SIZE, 1), dtype=float)

  for i in range(LOOP):
      for j in range(SIZE):
          tmp = 0
          for k in range(j):  # **既に計算した新しい値を使用**
              tmp += base_matrix[j][k] * x[k]  
          for k in range(j+1, SIZE):  # **まだ更新されていない古い値を使用**
              tmp += base_matrix[j][k] * x[k]  
          tmp += D_inv_b[j]
          x[j] = tmp  # **上書き更新**

  print(x)
  return x

実行してみた。

solution = 0

x = np.zeros((SIZE, 1), dtype=float)
if(solution == 0):
  x = Jacobi()
elif(solution == 1):
  x = GaussSeidel()


print("---- x ----")
print(x)

print("---- Ax ----")
print(matrix @ x)

print("---- b ----")
print(b)

出力がこちら↓

solution = 0

x = np.zeros((SIZE, 1), dtype=float)
if(solution == 0):
  x = Jacobi()
elif(solution == 1):
  x = GaussSeidel()


print("---- x ----")
print(x)

print("---- Ax ----")
print(matrix @ x)

print("---- b ----")
print(b)
[[nan]
 [nan]
 [nan]
 [nan]
 [nan]]
---- x ----
[[nan]
 [nan]
 [nan]
 [nan]
 [nan]]
---- Ax ----
[[nan]
 [nan]
 [nan]
 [nan]
 [nan]]
---- b ----
[[2.]
 [5.]
 [9.]
 [2.]
 [5.]]

はい！ごめんなさい！

ということで、もはやネタ記事枠ですが、以上です。最後までお読みいただき、ありがとうございます。