Objective of final fit lower than initialization

[KaiqianZhang]

While working on non-equispaced wavelet-based susie problem, we came across the following example. The objective of the final fit is lower than that of initialization in this example.

The following chunk just creates an interpolated wavelet basis matrix RW_t, which is used as an input X in SuSiE.

library(wavethresh)
#' @param x is an n-vector of data
#' @return R an n by K interpolation matrix
create_interpolation_matrix = function(x){
  n = length(x)
  K = 2^(ceiling(log2(n)))
  R = matrix(0, n, K)
  for (i in 1:n){
    for (j in 1:K){
      if (j == 1 & x[i] <= 1/K){
        R[i,j] = 1
      } else if (j == floor(K*x[i]) & x[i] > 1/K & x[i] <=1){
        R[i,j] = (j+1) - K*x[i]
      } else if (j == ceiling(K*x[i]) & x[i] > 1/K & x[i] <=1){
        R[i,j] = K*x[i] - (j-1)
      } else R[i,j] = 0
    }
  }
  return(R)
}
n = 100
K = 2^(ceiling(log2(n)))
set.seed(1)
x = sort(runif(n, 0,1))
W <- t(GenW(n=K, filter.number=1, family="DaubExPhase"))
R = create_interpolation_matrix(x)
# tcrossprod(R,W) computes R %*% t(W)
RW_t = tcrossprod(R, W)

Then we run SuSiE in our example. resid.RDS.zip

resid = readRDS("resid.RDS")
s.fix = susie(RW_t,resid,estimate_prior_variance = FALSE,L=1) 
s.est = susie(RW_t,resid,estimate_prior_variance = TRUE,L=1,s_init = s.fix)
fit_objective = susie_get_objective(s.est)
init_objective = susie_get_objective(s.fix)
fit_objective
init_objective

The results show that

[1] -72.86187
[1] -65.57433

Here objective of initialization -65.57433 is larger than that of final fit -72.86187. We hypothesize that we have problems in estimating prior variance. Specifically, it is likely to be with numerical optimization of the prior variance in the single effect regression, where we use uniroot().