Feature-sign search algorithm (R)

FeatureSignSearch.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))
}

fss <- function(A, y, lambda, Gram, intercept=TRUE){

	obj_vec <- NULL

	if(intercept){
		intercept <- mean(y)
		y <- y - intercept
	}

	Aty <- t(A) %*% y
	if(missing(Gram)){
		AtA <- t(A) %*% A
	}

	## step1 - - - - - 

	M <- ncol(A)
	N <- nrow(A)

	x <- numeric(M)
	theta <- integer(M)
	active <- NULL
	act_indx0 <- c(1:M)
	obj_value <- 0.5*sum(y^2)
	
	## step2 - - - - - 

	grad <- AtA %*% x - Aty
	repeat{

	max_abs <- max(abs(grad[act_indx0]))
	max_indx <- which.max(abs(grad[act_indx0]))

	if(max_abs > lambda){
		active <- c(active,act_indx0[max_indx])
		theta[act_indx0[max_indx]] <- -sign(grad[act_indx0[max_indx]])
	}else{
		return(x)
	}

	repeat{

		## step3 (Feature-sign step) - - - - -

		x_act <- x[active]
		
		x_new <- ginv(AtA[active,active]) %*% (Aty[active] - lambda * theta[active])
		x[active] <- x_new

		# FeatureSignProcess
		sign_check <- active[ (sign(x_act) != sign(x_new))]
		if(length(sign_check) > 0){

			ma <- x_new - x_act
			step_all <- as.vector( x_act / (-ma) )

			check <- step_all > 0 & step_all < 1
			step_vec <- step_all[check]

			if(length(step_vec) > 0){
				obj_values <- numeric(length(step_vec))
				over_step_value <- obj_func(y,A,x,lambda)
				for(j in 1:length(step_vec)){
					x_test <- x_act + step_vec[j] * ma
					obj_values[j] <- obj_func(y,A[,active],x_test,lambda)
				}
				
				min_value <- min(obj_values)
				min_indx <- which.min(obj_values)

				if(over_step_value > min_value){
					x_new <- x_act + step_vec[min_indx] * ma
					x[active] <- x_new
					x[sign_check[min_indx]] <- 0
					
					active <- setdiff(active,sign_check[min_indx])
				}
			}
		}
		theta <- sign(x)

		obj_check <- obj_func(y, A, x, lambda)
		if(obj_check >= obj_value){
			return(x)
		}else{
			obj_value <- obj_check
		}

		# step4 (Check the optimality conditions) - - - - -

		grad <- AtA %*% x - Aty
		act_indx0 <- (1:M)[-active]

		condition_a <- !sum(abs(grad[active] + lambda * sign(x[active])) > 1e-10 )
		
		if(condition_a){
			condition_b <- !sum(abs(grad[-active]) > lambda)
			if(!condition_b){
				break
			}else{
				return(x)
			}
		}
	}

	}
	return(x)
}

## <<References>>
## [1] Efficient sparse coding algorithms
## 	Honglak Lee, Alexis Battle, Rajat Raina, and Andrew Y. Ng.
## 	http://ai.stanford.edu/~hllee/softwares/nips06-sparsecoding.htm