An Introduction to cxreg

Introduction

is a package that fits complex-valued penalized regressions (Lasso) and Gaussian likelihoods (graphical Lasso) using an exact pathwise coordinate descent method. Similar to the well-known package (2010) for Lasso in real value settings, computes the regularization path for complex-valued linear regression with a lasso penalty, referred to as , over a grid of values for the regularization parameter lambda. Likewise, to fit complex-valued Graphical Lasso models (2008), the regularization path is computed via . The package includes methods for prediction, cross-validation functions, printing and plotting utilities, as well as fitting for both model types.

The authors of are Navonil Deb and Sumanta Basu, with contributions from Younghoon Kim. The R package is maintained by Younghoon Kim.

This vignette describes basic usage of functions related to and models in package in R. There is an original description of the algorithm that should be useful:

• “Regularized estimation of sparse spectral precision matrices”

solves the following problem: For \(Y\in\mathbb{C}^n\) and \(X\in\mathbb{C}^{n \times p}\),

\[ \min_{\beta \in \mathbb{C}^p} \frac{1}{2n} \|Y - X \beta\|_2^2 + \lambda\|\beta\|_1, \]

over a grid of values of \(\lambda\) covering the entire range of possible solutions. Here, since the absolute of complex-values means a complex modulus, \(\beta\in\mathbb{C}^p\) so that the \(\ell_1\)-penalty can be viewed as

\[ \|\beta\|_1 = \sum_{j=1}^p |\beta_j| = \sum_{j=1}^p\|\mathrm{Re}(\beta_j) + \mathrm{Im}(\beta_j)\|_2, \]

which is a group Lasso with \(p\) groups, each of size 2.

Similarly, seeks the minimizer of Whittle’s approximate likelihood (1951). For a \(p\)-dimensional Gaussian time series \({X_t}_{t=1}^n\), one can compute the \(p\)-dimensional complex-valued discrete Fourier transforms (DFT) at the Fourier frequencies \(\omega_j = 2\pi j / n\), where \(j \in F_n := \left\{- \left[\frac{n-1}{2}\right], \ldots, \left[\frac{n}{2}\right] \right\}\):

\[ d_j := d(\omega_j) = \frac{1}{\sqrt{n}}\sum_{t=1}^n X_t \exp(-i t \omega_j). \]

Note that \(d_j \sim \mathcal{N}_{\mathbb{C}}(0, f(\omega_j))\), where \(f(\omega_j)\) is the spectral density, commonly estimated by averaging the periodogram \(I(\omega_j) = d(\omega_j)d^{\dagger}(\omega_j)\) over a bandwidth of \(2m+1\) frequencies:

\[ \hat{f}(\omega_j) = \frac{1}{2\pi(2m+1)}\sum_{|k|\leq m}I(\omega_{j+k}), \quad j\in F_n. \tag{1} \]

The approximation to the negative log-likelihood then takes the form:

\[ \sum_{k=j-m}^{j+m} \log\det f^{-1}(\omega_j) - \sum_{k=j-m}^{j+m} d_{k}^{\dagger} f^{-1}(\omega_j)d_k. \] Finally, by rearranging the components used to construct \(\hat{f}(\omega_j)\) and adding an \(\ell_1\)-penalty on the off-diagonal entries of the inverse spectral density, \(\|\Theta\|_{1,\mathrm{off}} = \sum_{k \neq \ell} |\Theta_{k,\ell}|=\sum_{k \neq \ell}\|\mathrm{Re}(\Theta_{k,\ell}) + \mathrm{Im}(\Theta_{k,\ell})\|_2\), one obtains the estimator:

\[ \hat{\Theta}_j := \hat{\Theta}(\omega_j) = \arg\min_{\Theta \in \mathcal{H}_{++}^p} \left\{ \mathrm{trace}(\Theta\hat{f}(\omega_j)) - \log\det\Theta + \lambda \|\Theta\|_{1,\mathrm{off}} \right\}, \] where \(\mathcal{H}_{++}^p\) represents the set of \(p \times p\) symmetric positive definite matrices. Throughout the illustrative example in this document, we drop the subscript \(j\) and focus on deriving \(\hat{\Theta}\) for a given \(\hat{f}\).

From a numerical standpoint, there are two formulations for complex-valued graphical Lasso (see Section 5 in (2024)), where a similar argument can be made in standard graphical Lasso (e.g., (2018)).

: The first formulation is to solve the complex-valued graphical Lasso with scaled spectral density matrix (called spectral coherence), \(\hat{R} = D^{-1}\hat{f}D^{-1}\) with \(D^2 = \mathrm{diag}(\hat{f}_{1,1},\ldots,\hat{f}_{p,p})\), and scale back once the estimates are obtained. The corresponding optimization problem is

\[ \hat{\Theta} = D^{-1} \hat{K} D^{-1}, \quad \hat{K} = \arg\min_{\Theta\in \mathcal{H}_{++}^p} \left\{ \mathrm{trace}(\hat{R} \Theta) - \log\det \Theta + \lambda \|\Theta\|_{\mathrm{1,off}}\right\}. \tag{CGLASSO-sc1} \]

This is equivalent to a equation (9) in (2018) for the graphical lasso with weighted penalty. The entry-wise weights determined by \(D\) as \(\lambda_{k,\ell} \propto D_{k,k}D_{\ell,\ell} = \sqrt{ \hat{f}_{k,k} \hat{f}_{\ell,\ell}}\).

: The second formulation takes the spectral density as the original input. In each iteration of within Algorithm 3 in (2024), we scale the partitioned \(W\), working inverse spectral density, with the corresponding diagonal entries of the estimated spectral density. The Lasso outputs are scaled back to obtain rows and columns of \(W\). Specifically, we implement the following update for \(k\)th row and column:

\[ \begin{aligned} W_{11}^{\mathrm{scl}} &\leftarrow D_{11}^{-1} W_{11} D_{11}^{-1}, \quad D_{11}^2 = \mathrm{diag}(W_{11}) = \mathrm{diag}(P_{11}), \quad \mathbf{p}_{12}^{\mathrm{scl}} \leftarrow D_{11}^{-1} \mathbf{p}_{12}, \\ \hat{\beta}^{\mathrm{scl}} &= \mathrm{CLASSO.COV}(W_{11}^{\mathrm{scl}}, \mathbf{p}_{12}^{\mathrm{scl}}, \mathcal{B}^{(0)}_{\cdot, k}, \lambda), \\ \mathcal{B}^{(0)}_{\cdot,k} &\leftarrow \hat{\beta}^{\mathrm{scl}}, \quad \hat{\beta} \leftarrow D_{11}^{-1} \hat{\beta}^{\mathrm{scl}}, \quad \mathbf{w}_{12} \leftarrow W_{11} \hat{\beta}, \end{aligned} \tag{CGLASSO-sc2} \]

where the scaling matrix \(D_{11}\) is similar to \(D\) is CGLASSO-sc1, except that the last diagonal entry \(D_{p,p}\) is disregarded in every update due to absence of the last row and column.

Installation

Like other R packages on Github, it can be downloaded by the command:

library(devtools)
devtools::install_github("yk748/cxreg")

Example: classo

The purpose of this section is to give users a general sense of the package regarding classo. We will briefly go over the main functions, basic operations and outputs. can be loaded using the library command:

library(cxreg)

We load a set of data created beforehand for illustration:

data(classo_example)
x <- classo_example$x
y <- classo_example$y

Note that “x” is already standardized, which makes the columns in \(X\) orthogonal. The provides the orthogonalization, and is the default. The necessity of the standardization is described in the original paper (2024).

We fit the model using the most basic call to . In addition to standardization, the is used as a default, which means the complex-valued constant is not considered. This is another difference from the package.

fit <- classo(x,y)

is an object of class that contains all the relevant information of the fitted model for further use. Various methods are provided for the object such as , , and , just like package.

We can visualize the coefficients by executing the ‘plot’ method. First, we plot against the log \(\lambda\) values with labels:

plot(fit, xvar="lambda", label=TRUE)

There are more types of \(x\)-axis, similar to those in . If “norm” is chosen, the values of the \(x\)-axis become \(\ell_1\) (modulus) of the coefficients:

plot(fit, xvar="norm", label=TRUE)

When “dev” is chosen, the percentage deviance explained by the model is used. Similar to those in , we can observe the blowing up phenomenon at the end of the curves:

plot(fit, xvar="dev", label=TRUE)

Unlike , the real (Re) and imaginary (Im) parts of the coefficients are displayed separately. It shows the path of its coefficients in the two parts separately against the \(\ell_1\)-norm of the whole coefficient vector as \(\lambda\) varies. Users may also wish to annotate the curves: this can be done by setting “label = TRUE” in the plot command.

We can obtain the model coefficients at one or more \(\lambda\)’s within the range of the sequence. Similar to , if \(s\) (lambda value) is not included in the sequence, only approximate values are provided:

any(fit$lambda == 0.1)

## [1] FALSE

coef(fit, s=0.1, exact=FALSE)

##                            s1
## V1   0.40105529-0.2972178861i
## V2   0.98902875+1.9546441978i
## V3   0.10351678-0.0312120969i
## V4   0.00000000+0.0000000000i
## V5   0.00000000+0.0000000000i
## V6   0.02281945-0.0699368612i
## V7   0.00000000+0.0000000000i
## V8  -0.11289043+0.0262567645i
## V9   0.00000000+0.0000000000i
## V10  0.08413426-0.0002130587i
## V11 -0.03891600+0.0457059172i
## V12  0.07753989-0.0101036782i
## V13 -0.05600043-0.0055142505i
## V14  0.00000000+0.0000000000i
## V15 -0.08646406+0.0062733357i
## V16  0.00000000+0.0000000000i
## V17  0.00000000+0.0000000000i
## V18  0.00000000+0.0000000000i
## V19 -0.09367033-0.1051484558i
## V20 -0.03736639-0.0051862273i

However, when “exact = TRUE”, the data should be provided to create the original fit, as it merged to the original fit. In this case, both \(x\) and \(y\) need to be supplied as named arguments:

coef(fit, s=0.1, exact=TRUE, x=x, y=y)

##                            s1
## V1   0.40105217-0.2972192512i
## V2   0.98902283+1.9546408619i
## V3   0.10351114-0.0312224080i
## V4   0.00000000+0.0000000000i
## V5   0.00000000+0.0000000000i
## V6   0.02283305-0.0699231880i
## V7   0.00000000+0.0000000000i
## V8  -0.11288666+0.0262662482i
## V9   0.00000000+0.0000000000i
## V10  0.08412604-0.0001639598i
## V11 -0.03891491+0.0456985714i
## V12  0.07753673-0.0101048978i
## V13 -0.05599224-0.0055307140i
## V14  0.00000000+0.0000000000i
## V15 -0.08645909+0.0062757702i
## V16  0.00000000+0.0000000000i
## V17  0.00000000+0.0000000000i
## V18  0.00000000+0.0000000000i
## V19 -0.09367530-0.1051358624i
## V20 -0.03736433-0.0051727268i

Users can also make predictions at specific \(\lambda\)’s with new input data: Note that the third line in the following chunk is about standardization, whose purpose is mentioned above. Here, ‘response’ means that predicted values of \(y\).

set.seed(29)
nx <- array(rnorm(5*20), c(5,20)) + (1+1i) * array(rnorm(5*20), c(5,20))
for (j in 1:20) {
  nx[,j] <- nx[,j] / sqrt(mean(Mod(nx[,j])^2))
}
predict(fit, newx = nx, s = c(0.1, 0.05), type="response")

##                         s1                    s2
## [1,] -0.32481930+3.645038i -0.39629864+3.704465i
## [2,] -0.57391307+1.795762i -0.58827334+2.036638i
## [3,] -1.17694639-1.152282i -1.27970611-1.289929i
## [4,]  0.32383314-1.768570i  0.43415839-1.844754i
## [5,]  0.06312868+2.285649i  0.04910468+2.183336i

The other two output types are and , which are analogous to those in :

predict(fit, newx = nx, s = c(0.1, 0.05), type="coefficients")

##                            s1                        s2
## V1   0.40105529-0.2972178861i  0.428214371-0.318935563i
## V2   0.98902875+1.9546441978i  0.999660784+2.016798885i
## V3   0.10351678-0.0312120969i  0.139385420-0.040227509i
## V4   0.00000000+0.0000000000i -0.003254420+0.001500173i
## V5   0.00000000+0.0000000000i  0.015771136-0.014275686i
## V6   0.02281945-0.0699368612i  0.021275606-0.106942413i
## V7   0.00000000+0.0000000000i  0.006034593-0.056865053i
## V8  -0.11289043+0.0262567645i -0.168151788+0.040938458i
## V9   0.00000000+0.0000000000i  0.010278116-0.043948609i
## V10  0.08413426-0.0002130587i  0.164869520-0.030198139i
## V11 -0.03891600+0.0457059172i -0.072647633+0.085904684i
## V12  0.07753989-0.0101036782i  0.121524579-0.015692437i
## V13 -0.05600043-0.0055142505i -0.155899766-0.006893521i
## V14  0.00000000+0.0000000000i  0.000000000+0.000000000i
## V15 -0.08646406+0.0062733357i -0.119223730+0.011067907i
## V16  0.00000000+0.0000000000i -0.028189795-0.027876879i
## V17  0.00000000+0.0000000000i -0.001475252+0.023214736i
## V18  0.00000000+0.0000000000i  0.007287146-0.030802916i
## V19 -0.09367033-0.1051484558i -0.112400477-0.120328696i
## V20 -0.03736639-0.0051862273i -0.089875425-0.017883726i

predict(fit, newx = nx, s = c(0.1, 0.05), type="nonzero")

## $s1
##  [1]  1  2  3  6  8 10 11 12 13 15 19 20
## 
## $s2
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 15 16 17 18 19 20

The function returns a sequence of models for the users to choose from. In many cases, users may prefer the software to select one of them. Cross-validation is perhaps the simplest and most widely used method for that task. is the main function to do cross-validation here, along with various supporting methods such as plotting and prediction.

cvfit <- cv.classo(x,y,trace.it = 1)

## Training
## Fold: 1/10
## Fold: 2/10
## Fold: 3/10
## Fold: 4/10
## Fold: 5/10
## Fold: 6/10
## Fold: 7/10
## Fold: 8/10
## Fold: 9/10
## Fold: 10/10

print(cvfit)

## 
## Call:  cv.classo(x = x, y = y, trace.it = 1) 
## 
## Measure: Mean-Squared Error 
## 
##     Lambda Index Measure     SE Nonzero
## min 0.2748    24   3.054 0.4036       2
## 1se 0.5270    17   3.374 0.3948       2

returns a ‘cv.classo’ object, a list with all the ingredients of the cross-validated fit.

plot(cvfit)

This plots the cross-validation curve (red dotted line) along with upper and lower standard deviation curves along the \(\lambda\) sequence (error bars). Two special values along the \(\lambda\) sequence are indicated by the vertical dotted lines. “lambda.min” is the value of \(\lambda\) that gives minimum mean cross-validated error, while “lambda.lse” is the value of \(\lambda\) that gives the most regularized model such that the cross-validated error is within one standard error of the minimum.

We can use the following code to get the value of “lambda.min” and the model coefficients at that value of \(\lambda\):

cvfit$lambda.min

## [1] 0.2747633

coef(cvfit, s = "lambda.min")

##               lambda.min
## V1  0.2959829-0.2280322i
## V2  0.9421236+1.8129565i
## V3  0.0000000+0.0000000i
## V4  0.0000000+0.0000000i
## V5  0.0000000+0.0000000i
## V6  0.0000000+0.0000000i
## V7  0.0000000+0.0000000i
## V8  0.0000000+0.0000000i
## V9  0.0000000+0.0000000i
## V10 0.0000000+0.0000000i
## V11 0.0000000+0.0000000i
## V12 0.0000000+0.0000000i
## V13 0.0000000+0.0000000i
## V14 0.0000000+0.0000000i
## V15 0.0000000+0.0000000i
## V16 0.0000000+0.0000000i
## V17 0.0000000+0.0000000i
## V18 0.0000000+0.0000000i
## V19 0.0000000+0.0000000i
## V20 0.0000000+0.0000000i

To get the corresponding values at “lambda.1se”, simply replace “lambda.min” with “lambda.1se” above, or omit the “s” argument, since “lambda.1se” is the default.

Note that unlike package, the coefficients are not represented in sparse matrix format, rather they are in the dense format. This is because of the traits of complex values in the function.

Predictions can be made based on the fitted object as well. The code below gives predictions for the new input matrix “newx” at “lambda.min”:

predict(cvfit, newx = x[1:5,], s = "lambda.min")

##                  lambda.min
## [1,]  0.2364879+0.07173343i
## [2,] -0.3353579-0.52554252i
## [3,]  0.5016216-1.41923947i
## [4,]  0.1210804+0.71063784i
## [5,] -0.4024407+3.42377148i

This concludes the basic usage of .

Example: cglasso

In this section, we illustrate our function for the two variants described in Section above. Again, we load a set of data created beforehand for illustration:

data(cglasso_example)
f_hat <- cglasso_example$f_hat
n     <- cglasso_example$n
m     <- floor(sqrt(n))   # half-bandwidth used to compute f_hat

where the number of variables and sample size used in this example are \(p=30\) and \(n=500\), respectively, and the covariance matrix of the white noise process \(\{X_t\}_{t=1}^n\) is

\[ \Sigma = C^{-1}=(C_{k \ell})^{-1}, \quad C_{kk} = 0.7, \ C_{k,k-1}=C_{k-1,k} = 0.3. \] Then the estimated spectral density \(\hat{f}\) is obtained by DFT in (1).

First, consider CGLASSO-sc1: It is called by using “type=”I”“.

fit_cglasso_I <- cglasso(S=f_hat, m=m, type="I")

## The algorithm was terminated at 17 th lambda

Now, let’s look at the example of CGLASSO-sc2:

fit_cglasso_II <- cglasso(S=f_hat, m=m, type="II", nlambda=30, stop_criterion = "AIC")

## The algorithm was terminated at 30 th lambda

Here, we change the stopping criterion from “EBIC” to “AIC”. Although the maximum number of paramters “nlambda” is set 30, due to the early termination, not all 30 lambdas created inside the function was not explored. However, if “stopping_rule=FALSE”, cases of all 50 lambdas (by default) are explored:

fit_cglasso_another <- cglasso(S=f_hat, m=m, type="II", stopping_rule = FALSE)
fit_cglasso_another$lambda_grid

##  [1] 0.596381490 0.542884984 0.494187211 0.449857718 0.409504661 0.372771347
##  [7] 0.339333078 0.308894283 0.281185903 0.255963015 0.233002666 0.212101902
## [13] 0.193075975 0.175756708 0.159991012 0.145639526 0.132575395 0.120683141
## [19] 0.109857643 0.100003213 0.091032742 0.082866939 0.075433624 0.068667090
## [25] 0.062507527 0.056900488 0.051796411 0.047150178 0.042920721 0.039070654
## [31] 0.035565946 0.032375615 0.029471464 0.026827820 0.024421315 0.022230679
## [37] 0.020236547 0.018421291 0.016768868 0.015264670 0.013895401 0.012648959
## [43] 0.011514324 0.010481468 0.009541261 0.008685393 0.007906297 0.007197088
## [49] 0.006551496 0.005963815

Finally, we can plot the heatmap using the list of estimated spectral precision matrices and specify the index to select a particular matrix from the list. The argument determines whether the display shows the real part (), imaginary part (), both (using . in parallel), or the modulus ().

plot(fit_cglasso_I, index=fit_cglasso_I$min_index, type="mod",       label=TRUE)

plot(fit_cglasso_I, index=fit_cglasso_I$min_index, type="both",      label=TRUE)

plot(fit_cglasso_II, index=fit_cglasso_II$min_index, type="real",    label=FALSE)

plot(fit_cglasso_II, index=fit_cglasso_II$min_index, type="imaginary", label=FALSE)

This concludes the basic usage of .

Example: Sparse Spectral Precision Matrix Estimation

As the real-data example from the time series, the workflow using this package should cover these steps:

library(mvtnorm)
set.seed(1010)
p <- 10    # number of variables
n <- 500   # sample size

# sparse precision matrix (tridiagonal)
C <- diag(0.7, p)
C[row(C) == col(C) + 1] <- 0.3
C[row(C) == col(C) - 1] <- 0.3
Sigma <- solve(C)

# generate multivariate Gaussian time series
X_t <- rmvnorm(n = n, mean = rep(0, p), sigma = Sigma)

# compute DFT at a specific Fourier frequency j
m <- floor(sqrt(n))   # bandwidth
j <- 1                # frequency index
d_j <- dft.j(X_t, j, m)

# estimate spectral density via smoothed periodogram
f_j_hat <- t(d_j) %*% Conj(d_j) / (2*m + 1)

fit_sim <- cglasso(S = f_j_hat, m = m, type = "I")

## The algorithm was terminated at 17 th lambda

plot(fit_sim,
     index = fit_sim$min_index,
     type  = "mod",
     label = TRUE)

# selected lambda
fit_sim$lambda_grid[fit_sim$min_index]

## [1] 0.2516412

# number of nonzero off-diagonal entries
Theta_est <- fit_sim$Theta_list[[fit_sim$min_index]]
sum(Mod(Theta_est[row(Theta_est) != col(Theta_est)]) > 1e-6)

## [1] 30