holq {tensr}

Calculate the incredible higher-order 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

holq(X, tol = 10^-9, itermax = 1000, mode_rep = NULL, mode_diag = NULL,
  mode_ldu = NULL, print_diff = TRUE, start_vals = "identity",
  use_sig = TRUE)

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 use_sig).

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 sig from 1 (TRUE).

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 (2016)). 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 (2016).

Value

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 deviation" form.

Author(s)

David Gerard.

References

Gerard, D., & Hoff, P. (2016). A higher-order LQ decomposition for separable covariance models. Linear Algebra and its Applications, 505, 57-84. https://doi.org/10.1016/j.laa.2016.04.033 http://arxiv.org/pdf/1410.1094v1.pdf

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

#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)))