Block principal pivoting algorithm for the L1-regularized linear regression (R)

BlockPrincipalPivoting.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