double descent demystified

# はじめに
今回は機械学習でよく出てくるdouble descentの問題がどのような状況で起きるのかを定性的に理解しようとする[論文](https://arxiv.org/pdf/2303.14151)を紹介しようと思います。この論文では最小二乗法や特異値分解のような工学における基本的な線形代数の知識のみで話が進められてわかりやすいと思います.より詳しい内容についてはこの論文の参考文献になっている論文を読んでください.([例えば](https://arxiv.org/pdf/1903.08560))

特異値分解と最小二乗法についてはすでに資料は豊富にあるのでここでは説明しません.

要約
- 特定の予測を行うタスクにおいて、あるパラメータの数またはある訓練データ数においては学習が致命的に失敗してしまう場合があり、これをdouble descentという
- 線形回帰においてもdouble descentの現象は起こり、パラメータの数で2つに場合分けして定式化することで定性的に理解することができる
- Double Descentは以下の3つの条件のどれかが成り立っていれば起こりづらい
    - 訓練データにほとんど重要じゃない特徴ベクトルがない
    - テストデータが訓練データの重要でない特徴量ベクトルの直交補空間あたりに存在する
    - 観測においてノイズが存在しない


# double descent

double descentを説明するにあたり具体的な例を出しましょう.

$x \in [-1, 1]$の範囲で未知の関数$y$を多項式近似によって予測するタスクを考えます. ここで学習に用いるサンプルデータを$n$個とします. またパラメータ数を$p$の多項式近似とは、$y \sim f(x) = a_{p-1} x^{p-1} + ... + a_0$となるような$a_{p-1},...a_0$を探してくるということです.

![スクリーンショット 2024-12-14 午後6.15.02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3630516/014b334b-cd4e-38d0-3895-c1849d61ab32.png)

以上の図は未知の関数は実は$y = 2 \times x + \cos(50x)$であった場合の結果を表しています. 上の大きい図の横軸はサンプルデータの数/パラメータ数を表していて、縦軸は誤差を表しています. 青線は学習に使っていないデータにおける誤差を表していて、赤線は学習に使ったデータにおける誤差を表しています. 下の小さい図は大きくunderestimated, interpolation threshold, overparametrizedの3つの場合に分けて、それぞれ予測した関数と実際の関数を比べています.

全てのサンプルデータの点$(x,y)$を通るのには必要なパラメータ数は$n$です.つまりパラメータ数が$n$以上であればサンプルデータでの誤差は0にできるということです(上の大きい図の赤線を参照).しかし学習に用いていないデータを使うと本来の関数$y$と$f(x)$は異なるので誤差は0になることはありません. ここでinterpolation thresholdあたりで学習に使っていないデータにおける誤差が発散していることがわかります. このような現象をDouble Descentといいます.[^1]どうしてこのようなことが起きるのでしょうか?

実はdouble descentの現象はニューラルネットワークや線形回帰など様々な最適化手法で観測されることが知られています。今回は線形回帰の場合に限定して定性的な理解をしようと思います.

[^1]: このグラフからはなぜdouble descentという名前がついているのかよくわからないかもしれませんが, ニューラルネットワークなどだとパラメータ数を増やしていくと誤差が減少して増加してまた減少するという挙動を見ることができます.

# 記法の導入

### 線形回帰の教師あり学習の枠組み

学習に用いるデータ数を$N$とする.
データ$\mathcal{D}$を${(\vec{x_n},y_n)}^N_{n=1}$とする. ここで$\vec{x_n} \in \mathbb{R}^d, y_n \in \mathbb{R}$とする.また訓練データ全体を以下のように行列表示する.
```math
X := \begin{pmatrix}
\vec{x_1} \\
\vec{x_2} \\
\vdots \\
\vec{x_n}
\end{pmatrix} \in \mathbb{R}^{N \times D},
Y := \begin{pmatrix}
y_1 \\
y_2 \\
\vdots \\
y_n
\end{pmatrix} \in \mathbb{R}^{N \times 1}
```
教師あり学習の枠組みは以下のようになります.

|目標: $\mathcal{D}$を用いて$f(x) \sim y$となるような写像$f: \mathbb{R}^D \to \mathbb{R}$を探す|
|:-|

今回は線形回帰に限定するので目標を以下のように書き直します.
|目標': $f(\vec{x}) = \vec{x} \cdot \vec{\beta}$とする.この時$\vec{x}\cdot\vec{\beta} \sim y$となるような最適解$\hat{\vec{\beta}} \in \mathbb{R}^D$を探す|
|:-|

また仮定として$rank X = min (N, D)$を付け加えます. この仮定はサンプルされたデータのベクトルがすべて独立していることを意味します.

### Double Descentを理解するための3つのパラメータ
以下の3つの記号を導入します.
- $P$: モデルパラメータの数
- $N$: 訓練データの数
- $D$: データが住んでいると仮定している次元
先ほど分類した3つの場合をこの記号を使って定式化し直します.
- $N < P$のとき, overparametrized(以後overと略記)
- $N > P$のとき, underparametrized(以後underと略記)
- $N = P$のとき, interpolation threshold

この論文ではこれらの記法を使ってdouble descentを定性的に理解することが目標となります

# 線形回帰の場合を数学的に理解
線形回帰の場合は$\vec{x}, \vec{\beta} \in \mathbb{R}^D$より$P$(パラメータ)$=D$(次元)となります. $P$を固定して$N$(訓練データの数)を変化させて挙動の変化を調べます.[^2]

[^2]: ニューラルネットワークなどだとdouble descentはパラメータを増やすことで観測されることが多いのですが、定式化という意味では訓練データの数を変化させた方が考えやすいです.

underの場合とoverの場合で問題設定を少し変えます.
- under (N < P)

よくある最小二乗問題を解きます.

```math
\begin{align}
\hat{\vec{\beta}}_{under} := &argmin_{\vec{\beta}}\frac{1}{N} \sum_{n}||\vec{x_n}\cdot y_n||_{2}^{2} \\
= &argmin_{\vec{\beta}}||X\vec{\beta}-Y||_{2}^{2}
\end{align}
```
この問題の解は$\vec{\beta}_{under} = (X^TX)^{-1}X^TY$となる.[^3]

[^3]: 線形回帰の教師あり学習の枠組みの最後に置いた仮定が$X^TX$の逆行列の存在を保証しています.

- over (P < N)

この場合$X\vec{\beta}=Y$を満たす解は無数に存在するので以下の制約を加えた問題を考えます.
```math
\begin{align}
\vec{\beta}_{over} := argmin_{\vec{\beta}}||\vec{\beta}||_{2}^{2} \\
s.t. \forall n \in [1, N], \vec{x_n} \cdot \vec{\beta} = y_n
\end{align}
```

この制約は最適化で用いられるよくgradient descentは暗に$\|\vec{\beta}\|_{2}^{2}$を最小化している操作を行なっているという点で正当化されます(補足として後述します).

ラグランジュの未定乗数法を用いるとこの問題は解きます.$\mathcal{L}(\vec{\beta}, \vec{\lambda}) := ||\beta||_{2}^2+\vec{\lambda}^T(Y-X\vec{\beta})$とすると

```math
\begin{aligned}
\nabla_{\vec{\beta}} \mathcal{L}(\vec{\beta}, \vec{\lambda})=\overrightarrow{0}=2 \hat{\vec{\beta}}-X^T \vec{\lambda} & \Rightarrow \hat{\vec{\beta}}_{\text {over }}=\frac{1}{2} X^T \vec{\lambda} \\
\nabla_{\vec{\lambda}} \mathcal{L}(\beta, \lambda)=\overrightarrow{0}=Y-X \hat{\vec{\beta}}_{\text {over }} & \Rightarrow Y=\frac{1}{2} X X^T \vec{\lambda} \\
& \Rightarrow \vec{\lambda}=2\left(X X^T\right)^{-1} Y \\
& \Rightarrow \hat{\vec{\beta}}_{\text {over }}=X^T\left(X X^T\right)^{-1} Y
\end{aligned}
```
[^4]

[^4]: 脚注3と同様に線形回帰の教師あり学習の枠組みの最後に置いた仮定が$XX^T$の逆行列の存在を保証しています.

よってまとめると, テストで用いるデータを$\vec{x}_{test}$とすると

```math
\begin{aligned}
\hat{y}_{\text {test }, \text { under }} & =\vec{x}_{\text {test }} \cdot \vec{\beta}_{\text {under }}=\vec{x}_{\text {test }} \cdot\left(X^T X\right)^{-1} X^T Y \\
\hat{y}_{\text {test }, \text { over }} & =\vec{x}_{\text {test }} \cdot \hat{\vec{\beta}}_{\text {over }}=\vec{x}_{\text {test }} \cdot X^T\left(X X^T\right)^{-1} Y
\end{aligned}
```

誤差を最小にするパラメータを $\vec{\beta}^* \in \mathbb{R}^P = \mathbb{R}^D$ とする．またノイズ $e_n$ を加えると $y_n = \vec{x_n} \cdot \vec{\beta}^* + e_n$ となります．行列形式では$Y = X\vec{\beta}^* + E$となります.$(E \in \mathbb{R}^{N \times 1})$

これを用いてモデルが予測する値をもう一度書き直してみます.
- underの場合

```math
\begin{aligned}
\hat{y}_{\text {test }, \text { under }} & =\vec{x}_{\text {test }} \cdot\left(X^T X\right)^{-1} X^T Y \\
& =\vec{x}_{t e s t} \cdot\left(X^T X\right)^{-1} X^T\left(X \beta^*+E\right) \\
& =\vec{x}_{\text {test }} \cdot\left(X^T X\right)^{-1} X^T X \beta^*+\vec{x}_{t e s t} \cdot\left(X^T X\right)^{-1} X^T E \\
& =\vec{x}_{t e s t}^* \cdot \beta^* \cdot\left(X^T X\right)^{-1} X^T E \\
\hat{y}_{\text {test }, \text { under }}-y_{\text {test }}^* & =\vec{x}_{\text {test }} \cdot\left(X^T X\right)^{-1} X^T E
\end{aligned}
```

- overの場合

```math
\begin{aligned}
\hat{y}_{\text {test }, \text { over }} & =\vec{x}_{\text {test }} \cdot X^T\left(X X^T\right)^{-1} Y \\
& =\vec{x}_{\text {test }} \cdot X^T\left(X X^T\right)^{-1}\left(X \beta^*+E\right) \\
& =\vec{x}_{\text {test }} \cdot X^T\left(X X^T\right)^{-1} X \beta^*+\vec{x}_{\text {test }} \cdot X^T\left(X X^T\right)^{-1} E \\
\hat{y}_{\text {test }, \text { over }}-\underbrace{\vec{x}_{\text {test }} \cdot \beta^*}_{\stackrel{\text { def }}{=} y_{\text {test }}^*} & =\vec{x}_{\text {test }} \cdot X^T\left(X X^T\right)^{-1} X \beta^*-\vec{x}_{\text {test }} \cdot I_D \beta^*+\vec{x}_{\text {test }} \cdot\left(X^T X\right)^{-1} X^T E \\
\hat{y}_{\text {test }, \text { over }}-y_{\text {test }}^* & =\vec{x}_{\text {test }} \cdot\left(X^T\left(X X^T\right)^{-1} X-I_D\right) \beta^*+\vec{x}_{\text {test }} \cdot\left(X^T X\right)^{-1} X^T E
\end{aligned}
```

2つを比べてみます. すると同じ項$\vec{x}_{test} \cdot (X^TX)^{-1}X^TE$が出てきていることがわかります．この同じ項はdouble descentがoverの場合でもunderの場合でも発散するということを考慮すると現象を理解する上で重要な項であることがわかります. しかしこの項だけでは意味が分かりずらいので後ほど特異値分解を用いてさらに式を分解することにします. 

# 式の解釈その1

まずoverの場合のみで出てくる項 

$\vec{x}_{\mathrm{test}} \cdot (X^T(XX^T)^{-1}X-I_D) \vec{\beta}^*$ 

について考察します．データ数は $N$ しかないため，モデルはパラメータ空間である $P$ 次元よりも小さい $N$ 次元の情報しか見ることができず，残りの $P-N$ 次元において見えない状態になります．この結果として最適な線形関係 $\vec{\beta}^*$ に関する情報の一部が失われてしまうことを意味しています．

# 式の解釈2
まず$X$について特異値分解を行い, $X = U\sigma V^T$とします. 仮定より$N = rank(X)$より$X$の特異値を大きい順に並べて$\sigma_1 \geq \cdots \sigma_n \geq 0$とします. ここで
```math
\begin{align}
(X^TX)^{-1}X^T &= (V\Sigma^TU^TU\Sigma V^T)^{-1}V\Sigma^TU^T \\
&= (V^T)^{-1}(\Sigma^T\Sigma)^{-1}V^{-1}V\Sigma^TU^T \\
&= V(\Sigma^T\Sigma)^{-1}\Sigma^TU^T
\end{align}
```
成分ごとに分けて書くと以下のようになります.
```math
\begin{align}
\hat{y}_{test,under}-y_{test}^* &= \vec{x_{test}}\cdot V\Sigma^*V^TE
&= \sum_{r=1}^N \frac{1}{\sigma_r}(\vec{x}_{test}\cdot \vec{v_r})(\vec{u_r}\cdot E)
\end{align}
```
この項を解釈していきます.

(1)
$\frac{1}{\sigma_r}$
訓練データの特徴量の重要度の逆数

(2)
$\vec{x}_{test} \cdot \vec{v_r}$
訓練データの特徴ベクトルと入力ベクトルの内積

(1)(2)を合わせて考えます
(1)によるとdouble descentは重要度が小さい特徴量により左右されることがわかります. またテストデータに重要じゃない特徴ベクトルの要素がより含まれているとよりdouble descentでの汎化誤差の増大につながることがわかります.

(3)
$\vec{u_r} \cdot E$
訓練データの特徴ベクトルとノイズの内積

(3)の意味は特徴ベクトルとノイズの区別がつきにくいとき, 汎化誤差の増大につながることがわかります.

# まとめ

これらの解釈を元にするとDouble Descentは以下の3つの条件のどれかが成り立っていれば起こりづらいことがわかります.
- 訓練データにほとんど重要じゃない特徴ベクトルがない
- テストデータが訓練データの重要でない特徴量ベクトルとはあまり関係しない.
- 観測においてノイズが存在しない

# 補足 どうしてgradient descentは暗に正則化が行われるのか

$y \in \mathbb{R}^n$を予測するモデル$X\omega$($\omega$はパラメータ)を持っているとし,
$L(\omega) = ||X \omega - y||^2_2=||e||^2_2$を最小化することを考えます. ここで$X \in \mathbb{R}^{n \times d_1}$とします.($n$はサンプル数, $d_1$はデータが住んでいると仮定している次元)

勾配学習のルールにより
```math
\begin{align}
\omega(t+1) &= \omega(t)-\eta X^T e(t) \\
\Rightarrow X \omega (t+1) - y &= (X \omega (t) - y) - \eta X X^T e(t) \\
e(t+1) &= e(t) - \eta XX^Te(t) \\
&= (I-\eta XX^T)e(t)
\end{align}
```

gradient descentを連続化するとパラメータの更新は以下のようになります.
```math 
\begin{align}
\dot{\omega} &= -\eta X^Te \\
\omega(0) &= 0
\end{align}
```

ここで$\dot{\omega}$は$X^T$の一次結合でかけることから$\omega$は$X^T$の一次結合で書けます.
gradient descentの収束先の$\omega$を$\omega^{*}$

とし$\omega^{*}=X^T\alpha$であるとします. overparametrizedの場合, $X\omega$は常に$y
となります.よって
```math
\begin{align}
y = X\omega^* = XX^T\alpha \\
\Rightarrow\alpha = (XX^T)^{-1}y
\Rightarrow\omega^* = X^T(XX^T)^{-1}y
\end{align}
```
この解は$L(\omega)$の[最小ノルム解](https://manabitimes.jp/math/2747)そのものになります.

# 参考文献
- [Double Descent Demystified](https://arxiv.org/pdf/2303.14151v1)
- [最小二乗法と特異値分解について詳しい講義ノート](https://webpark2049.sakura.ne.jp/shmdlb/wp-content/uploads/2024/10/MMVIIInotes2024.pdf)
- [最小ノルム解](https://manabitimes.jp/math/2747)