二重降下現象

最近，甘利さんの統計神経力学の本を研究室内で輪読していたが，自然勾配法の記事に引き続き，自分が二重降下の章を担当することになった．せっかく資料を作ったので一部をブログに載せておく．

ちなみに，以前にもグラフにするだけしてみていたりはする．深層学習まわりの理論の話では個人的にはこのトピックが最も好き．

potaxyz.hatenablog.jp

導入

二重降下現象（Double Descent）とは，統計学や機械学習において，モデルのパラメータ数を増やしていくにつれて一度は複雑性誤差が増加するが，その後再び誤差が減少していくという現象である．古典的な統計モデルでは、複雑性誤差が一度増加する手前でバイアスとバリアンスのトレードオフをとるようにモデルを調整することが一般的だったが，現在の深層学習に代表されるような過剰パラメータモデルではパラメータを増やせば増やすほど誤差を減らすことができる場合があるというお話．古典的な学習理論の話から始めて深層学習での二重降下現象の話をする．

統計的学習理論と情報量基準

ここでは，すこし遠くから始めるが統計的学習理論の枠組みに出てくる量をいくつか紹介しつつ，汎化誤差を訓練誤差を使った推定量で推定することを考える（というのも本がそういう流れっぽかったので）．

機械学習の目標は，観測データから適切な仮説（予測モデル） を学習するアルゴリズムを設計することである．

• 仮説: 入力空間 $\mathcal{X}$ から出力集合 $\mathcal{Y}$ への関数 $$ h:\mathcal{X}\to\mathcal{Y} $$

学習とは，真の入出力関係を精度よく近似できるように、データを使って $h$ を調整する過程のこと．

$h$ による予測のよさは損失関数で測る．

• 損失関数: 出力値と予測結果の間の誤差を測る関数 $$ l:\mathcal{Y}\times\mathcal{Y}\to[0,\infty) $$

どのように学習するか？といえば， 期待リスクをできるだけ小さくするように学習したい．

• 期待リスク（汎化誤差）: テストデータの分布における予測値の損失の期待値 $$ L(h)=E _ {(X,Y)}[l(Y,h(X))] $$

ただし，実際には $(X,Y)$ の分布はわからないので，観測データによって $L(h)$ を代替した指標として経験リスクを使う．

• 経験リスク（訓練誤差）: 観測データ $(X _ 1,Y _ 1),\ldots,(X _ n,Y _ n)$ に対する予測値 $h(X _ i)$ と観測値 $Y _ i$ の損失の平均 $$ L _ n(h)=\frac{1}{n}\sum _ {i=1} ^ nl(Y _ i,h(X _ i)) $$

汎化性能の良さは真に達成可能な期待リスク（ベイズリスク）との差で決まり，これとの差（残余誤差）が小さいほどよい．

• ベイズリスク: 損失関数 $l$ を定めたときに任意の可測関数 $h$ のもとで達成できる期待リスクの下限 $$ L ^ \ast:= L(h ^ \ast)=\inf _ {h:\mathbb{R} ^ d\to\mathbb{R}} L(h) $$
• ベイズ規則: ベイズリスクを達成する仮説 $$ h ^ *=\arg\inf _ {h:\mathbb{R} ^ d\to\mathbb{R}} L(h) $$

学習データを使って，できるだけ $h ^ \ast$ に近い $h$ をつくりたい（ちなみに，学習データ $\mathcal{D} _ n$ とすると， $E\ _ {\mathcal{D}\ _ n}[L(h _ n)]\to L ^ *,(n\to\infty)$ になる一致性をもつことが望ましい）．

ここでは $\theta$ でパラメトライズされた $h _ \theta$ を以下では考えることにする．ベイズリスクを達成するパラメータを $\theta ^ *$，経験リスク最小化で達成されるパラメータを $\hat{\theta}$ とかくと，ここまで紹介してきた量は

• 大数の法則から $L _ n(h _ {\theta ^ \ast})\to L(h _ {\theta ^ \ast})\,(n\to\infty,\mathrm{a.s.})$
• ベイズリスクと汎化誤差の関係から $L(h _ {\theta ^ \ast})\leq L(h _ {\hat{\theta}})$
• 経験リスク最小化で $\hat{\theta}$ を得ていたので $L _ n(h _ {\hat{\theta}})\leq L _ n(h _ {\theta ^ *})$

の関係がある．これらをまとめると，十分大きな $n$ で

$$ L _ n(h _ {\hat{\theta}})\leq L _ n(h _ {\theta ^ \ast})\approx L(h _ {\theta ^ \ast})\leq L(h _ {\hat{\theta}}) $$

が成り立つ．本当に実現したいのは $\theta ^ \ast$ だが，実際に実現できるのは $\hat{\theta}$なので，ここでは $L _ n(h _ {\hat{\theta}})$ で $L(h _ {\hat{\theta}})$ を推定することを考える．ここで $E[L _ n(h _ {\hat{\theta}})-L(h _ {\hat{\theta}})]=O(n ^ {-1})$ が示せるため， $n\to\infty$ では $L _ n(h _ {\hat{\theta}})$ は $L(h _ {\hat{\theta}})$ の漸近不偏推定量になっていることがいえる．しかし上の不等式を見れば明らかに過小評価しており，有限の $n$ ではこのオーダーの誤差が無視できないことがある．そこで，その分補正した推定量 $\hat{R} _ n$ として $E[\hat{R} _ n-L(h _ {\hat{\theta}})]=o(n ^ {-1})$ をみたすものを考え，この $\hat{R} _ n$ によって $L(h _ {\hat\theta})$ を推定することを考える（これは理論的にはリスク $R:= E[L(h _ {\hat{\theta}})]=E _ {\mathcal{D}}[L(h _ {\hat{\theta}(\mathcal{D})})]$ の推定になる．訓練データによらず $\hat{\theta}$ の作るプロセスが平均的に良いかを示す指標とみれる．）

ちなみに，この考え方を使っているのが情報量基準の導出になっている．この導出を追っていくと経験リスクと期待リスクの間の関係式は

$$ L(h _ {\hat\theta})\approx L _ n(h _ {\hat\theta})+\frac{p}{n} $$

となっていることがわかる．この式から

• パラメータ数 $p$ を増やすと訓練誤差 $L _ n$ は減少し， $n=p$ で $0$ となる．
• パラメータ数 $p$ を増やすと汎化誤差 $L$ は途中までは減少するが，途中から増加に転じ，訓練データに過剰適応する．結果として新たなデータに適応する汎化能力を失う（過学習）．

ことがわかる．したがって，この関係式で右辺が一番小さくなるモデルが，最もデータをよく説明するモデルになりそうだといえる（情報量規準を使ったモデル選択）．

図はReconciling modern machine-learning practice andthe classical bias–variance trade-offより

二重降下現象

実は，ここまでに紹介した描像は $p\leq n$ に限った話（古典的レジーム）だったことが最近わかってきた．近年のニューラルネットワークを使ったパラメータ設定では，多くの場合 $p\gt n$ を採用しており，この場合には異なる描像が出現する．イメージとしては次の図のような感じである．

図はReconciling modern machine-learning practice andthe classical bias–variance trade-offより

この現象をかなり網羅的に実験しているのが Deep Double Descent: Where Bigger Models and More Data Hurtである．この論文から得られる教訓は，必ずしも

• モデルが大きすぎると良くない（古典的な統計学）

• モデルが大きいほど良い（現代（当時？）の機械学習）

• データが多いほど良い

わけではないことである*1．

モデルサイズを大きくすると性能に悪影響を与えるモデルサイズの領域がある（There is a regime where bigger models are worse.）
図はDeep Double Descent: Where Bigger Models and More Data Hurt の Figure1 Left

サンプルを増加させると（同じモデルサイズでも）性能に悪影響を与える領域がある（There is a regime where more samples hurts.）
Deep Double Descent: Where Bigger Models and More Data Hurt の Figure3

長く訓練するとオーバーフィットから戻ってくる場合がある（There is a regime where training longer reverses overfitting.）
Deep Double Descent: Where Bigger Models and More Data Hurt の Figure2 Left + Figure9

この現象はもっと単純なモデルでもみることができる．2 層のニューラルネットで，入力層と中間層の間は適当にサンプリングして固定し，中間層と出力層だけ学習する（今回は適当なリッジ回帰を行なっている），というモデルを考える（Random Feature 回帰）．データとしては2次元で $(\cos\theta,\sin\theta)$ を20個与え，出力として $\sin\theta$ を出してもらうように学習を行う．中間層のパラメータの数によって次のように変化する．

パラメータが少ないとそこそこ合う，パラメータが20くらいになると過学習の典型的な見た目をするが，パラメータ数をものすごく増やすと安定する．

インタラクティブにプロットすると以下のようになる．下のスライダーを動かすとパラメータ数を変えたときにどのような関数が学習で得られたかをみることができる．

線形回帰モデルでの二重降下

深層学習モデルでみられたこの現象だが，なんと単純な線形回帰モデルにおいても観察することができる．ということで線形回帰モデルでこの現象を導出してみようというのがここからの話である．以下では次の問題設定で話を進めていく．

問題設定 特定モデルの状況で考える（観測ノイズを除いてデータの関係を説明できる状況）．

• 訓練データ: $\mathcal{D}=\lbrace(\boldsymbol{x} _ i,y _ i)\rbrace _ {i=1} ^ n\subset\mathbb{R} ^ p\times\mathbb{R}$, $\boldsymbol{x} _ i\overset{\mathrm{i.i.d}}{\sim}\mathcal{N}(0,I _ p)$
• 真のモデル: 線形モデル $$ y _ i = \boldsymbol{x} _ i ^ \top\boldsymbol{\theta}+ε _ i,\quad ε _ i\overset{\mathrm{i.i.d}}{\sim}\mathcal{N}(0,σ ^ 2) $$ 行列形式で書くと，$X=(\boldsymbol{x} _ 1,\ldots,\boldsymbol{x} _ n) ^ \top\in\mathbb{R} ^ {n\times p}, \boldsymbol{y}=(y _ 1,\ldots,y _ n) ^ \top, \boldsymbol{ε}=(ε _ 1,\ldots,ε _ n) ^ \top$ $$ \boldsymbol{y}=X\boldsymbol{\theta}+\boldsymbol{ε} $$
• パラメータの推定: $$ \hat{\boldsymbol{\theta}}=\arg\min _ {\boldsymbol{\theta}}\lVert \boldsymbol{y}-X\boldsymbol{\theta} \rVert _ 2 ^ 2 $$ これを勾配降下法で求めることを考える

目標はパラメータ数（もしくはデータ数/パラメータ数）を変数として汎化誤差の典型評価: $E[(y-\hat{y}) ^ 2]$ の式を求める（ $\boldsymbol{x}$ に対して $\hat{y}=\boldsymbol{x} ^ \top\hat{\boldsymbol{\theta}}$ とする）こと．

非漸近的な解析と，漸近的な解析の 2 つをここでは試してみる．前者はわりと線形回帰に特有な解析手法だけどわかりやすくて，後者の方が汎用的になっていると思う．

勾配降下法の陰的バイアス

汎化誤差の変形を始める前に，使用する勾配降下で求まる解の性質について言及しておこう．勾配降下法

$$ \hat{\boldsymbol{\theta}} _ {t+1}=\hat{\boldsymbol{\theta}} _ {t} - ηX ^ \top(X\hat{\boldsymbol{\theta}} _ t-\boldsymbol{y}) $$

で解を求めるとき， $p$ と $n$ の大小関係で得られる解が変化することが知られている．具体的には

• Underparametrized(優決定系): $p<n$ の場合には最小二乗解 $$ \hat{\boldsymbol{\theta}}=(X ^ \top X) ^ {-1}X ^ \top\boldsymbol{y} $$ が得られ，これは $$ \hat{\boldsymbol{\theta}}=\arg\min _ {\boldsymbol{\theta}}\lVert \boldsymbol{y}-X\boldsymbol{\theta} \rVert _ 2 ^ 2 $$ の解である．

• Overparametrized(劣決定系): $p>n$ の場合は最小ノルム解 $$ \hat{\boldsymbol{\theta}}=X ^ \top(XX ^ \top) ^ {-1}\boldsymbol{y} $$ が得られ，これは $$ \hat{\boldsymbol{\theta}}=\arg\min _ {\boldsymbol{\theta}}\lVert \boldsymbol{\theta} \rVert _ 2 ^ 2\,\mathrm{s.t.}\,\boldsymbol{y}=X\boldsymbol{\theta} $$ の解である．

後者で最小ノルム解が得られるのはある種非自明かもしれない．このように，選択したアルゴリズム（今回は勾配降下法）などの影響によって，陽にバイアスを加えたわけではないのに陰にバイアスが加わることを陰的バイアスと言ったりする．このことを証明できる形でちゃんと書くと次のような形になる．

勾配降下法の陰的バイアス $p>n$ で $\operatorname{rank} X=n$ ，固定学習率 $η\in (0,2/λ _ {\mathrm{max}}(X))$ の初期値 $\boldsymbol{\theta} _ 0 \in\mathrm{Im}(X ^ \top)$ から始める勾配降下法で $\hat{\boldsymbol{\theta}}=\arg\min _ {\boldsymbol{\theta}}\lVert \boldsymbol{y}-X\boldsymbol{\theta} \rVert _ 2 ^ 2$ の解を求めると，最小ノルム解に収束する： $$ \hat{\boldsymbol{\theta}}=\arg\min _ {\boldsymbol{\theta}}\lVert \boldsymbol{\theta} \rVert _ 2 ^ 2\,\mathrm{s.t.}\,\boldsymbol{y}=X\boldsymbol{\theta} $$

このステートメントは直感的には次のように説明できる．

• 毎ステップ $X ^ \top$ がかけられたベクトルが足されていく．
• ある $\boldsymbol{\alpha}$ があって $\lim _ {t\to\infty}\hat{\boldsymbol{\theta}} _ t=X ^ \top\boldsymbol{\alpha}$ のようになっているはず．
• $\boldsymbol{y}=XX ^ \top\boldsymbol{\alpha}$ をとくと $\boldsymbol{\alpha}=(XX ^ \top) ^ {-1}\boldsymbol{y}$．$\therefore \lim _ {t\to\infty}\hat{\boldsymbol{\theta}} _ t=X ^ \top(XX ^ \top) ^ {-1}\boldsymbol{y}$

しっかり証明するなら次のようになる．

$X$ を特異値分解する（$U \in \mathbb{R} ^ {n \times n}$ と $V \in \mathbb{R} ^ {p \times p}$ は直交行列，$Σ \in \mathbb{R} ^ {n \times p}$ は長方形の対角行列，$Σ _ {1} \in \mathbb{R} ^ {n \times n}$ は対角行列）：

$$ X=UΣ V ^ {\top}=U \begin{pmatrix} Σ _ {1} & 0 \end{pmatrix} \begin{pmatrix} V _ {1} ^ {\top} \cr V _ {2} ^ {\top} \end{pmatrix} $$

このとき，最小ノルム解 $\boldsymbol{\theta} ^ {*}$ は次のように書き換えることができる：

$$ \boldsymbol{\theta} ^ {*}=X ^ {\top}(X X ^ {\top}) ^ {-1}\boldsymbol{y}=V _ {1}Σ _ {1} ^ {-1}U ^ {\top}\boldsymbol{y} $$

勾配降下法の更新ルールは次の通り（ここで $η > 0$ はステップサイズ）：

$$ \boldsymbol{\theta} _ {t+1} = \boldsymbol{\theta} _ {t} - η \nabla L(\boldsymbol{\theta} _ t) = \boldsymbol{\theta} _ {t} - η X ^ {\top}(X \boldsymbol{\theta} _ {t} - \boldsymbol{y}) = (I - η X ^ {\top}X)\boldsymbol{\theta} _ {t} + η X ^ {\top}\boldsymbol{y} $$

これを繰り返すと，次のようになる：

$$ \boldsymbol{\theta} _ {t} = (I - η X ^ {\top}X) ^ t \boldsymbol{\theta} _ {0} + η \sum _ {s=0} ^ {t-1} (I - η X ^ {\top}X) ^ s X ^ {\top} \boldsymbol{y} $$

$X$ の特異値分解を用いると $X ^ {\top}X = V Σ ^ {\top} Σ V ^ {\top}$ であり，$V$ は直交行列なので $V ^ {\top}V = I$ ．これらに注意すると

$$ \begin{split} \boldsymbol{\theta} _ {t} &= V(I - η Σ ^ {\top} Σ) ^ t V ^ {\top} \boldsymbol{\theta} _ {0} + η V \left( \sum _ {s=0} ^ {t-1} (I - η Σ ^ {\top} Σ) ^ s \right) Σ ^ {\top} U ^ {\top} \boldsymbol{y} \cr &= V \begin{pmatrix} (I - η Σ _ {1} ^ 2) ^ t & 0 \cr 0 & I \end{pmatrix} V ^ {\top} \boldsymbol{\theta} _ {0} + η V \left( \sum _ {s=0} ^ {t-1} \begin{pmatrix} (I - η Σ _ {1} ^ 2) ^ s & 0 \cr 0 & I \end{pmatrix} \right) Σ ^ {\top} U ^ {\top} \boldsymbol{y} \end{split} $$

となる．

$0 \lt η \lt 2 / σ _ {\max} (Σ _ 1)$ を選ぶと，$σ _ {\max} (Σ _ 1)$ は $Σ _ 1$ の最大特異値であり，$I - η Σ ^ \top Σ$ の固有値の絶対値がすべて 1 未満であることが保証される．したがって $t\to\infty$ で次のようになる*2：

• 1 項目

  $$ V \begin{pmatrix} (I - η Σ _ 1 ^ 2) ^ t & 0 \cr 0 & I \end{pmatrix} V ^ {\top} \boldsymbol{\theta} _ {0} \rightarrow V \begin{pmatrix} 0 & 0 \cr 0 & I \end{pmatrix} V ^ {\top} \boldsymbol{\theta} _ {0} = V _ {2} V _ {2} ^ {\top} \boldsymbol{\theta} _ {0} $$

• 2 項目 $$ η\sum _ {s=0} ^ {t-1} \begin{pmatrix} (I - η Σ _ 1 ^ 2) ^ s Σ _ 1 \cr 0 \end{pmatrix} \rightarrow η \begin{pmatrix} \sum _ {s= 0} ^ {\infty} (I - η Σ _ 1 ^ 2) ^ s Σ _ 1 \cr 0 \end{pmatrix} = \begin{pmatrix} η(I-I+ηΣ _ 1 ^ 2) ^ {-1}Σ _ 1 \cr 0 \end{pmatrix} = \begin{pmatrix} Σ _ 1 ^ {-1} \cr 0 \end{pmatrix} $$

よって $\boldsymbol{\theta} _ {\infty}:=\lim _ {t\to\infty}\boldsymbol{\theta} _ t$ として，次のようになる：

$$ \boldsymbol{\theta} _ {\infty} = V _ 2 V _ 2 ^ {\top}\boldsymbol{\theta} _ {0} + V _ 1 Σ _ 1 ^ {-1} U ^ {\top} \boldsymbol{y} = V _ 2 V _ 2 ^ {\top} \boldsymbol{\theta} _ {0} + X ^ {\top}(X X ^ {\top}) ^ {-1} \boldsymbol{y} = V _ {2} V _ {2} ^ {\top} \boldsymbol{\theta} _ {0} + \boldsymbol{\theta} ^ {*} $$

ここで $\boldsymbol{\theta} _ {0}$ は $X ^ {\top}$ での写像先にあるので，$\boldsymbol{\theta} _ {0} = X ^ {\top} \boldsymbol{z}$ と書ける（ただし $\boldsymbol{z} \in \mathbb{R} ^ n$）．

$$ V _ 2 V _ 2 ^ {\top} \boldsymbol{\theta} _ 0 = V \begin{pmatrix} 0 & 0 \cr 0 & I \end{pmatrix} V ^ {\top} X ^ {\top} \boldsymbol{z} = V \begin{pmatrix} 0 & 0 \cr 0 & I \end{pmatrix} V ^ \top V Σ U ^ {\top} \boldsymbol{z} = V \begin{pmatrix} 0 & 0 \cr 0 & I \end{pmatrix} \begin{pmatrix} Σ _ 1 \cr 0 \end{pmatrix} U ^ {\top} \boldsymbol{z} = 0 $$

したがって， $\boldsymbol{\theta} _ {\infty}=\boldsymbol{\theta} ^ *$ となり，勾配降下法で求めた解は最小ノルム解に収束する．

非漸近的な導出

得られる解の形がわかったので，ここからは汎化誤差を変形して具体的な式を求めていく．汎化誤差 $E[(y-\hat{y}) ^ 2]$ は

$$ \begin{split} E[(y-\hat{y}) ^ 2] &=E[(\boldsymbol{x} ^ \top\boldsymbol{\theta} + ε - \boldsymbol{x} ^ \top\hat{\boldsymbol{\theta}}) ^ 2] \cr &=E[(\boldsymbol{x} ^ \top\boldsymbol{\theta} - \boldsymbol{x} ^ \top\hat{\boldsymbol{\theta}}) ^ 2] + E[ε ^ 2] \cr &=E[\boldsymbol{x} ^ \top(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}})(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}}) ^ \top \boldsymbol{x}] + E[ε ^ 2] \cr &=E[\mathrm{Tr}( (\boldsymbol{\theta} - \hat{\boldsymbol{\theta}})(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}}) ^ \top \boldsymbol{x}\boldsymbol{x} ^ \top)] + E[ε ^ 2] \cr &=\mathrm{Tr}(E[(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}})(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}}) ^ \top)] E[\boldsymbol{x}\boldsymbol{x} ^ \top] + E[ε ^ 2] \cr &=\mathrm{Tr}(E[(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}})(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}}) ^ \top]) + E[ε ^ 2] \cr &=E[\lVert \boldsymbol{\theta} - \hat{\boldsymbol{\theta}} \rVert ^ 2] + σ ^ 2 \end{split} $$

となるので $E[\lVert \boldsymbol{\theta} - \hat{\boldsymbol{\theta}} \rVert ^ 2]$ について考えたい．

先ほどもみたように，$p$ と $n$ の大小関係によって得られる解が変化するので，$p\leq n$ と $p\gt n$で場合分けを行う（当然，得られる解が変わるので汎化誤差の表式も変わる）．

まずは $p\leq n$ の場合を考える．

$$ \begin{split} \boldsymbol{\theta} - \hat{\boldsymbol{\theta}} &=\boldsymbol{\theta} - (X ^ \top X) ^ {-1}X ^ \top\boldsymbol{y}\cr &=\boldsymbol{\theta} - (X ^ \top X) ^ {-1}X ^ \top(X\boldsymbol{\theta}+\boldsymbol{ε})\cr &=\boldsymbol{\theta} - (X ^ \top X) ^ {-1}X ^ \top X\boldsymbol{\theta} - (X ^ \top X) ^ {-1}X ^ \top \boldsymbol{ε}\cr &=-(X ^ \top X) ^ {-1}X ^ \top\boldsymbol{ε} \end{split} $$

ここで $\boldsymbol{\theta} - (X ^ \top X) ^ {-1}X ^ \top X\boldsymbol{\theta}$ は $(X ^ \top X) ^ {-1}X ^ \top X=I$ より $0$ になることに注意する．これより

$$ \begin{split} E[\lVert \boldsymbol{\theta} - \hat{\boldsymbol{\theta}} \rVert ^ 2] &=E[\lVert -(X ^ \top X) ^ {-1}X ^ \top\boldsymbol{ε} \rVert ^ 2]\cr &=E[\boldsymbol{ε} ^ \top X(X ^ \top X) ^ {-\top}(X ^ \top X) ^ {-1}X ^ \top\boldsymbol{ε}] \cr &=E[\mathrm{Tr}(X(X ^ \top X) ^ {-\top}(X ^ \top X) ^ {-1}X ^ \top\boldsymbol{ε}\boldsymbol{ε} ^ \top)] \cr &=\mathrm{Tr}(E[X(X ^ \top X) ^ {-\top}(X ^ \top X) ^ {-1}X ^ \top]E[\boldsymbol{ε}\boldsymbol{ε} ^ \top]) \cr &=σ ^ 2\mathrm{Tr}(E[X(X ^ \top X) ^ {-\top}(X ^ \top X) ^ {-1}X ^ \top]) \cr &=σ ^ 2 E[\mathrm{Tr}( (X ^ \top X) ^ {-\top}(X ^ \top X) ^ {-1}X ^ \top X)] \cr &=σ ^ 2 \mathrm{Tr}(E[(X ^ \top X) ^ {-1}]) \cr &=σ ^ 2 \frac{p}{n-p-1} \end{split} $$

となる．最後の等式が成り立つことを言うには少し複雑なので補足として後ろに載せておく．inverse-Wishart分布とかで調べると証明は出てくる．以上から

$$ \therefore E[(y-\hat{y}) ^ 2] = σ ^ 2 \left(1 + \frac{p}{n-p-1}\right) $$

と導かれる．

次に $p\geq n$ の場合は

$$ \begin{split} \boldsymbol{\theta} - \hat{\boldsymbol{\theta}} &=\boldsymbol{\theta} - X ^ \top(XX ^ \top) ^ {-1}\boldsymbol{y}\cr &=\boldsymbol{\theta} - X ^ \top(XX ^ \top ) ^ {-1}(X\boldsymbol{\theta}+\boldsymbol{ε})\cr &=\boldsymbol{\theta} - X ^ \top(XX ^ \top) ^ {-1}X\boldsymbol{\theta} - X ^ \top(XX ^ \top ) ^ {-1} \boldsymbol{ε}\cr &=(I-X ^ \top(XX ^ \top) ^ {-1}X)\boldsymbol{\theta}-X ^ \top(XX ^ \top ) ^ {-1}\boldsymbol{ε} \end{split} $$

$$ \begin{split} E[\lVert \boldsymbol{\theta} - \hat{\boldsymbol{\theta}} \rVert ^ 2] &=E[\lVert (I-X ^ \top(XX ^ \top) ^ {-1}X)\boldsymbol{\theta}-X ^ \top(XX ^ \top ) ^ {-1}\boldsymbol{ε} \rVert ^ 2]\cr &=E[\lVert (I-X ^ \top(XX ^ \top) ^ {-1}X)\boldsymbol{\theta} \rVert ^ 2]+E[\lVert X ^ \top(XX ^ \top ) ^ {-1}\boldsymbol{ε} \rVert ^ 2] \end{split} $$

ここで

• $E[\lVert (I-X ^ \top(XX ^ \top) ^ {-1}X)\boldsymbol{\theta} \rVert ^ 2]$ をシグナル（バイアス）
• $E[\lVert X ^ \top(XX ^ \top ) ^ {-1}\boldsymbol{ε} \rVert ^ 2]$ をノイズ（バリアンス）

とよぶことにしよう．バイアス/バリアンスといっているものの方が多そう（というか元論文はそうしている）だが，ここではわざと別の言葉を使うことにする（伏線は後で回収される）．

• 1 項目: シグナル 実は $X ^ \top(XX ^ \top) ^ {-1}X$ は射影行列になっているので ピタゴラスの定理が使える:

  $$ \lVert (I-X ^ \top(XX ^ \top) ^ {-1}X)\boldsymbol{\theta} \rVert ^ 2 = \lVert \boldsymbol{\theta} \rVert ^ 2 - \lVert X ^ \top(XX ^ \top) ^ {-1}X\boldsymbol{\theta} \rVert ^ 2 $$

  $X ^ \top(XX ^ \top) ^ {-1}X$ は $n(<p)$ 次元方向への射影なので

  $$ E[\lVert X ^ \top(XX ^ \top) ^ {-1}X\boldsymbol{\theta} \rVert ^ 2]=\frac{n}{p}\lVert \boldsymbol{\theta} \rVert ^ 2 $$

  あわせて

  $$ E[\lVert (I-X ^ \top(XX ^ \top) ^ {-1}X)\boldsymbol{\theta} \rVert ^ 2] = \lVert \boldsymbol{\theta} \rVert ^ 2\left(1-\frac{n}{p}\right) $$

• 2 項目: ノイズ

  $$ \begin{split} E[\lVert X ^ \top(XX ^ \top) ^ {-1}\boldsymbol{ε} \rVert ^ 2] &=E[\boldsymbol{ε} ^ \top (XX ^ \top) ^ {-\top} XX ^ \top(XX ^ \top) ^ {-1}\boldsymbol{ε}]\cr &=E[\operatorname{Tr}( (XX ^ \top) ^ {-1}\boldsymbol{ε}\boldsymbol{ε} ^ \top)]\cr &=\operatorname{Tr}(E[(XX ^ \top) ^ {-1}]E[\boldsymbol{ε}\boldsymbol{ε} ^ \top])\cr &=σ ^ 2\operatorname{Tr}(E[(XX ^ \top) ^ {-1}])\cr &=σ ^ 2\frac{n}{p-n-1} \end{split} $$

したがって

$$ E[\lVert \boldsymbol{\theta} - \hat{\boldsymbol{\theta}} \rVert ^ 2]=\lVert \boldsymbol{\theta} \rVert ^ 2\left(1-\frac{n}{p}\right)+σ ^ 2\frac{n}{p-n-1} $$

となり，全体としては

$$ \therefore E[(y-\hat{y}) ^ 2] = \lVert \boldsymbol{\theta} \rVert ^ 2 \left(1 - \frac{n}{p}\right)+σ ^ 2\left(1+\frac{n}{p-n-1}\right) $$

となる．

$p\leq n$ と $p\geq n$ の場合をまとめると汎化誤差は

$$ E[(y-\hat{y}) ^ 2] = \begin{cases} σ ^ 2 \left(1 + \frac{p}{n-p-1}\right) & (p<n-1)\cr \infty & (p=n-1,n,n+1)\cr \lVert \boldsymbol{\theta} \rVert ^ 2 \left(1 - \frac{n}{p}\right)+σ ^ 2\left(1+\frac{n}{p-n-1}\right) & (p>n+1) \end{cases} $$

と導ける．これを $p$ を横軸にしてプロットすると下の図のようになる．$p < n$ では単調に増加していて，途中で減少して増加するのような描像はみられなかったが， $p=n$ でピークを持つ山のような形になっているというのは観察できている．

$n=40$, $\lVert \boldsymbol{\theta} \rVert ^ 2=1$, $σ=1/5$ での例

補足

補足1

inverse-Wishart分布の平均 $x _ {ij}\overset{\mathrm{i.i.d}}{\sim}\mathcal{N}(0,1)$ とし，これを使って行列 $X=(x _ {ij}) _ {1\leq i\leq n;1\leq j\leq d}=(x _ 1,\ldots,x _ n) ^ \top$ と $A:= X ^ \top X=\sum _ {i=1} ^ n x _ ix _ i ^ \top$ を定める． このとき $$ E[\operatorname{Tr}(A ^ {-1})]=\frac{1}{n-d-1} $$ が成り立つ．

$$ \begin{split} E[x _ i ^ \top A ^ {-1}x _ i] &\overset{\text{1.}}{=}E\left\lbrack \sum _ {j=1} ^ d\frac{\partial}{\partial x _ {ij}}(e _ j ^ \top A ^ {-1}x _ i)\right\rbrack \cr &=\sum _ {j=1} ^ dE\left\lbrack e _ j ^ \top\left(\frac{\partial}{\partial x _ {ij}}A ^ {-1}\right)x _ i + e _ j ^ \top A ^ {-1}\frac{\partial}{\partial x _ {ij}} x _ i\right\rbrack \cr &\overset{\text{2.}}{=}\sum _ {j=1} ^ dE\left\lbrack e _ j ^ \top \left(-A ^ {-1}\left(\frac{\partial}{\partial x _ {ij}}A\right) A ^ {-1}\right)x _ i + e _ j ^ \top A ^ {-1} e _ j\right\rbrack \cr &=\sum _ {j=1} ^ dE\left\lbrack e _ j ^ \top \left(-A ^ {-1}\left(e _ j e _ i ^ \top X + X ^ \top e _ i e _ j ^ \top\right) A ^ {-1}\right)x _ i+e _ j ^ \top A ^ {-1} e _ j\right\rbrack \cr &\overset{\text{3.,4.}}{=}E[-\operatorname{Tr}(A ^ {-1})x _ i ^ \top A ^ {-1}x _ i-\sum _ {j=1} ^ d(e _ j ^ \top A ^ {-1}x _ i) ^ 2+\operatorname{Tr}(A ^ {-1})] \cr &=E[-\operatorname{Tr}(A ^ {-1})x _ i ^ \top A ^ {-1}x _ i-x _ i ^ \top A ^ {-2} x _ i + \operatorname{Tr}(A ^ {-1})] \end{split} $$

ここで

1. $x\sim\mathcal{N}(0,I _ d)$ と 有界連続な関数 $f:\mathbb{R} ^ d\to\mathbb{R}$ に対して $E[xf(x)]=E[\nabla f(x)]$ となること（Stein's Lemma）を使う．
2. $\frac{d}{dA _ {ij}}A ^ {-1}=-A ^ {-1}\left(\frac{d}{dA _ {ij}}A\right)A ^ {-1}$
3. $e _ j ^ \top A ^ {-1}e _ je _ i ^ \top XA ^ {-1}x _ i=e _ i ^ \top A ^ {-1}e _ j x _ i ^ \top A ^ {-1}x _ i$
4. $e _ j ^ \top A ^ {-1}X ^ \top e _ i e _ j ^ \top A ^ {-1}x _ i = e _ i ^ \top X A ^ {-1} e _ j e _ i ^ \top A ^ {-1} x _ i = (x _ i ^ \top A ^ {-1} e _ j)(e _ j ^ \top A ^ {-1} x _ i)$

である．これを使うと

$$ \begin{split} E[\operatorname{Tr}(XA ^ {-1}X)]&=E[-\operatorname{Tr}(A ^ {-1})\operatorname{Tr}(XA ^ {-1}X ^ \top)-\operatorname{Tr}(XA ^ {-2}X ^ \top)+n\operatorname{Tr}(A ^ {-1})] \cr d&=-E[\operatorname{Tr}(A ^ {-1})]\cdot d-E[\operatorname{Tr}(A ^ {-1})]+n\operatorname{Tr}(A ^ {-1}) \end{split} $$

よって

$$ E[\operatorname{Tr}(A ^ {-1})]=\frac{1}{n-d-1} $$

が導ける．

補足2

擬似逆行列を使うと，解を統一的に記述できる．

$$ \hat{\boldsymbol{\theta}}=X ^ +\boldsymbol{y}=X ^ +(X\boldsymbol{\theta}+\boldsymbol{ε})=X ^ +X\boldsymbol{\theta}+X ^ +\boldsymbol{ε} $$

Moore-Penrose pseudo-inverse 行列 $X$ が与えられたとき，次の 4 条件をみたすもののことで， $X ^ +$ と表す．

1. $XX ^ +X = X$
2. $X ^ +XX ^ + = X ^ +$
3. $(X X ^ + ) ^ \top = X X ^ +$
4. $(X ^ + X) ^ \top = X ^ +X$

特に，$X\in\mathbb{R} ^ {n\times p}$ に対して

• $p\lt n$ ならば $\mathrm{rank}X=p$ で $X ^ + = (X ^ \top X) ^ {-1} X ^ \top$
• $p=n$ ならば $X ^ +=X ^ {-1}$
• $n\lt p$ ならば $\mathrm{rank}X=n$ で $X ^ + = X ^ \top(XX ^ \top) ^ {-1}$

後のトピックでは次の 2 つの事実を使うことになる（はじめから擬似逆行列を使って記述しておいてもよかったかもしれない）．

• $X ^ +X=(X ^ +X) ^ \top$ であることから

  $$ X ^ +X=\begin{cases} I & (p\lt n)\cr X ^ \top(XX ^ \top) ^ {-1}X & (n\lt p) \end{cases} $$ ちなみに後者は $n$ 次元の $X$ の行方向への射影である．

• トレースについて $$ \begin{split} \operatorname{Tr}( (X ^ +) ^ \top X ^ +) &=\begin{cases} \operatorname{Tr}(X(X ^ \top X) ^ {-\top}(X ^ \top X) ^ {-1}X ^ \top)=\operatorname{Tr}( (X ^ \top X) ^ {-\top}(X ^ \top X) ^ {-1}X ^ \top X) & (p\lt n)\cr \operatorname{Tr}( (XX ^ \top) ^ {-\top}XX ^ \top(XX ^ \top) ^ {-1}) & (n\lt p) \end{cases}\cr &=\begin{cases} \operatorname{Tr}( (X ^ \top X) ^ {-1}) & (p\lt n)\cr \operatorname{Tr}( (XX ^ \top) ^ {-1}) & (n\lt p) \end{cases} \end{split} $$

漸近的導出

データ数 $n$ ，モデルのパラメータ数（＝データの次元数） $p$ として， $n,p\to\infty$ かつ $p/n\toγ\in(0,\infty)$ の状況で導出をしてみる．ここでは一度 Ridge 回帰を考えて $λ\to+0$ によって解を得ることを考える．

$$ \hat{\boldsymbol{\theta}} _ λ:=\arg\min _ {\boldsymbol{\theta}}\frac{1}{n}\lVert \boldsymbol{y}-X\boldsymbol{\theta} \rVert ^ 2+λ\lVert \boldsymbol{\theta} \rVert ^ 2=(X ^ \top X+λ nI) ^ {-1}X ^ \top\boldsymbol{y},\quad\hat{\boldsymbol{\theta}}:=\lim _ {λ\to+0}\hat{\boldsymbol{\theta}} _ λ $$

このとき $E[(y-\hat{y}) ^ 2]$ から $σ ^ 2$ を除いたものは $E _ {\boldsymbol{x}}[\lVert \boldsymbol{x} ^ \top\boldsymbol{\theta}-\boldsymbol{x} ^ \top\hat{\boldsymbol{\theta}} _ λ \rVert ^ 2]=\lVert \boldsymbol{\theta}-\hat{\boldsymbol{\theta}} _ λ \rVert ^ 2$ ．ここではデータについての平均までとり，極限での汎化誤差の典型値

$$ R(γ):=\lim _ {λ\to+0}\lim _ {n,p\to\infty, p/n\toγ}E[\lVert \boldsymbol{\theta}-\hat{\boldsymbol{\theta}} _ λ \rVert ^ 2] $$

を求めることを目指す．

（典型）汎化誤差の変形

まず典型汎化誤差を分解する．

$$ \begin{split} E[\lVert \boldsymbol{\theta}-\hat{\boldsymbol{\theta}} _ λ \rVert ^ 2] &=E[\lVert \boldsymbol{\theta} - (X ^ \top X+λ nI) ^ {-1}X ^ \top(X\boldsymbol{\theta}+\boldsymbol{ε}) \rVert ^ 2]\cr &=E[\lVert \boldsymbol{\theta} - (X ^ \top X+λ nI) ^ {-1}X ^ \top X\boldsymbol{\theta}+(X ^ \top X+λ nI) ^ {-1}X ^ \top \boldsymbol{ε} \rVert ^ 2]\cr &=E[\lVert \boldsymbol{\theta} - (X ^ \top X+λ nI) ^ {-1}X ^ \top X\boldsymbol{\theta} \rVert ^ 2]+E[\lVert (X ^ \top X+λ nI) ^ {-1}X ^ \top \boldsymbol{ε} \rVert ^ 2]\cr &=E _ X\left\lbrack \lVert \boldsymbol{\theta}-E _ {\boldsymbol{ε}}[\hat{\boldsymbol{\theta}}] \rVert ^ 2+E _ {\boldsymbol{ε}}[\lVert \hat{\boldsymbol{\theta}}-E _ {\boldsymbol{ε}}[\hat{\boldsymbol{\theta}}] \rVert ^ 2]\right\rbrack \end{split} $$

ここで $B _ {X}(λ):=\lVert \boldsymbol{\theta}-E _ {\boldsymbol{ε}}[\hat{\boldsymbol{\theta}}] \rVert ^ 2$ と $V _ {X}(λ):=E _ {\boldsymbol{ε}}[\lVert \hat{\boldsymbol{\theta}}-E _ {\boldsymbol{ε}}[\hat{\boldsymbol{\theta}}] \rVert ^ 2$ とおくことにすると

$$ \begin{split} B _ {X}(λ) &=\lVert (I - (X ^ \top X+λ nI) ^ {-1}X ^ \top X)\boldsymbol{\theta} \rVert ^ 2\cr &=\lVert (I - (X ^ \top X+λ nI) ^ {-1}(X ^ \top X+λ nI-λ nI) )\boldsymbol{\theta} \rVert ^ 2\cr &=\lVert λ n(X ^ \top X+λ nI) ^ {-1}\boldsymbol{\theta} \rVert ^ 2\cr &=\lVert \boldsymbol{\theta} \rVert ^ 2λ ^ 2\frac{1}{p}\operatorname{Tr}\left(\left(\frac{1}{n}X ^ \top X+λ I\right) ^ {-2}\right)\cr &=\lVert \boldsymbol{\theta} \rVert ^ 2λ ^ 2\int _ {(0,\infty)}\frac{1}{(s+λ) ^ 2}d\hat{μ}(s) \end{split} $$

$$ \begin{split} V _ {X}(λ) &=E _ {\boldsymbol{ε}}[\lVert (X ^ \top X+λ nI) ^ {-1}X ^ \top \boldsymbol{ε} \rVert ^ 2]\cr &=σ ^ 2\operatorname{Tr}( (X ^ \top X+λ nI) ^ {-2}X ^ \top X)\cr &=σ ^ 2\frac{p}{n}\frac{1}{p}\operatorname{Tr}\left(\left(\frac{1}{n}X ^ \top X+λ I\right) ^ {-2}\frac{1}{n}X ^ \top X\right)\cr &=σ ^ 2 \frac{p}{n}\int _ {(0,\infty)}\frac{s}{(s+λ) ^ 2}\hat{μ}(ds) \end{split} $$

と積分を利用した形に書き換えることができる．ここで，$(X ^ \top X+λ nI) ^ {-2}$ はどの成分についても等しい分布なので，適当な直交行列を使って $\boldsymbol{\theta}$ を 1 つの成分だけが非零になるように変換して，トレースで成分の数だけ割る，といったことをしている．また $\hat{μ}(s)=\frac{1}{p}\sum _ {j=1} ^ pδ _ {λ _ j}(s)$ は$X ^ \top X/n$ の固有値計数測度．

ここからは $B:=\lim _ {n,p\to\infty,\,p/n\toγ}B _ {X}$, $V:=\lim _ {n,p\to\infty,\,p/n\toγ}V _ {X}$ として，これらを求めることを考える．

わざわざこんな変形をしたのは，ランダム行列理論という強力な道具を使えるからである．

ランダム行列理論からの道具

Marchenko-Pastur law

マルチェンコ・パスツール分布（Marchenko-Pastur distribution, Marchenko-Pastur law）は，次の条件をみたすようなランダム行列の固有値が（漸近的に）従う分布である．

Theorem (Marchenko-Pastur) $\lbrace x _ i\rbrace _ {i=1} ^ n\subset\mathbb{R} ^ p$ を独立同分布で $E[x _ i]=0, E[x _ ix _ i ^ \top]=σ ^ 2I\,(σ ^ 2<+\infty)$ をみたす確率変数とする．$X=(x _ 1,\ldots,x _ n)$ とし $XX ^ \top /n\in\mathbb{R} ^ {p\times p}$ の固有値を $λ _ i$ とする．$p,n\to\infty$, $p/n\toγ\in (0,\infty)$ のとき，任意の有界連続関数 $f$ に対して $$ \frac{1}{p}\sum _ {j=1} ^ p f(λ _ j)=\int f(s)\hat{μ}(ds)\longrightarrow\int f(s)μ(ds)\quad(\mathrm{a.s.}) $$ が成り立つ（分布収束）．ただし $μ$ は次をみたす**（Marchenko-Pastur distribution, Marchenko-Pastur law）**． $$ \frac{dμ(s)}{ds}=\frac{1}{2\pi}\frac{\sqrt{(λ _ +-s)(s-λ _ -)}}{γ s}1 _ {s\in [λ _ -,λ _ +]}+\max\left(0,1-\frac{1}{γ}\right)δ(s),\quad λ _ \pm=(1\pm\sqrt{γ}) ^ 2 $$

この分布は，自由確率論において自由ポアソン分布ともいわれることもある．古典的な確率論におけるポアソン少数の法則において，独立性の条件を自由独立性に変えるとこの分布が出てくる．

また，これを利用することから，実際には $X$ の期待値をとらなくても続く結果が成り立つことがわかる．

Stieltjes Transform

確率論における特性関数を使った変換のようなものは，ランダム行列理論においてはスティルチェス変換によって行われる．

Stieltjes Transform $\mathbb{R}$ 上の測度 $μ$ に対して次で定義される $S:\mathbb{C}\backslash\mathrm{supp}(μ)\to\mathbb{C}$ を与える変換 $μ\mapsto S$ のこと． $$ S(z)=\int \frac{1}{x-z}μ(dx) $$

いくつか今回使う性質にかぎって紹介すれば

• モーメントとの関係： $$ S(z)=-\sum _ {k=0} ^ \infty\frac{m _ k}{z ^ {k+1}},\quad m _ k:=\int x ^ k μ(dx) $$
• Stieltjes Inverse Transform: $μ(dx)=\rho(x)dx$ のとき $$ \rho(x)=\lim _ {ε\to+0}\frac{1}{\pi}\operatorname{Im}S(x+iε) $$
• $μ$ として $n\times n$ 行列 $X$ の固有値計数測度 $\hat{μ}(x)=\frac{1}{n}\sum _ {j=1} ^ nδ _ {λ _ j}(x)$ をとれば $$ S _ n(z)=\frac{1}{n}\operatorname{Tr}(X-zI) ^ {-1} $$

といった関係がある．上で紹介した Marchenko-Pastur 分布に対しては次のような変換式を導出できる．

Marchenko-Pastur law に対する Stieltjes Transform $$ S(z)=\frac{(1-γ-z)-\sqrt{(1-γ-z) ^ 2-4γ z}}{2γ z} $$

ちなみに，Cauchy 変換というのもあり，これは Stieltjes 変換の逆符号で定義される．ただし，Cauchy 変換と Stieltjes 変換を同じものとして定義しているものもあるので注意．

補足: Marchenko-Pastur law の導出

$X=(x _ 1,\ldots,x _ n)\in\mathbb{R} ^ {p\times n}$, $x _ i \sim\mathcal{N}(0,I _ p)$ とし， $XX ^ \top /n\in\mathbb{R} ^ {p\times p}$ とおく．

$Σ=\begin{pmatrix}Σ _ {11} & Σ _ {\cdot,1} ^ \top \cr Σ _ {\cdot,1} & Σ _ {p-1} \end{pmatrix}$ と書く． $Σ _ {11}\in\mathbb{R}$, $Σ _ {\cdot,1}\in\mathbb{R} ^ {p-1}$, $Σ _ {p-1}\in\mathbb{R} ^ {(p-1)\times(p-1)}$ ．

対称性から

$$ E[\left( (Σ-zI) ^ {-1}\right) _ {11}]=\frac{1}{n}E[\operatorname{Tr}(Σ-zI) ^ {-1}] $$

になるので，$\left( (Σ-zI) ^ {-1}\right) _ {11}$ について考えればよく，また

$$ \frac{1}{\left( (Σ-zI) ^ {-1}\right) _ {11}}=E\left\lbrack \frac{1}{\left( (Σ-zI) ^ {-1}\right) _ {11}}\right\rbrack +O(n ^ {-1/2})=\frac{1}{E[\left( (Σ-zI) ^ {-1}\right) _ {11}]}+O(n ^ {-1/2}) $$

になることが示せる．

Schur complement formula から

$$ \frac{1}{\left( (Σ-zI) ^ {-1}\right) _ {11}}=(z-Σ _ {11})-Σ _ {\cdot,1} ^ \top(Σ _ {p-1}-zI)Σ _ {\cdot,1} $$

が成り立つ．また

$$ \begin{split} E\left\lbrack \frac{1}{\left( (Σ-zI) ^ {-1}\right) _ {11}}\right\rbrack &=E\lbrack(z - Σ _ {11})-Σ _ {\cdot,1} ^ \top(Σ _ {p-1} - z I )Σ _ {\cdot,1}\rbrack \cr &= z - 1 - \frac{1}{n ^ 2} E \left\lbrack \sum _ {s , t = 1} ^ n \sum _ {j , k = 2} ^ p X _ {1 t} X _ {j t} (Σ _ {p - 1} - z I) _ {jk} ^ {-1} X _ {ks} X _ {1s} \right\rbrack \cr &= z - 1 - \frac{1}{n} E \left\lbrack \operatorname{Tr}Σ _ {p-1}(Σ _ {p-1}-zI) ^ {-1}\right\rbrack \cr &= z - 1 + \frac{p-1}{n} -z \frac{1}{n} E \left\lbrack \operatorname{Tr}(Σ _ {p-1} - zI) ^ {-1}\right\rbrack \cr &= \left(\frac{1}{n}E\lbrack \operatorname{Tr}(Σ - zI) ^ {-1}\rbrack\right) ^ {-1} + O(n ^ {-1/2}) \end{split} $$

ここで， $n,p\to\infty$, $p/n\toγ\in (0,\infty)$ を考えれば， $S(z)=\lim \frac{1}{n}E[\operatorname{Tr}(Σ-zI) ^ {-1}]$ として

$$ \frac{1}{S(z)}=z-1+γ-γ z S(z) $$

これを，固有値分布の性質を満たす解を選ぶことで（$S(z)\sim 1/z, |z|\to\infty$）Marchenko-Pastur law に対する Stieltjes Transform

$$ \begin{split} S(z)&=\frac{(1-γ-z)-\sqrt{(1-γ-z) ^ 2-4γ z}}{2γ z}\cr &=\frac{z-(1-γ)-\sqrt{(z-λ _ +)(z-λ _ -)}}{2γ z},\quadλ _ {\pm}=(1\pm\sqrt{γ}) ^ 2 \end{split} $$

を得ることができる．逆変換をすれば

$$ \rho(x)=\frac{1}{\pi}\lim _ {ε\to+0}\mathrm{Im}S(x+iε)=\frac{\sqrt{(λ _ +-x)(x-λ _ -)}}{2\piγ x}1 _ {x\in[λ _ -,λ _ +]} $$

もし $γ>1$ であれば $(1-1/γ)δ(x)$ が加わり，結局

$$ \rho(x)=\frac{\sqrt{(λ _ +-x)(x-λ _ -)}}{2\piγ x}1 _ {x\in[λ _ -,λ _ +]}+\max\left(0,1-\frac{1}{γ}\right)δ(x),\quad λ _ \pm=(1\pm\sqrt{γ}) ^ 2 $$

が得られる．

B と V の Stieltjes Transform による表示

簡単に使う道具について紹介したところで，$n,p\to\infty$ として， $B$ と $V$ の具体的な式を求める．以下で示すように，$B$ と $V$ は $S$ によって表現できて，これから求めることができる． $S'(-λ) = \frac{d}{dλ}S(-λ)$ とかくと

$$ B(λ)=\lVert \boldsymbol{\theta} \rVert ^ 2λ ^ 2\int _ {(0,\infty)}\frac{1}{(s+λ) ^ 2}dμ(s)=-\lVert \boldsymbol{\theta} \rVert ^ 2λ ^ 2S'(-λ) $$

$$ \begin{split} V(λ) &=σ ^ 2γ \int _ {(0,\infty)}\frac{s}{(s+λ) ^ 2}dμ(s)\cr &=σ ^ 2γ\left( \int _ {(0,\infty)}\frac{1}{s+λ}dμ(s)- \int _ {(0,\infty)}\frac{λ}{(s+λ) ^ 2}dμ(s)\right)\cr &=σ ^ 2γ(S(-λ)+λ S'(-λ)) \end{split} $$

ただし

$$ S(-λ)=\frac{-(1-γ+λ)+\sqrt{(1-γ+λ) ^ 2+4γ λ}}{2γ λ} $$

$$ \begin{split} S'(-λ) &=\frac{\frac{2(-γ+λ+1)+4γ}{2\sqrt{(-γ+λ+1) ^ 2+4γλ}}-1}{2γλ}-\frac{-(1-γ+λ)+\sqrt{(1-γ+λ) ^ 2+4γ λ}}{2γ λ ^ 2}\cr &=\frac{\frac{2(-γ+λ+1)+4γ}{2\sqrt{(-γ+λ+1) ^ 2+4γλ}}-1}{2γλ}-\frac{1}{λ}S(-λ) \end{split} $$

$λ\to +0$で

$$ \begin{split} -λ ^ 2S'(-λ) &=-\left(\frac{\frac{2(-γ+λ+1)+4γ}{2\sqrt{(-γ+λ+1) ^ 2+4γλ}}-1}{2γ}λ\right)+\frac{-(1-γ+λ)+\sqrt{(1-γ+λ) ^ 2+4γ λ}}{2γ}\cr &\to 0 + \frac{-(1-γ)+|1-γ|}{2γ}\cr &=\begin{cases} 0 & (γ\leq 1)\cr 1-\frac{1}{γ} &(γ>1) \end{cases} \end{split} $$

$$ \begin{split} S(-λ)+λ S'(-λ) &=\frac{\frac{2(-γ+λ+1)+4γ}{2\sqrt{(-γ+λ+1) ^ 2+4γλ}}-1}{2γ}\cr &\to\frac{1}{2γ}\frac{γ+1-|-γ+1|}{|-γ+1|}\cr &=\begin{cases} \frac{1}{1-γ} & (γ\leq 1)\cr \frac{1}{γ}\frac{1}{γ-1} &(γ>1) \end{cases} \end{split} $$

以上から

$$ \therefore B(0)=\lim _ {λ\to +0}-\lVert \boldsymbol{\theta} \rVert ^ 2λ ^ 2S'(-λ)=\begin{cases} 0 & (γ\leq 1)\cr \lVert \boldsymbol{\theta} \rVert ^ 2\left(1-\frac{1}{γ}\right) &(γ>1) \end{cases}, $$

$$ \begin{split} V(0)=\lim _ {λ\to +0}σ ^ 2γ (S(-λ)+λ S'(-λ))=\begin{cases} σ ^ 2\frac{γ}{1-γ} & (γ\leq 1)\cr σ ^ 2\frac{1}{γ-1} & (γ>1) \end{cases} \end{split} $$

と求めることができた．したがって

$$ R(γ) = B(0)+V(0) = \begin{cases} σ ^ 2 \frac{γ}{1-γ} & (γ\leq 1)\cr \lVert \boldsymbol{\theta} \rVert ^ 2 \left(1 - \frac{1}{γ}\right)+σ ^ 2\frac{1}{γ-1} & (γ>1) \end{cases} $$

となる．非漸近的導出の際に求めた $E[\lVert \boldsymbol{\theta}-\hat{\boldsymbol{\theta}} \rVert ^ 2]$ と比較すると，極限をとればちゃんと一致していることが確認できる．

バイアス・バリアンス分解による解釈

ここまでで導出してきた汎化誤差が，どの変数の影響でピークが現れることになるかを解釈したい．それのヒントとなるのがバイアス・バリアンス分解である．

ランダムデザイン・固定デザイン・より細かい分解

一般に $y(\boldsymbol{x})=f(\boldsymbol{x})+ε$, $ε\sim\mathcal{N}(0,σ ^ 2)$ に対して $f$ を $\hat{f}$ で推定するときのバイアス・バリアンス分解は次のように与えられる．今回のような線形回帰の場合には $f=\boldsymbol{\theta} ^ \top$, $\hat{f}=\hat{\boldsymbol{\theta}} ^ \top$ に対応している）．

• ランダムデザイン（Random Design）: 計画行列 $X$ の期待値をとる

  $$ E _ {\boldsymbol{x}}E _ XE _ {\boldsymbol{ε}}[(y(\boldsymbol{x})-\hat{f}(\boldsymbol{x})) ^ 2]=E _ {\boldsymbol{x}}[(f(\boldsymbol{x}) - E _ {X,\boldsymbol{ε}}[\hat{f}(\boldsymbol{x})]) ^ 2] + E _ {\boldsymbol{x}}E _ {X,\boldsymbol{ε}}[(\hat{f}(\boldsymbol{x}) - E _ {X,\boldsymbol{ε}}[\hat{f}(\boldsymbol{x})]) ^ 2] + σ ^ 2 $$

例えば PRML などを参照するとこれがバイアス・バリアンス分解として紹介されるが，しばしば以下のようなバイアス・バリアンス分解が採用される場合もある．

• 固定デザイン（Fixed Design）: 計画行列 $X$ の条件下で考える

  $$ E _ {\boldsymbol{x}}E _ {\boldsymbol{ε}}[(y(\boldsymbol{x})-\hat{f}(\boldsymbol{x})) ^ 2]=E _ {\boldsymbol{x}}[(f(\boldsymbol{x}) - E _ {\boldsymbol{ε}}[\hat{f}(\boldsymbol{x})]) ^ 2] + E _ {\boldsymbol{x}}E _ {\boldsymbol{ε}}[(\hat{f}(\boldsymbol{x}) - E _ {\boldsymbol{ε}}[\hat{f}(\boldsymbol{x})]) ^ 2] + σ ^ 2 $$

ここまでの解析で見てきたのはデータ $X$ で条件づけた下でのラベルノイズ $\boldsymbol{ε}$ についての分解を考え，後で $X$ について平均をとるものだった．これをUnderstanding Double Descent Requires A Fine-Grained Bias-Variance Decompositionの論文の言葉を借りて "Semi-classical Approach" と呼ぶことにする*3．

• Semi-classical Approach: $n,d\to\infty$ の極限では測度の集中の効果で $X$ の期待値がとれ，Fixed Design に一致する

  $$ E _ {\boldsymbol{x}}E _ XE _ {\boldsymbol{ε}}[(y(\boldsymbol{x})-\hat{f}(\boldsymbol{x})) ^ 2]=E _ {\boldsymbol{x}}E _ X[(f(\boldsymbol{x}) - E _ {\boldsymbol{ε}}[\hat{f}(\boldsymbol{x})]) ^ 2] + E _ {\boldsymbol{x}}E _ XE _ {\boldsymbol{ε}}[(\hat{f}(\boldsymbol{x}) - E _ {\boldsymbol{ε}}[\hat{f}(\boldsymbol{x})]) ^ 2] + σ ^ 2 $$

ここで考えたいのは，Random Design のバイアス・バリアンス分解と， Semi-classical / Fixed Design の分解の関係性はどうなっているか，また，どの変数がピークに寄与しているか解釈できるか？という点である．

以下の解析はUnderstanding Double Descent Requires A Fine-Grained Bias-Variance Decompositionを線形回帰の場合に焼き直したもので，Scaling and Renormalization in High-Dimensional Regressionも参考にしている．なお，こんな感じの細かい分解をしてはどうか？みたいな話はRethinking Bias-Variance Trade-off for Generalization of Neural Networksが発端っぽい？

線形回帰におけるより細かいバイアス・バリアンス分解

ここまでみてきた線形回帰の場合で，Random Design での分解をしてみる．先に $E _ {\boldsymbol{x}}$ をとっておけば，

$$ \begin{split} E _ {\boldsymbol{x}}E _ {X,\boldsymbol{ε}}[(y-\hat{y}) ^ 2] &= E _ {X,\boldsymbol{ε}}[\lVert \boldsymbol{\theta} - \hat{\boldsymbol{\theta}} \rVert ^ 2] + σ ^ 2 \cr &= \lVert \boldsymbol{\theta}-E _ {X,\boldsymbol{ε}}[\hat{\boldsymbol{\theta}}] \rVert ^ 2+E _ {X,\boldsymbol{ε}}[\lVert \hat{\boldsymbol{\theta}}-E _ {X,\boldsymbol{ε}}[\hat{\boldsymbol{\theta}}] \rVert ^ 2] +σ ^ 2 \cr &= B(\hat{\boldsymbol{\theta}}) + V _ {X,\boldsymbol{ε}}[\hat{\boldsymbol{\theta}}] +σ ^ 2 \end{split} $$

ここで

• バイアス（の 2 乗）

  $$ B(\hat{\boldsymbol{\theta}})=\lVert \boldsymbol{\theta}-E _ {X,\boldsymbol{ε}}[\hat{\boldsymbol{\theta}}] \rVert ^ 2 $$

• バリアンス

  $$ V _ {X,\boldsymbol{ε}}[\hat{\boldsymbol{\theta}}]=E _ {X,\boldsymbol{ε}}[\lVert \hat{\boldsymbol{\theta}}-E _ {X,\boldsymbol{ε}}[\hat{\boldsymbol{\theta}}] \rVert ^ 2] $$

• 除去できないノイズ $$ σ ^ 2 $$

とわかれる．各項を計算すると，$\hat{\boldsymbol{\theta}}=X ^ +X\boldsymbol{\theta}+X ^ +\boldsymbol{ε}$ と $E _ {X,\boldsymbol{ε}}[\hat{\boldsymbol{\theta}}]=E _ X[X ^ +X]\boldsymbol{\theta}+E _ X[X ^ +]E _ {\boldsymbol{ε}}[\boldsymbol{ε}]=E _ X[X ^ +X]\boldsymbol{\theta}$ を使うことで

$$ \begin{split} V _ {X,\boldsymbol{ε}}[\hat{\boldsymbol{\theta}}] &=E _ {X,\boldsymbol{ε}}[\lVert \hat{\boldsymbol{\theta}}-E _ {X,\boldsymbol{ε}}[\hat{\boldsymbol{\theta}}] \rVert ^ 2] \cr &=E _ {X,\boldsymbol{ε}}[\lVert (X ^ +X-E _ X[X ^ +X])\boldsymbol{\theta}+X ^ +\boldsymbol{ε} \rVert ^ 2] \cr &=E _ {X}[\lVert (X ^ +X-E _ X[X ^ +X])\boldsymbol{\theta} \rVert ^ 2]+E _ {X,\boldsymbol{ε}}[\lVert X ^ +\boldsymbol{ε} \rVert ^ 2] \cr &=E _ {X}[\lVert X ^ +X \boldsymbol{\theta} \rVert ^ 2]-\lVert E _ X[X ^ +X\boldsymbol{\theta}] \rVert ^ 2+σ ^ 2\operatorname{Tr}(E _ {X}[(X ^ +) ^ \top X ^ +]) \cr &=\begin{cases} σ ^ 2\frac{p}{n-p-1} & (p<n-1)\cr \infty &(p=n-1,n,n+1)\cr \lVert \boldsymbol{\theta} \rVert ^ 2\frac{n}{p}\left(1-\frac{n}{p}\right)+σ ^ 2\frac{n}{p-n-1} & (p>n+1) \end{cases} \end{split} $$

$$ \begin{split} B &:= \lVert \boldsymbol{\theta}-E _ {X,\boldsymbol{ε}}[\hat{\boldsymbol{\theta}}] \rVert ^ 2 \cr &= E _ {X,\boldsymbol{ε}}[\lVert \boldsymbol{\theta} - \hat{\boldsymbol{\theta}} \rVert ^ 2] - V _ {X,\boldsymbol{ε}}[\hat{\boldsymbol{\theta}}] \cr &=\begin{cases} 0& (p\leq n)\cr \lVert \boldsymbol{\theta} \rVert ^ 2\left(1-\frac{n}{p}\right) ^ 2 & (p>n) \end{cases} \end{split} $$

と計算できる．ここで

• $p\leq n$ の場合

  • $E _ {X}[\lVert X ^ +X \boldsymbol{\theta} \rVert ^ 2]-\lVert E _ X[X ^ +X\boldsymbol{\theta}] \rVert ^ 2=E _ {X}[\lVert \boldsymbol{\theta} \rVert ^ 2]-\lVert E _ X[\boldsymbol{\theta}] \rVert ^ 2=0$

  • $σ ^ 2\operatorname{Tr}( (X ^ \top X) ^ {-1})=σ ^ 2\frac{p}{n-p-1}$

• $p\geq n$ の場合

  • $E _ {X}[\lVert X ^ +X \boldsymbol{\theta} \rVert ^ 2]=E _ X[\boldsymbol{\theta} ^ \top(X ^ +X) ^ \top X ^ +X\boldsymbol{\theta}]=E[\boldsymbol{\theta} ^ \top X ^ +X \boldsymbol{\theta}]=\frac{1}{p}\operatorname{Tr}E[X ^ +X]\lVert \boldsymbol{\theta} \rVert ^ 2=\frac{n}{p}\lVert \boldsymbol{\theta} \rVert ^ 2$

  • $E _ {X}[X ^ +X \boldsymbol{\theta}]=\frac{n}{p}\boldsymbol{\theta}$ より $\lVert E _ X[X ^ +X\boldsymbol{\theta}] \rVert ^ 2=\frac{n ^ 2}{p ^ 2}\lVert \boldsymbol{\theta} \rVert ^ 2$

  • $σ ^ 2\operatorname{Tr}( (XX ^ \top) ^ {-1})=σ ^ 2\frac{n}{p-n-1}$

となることを利用した．

ここからはより細かくバリアンスの分解を行うことを考える．つまり，$X$ だけ， $\boldsymbol{ε}$ だけ，$X$ と $\boldsymbol{ε}$ の相互作用による寄与の項に分解する．

• $X$ だけからの寄与: $\boldsymbol{ε}$ について周辺化してから $X$ についての分散を考える

  $$ \begin{split} V _ X &:= V _ XE _ {\boldsymbol{ε}}[\hat{\boldsymbol{\theta}}] \cr &=V _ X[X ^ +X\boldsymbol{\theta}] \cr &=E _ {X}[\lVert (X ^ +X-E _ X[X ^ +X])\boldsymbol{\theta} \rVert ^ 2] \cr &=\begin{cases} 0 & (p\leq n)\cr \frac{n}{p}\left(1-\frac{n}{p}\right)\lVert \boldsymbol{\theta} \rVert ^ 2 & (p>n) \end{cases} \end{split} $$

• $\boldsymbol{ε}$ だけからの寄与: $X$ について周辺化してから $\boldsymbol{ε}$ についての分散を考える （素直にやると計算できないので，元の問題で $X=0$ だったと思うことにする）

  $$ V _ {\boldsymbol{ε}}:= V _ {\boldsymbol{ε}}E _ {X}[\hat{\boldsymbol{\theta}}]=0 $$

• $X$ と $\boldsymbol{ε}$ の相互作用による寄与: 全体から上の 2 つを引く

  $$ \begin{split} V _ {X,\boldsymbol{ε}}&:= V _ {X,\boldsymbol{ε}}[\hat{\boldsymbol{\theta}}] - V _ X - V _ {\boldsymbol{ε}}\cr &=E _ {X,\boldsymbol{ε}}[\lVert X ^ +\boldsymbol{ε} \rVert ^ 2] \cr &=σ ^ 2\operatorname{Tr}(E _ {X}[(X ^ +) ^ \top X ^ +]) \cr &=\begin{cases} σ ^ 2\frac{p}{n-p-1} & (p<n-1)\cr \infty & (p=n-1,n,n+1) \cr σ ^ 2\frac{n}{p-n-1} & (p>n+1) \end{cases} \end{split} $$

というわけで，以上により，汎化誤差は次のベン図のように分解できることがわかった．この結果をみると線形回帰における二重降下はデータと観測値のノイズの相互作用によって発生することがわかる．

また，"Semi-Classical Approach" では$B+V _ X$ をシグナル（バイアス）， $V _ {X,\boldsymbol{ε}}+V _ {\boldsymbol{ε}}$ をノイズ（バリアンス）と呼んでいたこともわかる（こういった背景があるのでわざとシグナルとノイズという言葉を使った）．

$$ \begin{array}{c||c:c:c:c} & B & V _ X & V _ {X,\boldsymbol{ε}} & V _ {\boldsymbol{ε}} \cr \hline γ<1 \vphantom{\displaystyle\lVert \boldsymbol{\theta} \rVert ^ 2\left(1-\frac{1}{γ}\right) ^ 2}& 0 & 0 & \displaystyleσ ^ 2\frac{γ}{1-γ} & 0 \cr γ>1 & \displaystyle\lVert \boldsymbol{\theta} \rVert ^ 2\left(1-\frac{1}{γ}\right) ^ 2 & \displaystyle\lVert \boldsymbol{\theta} \rVert ^ 2\frac{1}{γ}\left(1-\frac{1}{γ}\right) & \displaystyleσ ^ 2\frac{1}{γ-1} & 0 \end {array}\quad\left(\frac{p}{n}\toγ,\,n,p\to\infty\right) $$

Fine-Grained Bias-Variance Decomposition

汎化誤差のバイアス・バリアンス分解（$p$の関数として）

汎化誤差のバイアス・バリアンス分解（$γ$の関数として）

参考文献

• M.Lafon and A. Thomas, “Understanding the Double Descent Phenomenon in Deep Learning”, arXiv preprint arXiv:2403.10459, 2021.
• P.Nakkiran, G. Kaplun, Y. Bansal, T. Yang, B. Barak, and I. Sutskever, “Deep Double Descent: Where Bigger Models and More Data Hurt”, in International Conference on Learning Representations, 2020.
• OpenAI, “Deep Double Descent”.
• O.Shamir, “Deep Double Descent”.
• B.Adlam and J. Pennington, “Two Models of Double Descent for Weak Features”, arXiv preprint arXiv:1903.07571, 2019.
• M.Belkin, D. Hsu, S. Ma, and S. Mandal, “Surprises in High-Dimensional Ridgeless Least Squares Interpolation”, arXiv preprint arXiv:1903.08560, 2019.
• T.Hastie, A. Montanari, S. Rosset, and R. J. Tibshirani, “More Data Can Hurt for Linear Regression: Sample-wise Double Descent”, arXiv preprint arXiv:1912.07242, 2019.
• Z.Yang, Z. Zhang, M. Li, and J. Wang, “Rethinking Bias-Variance Trade-off for Generalization of Neural Networks”, arXiv preprint arXiv:2002.11328, 2020.
• S.D'Ascoli, L. Sagun, and G. Biroli, “Understanding Double Descent Requires a Fine-Grained Bias-Variance Decomposition”, in Advances in Neural Information Processing Systems, 2020.
• S.Mei and A. Montanari, “Scaling and Renormalization in High-Dimensional Regression”, arXiv preprint arXiv:2405.00592, 2024.
• S.Mei, “Lecture 16: Double Descent and Generalized Linear Models”.
• G.Livan, M. Novaes, and P. Vivo, A First Course in Random Matrix Theory: For Physicists, Engineers and Data Scientists. Cambridge University Press, 2018.
• Z.Bai and J. W. Silverstein, Spectral Analysis of Large Dimensional Random Matrices. Springer, 2010.
• 清水泰隆. (2023). 統計学への漸近論，その先は 現代の統計リテラシーから確率過程の統計学へ. 内田老鶴圃.

www.saiensu.co.jp

*1:ちなみに Regularization-Wise Double Descent なんてのもあったりする

*2:収束の速度が $λ _ {\mathrm{min}}(X)$ で決まることもわかる．

*3:ここでは線形回帰の場合を見越して $X$ だが，モデル設定にかかわる任意のランダム性のようなものを考えている．具体的には初期パラメータ，訓練のミニバッチサイズの選択などである