Skip to content

Instantly share code, notes, and snippets.

@toshi-k
Created August 14, 2015 09:12
Show Gist options
  • Save toshi-k/3289981840a99c081e60 to your computer and use it in GitHub Desktop.
Save toshi-k/3289981840a99c081e60 to your computer and use it in GitHub Desktop.
Block principal pivoting algorithm for the L1-regularized linear regression (R)
# = = = = = include = = = = = #
library(MASS)
# = = = = = function = = = = = #
obj_func <- function(y, A, x, lambda){
obj_value <- 0.5 * sum((y - A%*%x)^2) + lambda * sum(abs(x))
}
bpp <- function(A, y, lambda, Gram, intercept=TRUE){
# parameter
K_max <- 3
if(intercept){
intercept <- mean(y)
y <- y - intercept
}
Aty <- t(A) %*% y
if(missing(Gram)){
AtA <- t(A) %*% A
}
## step1 - - - - -
p <- ncol(A)
H <- 1:p
F_p <- c()
F_m <- c()
beta <- rep(0,p)
d <- - Aty
k <- K_max
t <- p + 1
## step2 (Block Principal Pivoting) - - - - -
while(1){
F <- union(F_p, F_m)
if(length(F)>0){
d[F_p] <- lambda
d[F_m] <- -lambda
beta[F] <- ginv( AtA[F,F] ) %*% ( Aty[F] - d[F] )
}
d[H] <- Aty[H] - AtA[H,F] %*% beta[F]
beta[H] <- 0
J1 <- H[ d[H] > lambda ]
J2 <- H[ d[H] < -lambda ]
J3 <- F_p[ beta[F_p] < 0]
J4 <- F_m[ beta[F_m] > 0]
G <- c(J1,J2,J3,J4)
if(length(G)==0){
break
}else if(length(G) < t){
t <- length(G)
k <- K_max
G_hat <- G
}else{
if(k >= 1){
k <- k-1
G_hat <- G
}else{
k <- 0
G_hat <- G[length(G)]
}
}
J1_hat <- intersect(G_hat, J1)
J2_hat <- intersect(G_hat, J2)
J3_hat <- intersect(G_hat, J3)
J4_hat <- intersect(G_hat, J4)
H <- union(setdiff(H,union(J1_hat,J2_hat)), union(J3_hat,J4_hat))
F_p <- union(setdiff(F_p,J3_hat), J1_hat)
F_m <- union(setdiff(F_m,J4_hat), J2_hat)
}
return(beta)
}
## <<References>>
## [1] Fast Active-set-type Algorithms for L1-regularized Linear Regression
## Jingu Kim, Haesun Park
## http://jmlr.csail.mit.edu/proceedings/papers/v9/kim10a/kim10a.pdf
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment