Skip to content

fix(svd): projection deflation for constrained factors (no repeated factors) #84

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
zdebruine merged 1 commit into from
Oct 19, 2026
Merged

fix(svd): projection deflation for constrained factors (no repeated factors) #84

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions NEWS.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
98 changes: 79 additions & 19 deletions inst/include/FactorNet/svd/deflation.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -626,6 +626,60 @@ SVDResult<Scalar> 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<Scalar> 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<Scalar>& x, const DenseMatrix<Scalar>& Q) {
if (n_proj == 0) return;
auto Qk = Q.leftCols(n_proj);
for (int pass = 0; pass < 2; ++pass) {
DenseVector<Scalar> c = Qk.transpose() * x;
x.noalias() -= Qk * c;
}
};
// Deflated products: out = R_k v and out = R_k' u
auto apply_R = [&](const DenseVector<Scalar>& v_in, DenseVector<Scalar>& 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<Scalar> 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<Scalar>& u_in, DenseVector<Scalar>& 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<Scalar> 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();

Expand All @@ -652,13 +706,11 @@ SVDResult<Scalar> deflation_svd(const MatrixType& A,
Scalar u_norm_init = u.norm();
if (u_norm_init > std::numeric_limits<Scalar>::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);
Expand Down Expand Up @@ -701,14 +753,12 @@ SVDResult<Scalar> 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);
}
Expand All @@ -717,8 +767,7 @@ SVDResult<Scalar> 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<Scalar> 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) {
Expand All @@ -729,8 +778,7 @@ SVDResult<Scalar> 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) {
Expand Down Expand Up @@ -762,8 +810,7 @@ SVDResult<Scalar> 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<Scalar> 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) {
Expand All @@ -774,8 +821,7 @@ SVDResult<Scalar> 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) {
Expand Down Expand Up @@ -839,7 +885,7 @@ SVDResult<Scalar> 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);

Expand Down Expand Up @@ -872,8 +918,7 @@ SVDResult<Scalar> 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
}
Expand All @@ -883,6 +928,21 @@ SVDResult<Scalar> deflation_svd(const MatrixType& A,
d_all(k) = sigma;
V_all.col(k) = v;

if (use_projection) {
const Scalar eps_q = std::numeric_limits<Scalar>::epsilon() * 100;
if (u_constrained) {
DenseVector<Scalar> 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<Scalar> 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 ===
Expand Down
16 changes: 16 additions & 0 deletions tests/testthat/test_svd.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Loading