【グラフでわかる】正則化がどのように多重共線性を持ったデータに対する学習を安定させているか【コード解説付き】

機械学習・データサイエンス

※ 本記事はQiitaに投稿した記事に加筆をしたものです。

はじめに

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

多重共線性のある擬似データの生成

今回、擬似データとして以下のような二次元の特徴量と、それを線形結合してノイズを加えた目的変数のデータを生成しています。(データ生成も含めたコードは全て最後に載せています。)

$$
\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)を採用しています。

$$
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\)の極小点が定まりにくい状態になっています。これによって、モデルの重み推定が不安定になり、場合によっては、最適化の結果極端に大きな係数が計算される可能性があります。

損失関数の形状(正則化あり)

次に、以下のように損失関数にL2正則化項を加えて計算をしてみます。

$$
L = \frac{1}{N} \sum_{i=1}^{N} (y_i – \hat{y}_{i})^2 +\alpha(a^2 + b^2)
$$

\(\alpha\)は正則化項の強さを決めるパラメータです。
これを同様に\(a,b\)を変化させて計算すると以下のようになります。正則化項の存在によって、重みが大きな絶対値を取ることができなくなり、きちんと極小点が見られるようになりました。これによって、重み推定の不安定さは解消することができることが分かると思います。

正則化で解決できないこと

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

使用したコード

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()

コード解説

使用コードの中で、

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)

部分の処理の流れが分かりづらいと思うので解説したいと思います。
具体的に以下のコードを例に挙げてみます(ndarrrayのどの軸がどの値に対応しているかわかりやすくするために、わざと形状を非対称にしています)。

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

# 係数の範囲を設定
coef_A_range = np.linspace(0, 3, 3)
coef_B_range = np.linspace(0, 3, 4)
A_mesh, B_mesh = np.meshgrid(coef_A_range, coef_B_range)

predictions = A_mesh[..., np.newaxis] * X[:, 0] + B_mesh[..., np.newaxis] * X[:, 1]
mse = ((predictions - y) ** 2).mean(axis=2)

np.meshgrid(coef_A_range, coef_B_range)

最初に、以下のコードの説明です。

coef_A_range = np.linspace(0, 3, 3)
coef_B_range = np.linspace(0, 3, 4)
A_mesh, B_mesh = np.meshgrid(coef_A_range, coef_B_range)

np.meshgridは、与えられた配列から各要素の組み合わせを網羅的に探索したい場合に使うメソッドで、今回の場合、A_mesh, B_meshは以下のような配列になります。


また、この二つの配列は以下のように、coef_A_range, coef_B_rangeに含まれるすべての要素の組み合わせを網羅していて、例として図中のA_mesh[2][1], B_mesh[2][1]のようにすることで各点にアクセスができます。

A_mesh[..., np.newaxis] * X[:, 0] + B_mesh[..., np.newaxis] * X[:, 1]

次に、今回のコードの中で一番わかりづらい箇所だと思いますが、predictionsを計算する箇所について説明します。
まず、A_mesh[..., np.newaxis]は、A_meshに新しい軸を追加する処理なのですが、これによって配列が以下のように変形します。(axis=0方向は全て同じ配列[a0, a1, a2]です。)


次に、A_mesh[..., np.newaxis] * X[:, 0]の処理では、ブロードキャストが行われることで以下のように配列がそれぞれ変化します。


同様に、B_meshB_mesh[..., np.newaxis]によって以下のように変換され、

B_mesh[..., np.newaxis] * X[:, 1]は以下のようにブロードキャストされて計算されます。

結果としてA_mesh[..., np.newaxis] * X[:, 0] + B_mesh[..., np.newaxis] * X[:, 1]は以下のような計算結果となります。

mse = ((predictions - y) ** 2).mean(axis=2)

最後に、損失関数であるMSEの計算です。
predictions = A_mesh[..., np.newaxis] * X[:, 0] + B_mesh[..., np.newaxis] * X[:, 1]
によって計算されたpredictionsから、(predictions - y) ** 2)を計算すると以下のようになります(yは計算の際、ブロードキャストされます)。

あるパラメータa, bに対してXから計算された予測値と、実際のyの値の二乗誤差がaxis=2方向に展開されているので、((predictions - y) ** 2).mean(axis=2)によって、axis=2方向に平均を取ることで各パラメータa, bに対するMSEの値が計算されます。
以上のようにして、パラメータの各組み合わせとそれに対する損失関数(MSE)の値が無事に計算されました。

書籍紹介

最近はこの本(見て試してわかる機械学習アルゴリズムの仕組み)を読んで勉強しています。

LightGBMなど新しめのモデルはないですが、機械学習をする上で誰もが通るモデルは一通り網羅されています。

教師なし学習の手法についても書かれているので、初学者の方にはおすすめです。

コメント

タイトルとURLをコピーしました