AC/DC search (Python)

A Python translation of the MATLAB reference implementation of AC/DC search — the Alternating Continuous and Discrete Combinatorial optimization behind the winning solution to the FlyWire Ventral Nerve Cord Matching Challenge.

D. D. Lee, A. Matsliah & L. K. Saul, "AC/DC search: behind the winning solution to the FlyWire graph-matching challenge", Transactions on Machine Learning Research (01/2026). OpenReview

What it does

Given two weighted directed graphs A and B (same number of nodes), AC/DC searches for the permutation P (a one-to-one node correspondence) that maximizes the overlapping edge weight

J(P) = sum_ij  min( A_ij , B_{p(i), p(j)} )

where p(i) is the node of B matched to node i of A. It alternates:

• a continuous phase — Frank–Wolfe optimization over the doubly stochastic relaxation (Birkhoff polytope), with an exact closed-form line search and a linear-assignment subproblem at each step; and
• a discrete phase — greedy pairwise-swap search over permutations, applying the exact-gain swaps that improve the score.

Each phase warm-starts the other.

Install

pip install acdc-search

The distribution is named acdc-search (the names acdc and acdc-py are taken on PyPI by unrelated projects); the import name is acdc. It is pure Python and needs only NumPy and SciPy.

For a development install from a checkout:

pip install -e ".[test]"        # numpy, scipy + pytest

Usage

import numpy as np
import scipy.sparse as sp
from acdc import acdc_match, score, matrix_to_perm

# --- build two graphs A, B as scipy.sparse n x n matrices ---
rng = np.random.default_rng(0)
A = sp.random(200, 200, density=0.05, random_state=0).tocsr()

# (toy example) B is A relabeled by a hidden permutation we will recover
perm = rng.permutation(200)
P = sp.csc_matrix((np.ones(200), (np.arange(200), perm)), shape=(200, 200))
B = (P.T @ A @ P).tocsr()

# --- run the full AC/DC algorithm ---
M = acdc_match(A, B, max_iter=5, num_frank_wolfe=10)

print("score:", score(M, A, B), " (optimum:", A.sum(), ")")
matching = matrix_to_perm(M)          # 0-based: node i of A -> matching[i] of B

High-level functions

All three return a sparse permutation matrix.

Progress output. The phases emit their MATLAB-style tables through the acdc logger at INFO. verbose=True (the default) attaches a stdout handler for the duration of the call, so they just appear; verbose=False leaves your logging configuration alone, which means an application that has configured logging at INFO still receives them and one that has not stays silent:

import logging
logging.basicConfig(level=logging.INFO)         # or logging.getLogger("acdc")
M = acdc_match(A, B, verbose=False)             # tables go to your handlers

Inputs. A and B may be a SciPy sparse matrix, a dense numpy.ndarray, or an edge-list tuple (rows, cols, weights) / (rows, cols, weights, n) with 0-based node indices. P0 may be a permutation vector, a sparse permutation matrix, or None (identity start).

Public API

acdc exports exactly six names: the three entry points above, plus score to evaluate a matching and perm_to_matrix / matrix_to_perm to convert between a permutation matrix and a 0-based permutation vector.

The individual algorithm phases live in the submodules that mirror the MATLAB files — acdc.objective (compute_gradient), acdc.frank_wolfe (do_frank_wolfe), acdc.swaps (evaluate_swaps, make_swaps, greedy_search), acdc.matching (permutation_match, as_graph, …) — and can be imported from there. They are internal: their signatures follow the call graph rather than any user-facing contract and may change without a major version bump.

File map (MATLAB → Python)

Notes on fidelity & performance

• solver='dense' (default) uses SciPy's exact Hungarian/LAPJV solver (linear_sum_assignment) — the faithful analogue of the MATLAB perfectMatching. solver='sparse' uses min_weight_full_bipartite_matching, which is much faster and lighter at scale but only considers stored edges (a perfect matching must exist on the sparsity pattern).
• The MATLAB warm-start preconditioner for the assignment solver is omitted on purpose: its transform only adds per-row/column constants and does not change the optimal matching — it existed solely to speed MATLAB's sparse solver.
• Like the MATLAB original, the gradient and swap-gain matrices are dense n × n, so memory is O(n²). At challenge scale (n = 18524, see below) that needs ~16 GB+ of RAM and is best run with solver='dense' on a machine with enough memory; small and medium graphs run comfortably anywhere.
• The quadratic part of the swap-gain matrix has two implementations. Six of its eight terms are min(x, y) with one argument taken from A or P B P', so for nonnegative edge weights they vanish off those supports and can be scattered onto the union of the supports instead of swept densely. On sparse graphs that holds 2 dense n × n arrays instead of 5 (≈5.5 GB rather than ≈13.7 GB at challenge scale). Scattered writes are slower per element, so dense or negatively-weighted graphs use the literal dense expression; the two paths agree bit for bit and the choice is automatic.
• The greedy phase recomputes the full score after each trial swap (faithful to the MATLAB original); incremental scoring is a possible future optimization.
• Both loops stop early at a fixed point: a Frank–Wolfe update that leaves the iterate unchanged, or an AC/DC alternation that leaves the matching unchanged, makes every remaining iteration a deterministic replay. max_iter and num_frank_wolfe are therefore upper bounds, not exact counts. The result is unchanged — only the wasted work is skipped.

TODO

Pinned (fixed) matches

Neither the MATLAB original nor this translation supports pinning a subset of matches so they are never swapped out. Every pair is free to move at every stage: the discrete phase picks its swap from a max over the whole gain matrix (make_swaps.m → swaps.py), and each Frank–Wolfe step solves an unconstrained linear assignment problem over the full gradient (permutation_match.m → matching.py) — the P0 argument there is only a warm-start accelerator and does not restrict the solution. This would be useful whenever part of the correspondence is known a priori, e.g. neurons matched with high confidence by cell type.

Adding it means two hooks, and would surface as a pinned= argument on the three entry points:

• Discrete phase — mask the swap-gain matrix D: zero the row and column of each pinned index where D is assembled in evaluate_swaps, so that dMax also respects the pins and the while dMax > 0 loop still terminates.
• Continuous phase — constrain the assignment problem. Prefer reducing it (drop the pinned rows/columns, solve on the free submatrix, reinsert the pinned pairs) over penalising forbidden entries with a large negative cost: the reduced LAP is strictly smaller and avoids -inf handling in linear_sum_assignment entirely. Only the two permutation_match calls need constraining — the interpolation P + step*(Q - P) with step in (0, 1] preserves any entry that is 1 in both P and Q, so pinned entries carry through the continuous iterate on their own.

Benchmark

benchmark.py runs AC/DC on synthetic planted instances of increasing size to give a feel for scaling. Each instance relabels a random graph by a hidden permutation, corrupts that alignment to make a warm start (40 % of nodes scrambled by default), then measures runtime and how much of the planted optimum is recovered.

python benchmark.py                                  # default sweep, ~1 min
python benchmark.py --sizes 100 200 400 800          # quick
python benchmark.py --sizes 1000 2000 --density 0.02 --solver sparse

Example run (laptop, solver='dense', max_iter=3, num_frank_wolfe=10):

      n     edges       mem     grad     eval      acdc     cont     disc  start%  final%  nodes%
-------------------------------------------------------------------------------------------------
    100       481    78.1KB    0.001    0.002     0.014    0.010    0.004    38.9   100.0   100.0
    200      1935   312.5KB    0.003    0.003     0.021    0.013    0.008    36.0   100.0   100.0
    400      7771     1.2MB    0.007    0.009     0.060    0.038    0.022    39.3   100.0   100.0
    800     31147     4.9MB    0.028    0.035     0.219    0.139    0.080    38.1   100.0   100.0
   2000    194913    30.5MB    0.522    0.572     3.483    2.255    1.227    38.2   100.0   100.0
   5000   1219064   190.7MB    8.149    8.837    50.284   32.975   17.308    38.2   100.0   100.0
- - - - - - - - - - - - - - - - - - - - - - - -
  extrapolated from the n=800..5000 rows (acdc ~ n^2.97, R2=1.000 in log-log) -- NOT measured
 ~10000  4.88e+06   762.9MB     71.6     72.2       399      266      134       -       -       -
 ~50000  1.22e+08    18.6GB 1.04e+04 9.29e+03  4.73e+04 3.25e+04 1.52e+04       -       -       -
~100000  4.88e+08    74.5GB 8.89e+04 7.52e+04   3.7e+05 2.57e+05 1.16e+05       -       -       -

• grad/eval are per-call times (in seconds) for compute_gradient / evaluate_swaps (best of as many repeats as fit in a short budget, so small sizes are not dominated by noise). acdc is the full run, split into the continuous (cont) and discrete (disc) phases — roughly a 2:1 split throughout. mem is the size of one dense n × n float64 matrix; the n = 5000 row peaks at about 1.1 GB resident.
• Rows prefixed ~ are extrapolations, not measurements. Each timing column is fitted as a power law over the asymptotic tail of the measured rows (n ≥ 800) and evaluated at the larger size. Fitting the small rows instead gives an exponent near 1.3 and under-predicts n = 5000 by ~25×, because those sizes are dominated by fixed overhead — the benchmark warns when the fitted exponent lands below 2.5. Held out against a real measurement, a fit over n = 400..2000 predicts 33 s at n = 5000 against 42 s measured.
• Scaling is O(n³) at fixed density, not O(n²): the gradient loops over the n matched nodes and does work quadratic in the degree, which at fixed density is itself proportional to n. The measured exponent is 2.97 with R² = 1.000. At fixed average degree — the realistic regime for a connectome — the same argument gives O(n²), dominated by the dense matrices. Which of the two applies to your graph is the single biggest factor in what it will cost.
• The extrapolated rows hold density fixed at 0.05, so they are an upper bound for a graph of that size rather than a prediction for a real one. The challenge connectomes are n = 18524 at density 0.012 (male) and 0.006 (female) — 4–8× sparser than this sweep. At these sizes memory binds well before time does: 74.5 GB for a single array at n = 100 000, and the greedy phase holds two to three of them.
• From a corrupted-but-informative warm start the planted alignment is recovered essentially perfectly (final% ≈ 100). Quality from a cold (identity) start is much lower — AC/DC is a local search and, like the paper, relies on a reasonable warm start.

The challenge data

The challenge was to align the connectomes of the ventral nerve cords (VNCs) of a male and a female fruit fly — not the FlyWire brain connectome, despite the challenge carrying the FlyWire name. Per the paper, each VNC connectome is a directed weighted graph with n = 18524 nodes (neurons) and millions of edges (synapse counts); the organizers also supplied a cell-type-derived baseline match, scoring 5 154 247, which teams could use as a warm start.

The two graphs are that size in the data shipped here as well: the male edge list spans node ids m1..m18524 and the female f1..f18524 (the female list has 18 523 nodes with at least one edge — f9574 is isolated, which is why _coerce_pair pads to a common n). Their densities are 0.012 and 0.006 respectively, so both are far sparser than the benchmark's default 0.05.

CSV / challenge I/O

This package translates the core algorithm only. The challenge-specific CSV readers/writers (read_connectome.m, read_solution.m, save_solution.m) and the figure scripts are intentionally out of scope; build A, B as sparse matrices (or pass edge lists) as shown above.

Tests

pytest -q

Covers: a finite-difference gradient check, the score identity, self-match optimality, recovery of a known permutation, exactness of the swap-gain matrix, dense/sparse backend agreement, and end-to-end smoke tests.

Attribution & licence

This package is a translation of the MATLAB reference implementation that Lee, Matsliah & Saul published as supplementary material to their TMLR paper:

D. D. Lee, A. Matsliah & L. K. Saul, "AC/DC search: behind the winning solution to the FlyWire graph-matching challenge", Transactions on Machine Learning Research (01/2026). https://openreview.net/forum?id=8MjCOMyaDf

That original MATLAB code is MIT-licensed, Copyright (c) 2025 Daniel Lee and Lawrence Saul. This Python translation is therefore also offered under the MIT licence; LICENSE carries both notices, with the upstream licence reproduced in full as it requires.

NOTICE summarises how the translation departs from the original — 0-based indexing, SciPy in place of MATLAB's internal perfectMatching, the omitted assignment preconditioner, the untranslated challenge I/O, and the added tests and benchmark. The "Notes on fidelity & performance" section above covers the same ground in more detail.

If you use this package, please cite the paper above — the algorithm is the authors' work; this is only a translation of it.