ebpmf_rank1_demo

Model: EBPMF-rank1

\[ \begin{align} & X_{ij} \sim Pois(l_i f_j)\\ & l_i \sim g_L(.), g_L \in \mathcal{G}\\ & f_j \sim g_F(.), g_F \in \mathcal{G} \end{align} \]

Algorithm

Described in https://www.overleaf.com/project/5bd084d90a33772e7a7f99a2

I start implememting it for mixture of exponential as \(\mathcal{G}\).

Seems that one iteration is enough. Is it the same case as in MLE for pmf?

# X: a matrix/array of shape n by p 
ebpmf_rank1_exponential <- function(X, init, m = 2, maxiter = 1){
  n = nrow(X)
  p = ncol(X)
  #El = init$ql$mean
  ql = init$ql
  #E_f = get_exp_F(init)
  for(i in 1:maxiter){
    ## update q(f), g(f)
    sum_El = sum(ql$mean)
    tmp = ebpm_exponential_mixture(x = colSums(X), 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 = rowSums(X), 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)
}


compute_elbo <- function(X, qg){
  ql = qg$ql
  gl = qg$gl
  qf = qg$qf
  gf = qg$gf
  ll_f = qg$ll_f
  ll_l = qg$ll_l
  ## compute Eq(logp(X | l, f))
  term1 = sum(- outer(ql$mean, qf$mean, "*") + X*outer(ql$mean_log, qf$mean_log, "+"))
  print(sprintf("term1: %f", term1))
  ## compute Eq(log(gL(l)/qL(l)))
  term2 = ll_l - sum(sum(qf$mean)*ql$mean + rowSums(X)*ql$mean_log) - sum(lgamma(rowSums(X + 1)))
  print(sprintf("term2: %f", term2))
  ## compute Eq(log(gF(f)/qF(f)))
  term3 = ll_f - sum(sum(ql$mean)*qf$mean + colSums(X)*qf$mean_log) - sum(lgamma(colSums(X + 1)))
  print(sprintf("term3: %f", term3))
  return(term1 + term2 + term3)
}

## ===========================================================================
## ==========================experiment setup=================================
## ===========================================================================
## sample from mixture of gamm distribution
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, seed = 123){
  set.seed(seed)
  ## simulate l
  a = replicate(dl,1)
  b = 0.1*runif(dl)
  pi <- rdirichlet(1,rep(1/dl, dl))
  gl = list(pi = pi, a = a, b= b)
  l = matrix(replicate(n, sim_mgamma(gl)), ncol = 1)
  ## simulate f
  a = replicate(df,1)
  b = 0.1*runif(df)
  pi <- rdirichlet(1,rep(1/df, df))
  gf = list(pi = pi, a = a, b= b)
  f = t(matrix(replicate(p, sim_mgamma(gf)), nrow = 1))
  ## 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)
}

## ===========================================================================
## ==========================helper functions ================================
## ===========================================================================
## sample from mixture of gamm distribution

rmse <- function(x,y){
  return(sqrt(mean((x-y)^2)))
}

compute_ll <- function(X, lam){
  return(sum(dpois(X, lam, log = T)))
}

n = 500
p = 1000
# n = 1000
# p = 2000
dl = 3
df = 5  
sim = simulate_pm(n, p, dl, df)

# ## init
# tmp = nnmf(sim$X, k = 1, loss  = "mkl",max.iter = 10)
# ql = list(mean = tmp$W[,1])
# qf = list(mean = tmp$H[1,])
# init = list(ql = ql, qf = qf)

ql = list(mean = runif(n, 0, 1))
qf = list(mean = runif(p, 0, 1))
init = list(ql = ql, qf = qf)

start =  proc.time()
out_ebpmf = ebpmf_rank1_exponential(sim$X, init,maxiter = 10)

[1] "term1: 4312955820.709294"
[1] "term2: -8708921514.333838"
[1] "term3: -13543922943.760603"
[1] "ELBO: -17939888637.385147"
[1] "term1: 4312955820.709294"
[1] "term2: -8708921516.426716"
[1] "term3: -13543922941.667723"
[1] "ELBO: -17939888637.385143"
[1] "term1: 4312955820.709294"
[1] "term2: -8708921518.519596"
[1] "term3: -13543922939.574842"
[1] "ELBO: -17939888637.385143"
[1] "term1: 4312955820.709294"
[1] "term2: -8708921520.612476"
[1] "term3: -13543922937.481964"
[1] "ELBO: -17939888637.385147"
[1] "term1: 4312955820.709294"
[1] "term2: -8708921522.705355"
[1] "term3: -13543922935.389084"
[1] "ELBO: -17939888637.385143"
[1] "term1: 4312955820.709294"
[1] "term2: -8708921524.798235"
[1] "term3: -13543922933.296204"
[1] "ELBO: -17939888637.385143"
[1] "term1: 4312955820.709294"
[1] "term2: -8708921526.891115"
[1] "term3: -13543922931.203323"
[1] "ELBO: -17939888637.385143"
[1] "term1: 4312955820.709294"
[1] "term2: -8708921528.983995"
[1] "term3: -13543922929.110443"
[1] "ELBO: -17939888637.385143"
[1] "term1: 4312955820.709294"
[1] "term2: -8708921531.076876"
[1] "term3: -13543922927.017563"
[1] "ELBO: -17939888637.385143"
[1] "term1: 4312955820.709294"
[1] "term2: -8708921533.169754"
[1] "term3: -13543922924.924685"
[1] "ELBO: -17939888637.385143"

runtime = (proc.time() - start)[[3]]
out_ebpmf[["runtime"]] = runtime

[1] "ebpmf_rank1_exponential fit with 0.650000 seconds"

[1] "ll_train using posterior mean: -2168601.604257"

[1] "ll_val   using posterior mean: -2169032.702840"

• ELBO stops increasing after the first iteration … though there are small updates going on still. Is it because one iteration can get to optimum, or is there a bug?
• Note it can be very slow if we choose m (multiple when selecting grid) to be small (like 1.1)

# plot(sim$l,out_ebpmf$ql$mean, xlab =  "l_sim", ylab = "l_fit", main = "l_sim vs l_fit")
# plot(sim$f,out_ebpmf$qf$mean, xlab =  "f_sim", ylab = "f_fit", main = "f_sim vs f_fit")

Let’s see how nmf does on this dataset

start =  proc.time()
tmp = nnmf(sim$X, k = 1, loss  = "mkl",method = "lee",max.iter = 1, rel.tol = -1, verbose = 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)[[3]]
out_nnmf = list(l = tmp$W[,1], f = tmp$H[1,], runtime = runtime)

[1] "nnmf fit with 0.076000 seconds"

[1] "ll_train using MLE         : -2168601.594216"

[1] "ll_val   using MLE         : -2169032.642280"

# plot(sim$l, out_nnmf$l, xlab =  "l_sim", ylab =  "l_fit", main = "l_sim vs l_fit")
# plot(sim$f, out_nnmf$f, xlab =  "f_sim", ylab =  "f_fit", main = "f_sim vs f_fit")

Note that we only need to run nnmf with lee’s update one iteration to get optimal (up to scaling between L,F), as we have shown before (there is analytic solution, and EM, which Lee’s is, gets to that solution in one step). However, if we use “scd”, one iteration is not enough!

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      ggplot2_3.2.1   mixsqp_0.1-120  gtools_3.8.1   
[5] ebpm_0.0.0.9000

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.2       compiler_3.5.1   pillar_1.4.2     git2r_0.25.2    
 [5] workflowr_1.4.0  tools_3.5.1      digest_0.6.21    evaluate_0.14   
 [9] tibble_2.1.3     gtable_0.3.0     pkgconfig_2.0.3  rlang_0.4.0     
[13] yaml_2.2.0       xfun_0.8         withr_2.1.2      stringr_1.4.0   
[17] dplyr_0.8.1      knitr_1.25       fs_1.3.1         rprojroot_1.3-2 
[21] grid_3.5.1       tidyselect_0.2.5 glue_1.3.1       R6_2.4.0        
[25] rmarkdown_1.13   purrr_0.3.2      magrittr_1.5     whisker_0.3-2   
[29] backports_1.1.5  scales_1.0.0     htmltools_0.3.6  assertthat_0.2.1
[33] colorspace_1.4-1 stringi_1.4.3    lazyeval_0.2.2   munsell_0.5.0   
[37] crayon_1.3.4