1/48

近似確率伝搬法による⼀般化線形モデルの予測誤差表現について

統計数理研究所
坂⽥

数理・推論研究系
綾⾹

Special thanks to Yukito Iba

2/48

今⽇の内容
l 予測誤差(の推定量)を評価する
動機• 予測誤差を何のために評価するのか

⽬標• 予測誤差の推定量の解釈と⽐較• 効率的な計算法の開発

⼿法• ファクターグラフ表現を⽤いた近似推論法▫ ファクターグラフ表現…確率分布をグラフで表す▫ 近似推論法…確率伝搬法

結果• 確率伝搬法を使うことで出来ること，もたらされる解釈について説明

3/48

予測誤差とは予測に基づく統計的モデリング4/48

統計学におけるモデリングn 与えられたデータを説明するモデルを作る
l モデル：確率分布l ⼀般に，データが⽣成される過程は未知Truth Realizationn 可能な確率分布のうち，どれが適切なのか？
“射影”

l ⽤意したモデル候補から適切なものを選ぶ
l これをモデル選択というモデル空間

n データは⼀つの実現値である
l 実現値をもとに，確率分布の性質を知ることが⽬標”The Elements ofStatistical Learning”
より

5/48

モデル選択
(例) データ ! = {$!, … , $" }を得たとする
l パラメータ(をもつ確率分布) * ( を導⼊する
l データ!のもとで，パラメータ(を推定する → (+ ! とする• (例) 最尤推定"(+#$ ! = max / (; ! ,%9! !
l モデル候補ℳ = )! *|6
適切なモデルを選ぶ/((; !) = 3 ln ) * = $& 6&'!

9( !
, )( *|69) !
, … , )) *|6の中から

6/48

モデル選択の例
2次関数

3次関数

どのモデルでデータを記述することが適切か？
客観的な指標が必要：予測性能

5次関数

7/48

「予測」に基づくモデリングl 与えられたデータに対するあてはまりの良さは必要
l ⼀⽅で，与えられたデータに過適合してしまうことも問題• 将来の予測ができなくなる可能性がある「予測」に基づくモデリングの⽴場
n “真のモデル”から⽣成される 新しいデータについて⾼い予測能⼒を持つモデルを良いモデルとする
p “真のモデル”を⾒つけることが⽬的ではない．
• 結果的に“真のモデル”が得られるのであれば，それはそれで良いこと．

p 予測能⼒を予測誤差で定量化

8/48

予測誤差とはIn-sample error

Extra-sample error：Training data
：Estimatedrelationship

OutputOutput：Test data

Input

Input

訓練データと異なる⼊出⼒関係をテストデータとする

訓練データと同じ⼊⼒に対する異なる出⼒をテストデータとする

推定量：交差検証誤差

推定量：情報量規準, Cp規準
漸近等価 [Stone (1977)]

9/48

予測誤差の定義と汎化ギャップ
以下では回帰問題を考える．
• 仮定する確率分布を!(⋅ |%)とする．% !, # =!
"()はパラメータ．

• 回帰係数：) ∈ ℝ"
• 説明変数：( ∈ ℝ#×"
• 出⼒：, (-次元ベクトル)

データセット! = :, ;l データ!のもとでの推定値を<= . と表現する．
1< !
• Training error： err*+,-. ! = − ln C :|6 ;, =B

予測誤差と訓練誤差の差を汎化ギャップと呼ぶ• Ex-sample
02< !
err/+0 ! = −D3,5!"# ln C *|6 F.06 , =prediction error：

1• In-sample-.
< !
err/+0 ! = − D1 ln C E|6 ;, =prediction error：B10/48

予測誤差の推定量l 予測誤差を定義通りに評価するには，真の確率分布が必要
l 真の確率分布は未知なので厳密評価は不可能→ 予測誤差の推定量を構成する• 確率変数の実現値を⼊れると，それに応じた値を返す関数

推定量の例• 交差検証誤差• 情報量規準• Cp規準

11/48

交差検証誤差(cross validation error)• ここではleave-one-out CV (LOOCV) errorを考える

l G番⽬のサンプルを抜いたデータを!∖& とする．
l !∖& のもとでの推定値をIH .∖& とする．
n LOOCV errorの定義は次の通り
"11
< !∖&
err$889: ! = − 3 ln C J& |
F& =B
K&'!

n 利点：さまざまなモデルに適⽤可能
n ⽋点：計算量が多い• Linear estimatorの場合は解析的に表現できる

12/48Linear estimatorl Linear estimator< = L:として与えられる / ∈ ℝ#×#
モデルによる出⼒の表現<:が :(例) 線形回帰問題における最尤推定 (- > 1)

< = argmin : −=;1 ' (! 1 '2=)
( (( ,,11)!

*"• 特に，⾏列 5 = ('( ((!

1K(;=(11 ' (! '2=,(2)= (( (( ,1111(' をHat matrixと呼ぶ [Hoaglin & Welsch (1978)]．

13/48Linear estimatorの場合のLOOCV error
l Linear estimatorの場合，leave-one-out sample !∖& のもとでの出⼒表現は次のように与えられるJH& !∖&

JH& ! − O&& J&=
1 − O&&[Cook (1977), Seber & Lee (2003)]

l 特にガウス分布の場合，LOOCV errorは次の通り"

J& − JH& !
1err$889: ! =32B1 − O&& !
&'!

n フルデータ!のもと１度だけ推定すれば良い
n Linear estimatorのみ適⽤可能(

14/48

Importance sampling cross validation errorl 正則化つき最尤法<< ! = argmax Q< !, = ,=Q< !, = = ln C :|6 ;, =;− ℎ =; Sl ベイズ推定(最⼤事後確率推定)で表す• 事後分布 6+,-. ) .; 8 ∝ ! / ,|% (, )exp 8ℎ ); >• 推定値の表現 @?0,2 . = lim DE@ @2 6+,-. ) .; 8/→4

l 対応するLOOCV error#1
err56678 .; 8 = − H ln DE@ ! / J& K& L& , ) 6+,-. )|.∖& ; 8&9!

正則化15/48

Importance sampling cross validation error#

1err56678 .; 8 = − H ln DE@ ! / J& K& L& , ) 6+,-. )|.∖& ; 8&9!
#

1= − H ln DE@ ! J& %& L& , )&9!

∏:;& ! / J: K: L: , ) exp 8ℎ ); >∫ E@ ∏:;& ! / J: K: L: , ) exp 8ℎ ); >#∫ E@ ∏:;& ! / J: K: L: , ) exp 8ℎ ); >1
= H ln/∫ E@ ∏#:9! ! J: K: L: , ) exp 8ℎ ); >&9!
#

1
1= H ln D E@6+,-. )|.; 8 /! J& K& L& , )&9!

この表現をImportance Sampling Cross Validation error と呼ぶ．
[Vehtari & Lampinen, Neural Computation (2002) ]16/48

Widely Applicable Information Criterion#1</J& K& L& , )• 次の関数を定義： O P = H ln DE@6+,-. )|.; 8 !
&9!

→ −O 1 ：Training error，O −1 ：ISCVl T ≪ 1での表現を適⽤した汎化ギャップの推定量 (Functional Variance)"

9 ! =ΔW
3 X/=>* ln C J& |(& C& , I
B&'!

Q+,-. ⋅ は事後分布における分散• 予測誤差の推定量「訓練誤差+FV」をWidely applicable information criterionとよぶ．
[Watanabe, JMLR (2010) & JJSDS (2021)]

n 利点：計算量が少ない
n ⽋点：収束半径を考慮する必要がある17/48Linear estimatorのin-sample error?@err+=>

12 .
. = A TB U − ,2S -12 .
err.=C?@ . = A,−,
2S -A
A

#Δ ?@AA

1. = A H CovD J& , J?& .
S &9!

• Linear estimatorの場合： Cov? J& , JH& !

[Efron (2004)]= Cov? J& , ∑@ O&@ J@• さらにJ ∼ ^(G, _ (`" )のときCov? J& , ∑@ O&@ J@ = O&& なのでΔ -.

1! = ( Tr O_ B

…Degrees of freedomと呼ばれる線形回帰問題では，Hat matrixの性質よりTr b = K =モデルパラメータ数→ Akaike Information criterionと⼀致する

18/48Beyond Linear estimator?@err+=>

12 .
. = A TB U − ,2S -12 .
err.=C?@ . = A,−,2S -A
A#Δ ?@AA

1. = A H CovD J& , J?& .
S &9!

[Efron (2004)]• , = Z + S A \, \ ∼ ^(0, `# )のとき

1
Steinʼs Lemma： ( Cov? J& , JH& !
_c= D?
JH& !
cd&#a ?@ . =Δ

1bHJ?& .
bc&

Generalized degrees of freedomと呼ばれる&9!

n 利点：⼀般の推定量に適⽤可能
n ⽋点：数値的微分が必要．ガウスノイズの場合のみ不偏19/48

いろいろな予測誤差の推定量がある必要とされること(1) 効率的な計算⽅法
l 近似評価法が必要(2) 推定量の間の関係性の理解
l 特に⾮漸近領域において• ファクターグラフ表現• 近似確率伝搬法を使って議論する(3) 実⽤上の指針
l 近似評価法の精度，適⽤限界

20/48

問題設定

21/48⼀般化線形モデル(GLM)• 回帰係数：= ∈ ℝA
• 説明変数：; ∈ ℝ"×A
• 出⼒：: (B次元ベクトル)

データセット! = :, ;l 対数尤度 ln C :|6 = 6C: − g 6 + i :• K& =

1"#H d&2 @2 ,1 29!

#e % = H e K& ,f , = H f J&&9!

be(K& ),& =bK&• TJ

T J& − T J&!

!

AA&9!
Abe K&=bK&A

p 線形回帰：e K = K A , f J = − J A
p ロジスティック回帰： e K = ln 1 + g E , f J = 0

22/48

制約付き最尤推定<< ! = argmax Q< !, = ,=;Q< !, = = ln C :|6 ;, =l 正則化パラメータSを選ぶことがモデル選択に対応する
2! . , !(i|% (, )2A . , …
• モデル候補：ℳ = !(i|% (, )l 定式化：最⼤事後確率を⽤いた推定として表現• 事後分布 6+,-. ) .; 8 ∝ ! / ,|% (, )exp 8ℎ ); >• 推定値の表現 @?0,2 . = lim DE@ @2 6+,-. ) .; 8/→4n 尤度⟷ハミルトニアン，事前分布⟷外場
n 最⼤事後確率推定⟷基底状態探索− ℎ =; S正則化

23/48

近似確率伝搬法
Generalized Approximate Message Passingとは

24/48Generalized Approximate Message Passingとはマルコフネットワーク(特にファクターグラフ)上で定義されるアルゴリズム
l⽬的：変数の周辺化（分配関数や周辺化分布の評価）
n Message Passing• グラフ構造が規定する条件付き独⽴性のもと，効率的に周辺化• “メッセージ”の更新により計算

n Approximate Message Passing• メッセージの⼀部をガウス近似したmessage passingn Generalized Approximate Message Passing• ⼀般の尤度と事前分布に適⽤可能なAMP

25/48ファクターグラフ表現の例•2-body spin-glass model$ % ∝ exp * + ,!# -! -#!,#≡ / 0!# -! , -# ; ,!#!,#_D,!#•3-body spin-glass modelwith external field$ % ∝ exp * + ,!#$ -! -# -$!,#,$≡ / 0!#$ -! , -# , -$ ; ,!#$!,#,$_El Message Passingの⼿続き_D,!#$_E_F• ファクターkをかける• 変数の和を実⾏する: 変数ノード: ファクターノード

26/46

ファクターグラフ上でのSum-Productl 例: n lG o = ∑J$ ⋯ ∑J% n q|o を計算する
input messagelG

32 (→' 4'

32 %→' 4'm!

m(lIlHl(n q|o ∝ k! l!, l(, lG; m! k( lH, lI, lG; m(n lG o = 3 k( lH, lI, lG; m(ts (→G lG

J& ,J$l!

3 k! l!, l(, lG; m!
J% ,J'n 周辺化分布はinput messageの積により与えられる．
詳しくは→ts !→G lG

27/46

ファクターグラフ上でのSum-Productl 例: n lI o = ∑J& ∑J' ∑J% n q|o を計算する3)→) 4)

lImH

lHm!

l!

m(l(n lI o = 3 kH lH, lI; mH 3 k( l(, lH; m( 3 k! l!, lH; m!
J&J'

J%t
s (→H lHts !→H lH
tH→H lHts H→I lI• 3!→* 4! : 5番⽬の変数ノードから6番⽬のファクターノードへのoutput message28/46

Message PassingのGLMへの適⽤• 事後分布 6+,-. ) .; 8 ∝ ! / ,|% (, )exp 8ℎ ); >• 推定値の表現 @?0,2 . = lim DE@ @2 6+,-. ) .; 8/→4

J!

J(JHJI• グラフの構造はI!
説明変数から決まる

I(IHIIIG

ILu!

u(uHuIuG

uL

ファクター＝尤度! / J! K! L! , )Montanari, “Graphical Models Concepts in Compressed Sensing”
arXiv:1011.4328

29/46

Output messageの意味• ファクターグラフがツリーであるとする．

7)• 変数ノードkとファクターノードe間のエッジを切り離したとする．

7(• このとき変数ノードkが含まれるツリーをl k ∖ e とする．

0+

Output message tD→M ID は7%v x ∖ g におけるID の厳密な周辺化分布である．

7.

3%→+ 7%0,0-

Output message tD→M ID は

7'

7/LOO sample !∖M のもとでのID の厳密な周辺化分布である．

v 1∖g70• @?2→F ≡ lim ∫ E@2 @2 o2→F @2 は.∖F の下での推定値/→4[Mezard & Montanari, (2009)]

30/46

ファクターグラフによるLOOサンプル下での推定の表現@!
@G@A
Leave-one-out
@!
@G

J&@A

メッセージを利⽤すれば，LOO sampleのもとでM回の推定を⾏わずに，LOOCV errorの近似値を得ることが出来る

31/46

ツリーではない場合
l ツリーでないグラフ（ループを含むグラフ）にMessage Passingを適⽤することも可能
その場合，グラフがツリーであると仮定したうえで，真の分布とのKullback-Leibler divergenceを最⼩化するような“試⾏”周辺化分布を求めている と解釈できる• Bethe-Peierls近似に相当

l 厳密性の保証はないが，近似計算法として妥当な精度が得られる場合がある．
• ターボ符号，低密度パリティ検査符号など [Gallager (1962), MacKay (1999)]

32/46

Message PassingからApproximate Message Passingへ
l Input message
op &→2 @2 ∝ DE)-H,+> & ∖2 ! / J& |K& L& , )oI→& @2I∈-H,+> & ∖2

ガウス近似

8Aop &→2 @2 ∝ exp −@2 − @u&→22t&→2
̃l Output message
o2→& @2 ∝ r2 @2 ; >qq<∈ℱ 2 ∖&op <→2 @2※ @u&→2 は平均， 8 (! t&→2̃ は分散
8A∝ r2 @2 ; > exp − H@2 − @u<→22t̃<→2<∈ℱ 2

平均と分散を評価@?2→& = DE@2 o2→& @2 ,t2→& = 8 DE@2 o2→& @2 @2 − @?2→&

A@u&→2 とt&→2
̃ ，@?2→& とt2→& の更新アルゴリズムをApproximate Message Passingと呼ぶ．

33/46

AMPの使⽤が妥当な場合
次の性質が満たされるとき，ガウス近似は妥当である．
i.

変数の数Kが⼗分⼤きいii. 説明変数の各成分がy(K N!/()• 説明変数が標準化(standardization)されている

iii. 説明変数間の相関が無視できるほど⼩さい• 多重共線性がないように説明変数を設計

iv. 説明変数がスパースではない
i〜ivの仮定のもと，中⼼極限定理からAMPの妥当性は保証される．

34/46

これから⽰す結果• GAMPの適⽤可能な範囲で，⼀般化線形モデルにおける次の結果を紹介
n ⼀般化⾃由度は熱揺らぎから評価可能
n Linear estimation ruleにおけるLOOCV error表現の拡張を実現する
n WAIC導出の根拠となる，ISCVの重みに関する級数展開は線形ガウスモデルにおいて スピングラス転移点で発散する• ⼀般化線形モデルで収束半径を導出することも可能

l GAMPの根拠となる仮定が破れているときについては最後に少し説明

35/46

GAMPによる予測誤差の表現

36/46

ファクターグラフによるLOOサンプル下での推定の表現@!
@G@A
Leave-one-out
@!
@G

J&@A

J& がないグラフでの@! の期待値：@?!→&(メッセージo!→&の期待値)

37/46

LOOCV error• ⼀般化線形モデルにおける!
∖&LOO sampleのもとでの推定値(+& = A ∑A
HD→& とD'! z&D I!
+フルデータのもとでの推定値(& = A ∑A
HD の関係D'! z&D I∖&
2(+& − (+&1= J& − gP (+& , 8* = + <*!( =!→*;{&!1%•代⼊#1∖&
err56678 . = − H ln ! J& |Kv&&9!

特にガウス分布の場合> = @?での表現
Linear estimator ?

∖&J& − J?& = J& − J?& 1 + w&

J& − J?& .∖&J& − J?& .
=1 − x&&

Linear estimation ruleにおける表現の⼀般化{&O&& ⟷= lim W
1 + {& Q→R(&( %((− (& %(※ ⋅ 3! ：Scope 6 の局所平均38/46

ファクターグラフによるin-sample errorの表現
@!
@G

J&@A

新しいデータを取得
@!
@G

J&P
@A

39/46

GAMPによる⼀般化⾃由度の推定量
E* → E*4 の変化に伴う@!
@G

汎化ギャップ∝ Scope[G]の “熱揺らぎ”

J&@!

@A

@G

J&汎化ギャップ

@A
@!

@G

J&P

9 -. = lim
Δ@A

Q→R"(_ W

B3&'!

(&( %((− (& %(※ ⋅ 3! ：Scope 6 の局所平均

a ?@ = Tr(x)の拡張• Linear estimatorにおけるΔ「予測における揺らぎと応答」 [Iba (private communication), Watanabe JJSDS (2021)]

40/46

ファクターグラフによるISCVの表現
T=1

WAICの導出：@!
@GT ≪ 1での汎化ギャップを|T| = 1に適⽤

J&@AT ≪1@!
@GJ&@A

T = −1@!
@G@A41/46

GAMPによるWAICの表現と適⽤限界
T=1@!
@G

J&

汎化ギャップは， P ≪ 1のときScope[}]の “熱揺らぎ”と残差から評価可能

@A@!
@G

T ≪1@GJ&@!
@A

T = −1

J&@A#8A
H J& − J?&/→4 -a = lim
Δ&9!

@!
@G
@A

K&A

E!

− K&AE!

• S A を推定量(残差2)で表した⼀般化⾃由度※ - ( を推定量で置き換える⽅法の提案 [Efron (2004)]

42/46

GAMPによるWAICの表現と適⽤限界
T=1@!
@GJ&@AT ≪1@GJ&@!
@A

T = −1@!
@G@A

尤度がガウス分布の場合，ISCVの級数展開は次の条件が成⽴するとき収束する
511P <,
8=+ 8*w
G*1%
1+ww
<1
H = 1で発散しないためには1+w
これはスピングラス転移点に対応する．
意味：
⼀般化⾃由度<データ数の場合に適⽤可能※ 正則化による推定値の縮⼩がある場合は必ずしもモデルパラメータ>データ数で発散するわけではない43/46

実際の値とGAMPの⽐較：Ridge回帰の場合
LOOCV errorの汎化ギャップ
S = 0.01, X = 1

In-sample errorの汎化ギャップ
GAMPによるLOOCV errorの汎化ギャップ
GAMPによるIn-sample errorの汎化ギャップ＊GAMPによるFV• GAMPは実際の値を精度良く評価できている．
• 漸近極限ではすべての推定量は⼀致するが，データサンプル数が⼩さくなるにつれてずれる様⼦が⾒える．
データサンプル数/パラメータ数• Logistic回帰などでも同様．

44/46GAMPにより得られたこと(1) 汎化ギャップの効率的な計算
l 汎化ギャップを評価するための追加の計算が必要ない
n GAMPによる表現は，Linear estimatorに対する表現の拡張と理解できる(2) 「揺らぎと応答」としての汎化ギャップの表現
l J& → J&L による汎化ギャップ ∝ Scope } に対する事後分布の分散n 未知データへの応答は，事後分布の性質から把握できる(3) ISCVの吸数展開の収束半径l Gaussian likelihoodの場合はスピングラス転移が起こらなければWAICは妥当
n そのほかのlikelihoodの場合にも収束半径は導出できる
Sakata, J Phys A (2023), (arXiv:2206.12832)

45/46

実データへの適⽤について• 実データにおいては，GAMPの仮定は必ずしも満たされない
TIa supernovae dataset, a=4 相転移点

l 説明変数に相関がある場合の補正
0.2• Obuchi & Kabashima, JSTAT (2016)

n GAMPを利⽤した“相転移点”の発⾒とモデル候補からの除外• Obuchi & Sakata, JSTAT (2019)

CV errors

n 分散共分散⾏列を正確に評価する⽅法

Approximate
LiteralInstability
CVE min

0.15
0.1
0.05
010-1

LOOCV errorの最⼩値

100

正則化
パラメータ

最も良いモデル46/46

おわりに
n 統計学・機械学習におけるモデリングと物理学におけるモデリングの考え⽅の違い• Ground truthがない問題に折り合いをつけるための客観的⽅法
n その⽅法の中にも“物理”がある• モデルを作る＝ハミルトニアンの設計• 正則化つき最尤推定＝外場のある問題での基底状態探索• 推定アルゴリズム＝ダイナミクス• 汎化ギャップ＝相互作⽤，温度に関する摂動への応答

n 推論⽅法の開発，計算の⾼度化 → 統計学・機械学習の深化へ• 特に統計⼒学は⾮漸近領域の記述が得意