Graphical lasso の事前の標準化について

先日 Twitter で『異常検知と変化検知』の著者である井手剛さんが，以下のようなツイートをしていました:

<blockquote class="twitter-tweet" data-lang="en"><p lang="ja" dir="ltr">「一般性を失わずに平均0分散1にあらかじめ標準化したと仮定します」，なんて自分の本にも書いちゃったけど，標準化って一般性を失いまくりなんだよな．たぶんコードが動かないという苦情の3割くらいはそれ系．</p>&mdash; Tsuyoshi Ide (井手 剛) (@Idesan) <a href="https://twitter.com/Idesan/status/926615938493009921?ref_src=twsrc%5Etfw">November 4, 2017</a></blockquote>

『異常検知と変化検知』には，標準化してからgraphical lassoをかける下りがあるので，その部分のことを指していらっしゃるのだと思います．私はこれまでこの手順にしたがってgraphical lassoを使っていたため，標準化してしまうとどうなるのか確認しようと思い，計算してみました．

## Graphical lasso
Graphical lassoは，$D$次元の独立同分布標本

$$
\mathbf{X} \equiv [\boldsymbol{x}_1, \ldots, \boldsymbol{x} _N]^\top
$$

が得られたとき，尤度として多変量Gauss分布

$$
p(\boldsymbol{x}\mid \mathbf\Theta) \equiv \mathcal{N}(\boldsymbol{x}; \boldsymbol\mu, \mathbf\Theta^{-1})
$$

を，精度行列 $\mathbf\Theta$ に関する事前分布として

$$
p(\mathbf\Theta) \propto \exp\left(-\rho \lVert\mathbf\Theta\rVert_1\right), \quad\text{where}~~ \lVert\mathbf\Theta\rVert_1 \equiv \sum_{i, j} \lvert\Theta_{ij}\rvert
$$

を仮定し，精度行列 $\mathbf\Theta$ を最大事後確率推定する手法です．

このとき事後確率を最大にする $\mathbf\Theta$ は，$\mathbf{S}$ を標本共分散行列として，

```math
\begin{align}
\hat{\mathbf\Theta}
& \equiv \arg\max_{\mathbf\Theta}\left\{\ln{p(\boldsymbol\Theta)}\prod_{n=1}^N\mathcal{N}(\boldsymbol{x}_n; \boldsymbol\mu, \mathbf\Theta^{-1})\right\} \nonumber \\
& = \arg\max_{\boldsymbol\Theta} \left\{\frac{N}{2}\ln\det\boldsymbol{\Theta} -\frac{1}{2}\sum_{n=1}^N(\boldsymbol{x}_n-\boldsymbol\mu)^\top \boldsymbol\Theta(\boldsymbol{x}_n-\boldsymbol\mu)-\rho||\boldsymbol\Theta||_1 \right\} \nonumber \\
& = \arg\max_{\mathbf\Theta}\left\{\ln\det\mathbf\Theta - \mathrm{tr}(\mathbf{S\Theta}) - \rho^\prime\lVert\mathbf\Theta\rVert_1\right\}.
\end{align}
```

ただし，$\rho^\prime\equiv 2\rho/N$です．この最適化問題は，座標降下法を用いて解くことができます．(正確には、この最適化を高速化したのがgraphical lassoです。)

## Graphical lasso の事前の標準化
先ほどの最適化の式は平均$\boldsymbol\mu$を含まないので，全観測から平均$\boldsymbol\mu$を引いてからgraphical lassoをかけても結果に影響しません(正確には共分散行列に含まれていますが，平均を引いても共分散行列は変わりません)．したがって，以降では

$$
\mathbf{Z} \equiv \boldsymbol{X} - \boldsymbol{1}_N \boldsymbol\mu^\top
$$

として，$\mathbf{Z}$の各次元を標準偏差で割ってからgraphical lassoをかけるのが，$\mathbf{Z}$にそのままgraphical lassoをかけるのとどう違うかを確かめます．

### 記号
以降，以下の表記を用います．

| 定義 | 説明 |
|---|---|
| $\sigma_i \equiv \sqrt{\frac{1}{N}\sum_{n=1}^N z_{ni}^2}$ | 各次元の標準偏差 |
| $\boldsymbol{t} \equiv [\sigma_1, \ldots, \sigma_D]^\top$ | 標準偏差を並べたベクトル |
| $\mathbf{T} \equiv \mathrm{diag}(\boldsymbol{t})$ | 標準偏差を並べた対角行列 |
| $\boldsymbol{b} \equiv [1/\sigma_1, \ldots, 1/\sigma_D]^\top$ | 精度の平方根を並べたベクトル |
| $\mathbf{B} \equiv \mathrm{diag}(\boldsymbol{b})$ | 精度の平方根を並べた対角行列 |

このとき，$\mathbf{T} = \mathbf{B}^{-1}$となります．

### 事前に標準化を行って最適化
標準化した観測は$\mathbf{Z}^\prime \equiv \mathbf{ZB}$と表せるので，標準化したときの共分散行列は，

$$
\mathbf{S}^\prime 
\equiv \frac{1}{N} \mathbf{Z^{\prime\top} Z^\prime}
= \frac{1}{N} \mathbf{(ZB)^\top ZB}
= \frac{1}{N} \mathbf{B^\top Z^\top ZB}
= \mathbf{BSB}.
$$

したがって，標準化した観測の下での最適化の式は，

```math
\begin{align}
\hat{\mathbf\Theta}
& \equiv \arg\max_{\mathbf\Theta}\left\{\ln\det\mathbf\Theta - \mathrm{tr}(\mathbf{S^\prime\Theta}) - \rho^\prime\lVert\mathbf\Theta\rVert_1\right\} \nonumber\\
& = \arg\max_{\mathbf\Theta}\left\{\ln\det\mathbf\Theta - \mathrm{tr}(\mathbf{BSB\Theta}) - \rho^\prime\lVert\mathbf\Theta\rVert_1\right\} \nonumber\\
& = \arg\max_{\mathbf\Theta}\left\{\ln\det\mathbf\Theta - \mathrm{tr}(\mathbf{SB\Theta B}) - \rho^\prime\lVert\mathbf\Theta\rVert_1\right\}.
\end{align}
```

ここで，$\mathbf\Theta^\prime \equiv \mathbf{B\Theta B}$と置くと，$\mathbf{T} = \mathbf{B}^{-1}$より，$\mathbf\Theta = \mathbf{T\Theta^\prime T}$です．これを代入し，改めて$\mathbf\Theta^\prime$について最適化すると，

```math
\begin{align}
\hat{\mathbf\Theta}^\prime
& = \arg\max_{\mathbf\Theta^\prime}\left\{\ln\det(\mathbf{T\Theta^\prime T}) - \mathrm{tr}(\mathbf{S\Theta^\prime}) - \rho^\prime\lVert\mathbf{T\Theta^\prime T}\rVert_1\right\} \nonumber\\
& = \arg\max_{\mathbf\Theta^\prime}\left\{\ln\det\mathbf{\Theta^\prime} - \mathrm{tr}(\mathbf{S\Theta^\prime}) - \rho^\prime\lVert\mathbf{T\Theta^\prime T}\rVert_1\right\} \nonumber\\
& = \arg\max_{\mathbf\Theta^\prime}\left\{\ln\det\mathbf{\Theta^\prime} - \mathrm{tr}(\mathbf{S\Theta^\prime}) - \rho^\prime\lVert(\boldsymbol{tt^\top}) * \mathbf{\Theta^\prime}\rVert_1\right\}.
\end{align}
```

ただし，$*$は行列の要素積を表します．

この式から，観測をあらかじめ標準化してからgraphical lassoをかけるということは，正則化項として$\rho^\prime\lVert\mathbf\Theta\rVert_1$を用いる代わりに，$\rho^\prime\lVert(\boldsymbol{tt^\top}) * \mathbf{\Theta^\prime}\rVert_1$を用いることに対応することが分かります．また，結果として求まる精度行列$\hat{\mathbf\Theta}$は，標準化を行わずに正則化項だけ上記のものを用いた場合に得られる精度行列$\hat{\mathbf\Theta}^\prime$の行と列に，それぞれ標準偏差をかけた行列になります．

この結果についてですが，正則化項に行列との要素積を用いることは，graphical lassoを提案したFriedmanの論文[^1]でも，次の箇所:

> ![Friedmanの論文から引用．](https://qiita-image-store.s3.amazonaws.com/0/117994/42e7bd68-c8b1-694b-672e-ff6f5ef62ff9.png)

のように触れられています．したがって，少なくともgraphical lassoに関しては，あらかじめ標準化を行っても問題ないのではないかと思いました．（ただし，正則化のかかり方が意図したものと違う可能性があることには注意が必要だと思います．）どこかに間違いがあるかもしれませんので，何か気づいた点などありましたら，ご指摘お願いします．


---

[^1]: Friedman et al., Sparse inverse covariance estimation with the graphical lasso, http://statweb.stanford.edu/~tibs/ftp/graph.pdf, 2008.