logo

機械学習における線形回帰モデルの最大事後確率推定 📂機械学習

機械学習における線形回帰モデルの最大事後確率推定

定理

データxi∈Rn\mathbf{x}_{i} \in \mathbb{R}^{n}とそのラベルyi∈Ry_{i} \in \mathbb{R}の関係が以下の線形モデルであると仮定する。

yi=wTxi+ϵi,i=1,…,K(1) y_{i} = \mathbf{w}^{\mathsf{T}} \mathbf{x}_{i} + \epsilon_{i}, \qquad i = 1, \ldots, K \tag{1}

事後確率が最大となるパラメータwMAP\mathbf{w}_{\text{MAP}}は次の通りである。y=[y1⋯yK]T\mathbf{y} = \begin{bmatrix} y_{1} & \cdots & y_{K} \end{bmatrix}^{\mathsf{T}}とX=[x1⋯xK]T∈Rn×K\mathbf{X} = \begin{bmatrix} \mathbf{x}_{1} & \cdots & \mathbf{x}_{K} \end{bmatrix}^{T} \in \mathbb{R}^{n \times K}について、

  • 事前分布が正規分布の場合: wMAP=(1σ2XTX+Σ−1)−1(1σ2XTy+Σ−1μ) \mathbf{w}_{\text{MAP}} = \left(\dfrac{1}{\sigma^{2}} \mathbf{X}^{T}\mathbf{X} + \boldsymbol{\Sigma}^{-1} \right)^{-1} \left(\dfrac{1}{\sigma^{2}} \mathbf{X}^{T} \mathbf{y} + \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu} \right)

    このとき、μ\boldsymbol{\mu}とΣ\boldsymbol{\Sigma}はそれぞれw\mathbf{w}の事前分布の平均ベクトルと共分散行列である。

  • 事前分布がラプラス分布の場合:

    最適解の明示的な形は存在しない。

説明

以下の(3)(3)を見れば、w\mathbf{w}の事前分布として標準正規分布N(0,I)N(\mathbf{0}, \mathbf{I})を仮定すると、リッジ回帰と同じ問題であることが分かる。

arg min⁡w[12σ2∥y−Xw∥22+12∥w∥22] \argmin_{\mathbf{w}} \left[\dfrac{1}{2\sigma^{2}} \| \mathbf{y} - \mathbf{X}\mathbf{w} \|_{2}^{2} + \dfrac{1}{2} \| \mathbf{w} \|_{2}^{2} \right]

一方、事前分布として🔒(25/05/16)ラプラス分布Laplace⁡(0,b)\operatorname{Laplace}(0, b)を仮定すると、ラッソ回帰と同じ問題である。

arg min⁡w[12σ2∥y−Xw∥22+1b∥w∥1] \argmin_{\mathbf{w}} \left[\dfrac{1}{2\sigma^{2}} \| \mathbf{y} - \mathbf{X}\mathbf{w} \|_{2}^{2} + \dfrac{1}{b} \| \mathbf{w} \|_{1} \right]

事前分布=正規分布

(1)(1)でw∈Rn\mathbf{w} \in \mathbb{R}^{n}は母数パラメータであり、ϵi∼N(0,σ2)\epsilon_{i} \sim N(0, \sigma^{2})をガウシアンノイズと仮定する。ϵi\epsilon_{i}がN(0,σ2)N(0, \sigma^{2})に従うと仮定したので、yi=wTxi+ϵiy_{i} = \mathbf{w}^{\mathsf{T}} \mathbf{x}_{i} + \epsilon_{i}はN(wTxi,σ2)N(\mathbf{w}^{\mathsf{T}} \mathbf{x}_{i}, \sigma^{2})に従う。

yi∼N(wTxi,σ2) y_{i} \sim N(\mathbf{w}^{\mathsf{T}} \mathbf{x}_{i}, \sigma^{2})

最大事後確率推定は次を満足するwMAP\mathbf{w}_{\text{MAP}}を求めるものである。

wMAP=arg max⁡wp(y∣w,X)p(w)(2) \mathbf{w}_{\text{MAP}} = \argmax_{\mathbf{w}} p(\mathbf{y} | \mathbf{w}, \mathbf{X}) p(\mathbf{w}) \tag{2}

p(y∣w,X)p(\mathbf{y} | \mathbf{w}, \mathbf{X})は尤度であり、p(w)p(\mathbf{w})は事前確率である。尤度関数は次の通りである。

p(y∣w,X)=1(2πσ2)K/2exp⁡[−12σ2∥y−Xw∥22] p(\mathbf{y} | \mathbf{w}, \mathbf{X}) = \dfrac{1}{(2\pi \sigma^{2})^{K/2}} \exp \left[ -\dfrac{1}{2\sigma^{2}} \| \mathbf{y} - \mathbf{X}\mathbf{w} \|_{2}^{2} \right]

また、w\mathbf{w}の事前分布が以下のような多変量正規分布に従うと仮定する。

w∼N(μ,Σ),p(w)=1(2π)ndet⁡Σexp⁡[−12(w−μ)TΣ−1(w−μ)] \mathbf{w} \sim N(\boldsymbol{\mu}, \boldsymbol{\Sigma}), \qquad p(\mathbf{w}) = \dfrac{1}{\sqrt{(2\pi)^{n} \det \boldsymbol{\Sigma}}} \exp \left[ -\dfrac{1}{2} (\mathbf{w} - \boldsymbol{\mu})^{\mathsf{T}} \boldsymbol{\Sigma}^{-1} (\mathbf{w} - \boldsymbol{\mu}) \right]

事後確率が指数関数で表現されるため、対数尤度を考慮することが計算上便利である。

wMAP=arg max⁡wlog⁡(p(y∣w,X)p(w))=arg max⁡w[−12σ2∥y−Xw∥22−12(w−μ)TΣ−1(w−μ)]=arg min⁡w[12σ2∥y−Xw∥22+12(w−μ)TΣ−1(w−μ)] \begin{align*} \mathbf{w}_{\text{MAP}} &= \argmax_{\mathbf{w}} \log (p(\mathbf{y} | \mathbf{w}, \mathbf{X}) p(\mathbf{w})) \\ &= \argmax_{\mathbf{w}} \left[-\dfrac{1}{2\sigma^{2}} \| \mathbf{y} - \mathbf{X}\mathbf{w} \|_{2}^{2} -\dfrac{1}{2} (\mathbf{w} - \boldsymbol{\mu})^{\mathsf{T}} \boldsymbol{\Sigma}^{-1} (\mathbf{w} - \boldsymbol{\mu}) \right] \\ &= \argmin_{\mathbf{w}} \left[\dfrac{1}{2\sigma^{2}} \| \mathbf{y} - \mathbf{X}\mathbf{w} \|_{2}^{2} + \dfrac{1}{2} (\mathbf{w} - \boldsymbol{\mu})^{\mathsf{T}} \boldsymbol{\Sigma}^{-1} (\mathbf{w} - \boldsymbol{\mu}) \right] \tag{3}\\ \end{align*}

したがって、上式を微分して、0\mathbf{0}となるようなw\mathbf{w}がwMAP\mathbf{w}_{\text{MAP}}である。勾配を計算してみよう。

∇w[12σ2∥y−Xw∥22+12(w−μ)TΣ−1(w−μ)]=∇w[12σ2(y−Xw)T(y−Xw)+12(w−μ)TΣ−1(w−μ)]=∇w[12σ2(y−Xw)T(y−Xw)]+∇w[12(w−μ)TΣ−1(w−μ)] \begin{align*} & \nabla_{\mathbf{w}} \left[\dfrac{1}{2\sigma^{2}} \| \mathbf{y} - \mathbf{X}\mathbf{w} \|_{2}^{2} +\dfrac{1}{2} (\mathbf{w} - \boldsymbol{\mu})^{\mathsf{T}} \boldsymbol{\Sigma}^{-1} (\mathbf{w} - \boldsymbol{\mu}) \right] \\ &= \nabla_{\mathbf{w}} \left[\dfrac{1}{2\sigma^{2}} (\mathbf{y} - \mathbf{X}\mathbf{w})^{T}(\mathbf{y} - \mathbf{X}\mathbf{w}) +\dfrac{1}{2} (\mathbf{w} - \boldsymbol{\mu})^{\mathsf{T}} \boldsymbol{\Sigma}^{-1} (\mathbf{w} - \boldsymbol{\mu}) \right] \\ &= \nabla_{\mathbf{w}} \left[\dfrac{1}{2\sigma^{2}} (\mathbf{y} - \mathbf{X}\mathbf{w})^{T}(\mathbf{y} - \mathbf{X}\mathbf{w})\right] + \nabla_{\mathbf{w}} \left[\dfrac{1}{2} (\mathbf{w} - \boldsymbol{\mu})^{\mathsf{T}} \boldsymbol{\Sigma}^{-1} (\mathbf{w} - \boldsymbol{\mu}) \right] \\ \end{align*}

微分の規則は以下を参照しよう。

ベクトルと行列の導関数

  • 内積 ∂f(x)∂x=∂(wTx)∂x=∂(xTw)∂x=w \frac{ \partial f(\mathbf{x})}{ \partial \mathbf{x} } = \frac{ \partial (\mathbf{w}^{T}\mathbf{x})}{ \partial \mathbf{x} } = \frac{ \partial (\mathbf{x}^{T}\mathbf{w})}{ \partial \mathbf{x} } = \mathbf{w}

  • ノルム ∇f(x)=∂∥x∥2∂x=∂(xTx)∂x=2x \nabla f(\mathbf{x}) = \dfrac{\partial \left\| \mathbf{x} \right\|^{2}}{\partial \mathbf{x}} = \dfrac{\partial (\mathbf{x}^{T}\mathbf{x})}{\partial \mathbf{x}} = 2\mathbf{x}

  • 二次形式 ∂f(x)∂x=∂(xRx)∂x=(R+RT)x \dfrac{\partial f(\mathbf{x})}{\partial \mathbf{x}} = \dfrac{\partial (\mathbf{x}\mathbf{R}\mathbf{x})}{\partial \mathbf{x}} = (\mathbf{R} + \mathbf{R}^{T})\mathbf{x} R\mathbf{R}が対称行列であれば、 ∂f(x)∂x=2Rx \dfrac{\partial f(\mathbf{x})}{\partial \mathbf{x}} = 2\mathbf{R}\mathbf{x}

最初の項から微分を計算すると次のようになる。

∇w[12σ2(y−Xw)T(y−Xw)]=∇w[12σ2(yTy−2yTXw+wTXTXw)]=12σ2[−2XTy+2XTXw]=1σ2[−XTy+XTXw]=−1σ2XT(y−Xw) \begin{align*} &\nabla_{\mathbf{w}} \left[\dfrac{1}{2\sigma^{2}} (\mathbf{y} - \mathbf{X}\mathbf{w})^{T}(\mathbf{y} - \mathbf{X}\mathbf{w})\right] \\ &=\nabla_{\mathbf{w}} \left[\dfrac{1}{2\sigma^{2}} \left( \mathbf{y}^{T}\mathbf{y} - 2\mathbf{y}^{T}\mathbf{X}\mathbf{w} + \mathbf{w}^{T}\mathbf{X}^{T}\mathbf{X}\mathbf{w} \right) \right] \\ &= \dfrac{1}{2\sigma^{2}}\left[- 2\mathbf{X}^{T}\mathbf{y} + 2\mathbf{X}^{T}\mathbf{X}\mathbf{w} \right] \\ &= \dfrac{1}{\sigma^{2}}\left[-\mathbf{X}^{T}\mathbf{y} + \mathbf{X}^{T}\mathbf{X}\mathbf{w} \right] \\ &= -\dfrac{1}{\sigma^{2}} \mathbf{X}^{T} \left(\mathbf{y} - \mathbf{X}\mathbf{w} \right) \\ \end{align*}

次に第二項を計算すると以下のようになる。共分散行列は対称行列であり、対称行列の逆行列も対称行列であることから

∇w[12(w−μ)TΣ−1(w−μ)]=Σ−1(w−μ) \nabla_{\mathbf{w}} \left[\dfrac{1}{2} (\mathbf{w} - \boldsymbol{\mu})^{\mathsf{T}} \boldsymbol{\Sigma}^{-1} (\mathbf{w} - \boldsymbol{\mu}) \right] \\ = \boldsymbol{\Sigma}^{-1}(\mathbf{w} - \boldsymbol{\mu})

したがって次を得る。

−1σ2XT(y−XwMAP)+Σ−1(wMAP−μ)=0 -\dfrac{1}{\sigma^{2}} \mathbf{X}^{T} \left(\mathbf{y} - \mathbf{X}\mathbf{w}_{\text{MAP}} \right) + \boldsymbol{\Sigma}^{-1}(\mathbf{w}_{\text{MAP}} - \boldsymbol{\mu}) = \mathbf{0}

wMAP\mathbf{w}_{\text{MAP}}を求めるために上式を解くと次のようになる。

−1σ2XTy+1σ2XTXwMAP+Σ−1wMAP−Σ−1μ=0  ⟹  (1σ2XTX+Σ−1)wMAP=1σ2XTy+Σ−1μ  ⟹  wMAP=(1σ2XTX+Σ−1)−1(1σ2XTy+Σ−1μ) \begin{align*} && - \dfrac{1}{\sigma^{2}} \mathbf{X}^{T} \mathbf{y} + \dfrac{1}{\sigma^{2}} \mathbf{X}^{T}\mathbf{X}\mathbf{w}_{\text{MAP}} + \boldsymbol{\Sigma}^{-1} \mathbf{w}_{\text{MAP}} - \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu} = \mathbf{0} \\ \implies && \left(\dfrac{1}{\sigma^{2}} \mathbf{X}^{T}\mathbf{X} + \boldsymbol{\Sigma}^{-1} \right) \mathbf{w}_{\text{MAP}} = \dfrac{1}{\sigma^{2}} \mathbf{X}^{T} \mathbf{y} + \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu} \\ \implies && \mathbf{w}_{\text{MAP}} = \left(\dfrac{1}{\sigma^{2}} \mathbf{X}^{T}\mathbf{X} + \boldsymbol{\Sigma}^{-1} \right)^{-1} \left(\dfrac{1}{\sigma^{2}} \mathbf{X}^{T} \mathbf{y} + \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu} \right) \end{align*}

事前分布=ラプラス分布

いくつかの仮定や計算は、事前分布がラプラス分布であることを除けば、上と同じである。最大事後確率推定は次を満足するwMAP\mathbf{w}_{\text{MAP}}を求めるものである。

wMAP=arg max⁡wp(y∣w,X)p(w) \mathbf{w}_{\text{MAP}} = \argmax_{\mathbf{w}} p(\mathbf{y} | \mathbf{w}, \mathbf{X}) p(\mathbf{w})

p(y∣w,X)p(\mathbf{y} | \mathbf{w}, \mathbf{X})は尤度であり、p(w)p(\mathbf{w})は事前確率である。尤度関数は次の通りである。

p(y∣w,X)=1(2πσ2)K/2exp⁡[−12σ2∥y−Xw∥22] p(\mathbf{y} | \mathbf{w}, \mathbf{X}) = \dfrac{1}{(2\pi \sigma^{2})^{K/2}} \exp \left[ -\dfrac{1}{2\sigma^{2}} \| \mathbf{y} - \mathbf{X}\mathbf{w} \|_{2}^{2} \right]

また、w\mathbf{w}の事前分布に関して、それぞれのwiw_{i}が独立に🔒(25/05/16)ラプラス分布Laplace⁡(μ,b)\operatorname{Laplace}(\mu, b)に従うと仮定しよう。

p(w)=∏i=1n12bexp⁡(−∣wi−μ∣b) p(\mathbf{w}) = \prod_{i=1}^{n} \dfrac{1}{2b} \exp \left( -\dfrac{|w_{i} - \mu|}{b} \right)   ⟹  log⁡p(w)=−nlog⁡2b−∑i=1n[∣wi−μ∣b] \implies \log p(\mathbf{w}) = -n\log 2b - \sum_{i=1}^{n} \left[ \dfrac{|w_{i} - \mu|}{b} \right]

事後確率が指数関数で表現されるため、対数尤度を考慮することが計算上便利である。

wMAP=arg max⁡wlog⁡(p(y∣w,X)p(w))=arg max⁡w[−12σ2∥y−Xw∥22−∑i=1n∣wi−μ∣b]=arg min⁡w[12σ2∥y−Xw∥22+1b∥w−μ1∥1] \begin{align*} \mathbf{w}_{\text{MAP}} &= \argmax_{\mathbf{w}} \log (p(\mathbf{y} | \mathbf{w}, \mathbf{X}) p(\mathbf{w})) \\ &= \argmax_{\mathbf{w}} \left[-\dfrac{1}{2\sigma^{2}} \| \mathbf{y} - \mathbf{X}\mathbf{w} \|_{2}^{2} - \sum_{i=1}^{n} \dfrac{|w_{i} - \mu|}{b} \right] \\ &= \argmin_{\mathbf{w}} \left[\dfrac{1}{2\sigma^{2}} \| \mathbf{y} - \mathbf{X}\mathbf{w} \|_{2}^{2} + \dfrac{1}{b} \| \mathbf{w} - \mu\mathbf{1} \|_{1} \right] \\ \end{align*}

ここで、1∈Rn\mathbf{1} \in \mathbb{R}^{n}はすべての成分が1のベクトルである。ここでμ=0\mu = 0とすれば、ラッソ回帰と同じ形態になる。

arg min⁡w[12σ2∥y−Xw∥22+1b∥w∥1] \argmin_{\mathbf{w}} \left[\dfrac{1}{2\sigma^{2}} \| \mathbf{y} - \mathbf{X}\mathbf{w} \|_{2}^{2} + \dfrac{1}{b} \| \mathbf{w} \|_{1} \right]

残念ながら、この場合、最適解の閉じた形closed formがないことが知られている。

関連情報