Max Flow via Electrical Flows, in depth: multiplicative weights over Laplacian solves, the energy certificate and a numpy implementation

The problem and the congestion view

The input is an undirected graph with capacities u_e, a source s and a sink t. A flow of value F sends F units out of s and into t, conserves flow at every other node and keeps |f_e| at most u_e on every edge. The max flow introduction covers the exact combinatorial algorithms. The goal here is weaker but scales differently: a flow of value at least (1 - eps) times the maximum, for a chosen eps.

Define the congestion of an edge as |f_e| / u_e. A flow is feasible exactly when the maximum congestion is at most 1. That reformulation matters, because max flow is now about controlling a maximum, the L-infinity norm of the congestion vector, while an electrical flow controls a sum of squares.

One electrical flow is not enough

Give each edge a resistance r_e. Among all flows of value F, the electrical flow is the unique one that minimises the energy, the sum of r_e f_e^2. It is found by solving L phi = F (e_s - e_t) for potentials phi, where L is the weighted graph Laplacian, and setting f_e = (phi_u - phi_v) / r_e. One linear solve gives the flow.

Take a diamond: edges s-a (capacity 3), s-b (2), a-b (1), a-t (2) and b-t (3), with the maximum flow equal to 5. With unit resistors and F = 5, symmetry sends 2.5 down each side and nothing across a-b; potentials are 5, 2.5, 2.5 and 0, and the energy is 25. That ignores capacities, so edges s-b and a-t are overloaded.

The natural fix sets r_e = 1 / u_e^2, which makes the energy equal to the sum of squared congestions. The solve then routes [2.647, 2.353, 0.294, 2.353, 2.647] across the five edges. Edges s-b and a-t still carry congestion 1.176. Minimising the sum of squares spreads load but cannot promise that no single edge exceeds 1. No fixed choice of resistances does that in general, which is why the algorithm iterates.

Multiplicative weights over resistances

The algorithm keeps a weight w_e per edge, starting at 1, and repeats. It sets r_e = (w_e + eps W / 3m) / u_e^2, where W is the total weight and m the number of edges, and computes the electrical flow of value F. Edges with congestion above 1 have their weight multiplied by (1 + (eps / rho) cong_e), so their resistance grows and the next flow avoids them. After N rounds it outputs the average of the flows. This is the multiplicative-weights framework: each electrical flow is a cheap oracle answer that is good on average, and the weights steer the oracle until the average is good everywhere.

The energy test gives the algorithm a certificate. If a flow of value F with congestion at most 1 exists, its energy under the current resistances is at most the sum of (w_e + eps W / 3m), which is (1 + eps / 3) W. The electrical flow has the smallest energy of all flows of value F, so its energy can only be lower. If the solve returns energy above (1 + eps) W, no feasible flow of value F exists, and the algorithm can stop and report that.

rho is the width, a bound on the congestion any single oracle answer can have. The paper shows that the eps W / 3m term keeps every resistance large enough that rho is O(sqrt(m / eps)), which bounds the number of rounds at O(rho log m / eps^2). The constants used below, rho = 3 sqrt(m / eps) and N = 2 rho ln m / eps^2, are one standard way of writing those bounds; treat them as a reading of the analysis rather than tuned values.

A numpy implementation

import numpy as np

def electrical_flow(n, edges, r, s, t, F):
    m = len(edges)
    B = np.zeros((m, n))                   # signed edge-vertex incidence matrix
    for e, (u, v) in enumerate(edges):
        B[e, u], B[e, v] = 1.0, -1.0
    L = B.T @ (B / r[:, None])             # weighted Laplacian B^T R^-1 B
    b = np.zeros(n); b[s], b[t] = F, -F
    keep = [i for i in range(n) if i != t] # ground the sink: phi_t = 0
    phi = np.zeros(n)
    phi[keep] = np.linalg.solve(L[np.ix_(keep, keep)], b[keep])  # dense: demo only
    f = (B @ phi) / r
    return f, float(np.sum(r * f * f))

def mwu_flow(n, edges, cap, s, t, F, eps=0.1, iters=None, rho=None):
    m = len(edges)
    rho = rho or 3 * np.sqrt(m / eps)
    iters = iters or int(np.ceil(2 * rho * np.log(m) / eps ** 2))
    w, total = np.ones(m), np.zeros(m)
    for k in range(iters):
        W = w.sum()
        r = (w + eps * W / (3 * m)) / cap ** 2
        f, energy = electrical_flow(n, edges, r, s, t, F)
        if energy > (1 + eps) * W:         # certificate: no flow of value F fits
            return None, k
        w *= 1 + (eps / rho) * (np.abs(f) / cap)
        total += f
    avg = total / iters
    return avg / max(1.0, np.max(np.abs(avg) / cap)), iters   # scale to feasibility

Two properties hold by construction. Every electrical flow conserves flow, so the average does too; in the 294-edge run the largest conservation error at an internal node was 2.6e-14. And dividing by the maximum congestion turns any conserving flow into a feasible one of proportionally smaller value, so the output is always a valid flow and the only question is how close its value is to F. The dense np.linalg.solve is O(n^3) per round and is there to keep the demo readable; the speed of the real method comes entirely from replacing it with a fast Laplacian solver.

Worked runs

On the diamond with F = 5 and eps = 0.1, the theoretical parameters give 6,829 rounds. The raw average is [2.983, 2.017, 0.967, 2.017, 2.983], close to the exact routing [3, 2, 1, 2, 3], with maximum congestion 1.0083; the function returns it scaled to [2.959, 2.000, 0.959, 2.000, 2.959], a feasible flow of value about 4.96. Asking for F = 5.5 instead trips the energy certificate in the very first round, correctly proving that 5.5 is infeasible.

The theoretical round count is conservative. A heuristic variant with rho = 2 and a fixed number of rounds falls outside the proof but shows the behaviour. On a 60-node, 294-edge random graph with capacities 1 to 10 and exact maximum flow 58, asking for F = 58 gave:

Asking for F = 63.8, ten percent above the optimum, with the same rho = 2 settings, produced the infeasibility certificate after 126 rounds; the round count depends on rho. Combined with binary search on F, the certificate and the scaled flow bracket the true maximum from both sides: every infeasible answer lowers the upper bound and every scaled flow raises the lower one.

Where the speed comes from

Each round is one Laplacian solve. Spielman and Teng showed in 2004 that symmetric diagonally dominant systems can be solved to good accuracy in near-linear time, and later solvers made that practical. With an O~(m) solver and O(sqrt(m) / eps^(5/2)) rounds, the simple algorithm above runs in O~(m^(3/2) eps^(-5/2)). The paper then removes the few edges that cause high congestion, which shrinks the width, and its abstract states the final bound as O~(m n^(1/3) eps^(-11/3)) for a (1 - eps)-approximate maximum s-t flow.

The idea kept going. Sherman (2013) and Kelner, Lee, Orecchia and Sidford (2014) reached almost-linear time for approximate undirected max flow using non-Euclidean gradient methods and oblivious routing. Madry (2013) used electrical flows inside an interior-point method to get O~(m^(10/7)) exact max flow on unit-capacity directed graphs. In 2022 Chen, Kyng, Liu, Peng, Probst Gutenberg and Sachdeva gave an almost-linear m^(1+o(1)) algorithm for exact max flow and min-cost flow on directed graphs, with an interior-point method driven by dynamic data structures for approximate min-ratio cycles rather than a Laplacian solve per step. Those results are theoretical milestones; none of them is what a production system should reach for first.

Operational guidance

• Know when not to use it. For exact answers on graphs that fit in memory, push-relabel or the algorithms in the ISAP article are simpler and faster in practice. The electrical approach is a fit for very large undirected graphs where an approximate value with a certificate is enough, or where sparse linear algebra on a GPU is the available hardware.
• Use a sparse solver. Replace the dense solve with conjugate gradients plus a preconditioner such as algebraic multigrid, and warm-start each round from the previous potentials, since weights change slowly.
• Rescale weights. Weights grow geometrically; divide them by their sum every round, since only ratios affect the resistances.
• Repair inexact solves. An approximate solve breaks conservation slightly. Measure the residual and fix the imbalance with a spanning-tree correction before scaling, or the output is not a flow.
• Get the cut separately. The method returns a flow value, not a minimum cut. For a cut, round the approximate flow with a few exact augmentations, or sweep thresholds over the final potentials; the max-flow min-cut article shows how an exact cut is read off the residual graph of a maximum flow.

Failure modes

• Directed graphs. Electrical flow is symmetric: current flows either way along a resistor. Applying this loop to a directed instance silently uses arcs backwards. Directed problems need the interior-point machinery instead.
• Extra components. The grounded Laplacian is singular whenever some component does not contain t, even one isolated vertex, so the solve fails although s and t are connected. Restrict the graph to the component containing s; if t is not in it, the answer is 0.
• Zero or tiny capacities. r_e = 1 / u_e^2 explodes. Delete zero-capacity edges and expect ill-conditioning when capacities span many orders of magnitude.
• Trusting the heuristic. A fixed small round count gives no guarantee; report the scaled value as a lower bound, never as the maximum.
• Misreading the certificate. The energy test proves infeasibility only when the solve is accurate. With a loose solver tolerance, leave a margin.

Trade-offs

What to do next

1. Implement the electrical flow solve for a small graph and check it against Ohm's and Kirchhoff's laws by hand.
2. Add the multiplicative-weights loop, run it on the diamond, and compare the average with the exact max flow from an augmenting-path solver.
3. Wrap it in binary search on F, using the energy certificate for the upper bound and the scaled flow for the lower bound.
4. Swap the dense solve for sparse conjugate gradients and measure rounds and time as the graph grows.
5. Use the reductions in the flow network reductions article to map your real problem to undirected s-t flow, and check it really is undirected.
6. Read the CKMST paper's Section 3 for the simple algorithm before moving on to the width-reduction refinements.