From b9ec1a61435bb9c0f84b79063f89fbe48d47950b Mon Sep 17 00:00:00 2001 From: Hao Zhu Date: Sun, 8 Mar 2026 09:48:22 +0100 Subject: [PATCH 1/3] fix sum_squares canonicalization for affine expressions in objective --- src/canon/canonicalizer.rs | 44 +++++++++++++++++---- tests/sum_squares_tests.rs | 80 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 116 insertions(+), 8 deletions(-) create mode 100644 tests/sum_squares_tests.rs diff --git a/src/canon/canonicalizer.rs b/src/canon/canonicalizer.rs index 6bea38b..718e13b 100644 --- a/src/canon/canonicalizer.rs +++ b/src/canon/canonicalizer.rs @@ -12,7 +12,7 @@ use nalgebra_sparse::CscMatrix; use super::lin_expr::{LinExpr, QuadExpr}; use crate::expr::{Array, Expr, ExprId, IndexSpec, Shape, VariableBuilder}; -use crate::sparse::{csc_repeat_rows, csc_to_dense, csc_vstack, dense_to_csc}; +use crate::sparse::{csc_add, csc_repeat_rows, csc_to_dense, csc_vstack, dense_to_csc}; /// A cone constraint in standard form: Ax + b in K. #[derive(Debug, Clone)] @@ -903,15 +903,43 @@ impl CanonContext { fn canonicalize_sum_squares_lin(&mut self, x: &LinExpr, for_objective: bool) -> CanonExpr { if for_objective { - // For objective, use native QP: ||x||^2 = x' I x - // The (1/2) factor for Clarabel is handled in stuffing.rs + // For objective, use native QP: ||Ax + c||^2 = x'(A'A)x + 2c'Ax + ||c||^2 + // stuffing.rs doubles P to account for Clarabel's (1/2)x'Px convention. let vars = x.variables(); - if vars.len() == 1 && x.constant.iter().all(|&v| v == 0.0) { - let var_id = vars[0]; - let size = x.size(); - let identity = CscMatrix::identity(size); - return CanonExpr::Quadratic(QuadExpr::quadratic(var_id, identity)); + let c = &x.constant; // dense (m, 1) constant + + let mut quad_coeffs = std::collections::HashMap::new(); + let mut linear_coeffs = std::collections::HashMap::new(); + + for &var_i in &vars { + let ai = csc_to_dense(&x.coeffs[&var_i]); // (m, ni) + for &var_j in &vars { + let aj = csc_to_dense(&x.coeffs[&var_j]); // (m, nj) + // A_i' * A_j: (ni, m) * (m, nj) = (ni, nj) + let ai_t_aj = dense_to_csc(&(ai.transpose() * &aj)); + quad_coeffs + .entry((var_i, var_j)) + .and_modify(|existing| *existing = csc_add(existing, &ai_t_aj)) + .or_insert(ai_t_aj); + } + // Linear term: 2 * c' * A_i → (1, ni) row coefficient + let q_col = ai.transpose() * c; // (ni, 1) + let q_row = dense_to_csc(&(q_col * 2.0).transpose()); // (1, ni) + linear_coeffs.insert(var_i, q_row); } + + // Constant: ||c||^2 + let constant: f64 = c.iter().map(|v| v * v).sum(); + + return CanonExpr::Quadratic(QuadExpr { + quad_coeffs, + linear: LinExpr { + coeffs: linear_coeffs, + constant: DMatrix::zeros(1, 1), + shape: Shape::scalar(), + }, + constant, + }); } // SOC reformulation: ||x||^2 <= t iff SOC(sqrt(t), x) diff --git a/tests/sum_squares_tests.rs b/tests/sum_squares_tests.rs new file mode 100644 index 0000000..7501986 --- /dev/null +++ b/tests/sum_squares_tests.rs @@ -0,0 +1,80 @@ +use cvxrust::prelude::*; + +/// sum_squares(Ax - b) must be canonicalized as ||r||^2, not ||r||_2. +/// The old bug used a plain SOC `||r||_2 <= t`, which minimized the L2 norm +/// instead of its square, causing solution.value to be off by a square root. + +#[test] +fn test_sum_squares_scalar_constrained() { + // minimize (x - 3)^2 s.t. x <= 1 + // True optimum: x* = 1, obj* = (1 - 3)^2 = 4 + // Old bug would report: |1 - 3| = 2 + let x = variable(()); + let residual = x.clone() - constant(3.0); + + let sol = Problem::minimize(sum_squares(&residual)) + .constraint(constraint!(x <= 1.0)) + .solve() + .unwrap(); + + let reported = sol.value.unwrap(); + let eval_sq = sum_squares(&residual).value(&sol).as_scalar().unwrap(); + let eval_norm = norm2(&residual).value(&sol).as_scalar().unwrap(); + + assert!((reported - 4.0).abs() < 1e-4, "objective should be 4.0, got {reported}"); + assert!((reported - eval_sq).abs() < 1e-4, "reported obj must equal sum_squares evaluated at solution"); + // Sanity check: the L2 norm (sqrt(4) = 2) is clearly different from the correct answer + assert!((eval_norm - 2.0).abs() < 1e-4); + assert!((reported - eval_norm).abs() > 0.5, "objective must not equal the L2 norm (old bug)"); +} + +#[test] +fn test_sum_squares_vector_constrained() { + // minimize ||x - [3, 4]||^2 s.t. x <= 2 + // True optimum: x* = [2, 2], obj* = (2-3)^2 + (2-4)^2 = 5 + // Old bug would report: sqrt(5) ≈ 2.236 + let x = variable(2); + let residual = x.clone() - constant_vec(vec![3.0, 4.0]); + + let sol = Problem::minimize(sum_squares(&residual)) + .constraint(constraint!(x <= 2.0)) + .solve() + .unwrap(); + + let reported = sol.value.unwrap(); + let eval_sq = sum_squares(&residual).value(&sol).as_scalar().unwrap(); + let eval_norm = norm2(&residual).value(&sol).as_scalar().unwrap(); + + assert!((reported - 5.0).abs() < 1e-4, "objective should be 5.0, got {reported}"); + assert!((reported - eval_sq).abs() < 1e-4, "reported obj must equal sum_squares evaluated at solution"); + // Sanity check: the L2 norm (sqrt(5) ≈ 2.236) is clearly different + assert!((eval_norm - 5f64.sqrt()).abs() < 1e-4); + assert!((reported - eval_norm).abs() > 0.5, "objective must not equal the L2 norm (old bug)"); +} + +#[test] +fn test_sum_squares_matmul_constrained() { + // minimize ||Ax - b||^2 s.t. x <= 1 + // A = [[1], [1]] (2x1), b = [2, 4], x scalar + // Unconstrained LS: x* = 3, obj* = (3-2)^2 + (3-4)^2 = 2 + // With x <= 1: x* = 1, residual = [-1, -3], obj* = 1 + 9 = 10 + // Old bug would report: sqrt(10) ≈ 3.162 + let a = constant_matrix(vec![1.0, 1.0], 2, 1); + let b = constant_vec(vec![2.0, 4.0]); + let x = variable(()); + let residual = matmul(&a, &x) - &b; + + let sol = Problem::minimize(sum_squares(&residual)) + .constraint(constraint!(x <= 1.0)) + .solve() + .unwrap(); + + let reported = sol.value.unwrap(); + let eval_sq = sum_squares(&residual).value(&sol).as_scalar().unwrap(); + let eval_norm = norm2(&residual).value(&sol).as_scalar().unwrap(); + + assert!((reported - 10.0).abs() < 1e-4, "objective should be 10.0, got {reported}"); + assert!((reported - eval_sq).abs() < 1e-4, "reported obj must equal sum_squares evaluated at solution"); + assert!((eval_norm - 10f64.sqrt()).abs() < 1e-4); + assert!((reported - eval_norm).abs() > 0.5, "objective must not equal the L2 norm (old bug)"); +} From 2c672bfee1b66be0c96f24eb52d6101eb0a3ace3 Mon Sep 17 00:00:00 2001 From: Hao Zhu Date: Mon, 9 Mar 2026 00:06:20 +0100 Subject: [PATCH 2/3] Update the name of MSRV check to version 1.71 in CI configuration --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index af58674..d789ee9 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -105,7 +105,7 @@ jobs: run: cargo run --example portfolio msrv: - name: Check MSRV (1.70) + name: Check MSRV (1.71) runs-on: ubuntu-latest steps: From e038daf80197840ac7cbce1418dbed60d77eed90 Mon Sep 17 00:00:00 2001 From: Hao Zhu Date: Mon, 9 Mar 2026 00:13:57 +0100 Subject: [PATCH 3/3] Fix lint --- src/canon/canonicalizer.rs | 2 +- tests/sum_squares_tests.rs | 45 ++++++++++++++++++++++++++++++-------- 2 files changed, 37 insertions(+), 10 deletions(-) diff --git a/src/canon/canonicalizer.rs b/src/canon/canonicalizer.rs index 718e13b..bcecb15 100644 --- a/src/canon/canonicalizer.rs +++ b/src/canon/canonicalizer.rs @@ -915,7 +915,7 @@ impl CanonContext { let ai = csc_to_dense(&x.coeffs[&var_i]); // (m, ni) for &var_j in &vars { let aj = csc_to_dense(&x.coeffs[&var_j]); // (m, nj) - // A_i' * A_j: (ni, m) * (m, nj) = (ni, nj) + // A_i' * A_j: (ni, m) * (m, nj) = (ni, nj) let ai_t_aj = dense_to_csc(&(ai.transpose() * &aj)); quad_coeffs .entry((var_i, var_j)) diff --git a/tests/sum_squares_tests.rs b/tests/sum_squares_tests.rs index 7501986..c1b59e1 100644 --- a/tests/sum_squares_tests.rs +++ b/tests/sum_squares_tests.rs @@ -21,11 +21,20 @@ fn test_sum_squares_scalar_constrained() { let eval_sq = sum_squares(&residual).value(&sol).as_scalar().unwrap(); let eval_norm = norm2(&residual).value(&sol).as_scalar().unwrap(); - assert!((reported - 4.0).abs() < 1e-4, "objective should be 4.0, got {reported}"); - assert!((reported - eval_sq).abs() < 1e-4, "reported obj must equal sum_squares evaluated at solution"); + assert!( + (reported - 4.0).abs() < 1e-4, + "objective should be 4.0, got {reported}" + ); + assert!( + (reported - eval_sq).abs() < 1e-4, + "reported obj must equal sum_squares evaluated at solution" + ); // Sanity check: the L2 norm (sqrt(4) = 2) is clearly different from the correct answer assert!((eval_norm - 2.0).abs() < 1e-4); - assert!((reported - eval_norm).abs() > 0.5, "objective must not equal the L2 norm (old bug)"); + assert!( + (reported - eval_norm).abs() > 0.5, + "objective must not equal the L2 norm (old bug)" + ); } #[test] @@ -45,11 +54,20 @@ fn test_sum_squares_vector_constrained() { let eval_sq = sum_squares(&residual).value(&sol).as_scalar().unwrap(); let eval_norm = norm2(&residual).value(&sol).as_scalar().unwrap(); - assert!((reported - 5.0).abs() < 1e-4, "objective should be 5.0, got {reported}"); - assert!((reported - eval_sq).abs() < 1e-4, "reported obj must equal sum_squares evaluated at solution"); + assert!( + (reported - 5.0).abs() < 1e-4, + "objective should be 5.0, got {reported}" + ); + assert!( + (reported - eval_sq).abs() < 1e-4, + "reported obj must equal sum_squares evaluated at solution" + ); // Sanity check: the L2 norm (sqrt(5) ≈ 2.236) is clearly different assert!((eval_norm - 5f64.sqrt()).abs() < 1e-4); - assert!((reported - eval_norm).abs() > 0.5, "objective must not equal the L2 norm (old bug)"); + assert!( + (reported - eval_norm).abs() > 0.5, + "objective must not equal the L2 norm (old bug)" + ); } #[test] @@ -73,8 +91,17 @@ fn test_sum_squares_matmul_constrained() { let eval_sq = sum_squares(&residual).value(&sol).as_scalar().unwrap(); let eval_norm = norm2(&residual).value(&sol).as_scalar().unwrap(); - assert!((reported - 10.0).abs() < 1e-4, "objective should be 10.0, got {reported}"); - assert!((reported - eval_sq).abs() < 1e-4, "reported obj must equal sum_squares evaluated at solution"); + assert!( + (reported - 10.0).abs() < 1e-4, + "objective should be 10.0, got {reported}" + ); + assert!( + (reported - eval_sq).abs() < 1e-4, + "reported obj must equal sum_squares evaluated at solution" + ); assert!((eval_norm - 10f64.sqrt()).abs() < 1e-4); - assert!((reported - eval_norm).abs() > 0.5, "objective must not equal the L2 norm (old bug)"); + assert!( + (reported - eval_norm).abs() > 0.5, + "objective must not equal the L2 norm (old bug)" + ); }