ebpmf_rankk

Last updated: 2019-10-05

Checks: 6 1

Knit directory: ebpmf_demo/

This reproducible R Markdown analysis was created with workflowr (version 1.4.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: uncommitted changes

The R Markdown file has unstaged changes. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish to commit the R Markdown file and build the HTML.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20190923)

The command set.seed(20190923) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: f1a17d7

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/

Untracked files:
    Untracked:  analysis/.ipynb_checkpoints/
    Untracked:  analysis/ebpmf_demo.Rmd
    Untracked:  analysis/ebpmf_rankk_demo2.Rmd
    Untracked:  analysis/softmax_experiments.ipynb
    Untracked:  docs/figure/test.Rmd/

Unstaged changes:
    Modified:   analysis/ebpmf_rankk_demo.Rmd
    Modified:   analysis/softmax_experiments.Rmd

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the R Markdown and HTML files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view them.

File	Version	Author	Date	Message
html	076df3e	zihao12	2019-10-03	Build site.
Rmd	4528871	zihao12	2019-10-03	update ebpmf_rankk_demo
html	72909fb	zihao12	2019-10-03	Build site.
Rmd	c0cced7	zihao12	2019-10-03	update ebpmf_rankk_demo
html	92a80e3	zihao12	2019-10-03	Build site.
Rmd	736baa2	zihao12	2019-10-03	update ebpmf_rankk_demo
html	28a19e5	zihao12	2019-10-03	Build site.
Rmd	c635194	zihao12	2019-10-03	update ebpmf_rankk_demo
html	6de69a2	zihao12	2019-10-03	Build site.
Rmd	8c5114a	zihao12	2019-10-03	very naive implementation
html	fa5fa67	zihao12	2019-10-02	Build site.
Rmd	2851c1b	zihao12	2019-10-02	start developing ebpmf rank k

rm(list = ls())
library(ebpm)
library(matrixStats)
library(Matrix)

Warning: package 'Matrix' was built under R version 3.5.2

library(gtools)
library(NNLM)

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)
  start = proc.time()
  qg = initialize_qg(X, K)
  runtime_init = (proc.time() - start)[[3]]
  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) 
      runtime_ez = (proc.time() - start)[[3]]
      print(sprintf("compute Ez takes %f seconds", runtime_ez))
      ## 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){
  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
## 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)

[1] "init takes: 0.080000 seconds"
[1] "compute Ez takes 0.310000 seconds"
[1] "compute Ez takes 0.243000 seconds"
[1] "compute Ez takes 0.250000 seconds"
[1] "compute Ez takes 0.237000 seconds"
[1] "compute Ez takes 0.241000 seconds"
[1] "compute Ez takes 0.240000 seconds"
[1] "compute Ez takes 0.237000 seconds"
[1] "compute Ez takes 0.239000 seconds"
[1] "compute Ez takes 0.239000 seconds"
[1] "compute Ez takes 0.240000 seconds"
[1] "compute Ez takes 0.242000 seconds"
[1] "compute Ez takes 0.243000 seconds"
[1] "compute Ez takes 0.234000 seconds"
[1] "compute Ez takes 0.240000 seconds"
[1] "compute Ez takes 0.245000 seconds"
[1] "compute Ez takes 0.249000 seconds"
[1] "compute Ez takes 0.236000 seconds"
[1] "compute Ez takes 0.294000 seconds"
[1] "compute Ez takes 0.442000 seconds"
[1] "compute Ez takes 0.270000 seconds"
[1] "compute Ez takes 0.232000 seconds"
[1] "compute Ez takes 0.239000 seconds"
[1] "compute Ez takes 0.239000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.228000 seconds"
[1] "compute Ez takes 0.228000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.229000 seconds"
[1] "compute Ez takes 0.219000 seconds"
[1] "compute Ez takes 0.225000 seconds"
[1] "compute Ez takes 0.226000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.228000 seconds"
[1] "compute Ez takes 0.232000 seconds"
[1] "compute Ez takes 0.249000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.225000 seconds"
[1] "compute Ez takes 0.232000 seconds"
[1] "compute Ez takes 0.247000 seconds"
[1] "compute Ez takes 0.228000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.227000 seconds"
[1] "compute Ez takes 0.227000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.227000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.241000 seconds"
[1] "compute Ez takes 0.228000 seconds"
[1] "compute Ez takes 0.245000 seconds"
[1] "compute Ez takes 0.226000 seconds"
[1] "compute Ez takes 0.229000 seconds"
[1] "compute Ez takes 0.221000 seconds"
[1] "compute Ez takes 0.225000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.227000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.226000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.234000 seconds"
[1] "compute Ez takes 0.231000 seconds"
[1] "compute Ez takes 0.228000 seconds"
[1] "compute Ez takes 0.225000 seconds"
[1] "compute Ez takes 0.221000 seconds"
[1] "compute Ez takes 0.221000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.229000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.227000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.230000 seconds"
[1] "compute Ez takes 0.218000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.219000 seconds"
[1] "compute Ez takes 0.221000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.221000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.221000 seconds"
[1] "compute Ez takes 0.218000 seconds"
[1] "compute Ez takes 0.312000 seconds"
[1] "compute Ez takes 0.229000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.221000 seconds"
[1] "compute Ez takes 0.225000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.219000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.219000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.233000 seconds"
[1] "compute Ez takes 0.226000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.225000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.219000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.221000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.232000 seconds"
[1] "compute Ez takes 0.229000 seconds"
[1] "compute Ez takes 0.228000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.219000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.225000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.227000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.228000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.221000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.221000 seconds"
[1] "compute Ez takes 0.225000 seconds"
[1] "compute Ez takes 0.235000 seconds"
[1] "compute Ez takes 0.218000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.229000 seconds"
[1] "compute Ez takes 0.218000 seconds"
[1] "compute Ez takes 0.227000 seconds"
[1] "compute Ez takes 0.221000 seconds"
[1] "compute Ez takes 0.221000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.221000 seconds"
[1] "compute Ez takes 0.217000 seconds"
[1] "compute Ez takes 0.221000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.230000 seconds"
[1] "compute Ez takes 0.227000 seconds"
[1] "compute Ez takes 0.218000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.217000 seconds"
[1] "compute Ez takes 0.221000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.219000 seconds"
[1] "compute Ez takes 0.320000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.218000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.221000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.227000 seconds"
[1] "compute Ez takes 0.232000 seconds"
[1] "compute Ez takes 0.225000 seconds"
[1] "compute Ez takes 0.225000 seconds"
[1] "compute Ez takes 0.228000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.222000 seconds"
[1] "compute Ez takes 0.231000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.233000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.228000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.228000 seconds"
[1] "compute Ez takes 0.229000 seconds"
[1] "compute Ez takes 0.233000 seconds"
[1] "compute Ez takes 0.220000 seconds"
[1] "compute Ez takes 0.226000 seconds"
[1] "compute Ez takes 0.223000 seconds"
[1] "compute Ez takes 0.227000 seconds"
[1] "compute Ez takes 0.224000 seconds"
[1] "compute Ez takes 0.221000 seconds"
[1] "compute Ez takes 0.228000 seconds"

runtime = proc.time() -  start

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

[1] "runtime: 47.749000 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.167000 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.179000 seconds"

[1] "ll train = -4448.016719"

[1] "ll val   = -6670.583905"

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

ebpmf_rankk_demo

zihao12

2019-10-01

Algorithm

Experiment Setup

Run ebpmf_rankk_exponential

Run nnmf with random initialization

Run nnmf with initialization from ebpmf result

Run `ebpmf_rankk_exponential`

Run `nnmf` with random initialization

Run `nnmf` with initialization from ebpmf result