Orthogonal Matching Pursuit (OMP)

最終更新 2020-07-01 読了時間 4 分

まとめ

Orthogonal Matching Pursuit (OMP) は、残差と最も相関の高い特徴量を 1 本ずつ選び、選んだ列だけで最小二乗法をやり直す貪欲法です。ラッソのように罰則の強さ \(\alpha\) を調整するのではなく、n_nonzero_coefs で「残す特徴量の本数」を直接指定するので、選ばれなかった列の係数は厳密に 0 になります。このページでは 40 本の特徴量のうち 4 本だけが効く人工データを作り、OrthogonalMatchingPursuit(n_nonzero_coefs=4) と Lasso(alpha=0.05) の係数を真の係数と並べて比べます。読み終えると、OMP がどの順番で列を拾い、なぜラッソの係数は小さな値が残るのかを説明できるようになります。

先にリッジ回帰とラッソ回帰で L1 罰則が係数を 0 にする仕組みを見ておくと、OMP が「罰則ではなく本数で疎にする」別の道だと分かります。

直感 #

40 本の特徴量があって、本当に効いているのは 4 本だけだと分かっているとします。OMP はまず \(y\) そのものを残差として、40 本の中で残差と最も相関の高い列を 1 本選び、その 1 本だけで最小二乗法を解いて残差を更新します。残った 39 本から再び相関最大の列を選びます。これを n_nonzero_coefs=4 回繰り返すと 4 本が選ばれ、残り 36 本の係数は一度も触られないので 0 のままです。「選ぶ → 解き直す → 残差を更新する」を 1 本ずつ進めるだけなので、各ステップで何が起きたかを追いやすいのが特徴です。

式で書くと各ステップで何をしているか #

初期残差を \(\mathbf{r}^{(0)} = \mathbf{y}\) とし、各ステップ \(t\) で次の 4 つの操作を繰り返します。まず、各特徴量 \(\mathbf{x}_j\) と残差 \(\mathbf{r}^{(t-1)}\) の内積を計算し、絶対値が最大の特徴量 \(j\) を選びます。この内積が「残差を最も説明できる列」の指標です。次に、選択済み特徴量集合 \(\mathcal{A}_t\) に \(j\) を追加します。続いて \(\mathcal{A}t\) の特徴量だけを使って最小二乗解 \(\hat{\boldsymbol\beta}{\mathcal{A}t}\) を求め直します。この「解き直し」が Orthogonal の由来で、新しい残差は選択済みの列すべてと直交します。最後に新しい残差 \(\mathbf{r}^{(t)} = \mathbf{y} - \mathbf{X}{\mathcal{A}t} \hat{\boldsymbol\beta}{\mathcal{A}_t}\) を計算します。

残差 r = y 残差と相関が最大の列 j を選ぶ選択済み集合に j を追加選択済みの列だけで最小二乗を解く残差を更新本数が n_nonzero_coefs に達した? 係数を出力

n_nonzero_coefs で指定した本数に達するか、tol で指定した残差の大きさを下回るまでこの操作を続けます。

Pythonを用いた実験や説明 #

次のコードは、200 サンプル × 40 特徴量のうち 4 本だけに大きさ 3 程度の係数を置いた人工データを作り、OrthogonalMatchingPursuit(n_nonzero_coefs=4) と Lasso(alpha=0.05) を学習させます。真の係数・OMP・ラッソの 3 本の棒を特徴量ごとに並べます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
from __future__ import annotations

import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import Lasso, OrthogonalMatchingPursuit
from sklearn.metrics import mean_squared_error

def run_omp_vs_lasso(
    n_samples: int = 200,
    n_features: int = 40,
    sparsity: int = 4,
    noise_scale: float = 0.5,
    xlabel: str = "feature index",
    ylabel: str = "coefficient",
    label_true: str = "true",
    label_omp: str = "OMP",
    label_lasso: str = "Lasso",
    title: str | None = None,
) -> dict[str, object]:
    """Compare OMP and lasso on synthetic sparse regression data.

    Args:
        n_samples: Number of training samples to generate.
        n_features: Total number of features in the dictionary.
        sparsity: Count of non-zero coefficients in the ground truth.
        noise_scale: Standard deviation of Gaussian noise added to targets.
        xlabel: Label for the coefficient plot x-axis.
        ylabel: Label for the coefficient plot y-axis.
        label_true: Legend label for the ground-truth bars.
        label_omp: Legend label for the OMP bars.
        label_lasso: Legend label for the lasso bars.
        title: Optional title for the bar chart.

    Returns:
        Dictionary containing recovered supports and MSE values.
    """
    rng = np.random.default_rng(0)

    X = rng.normal(size=(n_samples, n_features))
    true_coef = np.zeros(n_features)
    true_support = rng.choice(n_features, size=sparsity, replace=False)
    true_coef[true_support] = rng.normal(loc=0.0, scale=3.0, size=sparsity)
    y = X @ true_coef + rng.normal(scale=noise_scale, size=n_samples)

    omp = OrthogonalMatchingPursuit(n_nonzero_coefs=sparsity)
    omp.fit(X, y)
    lasso = Lasso(alpha=0.05)
    lasso.fit(X, y)

    omp_pred = omp.predict(X)
    lasso_pred = lasso.predict(X)

    fig, ax = plt.subplots(figsize=(10, 5))
    indices = np.arange(n_features)
    ax.bar(indices - 0.3, true_coef, width=0.2, label=label_true, color="#2ca02c")
    ax.bar(indices, omp.coef_, width=0.2, label=label_omp, color="#1f77b4")
    ax.bar(indices + 0.3, lasso.coef_, width=0.2, label=label_lasso, color="#d62728")
    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel)
    if title:
        ax.set_title(title)
    ax.legend()
    fig.tight_layout()
    plt.show()

    return {
        "true_support": np.flatnonzero(true_coef),
        "omp_support": np.flatnonzero(omp.coef_),
        "lasso_support": np.flatnonzero(np.abs(lasso.coef_) > 1e-6),
        "omp_mse": float(mean_squared_error(y, omp_pred)),
        "lasso_mse": float(mean_squared_error(y, lasso_pred)),
    }

metrics = run_omp_vs_lasso(
    xlabel="特徴量インデックス",
    ylabel="係数",
    label_true="真の係数",
    label_omp="OMP",
    label_lasso="Lasso",
    title="OMP と Lasso の係数比較",
)
print("真のサポート:", metrics['true_support'])
print("OMP のサポート:", metrics['omp_support'])
print("Lasso のサポート:", metrics['lasso_support'])
print(f"OMP の MSE: {metrics['omp_mse']:.4f}")
print(f"Lasso の MSE: {metrics['lasso_mse']:.4f}")

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

39-43 行目40 本のうち 4 本だけが効くデータを作る
true_coef を全部 0 で初期化し、rng.choice で選んだ 4 本（sparsity）だけに標準偏差 3.0 の係数を入れています。これが図の緑の棒です。X は 40 本が互いに無相関な標準正規乱数なので、OMP にとっては選びやすい条件です。noise_scale=0.5 を 3.0 くらいまで上げると、OMP が間違った列を拾い始める様子を見られます。
45-46 行目残す本数を直接指定する
n_nonzero_coefs=sparsity が OMP の唯一の主要パラメータで、ここでは真の本数 4 をそのまま渡しています。実データでは真の本数は分からないので、OrthogonalMatchingPursuitCV に交差検証で選ばせます。本数を 6 にすると、余った 2 本はノイズに最も相関した列が拾われます。
47-48 行目比較相手のラッソは罰則で疎にする
Lasso(alpha=0.05) は L1 罰則の強さで係数を 0 に落とすので、「何本残るか」は \(\alpha\) の結果として決まります。\(\alpha = 0.05\) は小さめの値で、4 本以外にも絶対値の小さな係数が残ります。この違いが赤い棒と青い棒の差になります。
55-57 行目3 本の棒をずらして並べる
同じ特徴量インデックスに対して、緑（真の係数）を 0.3 左、赤（ラッソ）を 0.3 右にずらして描いています。幅 0.2 の棒を 3 本並べるためのずらし幅で、重ねて描くと青が緑を隠してしまいます。4 本の位置で 3 色の高さがそろっていれば、両手法とも正しい列を見つけたと読めます。
66-72 行目ラッソだけ 1e-6 の閾値で数える
omp_support は np.flatnonzero(omp.coef_) で足りますが、ラッソの係数は厳密な 0 ではなく 1e-10 程度の値が残ることがあるため、np.abs(lasso.coef_) > 1e-6 で数えています。OMP は選ばなかった列に一度も値を入れないので、閾値なしで厳密に 4 本になります。

実行結果の読み方 #

図の 4 か所で緑・青・赤の棒がそろって立っていれば、両手法が正しいサポートを復元しています。青（OMP）は選ばれた 4 本以外が厳密に 0 で、omp_support の出力は真のサポートと同じ 4 つの番号になります。赤（ラッソ）は 4 本以外にも小さな棒が残り、lasso_support は 4 より多い番号を含みます。ラッソの 4 本は L1 罰則で真の値より縮んで低く描かれますが、OMP は選んだ列で最小二乗を解き直すので緑とほぼ同じ高さです。

使いどころと落とし穴 #

「残す本数」を仕様として決めたい場面、たとえばセンサーを 5 個だけ選ぶといった制約がある場合に OMP が向きます。落とし穴は特徴量どうしの相関です。相関の強い列が 2 本あると、最初の 1 本を選んだ時点で残差からもう 1 本の情報が消え、2 本目が選ばれないことがあります。この状況では Elastic Net の方が安定します。もう 1 つは尺度で、列の分散が違うと内積の大きさがスケールに左右されるため、StandardScaler で揃えてから fit します。

参考文献 #

Pati, Y. C., Rezaiifar, R., & Krishnaprasad, P. S. (1993). Orthogonal Matching Pursuit: Recursive Function Approximation with Applications to Wavelet Decomposition. In Conference Record of the Twenty-Seventh Asilomar Conference on Signals, Systems and Computers.
Tropp, J. A. (2004). Greed is Good: Algorithmic Results for Sparse Approximation. IEEE Transactions on Information Theory, 50(10), 2231–2242.

リッジ・ラッソ回帰 — 正則化による変数選択
Elastic Net — L1+L2正則化
特徴選択の基本 — 特徴選択の体系的解説

係数パスと特徴選択 #

スライダーで正則化係数 α を大きくすると、Lasso の係数が 1 本ずつ 0 に落ちていく順番を確認できます。OMP の n_nonzero_coefs は、この「残す本数」を α ではなく整数で直接指定するパラメータです。