Blossom minimum-weight perfect matching
# LP formulation of mwpm
# notation x(A) := ∑ i ∈ A, (x i)
# Let S ⊆ G.V
# δ: edges that start in S and end outside
def δ (S) : {e ∈ G.E | e.v₁ ∈ S ∧ e.v₂ ∉ S }
# let O be the set of all subsets of V
# of odd cardinality with three or more nodes.
# Potential blossoms
Primal problem
min {∑ e ∈ G.E, e.w xₑ |
∀ e ∈ G.E, xₑ ≥ 0 ∧
∀ v ∈ G.V, x(δ v) = 1 ∧
∀ S ∈ O, x(δ S) ≥ 1 }
Dual problem
max {∑ v ∈ G.V, v.y + ∑ S ∈ O, ∑ v ∈ S, v.y |
∀ e ∈ G.E, slack e ≥ 0 ∧
∑ S ∈ O, ∑ v ∈ S, v.y ≥ 0}
-
Complexity
-
“stage: a sequence of ops between two successive augmentations.”
-
There are at most n/2 stages.
-
Each augmentation reduces the number of unpaired v’s by two.
-
This can happen at most n/2 times before the set of unpaired v’s is exhausted.
-
-
Each stage can perform at most O(n) basic ops GROW, SHRINK, EXPAND.
-
Assuming (i) each dual update results in at least one basic op and
-
(ii) each GROW, SHRINK, EXPAND, and a dual update take O(m) time,
-
we arrive at a straightforward bound O(n^2*m) for the overall complexity.
-
Note that assumption (i) is valid for the first two approaches.
-
However, in the third approach values of delta T may be determined by the constraint (4e)
-
rather than (4a), (4b), (4c), (4d).
-
Thus the O(n^2 * m) complexity is not guaranteed for the variable delta approach.”
-
Unweighted perfect matching: All the weights are equal to 1.
“Weighted matching consists of Unweighted perfect matching together with small modifications suggested by theorem (M).”
“Let C be the complex polyhedron of vectors (x : ℕ → ℝ) formed by the intersection of all the half-spaces given by the following inequalities:
1- ∀ e ∈ G.E, e.matched ≥ 0
2- ∀ v ∈ G.V, (∑ e ∈ v.edges, e.matched) ≤ 1
3- For every subset S of 2r+1 nodes, ∑ x ≤ r (summed over x in R), where R is the set of variables corresponding to the edges of G with both ends in S.”
variables
(a : ℕ → ℕ → ℝ) -- adjacency matrix?
(b : ℕ → ℝ)
(c : ℕ → ℝ) -- edge weights
(x : ℕ → ℝ) -- edge.matched
[∀ i, x i ≥ 0] [∀ j, ∑ i, (a i j) * (x i) ≤ (b j)]
(y : ℕ → ℝ) -- vertex weights
[∀ j, y j ≥ 0] [∀ i, ∑ j, (a i j) * (y j) ≤ (c i)]
/-
"The duality theorem of LP:
max ∑ i, (c i) * (x i) = min ∑ i, (b i) * (y i)
when these extrema exist for vectors xi and yi."
-/
To get the linear program dual to maximizing W of (4) in the polyhedron C, intro a new variable corresponding to each inequality of (2) and (3).
Intro a variable z for each set S of 2r + 1 nodes in G.V.
(4) W = ∑ e ∈ G.E, if e.matched = 1 then e.w else 0
Let (y,z) denote the vector of all y’s and z’s.
The duality theorem says that max W = min of
(5) U = ∑ v ∈ G.V, v.y + ∑ S, r z. The vector (y,z) satisfies the nonnegativity ineq (6) and the ineq (7) obtained from the transpose of the matrix of coefficients of ineq (2) and (3).
(6) Let ∀ v ∈ G.V, v.y ≥ 0; Let ∀ S, z ≥ 0.
(7) ∀ e ∈ G.E, e.v₁.y + e.v₂.y + (∑ S, if e.v₁ ∈ S ∧ e.v₂ ∈ S then S.z else 0) ≥ e.w.
(8) unpaired v → v.y = 0
(9) e.matched = 1 → e.v₁.y + e.v₂.y + ∑ z = e.w
(10) for each z > 0, the set S contains both endpoints of r edges in M.”
def satisfies_vertex_constraints (G : Graph) : Prop :=
(v_weight_nonneg : ∀ v ∈ G.V, v.y ≥ 0)
(e_weight_leq_sum_ends : ∀ e ∈ G.E, e.v₁.y + e.v₂.y ≥ e.w)
(exists_blossom_if_not_done : ) -- For i=0,..,n-1, there is a length 2ai + 1 blossom Bi in Gi with ai matched edges.
(e_weight_eq_sum_ends_in_blossom : ∀ e ∈ B.E, e.v₁.y + e.v₂.y = e.w)
(smallest_v_in_blossom_of_unpaired : ∀ v ∈ B.V, unpaired v → (v.y is smallest among B.V))
(next_G_wghts_same_of_not_blossom : ) -- (g) Weights in G(i+1) are the same as corresponding weights in Gi. Except at u(i+1) and u(i+1).E.
(blossom_weight_leq_smallest_v : ) -- Let q ∈ B.V be the vertex with the minimum y. Then B.y ≤ q.y
(next_G_blossom_weights : ) -- ∀ (ei+1) ∈ B(i+1).E, meeting v of Bi. Then (ei+1).w = (ei).w - (vi).y + q.y.
(e_weight_eq_sum_ends_of_matched : ) -- ∀ e ∈ Gₙ.E, e.matched = 1 → e.v₁.y + e.v₂.y = e.w
def vertex_weight_zero_of_unpaired (G : Graph) : Prop :=
∀ v ∈ G.V, unpaired v → v.y = 0
-- some vertices are not connected by an edge at all?
theorem M: Matching M is perfect ↔ ∃ (Gs : ℕ → G),
satisfies_vertex_constraints (Gs 0)
∧ each graph is obtained from the previous graph following the algorithm
∧ vertex_weight_zero_of_unpaired (Gs n) :=
begin
sorry,
end
-- Conditions (a) through (j) are preserved by the steps of the algorithm."
dual updates update the weight of the nodes. And increase the value of the objective function whose max we seek.
-
the new node weights must satisfy the constraints of the dual program
-
slack(e) ≥ δT
-
slack(e) ≥ δT + δT’
-
B.y ≥ δT
-
slack(e) ≥ 2 δT
slack(e) = e.weight - e.v1.y - e.v2.y - ∑ over odd cycles that have e as a boundary edge, S.y
Primal updates
Look at the neighbors: is the edge to the neighbor tight? Then continue as in Unweighted blossom.
Initialize blossoms with weight 0. A blossom should have weight 0 before it opens
Dual updates: update the node weights
Blossom forming edges are tight because otherwise we are not allowed to walk them to make a cycle that makes the blossom
Application: Surface code
. ⇅ . ⇅ . ⇅ . ⇅ . ⇅ .
⇅ ⇅ ⇅ ⇅ ⇅ ⇅
. ⇅ . ⇅ . ⇅ . ⇅ . ⇅ .
⇅ ⇅ ⇅ ⇅ ⇅ ⇅
. ⇅ . ⇅ . ⇅ . ⇅ . ⇅ .
⇅ ⇅ ⇅ ⇅ ⇅ ⇅
. ⇅ . ⇅ . ⇅ . ⇅ . ⇅ .
⇅ ⇅ ⇅ ⇅ ⇅ ⇅
. ⇅ . ⇅ . ⇅ . ⇅ . ⇅ .
. ⇅ . ⇅ . ⇅ . ⇅ . ⇅ .
⇅ ⇅ ⇅ ⇅ ⇅ ⇅
. ⇅ . ⇅ . ⇅ Z ⇅ . ⇅ .
⇅ ⇅ ⇅ ⇅ ⇅ ⇅
. ⇅ . ⇅ . ⇅ . ⇅ . ⇅ .
⇅ ⇅ ⇅ ⇅ ⇅ ⇅
. ⇅ . ⇅ . ⇅ . ⇅ . ⇅ .
⇅ ⇅ ⇅ ⇅ ⇅ ⇅
. ⇅ . ⇅ . ⇅ . ⇅ . ⇅ .
Z ⇅ Z ⇅ Z ⇅ Z ⇅ Z ⇅ Z
⇅ ⇅ ⇅ ⇅ ⇅ ⇅
Z ⇅ Z ⇅ Z ⇅ Z ⇅ Z ⇅ Z
⇅ ⇅ ⇅ ⇅ ⇅ ⇅
Z ⇅ Z ⇅ Z ⇅ Z ⇅ Z ⇅ Z
⇅ ⇅ ⇅ ⇅ ⇅ ⇅
Z ⇅ Z ⇅ Z ⇅ Z ⇅ Z ⇅ Z
⇅ ⇅ ⇅ ⇅ ⇅ ⇅
Z ⇅ Z ⇅ Z ⇅ Z ⇅ Z ⇅ Z
. ⇅ ◯ ⇅ . ⇅ . ⇅ . ⇅ .
⇅ ⇅ ⇅ ⇅ ⇅ ⇅
. ⇅ . ⇅ . ⇅ . ⇅ ◯ ⇅ .
⇅ ⇅ ⇅ ⇅ ⇅ ⇅
. ⇅ . ⇅ ◯ ⇅ . ⇅ . ⇅ .
⇅ ⇅ ⇅ ⇅ ⇅ ⇅
. ⇅ ◯ ⇅ . ⇅ . ⇅ ◯ ⇅ .
⇅ ⇅ ⇅ ⇅ ⇅ ⇅
. ⇅ . ⇅ ◯ ⇅ . ⇅ . ⇅ .
. ◯ . . . . . . . . ◯ . . . ◯ . . . . ◯ . . ◯ . . . ◯ . . .
We want a path to an unpaired node
s
●---.---.---o
| |
o---.---.---o
s
●xxx.xxx.xxxo
| |
o---.---.---o
s s
●xxx.xxx.xxxo ●---.---.---o
| | x |
o---.---.---o o---.---.---o
A new rule disallows walking to the faraway node:
s s
●xxx.xxx.xxxo ●---.---.---o
| | x |
o---.---.---o o---.---.---o
A new rule disallows walking to the faraway node:
s s
●xxx.xxx.xxxo ●---.---.---o
| | x |
o---.---.---o o---.---.---o
Walk only over "tight" edges
Initialize each node with weight w/2
o---.---.---o
| |
o---.---.---o
w is the weight of the shortest incident edge
Initialize each node with weight w/2
1/2 1/2
o---.---.---o
| |
o---.---.---o
1/2 1/2
w is the weight of the shortest incident edge
the short edges are tight
1/2 1/2
o---.---.---o
| |
o---.---.---o
1/2 1/2
Node unpaired? Increase weight
1 2
..---●xxx.xxx●---.---o ..---●xxx.xxx●---.---o
.............. ..............
until an incident edge becomes tight
Node is a leaf? Increase its weight
2 1 1 2
..---●xxx.xxx.xxx●---.---o ..---●xxx.xxx.xxx●---.---o
................. ..................
Decrease mate's weight by the same amount
New blossom? Assign weight 0
: :
/ /
. s ...... / /
: ●-------● : /
: | x : ❀ 0
: ---●xxxx : |
..... | ..... |
| |
: :
from Unweighted_perfect_matching import Node, Edge, Blossom, Tree
class WeightedEdge(Edge):
def __init__(self, weight):
super(self, Edge).__init__()
self.weight = weight
def tight(self):
return self.weight - self.node1.weight - self.node2.weight == 0
class WeightedNode(Node):
def __init__(self, weight):
super(self, Node).__init__()
self.weight = weight
class WeightedBlossom(Blossom):
def __init__(self, weight):
super(self, Blossom).__init__()
self.weight = 0
class Tree:
def __init__(self, seed):
self.seed = seed
self.seed.tree = self
self.leaves = []
self.active_node = self.seed
def grow(self):
found_unpaired_node = False
while not found_unpaired_node:
for neighbor in self.active_node.neighbors():
if neighbor.unpaired():
found_unpaired_neighbor = True
self.handle_unpaired(neighbor)
return neighbor
else:
if neighbor.tree == self:
cycle = self.handle_cycle(neighbor)
if len(cycle) % 2 == 1:
self.active_node = Blossom(cycle)
else:
self.handle_paired(neighbor)
self.active_node = self.leaves.pop(0)
def handle_unpaired(self, node):
self.active_node.add_child(node)
self.walk_back(node)
def handle_paired(self, node):
self.active_node.add_child(node)
mate = node.mate()
node.add_child(mate)
self.leaves.append(mate)
def handle_cycle(self, neighbor):
# branch from common node to neighbor
path1 = [neighbor]
# branch from common node to active_node
path2 = [active_node]
while path1[-1] != path2[-1]:
while len(path1[-1].children) == 1:
path1.append(path1[-1].parent)
while len(path2[-1].children) == 1:
path2.append(path2[-1].parent)
# common starting point found
cycle = path1 + path2[:-2]
return cycle
# walk back to the seed and invert edges
def walk_back(self, node):
pos = node
edge_state = False
# open blossoms in the neighborhood
while pos.parent.is_blossom:
pos.parent.open()
while pos != self.seed:
pos.get_edge_to(pos.parent).invert()
pos = pos.parent
while len(nodes) > 0:
t = Tree(seed=nodes.pop())
nodes.remove(t.grow())
for e in edges:
if e.matched:
print("paired nodes", e.node1, e.node2)
Modify Tree.grow() method in Unweighted perfect matching so only tight edges are used for growing.