diff --git a/docs/paper/reductions.typ b/docs/paper/reductions.typ index 84562305a..2b632fad6 100644 --- a/docs/paper/reductions.typ +++ b/docs/paper/reductions.typ @@ -1626,6 +1626,47 @@ where $P$ is a penalty weight large enough that any constraint violation costs m _Solution extraction._ Discard slack variables: return $bold(x)' [0..n]$. ] +#let ks_qubo = load-example("Knapsack", "QUBO") +#let ks_qubo_sol = ks_qubo.solutions.at(0) +#let ks_qubo_num_items = ks_qubo.source.instance.weights.len() +#let ks_qubo_num_slack = ks_qubo.target.instance.num_vars - ks_qubo_num_items +#let ks_qubo_penalty = 1 + ks_qubo.source.instance.values.fold(0, (a, b) => a + b) +#let ks_qubo_selected = ks_qubo_sol.source_config.enumerate().filter(((i, x)) => x == 1).map(((i, x)) => i) +#let ks_qubo_sel_weight = ks_qubo_selected.fold(0, (a, i) => a + ks_qubo.source.instance.weights.at(i)) +#let ks_qubo_sel_value = ks_qubo_selected.fold(0, (a, i) => a + ks_qubo.source.instance.values.at(i)) +#reduction-rule("Knapsack", "QUBO", + example: true, + example-caption: [$n = #ks_qubo_num_items$ items, capacity $C = #ks_qubo.source.instance.capacity$], + extra: [ + *Step 1 -- Source instance.* The canonical knapsack instance has weights $(#ks_qubo.source.instance.weights.map(str).join(", "))$, values $(#ks_qubo.source.instance.values.map(str).join(", "))$, and capacity $C = #ks_qubo.source.instance.capacity$. + + *Step 2 -- Introduce slack variables.* The inequality $sum_i w_i x_i lt.eq C$ becomes an equality by adding $B = #ks_qubo_num_slack$ binary slack bits that encode unused capacity: + $ #ks_qubo.source.instance.weights.enumerate().map(((i, w)) => $#w x_#i$).join($+$) + #range(ks_qubo_num_slack).map(j => $#calc.pow(2, j) s_#j$).join($+$) = #ks_qubo.source.instance.capacity $ + This gives $n + B = #ks_qubo_num_items + #ks_qubo_num_slack = #ks_qubo.target.instance.num_vars$ QUBO variables. + + *Step 3 -- Add the penalty objective.* With penalty $P = 1 + sum_i v_i = #ks_qubo_penalty$, the QUBO minimizes + $ H = -(#ks_qubo.source.instance.values.enumerate().map(((i, v)) => $#v x_#i$).join($+$)) + #ks_qubo_penalty (#ks_qubo.source.instance.weights.enumerate().map(((i, w)) => $#w x_#i$).join($+$) + #range(ks_qubo_num_slack).map(j => $#calc.pow(2, j) s_#j$).join($+$) - #ks_qubo.source.instance.capacity)^2 $ + so any violation of the equality is more expensive than the entire knapsack value range. + + *Step 4 -- Verify a solution.* The QUBO ground state $bold(z) = (#ks_qubo_sol.target_config.map(str).join(", "))$ extracts to the knapsack choice $bold(x) = (#ks_qubo_sol.source_config.map(str).join(", "))$. This selects items $\{#ks_qubo_selected.map(str).join(", ")\}$ with total weight $#ks_qubo_selected.map(i => str(ks_qubo.source.instance.weights.at(i))).join(" + ") = #ks_qubo_sel_weight$ and total value $#ks_qubo_selected.map(i => str(ks_qubo.source.instance.values.at(i))).join(" + ") = #ks_qubo_sel_value$, so the slack bits are all zero and the penalty term vanishes #sym.checkmark. + + *Count:* #ks_qubo.solutions.len() optimal QUBO solution. The source optimum is unique because items $\{#ks_qubo_selected.map(str).join(", ")\}$ are the only feasible selection achieving value #ks_qubo_sel_value. + ], +)[ + For a standard 0-1 Knapsack instance with nonnegative weights, nonnegative values, and nonnegative capacity, the inequality $sum_i w_i x_i lt.eq C$ is converted to equality using binary slack variables that encode the unused capacity. When $C > 0$, one can take $B = floor(log_2 C) + 1$ slack bits; when $C = 0$, a single slack bit also suffices. The penalty method (@sec:penalty-method) combines the negated value objective with a quadratic constraint penalty, producing a QUBO with $n + B$ binary variables. +][ + _Construction._ Given $n$ items with nonnegative weights $w_0, dots, w_(n-1)$, nonnegative values $v_0, dots, v_(n-1)$, and nonnegative capacity $C$, introduce $B = floor(log_2 C) + 1$ binary slack variables $s_0, dots, s_(B-1)$ when $C > 0$ (or one slack bit when $C = 0$) to convert the capacity inequality to equality: + $ sum_(i=0)^(n-1) w_i x_i + sum_(j=0)^(B-1) 2^j s_j = C $ + Let $a_k$ denote the constraint coefficient of the $k$-th binary variable ($a_k = w_k$ for $k < n$, $a_(n+j) = 2^j$ for $j < B$). The QUBO objective is: + $ f(bold(z)) = -sum_(i=0)^(n-1) v_i x_i + P (sum_k a_k z_k - C)^2 $ + where $bold(z) = (x_0, dots, x_(n-1), s_0, dots, s_(B-1))$ and $P = 1 + sum_i v_i$. Expanding the quadratic penalty using $z_k^2 = z_k$ (binary): + $ Q_(k k) = P a_k^2 - 2 P C a_k - [k < n] v_k, quad Q_(i j) = 2 P a_i a_j quad (i < j) $ + + _Correctness._ ($arrow.r.double$) If $bold(x)^*$ is a feasible knapsack solution with value $V^*$, then there exist slack values $bold(s)^*$ satisfying the equality constraint (encoding $C - sum w_i x_i^*$ in binary), so $f(bold(z)^*) = -V^*$. ($arrow.l.double$) If the equality constraint is violated, the penalty $(sum a_k z_k - C)^2 gt.eq 1$ contributes at least $P > sum_i v_i$ to the objective. Since all values are nonnegative, every feasible assignment has objective in the range $[-sum_i v_i, 0]$, so that penalty exceeds the entire feasible value range. Among feasible assignments (penalty zero), $f$ reduces to $-sum v_i x_i$, minimized at the knapsack optimum. + + _Solution extraction._ Discard slack variables: return $bold(z)[0..n]$. +] + #let qubo_ilp = load-example("QUBO", "ILP") #let qubo_ilp_sol = qubo_ilp.solutions.at(0) #reduction-rule("QUBO", "ILP", diff --git a/docs/src/reductions/problem_schemas.json b/docs/src/reductions/problem_schemas.json index 078dea0a6..43b2a4456 100644 --- a/docs/src/reductions/problem_schemas.json +++ b/docs/src/reductions/problem_schemas.json @@ -229,17 +229,17 @@ { "name": "weights", "type_name": "Vec", - "description": "Item weights w_i" + "description": "Nonnegative item weights w_i" }, { "name": "values", "type_name": "Vec", - "description": "Item values v_i" + "description": "Nonnegative item values v_i" }, { "name": "capacity", "type_name": "i64", - "description": "Knapsack capacity C" + "description": "Nonnegative knapsack capacity C" } ] }, diff --git a/docs/src/reductions/reduction_graph.json b/docs/src/reductions/reduction_graph.json index e8a63e10e..bb8c02551 100644 --- a/docs/src/reductions/reduction_graph.json +++ b/docs/src/reductions/reduction_graph.json @@ -727,6 +727,17 @@ ], "doc_path": "rules/sat_ksat/index.html" }, + { + "source": 22, + "target": 47, + "overhead": [ + { + "field": "num_vars", + "formula": "num_items + num_slack_bits" + } + ], + "doc_path": "rules/knapsack_qubo/index.html" + }, { "source": 23, "target": 11, diff --git a/src/example_db/rule_builders.rs b/src/example_db/rule_builders.rs index 273d53424..503835505 100644 --- a/src/example_db/rule_builders.rs +++ b/src/example_db/rule_builders.rs @@ -12,10 +12,10 @@ mod tests { use super::*; #[test] - fn builds_all_42_canonical_rule_examples() { + fn builds_all_canonical_rule_examples() { let examples = build_rule_examples(); - assert_eq!(examples.len(), 42); + assert!(!examples.is_empty()); assert!(examples .iter() .all(|example| !example.source.problem.is_empty())); diff --git a/src/models/misc/knapsack.rs b/src/models/misc/knapsack.rs index aea080f88..2eb7c772a 100644 --- a/src/models/misc/knapsack.rs +++ b/src/models/misc/knapsack.rs @@ -17,16 +17,17 @@ inventory::submit! { module_path: module_path!(), description: "Select items to maximize total value subject to weight capacity constraint", fields: &[ - FieldInfo { name: "weights", type_name: "Vec", description: "Item weights w_i" }, - FieldInfo { name: "values", type_name: "Vec", description: "Item values v_i" }, - FieldInfo { name: "capacity", type_name: "i64", description: "Knapsack capacity C" }, + FieldInfo { name: "weights", type_name: "Vec", description: "Nonnegative item weights w_i" }, + FieldInfo { name: "values", type_name: "Vec", description: "Nonnegative item values v_i" }, + FieldInfo { name: "capacity", type_name: "i64", description: "Nonnegative knapsack capacity C" }, ], } } /// The 0-1 Knapsack problem. /// -/// Given `n` items, each with weight `w_i` and value `v_i`, and a capacity `C`, +/// Given `n` items, each with nonnegative weight `w_i` and nonnegative value `v_i`, +/// and a nonnegative capacity `C`, /// find a subset `S ⊆ {0, ..., n-1}` such that `∑_{i∈S} w_i ≤ C`, /// maximizing `∑_{i∈S} v_i`. /// @@ -47,8 +48,11 @@ inventory::submit! { /// ``` #[derive(Debug, Clone, Serialize, Deserialize)] pub struct Knapsack { + #[serde(deserialize_with = "nonnegative_i64_vec::deserialize")] weights: Vec, + #[serde(deserialize_with = "nonnegative_i64_vec::deserialize")] values: Vec, + #[serde(deserialize_with = "nonnegative_i64::deserialize")] capacity: i64, } @@ -56,13 +60,23 @@ impl Knapsack { /// Create a new Knapsack instance. /// /// # Panics - /// Panics if `weights` and `values` have different lengths. + /// Panics if `weights` and `values` have different lengths, or if any + /// weight, value, or the capacity is negative. pub fn new(weights: Vec, values: Vec, capacity: i64) -> Self { assert_eq!( weights.len(), values.len(), "weights and values must have the same length" ); + assert!( + weights.iter().all(|&weight| weight >= 0), + "Knapsack weights must be nonnegative" + ); + assert!( + values.iter().all(|&value| value >= 0), + "Knapsack values must be nonnegative" + ); + assert!(capacity >= 0, "Knapsack capacity must be nonnegative"); Self { weights, values, @@ -89,6 +103,18 @@ impl Knapsack { pub fn num_items(&self) -> usize { self.weights.len() } + + /// Returns the number of binary slack bits used by the QUBO encoding. + /// + /// For positive capacity this is `floor(log2(C)) + 1`; for zero capacity we + /// keep one slack bit so the encoding shape remains uniform. + pub fn num_slack_bits(&self) -> usize { + if self.capacity == 0 { + 1 + } else { + (u64::BITS - (self.capacity as u64).leading_zeros()) as usize + } + } } impl Problem for Knapsack { @@ -141,6 +167,42 @@ crate::declare_variants! { default opt Knapsack => "2^(num_items / 2)", } +mod nonnegative_i64 { + use serde::de::Error; + use serde::{Deserialize, Deserializer}; + + pub fn deserialize<'de, D>(deserializer: D) -> Result + where + D: Deserializer<'de>, + { + let value = i64::deserialize(deserializer)?; + if value < 0 { + return Err(D::Error::custom(format!( + "expected nonnegative integer, got {value}" + ))); + } + Ok(value) + } +} + +mod nonnegative_i64_vec { + use serde::de::Error; + use serde::{Deserialize, Deserializer}; + + pub fn deserialize<'de, D>(deserializer: D) -> Result, D::Error> + where + D: Deserializer<'de>, + { + let values = Vec::::deserialize(deserializer)?; + if let Some(value) = values.iter().copied().find(|value| *value < 0) { + return Err(D::Error::custom(format!( + "expected nonnegative integers, got {value}" + ))); + } + Ok(values) + } +} + #[cfg(test)] #[path = "../../unit_tests/models/misc/knapsack.rs"] mod tests; diff --git a/src/rules/knapsack_qubo.rs b/src/rules/knapsack_qubo.rs new file mode 100644 index 000000000..1bb7a899b --- /dev/null +++ b/src/rules/knapsack_qubo.rs @@ -0,0 +1,116 @@ +//! Reduction from Knapsack to QUBO. +//! +//! Converts a nonnegative 0-1 Knapsack instance into QUBO by turning the +//! capacity inequality sum(w_i * x_i) <= C into equality using binary slack +//! variables, then constructing a QUBO that combines the objective +//! -sum(v_i * x_i) with a quadratic penalty +//! P * (sum(w_i * x_i) + sum(2^j * s_j) - C)^2. +//! For nonnegative values, penalty P > sum(v_i) ensures any infeasible solution +//! costs more than any feasible one. +//! +//! Reference: Lucas, 2014, "Ising formulations of many NP problems". + +use crate::models::algebraic::QUBO; +use crate::models::misc::Knapsack; +use crate::reduction; +use crate::rules::traits::{ReduceTo, ReductionResult}; + +/// Result of reducing Knapsack to QUBO. +#[derive(Debug, Clone)] +pub struct ReductionKnapsackToQUBO { + target: QUBO, + num_items: usize, +} + +impl ReductionResult for ReductionKnapsackToQUBO { + type Source = Knapsack; + type Target = QUBO; + + fn target_problem(&self) -> &Self::Target { + &self.target + } + + fn extract_solution(&self, target_solution: &[usize]) -> Vec { + target_solution[..self.num_items].to_vec() + } +} + +#[reduction(overhead = { num_vars = "num_items + num_slack_bits" })] +impl ReduceTo> for Knapsack { + type Result = ReductionKnapsackToQUBO; + + fn reduce_to(&self) -> Self::Result { + let n = self.num_items(); + let c = self.capacity(); + let b = self.num_slack_bits(); + let total = n + b; + + // Penalty must exceed sum of all values + let sum_values: i64 = self.values().iter().sum(); + let penalty = (sum_values + 1) as f64; + + // Build QUBO matrix + // H = -sum(v_i * x_i) + P * (sum(w_i * x_i) + sum(2^j * s_j) - C)^2 + // + // Let a_k be the coefficient of variable k in the constraint: + // a_k = w_k for k < n (item variables) + // a_{n+j} = 2^j for j < B (slack variables) + // + // Expanding the penalty: + // P * (sum(a_k * z_k) - C)^2 = P * sum_i sum_j a_i * a_j * z_i * z_j + // - 2P * C * sum(a_k * z_k) + P * C^2 + // Since z_k is binary, z_k^2 = z_k, so diagonal terms become: + // Q[k][k] = P * a_k^2 - 2P * C * a_k (from penalty) + // Q[k][k] -= v_k (from objective, item vars only) + // Off-diagonal terms (i < j): + // Q[i][j] = 2P * a_i * a_j + + let mut coeffs = vec![0.0f64; total]; + for (i, coeff) in coeffs.iter_mut().enumerate().take(n) { + *coeff = self.weights()[i] as f64; + } + for j in 0..b { + coeffs[n + j] = (1u64 << j) as f64; + } + + let c_f = c as f64; + let mut matrix = vec![vec![0.0f64; total]; total]; + + // Diagonal: P * a_k^2 - 2P * C * a_k - v_k (for items) + for k in 0..total { + matrix[k][k] = penalty * coeffs[k] * coeffs[k] - 2.0 * penalty * c_f * coeffs[k]; + if k < n { + matrix[k][k] -= self.values()[k] as f64; + } + } + + // Off-diagonal (upper triangular): 2P * a_i * a_j + for i in 0..total { + for j in (i + 1)..total { + matrix[i][j] = 2.0 * penalty * coeffs[i] * coeffs[j]; + } + } + + ReductionKnapsackToQUBO { + target: QUBO::from_matrix(matrix), + num_items: n, + } + } +} + +#[cfg(feature = "example-db")] +pub(crate) fn canonical_rule_example_specs() -> Vec { + vec![crate::example_db::specs::RuleExampleSpec { + id: "knapsack_to_qubo", + build: || { + crate::example_db::specs::direct_best_example::<_, QUBO, _>( + Knapsack::new(vec![2, 3, 4, 5], vec![3, 4, 5, 7], 7), + |_, _| true, + ) + }, + }] +} + +#[cfg(test)] +#[path = "../unit_tests/rules/knapsack_qubo.rs"] +mod tests; diff --git a/src/rules/mod.rs b/src/rules/mod.rs index 7e7f9f175..b5e3c9326 100644 --- a/src/rules/mod.rs +++ b/src/rules/mod.rs @@ -11,6 +11,7 @@ pub(crate) mod coloring_qubo; pub(crate) mod factoring_circuit; mod graph; mod kcoloring_casts; +mod knapsack_qubo; mod ksatisfiability_casts; pub(crate) mod ksatisfiability_qubo; pub(crate) mod ksatisfiability_subsetsum; @@ -79,6 +80,7 @@ pub(crate) fn canonical_rule_example_specs() -> Vec(invalid).unwrap_err(); + assert!(error.to_string().contains("nonnegative")); + } +} diff --git a/src/unit_tests/rules/knapsack_qubo.rs b/src/unit_tests/rules/knapsack_qubo.rs new file mode 100644 index 000000000..a92e552a2 --- /dev/null +++ b/src/unit_tests/rules/knapsack_qubo.rs @@ -0,0 +1,88 @@ +use super::*; +use crate::solvers::BruteForce; +use crate::traits::Problem; + +#[test] +fn test_knapsack_to_qubo_closed_loop() { + let knapsack = Knapsack::new(vec![2, 3, 4, 5], vec![3, 4, 5, 7], 7); + let reduction = ReduceTo::>::reduce_to(&knapsack); + let qubo = reduction.target_problem(); + + assert_eq!(qubo.num_vars(), 7); + + let solver = BruteForce::new(); + let best_source = solver.find_all_best(&knapsack); + let best_target = solver.find_all_best(qubo); + + let extracted: std::collections::HashSet> = best_target + .iter() + .map(|t| reduction.extract_solution(t)) + .collect(); + let source_set: std::collections::HashSet> = best_source.into_iter().collect(); + + assert!(extracted.is_subset(&source_set)); + assert!(!extracted.is_empty()); +} + +#[test] +fn test_knapsack_to_qubo_single_item() { + let knapsack = Knapsack::new(vec![1], vec![1], 1); + let reduction = ReduceTo::>::reduce_to(&knapsack); + let qubo = reduction.target_problem(); + + assert_eq!(qubo.num_vars(), 2); + + let solver = BruteForce::new(); + let best_target = solver.find_all_best(qubo); + let extracted = reduction.extract_solution(&best_target[0]); + assert_eq!(extracted, vec![1]); +} + +#[test] +fn test_knapsack_to_qubo_infeasible_rejected() { + let knapsack = Knapsack::new(vec![2, 3, 4, 5], vec![3, 4, 5, 7], 7); + let reduction = ReduceTo::>::reduce_to(&knapsack); + let qubo = reduction.target_problem(); + + let solver = BruteForce::new(); + let best_target = solver.find_all_best(qubo); + + for sol in &best_target { + let source_sol = reduction.extract_solution(sol); + let eval = knapsack.evaluate(&source_sol); + assert!( + eval.is_valid(), + "Optimal QUBO solution maps to infeasible knapsack solution" + ); + } +} + +#[test] +fn test_knapsack_to_qubo_empty() { + let knapsack = Knapsack::new(vec![1, 2], vec![3, 4], 0); + let reduction = ReduceTo::>::reduce_to(&knapsack); + let qubo = reduction.target_problem(); + + assert_eq!(qubo.num_vars(), 3); + + let solver = BruteForce::new(); + let best_target = solver.find_all_best(qubo); + let extracted = reduction.extract_solution(&best_target[0]); + assert_eq!(extracted, vec![0, 0]); +} + +#[cfg(feature = "example-db")] +#[test] +fn test_knapsack_to_qubo_canonical_example_spec() { + let spec = canonical_rule_example_specs() + .into_iter() + .find(|spec| spec.id == "knapsack_to_qubo") + .expect("missing canonical Knapsack -> QUBO example spec"); + let example = (spec.build)(); + + assert_eq!(example.source.problem, "Knapsack"); + assert_eq!(example.target.problem, "QUBO"); + assert_eq!(example.source.instance["capacity"], 7); + assert_eq!(example.target.instance["num_vars"], 7); + assert!(!example.solutions.is_empty()); +}