diff --git a/tensorPredictors/NAMESPACE b/tensorPredictors/NAMESPACE
index 221aacf..b623699 100644
--- a/tensorPredictors/NAMESPACE
+++ b/tensorPredictors/NAMESPACE
@@ -24,6 +24,7 @@ export(matrixImage)
 export(mcrossprod)
 export(reduce)
 export(rowKronecker)
+export(rtensornorm)
 export(tensor_predictor)
 export(ttm)
 import(stats)
diff --git a/tensorPredictors/R/dist_subspace.R b/tensorPredictors/R/dist_subspace.R
index 19d23f8..fdcfc81 100644
--- a/tensorPredictors/R/dist_subspace.R
+++ b/tensorPredictors/R/dist_subspace.R
@@ -40,7 +40,7 @@ dist.subspace <- function (A, B, is.ortho = FALSE, normalize = FALSE,
 
     if (normalize) {
         rankSum <- ncol(A) + ncol(B)
-        c <- 1 / sqrt(min(rankSum, 2 * nrow(A) - rankSum))
+        c <- 1 / sqrt(max(1, min(rankSum, 2 * nrow(A) - rankSum)))
     } else {
         c <- sqrt(2)
     }
diff --git a/tensorPredictors/R/kpir_ls.R b/tensorPredictors/R/kpir_ls.R
new file mode 100644
index 0000000..871e717
--- /dev/null
+++ b/tensorPredictors/R/kpir_ls.R
@@ -0,0 +1,53 @@
+#' 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
+    }
+
+}
+
+
diff --git a/tensorPredictors/R/mlm.R b/tensorPredictors/R/mlm.R
new file mode 100644
index 0000000..c27ef59
--- /dev/null
+++ b/tensorPredictors/R/mlm.R
@@ -0,0 +1,71 @@
+#' Multi Linear Multiplication
+#'
+#'     C = A x { B1, ..., Br }
+#'
+#' @param A tensor (multi-linear array)
+#' @param B matrix or list of matrices
+#' @param ... further matrices, concatenated with \code{B}
+#' @param modes integer sequence of the same length as number of matrices
+#'  supplied (in \code{B} and \code{...})
+#'
+#' @examples
+#' # general usage
+#' dimA <- c(3, 17, 19, 2)
+#' dimC <- c(7, 11, 13, 5)
+#' A <- array(rnorm(prod(dimA)), dim = dimA)
+#' B <- Map(function(p, q) matrix(rnorm(p * q), p, q), dimC, dimA)
+#' C1 <- mlm(A, B)
+#' C2 <- mlm(A, B[[1]], B[[2]], B[[3]], B[[4]])
+#' C3 <- mlm(A, B[[3]], B[[1]], B[[2]], B[[4]], modes = c(3, 1, 2, 4))
+#' C4 <- mlm(A, B[1:3], B[[4]])
+#' stopifnot(all.equal(C1, C2))
+#' stopifnot(all.equal(C1, C3))
+#' stopifnot(all.equal(C1, C4))
+#'
+#' # selected modes
+#' C1 <- mlm(A, B, modes = 2:3)
+#' C2 <- mlm(A, B[[2]], B[[3]], modes = 2:3)
+#' C3 <- ttm(ttm(A, B[[2]], 2), B[[3]], 3)
+#' stopifnot(all.equal(C1, C2))
+#' stopifnot(all.equal(C1, C3))
+#'
+#' # analog to matrix multiplication
+#' A <- matrix(rnorm( 6), 2, 3)
+#' B <- matrix(rnorm(12), 3, 4)
+#' C <- matrix(rnorm(20), 4, 5)
+#' stopifnot(all.equal(
+#'     A %*% B %*% t(C),
+#'     mlm(B, list(A, C))
+#' ))
+#'
+#' # usage with repeated modes (non commutative)
+#' dimA <- c(3, 17, 19, 2)
+#' A <- array(rnorm(prod(dimA)), dim = dimA)
+#' B1 <- matrix(rnorm(9), 3, 3)
+#' B2 <- matrix(rnorm(9), 3, 3)
+#' C <- matrix(rnorm(4), 2, 2)
+#' # same modes do NOT commute
+#' all.equal(
+#'     mlm(A, B1, B2, C, modes = c(1, 1, 4)),  # NOT equal!
+#'     mlm(A, B2, B1, C, modes = c(1, 1, 4))
+#' )
+#' # but different modes do commute
+#' P1 <- mlm(A, C, B1, B2, modes = c(4, 1, 1))
+#' P2 <- mlm(A, B1, C, B2, modes = c(1, 4, 1))
+#' P3 <- mlm(A, B1, B2, C, modes = c(1, 1, 4))
+#' stopifnot(all.equal(P1, P2))
+#' stopifnot(all.equal(P1, P3))
+#'
+#' @export
+mlm <- function(A, B, ..., modes = seq_along(B)) {
+    # Collect all matrices in `B`
+    B <- c(if (is.matrix(B)) list(B) else B, list(...))
+
+    # iteratively apply Tensor Times Matrix multiplication over modes
+    for (i in seq_along(modes)) {
+        A <- ttm(A, B[[i]], modes[i])
+    }
+
+    # return result tensor
+    A
+}