In stephenslab/mvebnm: Ultimate Deconvolution for Multivariate Normal Means

set.seed(1)
n <- 500
V <- rbind(c(0.8,0.2),
           c(0.2,1.5))
U <- list(none   = rbind(c(0,0),
                         c(0,0)),
          shared = rbind(c(1.0,0.9),
                         c(0.9,1.0)),
          only1  = rbind(c(1,0),
                         c(0,0)),
          only2  = rbind(c(0,0),
                         c(0,1)))
w <- c(0.8,0.1,0.075,0.025)
X <- simulate_ud_data(n,w,U,V)

Compare unconstrained updates for common V vs. general V

Here, I specify a common V and extend it to a list of Vs. fit1 and fit2 both use unconstrained.update = "ed". In this case, they are doing the same update with the same Vs. So results are expected to be exactly the same.

set.seed(1)
fit1 <- ud_init(X,V = rep(list(V),n))
fit1 <- ud_fit(fit1,
               control = list(unconstrained.update = "ed",maxiter = 10,
                              resid.update = "none",scaled.update = "none",
                              rank1.update = "none"))

set.seed(1)
fit2 <- ud_init(X,V = V)
fit2 <- ud_fit(fit2,
               control = list(unconstrained.update = "ed",maxiter = 10,
                              resid.update = "none",scaled.update = "none",
                              rank1.update = "none"))

plot(fit1$progress$loglik,type = "l",col = "darkblue",lwd = 2)
lines(fit2$progress$loglik,col = "darkorange",lwd = 2,lty = "dashed")
range(simplify2array(fit1$U) - simplify2array(fit2$U))

Compare rank-1 updates for common V vs. general V

If rank1.update = "ed", we are using David's algorithm. If rank1.update = 'teem', we are using the TEEM rank-1 update.

Those results are not expected to be the same.

set.seed(1)
fit1 <- ud_init(X,V = rep(list(V),n))
fit1 <- ud_fit(fit1,control = list(unconstrained.update = "none",
                                   resid.update = "none",
                                   scaled.update = "none",
                                   rank1.update = "ed",
                                   minval = 0,maxiter = 40))

set.seed(1)
fit2 <- ud_init(X,V = V)
fit2 <- ud_fit(fit2,control = list(unconstrained.update = "none",
                                   resid.update = "none",
                                   scaled.update = "none",
                                   rank1.update = "teem",
                                   minval = 0,maxiter = 40))

set.seed(1)
fit3 <- ud_init(X,V = rep(list(V),n))
attr(fit3$U$rank1_1,"covtype") <- "unconstrained"
attr(fit3$U$rank1_2,"covtype") <- "unconstrained"
attr(fit3$U$rank1_3,"covtype") <- "unconstrained"
attr(fit3$U$rank1_4,"covtype") <- "unconstrained"
attr(fit3$U$unconstrained1,"covtype") <- "scaled"
attr(fit3$U$unconstrained2,"covtype") <- "scaled"
attr(fit3$U$unconstrained3,"covtype") <- "scaled"
attr(fit3$U$unconstrained4,"covtype") <- "scaled"
fit3 <- ud_fit(fit3,control = list(unconstrained.update = "ed",
                                   resid.update = "none",
                                   scaled.update = "none",
                                   rank1.update = "none",
                                   minval = 0,maxiter = 100))

plot(fit1$progress$loglik,type = "l",col = "darkblue",lwd = 2)
lines(fit2$progress$loglik,col = "darkorange",lwd = 2,lty = "dashed")
lines(fit3$progress[1:40,"loglik"],col = "limegreen",lwd = 1,lty = "dashed")
print(fit1$loglik,digits = 8)
print(fit2$loglik,digits = 8)
print(fit3$loglik,digits = 8)