R/ggm.R

In Libra: Linearized Bregman Algorithms for Generalized Linear Models

#' Linearized Bregman solver for composite conditionally likelihood of Gaussian Graphical
#'  model with lasso penalty.
#' 
#' Solver for the entire solution path of coefficients.
#' 
#' The data matrix X is assumed to follow the Gaussian Graohical model which is described as following:\cr
#' \deqn{X \sim N(\mu, \Theta^{-1})}\cr
#' where \eqn{\Theta} is sparse p-by-p symmetric matrix. Then conditional on \eqn{x_{-j}}\cr
#' \deqn{x_j \sim N(\mu_j - \sum_{k\neq j}\Theta_{jk}/\Theta_{jj}(x_k-\mu_k),1/\Theta_{jj}) }\cr
#' then the composite conditional likelihood is like this:\cr
#' \deqn{- \sum_{j} condloglik(X_j | X_{-j})}\cr
#' or in detail:\cr
#' \deqn{\sum_{j} \Theta_{j}^TS\Theta_{j}/2\Theta_{jj} - ln(\Theta_{jj})/2}\cr
#' where \eqn{S} is covariance matrix of data. It is easy to prove that this loss function
#' is convex.
#' 
#' 
#' @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 S The covariance matrix can be provided directly if data matrix X is missing.
#' @param c Normalized step-length. If alpha is missing, alpha is automatically generated by 
#' \code{alpha=c*n/(kappa*||X^T*X||_2)}. Default is 2. It should be in (0,4).
#' If beyond this range the path may be oscillated at large t values.
#' @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.
#' @return A "ggm" class object is returned. The list contains the call, 
#'  the path, value for alpha, kappa, t.
#' @author Jiechao Xiong
#' @keywords regression
#' @examples
#' 
#' library(MASS)
#' p = 20
#' Omega = diag(1,p,p)
#' Omega[0:(p-2)*(p+1)+2] = 1/3
#' Omega[1:(p-1)*(p+1)] = 1/3
#' S = solve(Omega)
#' X = mvrnorm(n=500,rep(0,p),S)
#' obj = ggm(X,10,trate=10)
#' obj$path[,,50]

ggm <- function(X, kappa, alpha, S=NA, c = 2, tlist,nt = 100,trate = 100,print=FALSE) 
{
  if (!missing(X)){
    n <- nrow(X)
    S <- cov(X)*(1-1/n)
  }else if(any(is.na(S))){stop('A data matrix or covariance matrix must be provided!')}
  p <- nrow(S)
  
  if (missing(tlist)) tlist<-rep(-1.0,nt)
  else nt<-length(tlist)
  
  if (missing(alpha)){
    sigma <- norm(S,"2")
    alpha <- c/kappa/sigma^2
  }
  
	result_r <- vector(length = nt*p^2)
	solution <- .C("ggm_C",
		as.numeric(S),
		as.integer(p),
		as.numeric(kappa),
		as.numeric(alpha),
		as.numeric(result_r),
		as.numeric(tlist),
		as.integer(nt),
		as.numeric(trate),
		as.integer(print!=0)
	)
	path <- sapply(0:(nt- 1), function(x)
	  matrix(solution[[5]][(1+x*p^2):((x+1)*p^2)], p, p),simplify = "array")

	object <- list(call = call,family="ggm", kappa = kappa, alpha = alpha, path = path, nt = nt,t=solution[[6]])
	class(object) <- "ggm"
	object
}