正則化による多重共線性を持ったデータの損失関数の変化を可視化してみた

## 目次
[1. はじめに](#はじめに)
[2. 多重共線性のある擬似データの生成](#多重共線性のある擬似データの生成)
[3. 損失関数の形状（正則化なし）](#損失関数の形状（正則化なし）)
[4. 損失関数の形状（正則化あり）](#損失関数の形状正則化あり)
[5. 正則化で解決できないこと](#正則化で解決できないこと)
[6. 使用したコード](#使用したコード)

## はじめに
特徴量間に相関の高い組み合わせが存在することを多重共線性といいますが、機械学習・統計の領域において特徴量の種類が多いデータを扱う場合などにはよく生じる問題です。
多重共線性があることで、モデルの重み推定が不安定になることや、過学習してしまう可能性があるなどの問題点があり、その解決方法の一つとして正則化項を損失関数に加えることが挙げられます。
正則化項を加えることで、どのように多重共線性が抱える問題が解決されるかをまとめました。

## 多重共線性のある擬似データの生成
今回、擬似データとして以下のような二次元の特徴量と、それを線形結合してノイズを加えた目的変数のデータを生成しています。（データ生成も含めたコードは全て最後に載せています。）
```math
\begin{align*}
x_{1} &= \text{rand}(100) \\
x_{2} &= 2x_{1} + 0.1\epsilon, & \epsilon \sim N(0, 1) \\
y &= x_{1} + 2x_{2} + \epsilon', & \epsilon' \sim N(0, 1)
\end{align*}
```
変数$x_{1}$, $x_{2}$間には強い相関があり、$x_{1}$, $x_{2}$の線形結合によって、$y$が表されているようなデータです。
## 損失関数の形状（正則化なし）
損失関数としてMean Squared Error(MSE)を採用しています。
```math
L = \frac{1}{N} \sum_{i=1}^{N} (y_i - \hat{y}_{i})^2
```

ここで、$ \hat{y}\_{i} = ax_{i,1} + bx_{i,2}$です。（$a,b$は係数）
$a,b$を変化させて、損失関数$L$を計算すると以下のような谷型の形状になり、$L$の極小点が定まりにくい状態になっています。
これによって、モデルの重み推定が不安定になり、場合によっては、最適化の結果極端に大きな係数が計算される可能性があります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/296599/1fdbd414-72a3-fd47-3f64-68c6a2230d74.png)

## 損失関数の形状（正則化あり）
次に、以下のように損失関数にL2正則化項を加えて計算をしてみます。
```math
L = \frac{1}{N} \sum_{i=1}^{N} (y_i - \hat{y}_{i})^2 +\alpha(a^2 + b^2)
```

$\alpha$は正則化項の強さを決めるパラメータです。
これを同様に$a,b$を変化させて計算すると以下のようになります。正則化項の存在によって、重みが大きな絶対値を取ることができなくなり、きちんと極小点が見られるようになりました。これによって、重み推定の不安定さは解消することができることが分かると思います。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/296599/23539013-ca2d-d7a2-d6e7-7de2c5bb396d.png)

## 正則化で解決できないこと
上で見たように、正則化には重みの推定値を安定化させる効果や、それに伴ってモデルが訓練データに過学習する効果を抑えられる効果があります。
一方で、正則化自体には特徴量間の存在自体を解消する効果はないので、モデルの予測精度向上には役立つ可能性がありますが、そこから得られる知見などに関して過度に期待してしまうと問題になる可能性があります。

## 使用したコード

```python
import numpy as np
import plotly.graph_objects as go
from sklearn.metrics import mean_squared_error


# 損失関数の計算（L2正則化項を含む）
def calculate_loss_surface(X, y, A_range, B_range, alpha=0):
    A_mesh, B_mesh = np.meshgrid(A_range, B_range)
    predictions = A_mesh[..., np.newaxis] * X[:, 0] + B_mesh[..., np.newaxis] * X[:, 1]
    # 平均二乗誤差を計算
    mse = ((predictions - y) ** 2).mean(axis=2)
    # 正則化項を計算
    regularisation = alpha * (A_mesh**2 + B_mesh**2)
    # 最終的な損失関数のサーフェス
    loss_surface = mse + regularisation
    print(np.min(loss_surface))
    return A_mesh, B_mesh, loss_surface


# Plotlyでサーフェスプロットを作成する関数
def create_surface_plot(A_mesh, B_mesh, loss_surface, title):
    fig = go.Figure(data=[go.Surface(z=loss_surface, x=A_mesh, y=B_mesh, colorscale='Viridis')])
    fig.update_layout(
        title=title,
        scene=dict(
            xaxis_title='a',
            yaxis_title='b',
            zaxis_title='Loss'
        ),
        autosize=False,
        width=800,
        height=800,
        margin=dict(l=65, r=50, b=65, t=90)
    )
    return fig


# データセットの生成
np.random.seed(0)  # 結果の再現性のためにシードを設定
featureA = np.random.rand(100)
featureB = 2*featureA + 0.1*np.random.randn(100)
output = 1*featureA + 2*featureB + np.random.randn(100)
X = np.column_stack((featureA, featureB))
y = output

# 係数の範囲を設定
coef_A_range = np.linspace(-3, 5, 100)
coef_B_range = np.linspace(-3, 5, 100)

# 損失関数のサーフェスを計算（L2正則化なし）
A_mesh, B_mesh, loss_surface_no_reg = calculate_loss_surface(X, y, coef_A_range, coef_B_range)

# 損失関数のサーフェスを計算（L2正則化あり）
alpha = 3  # L2正則化の強さを設定
A_mesh, B_mesh, loss_surface_with_reg = calculate_loss_surface(X, y, coef_A_range, coef_B_range, alpha=alpha)

# L2正則化なしのサーフェスプロットを作成
fig_no_reg = create_surface_plot(A_mesh, B_mesh, loss_surface_no_reg, 'Loss Surface without Regularization')

# L2正則化ありのサーフェスプロットを作成
fig_with_reg = create_surface_plot(A_mesh, B_mesh, loss_surface_with_reg, 'Loss Surface with L2 Regularization')

fig_no_reg.show()
fig_with_reg.show()

```