PLS回帰 (PLS1) を理解する

<script type="text/javascript" async 
  src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.7/MathJax.js?config=TeX-MML-AM_CHTML">
</script>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
  tex2jax: {
    inlineMath: [['$', '$']],
    displayMath: [['$$', '$$'], ["\\[", "\\]"]]
  }
});
</script>

# はじめに
　[スモールデータと機械学習](https://www.ohmsha.co.jp/book/9784274227783/)より、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$)  を予測するモデル 
```math
\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潜在変数ベクトル) が大きくなるようにする
```math
\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$) をかけて作成されるとします。
 ```math
 t_1 = X w_1
 ```
 $y$と$t_1$ の共分散 ($\sigma_{y, t_1}$) を最大化します。$\frac{1}{N-1}$ は固定なので記載を省略します。
 ```math
 \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$ に制約を加えます。
 ```math
 w_1^\mathsf{T} w_1 = ||w_1|| = 1
 ```
$||w_1|| = 1$ の制約下で$\sigma_{y, t_1}$ の極値を探すため、ラグランジュの未定乗数法で解きます。ラグランジュ乗数を $\lambda$ とすると、目的関数 $J_1$は
```math
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$ で偏微分することで
```math
\frac{\delta J_1}{\delta w_1} = \underbrace{X^\mathsf{T} y}_{w_1 と同じ列ベクトルの形に変換} - 2\lambda w_1 = 0 
```
```math
X^\mathsf{T} y = 2\lambda w_1
```
となります。そのため、$\sigma_{y,t_1}$ は
```math
\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}$ から
```math
\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$ は
```math
w_1 = \frac{X^\mathsf{T} y}{||X^\mathsf{T} y||}
```
となり、$t_1 = Xw_1$ を計算できるようになります。
$y$ を予測するため、$t_1$ の回帰係数 $d_1$ を考えます。
```math
y_{pred} = d_1 t_1
```
最小二乗法によって回帰係数 $d_1$ は以下のように求まります。
```math
\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}
```
<details>
<summary>最小二乗法による回帰係数の導出</summary>
最小二乗法における重回帰の回帰係数(β) は予測値 と実測値(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}$$
</details>

よって、第一潜在変数 $t_1$ を用いた予測値 $y_{pred}$ は以下のようにあらわせます。
```math
\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$ は以下のように導出されます。
```math
\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$ )を抽出します。この操作はデフレーションと呼ばれます。
```math
\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`](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_diabetes.html#sklearn.datasets.load_diabetes) を使用します。特徴量の数は 10、データ数は 442 です。第 4 潜在変数まで計算してみます。
```py
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)
```
<img src = "https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3230781/c8f61d78-60ed-4a8d-a255-b5a6127400de.png" width ="500">

第 3 潜在変数までで回帰精度が頭打ちになっていそうです。

## scikit-learn ライブラリとの比較
[`sklearn.cross_decomposition.PLSRegression`](https://scikit-learn.org/stable/modules/generated/sklearn.cross_decomposition.PLSRegression.html) の結果と比較してみます。
```py
# 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()
```
 <img src = "https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3230781/162c168f-f990-4b16-9fa9-c91eec5dd0ff.png
" width ="500">
sklearn との結果とも重なったので、実装面でも問題なさそうです。
# おわりに
　PLS1 を理解するために、自身の学びを出力しました。違っている点など見つかれば適宜修正します。ご指摘いただけると感謝です。