[thinkjrs]

Removes nalgebra's SVD from cutup. On rank-deficient and heavily shrunk covariances it could return factors that don't reconstruct the matrix, so min variance, mean-variance, max Sharpe and long-only min variance returned SingularCovariance, or slightly wrong weights, on valid input. Equality-constrained problems are now solved by the null-space method, and the decompositions that need a pseudo-inverse or a null space use a Jacobi eigensolver written for cutup. Weights are unchanged up to rounding, including the minimum-norm weights on a singular covariance. No API changes.

The problem

PR #11 moved the SVD to nalgebra's default tolerance and added iterative refinement. That fixed the shrunk-covariance case, but not rank-deficient ones:

• On positive semidefinite matrices, try_svd at 5ε still reconstructed with errors as large as 1.13, on 665 of 200k random matrices. At ε the worst error was 3.5, on 4,594 of them.
• nalgebra's symmetric_eigen does better, but still reaches errors of 1e-9, on 61 of 200k.
• ColPivQR pivots on the largest single entry rather than the largest column norm, so it doesn't reveal rank and can't stand in.

On a stress run of 60k random singular problems, about 4 failed with SingularCovariance. One failure is plain min variance on 3 observations of 4 assets, which is now a regression test.

Upgrading nalgebra isn't an option in a patch release, since nalgebra types are in cutup's public API.

Jacobi eigensolver (src/eigen.rs)

SymmetricEigen::new uses cyclic Jacobi (Golub and Van Loan, Algorithm 8.5.3). Jacobi rotations are backward stable: the factors reconstruct the matrix to within rounding of its norm.

• Warm start. Pure Jacobi took 3.8 s on a 300-asset long-only problem. The rotations now apply to VᵀAV, where V holds nalgebra's eigenvectors made exactly orthonormal by a QR decomposition. Any orthonormal V gives the same eigenvalues, and a good one leaves only a sweep or two of work. nalgebra runs at 5ε with a cap of 30·n QR iterations, as in LAPACK; if it doesn't converge, V is the identity.
• Scaling. A is divided by its largest entry first, so the stopping test can't overflow or underflow at extreme scales.
• Stopping. Sweeps stop when the off-diagonal norm is within ε‖A‖. Non-finite input, or 100 sweeps without converging, returns None, which callers report as SingularCovariance, the error the SVD path returned.
• Split. split(cutoff) separates the range from the null space for positive semidefinite matrices. Eigenvalues at or below the cutoff, negative rounding included, go to the null space. The KKT solve, the zero-curvature direction and max Sharpe all use it.
• Large τ. When τ² would overflow, the rotation uses t ≈ 1/2τ.

Null-space KKT solve (src/qp.rs)

solve_kkt solves min ½xᵀKx − cᵀx subject to Eᵀx = b by the null-space method (Nocedal and Wright, §16.2), in place of a pseudo-inverse of the whole KKT matrix:

1. Scale. K is divided by its average diagonal entry, and each constraint column by its largest entry.
2. Constraints. A pivoted Gram–Schmidt with reorthogonalization gives E = YT, where Y is an orthonormal basis of E's range. Dependent columns are dropped, which still detects DegenerateExpectedReturns. Householder QR of Y completes the null-space basis Z.
3. Feasible point. x₀ = Yu is the smallest point satisfying the constraints, found by forward substitution on T. Inconsistent constraints fail the residual check below.
4. Reduced system. ZᵀKZ is decomposed by the Jacobi solver. Its pseudo-inverse, over eigenvalues above the cutoff max(|λ|, |K|)·(n+1)·ε, gives v, and x = x₀ + Zv is the minimum-norm optimum. That's the same answer the pseudo-inverse of the KKT matrix gave.
5. Multipliers. λ comes from back substitution on T against c − Kx, with zeros for dependent constraints.
6. Check. If the stationarity and feasibility residuals together exceed 1e-8 of the problem's size, there is no optimum and the solve returns None, as before. A noise floor lets the active-set solver pass its own rounding.

The SVD tolerance constant and the iterative refinement are gone. The active-set solver's zero-curvature direction, used when Q is singular, is built from the same Z and the null eigenvectors of the reduced system.

Max Sharpe (src/lib.rs)

Max Sharpe splits Σ into its range and null space with SymmetricEigen::split instead of the SVD. Σ is positive semidefinite, so its eigenvalues are its singular values and the cutoff is unchanged.

Tests

• Eigensolver.
  • A known spectrum, also at scales of 1e300 and 1e-300.
  • Negative rounding going to the null space in split.
  • Diagonal, zero and empty matrices.
  • Non-finite input.
  • A rank-one matrix and an indefinite KKT matrix.
  • A property test on random symmetric matrices, many of them singular or indefinite, at scales from 1e-8 to 1e8. It checks reconstruction to 1e-14·‖A‖·n and that V is orthonormal.
• test_min_variance_on_rank_deficient_kkt. Min variance on 3 observations of 4 assets matches NumPy's minimum-norm answer to 1e-9, and long-only min variance reaches zero variance. Plain min variance fails on this input on develop.
• Stress. 0 errors and 0 KKT violations in 60k random rank-deficient problems, down from about 4 errors.
• Properties. Three runs at PROPTEST_CASES=50000 pass. That includes max_sharpe_ignores_column_order_and_joint_scale, which used to fail about one run in three.
• Reference fixtures. All 138 cases pass, and the drift check reports 0 of 138 drifted.
• Suite. The full suite, clippy, docs and fmt pass on stable and on MSRV 1.78.

Performance

Problem	Before	After
300 assets, long-only min variance	about 1.3 s	1.06 s
300 assets, min variance		37 ms
100 assets, long-only min variance		17 ms

Long-only problems still refactor the KKT system at every active-set step. Updating the factorization instead is review finding 7 from #11, left for later.

[thinkjrs]

Duplicate decomposition. When solve_noisy_kkt returns None in the active-set loop, descent_ray is called on the same q_free/e_free and rebuilds everything: Scaled::new, Gram–Schmidt, the full Householder QR, ZᵀQZ, the nalgebra eigen warm start and the Jacobi sweeps. That is O(n³) work done twice per zero-curvature step. Return (z, eigen, cutoff) from the failed solve, or factor a shared struct, so the ray reuses it.

[thinkjrs]

Deferred. This is the same work as finding 7 on #11, reusing factorizations across active-set steps, so it goes with that. It only runs on zero-curvature steps, which plain min variance never takes.