feat(strategy): rsvd captures 3 subspaces, not 4 (drop redundant col(B) sketch)

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():

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.