I've published some R code on CRAN that can be used to calculate the incredible higher-order LQ decomposition (incredible HOLQ), the incredible singular value decomposition (ISVD), and the incredible higher-order polar decomposition (IHOP). The HOLQ can be used to (1) derive the maximum likelihood estimates of the covariance matrices under the array normal model, (2) run a likelihood ratio test in separable covariance models, and (3) calculate AIC's and BIC's for separable covariance models. There are also a few useful array and matrix functions. A description of the methods can be found in:

Gerard, D. , & Hoff, P. (2016). A higher-order LQ decomposition for separable covariance models. Linear Algebra and its Applications, 505, 57-84. [Link to LAA] [Link to arXiv] [bib]

I also provide a brief vignette on how to use the functions available in the R code. To download, simply type in R:

install.packages("tensr")

Alternatively, you can download it directly from CRAN.

Please provide us with questions or comments: David Gerard (dcgerard) or Peter Hoff (pdhoff) -- @uchicago.edu and @uw.edu, respectively.