From ea10291cfa7a0e08d5bac6a0d0c9886759fb2063 Mon Sep 17 00:00:00 2001 From: "Jason R. Stevens, CFA" Date: Mon, 26 Oct 2026 20:38:55 -0500 Subject: [PATCH 1/4] =?UTF-8?q?=F0=9F=90=9B=20Replace=20nalgebra's=20SVD?= =?UTF-8?q?=20with=20a=20null-space=20KKT=20solve=20and=20Jacobi=20eigende?= =?UTF-8?q?composition?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/eigen.rs | 210 ++++++++++++++++++++++++++++++++++++++ src/lib.rs | 56 +++++++++-- src/qp.rs | 280 ++++++++++++++++++++++++++++++++++++--------------- 3 files changed, 456 insertions(+), 90 deletions(-) create mode 100644 src/eigen.rs diff --git a/src/eigen.rs b/src/eigen.rs new file mode 100644 index 0000000..c08b5b5 --- /dev/null +++ b/src/eigen.rs @@ -0,0 +1,210 @@ +//! 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 accurate to rounding (Demmel and Veselić, 1992), 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 crate::CutupError; +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; + +/// `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, + pub(crate) vectors: DMatrix, +} + +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. + /// + /// # Errors + /// + /// * [`CutupError::NotConverged`] if `a` has NaN or infinite entries, or the sweeps run + /// out. + pub(crate) fn new(a: &DMatrix) -> Result { + let n = a.nrows(); + if a.iter().any(|x| !x.is_finite()) { + return Err(CutupError::NotConverged); + } + // nalgebra panics on an empty matrix. + let start = (n > 0) + .then(|| a.clone().try_symmetric_eigen(f64::EPSILON, 10_000)) + .flatten(); + 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::(); + if off_diagonal.is_nan() { + break; + } + if (2.0 * off_diagonal).sqrt() <= tolerance { + return Ok(SymmetricEigen { + values: a.diagonal(), + vectors: v, + }); + } + for p in 0..n { + for q in p + 1..n { + rotate(&mut a, &mut v, p, q); + } + } + } + Err(CutupError::NotConverged) + } +} + +/// Applies the Jacobi rotation that zeroes `a[(p, q)]`: `A ← JᵀAJ` and `V ← VJ`. +fn rotate(a: &mut DMatrix, v: &mut DMatrix, 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) { + 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::::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 = 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!(matches!( + SymmetricEigen::new(&a), + Err(CutupError::NotConverged) + )); + } + + #[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); + } + } +} diff --git a/src/lib.rs b/src/lib.rs index d71feee..38de569 100644 --- a/src/lib.rs +++ b/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,18 +1191,14 @@ 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; + // Σ is symmetric, so its eigenvalues' magnitudes are its singular values. + let eigen = SymmetricEigen::new(&(cov / scale))?; + let cutoff = eigen.values.amax() * 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)); + for (k, &value) in eigen.values.iter().enumerate() { + let direction = eigen.vectors.column(k).into_owned(); + if value.abs() > cutoff { + range.push((value, direction)); } else { null.push(direction); } @@ -2796,6 +2795,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 diff --git a/src/qp.rs b/src/qp.rs index 1dd78b8..08b87ba 100644 --- a/src/qp.rs +++ b/src/qp.rs @@ -1,25 +1,26 @@ -//! Convex quadratic programs: equality-constrained ones solved through their KKT system, and +//! Convex quadratic programs: equality-constrained ones solved by the null-space method, and //! ones with nonnegative variables solved by a primal active-set method on top of it. +use crate::eigen::SymmetricEigen; use crate::CutupError; use nalgebra::{DMatrix, DVector}; -/// Convergence tolerance for nalgebra's SVD: its own default. With exactly `ε` the iteration -/// can stop with singular vectors that reconstruct the matrix only to about 1e-3, and even at -/// this tolerance rare matrices reconstruct only to about 1e-10. -pub(crate) const SVD_EPSILON: f64 = 5.0 * f64::EPSILON; - /// A solution `x` and its constraint multipliers `λ`. type Solution = (DVector, DVector); -/// Solves `minimize ½ xᵀKx − cᵀx subject to Eᵀx = b` through its KKT system -/// `[K E; Eᵀ 0] [x; λ] = [c; b]`, returning `x` and the multipliers `λ`. +/// Solves `minimize ½ xᵀKx − cᵀx subject to Eᵀx = b` for a positive semidefinite `K`, +/// returning `x` and the multipliers `λ` of the KKT system `Kx + Eλ = c`, `Eᵀx = b`. +/// +/// This is the null-space method (Nocedal and Wright, 2006, §16.2). With `Y` an orthonormal +/// basis for the range of `E` and `Z` one for the null space of `Eᵀ`, the smallest feasible +/// point is `x₀ = Yu`, and every feasible point is `x₀ + Zv`. The objective is then a +/// quadratic in `v` with the positive semidefinite Hessian `ZᵀKZ`, solved with its +/// pseudo-inverse (eigenvalues below `size · ε` of the largest count as zero, matching +/// `numpy.linalg.pinv`). Since `Y` and `Z` are orthogonal, a singular `K` still gives an exact +/// optimum, the one with the smallest norm. /// -/// The system is solved with the pseudo-inverse (singular values below `size · ε · σ_max` -/// count as zero, matching `numpy.linalg.pinv`), so a singular `K` still gives an exact -/// optimum: the one with the smallest norm. Two rounds of iterative refinement follow. The -/// solution is then checked against the system. If it doesn't satisfy it, the problem has no -/// optimum, either because the objective is unbounded or because the constraints are +/// The solution is then checked against the KKT system. If it doesn't satisfy it, the problem +/// has no optimum, either because the objective is unbounded or because the constraints are /// infeasible, and `Ok(None)` is returned. pub(crate) fn solve_kkt( k: &DMatrix, @@ -39,73 +40,203 @@ fn solve_noisy_kkt( b: &DVector, noise: f64, ) -> Result, CutupError> { - // nalgebra's SVD never converges on NaN, so check first and cap the iterations below. if k.iter().chain(c.iter()).any(|x| !x.is_finite()) { return Err(CutupError::NonFiniteCovariance); } - let (n, m) = (k.nrows(), e.ncols()); - - // Put every block on a comparable scale so the pseudo-inverse cutoff treats them alike. - // Scaling K and c together, or a constraint column of E with its entry of b, leaves the - // solution unchanged and scales the multipliers, which are scaled back at the end. - let k_scale = k_scale(k)?; - let column_scales: Vec = (0..m).map(|j| column_scale(e, j)).collect(); - let mut kkt = DMatrix::zeros(n + m, n + m); - let mut rhs = DVector::zeros(n + m); - kkt.view_mut((0, 0), (n, n)).copy_from(&(k / k_scale)); - rhs.rows_mut(0, n).copy_from(&(c / k_scale)); - for (j, &scale) in column_scales.iter().enumerate() { - let column = e.column(j) / scale; - kkt.view_mut((0, n + j), (n, 1)).copy_from(&column); - kkt.view_mut((n + j, 0), (1, n)) - .copy_from(&column.transpose()); - rhs[n + j] = b[j] / scale; - } + let problem = Scaled::new(k, e)?; + let (k, e) = (&problem.k, &problem.e); + let c = c / problem.k_scale; + let b = DVector::from_fn(b.len(), |j, _| b[j] / problem.column_scales[j]); - let size = n + m; - let svd = kkt - .clone() - .try_svd(true, true, SVD_EPSILON, 10_000) - .ok_or(CutupError::SingularCovariance)?; - let tolerance = svd.singular_values.max() * size as f64 * f64::EPSILON; - let pinv = svd - .pseudo_inverse(tolerance) - .map_err(|_| CutupError::SingularCovariance)?; - // Iterative refinement against the exact system removes the error of an inaccurate SVD. - let mut solution = &pinv * &rhs; - for _ in 0..2 { - solution += &pinv * (&rhs - &kkt * &solution); + let constraints = Constraints::new(e); + let x0 = constraints.smallest_feasible(&b); + let z = &constraints.null; + let gradient = z.transpose() * (&c - k * &x0); + let (eigen, cutoff) = reduced_hessian(k, z)?; + let mut v = DVector::zeros(z.ncols()); + for (i, &value) in eigen.values.iter().enumerate() { + if value.abs() > cutoff { + let w = eigen.vectors.column(i); + v += w * (w.dot(&gradient) / value); + } } + let x = x0 + z * v; + let multipliers = constraints.multipliers(&(&c - k * &x)); - let residual = (&kkt * &solution - &rhs).norm(); - if residual > 1e-8 * (kkt.norm() * solution.norm() + rhs.norm() + noise / k_scale) { + let stationarity = (k * &x + e * &multipliers - &c).norm(); + let feasibility = (e.transpose() * &x - &b).norm(); + let size = k.norm() * x.norm() + + e.norm() * (x.norm() + multipliers.norm()) + + c.norm() + + b.norm() + + noise / problem.k_scale; + if stationarity + feasibility > 1e-8 * size { return Ok(None); } - let x = solution.rows(0, n).into_owned(); if x.iter().any(|w| !w.is_finite()) { return Err(CutupError::SingularCovariance); } - let multipliers = DVector::from_fn(m, |j, _| solution[n + j] * k_scale / column_scales[j]); + let multipliers = DVector::from_fn(multipliers.len(), |j, _| { + multipliers[j] * problem.k_scale / problem.column_scales[j] + }); Ok(Some((x, multipliers))) } -/// The average diagonal entry of `k`, or 1 when `k` is zero. -fn k_scale(k: &DMatrix) -> Result { - let scale = k.trace() / k.nrows() as f64; - if scale == 0.0 { - Ok(1.0) - } else if scale.is_finite() && scale > 0.0 { - Ok(scale) - } else { - Err(CutupError::SingularCovariance) +/// `K` and `E` put on a comparable scale, so that the cutoffs treat them alike: `K` divided by +/// its average diagonal entry, and each column of `E` by its largest entry. Scaling `K` and `c` +/// together, or a column of `E` with its entry of `b`, leaves the solution unchanged and +/// scales the multipliers. +struct Scaled { + k: DMatrix, + e: DMatrix, + k_scale: f64, + column_scales: Vec, +} + +impl Scaled { + fn new(k: &DMatrix, e: &DMatrix) -> Result { + // A zero K, as the free block of a linear objective is, keeps scale 1. + let average = k.trace() / k.nrows() as f64; + let k_scale = if average == 0.0 { + 1.0 + } else if average.is_finite() && average > 0.0 { + average + } else { + return Err(CutupError::SingularCovariance); + }; + let column_scales: Vec = e + .column_iter() + .map(|column| match column.amax() { + s if s > 0.0 => s, + _ => 1.0, + }) + .collect(); + let e = DMatrix::from_fn(e.nrows(), e.ncols(), |i, j| e[(i, j)] / column_scales[j]); + Ok(Scaled { + k: k / k_scale, + e, + k_scale, + column_scales, + }) } } -/// The largest absolute entry of column `j` of `e`, or 1 when the column is zero. -fn column_scale(e: &DMatrix, j: usize) -> f64 { - match e.column(j).amax() { - s if s > 0.0 => s, - _ => 1.0, +/// The eigendecomposition of the reduced Hessian `ZᵀKZ`, and the cutoff below which its +/// eigenvalues count as zero. +fn reduced_hessian( + k: &DMatrix, + z: &DMatrix, +) -> Result<(SymmetricEigen, f64), CutupError> { + let eigen = SymmetricEigen::new(&(z.transpose() * k * z))?; + // Relative to K as well, so that a reduced Hessian that is all rounding error is zero. + let largest = eigen.values.amax().max(k.amax()); + let cutoff = largest * (k.nrows() + 1) as f64 * f64::EPSILON; + Ok((eigen, cutoff)) +} + +/// The constraint matrix `E` as `E = YT`, with `Y` an orthonormal basis for its range, and +/// `Z` an orthonormal basis for the null space of `Eᵀ`, which completes `Y`. +/// +/// `Y` comes from Gram–Schmidt with column pivoting and reorthogonalization: each step takes +/// the column of `E` with the most left over after removing the earlier basis vectors, until +/// what is left is rounding. A column that never gets taken depends on the others. `Z` comes +/// from a Householder QR of `Y`. +struct Constraints { + range: DMatrix, + null: DMatrix, + /// `T`, with a row per basis vector and a column per column of `E`. + coefficients: DMatrix, + /// The column of `E` each basis vector was built from, in order. + pivots: Vec, +} + +impl Constraints { + fn new(e: &DMatrix) -> Self { + let (n, m) = e.shape(); + let tolerance = + n.max(m) as f64 * f64::EPSILON * e.column_iter().map(|c| c.norm()).fold(0.0, f64::max); + let mut left = e.clone(); + let mut coefficients = DMatrix::zeros(n.min(m), m); + let mut basis: Vec> = Vec::new(); + let mut pivots = Vec::new(); + let mut remaining: Vec = (0..m).collect(); + + while basis.len() < n { + let Some((position, norm)) = remaining + .iter() + .map(|&j| left.column(j).norm()) + .enumerate() + .max_by(|a, b| a.1.total_cmp(&b.1)) + else { + break; + }; + if norm <= tolerance { + break; + } + let j = remaining.remove(position); + let row = basis.len(); + let y = left.column(j) / norm; + coefficients[(row, j)] = norm; + for &l in &remaining { + // Twice, so that rounding in the first pass doesn't leave y in the column. + for _ in 0..2 { + let t = y.dot(&left.column(l)); + coefficients[(row, l)] += t; + let updated = left.column(l) - &y * t; + left.set_column(l, &updated); + } + } + basis.push(y); + pivots.push(j); + } + + let rank = basis.len(); + let range = if rank == 0 { + DMatrix::zeros(n, 0) + } else { + DMatrix::from_columns(&basis) + }; + // The last n − rank columns of a full Q, with Q's first columns spanning range(Y). + let null = if rank == 0 { + DMatrix::identity(n, n) + } else { + let mut q_transpose = DMatrix::identity(n, n); + range.clone().qr().q_tr_mul(&mut q_transpose); + q_transpose.rows(rank, n - rank).transpose() + }; + Constraints { + range, + null, + coefficients: coefficients.rows(0, rank).into_owned(), + pivots, + } + } + + /// The smallest `x` with `Eᵀx = b` on the pivot columns: `x = Yu` with `Tᵀu = b`, solved + /// by forward substitution, since column `pᵢ` of `T` is zero below row `i`. A dependent + /// column's equation is left to the caller to check. + fn smallest_feasible(&self, b: &DVector) -> DVector { + let mut u = DVector::zeros(self.pivots.len()); + for (i, &j) in self.pivots.iter().enumerate() { + let known: f64 = (0..i).map(|k| self.coefficients[(k, j)] * u[k]).sum(); + u[i] = (b[j] - known) / self.coefficients[(i, j)]; + } + &self.range * u + } + + /// The `λ` that best solves `Eλ = g`: `Tλ = Yᵀg` on the pivot columns, solved by back + /// substitution, with zero for dependent columns. + fn multipliers(&self, g: &DVector) -> DVector { + let h = self.range.transpose() * g; + let mut lambda = DVector::zeros(self.coefficients.ncols()); + for (i, &j) in self.pivots.iter().enumerate().rev() { + let known: f64 = self.pivots[i + 1..] + .iter() + .map(|&l| self.coefficients[(i, l)] * lambda[l]) + .sum(); + lambda[j] = (h[i] - known) / self.coefficients[(i, j)]; + } + lambda } } @@ -278,22 +409,13 @@ fn descent_ray( c: &DVector, e: &DMatrix, ) -> Result, CutupError> { - let (n, m) = (q.nrows(), e.ncols()); - let mut stacked = DMatrix::zeros(n + m, n); - stacked.rows_mut(0, n).copy_from(&(q / k_scale(q)?)); - for j in 0..m { - let column = e.column(j) / column_scale(e, j); - stacked.row_mut(n + j).copy_from(&column.transpose()); - } - let svd = stacked - .try_svd(false, true, 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 + m) as f64 * f64::EPSILON; - let mut direction = DVector::zeros(n); - for (k, &sigma) in svd.singular_values.iter().enumerate() { - if sigma <= cutoff { - let v = v_t.row(k).transpose(); + let problem = Scaled::new(q, e)?; + let z = Constraints::new(&problem.e).null; + let (eigen, cutoff) = reduced_hessian(&problem.k, &z)?; + let mut direction = DVector::zeros(q.nrows()); + for (i, &value) in eigen.values.iter().enumerate() { + if value.abs() <= cutoff { + let v = &z * eigen.vectors.column(i); direction += &v * v.dot(c); } } From 1e24386b44e4779082d76de68bb9ba28c61fd45d Mon Sep 17 00:00:00 2001 From: "Jason R. Stevens, CFA" Date: Mon, 26 Oct 2026 20:38:56 -0500 Subject: [PATCH 2/4] =?UTF-8?q?=F0=9F=93=9D=20Describe=20the=20decompositi?= =?UTF-8?q?on=20fix=20in=20the=20changelog?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CHANGELOG.md | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0e90ea3..574abb1 100644 --- a/CHANGELOG.md +++ b/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 are accurate to rounding. Results are + unchanged up to rounding, including the minimum-norm weights on a singular covariance. ## [0.2.3] - 2026-10-23 From e0a55bd3f37c00190c13e54c3d0014e639f19585 Mon Sep 17 00:00:00 2001 From: "Jason R. Stevens, CFA" Date: Mon, 26 Oct 2026 23:04:31 -0500 Subject: [PATCH 3/4] =?UTF-8?q?=F0=9F=90=9B=20Treat=20negative=20rounding?= =?UTF-8?q?=20as=20null=20space,=20scale=20Jacobi=20input=20and=20keep=20S?= =?UTF-8?q?ingularCovariance=20on=20decomposition=20failure?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/eigen.rs | 95 ++++++++++++++++++++++++++++++++++++++++------------ src/lib.rs | 15 ++------- src/qp.rs | 29 ++++++++-------- 3 files changed, 90 insertions(+), 49 deletions(-) diff --git a/src/eigen.rs b/src/eigen.rs index c08b5b5..51e9fbe 100644 --- a/src/eigen.rs +++ b/src/eigen.rs @@ -2,17 +2,24 @@ //! //! 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 accurate to rounding (Demmel and Veselić, 1992), 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. +//! 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 crate::CutupError; 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)>, Vec>); + +/// 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 { @@ -27,21 +34,26 @@ impl SymmetricEigen { /// /// 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. + /// 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. /// - /// # Errors - /// - /// * [`CutupError::NotConverged`] if `a` has NaN or infinite entries, or the sweeps run - /// out. - pub(crate) fn new(a: &DMatrix) -> Result { + /// Returns `None` if `a` has NaN or infinite entries, or the sweeps run out. + pub(crate) fn new(a: &DMatrix) -> Option { let n = a.nrows(); if a.iter().any(|x| !x.is_finite()) { - return Err(CutupError::NotConverged); + return None; + } + let scale = a.amax(); + if scale == 0.0 { + return Some(SymmetricEigen { + values: DVector::zeros(n), + vectors: DMatrix::identity(n, n), + }); } - // nalgebra panics on an empty matrix. - let start = (n > 0) - .then(|| a.clone().try_symmetric_eigen(f64::EPSILON, 10_000)) - .flatten(); + 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), @@ -59,8 +71,8 @@ impl SymmetricEigen { break; } if (2.0 * off_diagonal).sqrt() <= tolerance { - return Ok(SymmetricEigen { - values: a.diagonal(), + return Some(SymmetricEigen { + values: a.diagonal() * scale, vectors: v, }); } @@ -70,7 +82,24 @@ impl SymmetricEigen { } } } - Err(CutupError::NotConverged) + 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) } } @@ -164,10 +193,32 @@ mod tests { #[test] fn rejects_non_finite_matrices() { let a = dmatrix![1.0, f64::NAN; f64::NAN, 1.0]; - assert!(matches!( - SymmetricEigen::new(&a), - Err(CutupError::NotConverged) - )); + 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 = 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] diff --git a/src/lib.rs b/src/lib.rs index 38de569..fe1cd65 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1191,18 +1191,9 @@ impl PortfolioAllocator { if !(scale.is_finite() && scale > 0.0) { return Err(CutupError::SingularCovariance); } - // Σ is symmetric, so its eigenvalues' magnitudes are its singular values. - let eigen = SymmetricEigen::new(&(cov / scale))?; - let cutoff = eigen.values.amax() * n as f64 * f64::EPSILON; - let (mut range, mut null) = (Vec::new(), Vec::new()); - for (k, &value) in eigen.values.iter().enumerate() { - let direction = eigen.vectors.column(k).into_owned(); - if value.abs() > cutoff { - range.push((value, 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. diff --git a/src/qp.rs b/src/qp.rs index 08b87ba..833d4e2 100644 --- a/src/qp.rs +++ b/src/qp.rs @@ -53,13 +53,12 @@ fn solve_noisy_kkt( let z = &constraints.null; let gradient = z.transpose() * (&c - k * &x0); let (eigen, cutoff) = reduced_hessian(k, z)?; - let mut v = DVector::zeros(z.ncols()); - for (i, &value) in eigen.values.iter().enumerate() { - if value.abs() > cutoff { - let w = eigen.vectors.column(i); - v += w * (w.dot(&gradient) / value); - } - } + let (range, _) = eigen.split(cutoff); + let v = range + .iter() + .fold(DVector::zeros(z.ncols()), |v, (value, w)| { + v + w * (w.dot(&gradient) / value) + }); let x = x0 + z * v; let multipliers = constraints.multipliers(&(&c - k * &x)); @@ -127,7 +126,8 @@ fn reduced_hessian( k: &DMatrix, z: &DMatrix, ) -> Result<(SymmetricEigen, f64), CutupError> { - let eigen = SymmetricEigen::new(&(z.transpose() * k * z))?; + let eigen = + SymmetricEigen::new(&(z.transpose() * k * z)).ok_or(CutupError::SingularCovariance)?; // Relative to K as well, so that a reduced Hessian that is all rounding error is zero. let largest = eigen.values.amax().max(k.amax()); let cutoff = largest * (k.nrows() + 1) as f64 * f64::EPSILON; @@ -412,13 +412,12 @@ fn descent_ray( let problem = Scaled::new(q, e)?; let z = Constraints::new(&problem.e).null; let (eigen, cutoff) = reduced_hessian(&problem.k, &z)?; - let mut direction = DVector::zeros(q.nrows()); - for (i, &value) in eigen.values.iter().enumerate() { - if value.abs() <= cutoff { - let v = &z * eigen.vectors.column(i); - direction += &v * v.dot(c); - } - } + let (_, null) = eigen.split(cutoff); + let direction = null.iter().fold(DVector::zeros(q.nrows()), |direction, w| { + let v = &z * w; + let along = v.dot(c); + direction + v * along + }); // The KKT system was inconsistent, so such a direction exists up to rounding; without one // the system was too ill-conditioned to tell. if direction.norm() <= 1e-12 * c.norm() || direction.norm() == 0.0 { From 2b23518250facbdef87f1cb6e5365ec22b806630 Mon Sep 17 00:00:00 2001 From: "Jason R. Stevens, CFA" Date: Mon, 26 Oct 2026 23:04:31 -0500 Subject: [PATCH 4/4] =?UTF-8?q?=F0=9F=93=9D=20Say=20the=20Jacobi=20factors?= =?UTF-8?q?=20are=20accurate=20to=20rounding=20of=20the=20norm?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 574abb1..7e2e747 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,7 +22,7 @@ All notable changes to this project are documented here. The format follows 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 are accurate to rounding. Results are + 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