R/potts.R

In Libra: Linearized Bregman Algorithms for Generalized Linear Models

#' Linearized Bregman solver for composite conditionally likelihood of Potts model 
#'  with lasso penalty and block-group penalty.
#' 
#' Solver for the entire solution path of coefficients.
#' 
#' The data matrix X is transformed into a 0-1 indicator matrix D with each column
#' \eqn{D_{jk}} means \eqn{1(X_j)==k}. The Potts model here used is described as following:
#' \deqn{P(x) \sim \exp(\sum_{jk} a_{0,jk}1(x_i=1) + d^T \Theta d/2)}
#' where \eqn{\Theta} is p-by-p symmetric and 0 on diagnal. Then conditional on \eqn{x_{-j}}\cr
#' \deqn{P(x_j=k) \sim exp(\sum_{k} a_{0,jk} + \sum_{i\neq j,r}\theta_{jk,ir}d_{ir})}\cr
#' then the composite conditional likelihood is like this:\cr
#' \deqn{- \sum_{j} condloglik(X_j | X_{-j})}
#' 
#' @param X An n-by-p matrix of variables.
#' @param kappa The damping factor of the Linearized Bregman Algorithm that is
#'  defined in the reference paper. See details. 
#' @param alpha Parameter in Linearized Bregman algorithm which controls the 
#' step-length of the discretized solver for the Bregman Inverse Scale Space. 
#' See details. 
#' @param c Normalized step-length. If alpha is missing, alpha is automatically generated by 
#' \code{alpha=c*n/(kappa*||XX^T*XX||_2)}, where XX is 0-1 indicator matrix 
#' induced by the class of each Xi. Default is 1. It should be in (0,2).
#' If beyond this range the path may be oscillated at large t values.
#' @param intercept if TRUE, an intercept is included in the model (and not 
#' penalized), otherwise no intercept is included. Default is TRUE.
#' @param tlist Parameters t along the path.
#' @param nt Number of t. Used only if tlist is missing. Default is 100.
#' @param trate tmax/tmin. Used only if tlist is missing. Default is 100.
#' @param print If TRUE, the percentage of finished computation is printed.
#' @param group Whether to use a block-wise group penalty, Default is FALSE
#' @return A "potts" class object is returned. The list contains the call, 
#'  the path, the intercept term a0 and value for alpha, kappa, t.
#' @author Jiechao Xiong
#' @keywords regression
#' @examples
#'
#' X = matrix(floor(runif(200*10)*3),200,10)
#' obj = potts(X,10,nt=100,trate=10,group=TRUE)
#'


potts <- function(X, kappa, alpha,c=1, tlist,nt = 100,trate = 100,group=FALSE, intercept = TRUE,print=FALSE) 
{
  if (!is.matrix(X)) stop("X must be a matrix!")
  np <- dim(X)
  n <- np[1]
  p0 <- np[2]
  XX <- matrix(nrow=n,ncol=0)
  num <- rep(1,p0)
  for (i in 1:p0){
    x_unique <- unique(X[,i])
    num[i] <- length(x_unique)
    XX <- cbind(XX,sapply(1:length(x_unique),function(x) as.numeric(X[,i]==x_unique[x])))
  }
  group_split <- c(0, cumsum(num))
  group_split_length <- length(group_split)
  p <- group_split[group_split_length]
  
  if (missing(tlist)) tlist<-rep(-1.0,nt)
  else nt=length(tlist)
  
  if (missing(alpha)){
    meanx <- drop(rep(1,n) %*% XX)/n
    tempX <- scale(XX, meanx, FALSE)
    sigma <- norm(tempX,"2")
    alpha <- c*n/kappa/sigma^2
  }
  group <- as.integer(group != 0)
	intercept <- as.integer(intercept != 0)
	result_r <- vector(length = nt * p*(p + intercept))
	solution <- .C("potts_C",
	               as.numeric(XX),
	               as.integer(n),
	               as.integer(p),
	               as.numeric(kappa),
	               as.numeric(alpha),
	               as.numeric(result_r),
	               as.integer(group_split),
	               as.integer(group_split_length),         
	               as.integer(intercept),
                 as.numeric(tlist),
	               as.integer(nt),
	               as.numeric(trate),
	               as.integer(group),
	               as.integer(print!=0)
	               )
	path <- sapply(0:(nt- 1), function(x)
	  matrix(solution[[6]][(1+x*p*(p+intercept)):((x+1)*p*(p+intercept))], p, p+intercept),simplify = "array")
	
	if (intercept){
	  a0 <- path[,p+1,,drop=TRUE]
	  path <- path[,-(p+1),,drop=FALSE]
	}else{
	  a0 <- matrix(0,p,length(nt))
	}
	object <- list(call = call, family="potts", kappa = kappa, alpha = alpha, path = path, nt = nt,t=solution[[10]],a0=a0)
	class(object) <- "potts"
	object
}