From c6d9bda1b66d1add30cda63cc711777d836c030a Mon Sep 17 00:00:00 2001 From: Zachary DeBruine Date: Mon, 19 Oct 2026 13:23:49 -0400 Subject: [PATCH] fix(svd): projection deflation for constrained factors (no repeated factors) With nonneg/L1/upper_bound constraints, deflation subtracted each fitted rank-1 term (Hotelling deflation). Constrained factors are not singular vectors, so the residual kept energy along earlier factors and later factors converged to near-copies of them. On pbmc3k (all 13,714 genes, k = 16, nonneg loadings) the 16 returned factors spanned only 8 distinct directions (max |cos| between factors 0.99). When either side is constrained, use generalized (projection) deflation, as in nsprcomp: R_k = P_u A P_v, with P = I - QQ' over an orthonormal basis of the earlier constrained factors. Both deflated products go through one helper; post-factor Gram-Schmidt is skipped in this mode. Unconstrained SVD/PCA keeps the Hotelling path unchanged. pbmc3k, k = 16: 16 distinct directions (max |cos| 0.43), 0 negative loadings, best-fit relative residual 0.9478 vs 0.9480 for nsprcomp. New test fails on the previous code (max |cos| 1.000) and passes here. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014zFkHdSUAE1RSwGtTqz9Fi --- NEWS.md | 7 ++ inst/include/FactorNet/svd/deflation.hpp | 98 +++++++++++++++++++----- tests/testthat/test_svd.R | 16 ++++ 3 files changed, 102 insertions(+), 19 deletions(-) diff --git a/NEWS.md b/NEWS.md index 703820e..63136b6 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,6 +2,13 @@ ## Bug fixes +- `svd()`/`pca()` with `method = "deflation"` and `nonneg`/`L1`/`upper_bound` + now use projection (generalized) deflation for constrained factors, as in + nsprcomp: each factor is fit to the data with earlier constrained directions + projected out. Hotelling deflation left energy along earlier constrained + factors, so at larger k later factors repeated earlier ones (e.g. pbmc3k, + k = 16: 16 factors spanned only 8 directions). Unconstrained SVD/PCA is + unchanged. CPU path only; the GPU deflation path is not yet updated. - `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 diff --git a/inst/include/FactorNet/svd/deflation.hpp b/inst/include/FactorNet/svd/deflation.hpp index 7ead197..61bc888 100644 --- a/inst/include/FactorNet/svd/deflation.hpp +++ b/inst/include/FactorNet/svd/deflation.hpp @@ -626,6 +626,60 @@ SVDResult deflation_svd(const MatrixType& A, const bool v_constrained = config.nonneg_v || config.L1_v > 0 || config.upper_bound_v > 0; + // ------------------------------------------------------------------ + // Projection deflation for constrained factors. + // + // Hotelling deflation (R_k = A - sum_r sigma_r u_r v_r') only removes a + // direction exactly when u_r, v_r are singular vectors. Constrained factors + // are not, so the residual keeps energy along u_r and later factors can + // re-find (near-)copies of earlier ones. When either side is constrained we + // instead use generalized (projection) deflation, as in nsprcomp: + // R_k = P_u A P_v, P_u = I - Qu Qu', P_v = I - Qv Qv', + // where Qu (Qv) is an orthonormal basis for the earlier constrained u (v) + // factors. Each new factor can then only explain variance not already + // captured. Unconstrained SVD/PCA keeps the Hotelling path unchanged. + // ------------------------------------------------------------------ + const bool use_projection = u_constrained || v_constrained; + DenseMatrix Qu, Qv; + if (u_constrained) { Qu.resize(m, k_max); Qu.setZero(); } + if (v_constrained) { Qv.resize(n, k_max); Qv.setZero(); } + int n_proj = 0; + auto project_out = [&](DenseVector& x, const DenseMatrix& Q) { + if (n_proj == 0) return; + auto Qk = Q.leftCols(n_proj); + for (int pass = 0; pass < 2; ++pass) { + DenseVector c = Qk.transpose() * x; + x.noalias() -= Qk * c; + } + }; + // Deflated products: out = R_k v and out = R_k' u + auto apply_R = [&](const DenseVector& v_in, DenseVector& out, int k_cur) { + if (!use_projection) { + spmv_f(v_in, out); + deflation_correct(out, v_in, U_all, d_all, V_all, k_cur); + return; + } + if (v_constrained && n_proj > 0) { + DenseVector vp = v_in; project_out(vp, Qv); spmv_f(vp, out); + } else { + spmv_f(v_in, out); + } + if (u_constrained) project_out(out, Qu); + }; + auto apply_Rt = [&](const DenseVector& u_in, DenseVector& out, int k_cur) { + if (!use_projection) { + spmv_t(u_in, out); + deflation_correct(out, u_in, V_all, d_all, U_all, k_cur); + return; + } + if (u_constrained && n_proj > 0) { + DenseVector up = u_in; project_out(up, Qu); spmv_t(up, out); + } else { + spmv_t(u_in, out); + } + if (v_constrained) project_out(out, Qv); + }; + for (int k = 0; k < k_max; ++k) { SVD_CHECK_INTERRUPT(); @@ -652,13 +706,11 @@ SVDResult deflation_svd(const MatrixType& A, Scalar u_norm_init = u.norm(); if (u_norm_init > std::numeric_limits::epsilon() * 100) { u /= u_norm_init; - spmv_t(u, spmv_result); - deflation_correct(spmv_result, u, V_all, d_all, U_all, k); + apply_Rt(u, spmv_result, k); v = spmv_result; Scalar v_norm_tmp = v.norm(); if (v_norm_tmp > 0) v /= v_norm_tmp; - spmv_f(v, spmv_result); - deflation_correct(spmv_result, v, U_all, d_all, V_all, k); + apply_R(v, spmv_result, k); u = spmv_result; for (int r = 0; r < k; ++r) { u -= u.dot(U_all.col(r)) * U_all.col(r); @@ -701,14 +753,12 @@ SVDResult deflation_svd(const MatrixType& A, // IRLS: compute row weights from current u, v, sigma // residual_i = (Av)_i - sigma * u_i // ============================================================ - spmv_f(v, spmv_result_unweighted); - deflation_correct(spmv_result_unweighted, v, U_all, d_all, V_all, k); + apply_R(v, spmv_result_unweighted, k); compute_row_weights(spmv_result_unweighted, u, sigma, config.robust_delta, row_weights, eff_threads); // Column weights from A'u residuals - spmv_t(u, spmv_result_unweighted); - deflation_correct(spmv_result_unweighted, u, V_all, d_all, U_all, k); + apply_Rt(u, spmv_result_unweighted, k); compute_col_weights(spmv_result_unweighted, v, sigma, config.robust_delta, col_weights, eff_threads); } @@ -717,8 +767,7 @@ SVDResult deflation_svd(const MatrixType& A, if (do_robust && iter > 0) { // Weighted v-update: v = A' diag(w_u) u_hat / (u_hat' diag(w_u) u_hat) DenseVector wu = u_hat.cwiseProduct(row_weights); - spmv_t(wu, spmv_result); - deflation_correct(spmv_result, wu, V_all, d_all, U_all, k); + apply_Rt(wu, spmv_result, k); Scalar u_sq = wu.dot(u_hat); // u_hat' diag(w) u_hat if (do_cv) u_sq *= cv_denom_correction; if (u_sq > 0) { @@ -729,8 +778,7 @@ SVDResult deflation_svd(const MatrixType& A, } } else { // Standard v-update - spmv_t(u_hat, spmv_result); - deflation_correct(spmv_result, u_hat, V_all, d_all, U_all, k); + apply_Rt(u_hat, spmv_result, k); Scalar u_sq = u_hat.squaredNorm(); if (do_cv) u_sq *= cv_denom_correction; if (u_sq > 0) { @@ -762,8 +810,7 @@ SVDResult deflation_svd(const MatrixType& A, if (do_robust && iter > 0) { // Weighted u-update: u = A diag(w_v) v / (v' diag(w_v) v) DenseVector wv = v.cwiseProduct(col_weights); - spmv_f(wv, spmv_result); - deflation_correct(spmv_result, wv, U_all, d_all, V_all, k); + apply_R(wv, spmv_result, k); Scalar v_sq = wv.dot(v); // v' diag(w) v if (do_cv) v_sq *= cv_denom_correction; if (v_sq > 0) { @@ -774,8 +821,7 @@ SVDResult deflation_svd(const MatrixType& A, } } else { // Standard u-update - spmv_f(v, spmv_result); - deflation_correct(spmv_result, v, U_all, d_all, V_all, k); + apply_R(v, spmv_result, k); Scalar v_sq = v.squaredNorm(); if (do_cv) v_sq *= cv_denom_correction; if (v_sq > 0) { @@ -839,7 +885,7 @@ SVDResult deflation_svd(const MatrixType& A, // 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) { + if (k > 0 && !use_projection) { auto U_prev = U_all.leftCols(k); auto V_prev = V_all.leftCols(k); @@ -872,8 +918,7 @@ SVDResult deflation_svd(const MatrixType& A, // The ALS-loop sigma was computed before GS modified u and v, // so it's stale. Rayleigh quotient: sigma = u' * (A*v corrected). { - spmv_f(v, spmv_result); - deflation_correct(spmv_result, v, U_all, d_all, V_all, k); + apply_R(v, spmv_result, k); sigma = u.dot(spmv_result); if (sigma < 0) sigma = -sigma; // singular values are non-negative } @@ -883,6 +928,21 @@ SVDResult deflation_svd(const MatrixType& A, d_all(k) = sigma; V_all.col(k) = v; + if (use_projection) { + const Scalar eps_q = std::numeric_limits::epsilon() * 100; + if (u_constrained) { + DenseVector q = u; project_out(q, Qu); + Scalar qn = q.norm(); + if (qn > eps_q) Qu.col(n_proj) = q / qn; else Qu.col(n_proj).setZero(); + } + if (v_constrained) { + DenseVector q = v; project_out(q, Qv); + Scalar qn = q.norm(); + if (qn > eps_q) Qv.col(n_proj) = q / qn; else Qv.col(n_proj).setZero(); + } + ++n_proj; + } + result.iters_per_factor.push_back(iter); // === Test loss evaluation === diff --git a/tests/testthat/test_svd.R b/tests/testthat/test_svd.R index eb24d1b..a7106a1 100644 --- a/tests/testthat/test_svd.R +++ b/tests/testthat/test_svd.R @@ -249,6 +249,22 @@ test_that("deflation L1 sparsity is preserved on later factors", { info = "each factor keeps exact zeros from L1 soft-thresholding") }) +test_that("deflation nonneg PCA returns distinct factors at larger k", { + # Hotelling deflation let constrained factors re-find near-copies of earlier + # ones; projection deflation must yield k distinct non-negative directions. + set.seed(11) + A <- matrix(rexp(120 * 400), 120, 400) + A[1:30, 1:200] <- A[1:30, 1:200] + 2 + A[31:60, 150:350] <- A[31:60, 150:350] + 2 + res <- RcppML::svd(A, k = 10, center = TRUE, nonneg = c(TRUE, FALSE), + method = "deflation", seed = 1, test_fraction = 0) + W <- sweep(res@u, 2, sqrt(colSums(res@u^2)), "/") + C <- abs(crossprod(W)); diag(C) <- 0 + expect_true(all(res@u >= 0)) + expect_lt(max(C), 0.95) + expect_equal(qr(W, tol = 1e-3)$rank, 10L) +}) + test_that("auto method uses deflation for nonneg/L1 constraints at any k", { set.seed(9) A <- matrix(abs(rnorm(200 * 60)), 200, 60)