knitr::opts_chunk$set(echo = TRUE)

gclm

This package contains the implementation of the algorithm in

Varando G, Hansen NR (2020) Graphical continuous Lyapunov models

gclm contains methods to estimate a sparse parametrization of covariance matrix as solution of a continuous time Lyapunov equation (CLE):

[ B\Sigma + \Sigma B^t + C = 0 ]

Solving the following (\ell_1) penalized loss minimization problem:

[ \arg\min L(\Sigma(B,C)) + \lambda \rho_1(B) + \lambda_C ||C - C_0||^2_F ]

subject to (B) stable and (C) diagonal, where (\rho_1(B)) is the (\ell_1) norm of the off-diagonal elements of (B) and (||C - C_0||^2_F ) is the squared frobenius norm of the difference between (C) and a fixed diagonal matrix ( C_0 ) (usually the identity).

Installation

## version on CRAN
install.packages("gclm")
## development version from github
devtools::install_github("gherardovarando/gclm")

Usage

library(gclm)

### define coefficient matrices
B <- matrix(nrow = 4, c(-4,  2,  0,   0, 
                         0, -3,  1,   0,
                         0,  0, -2, 0.5,
                         0,  0,  0,  -1), byrow = TRUE)
C <- diag(c(1,4,1,4))

### solve continuous Lyapunov equation 
### to obtain real covariance matrix
Sigma <- clyap(B, C) 

### obtain observations 
sample <- MASS::mvrnorm(n = 1000, mu = rep(0,4), Sigma =  Sigma)


### Solve minimization

res <- gclm(cov(sample), lambda = 0.4, lambdac = 0.01)

res$B

res$C

The CLE can be freely multiplied by a scalar and thus the (B,C) parametrization can be rescaled. As an example we can impose ( C_{11} = 1) as in the true ( C ) matrix, obtaining the estimators:

C1 <- res$C / res$C[1]
B1 <- res$B / res$C[1]

B1 
C1

Solutions along a regularization path

path <- gclm.path(cov(sample), lambdac = 0.01, 
                  lambdas = 10^seq(0, -3, length = 10))
t(sapply(path, function(res) c(lambda = res$lambda, 
                                npar = sum(res$B!=0),
                                fp = sum(res$B!=0 & B==0),
                                tp = sum(res$B!=0 & B!=0) ,
                                fn = sum(res$B==0 & B!=0),
                                tn = sum(res$B==0 & B==0),
                                errs = sum(res$B!=0 & B==0) + 
                                  sum(res$B==0 & B!=0))))

Related code



gherardovarando/clggm documentation built on April 17, 2023, 10:04 a.m.