Type: Package
Title: Graphical Continuous Lyapunov Models
Version: 0.0.1
Description: Estimation of covariance matrices as solutions of continuous time Lyapunov equations. Sparse coefficient matrix and diagonal noise are estimated with a proximal gradient method for an l1-penalized loss minimization problem. Varando G, Hansen NR (2020) <doi:10.48550/arXiv.2005.10483>.
License: MIT + file LICENSE
Encoding: UTF-8
LazyData: true
RoxygenNote: 7.1.0
URL: https://github.com/gherardovarando/gclm
BugReports: https://github.com/gherardovarando/gclm/issues
Suggests: testthat
NeedsCompilation: yes
Packaged: 2020-05-25 06:55:33 UTC; gherardo
Author: Gherardo Varando ORCID iD [aut, cre, cph], Niels Richard Hansen ORCID iD [aut]
Maintainer: Gherardo Varando <gherardo.varando@gmail.com>
Repository: CRAN
Date/Publication: 2020-06-04 08:40:07 UTC

Generate a naive stable matrix

Description

Generate a naive stable matrix

Usage

B0(p)

Arguments

p

dimension of the matrix

Value

a stable matrix with off-diagonal entries equal to 1 and diagonal entries equal to -p


Anti transpose (internal)

Description

Anti transpose (internal)

Usage

anti_t(m)

Arguments

m

square matrix

Value

the anti transpose of m


Solve continuous-time Lyapunov equations

Description

clyap solve the continuous-time Lyapunov equations

BX + XB' + C=0

Using the Bartels-Stewart algorithm with Hessenberg–Schur decomposition. Optionally the Hessenberg-Schur decomposition can be returned.

Usage

clyap(B, C, Q = NULL, all = FALSE)

Arguments

B

Square matrix

C

Square matrix

Q

Square matrix, the orthogonal matrix used to transform the original equation

all

logical

Details

If the matrix Q is set then the matrix B is assumed to be in upper quasi-triangular form (Hessenberg-Schur canonical form), as required by LAPACK subroutine DTRSYL and Q is the orthogonal matrix associated with the Hessenberg-Schur form of B. Usually the matrix Q and the appropriate form of B are obtained by a first call to clyap(B, C, all = TRUE)

clyap uses lapack subroutines:

Value

The solution matrix X if all = FALSE. If all = TRUE a list with components X, B and Q. Where B and Q are the Hessenberg-Schur form of the original matrix B and the orthogonal matrix that performed the transformation.

Examples

B <- matrix(data = rnorm(9), nrow = 3)
## make B negative diagonally dominant, thus stable:
diag(B) <- - 3 * max(B) 
C <- diag(runif(3))
X <- clyap(B, C)
## check X is a solution:
max(abs(B %*% X + X %*% t(B) + C)) 

l1 penalized loss estimation for GCLM

Description

Estimates a sparse continuous time Lyapunov parametrization of a covariance matrix using a lasso (L1) penalty.

Usage

gclm(
  Sigma,
  B = -0.5 * diag(ncol(Sigma)),
  C = rep(1, ncol(Sigma)),
  C0 = rep(1, ncol(Sigma)),
  loss = "loglik",
  eps = 0.01,
  alpha = 0.5,
  maxIter = 100,
  lambda = 0,
  lambdac = 0,
  job = 0
)

gclm.path(
  Sigma,
  lambdas = NULL,
  B = -0.5 * diag(ncol(Sigma)),
  C = rep(1, ncol(Sigma)),
  ...
)

Arguments

Sigma

covariance matrix

B

initial B matrix

C

diagonal of initial C matrix

C0

diagonal of penalization matrix

loss

one of "loglik" (default) or "frobenius"

eps

convergence threshold

alpha

parameter line search

maxIter

maximum number of iterations

lambda

penalization coefficient for B

lambdac

penalization coefficient for C

job

integer 0,1,10 or 11

lambdas

sequence of lambda

...

additional arguments passed to gclm

Details

gclm performs proximal gradient descent for the optimization problem

argmin L(\Sigma(B,C)) + \lambda \rho(B) + \lambda_C ||C - C0||_F^2

subject to B stable and C diagonal, where \rho(B) is the l1 norm of the off-diagonal element of B.

gclm.path simply calls iteratively gclm with different lambda values. Warm start is used, that is in the i-th call to gclm the B and C matrices are initialized as the one obtained in the (i-1)th call.

Value

for gclm: a list with the result of the optimization

for gclm.path: a list of the same length of lambdas with the results of the optimization for the different lambda values

Examples

x <- matrix(rnorm(50*20),ncol=20)
S <- cov(x)

## l1 penalized log-likelihood
res <- gclm(S, eps = 0, lambda = 0.1, lambdac = 0.01)

## l1 penalized log-likelihood with fixed C
res <- gclm(S, eps = 0, lambda = 0.1, lambdac = -1)

## l1 penalized frobenius loss
res <- gclm(S, eps = 0, lambda = 0.1, loss = "frobenius")

Recover lower triangular GCLM

Description

Recover the only lower triangular stable matrix B such that Sigma (\Sigma) is the solution of the associated continuous Lyapunov equation:

B\Sigma + \Sigma B' + C = 0

Usage

gclm.lowertri(Sigma, P = solve(Sigma), C = diag(nrow = nrow(Sigma)))

Arguments

Sigma

covariance matrix

P

the inverse of the covariance matrix

C

symmetric positive definite matrix

Value

A stable lower triangular matrix


List all the inverse matrices of the leading sub-matrix of P=Sigma^-1 (INTERNAL)

Description

List all the inverse matrices of the leading sub-matrix of P=Sigma^-1 (INTERNAL)

Usage

listInverseBlocks(Sigma)

Arguments

Sigma

an invertible matrix

Value

a list with the inverse of the leading sub-matrices of Sigma^(-1)