mixingmatrix

Compute the edge weights that make averaging, gossip or a random walk converge as fast as possible on a given graph.

This is the fastest-mixing Markov chain problem (Boyd, Diaconis & Xiao 2004) and its free-weight variant FDLA (Xiao & Boyd 2004).

pip install mixingmatrix

Quick start

import networkx as nx
import mixingmatrix

G = nx.random_geometric_graph(60, 0.25, seed=3)
sol = mixingmatrix.solve(G)

sol.weights          # scipy.sparse.csr_matrix, the mixing matrix
sol.to_dict()        # {(u, v): weight}, in your own node labels
sol.slem             # 0.946 — the mixing rate; lower is faster
sol.certified_gap    # 4.1e-04 — the true optimum is within this of sol.slem

print(sol.summary())

Two fields are worth reading every time:

• sol.status — "optimal" means the solver's stopping test was met. Anything else ("max_iter", "time_limit") means it stopped early. The matrix is still valid and usable.
• sol.certified_gap — how far from optimal you might be. Stays valid even when the solve stopped early.

Is it worth it on your graph?

mixingmatrix.compare(G)

       method    slem  spectral_gap  consensus_rounds  speedup  runtime
      optimal  0.9461       0.05392             249.3    2.29x    10.49
   metropolis  0.9761       0.02395             569.9     1.00     0.001
   max_degree  0.9817       0.01834             746.2     0.76     0.001
best_constant  0.9817       0.01834             746.2     0.76     0.001

speedup is the reduction in rounds to consensus against Metropolis-Hastings, which is what most gossip code already uses. Typical gains across graph families are 1.4x–2.3x, occasionally more on graphs with a bottleneck.

What to call, when

Situation	Call
Not sure it's worth doing	compare(G)
Up to ~500 nodes	solve(G)
Above ~1000 nodes, or an unfamiliar graph	solve(G, time_limit=60)
A few edges changed since last time	update(sol, changes=[...])
Need full precision, under ~100 nodes	solve(G, method="cvxpy")
Distributed averaging, not a Markov chain	solve(G, allow_negative=True)
Non-uniform target distribution	solve(G, stationary=pi)
Above 10,000 nodes	metropolis_hastings(G)

Set time_limit on anything large. Every iterate the solver produces is a valid mixing matrix, so a time limit returns the best one found rather than nothing.

Topology changes

If a few edges are added or removed, don't call solve again:

sol2 = mixingmatrix.update(sol, changes=[("remove", 3, 7), ("add", 1, 9)])
sol2.update_kind      # "certified_free" | "warm" | "cold"

Many single-edge changes provably don't move the optimum, and update skips the solve entirely when that's the case. On a 60-node graph, removing a zero-weight edge takes 0 solver iterations instead of 7745.

See docs/incremental.md.

The optimum is sparse

The optimal weights are often zero on a large fraction of edges, and those edges can be deleted without changing the mixing rate at all:

sol.zero_edges              # edges with weight 0
sol.support_fraction        # fraction of edges actually carrying weight
mixingmatrix.can_remove(sol, u, v).is_free    # True = provably safe to delete

On the 60-node example above, 98 of 260 edges get weight zero. If you are designing a network rather than weighting an existing one, this is often the more useful output.

API

# Solve
mixingmatrix.solve(G, method="auto", tol=1e-6, time_limit=None, ...)
mixingmatrix.update(sol, graph=None, changes=None)
mixingmatrix.compare(G, methods=(...))

# Other weightings, for comparison
mixingmatrix.metropolis_hastings(G)
mixingmatrix.max_degree(G)
mixingmatrix.best_constant(G)
mixingmatrix.lazy_random_walk(G, p=0.5)

# Measure any mixing matrix
mixingmatrix.slem(B)
mixingmatrix.spectral_gap(B)
mixingmatrix.consensus_rounds(B, eps=1e-6)
mixingmatrix.mixing_time(B, eps=1e-4)
mixingmatrix.is_valid(B, G)

# Run the chain
mixingmatrix.gossip(B, x0, rounds=50)
mixingmatrix.rounds_to_consensus(B, x0, eps=1e-6)

# Prove an edge change is safe
mixingmatrix.can_remove(sol, u, v)
mixingmatrix.can_improve(sol, u, v)

Command line:

mixingmatrix solve graph.edgelist --out weights.npz
mixingmatrix compare graph.edgelist
mixingmatrix info graph.edgelist

Full reference: docs/api.md.

Performance

One core, tol=1e-6, random 4-regular graphs:

Nodes	Time	Result
50	1.7 s	converged, gap 7e-05
100	3.2 s	converged, gap 2e-04
200	17 s	converged, gap 2e-03
400	77 s	hit iteration cap, gap 5e-03
800	41 s	hit iteration cap, gap 2e-02

Above 1000 nodes, use time_limit and read certified_gap. The solver runs up to 10,000 nodes in a few hundred MB, but does not fully converge at that size — you get a valid matrix better than the baselines, with a large but reliable gap estimate.

Graph structure matters more than size at the top end: the same 10,000-node solve is far cheaper on a grid than on a random regular graph.

If solves are unexpectedly slow, try OMP_NUM_THREADS=1. Multi-threaded BLAS can be dramatically slower on the small matrices this solver uses.

Limitations

• Graph must be connected and undirected.
• tol is a residual tolerance, not an error bound. Use certified_gap for accuracy.
• method="smoothing" cannot run above 2000 nodes.
• Non-reversible chains (symmetric=False) are not implemented.
• A better spectral gap does not always mean better decentralised-SGD performance (Vogels et al., NeurIPS 2022).

More detail: docs/limits.md.

Documentation

• Quickstart — a longer walkthrough
• Cookbook — recipes by situation
• API reference
• Incremental updates
• Limitations
• Internals — how it works, if you're curious

Runnable examples: examples/quickstart.py, examples/incremental.py.

Contributing

See CONTRIBUTING.md. Bug reports and feature requests go to the issue tracker. Participants are expected to follow the Code of Conduct.

Citing

See CITATION.cff. The papers this implements:

• Boyd, Diaconis & Xiao, Fastest Mixing Markov Chain on a Graph, SIAM Review 46(4), 2004.
• Xiao & Boyd, Fast linear iterations for distributed averaging, Systems & Control Letters 53(1), 2004.

MIT licensed.