正則化パラメータλをlogスケールで最適化する理由 - ProphetをOptunaで最適化する際に出た疑問 -

結論

正則化パラメータ(Prophetのprior_scale系のパラメータなど)はlogスケールでチューニングしよう！

はじめに

ProphetはFacebookが公開しているオープンソースの時系列予測ライブラリであり、シンプルな設計と軽量さ、また実用的な精度を兼ね備えたツールとして広く使われています。
https://facebook.github.io/prophet/

そんな便利なProphetですが、こちらの記事で示されているように、ハイパーパラメータチューニングを行うことでデフォルトよりも大きく性能改善することが多く、Prophetを活用する場合はハイパーパラメータチューニングが必須と言えます。

Prophetの公式ドキュメントには、ハイパーパラメータチューニング可能なパラメータとして以下のパラメータが記載されています。

• changepoint_prior_scale：トレンドの変化点の影響度を調整する正則化パラメータ
• seasonality_prior_scale：季節性の影響度を調整する正則化パラメータ
• holidays_prior_scale：休日効果の影響度を調整する正則化パラメータ
• seasonality_mode：季節性を加法(additive)的に取り扱うか、乗法(multiplicative)的に取り扱うかを調整するパラメータ

4つのうち上3つのパラメータはどれもprior_scaleであり、例えば、changepoint_prior_scaleの説明には、

Parameters like this (regularization penalties; this is effectively a lasso penalty) are often tuned on a log scale.

(訳)この種のパラメータ（正則化項。実質的にはlassoペナルティ）は、対数スケールで調整されることがよくあります。

と書かれています。

これを考慮してOptunaの最適化関数を書くとすると

import optuna

def objective(trial):
    # ...(略)...
    changepoint_prior_scale = trial.suggest_float(
        "changepoint_prior_scale", 0.001, 0.5, log=True
    )
    # ...(略)...

となるわけですが、デフォルトではlog=Falseであり、何も考えないと線形スケールでの最適化をすることになります。

正則化のパラメータだからlogスケールで最適化すると言えばそれまでですが、なんでそうしなきゃいけないのかがパッと説明できなかったので、色々と実験しながら整理し、自分なりに納得した内容を共有いたします。

事前分布と正則化の関係性

結論、prior_scale は正則化に関係するパラメータなのですが、これを理解するには、パラメータに事前分布を導入したベイズ推定の式と正則化項を導入した目的関数との数学的な対応を確認する必要があります。
以下、線形回帰モデルを前提として、細かい数式の展開を省いて説明します。

モデルの前提

以下のように変数を定義します。

• 説明変数：$\mathbf{X} \in \R^{n \times d}$
• 目的変数：$\mathbf{y} \in \R^n$
• モデルパラメータ：$\boldsymbol{\theta} \in \R^{d}$
• 観測ノイズ：$\boldsymbol{\epsilon} \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$
• $n$：サンプル数（データの件数）
• $d$：特徴量の次元数（説明変数の数）

なお、正規分布の確率密度関数は以下のように表されます(単純のため1次元で記載)。

$$\mathcal{N}(\mu, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}}\exp \bigg(-\frac{(x-\mu)^2}{2 \sigma^2} \bigg)$$

上の変数を使い、線形回帰モデルは、以下のように定義します。

$$\mathbf{y} = \mathbf{X}\boldsymbol{\theta} + \boldsymbol{\epsilon}$$

また、パラメータ$\boldsymbol{\theta}$に対しては、以下のような事前分布(prior)を仮定します。

$$p(\boldsymbol{\theta}) = \mathcal{N}(\mathbf{0}, \tau^2\mathbf{I})$$

この$\tau$がprior_scaleとなります。

MAP推定

事前分布を仮定した際は、データが得られた後のパラメータの事後確率を最大化するMAP(Maximum A Posteriori)推定を用いてパラメータを推定します。
ベイズの定理を使うと

$$\begin{aligned} \hat{\boldsymbol{\theta}}_{MAP} &= {\rm arg}\max_{\boldsymbol{\theta}} p(\boldsymbol{\theta} | \mathbf{X},\mathbf{y}) \\ &= {\rm arg}\max_{\boldsymbol{\theta}} p(\mathbf{y} | \mathbf{X} ; \boldsymbol{\theta}) p(\boldsymbol{\theta}) \end{aligned}$$

対数をとっても最適解は変わらないため、対数を取ると以下のように変形できます。

$$\hat{\boldsymbol{\theta}}_{MAP} = {\rm arg} \max_{\boldsymbol{\theta}} \left[ \log p(\mathbf{y} | \mathbf{X} ; \boldsymbol{\theta}) + \log p(\boldsymbol{\theta}) \right]$$

また、各項を計算すると、第1項目が、

$$\log p(\mathbf{y} | \mathbf{X} ; \boldsymbol{\theta}) = -\frac{1}{2} \| \mathbf{y} - \mathbf{X} \boldsymbol{\theta} \|^2 + \text{const}$$

第2項目が

$$\log p(\boldsymbol{\theta}) = -\frac{1}{2\tau^2} \| \boldsymbol{\theta} \|^2 + \text{const}$$

となり、MAP推定の最適化の式は以下のように書き換えられます。

$$\hat{\boldsymbol{\theta}}_{MAP} = {\rm arg} \max_{\boldsymbol{\theta}} \left[ - \| \mathbf{y} - \mathbf{X} \boldsymbol{\theta} \|^2 -\frac{1}{\tau^2} \| \boldsymbol{\theta} \|^2 \right]$$

定数倍は最適解に影響を与えないため落としています。
おや、この式に似た式をどこかで見たことありますね。

Ridge回帰の目的関数

Ridge回帰の目的関数は以下で定義されます。

$$\hat{\boldsymbol{\theta}}_{\text{Ridge}} = \arg\min_{\boldsymbol{\theta}} \left( \| \mathbf{y} - \mathbf{X} \boldsymbol{\theta} \|^2 + \lambda \| \boldsymbol{\theta} \|^2 \right)$$

ここで $\lambda$ は正則化係数であり、モデルの複雑さ（パラメータの大きさ）に対してどの程度のペナルティを与えるかを調整します。

先ほど導出したMAP推定の式において、最適化の方向を最大化から最小化に変えて

$$\lambda = \frac{1}{\tau^2}$$

とおくと、Ridge回帰は、平均0・分散$\tau^2$の正規分布をパラメータの事前分布と仮定したときのMAP推定と数式的に等価になります。

ここで、prior_scaleが$\tau$であるわけなので、

$$\lambda = \frac{1}{\texttt{prior\_scale}^2}$$

という関係が成り立ちます。

したがって、prior_scaleを小さく設定するほど、$\lambda$は大きくなり、モデルのパラメータはゼロ方向へ強く引き寄せられます（強い正則化）。
逆に、prior_scaleを大きくすると、正則化は弱まり、より自由なパラメータが許されるようになります。

正則化パラメータのスケールとモデルパラメータの変化の関係性

prior_scaleが正則化に関係するパラメータであることがわかったところで、なぜ正則化のパラメータはlogスケールで調整する必要があるのかを確認します。
この項以降は簡単のためprior_scaleではなく正則化パラメータ$\lambda$を使って説明します。

正則化の役割(簡易的な説明)

よくある絵を使ってRidge回帰の例に正則化の役割について説明します。

Ridge回帰の正則化のイメージ(キカガクより引用)
この図は、損失関数の等高線（楕円）と、正則化項による制約領域（円）の関係を示しています。

• 楕円：誤差関数（例えば $\|\mathbf{y} - \mathbf{X} \boldsymbol{\theta}\|^2$）の等高線。中心に近いほど誤差が小さい
• 円：正則化項（$\|\boldsymbol{\theta}\|^2$)により作られるモデルの複雑さに対する制限

Ridge回帰では、これらの2つが交わる点が最適解となります。つまり、「誤差を最小にしたい」という力学と「パラメータの大きさを抑えたい」という力学のバランスをとった解がRidge回帰では求められるというわけです。

このバランスを調整するのが正則化のパラメータ $\lambda$ であり、$\lambda$の値の大きさによってこの円の大きさが変化します。

モデルパラメータは正則化パラメータによってどう変化するのか(線形回帰の場合)

次に、正則化を加味した目的関数を最適化した場合に、$\lambda$がモデルパラメータ$\theta$にどう変化をもたらすのかを非常に簡易なモデルを使って数式で確認してみます。

ここでは説明変数$x$と目的変数$y$が1次元の単回帰モデルを定義し、観測データ$(x, y)$が1点のみあると仮定して、回帰係数$\theta$を推定します。

このとき、Ridge回帰の目的関数は次のようになります。

$$\mathcal{L}(\theta) = (y - \theta x)^2 + \lambda \theta^2$$

この目的関数を$\theta$について最小化する(微分して0となる$\theta$を求める)と、解析的に以下のような最適解が得られます。

$$\hat{\theta} = \frac{xy}{x^2 + \lambda}$$

この式から$\lambda$と$\theta$の関係性は非線形だということがわかります。
非常に簡単なモデルにおいても正則化のパラメータはモデルのパラメータに対して非線形に影響を及ぼすため、高次元のモデルでも同様に非線形に影響を及ぼすだろうことが推察されます。

この結果から、正則化のパラメータを変化させたときにモデルのパラメータが非線形に変化するならば、正則化のパラメータを線形に調整するとよくないのではないかという仮説が立てられます。

正則化パラメータをlogスケール、線形スケールで探索させた時の探索結果の比較

では実際に正則化パラメータをlogスケール、線形スケールで探索した際に探索結果がどう変わるのかみてみましょう。

今回は非常にシンプルな設定であるため、モデルの評価結果との関係性については考慮していません。本可視化の主目的は、$\lambda$ のスケールの選び方によって、最適パラメータの探索がどのように変化するかを確認することにあります。

条件設計

今回は簡単のため、2つのモデルパラメータ$w_1$, $w_2$がある線形回帰モデルで検証しました。
2つのモデルパラメータに対し、真のパラメータを所与とした上でデータをノイズを加えて100サンプル作成し、そのデータに対してRidge回帰を行います。
また、正則化パラメータ $\lambda$を logspace、linspace のそれぞれで等間隔に100点サンプリングします。

import numpy as np

# data creation
np_rng = np.random.RandomState(42)
n_samples = 100
w_true = np.array([10.4, 6.4])
X = np_rng.randn(n_samples, 2)
noise = np_rng.normal(0, 1, size=n_samples)
y = X @ w_true + noise

# lambda space
lambdas_log = np.logspace(-3, 3, 100)
lambdas_lin = np.linspace(0.001, 1000, 100)

パラメータ最適化と可視化

まずは、最小二乗法でパラメータを求めてみます。

from sklearn.linear_model import LinearRegression

# Simple Linear Regression check
model_ols = LinearRegression(fit_intercept=False)
model_ols.fit(X, y)
w_ols = model_ols.coef_
w_ols

array([10.57610297,  6.23139683])

データにノイズを与えているので最適解とは少しずれていることがわかります。

次に、上で定義した$\lambda$に対して、Ridge回帰でパラメータ推定を行い、探索点と最適パラメータを保存します。

import numpy as np
from sklearn.linear_model import Ridge

# weight calculation loop
def solver_loop(lambdas, X, y):
  w_list = []
  for l in lambdas:
    ridge = Ridge(alpha=l, fit_intercept=False)
    ridge.fit(X, y)
    w_list.append(ridge.coef_)
  return np.array(w_list)

w_ridge_log = solver_loop(lambdas_log, X, y)
w_ridge_lin = solver_loop(lambdas_lin, X, y)

最後に、$\lambda$の探索点と最適パラメータの変化を可視化します。

import numpy as np
import matplotlib.pyplot as plt

# visualize parameters while changing lambda (logspace vs linspace)
fig, ax = plt.subplots(1, 2, figsize=(16, 6), gridspec_kw={'width_ratios': [1, 1], 'wspace': 0.25})

# logspace
sc1 = ax[0].scatter(w_ridge_log[:, 0], w_ridge_log[:, 1],
                    c=np.log10(lambdas_log))
ax[0].plot(w_ridge_log[:, 0], w_ridge_log[:, 1], color='gray', linewidth=1, alpha=0.5)
ax[0].plot(w_true[0], w_true[1], 'ro', label='True parameters', markersize=8)
ax[0].plot(w_ols[0], w_ols[1], 'kx', label='OLS solution', markersize=8)
ax[0].set_title("Log-spaced λ (1e-3 to 1e3)")
ax[0].set_xlabel(r"$w_1$")
ax[0].set_ylabel(r"$w_2$")
ax[0].grid(True)
ax[0].legend()

# linspace
sc2 = ax[1].scatter(w_ridge_lin[:, 0], w_ridge_lin[:, 1],
                    c=np.log10(lambdas_lin))
ax[1].plot(w_ridge_lin[:, 0], w_ridge_lin[:, 1], color='gray', linewidth=1, alpha=0.5)
ax[1].plot(w_true[0], w_true[1], 'ro', label='True parameters', markersize=8)
ax[1].plot(w_ols[0], w_ols[1], 'kx', label='OLS solution', markersize=8)
ax[1].set_title("Linearly-spaced λ (0.001 to 1000)")
ax[1].set_xlabel(r"$w_1$")
ax[0].set_ylabel(r"$w_2$")
ax[1].grid(True)
ax[1].legend()

# color bar
fig.subplots_adjust(right=0.80)  
cbar_ax = fig.add_axes([0.83, 0.15, 0.02, 0.7]) 
cbar = fig.colorbar(sc2, cax=cbar_ax)
cbar.set_label(r'$\log_{10}(\lambda)$')

plt.suptitle("Comparison of Ridge Coefficient Trajectories (Log vs Linear λ)", y=1.02)
plt.show()

左側がlogスケールで正則化パラメータ$\lambda$を変化させた際の探索点とパラメータの値、右側が線形スケールです。
logスケールでは、$\lambda$の小さい領域に多くのサンプル点が割り当てられており、解がどのように真のパラメータから縮小されるかを解像度高く確認することができています。
一般的に$\lambda$が大きすぎるとモデルとしての複雑性が失われすぎて精度が下がる傾向にあるため、最適なパラメータ付近で多く探索点を持てていることで、精度と正則化のバランスを取った良いモデルを見つけ出すことができるのではないかと考えられます。
一方で、線形スケールの方は、最適なパラメータ付近での探索点が少なく、$\lambda$が大きいところで多くの探索点を持つ挙動となっています。
誤差関数の最適解から遠い部分に対して細かく調べることが最適なモデルパラメータの探索に対して効果的とは言えず、パラメータ探索の方策としては適切ではないと考えられます。

この結果から、正則化パラメータの最適化においてモデルの挙動変化を満遍なく確認するためには、logスケールでの探索が適切であることがわかりました。

Optunaなどの自動最適化における注意点

Optunaのようなベイズ最適化ライブラリを用いると、探索空間の形に合わせて効率的にパラメータをサンプリングしてくれるため、ランダムサーチやグリッドサーチよりも賢く探索してくれることが期待されます。
ですが、Optunaのベースサンプリングはスケールに合わせて一様分布でサンプリングするため、logスケールの指定を明示しなければベイズ最適化で適切な箇所を見つけるまでの初期のサンプリングでは線形スケールでの探索となり最適な領域を見逃す可能性があります。
そのため、正則化のような非線形影響をもつパラメータに対してはlog=Trueを明示することが重要だと言えそうです。

まとめ

Prophetにおける prior_scale 系パラメータが正則化に関連するパラメータであり、logスケールでチューニングすべきであることを解析的、実験的に理由づけしてみました。
生成AIがコーディングもしてくれる便利な世の中になりましたが、生成AIが出すコードが適切に記載されているかを確認するために、こういった基礎をしっかり身につけていきたいなと思いました。

ご覧いただきありがとうございました！

参考文献

Prophet Hyperparameter Tuning
https://facebook.github.io/prophet/docs/diagnostics.html
https://qiita.com/dcm_shimadas/items/8f049c4a258f0fbbfe23

Optuna
https://optuna.readthedocs.io/en/stable/reference/generated/optuna.trial.Trial.html#optuna.trial.Trial.suggest_float

ベイズ推定
https://zenn.dev/misya11p/articles/769a6cc2932262