diff --git a/CHANGELOG.md b/CHANGELOG.md index 0e90ea3..7e2e747 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 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 diff --git a/src/eigen.rs b/src/eigen.rs new file mode 100644 index 0000000..51e9fbe --- /dev/null +++ b/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)>, 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 { + 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. `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) -> Option { + 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::(); + 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, 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!(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] + 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..fe1cd65 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,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 diff --git a/src/qp.rs b/src/qp.rs index 1dd78b8..833d4e2 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 (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)); - 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)).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; + 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,25 +409,15 @@ 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(); - direction += &v * v.dot(c); - } - } + let problem = Scaled::new(q, e)?; + let z = Constraints::new(&problem.e).null; + let (eigen, cutoff) = reduced_hessian(&problem.k, &z)?; + 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 {