Paper deep dive
Support Discovery With Iteratively Reweighted Least Squares for Fixed-Charge Network Flow
Sindura Saraswathi, Christian Kümmerle
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 9/10/2026, 2:37:22 AM
Summary
This paper proposes a scalable continuous-optimization algorithm for the Fixed-Charge Network Flow Problem (FCNFP) using an Iteratively Reweighted Least Squares (IRLS) framework. The method replaces the discontinuous fixed-charge objective with a smooth nonconvex Lasry-Lions surrogate, solving a sequence of weighted quadratic flow subproblems via a warm-started dual semismooth Newton method. To enhance support discovery, the algorithm incorporates objective-driven perturbation restarts and an anchor-union restricted search. Experiments on 410 instances demonstrate that this approach achieves superior objective quality compared to other scalable heuristics, with a mean gap of 1.316% to MILP references.
Entities (8)
Relation Signals (6)
Iteratively Reweighted Least Squares → uses → Lasry-Lions Surrogate
confidence 96% · The method replaces the discontinuous fixed-charge and linear arc cost objective with a smooth nonconvex Lasry--Lions surrogate
Fixed-Charge Network Flow Problem → solvedby → Iteratively Reweighted Least Squares
confidence 95% · We propose a scalable continuous-optimization algorithm for large-scale single-commodity FCNFP based on an iteratively reweighted least-squares (IRLS) framework.
Iteratively Reweighted Least Squares → solvessubproblemswith → Dual Semismooth Newton Method
confidence 94% · Each subproblem is solved by a warm-started dual semismooth Newton method
Dual Semismooth Newton Method → utilizesstructure → Graph Laplacian
confidence 93% · whose Newton systems have weighted graph-Laplacian structure, enabling the use of modern Laplacian solvers.
Fixed-Charge Network Flow Problem → comparedto → Mixed Integer Linear Programming
confidence 90% · Exact mixed-integer linear programming formulations capture the fixed-charge structure faithfully, but often become difficult to solve on large networks.
Fixed-Charge Network Flow Problem → relatedto → Dynamic Slope Scaling Procedure
confidence 85% · The dynamic slope scaling procedure (DSSP) [19] replaces the fixed-charge objective with a sequence of linear min-cost-flow (MCF) subproblems
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:The fixed-charge network flow problem (FCNFP) couples continuous flow allocation with discrete arc-activation decisions, making it a canonical but computationally challenging model for a variety of network design and resource allocation problems. Exact mixed-integer linear programming formulations capture the fixed-charge structure faithfully, but often become difficult to solve on large networks. We propose a scalable continuous-optimization algorithm for large-scale single-commodity FCNFP based on an iteratively reweighted least-squares (IRLS) framework. The method replaces the discontinuous fixed-charge and linear arc cost objective with a smooth nonconvex Lasry--Lions surrogate and solves a sequence of weighted quadratic flow subproblems. Each subproblem is solved by a warm-started dual semismooth Newton method whose Newton systems have weighted graph-Laplacian structure, enabling the use of modern Laplacian solvers. To further improve the discovered arc supports of the challenging underlying combinatorial problem, we also develop an algorithmic variant that incorporates objective-driven perturbation restarts and an anchor-union restricted search that jointly leverages supports discovered by IRLS and by complementary FCNFP heuristics. Computational experiments on 410 benchmark, synthetic, and large-scale instances show that our method obtains the best objective quality among the evaluated scalable FCNFP algorithms, with a mean gap of $1.316\%$ to a time-limited MILP reference and a win-or-tie rate of $90.0\%$ among the non-MILP methods. The results indicate that combining smooth continuous optimization with support-level search is an effective strategy for producing high-quality feasible solutions to large-scale FCNFP.
Tags
Links
- Source: https://arxiv.org/abs/2609.09295v1
- Canonical: https://arxiv.org/abs/2609.09295v1
Trouble viewing inline? Open PDF directly →
Full Text
118,798 characters extracted from source content.
Expand or collapse full text
Support Discovery With Iteratively Reweighted Least Squares for Fixed-Charge Network Flow Sindura Saraswathi Affiliation: University of Central Florida Christian Kümmerle Affiliation: University of Central Florida Affiliation: sindura@ucf.edu Affiliation: kuemmerle@ucf.edu Abstract The fixed-charge network flow problem (FCNFP) couples continuous flow allocation with discrete arc-activation decisions, making it a canonical but computationally challenging model for a variety of network design and resource allocation problems. Exact mixed-integer linear programming formulations capture the fixed-charge structure faithfully, but often become difficult to solve on large networks. We propose a scalable continuous-optimization algorithm for large-scale single-commodity FCNFP based on an iteratively reweighted least-squares (IRLS) framework. The method replaces the discontinuous fixed-charge and linear arc cost objective with a smooth nonconvex Lasry–Lions surrogate and solves a sequence of weighted quadratic flow subproblems. Each subproblem is solved by a warm-started dual semismooth Newton method whose Newton systems have weighted graph-Laplacian structure, enabling the use of modern Laplacian solvers. To further improve the discovered arc supports of the challenging underlying combinatorial problem, we also develop an algorithmic variant that incorporates objective-driven perturbation restarts and an anchor-union restricted search that jointly leverages supports discovered by IRLS and by complementary FCNFP heuristics. Computational experiments on 410410 benchmark, synthetic, and large-scale instances show that our method obtains the best objective quality among the evaluated scalable FCNFP algorithms, with a mean gap of 1.316%1.316\% to a time-limited MILP reference and a win-or-tie rate of 90.0%90.0\% among the non-MILP methods. The results indicate that combining smooth continuous optimization with support-level search is an effective strategy for producing high-quality feasible solutions to large-scale FCNFP. I Introduction Many problems in AI and operations research couple a discrete activation decision with a continuous allocation: which resources to switch on, and how much to route through them. Large-scale LLM inference is a recent instance, where deploying a serving configuration incurs a setup cost while routing requests through it incurs load-dependent operating cost [1]; the same fixed-plus-variable structure underlies classical network design and resource allocation. What makes such problems hard is not the continuous part, which is often polynomial once the active set is fixed, but the discrete structure that must be discovered. The fixed-charge network flow problem (FCNFP) provides a canonical model for this structure. Given a directed network, the goal is to route flow satisfying supply-demand balance and capacity constraints while minimizing both proportional flow costs and fixed charges incurred when arcs are activated. FCNFP arises in transportation [2, 3, 4], supply-chain design [5], and infrastructure planning [6]. FCNFP is NP-hard [7], with its computational difficulty arising from the discrete arc-activation structure. If the set of active arcs were known in advance, the remaining problem would reduce to a linear minimum cost network flow problem, which can be solved efficiently by linear programming methods or specialized solvers in O(nmlog(n2/m)log(nΓ))O(nm (n^2/m) (n )), where n and m denote the numbers of nodes and arcs and Γ bounds the integer arc costs [8]. Exact mixed-integer linear programming (MILP) [9, 10, 11, 12] formulations capture this structure faithfully, but they can become prohibitively expensive on large instances. In our experiments, a time-limited OR-Tools [13] MILP baseline solved all instances up to 6060 nodes and 400400 arcs, but began failing to prove optimality at 120120 nodes and 15001500 arcs and failed to prove optimality on all 10001000-node, 20,00020,000-arc instances. This motivates a principled continuous optimization approach that preserves the computational advantages of network-flow structure while inducing sparse arc support. This perspective comes from sparse linear inverse problems, where the combinatorial task of support selection is often handled through sparsity-promoting continuous penalties and iterative reweighting schemes [14, 15]. In this paper, we propose a continuous optimization strategy for large-scale single-commodity FCNFP. The method is based on an iteratively reweighted least-squares (IRLS) [15, 16] framework that replaces the discontinuous fixed-plus-linear arc cost with a nonconvex surrogate smoothed by Lasry–Lions regularization [17]. The smoothed surrogate [18] helps stabilize the nonconvex optimization landscape induced by fixed-charge penalties. At each iteration, the method solves a weighted continuous flow subproblem and updates the weights so that arcs carrying small flow are increasingly penalized. This reweighting mechanism encourages small-flow arcs to vanish and guides the continuous flow towards a sparse, high-quality feasible solution without explicitly optimizing over binary activation variables. Despite the Lasry–Lions smoothing strategy avoiding many suboptimal local minima, the underlying problem is still fundamentally nonconvex, and there are no guarantees available that the discovered support set is optimal. To find even better support sets of the network flow, we propose to augment the base IRLS method with two support-improvement layers. First, objective-driven perturbation restarts modify the incumbent support by suppressing expensive active arcs and opening promising inactive arcs, followed by projection back onto the feasible capacitated-flow polytope. Second, an anchor-union restricted search combines IRLS supports with supports discovered by external FCNFP heuristics, and then reruns IRLS on the resulting restricted subgraphs. The main contribution of this work is a scalable primal algorithm that bridges smooth continuous optimization and support-level combinatorial search. Unlike exact MILP approaches, the proposed method does not attempt to certify global optimality. Instead, it is designed to produce high-quality feasible solutions on instances where exact methods are too slow or terminate with only feasible, but suboptimal flow vectors. Computational experiments on 410410 benchmark, synthetic, and large-scale instances show that our method achieves the best objective quality among the scalable evaluated FCNFP heuristics. It achieves a mean gap of 1.316%1.316\% to the MILP reference and, among the non-MILP methods, wins or ties on 90.0%90.0\% of the instances. I Related Work Existing FCNFP Solvers. FCNFP is a representative problem in the broader class of nonconvex minimum-cost flow problems and is NP-hard. Exact MILP and branch-and-cut methods can model the FCNFP structure directly, but their scalability is limited by the binary arc-activation variables. This has motivated a large body of approximation and reformulation methods that exploit the underlying network-flow structure. The dynamic slope scaling procedure (DSSP) [19] replaces the fixed-charge objective with a sequence of linear min-cost-flow (MCF) subproblems by converting each arc’s fixed cost into a flow-dependent per-unit slope, producing an LP with network-flow constraints. These subproblems can be solved by generic LP solvers or by specialized min-cost-flow algorithms, with cost-scaling variants admitting the same O(nmlog(n2/m)log(nΓ))O(nm (n^2/m) (n )) bound, but the linearization only locally approximates the discontinuous term, making DSSP sensitive to initialization and prone to poor early support decisions. Related dynamic updating methods refine this idea. ADCUP [20] approximates FCNFP by a concave piecewise-linear network-flow problem and updates it dynamically, while [21] give a continuous bilinear formulation that represents the fixed-charge structure without binary variables, at the cost of nonconvexity; [22] build on this with a smoothing parameter and coupled linear programming subproblems. More recently, [23] proposed the adaptive dynamic slope scaling procedure (ADSSP) with edge deletion, emphasizing the role of support identification and arc elimination in large-scale FCNFP. Our method is related to these approaches in that it avoids full mixed-integer search and instead solves a sequence of continuous network-flow subproblems. However, our methodology fundamentally differs in several key aspects. First, rather than using slope scaling, which approximates the fixed-charge term with a sequence of linear surrogates, or an exact but nonconvex bilinear reformulation, we implicitly optimize a family of tailored, continuous surrogate objectives derived from a Lasry–Lions smoothing strategy [24, 25, 17] that converge to the original FCNFP objective in the limit. By adjusting the smoothing parameters of this surrogate family iteratively, our method increases the degree of non-convexity over time, allowing for the discovery of high-quality supports. Furthermore, the efficiently solvable subproblems we use are not linear or bilinear, but quadratic, allowing for an appropriate balance between solver efficiency and modeling accuracy. Smoothing Strategies and IRLS Framework. Smoothing has been developed successfully for a variety of nonsmooth optimization problems, especially nonsmooth convex composite optimization [26, 27, 28]. Seen primarily as an algorithmic acceleration device for such problems, smoothing methods can also serve, in nonconvex settings, as a landscape regularization by exposing useful descent directions and delaying premature commitment to poor local supports [18, 29, 30, 31]. Specifically, the Lasry–Lions smoothing strategy [24, 25, 17] that we adapt in this work is tailored to high-dimensional, nonconvex problems due to its flexibility involving two smoothing parameters, but does not directly lead to a computationally tractable optimization method [31]. We obtain a tractable algorithm by majorizing the Lasry–Lions surrogates by an adaptive quadratic model, which makes the proposed method an instance of the iteratively reweighted least-squares (IRLS) framework [32, 15, 16]. IRLS methods are well-known tools to solve hard instances of estimation problems with sparse combinatorial structure in computer vision [33], signal processing [15, 34], and machine learning [35], but have not been well-studied so far for network flow problems. Efficient Laplacian Solvers. The scalability of our method is also tied to recent progress on Laplacian linear-system solvers. Nearly-linear-time algorithms for graph Laplacian and symmetric diagonally dominant systems have become a central tool in theoretical computer science [36, 37, 38, 39, 40], with best available runtime bounds of O(m(loglogn)O(1)log(1/ϵ))O(m( n)^O(1) (1/ε)) for finding an ϵε-accurate solution [41]. While provably fast approximation algorithms have been developed for several network-flow problems [42, 37, 43], implementable efficient algorithms have been so far only developed for Laplacian system solvers [44, 45, 46], which we leverage in this work. I Problem Formulation Let G=(V,E)G=(V,E) be a directed network with node set V and edge set E. Let n=|V|n=|V| and m=|E|m=|E|. For each edge e∈Ee∈ E, let fef_e denote the flow on edge e, cec_e its capacity, rer_e its proportional cost, and beb_e its fixed cost. Let C∈ℝn×mC ^n× m denote the incidence matrix and let d∈ℝnd ^n denote the demand vector, where dv>0d_v>0 at sink nodes, dv<0d_v<0 at source nodes, dv=0d_v=0 at transshipment nodes, with ∑vdv=0 _vd_v=0 so that total supply equals total demand. With these definitions, the fixed-charge network flow problem (FCNFP) can be formulated as minf∑e∈E(refe+be|fe|>0)=:Jtrue(f)s.t.Cf=d,0≤fe≤ce,∀e∈E. split _f & _e∈ E (r_ef_e+b_e1\|f_e|>0\ )=:J_true(f)\\ s.t. &Cf=d,0≤ f_e≤ c_e, ∀ e∈ E. split (1) The key difficulty of FCNFP lies in the combinatorial nature of the fixed-charge cost be|fe|>0b_e1\|f_e|>0\, which is active if and only if the binary variable ze:=|fe|>0z_e:=1\|f_e|>0\ is non-zero, corresponding to the opening of the edge e for a non-zero flow. The set of feasible flow vectors of (1) is ℱ=f∈ℝm:Cf=d, 0≤f≤cF= \f ^m:Cf=d,\;0≤ f≤ c \. I-A Lasry–Lions Double Envelope Combined Surrogate The fixed-charge objective is discontinuous at zero, since an arbitrarily small positive flow activates the fixed charge, and it promotes sparsity, since once an arc carries flow, increasing that flow incurs no further charge. To obtain a tractable continuous formulation we replace the fixed-plus-linear arc cost with a smooth nonconvex surrogate based on the Lasry–Lions double envelope [17], which removes the discontinuity while preserving the support-promoting structure. The surrogate has three regimes and is governed by two smoothing parameters λ>μ>0λ>μ>0 shared across all edges. The parameter λ controls the overall smoothing level, while μ controls the width of the transition region. Near zero, the surrogate is quadratic, which keeps the IRLS weights finite and the subproblem well defined. The transition region connects the quadratic and active-flow regimes smoothly. For larger flows, the exact envelope’s linear tail has derivative rer_e, so beb_e would enter only as an additive constant and vanish from the reweighting rule. We therefore replace the linear tail by a logarithmic term with coefficient αbeα b_e, where α=λ/(λ+sλ0)α=λ/(λ+s _0), λ0 _0 is the initial smoothing level and s=10−2s=10^-2. For each edge e, the first transition point is x1,e=Te(λ−μ)/λx_1,e=T_e(λ-μ)/λ, where Te=λre+2λbeT_e=λ r_e+ 2λ b_e. The second transition point is x2,e=(Se+Se2−4μαbe)/2x_2,e=(S_e+ S_e^2-4μα b_e)/2, where Se=Te−μreS_e=T_e-μ r_e. We restrict 0<μ<λ0<μ<λ by the edge-wise cap μ≤mineμmax,e(λ)μ≤ _e _ ,e(λ), which makes these transitions real and ordered; the cap and ordering proof are given in supplementary material, Section A-A, Section A-B. The resulting per-edge Lasry–Lions surrogate is LLeλ,μ(x)=x22(λ−μ),|x|≤x1,e,Te22λ−(|x|−Te)22μ,x1,e<|x|<x2,e,re(|x|−x2,e)+αbelog(|x|x2,e)+ψ2,e(x2,e),|x|≥x2,e;L_e^λ,μ(x)= cases x^22(λ-μ), |x|≤ x_1,e,\\[10.0pt] T_e^22λ- (|x|-T_e)^22μ, x_1,e<|x|<x_2,e,\\[10.0pt] aligned &r_e(|x|-x_2,e)+α b_e ( |x|x_2,e )\\ & + _2,e(x_2,e), |x|≥ x_2,e; aligned cases (2) where ψ2,e(x2,e)=Te2/(2λ)−(x2,e−Te)2/(2μ) _2,e(x_2,e)=T_e^2/(2λ)-(x_2,e-T_e)^2/(2μ) makes the surrogate continuous at x2,ex_2,e. The transition points are chosen such that LLeλ,μL_e^λ,μ is continuously differentiable across the three regimes (see supplementary material, Section A-A for derivations). The full smoothed objective is JLLλ,μ(f)=∑e∈ELLeλ,μ(fe).J_L^λ,μ(f)= _e∈ EL_e^λ,μ(f_e). (3) Figure 1 illustrates this for an edge at three values of λ with μ=λ/2μ=λ/2. For large λ, the quadratic and transition regions are wider and the discontinuity is smoothed more strongly. As λ decreases, these regions contract toward the origin and the surrogate more closely tracks the fixed-plus-linear objective. Fig. 1: Lasry–Lions surrogate LLeλ,μL^λ,μ_e for an edge with re=0.02r_e=0.02, be=1b_e=1, at three smoothing levels. x1,ex_1,e and x2,ex_2,e represent the transition points. I-B Iteratively Reweighted Least Squares We minimize the smoothed objective in (3) through an iteratively reweighted least-squares (IRLS) scheme [16]. The idea is to replace the nonconvex surrogate at each iteration by a convex weighted quadratic model whose weights are computed from the current flow. At iteration k, given the current flow f(k)f^(k), IRLS constructs edge weights we(k)w_e^(k) and sets Wk=diag(w1(k),…,wm(k))W_k=diag (w_1^(k),…,w_m^(k) ). With the weighted inner product ⟨u,v⟩Wk:=u⊤Wkv u,v _W_k:=u W_kv, define the shifted quadratic model Qλ,μ(f∣f(k))=JLLλ,μ(f(k))+12⟨f,f⟩Wk−12⟨f(k),f(k)⟩Wk.Q^λ,μ(f f^(k))=J_L^λ,μ(f^(k))+ 12 f,f _W_k- 12 f^(k),f^(k) _W_k. (4) The terms independent of f do not affect the minimizer. Hence the IRLS update and its equivalent weighted quadratic formulation can be written compactly as f(k+1) f^(k+1) =argminf∈ℱQλ,μ(f∣f(k))=argminf∈ℱ12⟨f,f⟩Wk. = *argmin_f Q^λ,μ(f f^(k))= *argmin_f 12 f,f _W_k. (5) The IRLS weights are obtained from the local curvature of the Lasry–Lions surrogate. For a nonzero flow value, we define the weight by differentiating the per-edge surrogate in (2) and normalizing by the flow magnitude as we(x)=(LLeλ,μ)′(|x|)|x|w_e(x)= (L_e^λ,μ ) (|x|)|x|. This choice couples the quadratic subproblem to the local geometry of the Lasry–Lions surrogate. For parameters in the capped domain and exact analytical weights, Qλ,μ(⋅∣f(k))Q^λ,μ(\,· f^(k)) globally majorizes JLLλ,μJ_L^λ,μ and is tangent to it at the current iterate. Since small flows receive large weights, arcs carrying negligible flow become increasingly expensive in subsequent subproblems and are encouraged to vanish. Larger active flows receive smaller weights, allowing the algorithm to concentrate flow on a sparse support. The weights are given by we(fe)=1λ−μ,|fe|≤x1,e,Te−|fe|μ|fe|,x1,e<|fe|<x2,e,re|fe|+αbe|fe|2,|fe|≥x2,e.w_e(f_e)= cases 1λ-μ,&|f_e|≤ x_1,e,\\[10.0pt] T_e-|f_e|μ|f_e|,&x_1,e<|f_e|<x_2,e,\\[12.0pt] r_e|f_e|+ α b_e|f_e|^2,&|f_e|≥ x_2,e. cases (6) For fe=0f_e=0, the first case is used, so the weight remains finite. In the active-flow region, the term re/|fe|r_e/|f_e| represents the proportional cost contribution, while the term αbe/|fe|2α b_e/|f_e|^2 acts as an ℓ0 _0-like fixed-charge penalty. Under the capped parameter condition, Proposition 1 (Supplementary Material, Section A-B) shows that the exact weight and exact update in (5) constitute a majorization–minimization scheme [47] for the fixed-parameter surrogate JLLλ,μJ_L^λ,μ. At iteration k, the analytical weights are evaluated at the current flow f(k)f^(k), giving we(k)=we(fe(k)),e∈Ew_e^(k)=w_e(f_e^(k)),\ e∈ E. These weights are then fixed, and the next iterate is obtained by solving the weighted quadratic flow subproblem in (5), which can be solved efficiently (see Section I-C). Each subproblem preserves the original flow-conservation and capacity constraints. The numerical implementation clips weights and solves QPs inexactly, so it verifies fixed-parameter surrogate descent directly rather than relying on exact majorization alone. Overall, the weight update and (5) provide a continuous mechanism for seeking low-cost FCNFP flows without explicitly introducing binary activation variables or tree-based search. I-C Weighted Quadratic Program Solver Each IRLS iteration requires solving the weighted quadratic flow subproblem in (5), which we solve by a warm-started dual semismooth Newton (SSN) method [48]. Introducing a dual variable y∈ℝny ^n for the flow-conservation constraint Cf=dCf=d, the box-constrained primal minimizer for a fixed y is available in closed form as f(y)=clip(W−1C⊤y, 0,c)f(y)=clip\! (W^-1C y;\,0,\,c ), where the clipping is applied component-wise. The corresponding dual residual is R(y)=d−Cf(y)R(y)=d-Cf(y). Thus, solving the subproblem is equivalent to finding a root of R(y)=0R(y)=0. An SSN step solves (CD(y)C⊤)Δy=d−Cf(y), (CD(y)C ) y=d-Cf(y), (7) where D(y)=diag(we−10<(C⊤y)e<wece)D(y)=diag\! (w_e^-11\0<(C y)_e<w_ec_e\ ). Only free arcs, those clipped at neither bound, contribute to this Newton system, so CD(y)C⊤CD(y)C is a weighted graph Laplacian on the currently free subgraph, singular on each connected component’s constant potential; we remove this gauge freedom by mean-centering. Because D(y)D(y) vanishes on arcs at a bound, a literal Newton step from y=0y=0 is degenerate; each solver call is instead seeded by a predictor read off the current primal flow, with a regularized, line-searched Newton correction when the predictor is rejected (see supplementary material, Section A-C). Each Newton system is solved by preconditioned conjugate gradient (PCG) using a freshly computed approximate-Cholesky factorization [49] of the free-subgraph Laplacian as preconditioner, with an escalating iteration budget if the residual has not converged. The solver is warm-started across IRLS iterations from the previous dual potential y and the previous lower and upper active sets. Each Newton system costs O~(mlog(1/ϵ)) O(m (1/ε)) [44], so a run of K outer iterations averaging s¯ s Newton steps per subproblem performs O~(Ks¯mlog(1/ϵ)) O\! (K\, s\,m (1/ε) ) work. In our runs, s¯ s stays below 44 across a 200200-fold range in m (see supplementary material, Section A-C). The IRLS weights are clipped to [10−12,1012][10^-12,10^12], and each dual-SSN QP call is capped at 4040 Newton iterations. If the SSN does not converge and leaves a flow-conservation residual above 10−610^-6, the algorithm falls back to the primal-dual interior-point method (see supplementary material, Section A-C), which exploits the same weighted graph-Laplacian Newton systems. IV Algorithm We outline a core IRLS method for FCNFP in Algorithm 1: each iteration evaluates the weights (6) at the current flow and solves the weighted quadratic subproblem over ℱF with the dual semismooth Newton method of Section I-C. A trial iterate is accepted only if it does not increase JLLλ,μJ_L^λ,μ beyond a numerical tolerance; otherwise the algorithm backtracks along the segment joining the current and the trial flows; if no tested step is accepted, the current flow is retained. The reported solution is the best feasible candidate encountered during the run, scored under JtrueJ_true. The run is organized into stages. Within a stage, λ is updated non-increasingly based upon the flow gap quantity λgap _gap (see Section V). A new algorithmic stage is triggered if the relative flow change falls below 10−310^-3, if a line search rejects the step or after six weighted-QP solves, in which case λ is decreased more aggressively. Let |f|(1)≥⋯≥|f|(m)|f|_(1)≥·s≥|f|_(m) be the sorted flow magnitudes, let Kgap=argmax1≤j<m|f|(j)/|f|(j+1)K_gap= _1≤ j<m|f|_(j)/|f|_(j+1), and set f↑=|f|(Kgap)f_ =|f|_(K_gap), f↓=|f|(Kgap+1)f_ =|f|_(K_gap+1), and ρgap=f↑/f↓ _gap=f_ /f_ . The pair (f↑,f↓)(f_ ,f_ ) brackets the widest gap in the flow magnitudes, so KgapK_gap is the support size that gap implies, f↓f_ the largest magnitude outside that support, and ρgap _gap a measure of how cleanly the two sides separate. This signal is used both by a λ-update rule and the support-refit procedure below. In practice, relevant arcs and a support estimate could be identified via an activation threshold τ (e.g., of 10−610^-6), below which flows are set to zero; however, such a thresholding will yield generally infeasible flow vectors. To improve on this, we solve a linear min-cost flow restricted to the support set estimate S, which is updated whenever ρgap≥10 _gap≥ 10 and 0<Kgap≤max4n,10240<K_gap≤ \4n,1024\ by setting S=e:|fe|>f↑f↓S=\e:|f_e|> f_ f_ \. While this is not an exact solver for the FCNFP problem on the S-restricted subgraph defined by S, this support refit procedure Refit provides a mechanism to meaningfully score candidate supports. Each stage boundary additionally produces a full-graph candidate, which may open arcs outside S. Candidates enter f⋆f only when feasible and strictly improving. Algorithm 1 terminates if ρgap≥100 _gap≥ 100, KgapK_gap has been stable for four consecutive iterations, and f⋆f has not improved for eight iterations. Alternatively, termination is triggered if the objective value JLLλ,μ(f(k))J_L^λ,μ(f^(k)) stagnates or if λ reaches a numerical precision floor. More details on support refit, step acceptance, stagnation, and candidate evaluation are given in the supplementary material, Section A-D. Algorithm 1 IRLS for FCNFP 1: Network (C,d,c,r,b)(C,d,c,r,b), f(0)∈ℱf^(0) , K=200K=200 2: (λ0,μ0)←InitializeSmoothing(f(0))( _0, _0)← InitializeSmoothing(f^(0)); f⋆←f(0)f ← f^(0) 3: for k=0,…,K−1k=0,…,K-1 do 4: w(k)←Weights(f(k),λk,μk)w^(k)← Weights(f^(k), _k, _k), ⊳ Eq. (6) 5: (f~,y)←DualSSN(w(k),ℱ,y)( f,y)← DualSSN(w^(k),F;y) ⊳ Eq. (5) 6: (f(k+1),y)←Accept(f(k),f~,λk,μk)(f^(k+1),y)← Accept(f^(k), f, _k, _k) ⊳ backtrack 7: g(k+1)=(Kgap,ρgap,f↑,f↓)←FlowGap(f(k+1))g^(k+1)=(K_gap, _gap,f_ ,f_ )← FlowGap(f^(k+1)) 8: f⋆←Best(f⋆,f(k+1),Refit(f(k+1),g(k+1)))f ← Best (f ,f^(k+1), Refit(f^(k+1),g^(k+1)) ) 9: (λk+1,μk+1)←Continuation(λk,μk,g(k+1))( _k+1, _k+1)← Continuation( _k, _k,g^(k+1)) 10: if Stop then break 11: end if⊳λ≤λmin \ λ≤ _ , stable gap, stagnation 12: end for 13: return f⋆f V Smoothing Continuation Scheme For the analytical surrogate with unmodified costs, JLLλ,μ(λ)J_L^λ,μ(λ) converges pointwise to JtrueJ_true on ℱF as λ↓0λ 0 along any admissible path μ(λ)μ(λ); in particular, this holds for μ(λ)=λ/2μ(λ)=λ/2 (Supplementary Material, Section A-A). This consistency motivates starting from a broadly smoothed model and progressively sharpening it toward the fixed-charge objective [50]. Thus, we decrease λ by a stage-wise continuation rule and set μ=minλ/2,mine∈Eμmax,e(λ)μ= \λ/2,\ _e∈ E _ ,e(λ)\ at every parameter update. We propose a specific continuation schedule that balances rapid sharpening against premature support commitment. Specifically, since Te=λre+2λbeT_e=λ r_e+ 2λ b_e has the scale of a flow magnitude, λ has units of squared flow per unit cost, which motivates fixing λ0 _0 at the scale of an initial feasible flow f(0)f^(0). With fmedf_med the median initial flow magnitude above one and bmedb_med the median positive fixed cost, sparse flow initializations use λ0=10−2maxfmed2/(2bmed),1 _0=10^-2 \f_med^2/(2b_med),1\, which places TeT_e ten times below fmedf_med so that the active arcs of the start begin in the logarithmic tail. Dense flow starts are instead over-smoothed with λ0=102maxmaxe|fe(0)|,1 _0=10^2 \ _e|f^(0)_e|,1\. We set μ0=minθ0λ0,mine∈Eμmax,e(λ0) _0= \ _0 _0,\ _e∈ E _ ,e( _0)\ with θ0=0.3 _0=0.3 and 0.70.7, respectively. Both initial parameter pairs are admissible under Proposition 1 (Section A-B) by construction. The update candidate for λk+1 _k+1 is the gap-based value λgap=f↓ 2/(2brep) _gap=f_ ^\,2/(2b_rep), where brepb_rep is a representative fixed cost of the arcs near the gap defined by KgapK_gap and is motivated by the scale match of Te≈2λgapbrep=f↓T_e≈ 2 _gapb_rep=f_ for negligible rer_e. Within a stage, we set λk+1=maxminλk,λgap,λk/5 _k+1= \ \ _k, _gap\, _k/5\, whereas we set λk+1=maxminλk/4,λgap,λk/100,10−9 _k+1= \ \ _k/4, _gap\, _k/100,10^-9\ across stage boundaries. We use this acceleration only after the detected gap is sufficiently pronounced and stable; the trust conditions and remaining schedule details are given in Section A-E. For exact weights and an exact QP solve, Proposition 1 (Section A-B) gives JLLλk,μk(f(k+1))≤JLLλk,μk(f(k))J_L _k, _k(f^(k+1))≤ J_L _k, _k(f^(k)); acceptance/backtracking, with rejection as a fallback, enforces the same current-pair nonincrease numerically. Continuation changes the surrogate, so no monotonic comparison across parameter pairs is claimed. Thus, fixed-pair descent supplies stability while continuation supplies sharpening. The best feasible incumbent remains nonincreasing under JtrueJ_true. VI Algorithmic Layers The IRLS procedure provides a continuous mechanism for refining a feasible flow under the smoothed fixed-charge surrogate. However, because the surrogate is nonconvex, the final support can strongly depend on the starting point. We therefore use a layered design that combines continuous reweighting with structured support exploration. Each layer maintains the best feasible solution, scored by the true fixed-charge objective JtrueJ_true. Layer 1: Multi-start IRLS. The first layer runs the IRLS core method of Algorithm 1 from five complementary feasible initializations: a linear-cost-weighted quadratic-flow solution, a quadratic-flow solution weighted by the combined linear and capacity-amortized fixed cost, a linearized MCF solution, a capacity-amortized MCF solution, and a DSSP slope-scaling iterate; implementation details are given in supplementary material, Section A-F. Duplicate supports are removed, and IRLS is run from each remaining start. This layer gives the algorithm access to several plausible support regions without introducing binary activation variables. Layer 2: Objective-driven perturbation restarts. The second layer explores alternatives around each support discovered by Layer 1. Given an incumbent flow f, active arcs are scored by their realized average cost re+be/maxfe,1r_e+b_e/ \f_e,1\, while inactive arcs are scored by their capacity-amortized opening cost re+be/maxce,1r_e+b_e/ \c_e,1\. These scores identify expensive active arcs that may be worth suppressing and cheap inactive arcs worth opening. Using these scores, we form a perturbed target f~ f by cutting flow on the expensive active arcs and putting flow on the cheap inactive ones. The resulting perturbed targets are not required to be feasible; instead, each target is projected back onto the capacitated flow polytope ℱ=g∈ℝ|E|:Cg=d, 0≤g≤cF=\g ^|E|:Cg=d,\;0≤ g≤ c\ via Πℱ(f~)=argming∈ℱ‖g−f~‖22 _F( f)= _g \|g- f\|_2^2. Projections that fail or duplicate a previously generated support are discarded. The surviving projections are ranked using their true objective values, support changes, and numbers of active arcs, and the highest-ranked candidates are used as new IRLS starts. This layer acts as a structured escape mechanism from locally stable but suboptimal supports. Details of the perturbation rules, candidate selection, and projection step are given in supplementary material, Section A-F. Layer 3: Anchor-union restricted IRLS. The third layer expands the search beyond the supports discovered by IRLS and its perturbation restarts. This is important because some useful arcs may never enter an IRLS trajectory once the surrogate becomes sharp, especially when opening them requires coordinated changes across an entire path. We therefore augment the IRLS supports with supports obtained from external FCNFP algorithms. Specifically, we use ADSSP solutions generated from the max_flow, unitflow, trans, and multi_unit initializations as support anchors. These solutions often provide sparse feasible supports containing useful combinatorial structure, even when they do not attain the best objective value. Let IRLSX_IRLS denote the candidate flows from Layer 2, and let ancX_anc denote anchor flows from the external solver. For any flow x, define its support S(x)=e∈E:xe>τS(x)=\e∈ E:x_e>τ\, where τ is the arc activation threshold. Each candidate union support has the form E′=S(f)∪S(a)E =S(f)∪ S(a), where f∈IRLSf _IRLS and S(a)S(a) is empty, a single anchor support, or the union of all anchor supports. Duplicate, empty, and infeasible supports are discarded. For each retained union E′E , we compute feasible refit flows by solving restricted linear flow problems with costs consisting of variable cost and capacity-amortized fixed cost. Several fixed-cost weights are tried, together with the corresponding Layer 2 flow, and the best feasible starts are retained. We then run IRLS on the subproblem obtained by fixing all arcs outside E′E to zero. A lightweight within-union refinement then tests deleting expensive active arcs and reintroducing promising inactive ones, keeping only feasible improvements. All candidates are scored under the original fixed-charge objective JtrueJ_true, and the best feasible solution is retained. Further implementation details are given in supplementary material, Section A-F. Algorithm 2 gives this complete anchor-union restricted IRLS procedure. Algorithm 2 Layer 3: Anchor-Union Restricted IRLS 1: Network data (C,d,c,r,b)(C,d,c,r,b), IRLS starts ℐI, ADSSP starts A 2: ←IRLS(f(0)):f(0)∈ℐX←\ IRLS(f^(0)):f^(0) \ ⊳ Layer 1 3: ←Πℱ(f~):f~∈Perturb()P←\ _F( f): f∈ Perturb(X)\ 4: IRLS←∪IRLS(p):p∈,p feasibleX_IRLS ∪\ IRLS(p):p ,\ p feasible\ 5: x⋆←argminf∈IRLSJtrue(f)x ← _f _IRLSJ_true(f); J⋆←Jtrue(x⋆)J ← J_true(x ) ⊳ Layer 2 6: anc←ADSSP(i):i∈X_anc←\ ADSSP(i):i \ 7: +←∅∪S(a):a∈anc∪⋃a∈ancS(a)A^+←\ \∪\S(a):a _anc\∪\ _a _ancS(a)\ 8: ←S(f)∪A:f∈IRLS,A∈+U←\\,S(f)∪ A:f _IRLS,\ A ^+\,\ 9: for E′∈E do 10: (E′)←f|E′:f∈IRLS∪Refit(E′)C(E )←\f|_E :f _IRLS\∪\ Refit(E )\ 11: x¯←IRLS(E′,(E′)) x← IRLS(E ;\,C(E )) 12: x^←Refine(x¯,E′) x← Refine( x,E ) ⊳ local add/drop within E′E 13: if Jtrue(x^)<J⋆J_true( x)<J then (x⋆,J⋆)←(x^,Jtrue(x^))(x ,J )←( x,\,J_true( x)) 14: end if 15: end for 16: return x⋆x VII Experiments We evaluate the proposed method on a collection of 410410 single-commodity FCNFP instances drawn from benchmark, synthetic, and large-scale families. Benchmark sets support comparison to prior work; scale families probe MILP-intractable sizes. The set includes FCNetLib_fc from [51], benchmark from [52], synthetically generated synthetic instances, and nine scale families from [23]. All methods are evaluated using the true fixed-charge objective in (1), with arc activation threshold τ=10−6τ=10^-6. A returned flow is treated as feasible only if the flow-conservation residual, capacity violation, and nonnegativity violation are each at most 10−610^-6. We compare the proposed method against three FCNFP baselines and a MILP reference. The baselines are linearized MCF, DSSP, and ADSSP. Linearized MCF solves a min-cost-flow problem using arc costs that approximate the combined fixed and proportional costs. DSSP [19] and ADSSP [23] are described in Section I. When the MILP proves optimality, its objective provides a certified optimum; otherwise, its feasible incumbent is used only as a time-limited reference. In our experiments, L1, L2 and L3 denote Layers 1–3 as described in Section VI. All heuristic methods use multiple feasible starts when applicable, and the best feasible solution under JtrueJ_true is reported. Each method is given a solver time limit of 300300 seconds per instance. Reported times are wall-clock and include model construction and solution extraction, so they can exceed this limit. Each individual IRLS run is limited to 200200 outer reweighting iterations. We report two objective-gap metrics. For a method z on instance i with objective Jz(i)J_z(i), the gap relative to MILP is GapMILP(z,i)=100⋅Jz(i)−JMILP(i)|JMILP(i)|,Gap_MILP(z,i)=100· J_z(i)-J_MILP(i)|J_MILP(i)|, where JMILP(i)J_MILP(i) is either a certified optimum or a time-limited MILP incumbent. The gap relative to the best non-MILP method is Gapbest(z,i)=100⋅Jz(i)−Jbest(i)|Jbest(i)|,Gap_best(z,i)=100· J_z(i)-J_best(i)|J_best(i)|, where Jbest(i)J_best(i) is the best feasible objective among the non-MILP methods on instance i. Full hardware, software versions, random seeds, and instance-generation details are given in supplementary material, Section A-G. Results. Table I reports the aggregate comparison over all 410410 instances with gaps taken relative to the MILP objective—a certified optimum where the MILP proved optimality, and a time-limited incumbent otherwise. MILP proved optimality on 235235 of the 410410 instances. L3 achieves the best objective quality among all heuristics, with a mean gap of 1.316%1.316\% and a median gap of 0.131%0.131\%. This improves over L2, L1, ADSSP, DSSP, and MCF. Here, L3- denotes L3 with anchors disabled; it still outperforms ADSSP, showing that L3’s gain does not depend on its anchor supports alone (see Supplementary Material Section A-H). TABLE I: Aggregate performance on the 410410 instances with a matched feasible MILP baseline. Gaps are relative to the MILP objective; times are seconds. Gap to MILP baseline (%) Time (s) Method mean p25 median p75 mean median MILP 0.000 0.000 0.000 0.000 144.939 64.516 L3 1.316 0.000 0.131 1.278 17.943 8.074 L3- 2.073 0.013 0.560 2.060 6.523 1.775 L2 2.334 0.037 0.749 2.424 5.999 2.632 L1 3.420 0.077 1.412 3.893 1.707 0.723 ADSSP 2.931 0.027 0.695 3.792 0.370 0.151 DSSP 4.084 0.045 2.164 5.903 1.391 0.369 MCF 8.304 0.118 4.199 11.577 0.089 0.043 Table I reports the non-MILP methods only comparison on all 410410 instances. This table isolates the contribution of the algorithmic layers without relying on a MILP reference. The base IRLS layer has a mean gap of 2.163%2.163\% to the best displayed heuristic. L2 reduces this gap to 1.106%1.106\%, while L3 reduces it further to 0.140%0.140\%. L3 is also the best or tied-best method on 90.0%90.0\% of instances. TABLE I: Heuristic layer-ablation performance over all 410410 instances. Gaps are relative to the best displayed heuristic on each instance. Method Mean gap (%) Median gap (%) Win/tie rate (%) L3 0.140 0.000 90.0 L2 1.106 0.441 22.9 L1 2.163 1.100 12.0 ADSSP 1.709 0.402 16.6 DSSP 2.848 1.468 12.2 MCF 6.903 3.765 6.1 Family-wise and MILP-status breakdowns further show that L3’s advantage persists across instance scales and MILP certification status (see Supplementary Material Section A-H). Discussion. The experiments indicate that support quality is the primary determinant of performance in large-scale FCNFP. IRLS (Algorithm 1) provides a suitable continuous mechanism for identifying flows with locally minimal fixed-charge cost objectives, with its final support still related to the initial flow vector due to the inherent combinatorial problem structure. Perturbation restarts (introduced via L2) empirically improve final flow costs by exploring nearby support with the anchor-union layer of the L3 pipeline of Algorithm 2 leading to further performance increase to achieve the state-of-the-art among non-MILP methods. The intuition is that a strictly better support can be found by jointly considering IRLS-discovered supports with support candidates identified from other methodologies, before reoptimizing via IRLS continuously over the resulting restricted subgraphs. The ADSSP anchor supports, which are part of the L3 pipeline, help to improve performance over L1 or L2, with L3, however, improving performance beyond a pure ADSSP method through its IRLS-enabled support discovery. The runtime results show the expected quality–time tradeoff. MCF and ADSSP are attractive when solutions are needed almost immediately. L3 is more expensive because it performs multiple IRLS solves, restricted refits, and support refinements. However, it remains far cheaper than MILP on average while consistently providing the strongest heuristic objective values. Therefore, the method is best viewed as a high-quality primal algorithm for settings where a strong feasible solution is more important than a certificate of global optimality. VIII Conclusion We proposed a scalable continuous-optimization algorithm based on IRLS for the fixed-charge network flow problem. The method replaces the discontinuous fixed-plus-linear cost objective with a smooth Lasry–Lions surrogate, solves the resulting weighted quadratic flow subproblems using Laplacian-structured linear algebra, and augments the continuous reweighting process with perturbation restarts and anchor-union restricted search. Across 410410 instances, our method achieves the best objective quality among the evaluated FCNFP baselines, including on large-scale instances where the MILP cannot certify optimality. These results show that combining smooth continuous optimization with structured support exploration is an effective strategy for large-scale FCNFP. Acknowledgments This work was supported in part by the NSF grant CCF-2549926. References [1] J. Cheng and D. T. Nguyen (2026) Fast heterogeneous serving: scalable mixed-scale llm allocation for slo-constrained inference. arXiv preprint arXiv:2604.07472. Cited by: §I. [2] S. Sadeghi-Moghaddam, M. Hajiaghaei-Keshteli, and M. Mahmoodjanloo (2019) New approaches in metaheuristics to solve the fixed charge transportation problem in a fuzzy environment. Neural computing and applications 31 (Suppl 1), p. 477–497. Cited by: §I. [3] M. Hajiaghaei-Keshteli, M. Aminnayeri, and S. F. Ghomi (2014) Integrated scheduling of production and rail transportation. Computers & Industrial Engineering 74, p. 240–256. Cited by: §I. [4] S. Golmohamadi, R. Tavakkoli-Moghaddam, and M. Hajiaghaei-Keshteli (2017) Solving a fuzzy fixed charge solid transportation problem using batch transferring by new approaches in meta-heuristic. Electronic Notes in Discrete Mathematics 58, p. 143–150. Cited by: §I. [5] J. Willet, K. Wetser, J. E. Dykstra, A. B. Bianchi, G. H. Oude Essink, and H. H. Rijnaarts (2021) WaterROUTE: a model for cost optimization of industrial water supply networks when using water resources with varying salinity. Water research 202, p. 117390. Cited by: §I. [6] S. Bogs, A. Abdelshafy, and G. Walther (2025) Planning minimum regret co2 pipeline networks. IISE Transactions, p. 1–14. Cited by: §I. [7] G. M. Guisewite and P. M. Pardalos (1990) Minimum concave-cost network flow problems: applications, complexity, and algorithms. Annals of Operations Research 25 (1), p. 75–99. Cited by: §I. [8] A. V. Goldberg and R. E. Tarjan (1990) Finding minimum-cost circulations by successive approximation. Mathematics of Operations Research 15 (3), p. 430–466. Cited by: §I. [9] K. G. Murty (1968) Solving the fixed charge problem by ranking the extreme points. Operations research 16 (2), p. 268–279. Cited by: §I. [10] B. Gendron and M. Larose (2014) Branch-and-price-and-cut for large-scale multicommodity capacitated fixed-charge network design. EURO Journal on Computational Optimization 2 (1-2), p. 55–75. Cited by: §I. [11] F. Ortega and L. A. Wolsey (2003) A branch-and-cut algorithm for the single-commodity, uncapacitated, fixed-charge network flow problem. Networks: An International Journal 41 (3), p. 143–158. Cited by: §I. [12] U. S. Palekar, M. H. Karwan, and S. Zionts (1990) A branch-and-bound method for the fixed charge transportation problem. Management science 36 (9), p. 1092–1105. Cited by: §I. [13] L. Perron and V. Furnon (2025) Google OR-Tools, version 9.12. Note: https://developers.google.com/optimization/Released February 17, 2025 Cited by: §I. [14] I. Daubechies, M. Defrise, and C. De Mol (2004) An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences 57 (11), p. 1413–1457. Cited by: §I. [15] I. Daubechies, R. DeVore, M. Fornasier, and C. S. Güntürk (2010) Iteratively reweighted least squares minimization for sparse recovery. Comm. Pure Appl. Math. 63, p. 1–38. Cited by: §I, §I, §I. [16] C. Kümmerle, C. Mayrink Verdun, and D. Stöger (2021) Iteratively reweighted least squares for basis pursuit with global linear convergence rate. Advances in Neural Information Processing Systems 34, p. 2873–2886. Cited by: §I, §I, §I-B. [17] M. Simões, A. Themelis, and P. Patrinos (2021) Lasry-lions envelopes and nonconvex optimization: a homotopy approach. In 2021 29th European Signal Processing Conference (EUSIPCO), p. 2089–2093. Cited by: §I, §I, §I, §I-A. [18] X. Chen (2012) Smoothing methods for nonsmooth, nonconvex minimization. Mathematical Programming 134 (1), p. 71–99. Cited by: §I, §I. [19] D. Kim and P. M. Pardalos (1999) A solution approach to the fixed charge network flow problem using a dynamic slope scaling procedure. Operations Research Letters 24 (4), p. 195–203. Cited by: §I, §VII. [20] A. Nahapetyan and P. Pardalos (2008) Adaptive dynamic cost updating procedure for solving fixed charge network flow problems. Computational Optimization and Applications 39 (1), p. 37–50. Cited by: §I. [21] S. Rebennack, A. Nahapetyan, and P. M. Pardalos (2009) Bilinear modeling solution approach for fixed charge network flow problems. Optimization Letters 3 (3), p. 347–355. Cited by: §I. [22] Z. Nie and S. Wang (2020) A dynamic method to solve the fixed charge network flow problem. IFAC-PapersOnLine 53 (2), p. 11231–11236. Cited by: §I. [23] L. Yang and Z. Yang (2024) A sequential reduction algorithm for the large-scale fixed-charge network flow problems. Optimization Letters 18 (5), p. 1131–1149. Cited by: §A-G, §I, §VII, §VII. [24] J. Lasry and P. Lions (1986) A remark on regularization in hilbert spaces. Israel Journal of Mathematics 55 (3), p. 257–266. Cited by: §I, §I. [25] H. Attouch and D. Aze (1993) Approximation and regularization of arbitrary functions in hilbert spaces by the lasry-lions method. In Annales de l’Institut Henri Poincaré C, Analyse non linéaire, Vol. 10, p. 289–312. Cited by: §I, §I. [26] A. Ben-Tal and M. Teboulle (2006) A smoothing technique for nondifferentiable optimization problems. In Optimization: Proceedings of the Fifth French-German Conference Held in Castel-Novel (Varetz), France, Oct. 3–8, 1988, p. 1–11. Cited by: §I. [27] Y. Nesterov (2005) Smooth minimization of non-smooth functions. Mathematical programming 103 (1), p. 127–152. Cited by: §I. [28] A. Beck and M. Teboulle (2012) Smoothing and first order methods: a unified framework. SIAM Journal on Optimization 22 (2), p. 557–580. Cited by: §I. [29] W. Bian and X. Chen (2013) Worst-case complexity of smoothing quadratic regularization methods for non-lipschitzian optimization. SIAM Journal on Optimization 23 (3), p. 1718–1741. Cited by: §I. [30] H. Yang, P. Antonante, V. Tzoumas, and L. Carlone (2020) Graduated non-convexity for robust spatial perception: from non-minimal solvers to global outlier rejection. IEEE Robotics and Automation Letters 5 (2), p. 1127–1134. Cited by: §I. [31] G. Kornowski and O. Shamir (2021) Oracle complexity in nonsmooth nonconvex optimization. Advances in Neural Information Processing Systems 34, p. 324–334. Cited by: §I. [32] I. F. Gorodnitsky and B. D. Rao (1997) Sparse signal reconstruction from limited data using focuss: a re-weighted minimum norm algorithm. IEEE Trans. Signal Process. 45 (3), p. 600–616. Cited by: §I. [33] P. Ochs, A. Dosovitskiy, T. Brox, and T. Pock (2015) On iteratively reweighted algorithms for nonsmooth nonconvex optimization in computer vision. SIAM J. Imaging Sci. 8 (1), p. 331–372. Cited by: §I. [34] D. Wipf and S. Nagarajan (2010) Iterative reweighted ℓ1 _1 and ℓ2 _2 methods for finding sparse solutions. IEEE Journal of Selected Topics in Signal Processing 4 (2), p. 317–329. Cited by: §I. [35] C. Kümmerle and C. M. Verdun (2021) A scalable second order method for ill-conditioned matrix completion from few samples. In International Conference on Machine Learning, p. 5872–5883. Cited by: §I. [36] D. A. Spielman and S. Teng (2004) Nearly-linear time algorithms for graph partitioning, graph sparsification, and solving linear systems. In Proceedings of the thirty-sixth annual ACM symposium on Theory of computing, p. 81–90. Cited by: §I. [37] S. I. Daitch and D. A. Spielman (2008) Faster approximate lossy generalized flow via interior point algorithms. In Proceedings of the Fortieth Annual ACM Symposium on Theory of Computing, p. 451–460. Cited by: §I. [38] S. Teng (2010) The laplacian paradigm: emerging algorithms for massive graphs. In Theory and Applications of Models of Computation: 7th Annual Conference, TAMC 2010, Prague, Czech Republic, June 7-11, 2010. Proceedings 7, p. 2–14. Cited by: §I. [39] D. A. Spielman and S. Teng (2014) Nearly linear time algorithms for preconditioning and solving symmetric, diagonally dominant linear systems. SIAM Journal on Matrix Analysis and Applications 35 (3), p. 835–885. Cited by: §I. [40] A. Madry (2018) Gradients and flows: continuous optimization approaches to the maximum flow problem. In Proceedings of the International Congress of Mathematicians: Rio de Janeiro 2018, p. 3361–3387. Cited by: §I. [41] A. Jambulapati and A. Sidford (2025) Ultrasparse ultrasparsifiers and faster laplacian system solvers. ACM Transactions on Algorithms 21 (3), p. 1–49. Cited by: §I. [42] P. Christiano, J. A. Kelner, A. Madry, D. A. Spielman, and S. Teng (2011) Electrical flows, laplacian systems, and faster approximation of maximum flow in undirected graphs. In Proceedings of the forty-third annual ACM symposium on Theory of computing, p. 273–282. Cited by: §I. [43] G. Karakostas (2008) Faster approximation schemes for fractional multicommodity flow problems. ACM Transactions on Algorithms (TALG) 4 (1), p. 1–17. Cited by: §I. [44] R. Kyng and S. Sachdeva (2016) Approximate gaussian elimination for laplacians-fast, sparse, and simple. In 2016 IEEE 57th Annual Symposium on Foundations of Computer Science (FOCS), p. 573–582. Cited by: §I, §I-C. [45] Y. Gao, R. Kyng, and D. A. Spielman (2023) Robust and practical solution of laplacian equations by approximate elimination. arXiv preprint arXiv:2303.00709. Cited by: §I. [46] Y. Gao, R. Kyng, and D. A. Spielman (2026) AC(k): robust solution of laplacian equations by randomized approximate cholesky factorization. SIAM Journal on Scientific Computing 48 (3), p. A1643–A1684. Cited by: §I. [47] Y. Sun, P. Babu, and D. P. Palomar (2017) Majorization-Minimization algorithms in signal processing, communications, and machine learning. IEEE Transactions on Signal Processing 65 (3), p. 794–816. External Links: Document Cited by: §I-B. [48] M. Hintermüller, K. Ito, and K. Kunisch (2002) The primal-dual active set strategy as a semismooth Newton method. SIAM Journal on Optimization 13 (3), p. 865–888. Cited by: §I-C. [49] aai-institute (2026) approx-chol: approximate cholesky factorization for graph laplacians. Note: https://github.com/aai-institute/approx-cholAccessed July 1, 2026 Cited by: §I-C. [50] L. Peng, C. Kümmerle, and R. Vidal (2023) On the convergence of irls and its variants in outlier-robust estimation. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, p. 17808–17818. Cited by: §V. [51] A. Atamtürk (2001) Flow pack facets of the single node fixed-charge flow polytope. Operations Research Letters 29 (3), p. 107–114. Cited by: §VII. [52] A. Gleixner, G. Hendel, G. Gamrath, T. Achterberg, M. Bastubbe, T. Berthold, P. Christophel, K. Jarck, T. Koch, J. Linderoth, et al. (2021) MIPLIB 2017: data-driven compilation of the 6th mixed-integer programming library. Mathematical Programming Computation 13 (3), p. 443–490. Cited by: §VII. [53] J. Moreau (1965) Proximité et dualité dans un espace hilbertien. Bulletin de la Société mathématique de France 93, p. 273–299. Cited by: §A-A. [54] M. Razaviyayn, M. Hong, and Z. Luo (2013) A unified convergence analysis of block successive minimization methods for nonsmooth optimization. SIAM Journal on Optimization 23 (2), p. 1126–1153. Cited by: §A-B. [55] J. Nocedal and S. J. Wright (2006) Numerical optimization. Springer. Cited by: §A-C. Appendix A Supplementary Material A-A Derivation of the Lasry–Lions Surrogate This section derives the Lasry–Lions double envelope surrogate used in the main paper (Section I-A): the first Moreau envelope, the exact double envelope, the logarithmic active-flow modification that yields the implemented surrogate, the admissible range for μ, continuous differentiability of the resulting surrogate, and the proof that it converges pointwise to the true fixed-charge objective as λ↓0λ 0. Throughout this subsection we fix a single, arbitrary edge and suppress its index: we write r and b for its proportional and fixed cost, and likewise J, T, S, x1x_1, x2x_2, ψ2 _2, μmax _ and Lλ,μL^λ,μ for the corresponding per-edge quantities. We assume the nonnegative-cost setting re≥0r_e≥ 0, be≥0b_e≥ 0 and re+be>0r_e+b_e>0 for every edge e∈Ee∈ E; for the fixed edge this reads r,b≥0r,b≥ 0 and r+b>0r+b>0. The per-edge objective is J(f)=r|f|+b 1|f|>0.J(f)=r|f|+b\,1\|f|>0\. For the directed FCNFP, f≥0f≥ 0, but the derivation is written in terms of t=|f|t=|f| so that it also applies to a signed-flow formulation. First Moreau envelope. For λ>0λ>0, define the Moreau envelope [53] MλJ(t)=infu≥0J(u)+(t−u)22λ.M^λJ(t)= _u≥ 0 \J(u)+ (t-u)^22λ \. There are two relevant candidates. Setting u=0u=0 gives q0(t)=t22λ.q_0(t)= t^22λ. For u>0u>0, the objective is b+ru+(t−u)22λ.b+ru+ (t-u)^22λ. Its stationary point is u=t−λru=t-λ r, which gives q1(t)=b+rt−λr22.q_1(t)=b+rt- λ r^22. The two candidates coincide when t22λ=b+rt−λr22, t^22λ=b+rt- λ r^22, or equivalently, (t−λr)2=2λb.(t-λ r)^2=2λ b. Hence, the transition point is T=λr+2λb,T=λ r+ 2λ b, and the first Moreau envelope is MλJ(t)=t22λ,0≤t≤T,b+rt−λr22,t>T.M^λJ(t)= cases t^22λ,&0≤ t≤ T,\\[8.0pt] b+rt- λ r^22,&t>T. cases (8) Lasry–Lions double envelope. For 0<μ<λ0<μ<λ, the Lasry–Lions double envelope can be written as Lλ,μJ(t)=supu≥0MλJ(u)−(t−u)22μ.L^λ,μJ(t)= _u≥ 0 \M^λJ(u)- (t-u)^22μ \. Consider first the quadratic branch of (8). The maximization problem becomes sup0≤u≤Tu22λ−(t−u)22μ. _0≤ u≤ T \ u^22λ- (t-u)^22μ \. The unconstrained maximizer is u⋆=λtλ−μ.u = λ tλ-μ. It remains in the quadratic branch when u⋆≤Tu ≤ T, or t≤x1,x1=λ−μλT.t≤ x_1, x_1= λ-μλT. Substitution gives Lλ,μJ(t)=t22(λ−μ),0≤t≤x1.L^λ,μJ(t)= t^22(λ-μ), 0≤ t≤ x_1. When t>x1t>x_1, the quadratic-branch maximizer would lie beyond T. For an intermediate range, the maximum is therefore attained at the junction u=Tu=T, giving Lλ,μJ(t)=T22λ−(t−T)22μ.L^λ,μJ(t)= T^22λ- (t-T)^22μ. The derivative of this middle branch is T−tμ. T-tμ. For the linear branch of (8), the maximizing point is u⋆=t+μru =t+μ r. It lies in the linear branch when t+μr≥Tt+μ r≥ T, or equivalently t≥S:=T−μrt≥ S:=T-μ r. Substitution gives the exact Lasry–Lions tail Lλ,μJ(t)=b+rt−(λ−μ)r22,t≥S.L^λ,μJ(t)=b+rt- (λ-μ)r^22, t≥ S. Thus, the exact Lasry–Lions double envelope is Lλ,μJ(t)=t22(λ−μ),0≤t≤x1,T22λ−(t−T)22μ,x1<t<S,b+rt−(λ−μ)r22,t≥S,L^λ,μJ(t)= cases t^22(λ-μ),&0≤ t≤ x_1,\\[9.0pt] T^22λ- (t-T)^22μ,&x_1<t<S,\\[9.0pt] b+rt- (λ-μ)r^22,&t≥ S, cases (9) which is continuously differentiable at both transition points. Logarithmic active-flow modification. The exact Lasry–Lions tail in (9) has derivative r. Consequently, an IRLS weight based on L′(t)/tL (t)/t would equal r/tr/t in the active-flow region. The fixed charge b would appear only as an additive constant and would therefore disappear from the reweighting rule. To retain a fixed-charge contribution in the active-flow weights, we replace the exact linear tail by the logarithmic continuation ψ3(t)=ψ2(x2)+r(t−x2)+αblog(tx2),t≥x2, _3(t)= _2(x_2)+r(t-x_2)+α b ( tx_2 ), t≥ x_2, (10) where ψ2(t)=T22λ−(t−T)22μ _2(t)= T^22λ- (t-T)^22μ is the middle Lasry–Lions branch. The coefficient α=λ+sλ0,s=10−2,α= λ+s _0, s=10^-2, is shared across all edges and controls the strength of the logarithmic term. Value continuity at x2x_2 holds automatically because of the additive constant ψ2(x2) _2(x_2). To obtain derivative continuity, we require ψ2′(x2)=ψ3′(x2). _2 (x_2)= _3 (x_2). Since ψ2′(t)=T−tμandψ3′(t)=r+αbt, _2 (t)= T-tμ _3 (t)=r+ α bt, the matching condition is T−x2μ=r+αbx2. T-x_2μ=r+ α bx_2. Multiplication by μx2μ x_2 gives x2(T−x2)=μrx2+μαb,x_2(T-x_2)=μ rx_2+μα b, and therefore x22−(T−μr)x2+μαb=0.x_2^2-(T-μ r)x_2+μα b=0. Defining S=T−μr,S=T-μ r, the two roots are x2±=S±S2−4μαb2.x_2^±= S± S^2-4μα b2. We select the larger root, x2=S+S2−4μαb2,x_2= S+ S^2-4μα b2, (11) because x2→Sx_2→ S as α→0α→ 0, thereby recovering the transition point of the exact Lasry–Lions envelope. The smaller root instead converges to zero. The resulting implemented surrogate is Lλ,μ(x)=|x|22(λ−μ),|x|≤x1,T22λ−(|x|−T)22μ,x1<|x|<x2,ψ2(x2)+r(|x|−x2)+αblog(|x|x2),|x|≥x2.L^λ,μ(x)= cases |x|^22(λ-μ),&|x|≤ x_1,\\[10.0pt] T^22λ- (|x|-T)^22μ,&x_1<|x|<x_2,\\[10.0pt] _2(x_2)+r(|x|-x_2)+\\ α b ( |x|x_2 ),&|x|≥ x_2. cases (12) It preserves the two Lasry–Lions parabolic regimes near the origin while introducing a logarithmic active-flow tail that retains the fixed-charge information required by IRLS. Admissible values of μ. For x2x_2 to be real, its discriminant must be nonnegative: (T−μr)2−4μαb≥0.(T-μ r)^2-4μα b≥ 0. Expanding gives r2μ2−(2Tr+4αb)μ+T2≥0.r^2μ^2-(2Tr+4α b)μ+T^2≥ 0. Let A=2Tr+4αb.A=2Tr+4α b. For r>0r>0, the admissible small-μ branch is μ≤A−A2−4r2T22r2.μ≤ A- A^2-4r^2T^22r^2. Rationalizing this expression gives the numerically stable bound μmax=2T2A+A2−4r2T2. _ = 2T^2A+ A^2-4r^2T^2. (13) Therefore, μmax=2T22Tr+4αb+(2Tr+4αb)2−4r2T2. _ = 2T^22Tr+4α b+ (2Tr+4α b)^2-4r^2T^2. (14) In the implementation, μ is capped so that μ<λμ<λ and μ≤mine∈Eμmax,eμ≤ _e∈ E _ ,e, where μmax,e _ ,e denotes (14) evaluated at the costs of edge e, so that the bound holds simultaneously on every edge. When r is zero, this bound reduces to λ/(2α)λ/(2α). Continuous Differentiability of the Lasry–Lions Surrogate. We verify that Lλ,μL^λ,μ in (12) is continuously differentiable on ℝR by differentiating each branch directly and checking that values and derivatives match at x1x_1 and x2x_2. Since every branch depends on x only through |x||x|, it suffices to check x≥0x≥ 0; the negative side follows by evenness. At x=0x=0, the quadratic branch x2/(2(λ−μ))x^2/(2(λ-μ)) is the same expression on both sides of the origin, with derivative x/(λ−μ)→0x/(λ-μ)→ 0 as x→0±x→ 0^±, so differentiability at 00 is immediate. At x1x_1. Differentiating the quadratic and middle branches gives dxx22(λ−μ) ddx x^22(λ-μ) =xλ−μ, = xλ-μ, dx[T22λ−(x−T)22μ] ddx [ T^22λ- (x-T)^22μ ] =T−xμ. = T-xμ. Using x1=(λ−μ)T/λx_1=(λ-μ)T/λ, x1λ−μ x_1λ-μ =(λ−μ)T/λ−μ=Tλ, = (λ-μ)T/λ-μ= Tλ, T−x1μ T-x_1μ =T−(λ−μ)T/λμ = T-(λ-μ)T/λμ =Tμ/λμ=Tλ, = Tμ/λμ= Tλ, so both derivatives equal T/λT/λ at x1x_1. The values match as well: since x1−T=T[(λ−μ)−λ]/λ=−μT/λx_1-T=T [(λ-μ)-λ ]/λ=-μ T/λ, T22λ−(x1−T)22μ T^22λ- (x_1-T)^22μ =T22λ−μT22λ2 = T^22λ- μ T^22λ^2 =T22λ⋅λ−μλ = T^22λ· λ-μλ =(λ−μ)T22λ2=x122(λ−μ), = (λ-μ)T^22λ^2= x_1^22(λ-μ), the last step using x12=(λ−μ)2T2/λ2x_1^2=(λ-μ)^2T^2/λ^2. This matches the quadratic branch evaluated at x1x_1, so the two branches meet with equal value and equal slope. At x2x_2. The middle branch has derivative (T−x)/μ(T-x)/μ as above; differentiating the log-tail branch ψ2(x2)+r(x−x2)+αblog(x/x2) _2(x_2)+r(x-x_2)+α b (x/x_2) gives dx[r(x−x2)+αblog(x/x2)]=r+αbx. ddx [r(x-x_2)+α b (x/x_2) ]=r+ α bx. Equality of the two derivatives at x=x2x=x_2, namely (T−x2)/μ=r+αb/x2(T-x_2)/μ=r+α b/x_2, holds if and only if, after multiplying both sides by μx2μ x_2, x2(T−μr)−x22−μαb x_2(T-μ r)-x_2^2-μα b =0 =0 ⟺x22−Sx2+μαb x_2^2-Sx_2+μα b =0, =0, using S=T−μrS=T-μ r. This is exactly the defining quadratic solved for x2x_2 in (11), so the derivative equality holds exactly, not merely in the limit. Value continuity at x2x_2 is immediate: the log-tail branch evaluates to ψ2(x2)+r(x2−x2)+αblog(x2/x2)=ψ2(x2) _2(x_2)+r(x_2-x_2)+α b (x_2/x_2)= _2(x_2), matching the middle branch value at x2x_2 by definition of ψ2(x2) _2(x_2). Degenerate case x1=x2x_1=x_2. When b=0b=0, the middle interval collapses to a point, x1=x2=(λ−μ)rx_1=x_2=(λ-μ)r (Section A-B). Neither derivative computation above assumed x1<x2x_1<x_2, so both identities still hold verbatim at this shared point: the quadratic, middle, and log-tail derivatives all equal T/λ=rT/λ=r there, and the quadratic and log-tail branches meet directly with equal value and slope even though the middle branch has zero width. Hence Lλ,μL^λ,μ is continuously differentiable everywhere, with (Lλ,μ)′(x)=sign(x)w(x)|x|,x≠0, (L^λ,μ ) (x)=sign(x)\,w(x)\,|x|, x≠ 0, where w(x)=(Lλ,μ)′(|x|)/|x|w(x)=(L^λ,μ) (|x|)/|x| is the weight function of (6) of the main paper, taken at x=0x=0 by continuous extension. Proof of Pointwise Convergence. Fix a path μ(λ)μ(λ) such that, for every λ>0λ>0, 0<μ(λ)<λ,μ(λ)≤mine∈Eμmax,e(λ),0<μ(λ)<λ, μ(λ)≤ _e∈ E _ ,e(λ), where μmax,e(λ) _ ,e(λ) denotes (14) evaluated at the costs of edge e. This is the admissible class (15); it is nonempty, since Te>0T_e>0 implies μmax,e(λ)>0 _ ,e(λ)>0, and it covers the implemented rule μ(λ)=minλ/2,mine∈Eμmax,e(λ)μ(λ)= \λ/2,\ _e∈ E _ ,e(λ)\ of Section V as a special case. Since 0<μ(λ)<λ0<μ(λ)<λ, it follows that μ(λ)→0μ(λ)→ 0 as λ↓0λ 0. Let further α(λ)=λ+sλ0,α(λ)= λ+s _0, where s>0s>0 and λ0>0 _0>0 are fixed as in Section I-A. Recall the standing assumptions r,b≥0r,b≥ 0 and r+b>0r+b>0. We show that Lλ,μ(λ)(x)→J(x)=r|x|+b 1|x|>0L^λ,μ(λ)(x)→ J(x)=r|x|+b\,1\|x|>0\ for every fixed x∈ℝx as λ↓0λ 0. Summing over edges then gives JLLλ,μ(λ)(f)→Jtrue(f)J_L^λ,μ(λ)(f)→ J_true(f) for every fixed f∈ℱf . Case x=0x=0. The origin always lies in the quadratic branch, and Lλ,μ(0)=0=J(0)L^λ,μ(0)=0=J(0) for every λ, so the limit is immediate. Transition points vanish. Since r,b≥0r,b≥ 0, T(λ) T(λ) =λr+2λb⟶0, =λ r+ 2λ b 0, S(λ) S(λ) =(λ−μ(λ))r+2λb. =(λ-μ(λ))r+ 2λ b. Since μ(λ)r≥0μ(λ)r≥ 0 and λ−μ(λ)>0λ-μ(λ)>0 together with r+b>0r+b>0, we have 0<S(λ)≤T(λ)0<S(λ)≤ T(λ). Hence S(λ)→0S(λ)→ 0 by the squeeze theorem. The admissibility condition guarantees that S2−4μαb≥0S^2-4μα b≥ 0. Since μαb≥0μα b≥ 0, we also have 0 0 ≤S2−4μαb≤S2, ≤ S^2-4μα b≤ S^2, 0 0 ≤S2−4μαb≤S. ≤ S^2-4μα b≤ S. It follows directly from the definition of the larger root that S2≤x2=S+S2−4μαb2≤S. S2≤ x_2= S+ S^2-4μα b2≤ S. Therefore x2(λ)→0x_2(λ)→ 0 by the squeeze theorem. Finally, 0≤x1≤x20≤ x_1≤ x_2, so x1(λ)→0x_1(λ)→ 0 as well. Case x≠0x≠ 0, b>0b>0. Because x2(λ)→0x_2(λ)→ 0, there is λ∗>0λ >0 with x2(λ)<|x|x_2(λ)<|x| for all 0<λ<λ∗0<λ<λ , so x lies in the log-tail branch throughout this range: Lλ,μ(x)=ψ2(x2)+r(|x|−x2)+αblog(|x|/x2).L^λ,μ(x)= _2(x_2)+r(|x|-x_2)+α b (|x|/x_2). We take each term to its limit. The linear term satisfies r(|x|−x2)→r|x|r(|x|-x_2)→ r|x| since x2→0x_2→ 0. For the logarithmic term, S≥2λbS≥ 2λ b and x2≥S/2x_2≥ S/2 give x2≥λb/2x_2≥ λ b/2. Moreover, 0<α≤λ/(sλ0)0<α≤λ/(s _0). Hence, for all sufficiently small λ, 0 0 ≤αblog(|x|x2) ≤α b \! ( |x|x_2 ) ≤λbsλ0log(|x|2λb) ≤ λ bs _0 \! (|x| 2λ b ) =bsλ0[λlog(|x|2b)+λ2log(1λ)]⟶0, = bs _0 [λ \! (|x| 2b )+ λ2 \! ( 1λ ) ] 0, where the last limit uses λlog(1/λ)→0λ (1/λ)→ 0. The squeeze theorem therefore gives αblog(|x|/x2)→0α b (|x|/x_2)→ 0. It remains to determine the limit of the junction value ψ2(x2) _2(x_2). The defining relation x22−Sx2+μαb=0x_2^2-Sx_2+μα b=0 gives x2(S−x2)=μαbx_2(S-x_2)=μα b. Since x2>0x_2>0, T−x2 T-x_2 =(T−S)+(S−x2) =(T-S)+(S-x_2) =μr+μαbx2 =μ r+ μα bx_2 =μ(r+αbx2). =μ (r+ α bx_2 ). The preceding lower bound on x2x_2 yields the explicit estimate 0≤αbx2≤2αb2λb=2bλ+sλ0≤2bλsλ0⟶0.0≤ α bx_2≤ 2α b 2λ b= 2bλ+s _0≤ 2bλs _0 0. Thus r+αb/x2→r+α b/x_2→ r and, since μ→0μ→ 0, (x2−T)22μ=μ2(r+αbx2)2⟶0. (x_2-T)^22μ= μ2 (r+ α bx_2 )^2 0. On the other hand, the exact identity T22λ=λr22+r2λb+b T^22λ= λ r^22+r 2λ b+b shows that T2/(2λ)→bT^2/(2λ)→ b. Hence ψ2(x2)=T22λ−(x2−T)22μ⟶b. _2(x_2)= T^22λ- (x_2-T)^22μ\ \ b. Adding the three limits gives Lλ,μ(λ)(x)→b+r|x|=J(x)L^λ,μ(λ)(x)→ b+r|x|=J(x). Case x≠0x≠ 0, b=0b=0. Then r>0r>0, and x1=x2=(λ−μ)r→0x_1=x_2=(λ-μ)r→ 0 since 0<λ−μ<λ→00<λ-μ<λ→ 0, so x again eventually lies in the log-tail branch, where αb=0α b=0 exactly. Here S=x2S=x_2 and T=S+μrT=S+μ r, so ψ2(x2) _2(x_2) =T22λ−(μr)22μ = T^22λ- (μ r)^22μ =(λ−μ)r22→0=b, = (λ-μ)r^22→ 0=b, and r(|x|−x2)→r|x|r(|x|-x_2)→ r|x|, so again Lλ,μ(λ)(x)→r|x|=J(x)L^λ,μ(λ)(x)→ r|x|=J(x). In every case Lλ,μ(λ)(x)→J(x)L^λ,μ(λ)(x)→ J(x) as λ↓0λ 0, which proves the claim. A-B Quadratic Majorization at Fixed Smoothing Parameters We next justify the majorization–minimization interpretation of the quadratic IRLS model. Throughout this subsection, the smoothing parameters (λ,μ)(λ,μ) are fixed. We assume the standard nonnegative-cost setting re,be≥0r_e,b_e≥ 0 and re+be>0r_e+b_e>0 for every e∈Ee∈ E, and note that the implemented coefficient satisfies 0<α≤10<α≤ 1. On instance families where some arcs have re=be=0r_e=b_e=0 exactly, this assumption is violated as stated; the implementation floors rer_e at 10−1210^-12 before evaluating the surrogate, its weights, and μmax,e _ ,e from (14), which restores re+be>0r_e+b_e>0 for every edge and brings such arcs within the scope of the proposition below. We take 0<μ<λ,μ≤mine∈Eμmax,e(λ).0<μ<λ, μ≤ _e∈ E _ ,e(λ). (15) The μmax,e _ ,e cap makes the discriminant defining x2,ex_2,e nonnegative. To show that the transition points in Section I-A of the main paper are also ordered, set Se=Te−μre=(λ−μ)re+2λbe>0S_e=T_e-μ r_e=(λ-μ)r_e+ 2λ b_e>0 and pe(x)=x2−(Te−μre)x+μαbe.p_e(x)=x^2-(T_e-μ r_e)x+μα b_e. Since the discriminant is nonnegative, the root sum is Se>0S_e>0, and the root product is μαbe≥0μα b_e≥ 0, both roots are real and nonnegative, and x2,ex_2,e is the larger one. If be=0b_e=0, then re>0r_e>0 and pe(x)=x(x−(λ−μ)re)p_e(x)=x (x-(λ-μ)r_e ), so x1,e=x2,e=(λ−μ)rex_1,e=x_2,e=(λ-μ)r_e. Suppose now that be>0b_e>0, and set qe=2be/λq_e= 2b_e/λ. Direct substitution gives pe(x1,e)=μ(αbe−qex1,e).p_e(x_1,e)=μ(α b_e-q_ex_1,e). If pe(x1,e)≤0p_e(x_1,e)≤ 0, then x1,ex_1,e lies between the roots and hence x1,e≤x2,ex_1,e≤ x_2,e. If pe(x1,e)>0p_e(x_1,e)>0, then qex1,e<αbe≤be=λqe22,sox1,e<λqe2≤Se2.q_ex_1,e<α b_e≤ b_e= λ q_e^22, x_1,e< λ q_e2≤ S_e2. Because Se/2S_e/2 is the vertex of the convex quadratic pep_e, this places x1,ex_1,e to the left of the smaller root (or of the double root), and again x1,e<x2,ex_1,e<x_2,e. Finally, if Δe=Se2−4μαbe≥0 _e=S_e^2-4μα b_e≥ 0, then x2,e=Se+Δe2≤Se≤Te.x_2,e= S_e+ _e2≤ S_e≤ T_e. Consequently, the capped domain gives 0<x1,e≤x2,e≤Te0<x_1,e≤ x_2,e≤ T_e for every edge. Proposition 1 (Fixed-parameter quadratic majorization). Suppose that the stated cost assumptions and (15) hold. Let f(k)∈ℱf^(k) , let we(k)=we(fe(k))w_e^(k)=w_e(f_e^(k)) be the exact weight given by (6) of the main paper, and let f(k+1)f^(k+1) be the exact update (c.f. (5) of the main paper). Then Wk≻0W_k 0, so ⟨⋅,⋅⟩Wk ·,· _W_k is an inner product, and the model Qλ,μ(⋅∣f(k))Q^λ,μ(\,· f^(k)) defined in (4) of the main paper is a global majorizer of the fixed-parameter smoothed objective: JLLλ,μ(f) J_L^λ,μ(f) ≤Qλ,μ(f∣f(k))for every f∈ℝm, ≤ Q^λ,μ(f f^(k)) every f ^m, (16) Qλ,μ(f(k)∣f(k)) Q^λ,μ(f^(k) f^(k)) =JLLλ,μ(f(k)). =J_L^λ,μ(f^(k)). (17) The two functions also have the same gradient at f(k)f^(k). Consequently, JLLλ,μ(f(k+1)) J_L^λ,μ(f^(k+1)) ≤Qλ,μ(f(k+1)∣f(k)) ≤ Q^λ,μ(f^(k+1) f^(k)) (18) ≤Qλ,μ(f(k)∣f(k))=JLLλ,μ(f(k)). ≤ Q^λ,μ(f^(k) f^(k))=J_L^λ,μ(f^(k)). More precisely, JLLλ,μ(f(k))−JLLλ,μ(f(k+1)) J_L^λ,μ(f^(k))-J_L^λ,μ(f^(k+1)) (19) ≥12⟨f(k+1)−f(k),f(k+1)−f(k)⟩Wk. ≥ 12 f^(k+1)-f^(k),f^(k+1)-f^(k) _W_k. Proof. The ordering established above implies positivity of every exact analytical weight. In the quadratic regime, we(x)=1/(λ−μ)>0w_e(x)=1/(λ-μ)>0. In a nonempty transition regime, 0<x<x2,e≤Te0<x<x_2,e≤ T_e, so we(x)=(Te−x)/(μx)>0w_e(x)=(T_e-x)/(μ x)>0. In the active-flow regime, we(x)=re/x+αbe/x2>0w_e(x)=r_e/x+α b_e/x^2>0 because re+be>0r_e+b_e>0 and α>0α>0. The weight at zero is defined by the quadratic branch. Hence Wk≻0W_k 0. For each edge, define the scalar function ge(t)=LLeλ,μ(t),t≥0.g_e(t)=L_e^λ,μ( t), t≥ 0. (20) For t>0t>0, differentiation of (2) of the main paper gives ge′(t)=12(λ−μ),0<t<x1,e2,Te−t2μt,x1,e2<t<x2,e2,re2t+αbe2t,t>x2,e2.g_e (t)= cases 12(λ-μ),&0<t<x_1,e^2,\\[8.0pt] T_e- t2μ t,&x_1,e^2<t<x_2,e^2,\\[10.0pt] r_e2 t+ α b_e2t,&t>x_2,e^2. cases (21) The values at the junctions are understood by continuous extension. Inside the three respective regions, ge′(t)=0,0<t<x1,e2,−Te4μt3/2,x1,e2<t<x2,e2,−re4t3/2−αbe2t2,t>x2,e2.g_e (t)= cases0,&0<t<x_1,e^2,\\[3.0pt] - T_e4μ t^3/2,&x_1,e^2<t<x_2,e^2,\\[10.0pt] - r_e4t^3/2- α b_e2t^2,&t>x_2,e^2. cases (22) Since re,be≥0r_e,b_e≥ 0, each expression in (22) is nonpositive. Moreover, the one-sided derivatives agree at both transition points. Indeed, Te−x1,eμx1,e=1λ−μ, T_e-x_1,eμ x_1,e= 1λ-μ, (23) while the defining equation x2,e2−(Te−μre)x2,e+μαbe=0x_2,e^2-(T_e-μ r_e)x_2,e+μα b_e=0 (24) implies Te−x2,eμx2,e=rex2,e+αbex2,e2. T_e-x_2,eμ x_2,e= r_ex_2,e+ α b_ex_2,e^2. (25) If x1,e=x2,ex_1,e=x_2,e, the middle interval is empty, and these identities show that the derivatives of the quadratic and active-flow branches agree directly at the common transition point. Value continuity follows from the matched middle-branch value used to define the active-flow constant. The right derivative at zero is ge′(0)=1/[2(λ−μ)]g_e (0)=1/[2(λ-μ)]. Hence geg_e is continuously differentiable and concave on [0,∞)[0,∞). The tangent inequality for a concave function now yields, for every fe∈ℝf_e , LLeλ,μ(fe)=ge(fe2) _e^λ,μ(f_e)=g_e(f_e^2) ≤ge((fe(k))2) ≤ g_e\! ((f_e^(k))^2 ) +ge′((fe(k))2)(fe2−(fe(k))2). +g_e \! ((f_e^(k))^2 ) (f_e^2-(f_e^(k))^2 ). (26) By (21) and (6) of the main paper, 2ge′((fe(k))2)=we(k),2g_e \! ((f_e^(k))^2 )=w_e^(k), (27) where the first branch gives the same identity when fe(k)=0f_e^(k)=0. Summing (26) over the edges proves (16). Substitution of f=f(k)f=f^(k) gives (17), and ∇fQλ,μ(f(k)∣f(k))=Wkf(k)=∇JLLλ,μ(f(k)), _fQ^λ,μ(f^(k) f^(k))=W_kf^(k)=∇ J_L^λ,μ(f^(k)), which establishes first-order tangency. It remains to prove the last claim. The terms in (4) of the main paper that are independent of f do not affect the minimizer, so (5) of the main paper is precisely the IRLS weighted quadratic update. This gives the descent chain (18). Since ℱF is convex, optimality of f(k+1)f^(k+1) for Qλ,μ(⋅∣f(k))Q^λ,μ(\,· f^(k)) also gives ⟨f(k+1),f(k)−f(k+1)⟩Wk≥0. f^(k+1),f^(k)-f^(k+1) _W_k≥ 0. Expanding the quadratic difference therefore shows that Qλ,μ(f(k)∣f(k))−Qλ,μ(f(k+1)∣f(k)) Q^λ,μ(f^(k) f^(k))-Q^λ,μ(f^(k+1) f^(k)) ≥12⟨f(k+1)−f(k),f(k+1)−f(k)⟩Wk. ≥ 12 f^(k+1)-f^(k),f^(k+1)-f^(k) _W_k. Combining this inequality with the majorization bound proves (19). ∎ Idealized fixed-parameter convergence. If the exact update is iterated indefinitely at one fixed admissible pair (λ,μ)(λ,μ) satisfying (15), then compactness of ℱF and continuity of JLLλ,μJ_L^λ,μ make the objective bounded below. Summing (19) gives ∑k‖f(k+1)−f(k)‖Wk2<∞ _k\|f^(k+1)-f^(k)\|^2_W_k<∞, so f(k+1)−f(k)→0f^(k+1)-f^(k)→ 0 in the WkW_k-weighted norm, and standard M arguments imply that every accumulation point is stationary for JLLλ,μJ_L^λ,μ over ℱF [54]. This idealized statement does not assert convergence of a finite continuation stage or of the full parameter-varying algorithm. Numerical qualification. Proposition 1 concerns a fixed admissible pair, the exact analytical weights in (6) of the main paper, and an exact solution of the weighted quadratic subproblem. If a flow magnitude is floored or a weight is clipped for numerical protection, the resulting quadratic need not retain either global majorization or first-order tangency. With exact weights, an inexact subproblem solve retains the descent conclusion whenever Qλ,μ(f(k+1)∣f(k))≤Qλ,μ(f(k)∣f(k))Q^λ,μ(f^(k+1) f^(k))≤ Q^λ,μ(f^(k) f^(k)). Because the implementation clips weights and solves the QP numerically, it does not infer descent from the proposition. Instead it evaluates JLLλ,μJ_L^λ,μ directly at the fixed smoothing pair after a trial step and backtracks or rejects any step that fails the strict acceptance rule specified in Section A-D. Scope with respect to continuation. The proposition is step-wise: it compares consecutive iterates using the same fixed pair (λ,μ)(λ,μ). It does not imply a comparison after a continuation update. With the parameter dependence displayed explicitly, the jump JLLλk+1,μk+1(f(k+1))−JLLλk,μk(f(k+1))J_L _k+1, _k+1(f^(k+1))-J_L _k, _k(f^(k+1)) has no fixed sign in general. Non-positivity of this quantity is desirable, but cannot be guaranteed as downstream effect of the particular shape of the Lasry–Lions functional Lλ,μL^λ,μ and the continuation scheme for λ and μ. A-C Weighted-QP Solver This section documents the weighted-QP solver of Section I-C in full: the dual semismooth Newton (SSN) scheme implemented, the empirical scaling of its inner iteration counts across problem sizes, and the primal-dual interior-point method (IPM) used as a fallback when SSN fails to converge. Dual semismooth Newton. Each weighted-QP solve is initialized with the dual potential and active-set masks from the previous accepted outer iterate. On the first call, no such state is available, so the input dual potential is set to y=0y=0 and the masks are left empty. The solver does not form a Newton system at this initial value: since D(0)=0D(0)=0, a literal Newton step from y=0y=0 would make (7) of the main paper degenerate. Instead, it first constructs a dual predictor from the current primal flow f, updates the corresponding active-set masks, and then applies the semismooth Newton correction. Predictor. The free and bound-active arcs are first estimated directly from the current primal flow (f) by thresholding it against the lower and upper bounds, rather than from the dual potential (y). If the resulting free subgraph is component-wise balanced with respect to the residual demand induced by the bound-active arcs, the solver performs one Laplacian solve on that subgraph. The resulting primal–dual candidate is accepted without any semismooth Newton iterations if it is primal feasible and the recovered dual variables satisfy the sign conditions associated with the estimated free and bound-active sets. Regularized correction. If the predictor is rejected or infeasible, the solver returns to the warm-started dual potential y and applies semismooth Newton corrections. At each step, arcs are classified as free, bound-active, or near a kink, with the latter lying within a numerical tolerance of either bound. On connected components for which the free and bound-active partition is component-wise balanced, the solver uses the exact generalized derivative D(y)D(y). Otherwise, it regularizes the derivative by adding a penalty η/weη/w_e on near-kink and bound-active arcs. The regularization parameter is initialized at η=10−2η=10^-2 and increased adaptively if the linear system cannot be solved accurately, the computed direction does not reduce the dual residual, or the free-subgraph solve remains poorly converged. The regularized Newton matrix is again a graph Laplacian and therefore remains singular up to an additive constant on each connected component. For a component-wise compatible right-hand side, the Newton equation is understood modulo these gauges (equivalently, through a Laplacian pseudoinverse); the implementation mean-centers the right-hand side and returned potential, which fixes the gauge when the regularized conductance graph is connected. The regularization still provides an inexact direction when the unregularized free-subgraph system is degenerate. Line search. Let ϕ(β)φ(β) denote the concave dual objective restricted to the search path y+βΔy+β y. Its directional derivative is ϕ′(β)=Δy⊤R(y+βΔy).φ (β)= y R(y+β y). This scalar function is continuous, nonincreasing, and piecewise linear in β. The solver accepts the full step if ϕ′(1)≥0φ (1)≥ 0; otherwise, it locates a zero of ϕ′φ in (0,1)(0,1) before updating the dual potential, thereby avoiding overshooting an intervening kink. Inner iteration counts. Each outer IRLS iteration issues one weighted-QP solve, and each dual semismooth Newton step inside that solve requires one PCG solve of a free-subgraph Laplacian system, using a freshly computed approximate-Cholesky factorization as preconditioner, at cost O~(mlog(1/ϵ)) O(m (1/ε)). Writing K for the number of outer iterations and s¯ s for the average number of Newton steps per subproblem, the dominant cost of a single IRLS run is O~(Ks¯mlog(1/ϵ)) O(K\, s\,m (1/ε)). Since K is capped at 200200 by construction, the practical question is how s¯ s behaves as m grows. Table I reports these counters for one Layer-1 run on an instance from each of the nine scale families, which span a 200200-fold range in m and account for 270270 of the 410410 instances. All nine runs terminate normally, so the counts are complete. Over this range s¯ s grows slowly, from 1.001.00 at m=100m=100 to 2.912.91 at m=2×104m=2× 10^4, and never exceeds 3.593.59; a least-squares fit gives s¯≈−1.03+0.41logm s≈-1.03+0.41 m, so a 200200-fold increase in m costs less than a threefold increase in Newton steps per subproblem. Since s¯ s itself grows only logarithmically in m, the per-run cost over these instances is therefore close to O~(Km) O(K\,m), with the log(1/ϵ) (1/ε) and logm m factors both absorbed into the O~ O. The number of outer iterations is likewise small, between 1212 and 2626, so no run approaches the iteration cap. On every instance the number of Laplacian solves exceeds the number of Newton steps by exactly K: the warm-start predictor contributes one additional solve per outer iteration, and no subproblem requires a restart beyond it. These figures are for a single instance per family and for Layer 1 alone; Layers 2 and 3 invoke the same solver repeatedly, so their totals scale with the number of restarts and restricted supports rather than with s¯ s. TABLE I: Inner iteration counts for one Layer-1 IRLS run per scale family. K is the number of outer iterations, which equals the number of weighted-QP solves; s¯ s is the average number of Newton steps per solve; Lapl. counts Laplacian solver calls. n m K Newton s¯ s Lapl. 2020 100100 12 12 1.00 24 6060 400400 12 17 1.42 29 120120 1,5001,500 20 23 1.15 43 150150 2,5002,500 15 31 2.07 46 200200 4,0004,000 17 61 3.59 78 250250 7,5007,500 21 53 2.52 74 300300 9,0009,000 26 57 2.19 83 500500 10,00010,000 21 65 3.10 86 1,0001,000 20,00020,000 22 64 2.91 86 Interior-point method formulation. When the dual SSN solver described above fails to converge to a flow-conservation residual below 10−610^-6, the same weighted quadratic subproblem is instead solved by a primal-dual interior-point method (IPM), which shares the weighted graph-Laplacian structure of the SSN Newton systems. The weighted quadratic subproblem can be written as minf _f 12f⊤Wf 12f Wf (28) s.t. .t. Cf=d, Cf=d, 0≤f≤c. 0≤ f≤ c. To solve (28), we use a primal-dual interior-point method [55]. Introducing dual variables y for the equality constraints and nonnegative dual variables s and q for the lower and upper bounds, respectively, the Lagrangian is ℒ(f,y,s,q)=12f⊤Wf−y⊤(Cf−d)−s⊤f−q⊤(c−f).L(f,y,s,q)= 12f Wf-y (Cf-d)-s f-q (c-f). (29) The KKT residuals are rb=Cf−d,r_b=Cf-d, (30) and rc=−C⊤y+Wf−s+q.r_c=-C y+Wf-s+q. (31) The complementarity conditions are fese=0,f_es_e=0, (32) and (ce−fe)qe=0.(c_e-f_e)q_e=0. (33) In the interior-point method, (32) and (33) are relaxed as fese≈β,f_es_e≈β, (34) and (ce−fe)qe≈β,(c_e-f_e)q_e≈β, (35) where β>0β>0 is the barrier target. Then, we define the complementarity residuals as rc1 r_c1 =f⊙s−β, =f s- 1, (36) rc2 r_c2 =(c−f)⊙q−β. =(c-f) q- 1. (37) Linearizing the perturbed KKT system and eliminating the slack directions yields the reduced system (CΘC⊤)Δy=−C(Θζ)−rb, (C C ) y=-C( ζ)-r_b, (38) where ζ=−rc−F−1rc1+(U−F)−1rc2ζ=-r_c-F^-1r_c1+(U-F)^-1r_c2 (39) and Θ=H−1. =H^-1. (40) The Newton system involves the diagonal matrix H=W+F−1S+(U−F)−1Q,H=W+F^-1S+(U-F)^-1Q, (41) where F=diag(f),S=diag(s),Q=diag(q),U=diag(c).F=diag(f),S=diag(s),Q=diag(q),U=diag(c). (42) CΘC⊤C C is a weighted graph Laplacian, so the method can use Laplacian solvers to accelerate the subproblem solve. Once Δy y is computed, the remaining directions are recovered from Δf f =Θ(ζ+C⊤Δy), = (ζ+C y ), (43) Δs s =F−1(−rc1−SΔf), =F^-1 (-r_c1-S f ), (44) Δq q =(U−F)−1(−rc2+QΔf). =(U-F)^-1 (-r_c2+Q f ). (45) A step length is then selected to preserve strict positivity of the primal-dual variables, and the iterates are updated accordingly. The overall method has a nested structure. The outer loop corresponds to the IRLS reweighting step, while the inner loop solves the resulting convex quadratic subproblem. At outer iteration k, given the current flow f(k)f^(k), we compute the edge weights we(k)=wλk,μk(fe(k)),∀e∈E,w_e^(k)=w_ _k, _k(f_e^(k)), ∀ e∈ E, (46) using (6) of the main paper. These weights define the quadratic model (28), which is then solved by the interior-point method to obtain f(k+1)f^(k+1). The surrogate parameters are subsequently updated, producing a new smoothed model for the next outer iteration. Hence, the outer loop is responsible for approximating the original non-convex objective, while the inner loop computes the minimizer of the current quadratic model. For exact analytical weights in the capped domain, this model is the majorizer of Proposition 1; the clipped numerical model is safeguarded by direct surrogate-value acceptance. A-D Algorithmic Details This section completes the specification of the IRLS core method of Section IV: the support-refit mechanism, termination conditions, numerical clamping, the step-acceptance rule, stagnation tests, and candidate feasibility evaluation. Support refit. The refit is attempted when ρgap≥10 _gap≥ 10 and 0<Kgap≤max4n,10240<K_gap≤ \4n,1024\. Its threshold f↑f↓ f_ f_ is the midpoint of the gap on a logarithmic scale. An attempt is skipped when the support repeats the previous one, and a cost model is discarded when the restricted network cannot route the demand. We solve it under two cost models, rer_e and re+be/maxfe,1r_e+b_e/ \f_e,1\, and keep the smaller JtrueJ_true: the former minimizes variable cost on S, while the latter rewards closing weakly loaded arcs and remains informative when r≡0r≡ 0. Both models use the same integer-cost network-simplex routine as the stage-boundary candidate. The result enters the best feasible solution only if it improves the incumbent by more than max10−9,10−6|Jtrue⋆| \10^-9,10^-6|J _true|\, so the refit cannot degrade the reported objective. The final refit uses the visible support e:|fe|>τ\e:|f_e|>τ\ of the last iterate. Each stage boundary additionally produces a full-graph feasible candidate, obtained by solving a min-cost flow with unit costs re+be/maxfe,1r_e+b_e/ \f_e,1\ on the active arcs and re+be/maxce,1r_e+b_e/ \c_e,1\ elsewhere. This candidate is kept only if it is feasible and improves JtrueJ_true. Termination. The run exits once ρgap≥100 _gap≥ 100, 0<Kgap≤max4n,10240<K_gap≤ \4n,1024\, KgapK_gap has taken the same value on four consecutive iterations, and the best feasible solution has not improved for eight iterations. The separation required here is two decades rather than the factor of five needed to trust λgap _gap, because exit is irreversible while an inaccurate target only perturbs the rate of sharpening. The 10241024 floor keeps this trigger permissive on small networks, where 4n4n alone would be too restrictive; for larger graphs (n>256n>256) the bound scales as 4n4n, a small multiple of the n−1n-1 arcs of a basic feasible flow. The run also stops when the surrogate fails to decrease for 1515 iterations with relative flow change below 10−410^-4, when λ reaches λmin _ , or at the limit of 200200 reweighting iterations or the time budget. Clamping. The IRLS weights of (6) of the main paper are clipped to [10−12,1012][10^-12,10^12] before each solve, and the variable costs entering the surrogate, the weights, and μmax,e _ ,e are floored at re←maxre,10−12r_e← \r_e,10^-12\. The floor matters on the pure fixed-charge families, for which r≡0r≡ 0 and several surrogate quantities would otherwise be degenerate. Objectives are evaluated with the unmodified costs, so the reported JtrueJ_true is unaffected. As noted after Proposition 1, a clipped weight may forfeit global majorization, which is why acceptance of a step is tested directly rather than assumed. Step acceptance. A trial iterate f~ f is evaluated using JLLλk,μk(⋅)J_L _k, _k(·) at the current parameter pair and accepted only if JLLλk,μk(f~) J_L _k, _k( f) ≤JLLλk,μk(f(k)) ≤ J_L _k, _k(f^(k)) (the numerical tolerance mentioned in the main paper is set to 00 in the implementation). Otherwise the step is backtracked along f(k)+β(f~−f(k))f^(k)+β( f-f^(k)) for β∈12,14,18,116,132β∈\ 12, 14, 18, 116, 132\, and the first β satisfying the same test is taken. If none does, the step is rejected and f(k+1)=f(k)f^(k+1)=f^(k). Whenever β<1β<1, the dual warm start is discarded and the semismooth Newton state is reinitialized to y=0y=0 and s=q=s=q=1, because the retained active set no longer matches the accepted iterate. Stagnation tests. Both counters use a relative threshold of the same form. The surrogate counter advances unless JLLλ,μJ_L^λ,μ decreases by more than max10−9,10−6|JLLλ,μ(f(k))| \10^-9,10^-6|J_L^λ,μ(f^(k))|\. The true-objective counter advances unless a feasible candidate improves the best value by more than max10−9,10−6|Jtrue⋆| \10^-9,10^-6|J _true|\. Candidate evaluation. All candidates are scored on the visible flow obtained by setting to zero every arc with fe≤τ=10−6f_e≤τ=10^-6. A candidate is feasible when the visible flow satisfies ‖Cf−d‖∞≤10−6\|Cf-d\|_∞≤ 10^-6, together with capacity and nonnegativity violations of at most 10−610^-6, and only feasible candidates are retained. Applying the threshold before the feasibility test is what makes an interior iterate ineligible, and is the reason the stage-boundary candidates are needed. Those candidates are passed to an integer-cost network-simplex routine, so arc costs are scaled by 10610^6 and rounded. A candidate is skipped when the node balances are not integral to within 10−610^-6 or do not sum to zero, which does not occur on the families reported here. A-E Smoothing Continuation: Implementation Details This section completes the specification of Section V. Table IV collects the constants of the schedule, which were held fixed across all reported runs. These constants were fixed once via manual tuning on a small subset of instances prior to the full 410410-instance run, rather than via a systematic search over a stated range of values. Initialization. Let fmedf_med be the median magnitude among initial flows exceeding one unit, and let bmedb_med be the median positive fixed cost. The choice λ=fmed2/(2bmed)λ=f_med^2/(2b_med) places TeT_e at fmedf_med when rer_e is negligible. Starts that already carry combinatorial support information, namely those derived from MCF, DSSP, and ADSSP solutions and from restart or cached supports, begin one decade below this level in flow units, and therefore two decades below in λ because Te∝λT_e λ: λ0=10−2maxfmed22bmed, 1. _0=10^-2 \ f_med^22b_med,\;1 \. The active flows of the start then lie in the logarithmic tail, so that sharpening acts on the incumbent support rather than dissolving it. If fewer than two arcs carry flow, this scale is undefined and we use λ0=10−2max3J0,1 _0=10^-2 \3J_0,1\ instead, where J0J_0 is the realized cost of the start. Diffuse quadratic-programming starts have no support worth preserving, and are deliberately over-smoothed with λ0=102maxmaxe|fe(0)|,1 _0=10^2 \ _e|f^(0)_e|,1\. This places every arc in the quadratic head, so the first subproblem is an almost uniformly weighted projection onto ℱF. The initial ratio is μ0=minθ0λ0,mineμmax,e(λ0) _0= \ _0 _0, _e _ ,e( _0)\ with θ0=0.3 _0=0.3 and θ0=0.7 _0=0.7 in the two cases. The cap is inactive for θ0=0.3 _0=0.3 but can reduce the diffuse-start ratio below 0.70.7; both choices satisfy the capped domain of Proposition 1. The special θ0 _0 applies to the initial weighted-QP solve, after which every continuation update uses μ=λ/2μ=λ/2. The log-tail coefficient is anchored to the start through α=λ/(λ+sλ0)α=λ/(λ+s _0) with s=10−2s=10^-2, so that α≈0.99α≈ 0.99 initially and α=1/2α=1/2 after two decades of continuation, past which the tail relaxes toward the exact Lasry–Lions envelope recovered as α→0α→ 0. Gap statistic. The sorted flow magnitudes are floored at λmin _ before the ratios are formed, and we write f↑=|f|(Kgap)f_ =|f|_(K_gap) and f↓=|f|(Kgap+1)f_ =|f|_(K_gap+1) for the magnitudes bracketing the gap. The representative fixed cost brepb_rep is the flow-weighted average of beb_e over the arcs within one octave of the gap, that is |fe|∈[f↓/2, 2f↑]|f_e|∈[f_ /2,\,2f_ ] with be>10−3b_e>10^-3, and the median positive beb_e when that set is empty. The gap is trusted, and λgap _gap used, only when ρgap≥5 _gap≥ 5, Kgap≤m/2K_gap≤ m/2, and KgapK_gap has taken the same value on three consecutive iterations. The last two conditions reject a ratio that lies in the interior of a broad flow distribution. Stage rule. A stage ends when the relative flow change falls below 10−310^-3, when the line search rejects the step, or after six weighted-QP solves. At a stage boundary, λ+=maxminλ/4,λgap,λ/100,λ^+= \ \λ/4,\; _gap\,\;λ/100 \, while between boundaries a trusted target may cap the level by at most a factor of five, that is λ+=maxminλ,λgap,λ/5λ^+= \ \λ, _gap\,λ/5\. After flooring at λmin _ , we set μ+=minλ+/2,mineμmax,e(λ+)μ^+= \λ^+/2, _e _ ,e(λ^+)\. The second bound is never active, since μmax,e(λ)≥λ/2 _ ,e(λ)≥λ/2 holds for every edge whenever α≤1α≤ 1. TABLE IV: Constants of the smoothing continuation schedule. Quantity Value Post-initialization θ=μ/λθ=μ/λ 1/21/2 Initial-QP requested θ0 _0 0.30.3 / 0.70.7 λ0 _0, sparse start 10−2maxfmed2/(2bmed),110^-2 \f_med^2/(2b_med),1\ λ0 _0, dense start 102maxmaxe|fe(0)|,110^2 \ _e|f^(0)_e|,1\ s (in α) 10−210^-2 λmin _ 10−910^-9 Stage decay λ/4λ/4 Stage clamp λ/100λ/100 Intra-stage clamp λ/5λ/5 Stage QP cap 66 solves Stage-close flow change 10−310^-3 Trust: ρgap _gap, KgapK_gap, stability 55, m/2m/2, 33 iters Refit trigger: ρgap _gap, KgapK_gap 1010, max4n,1024 \4n,1024\ Refit support threshold f↑f↓ f_ f_ Refit cost models rer_e and re+be/maxfe,1r_e+b_e/ \f_e,1\ Exit: ρgap _gap, KgapK_gap, stability 100100, max4n,1024 \4n,1024\, 44 iters Exit patience 88 iters Surrogate stagnation 1515 iters, 10−410^-4 brepb_rep window [f↓/2, 2f↑][f_ /2,\,2f_ ], be>10−3b_e>10^-3 Weight clip [10−12,1012][10^-12,10^12] Variable-cost floor 10−1210^-12 Backtracking factors 12,14,18,116,132 12, 14, 18, 116, 132 Activation threshold τ 10−610^-6 Feasibility tolerance 10−610^-6 Refit cost scaling 10610^6 Reweighting iteration limit 200200 A-F Algorithmic Layers This section completes the specification of the three algorithmic layers of Section VI: initialization strategies for multi-start IRLS, the objective-driven perturbation-restart mechanism, and the implementation details of the anchor-union restricted search. A-F1 Initialization Strategies The multi-start IRLS layer uses several feasible initial flows. Quadratic warm starts solve weighted least-squares flow problems over ℱF. The linear-cost-weighted quadratic flow initialization solves min∑e∈Ef∈ℱrefe2. _f _e∈ Er_ef_e^2. The quadratic initialization weighted by combined linear and capacity-amortized fixed cost solves min∑e∈Ef∈ℱ(re+bemaxce,1)fe2. _f _e∈ E (r_e+ b_e \c_e,1\ )f_e^2. The min-cost-flow warm starts solve linearized flow problems with combined per-unit costs. One start uses total-flow-amortized fixed charges, while the capacity-aware start uses costs re+bemaxce,1,r_e+ b_e \c_e,1\, which favors high-capacity arcs with low amortized fixed charge. We also include a DSSP-based start obtained from the slope-scaling iterate. All initializations are checked for feasibility, and duplicate supports are removed before launching IRLS. A-F2 Objective-Driven Perturbation Restarts Although the multi-start strategy improves robustness, each IRLS run still evolves locally from its initial support. Once the smoothing parameters become small, the surrogate strongly penalizes opening new arcs, and the method may continue improving the flow values without substantially changing the active arc set. To encourage larger support changes, we introduce an objective-driven perturbation restart layer. Let f denote an incumbent feasible flow returned by IRLS. We define the active and inactive arc sets as A(f)=e∈E:fe>τ,I(f)=e∈E:fe≤τ,A(f)=\e∈ E:f_e>τ\, I(f)=\e∈ E:f_e≤τ\, where τ>0τ>0 is the activation threshold used for scoring. Active arcs are ranked by their realized average cost σeact=re+bemaxfe,1,e∈A(f), _e^act=r_e+ b_e \f_e,1\, e∈ A(f), so arcs carrying little flow but incurring large fixed cost are treated as expensive. Inactive arcs are ranked by their capacity-amortized opening cost σeinact=re+bemaxce,1,e∈I(f), _e^inact=r_e+ b_e \c_e,1\, e∈ I(f), so arcs with low variable cost and low amortized fixed cost are treated as promising alternatives. Using these scores, we generate perturbed target vectors that modify the incumbent support in several complementary ways. Weak-arc deletions reduce flow on active arcs carrying very small amounts. Expensive-arc deletions reduce or remove active arcs with large σeact _e^act. Cheap-arc openings assign trial flow to inactive arcs with small σeinact _e^inact. Rerouting perturbations simultaneously suppress expensive active arcs and open cheap inactive arcs. Support swaps replace weak active arcs with a small set of cheap inactive alternatives. A perturbed target f~ f is not required to satisfy flow conservation or capacity constraints. For example, opening a cheap inactive arc does not guarantee that a complete feasible route exists through that arc. Therefore, the target is not accepted directly. Instead, it is projected onto the feasible capacitated-flow polytope ℱ=g∈ℝ|E|:Cg=d, 0≤g≤c.F=\g ^|E|:Cg=d,\;0≤ g≤ c\. The projected restart point is obtained from Πℱ(f~)=argming∈ℱ‖g−f~‖22. _F( f)= _g \|g- f\|_2^2. The restart initialization is then frestart(0)=Πℱ(f~)f_restart^(0)= _F( f). This projection allows the perturbation step to be aggressive at the support level while restoring feasibility through a controlled continuous optimization problem. If the projection fails, produces only a negligible change from the incumbent, or duplicates a support already generated by another perturbation, the candidate is discarded. The perturbed targets are first ranked by their realized objective value and support size. After projection and duplicate filtering, IRLS is rerun from the surviving restart candidates, and the best resulting flow is selected lexicographically by objective value, feasibility residual, and runtime. A-F3 Layer-3 Implementation Details This section details the anchor-union restricted search of Section VI: how anchor supports are generated, how candidate unions are formed and filtered, the restricted refit portfolio, restricted IRLS runs, runtime control, within-union refinement, and final candidate selection. Anchor generation. The anchor set ancX_anc is generated using ADSSP runs with four initializations: max_flow, unitflow, trans, multi_unit. Each anchor run is limited to 120120 iterations. Only feasible returned flows are retained as anchors. Union supports. For each restart-layer candidate f, the implementation constructs three types of restricted supports: S(f)S(f), S(f)∪S(a)S(f)∪ S(a) for each a∈anca _anc, and S(f)∪⋃a∈ancS(a)S(f)∪ _a _ancS(a). Duplicate and empty supports are removed. A support is also discarded if the restricted directed graph fails a reachability screening test from supply nodes to demand nodes. Restricted refit portfolio. For each retained support E′E , feasible warm starts are obtained from two sources: the corresponding Layer 2 candidate flow restricted to E′E , and restricted linear flow problems ming _g ∑e∈E′(re+γbemaxce,1)ge, _e∈ E (r_e+γ b_e \c_e,1\ )g_e, (47) s.t. .t. CE′g=d, C_E g=d, 0≤ge≤ce,e∈E′, 0≤ g_e≤ c_e, e∈ E , with γ∈0, 0.05, 0.25, 1, 4.γ∈\0,\;0.05,\;0.25,\;1,\;4\. The case γ=0γ=0 gives a pure variable-cost refit, while positive values penalize arcs with large fixed cost relative to capacity. These objectives are used only to generate feasible starts; all candidates, including the restricted Layer 2 flows, are scored by JtrueJ_true. If no feasible refit is found, the support is discarded. Among all feasible starts, duplicate numerical supports are removed and at most the three best starts are retained, ordered by JtrueJ_true and then by support size. Restricted IRLS runs. Restricted IRLS uses the same smoothing schedule, weighted-QP solver, support-refit rule, feasibility tolerance, and stopping criteria as the main IRLS method. The returned restricted flow is lifted to the original graph by assigning zero flow to arcs outside E′E . Support selection and runtime control. Candidate supports are ranked by the true objective value of their best feasible refit, with support size used as a secondary criterion. To control runtime, only a bounded number of supports are processed: Nunion=2,n≥4000orm≥200000,3,n≥1500orm≥75000,6,otherwise.N_union= cases2,&n≥ 4000\ or\ m≥ 200000,\\ 3,&n≥ 1500\ or\ m≥ 75000,\\ 6,&otherwise. cases The remaining Layer-3 time budget is divided among the selected supports. Within-union refinement. After restricted IRLS, the algorithm performs a monotone support refinement inside E′E . Active arcs are ranked for deletion using qedel=re+bemaxfe,10−6+12bemaxce,1.q_e^del=r_e+ b_e \f_e,10^-6\+ 12 b_e \c_e,1\. (48) Deletion batches of approximately 12%12\% and 6%6\% of the active support are tested first, followed by individual trials among the 2525 highest-ranked active arcs. At most five successful deletion rounds are performed. Unused arcs in E′E are ranked for reopening using qeadd=re+bemaxce,1.q_e^add=r_e+ b_e \c_e,1\. (49) The algorithm tests bundles of up to 88, 44, and 22 arcs, followed by individual trials among the 2020 highest-ranked unused arcs. At most three successful reopening rounds are performed. Each deletion or reopening trial is followed by the restricted refit portfolio in (47). A modification is accepted only if the refit is feasible and strictly decreases JtrueJ_true. Thus, this refinement phase is monotone with respect to the reported fixed-charge objective. Final selection. Layer 3 returns the feasible candidate with the smallest JtrueJ_true among the restart-layer candidates, restricted refits, restricted IRLS outputs, and within-union refinement outputs. A-G Experiments This section gives the hardware and software environment, random seeds, and instance-generation parameters underlying the results in Section VII. Hardware and Software. All experiments were run on a MacBook Pro (Apple M4 Pro, 48 GB RAM) under macOS Sequoia 15.7.5. Methods were implemented in Python 3.14.3. Min-cost-flow and MILP subproblems use OR-Tools 9.15.6755. The MILP baseline is formulated through the OR-Tools pywraplp interface and solved with its SAT backend. For consistency across methods, reported objective values are obtained by re-evaluating the returned flow under the original fixed-charge objective. Reported times are wall-clock, single-run, on an otherwise idle machine. Random seeds. All reported heuristic runs use a master seed of 00 (default). Per-initialization seeds within a run are derived deterministically from the master seed (e.g. offsets of the form seed+idxseed+idx for multi-start IRLS and ADSSP initializations), so a fixed master seed reproduces the full 410-instance run exactly. Instance generation. The synthetic family is generated by a parametric single-commodity generator with per-family specifications of node count, arc density, supply/demand node counts, capacity factor, and fixed-cost ratio, drawn with numpy’s default_rng(seed) for reproducibility. The ten synthetic subfamilies vary these parameters (e.g. high_fixed_cost, low_fixed_cost, tight_capacity) to probe different regimes of the fixed-charge-to-variable-cost ratio and capacity tightness, with 1010 instances drawn per subfamily. The scale_* families follow the generation procedure of [23]. FCNetLib_fc and Benchmark are fixed, externally sourced instance sets and are not randomly generated. A-H Additional Experimental Results This section supports Section VII with the full per-family breakdown, the split by MILP termination status, and the anchor-only ablation referenced there. TABLE V: Mean gap (%) to the MILP baseline and mean time (s) by family. The best heuristic mean gap per family is bold. A negative entry indicates L3 improves on a time-limited MILP incumbent. Mean gap (%) Time (s) Family # instances L3 ADSSP DSSP L2 L1 MCF MILP L3 L2 L1 ADSSP DSSP MCF FCNetLib_fc 20 6.972 12.881 12.947 12.227 17.312 39.067 2.33 25.53 19.11 5.10 0.03 0.04 0.01 Benchmark 20 7.642 9.462 10.270 11.934 13.363 21.392 155.96 10.15 5.54 3.74 0.15 0.63 0.03 synthetic 100 1.706 4.766 6.089 2.688 4.359 10.979 75.85 7.30 3.17 0.89 0.15 0.78 0.04 sc.20_100 30 0.014 0.378 0.550 0.132 0.406 0.627 0.05 0.30 0.12 0.04 0.01 0.01 0.00 sc.60_400 30 0.348 1.486 1.661 0.776 1.262 3.678 1.97 1.40 0.45 0.14 0.03 0.03 0.01 sc.120_1500 30 0.555 1.462 2.696 1.049 1.785 4.546 109.24 4.20 1.63 0.43 0.11 0.20 0.03 sc.150_2500 30 0.809 1.816 2.850 1.291 2.050 5.148 173.13 6.28 1.99 0.52 0.18 0.41 0.05 sc.200_4000 30 0.487 1.285 2.548 0.939 1.600 4.245 209.38 10.47 3.33 0.89 0.32 0.89 0.08 sc.250_7500 30 0.300 0.827 2.591 0.938 1.421 4.375 244.92 17.69 6.87 1.80 0.53 2.08 0.14 sc.300_9000 30 0.129 0.477 2.090 0.609 0.999 3.667 257.96 22.04 8.62 2.31 0.65 2.69 0.17 sc.500_10000 30 0.199 1.020 2.613 0.778 1.325 4.949 324.28 34.78 10.02 2.50 0.84 2.82 0.19 sc.1000_20000 30 -0.280 0.519 2.446 0.317 0.912 5.356 301.56 99.92 21.96 5.84 1.78 6.82 0.39 TABLE VI: Mean/median gap to the MILP baseline (%) and mean time (s) by MILP termination status. optimal (n=235n=235) time-limited incumbent (n=175n=175) Method mean med. time mean med. time MILP 0.00 0.000 22.76 0.00 0.000 309.00 L3 1.48 0.014 6.76 1.10 0.329 32.96 L2 2.46 0.139 3.42 2.17 1.069 9.47 L1 3.69 0.613 1.04 3.06 1.683 2.61 ADSSP 3.50 0.386 0.12 2.17 0.883 0.71 DSSP 3.89 0.694 0.24 4.35 3.347 2.94 MCF 8.72 2.332 0.03 7.75 6.908 0.17 TABLE VII: Layer-ablation performance over all 410410 instances. Gaps are relative to the MILP baseline on each instance; win/tie rates are relative to the best displayed heuristic. Method Mean gap (%) Median gap (%) Win/tie rate (%) Mean time (s) L3 1.316 0.131 87.6 18.001 L3- 2.073 0.560 29.0 6.523 L2 2.334 0.749 23.7 6.052 L1 3.420 1.412 14.4 1.781 ADSSP anchors 2.113 0.384 30.7 0.620 ADSSP 2.891 0.695 18.3 0.346 L3 attains the best objective quality consistently across dataset families, including the larger-scale instances where the MILP baseline often terminates without proving optimality. A family-wise breakdown is given in Table V. The MILP status split in Table VI further supports the conclusion that the advantage of L3 is not limited to instances where the MILP solves easily. L3 remains the strongest heuristic both on instances where the MILP proves optimality and on instances where the MILP terminates with only a feasible incumbent. On the 235235 instances where the MILP proves optimality, L3 has mean and median gaps of 1.48%1.48\% and 0.014%0.014\%, respectively. On the 175175 instances where the MILP terminates with a time-limited incumbent, L3 has mean and median gaps of 1.10%1.10\% and 0.329%0.329\%. Thus, the proposed method performs well both when the MILP reference is certified and when exact search is unable to close the gap within the time limit. Table VII isolates the anchor-union contribution with two ablations. Note that the win/tie rates in this table are computed over a different method pool than Table I of the main paper (L3, L3-, L2, L1, ADSSP anchors, and ADSSP here, versus L3, L2, L1, ADSSP, DSSP, and MCF there), so L3’s 87.6%87.6\% rate here is not directly comparable to the 90.0%90.0\% reported in the main paper. Running ADSSP restricted to the same anchor supports that L3 uses reaches a mean gap of 2.113%2.113\% with a 30.7%30.7\% win/tie rate, well short of L3, showing that L3 is not simply returning a good ADSSP solution. L3-, which corresponds to L3 with the anchor set restricted to ∅ , i.e., L3 without using any ADSSP anchors, but keeping the union search, restricted refits, and add/drop refinement, reaches a mean gap of 2.073%2.073\%, a 0.7570.757-percentage-point degradation from full L3. Since +A^+ in Algorithm 2 (main paper) already contains ∅ , L3 is never worse than this restricted variant by construction, and is strictly better on 278/410278/410 instances. Together, the two ablations show that neither the ADSSP anchor supports nor the IRLS-based union search alone accounts for L3’s full gain: the two mechanisms are complementary and contribute comparably. Notably, L3- already outperforms the standard ADSSP baseline itself, winning or tying head-to-head on 280/410280/410 instances and reducing the mean gap by 0.8180.818 percentage points, so L3’s advantage does not depend on borrowing ADSSP’s solution quality even indirectly.