#' Per mode (axis) alternating least squares estimate
#'
#' @param sample.mode index of the sample mode, a.k.a. observation axis index
#'
#' @export
kpir.ls <- function(X, Fy, max.iter = 20L, sample.mode = 1L,
    eps = .Machine$double.eps, logger = NULL
) {
    # Check if X and Fy have same number of observations
    if (!is.array(Fy)) {
        # scalar response case (add new axis of size 1)
        dim(Fy) <- local({
            dims <- rep(1, length(dim(X)))
            dims[sample.mode] <- length(Fy)
            dims
        })
    } else {
        stopifnot(dim(X)[sample.mode] == dim(Fy)[sample.mode])
    }
    # and check shape
    stopifnot(length(X) == length(Fy))

    # mode index sequence (exclude sample mode, a.k.a. observation axis)
    modes <- seq_along(dim(X))[-sample.mode]

    ### Step 1: initial per mode estimates
    alphas <- Map(function(mode, ncol) {
        La.svd(mcrossprod(X, mode), ncol)$u
    }, modes, dim(Fy)[modes])

    # Call history callback (logger) before the first iteration
    if (is.function(logger)) { logger(0L, alphas) }


    ### Step 2: iterate per mode (axis) least squares estimates
    for (iter in seq_len(max.iter)) {
        # cyclic iterate over modes
        for (j in seq_along(modes)) {
            # least squares solution for `alpha_j | alpha_i, i != j`
            Z <- mlm(Fy, alphas[-j], modes = modes[-j])
            alphas[[j]] <- t(solve(mcrossprod(Z, j), tcrossprod(mat(Z, j), mat(X, j))))
            # TODO: alphas[[j]] <- t(solve(mcrossprod(Z, j), mcrossprod(Z, X, j)))
        }

        # Call logger (invoke history callback)
        if (is.function(logger)) { logger(iter, alphas) }

        # TODO: add some kind of break condition
    }

}