linalg: add matrix_rank and cond

[axiom-of-choice]

Closes #125.

What

Two new functions in aneforge/linalg.py, composed from the existing on-ANE SVD:

• matrix_rank(A, tol=None) — numerical rank by counting singular values above tol. Default follows numpy: max(m, n) * eps(fp32) * sigma_max.
• cond(A) — 2-norm condition number sigma_max / sigma_min. Returns inf for singular matrices. Works for rectangular A.

Checks

• ruff check — clean
• pylint 2-space — 10.00/10
• pyright — 0 errors
• compileall — clean
• pytest -m "not requires_ane" — all pass (off-device)
• On-device: 11 new tests pass on M2 Pro / macOS 26.5.2

Notes

• fp16 SVD floor: cond returns inf when condition number exceeds ~1e2 (sigma_min underflows to 0 in fp16). Tests use cond ≤ 5 for oracle comparisons.
• matrix_rank default tolerance follows numpy convention.

[sbryngelson]

This one needs a change before merge: test_matrix_rank_deficient fails on device. CI is green
because the ANE tests are all skipped there, so the checkmarks don't cover this.

E   assert 5 == 1
E    +  where 5 = matrix_rank(array([[ 0.01505, 0.08887, ...]], dtype=float16))

The cause is that both functions read singular values from svd() (Gram matrix + fp16 cyclic-Jacobi
eigh), which is weakest exactly where rank and conditioning live -- in the smallest singular values.
On input that is exactly rank-deficient in fp16 (duplicated rows, so the true rank is unambiguous at
every precision):

n=6, true rank 2    numpy fp64  [1.0e+00 7.8e-01 8.2e-17 1.9e-17 1.1e-33 2.6e-34]
                    svd()       [1.0e+00 7.8e-01 1.6e-02 4.2e-03 3.1e-03 0.0e+00]
n=6, true rank 3    numpy fp64  [1.0e+00 6.4e-01 2.8e-01 6.8e-17 1.7e-17 4.2e-18]
                    svd()       [1.0e+00 6.4e-01 2.5e-01 1.1e-01 5.7e-02 0.0e+00]

In that second case the null space sits at 1.1e-1 and 5.7e-2 of sigma_max while the smallest real
singular value is 2.8e-1. No tolerance separates them, so this isn't fixable by changing the default
tol. More Jacobi sweeps make it worse, not better (n=6 rank-1, spurious/sigma_max: 2.4e-2 at 8
sweeps, 2.1e-1 at 16, 1.0e-1 at 24) -- the fp16 iteration wanders in its own noise rather than
converging.

The same weakness makes cond silently wrong, which I think is the more serious half:

target    numpy fp64    cond via svd()    cond via randomized_svd()
1e+02          99.79              2.03                        99.81
1e+03         992.49               inf                       992.49

Returning 2.03 for a matrix whose condition number is 99.79 is a 50x underestimate with no signal to
the caller. The current test_cond_matches_numpy works around this by restricting the oracle
comparison to cond 2.0 and 5.0, with a comment explaining the floor -- worth treating that as the
bug report it is rather than the documented limit.

The fix is to source both from randomized_svd instead. Its matmuls stay on the ANE and only the QR
and the small SVD go to the host in float64, which is the split this module already documents in its
header, and it puts the null space at ~1e-5 instead of ~1e-1. It's also 13x faster here: the
matrix_rank/cond tests go from 145s to 11s.

Separately, test_matrix_rank_deficient can't pass as written even with a perfect SVD. Rounding
np.outer(v, v) to fp16 makes it genuinely full rank -- np.linalg.matrix_rank on that exact array
returns 8, not 1. I replaced it with duplicated-row input, which is exactly rank-deficient at any
precision, and added an ill-conditioned full-rank case so the default tolerance can't drift into
eating sigma_min.

With this patch the matrix_rank/cond tests go 11 passed + 1 failed -> 18 passed, and the full linalg
suite is 104 passed on device. Reverting just the source change makes the 5 new tests fail, so they
do pin the behaviour.

diff --git a/aneforge/linalg.py b/aneforge/linalg.py
index 28e81e6..28a5eb3 100644
--- a/aneforge/linalg.py
+++ b/aneforge/linalg.py
@@ -774,30 +774,49 @@ def matrix_power(A, n: int):
   return np.asarray(acc, np.float32)
 
 
+# Rank and conditioning both hinge on the SMALLEST singular values, which is exactly where the
+# pure-ANE svd() (Gram matrix + fp16 cyclic-Jacobi eigh) is weakest: on an exactly rank-deficient
+# matrix it leaves the null space at 1e-2..1e-1 of sigma_max, indistinguishable from a real
+# singular value, and no tolerance or sweep count separates them (more sweeps make it worse).
+# The sketch path keeps its matmuls on the ANE but does the QR and the small SVD on the host in
+# float64, which puts that null space at ~1e-5 instead -- a usable gap.  Measured, n=6..8:
+#   exactly rank-3 input   svd() trailing 1.1e-1 / 5.7e-2      randomized_svd() trailing 1.2e-5
+#   cond(A) = 99.8         svd() says 2.03 (silently wrong)    randomized_svd() says 99.81
+_RANK_FLOOR = 1e-5  # measured trailing-sv floor of the sketch path on fp16 input; NOT a machine eps
+
+
+def _svals(A16):
+  """Singular values of A (descending) via the sketch path: ANE matmuls, host float64 QR/SVD."""
+  return randomized_svd(A16, k=min(A16.shape), oversample=5, power_iters=2)[1]
+
+
 def matrix_rank(A, tol=None):
   """Numerical rank of A by counting singular values above `tol`, on the ANE.
 
-  Default `tol` follows numpy: max(m, n) * float32_eps * sigma_max.  Returns 0 for
-  a zero matrix (no positive singular values).  Oracle: np.linalg.matrix_rank."""
+  Default `tol` is max(m, n) * 1e-5 * sigma_max, mirroring numpy's max(m, n) * eps * sigma_max but
+  with the measured floor of this backend in place of a machine epsilon -- the input is fp16, so
+  nothing here resolves a singular value below ~1e-5 of sigma_max.  Pass `tol` explicitly if you
+  know your spectrum.  Returns 0 for a zero matrix.  Oracle: np.linalg.matrix_rank."""
   A16 = np.asarray(A, f16)
   if A16.ndim != 2: raise ValueError(f"matrix_rank: expected 2-D; got shape {A16.shape}")
   m, n = A16.shape
-  S = svd(A16)                                        # descending, float32
+  S = _svals(A16)
   if S.size == 0: return 0
   if tol is None:
-    tol = max(m, n) * np.finfo(np.float32).eps * float(S[0])
+    tol = max(m, n) * _RANK_FLOOR * float(S[0])
   return int(np.sum(S > tol))
 
 
 def cond(A):
   """2-norm condition number sigma_max / sigma_min via SVD, on the ANE.
 
-  Works for square and rectangular A.  Returns inf for a zero matrix.
+  Works for square and rectangular A.  Accurate to cond ~1e3, past which the fp16 input itself
+  stops resolving sigma_min.  Returns inf for a zero or exactly singular matrix.
   Oracle: np.linalg.cond(A)."""
   A16 = np.asarray(A, f16)
   if A16.ndim != 2: raise ValueError(f"cond: expected 2-D; got shape {A16.shape}")
-  S = svd(A16)
-  if S.size == 0 or float(S[-1]) == 0.0: return float("inf")
+  S = _svals(A16)
+  if S.size == 0 or float(S[-1]) <= 0.0: return float("inf")
   return float(S[0]) / float(S[-1])
 
 
diff --git a/tests/test_linalg.py b/tests/test_linalg.py
index e370071..50ff71c 100644
--- a/tests/test_linalg.py
+++ b/tests/test_linalg.py
@@ -484,18 +484,28 @@ def test_matrix_rank_full_rank():
   assert L.matrix_rank(A) == 8
 
 
-def test_matrix_rank_deficient():
-  # rank-1 matrix (outer product) -> clear gap in singular values
-  r = np.random.default_rng(91)
-  v = r.standard_normal(8).astype(f16)
-  A = np.outer(v, v)  # rank 1, sigma_1 >> sigma_2..8 = 0 (in fp16 noise)
-  assert L.matrix_rank(A) == 1
+@pytest.mark.parametrize("n,rank", [(6, 2), (6, 3), (8, 3)])
+def test_matrix_rank_deficient(n, rank):
+  """Duplicated rows are exactly rank-deficient in fp16, so the true rank is unambiguous.
+
+  np.outer(v, v) is NOT: rounding the outer product to fp16 makes it genuinely full rank, and
+  np.linalg.matrix_rank on the same array agrees (8, not 1)."""
+  base = np.random.default_rng(7).standard_normal((rank, n)).astype(f16)
+  A = np.asarray(np.vstack([base] + [base[i % rank] for i in range(n - rank)]), f16)
+  assert np.linalg.matrix_rank(A.astype(np.float64)) == rank, "test input is not exactly rank-deficient"
+  assert L.matrix_rank(A) == rank
 
 
 def test_matrix_rank_zero():
   assert L.matrix_rank(np.zeros((5, 5), f16)) == 0
 
 
+def test_matrix_rank_full_rank_ill_conditioned():
+  """cond=1e3 is full rank, not rank-deficient: the default tol must not eat sigma_min."""
+  A = f16(_square(8, 1e3, 96))
+  assert L.matrix_rank(A) == 8
+
+
 def test_matrix_rank_matches_numpy():
   A = f16(_square(8, 1e1, 92))
   assert L.matrix_rank(A) == np.linalg.matrix_rank(A.astype(np.float64))
@@ -521,14 +531,13 @@ def test_cond_well_conditioned():
   assert abs(c - ref) / ref <= 5e-2
 
 
-def test_cond_matches_numpy():
-  # cond~1e2 is at the fp16 SVD floor: sigma_min can round to 0, giving inf.
-  # Use cond~5 (well within fp16 range) for the oracle comparison.
-  for cond_target in (2.0, 5.0):
-    A = f16(_square(8, cond_target, int(cond_target * 10) + 95))
-    c = L.cond(A)
-    ref = np.linalg.cond(A.astype(np.float64))
-    assert abs(c - ref) / ref <= 5e-2, f"cond={c} vs numpy={ref}"
+@pytest.mark.parametrize("cond_target", [2.0, 5.0, 1e1, 1e2, 1e3])
+def test_cond_matches_numpy(cond_target):
+  """Through cond~1e3. The Gram+Jacobi svd() path silently returned 2.03 here for a true 99.79."""
+  A = f16(_square(8, cond_target, int(cond_target * 10) + 95))
+  c = L.cond(A)
+  ref = np.linalg.cond(A.astype(np.float64))
+  assert abs(c - ref) / ref <= 5e-2, f"cond={c} vs numpy={ref}"
 
 
 def test_cond_singular():

Scheduling note: this, #259 and #262 all append to the same __all__ block and the same tail of
linalg.py, so whichever lands first will force a rebase on the other two.