| Title: | Estimation of Low Rank Plus Sparse Structured Vector Auto-Regressive (VAR) Model |
|---|---|
| Description: | Implementations of estimation algorithm of low rank plus sparse structured VAR model by using Fast Iterative Shrinkage-Thresholding Algorithm (FISTA). It relates to the algorithm in Sumanta, Li, and Michailidis (2019) <doi:10.1109/TSP.2018.2887401>. |
| Authors: | Peiliang Bai [aut, cre] |
| Maintainer: | Peiliang Bai <[email protected]> |
| License: | GPL-2 |
| Version: | 1.2 |
| Built: | 2024-11-17 03:54:42 UTC |
| Source: | https://github.com/peiliangbai92/lsvar |
Main loss function
f.func(x, A, b)f.func(x, A, b)
x |
Model parameters |
A |
Design matrix with size of n by p |
b |
Correspond vector or matrix |
Value of objective function
A function to solve low rank plus sparse model estimation
fista.LpS( data, lambda, mu, alpha_L = 0.25, niter = 100, backtracking = TRUE, x.true = NULL )fista.LpS( data, lambda, mu, alpha_L = 0.25, niter = 100, backtracking = TRUE, x.true = NULL )
data |
A numeric dataset with size of n by p |
lambda |
A positive numeric value, indicating the tuning parameter for sparse component |
mu |
A positive numeric value, indicating the tuning parameter for low rank component |
alpha_L |
The constraint coefficient of low rank component, default is 0.25 |
niter |
The maximum number of iterations required for FISTA |
backtracking |
A boolean argument, indicating that use backtracking in the FISTA |
x.true |
A p by p matrix, the true model parameter. Only available for simulation. |
A S3 object of class LSVAR, including
estimated model parameter
Estimated sparse component
Estimated low-rank component
Values of objective function
Relative errors compared with the true model parameters if available
n <- 300 p <- 20 try <- testVAR(n, p, struct = "LS", signal = 0.75, rank = 2, singular_vals = c(1, 0.8)) data <- as.matrix(try$series) lambda <- 0.1; mu <- 1 fit <- fista.LpS(data, lambda = lambda, mu = mu, x.true = try$model_param) summary(fit, threshold = 0.2)n <- 300 p <- 20 try <- testVAR(n, p, struct = "LS", signal = 0.75, rank = 2, singular_vals = c(1, 0.8)) data <- as.matrix(try$series) lambda <- 0.1; mu <- 1 fit <- fista.LpS(data, lambda = lambda, mu = mu, x.true = try$model_param) summary(fit, threshold = 0.2)
Gradient function of quardratic loss
gradf.func(x, AtA, Atb)gradf.func(x, AtA, Atb)
x |
A vector, or matrix, indicating the model parameter |
AtA |
A p by p Gram matrix for corresponding design matrix A |
Atb |
An inner product for design matrix A and corresponding matrix (vector) b |
Value of gradients
Nuclear norm penalty for low-rank component
nuclear.pen(x, lambda)nuclear.pen(x, lambda)
x |
Model parameter |
lambda |
Tuning parameter |
Value of nuclear norm penalty term
objective function, main loss function and penalties
obj.func(x.lr, x.sparse, A, b, lambda, mu)obj.func(x.lr, x.sparse, A, b, lambda, mu)
x.lr |
low-rank component |
x.sparse |
sparse component |
A |
design matrix |
b |
correspond vector |
lambda |
a tuning parameter for sparse component |
mu |
a tuning parameter for low-rank component |
value of objective function
Plot a network to illustrate the estimated sparse component
plot_network(mat, threshold = 0.1)plot_network(mat, threshold = 0.1)
mat |
a p by p matrix, indicating the sparse component |
threshold |
the threshold for presenting the edges in the network |
A network plot for the sparse component
set.seed(1) est_mats <- matrix(rnorm(400, 0, 1), 20, 20) plot_network(est_mats, threshold = 0.1)set.seed(1) est_mats <- matrix(rnorm(400, 0, 1), 20, 20) plot_network(est_mats, threshold = 0.1)
Proximal function with nuclear norm
prox.nuclear.func(w1, y, A, b, L, lambda, AtA, Atb)prox.nuclear.func(w1, y, A, b, L, lambda, AtA, Atb)
w1 |
previously updated model parameter |
y |
updated model parameter |
A |
design matrix |
b |
correspond vector, or matrix |
L |
learning rate |
lambda |
tuning parameter for low-rank component |
AtA |
Gram matrix of design matrix A |
Atb |
inner product of design matrix A and correspond vector b |
Value of proximal function with nuclear norm penalty
Proximal function with l1-norm
prox.sparse.func(w1, y, A, b, L, lambda, AtA, Atb)prox.sparse.func(w1, y, A, b, L, lambda, AtA, Atb)
w1 |
previously updated model parameter |
y |
updated model parameter |
A |
design matrix |
b |
correspond vector, or matrix |
L |
learning rate |
lambda |
tuning parameter for sparse component |
AtA |
Gram matrix of design matrix A |
Atb |
inner product of design matrix A and correspond vector b |
Value of proximal function with l1-norm penalty
Auxiliary function for FISTA implementation
Q.func(x, y, A, b, L, AtA, Atb)Q.func(x, y, A, b, L, AtA, Atb)
x |
Model parameter for previous update |
y |
Model parameter for updating |
A |
An n by p design matrix |
b |
A correspond vector, or matrix with size of n by 1 or n by p |
L |
Learning rate |
AtA |
Gram matrix for design matrix A |
Atb |
Inner product for design matrix A and correspond vector b |
Value of function Q
Shrinkage function for sparse soft-thresholding
shrinkage(y, tau)shrinkage(y, tau)
y |
A matrix, or a vector for thresholding |
tau |
A positive number, threshold |
A thresholded matrix, or vector
Shrinkage function for low-rank soft-thresholding
shrinkage.lr(y, tau)shrinkage.lr(y, tau)
y |
A matrix, or a vector for thresholding |
tau |
A positive number, threshold |
A thresholded matrix, or vector
L1-norm penalty for sparse component
sparse.pen(x, lambda)sparse.pen(x, lambda)
x |
Model parameter |
lambda |
Tuning parameter |
Value of l1-norm penalty term
summary function for S3 class for the fitting result
## S3 method for class 'LSVAR' summary(object, threshold = 0.2, ...)## S3 method for class 'LSVAR' summary(object, threshold = 0.2, ...)
object |
the S3 class object of |
threshold |
the threshold for sparse component visualization |
... |
not in use |
A series of summary for the S3 result
n <- 300 p <- 20 try <- testVAR(n, p, struct = "LS", signal = 0.75, rank = 2, singular_vals = c(1, 0.8)) data <- as.matrix(try$series) lambda <- 0.1; mu <- 1 fit <- fista.LpS(data, lambda = lambda, mu = mu, x.true = try$model_param) summary(fit, threshold = 0.2)n <- 300 p <- 20 try <- testVAR(n, p, struct = "LS", signal = 0.75, rank = 2, singular_vals = c(1, 0.8)) data <- as.matrix(try$series) lambda <- 0.1; mu <- 1 fit <- fista.LpS(data, lambda = lambda, mu = mu, x.true = try$model_param) summary(fit, threshold = 0.2)
A function to generate synthetic time series process based on the given structure
testVAR( n, p, struct = c("sparse", "low rank", "LS")[1], sp_density = 0.1, signal = NULL, rank = NULL, singular_vals, spectral_radius = 0.9, sigma = NULL, skip = 50, seed = 1 )testVAR( n, p, struct = c("sparse", "low rank", "LS")[1], sp_density = 0.1, signal = NULL, rank = NULL, singular_vals, spectral_radius = 0.9, sigma = NULL, skip = 50, seed = 1 )
n |
the length of time series |
p |
the number of multivariate time series |
struct |
a character string indicating the structure of the transition matrix, here are three options: sparse, low rank and LS (low rank plus sparse) |
sp_density |
a numeric value, indicating the sparsity density of sparse components, default is 0.1 |
signal |
a numeric value, indicating the magnitude of transition matrix |
rank |
a positive integer, the rank for low rank component |
singular_vals |
a numeric vector, indicating the singular values for the low rank component, the length of singular value must equal to the rank |
spectral_radius |
a numeric value, controlling the stability of the process, default is 0.9 |
sigma |
a numeric matrix, indicating the covariance matrix of noise term |
skip |
a numeric value, indicating the number of skipped time points in the beginning of the process |
seed |
an integer, indicating the seed for random seed. |
A list object, including
the generated time series
the noise term
true transition matrix
n <- 300; p <- 15 signal <- 0.75 rank <- 3 singular_vals <- c(1, 0.75, 0.5) try <- testVAR(n, p, struct = "LS", signal = signal, rank = rank, singular_vals = singular_vals) data <- as.matrix(try$series)n <- 300; p <- 15 signal <- 0.75 rank <- 3 singular_vals <- c(1, 0.75, 0.5) try <- testVAR(n, p, struct = "LS", signal = signal, rank = rank, singular_vals = singular_vals) data <- as.matrix(try$series)