Create svd_wrapper and copy ginv to avoid LAPACK errors

DESCRIPTION

@@ -34,7 +34,6 @@ Imports:
     ggrepel,
     graphics,
     gridExtra,
-    MASS,
     matrixStats,
     methods,
     parallel,
@@ -45,6 +44,7 @@ Suggests:
     devtools,
     FactoMineR,
     knitr,
+    MASS,
     pander,
     rmarkdown,
     rticles,

NAMESPACE

@@ -41,7 +41,6 @@ export(rgcca_predict)
 export(rgcca_stability)
 export(rgcca_transform)
 importFrom(Deriv,Deriv)
-importFrom(MASS,ginv)
 importFrom(caret,confusionMatrix)
 importFrom(caret,multiClassSummary)
 importFrom(caret,postResample)

R/block_init.R

@@ -1,5 +1,3 @@
-#' @importFrom MASS ginv
-
 block_init <- function(x, init = "svd") {
   UseMethod("block_init")
 }

R/block_project.R

@@ -15,7 +15,12 @@ block_project.block <- function(x) {
 #' @export
 block_project.dual_block <- function(x) {
   if (any(x$alpha != 0)) {
-    x$alpha <- x$alpha / drop(sqrt(t(x$alpha) %*% x$K %*% x$alpha))
+    alpha_norm <- t(x$alpha) %*% x$K %*% x$alpha
+    if (alpha_norm > 0) {
+      x$alpha <- x$alpha / drop(sqrt(alpha_norm))
+    } else {
+      x$alpha <- x$alpha * 0
+    }
   }
   x$a <- pm(t(x$x), x$alpha, na.rm = x$na.rm)
 
@@ -26,7 +31,12 @@ block_project.dual_block <- function(x) {
 #' @export
 block_project.primal_regularized_block <- function(x) {
   if (any(x$a != 0)) {
-    x$a <- x$M %*% x$a / drop(sqrt(t(x$a) %*% x$M %*% x$a))
+    a_norm <- t(x$a) %*% x$M %*% x$a
+    if (a_norm > 0) {
+      x$a <- x$M %*% x$a / drop(sqrt(a_norm))
+    } else {
+      x$a <- x$a * 0
+    }
   }
 
   x$Y <- pm(x$x, x$a, na.rm = x$na.rm)
@@ -36,9 +46,12 @@ block_project.primal_regularized_block <- function(x) {
 #' @export
 block_project.dual_regularized_block <- function(x) {
   if (any(x$alpha != 0)) {
-    x$alpha <- x$M %*% x$alpha / drop(sqrt(
-      t(x$alpha) %*% x$M %*% x$K %*% x$alpha
-    ))
+    alpha_norm <- t(x$alpha) %*% x$M %*% x$K %*% x$alpha
+    if (alpha_norm > 0) {
+      x$alpha <- x$M %*% x$alpha / drop(sqrt(alpha_norm))
+    } else {
+      x$alpha <- x$alpha * 0
+    }
   }
   x$a <- pm(t(x$x), x$alpha, na.rm = x$na.rm)
 

R/deflation.R

@@ -1,11 +1,15 @@
 deflation <- function(X, y, na.rm = TRUE, left = TRUE) {
   # Computation of the residual matrix R
   # Computation of the vector p.
+  y_norm <- drop(crossprod(y))
+  if (y_norm == 0) {
+    y_norm <- 1
+  }
   if (left) {
-    p <- pm(t(X), y, na.rm = na.rm) / as.vector(crossprod(y))
+    p <- pm(t(X), y, na.rm = na.rm) / y_norm
     R <- X - tcrossprod(y, p)
   } else {
-    p <- pm(X, y, na.rm = na.rm) / as.vector(crossprod(y))
+    p <- pm(X, y, na.rm = na.rm) / y_norm
     R <- X - tcrossprod(p, y)
   }
   return(list(p = p, R = R))

R/ginv.R

@@ -0,0 +1,20 @@
+# Copy of the ginv function from the MASS package.
+# We replace calls to svd by our wrapper to avoid LAPACK errors.
+ginv <- function (X, tol = sqrt(.Machine$double.eps))
+{
+  if (length(dim(X)) > 2L || !(is.numeric(X) || is.complex(X)))
+    stop("'X' must be a numeric or complex matrix")
+  if (!is.matrix(X))
+    X <- as.matrix(X)
+  Xsvd <- svd_wrapper(X)
+  if (is.complex(X))
+    Xsvd$u <- Conj(Xsvd$u)
+  Positive <- Xsvd$d > max(tol * Xsvd$d[1L], 0)
+  if (all(Positive))
+    Xsvd$v %*% (1/Xsvd$d * t(Xsvd$u))
+  else if (!any(Positive))
+    array(0, dim(X)[2L:1L])
+  else Xsvd$v[, Positive, drop = FALSE] %*% (
+    (1/Xsvd$d[Positive]) * t(Xsvd$u[, Positive, drop = FALSE])
+  )
+}

R/initsvd.R

@@ -21,10 +21,10 @@ initsvd <- function(X, dual = TRUE) {
 
   if (dual) {
     ifelse(n >= p,
-      return(svd(X, nu = 0, nv = 1)$v),
-      return(svd(X, nu = 1, nv = 0)$u)
+      return(svd_wrapper(X, nu = 0, nv = 1)$v),
+      return(svd_wrapper(X, nu = 1, nv = 0)$u)
     )
   } else {
-    return(svd(X, nu = 0, nv = 1)$v)
+    return(svd_wrapper(X, nu = 0, nv = 1)$v)
   }
 }

R/rgcca_predict.R

@@ -105,7 +105,9 @@ rgcca_predict <- function(rgcca_res,
 
   ### Get projected train and test data
   projection <- rgcca_transform(rgcca_res, blocks_test[-test_idx])
-  X_train <- rgcca_res$Y[names(projection)]
+  X_train <- rgcca_transform(
+    rgcca_res, rgcca_res$call$blocks[names(projection)]
+  )
   X_train <- reformat_projection(X_train)
   X_test <- reformat_projection(projection)
 

R/rgcca_transform.R

@@ -86,6 +86,10 @@ rgcca_transform <- function(rgcca_res, blocks_test = rgcca_res$call$blocks) {
   # Otherwise we directly use astar to project the individual blocks
   } else {
     astar <- rgcca_res$astar[names(X_train)]
+    # Remove zero columns of astar
+    astar <- lapply(astar, function(x) {
+      x[, which(apply(x, 2, function(y) sum(abs(y)) > 0))]
+    })
     projection <- lapply(seq_along(blocks_test), function(j) {
       x <- pm(as.matrix(blocks_test[[j]]), astar[[j]])
       rownames(x) <- rownames(blocks_test[[j]])

R/svd_wrapper.R

@@ -0,0 +1,25 @@
+# workaround by Art Owen to avoid LAPACK errors
+# See https://stat.ethz.ch/pipermail/r-help/2007-October/143508.html
+svd_wrapper <- function(x, nu = min(n, p), nv = min(n, p), ...) {
+  success <- FALSE
+  n <- NROW(x)
+  p <- NCOL(x)
+  try({
+    svd_x <- base::svd(x, nu, nv, ...)
+    success <- TRUE
+  }, silent = TRUE)
+  if( success ) {
+    return(svd_x)
+  }
+  try( {
+    svd_tx <- base::svd(t(x), nv, nu, ...)
+    success <- TRUE
+  }, silent = TRUE )
+  if( !success ) {
+    stop("Error: svd(x) and svd(t(x)) both failed to converge.")
+  }
+  temp <- svd_tx$u
+  svd_tx$u <- svd_tx$v
+  svd_tx$v <- temp
+  return(svd_tx)
+}

codemeta.json

@@ -112,6 +112,11 @@
       },
       "sameAs": "https://CRAN.R-project.org/package=knitr"
     },
+    {
+      "@type": "SoftwareApplication",
+      "identifier": "MASS",
+      "name": "MASS"
+    },
     {
       "@type": "SoftwareApplication",
       "identifier": "pander",
@@ -247,18 +252,6 @@
       "sameAs": "https://CRAN.R-project.org/package=gridExtra"
     },
     "8": {
-      "@type": "SoftwareApplication",
-      "identifier": "MASS",
-      "name": "MASS",
-      "provider": {
-        "@id": "https://cran.r-project.org",
-        "@type": "Organization",
-        "name": "Comprehensive R Archive Network (CRAN)",
-        "url": "https://cran.r-project.org"
-      },
-      "sameAs": "https://CRAN.R-project.org/package=MASS"
-    },
-    "9": {
       "@type": "SoftwareApplication",
       "identifier": "matrixStats",
       "name": "matrixStats",

tests/testthat/test_defl_select.r

@@ -21,8 +21,17 @@ test_that("defl_select does not deflate block which reached ncomp", {
   expect_equal(res$resdefl[[1]], A[[1]])
 })
 
+test_that("defl_select does not deflate block when y is zero", {
+  yy0 <- yy
+  yy0[[1]] <- yy0[[1]] * 0
+  res <- defl_select(
+    yy = yy0, rr = A, nncomp = c(2, 2, 2), nn = 1, nbloc = 3
+  )
+  expect_equal(res$resdefl[[1]], A[[1]])
+})
+
 test_that("defl_select outputs coherent residuals and projections", {
-  res <- defl_select(yy = yy, rr = A, nncomp = c(1, 1, 1), nn = 1, nbloc = 3)
+  res <- defl_select(yy = yy, rr = A, nncomp = c(2, 2, 2), nn = 1, nbloc = 3)
   for (j in seq_along(A)) {
     expect_equal(A[[j]] - yy[[j]] %*% t(res$pdefl[[j]]), res$resdefl[[j]])
   }

tests/testthat/test_rgcca_transform.R

@@ -163,3 +163,17 @@ test_that("rgcca_transform creates projection with the right number of
     expect_true(all(dim(projection[[j]]) == c(nrow(projection[[j]]), ncomp[j])))
   }
 })
+
+#-------------------------------------------------------------------------
+# Checking rgcca_transform removes zero columns
+#-------------------------------------------------------------------------
+fit.rgcca_with_zeros <- fit.rgcca
+fit.rgcca_with_zeros$astar[[1]][, 3] <- 0
+fit.rgcca_with_zeros$astar[[3]][, c(3, 4)] <- 0
+
+projection <- rgcca_transform(fit.rgcca_with_zeros, A_test)
+
+test_that("rgcca_transform creates projection without zero columns", {
+  expect_equal(ncol(projection[[1]]), 2)
+  expect_equal(ncol(projection[[3]]), 2)
+})