Calculate the incredible higherorder LQ decomposition (HOLQ).
Description
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.
Usage
1 2 3 
Arguments
X 
An array of numerics. 
tol 
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

itermax 
An integer. The maximum number of iterations of the LQ decomposition to do before stopping. 
mode_rep 
A vector of integers. The optional mode(s) to be considered identity matrices. 
mode_diag 
A vector of integers. The optional mode(s) to be considered as independent but heteroscedastic. 
mode_ldu 
A vector of integers. The optional modes(s) to be considered as having unit diagonal. 
print_diff 
A logical. Should the updates be printed to the screen each iteration? 
start_vals 
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. 
use_sig 
A logical. What stopping criterion should we use? Frobenius
norm of difference of cores (FALSE) or absolute value of difference of
ratio of 
Details
Given an array X
, the default version of this function will calculate
(1) 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
such that X
is the same as sig * atrans(Z,L)
up to numeric
precision.
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
juniors.
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
implemented in lrt_stat
.
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).
Value
Z
The core array with scaled allorthonormality property.
A
A list of matrices.
sig
A numeric. The total variation parameter. This is the "standard
deviation" form.
Author(s)
David Gerard.
References
Gerard, D. C., & Hoff, P. D. (2014). A higherorder LQ decomposition for separable covariance models. arXiv preprint arXiv:1410.1094.
See Also
array_bic_aic
for using the output of holq
to
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 holq
to
calculate the maximum likelihood estimates of the component covariance
matrices under the array normal model.
Examples
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)))
