linalg: add polar decomposition

[axiom-of-choice]

Closes #168.

What

New polar(A) function in aneforge/linalg.py, composed from the existing on-ANE randomized SVD.

Returns (U, P) where A = U @ P, U has orthonormal columns, and P is symmetric positive-semidefinite. Oracle: scipy.linalg.polar(A, side='right').

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: 8 new tests pass on M2 Pro / macOS 26.5.2

Notes

• Uses randomized SVD (not full SVD), so P has ~1e-4 asymmetry from the approximation. Symmetrized in the implementation.
• Works for square and rectangular (tall) matrices.

[sbryngelson]

Verified on device (M5): the 8 polar tests pass and the full linalg suite stays green (97 passed).
I also ran the rectangular shapes the tests don't cover, and they work -- A[8,5] and A[5,8] both
reconstruct to ~7e-4 -- so this is correct as written. Three suggestions, none of them bugs.

The main() demo feeds polar _make_spd(...), and for symmetric positive-definite input polar is
the trivial case: U = I and P = A, so the demo never exercises the interesting part. Swapping in a
general matrix takes U[0,0] from 1.0 to 0.799.

The docstring's "orthogonal (or semi-orthogonal for wide A)" isn't quite right -- U is
semi-orthogonal for either rectangular shape, U^T U = I for tall A and U U^T = I for wide A.

And since rectangular does work, it's worth a test, so the docstring's claim is covered.

As a patch:

diff --git a/aneforge/linalg.py b/aneforge/linalg.py
index 72e113d..3c01051 100644
--- a/aneforge/linalg.py
+++ b/aneforge/linalg.py
@@ -534,6 +534,13 @@ __all__ = [
 
 # __main__ - self-test / validation vs numpy/scipy
 
+def _general_square(n, cond, seed):
+  """A general (non-symmetric) n x n matrix with the target condition number: U diag(s) V^T."""
+  rng = np.random.default_rng(seed)
+  U = np.linalg.qr(rng.standard_normal((n, n)))[0]; V = np.linalg.qr(rng.standard_normal((n, n)))[0]
+  return (U * np.geomspace(1.0, cond, n)) @ V.T
+
+
 def _make_spd(n, cond, seed):
   """SPD A with target condition number, fp16-stored (reference solves the fp16-rounded system)."""
   rng = np.random.default_rng(seed)
@@ -777,7 +784,8 @@ def matrix_power(A, n: int):
 def polar(A):
   """Polar decomposition A = U @ P, on the ANE.
 
-  U is orthogonal (or semi-orthogonal for wide A), P is symmetric positive-semidefinite.
+  U is orthogonal for square A and semi-orthogonal for either rectangular shape (U^T U = I for
+  tall A, U U^T = I for wide A); P is symmetric positive-semidefinite [n,n].
   Composed from the on-ANE SVD: U_ S V^T -> U = U_ V^T, P = V diag(S) V^T.
   Oracle: scipy.linalg.polar(A, side='right')."""
   A16 = np.asarray(A, f16)
@@ -911,7 +919,8 @@ def main():
   print("-" * 90)
   import scipy.linalg
   for n in (6, 8):
-    A_pol = _make_spd(n, 1e1, 80 + n)[0]
+    # a GENERAL matrix: for symmetric positive-definite input polar is the trivial case (U = I, P = A)
+    A_pol = f16(_general_square(n, 1e1, 80 + n))
     U, P = polar(A_pol)
     ref_P = scipy.linalg.polar(A_pol, side="right")[1]
     recon_err = _relerr(U @ P, A_pol)
diff --git a/tests/test_linalg.py b/tests/test_linalg.py
index a222169..661f58d 100644
--- a/tests/test_linalg.py
+++ b/tests/test_linalg.py
@@ -512,6 +512,17 @@ def test_polar_matches_scipy():
   assert relerr(P, ref_P) <= 5e-2, f"P vs scipy: relerr {relerr(P, ref_P)}"
 
 
+@pytest.mark.parametrize("m,n", [(8, 5), (5, 8)])
+def test_polar_rectangular(m, n):
+  """Both rectangular shapes: U is [m,n] semi-orthogonal, P is [n,n], and U P reconstructs A."""
+  A = np.asarray(np.random.default_rng(3).standard_normal((m, n)), f16)
+  U, P = L.polar(A)
+  assert U.shape == (m, n) and P.shape == (n, n)
+  assert relerr(U @ P, A) <= 5e-2, f"polar recon [{m},{n}]: relerr {relerr(U @ P, A)}"
+  I = U.T @ U if m >= n else U @ U.T
+  assert relerr(I, np.eye(I.shape[0])) <= 5e-2, f"U not semi-orthogonal: relerr {relerr(I, np.eye(I.shape[0]))}"
+
+
 def test_polar_rejects():
   with pytest.raises(ValueError): L.polar(np.zeros(4, f16))              # not 2-D
   with pytest.raises(ValueError): L.polar(np.zeros((2, 2, 2), f16))

Scheduling note: this, #260 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.