Calculate the incredible higher-order LQ decomposition (HOLQ).
This function calculates a generalization of the LQ decomposition to tensors. This decomposition has a close connection to likelihood inference in Kronecker structured covariande models.
1 2 3
An array of numerics.
A numeric. The maximum difference in frobenius norm between two
successive core arrays before stopping. Or maximum difference of the ratio
of sigs from 1 before stopping (which depends on the value of
An integer. The maximum number of iterations of the LQ decomposition to do before stopping.
A vector of integers. The optional mode(s) to be considered identity matrices.
A vector of integers. The optional mode(s) to be considered as independent but heteroscedastic.
A vector of integers. The optional modes(s) to be considered as having unit diagonal.
A logical. Should the updates be printed to the screen each iteration?
Determines how to start the optimization algorithm. Either 'identity' (default), or 'residuals', which results in using the cholesky square roots of the sample covariance matrices along each mode scaled to have unit determinant. You can also use your own start values.
A logical. What stopping criterion should we use? Frobenius
norm of difference of cores (FALSE) or absolute value of difference of
Given an array
X, the default version of this function will calculate
L a list of lower triangular matricies with positive diagonal
elements and unit determinant,
Z an array of the same dimensions as
X that has special orthogonal properties, and (3)
sig a numeric
X is the same as
sig * atrans(Z,L) up to numeric
This output (1) can be considered a generalization of the LQ decomposition to tensors, (2) solves an optimization problem which the matrix LQ decomposition solves, and (3) has a special connection to likelihood inference in the array normal model.
There are options to constrain the matrices in
L to either be
diagonal, lower triangular with unit diagonal, or the identity matrix. Each
of these correspond to submodels in Kronecker structured covariance models.
The core array corresponding to each of these options has different
properities (see Gerard and Hoff
(2014)). These more constrained tensor decompositions are called HOLQ
The MLE of the ith component covariance matrix under any
elliptically contoured Kronecker structured covariance model is given by
L[[i]] %*% t(L[[i]]). The MLE for the total variation pamarameter
will be different depending on the distribution of the array, but for the
array normal it is
sig ^ 2 / prod(p) (the "variance" form for the
total variation parameter).
The likelihood ratio test statistic depends only on
sig and can be
The algorithm used to fit the HOLQ iteratively repeats the LQ decomposition along each mode.
For more details on the incredible HOLQ, see Gerard and Hoff (2014).
Z The core array with scaled all-orthonormality property.
A A list of matrices.
sig A numeric. The total variation parameter. This is the "standard
Gerard, D. C., & Hoff, P. D. (2014). A higher-order LQ decomposition for separable covariance models. arXiv preprint arXiv:1410.1094.
array_bic_aic for using the output of
calculate AIC and BIC,
get_isvd for using the output of
holq to calculate a
tensor generalization of the singular value decomposition.
lq for the matrix LQ decomposition.
lrt_stat for using the output of
holq to calculate
likelihood ratio test statistics.
mle_from_holq for using the output of
calculate the maximum likelihood estimates of the component covariance
matrices under the array normal model.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
#Genrate random data. p <- c(2, 3, 4) X <- array(stats::rnorm(prod(p)), dim = p) #Calculate HOLQ with unit diagonal on 2nd mode, # and diagonal along 3rd mode. holq_x <- holq(X, mode_ldu = 2, mode_diag = 3) Z <- holq_x$Z A <- holq_x$A sig <- holq_x$sig #Reconstruct X trim(X - sig * atrans(Z, A), 10^-5) #Properties of core #First mode has orthonormal rows. trim(mat(Z, 1) %*% t(mat(Z, 1)), 10^-5) #Second mode has orthogonal rows. trim(mat(Z, 2) %*% t(mat(Z, 2)), 10^-5) #Third mode has unit diagonal (up to scale). diag(mat(Z, 3) %*% t(mat(Z, 3)))