network equilibrium, traffic assignment, maximum flow, algebraic multigrid, nonlinear multigrid, inexact Newton methods, continuation, graph Laplacian

1 Introduction↩︎

Flows on networks at equilibrium are everywhere: drivers distributing over a road map until no route is faster (traffic assignment [1]–[3]), packets routed to minimize delay in a communication network [4], [5], current in a resistor network [6], and—at the combinatorial extreme—the maximum flow through a capacitated graph [7]–[9]. These look like different problems with different algorithms: traffic assignment is solved by Frank–Wolfe or general convex optimization, max-flow by augmenting-path or push–relabel combinatorics. We make a unifying observation and build a single fast solver on it.

1.0.0.1 The observation

Each of these is the stationary point of a convex, edge-separable energy in the node potentials, and its Euler equation is one and the same nonlinear graph Laplacian, \[\label{eq:master} B\,\rho(B^{\top}\phi)=\alpha\,d ,\tag{1}\] where \(B\) is the node–edge incidence matrix, \(\phi\) the node potentials, \(d\) a source/sink demand, \(\alpha\) a load, and \(\rho_e\) a monotone edge law that alone encodes the physics: a saturating \(\rho_e\) gives a hard capacity (max flow); a gently rising \(\rho_e\) gives a congestion cost (traffic)³. The flow is \(f_e=\rho_e\big((B^{\top}\phi)_e\big)\), and every Newton linearization of 1 is an ordinary weighted graph Laplacian \(J=B\,\mathrm{diag}(\rho')B^{\top}\). For congestion and routing problems \(\alpha\) is a given input; for maximum flow, \(\alpha\) is an unknown to be maximized. To obtain a fast nonlinear solver, we use an inexact chord-Newton iteration with the linearization frozen across most steps, lazily refreshed. For maximum flow we weave a continuation in \(\alpha\) into the same Newton loop to track the solution through the fold. The nonlinear solve costs \(2\)–\(4\) linear Laplacian solves (§4.8). Any near-linear Laplacian solver can serve as the inner engine; we default to approximate Cholesky [10], the faster engine across the corpus at equal robustness (§4.7), with LAMG+ [11] interchangeable in robustness (its constant differs by graph class). One incidental convenience of LAMG+’s affinity-based aggregation is that it naturally respects the saturating min cut, though the framework does not depend on it (§3.1).

1.1 Contribution↩︎

We make three contributions, in order of practical weight.

1. A unifying resistor-network formulation (§2) that casts maximum flow, congestion (traffic) routing, and minimum-delay routing as one nonlinear Laplacian 1 , with the problem class fixed by a single monotone edge law and the easy/hard dichotomy (no-fold vs.fold; §2) fixed by whether that law saturates. The stationarity equation itself is not new—it is the optimality system of a convex, edge-separable network-flow program (monotropic programming [12]); our contribution is the unification, that one edge law parametrizes all three problems, and a single solver that addresses the whole class. Unlike the electrical-flow and IPM Laplacian paradigm for flow [13], [14], which drive an outer reweighting/barrier loop around linear solves to an \((1+\epsilon)\)-approximate optimum, here the nonlinearity sits in the edge law and the equilibrium reached is exact.

2. A linear-time solver for convex congestion flow (§4)—establishing two claims that serve as a building block for a larger congestion solver: (a) the convex Beckmann equilibrium is exactly 1 ; (b) a robust \(\mathcal{O}(m)\) near-linear Laplacian inner solve makes it linear-time as well. No combinatorial algorithm addresses this problem. NLF matches a state-of-the-art IPM to \(10^{-10}\) in a flat step count; on poorly-separable graphs NLF is the only \(\mathcal{O}(m)\) solver and wins past \(45\times\) (§4.4). Robustness is corpus-wide: \(2{,}003\) SuiteSparse graphs converge in time \(\propto m^{1.02}\) (§4.5, 4.8). A multicommodity extension (§6) carries \(K\) commodities through one hierarchy at \(\mathcal{O}(Km)\). We find that on real capacity data coupling inflates shared-corridor costs (\({\le}8.8\times\)) without rerouting.

3. An \(\mathcal{O}(m)\) maximum-flow solver (§5). The saturating instance recovers the exact combinatorial max-flow as a short, bounded sequence of Laplacian solves, with the cut respecting multigrid coarsening arising for free. Not faster than a tuned combinatorial solver, but returns the cut potential and full parametric flow curve.

Extensions are outlined in §7: directed multicommodity assignment, a full-multigrid variant, and related optimization problems: optimal transport and DC optimal flow.

1.2 Motivating applications↩︎

Computing an equilibrium flow on a large network is a recurring industrial task, and in each domain the computational bottleneck is the same Laplacian solve at the core of 1 . The production formulations of these tasks are typically directed (arc-based, with asymmetric per-direction costs); this paper studies the undirected relaxation as a feasibility study, leaving the directed models these domains deploy to future work (§7). The applications below are thus motivating problem classes; we don’t claim that NLF replaces an incumbent in its directed production setting.

Road traffic assignment. Metropolitan planning agencies compute user-equilibrium traffic assignment, i.e., how drivers distribute over a road map until no route is faster, to evaluate road projects and pricing [1]–[3]. Its undirected convex relaxation is the Beckmann program with the BPR congestion law, which is an NLF no-fold instance and the setting of our sharpest result (§4). The incumbent solvers are Frank–Wolfe and bush-based methods (§4.2).

Communication and data-center routing. Routing packets to minimize end-to-end delay under the Kleinrock \(M/M/1\) model [4], [5] is the same equilibrium with a capacity-saturating delay law—a fold instance handled by the continuation of §3.3.

Maximum flow (logistics, scheduling, vision). The capacitated maximum-flow problem [7], [8] is the saturating-law extreme of the same energy (§5). NLF recovers the exact value together with the min cut and the full flow-vs-load curve. Combinatorial solvers are faster to obtain the maximum flow value, but do not return the smooth interior flow or the parametric curve that an outer optimizer manipulates when max flow is an inner step.

Nonlinear resistor and analog-circuit networks. The most literal instance: monotone two-terminal devices (diodes, varistors, memristors) in DC steady state give 1 with \(\rho_e\) the device \(I\)–\(V\) law—the form of DC circuit analysis and in-memory-computing crossbars. As with OPF (§7.5) these topologies are typically structured, where direct solvers are efficient; NLF would have an edge only for large, irregular analog networks.

Across all applications, in addition to a scalar optimum (e.g., the maximum flow value) the equilibrium (flow), its binding-constraint structure (min cut), and its sensitivity to load are also required for design optimization, which is the gap NLF targets.

1.3 Prior work, and how NLF differs↩︎

Three lines of work meet at 1 ; NLF sits in the gap none of them fills.

The Laplacian paradigm for flow. Electrical-flow algorithms [13] and IPM for generalized flow [14] solve a sequence of linear weighted-Laplacian systems inside an outer multiplicative-weights or barrier loop, reaching an \((1+\epsilon)\)-approximate optimum; the flow is a linear function of the potential and all nonlinearity lives in the outer schedule. NLF instead places the nonlinearity in the constitutive law \(f=\rho(B^{\top}\phi)\), so a single energy encodes the capacity and the equilibrium reached is exact (§2).

Fast linear Laplacian solvers. Nearly-linear solvers for \(Lx=b\) such as combinatorial multigrid [15], support-tree/Spielman–Teng theory [16], approximate Cholesky [10], [17], and algebraic multigrid including LAMG [18] address the linear, unconstrained problem only. They are the engine, not the method: NLF calls a robust \(\mathcal{O}(m)\) Laplacian solver as its inner kernel but adds the nonlinear, capacity-constrained equilibrium layer on top. AMG has also been applied to the complex-valued non-separable AC power-flow equations [19]—a nonlinear system with no objective and a non-Laplacian Jacobian, solved by nonlinear FAS multigrid with a multiplicative coarse correction (no Newton, no fold); a different problem class from the real, edge-separable, SPD-Laplacian constrained flow program here.

Multilevel optimization. NLF is not multilevel optimization in the MG/OPT sense [20]–[22]: it does not build coarse optimization models or recurse the objective. While the equilibrium is nonlinear equation 1 , NLF employs a near-linear inner solve only as the linear inner solve inside a frozen-linearization chord-Newton continuation; the nonlinearity is carried by the edge law and the outer continuation, not by a coarse-grid objective. However, Sec. 7.5 outlines how NLF might be extended to a multilevel optimization paradigm when a cost functional is added.

Combinatorial and transport-specific solvers. Augmenting-path and push–relabel max-flow [8], [23] and the almost-linear-time exact algorithm [9]—a theoretical breakthrough with galactic constants, not a practical implementation—target directed max/min-cost flow; Frank–Wolfe [24] and bush-based assignment (Algorithm B [25], TAPAS [26]) target the smooth-cost traffic program. Each is specialized to one problem class. NLF’s contribution is orthogonal: a single solver, parametrized by one edge law, that spans all of them and inherits graph-class robustness from its multigrid core. We benchmark against representatives of each family where a runnable one exists (Ipopt [27] for congestion, Boykov–Kolmogorov for max flow) and position the rest honestly where it does not (Frank–Wolfe, bush methods; §4.2).

2 A unified resistor-network formulation↩︎

Let \(B\) be the \(n\times m\) signed node–edge incidence matrix of a connected graph (\(B_{ie}=+1\) if node \(i\) is the head of edge \(e\), \(-1\) if the tail). Write \(\phi\in\mathbb{R}^n\) for the node potentials, \(d\in\mathbb{R}^n\) for a balanced demand (\(\mathbf{1}^{\top}d=0\)), and \(\alpha\ge0\) for a scalar load. Each edge \(e\) carries a flow \(f_e\) driven by its potential difference \(g_e = (B^{\top}\phi)_e = \phi_i - \phi_j\) where \(i,j\) are the head and tail of \(e\). The physics of the medium is encoded in a smooth, strictly monotone edge law \(\rho_e:\mathbb{R}\to\mathbb{R}\): \[\label{eq:law} f_e = \rho_e(g_e).\tag{2}\] Different choices of \(\rho_e\) give different physical problems (Table 1): an Ohmic resistor has \(\rho_e(g)=g/R_e\); a saturating capacity has \(\rho_e(g)=c_e\tanh(g/c_e)\); a BPR congestion link has \(\rho_e\) equal to the inverse of the BPR travel-time function.

2.0.0.1 The equilibrium equation

Flow conservation at every node, \(Bf=\alpha d\), together with the constitutive relation 2 , gives the NLF master equation (1 ), a nonlinear system in the potentials \(\phi\) alone. This is the equation NLF solves.

2.1 Primal–dual structure and the co-energy↩︎

Because each \(\rho_e\) is monotone, it is the derivative of a strictly convex co-energy \(\Psi_e\): \(\rho_e=\Psi_e'\). Equation 1 is then the stationarity condition \(\nabla_\phi E=0\) of the unconstrained, convex dual energy \[\label{eq:dual} E(\phi)=\sum_e\Psi_e\!\big((B^{\top}\phi)_e\big)-\alpha\,d^{\top}\phi,\tag{3}\] Solving 1 is therefore equivalent to minimizing \(E\).

The dual energy \(E\) arises naturally from a primal convex program in the flow variables \(f\). Each edge has a primal energy (cost) \(\Phi_e(f_e)\), and the primal problem minimizes total cost subject to flow conservation: \[\label{eq:primal} \min_f\;\sum_e\Phi_e(f_e)\quad\text{s.t.}\quad Bf=\alpha d .\tag{4}\] The Lagrangian of 4 introduces node potentials \(\phi\) as multipliers on the conservation constraint; the KKT stationarity condition \(\Phi_e'(f_e)=(B^{\top}\phi)_e\) says the marginal cost equals the potential difference, i.e.\(t_e(f_e)=g_e\), or \(f_e=\rho_e(g_e)\). The marginal cost \(t_e=\Phi_e'=\rho_e^{-1}\) is thus the inverse edge law. Substituting into \(Bf=\alpha d\) recovers 1 .

The relationship between \(\Phi_e\) and \(\Psi_e\) is the Legendre–Fenchel duality [28]: \(\Psi_e=\Phi_e^*\), i.e. \(\Psi_e(g)=\sup_f\bigl[g\cdot f-\Phi_e(f)\bigr],\) and \(\Psi_e'=(\Phi_e^*)'=(\Phi_e')^{-1}=\rho_e\). In words: the co-energy \(\Psi_e\) is the convex conjugate of the primal cost \(\Phi_e\), and its derivative is the edge law \(\rho_e\). The potentials \(\phi\) are sometimes called co-potentials in the resistor-network literature [6] because they live in the dual (potential) space, conjugate to the primal (flow) space. Appendix 9 derives \(\Phi_e\) and \(\Psi_e\) for each application in Table 1.

2.1.0.1 The Newton linearization

Differentiating 1 with respect to \(\phi\) gives the Newton (Jacobian) matrix \[\label{eq:jac} J=B\,\mathrm{diag}\!\bigl(\rho'(B^{\top}\phi)\bigr)\,B^{\top},\qquad w_e=\rho'_e(g_e)=\Psi_e''(g_e)=1/t_e'(f_e)\ge0 ,\tag{5}\] a weighted graph Laplacian with solution-dependent conductances \(w_e\). Everything problem-specific lives in \(\rho_e\); the solver never sees anything but \(\rho_e\), \(\rho'_e\), and \(J\).

2.2 The fold dichotomy↩︎

The single feature that governs difficulty is whether the edge law saturates. Table 1 lists three standard instances.

• Soft congestion (rising \(\rho\), no saturation). The load \(\alpha\) is a given input. The cost rises with flow but imposes no hard cap, so \(\rho_e\) is unbounded and the law never saturates: \(\rho'_e>0\) at every finite flow (for BPR \(\rho'_e\) is large near free flow and decays as \(f\) grows, but never reaches \(0\)). Hence \(J\) is nonsingular at every equilibrium—no feasibility limit and no fold. Traffic assignment with BPR cost is of this type; on its bounded-flow equilibria \(J\) is well-conditioned, so 1 is solved directly at the target load, with no continuation (§4).

• Hard capacity (saturating \(\rho\), conductance \(\to0\)). \(\alpha\) is an unknown to be maximized. The flow cannot exceed a capacity \(c_e\), so \(\rho_e\) is bounded and \(\rho'_e\to0\) as the edge saturates. As \(\alpha\) rises to a feasibility limit, the saturating edges form a cut and \(J\) folds—a continuation limit point where \(\alpha\) caps at the min cut \(F^*\) and \(J\) becomes singular. Maximum flow and minimum-delay communication routing (Kleinrock delay) are of this type; they require arclength continuation through the fold (§3.3).

2.3 The max-flow law in detail↩︎

For the saturating instance used in §5, take \[\label{eq:rho} \rho_e(g)=\begin{cases} c^+_e\tanh(g/c^+_e), & g\ge0,\\[2pt] -c^-_e\tanh(g/c^-_e), & g<0,\end{cases} \qquad \Psi_e(g)=c_e^{2}\log\cosh(g/c_e),\tag{6}\] a soft-capacity barrier: \(\Psi_e(g)=\tfrac12 g^2\) for \(|g|\ll c_e\) (a unit resistor) and grows only linearly, \(\Psi_e\to c_e|g|\), for \(|g|\gg c_e\) (a saturated edge), so \(\rho_e=\Psi_e'\) clamps the flow smoothly to \(-c^-_e\le f_e\le c^+_e\). Here \(w_e=1-(f_e/c_e)^2\in[0,1]\) represents “edge openness,” \(1\) for a slack edge and \(0\) for a saturated one, and the forming min-cut is exactly the zero-conductance set, which is why AMG aggregation driven by a relaxation-based strong-connection measure refuses to coarsen across it (§5). The congestion (BPR) law is given where it is used, in §4.

2.4 Relation to prior Laplacian-paradigm formulations↩︎

Solving flow through graph-Laplacian systems is the core of a productive line of work, but NLF’s mechanism is fundamentally different. Electrical-flow methods [13] and their \(p\)-Laplacian descendants compute a sequence of linear electrical flows \(f=\mathrm{diag}(w)\,B^{\top}\phi\) —each a fixed weighted-Laplacian solve—and steer an outer multiplicative-weights loop to a \((1+\epsilon)\)-approximate max-flow; IPM for (generalized) flow [14] likewise wrap a Laplacian solve in a barrier path-following loop. NLF instead places the nonlinearity in a constitutive law: \(f=\rho(B^{\top}\phi)\) is a saturating function of the gradient, so the single energy of 1 encodes the capacity exactly and its feasibility limit is the exact min cut. The outer loop is then a short continuation through one scalar fold, not a reweighting or barrier schedule. The Beckmann congestion instance 4 is already known [1]; the novelty is the saturating-resistor reformulation of max-flow, and placing max-flow, congestion, and minimum-delay routing under one nonlinear Laplacian solved via existing linear Laplacian algorithms.

3 The NLF algorithm↩︎

The two problem classes need different outer loops. For congestion (no fold), \(\alpha\) is a given input, so 1 is solved directly at that load by a handful of inexact chord-Newton steps (Algorithm 2), warm-started from \(\phi=0\); no load ramp is needed—a single fixed-load solve converges across the corpus (§4.5). For max-flow (fold), \(\alpha\) is the unknown maximizer, and continuation from \(\alpha=0\) in pseudo-arclength parametrization tracks the solution curve through the fold to \(F^*\), each fixed-\((\phi,\alpha)\) subproblem again a handful of chord-Newton steps. The inner solver’s setup—its factor or hierarchy—is built once, frozen across steps and re-formed only when an observed-convergence monitor signals it has drifted too far. The numbers of steps and inner solves are both bounded, implying empirical-\(\mathcal{O}(m)\) complexity. We describe the linear engine (§3.1), the fixed-load solve (§3.2), and the continuation (§3.3).

3.1 The inner Laplacian solve↩︎

The outer loop treats the inner solver as a black box that is required to do three things: (i) it solves a graph-Laplacian system \(L\phi=b\) in near-linear \(\mathcal{O}(m)\) work; (ii) its accuracy is tunable, so it can match the \(\eta\|r\|\) forcing tolerance of the chord-Newton step (§3.2) rather than to machine precision; and (iii) its setup is separable from its solve, so the setup can be built once at a frozen linearization and reused across many continuation steps.

Two public near-\(\mathcal{O}(m)\) solvers meet all three with different internals. Approximate Cholesky [10] sets up a randomized incomplete sparse-Cholesky factor preconditioner and solves by preconditioned conjugate gradients to the requested tolerance. LAMG+ [11] sets up an aggregation-multigrid hierarchy (relaxation-based affinity aggregation, low-degree Schur elimination, a guarded coarse operator) and solves by a multigrid cycle with min-residual recombination. We use approximate Cholesky as the default as it is the faster default across the corpus at equal robustness (§4.7). LAMG+ is as robust, while its wall-clock constant differs by graph class.

While both inner solvers yield the same overall nonlinear convergence (§4.5), one convenient property of LAMG+’s relaxation-based aggregation is that it is cut-respecting: as the cut conductances vanish they decouple relaxation across the cut, so test vectors decorrelate and the affinity declines to aggregate across it. Its measured Asymptotic Convergence Factor (ACF) of \(\approx0.01\) thus holds through the continuation with no fold-aware tuning (§5). This is a specific feature of LAMG+, not a framework requirement: adequate refresh of the frozen operator near the fold suffices for any inner solver, and over-freezing degrades all of them equally. (\(\alpha\), the continuation, and the fold never enter the linear solve; the singular cut direction near \(F^*\) is deflated outside the engine, §3.3.)

3.2 The fixed-load solve: a frozen-setup inner solve↩︎

Eq. 1 with fixed \(\alpha\) is solved by inexact chord-Newton iteration (Algorithm 2). At each step we evaluate the residual \(r(\phi)=\alpha d-B\rho(B^{\top}\phi)\) and solve for the Newton direction \(\delta\approx J_H^{-1}r\) using the inner Laplacian solver to tolerance \(\eta\|r\|\), \(\eta=0.05\) [30]. The Jacobian \(J_H\) and the expensive setup hierarchy are only recalculated when chord Newton’s linear convergence rate rises above \(0.25\). This addresses the no-fold congestion class; we verified that \(\phi=0\) is always a good initial guess there.

3.2.0.1 The damped update

The line search (step 5) is standard backtracking on the residual norm [31]. Because \(\delta\) comes from a frozen linearization, a full step can overshoot right after a load increase where the law’s curvature is strongest; backtracking by halving and accepting the first non-increasing residual guarantees monotone progress at negligible cost. Near the solution the full step is always accepted (one trial, \(\tau=1\)).

3.2.0.2 When the setup is refreshed

Three conditions trigger: (i) no setup yet; (ii) per-step residual factor \(>0.25\) (frozen \(J\) has drifted); (iii) line-search stall on a stale \(H\)—rebuild and retry once; a fresh-\(H\) stall breaks the loop. The corpus median refresh count is small (often one). A lazy refresh is \(1.5\)–\(4.7\times\) faster than rebuilding every step.

3.2.0.3 Deep hierarchies (LAMG+ engine only)

On graphs with \({\ge}20\)-level hierarchies (country-scale roads, large FEM meshes) the dipole right-hand side excites the longest mode that the stock cycle under-treats. A grown cycle index (\(1.15\times\) per level, cap \(\gamma_l\tau_l\le0.95\), each cycle stays \(\mathcal{O}(m)\)) is used, restarted if the factor exceeds \(0.6\). This engages on \(11\) of \(2{,}003\) graphs (\(0.5\%\)) and restores textbook convergence.

3.3 Continuation↩︎

Only the saturating/fold class needs continuation. Here \(\alpha\) is the unknown maximizer. A fixed-\(\alpha\) solve is feasible iff \(\alpha<F^*\) (the min-cut value), and as \(\alpha\uparrow F^*\) the solution curve folds—\(\alpha(V)\) flattens onto the asymptote \(F^*\) while the potentials grow without bound (Fig. 4); at the fold \(\kappa(J)\sim(F^*-\alpha)^{-1}\) and the source–sink gap is \(V=|\phi_s-\phi_t|\sim-\log(F^*-\alpha)\), so \(dV/d\alpha\sim(F^*-\alpha)^{-1}\to\infty\). Thus \(\alpha\) is a bad continuation parameter: the feasible loads crowd at the fold and the step is forced to zero. The well-conditioned variable is the cut-mode amplitude \(\psi=\chi^{\top}\phi\sim-\log(F^*-\alpha)\),⁴ \(\chi\) the (a priori unknown) min-cut indicator. We realize this by pseudo-arclength continuation [32], [33] (Algorithm 3), which never forms \(\chi\).

The equilibria of 1 trace a smooth curve \(\big(\phi(s),\alpha(s)\big)\) as the load is varied; we parametrize it by its own arclength \(s\) and write a dot for the derivative along the curve, \(\dot{(\cdot)}\equiv d(\cdot)/ds\), so the tangent \((\dot{\phi},\dot{\alpha})\) is a unit vector, \(\|\dot{\phi}\|^2+\dot{\alpha}^2=1\). Differentiating 1 along \(s\) gives: \[\label{eq:tangent} J\dot{\phi}=\dot{\alpha}\,d\quad\Longrightarrow\quad \dot{\phi}=J^{+}d,\;\;\dot{\alpha}=1,\;\;\text{then rescale to }\|(\dot{\phi},\dot{\alpha})\|=1 .\tag{7}\] As \(\alpha\to F^*\), \(J^{+}d\propto\chi/\lambda_{\min}\) (the pseudoinverse \(J^{+}\) acts on \(\mathbf{1}^{\perp}\), so \(\lambda_{\min}\) is again the smallest nonzero eigenvalue), so the tangent rotates onto the cut mode by itself—stepping in \(s\) becomes stepping in \(\psi\), the well-conditioned direction.

Each continuation step is a predictor that advances a fixed arclength \(\Delta s\) along the tangent, \((\phi,\alpha)\mathrel{+}=\Delta s\,(\dot{\phi},\dot{\alpha})\), then a corrector that drives two residuals to zero: the equilibrium residual \(r=\alpha d-B\rho(B^{\top}\phi)\) of 1 (as in Algorithm 2) and the arclength residual \[\label{eq:arclength} N=\dot{\phi}^{\top}(\phi-\phi_0)+\dot{\alpha}\,(\alpha-\alpha_0)-\Delta s ,\tag{8}\] how far the point has advanced along the tangent past the previous solution \((\phi_0,\alpha_0)\), minus the target \(\Delta s\). Setting \(r=0\) holds the point on the equilibrium curve; \(N=0\) fixes how far along it (hence “pseudo”-arclength: the exact arclength sphere \(\|\delta\phi\|^2+\delta\alpha^2=\Delta s^2\) is linearized to its tangent plane \(N=0\).) One Newton step on \((r,N)=0\) is the bordered linear system \[\label{eq:bordered} \begin{bmatrix} J & -d\\ \dot{\phi}^{\top} & \dot{\alpha}\end{bmatrix} \begin{bmatrix}\delta\phi\\ \delta\alpha\end{bmatrix}=\begin{bmatrix}r\\ -N\end{bmatrix}.\tag{9}\] Block elimination of \(\delta\alpha\) gives \[\label{eq:deflate} \delta\alpha=-\,\frac{\dot{\phi}^{\top}w+N}{\dot{\alpha}+\dot{\phi}^{\top}u},\qquad \delta\phi=\delta\alpha\,u+w ,\tag{10}\] where \(u=J^{+}d\) and \(w=J^{+}r\) are computed by two applications of the inner Laplacian solver. Both \(u,w\in\mathrm{range}(J^{+})=\chi^{\perp}\): the near-null cut component of the solution (the \(\chi\) direction) is pinned by the scalar arclength row 8 , not by these solves. The border keeps the system nonsingular through the fold (\(\kappa=\mathcal{O}(10^3)\) while \(\kappa(J)\to\infty\)), so the inner Laplacian solver stays fast (approximate Cholesky PCG and LAMG+ cycle ACF both \(\approx0.01\)).

The refresh policy of §3.2 is woven into the corrector as a staleness safeguard (Algorithm 3). Both \(u\) and \(w\) computations reuse one frozen hierarchy; since the cut conductances drain exponentially in the cut amplitude, that hierarchy ages faster as the fold sharpens, and the monitor rebuilds it more often — about every 4 arclength steps at \(0.97 F^*\), every 2 at \(0.995 F^*\). The predictor’s tangent \(\dot{\phi} = J^{+}d\), the single direction most sensitive to \(J\)’s drift, is instead recomputed on a fresh hierarchy once per accepted step. The fixed-load solution remains machine-precise as close to the fold as \(\alpha = 0.995 F^*\) (a fixed-\(\alpha\) equilibrium is well-posed only for \(\alpha < F^*\)), while the arclength continuation tracks through the fold to the exact F* without stalling (Table 7).

4 Transportation: undirected congestion equilibrium↩︎

This is NLF’s sharpest result: on convex congestion flow—which has no combinatorial competitor and whose operator never folds—NLF is an empirically linear-time solver, verified against a state-of-the-art interior-point method on real road-network topologies.

4.0.0.1 Scope: an undirected, single-commodity

Throughout this section the congestion problem is the undirected, single-commodity Beckmann relaxation: flows are signed, the cost depends on \(|f_e|\), and Wardrop’s path-nonnegativity is dropped (App. 9). Our experiments establish two claims: (i) this undirected congestion equilibrium is exactly 1 for the inverse-BPR law, and (ii) a robust \(\mathcal{O}(m)\) linear graph Laplacian inner solve makes it linear-time. The real road networks supply real topologies, capacities, and BPR data exercised with a single source–sink demand; we cross-check NLF against an independent high-accuracy solver on the identical convex program.

4.1 The BPR congestion law↩︎

On a congested road the travel time rises with the flow it carries. The field-standard model is the Bureau of Public Roads (BPR) cost [2], \[\label{eq:bpr} t_e(f)=t^0_e\Big(1+b\,(f/c_e)^{p}\Big),\qquad b=0.15,\;p=4 ,\tag{11}\] with \(t^0_e\) the free-flow travel time and \(c_e\) the practical capacity. The undirected Beckmann congestion equilibrium (signed flows; Wardrop’s path-nonnegativity relaxed, App. 9) is the convex program 4 with \(\Phi_e(f)=\int_0^f t_e=t^0_e f+\tfrac{b\,t^0_e}{p+1}f^{p+1}/c_e^{p}\); its edge law is the inverse marginal cost \(\rho_e=t_e^{-1}\), with conductance \(\rho'_e=1/t_e'(f)=\big(t^0_e b\,p\,(f/c_e)^{p-1}/c_e\big)^{-1}>0\) at every finite flow: the cost rises but never caps the flow (the law does not saturate), so there is no feasibility limit and no fold.⁵ The Jacobian \(J=B\,\mathrm{diag}(\rho')B^{\top}\) is a well-conditioned Laplacian at every bounded-flow iterate, and NLF runs the fixed-load Algorithm 2.

4.2 No combinatorial competitor; an interior-point baseline↩︎

No combinatorial algorithm addresses the convex congestion program 4 . Frank–Wolfe [24] converges only sublinearly on the Beckmann program: its \(\mathcal{O}(1/k)\) rate (\(k\) = iteration count) drives the iteration count to even a \(10^{-4}\) duality gap up with problem size and cannot reach the high accuracy NLF attains—the classical slow-tail limitation of the method on traffic assignment [3]. The field’s high-accuracy standard, the bush-based family (Algorithm B [25], TAPAS [26]), solves the directed user equilibrium (nonnegative arc flows, one-sided cost)—a different convex program from the undirected, signed-flow Beckmann relaxation NLF solves here (App. 9), so a same-instance head-to-head is ill-posed. We therefore position bush methods qualitatively and benchmark quantitatively against Ipopt [27], a mature IPM solver on the identical program. This is conservative: Ipopt attains the same \(10^{-10}\), whereas Frank–Wolfe and bush-based methods do not address the max-flow instance (§5). Ipopt’s per-step cost is a sparse-direct KKT factorization: cheap on planar graphs but fill-in grows superlinearly as separators worsen—exactly the axis on which NLF is immune.

4.3 Real road data: NLF and Ipopt↩︎

We use the Transportation Networks for Research benchmark [34], the field-standard repository of real metropolitan road-network topologies with real capacities, free-flow times and BPR parameters. Across instances spanning three decades of size, NLF and Ipopt reach the identical equilibrium of the undirected single-commodity Beckmann program (Table 2): the Beckmann objective agrees to the solver tolerance and the equilibrium link flows to \(\|f_{\mathrm{NLF}}-f_{\mathrm{Ipopt}}\|/\|f_{\mathrm{Ipopt}}\|\sim10^{-10}\), in \(\approx9\) Newton steps with no size trend.

4.4 Linear scaling, and the direct-vs-iterative crossover↩︎

Figure 5 and Table 3 time both solvers on the real road networks (planar, good separators) and on synthetic Erdős–Rényi graphs carrying BPR parameters (poorly separable), to \(2.4\times10^5\) edges. NLF’s per-edge cost is flat on both families—empirically near-linear (the full-corpus fit of §4.5 gives \(t\propto m^{1.02}\)); neither the step count (\(\approx9\)) nor the inner linear solver ACF sees the graph’s separability, while Ipopt does. On the real road networks its sparse-direct core is near-linear and beats NLF by \(\approx3\times\); planar graphs are where direct methods excel. On poorly-separable graphs the KKT fill-in explodes—\(t\propto m^{2.2}\)—so the gap inverts and widens with size: NLF is \(5\times\) faster at \(m=3.0\times10^4\), \(45\times\) at \(m=1.2\times10^5\), and at \(m=2.4\times10^5\) Ipopt does not finish within its \(120\) s CPU budget while NLF finishes in \(3\) s.

4.5 Robustness across the graph corpus↩︎

The road networks above exercise NLF on real flow topologies; we now show it is robust across graph classes in full. We impose a single-commodity BPR congestion instance (random capacities and free-flow times, a far source–sink demand) on every graph of a \(2{,}003\) real-world SuiteSparse collection, each reduced to its largest connected component. This includes spanning structured grids, finite-element and circuit meshes, web, social, road, and citation networks with up to \(1.8\times10^7\) edges.

Under this protocol (one BPR instance per graph, random capacities, a single far source–sink demand, stock settings) NLF converged on all \(2{,}003\) graphs to \(10^{-9}\) residual tolerance, with a Newton step count median of \(6\) and no size trend (Fig. 6B). The wall-clock fits \(t\propto m^{1.02}\) over the \(1{,}788\) graphs above \(10^3\) edges (Fig. 6A), per-edge cost flat at \({\approx}0.86~\mu\)s/edge median (\(95\)th percentile \(3.7~\mu\)s/edge). Three giants confirm the scaling beyond the sweep’s cap, to hollywood-2009 at \(5.6\times10^7\) edges.

4.6 Comparison with other convex optimizers.↩︎

We compare three solver paradigms for the same convex BPR program: near-linear (NLF), sparse-direct (Ipopt), and matrix-free first-order (L-BFGS on the dual 3 ; on well-conditioned graphs L-BFGS also scales near-linearly and beats Ipopt, so that win is a fill-in effect any matrix-free method inherits). We ran all three on the \(1{,}669\) corpus graphs with \(n\le300{,}000\) (the subset where Ipopt’s sparse-direct KKT does not run out of memory on poorly-separable instances), each competitor capped at \(5\times\) NLF’s wall-clock (Table 4). NLF converged on all \(1{,}669\); first-order failed on \(481\) (ill-conditioning, \(\kappa\sim N\)) and interior-point on \(143\) (fill-in). On the \(90\) graphs hard on both axes, neither competitor finished within budget. Where they converge, NLF is median \(2.6\times\)/\(4.2\times\) faster than Ipopt/L-BFGS and fastest-or-only on \(85\%\). Timing protocol: all three solvers share one global precompilation pass, and NLF is then timed on a warm per-graph call; the competitors on their first per-graph call. This per-graph warm-up favours NLF by a small-graph-weighted margin that is negligible at scale, where NLF’s advantage is the genuine avoidance of direct fill-in (§4.4); the “within budget” counts use the \(5\times\)-NLF wall-clock cap and are a wall-clock criterion, not a claim the others cannot converge given more time.

4.7 Inner-solver ablation: which near-linear engine↩︎

The framework can use any near-linear inner solver (§3.1); we settle the default empirically. Holding NLF’s outer method fixed (load-continuation chord-Newton, lazy refresh), we swap only the inner Laplacian solve and rerun the congestion corpus of §4.5 (one BPR instance each; over \(1{,}800\) graphs except the very largest graphs). The candidates are the two public \(\mathcal{O}(m)\) solvers of §3.1, approximate Cholesky and LAMG+; we exclude sparse-direct, which is not near-linear—its separability crossover against an \(\mathcal{O}(m)\) inner is the scaling result of §4.4, not an inner we would ship.

Both converge on every graph, to the identical equilibrium: the dual energy agrees to floating-point precision and the link flows to \(\le2\times10^{-8}\) relative error. So robustness is therefore not a discriminator. Speed splits by graph class, both ways: approximate Cholesky is faster on \(81\%\) of graphs (median \(1.8\times\), up to \(43\times\)—its cheap factor wins on low-treewidth mesh and finite-element graphs) and LAMG+ on the other \(19\%\) (up to \(24\times\)—its hierarchy wins on the poorly-separable and structural graphs); the overall median is \(1.6\times\) in approximate Cholesky’s favour, with \(39\%\) of graphs within \(1.5\times\) either way (Table 5). We therefore adopt the faster approximate Cholesky as the default and keep LAMG+ as a robust alternative: it wins wall-clock on the poorly-separable and structural minority, and at the saturating fold its cut-respecting aggregation (§3.1) gives a structurally clean near-singular solve with no fold-aware tuning, though approximate Cholesky’s constant is often smaller even there (§5.1).

Residual tolerance criterion sensitivity is mild (Table 6): freezing the factor more aggressively trades a few extra outer steps for fewer setups until it over-freezes, and a looser inner forcing \(\eta\) costs outer steps; the stock \((\text{refresh}=0.25,\;\eta=0.05)\) sits in the flat interior.

4.8 Cost as a function of the target accuracy↩︎

The dependence on the requested solution accuracy \(\varepsilon\) is logarithmic. Each inner solve to relative residual \(\delta\) costs \(\mathcal{O}(m\log(1/\delta))\) (a bounded inner contraction gives \(\mathcal{O}(\log(1/\delta))\) inner iterations). The outer iteration is a frozen chord-Newton with inexact-Newton forcing \(\eta_k\) tracking the nonlinear residual \(\|r_k\|\) (§3.1); its per-step cycle counts form a series dominated by the final solve (where \(\|r_k\|\!\to\!\varepsilon\)), the earlier looser steps adding only a bounded geometric tail, while the continuation and globalization run at a fixed, \(\varepsilon\)-independent path tolerance. Each contribution is thus a bounded multiple of a single linear solve, so the total cost to accuracy \(\varepsilon\) is, empirically, \[\label{eq:eps} \mathcal{O}\!\big(m\,\log(1/\varepsilon)\big).\tag{12}\] The constant is the whole-nonlinear-solve-to-one-linear-solve ratio reported below (\(2\)–\(4\)), and it stays bounded because the frozen chord-Newton contraction stays bounded away from \(1\): on the no-fold congestion class directly (§4), and through the fold because the continuation deflates the diverging cut direction and keeps the bordered corrector uniformly conditioned (§3.3); empirically the contraction does not degrade at \(F^*\). We verify this by sweeping \(\varepsilon\) from \(10^{-2}\) to \(10^{-12}\) on one graph per class (web, social, road, finite-element): the total inner-solve count grows linearly in \(\log(1/\varepsilon)\), at \(0.24\)–\(0.60\) per decade of accuracy, confirming 12 .

4.8.0.1 The nonlinear/linear cost ratio is \(\mathcal{O}(1)\)

Brandt’s yardstick says: solving the nonlinear problem should cost only a small, size-independent multiple of solving one linear problem of the same kind [35]. We measure exactly this: the full congestion solve (load continuation, all steps, all cycles, lazy refresh) against a single inner setup-and-solve of the linearization at \(\alpha=0\) (the linear resistor network \(B\,\mathrm{diag}(\rho'(0))B^{\top}\phi=d\)) on the same graph. Across the classes the ratio is \(1.9\)–\(4.4\)—grid \(3.1\), finite-element \(2.5\), road \(1.9\), social \(3.8\), web (\(1.9\times10^6\) edges) \(4.4\)—flat in size and class.

4.9 How accurately to solve each corrector step at the fold↩︎

At the fold, continuation is mandatory (\(J\) singular at \(F^*\)), but each arclength corrector need only stay on the solution path [35]—a loose per-step tolerance reaches the same min cut in \(2\)–\(7\times\) fewer corrector solves and hierarchy rebuilds, while a single corrector step per arclength step overshoots \(F^*\).

5 Maximum flow↩︎

The saturating instance 6 recovers the exact combinatorial maximum flow. Read \(f=\rho(B^{\top}\phi)\) as an edge-flow vector: \(\rho\)’s range is the capacity box \(|f_e|\le c_e\), and \(Bf=\alpha d\) is conservation, so any \(\phi\) solving 1 is a capacity-feasible, conserved flow of value \(\alpha\). The recession slope of \(\Psi_e\) is \(c_e\), so \(E\) has a finite minimizer iff \(\alpha<F^*\); as \(\alpha\uparrow F^*\) the min-cut edges saturate (\(\rho'_e\to0\) there only) and \(\phi\) develops a cut-indicator step. By max-flow/min-cut duality the largest feasible \(\alpha\) is \(F^*\), the min-cut capacity, and \(\phi\) is the LP dual. Because \(f=\rho(B^{\top}\phi)\) is a nonlinear function of the gradient, it realizes the true max-flow; a linear resistor network computes only the sub-maximal electrical flow (\(0.02\)–\(0.20\,F^*\) on heterogeneous capacities, except single-bottleneck instances where the two coincide). The nonlinearity is essential.

5.0.0.1 Cut-respecting coarsening

As \(\alpha\uparrow F^*\) the linearization \(J\) becomes singular along a single near-null direction—the min-cut indicator—so \(\kappa(J)\to\infty\). We assume a simple fold: the min cut is unique, so the limiting near-null space is one-dimensional (a degenerate cut of multiplicity \(k\), e.g.in symmetric unweighted graphs, would make it \(k\)-dimensional and require a rank-\(k\) block border in place of the scalar one below). The continuation of §3.3 deflates that one direction to a scalar outside the engine (the bordered/pseudo-arclength system); what the inner solver sees is the deflated complement (Eq. 10 ). The deflation is the enabling property—it is what keeps the inner operator well-conditioned through the fold, and it is a property of the framework, not of any one inner solver (swapping multigrid for a sparse-direct or another near-linear Laplacian solve reaches the identical \(F^*\); §5.1). Multigrid’s relaxation-based affinity is a natural fit on the complement—since the forming cut is the zero-conductance set of \(J\) (§2), the affinity declines to aggregate across it, so the coarse hierarchy mirrors the cut without being told where it is—but it is helpful, not required. The measured multigrid factor on the deflated operator is \(\mu\approx0.01\), uniformly in \(\alpha/F^*\) down to \(0.999\,F^*\) (Fig. 7); the uniformity is a property of the deflation, not of aggregation alone—on the undeflated \(J\), \(\mu\) degrades with \(\kappa(J)\) as expected.

5.0.0.2 Exactness and scaling

NLF returns the exact max-flow (Table 7), checked against \(F^*\) from a linear program over free flows and from push–relabel (which agree to \(10^{-4}\)), on every graph, to solver tolerance. The pseudo-arclength step count is flat in size (Fig. 7B, mean \(6.7\)) and the total wall-clock is empirically \(\mathcal{O}(m)\): \(t\propto m^{1.00}\) over the max-flow sweep (to \({\sim}10^4\) edges, Fig. 7A). Calibration, not a contest: Boykov–Kolmogorov, on its target vision grids, runs at \(\approx0.17~\mu\)s/edge (Fig. 7A), far below NLF’s constant—NLF is not meant to be a faster combinatorial max-flow solver. Its value is being \(\mathcal{O}(m)\), returning the cut potential (LP dual), the full flow-vs-load curve, and a smooth interior flow—the objects an outer method manipulates when max-flow is an inner step—one instance of a framework that, on the congestion class, has no combinatorial competitor at all.

5.1 Inner solver ablation↩︎

The deflation, not the choice of inner solver, is what carries the continuation through the fold (§5); the inner Laplacian solve is therefore a replaceable module. We check this directly: holding NLF’s outer method fixed (chord-Newton with the deflated pseudo-arclength continuation to \(F^*\)), we swap only the inner solver—multigrid (LAMG+), approximate Cholesky [10], an exact sparse Cholesky factorization, and a diagonally preconditioned CG—on max-flow instances spanning FEM, structural, and circuit graphs (Table 8).

The result is modularity, not a winner: every inner solver that converges returns the identical \(F^*\) to solver tolerance, so the continuation is genuinely solver-agnostic. The distinction that matters is near-linear versus not: the exact factorization is fastest where it fits but exhausts memory past \({\sim}150\)k poorly-separable nodes (G2_circuit, scircuit), exactly the regime where a near-linear inner solve is mandatory, while the diagonal-PCG baseline—no coarse correction—is typically up to two orders of magnitude slower and the least scalable. Among the near-linear solvers, approximate Cholesky and LAMG+ are interchangeable in robustness; approximate Cholesky often carries the smaller constant. The single stiff case (bcsstk38) is instructive: there all iterative inner solves stall and only the exact factorization reaches \(F^*\)—the deflation removes the global cut singularity but not the residual ill-conditioning of the deflated system, for which an exact inner solve is still required.

6 Multicommodity congestion↩︎

§4 established the formulation and the \(\mathcal{O}(m)\) solver for a single commodity. This section extends both, end to end, to \(K\) simultaneously routed commodities: a vector edge law that couples the commodities through shared congestion (§6.1), a solver that reuses Algorithm 2 and one scalar inner setup (factor or hierarchy) for all \(K\) commodities (§6.2), validation against an exact-Newton baseline (full \(Kn\times Kn\) block solve, §6.3), and the same full-corpus robustness sweep as §4.5 (§6.5). The implementation is a modular layer over the single-commodity solver, which remains intact.

6.1 Formulation: a vector edge law↩︎

\(K\) commodities ship demands \(\alpha d^1,\dots,\alpha d^K\) (\(\mathbf{1}^{\top}d^k=0\)) over one network; commodity \(k\) has flow \(f^k\in\mathbb{R}^m\) and potential \(\phi^k\in\mathbb{R}^n\). Collect the per-edge flows into the vector \(\mathbf{f}_e=(f^1_e,\dots,f^K_e)\in\mathbb{R}^K\) and the per-edge potential gradients into \(\mathbf{g}_e=\big((B^{\top}\phi^1)_e,\dots,(B^{\top}\phi^K)_e\big)\). The modeling question is how the commodities couple through shared congestion. Three natural choices: the signed total \(\sum_k f^k_e\) (degenerate—collapses to single-commodity); the volume \(\sum_k|f^k_e|\) (classical but nonsmooth, forces a per-edge complementarity akin to the max-flow cut); and the Euclidean magnitude \(\|\mathbf{f}_e\|\) (smooth, strictly convex, reduces exactly to \(K=1\)). We take the third for this feasibility study; the volume-coupled and directed models remain future work (§7.1). The program is 4 with the congestion cost charged to the magnitude of the commodity-flow vector, \[\label{eq:mcprimal} \min_{f^1,\dots,f^K}\;\sum_e\Phi_e\big(\|\mathbf{f}_e\|\big)\quad\text{s.t.}\quad Bf^k=\alpha d^k,\quad k=1,\dots,K ,\tag{13}\] and stationarity gives, per edge, the vector law \[\label{eq:mclaw} \mathbf{f}_e=\rho_e\big(\|\mathbf{g}_e\|\big)\,\frac{\mathbf{g}_e}{\|\mathbf{g}_e\|}\,:\tag{14}\] the scalar law of 2 applied to the magnitude of the gradient vector and acting along it, so the \(K\) gradient components share one congested conductance. At \(K=1\), 14 is 2 . The system \(Bf^k=\alpha d^k\) with 14 is the stationarity of the convex energy \(E(\phi^1,\dots,\phi^K)=\sum_e\Psi_e(\|\mathbf{g}_e\|)-\alpha\sum_k(d^k)^{\top}\phi^k\), strictly convex modulo the \(K\) per-commodity constants, so the equilibrium is unique. Physically the coupling turns on exactly where congestion does. On real capacity data at rush-hour load, commodity flows superpose to within solver tolerance (\(\lesssim10^{-4}\) relative)—real networks place each edge in one of BPR’s two power-like regimes where path splits are load-independent—while the shared load inflates marginal costs by up to \(8.8\times\): the congestion externality is slowdown, not rerouting. Synthetic instances with random capacities do reroute (\(12\)–\(19\%\) flow deviation), so the corpus sweep (§6.5) exercises the solver where commodities genuinely interact.

The Newton linearization of 14 is the \(K\)-commodity block Laplacian whose \((k,l)\) block is \(B\,\mathrm{diag}\big(c_e\delta_{kl}+(\rho'_e-c_e)\,u^k_eu^l_e\big)B^{\top}\): per edge, the \(K\times K\) block \[\label{eq:mcjac} \frac{\partial\mathbf{f}_e}{\partial\mathbf{g}_e} = c_e\,I_K+(\rho'_e-c_e)\,\mathbf{u}_e\mathbf{u}_e^{\top},\qquad c_e=\frac{\rho_e(s_e)}{s_e},\quad \mathbf{u}_e=\frac{\mathbf{g}_e}{s_e},\quad s_e=\|\mathbf{g}_e\| ,\tag{15}\] with eigenvalues \(c_e\) (multiplicity \(K-1\), across the flow direction) and \(\rho'_e\) (along it): symmetric positive definite for any monotone law, with anisotropy ratio \(c_e/\rho'_e=t_e'(F)F/t_e(F)\in[1,p]\) bounded by the BPR exponent.

6.2 Algorithm: one scalar hierarchy carries all \(K\) commodities↩︎

Algorithm 2 carries over with three substitutions and no new machinery. (i) Stacked state. The iterate is \(\Phi=(\phi^1,\dots,\phi^K)\), the residual \(R=(r^1,\dots,r^K)\), \(r^k=\alpha d^k-Bf^k\), and the stopping test is on the stacked norm. (ii) A scalar frozen operator. Rather than freeze the \(K\times K\) block linearization, NLF freezes a scalar hierarchy \(H\) on the single conductance \(w_e=\sqrt{c_e\rho'_e}\)—the geometric mean of the block’s two eigenvalues—and computes the correction as \(K\) independent inner solves, \(\delta^k\) from \(H\) against \(r^k\). Thus one \(\mathcal{O}(m)\) inner setup serves all \(K\) commodities and all chord steps. Because BPR cost depends only on total edge flow, the \(K\times K\) coupling per edge is rank-one, with eigenvalues \(c_e\) (multiplicity \(K-1\)) and \(\rho'_e\); the \(w_e\)-preconditioned per-edge block then has eigenvalues \(\{\sqrt{c_e/\rho'_e},\sqrt{\rho'_e/c_e}\}\)—symmetric about \(1\) with ratio \(c_e/\rho'_e\le p\) independent of \(K\), \(m\), and the iterate. The per-edge conditioning is thus uniformly controlled; the realized contraction of the assembled iteration is governed by the graph-wide anisotropy and is measured at run time by the refresh monitor below, not bounded a priori. (iii) A self-calibrating refresh monitor. The block-diagonal frozen operator has an intrinsic rate floor set by the anisotropy that no rebuild can mitigate, so a refresh monitor measures the contraction on the step following each rebuild and treats only degradation beyond that baseline as staleness; otherwise a heavily congested instance would rebuild every step to no effect.

The cost accounting is the point of the construction. A chord step costs one vector-law evaluation—\(m\) scalar inversions, the same count as a single-commodity step, since only \(\|\mathbf{g}_e\|\) is inverted—plus \(K\) cycles and \(K\) residual products: \(\mathcal{O}(Km)\), on top of one shared \(\mathcal{O}(m)\) setup. The entire \(K\)-dependence is \(K\) right-hand sides through one hierarchy.

6.3 Validation↩︎

The frozen block Jacobian 15 matches the finite-difference derivative of the residual; at \(K=1\) both the multigrid and the direct paths reproduce the single-commodity solver’s equilibrium, and both agree with the exact-Newton baseline (full \(Kn\times Kn\) block solve); identical commodities receive identical potentials and a negated demand a negated potential (neither symmetry is imposed); and the joint-vs-alone coupling check above. Table 9 shows the reference comparison on real TNTP road networks at \(K=4\) (four well-separated demand dipoles, benchmark load \(\alpha=0.3\,\mathrm{median}(\{c_e\}_e)\) per commodity): NLF reaches the exact-Newton equilibrium to \(10^{-10}\) with a handful of chord steps and a single hierarchy setup—whereas the exact-Newton competitor must assemble and factor the full \(Kn\times Kn\) block linearization at every step. This is the relevant high-accuracy baseline for the coupled multicommodity program (general convex solvers such as Ipopt reduce to the same the same \(Kn\times Kn\) block linear algebra); NLF matches its equilibrium while replacing the \(K\)-fold-larger block solve with one shared scalar \(\mathcal{O}(m)\) hierarchy.

6.4 Applications across domains↩︎

We exercise the \(K=4\) solver on the multicommodity applications the formulation targets: real TNTP metropolitan road networks (topology, capacities, free-flow times, BPR parameters, and the four heaviest origin–destination pairs all from the benchmark data, jointly scaled to a loaded rush hour, peak \(\|\mathbf{f}_e\|/c_e\approx1.07\)), and corpus-protocol instances from other coupled-flow domains where flows share one network—internet AS topology, a power transmission grid, national road networks. Every instance converges to \(10^{-9}\) with a handful of chord steps and \(1\)–\(3\) hierarchy setups, from \(38\) edges to \(1.5\times10^6\) (e.g.roadNet-PA, \(1.5\times10^6\) edges, \(26\) steps, \(2\) setups).

6.5 Corpus robustness and the cost of \(K\)↩︎

The decisive test is the same full-corpus sweep as §4.5, repeated at \(K=4\): the identical \(2{,}003\) real-world SuiteSparse graphs and instance protocol, with four well-separated demand dipoles per graph, each at the benchmark load. Every one of the \(2{,}003\) graphs converges to \(10^{-9}\) relative residual (Fig. 8), to \(m=1.8\times10^{7}\) edges, at a flat median of \(18\) chord steps (range \(9\)–\(51\)) and a median of \(1\) hierarchy setup (maximum \(5\); more than two on only \(120\) of \(2{,}003\), triggered by the refresh conditions of §3.1): the shared scalar hierarchy carries the coupled four-commodity system essentially as cheaply as it carries one. The fitted total cost is \(t\propto m^{1.04}\) over four decades, median \(9.1~\mu\)s per edge. Joined per graph against the single-commodity sweep, the \(K=4\) solve costs a median \(3.8\times\) the single-commodity solve (interquartile range \([3.0,5.6]\))—the \(\approx\!K\) per-step work with the setup amortized, as the cost accounting of §6.2 predicts.

The cost of \(K\) is measured directly by sweeping \(K\in\{1,2,4,8\}\) with nested demand sets on FEM and road graphs to \(1.5\times10^6\) edges (olesnik0, srb1, roadNet-PA): the chord-step count and the hierarchy-setup count are roughly constant, and wall-clock grows close to linearly in \(K\) (\(K{=}8\) at \(4.2\)–\(8.5\times\) the single-commodity cost), i.e., \(\mathcal{O}(Km)\) in total.

7 Extensions and Future work↩︎

7.1 Directed multicommodity assignment↩︎

§6 routes \(K\) coupled commodities under a smooth magnitude coupling with signed flows; real traffic assignment is directed (nonnegative arc flows, costs coupling through the directed link flow), validated against published origin–destination equilibria (Sioux Falls, Anaheim, …) [25], [34]. Two steps separate §6 from that target: the volume coupling \(\sum_k|f^k_e|\) forces a per-edge complementarity (only commodities whose potential drop matches the common congested cost may flow)—a free boundary of the same character as the max-flow cut, attacked by the continuation-with-safeguard machinery of §3.3; and the scale of \(K\) (thousands of demand columns), made conceivable by the \(\mathcal{O}(Km)\) per-step cost with one shared hierarchy (§6.2) plus per-origin aggregation borrowed from the bush-method literature [25], [26], making the \(K\) scaling conceivable. Benchmarking against Frank–Wolfe and Algorithm B on the published equilibria would be required.

7.1.0.1 The payoff: design optimization

The forward equilibrium is the means, not the end. Its motivating problems—congestion pricing, origin–destination matrix calibration, capacity network design—are bilevel programs optimized over a design \(x\) (tolls, demand, capacities) subject to the equilibrium \(f(x)\), where the load-bearing object is the sensitivity \(\partial f/\partial x\). Because \(f\) is the stationarity of 1 , that sensitivity is a single adjoint solve with the same symmetric Laplacian Jacobian: the gradient of an outer objective with respect to all \(k\) design variables costs one additional \(\mathcal{O}(m)\) solve, independent of \(k\). Equilibrium sensitivity is itself classical—the adjoint/implicit-function calculus of Tobin and Friesz [37] and the bilevel network-design literature it seeds compute exactly this derivative—so the contribution is not the adjoint but that here it is the same \(\mathcal{O}(m)\) Laplacian solve: available natively, and near-linear rather than a direct KKT factorization on poorly-separable designs. A directed single-commodity prototype we developed validates this on optimal tolling: the adjoint gradient matches finite differences. The undirected of this paper is thus a stepping stone: the same \(\mathcal{O}(m)\) Laplacian machinery that solves 1 also differentiates it, which is exactly what bilevel network design needs.

7.2 One-shot FAS↩︎

The methodologically pure alternative to Algorithm 2 is a Full Approximation Scheme (FAS) multilevel cycle (local linearization) [35], [38]. In initial experiments, the LAMG+ caliber-1 interpolation misrepresented the flux feeding the coarse \(\tau\) correction: caliber-1 (piecewise-constant) prolongation implicitly assumes the fine-level conductances \(\rho'_e\) are near-uniform within an aggregate, but near saturation \(\rho'\) varies by orders of magnitude across a single aggregate, so the Galerkin coarse operator \(P^{\top}JP\) no longer represents the fine-level flux balance and the nonlinear coarse correction is inconsistent. Building a variationally faithful transfer for the saturating law—e.g.a flux-based interpolation [35]—remains an open question.

7.3 Full multigrid: \(\alpha\) update at the coarsest level↩︎

Brandt’s nested-iteration ideal [35] would update the global scalar \(\alpha\) at the coarsest level rather than the finest, prolongating \((\phi,\alpha)\) upward with per-level arclength finishing. A first version was \(1.6\)–\(2.1\times\) faster on small instances but \(2.2\times\) slower on a larger grid (per-level finishing dominates)—a non-uniform gain—so we keep the simple \(\mathcal{O}(m)\) finest-level architecture and leave a tuned FMG variant to future work.

7.4 Optimal transport as a further instance↩︎

The same edge-separable energy expresses optimal transport in Beckmann (flow) form [1]: minimizing a monotone edge-separable cost subject to \(Bf=\mu-\nu\) is exactly 1 . The \(W_p\) Wasserstein distance [39] between measures \(\mu\) and \(\nu\) on a graph corresponds to the \(\ell^p\) edge cost; the quadratic case (\(p=2\), \(H^{-1}\) cost) reduces to a single \(\mathcal{O}(m)\) Laplacian solve, while the \(W_1\) limit (\(p\to1\)) gives the non-smooth saturating member, a natural fold instance. The quadratically-regularized form, whose dual Hessian is an active-subgraph Laplacian solved in the literature by direct factorization [40], is where a near-linear inner could help; whether it is competitive is an open question, since the \(W_1\)/Beckmann limit is a min-cost flow with strong combinatorial and direct solvers.

7.5 DC optimal power flow↩︎

The DC approximation to optimal power flow (DC-OPF) minimizes generation cost subject to power balance and thermal line limits: \[\label{eq:dcopf} \begin{align} &\min_{\theta,\,p^g}\;\sum_i C_i(p^g_i) \quad\text{s.t.}\quad B\,\mathrm{diag}(b)\,B^{\top}\theta = p^g - p^d,\\ &\qquad\quad |b_e(B^{\top}\theta)_e|\le c_e,\quad p^{\min}\le p^g\le p^{\max}, \end{align}\tag{16}\] where \(b_e\) is the line susceptance, \(p^d\) fixed demand, and \(C_i\) a convex (typically quadratic or piecewise-linear) generation cost. The power-balance constraint is a susceptance-weighted graph Laplacian; the line limits are the box constraints of the saturating law 6 , so 16 is an instance of 1 extended by a nodal cost and generator box constraints.

We merely record DC-OPF as an instance of the framework but do not solve the constrained problem here: the generator box constraints and nodal cost \(\sum_i C_i(p^g_i)\) require an FAS extension [38] that carries \(p^g\) and its bounds into the coarse levels (§7.2), which we have not built. As a back-of-envelope projection only—not a measured result—current IPM solvers (MATPOWER, PowerModels.jl) run \(20\)–\(50\) Newton steps, each a fresh sparse-direct KKT factorization; on near-planar grids (\(118\) to \(30{,}000\) buses) we timed that factorization at \(2\)–\(6\times\) cheaper than one NLF cycle, so whether a continuation-Newton loop reusing one frozen hierarchy would win on total cost turns on the iteration-count trade and is an open question.

8 Conclusion↩︎

Placing the nonlinearity in the edge law collapses a family of convex, edge-separable network-flow equilibria—electrical, congestion, minimum-delay, maximum-flow—onto one nonlinear Laplacian \(B\,\rho(B^{\top}\phi)=\alpha d\), solved by a damped chord-Newton on a frozen linearization cycled against the nonlinear residual. Any near-linear Laplacian solver is shown to empirically translate into the same complexity on the whole nonlinear class.

8.0.0.1 Reproducibility

Every table is produced by released scripts, on public benchmarks—the SuiteSparse collection (each instance reduced to its largest connected component, §4.5), the Transportation Networks road data [34], and MATPOWER grids [41]—and on seeded synthetic generators (Erdős–Rényi topologies for the scaling study) where noted. The interior-point comparator is Ipopt 3.14.19 (MUMPS 5.9.0, \(120\) s CPU budget); the BPR law uses the strictly convex regularization of §4.1; capacity/free-flow distributions and demand seeds are fixed and recorded. The inner solvers are public [10], [11]; the NLF solver, instance-generation scripts, and a reproducibility notebook that re-runs the timings and checks them against these tables are at https://github.com/orenlivne/nlf. All timings were collected on a MacBook Pro with Apple M5 Pro (18-core: 6 efficiency + 12 performance) and 48 GB unified memory, running a single Julia process with no explicit parallelism.

9 Why each application is an instance of the framework↩︎

Per application: equilibrium principle \(\Leftrightarrow\) convex program 4 \(\Leftrightarrow\) KKT \(\Leftrightarrow\) \(B\rho(B^{\top}\phi)=\alpha d\). Each instance differs only in the edge law; directed complementarity is absorbed into that law or relaxed by the signed formulation.

9.0.0.1 The chain in general

Stationarity of 4 in \(f_e\) gives the nonlinear Ohm law \(t_e(f_e)=(B^{\top}\phi)_e\); inversion gives \(f_e=\rho_e((B^{\top}\phi)_e)\); and Kirchhoff (\(Bf=\alpha d\)) yields 1 . Summing along any source–sink path telescopes: every path has marginal cost \(\phi_s-\phi_t\), the scalar equilibrium cost. The signed formulation has no inequality constraints—a smooth equation—which is why Newton with a near-linear Laplacian solve applies with no active sets.

9.0.0.2 Electrical flow (linear law)

Ohm’s law is linear, \(f_e=g_e/r_e\); the energy is the dissipated power and 4 is Thomson’s principle (current minimizes dissipation) [6]. Equation 1 is the ordinary weighted Laplacian system—the framework’s fixed point under linearization, and the object the inner solver [11] is built for.

9.0.0.3 Congestion (traffic) equilibrium

Wardrop’s principle—per origin–destination pair, path flows \(h_p\ge0\) with \(\sum_p h_p=d\) obeying \(T_p-\pi\ge0\) and \(h_p(T_p-\pi)=0\)—is, by Beckmann’s theorem [1], the KKT system of \(\min_h\sum_e\int_0^{f_e}t_e\): differentiating gives \(\partial F/\partial h_p=\sum_{e\in p}t_e=T_p\), and \(\pi\) (the multiplier of \(\sum_p h_p=d\)) is both the marginal cost of demand and the equilibrium travel time. In link variables the demand multiplier unfolds into the node field \(\phi\) (\(\pi=\phi_s-\phi_t\)) and the program becomes 4 with \(\Phi_e'=t_e\), the BPR latency of §4. The signed relaxation keeps the first Wardrop clause exactly (all routes at potential-drop cost \(\pi\)) and discards \(h_p\ge0\): an unused route may carry counterflow—the undirected scope, with per-commodity nonnegativity the free-boundary extension of §7.

9.0.0.4 Minimum-delay routing

M/M/1 queueing with Little’s law gives the per-packet link delay \(t_e(f)=1/(c_e-f)\), and Kleinrock’s independence approximation makes total delay edge-separable [4]. Taking the delay as marginal cost, \(\Phi_e=\log\!\big(c_e/(c_e-f)\big)\), gives 4 directly; the inverse law \(\rho_e(g)=c_e-1/g\) (odd-extended) saturates at capacity with conductance \(\rho'_e=(c_e-f)^2\to0\), placing the problem in the fold class of Table 1: the capacity constraint \(f<c_e\) is never imposed—the law cannot produce it—and the feasibility limit, where the filling edges form a cut in the spare capacities, is the fold handled by §3.3.

9.0.0.5 Maximum flow

The LP \(\max\,\alpha\) s.t.\(Bf=\alpha d\), \(|f_e|\le c_e\) has the box constraint as its only nonsmoothness; its complementarity (an edge is either slack with zero reduced cost or saturated on the cut) is absorbed into the smooth saturating law 6 , whose energy \(\Psi_e=c_e^2\log\cosh(g/c_e)\) is quadratic for slack edges and exactly linear, \(c_e|g|\), for saturated ones. The active set is not tracked: it emerges as the zero-conductance set \(w_e\to0\), the forming min cut, and the LP’s degeneracy reappears as the one scalar fold at \(\alpha=F^{*}\) traversed by the continuation of §3.3. At the fold the saturated set is the exact min cut and \(\phi\) steps across it (§5); no \((1+\epsilon)\) relaxation is involved.

References↩︎

---

1. Pine Birch Analytics, 35 Kelinger Rd, Churchville, PA 18966-1033 (oren.livne@gmail.com, tel. 312-533-9130, pinebirchanalytics.com; ORCID: 0000-0001-6700-483X).↩︎

2. Submitted to the editors 2026-07-05.↩︎

3. Notation: \(\rho\) acts componentwise and \(\rho_e\) denotes the law on an individual edge \(e\).↩︎

4. Justification of both asymptotics. Near the fold \(\phi=V\chi+O(1)\), \(\chi\) the source-side indicator, so every cut edge carries gradient \(g_e=V+O(1)\). The saturating tail of 6 is \(\rho_e(g)=c_e\tanh(g/c_e)=c_e-2c_e\,e^{-2g/c_e}+O(e^{-4g/c_e})\), and all \(s\)–\(t\) flow crosses the cut, so \(\alpha=\sum_{e\in\mathrm{cut}}\rho_e(g_e)=F^*-2\sum_{e\in\mathrm{cut}}c_e\, e^{-2V/c_e}\,(1+o(1))\). Hence \(F^*-\alpha\asymp e^{-2V/\bar c}\), \(\bar c\) the largest cut capacity, i.e.\(V=\tfrac{\bar c}{2}\log\tfrac{1}{F^*-\alpha}+O(1)\) and \(\psi=\chi^{\top}\phi\propto V\). The same expansion gives the cut conductances \(\rho'_e(g_e)=\mathrm{sech}^2(g_e/c_e)=4e^{-2V/c_e}(1+o(1))\), whence \(\lambda_{\min}(J)\le\sum_{\mathrm{cut}}\rho'_e\propto F^*-\alpha\) (here and below \(\lambda_{\min}\) is the smallest nonzero eigenvalue—\(J\) being a graph Laplacian, its zero eigenvalue with constant eigenvector \(\mathbf{1}\) is excluded, and \(\kappa(J)\) is correspondingly the ratio over the nonzero spectrum; the bound is the Rayleigh quotient at \(\chi\), which loads only the cut edges); for a simple fold the next eigenvalue stays \(\Theta(1)\) (the deflated spectral gap) and \(\lambda_{\max}(J)=\Theta(1)\), so the bound is two-sided and \(\kappa(J)\sim(F^*-\alpha)^{-1}\).↩︎

5. We use a strictly convex, twice-differentiable form of 11 with the same free-flow and degree-\(p\) congestion structure and real \((t^0_e,c_e,b,p)\) from the data; this regularizes the free-flow threshold so the dual problem is well posed. “Practical capacity” sentinels coded as uncapacitated links are given a large finite capacity.↩︎