Skip to content

Bound the MILP FiveTuple rescaling by the encoded depth (#65) #106

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
jkhav wants to merge 4 commits into from
Closed

Bound the MILP FiveTuple rescaling by the encoded depth (#65) #106

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
9eec3f5
fix: bound the FiveTuple rescaling by the depth, not by the input's m…
ondrik Oct 13, 2026
18e94ba
fix: reject a non-integral FiveTuple k in vector pairs
ondrik Oct 13, 2026
f269ec7
fix: integral_k check before synthesis even begins -- crashed pysmt b…
jkhav Oct 17, 2026
0133d27
fix: bound the MILP FiveTuple rescaling by the encoded depth (#65)
jkhav Oct 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions src/gates.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,12 @@

self_adjoints = ['h', 'x', 'y', 'z', 'cx', 'cz', 'cy', 'ch', 'swap', 'ccx', 'ccz']

# FiveTuple: how much a gate raises k (the sqrt(2)^k denominator), the other gates keep it
k_increments = {'h': 1, 'ch': 1, 'sx': 2, 'sxdg': 2, 'csx': 2, 'sqrtswap': 2}

def max_k_increment(gate_set) -> int:
return max((k_increments.get(gate, 0) for gate in gate_set), default=0)

# OpenQASM 3 definitions of supported gates that are not part of stdgates.inc
nonstandard_gate_definitions = {
'sxdg': "gate sxdg a { inv @ sx a; }",
Expand Down
18 changes: 10 additions & 8 deletions src/generator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -752,8 +752,9 @@ def multiply_by_identity_unscaled(result_vec, source_vec):
return self.And(*eqs)

if self.logic == "QF_LIA" or self.logic == "QF_NIA":
if (d // 2) == 0:
d = 2
# the case split below covers k in [0, d); a larger k would leave r1 and r2 unconstrained, and any
# circuit would pass the equivalence check (#65)
self.add_assertion(self.LT(k, self.Int(d)))
for i in range(d):
pow2 = self.Int(2 ** i)
r1_multiply_by_m = self.And(multiply_by_m_scaled(r1, v1, pow2), multiply_by_identity_unscaled(r2, v2))
Expand All @@ -768,14 +769,15 @@ def multiply_by_identity_unscaled(result_vec, source_vec):
else:
raise ValueError("rescaling with wrong logic")

def add_milp_rescaling(self, r1, r2, v1, v2, pair_idx, d):
# first encode k as the result of the operation floor(abs(k1 - k2)/2)
def add_milp_rescaling(self, r1, r2, v1, v2, pair_idx, max_gap, value_bigM):
# first encode k as the result of the operation floor(abs(k1 - k2)/2), where abs(k1 - k2) <= max_gap
# value_bigM bounds |v - r| for every entry of v1, v2 and r1, r2
complex_representation = v1.element_representation
if self.mode == "milp":
k = self.declare_integer(f"k{pair_idx}", lb=0, ub=d)
d = max_gap // 2
k = self.declare_integer(f"k{pair_idx}", lb=0)
q = self.declare_bool(f"q{pair_idx}")

bigM = 2*d + 1
bigM = 2*max_gap + 1
sleq = self.declare_bool(f"sleq{pair_idx}")
self.add_assertion((v1.k - v2.k) - (2*k + q) <= bigM * sleq)
self.add_assertion((2*k + q) - (v1.k - v2.k) <= bigM * sleq)
Expand All @@ -789,7 +791,7 @@ def add_milp_rescaling(self, r1, r2, v1, v2, pair_idx, d):
self.add_assertion(k == self.Sum([constants[i] * i for i in range(d+1)]))
# now sleq, si, even encode all case splits needed for the rescaling
# encode all combinations to assign to r1, r2
bigM = (2**(d+1)) + 1
bigM = value_bigM
# for every si, I,M is multiplied by 2^i
for i in range(len(v1)):
for j, s in enumerate(constants):
Expand Down
50 changes: 36 additions & 14 deletions src/synth.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,14 +18,20 @@
from pysmt.shortcuts import Portfolio, Symbol, Real, And, Equals, Plus, GT, LT, get_env, Int, Or, Not, Implies, GE, LE, Ite, Minus, Div, Times, write_smtlib, Solver
from pysmt.typing import REAL, INT
from pareto import Pareto
from gates import GateSet, supported_gates, check_supported, Circuit
from gates import GateSet, check_supported, Circuit, k_increments, max_k_increment
from pulp import *
from gurobipy import GRB
from solvers import PortfolioSMTSolver, SMTSolver
import time
from logger import Logger
import numpy as np

def integral_k(k):
# a FiveTuple vector is divided by sqrt(2)^k, k has to be an integer (4.0 is accepted)
if k is None or int(k) != k:
raise ValueError(f"k of a FiveTuple vector must be an integer, got {k!r}")
return int(k)

def rus_recovery_pairs(output_vectors, vector_pairs, measured_qubits):
# (input, output) pairs for the RUS recovery operation: it maps the state after a failed
# attempt (measured qubits projected to 1) back to the original input state
Expand All @@ -52,6 +58,8 @@ def __init__(self) -> None:
self.q = None
self.d = None
self.max_k = 0
self.pair_k = []
self.k_increment = 0
self.encoding_method = None
self.curr_depth = 0
self.v = 0
Expand Down Expand Up @@ -175,7 +183,7 @@ def add_k_incr(sel, out_k, inp_k, increment):
propagate_identities[pos].append(bool_var)
propagate_identities[other].append(bool_var)
if self.complex_representation == FiveTuple:
add_k_incr(bool_var, out.k, inp.k, 1)
add_k_incr(bool_var, out.k, inp.k, k_increments[gate])
elif self.basis == "pauli":
raise NotImplementedError("Pauli basis not supported for H gate")
elif gate == 's':
Expand Down Expand Up @@ -259,7 +267,7 @@ def add_k_incr(sel, out_k, inp_k, increment):
propagate_identities[pos].append(bool_var)
propagate_identities[other].append(bool_var)
if self.complex_representation == FiveTuple:
add_k_incr(bool_var, out.k, inp.k, 2)
add_k_incr(bool_var, out.k, inp.k, k_increments[gate])
elif self.basis == "pauli":
raise NotImplementedError("Pauli basis not supported for SX gate")
elif gate == 'sxdg':
Expand All @@ -275,7 +283,7 @@ def add_k_incr(sel, out_k, inp_k, increment):
propagate_identities[pos].append(bool_var)
propagate_identities[other].append(bool_var)
if self.complex_representation == FiveTuple:
add_k_incr(bool_var, out.k, inp.k, 2)
add_k_incr(bool_var, out.k, inp.k, k_increments[gate])
elif self.basis == "pauli":
raise NotImplementedError("Pauli basis not supported for SXDG gate")
elif gate == 'y':
Expand Down Expand Up @@ -392,7 +400,7 @@ def add_k_incr(sel, out_k, inp_k, increment):
elif cannot_propagate:
saved_operations[pos].append((bool_var, inp[pos]))
if self.complex_representation == FiveTuple:
add_k_incr(bool_var, out.k, inp.k, 1)
add_k_incr(bool_var, out.k, inp.k, k_increments[gate])
elif self.basis == "pauli":
raise NotImplementedError("Pauli basis not supported for CH gate")
elif gate == 'csx':
Expand All @@ -416,7 +424,7 @@ def add_k_incr(sel, out_k, inp_k, increment):
elif cannot_propagate:
saved_operations[pos].append((bool_var, inp[pos]))
if self.complex_representation == FiveTuple:
add_k_incr(bool_var, out.k, inp.k, 2)
add_k_incr(bool_var, out.k, inp.k, k_increments[gate])
elif self.basis == "pauli":
raise NotImplementedError("Pauli basis not supported for CSX gate")
elif gate == 'cy':
Expand Down Expand Up @@ -553,7 +561,7 @@ def add_k_incr(sel, out_k, inp_k, increment):
elif cannot_propagate:
saved_operations[pos].append((bool_var, inp[pos]))
if self.complex_representation == FiveTuple:
add_k_incr(bool_var, out.k, inp.k, 2)
add_k_incr(bool_var, out.k, inp.k, k_increments[gate])
elif self.basis == "pauli":
raise NotImplementedError("Pauli basis not supported for SQRTSWAP gate")
elif len(rest) == 4: # three qubit
Expand Down Expand Up @@ -765,7 +773,7 @@ def encode_equivalence(self, vectors1, vectors2, depth, already_measured : bool
for i in range(8):
selectors.append(self.gen.declare_bool(f"global_phase_sel_{i}"))
self.gen.ExactlyOne(*selectors)
vec = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Multiplied_by_global_phase_{vec2.name}", bound=2**self.max_k)
vec = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Multiplied_by_global_phase_{vec2.name}", bound=vec2.bound)
vec_copy = vec2.copy()
for i in range(8):
# add selection shift
Expand All @@ -783,12 +791,19 @@ def encode_equivalence(self, vectors1, vectors2, depth, already_measured : bool
vec_target = vec

# fivetuple rescaling
rescaled1 = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Rescaled1_{vec1.name}", bound=2**self.max_k)
rescaled2 = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Rescaled2_{vec2.name}", bound=2**self.max_k)
# after `depth` layers, k lies in [k_in, k_in + self.k_increment*depth], which bounds |k - k_target|
# (self.max_k is fixed from the input and too small past its depth)
k_in, k_target = self.pair_k[pair_idx]
max_gap = max(k_in + self.k_increment * depth - k_target, k_target - k_in)
# the rescaling multiplies by at most 2^(max_gap // 2), and an odd gap adds two components
value_bound = max(vec1.bound, vec_target.bound)
rescaled_bound = 2 * value_bound * 2**(max_gap // 2)
rescaled1 = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Rescaled1_{vec1.name}", bound=rescaled_bound)
rescaled2 = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Rescaled2_{vec2.name}", bound=rescaled_bound)
if self.gen.mode == "milp" or self.gen.mode == "gurobi":
self.gen.add_milp_rescaling(rescaled1, rescaled2, vec1, vec_target, pair_idx, self.max_k)
self.gen.add_milp_rescaling(rescaled1, rescaled2, vec1, vec_target, pair_idx, max_gap, value_bound + rescaled_bound)
else:
self.gen.add_rescaling(rescaled1, rescaled2, vec1, vec_target, pair_idx, self.max_k)
self.gen.add_rescaling(rescaled1, rescaled2, vec1, vec_target, pair_idx, max_gap // 2 + 1)
for i in range(2**self.q):
self.gen.add_assertion(self.gen.Equals(rescaled1[i], rescaled2[i]))

Expand Down Expand Up @@ -943,6 +958,9 @@ def synthesis(self, qasm_file=None, matrix=None, vectors=None, vector_pairs=None
self.gate_set = gate_set
self.qubits_to_measure = []
self.v = len(vector_pairs)
if self.complex_representation != Complex:
for state in (state for pair in vector_pairs for state in pair):
state.k = integral_k(state.k)
self.curr_depth = self.d
max = 0
for pair_idx, (input_vector, output_vector) in enumerate(vector_pairs):
Expand Down Expand Up @@ -1043,6 +1061,10 @@ def _synth(self, vector_pairs, output_qasm="circuit.qasm") -> bool:
# encode vector bounds
inter_states = []
target_states = []
# k of every (input, target) pair and its largest increase per layer, used to bound the FiveTuple k
if self.complex_representation != Complex:
self.pair_k = [(input_state.k, output_state.k) for input_state, output_state in vector_pairs]
self.k_increment = max_k_increment(self.gate_set)
bounds = [1] * (self.d + 1)
if self.complex_representation != Complex:
for d in range(self.d + 1):
Expand All @@ -1068,7 +1090,7 @@ def _synth(self, vector_pairs, output_qasm="circuit.qasm") -> bool:
for d in range(self.d + 1):
Inter = None
if self.complex_representation == FiveTuple:
Inter = class_to_use(q=2**self.q, name=f"I_{pair_idx}_{d}", generator=self.gen, element_representation=self.complex_representation, k=input_state.k if d == 0 else 0, n = input_state.n, bound = bounds[d], k_bound = 2 * (d+1))
Inter = class_to_use(q=2**self.q, name=f"I_{pair_idx}_{d}", generator=self.gen, element_representation=self.complex_representation, k=input_state.k if d == 0 else 0, n = input_state.n, bound = bounds[d], k_bound = input_state.k + self.k_increment * d)
else:
Inter = class_to_use(q=2**self.q, name=f"I_{pair_idx}_{d}", generator=self.gen, element_representation=self.complex_representation, bound = bounds[d], k_bound = 2 * (d+1))

Expand Down Expand Up @@ -1169,7 +1191,7 @@ def _synth(self, vector_pairs, output_qasm="circuit.qasm") -> bool:
for pair_idx in range(len(vector_pairs)):
if self.curr_depth > self.d:
if self.complex_representation == FiveTuple:
state = class_to_use(q=2**self.q, name=f"I_{pair_idx}_{self.curr_depth}", generator=self.gen, element_representation=self.complex_representation, k=0, n = vector_pairs[pair_idx][0].n, bound = self.bound, k_bound = 2 * (self.curr_depth+1))
state = class_to_use(q=2**self.q, name=f"I_{pair_idx}_{self.curr_depth}", generator=self.gen, element_representation=self.complex_representation, k=0, n = vector_pairs[pair_idx][0].n, bound = self.bound, k_bound = vector_pairs[pair_idx][0].k + self.k_increment * self.curr_depth)
else:
state = class_to_use(q=2**self.q, name=f"I_{pair_idx}_{self.curr_depth}", generator=self.gen, element_representation=self.complex_representation, bound = self.bound, k_bound = 2 * (self.curr_depth+1))
inter_states[pair_idx].append(state)
Expand Down
Loading
Loading