行列分解の更新式

(6.11)の劣微分をとって、(6.12)になります。そして、(6.13)を繰り返し適用すると、(6.12)の等号が成立するようになります。$M=\sum_{i=1}^r d_i{u}_i{v}_i$ ($r$は$M$の階数)と書くとき、$[\Omega,Z,M]$の$n$個の特異値で$\lambda$以下のものは、$M$の$r$個の特異ベクトルに対応する場合、しない場合のいずれでも、$S_\lambda(\cdot)$の操作によって0になります。

前者の場合、$S_\lambda([\Omega,Z,M])$の対応する特異値が0になり、$M=S_\lambda([\Omega,Z,M])$であれば、階数$r$が$M$と一致せず、$M$の階数を減らして、(6.13)の更新をさらに続けます。後者の場合は、$M$の特異値とは直交する特異ベクトルで特異値が$\lambda$以下のもの($\sum_{k=r+1}^nd_k\tilde{u}_k\tilde{v}_k$, $d_k\leq \lambda$)が、(6.13)の更新を行っても、$[\Omega,Z,M]$に含まれたままになります。

したがって、(6.12)の最初の項が$M- [\Omega,Z,M]$とかけ、(6.13)が平衡状態に達した場合、(6.12)全体が0になります。

Graphical Lassoの解がなぜ正定値になるのか

スパース推定の 5.2節にあるgraphical lassoの解が非負定値になるのかというご質問をいただきました。研究者でも理解している人がほとんどいない問題です。難しいと感じた方は、遠慮なくご質問ください。

$W$は最初サンプル共分散行列$S$であり($W=S+\lambda I$で始める方法もある)、正定値行列の逆行列は正定値ですから、$W$が最後まで非負定値であることを示せればよいことになります。すなわち、毎回の$W$の更新が非負定値を非負定値に変換するものであれば良いことになります。まず、$$ \det \left[ \begin{array}{ll}W_{11}&w_{12}\\w_{21}&w_{22}\end{array}\right]=\det \left[ \begin{array}{ll}W_{11}&w_{12}\\0&w_{22}-w_{21}W_{11}^{-1}w_{12}\end{array}\right]$$となることに注意します。更新前に$W$が非負定値なら$,W_{11}$も非負定値であり、$W_{11}$は更新しないので、右辺の右下の成分が減少しなければ、更新後も非負定値です。以下では、$w_{21}W^{-1}w_{12}$が増加しないことを示します。(5.11)は、$$\frac{1}{2}\beta^TW_{11}^{-1}\beta-\beta^{T}s_{12}+\lambda\|\beta\|_1$$の最小化になります。$\beta=W^{^1}_{11}(s_{12}+\lambda\psi_{12})$を代入すると、$\beta_j\not=0, =0$によらず$|\beta_j|=-\lambda\psi_{12,j}\beta_j$となり($\psi_{12,j}$は$\psi_{12}$の第$j$成分)、$\|\beta\|_1=-\lambda\psi_{12}^T\beta$となる。したがって、
$$-\frac{1}{2}(s_{12}+\lambda\psi_{12})^TW_{11}^{-1}(s_{12}+\lambda\psi_{12})$$となるので、双対問題である$\|\gamma\|_\infty\leq \lambda$のもとでの
$$-\frac{1}{2}(s_{12}+\gamma)^TW_{11}^{-1}(s_{12}+\gamma)$$の最大化になり、その最適解$\gamma$を用いて、$\psi_{12}=\gamma/\lambda$および $w_{12}=s_{12}+\lambda\psi_{12}$((5.9)式)が得られます。そして、次回($p$回後)回ってきたときに$W_{11}$の値は変わっていますが、$W$と$W_{11}$の正定値性は維持されます。更新によって、$w_{12}$の値はより良い値となるので、$w_{21}W^{-1}_{11}w_{12}$の値はより小さくなり、正定値性が維持されます。

参考文献: R. Mazumder and T. Hastie, “The graphical lasso: New insights and alternatives”, Electronic Journal of Statistics. Vol 6, pp 2125-2149 (2012)

「カーネルの機械学習への応用」が登場

カーネルの機械学習への応用の原稿を共立出版に提出しました。誤植を直したり、まえがきや索引をつけるのに2ヶ月位かかるのが普通です。数学だけでは難しくてめげやすい分野です。今回もプログラムをたくさん入れました。Python版もまもなく出てくると思います。

機械学習の数理100 問シリーズ3 「カーネルの機械学習への応用100 問with R」鈴木 讓

第1章 正定値カーネル
1.1 行列の正定値性
1.2 カーネル
1.3 正定値カーネル
1.4 確率
1.5 Bochnerの定理
1.6 文字列、木、グラフのカーネル

第2章 Hilbert空間
2.1 距離空間と完備性
2.2 線形空間と内積空間
2.3 Hilbert 空間
2.4 射影定理
2.5 線形作用素
2.6 コンパクト作用素

第3章 再生核Hilbert空間
3.1 RKHS
3.2 Sobolev 空間
3.3 Mercer の定理

第4章 カーネル計算の実際
4.1 カーネルRidge 回帰
4.2 カーネル主成分分析
4.3 カーネルSVM
4.4 スプライン
4.5 Random Fourier Features
4.6 Nystrom 近似
4.7 不完全Cholesky 分解

第5章 MMDとHSIC
5.1 RKHSにおける確率変数
5.2 MMDと2標本問題
5.3 HSICと独立性検定
5.4 特性カーネルと普遍カーネル
5.5 経験過程入門

第6章 Gauss 過程と関数データ解析
6.1 回帰
6.2 分類
6.3 補助変数法
6.4 Karhunen-Loeve 展開
6.5 関数データ解析

スパース推定の命題17の証明

条件付き独立性などの関係が、以下の条件を満足するときに、無向グラフの分離性によってすべて表現できることが知られている(J. Pearl, “Probabilistic Reasoning in Intelligent Systems”, 1988)。$$X_A\perp X_B\mid X_C \Longleftrightarrow X_B\perp X_A|X_C\tag{1}$$$$X_A\perp (X_B\cup X_D)\mid X_C \Longrightarrow X_A\perp X_B|X_C, X_A\perp X_D\mid X_C\tag{2}$$ $$X_A\perp (X_C\cup X_D)\mid X_B, X_A\perp_E X_D\mid (X_B\cup X_C) \Longrightarrow X_A\perp (X_B\cup X_D)\mid X_C\tag{3}$$ $$X_A\perp X_B\mid X_C \Longrightarrow X_A\perp X_B\mid (X_C\cup X_D)\tag{4}$$ $$X_A\perp X_B\mid X_C \Longrightarrow X_A\perp X_k \mid X_C\ {\rm or}\ X_k \perp X_B|X_C\tag{5}\ ,\ k\not \in C$$ ただし、各方程式で$A,B,C,D$は、排他的な$V$の部分集合であるとする。このうち、条件付き独立性に関しては(4)(5)以外は、満足されることがわかる。たえとえば(3)は、仮定$$f(x_A,x_B,x_C,x_D)f(x_C,x_D)=f(x_A,x_C,x_D)f(x_B,x_C,x_D) \tag{6}$$ $$f(x_A,x_B,x_C,x_D)f(x_B,x_C)=f(x_A,x_B,x_C)f(x_B,y_C,x_D)\ ,$$から、$f(x_B,x_C)f(x_A,x_C,x_D)=f(x_A,x_B,x_C)f(x_C,x_D)$が得られ、これを周辺化して、$$f(x_C)f(x_A,x_C,x_D)=f(x_A,x_C)f(x_C,x_D)\tag{7}$$が得られる。そして、(6)(7)は$$f(x_A,x_B,x_c,x_D)f(x_C)=f(x_A,x_c)f(x_B,x_C,x_D)$$を意味し、さらに(3)を意味する。また、(4)は命題16より成立する。したがって、(5)が成立することを言えば十分である。まず、$k\in A$または$k\in B$であれば、自明に成立する。そして、正規分布であるので、$k\not\in A\cup B\cup C$のもとで、$X_A\perp X_k \mid X_C$と $X_k \perp X_B|X_C$の両方が成立しないのは、 $X_A\perp X_B|X_C$と矛盾する。したがって、(5)が成立する。

以上より、正規分布の場合、無向グラフの分離性によって、条件付き独立性のすべてが無向グラフで表現できる。

LARS

Lasso では、非ゼロ係数が1個以上ある$\lambda$で最大の値が一意に決まります。2乗誤差を$N$で割らない場合、そのような$\lambda$ は$\langle x^{(j)}, y\rangle$の絶対値の最大値($1\leq j\leq p$)になります。ただし、$x^{(j)}$は$X\in {\mathbb R}^{N\times p}$の$j$列目、$y\in {\mathbb R}^N$であるとします。そのような$\lambda$, $j$を今後$\lambda_0$, $j_0$と書きます。

LARSは、$\lambda\leq \lambda_0$に対して、$\beta(\lambda)$を以下のように定義します: それまでに$\lambda_0,\lambda_1,\ldots,\lambda_{k-1}$, $\beta_0=0,\beta_1,\ldots,\beta_{k-1}$, ${\cal S}=\{j_0,\ldots,j_{k-1}\}$が得られているとして($k\geq 1$)、以下を繰り返します。

1. 最初の$k$成分がゼロ、残り$p-k$成分が0のベクトル$\Delta_{k-1}$を定義
2. $\lambda\leq \lambda_{k-1}$について、 $\beta(\lambda):=\beta_{k-1}+(\lambda_{k-1}-\lambda)\Delta_{k-1}$および$r(\lambda):=y-X\beta(\lambda)$を定義。
3. $\langle x^{(j_k)},r(\lambda_k)\rangle$の絶対値が$\lambda_k$に一致する 最大の$\lambda_k$およびその$j_k$を見出し、$j_k$を$\cal S$に加える。
4. $\beta(\lambda)$, $r(\lambda)$の定義域を$\lambda_k\leq \lambda$に限定し、$\beta_k:=\beta(\lambda_k)$とおく。

ここで、$$\Delta_{k-1}:=\left[
\begin{array}{c}
(X_{\cal S}^TX_{\cal S})^{-1}X_{\cal S}^Tr(\lambda_{k-1})/\lambda_{k-1}\\
0
\end{array}
\right]\ ,
$$
と選べば、各$k=0,1,\ldots,p-1$について、$$\lambda\leq \lambda_k\Longrightarrow \langle x^{(j_k)}, r(\lambda)\rangle=\lambda$$
が成立します。ただし、$X_{\cal S}\in {\mathbb R}^{N\times k}$は、$x^{(j)}$, $j\in {\cal S}$を列に持つ$X$の部分行列とします。

ロジスティック回帰におけるグループLasso

$\beta_j,\gamma_j\in R^{K}$ ($j=1,\ldots,p$)に関して、
$$(\beta_j-\gamma_j)\{\sum_{i=1}^Nx_{i,j}x_{i,j’}W_i\}(\beta_{j’}-\gamma_{j’})^T\leq t\sum_{i=1}^Nx_{i,j}x_{i,j’}(\beta_j-\gamma_j)(\beta_{j’}-\gamma_{j’})^T$$
とできます。$W_i$は非負定値ではありませんので、$tI -W_i$, $i=1,\ldots,N$が非負定値となる$t$の中で最小のものを求めます。したがって、(3.16)の場合と同様、第3項を凸な上界に置き換えることができます。必ずしも、
$$\sum_{i=1}^Nx_{i,j}x_{i,j’}W_i\leq L I$$
(L: Lipschitz定数)というように、最大固有値で置き換えなくてもよいということです(この方が厳密な上界になります)。

これに近接勾配法を適用します。

スパースグループLasso

(3.10)(3.11)の更新式を以前に説明しました。スパースグループLassoは、$f=g+k+h$, $g,k,h$は凸、$g$は微分可能という一般的な場合を扱います。
\[\gamma\leftarrow \beta+\nu X^T(y-X\beta) \tag{3.10}\]
\[\beta \leftarrow (1-\frac{\nu \lambda}{\|\gamma\|_2})_+\gamma\tag{3.11}\]
の間に更新式をいれて、
$$\gamma\leftarrow \beta+\nu X^T(y-X\beta) $$
$$\delta \leftarrow {S}_{\lambda\alpha}(\gamma)$$
$$\beta \leftarrow (1-\frac{\nu \lambda(1-\alpha)}{\|\delta\|_2})_+\delta$$
としたのがスパースグループLassoです。特に、最後の2式をまとめたのが、
\[\beta\leftarrow (1-\frac{\nu\lambda(1-\alpha)}{\|S_{\lambda \alpha}(\gamma)\|})_+S_{\lambda\alpha}(\gamma) \tag{3.23}\]
になります。$k(\gamma)=\alpha\lambda \|\gamma\|_1$, $h(\delta)=(1-\alpha)\lambda\|\delta\|_2$が微分可能であれば
$$\gamma \leftarrow \beta-\nu \nabla g(\beta)$$
$$\delta \leftarrow \gamma-\lambda\alpha$$
$$ \beta \leftarrow \delta -\nu\lambda(1-\alpha)\frac{\delta}{\|\delta\|_2}$$
とできます。

グループLassoの更新式(3.10)(3.11)

\[\gamma\leftarrow \beta+\nu X^T(y-X\beta) \tag{3.10}\]
\[\beta \leftarrow (1-\frac{\nu \lambda}{\|\gamma\|_2})_+\gamma\tag{3.11}\]
を更新することによって、1グループ(複数変数)のグループLassoの解が得られます。これは、近接勾配法で、ある$\nu>0$に対して
$$\beta \leftarrow \beta-\nu (\nabla g+\partial h)$$
なる更新を繰り返すことと同じです。それによって、$f=g+h$の最小化がはかれます。ただし、$g,h$は凸, $g$は微分可能, $\partial h$で$h$の劣微分をあらわします。(3.10)(3.11)の両方に$\nu$があるのは、近接勾配法の漸化式で$\nabla g$, $\partial h$同じが$\nu$倍されていることに起因します。$h(\gamma)=\lambda \|\gamma\|_2$が微分できる場合、$\partial h(\gamma)=\lambda \frac{\gamma}{\|\gamma\|_2}$となり、(3.11)は
$$\beta \leftarrow \gamma-\nu \lambda\frac{\gamma}{\|\gamma\|_2}$$
とかけます。また$g(\beta)=\frac{1}{2}\|y-X\beta\|^2$を微分すると$\nabla g(\beta)=-X^T(y-X\beta)$となり、(3.10)(3.11)は、
$$\gamma \leftarrow \beta-\nu \nabla g$$
$$\beta \leftarrow \gamma-\nu \partial h$$
を交互に行っていることになります。

講義の学生がglmnet for Pythonを動かす環境を作成してくれました。

スパース推定のglmnetは、PythonだとCentOSRでないとが動かないので、悩んでいました。講義は、Rの学生とPythonの学生がいるので。dockerで動作する環境を、受講する学生の一人が作成してくれました。windowsでも動きます。

https://github.com/oriki101/sparse_estimation_docker