分からんこと多すぎ

何故か就職できたので技術磨く。一人前の多重線形代数使いを目指しつつ、機械学習、データマイニングをやる

ICML2013:疎行列の固有ベクトル計算(sPCA)読んで実装した

大昔に書いた下書きを発見したので公開


論文

Fast algorithms for sparse principal component analysis based on Rayleigh quotient iteration
Volodymyr Kuleshov著

概要

疎行列の固有ベクトルを、レイリー商を使う事により高速に解く。

著者によると、既存のSparse Principal Component Analysis(sPCA)と比較しても、勝るとも劣らない精度で、1桁か2桁くらい計算量が少ないそうだ。
理屈上は、行列の行数(=列数)をn、非ゼロ要素の数をkとして、イテレーション毎にO(nk+k^3)
ちなみに既存手法には、べき乗法(PM)とか、半正定値計画問題(SDP)とか、Greedyな探索によるものとか、色々あるらしい。

10000×10000の疎行列(半正定値対象。ランクは3。密度は3*100*100/10000*10000=0.03%)を対象として実験してみた。
結果としてRのeigen関数と、この論文の手法で固有値分解した場合の速度を比較すると、200倍ほど速い。(10秒:1800秒)
自分で書いたゴミコードなので、もっと頑張れると思う。

レイリー商(Rayleigh’s quotient)とは

これを読むとよく分かる。
もうちょっと詳しく知りたい人、例題とか欲しい人はこっち。

要するに、二次形式の問題から発生した、行列の固有値に関する評価指標の一つで、微分すると一般化固有値問題が出てくるものである。
本論文では、行列がスパースなため、こいつを利用し、固有ベクトルに関してL2とL0ノルムで制約しつつ求める。

特徴(飛ばして良いところ)

sPCAにおいて、最も一般的なのは、2010年に出たGPM(Generalized Power Method)と呼ばれるものである。
これは他の既存手法と比較して、同等以上の性能を誇る。
GPMは二次形式の問題を勾配法で解いている。
(さらに疎な行列の特性と、べき乗法を使って加速しているっぽい)
したがって原理としてとても簡単で直感的なものとなっているものの、対象となる行列が大きいと収束が遅い。

対して本論文の手法GRQI(Generalized Rayleigh Quotient Iteration)は、QR分解にも使われるレイリー商を、ニュートン法みたいな二次収束する最適化によって解く
したがって早くて当然という論法。

具体的な理屈

f:id:rishida:20130614220225p:plain

これを上から順に解説していく。

Algorithm1 GRQIでは、イテレーション(反復)毎にレイリー商とべき乗法による更新する。
そして結果得られる固有ベクトルを、Euclidean Projectionする。
EPの日訳がよく分からないが、要するにL0とL2の正則化が作る超平面(というか球)に対して、射影する。
詳しい理屈はこちら

  1. レイリー商の計算
    レイリー商R_{\Sigma}は、定義より以下のように求められる。
    R_{\Sigma} = \frac{x^T \Sigma x}{x^T x}

    余計な計算をしたくないので、固有ベクトルxの非ゼロ要素のみを抜き出す。
    ゼロの要素は、積をとってもゼロのままなので無意味。

  2. レイリー商の更新
    x_{new} = \left( \Sigma - \mu x\right)^{-1} x
    \mu = \frac{x^T \Sigma x}{x^T x}

    もし\mu固有値\lambda_iに近ければ、\left( \Sigma - \mu x\right)^{-1}のi番目の固有値は、無限大になる。

    どういうことかと言うと固有ベクトルを利用して、行列は対角化できる。
    P^{-1} A P = \Lambda
    この原理を利用する。

    A - \lambda_{i} I
    = A - \lambda_{i} P P^{-1}
    = P \left( P^{-1} A P - \lambda_{i} P^{-1} P P^{-1} P \right) P^{-1}
    = P \left( P^{-1} A P - \lambda_{i} I \right) P^{-1}

    したがって、 \lambda_{i}の要素だけがほぼ0になる。
    この逆行列を取るので、固有ベクトルx_{i}にかかる固有値(の逆比)が無限大になる。

    よって、固有ベクトルx_{i}との内積が大きい要素が、とんでもない勢いで拡大される。
    つまり、x_{new}固有ベクトルとほぼ平行になる。

  3. べき乗法の計算
    べき乗法とはこういうもので、
    要するに、行列を何回もかけて、高速に最大固有値固有ベクトルを求める方法である。
  4. Euclidean Projection
    要するに、制約条件を満たす範囲への射影である。
    今回の場合、固有ベクトルの要素の内、絶対値が大きいものからk個のみを採用して、残りをゼロにすれば良い。

以上の原理により、非ゼロ要素のみに注目して、固有ベクトルを高速に求めるのが、本手法である。

ソースコード

#generalized Rayleigh quotient iteration
#固有ベクトルを求める際にレイリー商を使う方法の一般化

#固有ベクトルの導出
grqi <- function(Sig,x.new=NULL,k=3,J=2,e=0.01,END=100){
  if(is.null(x.new)){
    x.new <- matrix(1,nrow=dim(Sig)[1],ncol=1)
  }

  for(i in 1:END){
    #keep old vec
    x.old <- x.new

    #compute rayleigh quotient and working set
    mu <- c( (t(x.new) %*% Sig %*% x.new) / (t(x.new) %*% x.new))
    w <- which(x.new != 0)
    if(length(w) == 0){
      break
    }

    #rayleigh quotient iteration updata
    sol.mat <- svd(Sig[w,w] - mu * diag(1,nrow=length(w), ncol=length(w)))
    sol.diag <- 1/sqrt(sol.mat$d)
    sol.diag[sol.diag > 10^10] <- 10^10
    ##これを作らないと1*1行列でバグる
    tmp.diag <- matrix(0,nrow=length(w),ncol=length(w))
    diag(tmp.diag) <- sol.diag
    sol.mat <- sol.mat$u %*% tmp.diag %*% t(sol.mat$v)
    x.new[w,1] <- matrix( sol.mat %*% as.matrix(x.new[w,1]), nrow=length(w), ncol=1)
    x.new <- x.new / sqrt(sum(t(x.new) %*% x.new))
    #power method update
    if(i < J){
      vec <- matrix(Sig %*% x.new,nrow=dim(Sig)[1],ncol=1)
      x.new <- vec / c(t(vec) %*% vec)
    }

    #update
    v <- rep(0,dim(Sig)[1])
    v[order(x.new,decreasing=TRUE)[1:k]] <- 1
    x.new <- x.new * v
    x.new <- matrix(x.new / sqrt(sum(x.new * x.new)),nrow=dim(Sig)[1],ncol=1)

    if(sum(is.nan(x.new)) > 0){
      break
    }
    if( sum((x.new - x.old)^2) < e){
      break
    }
  }
  return(x.old)
}

#sparse principal component analysis(疎行列PCA)
spca.grqi <- function(Sig,k=NULL,J=100,e=0.01,K=NULL,d=1,END=100,DIAG=FALSE){
  #Sig:データの行列(半正定値対象)
  #k:一つの固有ベクトルに現われる非ゼロ要素数
  #J:べき乗法の使用頻度(Jが大きいほど高速で雑)
  #e:許容誤差
  #K:求める固有ベクトルの数
  #d:固有ベクトルの重複を制御する(d=1で重複なし)
  #END:ひとつの固有ベクトルの最大更新回数
  #DIAG:固有値計算(データを複製するのでメモリを消費する)

  if(is.null(K)){
    K <- dim(Sig)[2]
  }
  if(is.null(k)){
    k <- dim(Sig)[2]
  }
  if(DIAG == TRUE){
    data <- Sig
  }

  egn <- matrix(0,nrow=dim(Sig)[1],ncol=K)

  for(j in 1:K){
    i <- which.max(colSums(Sig))
    x.0 <- matrix(Sig[,i],nrow=dim(Sig)[1],ncol=1)
    x.k <- grqi(Sig,x.0,k,J,e,END)
    egn[,j] <- x.k
    Sig <- Sig - d * c( t(x.k) %*% Sig %*% x.k) * (x.k %*% t(x.k))
  }

  ans <- list(egn)
  names(ans) <- c('vectors')
  if(DIAG == TRUE){
    values <- solve(egn) %*% data %*% solve(t(egn))
    ans[[2]] <- values
    names(ans[[2]]) <- c('values')
  }

  return(ans)
}


#使い方
if(FALSE){
  #普通の行列に対する精度の確認
  ##データの生成
  data <- matrix(1:25,nrow=5,ncol=5)
  data[lower.tri(data)] <- t(data)[lower.tri(data)]
  ##処理
  ans.spca <- spca.grqi(data, k=dim(data)[2], J=100, e=0.00001, K=5, d=1, END=1000, DIAG=TRUE)
  ##精度確認
  print(sum((ans.spca[[1]] %*% ans.spca[[2]] %*% t(ans.spca[[1]]) - data)^2))

  #疎行列に対する速度の確認
  ##データ生成
  data <- matrix(0,10000,10000)
  data[1:100,1:100] <- 1
  data[101:200,101:200] <- 10
  data[201:300,201:300] <-5
  ##処理
  system.time(ans.spca <- spca.grqi(data, k=dim(data)[1], J=100, e=0.001, K=3, END=1000, DIAG=FALSE))
  ##eigenとの比較:30分くらいかかるのでおすすめしません
  ##sysytem.time(ans.egn <- eigen(data))
  ##固有ベクトルの確認
  image(log(sqrt(ans.spca[[1]])))
}