Path Finding II : An \~O(m sqrt(n)) Algorithm for the Minimum Cost Flow Problem

1.1 Previous Work

• Previous work improved special cases and related settings, including undirected flow, exact directed flow, and unit-capacity graphs, but general directed maximum flow retained the Goldberg–Rao bound [10].Mądry [25] achieved eO(m^10/7) for uncapacitated directed graphs, while Even and Tarjan and Karzanov established the earlier unit-capacity barrier.Source: §1 · p. 1, §1.1 · pp. 2–3

1.4 Paper Organization

• The paper hides polylogarithmic factors in eO(·), restricts attention to weakly polynomial algorithms, and does not attempt to reduce those factors.Its discussion also notes that standard barrier analysis faces an Ω(m) self-concordance lower bound for parallel-edge instances.Source: §1 · p. 1, §1.1 · p. 2, §1.2 · p. 4

1 Introduction

• The algorithm solves directed maximum flow in eO(m√n logO(1)(U)) time, improving the best general bound unchanged since Goldberg and Rao [10].It also gives eO(m√n logO(1)(U)) work and eO(√n logO(1)(U)) depth using parallel Laplacian solvers.Source: §1 · pp. 1–2
• The method extends to minimum-cost flow and lossy generalized minimum-cost flow, with eO(m√n logO(1)(U/ϵ)) time for ϵ-approximate generalized solutions.The standard minimum-cost maximum-flow algorithm has matching eO(m√n logO(1)(U)) work and eO(√n logO(1)(U)) depth bounds.Source: §1 · p. 2, §1.3 · p. 5

1.2 Our Approach

• The paper’s central novelty is a general interior-point method whose convergence rate beats that predicted by the best self-concordant barrier for the feasible region.The analysis uses weighted path following and measures Newton steps with a mixed norm combining ℓ∞ and weighted ℓ2 behavior.Source: §1 · p. 2, §1.2 · p. 4, §1.3 · pp. 4–5

1 Introduction

• The framework reuses sampling, fast Laplacian solvers, and Daitch–Spielman reductions to make interior-point iterations inexpensive for flow problems.For maximum flow, the operational interpretation alternates cost reweighting, electric-flow computation, and effective-resistance estimation.Source: §1 · p. 2, §1.1 · p. 3

1.3 Our Contributions

• For flow linear programs, the method directly handles edge-flow bounds through coordinate barriers and weights, avoiding the rank increase caused by naive reformulations.The resulting iterations depend on rank(A), while each iteration requires only eO(1) linear-system solves of the form A^TDAx=d.Source: §1.2 · p. 4, §1.3 · p. 5

• Vectors use bold symbols, with unbolded symbols denoting their entries; [z] denotes the integers from 1 to z.Source: §2 · p. 6
• Multiplicative closeness is expressed through coordinate-wise bounds, matrix inequalities, or approximately equivalent log-scale distances for small errors.Source: §2 · p. 6
• Projection matrices, leverage scores, leverage matrices, and projection Laplacians are defined from a reweighted matrix A and its associated projection.Source: §2 · p. 6
• Convex sets, symmetric sets, scaled sets, and ℓ_p balls are defined geometrically through convex combinations, negation, scaling, and norm bounds.Source: §2 · pp. 6–7
• Gradients, Hessians, Jacobians, coordinate-wise vector operations, positive-definite norms, matrix order, Schur products, and diagonal matrices receive standard formal definitions.Source: §2 · pp. 6–7

3.1 The Problem

• 3.1 The Problem: The central task is minimizing a linear objective subject to A^T x = b and coordinate bounds l_i ≤ x_i ≤ u_i.The formulation assumes nondegenerate domains, full column rank, and a nonempty interior.Source: §3.1 · p. 7
• 3.1 The Problem: The formulation is a direct generalization of standard form, avoiding transformations that could increase rank(A).Source: §3.1 · p. 7

3.2 Coordinate Barrier Functions

• 3.2 Coordinate Barrier Functions: Each coordinate domain is assigned a 1-self-concordant barrier, with explicit logarithmic barriers for one-sided bounds.Source: §3.2 · pp. 7–8
• 3.2 Coordinate Barrier Functions: For finite lower and upper bounds, the paper uses a trigonometric barrier that diverges at both endpoints and is derived from a third-order differential equation.The authors state that they are unaware of this barrier being used previously.Source: §3.2 · p. 8
• 3.2 Coordinate Barrier Functions: Self-concordance controls how quickly the barrier Hessian changes and how much force the barrier can exert.Source: §3.2 · p. 8

4.1 The Weighted Central Path

• 4.1 The Weighted Central Path: For fixed positive weights, minimizers of the penalized objective as t increases form a weighted central path whose limit solves the linear program.Source: §4.1 · p. 9
• 4.1 The Weighted Central Path: The algorithm alternates between increasing t, updating x toward the path, and updating weights to select a better path.Source: §4.1 · p. 9

4.2 Measuring Centrality.

• 4.2 Measuring Centrality.: Centrality measures the distance from the current feasible point to the minimizer of the penalized objective, using a mixed norm to track relevant quantities.Source: §4.2 · p. 10
• 4.2 Measuring Centrality.: The projected Newton step is defined using barrier gradients, Hessians, and projections onto the constraint-kernel space.Source: §4.2 · pp. 10–11

4.3 The Weight Function

• 4.3 The Weight Function: The iterative scheme increases t when centrality is small, then applies a projected Newton step and moves the weights toward the weight function.Source: §4.3 · p. 11
• 4.3 The Weight Function: Weight functions are differentiable maps satisfying slack-sensitivity, step-consistency, and uniformity conditions.The paper states that such a weight function exists and uses it to control changes in centrality.Source: §4.3 · p. 11

5 Progressing Along Weighted Paths

• 5 Progressing Along Weighted Paths: The section develops an iteration scheme with approximately O(rank(A) log(U/ε)) weighted path-following iterations.Source: §5 · p. 12

5.1 Changing t

• 5.1 Changing t: Increasing t raises centrality at a rate governed by the norm parameter C_norm and the current weights.Source: §5.1 · p. 12

5.2 Changing⃗x

• 5.2 Changing⃗x: Under sufficient centrality and multiplicative weight accuracy, a projected Newton step decreases centrality quadratically: δ_t(x_1,w) ≤ 4(δ_t(x_0,w))^2.Source: §5.2 · pp. 12–13

5.3 Changing⃗w

• 5.3 Changing⃗w: The analysis bounds movement of the weight function as x changes and bounds the centrality cost of moving weights toward it.Source: §5.3 · pp. 14–15

5.4 Centering

• 5.4 Centering: Approximate weight computation is handled through a chasing-0 game that keeps weight error small while preserving progress and centrality invariants.Source: §5.4 · p. 16
• 5.4 Centering: The centering procedure preserves bounded centrality and maintains the stated potential bound after updating the point and weights.Source: §5.4 · pp. 17–18

6 Weight Function

• The weight function is a central component of the algorithm, extending the approach from Part I [22] to obtain the paper’s improved running times.The section develops its properties and convergence guarantees for the broader barrier-function framework.Source: §6 · p. 19
• The analysis chooses α ∈[1, 2), β ∈(0, 1), and β^(1−α) ≤2, with β small and α just slightly larger than 1.9.Source: §6 · p. 19
• The weight function is defined through the unique minimizer of a strongly convex function, with weights constrained between β and 1 + β.The Hessian is positive definite in the relevant region, ensuring uniqueness and well-definedness.Source: §6 · pp. 19–20
• Its derivative obeys G′(x) = −G(x)(G(x) + αΛ)^−1Λ, enabling the analysis of how weights change along the path.The section derives this relation using the defining optimality condition and the implicit function theorem.Source: §6 · pp. 19–20

6.1 Computing and Correcting Weight Function

• The computeWeight procedure obtains an approximate weight using gradient descent and approximate leverage scores, tolerating the resulting multiplicative noise.Leverage scores can be computed with high probability using Õ(1) linear-system solves.Source: §6.1 · pp. 22–23
• Theorem 17 gives weight computation and correction using Õ(log^3(1/K)/K^2) linear-system solves, while initialization achieves the same guarantee with constant probability.The stated initialization bound includes a rank(A) factor and logarithmic dependence on m/K.Source: §6.1 · p. 23

7 The Algorithm

• Theorem 18 states that LPSolve returns x with c^T x ≤ OPT + ϵ in the theorem’s stated running time.The algorithm starts from an interior point and alternates path-following and inexact centering steps.Source: §7 · pp. 23, 25
• The convergence proof relates progress along the weighted central path to optimality through Lemmas 19 and 20.The proof uses path centrality to control distance from the path and then bounds the resulting objective error.Source: §7 · pp. 23–24
• PathFollowing maintains bounded centrality and weight-error invariants while decreasing δ_t until δ_t ≤ ϵ.Its running time is expressed in terms of the cost T of solving the required linear system.Source: §7 · p. 25
• The stability analysis shows that approximate weight computations can be incorporated into centering while preserving the algorithm’s invariants.The first iteration requires O(log(U/ϵ)) calls to S_x,w, but the amortized number of calls is constant per iteration.Source: §7 · p. 25
• The algorithm uses an initial path-following phase followed by a second phase that drives the parameter to t_2 and achieves the target approximation.The implementation selects t_1, t_2, ϵ_1, and ϵ_2 as functions of m, U, and ϵ.Source: §7 · p. 26

8 Linear System Solver Requirements

• Approximate solvers with error ϵS = 1/m^d suffice for the algorithms, addressing accuracy and equality-constraint challenges from inexact linear algebra.The section’s solver implementation is based on showing that this accuracy supports Sections 8.1–8.4, while Section 8.5 establishes stability properties for efficient construction.Source: §8 · p. 27

8.1 The normal force A⃗η.

• The maintained normal force removes the potentially large equality-constraint component from Newton systems, enabling stable approximate steps.The normal force can be computed explicitly and maintained using Sx,w with an amortized constant number of calls per iteration, although the first iteration requires O(log(U/ϵ)) calls.Source: §8.1 · pp. 27–28

8.2 An efficient⃗x step

• The efficient approximate x-step preserves progress while controlling infeasibility, with W ≤ 3 and periodic correction reducing infeasibility more than it harms centrality.The approach relies on ϵS < C/m^2 and cannot enforce ATx(apx) = b exactly, so infeasibility is measured and corrected explicitly.Source: §8.2 · pp. 31–34

8.3 An efficient⃗w step

• Approximate weights require no linear solve for centering, while leverage-score estimation can use approximate solvers as described in Part I [22, ArXiv v3, Section D].The weight-changing step uses approximate weights, and only the computation of those weights requires approximate linear-system solves.Source: §8.3 · p. 34

8.4 The stable algorithm

• The stable algorithm implements LPSolve with approximate linear algebra while maintaining objective accuracy, approximate feasibility, and bounded variables; it also derives dual solutions from the maintained normal force.When infeasibility exceeds 1/m^2, a feasibility-improvement step reduces it without significantly harming progress, and the dual construction uses an additional ˜O(log(mU/ϵ)) iterations.Source: §8.4 · pp. 34–35

8.5 Well conditioned

• The algorithm’s systems are well conditioned: Lemma 32 bounds Dk between polynomially scaled copies of A^T A, enabling fast matrix solvers.The conditioning bound is poly(mU/ϵ), and the preceding lemmas also bound weighted-central-path slacks and changes in successive systems.Source: §8.5 · pp. 37–38
• The linear system matrix has the form A^T D A with D = W^-1Φ''(x)^-1.The weights satisfy the stated bounds, while Theorem 18 and Lemma 31 provide bounds on Φ''(x) and the slacks.Source: §8.5 · p. 38
• The condition number of D is bounded by poly(mU/ϵ).This follows from lower and upper bounds on Φ''(x), including the slack lower bound.Source: §8.5 · p. 38
• Changes in D are bounded by separately controlling changes in Φ''(x) and W.For W, centeringInexact ensures that changes in log(w) lie in (1+ϵ)U, with U defined by a weighted infinity-norm constraint.Source: §8.5 · p. 39

9 Generalized Minimum Cost Flow

• The generalized-flow result is an application of the same path-following algorithm used for maximum flow and minimum cost flow.The paper presents the generalized formulation directly because it subsumes the simpler cases.Source: §9 · pp. 39, 41
• The generalized formulation allows edge multipliers γ_e ≤ 1, with maximum flow and minimum cost flow recovered when γ_e = 1.For generalized flow, the paper targets approximation because efficient exact algorithms are difficult.Source: §9 · pp. 39–40
• The reduction constructs a graph with O(m) edges and O(n) vertices, together with an interior point and a well-conditioned linear program.The transformed system can be solved in nearly linear time under the stated diagonal-matrix condition-number bound.Source: §9 · p. 40
• The solver supports explicit lower and upper bounds on flow variables, avoiding the need to split each flow variable into positive and negative parts.The primal path-following formulation also adds bounded auxiliary variables without changing the optimum.Source: §9 · p. 41
• The algorithm computes ϵ-approximate generalized minimum cost maximum flow in Õ(m√n log^O(1)(U/ϵ)) work and Õ(√n log^O(1)(U/ϵ)) depth.It also computes exact standard minimum cost maximum flow with Õ(m√n log^O(1)(U)) work and Õ(√n log^O(1)(U)) depth.Source: §9 · p. 41
• The reduction combines Theorem 18 with nearly linear-work Laplacian solvers after converting each required system into Õ(log(U/ϵ)) Laplacian systems.This yields the parallel work-depth guarantees stated in Theorem 35.Source: §9 · p. 41

• The linear-program notation uses A, c, b, x, and w for matrices, costs, constraints, solutions, and constraint weights, with m variables and n constraints.Source: A Glossary · p. 44
• Matrix notation represents vector quantities with diagonal matrices S, W, and Φ corresponding to s, w, and φ.Source: A Glossary · p. 44

Barrier functions

• The paper specifies logarithmic barrier functions for lower-bounded, upper-bounded, and interval-constrained domains.Source: A Glossary · p. 44

Weight function

• The glossary records size and step-consistency properties for the weight function, with deviations controlled by polylogarithmic factors in m.Source: A Glossary · p. 44

Potential function

• The tracing potential is Φ_µ(x) = e^µx + e^-µx and is maintained at approximately polynomial magnitude in m.Source: A Glossary · p. 44

• The lemmas establish a logarithmic approximation when |log(a) − log(b)| is at most 1/2.Source: §B.1 · p. 45
• A projection-matrix lemma relates the diagonal of a projection matrix to vectors and weights.Source: §B.1 · p. 45
• The section states a general inequality for x, epsilon, and positive lambda, with the displayed bound involving (x + epsilon)tan^-1(lambda(x + epsilon)).Source: §B.1 · pp. 45–46

B.2 Projection on Mixed Norm Ball

• A change of variables extends the projection algorithm from the ball formulation to the related problem in (B.2).Source: §B.2 · pp. 46, 48
• The analysis assumes normalized inputs and descending coordinate ratios, yielding a structured solution form based on signed boundary coordinates.Source: §B.2 · p. 46
• The derivation uses concavity to obtain a closed-form solution on each interval and then selects the interval with the best value.Source: §B.2 · p. 47
• The algorithm normalizes the input, sorts coordinates by descending |a_i|/l_i, precomputes prefix quantities, optimizes over intervals, and outputs the best candidate.Source: §B.2 · pp. 46–47
• The exposition notes that some boundary cases are omitted for simplicity.Source: §B.2 · p. 47
• Theorem 40 gives a solution to the mixed-norm projection problem in total work O~(m) and depth O~(1).Source: §B.2 · p. 48