ebpmf_rankk_demo

Algorithm

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

(Note: in rank1 case, we only need row & column sum of X as input)

## ================== 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)
  qg = initialize_qg(X, K)
  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>
      Ez = get_Ez(X, qg, k) 
      ## update q, g
      tmp = ebpmf_rank1_exponential_helper(Ez$rsum,Ez$csum,NULL,m, maxiter.int) 
      qg = update_qg(tmp, qg, 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){
  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
## since <Z_ijk> != 0 only if X_ij != 0, we only need to loop over nonzero elements of X
get_Ez <- function(X, qg, k){
  #browser()
  rsum = replicate(nrow(X), 0)
  csum = replicate(ncol(X), 0)
  for(l in 1:length(X@i)){
    i = X@i[l] + 1
    j = X@j[l] + 1 ## well the index is  zero-based
    current = X[i,j] * softmax1d(qg$qls_mean_log[i,] + qg$qfs_mean_log[j,])[k] ## <Z_ijk> = X_ij * psi_ijk
    rsum[i] = rsum[i] + current
    csum[j] = csum[j] + current
  }
  return(list(rsum = rsum, csum = csum))
}

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

start = proc.time()
out_ebpmf = ebpmf_rankk_exponential(sim$X, K, maxiter.out = 100)
runtime = proc.time() -  start

It is very slow, and when the data gets much denser, it will be even much slower…

[1] "runtime: 49.857000 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.136000 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.165000 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 ebpm_0.0.0.9000   

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