README.Rmd

---
output: github_document
---

```{r setup, include=FALSE}
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](https://arxiv.org/abs/2005.10483) 

`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 

```{r, eval = FALSE}
## version on CRAN
install.packages("gclm")
## development version from github
devtools::install_github("gherardovarando/gclm")
```


## Usage

```{r}
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:

```{r}
C1 <- res$C / res$C[1]
B1 <- res$B / res$C[1]

B1 
C1
```


#### Solutions along a regularization path

```{r}
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

- Some inspiration is from the `lyapunov` package (https://github.com/gragusa/lyapunov).