ihop {tensr}

The incredible higher-order polar decomposition (IHOP).

Description

Mmm, pancakes.

Usage

ihop(X, itermax = 100, tol = 10^-9, print_diff = TRUE, mode_rep = NULL,
  use_sig = TRUE)

Arguments

X

An array of numerics.

itermax

An integer. The maximum number of iterations to perform during the optimization procedure.

tol

A numeric. The algorithm will stop when the Frobenius norm of the difference of core arrays between subsequent iterations is below tol (for use_sig = FALSE) or when the absolute difference between the ratio of subsequent values of sig and 1 is less than tol (for use_sig = TRUE).

print_diff

A logical. Should we print the updates of the algorithm?

mode_rep

A vector. Which component matrices should be set to be the identity?

use_sig

A logical. See tol.

Details

This function will calculate the higher-order polar decomposition, a generalization of the polar decomposition to tensors. It generalizes a minimization formulation of the polar decomposition.

Given an array X, ihop will output L a list of lower triangular matrices with positive diagonal elements and unit Frobenius norm, R a core array with certain orthogonality properties, and sig a total variation parameter. We have that X is equal to sig * atrans(R, L) up to numerical precision.

t(solve(L[[i]])) %*% mat(R, i) will have orthonormal rows for all i.

For more details on the IHOP, see Gerard and Hoff (2016).

Value

R A core array which, in combination with L, has certain orthogonality properties.

L A list of lower triangular matrices with unit Frobenius norm.

sig A numeric.

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

Examples

#Generate random data.
p <- c(2, 3, 4)
X <- array(stats::rnorm(prod(p)), dim = p)

#Calculate IHOP.
ihop_x <- ihop(X)
R <- ihop_x$R
L <- ihop_x$L
sig <- ihop_x$sig

#Reconstruct X
trim(X - sig * atrans(R, L))

#Orthogonality properties
ortho_1 <- t(solve(L[[1]])) %*% mat(R, 1)
trim(ortho_1 %*% t(ortho_1))

ortho_2 <- t(solve(L[[2]])) %*% mat(R, 2)
trim(ortho_2 %*% t(ortho_2))

ortho_3 <- t(solve(L[[3]])) %*% mat(R, 3)
trim(ortho_3 %*% t(ortho_3))