#' A multi-linear model
#'
#' @export
gmlm_mlm <- function(X, F, sample.axis) {
    # get (set) problem (observation) dimensions
    modes <- seq_along(dim(X))[-sample.axis]
    dimX <- dim(X)[modes]
    sample.size <- dim(X)[sample.axis]
    if (!is.array(F)) {
        dim(F) <- c(1L, sample.size)[(seq_along(dim(X)) == sample.axis) + 1L]
    }
    dimF <- dim(F)[modes]

    # vectorize the tensor valued data (columns are the vectorized samples)
    matX <- mat(X, modes)
    matF <- mat(F, modes)

    # center
    matX <- matX - (meanX <- rowMeans(matX))
    matF <- matF - rowMeans(matF)

    # solve vectorized linear model
    B <- tcrossprod(matX, matF) %*% pinv(tcrossprod(matF))

    # decompose linear model solution as Kronecker product
    betas <- approx.kron(B, dimX, dimF)

    # reshape centered vectorized `X` and `F` into tensors (now, sample axis
    # is last axis)
    X <- `dim<-`(matX, c(dimX, sample.size))
    F <- `dim<-`(matF, c(dimF, sample.size))

    # and estimate covariances (sample axis is last axis)
    Sigmas <- mcov(X - mlm(F, betas), sample.axis = length(dim(X)))

    # finaly, invert covariances to get the scatter matrices `Omegas`
    Omegas <- Map(solve, Sigmas)

    list(
        eta1 = array(meanX, dim = dim(X)[-sample.axis]),
        betas = betas,
        Omegas = Omegas
    )
}