SVD via Golub-Kahan-Reinsch with Givens rotations - 100% convergence, PCA-ready

[cyancirrus]

SVD (Golub–Kahan–Reinsch)

This adds a full SVD implementation using Golub–Kahan bidiagonalization followed by implicit-shift bulge chasing (the SVD analogue of Francis QR) — same family of algorithm as the eigenvalue PR, applied to the non-symmetric/rectangular case.

A few things worth calling out before you read the diff:

• As shown in the earlier eigenvalue PR, the singular values of X and the eigenvalues of the covariance matrix X'X are equivalent (derivation below, under "Why this replaces the eigen-approach"). This means you can get PCA directly from this SVD without ever forming X'X — which matters, because forming X'X squares the condition number of your data and is a real source of numerical error, especially on ill-conditioned inputs.
• This has 100% convergence in testing — 10,000 trials with zero orthogonality or reconstruction failures (see test_svd_reconstruct_trials, tightened to assert == 0 rather than a tolerance band, so the test actually backs up the claim). Iteration budget is 20 per eigenvalue plus a 20-iteration bank, and across everything I've run it has not diverged once. I think that's a structurally sound bound, not just an empirical accident, but flagging it as something worth stress-testing further on your end too if you want extra confidence before leaning on it hard.
• Structure mirrors the eigenvalue PR pretty closely (same full_*/bare split, same rotation-tracking pattern, same zero-alloc core with an allocating auto_* wrapper on top), so hopefully the shape of this is already familiar.

Why this replaces the eigen-approach for your PCA case

The eigen-solver would work here too — this isn't "that approach was wrong," it's that this one is better suited to what you're doing, for two concrete reasons:

1. Condition number. X'X has condition number κ(X)². SVD works on X directly, so you're never squaring your conditioning — meaningfully better numerical stability, especially as your data gets less well-conditioned.
2. No rotational trap. Because SVD applies independent left and right orthogonal transforms (U and V) rather than a single similarity transform, it doesn't get stuck in the kind of degenerate rotational configurations that can slow or destabilize a symmetric eigensolver on nearly-degenerate inputs.

The math connecting the two (so this isn't just assertion):

// center the data
X <- X - mean(X);

// Method 1: Schur decomposition of X'X
X'X = Y;
Schur(Y) => Y = VTV';

// Y is symmetric (since Y = X'X)
=> T is diagonal here
Y = VDV';

// Method 2: SVD of X directly
SVD(X) => X = USV';
X'X = (USV')'(USV');
X'X = VS²V';

// Since Y = X'X, comparing the two results:
VDV' = VS²V';
=> D = S²;

V ~ eigenvectors of X'X = principal components (columns of V)

So V from this SVD gives you your principal components directly, and S² / (n-1) gives you the variances — no eigendecomposition of X'X needed at all.

What you'll need to do on your end

This gets you 95% of the way to PCA, but there are three small pieces I've deliberately left for you to wire in rather than handing you a finished pca() — partly because they're genuinely trivial, and partly because I think it's worth actually touching this code rather than just importing a black box, especially after ~1700 lines went into the eigenvalue side of this. All three are short:

1. Center the incoming matrix by its column means before calling into the SVD. This is required, not optional — SVD has no notion of "variance," it just factors whatever matrix you hand it, so skipping this gets you components that point at your data's mean rather than its spread.
2. Sorting. The singular values come out very nearly sorted already, but not guaranteed-strict — two values can be close enough that order isn't guaranteed. If you want strict descending order, write a sort over the resulting diagonal and apply the same permutation as a column permutation to U and V (both are just column permutations, not a re-derivation of anything).
3. Sign. Singular values are inherently sign-ambiguous — (u_i, v_i) and (-u_i, -v_i) are both valid. If you care about consistent sign conventions (e.g. reproducibility across runs), negate the corresponding column of U or V to match; if you only care about magnitudes, you can ignore this entirely.

All three of the above are direct consequences of what's already computed — no new math, just plumbing.

Status

This port was done a bit quickly and should be up to spec, but it's probably due for a more in-depth review pass before you fully rely on it — flagging that honestly rather than presenting it as more battle-tested than it is.

Main branch highlights

• Golub–Kahan bidiagonalization
• Reinsch-style implicit shift, driven off the singular values
• The incoming matrix gets routed to either the upper- or lower-bidiagonal path depending on whether rows or cols is larger — both paths are needed (not just one generalized path) because of how the trailing value and iteration order work; there's a top-level svd/full_svd_decomposition that dispatches into whichever is correct for the input shape.
• Zero data movement, zero allocation in the core path. For very large matrices where cols >> rows (or vice versa) in a way that would blow out L1 cache, transposing the input first is worth doing — should be quick even at that scale, but it's a deliberate option, not something forced on you.

---

One last thing, unrelated to the code: sorry for how the earlier back-and-forth about the eigen-solver landed — I know it probably read as pushback for its own sake at the time. This is the reasoning behind it, and hopefully it's useful now that it's concrete instead of abstract.

[cyancirrus]

Performance vs. nalgebra

Benchmarked Autumn_SVD against nalgebra's SVD across sizes 8–64.

Size	Autumn_SVD	Nalgebra_SVD	Full_Autumn_SVD	Full_Nalgebra_SVD
8	4.77 µs	4.89 µs	6.69 µs	6.48 µs
16	18.77 µs	15.39 µs	34.69 µs	22.60 µs
32	77.76 µs	56.96 µs	187.72 µs	101.37 µs
64	364.34 µs	247.17 µs	1237.5 µs	567.24 µs

At size 8 the two implementations are essentially at parity. From 16 onward, nalgebra pulls ahead, and the gap widens with size — most noticeably in the Full_* variants (full U/Σ/V vs. values-only), where Full_Autumn_SVD is ~2.2x slower than Full_Nalgebra_SVD at size 64.

Available optimizations

(out of scope for this PR)

• Block the U/V rotations and batch them instead of applying one at a time
• Cache U'/V' to avoid recomputing transforms
• Pack the diagonal entries densely for better cache locality

Autumn SVD is this implementation