In stephenslab/flashr2: Empirical Bayes Matrix Factorization

knitr::opts_chunk$set(echo = TRUE)

A single example

First we will simulate some rank 0 data:

  library("flashr")
  set.seed(1)
  n=100
  p=1000
  Y = matrix(rnorm(n*p),nrow=n)

Run flash to fit a rank 1 model (Kmax=1). Here nullcheck=FALSE tells flash to keep the final value of the variational optimization without checking it against the rank 0 solution (since we are going to do that by hand here.)

  data = flash_set_data(Y)
  f1= flash(data,Kmax=1,nullcheck=FALSE)
  flash_get_objective(data,f1)

Now try deleting the factor - that gives us a rank 0 solution in this case. We can see here that deleting the factor actually increases the objective function! That is the optimization did not work very well (because deleting the factor is actually a special case of the model).

f0 = flash_zero_out_factor(f1,1)
f0 = flash_update_precision(data,f0)
flash_get_objective(data,f0)

This is actually quite a common occurance (see below). So we have implemented this check as the default in flash. Here is what happens if you run flash with default (nullcheck=TRUE): it automatically deletes the factor.

f1 = flash(data,Kmax=1,verbose=TRUE)

Small simulation study:

Let's see how often this happens:

set.seed(1)
n=100
p=1000
F1minusF0 = rep(0,100)
for(i in 1:100){
  Y = matrix(rnorm(n*p),nrow=n)
  data = flash_set_data(Y)
  f1= flash(data,Kmax=1,nullcheck=FALSE)
  f0 = flash_zero_out_factor(f1,1)
  f0 = flash_update_precision(data,f0)
  F1minusF0[i] = flash_get_objective(data,f1) - flash_get_objective(data,f0)
}
hist(F1minusF0)