#1 Disclaimer
論文が読んだ時のメモを公開するつもりで書きました．マサカリ投げていただけると有り難いです．
#2 はじめに
$D$次元正規分布は平均を０にする標準化を仮定すると，精度行列$\Lambda$，$D$次元ベクトル$\boldsymbol{x}$を用いて以下のようにかけます．この記事では異常検知などの機械学習に必要な精度行列$\Lambda$の推定について書いてみます．

$$p(\boldsymbol{x} ; \Lambda)=\frac{\sqrt{\det \Lambda}}{\left(\sqrt{2 \pi}\right)^D}\exp \left( -\frac{\boldsymbol{x}^T \Lambda\boldsymbol{x}}{2} \right)$$

##2.1 精度行列の最尤推定
平均０の多変量正規分布にしたがうデータから精度行列を最尤推定することを考えます．精度行列を$\Lambda$,分散共分散行列を$S$,対数尤度関数を$L(\Lambda)$とすると

\begin{align}
L(\Lambda) 
&= \sum_{i=1}^{N} \ln p(\boldsymbol{x}_i;\Lambda)\\ 
&= - \frac{ND}{2} \ln 2 \pi + \frac{N}{2} \ln \det \Lambda - \frac{N}{2} \mathrm{trace}(S\Lambda)
\end{align}

したがって正規分布の精度行列の最尤推定は以下になります.
$$
\mathrm{maximize}\quad \ln \det \Lambda -
\mathrm{trace}(S\Lambda)
$$
これの解は偏微分を行うことで
$$\Lambda = S^{-1}$$となります．

##2.2 Graphical Lasso
一般に多変量になり多重共線性があると，分散共分散$S$がランク落ちもしくはランク落ちに近い状態になります．そこで，精度行列の推定にL1正則化項を加え正則にし，さらに疎な精度行列を求めることで解釈性も与える手法があります．最適化問題として書くと以下のようになります．

\begin{align}
\begin{cases}  
\text{maximize} & \ln \det \Lambda -\mathrm{trace}(S\Lambda)-\rho \|\Lambda\|_1\\ 
\text{subject to} & \Lambda \succ O 
\end{cases}\\
\|\Lambda\|_1 = \sum_{ij}\|\Lambda_{ij}\|_1
\end{align}

ここでの正則化は次のように考えられます．分散共分散行列がランク落ちに近い時．固有値に0に近いものが存在します．精度行列を計算する時に,その0に近い固有値の逆数を取るので精度行列に非常に大きな要素が存在することになります．そこで精度行列に大きな要素が存在することに罰則が加わるような正則化項を加えます．Graphical Lassoについての詳細はこの記事がわかりやすいと思います．
##2.3 ROPE
Graphical LassoのRidge版が今回紹介するROPE(Ridge-type Operator for the Precision matrix Estimation)です[1]．

\begin{align}
\text{maximize} \ln \det \Lambda -\mathrm{trace}(S\Lambda)-\rho \|\Lambda\|_2 \tag{1}
\end{align}

式(1)は凸最適化問題であり,最適性条件から劣微分が,要素がすべて0の行列になるときが最適解です．つまり
$$
\Lambda^{-1} - S - 2\rho \Lambda = 0
$$
ここで右辺の0は要素がすべて0の行列を示しています．
$\Lambda$を右から掛けることで以下の式になります.
$$
D(\Lambda) = 2\rho \Lambda^2 + S\Lambda - I = 0
\tag{2}
$$
ここで式(2)の解が対称行列であることを仮定します．式(2)に\式(2)を転置したものを加えると，
$$
-S\Lambda - \Lambda S - \Lambda 4\rho I \Lambda +2I = 0
\tag{3}
$$
これは代数的Ricatti方程式の形です．
#3 代数的Ricatti方程式
[2]を参考にして．代数的Ricatti方程式について軽くまとめます．数式がややこしいので，この節は読み飛ばしてもらっても構いません．以下の式を代数的Ricatti方程式と呼びます，
$$
A^HX+XA-XRX+Q=0\tag{4}
$$
$R,Q$はエルミート行列，$A,R,Q, X\in \boldsymbol{C}^{n\times n}$．ここで$G$を以下のように定義します．

\begin{align}
  G = 
    \begin{bmatrix}
    A & -R\\
    -Q & -A^H
    \end{bmatrix}
\end{align}

補題3.1　

$\lambda$が$G$の固有値の時$\bar{\lambda}$も$G$の固有値となる．
証明は以下のようになります．

\begin{align}
 S \equiv 
\begin{bmatrix}
    0 & -I_n\\
    I_n & 0
    \end{bmatrix}
\end{align}

とすると，

\begin{align}
SG&= 
\begin{bmatrix}
 0 & -I_n\\
I_n & 0
\end{bmatrix}
\begin{bmatrix}
A & -R\\
-Q & -A^H
\end{bmatrix}
=
\begin{bmatrix}
Q&A^H\\
A &-R
\end{bmatrix}\\
-G^HS&=
-
\begin{bmatrix}
A^H & -Q^H\\
-R^H & -A
\end{bmatrix}
\begin{bmatrix}
0 & -I_n\\
I_n & 0
\end{bmatrix}
=
\begin{bmatrix}
Q&A^H\\
A &-R
\end{bmatrix}\\
\end{align}

となるから
$$
SGS^{-1} = -G^H
$$
ここで随伴行列の行列式は、もとの行列の行列式の複素共役であるから$p(\lambda ) = \det (\lambda I_{2n} -G) = 0$から$ \det (-\bar{\lambda} I_{2n} -G^H) = 0$
$$
p(\bar{\lambda}) = \det (\bar{\lambda} I_{2n} -G)
= \det (S) \det (\bar{\lambda} I_{2n} -G) \det(S^{-1})
= \det (\bar{\lambda} I_{2n} +G^H) = 0
$$

3.2 Jordan分解

ここで．$R$と$Q$は半正定値行列であると仮定すると$G$は純虚数の固有値を持たない[3]~~証明難しくて追えませんでした.~~よって，$G$が以下の形でJordan分解することができます．

\begin{align}
    G = U \begin{bmatrix}
    J_1 & 0\\
    0 & J_2 
    \end{bmatrix}
    U^{-1}
\end{align}

ここで$J_1,J_2 \in \boldsymbol{C}^{n\times n}$,$\sigma(J_1) \in \Pi_-,\sigma(J_2) \in \Pi_+$,$\sigma$は固有値の集合を表し,$\Pi_-$は複素平面での左半面,$\Pi_+$は右半面を表す．
$U$は以下のように$n\times n$の行列を用いて表せます．

\begin{align}
U = \begin{bmatrix}
U_{11}& U_{12}\\
U_{21} & U_{22}\\
\end{bmatrix}
\end{align}

ここで式(4)の解は以下のようになります．詳しくは[4]を参考にしてください．
$$
X = U_{21}U_{11}^{-1}
$$
試しに，$A^HX+XA-XRX+Q=0$に代入してみます．

\begin{align}
A^HX+XA-XRX+Q&=
\begin{bmatrix}
X &  -I_n
\end{bmatrix}
\begin{bmatrix}
A & -R\\
-Q & -A^H
\end{bmatrix}
\begin{bmatrix}
I_n\\
X
\end{bmatrix}\\
&=
\begin{bmatrix}
U_{21}U_{11}^{-1} &  -I_n
\end{bmatrix}
G
\begin{bmatrix}
I_n\\
U_{21}U_{11}^{-1}
\end{bmatrix}\\
&=
\begin{bmatrix}
U_{21}U_{11}^{-1} &  -I_n
\end{bmatrix}
G
\begin{bmatrix}
U_{11}\\
U_{21}
\end{bmatrix}
U_{11}^{-1}\\
&=
\begin{bmatrix}
U_{21}U_{11}^{-1} &  -I_n
\end{bmatrix}
\begin{bmatrix}
U_{11}\\
U_{21}
\end{bmatrix}
J_1
U_{11}^{-1}\\
&=
(
U_{21}-U_{21}
)
J_1
U_{11}^{-1}\\
&= 0
\end{align}

4 分散共分散行列の固有値との関連

3節で代数的Ricatti方程式の解について述べましたが，[1]によって式(3)の解が分散共分散行列$S$の固有値と固有ベクトルを用いて表せることが示されています．~~こっちの形の解を導出するべきだが，線形代数力が足らなかった．~~　$l_i$を$S$の固有値として

\begin{align}
    \hat{\lambda_i} &= \frac{2}{l_i + \sqrt{l^2_i + 8\rho}}\\
 V &= \mathrm{diag} (\hat{\lambda_1},......,\hat{\lambda_p}),
    \hat{\lambda_1} \geq\cdots  \geq\hat{\lambda_p} 
\end{align}

とすると
$$
\hat{\Lambda} =
M V M^T
$$と表せます．ここで$M$は$l_i$に対応する固有ベクトルを，$V$と同じ順序になるように並べて作った行列です．

#5 PCAとの関連

5.1 マハラノビス距離

マハラノビス距離とは正規分布に従う多変量の距離であり，ホテリング$T2$理論といった異常検知手法に用いられます．
具体的には,$\Lambda$を精度行列として，以下のように示せます．
$$
a(\boldsymbol{x}^\prime) = (\boldsymbol{x}^\prime - \hat{\boldsymbol{\mu}})^T \Lambda (\boldsymbol{x}^\prime - \hat{\boldsymbol{\mu}})
$$
多重共線性があると分散共分散行列はランク落ちしがちで，精度行列をうまく求められなくなるので，ホテリング$T2$理論と言った異常検知手法は，多重共線性があるときにうまく行きません．

5.2 PCA

PCAは，分散共分散行列の固有値が大きいもの順に選び，小さいものを捨てるといった次元削減の考えで前述の問題を解決しています．以下，標準化を仮定してPCAについて軽く述べます．
データ行列を$X\in \boldsymbol{R}^{n\times p}$とするとその特異値分解は
$U \in \boldsymbol{R}^{n\times n}$,$V \in \boldsymbol{R}^{p\times p}$,$\Sigma \in \boldsymbol{R}^{n\times p}$として
$$
\frac{1}{\sqrt{n}}X = U \Sigma V^T
$$
これによって，分散共分散行列$S$の固有値分解は以下のように与えることができます．
$$
S = \frac{1}{n}X^TX= V \Sigma^2 V^T
$$
PCAを用いたマハラノビス距離のようなもの($T^2$統計量と呼ばれる)は次のように書けます．
$$
T^2 = x^T P \Sigma_a^{-2} P^T x
$$
ここで$\Sigma_a$は$\Sigma$のうち大きいものから$a$個選んだもので$P$はそれらに対応する固有ベクトルからなる行列です．これが前述の固有値が小さいものを捨てるということです．

##5.3 ROPE
ROPEでは，精度行列を求める時に，精度行列の固有値を純粋に分散共分散の固有値$l_i$の逆数にするのではなく$l_i$よりも大きなものの逆数をとる．つまり
$$
\frac{1}{\hat{\lambda_i}} = \frac{l_i + \sqrt{l^2_i + 8\rho}}{2}
$$
とすることで$l_i$が0に近い場合でも$\hat{\lambda_i}$が大きくなることを防ぎます．PCAは共分散行列の固有値が小さいものを捨てる一方で．ROPEは小さい固有値を大きくしていると考えることができます．

#6 まとめ
ROPEはPCAと似てますが，PCAは分散共分散行列の固有値が小さいものをすてる一方．ROPEは固有値を大きくするという異なったアプローチをとっています．似てるなと思い，シミュレーションデータの異常検知において，PCAとROPEを比較してみましたが，PCAの方が強かったですね．なんでだろう．

コード

w,v = np.linalg.eigh(train_covariance)
def ridge(rho):
    lambda_hat = 2/(w+np.sqrt(w**2+8*rho))
    precision = v @ np.diag(lambda_hat) @ v.T
    return precision

#参考文献
[1]M. O. Kuismin, J. T. Kemppainen, M. J.Sillanpää. Precision matrix estimation with rope.Journalof Computational and Graphical Statistics, Vol. 26, No. 3, pp. 682-694, 2017.
[2]Harry Dym.Linear Algebra in Action (Graduate Studies in Mathematics). American MathematicalSociety, 2 edition, 12 2013.
[3]Chris Paige and Charles Van Loan. A schur decomposition for hamiltonian matrices.Linear Algebraand its Applications, Vol. 41, pp. 11 - 32, 1981.
[4]Alan J. Laub. Schur method for solving algebraic riccati equations.Laboratory for Info. and DecisionSystems Report No. LIDS-R-859, 10 1978