Fix #91: [Rule] ClosestVectorProblem to QUBO

docs/paper/reductions.typ

@@ -4336,6 +4336,65 @@ Each reduction is presented as a *Rule* (with linked problem names and overhead
   _Solution extraction._ Convert binary to spins: $s_i = 2x_i - 1$, i.e.\ $x_i = 1 arrow.r s_i = +1$, $x_i = 0 arrow.r s_i = -1$.
 ]
 
+#let cvp_qubo = load-example("ClosestVectorProblem", "QUBO")
+#let cvp_qubo_sol = cvp_qubo.solutions.at(0)
+#{
+  let basis = cvp_qubo.source.instance.basis
+  let bounds = cvp_qubo.source.instance.bounds
+  let target = cvp_qubo.source.instance.target
+  let offsets = cvp_qubo_sol.source_config
+  let coords = offsets.enumerate().map(((i, off)) => off + bounds.at(i).lower)
+  let matrix = cvp_qubo.target.instance.matrix
+  let bits = cvp_qubo_sol.target_config
+  let lo = bounds.map(b => b.lower)
+  let anchor = range(target.len()).map(d => lo.enumerate().fold(0.0, (acc, (i, x)) => acc + x * basis.at(i).at(d)))
+  let constant = range(target.len()).fold(0.0, (acc, d) => acc + calc.pow(anchor.at(d) - target.at(d), 2))
+  let qubo-value = range(bits.len()).fold(0.0, (acc, i) => acc + if bits.at(i) == 0 { 0.0 } else {
+    range(bits.len() - i).fold(0.0, (row-acc, delta) => row-acc + if bits.at(i + delta) == 0 { 0.0 } else { matrix.at(i).at(i + delta) })
+  })
+  let fmt-vec(v) = $paren.l #v.map(e => str(e)).join(", ") paren.r^top$
+  let rounded-constant = calc.round(constant, digits: 2)
+  let rounded-qubo = calc.round(qubo-value, digits: 1)
+  let rounded-distance-sq = calc.round(qubo-value + constant, digits: 2)
+  [
+    #reduction-rule("ClosestVectorProblem", "QUBO",
+      example: true,
+      example-caption: [2D bounded CVP with two 3-bit exact-range encodings],
+      extra: [
+        *Step 1 -- Source instance.* The canonical CVP example uses basis columns $bold(b)_1 = #fmt-vec(basis.at(0))$ and $bold(b)_2 = #fmt-vec(basis.at(1))$, target $bold(t) = #fmt-vec(target)$, and bounds $x_1, x_2 in [#bounds.at(0).lower, #bounds.at(0).upper]$.
+
+        *Step 2 -- Exact bounded encoding.* Each variable has #bounds.at(0).upper - bounds.at(0).lower + 1 admissible values, so the implementation uses the capped binary basis $(1, 2, 3)$ rather than $(1, 2, 4)$: the first two bits are powers of two, and the last weight is capped so every bit pattern reconstructs an offset in ${0, dots, 6}$. Thus
+        $ x_1 = #bounds.at(0).lower + z_0 + 2 z_1 + 3 z_2, quad x_2 = #bounds.at(1).lower + z_3 + 2 z_4 + 3 z_5 $
+        giving #cvp_qubo.target.instance.num_vars QUBO variables in total.
+
+        *Step 3 -- Build the QUBO.* For this instance, $G = A^top A = ((4, 2), (2, 5))$ and $h = A^top bold(t) = (5.6, 5.8)^top$. Expanding the shifted quadratic form yields the exported upper-triangular matrix with representative entries $Q_(0,0) = #matrix.at(0).at(0)$, $Q_(0,1) = #matrix.at(0).at(1)$, $Q_(0,2) = #matrix.at(0).at(2)$, $Q_(2,5) = #matrix.at(2).at(5)$, and $Q_(5,5) = #matrix.at(5).at(5)$.
+
+        *Step 4 -- Verify a solution.* The fixture stores the canonical witness $bold(z) = (#bits.map(str).join(", "))$, which extracts to source offsets $bold(c) = (#offsets.map(str).join(", "))$ and actual lattice coordinates $bold(x) = (#coords.map(str).join(", "))$. The QUBO value is $bold(z)^top Q bold(z) = #rounded-qubo$; adding back the dropped constant #rounded-constant yields the original squared distance #(rounded-distance-sq), so the extracted point is the closest lattice vector #sym.checkmark.
+
+        *Multiplicity.* Offset $3$ has two bit encodings ($(0, 0, 1)$ and $(1, 1, 0)$), so the fixture stores one canonical witness even though the QUBO has multiple optimal binary assignments representing the same CVP solution.
+      ],
+    )[
+      A bounded Closest Vector Problem instance already supplies a finite integer box $x_i in [ell_i, u_i]$ for each coefficient. Following the direct quadratic-form reduction of Canale, Qureshi, and Viola @canale2023qubo, encoding each offset $c_i = x_i - ell_i$ with an exact in-range binary basis turns the squared-distance objective into an unconstrained quadratic over binary variables. Unlike penalty-method encodings, no auxiliary feasibility penalty is needed: every bit pattern decodes to a legal coefficient vector by construction.
+    ][
+      _Construction._ Let $A in ZZ^(m times n)$ be the basis matrix with columns $bold(a)_1, dots, bold(a)_n$, let $bold(t) in RR^m$ be the target, and let $x_i in [ell_i, u_i]$ with range $r_i = u_i - ell_i$. Define $L_i = ceil(log_2(r_i + 1))$ when $r_i > 0$ and omit bits when $r_i = 0$. For each variable, introduce binary variables $z_(i,0), dots, z_(i,L_i-1)$ with exact-range weights
+      $ w_(i,p) = 2^p quad (0 <= p < L_i - 1), quad w_(i,L_i-1) = r_i + 1 - 2^(L_i - 1) $
+      so that every bit vector represents an offset in ${0, dots, r_i}$. Then
+      $ x_i = ell_i + sum_(p=0)^(L_i-1) w_(i,p) z_(i,p) $
+      and the total number of QUBO variables is $N = sum_i L_i$, exactly the exported overhead `num_vars = num_encoding_bits`.
+
+      Let $G = A^top A$ and $h = A^top bold(t)$. Writing $bold(x) = bold(ell) + B bold(z)$ for the encoding matrix $B in RR^(n times N)$ gives
+      $ norm(A bold(x) - bold(t))_2^2 = bold(z)^top (B^top G B) bold(z) + 2 bold(z)^top B^top (G bold(ell) - h) + "const" $
+      where the constant $norm(A bold(ell) - bold(t))_2^2$ is dropped. Therefore the QUBO coefficients are
+      $ Q_(u,u) = (B^top G B)_(u,u) + 2 (B^top (G bold(ell) - h))_u, quad Q_(u,v) = 2 (B^top G B)_(u,v) quad (u < v) $
+      using the usual upper-triangular convention.
+
+      _Correctness._ ($arrow.r.double$) Every binary vector $bold(z) in {0,1}^N$ decodes to a coefficient vector $bold(x)$ inside the prescribed bounds because each exact-range basis reaches only offsets in ${0, dots, r_i}$. Substituting this decoding into the CVP objective yields $bold(z)^top Q bold(z) + "const"$, so any QUBO minimizer maps to a bounded CVP minimizer. ($arrow.l.double$) Every bounded CVP solution $bold(x)$ has at least one bit encoding for each coordinate offset, hence at least one binary vector $bold(z)$ with the same objective value up to the dropped constant. Thus the minimizers correspond exactly, although several binary witnesses may decode to the same CVP solution.
+
+      _Solution extraction._ For each source variable, sum its selected encoding weights to recover the source configuration offset $c_i = x_i - ell_i$. This is exactly the configuration format expected by the `ClosestVectorProblem` model.
+    ]
+  ]
+}
+
 == Penalty-Method QUBO Reductions <sec:penalty-method>
 
 The _penalty method_ @glover2019 @lucas2014 converts a constrained optimization problem into an unconstrained QUBO by adding quadratic penalty terms. Given an objective $"obj"(bold(x))$ to minimize and constraints $g_k (bold(x)) = 0$, construct:
@@ -5250,6 +5309,7 @@ The following table shows concrete variable overhead for example instances, take
   (source: "SpinGlass", target: "MaxCut"),
   (source: "SpinGlass", target: "QUBO"),
   (source: "QUBO", target: "SpinGlass"),
+  (source: "ClosestVectorProblem", target: "QUBO"),
   (source: "KColoring", target: "QUBO"),
   (source: "MaximumSetPacking", target: "QUBO"),
   (

docs/paper/references.bib

@@ -507,6 +507,15 @@ @article{pan2025
   archivePrefix = {arXiv}
 }
 
+@article{canale2023qubo,
+  author  = {Eduardo Canale and Claudio Qureshi and Alfredo Viola},
+  title   = {Qubo model for the Closest Vector Problem},
+  journal = {arXiv preprint},
+  year    = {2023},
+  eprint  = {2304.03616},
+  archivePrefix = {arXiv}
+}
+
 @article{goemans1995,
   author  = {Michel X. Goemans and David P. Williamson},
   title   = {Improved Approximation Algorithms for Maximum Cut and Satisfiability Problems Using Semidefinite Programming},
@@ -793,17 +802,6 @@ @article{cygan2014
   doi     = {10.1137/140990255}
 }
 
-@article{dreyfuswagner1971,
-  author  = {Stuart E. Dreyfus and Robert A. Wagner},
-  title   = {The Steiner Problem in Graphs},
-  journal = {Networks},
-  volume  = {1},
-  number  = {3},
-  pages   = {195--207},
-  year    = {1971},
-  doi     = {10.1002/net.3230010302}
-}
-
 @inproceedings{bjorklund2007,
   author    = {Andreas Bj\"{o}rklund and Thore Husfeldt and Petteri Kaski and Mikko Koivisto},
   title     = {Fourier Meets M\"{o}bius: Fast Subset Convolution},

src/models/algebraic/closest_vector_problem.rs

@@ -98,6 +98,42 @@ impl VarBounds {
             _ => None,
         }
     }
+
+    /// Returns an exact bounded binary basis for offsets in this range.
+    ///
+    /// For a bounded variable with offsets `0..=hi-lo`, the returned weights
+    /// ensure that every bit-pattern reconstructs an in-range offset. Low-order
+    /// weights use powers of two; the final weight is capped so the maximum
+    /// reachable offset is exactly `hi-lo`.
+    pub(crate) fn exact_encoding_weights(&self) -> Vec<i64> {
+        let Some(num_values) = self.num_values() else {
+            panic!("CVP QUBO encoding requires finite variable bounds");
+        };
+        if num_values <= 1 {
+            return Vec::new();
+        }
+
+        let max_offset = (num_values - 1) as i64;
+        let num_bits = (usize::BITS - (num_values - 1).leading_zeros()) as usize;
+        let mut weights = Vec::with_capacity(num_bits);
+
+        for bit in 0..num_bits.saturating_sub(1) {
+            weights.push(1_i64 << bit);
+        }
+
+        let covered_by_lower_bits = if num_bits <= 1 {
+            0
+        } else {
+            (1_i64 << (num_bits - 1)) - 1
+        };
+        weights.push(max_offset - covered_by_lower_bits);
+        weights
+    }
+
+    /// Returns the number of encoding bits needed for the exact bounded basis.
+    pub(crate) fn num_encoding_bits(&self) -> usize {
+        self.exact_encoding_weights().len()
+    }
 }
 
 /// Closest Vector Problem (CVP).
@@ -175,6 +211,11 @@ impl<T> ClosestVectorProblem<T> {
         &self.bounds
     }
 
+    /// Returns the total number of bounded-encoding bits used by the QUBO form.
+    pub fn num_encoding_bits(&self) -> usize {
+        self.bounds.iter().map(VarBounds::num_encoding_bits).sum()
+    }
+
     /// Convert a configuration (offsets from lower bounds) to integer values.
     fn config_to_values(&self, config: &[usize]) -> Vec<i64> {
         config

src/rules/closestvectorproblem_qubo.rs

@@ -0,0 +1,197 @@
+//! Reduction from ClosestVectorProblem to QUBO.
+//!
+//! Encodes each bounded CVP coefficient with an exact in-range binary basis and
+//! expands the squared-distance objective into a QUBO over those bits.
+
+#[cfg(feature = "example-db")]
+use crate::export::SolutionPair;
+use crate::models::algebraic::{ClosestVectorProblem, QUBO};
+use crate::reduction;
+use crate::rules::traits::{ReduceTo, ReductionResult};
+
+#[derive(Debug, Clone)]
+struct EncodingSpan {
+    start: usize,
+    weights: Vec<usize>,
+}
+
+/// Result of reducing a bounded ClosestVectorProblem instance to QUBO.
+#[derive(Debug, Clone)]
+pub struct ReductionCVPToQUBO {
+    target: QUBO<f64>,
+    encodings: Vec<EncodingSpan>,
+}
+
+impl ReductionResult for ReductionCVPToQUBO {
+    type Source = ClosestVectorProblem<i32>;
+    type Target = QUBO<f64>;
+
+    fn target_problem(&self) -> &Self::Target {
+        &self.target
+    }
+
+    /// Reconstruct the source configuration offsets from the encoded QUBO bits.
+    fn extract_solution(&self, target_solution: &[usize]) -> Vec<usize> {
+        self.encodings
+            .iter()
+            .map(|encoding| {
+                encoding
+                    .weights
+                    .iter()
+                    .enumerate()
+                    .map(|(offset, weight)| {
+                        target_solution
+                            .get(encoding.start + offset)
+                            .copied()
+                            .unwrap_or(0)
+                            * weight
+                    })
+                    .sum()
+            })
+            .collect()
+    }
+}
+
+#[cfg(feature = "example-db")]
+fn canonical_cvp_instance() -> ClosestVectorProblem<i32> {
+    ClosestVectorProblem::new(
+        vec![vec![2, 0], vec![1, 2]],
+        vec![2.8, 1.5],
+        vec![
+            crate::models::algebraic::VarBounds::bounded(-2, 4),
+            crate::models::algebraic::VarBounds::bounded(-2, 4),
+        ],
+    )
+}
+
+fn encoding_spans(problem: &ClosestVectorProblem<i32>) -> Vec<EncodingSpan> {
+    let mut start = 0usize;
+    let mut spans = Vec::with_capacity(problem.num_basis_vectors());
+    for bounds in problem.bounds() {
+        let weights = bounds
+            .exact_encoding_weights()
+            .into_iter()
+            .map(|weight| usize::try_from(weight).expect("encoding weights must be nonnegative"))
+            .collect::<Vec<_>>();
+        spans.push(EncodingSpan { start, weights });
+        start += spans.last().expect("just pushed").weights.len();
+    }
+    spans
+}
+
+fn gram_matrix(problem: &ClosestVectorProblem<i32>) -> Vec<Vec<f64>> {
+    let basis = problem.basis();
+    let n = basis.len();
+    let mut gram = vec![vec![0.0; n]; n];
+    for i in 0..n {
+        for j in i..n {
+            let dot = basis[i]
+                .iter()
+                .zip(&basis[j])
+                .map(|(&lhs, &rhs)| lhs as f64 * rhs as f64)
+                .sum::<f64>();
+            gram[i][j] = dot;
+            gram[j][i] = dot;
+        }
+    }
+    gram
+}
+
+fn at_times_target(problem: &ClosestVectorProblem<i32>) -> Vec<f64> {
+    problem
+        .basis()
+        .iter()
+        .map(|column| {
+            column
+                .iter()
+                .zip(problem.target())
+                .map(|(&entry, &target)| entry as f64 * target)
+                .sum()
+        })
+        .collect()
+}
+
+#[reduction(overhead = { num_vars = "num_encoding_bits" })]
+impl ReduceTo<QUBO<f64>> for ClosestVectorProblem<i32> {
+    type Result = ReductionCVPToQUBO;
+
+    fn reduce_to(&self) -> Self::Result {
+        let encodings = encoding_spans(self);
+        let total_bits = encodings
+            .last()
+            .map(|encoding| encoding.start + encoding.weights.len())
+            .unwrap_or(0);
+        let mut matrix = vec![vec![0.0; total_bits]; total_bits];
+
+        if total_bits == 0 {
+            return ReductionCVPToQUBO {
+                target: QUBO::from_matrix(matrix),
+                encodings,
+            };
+        }
+
+        let gram = gram_matrix(self);
+        let h = at_times_target(self);
+        let lowers = self
+            .bounds()
+            .iter()
+            .map(|bounds| {
+                bounds
+                    .lower
+                    .expect("CVP QUBO reduction requires finite lower bounds")
+            })
+            .map(|lower| lower as f64)
+            .collect::<Vec<_>>();
+        let g_lo_minus_h = (0..self.num_basis_vectors())
+            .map(|i| {
+                (0..self.num_basis_vectors())
+                    .map(|j| gram[i][j] * lowers[j])
+                    .sum::<f64>()
+                    - h[i]
+            })
+            .collect::<Vec<_>>();
+
+        let mut bit_terms = Vec::with_capacity(total_bits);
+        for (var_index, encoding) in encodings.iter().enumerate() {
+            for &weight in &encoding.weights {
+                bit_terms.push((var_index, weight as f64));
+            }
+        }
+
+        for u in 0..total_bits {
+            let (var_u, weight_u) = bit_terms[u];
+            matrix[u][u] =
+                gram[var_u][var_u] * weight_u * weight_u + 2.0 * weight_u * g_lo_minus_h[var_u];
+
+            for v in (u + 1)..total_bits {
+                let (var_v, weight_v) = bit_terms[v];
+                matrix[u][v] = 2.0 * gram[var_u][var_v] * weight_u * weight_v;
+            }
+        }
+
+        ReductionCVPToQUBO {
+            target: QUBO::from_matrix(matrix),
+            encodings,
+        }
+    }
+}
+
+#[cfg(feature = "example-db")]
+pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
+    vec![crate::example_db::specs::RuleExampleSpec {
+        id: "closestvectorproblem_to_qubo",
+        build: || {
+            crate::example_db::specs::rule_example_with_witness::<_, QUBO<f64>>(
+                canonical_cvp_instance(),
+                SolutionPair {
+                    source_config: vec![3, 3],
+                    target_config: vec![0, 0, 1, 0, 0, 1],
+                },
+            )
+        },
+    }]
+}
+
+#[cfg(test)]
+#[path = "../unit_tests/rules/closestvectorproblem_qubo.rs"]
+mod tests;

src/rules/mod.rs

@@ -7,6 +7,7 @@ pub use cost::{CustomCost, Minimize, MinimizeSteps, PathCostFn};
 pub use registry::{ReductionEntry, ReductionOverhead};
 
 pub(crate) mod circuit_spinglass;
+mod closestvectorproblem_qubo;
 pub(crate) mod coloring_qubo;
 pub(crate) mod factoring_circuit;
 mod graph;
@@ -89,6 +90,7 @@ pub use traits::{ReduceTo, ReductionAutoCast, ReductionResult};
 pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
     let mut specs = Vec::new();
     specs.extend(circuit_spinglass::canonical_rule_example_specs());
+    specs.extend(closestvectorproblem_qubo::canonical_rule_example_specs());
     specs.extend(coloring_qubo::canonical_rule_example_specs());
     specs.extend(factoring_circuit::canonical_rule_example_specs());
     specs.extend(knapsack_qubo::canonical_rule_example_specs());