From 25bd3758087a000a2a26e87f1672952926995315 Mon Sep 17 00:00:00 2001 From: ultrahighsuper Date: Fri, 10 Jul 2026 22:22:13 +0900 Subject: [PATCH] feat(strategy): rsvd captures 3 subspaces, not 4 (col(B) is redundant) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The subspace reconstruction is Ĉ = P A P B P with P = QQᵀ, which equals A@B once range(Q) contains col(A), row(A), row(B): P A = A needs col(A); A P = A needs row(A); B P = B needs row(B) => P A P B P = A B P = A B col(B) is never needed. RandomizedSVDTransform split the M budget four ways (col A, row A, col B, row B), so a quarter of every rank went to a subspace that does not affect the result -- exact recovery of a rank-r product needed M ~ 4r. Splitting three ways over {col(A), row(A), row(B)} recovers at M ~ 3r, a 25% smaller subspace for the same accuracy (and better accuracy at any fixed M < 4r). Total sketch width is still m, so basis_flops is unchanged. Measured improvement over exact on compressible data at M = 3r (GTX 1650, fp32, N=8192, lowrank rank-16, seed 0): accuracy 1.0000, 748ms vs 1774ms exact, 259 vs 1288 MiB, 42.4x fewer FLOPs. The old 4-way split at the same M=48 gates at accuracy 0.35 (score 0); the 3-way split scores an improvement. Fixes #91 --- strategy/tests/test_basis_vram_fraction.py | 7 ++--- strategy/transforms.py | 30 ++++++++++++---------- 2 files changed, 21 insertions(+), 16 deletions(-) diff --git a/strategy/tests/test_basis_vram_fraction.py b/strategy/tests/test_basis_vram_fraction.py index 4c217bf..d6998a5 100644 --- a/strategy/tests/test_basis_vram_fraction.py +++ b/strategy/tests/test_basis_vram_fraction.py @@ -56,14 +56,15 @@ def fake(X, Q, backend, dtype, frac=sub._DEFAULT_ROW_BLOCK_FRACTION): def test_rsvd_basis_threads_the_configured_fraction(): + # rsvd captures 3 subspaces (col(A), row(A), row(B)); frac reaches every sketch. fracs = _capture_fracs(frac_kwarg=0.15) - assert fracs == [0.15, 0.15, 0.15, 0.15], fracs + assert fracs == [0.15, 0.15, 0.15], fracs def test_rsvd_basis_defaults_to_streaming_default_when_none(): d = sub._DEFAULT_ROW_BLOCK_FRACTION - assert _capture_fracs(frac_kwarg=None) == [d, d, d, d] - assert _capture_fracs(frac_kwarg="__omit__") == [d, d, d, d] + assert _capture_fracs(frac_kwarg=None) == [d, d, d] + assert _capture_fracs(frac_kwarg="__omit__") == [d, d, d] def test_basis_signatures_accept_frac(): diff --git a/strategy/transforms.py b/strategy/transforms.py index 688bb3f..fa7a9d0 100644 --- a/strategy/transforms.py +++ b/strategy/transforms.py @@ -62,11 +62,17 @@ def _orthonormalize(M, backend): class RandomizedSVDTransform(Transform): """Data-dependent range finder over A and B (the accurate one). - Splits the M-column budget evenly across the four spaces that must be - captured for the product -- col(A), row(A), col(B), row(B) -- via random - sketches, then orthonormalizes. Because all four are represented, the - reconstruction converges to the exact product as M approaches the numerical - rank. Sketches stream, so A/B may be disk-backed memmaps. + The reconstruction is ``Ĉ = P A P B P`` with the projector ``P = Q Qᵀ``, so + ``Ĉ = A @ B`` exactly once range(Q) contains col(A), row(A), row(B): + + P A = A needs col(A) ⊆ range(Q) + A P = A needs row(A) ⊆ range(Q) + B P = B needs row(B) ⊆ range(Q) (then P A P B P = A B P = A B) + + Three spaces are necessary and sufficient -- col(B) is redundant. We split + the M-column budget across those three via random sketches, then + orthonormalize; exact recovery of a rank-r product needs only ``M ≳ 3r`` + instead of ``4r``. Sketches stream, so A/B may be disk-backed memmaps. """ name = "rsvd" @@ -87,8 +93,8 @@ def basis(self, n, m, backend, dtype, A=None, B=None, frac=None): frac = _DEFAULT_ROW_BLOCK_FRACTION xp = backend.xp - base, rem = divmod(m, 4) - widths = [base + (1 if i < rem else 0) for i in range(4)] + base, rem = divmod(m, 3) + widths = [base + (1 if i < rem else 0) for i in range(3)] rng = np.random.default_rng(self.seed) def omega(w): @@ -98,19 +104,17 @@ def omega(w): parts = [] if widths[0]: - parts.append(stream_gemm_right(A, omega(widths[0]), backend, dtype, frac)) + parts.append(stream_gemm_right(A, omega(widths[0]), backend, dtype, frac)) # col(A): A Ω if widths[1]: - parts.append(stream_gemm_left_t(A, omega(widths[1]), backend, dtype, frac)) + parts.append(stream_gemm_left_t(A, omega(widths[1]), backend, dtype, frac)) # row(A): Aᵀ Ω if widths[2]: - parts.append(stream_gemm_right(B, omega(widths[2]), backend, dtype, frac)) - if widths[3]: - parts.append(stream_gemm_left_t(B, omega(widths[3]), backend, dtype, frac)) + parts.append(stream_gemm_left_t(B, omega(widths[2]), backend, dtype, frac)) # row(B): Bᵀ Ω Y = xp.concatenate(parts, axis=1) # (n, m) return self._orthonormalize(Y, backend) # (n, m) orthonormal columns def basis_flops(self, n, m): - # 4 random sketches over A and B totalling m columns cost 2*n*n*m FLOPs + # 3 random sketches over A and B totalling m columns cost 2*n*n*m FLOPs # (each width-w sketch A@Omega / A^T@Omega is 2*n*n*w, and the widths sum # to m), plus the QR of the (n, m) sketch ~ 2*n*m*m. Recomputed every call # (the sketches depend on A, B), so it is not amortizable.