はじめに
スモールデータと機械学習より、PLS回帰の理解を深めるために記事を作成します。今回は1つの目的変数を予測するPLS1 についてです。n番煎じです。この記事では基本的にベクトルは縦ベクトルとして扱います。
PLS 回帰の概要
- データセット ($X, X \in \mathbb{R}^{N \times M}$, $N$: データ数, $M$: 説明変数の次元) から潜在変数 ($T, T \in \mathbb{R}^{N \times R}$, $R$: 潜在変数の次元) を作成し、潜在変数を用いて目的変数 ($Y, Y \in \mathbb{R}^N$) を予測するモデル
\begin{align}
T
&= X\underbrace{w}_{重み} \\
y_{pred}
&= \underbrace{d}_{T に対する回帰係数} T
\end{align}
- 目的変数 ($Y$) をよく予測できるような潜在変数 ($T$) を作成する
➡$Y$ と $T$ の共分散 ($\sigma_{y, t_1}$, $t_1$: $T$ の第1潜在変数ベクトル) が大きくなるようにする
\sigma_{y, t_1} = \frac{1}{N-1} y^\mathsf{T} t_1
PLSを行って何がうれしいのか
- 説明変数間に多重共線性があってもモデルを作成できる
- 説明変数に対する次元削減の一助になる
導出
第1潜在変数ベクトル
$t_1$ が$X$ に適当な重みベクトル ($w_1, w_1 \in \mathbb{R}^M$) をかけて作成されるとします。
t_1 = X w_1
$y$と$t_1$ の共分散 ($\sigma_{y, t_1}$) を最大化します。$\frac{1}{N-1}$ は固定なので記載を省略します。
\begin{align}
\sigma_{y, t_1}
&= y^\mathsf{T} t_1 \\
&= y^\mathsf{T} X w_1
\end{align}
この時、$\sigma_{y, t_1}$ の大きさを制限するため、$w_1$ に制約を加えます。
w_1^\mathsf{T} w_1 = ||w_1|| = 1
$||w_1|| = 1$ の制約下で$\sigma_{y, t_1}$ の極値を探すため、ラグランジュの未定乗数法で解きます。ラグランジュ乗数を $\lambda$ とすると、目的関数 $J_1$は
J_1 = y^\mathsf{T} X w_1 - \lambda(w_1^\mathsf{T} w_1 - 1)
となります。$y^\mathsf{T} X w_1$ は共分散 (スカラー) であるため、その微分は $w_1$ と同じ列ベクトルの形になります。$J_1$ を $w_1$ で偏微分することで
\frac{\delta J_1}{\delta w_1} = \underbrace{X^\mathsf{T} y}_{w_1 と同じ列ベクトルの形に変換} - 2\lambda w_1 = 0
X^\mathsf{T} y = 2\lambda w_1
となります。そのため、$\sigma_{y,t_1}$ は
\begin{align}
\sigma_{y,t_1}
&= y^\mathsf{T} X w_1 \\
&= (X^\mathsf{T} y)^\mathsf{T} w_1 \\
&= 2\lambda w_1^\mathsf{T} w_1 \\
&= 2\lambda
\end{align}
であり、$w_1 = \frac{X^\mathsf{T} y}{2\lambda}$ から
\begin{align}
\sigma_{y, t_1}
&= (X^\mathsf{T} y)^\mathsf{T} w_1 \\
&= \frac{(X^\mathsf{T} y)^\mathsf{T} X^\mathsf{T} y}{2\lambda} \\
&= \frac{||X^\mathsf{T} y||^2}{2\lambda} \\
&= 2\lambda \\
(2\lambda)^2
&= ||X^\mathsf{T} y||^2 \\
\sigma_{y, t_1}
&= ||X^\mathsf{T} y|| = 2\lambda
\end{align}
となります。よって $w_1$ は
w_1 = \frac{X^\mathsf{T} y}{||X^\mathsf{T} y||}
となり、$t_1 = Xw_1$ を計算できるようになります。
$y$ を予測するため、$t_1$ の回帰係数 $d_1$ を考えます。
y_{pred} = d_1 t_1
最小二乗法によって回帰係数 $d_1$ は以下のように求まります。
\begin{align}
d_1
&= (t_1^\mathsf{T} t_1)^{-1} t_1^\mathsf{T} y \\
&= \underbrace{\frac{t_1^\mathsf{T} y}{||t_1||^2}}_{ベクトル t_1 の内積はスカラーになるため、(t_1^\mathsf{T} t_1)^{-1} = \frac{1}{||t_1||^2}}
\end{align}
最小二乗法による回帰係数の導出
最小二乗法における重回帰の回帰係数(β) は予測値 と実測値(y) の差(Q) を最小にするβとなります。 $$\begin{align} y_{pred} &= X \beta \\ Q &= (y - X \beta)^\mathsf{T}(y - X \beta) \\ &= y^\mathsf{T} y -2y^\mathsf{T} X\beta + (X\beta)^\mathsf{T} X\beta \\ \frac{\delta Q}{\delta \beta} &= 2X^\mathsf{T} X \beta - 2 X^\mathsf{T} y &= 0 \\ \beta &= (X^\mathsf{T} X)^{-1}X^\mathsf{T} y \end{align}$$よって、第一潜在変数 $t_1$ を用いた予測値 $y_{pred}$ は以下のようにあらわせます。
\begin{align}
y_{pred}
&= d_1 t_1 \\
&= d_1 X w_1 \\
\quad
w_1
&= \frac{X^\mathsf{T} y}{||X^\mathsf{T} y||} \\
t_1
&= Xw_1 \\
d_1
&= \frac{t_1^\mathsf{T} y}{||t_1||^2}\\
\end{align}
また、$ t_1 = X w_1$ であるため、$X = t_1 p_1^\mathsf{T}$ となるローディングベクトル $p_1$ は以下のように導出されます。
\begin{align}
X
&= t_1 p_1^\mathsf{T} \\
\underbrace{t_1^\mathsf{T}}_{両辺に左からt_1^\mathsf{T}をかける} X
&= t_1^\mathsf{T} t_1 p_1^\mathsf{T} \\
p_1^\mathsf{T}
&= \frac{t_1^\mathsf{T} X}{t_1^\mathsf{T} t_1}\\
p_1
&= \frac{X^\mathsf{T} t_1}{||t_1||^2}
\end{align}
その後、第一潜在変数 $t_1$ で表現ができなった $y$、$X$ の成分( $y_2, X_2$ )を抽出します。この操作はデフレーションと呼ばれます。
\begin{align}
y_2
&= y - d_1 t_1 \\
&= y - d_1 X w_1 \\
X_2
&= X - t_1p_1^\mathsf{T} \\
&= X - t_1 \frac{t_1^\mathsf{T} X}{||t_1||^2}
\end{align}
PLS1では、デフレーション後の$y_2, X_2$ に対して、同様に第2潜在変数 $t_2$ を求めることで回帰精度の向上を図ります。
実装
関数の作成
numpy で実装してみます。データは sklearn.datasets.load_diabetes を使用します。特徴量の数は 10、データ数は 442 です。第 4 潜在変数まで計算してみます。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_diabetes
from sklearn.cross_decomposition import PLSRegression
from sklearn.metrics import r2_score, mean_squared_error
def pls1_nipals(X, y, r = 4):
w_l = []
t_l = []
p_l = []
d_l = []
y_pred_l = []
y_data = y.reshape(-1, 1)
X_data = X
# fig, ax
plt.rcParams["font.size"] = 14
fig, ax = plt.subplots(2,2, figsize=(10, 10), tight_layout = True)
ax = ax.flatten()
for i in range(r):
w = X_data.T @ y_data /(np.linalg.norm(X_data.T @ y_data))
t = X_data @ w
p = X_data.T @ t / (np.linalg.norm(t)**2)
d = t.T @ y_data / (np.linalg.norm(t)**2)
y_pred = d * t
w_l.append(w)
t_l.append(t)
p_l.append(p)
d_l.append(d)
y_pred_l.append(y_pred)
# visualizaion
y_pred_vis = np.sum(y_pred_l, axis=0)
r2 = r2_score(y, y_pred_vis)
rmse = np.sqrt(mean_squared_error(y, y_pred_vis))
ax[i].scatter(y, y_pred_vis)
ax[i].set_ylabel("y_pred")
ax[i].set_xlabel("y")
text_str = f"R² Score: {r2:.3f}\nRMSE: {rmse:.3f}"
boxdic = {
"facecolor" : "white",
"edgecolor" : "black",
"boxstyle" : "Round",
"linewidth" : 2}
ax[i].text(max(y)*0.5, min(y_pred_vis)+np.mean(y_pred_vis)/2, text_str, fontsize=14, bbox=boxdic) #,
ax[i].set_title(f"PLS1 NIPALS: Component {i+1}")
# deflation
y_data = y_data - y_pred
X_data = X_data - t @ p.T
plt.show()
return t_l, p_l, w_l, d_l
diabetes = load_diabetes()
X, y = diabetes.data, diabetes.target
# autoscaling
X = (X - np.mean(X, axis = 0))/np.std(X, axis = 0)
y = (y - np.mean(y))/np.std(y)
t, p, w, d = pls1_nipals(X, y, 4)
第 3 潜在変数までで回帰精度が頭打ちになっていそうです。
scikit-learn ライブラリとの比較
sklearn.cross_decomposition.PLSRegression の結果と比較してみます。
# pls by numpy
y_pred = np.squeeze(t).T @ np.squeeze(d)
# pls by sklearn
pls = PLSRegression(n_components=4)
pls.fit(X, y)
y_pred_sk = pls.predict(X)
plt.scatter(y, y_pred, label = "numpy")
plt.scatter(y, y_pred_sk, label = "sklearn")
plt.ylabel("y_pred")
plt.xlabel("y")
plt.legend()
plt.show()
sklearn との結果とも重なったので、実装面でも問題なさそうです。
おわりに
PLS1 を理解するために、自身の学びを出力しました。違っている点など見つかれば適宜修正します。ご指摘いただけると感謝です。