diff --git a/NEWS.md b/NEWS.md index 885058b..703820e 100644 --- a/NEWS.md +++ b/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 diff --git a/R/svd.R b/R/svd.R index 8ce0968..7b373b9 100644 --- a/R/svd.R +++ b/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) { diff --git a/inst/include/FactorNet/svd/deflation.hpp b/inst/include/FactorNet/svd/deflation.hpp index 1f15198..7ead197 100644 --- a/inst/include/FactorNet/svd/deflation.hpp +++ b/inst/include/FactorNet/svd/deflation.hpp @@ -619,6 +619,13 @@ SVDResult 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 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(rng.uniform()); @@ -824,28 +834,37 @@ SVDResult 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 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::epsilon() * 100) { - u /= u_norm_gs; + if (!u_constrained) { + for (int pass = 0; pass < 2; ++pass) { + DenseVector 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::epsilon() * 100) { + u /= u_norm_gs; + } } // Reorthogonalize v against V_all[:, 0:k-1] (two passes) - for (int pass = 0; pass < 2; ++pass) { - DenseVector 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::epsilon() * 100) { - v /= v_norm_gs; + if (!v_constrained) { + for (int pass = 0; pass < 2; ++pass) { + DenseVector 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::epsilon() * 100) { + v /= v_norm_gs; + } } } diff --git a/inst/include/FactorNet/svd/deflation_gpu.cuh b/inst/include/FactorNet/svd/deflation_gpu.cuh index 74e3df6..d87b36c 100644 --- a/inst/include/FactorNet/svd/deflation_gpu.cuh +++ b/inst/include/FactorNet/svd/deflation_gpu.cuh @@ -636,6 +636,12 @@ SVDResult deflation_svd_gpu( int m, int n, int nnz, const SVDConfig& 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; using DenseVectorS = Eigen::Matrix; @@ -1253,27 +1259,31 @@ SVDResult 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( - ctx.cublas, d_U_factors.get(), d_u.get(), m, k, d_gs_h); - // Renormalize d_u - { - detail::cublas_nrm2_device( - ctx.cublas, m, d_u.get(), d_scalars.get()); - int blocks = (m + BLOCK_SIZE - 1) / BLOCK_SIZE; - detail::fused_normalize_kernel<<>>( - d_u.get(), d_u.get(), d_scalars.get(), m); + if (!u_constrained) { + detail::gpu_orthog_against( + ctx.cublas, d_U_factors.get(), d_u.get(), m, k, d_gs_h); + // Renormalize d_u + { + detail::cublas_nrm2_device( + ctx.cublas, m, d_u.get(), d_scalars.get()); + int blocks = (m + BLOCK_SIZE - 1) / BLOCK_SIZE; + detail::fused_normalize_kernel<<>>( + d_u.get(), d_u.get(), d_scalars.get(), m); + } } // Reorthogonalize d_v against d_V_factors[:, 0:k-1] - detail::gpu_orthog_against( - ctx.cublas, d_V_factors.get(), d_v.get(), n, k, d_gs_h); - // Renormalize d_v - { - detail::cublas_nrm2_device( - ctx.cublas, n, d_v.get(), d_scalars.get()); - int blocks = (n + BLOCK_SIZE - 1) / BLOCK_SIZE; - detail::fused_normalize_kernel<<>>( - d_v.get(), d_v.get(), d_scalars.get(), n); + if (!v_constrained) { + detail::gpu_orthog_against( + ctx.cublas, d_V_factors.get(), d_v.get(), n, k, d_gs_h); + // Renormalize d_v + { + detail::cublas_nrm2_device( + ctx.cublas, n, d_v.get(), d_scalars.get()); + int blocks = (n + BLOCK_SIZE - 1) / BLOCK_SIZE; + detail::fused_normalize_kernel<<>>( + d_v.get(), d_v.get(), d_scalars.get(), n); + } } } @@ -1432,6 +1442,12 @@ SVDResult deflation_svd_gpu_dense( int m, int n, const SVDConfig& 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; using DenseVectorS = Eigen::Matrix; @@ -1797,23 +1813,27 @@ SVDResult 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( - ctx.cublas, d_U_factors.get(), d_u.get(), m, k, d_gs_h); - { - detail::cublas_nrm2_device( - ctx.cublas, m, d_u.get(), d_scalars.get()); - int blocks = (m + BLOCK_SIZE - 1) / BLOCK_SIZE; - detail::fused_normalize_kernel<<>>( - d_u.get(), d_u.get(), d_scalars.get(), m); + if (!u_constrained) { + detail::gpu_orthog_against( + ctx.cublas, d_U_factors.get(), d_u.get(), m, k, d_gs_h); + { + detail::cublas_nrm2_device( + ctx.cublas, m, d_u.get(), d_scalars.get()); + int blocks = (m + BLOCK_SIZE - 1) / BLOCK_SIZE; + detail::fused_normalize_kernel<<>>( + d_u.get(), d_u.get(), d_scalars.get(), m); + } } - detail::gpu_orthog_against( - ctx.cublas, d_V_factors.get(), d_v.get(), n, k, d_gs_h); - { - detail::cublas_nrm2_device( - ctx.cublas, n, d_v.get(), d_scalars.get()); - int blocks = (n + BLOCK_SIZE - 1) / BLOCK_SIZE; - detail::fused_normalize_kernel<<>>( - d_v.get(), d_v.get(), d_scalars.get(), n); + if (!v_constrained) { + detail::gpu_orthog_against( + ctx.cublas, d_V_factors.get(), d_v.get(), n, k, d_gs_h); + { + detail::cublas_nrm2_device( + ctx.cublas, n, d_v.get(), d_scalars.get()); + int blocks = (n + BLOCK_SIZE - 1) / BLOCK_SIZE; + detail::fused_normalize_kernel<<>>( + d_v.get(), d_v.get(), d_scalars.get(), n); + } } } diff --git a/tests/testthat/test_svd.R b/tests/testthat/test_svd.R index 8132272..eb24d1b 100644 --- a/tests/testthat/test_svd.R +++ b/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 # ===========================================================================