ebpmf_rankk_demo2

Algorithm

See detail in https://www.overleaf.com/project/5bd084d90a33772e7a7f99a2

Implementation

• I did a demo of the same model in https://zihao12.github.io/ebpmf_demo/ebpmf_rankk_demo.html\

• However, the implementation is too slow. The computational bottle-neck is the step to compute \(E(Z_{ijk})\). In the previous implemementation, I iterate over nonzero elements in X to compute \(E(Z_{ijk})\), but even with very sparse dataset, its advantage is not justified because the for loop in R is too slow.
  

• Instead, here I use efficient algorithm to compute softmax3d for that step (see https://zihao12.github.io/ebpmf_demo/softmax_experiments.html). It turns out to be much faster (though still 20 times slower than nnmf).

Summary

Now the computational bottleneck becomes the ebpm algorithm.

## ================== main function ==================================
ebpmf_rankk_exponential <- function(X, K, m = 2, maxiter.out = 10, maxiter.int = 1, seed = 123){
  X = as(X, "dgTMatrix") ## triplet representation: i,j, x
  set.seed(123)
  start = proc.time()
  qg = initialize_qg(X, K)
  runtime_init = (proc.time() - start)[[3]]
  
  runtime_rank1 = 0
  runtime_ez = 0
  
  #print(sprintf("init takes: %f seconds", runtime_init))
  for(iter in 1:maxiter.out){
    #print(sprintf("iter: %d", iter))
    for(k in 1:K){
      #print(sprintf("k: %d", k))
      ## get row & column sum of <Z_ijk>
      start = proc.time()
      Ez = get_Ez(X, qg, k, K) 
      runtime_ez = runtime_ez + (proc.time() - start)[[3]]
      #print(sprintf("compute Ez takes %f seconds", runtime_ez))
      ## update q, g
      start = proc.time()
      tmp = ebpmf_rank1_exponential_helper(Ez$rsum,Ez$csum,NULL,m, maxiter.int) 
      runtime_rank1 = runtime_rank1 + (proc.time() - start)[[3]]
      qg = update_qg(tmp, qg, k)
    }
  }
  print("summary of  runtime:")
  print(sprintf("init           : %f", runtime_init))
  print(sprintf("Ez     per time: %f", runtime_ez/(iter*K)))
  print(sprintf("rank1  per time: %f", runtime_rank1/(iter*K)))
  return(qg)
}
## ================== helper functions ==================================
## for each pair of l, f, give them 1/k of the row & col sum
initialize_qg <- function(X, K, seed = 123){
  n = nrow(X)
  p = ncol(X)
  set.seed(seed)
  X_rsum = rowSums(X)
  X_csum = colSums(X)
  prob_r = replicate(n, rdirichlet(1,replicate(K, 1/K)))[1,,] ## K by n
  prob_c = replicate(p, rdirichlet(1,replicate(K, 1/K)))[1,,] ## K  by p
  rsums = matrix(replicate(K*n,0), nrow = K)
  csums = matrix(replicate(K*p,0), nrow = K)
  for(i in  1:n){
    if(X_rsum[i] == 0){rsums[,i] = replicate(K, 0)}
    else{rsums[,i] = rmultinom(1, X_rsum[i],prob_r[,i])}
  }
  for(j in  1:p){
    if(X_csum[j] == 0){csums[,j] = replicate(K, 0)}
    else{csums[,j] = rmultinom(1, X_csum[j],prob_c[,j])}
  }
  qg = list(qls_mean = matrix(replicate(n*K, 0), ncol =  K), qls_mean_log = matrix(replicate(n*K, 0), ncol =  K), gls = replicate(K, list(NaN)), 
            qfs_mean = matrix(replicate(p*K, 0), ncol =  K), qfs_mean_log = matrix(replicate(p*K, 0), ncol =  K), gfs = replicate(K, list(NaN))
            )
  for(k in 1:K){
    qg_ = ebpmf_rank1_exponential_helper(rsums[k,], csums[k, ], init = NULL, m = 2, maxiter = 1)
    qg   = update_qg(qg_, qg, k)
  }  
  return(qg)
}

## compute the row & col sum of <Z_ijk> for a given k
get_Ez <- function(X, qg, k, K){
  n = nrow(X)
  p = ncol(X)
  psi = array(dim = c(n, p, K))
  ## get <ln l_ik> + <ln f_jk>
  for(d in 1:K){
    psi[,,d] = outer(qg$qls_mean_log[,d], qg$qfs_mean_log[,d], "+")
  }
  ## do softmax
  #browser()
  psi = softmax3d(psi)
  Ez = as.vector(psi)*as.vector(X)
  dim(Ez) = dim(psi)
  return(list(rsum = rowSums(Ez[,,k]), csum = colSums(Ez[,,k])))
}

softmax3d <- function(x){
  score.exp <- exp(x)
  probs <-as.vector(score.exp)/as.vector(rowSums(score.exp,dims=2))
  dim(probs) <- dim(x)
  return(probs)
}

softmax1d <- function(x){
  return(exp(x - logSumExp(x)))
}

update_qg <- function(tmp, qg, k){
  qg$qls_mean[,k] = tmp$ql$mean
  qg$qls_mean_log[,k] = tmp$ql$mean_log
  qg$qfs_mean[,k] = tmp$qf$mean
  qg$qfs_mean_log[,k] = tmp$qf$mean_log
  qg$gls[[k]] = tmp$gl
  qg$gfs[[k]] = tmp$gf
  return(qg)
}

ebpmf_rank1_exponential_helper <- function(X_rowsum,X_colsum, init = NULL, m = 2, maxiter = 1){
  if(is.null(init)){init = list(mean = runif(length(X_rowsum), 0, 1))}
  ql = init
  for(i in 1:maxiter){
    ## update q(f), g(f)
    sum_El = sum(ql$mean)
    tmp = ebpm::ebpm_exponential_mixture(x = X_colsum, s = replicate(p,sum_El), m = m)
    qf = tmp$posterior
    gf = tmp$fitted_g
    ll_f = tmp$log_likelihood
    ## update q(l), g(l)
    sum_Ef = sum(qf$mean)
    tmp = ebpm_exponential_mixture(x = X_rowsum, s = replicate(n,sum_Ef), m = m)
    ql = tmp$posterior
    gl = tmp$fitted_g
    ll_l = tmp$log_likelihood
    qg = list(ql = ql, gl = gl, qf = qf, gf = gf, ll_f = ll_f, ll_l = ll_l)
  }
  return(qg)
}

Experiment Setup

I simulate all columns of \(L\) from the same exponential mixture, and all columns of \(F\) from another exponential mixture. Then I get \(X_{ij} \sim Pois(\sum_k l_{ik} f_{jk})\) as training, and \(Y_{ij} \sim Pois(\sum_k l_{ik} f_{jk})\) as validation.
In order to get a sparse matrix, I set the rate (scale_b) for the exponential to be large.

sim_mgamma <- function(dist){
  pi = dist$pi
  a = dist$a
  b = dist$b
  idx = which(rmultinom(1,1,pi) == 1)
  return(rgamma(1, shape = a[idx], rate =  b[idx]))
}

## simulate a poisson mean problem
## to do: 
simulate_pm  <-  function(n, p, dl, df, K,scale_b = 10, seed = 123){
  set.seed(seed)
  ## simulate L
  a = replicate(dl,1)
  b = 10*runif(dl)
  pi <- rdirichlet(1,rep(1/dl, dl))
  gl = list(pi = pi, a = a, b= b)
  L = matrix(replicate(n*K, sim_mgamma(gl)), ncol = K)
  ## simulate F
  a = replicate(df,1)
  b = 10*runif(df)
  pi <- rdirichlet(1,rep(1/df, df))
  gf = list(pi = pi, a = a, b= b)
  F = matrix(replicate(p*K, sim_mgamma(gf)), ncol = K)
  ## simulate X
  lam = L %*% t(F)
  X = matrix(rpois(n*p, lam), nrow = n)
  Y = matrix(rpois(n*p, lam), nrow = n)
  ## prepare output
  g = list(gl = gl, gf = gf)
  out = list(X = X, Y = Y, L = L, F = F, g = g)
  return(out)
}

I generate a very sparse, small matrix.

n = 100
p = 200
K = 2
dl = 10
df = 10 
scale_b = 5
sim = simulate_pm(n, p, dl, df, K, scale_b = scale_b)

A summary of the simulation:

[1] "nonzero ratio: 0.075400"

[1] "ll train = -4964.163979"

[1] "ll val   = -5263.255923"

Run ebpmf_rankk_exponential

I put the functions above into a function in package ebpmf.

start = proc.time()
##out_ebpmf = ebpmf_rankk_exponential(sim$X, K, maxiter.out = 100)
out_ebpmf = ebpmf::ebpmf_exponential_mixture(sim$X, K, maxiter.out = 100)

[1] "summary of  runtime:"
[1] "init           : 0.033000"
[1] "Ez     per time: 0.002955"
[1] "rank1  per time: 0.009700"

runtime = proc.time() -  start

[1] "runtime: 2.586000 seconds"

[1] "ll train = -4759.116930"

[1] "ll val   = -5615.310826"

Run nnmf with random initialization

start = proc.time()
out_nmf = nnmf(sim$X, K, loss = "mkl", method = "lee", max.iter = 100, rel.tol = -1)
runtime = proc.time() -  start

[1] "runtime: 0.140000 seconds"

[1] "ll train = -4716.149696"

[1] "ll val   = -Inf"

Run nnmf with initialization from ebpmf result

W0 = out_ebpmf$qls_mean
H0 = t(out_ebpmf$qfs_mean)

start = proc.time()
out_nmf_init = nnmf(sim$X, K,init = list(W0 = W0, H0 = H0), loss = "mkl", method = "lee", max.iter = 100, rel.tol = -1)
runtime = proc.time() -  start

[1] "runtime: 0.167000 seconds"

[1] "ll train = -4448.016719"

[1] "ll val   = -6670.583905"

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] NNLM_0.4.2         gtools_3.8.1       Matrix_1.2-17     
[4] matrixStats_0.54.0 ebpmf_0.1.0        ebpm_0.0.0.9000   

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.2      knitr_1.25      whisker_0.3-2   magrittr_1.5   
 [5] workflowr_1.4.0 lattice_0.20-38 stringr_1.4.0   tools_3.5.1    
 [9] grid_3.5.1      xfun_0.8        git2r_0.25.2    htmltools_0.3.6
[13] yaml_2.2.0      rprojroot_1.3-2 digest_0.6.21   mixsqp_0.1-120 
[17] fs_1.3.1        glue_1.3.1      evaluate_0.14   rmarkdown_1.13 
[21] stringi_1.4.3   compiler_3.5.1  backports_1.1.5