🐛 Replace nalgebra's SVD with a null-space KKT solve and Jacobi eigendecomposition

CHANGELOG.md

@@ -17,11 +17,13 @@ All notable changes to this project are documented here. The format follows
 
 ### Fixed
 
-- Min variance, mean-variance and max Sharpe could return `SingularCovariance` for a
-  well-conditioned covariance, such as one heavily shrunk. nalgebra's SVD with a convergence
-  tolerance of exactly machine epsilon can return a decomposition that doesn't reconstruct
-  the matrix. It now uses nalgebra's default tolerance, and KKT solutions are refined
-  against the exact system.
+- Min variance, mean-variance and max Sharpe could return `SingularCovariance`, or slightly
+  inaccurate weights, on valid input: a heavily shrunk covariance, or a rank-deficient one
+  from fewer observations than assets. nalgebra's SVD can return factors that don't
+  reconstruct the matrix, by up to 1e-2. cutup no longer uses it. Equality-constrained
+  problems are solved by the null-space method, with the reduced system and max Sharpe's
+  covariance decomposed by Jacobi rotations, which reconstruct it to within rounding. Results are
+  unchanged up to rounding, including the minimum-norm weights on a singular covariance.
 
 ## [0.2.3] - 2026-10-23
 

src/eigen.rs

@@ -0,0 +1,261 @@
+//! Symmetric eigendecomposition by the cyclic Jacobi method.
+//!
+//! nalgebra's SVD and symmetric eigendecomposition can both return factors that reconstruct
+//! the matrix only to 1e-3 or worse, on ordinary covariance and KKT matrices. Jacobi rotations
+//! are backward stable: the factors reconstruct the matrix to within rounding of its norm, so
+//! cutup uses them wherever it needs a pseudo-inverse or a null space. They start from
+//! nalgebra's eigenvectors, which are close enough that a sweep or two usually finishes the
+//! job.
+
+use nalgebra::{DMatrix, DVector};
+
+/// Sweeps allowed before giving up. Convergence is quadratic, so well under 20 are needed in
+/// practice.
+const MAX_SWEEPS: usize = 100;
+
+/// Eigenvectors spanning a range, each with its eigenvalue, and ones spanning a null space.
+type Split = (Vec<(f64, DVector<f64>)>, Vec<DVector<f64>>);
+
+/// QR iterations per row allowed for nalgebra's eigendecomposition, which only supplies the
+/// starting point, as in LAPACK. Past this, Jacobi starts from the identity instead.
+const START_ITERATIONS_PER_ROW: usize = 30;
+
+/// `A = V diag(λ) Vᵀ` for a symmetric `A`: eigenvalues `λ`, unordered, and orthonormal
+/// eigenvectors in the matching columns of `V`.
+pub(crate) struct SymmetricEigen {
+    pub(crate) values: DVector<f64>,
+    pub(crate) vectors: DMatrix<f64>,
+}
+
+impl SymmetricEigen {
+    /// Decomposes the symmetric matrix `a` by cyclic Jacobi rotations (Golub and Van Loan,
+    /// 2013, Algorithm 8.5.3). Each rotation zeroes one off-diagonal pair, and sweeps repeat
+    /// until the off-diagonal part is below rounding of the whole.
+    ///
+    /// The rotations apply to `VᵀAV`, where `V` is nalgebra's eigenvector matrix made exactly
+    /// orthonormal by a QR decomposition, or the identity if nalgebra doesn't converge. Any
+    /// orthonormal `V` gives the same eigenvalues; a good one leaves little to rotate. `A` is
+    /// first divided by its largest entry, so the stopping test can't overflow or underflow.
+    ///
+    /// Returns `None` if `a` has NaN or infinite entries, or the sweeps run out.
+    pub(crate) fn new(a: &DMatrix<f64>) -> Option<Self> {
+        let n = a.nrows();
+        if a.iter().any(|x| !x.is_finite()) {
+            return None;
+        }
+        let scale = a.amax();
+        if scale == 0.0 {
+            return Some(SymmetricEigen {
+                values: DVector::zeros(n),
+                vectors: DMatrix::identity(n, n),
+            });
+        }
+        let a = a / scale;
+        let start = a
+            .clone()
+            .try_symmetric_eigen(5.0 * f64::EPSILON, START_ITERATIONS_PER_ROW * n);
+        let mut v = match start {
+            Some(start) => start.eigenvectors.qr().q(),
+            None => DMatrix::identity(n, n),
+        };
+        let rotated = v.transpose() * a * &v;
+        let mut a = (&rotated + rotated.transpose()) * 0.5;
+        let tolerance = f64::EPSILON * a.norm();
+
+        for _ in 0..MAX_SWEEPS {
+            let off_diagonal = (0..n)
+                .flat_map(|p| (p + 1..n).map(move |q| (p, q)))
+                .map(|(p, q)| a[(p, q)] * a[(p, q)])
+                .sum::<f64>();
+            if off_diagonal.is_nan() {
+                break;
+            }
+            if (2.0 * off_diagonal).sqrt() <= tolerance {
+                return Some(SymmetricEigen {
+                    values: a.diagonal() * scale,
+                    vectors: v,
+                });
+            }
+            for p in 0..n {
+                for q in p + 1..n {
+                    rotate(&mut a, &mut v, p, q);
+                }
+            }
+        }
+        None
+    }
+
+    /// For a positive semidefinite matrix, splits the eigenvectors into those of eigenvalues
+    /// above `cutoff`, paired with their eigenvalue, which span the range, and the rest, which
+    /// span the null space. Eigenvalues at or below `cutoff`, negative ones included, are
+    /// rounding error.
+    pub(crate) fn split(&self, cutoff: f64) -> Split {
+        let (mut range, mut null) = (Vec::new(), Vec::new());
+        for (k, &value) in self.values.iter().enumerate() {
+            let vector = self.vectors.column(k).into_owned();
+            if value > cutoff {
+                range.push((value, vector));
+            } else {
+                null.push(vector);
+            }
+        }
+        (range, null)
+    }
+}
+
+/// Applies the Jacobi rotation that zeroes `a[(p, q)]`: `A ← JᵀAJ` and `V ← VJ`.
+fn rotate(a: &mut DMatrix<f64>, v: &mut DMatrix<f64>, p: usize, q: usize) {
+    let apq = a[(p, q)];
+    if apq == 0.0 {
+        return;
+    }
+    // t = tan θ is the smaller root of t² + 2τt − 1 = 0, which keeps the rotation below 45°.
+    // When τ is so large that τ² overflows, t ≈ 1 / 2τ.
+    let tau = (a[(q, q)] - a[(p, p)]) / (2.0 * apq);
+    let t = if tau.abs() > 1e150 {
+        0.5 / tau
+    } else {
+        tau.signum() / (tau.abs() + (tau * tau + 1.0).sqrt())
+    };
+    let c = 1.0 / (t * t + 1.0).sqrt();
+    let s = t * c;
+
+    let n = a.nrows();
+    for k in 0..n {
+        let (akp, akq) = (a[(k, p)], a[(k, q)]);
+        a[(k, p)] = c * akp - s * akq;
+        a[(k, q)] = s * akp + c * akq;
+    }
+    for k in 0..n {
+        let (apk, aqk) = (a[(p, k)], a[(q, k)]);
+        a[(p, k)] = c * apk - s * aqk;
+        a[(q, k)] = s * apk + c * aqk;
+    }
+    a[(p, q)] = 0.0;
+    a[(q, p)] = 0.0;
+    for k in 0..n {
+        let (vkp, vkq) = (v[(k, p)], v[(k, q)]);
+        v[(k, p)] = c * vkp - s * vkq;
+        v[(k, q)] = s * vkp + c * vkq;
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use nalgebra::dmatrix;
+    use proptest::prelude::*;
+
+    /// Checks that the decomposition reconstructs `a` and that `V` is orthonormal, both to
+    /// within rounding.
+    fn assert_decomposes(a: &DMatrix<f64>) {
+        let n = a.nrows();
+        let eigen = SymmetricEigen::new(a).unwrap();
+        let v = &eigen.vectors;
+        let rebuilt = v * DMatrix::from_diagonal(&eigen.values) * v.transpose();
+        let scale = a.norm().max(f64::MIN_POSITIVE);
+        assert!(
+            (&rebuilt - a).norm() <= 1e-14 * scale * n as f64,
+            "{a} rebuilt as {rebuilt}"
+        );
+        let identity = DMatrix::<f64>::identity(n, n);
+        assert!((v.transpose() * v - identity).amax() <= 1e-14 * n as f64);
+    }
+
+    #[test]
+    fn known_spectrum() {
+        let eigen = SymmetricEigen::new(&dmatrix![2.0, 1.0; 1.0, 2.0]).unwrap();
+        let mut values: Vec<f64> = eigen.values.iter().copied().collect();
+        values.sort_by(f64::total_cmp);
+        assert!((values[0] - 1.0).abs() < 1e-14 && (values[1] - 3.0).abs() < 1e-14);
+    }
+
+    #[test]
+    fn diagonal_and_empty_matrices() {
+        let eigen = SymmetricEigen::new(&DMatrix::from_diagonal_element(3, 3, 4.0)).unwrap();
+        assert!(eigen
+            .values
+            .iter()
+            .all(|&value| (value - 4.0).abs() < 1e-14));
+        assert_decomposes(&DMatrix::from_diagonal(&DVector::from_column_slice(&[
+            3.0, -1.0, 0.0, 2.0,
+        ])));
+        assert_eq!(
+            SymmetricEigen::new(&DMatrix::zeros(0, 0))
+                .unwrap()
+                .values
+                .len(),
+            0
+        );
+        assert_decomposes(&DMatrix::zeros(4, 4));
+    }
+
+    #[test]
+    fn rejects_non_finite_matrices() {
+        let a = dmatrix![1.0, f64::NAN; f64::NAN, 1.0];
+        assert!(SymmetricEigen::new(&a).is_none());
+    }
+
+    #[test]
+    fn extreme_scales() {
+        // At 1e300 the squared entries overflow, and at 1e-300 they underflow, unless the
+        // matrix is scaled first.
+        let a = dmatrix![2.0, 1.0; 1.0, 2.0];
+        for scale in [1e300, 1e-300] {
+            let eigen = SymmetricEigen::new(&(&a * scale)).unwrap();
+            let mut values: Vec<f64> = eigen.values.iter().map(|v| v / scale).collect();
+            values.sort_by(f64::total_cmp);
+            assert!((values[0] - 1.0).abs() < 1e-14 && (values[1] - 3.0).abs() < 1e-14);
+        }
+    }
+
+    #[test]
+    fn split_puts_negative_rounding_in_the_null_space() {
+        let eigen = SymmetricEigen {
+            values: DVector::from_column_slice(&[2.0, -1e-3, 1e-16]),
+            vectors: DMatrix::identity(3, 3),
+        };
+        let (range, null) = eigen.split(1e-15);
+        assert_eq!(range.len(), 1);
+        assert_eq!(range[0].0, 2.0);
+        assert_eq!(null.len(), 2);
+    }
+
+    #[test]
+    fn singular_and_indefinite_matrices() {
+        // Rank one: one eigenvalue of 14, the rest zero.
+        let x = DVector::from_column_slice(&[1.0, 2.0, 3.0]);
+        let rank_one = &x * x.transpose();
+        assert_decomposes(&rank_one);
+        let eigen = SymmetricEigen::new(&rank_one).unwrap();
+        assert!((eigen.values.amax() - 14.0).abs() < 1e-13);
+        assert!((eigen.values.sum() - 14.0).abs() < 1e-13);
+
+        // A KKT matrix, which has negative eigenvalues.
+        assert_decomposes(&dmatrix![
+            1.0, 0.5, 1.0;
+            0.5, 2.0, 1.0;
+            1.0, 1.0, 0.0
+        ]);
+    }
+
+    proptest! {
+        #[test]
+        fn decomposes_random_symmetric_matrices(
+            n in 1usize..12,
+            values in prop::collection::vec(-1.0f64..1.0, 144),
+            rank in 0usize..12,
+            exponent in -8i32..8,
+        ) {
+            // Gram matrices of `rank` random vectors, so many are singular, plus an
+            // indefinite symmetric part, at varied scales.
+            let x = DMatrix::from_fn(n, rank.min(n), |i, j| values[i * 12 + j]);
+            let mut a = &x * x.transpose();
+            if rank % 2 == 1 {
+                a += DMatrix::from_fn(n, n, |i, j| values[(i + j) % 144] * 0.1);
+            }
+            a *= 10f64.powi(exponent);
+            assert_decomposes(&a);
+        }
+    }
+}

src/lib.rs

@@ -2,10 +2,13 @@
 #![warn(missing_docs)]
 
 mod cluster;
+mod eigen;
 mod qp;
 
 pub use cluster::Linkage;
 
+use eigen::SymmetricEigen;
+
 use nalgebra::{DMatrix, DVector};
 use std::collections::HashMap;
 use std::fmt;
@@ -1188,22 +1191,9 @@ impl PortfolioAllocator {
         if !(scale.is_finite() && scale > 0.0) {
             return Err(CutupError::SingularCovariance);
         }
-        // nalgebra's SVD never converges on NaN, which the check above rules out; the
-        // iteration cap is a backstop.
-        let svd = (cov / scale)
-            .try_svd(false, true, qp::SVD_EPSILON, 10_000)
-            .ok_or(CutupError::SingularCovariance)?;
-        let v_t = svd.v_t.as_ref().ok_or(CutupError::SingularCovariance)?;
-        let cutoff = svd.singular_values.max() * n as f64 * f64::EPSILON;
-        let (mut range, mut null) = (Vec::new(), Vec::new());
-        for k in 0..n {
-            let direction = v_t.row(k).transpose();
-            if svd.singular_values[k] > cutoff {
-                range.push((svd.singular_values[k], direction));
-            } else {
-                null.push(direction);
-            }
-        }
+        // Σ is positive semidefinite, so its eigenvalues are its singular values.
+        let eigen = SymmetricEigen::new(&(cov / scale)).ok_or(CutupError::SingularCovariance)?;
+        let (range, null) = eigen.split(eigen.values.amax() * n as f64 * f64::EPSILON);
 
         if !null.is_empty() {
             // Coordinates of 1 and of the excess returns in the null space.
@@ -2796,6 +2786,41 @@ mod tests {
         }
     }
 
+    #[test]
+    fn test_min_variance_on_rank_deficient_kkt() {
+        // Three rows for four assets give Σ rank 2, and its KKT matrix a zero singular value.
+        // nalgebra's SVD left a residual of about 4e-7 here, so the solve was rejected as
+        // inconsistent; NumPy's solves it to rounding.
+        let returns = DMatrix::from_vec(
+            3,
+            4,
+            vec![
+                0.03923072979353582,
+                -0.027034431663761255,
+                -0.02017662174798146,
+                -0.030105448299475653,
+                0.02262405031858968,
+                0.01653614746213769,
+                -0.009248234151376263,
+                0.02095233950501841,
+                -0.049975234461268725,
+                -0.010174824524215684,
+                -0.043910804533156975,
+                0.04613336884944048,
+            ],
+        );
+        let allocator = PortfolioAllocator::from_returns(returns).unwrap();
+        // numpy.linalg.pinv's minimum-norm solution of the same KKT system.
+        let expected = Weights(vec![0.230123253, 0.270499312, 0.278914760, 0.220462675]);
+        assert_weights_close(
+            &allocator.min_variance_allocation().unwrap(),
+            &expected,
+            1e-9,
+        );
+        let weights = allocator.long_only_min_variance_allocation().unwrap();
+        assert!(portfolio_variance(&allocator, &weights) < 1e-18);
+    }
+
     #[test]
     fn test_long_only_min_variance_matches_unconstrained_when_long() {
         // A diagonal covariance gives all-positive minimum variance weights, so the