fix(svd): keep nonneg/L1 constraints on deflation factors after the first

NEWS.md

@@ -1,3 +1,21 @@
+# RcppML (development version)
+
+## Bug fixes
+
+- `svd()`/`pca()` with `method = "deflation"`: factors with projection
+  constraints (`nonneg`, `L1`, `upper_bound`) no longer go through the
+  post-factor Gram-Schmidt step, which re-introduced negative entries and
+  removed L1 zeros on every component after the first. Constrained factors are
+  now returned as fitted (non-negative/sparse, not mutually orthogonal);
+  unconstrained SVD/PCA is unchanged. Applies to the CPU and GPU deflation paths.
+- Constrained factors after the first now start from a non-negative random
+  vector instead of the orthogonalized power-step warm start, which could
+  collapse to zero after projection.
+- `method = "auto"` now selects `"deflation"` whenever `nonneg`, `L1`, or
+  `upper_bound` is set (previously `"krylov"` for `k >= 8`), and requesting
+  `method = "krylov"` with these constraints and `k > 1` warns that components
+  may be redundant or all-zero.
+
 # RcppML 1.0.0
 
 ## Breaking Changes

R/svd.R

@@ -285,6 +285,9 @@ svd <- function(A, k = 10, tol = 1e-5, maxit = 200,
 
   has_element_constraints <- any(L1 > 0) || any(L2 > 0) || any(nonneg) ||
                              any(upper_bound > 0) || any(L21 > 0)
+  # Constraints enforced by projection (sign / sparsity / bounds); these are
+  # incompatible with orthogonalizing factors after the fact.
+  has_projection_constraints <- any(nonneg) || any(L1 > 0) || any(upper_bound > 0)
   has_gram_constraints <- any(angular > 0) ||
                           (!is.null(graph_U) && graph_lambda[1] > 0) ||
                           (!is.null(graph_V) && graph_lambda[2] > 0)
@@ -304,6 +307,14 @@ svd <- function(A, k = 10, tol = 1e-5, maxit = 200,
         method))
     }
 
+    if (method == "krylov" && has_projection_constraints &&
+        (isTRUE(auto_rank) || (is.numeric(k) && k > 1))) {
+      warning("method = 'krylov' with nonneg/L1/upper_bound constraints and k > 1 ",
+              "can return redundant or all-zero components; ",
+              "method = 'deflation' (the default for constrained problems) is recommended.",
+              call. = FALSE)
+    }
+
     if (has_gram_constraints && !method %in% gram_methods) {
       stop(sprintf(
         "method '%s' does not support Gram-level regularization (angular/graph). Use 'krylov'.",
@@ -378,6 +389,11 @@ svd <- function(A, k = 10, tol = 1e-5, maxit = 200,
     has_any_constraints <- has_element_constraints || has_gram_constraints
     if (has_robust) {
       method <- "deflation"
+    } else if (has_projection_constraints) {
+      # Sign/sparsity constraints: Krylov's block solve-then-project can
+      # return redundant or all-zero components for k > 1; deflation keeps
+      # each factor's constraints and fits it against the residual.
+      method <- "deflation"
     } else if (has_any_constraints) {
       method <- if (k >= 8) "krylov" else "deflation"
     } else if (has_cv || auto_rank) {

inst/include/FactorNet/svd/deflation.hpp

@@ -619,6 +619,13 @@ SVDResult<Scalar> deflation_svd(const MatrixType& A,
         col_weights.setOnes();
     }
 
+    // Element-wise constraints on a side make that side's factors
+    // incompatible with Gram-Schmidt (see post-factor reorthogonalization).
+    const bool u_constrained = config.nonneg_u || config.L1_u > 0 ||
+                               config.upper_bound_u > 0;
+    const bool v_constrained = config.nonneg_v || config.L1_v > 0 ||
+                               config.upper_bound_v > 0;
+
     for (int k = 0; k < k_max; ++k) {
         SVD_CHECK_INTERRUPT();
 
@@ -627,8 +634,11 @@ SVDResult<Scalar> deflation_svd(const MatrixType& A,
         // Power-step warm start gives good initial direction
         // Nesterov momentum halves iteration count for clustered SVs
         // ============================================================
-
-        if (k == 0) {
+
+        // Constrained u: the orthogonalized power-step warm start can point
+        // into the negative orthant and collapse to zero after projection,
+        // so use a fresh non-negative random start instead.
+        if (k == 0 || u_constrained) {
             // Random initialization for first factor
             for (int i = 0; i < m; ++i) {
                 u(i) = static_cast<Scalar>(rng.uniform<double>());
@@ -824,28 +834,37 @@ SVDResult<Scalar> deflation_svd(const MatrixType& A,
         // Two-pass GS ensures orthogonality drift stays below eps
         // at high k. Cost: O(k*(m+n)) — negligible vs O(nnz) SpMV.
         // ============================================================
+        // Gram-Schmidt is skipped on any side with projection constraints
+        // (nonneg, L1, upper bound): subtracting projections onto prior
+        // factors re-introduces negative entries and destroys sparsity, so
+        // constrained factors are returned as fitted (non-orthogonal), as in
+        // other constrained/sparse PCA methods.
         if (k > 0) {
             auto U_prev = U_all.leftCols(k);
             auto V_prev = V_all.leftCols(k);
 
             // Reorthogonalize u against U_all[:, 0:k-1] (two passes)
-            for (int pass = 0; pass < 2; ++pass) {
-                DenseVector<Scalar> dots_u = U_prev.transpose() * u;
-                u.noalias() -= U_prev * dots_u;
-            }
-            Scalar u_norm_gs = u.norm();
-            if (u_norm_gs > std::numeric_limits<Scalar>::epsilon() * 100) {
-                u /= u_norm_gs;
+            if (!u_constrained) {
+                for (int pass = 0; pass < 2; ++pass) {
+                    DenseVector<Scalar> dots_u = U_prev.transpose() * u;
+                    u.noalias() -= U_prev * dots_u;
+                }
+                Scalar u_norm_gs = u.norm();
+                if (u_norm_gs > std::numeric_limits<Scalar>::epsilon() * 100) {
+                    u /= u_norm_gs;
+                }
             }
 
             // Reorthogonalize v against V_all[:, 0:k-1] (two passes)
-            for (int pass = 0; pass < 2; ++pass) {
-                DenseVector<Scalar> dots_v = V_prev.transpose() * v;
-                v.noalias() -= V_prev * dots_v;
-            }
-            Scalar v_norm_gs = v.norm();
-            if (v_norm_gs > std::numeric_limits<Scalar>::epsilon() * 100) {
-                v /= v_norm_gs;
+            if (!v_constrained) {
+                for (int pass = 0; pass < 2; ++pass) {
+                    DenseVector<Scalar> dots_v = V_prev.transpose() * v;
+                    v.noalias() -= V_prev * dots_v;
+                }
+                Scalar v_norm_gs = v.norm();
+                if (v_norm_gs > std::numeric_limits<Scalar>::epsilon() * 100) {
+                    v /= v_norm_gs;
+                }
             }
         }
 

inst/include/FactorNet/svd/deflation_gpu.cuh

@@ -636,6 +636,12 @@ SVDResult<Scalar> deflation_svd_gpu(
     int m, int n, int nnz,
     const SVDConfig<Scalar>& config)
 {
+    // Projection constraints (nonneg/L1/bounds) are incompatible with
+    // post-factor Gram-Schmidt; see deflation.hpp.
+    const bool u_constrained = config.nonneg_u || config.L1_u > 0 ||
+                               config.upper_bound_u > 0;
+    const bool v_constrained = config.nonneg_v || config.L1_v > 0 ||
+                               config.upper_bound_v > 0;
     using namespace FactorNet::gpu;
     using DenseMatrixS = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
     using DenseVectorS = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
@@ -1253,27 +1259,31 @@ SVDResult<Scalar> deflation_svd_gpu(
             Scalar* d_gs_h = d_defl_workspace.get();
 
             // Reorthogonalize d_u against d_U_factors[:, 0:k-1]
-            detail::gpu_orthog_against<Scalar>(
-                ctx.cublas, d_U_factors.get(), d_u.get(), m, k, d_gs_h);
-            // Renormalize d_u
-            {
-                detail::cublas_nrm2_device<Scalar>(
-                    ctx.cublas, m, d_u.get(), d_scalars.get());
-                int blocks = (m + BLOCK_SIZE - 1) / BLOCK_SIZE;
-                detail::fused_normalize_kernel<<<blocks, BLOCK_SIZE, 0, ctx.stream>>>(
-                    d_u.get(), d_u.get(), d_scalars.get(), m);
+            if (!u_constrained) {
+                detail::gpu_orthog_against<Scalar>(
+                    ctx.cublas, d_U_factors.get(), d_u.get(), m, k, d_gs_h);
+                // Renormalize d_u
+                {
+                    detail::cublas_nrm2_device<Scalar>(
+                        ctx.cublas, m, d_u.get(), d_scalars.get());
+                    int blocks = (m + BLOCK_SIZE - 1) / BLOCK_SIZE;
+                    detail::fused_normalize_kernel<<<blocks, BLOCK_SIZE, 0, ctx.stream>>>(
+                        d_u.get(), d_u.get(), d_scalars.get(), m);
+                }
             }
 
             // Reorthogonalize d_v against d_V_factors[:, 0:k-1]
-            detail::gpu_orthog_against<Scalar>(
-                ctx.cublas, d_V_factors.get(), d_v.get(), n, k, d_gs_h);
-            // Renormalize d_v
-            {
-                detail::cublas_nrm2_device<Scalar>(
-                    ctx.cublas, n, d_v.get(), d_scalars.get());
-                int blocks = (n + BLOCK_SIZE - 1) / BLOCK_SIZE;
-                detail::fused_normalize_kernel<<<blocks, BLOCK_SIZE, 0, ctx.stream>>>(
-                    d_v.get(), d_v.get(), d_scalars.get(), n);
+            if (!v_constrained) {
+                detail::gpu_orthog_against<Scalar>(
+                    ctx.cublas, d_V_factors.get(), d_v.get(), n, k, d_gs_h);
+                // Renormalize d_v
+                {
+                    detail::cublas_nrm2_device<Scalar>(
+                        ctx.cublas, n, d_v.get(), d_scalars.get());
+                    int blocks = (n + BLOCK_SIZE - 1) / BLOCK_SIZE;
+                    detail::fused_normalize_kernel<<<blocks, BLOCK_SIZE, 0, ctx.stream>>>(
+                        d_v.get(), d_v.get(), d_scalars.get(), n);
+                }
             }
         }
 
@@ -1432,6 +1442,12 @@ SVDResult<Scalar> deflation_svd_gpu_dense(
     int m, int n,
     const SVDConfig<Scalar>& config)
 {
+    // Projection constraints (nonneg/L1/bounds) are incompatible with
+    // post-factor Gram-Schmidt; see deflation.hpp.
+    const bool u_constrained = config.nonneg_u || config.L1_u > 0 ||
+                               config.upper_bound_u > 0;
+    const bool v_constrained = config.nonneg_v || config.L1_v > 0 ||
+                               config.upper_bound_v > 0;
     using namespace FactorNet::gpu;
     using DenseMatrixS = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
     using DenseVectorS = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
@@ -1797,23 +1813,27 @@ SVDResult<Scalar> deflation_svd_gpu_dense(
         // Post-factor Gram-Schmidt reorthogonalization (same as sparse path)
         if (k > 0) {
             Scalar* d_gs_h = d_defl_workspace.get();
-            detail::gpu_orthog_against<Scalar>(
-                ctx.cublas, d_U_factors.get(), d_u.get(), m, k, d_gs_h);
-            {
-                detail::cublas_nrm2_device<Scalar>(
-                    ctx.cublas, m, d_u.get(), d_scalars.get());
-                int blocks = (m + BLOCK_SIZE - 1) / BLOCK_SIZE;
-                detail::fused_normalize_kernel<<<blocks, BLOCK_SIZE, 0, ctx.stream>>>(
-                    d_u.get(), d_u.get(), d_scalars.get(), m);
+            if (!u_constrained) {
+                detail::gpu_orthog_against<Scalar>(
+                    ctx.cublas, d_U_factors.get(), d_u.get(), m, k, d_gs_h);
+                {
+                    detail::cublas_nrm2_device<Scalar>(
+                        ctx.cublas, m, d_u.get(), d_scalars.get());
+                    int blocks = (m + BLOCK_SIZE - 1) / BLOCK_SIZE;
+                    detail::fused_normalize_kernel<<<blocks, BLOCK_SIZE, 0, ctx.stream>>>(
+                        d_u.get(), d_u.get(), d_scalars.get(), m);
+                }
             }
-            detail::gpu_orthog_against<Scalar>(
-                ctx.cublas, d_V_factors.get(), d_v.get(), n, k, d_gs_h);
-            {
-                detail::cublas_nrm2_device<Scalar>(
-                    ctx.cublas, n, d_v.get(), d_scalars.get());
-                int blocks = (n + BLOCK_SIZE - 1) / BLOCK_SIZE;
-                detail::fused_normalize_kernel<<<blocks, BLOCK_SIZE, 0, ctx.stream>>>(
-                    d_v.get(), d_v.get(), d_scalars.get(), n);
+            if (!v_constrained) {
+                detail::gpu_orthog_against<Scalar>(
+                    ctx.cublas, d_V_factors.get(), d_v.get(), n, k, d_gs_h);
+                {
+                    detail::cublas_nrm2_device<Scalar>(
+                        ctx.cublas, n, d_v.get(), d_scalars.get());
+                    int blocks = (n + BLOCK_SIZE - 1) / BLOCK_SIZE;
+                    detail::fused_normalize_kernel<<<blocks, BLOCK_SIZE, 0, ctx.stream>>>(
+                        d_v.get(), d_v.get(), d_scalars.get(), n);
+                }
             }
         }
 

tests/testthat/test_svd.R

@@ -167,10 +167,10 @@ test_that("svd() with nonneg applies non-negativity constraint (deflation)", {
   nn <- RcppML::svd(A_nn, k = 3, nonneg = TRUE, method = "deflation",
                     seed = 1, maxit = 200, tol = 1e-6, test_fraction = 0)
   expect_s4_class(nn, "svd")
-  # Note: deflation applies Gram-Schmidt reorthogonalization AFTER nonneg
-  # projection, which can reintroduce negatives. This is a known limitation.
-  # For strict non-negativity, use method="krylov" instead.
-  # Here we just check reasonable reconstruction quality.
+  # Constrained factors skip post-factor Gram-Schmidt, so every factor
+  # (not just the first) stays non-negative.
+  expect_true(all(nn@u >= 0), info = "deflation nonneg: U non-negative")
+  expect_true(all(nn@v >= 0), info = "deflation nonneg: V non-negative")
   A_hat <- reconstruct(nn)
   residual_frac <- sum((A_nn - A_hat)^2) / sum(A_nn^2)
   expect_lt(residual_frac, 0.5)
@@ -179,8 +179,8 @@ test_that("svd() with nonneg applies non-negativity constraint (deflation)", {
 test_that("svd() with nonneg via krylov produces non-negative factors", {
   set.seed(77)
   A_nn <- matrix(abs(rnorm(80 * 40)), 80, 40)
-  nn <- RcppML::svd(A_nn, k = 3, nonneg = TRUE, method = "krylov",
-                    seed = 1, maxit = 200, tol = 1e-6, test_fraction = 0)
+  nn <- suppressWarnings(RcppML::svd(A_nn, k = 3, nonneg = TRUE, method = "krylov",
+                    seed = 1, maxit = 200, tol = 1e-6, test_fraction = 0))
   expect_s4_class(nn, "svd")
   expect_true(all(nn@v >= -1e-8), info = "krylov nonneg: V should be non-negative")
   expect_true(all(nn@u >= -1e-8), info = "krylov nonneg: U should be non-negative")
@@ -189,8 +189,8 @@ test_that("svd() with nonneg via krylov produces non-negative factors", {
 test_that("svd() with L1 via krylov applies sparsity", {
   set.seed(99)
   A_sp <- matrix(abs(rnorm(100 * 50)), 100, 50)
-  sp <- RcppML::svd(A_sp, k = 3, L1 = 1.0, method = "krylov",
-                    seed = 1, maxit = 200, tol = 1e-6, test_fraction = 0)
+  sp <- suppressWarnings(RcppML::svd(A_sp, k = 3, L1 = 1.0, method = "krylov",
+                    seed = 1, maxit = 200, tol = 1e-6, test_fraction = 0))
   expect_s4_class(sp, "svd")
   expect_equal(ncol(sp@v), 3L)
   # L1 should induce some zero entries in V
@@ -200,9 +200,9 @@ test_that("svd() with L1 via krylov applies sparsity", {
 test_that("svd() krylov+nonneg+L1 combined constraints", {
   set.seed(42)
   A_combo <- matrix(abs(rnorm(80 * 40)), 80, 40)
-  res <- RcppML::svd(A_combo, k = 3, nonneg = TRUE, L1 = c(0, 0.5),
+  res <- suppressWarnings(RcppML::svd(A_combo, k = 3, nonneg = TRUE, L1 = c(0, 0.5),
                      method = "krylov", seed = 1, maxit = 200,
-                     tol = 1e-6, test_fraction = 0)
+                     tol = 1e-6, test_fraction = 0))
   expect_s4_class(res, "svd")
   expect_true(all(res@u >= -1e-8), info = "krylov nonneg+L1: U non-negative")
   expect_true(all(res@v >= -1e-8), info = "krylov nonneg+L1: V non-negative")
@@ -215,8 +215,8 @@ test_that("krylov vs deflation nonneg quality is comparable", {
   A_nn <- matrix(abs(rnorm(80 * 40)), 80, 40)
   defl <- RcppML::svd(A_nn, k = 3, nonneg = TRUE, method = "deflation",
                       seed = 1, maxit = 200, tol = 1e-6, test_fraction = 0)
-  kry <- RcppML::svd(A_nn, k = 3, nonneg = TRUE, method = "krylov",
-                     seed = 1, maxit = 200, tol = 1e-6, test_fraction = 0)
+  kry <- suppressWarnings(RcppML::svd(A_nn, k = 3, nonneg = TRUE, method = "krylov",
+                     seed = 1, maxit = 200, tol = 1e-6, test_fraction = 0))
   # Both should reconstruct to similar quality
   A_hat_defl <- reconstruct(defl)
   A_hat_kry <- reconstruct(kry)
@@ -227,6 +227,45 @@ test_that("krylov vs deflation nonneg quality is comparable", {
   expect_lt(err_kry, 0.5)
 })
 
+test_that("deflation nonneg PCA keeps every component non-negative (mixed-sign data)", {
+  # Half the features are anti-correlated with the other half, so the
+  # unconstrained PC1 has mixed signs and later factors overlap earlier ones.
+  set.seed(3)
+  f <- rnorm(600)
+  X <- sapply(1:30, function(j) (if (j <= 15) f else -f) * runif(1, .3, 1) + rnorm(600))
+  res <- RcppML::svd(t(X), k = 3, center = TRUE, nonneg = c(TRUE, FALSE),
+                     method = "deflation", seed = 1, test_fraction = 0)
+  expect_true(all(res@u >= 0))
+  expect_true(all(colSums(res@u > 0) > 0), info = "no all-zero components")
+})
+
+test_that("deflation L1 sparsity is preserved on later factors", {
+  set.seed(5)
+  A <- matrix(rnorm(60 * 300), 60, 300)
+  A[1:20, ] <- A[1:20, ] + outer(rnorm(20, 2), rnorm(300))
+  res <- RcppML::svd(A, k = 3, center = TRUE, L1 = c(2, 0),
+                     method = "deflation", seed = 1, test_fraction = 0)
+  expect_true(all(colSums(res@u == 0) > 0),
+              info = "each factor keeps exact zeros from L1 soft-thresholding")
+})
+
+test_that("auto method uses deflation for nonneg/L1 constraints at any k", {
+  set.seed(9)
+  A <- matrix(abs(rnorm(200 * 60)), 200, 60)
+  res <- RcppML::svd(A, k = 8, nonneg = TRUE, seed = 1, test_fraction = 0)
+  expect_equal(res@misc$method, "deflation")
+})
+
+test_that("krylov with nonneg and k > 1 warns", {
+  set.seed(9)
+  A <- matrix(abs(rnorm(100 * 40)), 100, 40)
+  expect_warning(RcppML::svd(A, k = 3, nonneg = TRUE, method = "krylov",
+                             seed = 1, test_fraction = 0),
+                 "redundant or all-zero")
+  expect_no_warning(RcppML::svd(A, k = 1, nonneg = TRUE, method = "krylov",
+                                seed = 1, test_fraction = 0))
+})
+
 # ===========================================================================
 # S4 methods
 # ===========================================================================