Fix single-column inputs in tpls() and facewise products

R/tens.mproduct.R

@@ -263,18 +263,28 @@ dctii_m_transforms <- function(t, bpparam = NULL) {
     # for-loop algorithm -------------------------------------------------------
     out <- array(0, dim = c(na, pb, t))
     for (i in seq_len(t)) {
-      out[, , i] <- a[, , i] %*% b[, , i]
+      # matrix() keeps each face 2D when a dimension is 1, where `[` would drop
+      # it to a vector and %*% would either fail or, for an n x 1 by 1 x n
+      # product, silently return the inner product instead of the outer one
+      out[, , i] <- matrix(a[, , i], nrow = na, ncol = pa) %*%
+        matrix(b[, , i], nrow = nb, ncol = pb)
     }
     # --------------------------------------------------------------------------
   } else {
     # BiocParallel algorithm ---------------------------------------------------
     # bltodo: benchmark / investigate preallocation here?
-    out <- simplify2array(
-      BiocParallel::bplapply(
+    # build the array from the known dimensions: simplify2array() flattens
+    # the result to a plain vector when every face is 1 x 1
+    out <- array(
+      unlist(BiocParallel::bplapply(
         array(seq_len(t)),
-        FUN = function(i) a[, , i] %*% b[, , i],
+        FUN = function(i) {
+          matrix(a[, , i], nrow = na, ncol = pa) %*%
+            matrix(b[, , i], nrow = nb, ncol = pb)
+        },
         BPPARAM = bpparam
-      )
+      )),
+      dim = c(na, pb, t)
     )
     # --------------------------------------------------------------------------
   }

R/tens.tpca.R

@@ -54,10 +54,16 @@
 #' @keywords internal
 #' @noRd
 .extract_tensor_columns <- function(tensor, k_t_indices) {
-  return(apply(
+  out <- apply(
     k_t_indices, 2,
     FUN = function(k_t_index) tensor[, k_t_index[1], k_t_index[2]]
-  ))
+  )
+  # apply() simplifies to a vector when the tensor has a single row, so put
+  # the matrix shape back
+  if (is.null(dim(out))) {
+    out <- matrix(out, nrow = dim(tensor)[1])
+  }
+  return(out)
 }
 
 #' Helper function to convert compressed matrix form of the singular values into

R/tens.tpls.R

@@ -29,7 +29,7 @@
   # bltodo: if ultimately only use once just place it in function body directly
   return(
     .unravel_index(
-      order(as.vector(s_mat))[1:ncomp],
+      order(as.vector(s_mat), decreasing = TRUE)[1:ncomp],
       dim(s_mat)
     )
   )
@@ -308,7 +308,10 @@ tpls <- function(
       # xi, omega (Wegelin); t, u (Kim-Anh); scores (scikit-learn)
       # note that only one face of xhat and yhat is relevant per iteration
       curr_x_projected <- xhat[, , k_t_top[2]] %*% curr_x_loadings
-      curr_y_projected <- yhat[, , k_t_top[2]] %*% curr_y_loadings
+      # matrix() keeps the face n x q when q = 1, where `[` would drop it to a
+      # vector and %*% would reject the vector-by-vector product
+      curr_y_projected <- matrix(yhat[, , k_t_top[2]], nrow = n) %*%
+        curr_y_loadings
 
       # gamma (Wegelin); c (Kim-Anh); loadings (scikit-learn)
       # the calculation of the x regression coefficient and deflation step of

tests/testthat/test-tens.mproduct.R

@@ -77,6 +77,39 @@ test_that(
   }
 )
 
+test_that(
+  "facewise product handles faces with a single row or column",
+  code = {
+    t <- 3
+    set.seed(1)
+    # `[` drops these faces to vectors: n x 1 by 1 x 1 used to error, n x 1 by
+    # 1 x n (an outer product) used to return the inner product, and 1 x p by
+    # p x 1 gives 1 x 1 faces, which the BiocParallel branch used to flatten
+    for (dims in list(
+      list(c(4, 1), c(1, 1)), list(c(4, 1), c(1, 4)), list(c(1, 4), c(4, 1))
+    )) {
+      test_tensor1 <- array(rnorm(prod(dims[[1]]) * t), dim = c(dims[[1]], t))
+      test_tensor2 <- array(rnorm(prod(dims[[2]]) * t), dim = c(dims[[2]], t))
+      expected_result <- array(0, dim = c(dims[[1]][1], dims[[2]][2], t))
+      for (i in 1:t) {
+        expected_result[, , i] <-
+          matrix(test_tensor1[, , i], dims[[1]][1], dims[[1]][2]) %*%
+          matrix(test_tensor2[, , i], dims[[2]][1], dims[[2]][2])
+      }
+      expect_equal(test_tensor1 %fp% test_tensor2, expected_result)
+      expect_equal(
+        suppressWarnings(
+          facewise_product(
+            test_tensor1, test_tensor2,
+            bpparam = BiocParallel::MulticoreParam()
+          )
+        ),
+        expected_result
+      )
+    }
+  }
+)
+
 test_that(
   "facewise product resolves names correctly, and propogates them",
   code = {

tests/testthat/test-tens.tpls.R

@@ -326,3 +326,60 @@ test_that(
     )
   }
 )
+
+test_that(
+  "tsvdm-tpls matrix output keeps the largest singular values",
+  code = {
+    n <- 6
+    p <- 5
+    q <- 4
+    t <- 3
+    ncomp <- 3
+
+    set.seed(1)
+    test_x <- array(rnorm(n * p * t, mean = 0, sd = 5), dim = c(n, p, t))
+    test_y <- array(rnorm(n * q * t, mean = 0, sd = 3), dim = c(n, q, t))
+
+    transforms <- dctii_m_transforms(t)
+    m <- transforms$m
+
+    # center = FALSE so the singular values can be recomputed directly
+    tpls_tsvdm <- tpls(
+      test_x, test_y,
+      ncomp = ncomp,
+      mode = "tsvdm",
+      center = FALSE
+    )
+
+    svals <- tsvdm(
+      ft(m(test_x)) %fp% m(test_y),
+      transform = FALSE,
+      svals_matrix_form = TRUE
+    )$s
+    expect_equal(
+      svals[cbind(tpls_tsvdm$features, tpls_tsvdm$faces)],
+      sort(svals, decreasing = TRUE)[1:ncomp]
+    )
+  }
+)
+
+test_that(
+  "tpls accepts a single-column response",
+  code = {
+    n <- 12
+    p <- 5
+    t <- 4
+    ncomp <- 2
+
+    set.seed(1)
+    test_x <- array(rnorm(n * p * t), dim = c(n, p, t))
+    test_y <- array(rnorm(n * t), dim = c(n, 1, t))
+
+    for (mode in c("regression", "canonical", "tsvdm")) {
+      tpls_obj <- tpls(test_x, test_y, ncomp = ncomp, mode = mode)
+      expect_equal(dim(tpls_obj$x_projected), c(n, ncomp))
+      expect_equal(dim(tpls_obj$y_loadings), c(1, ncomp))
+      expect_true(all(is.finite(tpls_obj$y_projected)))
+    }
+  }
+)