Lassoの誤差関数を可視化してみる

# Lasso の誤差関数を可視化

Lasso Regression の誤差関数を自分なりにわかりやすく可視化してみました。

```Python
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
```

2次元の説明変数を適当に作る
```Python
np.random.seed(0)
x = np.random.rand(2, 10) - 0.5
```

変数を適当に決めて次のような y を作成する
$y = x[0] \times a[0] + x[1] \times a[1] + {\rm noise}$
```Python
a = [-8, 2]  # 係数を適当に決める
n = 0.2

# y = x[0] * a[0] + x[1] * a[1] + noise
y = (x.T @ a) + np.random.rand(*x[0].shape) * n

# lassoの係数
lambda_ = 0.5
```

プロットするために mesgrid で-10~10 の範囲で誤差関数を算出する
```Python
test_a = np.linspace(-10, 10, 1001)
test_a1, test_a2 = np.meshgrid(test_a, test_a)  # 係数を a_1, a_2 で -10~10 を 1001 分割で用意

avesq_e = np.empty_like(test_a1)
lasso_e = np.empty_like(test_a1)

for i, (a1, a2) in enumerate(zip(test_a1, test_a2)):
    # 最小二乗誤差
    differ = (x.T@np.array([a1, a2])).T - y  # 誤差
    avesq_e[i] = (differ**2).sum(axis=1)/x.shape[1]  # 2乗平均

    # 正則化項の誤差
    abs_coef = np.abs(np.array([a1, a2]))
    lasso_e[i] = lambda_ * np.sum(abs_coef, axis=0)

test_e = avesq_e + lasso_e  # Lasso の誤差
```

各誤差勾配の最適な係数を求める
```Python
max_point_a = np.unravel_index(np.argmax(-avesq_e), avesq_e.shape)
max_point_l = np.unravel_index(np.argmax(-lasso_e), lasso_e.shape)
max_point = np.unravel_index(np.argmax(-test_e), test_e.shape)

print("重回帰分析の係数", test_a1[max_point_a], test_a2[max_point_a])
print("L1ノルムの中心", test_a1[max_point_l], test_a2[max_point_l])
print("Lassoの最適解", test_a1[max_point], test_a2[max_point])
```

> 重回帰分析の係数 -7.640000000000001 2.08
> L1ノルムの中心 0.0 0.0
> Lassoの最適解 -1.2200000000000006 0.0

# 等高線の表示
上の結果をもとに等高線をプロットする

```Python
fig = plt.figure(figsize=(10, 5))

levels = np.linspace(-50, 0, 51)

ax_t = fig.add_subplot(1, 3, 1)
ax_t.contour(test_a1, test_a2, -avesq_e, levels=levels)
ax_t.scatter(test_a1[max_point_a], test_a2[max_point_a], color='k')
ax_t.set_title('Mean squared error')

ax_l = fig.add_subplot(1, 3, 2)
ax_l.contour(test_a1, test_a2, -lasso_e, levels=levels)
ax_l.scatter(test_a1[max_point_l], test_a2[max_point_l], color='k')
ax_l.set_title('Regularization (L1 Norm)')

ax_e = fig.add_subplot(1, 3, 3)
ax_e.contour(test_a1, test_a2, -test_e, levels=levels)
ax_e.scatter(test_a1[max_point], test_a2[max_point], color='k')
ax_e.set_title('MSE + L1')

plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/662712/a8e37083-3de4-0b4e-737b-91e20b70caf3.png)

+ Mean Squared Error で示される最適解は, 最初に設定した $a = (-8, 2)$ を示してる。
+ L1 Norm は $(0, 0)$ がもちろん最適解. 軸上 ($a_0 = 0$ or $a_1 = 0$) が尖っている
+ MSE + L1 ではこの2つの勾配が足し合わされ, L1 Normの影響が強いときは軸上に点がある(最適な係数が一部 $0$ になる→スパース性)

# 3D の表示

```Python
%matplotlib widget

# 3D で誤差関数をプロットする

plot_y = -test_e

fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.plot_surface(test_a1, test_a2, plot_y, cmap=cm.plasma, alpha=0.8)
ax.contourf(test_a1, test_a2, plot_y, levels=30, cmap=cm.Greys)
ax.scatter(test_a1[max_point], test_a2[max_point], plot_y[max_point], color='k')

ax.set_xlabel('a_0')
ax.set_ylabel('a_1')
ax.set_zlabel('y')
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/662712/16f0460b-46ed-6c4e-b795-3faebf5cb940.png)

3D でぐるぐる回して確認できる

# scikit-learn で解を求める
```Python
from sklearn.linear_model import Lasso  # scikit-learn の Lasso モジュール

alpha = lambda_*(1 / (2*x.shape[1])) * x.shape[1]  # sklearn の定義に合わせて正則化項の係数を算出
lasso_model = Lasso(alpha=alpha, fit_intercept=False)  # 上のグラフに合わせるため "fit_intercept=False" とする

lasso_model.fit(x.T, y)
print("係数:", *lasso_model.coef_)  # 上の 'Lassoの最適解' と同じ結果が得られる
```
> 係数: -1.214420120567555 0.0

scikit-learnで求めた結果は, 上の結果
> Lassoの最適解 -1.2200000000000006 0.0

とほぼ同じものが求められる

最初に設定した, a, lambda_ を変えても同じ結果が得られる

# まとめ
使われがちだけれどわかりづらい Lasso を可視化してみた。
変数をいろいろ変えてみることでより理解が深まった

# 履歴
初稿: 2023-11-06