Lasso の Post-Selection Inference

Post-Selection Inferenceの切断分布を求める処理を書いてみました。

J. Lee and Dennis L. Sun and Yuekai Sun and Jonathan E. Taylor, “Exact post-selection inference, with application to the lasso”, The Annals of Statistics, volume 44,
number=3, pages 907-927 (2016}

# X, yを中心化する関数
cent=function(z){
    if(is.matrix(z)){
        p=ncol(z)
        for(j in 1:p)z[,j]=z[,j]-mean(z[,j])
    }
    else z=z-mean(z)
  return(z)
}

# 線形回帰のLassoの係数を求める関数
lasso=function(X,y,lambda){
## X,yは、中心化されていると仮定 
    p=ncol(X)
    out=NULL
    for(j in 1:p){
        SD=sd(X[,j])
        X[,j]=X[,j]/SD
        beta=sum(X[,j]*y)
        if(beta>lambda)out=c(out,(beta-lambda)/SD)
        else if(beta< -lambda)out=c(out,(beta+lambda)/SD)
        else out=c(out,0)
    }
    return(out)
}

# 多面体で表現される区間(モデルと符号で条件付)
intervals=function(X,y,lambda,M,s,k){
    n=nrow(X)
    p=ncol(X)
    m=length(M)
    P=X[,M]%*%solve(t(X[,M])%*%X[,M])%*%t(X[,M])
    XX=X[,M]%*%solve(t(X[,M])%*%X[,M])
    A=rbind(1/lambda*t(X[,-M])%*%(diag(n)-P),
        -1/lambda*t(X[,-M])%*%(diag(n)-P),
        -diag(s)%*%solve(t(X[,M])%*%X[,M])%*%t(X[,M])
        )
    b=c(rep(1,p-m)-as.vector(t(X[,-M])%*%XX%*%s),
        as.vector(rep(1,p-m)+t(X[,-M])%*%XX%*%s),
        -as.vector(lambda*diag(s)%*%solve(t(X[,M])%*%X[,M])%*%s)
        )
    eta=(solve(t(X)%*%X)%*%(t(X)))[k,]
    cc=eta/as.vector(t(eta)%*%eta)
    z=(diag(n)-cc%*%t(eta))%*%as.matrix(y)
    Ac=as.vector(A%*%cc)
    Az=as.vector(A%*%z)
    nu.max=0
    nu.min=0
    nu.zero=0
    for(j in 1:p){
      if(Ac[j]>0)nu.max=max(nu.max,(b[j]-Az[j])/Ac[j])
      else if(Ac[j]<0)nu.min=min(nu.min,(b[j]-Az[j])/Ac[j])
      else nu.zero=min(nu.zero,b[j]-Az[j])
    }
    return(c(as.vector(t(eta)%*%as.matrix(y)),nu.min, nu.max,nu.zero))
}

# 10進数を2進数に変換する関数
binary=function(i,m){
  out=NULL
  for(j in 1:m){
    out=c(2*(i%%2)-1,out)
    i=i%/%2
  }
  return(out)
}

# モデルと符号に関する多面体の区間を、符号を動かして、合併させている。
bind.intervals=function(u,v){
  p=length(u); q=length(v)
  u=c(u,Inf); v=c(v,Inf)
  u.state=1; v.state=1
  w=NULL
  i=1; j=1
  while(i<=p||j<=q){
    if(u[i]<v[j]){
      if(v.state==1)w=c(w,u[i])
      i=i+1
      u.state=-u.state
    }
    else if(u[i]>v[j]){
      if(u.state==1)w=c(w,v[j])
      j=j+1
      v.state=-v.state
    }
    else {
      if(i!=p+1&&u.state==v.state)w=c(w,v[j])
      i=i+1; j=j+1
      u.state=-u.state; v.state=-v.state
    }
  }
  return(w)
}
bind.intervals(u,v)

##  実行例
# データ生成
n=100
p=5
X=matrix(rnorm(n*p),n,p)
y=X[,1]*2-X[,2]*3+rnorm(n)*0.4
# 中心化して、Lassoの係数を求め、アクティブ集合とそれぞれの符号を求める
X=cent(X)
y=cent(y)
lambda=40
beta=lasso(X,y,lambda)
M=NULL
ss=NULL
for(j in 1:p){
  if(beta[j]>0){M=c(M,j); ss=c(ss,1)}
  else if(beta[j]<0){M=c(M,j); ss=c(ss,-1)}
}
# 多面体の区間を求める
m=length(M)
L=2^m
print(intervals(X,y,lambda,M,ss,1)[1])
print(M)
print(ss)
S=NULL
for(i in 1:L){
  s=binary(i,m)
  u=intervals(X,y,lambda,M,s,1)
  print(s)
  print(u[2:3])
  S=bind.intervals(S,u[2:3])
}
S
pnorm(-1)

カーネルの苦手意識をいかに克服するか: 「機械学習のためのカーネル100問 with R」(共立出版)まえがきから抜粋

私は、機械学習の手法の中でもカーネルは、特に苦手でした。福水健次著「カーネル法入門」(朝倉書店)を読もうとして、何度も挫折していました。福水先生を、大阪大学の集中講義にお呼びして、学生と一緒に1週間講義を聞きましたが、本質は理解できませんでした。この書籍は、執筆を着手した当初は、自分自身の苦手意識を払拭することを目的としていました。しかし、本書が完成した現在、どうすれば読者にカーネルの苦手意識から脱却できるかを、お伝えできるようになりました。

ところで、機械学習の研究者ですら、カーネルを理解しないで、使っているだけの人がほとんどです。このページを開いている方は、苦手意識を克服したいという前向きな気持ちをお持ちだと思います。

そのための最短経路として最もおすすめしたいのが、数学を基礎から学ぶということです。カーネルは、その背後にある数学にしたがって動作します。理解するまで考え抜くことが重要です。カーネルの理解に必要な数学は、関数解析と呼ばれるものです(第2章)。線形代数や微分積分ならわかるという方でも、戸惑うことがあるかもしれません。ベクトルといえば有限次元ですが、関数の集合は無限次元で、線形代数として扱えます。完備化という概念が初めてという場合、時間をかけていただければと思います。しかし、この第2章を突破すれば、カーネルのすべてが理解できると思います。

本書は、機械学習の数理100問のシリーズの3巻目(全6巻)になります。書籍ですから、既存のカーネルの書籍が存在するのにどうして出版するのか(いわゆる大義)がないと、出版には至りません。本書の特徴として、以下の点をあげることができます。

  1. 数学的命題として証明し、正しい結論を述べているので、読者が本質までたどり着くことができる。
  2. 機械学習の数理100問シリーズの他書と同様、ソースプログラムと実行例を提示して理解を促している。数式だけであれば、特にカーネルの場合、読者が最後まで理解することは容易ではない。
  3. 関数解析の基本的事柄を理解(第2章)してから、それ以降の章での応用を検討していて、数学の予備知識を前提としていない。
  4. RKHSのカーネルと、ガウス過程のカーネルの両方を検討し、しかも両者の扱いを明確に区別している。本書では、それぞれ第5章、第6章で述べている。

国内外のカーネルの書籍を調査しましたが、上記で2個以上満足しているものはありませんでした。

本書の出版に至るまでに、いろいろな失敗を経験してきました。毎年、機械学習の各分野を数学とプログラミングの演習問題を100問解きながら学ぶ講義(大阪大学大学院)をしています。スパース推定(2018年)、グラフィカルモデル(2019年)では人気をはくし、2020年のカーネルも履修者が100名以上になりました。しかし、毎週2日以上講義の予習をしてのぞんだものの、苦手意識も手伝ってか、その講義はうまくいきませんでした。学生の授業アンケートを見ても明らかでした。ただ、そうした問題点の一つ一つを分析し、改良を加えて、本書が誕生しました。

読者の皆さんが,私と同じ道(試行錯誤で時間やエネルギーを消耗する)を辿らずに,効率的にカーネルを学ぶことができればという思いがあります。本書を読んだからといって直ちに論文が書けるわけではありませんが、確実な基礎が身につきます。難しそうに思えていたカーネルの論文がスムーズに読め、一段高いところからカーネルの全体が見えるようになります。また、機械学習の研究者の方でも楽しめるような内容になっています。皆さんが本書を活用し、それぞれの分野で成功を収めていただければ、幸いと考えています。

スパース推定でSCoTALASSの第2主成分以下を選択するとき

テキスト(7.7)の導出は、$L=u_k^{\top}Xv_k-\mu(u_k^\top u_k-1)$の最大化ではなく、
$$L’:=u_k^{\top}Xv_k-\mu(u_k^\top u_k-1)-\lambda\sum_{i=1}^{k-1}(u_i^\top u_k)^2$$の最大化によって得られます。また、(7.7)の$u_k$は大きさが1である$u_1,\ldots,u_{k-1}$と直交しています。

数学的帰納法で、$h=1,\ldots,k-1$の場合に成立していて、$k$について正しいことを示します。$L’$を$u_k$で微分して0とおくと、$$Xv_k-2\mu u_k-2\lambda\sum_{i=1}^{k-1}u_i=0$$となり、左から$P_{k-1}^\perp$をかけると、$$P_{k-1}^\perp Xv_k-2\mu(I- \sum_{j=1}^{k-1}u_ju_j^\top )u_k-2\lambda\sum_{i=1}^{k-1}(I- \sum_{j=1}^{k-1}u_ju_j^\top )u_i=0$$となります。また$L’$を$\lambda$で微分すると$u_ku_j^\top =0$, $1\leq j\leq k-1$となり、また、数学的帰納法の仮定から、$i\leq j\leq k-1$について $u_iu_j^\top =\delta_{i,j}$が成立します。したがって、前式は、$$P_{k-1}^\perp Xv_k-2\mu u_k=0$$となり、$L’$を$\mu$で微分して0とおいて得られる$\mu^\top \mu=1$より、(7.7)が得られます。

行列分解の更新式

(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定数)というように、最大固有値で置き換えなくてもよいということです(この方が厳密な上界になります)。

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