Objective occasionally decreases -- an example

[gaow]

Related to #26 I'm adding an example here, derived from this vignette.

Simulate some data

library(ashr)
library(mashr)
set.seed(1)
simdata = simple_sims(10000,5,1) # simulates data on 40k tests

# identify a subset of strong tests
m.1by1 = mash_1by1(mash_set_data(simdata$Bhat,simdata$Shat))
strong.subset = get_significant_results(m.1by1,0.05)

# identify a random subset of 5000 tests
random.subset = sample(1:nrow(simdata$Bhat),5000)
data.temp = mash_set_data(simdata$Bhat[random.subset,],simdata$Shat[random.subset,])
Vhat = estimate_null_correlation_simple(data.temp)
rm(data.temp)
data.random = mash_set_data(simdata$Bhat[random.subset,],simdata$Shat[random.subset,],V=Vhat)
data.strong = mash_set_data(simdata$Bhat[strong.subset,],simdata$Shat[strong.subset,], V=Vhat)

Configure flash with customized init function for factors (non-neg constraint)

library(flashr)
library(softImpute)
fd <- flash_set_data(as.matrix(data.strong$Bhat))
nonneg <- function(Y, K = 1) {
    # this is the flashr:::udv_si function
    suppressWarnings(ret <- softImpute(Y, rank.max = K, type = "als", lambda = 0))
    pos_sum = sum(ret$v[ret$v > 0])
    neg_sum = -sum(ret$v[ret$v < 0])
    if (neg_sum > pos_sum) {
      return(list(u = -ret$u, d = ret$d, v = -ret$v))
    } else
    return(ret)
}
ebnm_fn = "ebnm_ash"
ebnm_param = list(l = list(mixcompdist = "normal",
                               optmethod = "mixSQP"),
                           f = list(mixcompdist = "+uniform",
                               optmethod = "mixSQP"))

Run flash

U.f = flash(fd, ebnm_fn=ebnm_fn, ebnm_param=ebnm_param, init_fn = nonneg)

I get the objective decrease warning:

Fitting factor/loading 1 (stop when difference in obj. is < 1.00e-02):
  Iteration      Objective   Obj Diff
          1      -15458.33        Inf
          2      -15330.14   1.28e+02
          3      -15302.79   2.73e+01
          4      -15294.56   8.23e+00
          5      -15292.58   1.97e+00
          6      -15292.24   3.44e-01
          7      -15292.22   2.11e-02
          8      -15292.24  -1.86e-02
Performing nullcheck...
  Deleting factor 1 decreases objective by 6.51e+02. Factor retained.
  Nullcheck complete. Objective: -15292.24
Fitting factor/loading 2 (stop when difference in obj. is < 1.00e-02):
  Iteration      Objective   Obj Diff
          1      -15923.07        Inf
          2      -15588.78   3.34e+02
          3      -15477.17   1.12e+02
          4      -15426.12   5.10e+01
          5      -15392.37   3.37e+01
          6      -15360.21   3.22e+01
          7      -15337.61   2.26e+01
          8      -15327.95   9.66e+00
          9      -15316.61   1.13e+01
         10      -15307.40   9.21e+00
         11      -15298.36   9.04e+00
         12      -15295.33   3.03e+00
         13      -15295.07   2.60e-01
         14      -15294.37   6.99e-01
         15      -15292.93   1.44e+00
         16      -15293.24  -3.08e-01
Performing nullcheck...
  Deleting factor 2 increases objective by 1.00e+00. Factor zeroed out.
  Nullcheck complete. Objective: -15292.23
Warning messages:
1: In verbose_obj_decrease_warning() :
  An iteration decreased the objective. This happens occasionally, perhaps due to numeric reasons. You could ignore this warning, but you might like to check out https://github.com/stephenslab/flashr/issues/26 for more details.
2: In verbose_obj_decrease_warning() :
  An iteration decreased the objective. This happens occasionally, perhaps due to numeric reasons. You could ignore this warning, but you might like to check out https://github.com/stephenslab/flashr/issues/26 for more details.

The difference can be as "large" as -3.08e-01 .

BTW the flash interface does not provide an option to use other stopping rule eg factor and/or loading. Perhaps that's intentional?

[willwerscheid]

This happens because the ashr grid can change at each iteration so there's really no monotonicity guarantee. The warning should probably be updated with this information. The lack of stopping rule options isn't really "intentional" -- we just never got around to it. Other stopping rules are available in flashier.

[gaow]

Thanks @willwerscheid the question stems from this addition we are making to mashr:

stephenslab/mashr#88

Since the group seem to be focusing on getting flashier on CRAN I'm thinking of replacing that implementation with flashier based. If you could have a suggestion on the interface that'd be very much appreciated. I'll approach you on DM if that is easier.