# gclm

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

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 = 100, mu = rep(0,4),Sigma =  Sigma)

### Solve minimization

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

res$B  ## [,1] [,2] [,3] [,4] ## [1,] -1.012851 0.3646450 0.0000000 0.0000000 ## [2,] 0.000000 -0.5771159 0.0000000 0.0000000 ## [3,] 0.000000 0.0000000 -0.9002638 0.1671158 ## [4,] 0.000000 0.0000000 0.0000000 -0.3498257  res$C

## [1] 0.3169109 0.8481827 0.5380383 1.2108944


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

##           [,1]      [,2]      [,3]       [,4]
## [1,] -3.196011  1.150623  0.000000  0.0000000
## [2,]  0.000000 -1.821067  0.000000  0.0000000
## [3,]  0.000000  0.000000 -2.840747  0.5273274
## [4,]  0.000000  0.000000  0.000000 -1.1038614

C1

## [1] 1.000000 2.676407 1.697759 3.820930


#### Solutions along a regularization path

path <- gclm.path(cov(sample), lambdac = 0.01)
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))))

##           lambda npar fp tp fn tn errs
##  [1,] 0.85588233    4  0  4  3  9    3
##  [2,] 0.76078430    4  0  4  3  9    3
##  [3,] 0.66568626    5  0  5  2  9    2
##  [4,] 0.57058822    6  0  6  1  9    1
##  [5,] 0.47549019    6  0  6  1  9    1
##  [6,] 0.38039215    6  0  6  1  9    1
##  [7,] 0.28529411    6  0  6  1  9    1
##  [8,] 0.19019607    6  0  6  1  9    1
##  [9,] 0.09509804    7  1  6  1  8    2
## [10,] 0.00000000   16  9  7  0  0    9


## Try the gclm package in your browser

Any scripts or data that you put into this service are public.

gclm documentation built on July 1, 2020, 10:35 p.m.