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 needs NumPy, SciPy and numba.

For a development install from a checkout:

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

numba. The three inner loops that dominate the runtime are compiled (acdc/_kernels.py); see Notes on fidelity & performance. They are cached to disk, so the first call in a fresh install pays about a second of compilation and later runs start in ~0.1 s. If the package is installed somewhere unwritable, numba silently recompiles once per process instead — set NUMBA_CACHE_DIR to a writable path to avoid that.

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

Function	MATLAB equivalent	Description
acdc_match(A, B, P0=None, max_iter=5, num_frank_wolfe=10, max_swap=inf, solver='dense')	main_acdc.m	Full alternating algorithm.
frank_wolfe_search(A, B, P0=None, num_updates=40, solver='dense')	main_continuous.m	Continuous (Frank–Wolfe) phase only.
greedy_match(A, B, P0=None, max_swap=inf)	main_discrete.m	Discrete (greedy-swap) phase only.

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, gradient_entries), 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. acdc._kernels holds the compiled inner loops and is private.

File map (MATLAB → Python)

MATLAB (../src/)	Python (acdc/)
compute_gradient.m + score expression	objective.py
permutation_match.m	matching.py
do_frank_wolfe.m	frank_wolfe.py
evaluate_swaps.m, make_swaps.m, greedy_search.m	swaps.py
main_acdc.m, main_continuous.m, main_discrete.m	core.py
— (compiled inner loops, no MATLAB counterpart)	_kernels.py

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.
• Three inner loops are compiled with numba (acdc/_kernels.py): the gradient accumulation, the greedy pairwise-swap loop, and the support-restricted part of the swap-gain quadratic term. Each walks sparse adjacency structure one entry at a time over only a few dozen elements, so in NumPy they were dominated by per-call dispatch overhead rather than arithmetic. The kernels visit elements in the same order as the reference implementation, so results are bit-identical to the pure-NumPy version they replace, not merely close.
• 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 are applied cell by cell over 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). 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 scores each trial swap incrementally. Exchanging the images of i and j can only change cells in rows and columns i and j, so the exact score change follows from O(deg) work where make_swaps.m recomputes the whole O(nnz) score after every trial. The accept/reject decision is identical — only the cost of reaching it differs — and the score is re-derived in full once per greedy iteration, so nothing accumulates. This was the single largest saving: on one greedy pass at n = 1200 it took 148.6 s down to 0.65 s.
• The Frank–Wolfe line search never materializes the vertex gradient. It reads ∇J(Q) only on the support of Q − P, which has O(n) nonzeros, so those entries are evaluated pointwise instead of building a dense n × n array and discarding all but O(n) of it.
• A and B are needed in both CSR and CSC form by every gradient and every swap evaluation, but are fixed for a whole run; both orientations are built once and cached on the matrix rather than re-derived per call.
• 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.000    0.000     0.004    0.003    0.001    38.9   100.0   100.0
    200      1935   312.5KB    0.000    0.000     0.003    0.001    0.001    36.0   100.0   100.0
    400      7771     1.2MB    0.000    0.002     0.008    0.003    0.004    39.3   100.0   100.0
    800     31147     4.9MB    0.002    0.008     0.034    0.014    0.019    38.1   100.0   100.0
   2000    194913    30.5MB    0.108    0.157     0.695    0.371    0.323    38.2   100.0   100.0
   5000   1219064   190.7MB    2.722    3.094    14.461    8.295    6.165    38.2   100.0   100.0
- - - - - - - - - - - - - - - - - - - - - - - -
  extrapolated from the n=800..5000 rows (acdc ~ n^3.30, R2=1.000 in log-log) -- NOT measured
 ~10000  4.88e+06   762.9MB     45.5       29       142     94.5     53.4       -       -       -
 ~50000  1.22e+08    18.6GB  2.4e+04 5.31e+03  2.89e+04 2.53e+04 8.42e+03       -       -       -
~100000  4.88e+08    74.5GB 3.58e+05 5.01e+04  2.85e+05 2.81e+05 7.45e+04       -       -       -

extrapolated n	full acdc_match	one dense n × n array
10 000	~2.4 min	762.9 MB
50 000	~8 h	18.6 GB
100 000	~3.3 d	74.5 GB

• 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 3.30 with R² = 1.000. (It reads higher than the 2.97 measured before the inner loops were compiled, because removing a large constant per-element overhead makes the small rows — which are dominated by fixed costs — relatively cheaper, and steepens the fitted line. The asymptotic behaviour is unchanged.) 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. Two tests pin the compiled kernels to the definitions they stand in for: gradient_entries against the dense gradient, and make_swaps against a literal recompute-the-whole-score-per-trial reference.

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.