ebpmf_rankk_demo

workflowr

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

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

Unstaged changes:
    Modified:   analysis/softmax_experiments.Rmd

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

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

library(gtools)

## Note: in rank1 case, what we need is just row and  column sum of X

## TODO: 
## 1. think about how to store qg. They include:
# qls_mean_log = matrix(replicate(n*K, 0), ncol =  K)
# qfs_mean_log = matrix(replicate(p*K, 0), ncol =  K)
# qls_mean = matrix(replicate(n*K, 0), ncol =  K)
# qfs_mean = matrix(replicate(p*K, 0), ncol =  K)
# # ... gls, gfs

## 2. think about the importance of initialization, and how to
## well it is not important at all!! only need the sum of say sum_k <l_ik>, and elbo is "scale invariant""

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){
    for(k in 1: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) ## need to deal with (manipulate) sparse matrix
      qg = update_qg(tmp, qg, k)
    }
  }
  return(qg)
}

## 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)
}

# X: a matrix/array of shape n by p 
# different in that only row and column sum of X is provided
ebpmf_rank1_exponential_helper <- function(X_rowsum,X_colsum, init, m = 2, maxiter = 1){
  #El = init$ql$mean
  if(is.null(init)){init = list(mean = runif(length(X_rowsum), 0, 1))}
  ql = init
  #E_f = get_exp_F(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)
    # elbo = compute_elbo(X, qg)
    # print(sprintf("ELBO: %f", elbo))
  }
  return(qg)
}

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: 
## compute loglik for g (well, is it do-able?)
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.
Note that for some data matrix, we get an error when computing L in later iterations. Error in verify.likelihood.matrix(L) : Input argument "L" should be a numeric matrix with >= 2 columns, >= 1 rows, all its entries should be non-negative, finite and not NA, and some entries should be positive

n = 100
p = 50
K = 2
dl = 10
df = 10 
sim = simulate_pm(n, p, dl, df, K, scale_b = 8)
print(sprintf("nonzero ratio: %f", sum(sim$X != 0)/(n*p)))

[1] "nonzero ratio: 0.082800"

Run ebpmf_rankk_exponential

#out = initialize_qg(X, K)
start = proc.time()
out = 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: 13.664000 seconds"

[1] "ll train = -1228.332333"

[1] "ll val   = -1379.921692"

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

Warning in system.time(out <- .Call("NNLM_nnmf", A, as.integer(k),
init.mask$Wi, : Target tolerance not reached. Try a larger max.iter.

runtime = proc.time() -  start

[1] "runtime: 0.099000 seconds"

[1] "ll train = -1215.263872"

[1] "ll val   = -Inf"

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