行列分解法による線形回帰をPyTorchで解く

線形回帰はいろんな方法で行えますが、この記事では行列分解法（matrix decomposition method）を PyTorch で実装してみます。

行列分解法（matrix decomposition method）

行列分解法については、たとえばこちらを参照していただくと解説があります（Tensorflowによる実装もあります）。

PyTorch による行列分解法

必要なライブラリのインポート

import matplotlib.pyplot as plt
import numpy as np
import torch

PyTorchの基礎は、たとえばこちら。

データを作成します。

x_vals = np.linspace(0, 10, 100)
y_vals = x_vals + np.random.normal(0, 1, 100)

逆行列法に用いる行列を作成します。

# xデータの列と1の列から成る行列Aを作成
x_vals_column = np.transpose(np.matrix(x_vals))
ones_column = np.transpose(np.matrix(np.repeat(1, 100)))
A = np.column_stack((x_vals_column, ones_column))

# yデータの行列bを作成
b = np.transpose(np.matrix(y_vals))

テンソルに変換します

# テンソルを作成
A_tensor = torch.from_numpy(A).float()
b_tensor = torch.from_numpy(b).float()

ここまでは、逆行列法による線形回帰をPyTorchで解くと全く同じです。ここからは違います。

正方行列のコレスキー分解を実行します。

# これ好きー分解
tA_A = torch.matmul(torch.transpose(A_tensor, 0, 1), A_tensor)
L = torch.linalg.cholesky(tA_A)

# L * y = t(A) * b を解く
tA_b = torch.matmul(torch.transpose(A_tensor, 0, 1), b_tensor)
sol1 = torch.linalg.solve(L, tA_b)

# L' * y = sol1 を解く
solution = torch.linalg.solve(torch.transpose(L, 0, 1), sol1)

ここから先も、逆行列法による線形回帰をPyTorchで解くと一緒です。

# 係数を抽出
slope = solution[0][0]
y_intercept = solution[1][0]
slope, y_intercept

(tensor(1.0033), tensor(0.0988))

# 回帰直線
best_fit = []
for i in x_vals:
    best_fit.append(slope * i + y_intercept)

# 結果をプロット
plt.scatter(x_vals, y_vals, label="data")
plt.plot(x_vals, best_fit, label="fitted line")
plt.legend()
plt.show()

行列が大きい場合などに、行列分解法のほうが計算効率が良いそうです。