diff --git a/NAMESPACE b/NAMESPACE
index cd214ab6..c3316f5f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -48,6 +48,7 @@ export(rgcca_stability)
 export(rgcca_transform)
 importFrom(Deriv,Deriv)
 importFrom(MASS,ginv)
+importFrom(abind,abind)
 importFrom(caret,confusionMatrix)
 importFrom(caret,multiClassSummary)
 importFrom(caret,postResample)
@@ -66,7 +67,6 @@ importFrom(matrixStats,rowSds)
 importFrom(matrixStats,rowSums2)
 importFrom(methods,is)
 importFrom(rlang,.data)
-importFrom(stats,complete.cases)
 importFrom(stats,cor)
 importFrom(stats,median)
 importFrom(stats,model.matrix)
diff --git a/R/block.R b/R/block.R
index d2093416..03add9f1 100644
--- a/R/block.R
+++ b/R/block.R
@@ -45,19 +45,30 @@ new_sparse_block <- function(x, j, sparsity, tol = 1e-08, ...) {
   )
 }
 
-new_tensor_block <- function(x, j, ..., class = character()) {
+new_tensor_block <- function(x, j, rank, mode_orth, ..., class = character()) {
   new_block(
-    x, j, factors = NULL, weights = NULL, ..., class = c(class, "tensor_block")
+    x, j, rank = rank, mode_orth = mode_orth, factors = NULL,
+    weights = NULL, ..., class = c(class, "tensor_block")
+  )
+}
+
+new_regularized_tensor_block <- function(x, j, rank, mode_orth, tau, ...) {
+  new_tensor_block(
+    x, j, rank = rank, mode_orth = mode_orth, tau = tau,
+    M = NULL, ..., class = "tensor_regularized_block"
   )
 }
 
 ### Utility method to choose the adequate class
-create_block <- function(x, j, bias, na.rm, tau, sparsity, tol) {
+create_block <- function(x, j, bias, na.rm, tau, sparsity,
+                         tol, rank, mode_orth) {
   if (length(dim(x)) > 2) {         # TGCCA
     if (tau < 1) {
-      res <- new_regularized_tensor_block(x, j, tau, bias = bias, na.rm = na.rm)
+      res <- new_regularized_tensor_block(
+        x, j, rank, mode_orth, tau, bias = bias, na.rm = na.rm
+      )
     } else {
-      res <- new_tensor_block(x, j, bias = bias, na.rm = na.rm)
+      res <- new_tensor_block(x, j, rank, mode_orth, bias = bias, na.rm = na.rm)
     }
   } else {
     if (sparsity < 1) {             # SGCCA
diff --git a/R/block_init.R b/R/block_init.R
index 0ca52048..ba05ed71 100644
--- a/R/block_init.R
+++ b/R/block_init.R
@@ -44,14 +44,20 @@ block_init.dual_regularized_block <- function(x, init = "svd") {
 block_init.tensor_block <- function(x, init = "svd") {
   if (init == "svd") {
     x$factors <- lapply(seq_along(dim(x$x))[-1], function(m) {
-      svd(apply(x$x, m, c), nu = 0, nv = 1)$v
+      svd(apply(x$x, m, c), nu = 0, nv = x$rank)$v
     })
   } else {
     x$factors <- lapply(seq_along(dim(x$x))[-1], function(m) {
-      svd(matrix(rnorm(dim(x$x)[m]), dim(x$x)[m]), nu = 1, nv = 0)$u
+      if (m == x$mode_orth) {
+        svd(matrix(
+          rnorm(dim(x$x)[m] * x$rank), dim(x$x)[m]
+        ), nu = x$rank, nv = 0)$u
+      } else {
+        matrix(rnorm(dim(x$x)[m] * x$rank), dim(x$x)[m])
+      }
     })
   }
-  x$weights <- rep(1 / sqrt(1), 1)
+  x$weights <- rep(1 / sqrt(x$rank), x$rank)
 
   return(block_project(x))
 }
diff --git a/R/block_update.R b/R/block_update.R
index c456f484..90a3b4b1 100644
--- a/R/block_update.R
+++ b/R/block_update.R
@@ -29,8 +29,17 @@ block_update.tensor_block <- function(x, grad) {
     grad_m <- grad_m %*% khatri_rao(
       Reduce(khatri_rao, rev(x$factors[-seq_len(m)])), other_factors
     )
-    x$factors[[m]] <- grad_m %*% diag(x$weights, nrow = length(x$weights))
-    x$factors[[m]] <- x$factors[[m]] / norm(x$factors[[m]], type = "2")
+    if (m == x$mode_orth) {
+      SVD <- svd(
+        grad_m %*% diag(x$weights, nrow = x$rank), nu = x$rank, nv = x$rank
+      )
+      x$factors[[m]] <- SVD$u %*% t(SVD$v)
+    } else {
+      x$factors[[m]] <- grad_m %*% diag(x$weights, nrow = x$rank)
+      x$factors[[m]] <- apply(
+        x$factors[[m]], 2, function(y) y / norm(y, type = "2")
+      )
+    }
 
     other_factors <- khatri_rao(x$factors[[m]], other_factors)
   }
diff --git a/R/check_blocks.r b/R/check_blocks.r
index f1eb9f0e..0cd293bd 100644
--- a/R/check_blocks.r
+++ b/R/check_blocks.r
@@ -33,6 +33,7 @@
 #'   dimnames in the same order.
 #' }
 #' @inheritParams rgcca
+#' @importFrom abind abind
 #' @noRd
 check_blocks <- function(blocks, quiet = FALSE, response = NULL) {
   blocks <- check_blocks_is_list(blocks)
@@ -258,7 +259,7 @@ check_blocks_align <- function(blocks, m = 1) {
     y <- array(NA, dim = extra_dim)
     dimnames(y)[[m]] <- missing_names
     rownames(y) <- missing_names
-    return(abind::abind(x, y, along = m))
+    return(abind(x, y, along = m))
   })
 
   # Align blocks using names on dimension m
diff --git a/R/checks.R b/R/checks.R
index 1f7ff514..50b0c5cc 100644
--- a/R/checks.R
+++ b/R/checks.R
@@ -219,8 +219,12 @@ check_ncomp <- function(ncomp, blocks, min = 1, superblock = FALSE,
         "only one number of components must be specified (superblock)."
       )
     }
+    is_superblock_present <- "superblock" %in% names(blocks)
+    ncomp <- rep(max(ncomp), ifelse(
+      is_superblock_present, length(blocks), length(blocks) + 1
+    ))
     max_ncomp <- ifelse(
-      "superblock" %in% names(blocks),
+      is_superblock_present,
       NCOL(blocks[[length(blocks)]]),
       sum(vapply(blocks, NCOL, FUN.VALUE = integer(1)))
     )
@@ -240,6 +244,11 @@ check_ncomp <- function(ncomp, blocks, min = 1, superblock = FALSE,
 
   ncomp <- elongate_arg(ncomp, blocks)
   check_size_blocks(blocks, "ncomp", ncomp)
+
+  if (!is.null(response)) {
+    ncomp[response] <- max(ncomp[-response])
+  }
+
   ncomp <- vapply(
     seq_along(ncomp),
     function(x) {
@@ -319,7 +328,7 @@ check_size_blocks <- function(blocks, x, y = x, n_row = NULL,
 }
 
 check_penalty <- function(penalty, blocks, method = "rgcca", superblock = FALSE,
-                          ncomp = NULL) {
+                          ncomp = 1) {
   penalty <- elongate_arg(penalty, blocks)
   is_matrix <- is.matrix(penalty)
   DIM <- dim(penalty)
@@ -348,7 +357,10 @@ check_penalty <- function(penalty, blocks, method = "rgcca", superblock = FALSE,
     )
   }
 
-  if (is_matrix) penalty <- matrix(penalty, DIM[1], DIM[2])
+  penalty <- matrix(penalty, ncol = size)
+  if (nrow(penalty) < ncomp) {
+    penalty <- matrix(penalty, nrow = ncomp, ncol = size, byrow = TRUE)
+  }
 
   return(penalty)
 }
@@ -378,6 +390,49 @@ check_tau <- function(tau) {
   invisible(tau)
 }
 
+check_mode_orth <- function(mode_orth, blocks) {
+  mode_orth <- elongate_arg(mode_orth, blocks)
+  mode_orth <- check_integer("mode_orth", mode_orth, min = 1, type = "vector")
+  check_size_blocks(blocks, "mode_orth", mode_orth)
+  lapply(seq_along(blocks), function(j) {
+    if (length(dim(blocks[[j]])) > 2 &&
+        mode_orth[j] > length(dim(blocks[[j]])) - 1) {
+      stop_rgcca(
+        "mode_orth should be comprise between ", 1 ,
+        " and ", length(dim(blocks[[j]])) - 1, " (that is the number of
+          modes of block ", j, " minus 1)."
+      )
+    }
+  })
+  return(mode_orth)
+}
+
+check_rank <- function(rank, blocks, mode_orth, ncomp) {
+  rank <- elongate_arg(rank, blocks)
+  p <- ifelse(is.vector(rank), 1, nrow(rank))
+  J <- length(blocks)
+  rank <- vapply(seq_along(rank), function(x) {
+    y <- check_integer("rank", rank[x], min = 1)
+    j <- (x - 1) %% J + 1
+    if (length(dim(blocks[[j]])) > 2 &&
+        y > dim(blocks[[j]])[mode_orth[j] + 1]) {
+      stop_rgcca(
+        "rank[", x, "] should be comprise between ", 1 ,
+        " and ", dim(blocks[[j]])[mode_orth[j] + 1], " (that is the number of",
+        " variables of the mode bearing the orthogonality for block ", j, ")."
+      )
+    }
+    return(y)
+  }, FUN.VALUE = 1L)
+
+  rank <- matrix(rank, nrow = p)
+  check_size_blocks(blocks, "rank", rank, n_row = ncomp)
+  if (p < ncomp) {
+    rank <- matrix(rank, nrow = ncomp, ncol = J, byrow = TRUE)
+  }
+  return(rank)
+}
+
 check_scheme <- function(scheme) {
   if (mode(scheme) != "function") {
     scheme <- tolower(scheme)
diff --git a/R/format_output.R b/R/format_output.R
index e3bcfb4c..a7850d1a 100644
--- a/R/format_output.R
+++ b/R/format_output.R
@@ -43,11 +43,11 @@ format_output <- function(func_out, rgcca_args, opt, blocks) {
   names(func_out$AVE$AVE_X_cor) <- names_AVE
 
   ### Set names and shave
-  array_idx <- which(
-    vapply(blocks, function(x) length(dim(x)) > 2, FUN.VALUE = logical(1L))
+  array_idx <- vapply(
+    blocks, function(x) length(dim(x)) > 2, FUN.VALUE = logical(1L)
   )
 
-  for (j in array_idx) {
+  for (j in seq_along(blocks)[array_idx]) {
     func_out$factors[[j]] <- lapply(
       seq_along(func_out$factors[[j]]),
       function(m) {
@@ -60,7 +60,7 @@ format_output <- function(func_out, rgcca_args, opt, blocks) {
     rownames(func_out$a[[j]]) <- do.call(paste, c(grid, sep = "_"))
   }
 
-  for (j in seq_along(blocks)[-array_idx]) {
+  for (j in seq_along(blocks)[!array_idx]) {
     rownames(func_out$a[[j]]) <- colnames(blocks[[j]])
   }
 
diff --git a/R/get_rgcca_args.R b/R/get_rgcca_args.R
index 12ce4b41..1b9f61ad 100644
--- a/R/get_rgcca_args.R
+++ b/R/get_rgcca_args.R
@@ -23,6 +23,7 @@ get_rgcca_args <- function(object, default_args = list()) {
       tol = default_args$tol,
       init = tolower(default_args$init),
       bias = default_args$bias,
+      rank = default_args$rank,
       quiet = default_args$quiet,
       scale = default_args$scale,
       ncomp = default_args$ncomp,
@@ -34,6 +35,7 @@ get_rgcca_args <- function(object, default_args = list()) {
       response = default_args$response,
       NA_method = tolower(default_args$NA_method),
       comp_orth = default_args$comp_orth,
+      mode_orth = default_args$mode_orth,
       n_iter_max = default_args$n_iter_max,
       connection = default_args$connection,
       superblock = default_args$superblock,
diff --git a/R/plot.rgcca_cv.R b/R/plot.rgcca_cv.R
index 55abee8a..034b2894 100644
--- a/R/plot.rgcca_cv.R
+++ b/R/plot.rgcca_cv.R
@@ -30,7 +30,7 @@ plot.rgcca_cv <- function(x, type = c("sd", "quantile"),
   }
 
   best <- which(apply(
-    x$params, 1, function(z) identical(z, x$best_params)
+    x$params, 1, function(z) identical(z, drop(x$best_params))
   ))
 
   idx_order <- seq_len(nrow(df))
diff --git a/R/plot.rgcca_permutation.R b/R/plot.rgcca_permutation.R
index 0bfaabdd..52ee1c39 100644
--- a/R/plot.rgcca_permutation.R
+++ b/R/plot.rgcca_permutation.R
@@ -32,7 +32,7 @@ plot.rgcca_permutation <- function(x,
 
   # Mark the best parameter set
   best <- which(apply(
-    x$params[idx_order, ], 1, function(z) identical(z, x$best_params)
+    x$params[idx_order, ], 1, function(z) identical(z, drop(x$best_params))
   ))
   df$label[best] <- "Best parameter set"
 
diff --git a/R/rgcca.R b/R/rgcca.R
index ece5e34f..ce407444 100644
--- a/R/rgcca.R
+++ b/R/rgcca.R
@@ -174,6 +174,18 @@
 #' orthogonal block components or orthogonal block weight vectors.
 #' @param A Deprecated argument, please use blocks instead.
 #' @param C Deprecated argument, please use connection instead.
+#' @param rank Either an integer, an integer vector of
+#' size \eqn{J} or an integer matrix
+#' of dimension \eqn{\textrm{max}(\textrm{ncomp}) \times J} giving the rank
+#' of the decomposition sought for the canonical vectors in TGCCA.
+#' If block \eqn{j} is an array with at least three dimensions, rank must be
+#' comprised between 1 and the number of variables on the mode bearing the
+#' orthogonality constraint. See \textrm{mode_orth}.
+#' @param mode_orth Either an integer or an integer vector of size \eqn{J}
+#' designating the mode which associated set of factors will be orthogonal
+#' in the decomposition sought by TGCCA. If block \eqn{j} is an array with
+#' \eqn{d > 2} dimensions, \textrm{mode_orth} must be comprised between
+#' 1 and \eqn{d - 1}.
 #' @return A fitted rgcca object.
 #' @return \item{Y}{A list of \eqn{J} elements. The jth element
 #' of the list \eqn{Y}
@@ -433,6 +445,7 @@ rgcca <- function(blocks, connection = NULL, tau = 1, ncomp = 1,
                   superblock = FALSE,
                   NA_method = "na.ignore", quiet = TRUE,
                   n_iter_max = 1000, comp_orth = TRUE,
+                  rank = 1, mode_orth = 1,
                   A = NULL, C = NULL) {
   # Check for deprecated arguments
   if (!missing(A)) {
@@ -478,14 +491,12 @@ rgcca <- function(blocks, connection = NULL, tau = 1, ncomp = 1,
   }
 
   ### Call the gcca function
-  gcca_args <- rgcca_args[c(
-    "connection", "ncomp", "scheme", "init", "bias", "tol",
-    "verbose", "superblock", "response", "n_iter_max", "comp_orth"
-  )]
+  gcca_args <- rgcca_args[-which(names(rgcca_args) %in% c(
+    "quiet", "scale", "scale_block", "method", "NA_method"
+  ))]
   gcca_args[["na.rm"]] <- na.rm
   gcca_args[["blocks"]] <- blocks
   gcca_args[["disjunction"]] <- opt$disjunction
-  gcca_args[[opt$param]] <- rgcca_args[[opt$param]]
   func_out <- do.call(rgcca_outer_loop, gcca_args)
 
   ### Format the output
diff --git a/R/rgcca_bootstrap.R b/R/rgcca_bootstrap.R
index 3ba651e7..971db03e 100644
--- a/R/rgcca_bootstrap.R
+++ b/R/rgcca_bootstrap.R
@@ -105,8 +105,7 @@ rgcca_bootstrap <- function(rgcca_res, n_boot = 100,
     rgcca_res$call$blocks <- Map(
       function(x, y) x[, y, drop = FALSE], rgcca_res$call$blocks, keep_var
     )
-    rgcca_res$call$tau <-
-      rgcca_res$call$sparsity <- rep(1, length(rgcca_res$blocks))
+    rgcca_res$call$tau[] <- rgcca_res$call$sparsity[] <- 1
 
     rgcca_res <- rgcca(rgcca_res)
   }
diff --git a/R/rgcca_cv.r b/R/rgcca_cv.r
index 05e1607e..04134844 100644
--- a/R/rgcca_cv.r
+++ b/R/rgcca_cv.r
@@ -186,6 +186,8 @@ rgcca_cv <- function(blocks,
                      verbose = TRUE,
                      n_iter_max = 1000,
                      comp_orth = TRUE,
+                     rank = 1,
+                     mode_orth = 1,
                      ...) {
   ### Try to retrieve parameters from a rgcca object
   rgcca_args <- as.list(environment())
@@ -338,7 +340,7 @@ rgcca_cv <- function(blocks,
     par_type = param$par_type,
     params = param$par_value,
     validation = validation,
-    best_params = param$par_value[best_param_idx, ],
+    best_params = param$par_value[best_param_idx, , drop = FALSE],
     classification = model$classification,
     prediction_model = model$model_name
   )
diff --git a/R/rgcca_inner_loop.R b/R/rgcca_inner_loop.R
index 49924ccd..9d8af5e5 100644
--- a/R/rgcca_inner_loop.R
+++ b/R/rgcca_inner_loop.R
@@ -1,7 +1,9 @@
 rgcca_inner_loop <- function(A, C, g, dg, tau = rep(1, length(A)),
                              sparsity = rep(1, length(A)),
                              verbose = FALSE, init = "svd", bias = TRUE,
-                             tol = 1e-08, na.rm = TRUE, n_iter_max = 1000) {
+                             tol = 1e-08, na.rm = TRUE, n_iter_max = 1000,
+                             rank = rep(1, length(A)),
+                             mode_orth = rep(1, length(A))) {
   if (!is.numeric(tau)) {
     # From Schafer and Strimmer, 2005
     tau <- vapply(A, tau.estimate, na.rm = na.rm, FUN.VALUE = 1.0)
@@ -15,7 +17,9 @@ rgcca_inner_loop <- function(A, C, g, dg, tau = rep(1, length(A)),
 
   ### Initialization
   block_objects <- lapply(seq_along(A), function(j) {
-    create_block(A[[j]], j, bias, na.rm, tau[j], sparsity[j], tol)
+    create_block(
+      A[[j]], j, bias, na.rm, tau[j], sparsity[j], tol, rank[j], mode_orth[j]
+    )
   })
 
   block_objects <- lapply(block_objects, block_init, init = init)
diff --git a/R/rgcca_outer_loop.R b/R/rgcca_outer_loop.R
index 28f8dfcb..63b54c45 100644
--- a/R/rgcca_outer_loop.R
+++ b/R/rgcca_outer_loop.R
@@ -8,7 +8,8 @@ rgcca_outer_loop <- function(blocks, connection = 1 - diag(length(blocks)),
                              verbose = TRUE,
                              na.rm = TRUE, superblock = FALSE,
                              response = NULL, disjunction = NULL,
-                             n_iter_max = 1000, comp_orth = TRUE) {
+                             n_iter_max = 1000, comp_orth = TRUE,
+                             rank = 1, mode_orth = 1) {
   if (verbose) {
     scheme_str <- ifelse(is(scheme, "function"), "user-defined", scheme)
     cat(
@@ -60,22 +61,6 @@ rgcca_outer_loop <- function(blocks, connection = 1 - diag(length(blocks)),
     P <- lapply(seq(J), function(b) c())
   }
 
-  # Save computed shrinkage parameter in a new variable
-  computed_tau <- tau
-  if (is.vector(tau)) {
-    computed_tau <- matrix(
-      rep(tau, N + 1),
-      nrow = N + 1, J, byrow = TRUE
-    )
-  }
-
-  if (is.vector(sparsity)) {
-    sparsity <- matrix(
-      rep(sparsity, N + 1),
-      nrow = N + 1, J, byrow = TRUE
-    )
-  }
-
   # Whether primal or dual
   primal_dual <- matrix("primal", nrow = N + 1, ncol = J)
   primal_dual[which((sparsity == 1) & (nb_ind < matrix(
@@ -91,15 +76,16 @@ rgcca_outer_loop <- function(blocks, connection = 1 - diag(length(blocks)),
       ))
     }
     gcca_result <- rgcca_inner_loop(R, connection, g, dg,
-                                    tau = computed_tau[n, ],
+                                    tau = tau[n, ],
                                     sparsity = sparsity[n, ],
                                     init = init, bias = bias, tol = tol,
                                     verbose = verbose, na.rm = na.rm,
-                                    n_iter_max = n_iter_max
+                                    n_iter_max = n_iter_max,
+                                    rank = rank[n, ], mode_orth = mode_orth
     )
 
     # Store tau, crit
-    computed_tau[n, ] <- gcca_result$tau
+    tau[n, ] <- gcca_result$tau
     crit[[n]] <- gcca_result$crit
 
     # Store Y, a, factors and weights
@@ -137,9 +123,9 @@ rgcca_outer_loop <- function(blocks, connection = 1 - diag(length(blocks)),
   ##### Generation of the output #####
   if (N == 0) {
     crit <- unlist(crit)
-    computed_tau <- as.numeric(computed_tau)
+    tau <- as.numeric(tau)
   } else {
-    computed_tau <- apply(computed_tau, 2, as.numeric)
+    tau <- apply(tau, 2, as.numeric)
   }
 
   astar <- compute_astar(a, P, superblock, comp_orth, N)
@@ -150,7 +136,7 @@ rgcca_outer_loop <- function(blocks, connection = 1 - diag(length(blocks)),
     astar = astar,
     factors = factors,
     weights = weights,
-    tau = computed_tau,
+    tau = tau,
     crit = crit, primal_dual = primal_dual
   )
 
diff --git a/R/rgcca_permutation.R b/R/rgcca_permutation.R
index b004ed02..8cb5db2f 100644
--- a/R/rgcca_permutation.R
+++ b/R/rgcca_permutation.R
@@ -235,7 +235,7 @@ rgcca_permutation <- function(blocks, par_type = "tau", par_value = NULL,
                               response = NULL, superblock = FALSE,
                               NA_method = "na.ignore", rgcca_res = NULL,
                               verbose = TRUE, n_iter_max = 1000,
-                              comp_orth = TRUE) {
+                              comp_orth = TRUE, rank = 1, mode_orth = 1) {
   ### Try to retrieve parameters from a rgcca object
   rgcca_args <- as.list(environment())
   tmp <- get_rgcca_args(blocks, rgcca_args)
@@ -356,7 +356,8 @@ rgcca_permutation <- function(blocks, par_type = "tau", par_value = NULL,
 
   structure(list(
     opt = opt, call = rgcca_args, par_type = par_type,
-    n_perms = n_perms, best_params = param$par_value[which.max(zstat), ],
+    n_perms = n_perms,
+    best_params = param$par_value[which.max(zstat), , drop = FALSE],
     permcrit = permcrit, params = param$par_value, stats = stats
   ), class = "rgcca_permutation")
 }
diff --git a/R/rgcca_stability.R b/R/rgcca_stability.R
index fd460009..661ee05e 100644
--- a/R/rgcca_stability.R
+++ b/R/rgcca_stability.R
@@ -133,8 +133,7 @@ rgcca_stability <- function(rgcca_res,
   rgcca_res$call$blocks <- Map(
     function(x, y) x[, y, drop = FALSE], rgcca_res$call$blocks, keepVar
   )
-  rgcca_res$call$tau <-
-    rgcca_res$call$sparsity <- rep(1, length(rgcca_res$call$blocks))
+  rgcca_res$call$tau[] <- rgcca_res$call$sparsity[] <- 1
 
   rgcca_res <- rgcca(rgcca_res)
 
diff --git a/R/select_analysis.R b/R/select_analysis.R
index 4b9cc27b..2e45815c 100644
--- a/R/select_analysis.R
+++ b/R/select_analysis.R
@@ -21,6 +21,7 @@
 #' @noRd
 select_analysis <- function(rgcca_args, blocks) {
   tau <- rgcca_args$tau
+  rank <- rgcca_args$rank
   ncomp <- rgcca_args$ncomp
   quiet <- rgcca_args$quiet
   scheme <- rgcca_args$scheme
@@ -28,6 +29,7 @@ select_analysis <- function(rgcca_args, blocks) {
   response <- rgcca_args$response
   sparsity <- rgcca_args$sparsity
   comp_orth <- rgcca_args$comp_orth
+  mode_orth <- rgcca_args$mode_orth
   connection <- rgcca_args$connection
   superblock <- rgcca_args$superblock
   scale_block <- rgcca_args$scale_block
@@ -63,6 +65,7 @@ select_analysis <- function(rgcca_args, blocks) {
       penalty <- sparsity
     },
     "tgcca" = {
+      mode_orth <- check_mode_orth(mode_orth, blocks)
       param <- "tau"
       penalty <- tau
       superblock <- FALSE
@@ -410,6 +413,11 @@ select_analysis <- function(rgcca_args, blocks) {
     }
   }
 
+  ncomp <- check_ncomp(
+    ncomp, blocks,
+    superblock = superblock, response = response
+  )
+
   if (method %in% c("rgcca", "sgcca", "tgcca")) {
     scheme <- check_scheme(scheme)
     if (any(sparsity != 1)) {
@@ -419,50 +427,44 @@ select_analysis <- function(rgcca_args, blocks) {
     }
     if (!is.null(response)) {
       check_blockx("response", response, blocks)
-      ncomp[response] <- max(ncomp[-response])
       superblock <- FALSE
       connection <- connection_matrix(
         blocks, type = "response", response = response
       )
     }
-    penalty <- check_penalty(
-      penalty, blocks, method, superblock = superblock, ncomp = max(ncomp)
-    )
-    if (superblock) {
-      ncomp <- rep(max(ncomp), J + 1)
-      connection <- connection_matrix(blocks, type = "response", J = J + 1)
-      if (is.matrix(penalty)) {
-        if (ncol(penalty) < J + 1) {
-          pen <- 1
-        } else {
-          pen <- penalty[, J + 1]
-        }
-        penalty <- cbind(penalty[, seq(J)], pen)
-      } else {
-        pen <- ifelse(length(penalty) < J + 1, 1, penalty[J + 1])
-        penalty <- c(penalty[seq(J)], pen)
-      }
-    } else {
-      if (is.null(connection)) {
-        connection <- connection_matrix(blocks, type = "pair")
-      } else {
-        connection <- check_connection(connection, blocks)
-      }
+    rank <- check_rank(rank, blocks, mode_orth, ncomp = max(ncomp))
+    if (is.null(connection)) {
+      connection <- connection_matrix(blocks, type = "pair")
     }
   }
-  ncomp <- check_ncomp(
-    ncomp, blocks,
-    superblock = superblock, response = response
+
+  penalty <- check_penalty(
+    penalty, blocks, method, superblock = superblock, ncomp = max(ncomp)
   )
+  if (superblock) {
+    connection <- connection_matrix(blocks, type = "response", J = J + 1)
+    if (ncol(penalty) < J + 1) {
+      pen <- 1
+    } else {
+      pen <- penalty[, J + 1]
+    }
+    penalty <- cbind(penalty[, seq(J), drop = FALSE], pen)
+  } else {
+    connection <- check_connection(connection, blocks)
+  }
 
-  rgcca_args[[param]] <- penalty
+  other_param <- ifelse(param == "tau", "sparsity", "tau")
+  rgcca_args[[param]] <- rgcca_args[[other_param]] <- penalty
+  rgcca_args[[other_param]][] <- 1
 
   rgcca_args <- modifyList(rgcca_args, list(
+    rank = rank,
     ncomp = ncomp,
     scheme = scheme,
     method = method,
     response = response,
     comp_orth = comp_orth,
+    mode_orth = mode_orth,
     connection = connection,
     superblock = superblock,
     scale_block = scale_block
diff --git a/R/subset_rows.R b/R/subset_rows.R
index 49408036..1e17ee9c 100644
--- a/R/subset_rows.R
+++ b/R/subset_rows.R
@@ -2,7 +2,6 @@
 #' data.frame).
 #' @param x An object from which we want to extract rows
 #' @param rows A set of rows
-#' @importFrom stats complete.cases
 #' @noRd
 subset_rows <- function(x, rows) {
   is.x.data.frame <- is.data.frame(x)
diff --git a/R/summary.rgcca.R b/R/summary.rgcca.R
index c2f4af7c..64d9ed8a 100644
--- a/R/summary.rgcca.R
+++ b/R/summary.rgcca.R
@@ -109,7 +109,7 @@ summary.rgcca <- function(object, ...) {
   cat("\n")
   if (!tolower(object$call$method) %in% sparse_methods()) {
     param <- "regularization"
-    if (!is.matrix(object$call$tau)) {
+    if (!is.matrix(drop(object$call$tau))) {
       for (i in seq_len(NCOL(object$call$connection))) {
         tau <- object$call$tau[i]
         cat("The", param, "parameter used for", names(object$blocks)[i],
@@ -130,7 +130,7 @@ summary.rgcca <- function(object, ...) {
       function(a) apply(a, 2, function(l) sum(l != 0))
     )
     param <- "sparsity"
-    if (!is.matrix(object$call$sparsity)) {
+    if (!is.matrix(drop(object$call$sparsity))) {
       for (i in seq_len(NCOL(object$call$connection))[-response]) {
         sparsity <- object$call$sparsity[i]
 
diff --git a/R/summary.rgcca_cv.r b/R/summary.rgcca_cv.r
index 60395f1b..2718b6c2 100644
--- a/R/summary.rgcca_cv.r
+++ b/R/summary.rgcca_cv.r
@@ -42,7 +42,7 @@ summary.rgcca_cv <- function(object, type = c("sd", "quantile"), ...) {
   cat("\n")
 
   best <- which(apply(
-    object$params, 1, function(z) identical(z, object$best_params)
+    object$params, 1, function(z) identical(z, drop(object$best_params))
   ))
   optimal_y <- object$stats[best, "mean"]
 
diff --git a/R/summary.rgcca_permutation.R b/R/summary.rgcca_permutation.R
index 2f00a35b..c66a1a6b 100644
--- a/R/summary.rgcca_permutation.R
+++ b/R/summary.rgcca_permutation.R
@@ -22,7 +22,7 @@ summary.rgcca_permutation <- function(object, ...) {
   print(tab, quote = FALSE, ...)
 
   best <- which(apply(
-    object$params, 1, function(z) identical(z, object$best_params)
+    object$params, 1, function(z) identical(z, drop(object$best_params))
   ))
   cat(strwrap(paste0(
     "\nThe best combination is: ",
diff --git a/man/rgcca.Rd b/man/rgcca.Rd
index 4e3c3df4..3d93760b 100644
--- a/man/rgcca.Rd
+++ b/man/rgcca.Rd
@@ -24,6 +24,8 @@ rgcca(
   quiet = TRUE,
   n_iter_max = 1000,
   comp_orth = TRUE,
+  rank = 1,
+  mode_orth = 1,
   A = NULL,
   C = NULL
 )
@@ -165,6 +167,20 @@ iterations.}
 \item{comp_orth}{A logical value indicating if the deflation should lead to
 orthogonal block components or orthogonal block weight vectors.}
 
+\item{rank}{Either an integer, an integer vector of
+size \eqn{J} or an integer matrix
+of dimension \eqn{\textrm{max}(\textrm{ncomp}) \times J} giving the rank
+of the decomposition sought for the canonical vectors in TGCCA.
+If block \eqn{j} is an array with at least three dimensions, rank must be
+comprised between 1 and the number of variables on the mode bearing the
+orthogonality constraint. See \textrm{mode_orth}.}
+
+\item{mode_orth}{Either an integer or an integer vector of size \eqn{J}
+designating the mode which associated set of factors will be orthogonal
+in the decomposition sought by TGCCA. If block \eqn{j} is an array with
+\eqn{d > 2} dimensions, \textrm{mode_orth} must be comprised between
+1 and \eqn{d - 1}.}
+
 \item{A}{Deprecated argument, please use blocks instead.}
 
 \item{C}{Deprecated argument, please use connection instead.}
diff --git a/man/rgcca_cv.Rd b/man/rgcca_cv.Rd
index 2ae5c03c..dac111d9 100644
--- a/man/rgcca_cv.Rd
+++ b/man/rgcca_cv.Rd
@@ -33,6 +33,8 @@ rgcca_cv(
   verbose = TRUE,
   n_iter_max = 1000,
   comp_orth = TRUE,
+  rank = 1,
+  mode_orth = 1,
   ...
 )
 }
@@ -219,6 +221,20 @@ iterations.}
 \item{comp_orth}{A logical value indicating if the deflation should lead to
 orthogonal block components or orthogonal block weight vectors.}
 
+\item{rank}{Either an integer, an integer vector of
+size \eqn{J} or an integer matrix
+of dimension \eqn{\textrm{max}(\textrm{ncomp}) \times J} giving the rank
+of the decomposition sought for the canonical vectors in TGCCA.
+If block \eqn{j} is an array with at least three dimensions, rank must be
+comprised between 1 and the number of variables on the mode bearing the
+orthogonality constraint. See \textrm{mode_orth}.}
+
+\item{mode_orth}{Either an integer or an integer vector of size \eqn{J}
+designating the mode which associated set of factors will be orthogonal
+in the decomposition sought by TGCCA. If block \eqn{j} is an array with
+\eqn{d > 2} dimensions, \textrm{mode_orth} must be comprised between
+1 and \eqn{d - 1}.}
+
 \item{...}{Additional parameters to be passed to prediction_model.}
 }
 \value{
diff --git a/man/rgcca_permutation.Rd b/man/rgcca_permutation.Rd
index d8e528a1..53d4466d 100644
--- a/man/rgcca_permutation.Rd
+++ b/man/rgcca_permutation.Rd
@@ -29,7 +29,9 @@ rgcca_permutation(
   rgcca_res = NULL,
   verbose = TRUE,
   n_iter_max = 1000,
-  comp_orth = TRUE
+  comp_orth = TRUE,
+  rank = 1,
+  mode_orth = 1
 )
 }
 \arguments{
@@ -199,6 +201,20 @@ iterations.}
 
 \item{comp_orth}{A logical value indicating if the deflation should lead to
 orthogonal block components or orthogonal block weight vectors.}
+
+\item{rank}{Either an integer, an integer vector of
+size \eqn{J} or an integer matrix
+of dimension \eqn{\textrm{max}(\textrm{ncomp}) \times J} giving the rank
+of the decomposition sought for the canonical vectors in TGCCA.
+If block \eqn{j} is an array with at least three dimensions, rank must be
+comprised between 1 and the number of variables on the mode bearing the
+orthogonality constraint. See \textrm{mode_orth}.}
+
+\item{mode_orth}{Either an integer or an integer vector of size \eqn{J}
+designating the mode which associated set of factors will be orthogonal
+in the decomposition sought by TGCCA. If block \eqn{j} is an array with
+\eqn{d > 2} dimensions, \textrm{mode_orth} must be comprised between
+1 and \eqn{d - 1}.}
 }
 \value{
 A rgcca_permutation object that can be printed and plotted.
diff --git a/tests/testthat/test_checks.r b/tests/testthat/test_checks.r
index 9a07fb51..6b60adc3 100644
--- a/tests/testthat/test_checks.r
+++ b/tests/testthat/test_checks.r
@@ -393,27 +393,31 @@ test_that("check_penalty raises an error for invalid penalty", {
   expect_error(check_penalty("optimal", blocks, method = "sgcca"))
 })
 test_that("check_penalty passes and returns penalty when penalty is valid", {
-  expect_equal(check_penalty(1, blocks, method = "rgcca"), c(1, 1, 1))
+  expect_equal(check_penalty(1, blocks, method = "rgcca"), matrix(1, 1, 3))
   expect_equal(
-    check_penalty(c(0.8, 1, 0.5), blocks, method = "rgcca"), c(0.8, 1, 0.5)
+    check_penalty(c(0.8, 1, 0.5), blocks, method = "rgcca"),
+    matrix(c(0.8, 1, 0.5), 1, 3)
   )
   expect_equal(
     check_penalty("optimal", blocks, method = "rgcca"),
-    c("optimal", "optimal", "optimal")
+    matrix("optimal", 1, 3)
   )
   expect_equal(
-    check_penalty(matrix(1, 5, 3), blocks, method = "rgcca"), matrix(1, 5, 3)
+    check_penalty(matrix(1, 5, 3), blocks, method = "rgcca", ncomp = 5),
+    matrix(1, 5, 3)
   )
   expect_equal(
     check_penalty(rep(1, 4), blocks, method = "rgcca", superblock = TRUE),
-    rep(1, 4)
+    matrix(1, 1, 4)
   )
-  expect_equal(check_penalty(1, blocks, method = "sgcca"), c(1, 1, 1))
+  expect_equal(check_penalty(1, blocks, method = "sgcca"), matrix(1, 1, 3))
   expect_equal(
-    check_penalty(c(0.8, 1, 0.5), blocks, method = "sgcca"), c(0.8, 1, 0.5)
+    check_penalty(c(0.8, 1, 0.5), blocks, method = "sgcca"),
+    matrix(c(0.8, 1, 0.5), 1, 3)
   )
   expect_equal(
-    check_penalty(matrix(1, 5, 3), blocks, method = "sgcca"), matrix(1, 5, 3)
+    check_penalty(matrix(1, 5, 3), blocks, method = "sgcca", ncomp = 5),
+    matrix(1, 5, 3)
   )
 })
 
@@ -458,3 +462,37 @@ test_that("check_scheme passes when scheme is valid", {
   g <- function(x) x^2
   expect_error(check_scheme(g), NA)
 })
+
+# Test check_rank
+n <- nrow(blocks[[1]])
+array_blocks <- blocks
+array_blocks[[4]] <- array(seq(n * 10 * 5), dim = c(n, 10, 5))
+test_that("check_rank raises an error for invalid rank", {
+  expect_error(
+    check_rank("toto", array_blocks, rep(1, 4), 1),
+    "rank must be numeric.", fixed = TRUE
+  )
+  expect_error(
+    check_rank(rep(1, 5), array_blocks, rep(1, 4), 1),
+    paste0("rank must have the same number of columns (actually 5) ",
+           "as the number of blocks (4)."),
+    fixed = TRUE
+  )
+  expect_error(
+    check_rank(c(1, 1, 1, 7), array_blocks, c(1, 1, 1, 2), 1),
+    paste0("rank[4] should be comprise between 1 and 5 (that is the number of ",
+           "variables of the mode bearing the orthogonality for block 4)."),
+    fixed = TRUE
+  )
+  expect_error(
+    check_rank(matrix(rep(1, 8), nrow = 2), array_blocks, rep(1, 4), 1),
+    "rank must have 1 rows.", fixed = TRUE
+  )
+})
+test_that("check_rank passes when rank is valid", {
+  expect_equal(check_rank(1, array_blocks, rep(1, 4), 1), matrix(1, 1, 4))
+  r <- matrix(1, 2, 4)
+  expect_equal(check_rank(r, array_blocks, rep(1, 4), 2), r)
+  r[2, 4] <- 7
+  expect_equal(check_rank(r, array_blocks, rep(1, 4), 2), r)
+})
diff --git a/tests/testthat/test_select_analysis.R b/tests/testthat/test_select_analysis.R
index 82690d61..981ebfb3 100644
--- a/tests/testthat/test_select_analysis.R
+++ b/tests/testthat/test_select_analysis.R
@@ -13,12 +13,14 @@ run_selection <- function(method, quiet = TRUE, ...) {
 
   rgcca_args <- list(
     tau = rep(1, J),
+    rank = rep(1, J),
     ncomp = rep(1, J),
     quiet = quiet,
     scheme = "centroid",
     method = method,
     response = NULL,
     sparsity = rep(1, J),
+    mode_orth = rep(1, J),
     connection = 1 - diag(J),
     superblock = FALSE
   )
@@ -72,7 +74,7 @@ test_that("superblock methods sets all attributes of a superblock", {
   method <- "spca"
   res <- run_selection(method, superblock = TRUE, sparsity = 0.7)$res
   expect_equal(res$rgcca_args$method, "spca")
-  expect_equal(res$rgcca_args[[res$opt$param]], c(0.7, 0.7))
+  expect_equal(unname(res$rgcca_args[[res$opt$param]]), matrix(0.7, 1, 2))
   expect_equal(res$opt$param, "sparsity")
 })