Accelerating co-ordinate ascent updates for linear regression using DAAREM

Analysis settings

These variables specify how the data are generated: n is the number of simulated samples, p is the number of simulated predictors, na is the number of simulated predictors that have a nonzero effect, se is the variance of the residual.

n  <- 200 
p  <- 500
na <- 10
se <- 2

This specifies the prior on the regression coefficients: it is normal with zero mean and standard deviation s0.

s0 <- 1/se

Set up environment

Load some packages and function definitions used in the example below.

library(MASS)
library(daarem)
library(ggplot2)
library(cowplot)
source("../code/misc.R")
source("../code/datasim.R")
source("../code/ridge.R")

Initialize the sequence of pseudorandom numbers.

set.seed(1)

Simulate data

Simulate predictors with “decaying” correlations.

X <- simulate_predictors_decaying_corr(n,p,s = 0.5)
X <- scale(X,center = TRUE,scale = FALSE)

Generate additive effects for the markers so that exactly na of them have a nonzero effect on the trait.

i    <- sample(p,na)
b    <- rep(0,p)
b[i] <- rnorm(na)

Simulate the continuous outcomes, and center them.

y <- drop(X %*% b + se*rnorm(n))
y <- y - mean(y)

Run basic co-ordinate ascent updates

Set the initial estimate of the coefficients.

b0 <- rep(0,p)

Fit the ridge regression model by running 100 iterations of the basic co-ordinate ascent updates.

out <- system.time(fit1 <- ridge(X,y,b0,s0,numiter = 100))
f1  <- ridge.objective(X,y,fit1$b,s0)
cat(sprintf("Computation took %0.2f seconds.\n",out["elapsed"]))
cat(sprintf("Objective value at solution is %0.12f.\n",f1))
# Computation took 3.80 seconds.
# Objective value at solution is -20.573760535831.

Run accelerated co-ordinate ascent updates

Fitting ridge regression with accelerated co-ordinate ascent updates

out <- system.time(fit2 <- daarridge(X,y,b0,s0,numiter = 100))
f2  <- ridge.objective(X,y,fit2$b,s0)
cat(sprintf("Computation took %0.2f seconds.\n",out["elapsed"]))
cat(sprintf("Objective value at solution is %0.12f.\n",f2))
# Computation took 3.66 seconds.
# Objective value at solution is -20.332771749786.

Plot improvement in solution over time

bhat <- drop(solve(t(X) %*% X + diag(p)/s0^2,t(X) %*% y))
f    <- ridge.objective(X,y,bhat,s0)
pdat <-
  rbind(data.frame(iter = 1:100,dist = f - fit1$value,method = "basic"),
        data.frame(iter = 1:100,dist = f - fit2$value,method = "accelerated"))
p <- ggplot(pdat,aes(x = iter,y = dist,col = method)) +
  geom_line(size = 1) +
  scale_y_continuous(trans = "log10",breaks = 10^seq(-8,4)) +
  scale_color_manual(values = c("darkorange","dodgerblue")) +
  labs(x = "iteration",y = "distance from solution")
print(p)

Session information

sessionInfo()
# R version 3.4.3 (2017-11-30)
# Platform: x86_64-apple-darwin15.6.0 (64-bit)
# Running under: macOS High Sierra 10.13.6
# 
# Matrix products: default
# BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
# LAPACK: /Library/Frameworks/R.framework/Versions/3.4/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] cowplot_0.9.4        ggplot2_3.1.0        daarem_0.3          
# [4] MASS_7.3-48          workflowr_1.3.0.9000
# 
# loaded via a namespace (and not attached):
#  [1] Rcpp_1.0.1        compiler_3.4.3    pillar_1.3.1     
#  [4] git2r_0.25.2.9007 highr_0.6         plyr_1.8.4       
#  [7] tools_3.4.3       digest_0.6.17     evaluate_0.11    
# [10] tibble_2.1.1      gtable_0.2.0      debugme_1.1.0    
# [13] pkgconfig_2.0.2   rlang_0.3.1       yaml_2.2.0       
# [16] withr_2.1.2       stringr_1.3.1     dplyr_0.8.0.1    
# [19] knitr_1.20        fs_1.2.6          rprojroot_1.3-2  
# [22] grid_3.4.3        tidyselect_0.2.5  glue_1.3.0       
# [25] R6_2.2.2          rmarkdown_1.10    callr_2.0.3      
# [28] purrr_0.2.5       magrittr_1.5      whisker_0.3-2    
# [31] backports_1.1.2   scales_0.5.0      htmltools_0.3.6  
# [34] assertthat_0.2.0  colorspace_1.4-0  labeling_0.3     
# [37] stringi_1.2.4     lazyeval_0.2.1    munsell_0.4.3    
# [40] crayon_1.3.4