clime (version 0.5.0)

clime: solve for the inverse matrix

Description

Solve for a series of the inverse covariance matrix estimates at a grid of values for the constraint lambda.

Usage

clime(x, lambda=NULL, nlambda=ifelse(is.null(lambda),100,length(lambda)),
      lambda.max=0.8, lambda.min=ifelse(nrow(x)>ncol(x), 1e-4, 1e-2),
      sigma=FALSE, perturb=TRUE, standardize=TRUE, logspaced=TRUE,
      linsolver=c("primaldual", "simplex"), pdtol=1e-3, pdmaxiter=50)

Arguments

x

Input matrix of size n (observations) times p (variables). Each column is a variable of length n. Alternatively, the sample covariance matrix may be set here with the next option sigma set to be TRUE. When the input is the sample covariance matrix, cv.clime can not be used for this object.

lambda

Grid of non-negative values for the constraint parameter lambda. If missing, nlambda values from lambda.min to lambda.max will be generated.

standardize

Whether the variables will be standardized to have mean zero and unit standard deviation. Default TRUE.

nlambda

Number of values for program generated lambda. Default 100.

lambda.max

Maximum value of program generated lambda. Default 0.8.

lambda.min

Minimum value of program generated lambda. Default 1e-4(\(n > p\)) or 1e-2(\(n < p\)).

sigma

Whether x is the sample covariance matrix. Default FALSE.

perturb

Whether a perturbed Sigma should be used or the positive perturbation added if it is numerical. Default TRUE.

logspaced

Whether program generated lambda should be log-spaced or linear spaced. Default TRUE.

linsolver

Whether primaldual (default) or simplex method should be employed. Rule of thumb: primaldual for large p, simplex for small p.

pdtol

Tolerance for the duality gap, ignored if simplex is employed.

pdmaxiter

Maximum number of iterations for primaldual, ignored if simplex is employed.

Value

An object with S3 class "clime". You can also use it as a regular R list with the following fields:

Omega

List of estimated inverse covariance matrix for a grid of values for lambda.

lambda

Actual sequence of lambda used in the program

perturb

Actual perturbation used in the program.

standardize

Whether standardization is applied to the columns of x.

x

Actual x used in the program.

lpfun

Linear programming solver used.

Details

A constrained \(\ell_1\) minimization approach for sparse precision matrix estimation (details in references) is implemented here using linear programming (revised simplex or primal-dual interior point method). It solves a sequence of lambda values on the following objective function $$ \min | \Omega |_1 \quad \textrm{subject to: } || \Sigma_n \Omega - I ||_\infty \le \lambda $$ where \(\Sigma_n\) is the sample covariance matrix and \(\Omega\) is the inverse we want to estimate.

References

Cai, T.T., Liu, W., and Luo, X. (2011). A constrained \(\ell_1\) minimization approach for sparse precision matrix estimation. Journal of the American Statistical Association 106(494): 594-607.

Examples

# NOT RUN {
## trivial example
n <- 50
p <- 5
X <- matrix(rnorm(n*p), nrow=n)
re.clime <- clime(X)

## tridiagonal matrix example
bandMat <- function(p, k) {
  cM <- matrix(rep(1:p, each=p), nrow=p, ncol=p)
  return((abs(t(cM)-cM)<=k)*1)
}
## tridiagonal Omega with diagonal 1 and off-diagonal 0.5
Omega <- bandMat(p, 1)*0.5
diag(Omega) <- 1
Sigma <- solve(Omega)
X <- matrix(rnorm(n*p), nrow=n)%*%chol(Sigma)
re.clime <- clime(X, standardize=FALSE, linsolver="simplex")
re.cv <- cv.clime(re.clime)
re.clime.opt <- clime(X, standardize=FALSE, re.cv$lambdaopt)

## Compare Frobenius norm loss
## clime estimator
sqrt( sum( (Omega-re.clime.opt$Omegalist[[1]])^2 ) )
## Not run: 0.3438533
## Sample covariance matrix inversed
sqrt( sum( ( Omega-solve(cov(X)*(1-1/n)) )^2 ) )
## Not run: 0.874041
sqrt( sum( ( Omega-solve(cov(X)) )^2 ) )
## Not run: 0.8224296
# }