glamlassoRR {glamlasso} | R Documentation |
Penalized reduced rank regression in a GLAM
Description
Efficient design matrix free procedure for fitting large scale penalized reduced rank
regressions in a 3-dimensional generalized linear array model. To obtain a factorization of the parameter array,
the glamlassoRR
function performes a block relaxation scheme within the gdpg algorithm, see Lund and Hansen, 2018.
Usage
glamlassoRR(X,
Y,
Z = NULL,
family = "gaussian",
penalty = "lasso",
intercept = FALSE,
weights = NULL,
betainit = NULL,
alphainit = NULL,
nlambda = 100,
lambdaminratio = 1e-04,
lambda = NULL,
penaltyfactor = NULL,
penaltyfactoralpha = NULL,
reltolinner = 1e-07,
reltolouter = 1e-04,
reltolalt = 1e-04,
maxiter = 15000,
steps = 1,
maxiterinner = 3000,
maxiterouter = 25,
maxalt = 10,
btinnermax = 100,
btoutermax = 100,
iwls = "exact",
nu = 1)
Arguments
X |
A list containing the 3 tensor components of the tensor design matrix. These are matrices of sizes |
Y |
The response values, an array of size |
Z |
The non tensor structrured part of the design matrix. A matrix of size |
family |
A string specifying the model family (essentially the response distribution). Possible values
are |
penalty |
A string specifying the penalty. Possible values are |
intercept |
Logical variable indicating if the model includes an intercept. When |
weights |
Observation weights, an array of size |
betainit |
A list (length 2) containing the initial parameter values for each of the parameter factors. Default is NULL in which case all parameters are initialized at 0.01. |
alphainit |
A |
nlambda |
The number of |
lambdaminratio |
The smallest value for |
lambda |
The sequence of penalty parameters for the regularization path. |
penaltyfactor |
A list of length two containing an array of size |
penaltyfactoralpha |
A |
reltolinner |
The convergence tolerance for the inner loop |
reltolouter |
The convergence tolerance for the outer loop. |
reltolalt |
The convergence tolerance for the alternation loop over the two parameter blocks. |
maxiter |
The maximum number of inner iterations allowed for each |
steps |
The number of steps used in the multi-step adaptive lasso algorithm for non-convex penalties. Automatically set to 1 when |
maxiterinner |
The maximum number of inner iterations allowed for each outer iteration. |
maxiterouter |
The maximum number of outer iterations allowed for each lambda. |
maxalt |
The maximum number of alternations over parameter blocks. |
btinnermax |
Maximum number of backtracking steps allowed in each inner iteration. Default is |
btoutermax |
Maximum number of backtracking steps allowed in each outer iteration. Default is |
iwls |
A string indicating whether to use the exact iwls weight matrix or use a tensor structured approximation to it. |
nu |
A number between 0 and 1 that controls the step size |
Details
Given the setting from glamlasso
we place a reduced rank
restriction on the parameter array
given by
The glamlassoRR
function solves the PMLE problem by combining a
block relaxation scheme with the gdpg algorithm. This scheme alternates
between optimizing over the first parameter block
and the second block
while fixing the second resp.
first block.
Note that the individual parameter blocks are only identified up to a
multiplicative constant. Also note that the algorithm is sensitive to
inital values betainit
which can prevent convergence.
Value
An object with S3 Class "glamlasso".
spec |
A string indicating the model family and the penalty. |
coef12 |
A |
coef3 |
A |
alpha |
A |
lambda |
A vector containing the sequence of penalty values used in the estimation procedure. |
df |
The number of nonzero coefficients for each value of |
dimcoef |
A vector giving the dimension of the model coefficient array
|
dimobs |
A vector giving the dimension of the observation (response)
array |
Iter |
A list with 4 items: |
Author(s)
Adam Lund
Maintainer: Adam Lund, adam.lund@math.ku.dk
References
Lund, A., M. Vincent, and N. R. Hansen (2017). Penalized estimation in large-scale generalized linear array models. Journal of Computational and Graphical Statistics, 26, 3, 709-724. url = https://doi.org/10.1080/10618600.2017.1279548.
Lund, A. and N. R. Hansen (2019). Sparse Network Estimation for Dynamical Spatio-temporal Array Models. Journal of Multivariate Analysis, 174. url = https://doi.org/10.1016/j.jmva.2019.104532.
Examples
##size of example
n1 <- 65; n2 <- 26; n3 <- 13; p1 <- 12; p2 <- 6; p3 <- 4
##marginal design matrices (tensor components)
X1 <- matrix(rnorm(n1 * p1), n1, p1)
X2 <- matrix(rnorm(n2 * p2), n2, p2)
X3 <- matrix(rnorm(n3 * p3), n3, p3)
X <- list(X1, X2, X3)
Beta12 <- matrix(rnorm(p1 * p2), p1, p2) * matrix(rbinom(p1 * p2, 1, 0.5), p1, p2)
Beta3 <- matrix(rnorm(p3) * rbinom(p3, 1, 0.5), p3, 1)
Beta <- outer(Beta12, c(Beta3))
Mu <- RH(X3, RH(X2, RH(X1, Beta)))
Y <- array(rnorm(n1 * n2 * n3, Mu), dim = c(n1, n2, n3))
system.time(fit <- glamlassoRR(X, Y))
modelno <- length(fit$lambda)
oldmfrow <- par()$mfrow
par(mfrow = c(1, 3))
plot(c(Beta), type = "h")
points(c(Beta))
lines(c(outer(fit$coef12[, modelno], c(fit$coef3[, modelno]))), col = "red", type = "h")
plot(c(Beta12), ylim = range(Beta12, fit$coef12[, modelno]), type = "h")
points(c(Beta12))
lines(fit$coef12[, modelno], col = "red", type = "h")
plot(c(Beta3), ylim = range(Beta3, fit$coef3[, modelno]), type = "h")
points(c(Beta3))
lines(fit$coef3[, modelno], col = "red", type = "h")
par(mfrow = oldmfrow)
###with non tensor design component Z
q <- 5
alpha <- matrix(rnorm(q)) * rbinom(q, 1, 0.5)
Z <- matrix(rnorm(n1 * n2 * n3 * q), n1 * n2 * n3, q)
Y <- array(rnorm(n1 * n2 * n3, Mu + array(Z %*% alpha, c(n1, n2, n3))), c(n1, n2, n3))
system.time(fit <- glamlassoRR(X, Y, Z))
modelno <- length(fit$lambda)
oldmfrow <- par()$mfrow
par(mfrow = c(2, 2))
plot(c(Beta), type = "h")
points(c(Beta))
lines(c(outer(fit$coef12[, modelno], c(fit$coef3[, modelno]))), col = "red", type = "h")
plot(c(Beta12), ylim = range(Beta12,fit$coef12[, modelno]), type = "h")
points(c(Beta12))
lines(fit$coef12[, modelno], col = "red", type = "h")
plot(c(Beta3), ylim = range(Beta3, fit$coef3[, modelno]), type = "h")
points(c(Beta3))
lines(fit$coef3[, modelno], col = "red", type = "h")
plot(c(alpha), ylim = range(alpha, fit$alpha[, modelno]), type = "h")
points(c(alpha))
lines(fit$alpha[, modelno], col = "red", type = "h")
par(mfrow = oldmfrow)
################ poisson example
set.seed(7954) ## for this seed the algorithm fails to converge for default initial values!!
set.seed(42)
##size of example
n1 <- 65; n2 <- 26; n3 <- 13; p1 <- 12; p2 <- 6; p3 <- 4
##marginal design matrices (tensor components)
X1 <- matrix(rnorm(n1 * p1), n1, p1)
X2 <- matrix(rnorm(n2 * p2), n2, p2)
X3 <- matrix(rnorm(n3 * p3), n3, p3)
X <- list(X1, X2, X3)
Beta12 <- matrix(rnorm(p1 * p2, 0, 0.5) * rbinom(p1 * p2, 1, 0.1), p1, p2)
Beta3 <- matrix(rnorm(p3, 0, 0.5) * rbinom(p3, 1, 0.5), p3, 1)
Beta <- outer(Beta12, c(Beta3))
Mu <- RH(X3, RH(X2, RH(X1, Beta)))
Y <- array(rpois(n1 * n2 * n3, exp(Mu)), dim = c(n1, n2, n3))
system.time(fit <- glamlassoRR(X, Y ,family = "poisson"))
modelno <- length(fit$lambda)
oldmfrow <- par()$mfrow
par(mfrow = c(1, 3))
plot(c(Beta), type = "h")
points(c(Beta))
lines(c(outer(fit$coef12[, modelno], c(fit$coef3[, modelno]))), col = "red", type = "h")
plot(c(Beta12), ylim = range(Beta12, fit$coef12[, modelno]), type = "h")
points(c(Beta12))
lines(fit$coef12[, modelno], col = "red", type = "h")
plot(c(Beta3), ylim = range(Beta3, fit$coef3[, modelno]), type = "h")
points(c(Beta3))
lines(fit$coef3[, modelno], col = "red", type = "h")
par(mfrow = oldmfrow)