knitr::opts_chunk$set(
    message = FALSE,
    warning = FALSE,
    include = TRUE
)
library(Rcpp)
suppressMessages(library(tidyverse))
suppressMessages(library(R.matlab))
suppressMessages(library(glmnet))
library(rbenchmark)

AEFFP Testing

The purpose of this document is to verify my ADMM implementation matches that of Stephen Boyd https://stanford.edu/~boyd/papers/admm/logreg-l1/logreg.html

Validating agaist matlab implementation

Generate Data

n <- 200
p <- 5
set.seed(815)
beta <- c(rnorm(floor(p/2)), rep(0,ceiling(p/2)))
b0 <- rnorm(1) # random intercept
X <- matrix(rnorm(n*p), nrow = n, ncol = p)
pr <- 1/(1+exp(-X %*% beta + b0))
y = rbinom(n,1,pr)
y <- ifelse(y==1,1,-1)

qplot(as.numeric(y==1), geom = "histogram")

Now save the data in format matlab can read:

writeMat("input.mat", X=X, y=y)

Run with my implementation

Note: Boyd's implementation uses $m \times \mu$ in place of $\lambda$ for some reason, where m is the number of observations.

#sourceCpp("src/lasso_logisticC.cpp")
library(aeffp)
admmlasso_logC(X,y, 3, 0.2)

Runtime Comparison: aeffp versus glmnet

Simulate a dataset

n <- 200
p <- 50
set.seed(815)
beta <- c(rnorm(10), rep(0,40))
b0 <- rnorm(1) # random intercept
X <- matrix(rnorm(n*p), nrow = n, ncol = p)
pr <- 1/(1+exp(-X %*% beta + b0))
y = rbinom(n,1,pr)
y <- ifelse(y==1,1,-1)

qplot(as.numeric(y==1), geom = "histogram")

How long does it take to run glmnet?

benchmark("aeffp" = {
  mod1 <- admmlasso_logC(X,y,1)
  },
  "glmnet" = {
    glmnet(X, y, family = "binomial", alpha = 1)
  },
  replications = 10
)


theandyb/aeffp documentation built on May 8, 2019, 9:09 a.m.