Flows in Almost Linear Time via Adaptive Preconditioning

TL;DR

This paper introduces nearly linear time algorithms for solving a broad class of flow and regression problems on graphs, extending techniques from linear systems to convex objectives like _p-norm flows and total variation minimization.

Contribution

It develops adaptive non-linear preconditioning and graph decomposition methods, enabling almost-linear time solutions for complex convex optimization problems on graphs.

Findings

Achieved almost-linear time algorithms for _p-norm flow and semi-supervised learning.

Provided new approximation methods for max-flow and total variation minimization.

Extended linear system techniques to a wider class of convex objectives.

Abstract

We present algorithms for solving a large class of flow and regression problems on unit weighted graphs to $(1 + 1 / poly(n))$ accuracy in almost-linear time. These problems include $\ell_p$-norm minimizing flow for $p$ large ($p \in [\omega(1), o(\log^{2/3} n) ]$), and their duals, $\ell_p$-norm semi-supervised learning for $p$ close to $1$. As $p$ tends to infinity, $\ell_p$-norm flow and its dual tend to max-flow and min-cut respectively. Using this connection and our algorithms, we give an alternate approach for approximating undirected max-flow, and the first almost-linear time approximations of discretizations of total variation minimization objectives. This algorithm demonstrates that many tools previous viewed as limited to linear systems are in fact applicable to a much wider range of convex objectives. It is based on the the routing-based solver for Laplacian linear systems by Spielman and Teng (STOC '04, SIMAX '14), but require several new tools: adaptive non-linear preconditioning, tree-routing based ultra-sparsification for mixed $\ell_2$ and $\ell_p$ norm objectives, and decomposing graphs into uniform expanders.

Tables

Table 1. Table 1: Glossary of Notations in Algorithm and Analysis.

Notations in Ultrasparsify
$𝒢$	Input instance with $(V^{𝒢},E^{𝒢},{𝒈}^{𝒢},{𝒓}^{𝒢},s^{𝒢})$.
$T$	Low stretch spanning tree of $𝒢$ (stretch with respect to ${𝒓}^{𝒢}$).
$E^r$	All in $E^{𝒢}\setminus T$ whose resistance after rounding is $r$.
$E_i^r$	The remaining edges in $E^r$ at the $i$-th iteration of tree-portal routing $E^r$.
${\hat{G}}_i^r$	The image of edges in $E_i^r$ by the mapping ${\text{Move}}_{T,\hat{V}}$, i.e. moving off-tree edges along tree-portal path. The gradients of edges in ${\hat{G}}_i^r$ are set to preserve the linear flow dot gradient term under tree-portal routing.
$m_{r,i}$	The number of edges in $E_i^r$ (also the size of ${\hat{G}}_i^r$).
${\hat{G}}_{i,j}^r$	The $j$-th expander we get from decomposing ${\hat{G}}_i^r$. Edges keep their gradients from ${\hat{G}}_i^r$, and $𝒓,s$ are scaled.
${\hat{E}}_i^r$	The union of edges contained in the expanders ${\hat{G}}_{i,j}^r$ (i.e. over all $j$’s).
${\tilde{E}}_i^r$	The pre-image of edges in ${\hat{E}}_i^r$.
${\text{E}}_{last}^r$	The set of edges remaining in $E^r$ after the last iteration for $r$.
${\tau }_{r,i}$	The probability we use in SampleAndFixGradient to uniformly sample ${\hat{G}}_{i,j}^r$.
${\overline{H}}_{i,j}^r$	The sparsified graph of ${\hat{G}}_{i,j}^r$ computed by SampleAndFixGradient. $𝒈,𝒓,s$ on edges are computed by the subroutine.
Additional notations in the analysis
$\overline{G}$	The instance with the same edge set as $𝒢$. Note$†$  $E^{𝒢}=T+{\sum }_r{E}_{last}^r+{\sum }_{r,i}{\tilde{E}}_i^r$. Edges in $\overline{G}$ has the same $𝒈,𝒓,s$ as in $𝒢$ except for those in ${\sum }_{r,i}{\tilde{E}}_i^r$. For any resistance $r$ and iteration $i$, $e\in {\tilde{E}}_i^r$ has the same gradient as in $𝒢$, but ${𝒓}_e=r\kappa {\log }^{2}n$, $s_e={\tau }_{r,i}^{-p}{s}^{𝒢}$ are scaled.
${\overline{G}}_i^r$	The instance $\overline{G}$ restricted to the set of edges in ${\tilde{E}}_i^r$.
${\overline{G}}_{rest}$	The instance $\overline{G}$ restricted to the set of edges in ${\sum }_r{E}_{last}^r$.

Equations

$$\min_{\text{flow $\boldsymbol{\mathit{f}}$ with residue $\boldsymbol{\mathit{b}}$}}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|_{\odot}.$$

$$\displaystyle\min_{\text{flow $\boldsymbol{\mathit{f}}$ with residue $\boldsymbol{\mathit{b}}$}}\sum_{e}\boldsymbol{\mathit{g}}_{e}\boldsymbol{\mathit{f}}_{e}+\boldsymbol{\mathit{r}}_{e}\boldsymbol{\mathit{f}}_{e}^{2}+\mathopen{}\mathclose{{}\left|\boldsymbol{\mathit{f}}_{e}}\right|^{p},$$

$$\displaystyle\min_{\boldsymbol{\mathit{B}}^{\top}\boldsymbol{\mathit{f}}=\boldsymbol{\mathit{b}}}\sum_{e}\boldsymbol{\mathit{g}}_{e}\boldsymbol{\mathit{f}}_{e}+\boldsymbol{\mathit{r}}_{e}\boldsymbol{\mathit{f}}_{e}^{2}+\mathopen{}\mathclose{{}\left|\boldsymbol{\mathit{f}}_{e}}\right|^{p},$$

$$\text{val}(\boldsymbol{\widetilde{\mathit{f}}})-\text{OPT}\leq\frac{1}{{\textrm{poly}}(m)}\mathopen{}\mathclose{{}\left(\text{val}(\boldsymbol{\mathit{f}}^{(0)})-\text{OPT}}\right)+\frac{1}{{\textrm{poly}}(m)}$$

$$\mathopen{}\mathclose{{}\left\|\boldsymbol{\widetilde{f}}}\right\|^{p}_{p}\leq\mathopen{}\mathclose{{}\left(1+\frac{1}{{\textrm{poly}}(m)}}\right)\min_{\boldsymbol{\mathit{f}}:\boldsymbol{\mathit{B}}^{G\top}\boldsymbol{\mathit{f}}=\boldsymbol{\mathit{b}}}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|^{p}_{p}.$$

$$\min_{\boldsymbol{\mathit{x}}\in\Re^{V}\mid\boldsymbol{\mathit{x}}_{T}=\boldsymbol{\mathit{s}}_{T}}\sum_{u\sim v}\mathopen{}\mathclose{{}\left|\boldsymbol{\mathit{x}}_{u}-\boldsymbol{\mathit{x}}_{v}}\right|^{p}.$$

$$\sum_{e}\boldsymbol{\mathit{r}}_{e}\boldsymbol{\mathit{f}}_{e}^{2}+\frac{\epsilon\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{r}}}\right\|_{1}}{m}\mathopen{}\mathclose{{}\left|\boldsymbol{\mathit{f}}_{e}}\right|^{p}.$$

$$\sum_{e}\boldsymbol{\mathit{r}}_{e}\boldsymbol{\mathit{f}}_{e}^{2}+\frac{\epsilon\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{r}}}\right\|_{1}}{m}\mathopen{}\mathclose{{}\left|\boldsymbol{\mathit{f}}_{e}}\right|^{p}\leq\mathopen{}\mathclose{{}\left(1+\epsilon}\right)\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{r}}}\right\|_{1}.$$

$$\sum_{e}\boldsymbol{\mathit{r}}_{e}\boldsymbol{\mathit{f}}_{e}^{2}+\frac{\epsilon\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{r}}}\right\|_{1}}{m}\mathopen{}\mathclose{{}\left|\boldsymbol{\mathit{f}}_{e}}\right|^{p}\leq\mathopen{}\mathclose{{}\left(1+1.1\epsilon}\right)\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{r}}}\right\|_{1}.$$

$$\frac{\epsilon\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{r}}}\right\|_{1}}{m}\mathopen{}\mathclose{{}\left|\boldsymbol{\mathit{f}}_{e}}\right|^{p}\leq\mathopen{}\mathclose{{}\left(1+1.1\epsilon}\right)\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{r}}}\right\|_{1}\leq 2\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{r}}}\right\|_{1},$$

$$\mathopen{}\mathclose{{}\left|\boldsymbol{\mathit{f}}_{e}}\right|\leq\mathopen{}\mathclose{{}\left(m/\epsilon}\right)^{1/p}\leq m^{o\mathopen{}\mathclose{{}\left(1}\right)}$$

$$\sum_{u}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{x}}_{u}-\boldsymbol{\mathit{s}}_{u}}\right)^{2}+\sum_{i}\sqrt{\sum_{e\in E_{i}}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{x}}_{u}-\boldsymbol{\mathit{x}}_{v}}\right)^{2}}.$$

$$\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}_{E_{i}}}\right\|_{2}^{2}\leq 1.$$

$$h_{p}(r,s,x)=rx^{2}+s\mathopen{}\mathclose{{}\left|x}\right|^{p}.$$

$$h_{p}(\boldsymbol{\mathit{r}},s,\boldsymbol{\mathit{x}})=\sum_{i=1}^{m}h_{p}(\boldsymbol{\mathit{r}}_{i},s,\boldsymbol{\mathit{x}}_{i})=\sum_{i=1}^{m}(\boldsymbol{\mathit{r}}_{i}\boldsymbol{\mathit{x}}_{i}^{2}+s\mathopen{}\mathclose{{}\left|\boldsymbol{\mathit{x}}_{i}}\right|^{p}).$$

$$\mathcal{G}\stackrel{\mathrm{def}}{=}(V^{\mathcal{G}},E^{\mathcal{G}},g^{\mathcal{G}},r^{\mathcal{G}},s^{\mathcal{G}}),$$

$$\mathcal{E}^{\mathcal{G}}(\boldsymbol{\mathit{f}})=\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{g}}^{\mathcal{G}}}\right)^{\top}\boldsymbol{\mathit{f}}-h_{p}(\boldsymbol{\mathit{r}},s,\boldsymbol{\mathit{f}}).$$

$$\max_{\boldsymbol{\mathit{f}}:(\boldsymbol{\mathit{B}}^{\mathcal{G}})^{\top}\boldsymbol{\mathit{f}}=\boldsymbol{\mathit{b}}}\mathcal{E}^{\mathcal{G}}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{f}}}\right).$$

$$\frac{1}{\kappa}\mathcal{E}^{\mathcal{H}}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{f}}^{\mathcal{H}}}\right)\leq\mathcal{E}^{\mathcal{G}}\mathopen{}\mathclose{{}\left(\frac{1}{\kappa}\boldsymbol{\mathit{f}}^{\mathcal{G}}}\right).$$

$$\mathcal{G}=(V^{{\mathcal{G}}_1},E^{{\mathcal{G}}_1}\cup E^{{\mathcal{G}}_2},(g^{{\mathcal{G}}_1},g^{{\mathcal{G}}_2}),(r^{{\mathcal{G}}_1},r^{{\mathcal{G}}_2}),(s^{{\mathcal{G}}_1},s^{{\mathcal{G}}_2})).$$

$$\mathcal{M}_{\mathcal{H}\rightarrow\mathcal{G}}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{f}}^{\mathcal{H}}=(\boldsymbol{\mathit{f}}^{\mathcal{H}_{1}},\boldsymbol{\mathit{f}}^{\mathcal{H}_{2}})}\right)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathopen{}\mathclose{{}\left(\mathcal{M}_{\mathcal{H}_{1}\rightarrow\mathcal{G}_{1}}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{f}}^{\mathcal{H}_{1}}}\right),\mathcal{M}_{\mathcal{H}_{2}\rightarrow\mathcal{G}_{2}}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{f}}^{\mathcal{H}_{2}}}\right)}\right),$$

$$g^{\mathcal{G}}={\hat{g}}^{\mathcal{G}}+B^{\mathcal{G}}{\psi }^{\mathcal{G}},\text{ s.t. }{B^{\mathcal{G}}}^{\top }{\hat{g}}^{\mathcal{G}}=0.$$

$$\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}}\right\|_{2}^{2}=\min_{\boldsymbol{\mathit{x}}}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{g}}+\boldsymbol{\mathit{B}}\boldsymbol{\mathit{x}}}\right\|_{2}^{2}.$$

$$\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}}\right\|_{2}^{2}\leq\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}}\right\|_{2}^{2}\leq\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{g}}^{\mathcal{G}}}\right\|_{2}^{2}.$$

$$\displaystyle\min_{\boldsymbol{\mathit{B}}\boldsymbol{\mathit{f}}=\boldsymbol{\mathit{b}}}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|_{p}^{p}$$

$$f\leftarrow f+\mathbf{\delta },$$

$$2^{-p}\cdot h_p(r+|x{|}^{p-2},s,\mathbf{\delta })\le h_p(r,s,x+\mathbf{\delta })-h_p(r,s,x)-{\mathbf{\delta }}^{\top }{\nabla }_x{h}_p(r,s,x)\le 2^{2p}\cdot h_p(r+|x{|}^{p-2},s,\mathbf{\delta }).$$

$$\displaystyle\begin{array}[]{rl}\max_{\boldsymbol{\mathit{x}}}&\mathcal{E}_{1}(\boldsymbol{\mathit{x}})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\boldsymbol{\mathit{g}}^{\top}\boldsymbol{\mathit{x}}-h_{p}(\boldsymbol{\mathit{r}},s,\boldsymbol{\mathit{x}})\\ \text{s.t.}&\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{b}}\end{array}$$

$$\displaystyle\begin{array}[]{rl}\max_{\boldsymbol{\delta}}&\mathcal{E}_{2}(\boldsymbol{\delta})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}(\boldsymbol{\mathit{g}}^{\prime})^{\top}\boldsymbol{\delta}-h_{p}(\boldsymbol{\mathit{r}}^{\prime},s,\boldsymbol{\delta})\\ \text{s.t.}&\boldsymbol{\mathit{A}}\boldsymbol{\delta}=\boldsymbol{\mathbf{0}},\end{array}$$

Full text

As $p$ tends to infinity, $\ell_{p}$-norm flow and its dual tend to max-flow and min-cut respectively. Using this connection and our algorithms, we give an alternate approach for approximating undirected max-flow, and the first almost-linear time approximations of discretizations of total variation minimization objectives.

This algorithm demonstrates that many tools previous viewed as limited to linear systems are in fact applicable to a much wider range of convex objectives. It is based on the the routing-based solver for Laplacian linear systems by Spielman and Teng (STOC ’04, SIMAX ’14), but require several new tools: adaptive non-linear preconditioning, tree-routing based ultra-sparsification for mixed $\ell_{2}$ and $\ell_{p}$ norm objectives, and decomposing graphs into uniform expanders.

1. 1 Introduction

2. 1.1 Main Results and Applications

3. 1.1.1 $\ell_{p}$-Norm Flows

4. 1.1.2 Semi-Supervised Learning on Graphs.

5. 1.1.3 Use as Oracle in Conjunction with Multiplicative Weight Updates

6. 1.2 Related Work

7. 1.3 Overview

8. 1.4 Open Questions

9. 2 Preliminaries

10. 2.1 Smoothed $\ell_{p}$-norm functions

11. 2.2 Flow Problems and Approximation

12. 2.3 Approximating Smoothed $p$-norm Instances

13. 2.4 Orthogonal Decompositions of Flows

14. 3 Numerical Methods

15. 3.1 Iterative Refinement

16. 3.2 Vertex Elimination

17. 3.3 Recursive Preconditioning

18. 4 Graph Theoretic Preconditioners

19. 4.1 Tree-Portal Routing

20. 4.2 Partitioning Trees into Subtrees and Portals

21. 4.3 Graph Sparsification

22. 4.4 Ultra-sparsification Algorithm and Error Analysis

23. 5 Decomposing into Uniform Expanders

24. Acknowledgements

25. A Deferred Proofs from Prelims, Section 2.2

26. B Deferred Proofs for Numerical Methods from Section 3

27. C Elimination of Low-Degree Vertices, and Loops

28. D Sparsifying Uniform Expanders

29. E Using Approximate Projections

30. F $\ell_{p}$-norm Semi-Supervised Learning on Graphs.

1 Introduction

Graphs are among the most ubiquitous representations of data, and efficiently computing on graphs is a task central to operations research, machine learning, and network science. Among graph algorithms, network flows have been extensively studied  [EK72, Kar73, ET75, GT88, GR98, Sch02, Hoc08, CKM*+*11, HO13, Orl13, GT14], and have wide ranges of applications [KBR07, LSBG13, PZZ13]. Over the past decade, the ‘Laplacian paradigm’ of designing graph algorithms spurred a revolution in the best run-time upper bounds for many fundamental graph optimization problems. Many of these new graph algorithms incorporated numerical primitives: even for the $s$-$t$ shortest path problem in graphs with negative edge weights, the current best running times [CMTV17] are from invoking linear system solvers.

This incorporation of numerical routines [DS08, CKM*+*11] in turn led to a dependence on $\epsilon$, the approximation accuracy. While maximum flow and transshipment problems on undirected graphs can now be approximated in nearly-linear time [KLOS14, She13, BKKL17, Pen16, She17b, She17a] (and the distributed setting has also been studied [GKK*+*15, BKKL17]), these algorithms are low accuracy in that their running times have factors of $1/\epsilon$ or higher. This is in contrast to high accuracy solvers for linear systems and convex programs, which with $polylog(n)$ overhead give $1/poly(n)$-approximate solutions. Prior to our result, such high accuracy runtime bounds for problems beyond linear systems all utilize second order methods [DS08, Mad13, LS14, CMTV17, ALdOW17, BCLL18] from convex optimization.

The main contribution of this paper is giving almost-linear time, high accuracy solutions to a significantly wider range of graph optimization problems that can be viewed as interpolations between maximum flow, shortest paths, and graph-structured linear systems. Our unified formulation of these problems is based on the following unified formulation of flow/path problems as norm minimization over a demand vector $\boldsymbol{\mathit{b}}\in\mathbb{R}_{\geq 0}^{V}$:

[TABLE]

In particular, when \mathopen{}\mathclose{{}\left\|\cdot}\right\|_{\odot} is the $\ell_{\infty}$-norm, this formulation is equivalent finding the flow of minimum congestion, which is in turn equivalent to computing maximum flows and bipartite matchings in unit capacitated graphs [Mad11]. Our main result is that for any $p\geq 2,$ given weights $\boldsymbol{\mathit{r}}\in\mathbb{R}_{\geq 0}^{E},$ a “gradient” $\boldsymbol{\mathit{g}}\in\mathbb{R}^{E}$, and a demand vector $\boldsymbol{\mathit{b}}\in\mathbb{R}^{V}$ (with $\boldsymbol{\mathit{b}}^{\top}\boldsymbol{\mathbf{1}}=0$), we can solve

[TABLE]

to $\nicefrac{{1}}{{{\textrm{poly}}(n)}}$ additive error in time $2^{O(p^{3/2})}m^{1+O(\frac{1}{\sqrt{p}})}$. We will formally state this result as Theorem 1.1 at the start of Section 1.1, and discuss several of its applications in flows, semi-supervised learning, and total variation minimization.

We believe that our algorithm represents a new approach to designing high accuracy solvers for graph-structured optimization problems. A brief survey of relevant works is in Section 1.2: previous high accuracy algorithms treat linear systems as the separation between graph theoretic and numerical components: the outer loop adjusts the numerics, while the inner loop quickly solves the resulting linear systems using the underlying graph structures. Our result, in contrast, directly invoke analogs of linear system solving primitives to the non-linear (but still convex) objective functions, and no longer has this clear separation between graph theoretic and numerical components.

We will overview key components of our approach, as well as how they are combined, in Section 1.3. Discussions of possible avenues for addressing shortcomings of our result, namely the exponential dependence on $p$, the restriction to unweighted graphs, and gap between $\ell_{\sqrt{\log{n}}}$-norm flow and $\ell_{\infty}$ are in Section 1.4.

1.1 Main Results and Applications

The formal formulation of our problem relies on the following objects defined on a graph $G=(V,E)$ with $n$ vertices and $m$ edges:

edge-vertex incidence matrix $\boldsymbol{\mathit{B}}$, 2. 2.

a vector $\boldsymbol{\mathit{b}}$ indicating the required residues on vertices (satisfying $\boldsymbol{\mathbf{1}}^{T}\boldsymbol{\mathit{b}}=0$), and 3. 3.

a norm $p$ as finding a flow $\boldsymbol{\mathit{f}}$ with demands $\boldsymbol{\mathit{b}}$ that minimize a specified norm $\|\cdot\|$.

The normed flow problem that we solve can then be formulated as:

[TABLE]

Using \mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|_{2,\boldsymbol{\mathit{r}}}=\sqrt{\sum_{e}\boldsymbol{\mathit{r}}_{e}\boldsymbol{\mathit{f}}_{e}^{2}} to denote the $\boldsymbol{\mathit{r}}$-weighted 2-norm, the objective can also be viewed as \boldsymbol{\mathit{g}}^{\top}\boldsymbol{\mathit{f}}+\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|_{2,\boldsymbol{\mathit{r}}}^{2}+\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|_{p}^{p}. Let $\text{val}(\boldsymbol{\mathit{f}})$ denote value of a flow $\boldsymbol{\mathit{f}}$ according to the above objective, and let OPT denote value of the optimal solution to Problem (2). Our main technical result is the following statement which we prove as corollary of our main technical theorem in Section 3.3.

Theorem 1.1** (Smoothed $\ell_{p}$-norm flows).**

For any $p\geq 2,$ given weights $\boldsymbol{\mathit{r}}\in\mathbb{R}_{\geq 0}^{E},$ a “gradient” $\boldsymbol{\mathit{g}}\in\mathbb{R}^{E}$, a demand vector $\boldsymbol{\mathit{b}}\in\mathbb{R}^{V}$ (with $\boldsymbol{\mathit{b}}^{\top}\boldsymbol{\mathbf{1}}=0$), and an initial solution $\boldsymbol{\mathit{f}}^{(0)}$ such that all parameters are bounded by $2^{{\textrm{poly}}(\log n)},$ we can compute a flow $\boldsymbol{\widetilde{f}}$ satisfying demands $\boldsymbol{\mathit{b}},$ i.e., $\boldsymbol{\mathit{B}}^{G\top}\boldsymbol{\widetilde{f}}=\boldsymbol{\mathit{b}},$ such that

[TABLE]

in $2^{O(p^{\nicefrac{{3}}{{2}}})}m^{1+O(\frac{1}{\sqrt{p}})}$ time, where $m$ denotes the number of edges in $G.$

1.1.1 $\ell_{p}$-Norm Flows

From this, we also get a (slightly simpler) statement about $\ell_{p}$-norm flows.

Theorem 1.2** ($\ell_{p}$-norm flows).**

For any $p\geq 2,$ given an unweighted graph $G(V,E)$ and demands $\boldsymbol{\mathit{b}}$, using the routine $\textsc{pFlows}(\mathcal{G},\boldsymbol{\mathit{b}})$ (Algorithm 2) we can compute a flow $\boldsymbol{\widetilde{f}}$ satisfying satisfying $\boldsymbol{\mathit{b}},$ i.e., $\boldsymbol{\mathit{B}}^{G\top}\boldsymbol{\widetilde{f}}=\boldsymbol{\mathit{b}},$ such that

[TABLE]

in $2^{O(p^{\nicefrac{{3}}{{2}}})}m^{1+O(\frac{1}{\sqrt{p}})},$ time, where $m$ denotes the number of edges in $G.$

This corollary is also proven in Section 3.3.

Picking $\boldsymbol{\mathit{g}},\boldsymbol{\mathit{r}}=\boldsymbol{\mathbf{0}}$ gives us an $2^{O(p^{3/2})}m^{1+O(\frac{1}{\sqrt{p}})}$ time high-accuracy algorithm for $\ell_{p}$-norm minimizing flows on unit weighted undirected graphs ($p\geq 2$). For large $p$, e.g. $p=\sqrt{\log n}$ this is an $m^{1+o(1)}$ time algorithm, and to our knowledge the first almost linear time high-accuracy algorithm for a flow problem other than Laplacian solvers ($\ell_{2}$) or shortest-paths ($\ell_{1}$).

1.1.2 Semi-Supervised Learning on Graphs.

Semi-supervised learning on graphs in machine learning is often based on solving an optimization problem where voltages (labels) are fixed at some vertices in a graph the voltages at remaning nodes are chosen so that some overall objective is minimized (e.g. $\ell_{p}$-norm of the vector of voltage differences across edges) [AL11, KRSS15, EACR*+*16]. Formally, given a graph $G=(V,E)$ and a labelled subset of the nodes $T\subset V$ with labels $\boldsymbol{\mathit{s}}_{T}\in\mathbb{R}^{T}$, we can write the problem as

[TABLE]

By converting this problem to its dual, we get an almost linear time algorithm for solving it to high accuracy, provided the initial voltage problem uses $p$ close to 1: In this case, voltage solutions are “cut-like”. Given $p<2$, we get a solver that computes a $(1+1/{\textrm{poly}}(m))$ multiplicative accuracy solution in time $2^{O((\frac{1}{p-1})^{\nicefrac{{3}}{{2}}})}m^{1+O(\sqrt{p-1})}$. For $p=1+\frac{1}{\sqrt{\log n}}$, this is time is bounded by $m^{1+o(1)}$.

Converting the dual of Problem (3) into a form solvable by our algorithms requires a small transformation, which we describe in Appendix F.

1.1.3 Use as Oracle in Conjunction with Multiplicative Weight Updates

The mixed $\ell_{2}^{2}$ and $\ell_{p}^{p}$ objective in our Problem (2) is useful for building oracles to use in multiplicative weight update algorithms based on flows, as they appear in [CKM*+*11, AKPS19]. Assume we are looking to solve some problem to $(1+\epsilon)$-accuracy measured as multiplicative error, and let us assume $\frac{1}{{\textrm{poly}}(m)}<\epsilon<0.5$. Specifically we can solve for the following objective subject to certain linear constraints.

[TABLE]

This gives an oracle for several problems. Algorithms based on oracle solutions to this type of objective work by noting that any $\boldsymbol{\mathit{f}}$ with \mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|_{\infty}\leq 1 gives an objective value at most

[TABLE]

Since such a flow must exist in the context where the oracle is applied, the optimum flow must also meet this bound. Now, if we compute a $(1+0.01\epsilon)$ approximately optimal solution to this problem, it must satisfy

[TABLE]

By Cauchy-Schwarz, we get \sum_{e}\boldsymbol{\mathit{r}}_{e}\mathopen{}\mathclose{{}\left|\boldsymbol{\mathit{f}}_{e}}\right|\leq\sqrt{\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{r}}}\right\|_{1}\sum_{e}\boldsymbol{\mathit{r}}_{e}\boldsymbol{\mathit{f}}_{e}^{2}}\leq\mathopen{}\mathclose{{}\left(1+1.1\epsilon}\right)\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{r}}}\right\|_{1}, which tells us the oracle is “good-on-average” according to the weights $\boldsymbol{\mathit{r}}$. The objective value also implies for every edge that

[TABLE]

which simplifies to:

[TABLE]

when we set $p=\log^{0.1}n$. This is the width of the oracle, and together these conditions demonstrate that the oracle suffices for a multiplicative weights algorithm and bounds the number of calls to the oracle by m^{o\mathopen{}\mathclose{{}\left(1}\right)}{\textrm{poly}}(1/\epsilon).

This oracle has multiple uses:

Approximate undirected maximum flow.

Using the oracle, we can approximate maximum flow using [CKM*+*11], giving an algorithm for undirected maximum flow that is not based on oblivious routings unlike other fast algorithms for approximate maximum flow [She13, KLOS14, Pen16]. Our algorithm obtains almost-linear time, albeit only for unit weighted graphs.

Isotropic total variation denoising.

Using our algorithm, we can give the first almost linear time, low accuracy algorithm for total variation denoising on unit weighted graphs [ROF92, ZWC10]. While there has been significant advances in image processing since the introduction of this objective, it still remains a representative objective in pixel vision tasks. The total variation objectives can be viewed as variants of semi-supervised learning on graphs: Given a “signal” vector $\boldsymbol{\mathit{s}}$ which corresponds to noisy observations of pixels of an image, we want to find a denoised version of $\boldsymbol{\mathit{s}}$, which we refer to as $\boldsymbol{\mathit{x}}$. The denoised output $\boldsymbol{\mathit{x}}$ should minimize an objective that measures both the between pixels in $\boldsymbol{\mathit{x}}$ that are close to each other in the image (which should be small), and the difference between $\boldsymbol{\mathit{x}}$ and $\boldsymbol{\mathit{s}}$ (which should also be small). The most popular version of this problem, known as isotropic total variation denoising, allows the input to specify a collection of groups of pixels with connections inside each group $i$ given by a set of edges $E_{i}$, and asks that 1) the denoised pixels are close in an $\ell_{2}$ sense to the measured signal, 2) in each group, the standard deviation between denoised pixels is not too high. These goals are expressed in the objective

[TABLE]

The dual of this problem is grouped flows, which is finding $\boldsymbol{\mathit{f}}$ such that $\boldsymbol{\mathit{B}}^{\top}\boldsymbol{\mathit{f}}=\boldsymbol{\mathit{d}}$ and for edge sets $E_{i}$,

[TABLE]

Our oracle gives the first routine for approximate isotropic $TV$ denoising that runs in almost linear time. The previous best running time was about $m^{4/3}$ [CMMP13].

1.2 Related Work

Network flow problems have a long history of motivating broader developments in algorithms, including the introduction of strongly polynomial time as a benchmark of algorithmic efficiency [Edm65, EK72], the development of tree data structures [GN79, ST83, ST85], and randomized graph algorithms and graph approximations [KS96, BK96]. For general capacities, the best strongly polynomial time algorithms run in about quadratic time due to the flow decomposition barrier [EK72, GN79, GT88, HO13, Orl13], which says that the there exists graphs where the path decomposition of an optimum flow must have quadratic size.

The flow decomposition barrier suggest that sub-quadratic time algorithms for network flows should decompose solutions numericallly, and this has indeed been the case in the development of such algorithms [GR98, GT14]. These numerical approaches recently culminated in nearly-linear time algorithms for undirected maximum flow and transshipment (the $\ell_{1}$ case of Problem (1)), yielding nearly-linear time algorithms [CKM*+*11, She13, KLOS14, Pen16, She17b, She17a]. Much of these progress were motivated by the development of nearly-linear time high-accuracy solvers for Laplacian linear systems [ST14, KMP12, KOSZ13, LS13, KS16], whose duals, electrical flows are the $\ell_{\infty}$ case of Problem (1). Such solvers can in turn be used to give the current best high accuracy flow algorithms. For graphs with polynomially bounded capacities, the current best running time is $\widetilde{O}(m\sqrt{n})$ due to Lee and Sidford [LS14]. On sparse graphs, this bound still does not break the long-standing $O(n^{1.5})$ barrier dating back to the early 70s [HK73, Kar73, ET75]. Recently Madry [Mad13, Mad16] broke this barrier on unit capacitated graphs, obtaining $\widetilde{O}(m^{10/7})$ running time.

Our result has in common with all previous results on almost-linear time optimization problems on graphs [KLOS14, She13, BKKL17, Pen16, She17b, She17a] in that it is based on white-box modifications of a linear system solver. In particular, our high level algorithmic strategy in creating edge and vertex reductions is identical to the first nearly-linear time solver by Spielman and Teng [ST14]. Much of this similarity is due to the lack of understanding of more general versions of key components: some possibilities for simplifying the result will be discussed in Section 1.4. On the other hand, our algorithms differ from previous adaptations of solvers in that it obtains high accuracy 555 The nearly-linear time matrix scaling algorithm [CMTV17] has a linear dependence on the condition number $\kappa$, while convex optimization methods for matrix scaling have dependencies of $\log\kappa$ instead. . This requires us to tailor the scheme to the residual problems from the $p$-norm iterative methods, and results in us taking a more numerical approach, instead of the more routing and path embedding-based approaches utilized in similar adaptations of Spielman and Teng [ST14] to cuts [Mad10], flows [She13, KLOS14], and shortest paths [BKKL17].

The development of high-accuracy algorithms for $p$-norm minimization that are faster than interior point methods (IPMs) [NN94] was pioneered by the recent work of Bubeck et al. [BCLL18] which introduced the $\gamma$-functions that were also used in [AKPS19]. However, the methods in [BCLL18] are conceptually similar to interior point methods (IPMs) [NN94] (as in they are homotopy methods). Their runtime for large $p$ behaves essentially like IPMs, requiring about $m^{\nicefrac{{3}}{{2}}-o(1)}$ time for solving $p$-norm flow problems, whereas the limiting behavior of our result is about $m^{1+o(1)}$.

1.3 Overview

At a high level, our approach can be viewed as solving a graph optimization problem as a linear system. This is done by combining the numerical methods for $\ell_{p}$-norms by Adil et al. [AKPS19] with the recursive preconditioning of graph structured linear systems by Spielman and Teng [ST14]. Many conceptual obstacles arise in trying to realize this vision, preventing us from adopting later Laplacian linear solvers that have greatly simplified the result of Spielman and Teng. The main one is the lack of concentration theory for the smoothed $p$-norm objectives integral to our algorithms: these concentration arguments are at the core of all subsequent improvements to Laplacian solver algorithms [KMP11, KOSZ13, LPS15, KS16].

Our starting point is a recent result involving a subset of the authors [AKPS19] that significantly generalized the phenomenon of high-accuracy numerical methods. In particular, this method is applicable to general $\ell_{p}$-norm optimization problems, for all $p$ that are bounded away from $1$ and $\infty$. It also opens up a major question: can we develop an appropriate notion of preconditioning, the other central ingredient of fast solvers for linear systems, applicable to $\ell_{p}$-norms? We resolve this question in the affirmative, and develop a theory of preconditioning that works for a wide class of non-linear problems in Section 3. In particular, we show that the second and pth order terms from the main formulation in Equation 2 form a class of functions that’s closed under taking residual problems. We will formally define these as smoothed $\ell_{p}$-norms in Section 2.1.

The crux of our problem then becomes designing preconditioners that interact well with these smoothed $\ell_{p}$-norms. Here it’s worth noting that earlier works on preconditioning for non-linear (maximum) flow problems all relied on oblivious routing which gives rise to linear preconditioners. Such an approach encounters a significant obstacle with $\ell_{p}$ norms: consider changing a single coordinate from, say $1$, to $(1+\delta)$:

• •

If the update $\delta$ is much smaller than $1$ in absolute value, the change in the objective from $1^{p}$ to $(1+\delta)^{p}$ is dominated by terms that scale as $\delta$ and $\delta^{2}$.

• •

However, if the update is much larger than $1$, the change is dominated by a $\delta^{p}$ term.

This means that good preconditioning across small and large updates is inherently highly dependent on the current coordinate value.

This example captures the core difficulties of our preconditioned iterative methods for smoothed $\ell_{p}$-norm problems, which heavily rely on both the second and pth power terms the objective functions. It means our graph theoretic components must simultaneously control terms of different degrees (namely scaling as $\delta$, $\delta^{2}$, and $\delta^{p}$) related to the flows on graphs. Here our key idea is that unit-weighted graphs have “multi-objective low-stretch trees” that simultaneously preserve the $\delta^{2}$ and $\delta^{p}$ terms, while the linear (gradient) terms can be preserved exactly when routing along these trees. Here a major difficulty is that the tree depends on the second order derivatives of the current solution point, and thus must continuously change as the algorithm proceeds. Additionally, after rerouting graph edges along the tree, we need to sparsify the graph according to weights defined by the same second derivatives at the current solution, which makes the adaptive aspect of the algorithm even more important. We defer the construction of our adaptive preconditioner to Section 4, after first formally defining our objective functions in Section 2, and introducing numerical methods based on them in Section 3.

1.4 Open Questions

We expect that our algorithm can be greatly simplified and adapted to non-unit weight graphs in ways similar to the sampling based solvers for Laplacian linear systems [KMP14, KOSZ13, KS16]. The current understanding of concentration theory for $\ell_{p}$ norms rely heavily on tools from functional analysis [CP15]: generalizing these tools to smoothed $\ell_{p}$-norm objectives is beyond the scope of this paper.

A major limitation of our result is the restriction to unit capacitated graphs. We believe this limitation is inherent to our approach of constructing preconditioners from trees: for general weights, there are cases where no tree can simultaneously have small stretch w.r.t. $\ell_{2}$-norm and $\ell_{p}$-norm weights. We believe that by developing a more complete theory of elimination and sparsification for these objectives, it will be possible to sparsify non-unit weight graphs, and develop solvers following the patterns of non-tree based Laplacian solvers [PS14, LPS15, KS16].

We also believe that the overall algorithmic approach established here is applicable far beyond the class of objective functions studied in this paper. Here a direct question is whether the dependency on $p$ can be improved to handling $\ell_{m}$ flows, which in unit weighted graphs imply maximum flows. The exponential dependence on $p$ has already been shown to be improvable to about $\widetilde{O}(p^{2})$ [Sac19]. For even larger values of $p$, a natural approach is to use homotopy methods that modify the $p$ values gradually. Here it is also plausible that our techniques, or their possible generalizations to weighted cases, can be used as algorithmic building blocks.

2 Preliminaries

2.1 Smoothed $\ell_{p}$-norm functions

We consider $p$-norms smoothed by the addition of a quadratic term. First we define such a smoothed $p\textsuperscript{th}$-power on $\mathbb{R}.$

Definition 2.1** (Smoothed $p\textsuperscript{th}$-power).**

Given $r,x\in\mathbb{R},r\geq 0$ define the $r$-smoothed $s$-weighted $p\textsuperscript{th}$-power of $x$ to be

[TABLE]

This definition can be naturally extended to vectors to obtained smoothed $\ell_{p}$-norms.

Definition 2.2** (Smoothed $\ell_{p}$-norm).**

Given vectors $\boldsymbol{\mathit{x}}\in\mathbb{R}^{m},\boldsymbol{\mathit{r}}\in\mathbb{R}^{m}_{\geq 0},$ and a positive scalar $s\in\mathbb{R}_{\geq 0},$ define the $\boldsymbol{\mathit{r}}$-smooth $s$-weighted $p$-norm of $\boldsymbol{\mathit{x}}$ to be

[TABLE]

2.2 Flow Problems and Approximation

We will consider problems where we seek to find flows minimizing smoothed $p$-norms. We first define these problem instances.

Definition 2.3** (Smoothed $p$-norm instance).**

A smoothed $p$-norm instance is a tuple $\mathcal{G},$

[TABLE]

where $V^{\mathcal{G}}$ is a set of vertices, $E^{\mathcal{G}}$ is a set of undirected edges on $V^{\mathcal{G}},$ the edges are accompanied by a gradient, specified by $\boldsymbol{\mathit{g}}^{\mathcal{G}}\in\mathbb{R}^{E^{\mathcal{G}}},$ the edges have $\ell_{2}^{2}$-resistances given by $\boldsymbol{\mathit{r}}^{\mathcal{G}}\in\mathbb{R}^{E^{\mathcal{G}}}_{\geq 0},$ and $s\in\mathbb{R}_{\geq 0}$ gives the $p$-norm scaling.

Definition 2.4** (Flows, residues, and circulations).**

Given a smoothed $p$-norm instance $\mathcal{G},$ a vector $\boldsymbol{\mathit{f}}\in\mathbb{R}^{E^{\mathcal{G}}}$ is said to be a flow on $\mathcal{G}$. A flow vector $\boldsymbol{\mathit{f}}$ satisfies residues $\boldsymbol{\mathit{b}}\in\mathbb{R}^{V^{\mathcal{G}}}$ if \mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{B}}^{\mathcal{G}}}\right)^{\top}\boldsymbol{\mathit{f}}=\boldsymbol{\mathit{b}}, where $\boldsymbol{\mathit{B}}^{\mathcal{G}}\in\mathbb{R}^{E^{\mathcal{G}}\times V^{\mathcal{G}}}$ is the edge-vertex incidence matrix of the graph $(V^{\mathcal{G}},E^{\mathcal{G}}),$ i.e., \mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{B}}^{\mathcal{G}}}\right)^{\top}_{(u,v)}=\boldsymbol{\mathbf{1}}_{u}-\boldsymbol{\mathbf{1}}_{v}.

A flow $\boldsymbol{\mathit{f}}$ with residue $\boldsymbol{\mathbf{0}}$ is called a circulation on $\mathcal{G}$.

Note that our underlying instance and the edges are undirected. However, for every undirected edge $e=(u,v)\in E$, we assign an arbitrary fixed direction to the edge, say $u\to v,$ and interpret $\boldsymbol{\mathit{f}}_{e}\geq 0$ as flow in the direction of the edge from $u$ to $v,$ and $\boldsymbol{\mathit{f}}_{e}<0$ as flow in the reverse direction. For convenience, we assume that for any edge $(u,v)\in E,$ we have $\boldsymbol{\mathit{f}}_{(u,v)}=-\boldsymbol{\mathit{f}}_{(v,u)}.$

Definition 2.5** (Objective, $\mathcal{E}^{\mathcal{G}}$).**

Given a smoothed $p$-norm instance $\mathcal{G},$ and a flow $\boldsymbol{\mathit{f}}$ on $\mathcal{G},$ the associated objective function, or the energy, of $\boldsymbol{\mathit{f}}$ is given by

[TABLE]

Definition 2.6** (Smoothed $p$-norm flow / circulation problem).**

Given a smoothed $p$-norm instance $\mathcal{G}$ and a residue vector $\boldsymbol{\mathit{b}}\in\mathbb{R}^{E^{\mathcal{G}}},$ the smoothed $p$-norm flow problem $(\mathcal{G},\boldsymbol{\mathit{b}})$, finds a flow $\boldsymbol{\mathit{f}}\in\mathbb{R}^{E^{\mathcal{G}}}$ with residues $\boldsymbol{\mathit{b}}$ that maximizes $\mathcal{E}^{\mathcal{G}}(\boldsymbol{\mathit{f}}),$ i.e.,

[TABLE]

If $\boldsymbol{\mathit{b}}=\boldsymbol{\mathbf{0}},$ we call it a smoothed $p$-norm circulation problem.

Note that the optimal objective of a smoothed $p$-norm circulation problem is always non-negative, whereas for a smoothed $p$-norm flow problem, it could be negative.

2.3 Approximating Smoothed $p$-norm Instances

Since we work with objective functions that are non-standard (and not even homogeneous), we need to carefully define a new notion of approximation for these instances.

Definition 2.7** ($\mathcal{H}\preceq_{\kappa}\mathcal{G}$).**

For two smoothed $p$-norm instances, $\mathcal{G},\mathcal{H},$ we write $\mathcal{H}\preceq_{\kappa}\mathcal{G}$ if there is a linear map $\mathcal{M}_{\mathcal{H}\rightarrow\mathcal{G}}:\mathbb{R}^{E^{\mathcal{H}}}\rightarrow\mathbb{R}^{E^{\mathcal{G}}}$ such that for every flow $\boldsymbol{\mathit{f}}^{\mathcal{H}}$ on $\mathcal{H},$ we have that $\boldsymbol{\mathit{f}}^{\mathcal{G}}=\mathcal{M}_{\mathcal{H}\rightarrow\mathcal{G}}(\boldsymbol{\mathit{f}}^{\mathcal{H}})$ is a flow on $\mathcal{G}$ such that

$\boldsymbol{\mathit{f}}^{\mathcal{G}}$ has the same residues as $\boldsymbol{\mathit{f}}^{\mathcal{H}}$ i.e., $(\boldsymbol{\mathit{B}}^{\mathcal{G}})^{\top}\boldsymbol{\mathit{f}}^{\mathcal{G}}=(\boldsymbol{\mathit{B}}^{\mathcal{H}})^{\top}\boldsymbol{\mathit{f}}^{\mathcal{H}},$ and 2. 2.

has energy bounded by:

[TABLE]

For some of our transformations on graphs, we will be able to prove approximation guarantees only for circulations. Thus, we define the following notion restricted to circulations.

Definition 2.8** ($\mathcal{H}\preceq^{\text{cycle}}_{\kappa}\mathcal{G}$).**

For two smoothed $p$-norm instances, $\mathcal{G},\mathcal{H},$ we write $\mathcal{H}\preceq^{\text{cycle}}_{\kappa}\mathcal{G}$ if there is a linear map $\mathcal{M}_{\mathcal{H}\rightarrow\mathcal{G}}:\mathbb{R}^{E^{\mathcal{H}}}\rightarrow\mathbb{R}^{E^{\mathcal{G}}}$ such that for any circulation $\boldsymbol{\mathit{f}}^{\mathcal{H}}$ on $\mathcal{H}$, i.e., $(\boldsymbol{\mathit{B}}^{\mathcal{H}})^{\top}\boldsymbol{\mathit{f}}^{\mathcal{H}}=\boldsymbol{\mathbf{0}},$ the flow $\boldsymbol{\mathit{f}}^{\mathcal{G}}=\mathcal{M}_{\mathcal{H}\rightarrow\mathcal{G}}(\boldsymbol{\mathit{f}}^{\mathcal{H}})$ is a circulation, i.e., $(\boldsymbol{\mathit{B}}^{\mathcal{G}})^{\top}\boldsymbol{\mathit{f}}^{\mathcal{G}}=\boldsymbol{\mathbf{0}},$ and satisfies

[TABLE]

Observe that $\mathcal{H}\preceq_{\kappa}\mathcal{G}$ implies $\mathcal{H}\preceq^{\text{cycle}}_{\kappa}\mathcal{G}.$

These definitions satisfy most properties that we want from comparisons.

Lemma 2.9** (Reflexivity).**

For every smoothed $p$-norm instance $\mathcal{G},$ and every $\kappa\geq 1$, $\mathcal{G}\preceq_{\kappa}\mathcal{G}$ and $\mathcal{G}\preceq^{\text{cycle}}_{\kappa}\mathcal{G}$ with the identity map.

It behaves well under composition.

Lemma 2.10** (Composition).**

Given two smoothed $p$-norm instances, $\mathcal{G}_{1},\mathcal{G}_{2},$ such that $\mathcal{G}_{1}\preceq_{\kappa_{1}}\mathcal{G}_{2}$ with the map $\mathcal{M}_{\mathcal{G}_{1}\rightarrow\mathcal{G}_{2}}$ and $\mathcal{G}_{2}\preceq_{\kappa_{2}}\mathcal{G}_{3}$ with the map $\mathcal{M}_{\mathcal{G}_{2}\rightarrow\mathcal{G}_{3}}$, then $\mathcal{G}_{1}\preceq_{\kappa_{1}\kappa_{2}}\mathcal{G}_{3}$ with the map $\mathcal{M}_{\mathcal{G}_{1}\rightarrow\mathcal{G}_{3}}=\mathcal{M}_{\mathcal{G}_{2}\rightarrow\mathcal{G}_{3}}\circ\mathcal{M}_{\mathcal{G}_{1}\rightarrow\mathcal{G}_{2}}.$

Similarly, for any $\mathcal{G}_{1},\mathcal{G}_{2},$ if $\mathcal{G}_{1}\preceq^{\text{cycle}}_{\kappa_{1}}\mathcal{G}_{2}$ with the map $\mathcal{M}_{\mathcal{G}_{1}\rightarrow\mathcal{G}_{2}}$ and $\mathcal{G}_{2}\preceq^{\text{cycle}}_{\kappa_{2}}\mathcal{G}_{3}$ with the map $\mathcal{M}_{\mathcal{G}_{2}\rightarrow\mathcal{G}_{3}}$, then $\mathcal{G}_{1}\preceq^{\text{cycle}}_{\kappa_{1}\kappa_{2}}\mathcal{G}_{3}$ with the map $\mathcal{M}_{\mathcal{G}_{1}\rightarrow\mathcal{G}_{3}}=\mathcal{M}_{\mathcal{G}_{2}\rightarrow\mathcal{G}_{3}}\circ\mathcal{M}_{\mathcal{G}_{1}\rightarrow\mathcal{G}_{2}}.$

The most important property of this is that this notion of approximation is also additive, i.e., it works well with graph decompositions.

Definition 2.11** (Union of two instances).**

Consider smoothed $p$-norm instances, $\mathcal{G}_{1},\mathcal{G}_{2},$ with the same set of vertices, i.e. $V^{\mathcal{G}_{1}}=V^{\mathcal{G}_{2}}.$ Define $\mathcal{G}=\mathcal{G}_{1}\cup\mathcal{G}_{2}$ as the instance on the same set of vertices obtained by taking a disjoint union of the edges (potentially resulting in multi-edges). Formally,

[TABLE]

Lemma 2.12** ($\preceq_{\kappa}$ under union).**

Consider four smoothed $p$-norm instances, $\mathcal{G}_{1},\mathcal{G}_{2},\mathcal{H}_{1},\mathcal{H}_{2},$ on the same set of vertices, i.e. $V^{\mathcal{G}_{1}}=V^{\mathcal{G}_{2}}=V^{\mathcal{H}_{1}}=V^{\mathcal{H}_{2}},$ such that for $i=1,2,$ $\mathcal{H}_{i}\preceq_{\kappa}\mathcal{G}_{i}$ with the map $\mathcal{M}_{\mathcal{H}_{i}\rightarrow\mathcal{G}_{i}}.$ Let $\mathcal{G}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathcal{G}_{1}\cup\mathcal{G}_{2},$ and $\mathcal{H}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathcal{H}_{1}\cup\mathcal{H}_{2}.$ Then, $\mathcal{H}\preceq_{\kappa}\mathcal{G}$ with the map

[TABLE]

where $(\boldsymbol{\mathit{f}}^{H_{1}},\boldsymbol{\mathit{f}}^{H_{2}})$ is the decomposition of $\boldsymbol{\mathit{f}}^{H}$ onto the supports of $H_{1}$ and $H_{2}$.

This notion of approximation also behaves nicely with scaling of $\ell_{2}$ and $\ell_{p}$ resistances.

Lemma 2.13**.**

For all $\kappa\geq 1,$ and for all pairs of smoothed $p$-norm instances, $\mathcal{G},\mathcal{H},$ on the same underlying graphs, i.e., $(V^{\mathcal{G}},E^{\mathcal{G}})=(V^{\mathcal{H}},E^{\mathcal{H}}),$ such that,

the gradients are identical, $\boldsymbol{\mathit{g}}^{\mathcal{G}}=\boldsymbol{\mathit{g}}^{\mathcal{H}},$ 2. 2.

the $\ell_{2}^{2}$ resistances are off by at most $\kappa,$ i.e., $\boldsymbol{\mathit{r}}_{e}^{\mathcal{G}}\leq\kappa\boldsymbol{\mathit{r}}_{e}^{\mathcal{H}}$ for all edges $e,$ and 3. 3.

the $p$-norm scaling is off by at most $\kappa^{p-1},$ i.e., $s^{\mathcal{G}}\leq\kappa^{p-1}s^{\mathcal{H}},$

then $\mathcal{H}\preceq_{\kappa}\mathcal{G}$ with the identity map.

2.4 Orthogonal Decompositions of Flows

At the core of our graph decomposition and sparsification procedures is a decomposition of the gradient $\boldsymbol{\mathit{g}}$ of $\mathcal{G}$ into its cycle space and potential flow space. We denote such a splitting using

[TABLE]

Here $\boldsymbol{\widehat{\mathit{g}}}$ is a circulation, while $\boldsymbol{\mathit{B}}\boldsymbol{\mathit{\psi}}$ gives a potential induced edge value. We will omit the superscripts when the context is clear.

The following minimization based formulation of this splitting of $\boldsymbol{\mathit{g}}$ is critical to our method of bounding the overall progress of our algorithm

Fact 2.14**.**

The projection of $\boldsymbol{\mathit{g}}$ onto the cycle space is obtained by minimizing the energy added to a potential flow to $\boldsymbol{\mathit{g}}$. Specifically,

[TABLE]

Lemma 2.15**.**

Given a graph/gradient instance $\mathcal{G}$, consider $\mathcal{H}$ formed from a subset of its edges. The projections of $\boldsymbol{\mathit{g}}^{\mathcal{G}}$ and $\boldsymbol{\mathit{g}}^{\mathcal{H}}$ onto their respective cycle spaces, $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}$ and $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}$ satsify:

[TABLE]

3 Numerical Methods

The general idea of (preconditioned) numerical methods, which are at the core of solvers for graph-structured linear systems [ST14] is to repeatedly update a current solution in ways that multiplicative reduce the difference in objective value to optimum. In the setting of flows, suppose we currently have some tentative solution $\boldsymbol{\mathit{f}}$ to the minimization problem

[TABLE]

by performing the step

[TABLE]

with the goal of improving the objective value substantially.

The work of Adil et al. [AKPS19] proved that $\ell_{p}$-norm minimization problems could be iteratively refined. While that result hinted at a much more general theory of numerical iterative methods for minimizing convex objectives, this topic is very much under development. In this section, we will develop the tools necessary for preconditioning $\ell_{p}$-norm based functions, and formalize the requirements for preconditioners necessary for recursive preconditioning algorithms.

3.1 Iterative Refinement

The following key Lemma from [AKPS19] allows us to approximate the change in the smoothed $p$-norm of $\boldsymbol{\mathit{x}}+\boldsymbol{\delta}$ relative to the norm of $\boldsymbol{\mathit{x}},$ in terms of another smoothed $p$-norm of $\boldsymbol{\delta}.$

Lemma 3.1** ([AKPS19]).**

For all $\boldsymbol{\mathit{r}},\boldsymbol{\mathit{x}},\boldsymbol{\delta}\in\mathbb{R}^{m},$ with $\boldsymbol{\mathit{r}}\in\mathbb{R}_{\geq 0}^{m},$ and $s\geq 0,$ we have

[TABLE]

The above lemma gives us the following theorem about iteratively refining smoothed $\ell_{p}$-norm minimization problems. While the lemma was essentially proven in  [AKPS19], they used slightly different definitions, and for completeness we prove the lemma in Appendix B.

The following theorem also essentially appeared in [AKPS19], but again as slightly different definitions were used in that paper, we prove the theorem in Appendix B for completeness.

Theorem 3.2** ([AKPS19]).**

Given the following optimization problem,

[TABLE]

and an initial feasible solution $\boldsymbol{\mathit{x}}_{0},$ we can construct the following residual problem:

[TABLE]

where \boldsymbol{\mathit{g}}^{\prime}=2^{p}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{g}}-\nabla_{\boldsymbol{\mathit{x}}}h(\boldsymbol{\mathit{r}},s,\boldsymbol{\mathit{x}})|_{\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{x}}_{0}}}\right), and \boldsymbol{\mathit{r}}^{\prime}=\boldsymbol{\mathit{r}}+s\mathopen{}\mathclose{{}\left|\boldsymbol{\mathit{x}}_{0}}\right|^{p-2}.

There exists a feasible solution $\widetilde{\boldsymbol{\delta}}$ to the residual problem 13 that achieves an objective of \mathcal{E}_{2}\mathopen{}\mathclose{{}\left(\widetilde{\boldsymbol{\delta}}}\right)\geq 2^{p}(\mathcal{E}_{1}(\boldsymbol{\mathit{x}}^{\star})-\mathcal{E}_{1}(\boldsymbol{\mathit{x}}_{0})), where $\boldsymbol{\mathit{x}}^{\star}$ is an optimal solution to problem 10.

Moreover, given any feasible solution $\boldsymbol{\delta}$ to Program 13, the vector $\boldsymbol{\mathit{x}}_{1}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\boldsymbol{\mathit{x}}_{0}+{2^{-3p}}\boldsymbol{\delta}$ is a feasible solution to the Program 10 and obtains the objective

[TABLE]

Importantly, we can apply the above theorem to smoothed $p$-norm flow/circulation problems.

Corollary 3.3** **(Iterative refinement for smoothed $p$-norm

flow/circulation problems).

Given any smoothed $p$-norm flow problem $(\mathcal{G},\boldsymbol{\mathit{b}})$ with optimal objective $\mathcal{E}^{\star}(\mathcal{G}),$ and any initial circulation $\boldsymbol{\mathit{f}}_{0},$ we can build, in O(\mathopen{}\mathclose{{}\left|E^{\mathcal{G}}}\right|) time, a smoothed $p$-norm circulation problem $\mathcal{H}$ with the same underlying graph $(V^{\mathcal{H}},E^{\mathcal{H}})=(V^{\mathcal{G}},E^{\mathcal{G}}),$ such that $\mathcal{E}^{\star}(\mathcal{H})\geq 2^{p}(\mathcal{E}^{\star}(\mathcal{G})-\mathcal{E}^{\mathcal{G}}(\boldsymbol{\mathit{f}}_{0}))$ and for any circulation $\boldsymbol{\mathit{f}}^{\mathcal{H}}$ on $\mathcal{H},$ the flow $\boldsymbol{\mathit{f}}_{1}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\boldsymbol{\mathit{f}}_{0}+2^{-3p}\boldsymbol{\mathit{f}}^{\mathcal{H}}$ satisfies residues $\boldsymbol{\mathit{b}}$ on $\mathcal{G}$ and has an objective

[TABLE]

This means if we find even a crude approximate minimizer $\boldsymbol{\widetilde{\mathit{\delta}}}$ of this update problem, we can move to a new point $\boldsymbol{\mathit{f}}^{\prime}=\boldsymbol{\mathit{f}}+\boldsymbol{\widetilde{\mathit{\delta}}}$, so that the gap to the optimum in the original optimization problem (7) will decrease by a constant factor (depending only on $p$) from \mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|_{p}^{p}-\text{OPT} to \mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{\prime}}\right\|_{p}^{p}-\text{OPT}\leq(1-2^{-O(p)})(\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|_{p}^{p}-\text{OPT}). In other words, we have a kind of iterative refinement: crude solutions to an update problem directly give constant factor progress in the original objective.

Note that \mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|_{p}^{p}=\sum_{i}\boldsymbol{\mathit{f}}_{i}^{p}. This will help us understand the objective function of the update problem coordinate-wise. Our update problem objective function is motivated by the following observations. Our function differs slightly from the function used in [AKPS19], which in turn was based on functions from [BCLL18], but our function still uses a few special properties of the [BCLL18] functions. Suppose $p\geq 2$ is an even integer (only to avoid writing absolute values), then

[TABLE]

Of course, the exact expansion gives

[TABLE]

So essentially we can approximate this expansion using only the zeroth, first, second, and last term in the expansion. We use $\boldsymbol{\mathit{g}}(\boldsymbol{\mathit{f}})$ to denote the vector with $\boldsymbol{\mathit{g}}_{i}(\boldsymbol{\mathit{f}})=p\boldsymbol{\mathit{f}}_{i}^{p-1}$ (i.e. the gradient), and let $\boldsymbol{\mathit{f}}^{p-2}$ denote the entrywise powered vector, and define $h_{p}(\boldsymbol{\mathit{f}}^{p-2},\boldsymbol{\delta})=\sum_{i}h_{p}(\boldsymbol{\mathit{f}}_{i}^{p-2},\boldsymbol{\delta}_{i})$. Thus we have

[TABLE]

Note that for any scalar $0<\lambda<1$,

[TABLE]

Together, these observations are enough to ensure that if we have $\boldsymbol{\widetilde{\mathit{\delta}}}$ which is a constant factor approximate solution to the follow optimization problem, which we define as our update problem

[TABLE]

then we can find a $\lambda$ s.t. \mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}+\lambda\boldsymbol{\widetilde{\mathit{\delta}}}}\right\|_{p}^{p}-\text{OPT}\leq(1-2^{-O(p)})(\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|_{p}^{p}-\text{OPT}).

But what have we gained? Why is Problem (15) more tractable than the one we started with?

A key reason is that unlike the exact expansion of an update as given by Equation (14), all the higher order terms in the objective function of (15) are coordinate-wise even functions, i.e. flipping the sign of a coordinate does not change the value of the function. [AKPS19] used a different but still even function instead of our $h_{p}$. This symmetrization allowed them to develop a multiplicative weights update algorithm motivated by [CKM*+*11] for their version of Problem (15), reducing the problem to solving a sequence of linear equations.

For our choice of $h_{p}$, it is particularly simple to show another very important property: Consider solving Problem (15) by again applying iterative refinement to this problem. That is, at some intermediate step with $\delta$ being the current solution, we aim to find an update $\boldsymbol{\widehat{\mathit{\delta}}}$ s.t. $\boldsymbol{\mathit{B}}\boldsymbol{\widehat{\mathit{\delta}}}=\boldsymbol{\mathbf{0}}$ and $\boldsymbol{\mathit{g}}(\boldsymbol{\mathit{f}})^{\top}(\boldsymbol{\delta}+\boldsymbol{\widehat{\mathit{\delta}}})+h_{p}(\boldsymbol{\mathit{f}}^{p-2},\boldsymbol{\delta}+\boldsymbol{\widehat{\mathit{\delta}}})$ is smaller than $\boldsymbol{\mathit{g}}(\boldsymbol{\mathit{f}})^{\top}(\boldsymbol{\delta})+h_{p}(\boldsymbol{\mathit{f}}^{p-2},\boldsymbol{\delta})$. Then by expanding the two non-linear terms of $h_{p}(\boldsymbol{\mathit{f}}_{i}^{p-2},\boldsymbol{\delta}_{i})$, i.e. $(\boldsymbol{\delta}_{i}+\boldsymbol{\widehat{\mathit{\delta}}})^{p}$ and $(\boldsymbol{\delta}_{i}+\boldsymbol{\widehat{\mathit{\delta}}})^{2}$, similar to Equation (14), we get a sequence of terms with powers of $\boldsymbol{\delta}_{i}$ ranging from $2$ to $p$. If we approximate this sequence again using only the $\boldsymbol{\widehat{\mathit{\delta}}}_{i}^{2}$ and $\boldsymbol{\widehat{\mathit{\delta}}}_{i}^{p}$ terms, we get another update problem. This update problem is an instance of Problem (2). And in general, we can set up iterative refinement update problems for instances of Problem (2), and get back another problem of the that class (after our approximation based on dropping intermediate terms). Thus, the problem class (2) is closed under repeatedly creating iterative update problems. This observation is central because it allows us to develop recursive algorithms.

3.2 Vertex Elimination

Following the template of the Spielman-Teng Laplacian solver, we recursively solve a problem on $m$ edges by solving about $\kappa$ problems on graphs with $n-1+m/\kappa$ edges. These ultra-sparse graphs allow us to eliminate degree 1 and 2 vertices and obtain a smaller problem. Because our update problem (Problem (15)) corresponds to a flow-circulation problem with some objective, we are able to understand elimination on these objectives in a relatively simple way: the flow on degree 1 and 2 vertices is easily related to flow in a smaller graph created by elimination. Unlike Spielman-Teng, every recursive call must rely on a new graph sparsifier, because the “graph” we sparsify depends heavily on the current solution that we are seeking to update: we have to simultaneously preserve linear, quadratic and $p$-th order terms, whose weights depend on the current solution.

A critical component of this schema is the mapping of flows back and forth between the original graph and the new graph so a good solution on a smaller graph can be transformed into a good solution on the larger graph. These mappings are direct analogs of eliminating degrees $1$ and $2$ vertices. In Appendix C, we generalize these processes to smoothed $\ell_{p}$-norm objectives, proving the following statements:

Theorem 3.4** (Eliminating vertices with degree 1 and 2).**

Given a smoothed $p$-norm instance $\mathcal{G},$ the algorithm $\textsc{Eliminate}(\mathcal{G})$ returns another smoothed $p$-norm instance $\mathcal{G}^{\prime},$ along with the map $\mathcal{M}_{\mathcal{G}^{\prime}\rightarrow\mathcal{G}}$ in O(\mathopen{}\mathclose{{}\left|V^{\mathcal{G}}}\right|+\mathopen{}\mathclose{{}\left|E^{\mathcal{G}}}\right|) time, such that the graph $G^{\prime}=(V^{\mathcal{G}^{\prime}},E^{\mathcal{G}^{\prime}})$ is obtained from the graph $G=(V^{\mathcal{G}},E^{\mathcal{G}})$ by first repeatedly removing vertices with non-selfloop degree666By non-selfloop degree, we mean that self-loops do not count towards the degree of a vertex. 1 in $G$, and then replacing every path $u\leadsto v$ in $G$ where all internal path vertices have non-selfloop degree exactly 2 in $G,$ with a new edge $(u,v).$

Moreover,

[TABLE]

where n=\mathopen{}\mathclose{{}\left|V^{\mathcal{G}}}\right|, and the map $\mathcal{M}_{\mathcal{G}^{\prime}\rightarrow\mathcal{G}}$ can be applied in O(\mathopen{}\mathclose{{}\left|V^{\mathcal{G}}}\right|+\mathopen{}\mathclose{{}\left|E^{\mathcal{G}}}\right|) time.

Lemma 3.5** (Eliminating Self-loops).**

There is an algorithm RemoveLoops such that, given a smoothed $p$-norm instance $\mathcal{G}$ with self-loops in $E^{\mathcal{G}},$ in O(\mathopen{}\mathclose{{}\left|V^{\mathcal{G}}}\right|+\mathopen{}\mathclose{{}\left|E^{\mathcal{G}}}\right|) time, it returns instances $\mathcal{G}_{1},\mathcal{G}_{2},$ such that $\mathcal{G}=\mathcal{G}_{1}\cup\mathcal{G}_{2},$ where $\mathcal{G}_{1}$ is obtained from $\mathcal{G}$ by eliminating all self-loops from $E^{\mathcal{G}},$ and $\mathcal{G}_{2}$ is an instance consisting of just the self-loops from $\mathcal{G}.$ Thus, any flow $\boldsymbol{\mathit{f}}^{\mathcal{G}_{2}}$ on $\mathcal{G}_{2}$ is a circulation.

Moreover, there is an algorithm SolveLoops that, given $\mathcal{G}_{2},$ for any $\delta\leq\nicefrac{{1}}{{p}},$ in time $O(|E^{\mathcal{G}_{2}}|\log\nicefrac{{1}}{{\delta}}),$ finds a circulation $\boldsymbol{\widetilde{f}}{}^{\mathcal{G}_{2}}$ on $\mathcal{G}_{2}$ such that

[TABLE]

We remark that only the map from the smaller graph to the larger has to be constructive.

3.3 Recursive Preconditioning

We can now present our main recursive preconditioning algorithm, RecursivePreconditioning (Algorithm 1). Earlier work on preconditioning for non-linear (maximum) flow problems all relied on oblivious routing which gives rise to linear preconditioners. These inherently cannot work well for high-accuracy iterative refinement, and the issue is not merely linearity: Consider Problem (15): the optimal $\boldsymbol{\delta}$ is highly dependent on the current $\boldsymbol{\mathit{f}}$, and when a coordinate $\boldsymbol{\delta}_{i}$ is large compared to the current \mathopen{}\mathclose{{}\left|\boldsymbol{\mathit{f}}_{i}}\right|, the function depends on it as $\boldsymbol{\delta}_{i}^{p}$, while when $\boldsymbol{\delta}_{i}$ is small compared to \mathopen{}\mathclose{{}\left|\boldsymbol{\mathit{f}}_{i}}\right|, it behaves as $\boldsymbol{\delta}_{i}^{2}$. Thus the behavior is highly dependent on the current solution. This necessitates adaptive and non-linear preconditioners.

To develop adaptive preconditioners, we employ recursive chains of alternating calls to non-linear iterative refinement and a new type of (ultra-)sparsification that is more general and stronger, allowing us to simultanously preserve multiple different properties of our problem. And crucially, every time our solution is updated, our preconditioners change. The central theorem governing the combinatorial components of our algorithm, which is the main result proven in Section 4, is:

Theorem 3.6** (Ultra-Sparsification).**

Given any instance $\mathcal{G}=(V^{\mathcal{G}},E^{\mathcal{G}},\boldsymbol{\mathit{g}}^{\mathcal{G}},\boldsymbol{\mathit{r}}^{\mathcal{G}},s^{\mathcal{G}})$ with $n$ nodes, $m$ edges, and parameters $\kappa,\delta$ where $\log\frac{1}{\delta}$ and \log\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{r}}^{\mathcal{G}}}\right\|_{\infty} are both $O(\log^{c}n)$ for some constant $c$ and $\kappa<m$ , UltraSparsify computes in $\widetilde{O}(m)$ running time another instance $\mathcal{H}=(V^{\mathcal{H}},E^{\mathcal{H}},\boldsymbol{\mathit{g}}^{\mathcal{H}},\boldsymbol{\mathit{r}}^{\mathcal{H}},s^{\mathcal{H}}=s^{\mathcal{G}})$ along with flow mapping functions $\mathcal{M}_{\mathcal{H}\rightarrow\mathcal{G}},\mathcal{M}_{\mathcal{G}\rightarrow\mathcal{H}}$ such that $V^{\mathcal{H}}=V^{\mathcal{G}}$, and with high probability we have

$E^{H}$* consists of a spanning tree in the graph $(V^{\mathcal{G}},E^{\mathcal{G}})$, up to $m-n+1$ self-loops and at most $\widetilde{O}(\frac{m}{\kappa})$ other non self-loop edges. * 2. 2.

With $\kappa_{\mathcal{G}\to\mathcal{H}}=\widetilde{O}(\kappa m^{3/(p-1)})$ for any flow $\boldsymbol{\mathit{f}}^{\mathcal{G}}$ of $\mathcal{G}$ we have

[TABLE]

and with $\kappa_{\mathcal{H}\to\mathcal{G}}=\widetilde{O}(m^{2/(p-1)})$, for any flow solution $\boldsymbol{\mathit{f}}^{\mathcal{H}}$ of $\mathcal{H}$ we have

[TABLE]

The flow mappings $\mathcal{M}_{\mathcal{H}\rightarrow\mathcal{G}},\mathcal{M}_{\mathcal{G}\rightarrow\mathcal{H}}$ preserve residue of flow, and can be applied in $\widetilde{O}(m)$ time.

Our key theorem about the performance of the algorithm is then:

Theorem 3.7** (Recursive Preconditioning).**

For all $p\geq 2,$ say we are given a smoothed $p$-norm instance $\mathcal{G},$ residues $\boldsymbol{\mathit{b}},$ initial solution $\boldsymbol{\mathit{f}}^{(0)},$ and $\delta\leq 1$ such that \log\nicefrac{{1}}{{\delta}},\log\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{g}}^{\mathcal{G}}}\right\|,\log\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{r}}^{\mathcal{G}}}\right\|,\log\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{(0)}}\right\|\leq\widetilde{O}(1). We can pick $\kappa=\widetilde{\Theta}(m^{\frac{1}{\sqrt{p-1}}})$ so that the procedure $\textsc{RecursivePreconditioning}(\mathcal{G},\boldsymbol{\mathit{b}},\boldsymbol{\mathit{f}}^{(0)},\kappa,\delta)$ runs in time $2^{O(p^{\nicefrac{{3}}{{2}}})}m^{1+O(\frac{1}{\sqrt{p-1}})},$ and returns a flow $\boldsymbol{\mathit{f}}$ on $\mathcal{G}$ such that $\boldsymbol{\mathit{f}}$ satisfies residues $\boldsymbol{\mathit{b}},$ and

[TABLE]

Proof.

(of Theorem 3.7) By scaling $\boldsymbol{\mathit{g}}^{\mathcal{G}},\boldsymbol{\mathit{r}}^{\mathcal{G}},$ we can assume that $s^{\mathcal{G}}=1$ without loss of generality.

Let us consider iteration $t$ of the for loop in RecursivePreconditioning. First, let us prove guarantees on the optimal solutions of all the relevant instances. By the guarantees of Corollary 3.3, we know that $\mathcal{H}_{1}$ is a smoothed $p$-norm circulation problem with the same underlying graph $(V^{\mathcal{H}_{1}},E^{\mathcal{H}_{1}})=(V^{\mathcal{G}},E^{\mathcal{G}}),$ such that $\mathcal{E}^{\star}(\mathcal{H}_{1})\geq 2^{p}(\mathcal{E}^{\star}(\mathcal{G})-\mathcal{E}^{\mathcal{G}}(\boldsymbol{\mathit{f}}^{(t)})).$

From Theorem 3.6, we know that UltraSparsify returns a smoothed $p$-norm circulation instance $\mathcal{H}_{2}$ on the same set of vertices such that \mathcal{E}^{\star}(\mathcal{H}_{2})\geq{\kappa^{-1}_{{\mathcal{H}_{1}}\to\mathcal{H}_{2}}}\mathcal{E}^{\star}(\mathcal{H}_{1})-\delta^{\prime}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{g}}^{\mathcal{H}_{1}}}\right\|\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{\star}{}^{\mathcal{H}_{1}}}\right\|.

From Theorem 3.4, we know that the instance $\mathcal{H}_{3}$ returned by ${\textsc{Eliminate}}(\mathcal{H}_{2})$ satisfies $\mathcal{H}_{2}\preceq^{\text{cycle}}_{1}\mathcal{H}_{3},$ and hence $\mathcal{E}^{\star}(\mathcal{H}_{2})\leq\mathcal{E}^{\star}(\mathcal{H}_{3}).$ From Lemma 3.5, we know that $\mathcal{E}^{\star}(\mathcal{H}_{4})+\mathcal{E}^{\star}(\mathcal{H}_{\text{loops}})=\mathcal{E}^{\star}(\mathcal{H}_{3}).$ Combining these guarantees, we obtain,

[TABLE]

Now, we analyze the approximation guarantee provided by the solutions to these instances. From Lemma 3.5, $\textsc{SolveLoops}(\mathcal{H}_{\text{loops}})$ returns a $\boldsymbol{\mathit{\Delta}}^{\mathcal{H}_{\text{loops}}}$ that satisfies $\mathcal{E}^{\mathcal{H}_{\text{loops}}}(\boldsymbol{\mathit{\Delta}}^{\mathcal{H}_{\text{loops}}})\geq\frac{1}{2}\alpha^{\star}(\mathcal{H}_{\text{loops}}).$ By induction, $\textsc{RecursivePreconditioning}(\mathcal{H}_{4},\boldsymbol{\mathbf{0}},\kappa,\delta T^{-1}),$ upon starting with the initial solution $\boldsymbol{\mathbf{0}},$ returns a $\boldsymbol{\mathit{\Delta}}^{\mathcal{H}_{4}}$ that satisfies, $\mathcal{E}^{\mathcal{H}_{4}}(\boldsymbol{\mathit{\Delta}}^{\mathcal{H}_{4}})\geq\frac{1}{2}\mathcal{E}^{\star}(\mathcal{H}_{4})-\delta T^{-1}$ . Combining these guarantees, we have,

[TABLE]

From Theorem 3.4, we also have $\mathcal{H}_{3}\preceq^{\text{cycle}}_{\kappa_{\text{elim}}}\mathcal{H}_{2},$ for \kappa_{\text{elim}}=\mathopen{}\mathclose{{}\left|V^{\mathcal{H}_{2}}}\right|^{\frac{1}{p-1}}, and hence

[TABLE]

Finally, from Theorem 3.6, we have,

[TABLE]

Combining these guarantees, we obtain,

[TABLE]

\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{\star}{}^{\mathcal{H}_{1}}}\right\|\leq\sqrt{m}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{g}}^{\mathcal{H}_{1}}}\right\|^{\frac{1}{p-1}}, since $\mathcal{E}^{\mathcal{H}_{1}}(\boldsymbol{\mathit{f}}^{\star}{}^{\mathcal{H}_{1}})\geq 0$ implies that m^{1-\nicefrac{{p}}{{2}}}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{\star}{}^{\mathcal{H}_{1}}}\right\|_{2}^{p}\leq\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{\star}{}^{\mathcal{H}_{1}}}\right\|_{p}^{p}\leq\boldsymbol{\mathit{g}}^{\mathcal{H}_{1}\top}\boldsymbol{\mathit{f}}^{\star}{}^{\mathcal{H}_{1}}\leq\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{g}}^{\mathcal{H}_{1}}}\right\|\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{\star}{}^{\mathcal{H}_{1}}}\right\|. Similarly, \mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{\Delta}}^{\mathcal{H}_{1}}}\right\|\leq\sqrt{m}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{g}}^{\mathcal{H}_{1}}}\right\|^{\frac{1}{p-1}}, and \mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{\Delta}}^{\mathcal{H}_{2}}}\right\|\leq\sqrt{m}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{\Delta}}^{\mathcal{H}_{1}}}\right\|, by the reverse tree routing map.

Thus, by Theorem 3.3, $\boldsymbol{\mathit{f}}^{(t+1)}$ satisfies

[TABLE]

Thus, repeating for loop $\widetilde{O}(2^{3p}\kappa m^{\frac{6}{p-1}})$ times gives us a solution $\boldsymbol{\mathit{f}}^{(T)}$ such that

[TABLE]

Now, we analyze the running time. In a single iteration, the total cost of all the operations other than the recursive call to RecursivePreconditioning is $\widetilde{O}(m).$ Note that $\mathcal{H}_{1}$ has $m$ edges. Theorem 3.6 tells us that $\mathcal{H}_{2}$ consists of a tree ($n-1$ edges), at most $\widetilde{O}(m/\kappa)$ non-selfloop edges, plus many self-loops. After invoking $\textsc{Eliminate},$ and $\textsc{RemoveLoops},$ the instance $\mathcal{H}_{4}$ has no self-loops left, and after dropping vertices with degree 0, only has vertices with degree at least 3. Every edge removed in Eliminate decreases the number of edges and vertices with non-zero non-selfloop degree by 1. Suppose $n^{\prime}$ is the number of vertices in $\mathcal{H}_{4}$ with non-zero degree. Then, Eliminate must have removed $n-n^{\prime}$ edges from $\mathcal{H}_{2}.$ Since $\mathcal{H}_{4}$ must have at least $\frac{3}{2}n^{\prime}$ vertices, we have $n-1+\widetilde{O}(m/\kappa)-(n-n^{\prime})\geq\frac{3}{2}n^{\prime}.$ Hence, $n^{\prime}\leq\widetilde{O}(m/\kappa).$ Thus $\mathcal{H}_{4}$ is an instance with at most $\widetilde{O}(m/\kappa)$ vertices and edges.

Thus, the total running time recurrence is

[TABLE]

Note that $\kappa$ is fixed throughout the recursion. By picking $\kappa=\widetilde{\Theta}(m^{\frac{1}{\sqrt{p-1}}}),$ we can fix the depth of the recursion to be $O(\sqrt{p-1}).$ The total cost is dominated by the cost at the bottom level of the recursion, which adds up to a total running time of $2^{O(p^{\nicefrac{{3}}{{2}}})}m^{1+O(\frac{1}{\sqrt{p-1}})}.$

The above discussion does not take into account the reduction in $\delta$ as we go down the recursion. Observe that $\delta$ is lower bounded by

[TABLE]

Thus, we always satisfy $\log\delta=\widetilde{O}(1).$ ∎

We can now prove the central collaries regardling smoothed $\ell_{p}$-norm flows and $\ell_{p}$-norm flows.

See 1.1

Proof.

First note that the Problem (2) is a smoothed $p$-norm instance after flipping the sign of $\boldsymbol{\mathit{g}},$ and the sign of the objective function. We can solve this smoothed $p$-norm instance $\mathcal{G}$ by using Theorem 3.7 to compute the desired approximate solution to the residual problems. We start with $\boldsymbol{\mathit{f}}^{(0)}$ as our initial solution. At iteration $t,$ we invoke Theorem 3.7 using $\boldsymbol{\mathit{f}}^{(t-1)}$ as the initial solution,

[TABLE]

where $\kappa$ is given by Theorem 3.7.

We know that $\boldsymbol{\mathit{f}}^{(t)}$ satisfies,

[TABLE]

Iterating $O(\log n)$ times, we obtain,

[TABLE]

Finally, noting that we had flipped the sign of the objective function in Problem (2), we obtain our claim. ∎

See 1.2

Proof.

The pseudocode for our procedure $\textsc{pFlows}(\mathcal{G},\boldsymbol{\mathit{b}})$ for this problem is given in Algorithm 2. Our goal is to compute a flow $\boldsymbol{\widetilde{f}}$ satisfying $\boldsymbol{\mathit{B}}^{G\top}\boldsymbol{\widetilde{f}}=\boldsymbol{\mathit{b}},$ such that

[TABLE]

where $\boldsymbol{\mathit{f}}^{\star}$ is the flow minimizing the $\ell_{p}$-norm with residue $\boldsymbol{\mathit{b}}$. For concreteness, we take this to mean \mathopen{}\mathclose{{}\left\|\boldsymbol{\widetilde{f}}}\right\|^{p}_{p}\leq\mathopen{}\mathclose{{}\left(1+3m^{-c}}\right)\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{\star}}\right\|^{p}_{p}, for some constant $c$. We construct a smoothed $p$-norm instance $\mathcal{G}=(V,E,\boldsymbol{\mathbf{0}},\boldsymbol{\mathbf{0}},1).$ Note that the smoothed $p$-norm flow problem $(\mathcal{G},\boldsymbol{\mathit{b}})$ finds a flow satisfying residues $\boldsymbol{\mathit{b}},$ and maximizing \mathcal{E}^{\mathcal{G}}(\boldsymbol{\mathit{f}})=-\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|^{p}_{p}. We can solve this smoothed $p$-norm instance by iteratively refining using Corollary 3.3, and using Theorem 3.7 to compute the desired approximate solution to the residual problems.

Formally, we use Laplacian solvers to compute in $\widetilde{O}(m)$ time $\boldsymbol{\mathit{f}}^{(0)}$ as a 2-approximation to \min_{\boldsymbol{\mathit{f}}:\boldsymbol{\mathit{B}}^{\top}\boldsymbol{\mathit{f}}=\boldsymbol{\mathit{b}}}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|. We have,

[TABLE]

At each iteration $t,$ we construct the residual smoothed $p$-norm circulation problem $\mathcal{H}_{t}$ for $(\mathcal{G},\boldsymbol{\mathit{b}})$ with the current solution $\boldsymbol{\mathit{f}}^{(t)},$ given by Corollary 3.3. We then invoke

[TABLE]

where $\kappa=\widetilde{\Theta}(m^{\frac{1}{\sqrt{p-1}}})$ is given by Theorem 3.7. Let $\boldsymbol{\mathit{\Delta}}^{(t)}$ be the flow returned. We know

[TABLE]

We let $\boldsymbol{\mathit{f}}^{(t+1)}\leftarrow\boldsymbol{\mathit{f}}^{(t)}+2^{-3p}\boldsymbol{\mathit{\Delta}}^{(t)}.$ Thus, by Corollary 3.3, at every iteration, we have

[TABLE]

where we have used, \mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{(0)}}\right\|_{p}\leq 2m^{\frac{1}{2}-\frac{1}{p}}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{\star}}\right\|_{p},

[TABLE]

Thus

[TABLE]

Where first step follows by rearranging terms. To establish the second inequality, first consider the case when \mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{(t)}}\right\|_{p}^{p}-\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{\star}}\right\|_{p}^{p}}\right)\geq 2m^{-c}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{\star}}\right\|_{p}^{p} and hence \mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{(t)}}\right\|_{p}^{p}-\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{\star}}\right\|_{p}^{p}\leq(1-2^{-4p}/2)\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{(t)}}\right\|_{p}^{p}-\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{\star}}\right\|_{p}^{p}}\right), meanwhile, when \mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{(t)}}\right\|_{p}^{p}-\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{\star}}\right\|_{p}^{p}}\right)<2m^{-c}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{\star}}\right\|_{p}^{p}, the inequality is immediate.

Iterating $T=\Theta((c+p)2^{4p}\log m)$ times gives us

[TABLE]

∎

4 Graph Theoretic Preconditioners

In this section, we discuss at a high level of the construction of ultra-sparsifiers for a smooth $\ell_{p}$-norm instance. We start by restating the main theorem of our ultra-sparsifier. After establishing the necessary tools, we prove this theorem at the end of this section.

See 3.6

Our high-level approach is the same as Spielman-Teng [ST14], where we utilize a low-stretch spanning tree, and move off-tree edges to a small set of portal nodes. Once most off-tree edges are only between a small set of portal nodes, we sparsify the graph over the portal nodes to reduce the number of edges. As we need to map flow solutions between the original instance and the sparsified instance, our main concern is to carry out these step without incurring too much error on the objective function values. In our case we have $\ell_{p}^{p}$ resistances and gradients on edges in addition to $\ell_{2}^{2}$ resistances, and thus the main challenge is to simultaneously preserve their respective terms.

4.1 Tree-Portal Routing

Faced with a sparse graph or dense graph, we wish to move most edges onto a few vertices of a tree, so that many of the remaining vertices are low degree and can be eliminated. Our high-level approach is the same as Spielman-Teng [ST14], where we utilize a low-stretch spanning tree, and move off-tree edges to a small set of portal nodes. However, when we reroute flow on the graph where most edges are moved to be among a few vertices using the tree, we need to (approximately) preserve three different properties of the flow on the original graph, namely the inner product between gradients and flows $\sum_{e}\boldsymbol{\mathit{g}}_{e}\boldsymbol{\mathit{f}}_{e}$, the 2-norm energy $\boldsymbol{\mathit{r}}_{e}\boldsymbol{\mathit{f}}_{e}^{2}$, and the $\ell_{p}$-norm energy $\sum_{e}\boldsymbol{\mathit{f}}_{e}^{p}$. It turns out that we can move edges around on our graphs to produce a new graph while exactly preserving the linear term $\sum_{e}\boldsymbol{\mathit{g}}_{e}\boldsymbol{\mathit{f}}_{e}$ for flows mapped between one graph and the other. This means any tree is acceptable from the point of preserving the linear term. To move edges around and bound distortion of solutions w.r.t. the quadratic $\boldsymbol{\mathit{r}}_{e}\boldsymbol{\mathit{f}}_{e}^{2}$ term, we use a low stretch tree w.r.t. the $\boldsymbol{\mathit{r}}$ weights as resistances. This leaves us with little flexibility for the $\sum_{e}\boldsymbol{\mathit{f}}_{e}^{p}$ term. However, for large $p$, provided every $p$-th order term is weighted the same (i.e. we have $s\sum_{e}\boldsymbol{\mathit{f}}_{e}^{p}$ instead of $\sum_{e}\boldsymbol{\mathit{s}}_{e}\boldsymbol{\mathit{f}}_{e}^{p}$), it turns out that, moving edges along any tree will result in bounded distortion of the solution, provided we are careful about how we move those edges. Thus, we can move edges around carefully to be among a small subset of the portal nodes while simultaneously controlling all linear, $2$-nd order and $p$-th order terms. But, this only works if all the $p$-th order terms are weighted the same. This leads us to maintain uniform-weighted $p$-th order terms as an invariant throughout the algorithm. The iterative refinement steps naturally weigh all $p$-th order terms the same provided the original function does. However, our sparsification procedures do not immediately achieve this, but we show we can enforce this uniform-weight invariant with only a manageable additional distortion of our solutions. Elimination also does not naturally weigh all $p$-th order terms the same even in our case when the original function does, but we can bound the distortion incurred by explicitly making the weights uniform. Our tree-based edge re-routing naturally creates maps between solutions on the old and new graphs.

We first formalize what we mean by moving off-tree edges. Suppose we have a spanning tree $T$ of a graph $(V,E)$ and a subset set of nodes ${\widehat{{V}}}\subset V$ designated as portal nodes, for any off-tree edge $e=\{u,v\}\in E\setminus T$, there is a unique tree path ${\mathcal{{P}}}_{T}(u,v)$ in $T$ from $u$ to $v$. We define a tree-portal path ${\mathcal{{P}}}_{T,{\widehat{{V}}}}(u,v)$, which is not necessarily a simple path.

Definition 4.1** (Tree-portal path and edge moving).**

Given spanning tree $T$ and set of portal nodes ${\widehat{{V}}}$, let $e=\{u,v\}$ be any edge not in $T$, and $P_{T}(u,v)$ the unique tree path in $T$ from $u$ to $v$. We define $e$’s tree-portal path $P_{T,{\widehat{{V}}}}(u,v)$ and $e$’s image under tree-portal edge moving as follows

If $P_{T}(u,v)$ doesn’t go through any portal vertex. In this case, we replace $\{u,v\}$ with a distinct self-loop of $v$. We let $P_{T,{\widehat{{V}}}}(u,v)$ be the path $P_{T}(u,v)$ followed by the self-loop at $v$. 2. 2.

If $P_{T}(u,v)$ goes through exactly one portal vertex ${\widehat{{v}}}$. In this case, we replace $\{u,v\}$ with a distinct self-loop at ${\widehat{{v}}}$. We let $P_{T,{\widehat{{V}}}}(u,v)$ be the tree path $P_{T}(u,{\widehat{{v}}})$ followed by the self-loop at ${\widehat{{v}}}$ and then the tree path $P_{T}({\widehat{{v}}},v)$. 3. 3.

If $P_{uv}$ goes through at least two portal vertices. In this case, let ${\widehat{{u}}}$ (and ${\widehat{{v}}}$) be closest the portal vertex to $u$ (and $v$) on $P_{uv}$, we replace $\{u,v\}$ with a distinct edge777We will keep multi-edges explicitly between portal nodes. $\{{\widehat{{u}}},{\widehat{{v}}}\}$. We let $P_{T,{\widehat{{V}}}}(u,v)$ be the tree path $P_{T}(u,{\widehat{{u}}})$ followed by the new edge from ${\widehat{{u}}}$ to ${\widehat{{v}}}$ and then the tree path $P_{T}({\widehat{{v}}},v)$.

This maps any off-tree edge $e$ to a unique (edge or self-loop) ${\widehat{{e}}}$ given any $T,{\widehat{{V}}}$. We denote the tree-portal edge moving with the map ${\widehat{{e}}}=\mbox{{Move}}_{T,{\widehat{{V}}}}(e)$.

Although we will get self-loops in tree-portal routing, to keep the discussion simple, we ignore the possibility of getting self-loops. This still captures all the main ideas, and the algorithm/analysis extends to self-loops in a very straightforward but slightly tedious way. We discuss self-loops briefly at the end of the section.

Tree-portal routing is a mapping from flow solutions on the original off-tree edges to a flow solution (with the same residue) using the edges they are mapped to. Any flow along off-tree edge $(u,v)$ in the original graph is rerouted (again from $u$ to $v$) using the tree-portal path $P_{T,{\widehat{{V}}}}(u,v)$ instead. Rerouting the flow of any off-tree edge along its tree-portal path increases the congestion on tree edges, which in turn incurs error in the $\ell^{2}_{2}$ and $\ell^{p}_{p}$ terms in the objective function. We need to pick the tree and portal nodes carefully to bound the error.

Definition 4.2**.**

Given any graph $(V,E)$, resistance $\boldsymbol{\mathit{r}}$ on edges, a spanning tree $T$, and a set of portals ${\widehat{{V}}}\subset V$, for any $e=\{u,v\}\in E$, let ${\widehat{{e}}}=\mbox{{Move}}_{T,{\widehat{{V}}}}(e)$ and $P_{T,{\widehat{{V}}}}(u,v)$ be as specified above. The stretch of $e=\{u,v\}\in E\setminus T$ in the tree-portal routing is

[TABLE]

and the stretch of a tree edge $e\in T$ is \mathsf{Str}_{T,{\widehat{{V}}}}\mathopen{}\mathclose{{}\left(e}\right)=1. Note with our definition $\mathsf{Str}_{T,\emptyset}(e)$ gives the standard stretch.

The starting point is low stretch spanning trees [AN12], which provide good bounds on the total $\ell_{2}^{2}$ stretch.

Lemma 4.3** (Low-Stretch Trees [AN12]).**

Given any graph $G=(V,E)$ of $m$ edges and $n$ nodes, as well as resistance $\boldsymbol{\mathit{r}}$, $\textsc{LSST}(\boldsymbol{\mathit{r}})$ finds a spanning tree in O\mathopen{}\mathclose{{}\left(m\log{n}\log\log{n}}\right) time such that

[TABLE]

We will construct a low-stretch spanning tree $T$ of $(V^{\mathcal{G}},E^{\mathcal{G}},\boldsymbol{\mathit{r}}^{\mathcal{G}})$ using the above result. Still, the error will be too large if we only use tree edges to reroute the flow of all the off-tree edges, since the low average stretch doesn’t prevent one tree edge to be on the tree path for many off-tree edges. Thus, we need to add portal nodes so we can shortcut between them to reduce the extra congestion on tree edges.

4.2 Partitioning Trees into Subtrees and Portals

Next, we show how to find a small set of good portal nodes so that rerouting flow on off-tree edges using their tree-portal paths incurs small error in the objective function. Pseudocode of this routine is in Algorithm 3, and its guarantees are stated in Lemma 4.4 below.

Lemma 4.4**.**

There is a linear-time routine FindPortals that given any graph $G$, a spanning tree $T$, with ${\widehat{{m}}}$ off-tree edges, and a portal count ${\widehat{{n}}}\leq{\widehat{{m}}}$, returns a subset of ${\widehat{{V}}}$ of ${\widehat{{n}}}$ vertices so that for all edges ${\widehat{{e}}}\in T$, we have

[TABLE]

This lemma will be a fairly straightforward using the tree decomposition subroutine (page $881$ of journal version) from Spielman and Teng [ST14], which we include below for completeness.

Definition 4.5** ([ST14] Definition $10.2$).**

Given a tree $T$ that spans a set of vertices $V$, a $T$-decomposition is a decomposition of $V$ into sets $W_{1},\ldots,W_{h}$ such that $V=\bigcup W_{i}$, the graph induced by $T$ on each $W_{i}$ is a tree, possibly with just one vertex, and for all $i\neq j$, $|W_{i}\bigcap W_{j}|\leq 1$.

Given an additional set of edges $E$ on $V$, a $(T,E)$-decomposition is a pair $(\{W_{1},\ldots,W_{h}\},\rho)$ where $\{W_{1},\ldots,W_{h}\}$ is a $T$-decomposition and $\rho$ is a map that sends each edge of $E$ to a set or pair of sets in $\{W_{1},\ldots,W_{h}\}$ so that for each edge in $(u,v)\in E$,

if $\rho(u,v)=\{W_{i}\}$, then $\{u,v\}\subset W_{i}$, and 2. 2.

if $\rho(u,v)=\{W_{i},W_{j}\}$, then either $u\in W_{i}$ and $v\in W_{j}$ or vice versa.

Theorem 4.6** ([ST14] Theorem $10.3$).**

There exists a linear-time algorithm such that on input a set of edges $E$ on vertex set $V$, a spanning tree $T$ on $V$, a function $\eta:E\rightarrow\mathbb{R}^{+}$, and an integer $1<t\leq\sum_{e\in E}\eta(e)$, outputs a $(T,E)$-decomposition $(\{W_{1},\ldots,W_{h}\},\rho)$, such that

$h\leq t$** 2. 2.

for all $W_{i}$ such that $|W_{i}|>1$,

[TABLE]

We can use the above theorem to show Lemma 4.4.

Proof of Lemma 4.4.

We will apply Theorem 4.6 with $t={\widehat{{n}}}$, and the function $\eta$ will be

[TABLE]

Note by construction $\sum_{e}\eta(e)\leq 2\sum_{e^{\prime}\in E}\mathsf{Str}_{T,\emptyset}(e^{\prime})$. We get $\{W_{1},\ldots,W_{{\widehat{{n}}}}\}$ back, and let $T_{i}$ be the tree induced by $T$ on $W_{i}$. Note the $T_{i}$’s will be edge disjoint, and cover all tree edges of $T$. Our set of portals will be the set of nodes that are in more than one of the $W_{i}$’s, i.e. the nodes where different $T_{i}$’s overlap. Note the number of portals is at most the number of $T_{i}$’s by an inductive argument from any $T_{i}$ that is a sub-tree in $T$. Such $T_{i}$ have exactly one portal, and we can remove $T_{i}$ from $T$ and continue the argument until all that remain in $T$ is one sub-tree.

Consider any tree edge ${\widehat{{e}}}\in T$, suppose it in $T_{i}$ for some $i$. ${\widehat{{e}}}$ can only be on the tree-portal routing for some edge $\{u,v\}$ when $W_{i}\in\rho(u,v)$. Note as $T_{i}$ contains at least one tree edge, we know $|W_{i}|>1$, the second guarantee in Theorem 4.6 gives

[TABLE]

which directly gives the bounds we want in the lemma. ∎

4.3 Graph Sparsification

Once we are able to move most of the edges onto a small subset of vertices, we wish to sparsify the resulting dense graph over those vertices. This sparsification has to simultaneously preserve properties of $1$-st, $2$-nd and $p$-th order terms, as well as the interactions between them, which turns out to be challenging. We resort to expander decomposition which allows us to partition the vertex set s.t. the edges internal to each subset form an expander and not too many edges cross the partitions. Just having an expander graph is not enough to allow us to sample the graph due to the need of preserving the linear terms. Thus, we also require that on each expander the orthogonal projection of the gradient to the cycle space of the sub-graph has its maximum squared entry not much larger than the the average squared entry. We refer to this as a uniform (projected) gradient. We discuss how to obtain an expander decomposition that guarantees the projected gradients are uniform in the expanders later in this overview. Given the uniform projected gradient condition, we show that we can uniformly sample edges of these expanders to create sparsified versions of them. We construct maps between the flows on an original expander and its sampled version that work for any flow, not only a circulation. These maps preserve the linear term $\sum_{e}\boldsymbol{\mathit{g}}_{e}\boldsymbol{\mathit{f}}_{e}$ exactly, while bounding the cost of the $2$-norm and $\ell_{p}$-norm terms by relating them to the cost of the optimal routing of a flow with the same demands and same gradient inner product, and showing that optimal solutions are similar on the original expander and its sampled version. This strategy resembles the flow maps developed in [KLOS14], and like their maps, we route demands using electrical flows on individual expanders, but additionally we need create a flow in the cycle space that depends on projection of the gradient onto the cycle space.

Tree-portal routing will give us an instance where all the off-tree edges are between portal nodes. We can look at the sub-graph restricted to the portal nodes and the off-tree edges between them. This graph has many fewer nodes comparing to the original graph but roughly the same number of edges, and thus is much denser. We can then sparsify this graph to reduce the number of off-tree edges similar to the construction of spectral sparsifiers. The main technical difficulty is that in the sparsified graph, we still want the $\ell_{p}^{p}$ terms in our objective function to have a same scalar $s$ for every edge, but similar to the case of how resistances are scaled in spectral sparsification, to preserve the total value of the $\ell_{p}^{p}$ terms, we would naturally want to scale $s$ according to the probability we sample an edge $e$. Thus, to get a same scalar $s$ for all sampled edges, we are limited to uniform sampling. We know uniform sampling works in expanders (c.f. [ST14, SS11] and [KLOS14, She13] for $\ell_{2}$ and $\ell_{\infty}$ respectively), so the natural approach is to first decompose the graph into expanders, and sampling uniformly inside each expander. However, because of the presence of a gradient, we need to be a bit more careful than even expanderdecomposition-based sparsification steps. Thus, we work with uniform expanders.

Definition 4.7**.**

A graph 888We use an instance and its underlying graph interchangeably in our discussion. $G$ is a $\alpha$-uniform $\phi$-expander (or uniform expander when parameters not spelled out explicitly) if

$\boldsymbol{\mathit{r}}$ on all edges are the same. 2. 2.

$\boldsymbol{\mathit{s}}$ on all edges are the same. 3. 3.

$G$ has conductance999$\boldsymbol{\mathit{r}}$ are uniform, so conductance is defined as in unweighted graphs. We use the standard definition of conductance. For graph $G=(V,E)$, the conductance of any $\emptyset\neq S\subsetneq V$ is \phi(S)=\frac{\delta(S)}{\min\mathopen{}\mathclose{{}\left(vol(S),vol(V\setminus S)}\right)} where $\delta(S)$ is the number of edges on the cut $(S,V\setminus S)$ and $vol(S)$ is the sum of the degree of nodes in $S$. The conductance of a graph is $\phi_{G}=\min_{S\neq\emptyset,V}\phi(S)$. at least $\phi$. 4. 4.

The projection of $\boldsymbol{\mathit{g}}$ onto the cycle space of $G$, $\boldsymbol{\widehat{\mathit{g}}}^{G}=(I-\boldsymbol{\mathit{B}}\boldsymbol{\mathit{L}}^{{\dagger}}\boldsymbol{\mathit{B}}^{\top})\boldsymbol{\mathit{g}}$, is $\alpha$-uniform (see next definition), where $\boldsymbol{\mathit{B}}$ is the edge-vertex incidence matrix of $G$, and $\boldsymbol{\mathit{L}}=\boldsymbol{\mathit{B}}^{\top}\boldsymbol{\mathit{B}}$ is the Laplacian.

Definition 4.8**.**

A vector $\boldsymbol{\mathit{y}}\in\mathbb{R}^{m}$ is said to be $\alpha$-uniform if

[TABLE]

We abuse the notation to also let the all zero vector $\boldsymbol{\mathbf{0}}$ be $1$-uniform.

In Section 5 we show how to decompose the graph consisting of portals and the off-tree edges between them into vertex disjoint uniform expanders such that more than half of the edges are inside the expanders.101010Some of the expanders we find actually won’t satisfy the projected gradient being $\alpha$-uniform constraint (case $3(b)$ in Theorem 4.9). For those expanders, the projection of the gradient in the cycle space is tiny so we make it [math]. This leads to the additive error in Theorem 3.6.

Theorem 4.9** (Decomposition into Uniform Expanders).**

Given any graph/gradient/resistance instance $\mathcal{G}$ with $n$ vertices, $m$ edges, unit resistances, and gradient $\boldsymbol{\mathit{g}}^{\mathcal{G}}$, along with a parameter $\delta$, $\textsc{Decompose}(\mathcal{G},\delta)$ returns vertex disjoint subgraphs $\mathcal{G}_{1},\mathcal{G}_{2},\ldots$ in $O(m\log^{7}n\log^{2}(n/\delta))$ time such that at least $m/2$ edges are contained in these subgraphs, and each $\mathcal{G}_{i}$ satisfies (for some absolute constant $c_{partition}$):

The graph $(V^{\mathcal{G}_{i}},E^{\mathcal{G}_{i}})$ has conductance at least

[TABLE]

and degrees at least $\phi\cdot\frac{m}{3n}$, where $c_{partition}$ is an absolute constant. 2. 2.

The projection of its gradient $\boldsymbol{\mathit{g}}^{\mathcal{G}_{i}}$ into the cycle space of $\mathcal{G}_{i}$, $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}_{i}}$ satisfies one of:

1. (a)

$\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}_{i}}$* is $O(\log^{8}{n}\log^{3}(n/\delta))$-uniform,*

[TABLE]

Here $m_{i}$ is the number of edges in $\mathcal{G}^{\mathcal{G}_{i}}$. 2. (b)

The $\ell_{2}^{2}$ norm of $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}_{i}}$ is smaller by a factor of $\delta$ than the unprojected gradient:

[TABLE]

Moreover, the min degree of any node in the expanders is up to a polylog factor close to the average degree. For the off-tree edges not included in these uniform expanders, we work on the pre-image111111By pre-image of ${\widehat{{e}}}$ we mean the original off-tree edge $e$ that gets moved to ${\widehat{{e}}}$ in the tree-portal routing, i.e. $e=\textsc{Move}^{-1}_{T,{\widehat{{V}}}}({\widehat{{e}}})$. of them in the next iteration. That is, for any edge ${\widehat{{e}}}$ inside one of the expanders, we remove its pre-image from the instance $\mathcal{G}$, and work on the remaining off-tree edges in $\mathcal{G}$ in the next iteration. This iterative process terminates when the number of remaining off-tree edges is small enough (i.e. $\widetilde{O}(|E^{\mathcal{G}}|/\kappa)$). This takes $O(\log|E^{\mathcal{G}}|)$ iterations as a constant fraction of off-tree edges are moved to be inside the expanders each iteration.

Sparsify Uniform Expanders

If we append a column containing the gradient of edges to the edge-vertex incidence matrix $\boldsymbol{\mathit{B}}$, the conditions of a $\alpha$-uniform $\phi$-expander is equivalent to each row of $\boldsymbol{\mathit{B}}$ having leverage score at most $\frac{n\alpha\phi^{-1}}{m}$ where $n,m$ are number of nodes and edges. An underlying connection with the $\ell_{p}$-norm row sampling result by Cohen and Peng [CP15] is that this is also a setting under which $\ell_{q}$-norm functionals are preserved under uniform sampling. We refrain from developing a more complete picture of such machinery here, and will utilize ideas closer to routing on expanders [KM09, KLOS14] to show a cruder approximation in Section D.

Theorem 4.10** (Sampling Uniform Expanders).**

Given an $\alpha$-uniform $\phi$-expander $\mathcal{G}=(V^{\mathcal{G}}E^{\mathcal{G}},r^{\mathcal{G}},s^{\mathcal{G}},\boldsymbol{\mathit{g}}^{\mathcal{G}})$ with $m$ edges and vertex degrees at least $d_{\min}$, for any sampling probability $\tau$ satisfying

[TABLE]

where $c_{sample}$ is some absolute constant, $\textsc{SampleAndFixGradient}(\mathcal{G},\tau)$ w.h.p. returns a partial instance $\mathcal{H}=(H,r^{\mathcal{H}},s^{\mathcal{H}},\boldsymbol{\mathit{g}}^{\mathcal{H}})$ and maps $\mathcal{M}_{\mathcal{G}\rightarrow\mathcal{H}}$ and $\mathcal{M}_{\mathcal{H}\rightarrow\mathcal{G}}$. The graph $H$ has the same vertex set as $G$, and $H$ has at most $2\tau m$ edges. Furthermore, $r^{\mathcal{H}}=\tau\cdot r^{\mathcal{G}}$ and $s^{\mathcal{H}}=\tau^{p}\cdot s^{\mathcal{G}}$. The maps $\mathcal{M}_{\mathcal{G}\rightarrow\mathcal{H}}$ and $\mathcal{M}_{\mathcal{H}\rightarrow\mathcal{G}}$ certify

[TABLE]

*where $\kappa=m^{1/(p-1)}\phi^{-9}\log^{3}n$. *

4.4 Ultra-sparsification Algorithm and Error Analysis

Now we put all the pieces together. We need to show that adding together our individual sparsifiers results in a sparsifier of the overall graph. This is fairly immediate given the strong guarantees we established on the individual graphs. We also need to be able to repeatedly decompose and sparsify enough times that the overall graph becomes sparse. To address this issue, we use ideas from [KMP11] that suggest scaling up the tree from the tree routing section limits the error incurred during sampling. Here it again becomes important that because we rely on [SW18], we know exactly which edges belong to a sparsifier. This guarantee limits the interaction between sparsification of different expanders.

After constructing a low-stretch spanning tree $T$, we round each $\boldsymbol{\mathit{r}}^{\mathcal{G}}_{e}$ of off-tree edges $e\in E^{\mathcal{G}}\setminus T$ down to the nearest power of $2$ (can be less than $1$) if $\boldsymbol{\mathit{r}}^{\mathcal{G}}_{e}\geq\delta$, and round to [math] otherwise. This gives of \log\frac{\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{r}}^{\mathcal{G}}}\right\|_{\infty}}{\delta} bucket of edges with uniform resistances, and we work with one bucket of edges at a time, since the edges in a uniform expander need to have uniform $\boldsymbol{\mathit{r}}_{e}$. If $\mathcal{G}^{\prime}$ is the instance of $\mathcal{G}$ after rounding the resistance of off-tree edges, it is easy to see the following error guarantee.

Lemma 4.11**.**

$\mathcal{G}\preceq_{1}\mathcal{G}^{\prime}$* with the identity mapping, , and for any flow solution $\boldsymbol{\mathit{f}}^{\mathcal{G}^{\prime}}$ of $\mathcal{G}^{\prime}$, again using the identity mapping, we have*

[TABLE]

To avoid using too many symbols, we reuse $\mathcal{G}$ to refer to the original instance after the resistance rounding (i.e. the $\mathcal{G}^{\prime}$ above). Denote $E^{r}$ the subset of edges in $E^{\mathcal{G}}\setminus T$ containing edges with $\boldsymbol{\mathit{r}}_{e}=r$ for some particular $r$, note there are at most \log\frac{\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{r}}^{\mathcal{G}}}\right\|_{\infty}}{\delta} possible value of $r$. We work iteratively on the set $E^{r}$, starting with $E^{r}_{0}=E^{r}$. In the $i$-th iteration, we use $\textsc{FindPortal}(T,E^{r}_{i},m/\kappa)$ (Lemma 4.4) to find a set of $m/\kappa$ portal vertices for the edges remaining in $E^{r}_{i}$, note the low-stretch spanning tree is fixed through the process, but each iteration we find a new set of portals using FinalPortal as introduced in Section 4.2.

We then move edges in $E^{r}_{i}$ using the tree-portal routing. We let ${\widehat{{G}}}^{r}_{i}$ to be the graph of the $m/\kappa$ portal nodes and the off-tree edges between them. Note the number of edges in ${\widehat{{G}}}^{r}_{i}$ is $|E^{r}_{i}|$ and the number of nodes is $m/\kappa$.

So far we haven’t specified the $\boldsymbol{\mathit{r}},s,\boldsymbol{\mathit{g}}$ values on the edges in ${\widehat{{G}}}^{r}_{i}$, and these values will depend on the tree-portal routing as well as the average degree in ${\widehat{{G}}}^{r}_{i}$. For now we focus on discuss the edge set in our final sparsified instance, and assume we have $\boldsymbol{\mathit{r}},s,\boldsymbol{\mathit{g}}$ values for ${\widehat{{G}}}^{r}_{i}$. We will come back to specify these quantities later.

We use Decompose (Theorem 4.9) on the graph ${\widehat{{G}}}^{r}_{i}$ to compute a collection of vertex disjoint sub-graphs $\{{\widehat{{G}}}^{r}_{i,1},{\widehat{{G}}}^{r}_{i,2},\ldots\}$, and at least half of the edges in ${\widehat{{G}}}^{r}_{i}$ are inside these sub-graphs. We let ${\widehat{{E}}}^{r}_{i}$ to be the edges contained in these sub-graphs, and ${\widetilde{{E}}}^{r}_{i}$ be the set of pre-images of edges in ${\widehat{{E}}}^{r}_{i}$ (in terms of tree-portal off-tree edge moving). We remove ${\widetilde{{E}}}^{r}_{i}$ from $E^{r}_{i}$ and proceed to iteration $i+1$. If at the beginning of some iteration $i$, the size of $E^{r}_{i}$ is at most $\widetilde{O}(m/\kappa)$, we leave them as off-tree edges, and denote $E^{r}_{last}$ as the set containing them. Note for any $r$, the iterative process must finish in $O(\log\kappa)$ iterations as we start with $|E^{r}|<m$ edges. We do this for all $r$.

So far any edge in the original instance $\mathcal{G}$ we get either $(1)$ a tree edge in $T$, or $(2)$ an off-tree edge in a ${\widehat{{G}}}^{r}_{i,j}$ for some resistance value $r$, iteration $i$, and $j$-th expander computed in that iteration, or $(3)$ an off-tree edge remaining in $\mbox{{\bf E}}^{r}_{last}$ for some resistance value $r$. There are $n-1$ edges in group $(1)$, and $\widetilde{O}(m/\kappa)$ edges in group $(3)$, so we can keep all these edges in the ultra-sparsifier $\mathcal{H}$. For the off-tree edges in group $(2)$, we uniformly sample the edges in each ${\widehat{{G}}}^{r}_{i,j}$ to get a sparsified graph ${\overline{{H}}}^{r}_{i,j}$. Technically our sampling result only applies to an $\alpha$-uniform $\phi$-expander ${\widehat{{G}}}^{r}_{i,j}$ (i.e. case $3(a)$ in Theorem 4.9). If the ${\widehat{{G}}}^{r}_{i,j}$ we get back from Decompose is in case $3(b)$ of, we perturb the gradient on edges so that the projection of the gradient to the cycle space of the expander is [math], i.e. project the gradient to the space orthogonal to the cycle space. Then we have ${\widehat{{G}}}^{r}_{i,j}$ after perturbation is a $1$-uniform $\phi$-expander.

The edges in our final ultra-sparsifier $\mathcal{H}$ will be the tree edges in $T$, the off-tree edges in the $\mbox{{\bf E}}^{r}_{last}$’s over all resistance bucket value $r$, and the off-tree edges in the ${\overline{{H}}}^{r}_{i,j}$’s over all resistance $r$, iteration $i$ and expander $j$. We argued about the size of all but the edges in the ${\overline{{H}}}^{r}_{i,j}$’s, which we will do now.

Sampling Probability

We first specify the probability we sample each edge in ${\widehat{{G}}}^{r}_{i,j}$ to get ${\overline{{H}}}^{r}_{i,j}$ which we denote by $\tau_{r,i}$ (same across all the expanders, i.e. $j$’s, for any resistance $r$ and iteration $i$). By Theorem 4.10 we need the probability to be at least $c_{sample}\log{n}(d_{\min}^{-1}\phi^{-2}+\alpha m^{-1})$. Here $c_{sample}$ is a fixed constant across all $r,i,j$’s, and the guarantees on ${\widehat{{G}}}^{r}_{i,j}$ from Theorem 4.9 allow us to use some fixed polylog$n$ as $\phi^{-2}$ and $\alpha$ across all $r,i,j$’s. The only parameter that varies across different $r,i$’s is $d_{\min}$, a lower bound on the minimum vertex degree in ${\widehat{{G}}}^{r}_{i,j}$, which by Theorem 4.9 is within a (fixed) polylog factor of the average degree in ${\widehat{{G}}}^{r}_{i}$. As there are $m_{r,i}$ edges and $m/\kappa$ nodes in ${\widehat{{G}}}^{r}_{i}$, the average degree is $m_{r,i}\kappa/m$. Thus, we can write $\tau_{r,i}=\frac{c_{1}m\log^{c_{2}}n}{\kappa m_{r,i}}$ for some global constants $c_{1},c_{2}$, and since both $m_{r,i}$ and $\kappa$ is at most $m$, $\tau_{r,i}$ satisfies the requirement on $\tau$ in Theorem 4.10. With this particular choice of $\tau_{r,i}$ we can use SampleAndFixGradient to sample ${\widehat{{G}}}^{r}_{i,j}$ and the the guarantees from Theorem 4.10. Now we can prove the statement about the number of off-tree edges in $\mathcal{H}$.

Lemma 4.12**.**

The total number of edges over all ${\overline{{H}}}^{r}_{i,j}$’s is $\widetilde{O}(\frac{m}{\kappa})$ with high probability.

Proof.

Pick any $r,i$, recall when we call Decompose in that iteration, we have ${\widehat{{G}}}^{r}_{i}$ with uniform $\boldsymbol{\mathit{r}}$, $m_{r,i}=|E^{r}_{i}|$ edges and $n_{i}=m/\kappa$ nodes. From the previous discussion of the sampling probability, we know it is sufficient to call SampleAndFixGradient on ${\widehat{{G}}}^{r}_{i,j}$ with probability

[TABLE]

for some constants $c_{1},c_{2}$. By Theorem 4.10, the number of edges in ${\overline{{H}}}^{r}_{i,j}$’s over all $j$ is at most $\widetilde{\Theta}(\frac{m}{\kappa})$ with high probability since over all $j$ the ${\widehat{{G}}}^{r}_{i,j}$’s contain $\Theta(m_{r,i})$ edges.

Since for each $r$ the number of iterations is at most $i\leq\log\kappa$, and there are \log\frac{\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{r}}^{\mathcal{G}}}\right\|_{\infty}}{\delta} possible $r$ values, the final bound in the lemma follows from summing over all $r,i$. Note we can hide all $\log$ factors as $\log n$ factors by our assumption in Theorem 3.6 that \log\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{r}}^{\mathcal{G}}}\right\|_{\infty} and $\log\frac{1}{\delta}$ are both polylog in $n$. ∎

Now we discuss the $\boldsymbol{\mathit{g}},\boldsymbol{\mathit{r}},s$ values we put on the edges in all the steps. Note we need the final instance $\mathcal{H}$ to have a uniform scalar $s^{\mathcal{H}}$ for every $|f_{e}|^{p}$ term, so we can recursively optimize the instance. However, in the intermediate steps, we will divide the instance into sub-instances induced by the different subsets of edges, e.g. ${\widehat{{G}}}^{r}_{i,j}$’s, and later combine sub-instances induced by the sampled sub-graphs ${\overline{{H}}}^{r}_{i,j}$’s to get $\mathcal{H}$. Each of these sub-instances will have its own scalar, e.g. $s^{G}_{r,i},s^{H}_{r,i}$, but in general they won’t necessarily have the same value across different sub-instances. Notation-wise, in the following discussion, we assume each edge has its own scalar $\boldsymbol{\mathit{s}}_{e}$ associated with the term $|f_{e}|^{p}$ in the intermediate instances. Eventually, the different scaling we do to $s$ in the intermediate steps will cancel so that in $\mathcal{H}$ we have the scalar $s^{\mathcal{H}}$. The input $\mathcal{G}$ has a uniform scalar $s^{\mathcal{G}}$, and we will make $s^{\mathcal{H}}=s^{\mathcal{G}}$.

Now we specify the $\boldsymbol{\mathit{g}}$, $\boldsymbol{\mathit{r}}$ and $s$ values of the edges in the final instance $\mathcal{H}$ as well as in some of the key intermediate sub-instances we consider.

$e\in T$: The gradient, resistance and $\boldsymbol{\mathit{s}}_{e}$ on these edges in $\mathcal{H}$ remain the same as in $\mathcal{G}$, that is $\boldsymbol{\mathit{g}}^{\mathcal{H}}_{e}=\boldsymbol{\mathit{g}}^{\mathcal{G}}_{e}$, $\boldsymbol{\mathit{r}}^{\mathcal{H}}_{e}=\boldsymbol{\mathit{r}}^{\mathcal{G}}_{e}$, and $\boldsymbol{\mathit{s}}_{e}=s^{\mathcal{G}}$. 2. 2.

$e\in E^{r}_{last}$: These off-tree edges remain at the end for each bucket $E^{r}$. We keep their gradient, resistance, and $\boldsymbol{\mathit{s}}_{e}=s^{\mathcal{G}}$ as in the original instance. 3. 3.

${\widehat{{e}}}\in{\widehat{{G}}}^{r}_{i,j}$, the $j$-th expander computed in iteration $i$ for resistance $r$: In the intermediate sub-instance induced by ${\widehat{{G}}}^{r}_{i,j}$, we have $\boldsymbol{\mathit{r}}_{{\widehat{{e}}}}=r\kappa\log^{2}n$, $\boldsymbol{\mathit{s}}_{{\widehat{{e}}}}=\tau_{r,i}^{-p}s^{\mathcal{G}}$. For the gradient on ${\widehat{{e}}}$, recall ${\widehat{{e}}}$ is the image of some off-tree edge $e$ under the mapping $\textsc{Move}_{T,{\widehat{{V}}}^{r}_{i}}$ where ${\widehat{{V}}}^{r}_{i}$ is the set of portals in the $i$-th iteration for resistance $r$. Under the tree-portal routing, any flow along $e=(u,v)$ will be rerouted along the tree-portal path $P_{T,{\widehat{{V}}}^{r}_{i}}(u,v)$. We want the linear term (i.e. gradient times flow) in the objective function to remain the same under this rerouting, so routing $1$ unit of flow from $u$ to $v$ along $P_{T,{\widehat{{V}}}^{r}_{i}}(u,v)$ should give the same dot product with the gradients as routing $1$ unit of flow from $u$ to $v$ along $e$ in the original instance (i.e. $\boldsymbol{\mathit{g}}^{\mathcal{G}}_{e}$). As the only off-tree edge on the tree-portal path is ${\widehat{{e}}}$, and we are keeping the original gradients on all the tree edges, this uniquely determines $\boldsymbol{\mathit{g}}_{{\widehat{{e}}}}$. 4. 4.

${\overline{{e}}}\in{\overline{{H}}}^{r}_{i,j}$ for some $r,i,j$: As specified in Theorem 4.10, if the edge ${\overline{{e}}}$ is sampled (with uniform probability $\tau_{r,i}$), and $\boldsymbol{\mathit{r}}_{{\overline{{e}}}}$,$\boldsymbol{\mathit{s}}_{{\overline{{e}}}}$ are their corresponding values in ${\widehat{{G}}}^{r}_{i,j}$ scaled up by $\tau_{r,i}$ and $\tau_{r,i}^{p}$ respectively. In particular we get back $\boldsymbol{\mathit{s}}_{{\overline{{e}}}}=s^{\mathcal{G}}$ as the $\tau_{r,i}^{p}$ scaling cancels the $\tau_{r,i}^{-p}$ scaling in ${\widehat{{G}}}^{r}_{i,j}$.

Note all the edges in our $\mathcal{H}$ (i.e. group $1,2,4$ above) end up with the same scalar $s^{\mathcal{H}}=s^{\mathcal{G}}$.

Now we bound the approximation error. For simplicity, we carry out the analysis ignoring the additive errors in the bound, and defer the discussion of them to the end. In particular, additive errors come in at two cases. The first is when we round an original resistance to [math] when it is less than $\delta$, and the second is in Decompose, we may get an expander whose projected gradient is not $\widetilde{O}(1)$-uniform but has tiny norm (i.e. case $3(b)$), and we zero out its projection to the cycle space before sampling. For now we assume we don’t have these cases.

We summarize the notations in our algorithm and analysis in Table 1. We explicitly point out whenever we change the gradient, resistance or $s$ value on an edge. We will use instances and their underlying graphs interchangeably, and when we refer to a subgraph as an instance, it will be clear what are the $\boldsymbol{\mathit{g}},\boldsymbol{\mathit{r}},s$ values for the instance.

First we let ${\overline{{G}}}$ be the instance on the same nodes and edges as $\mathcal{G}$, but for any $e\in{\widetilde{{E}}}^{r}_{i}$ (i.e. $e$ will be mapped to some ${\widehat{{e}}}$ in ${\widehat{{G}}}^{r}_{i,j}$), we rescale the resistance and $s$ to be $\boldsymbol{\mathit{r}}_{e}=r\kappa\log^{2}n$, and $\boldsymbol{\mathit{s}}_{e}=\tau_{r,i}^{-p}s^{\mathcal{G}}$ . Note the gradient of $e$ in ${\overline{{G}}}$ stays the same as in $\mathcal{G}$. We first bound the approximation error between $\mathcal{G}$ and this rescaled instance ${\overline{{G}}}$.

Lemma 4.13**.**

$\mathcal{G}\preceq_{\widetilde{O}(m^{1/(p-1)}\kappa)}{\overline{{G}}}\preceq_{1}\mathcal{G}$* with the identity mapping in both directions.*

Proof.

For any edge $e$, we have $\boldsymbol{\mathit{g}}^{{\overline{{G}}}}_{e}=\boldsymbol{\mathit{g}}^{\mathcal{G}}_{e}$. As to the $\ell^{p}_{p}$ scalar, we have either $\boldsymbol{\mathit{s}}^{{\overline{{G}}}}_{e}=\boldsymbol{\mathit{s}}^{\mathcal{G}}_{e}$, or if $e$ is eventually moved to some ${\widehat{{G}}}^{r}_{i,j}$ then

[TABLE]

as $m_{i,r}\geq\widetilde{O}{m/\kappa}$ or otherwise we would have stopped for resistance value $r$, we can assume $m_{r,i}\kappa\geq c_{1}m\log^{c_{2}}n$ so

[TABLE]

where the second inequality is by $m_{r,i}\leq m$, and the third inequality is by $\kappa<m$. Similar calculation gives $\boldsymbol{\mathit{r}}^{\mathcal{G}}_{e}\leq\boldsymbol{\mathit{r}}^{{\overline{{G}}}}_{e}\leq\kappa\log^{2}n\cdot\boldsymbol{\mathit{r}}^{\mathcal{G}}_{e}$. Our result directly follow by Lemma 2.13. ∎

Now we break ${\overline{{G}}}$ into sub-instances induced on the disjoint edge sets. Let ${\overline{{E}}}^{r}_{i}$ be the instance of ${\overline{{G}}}$ restricted to edges in ${\widetilde{{E}}}^{r}_{i}$, $T$ the instance restricted to the tree edges, and ${\overline{{G}}}_{rest}$ the instance restricted to edges in any of the $E^{r}_{last}$’s. When use addition as union on sets when the sets are disjoint. The objective of the sum of two instances is simply the sum of the individual instances objectives.

Lemma 4.14**.**

For any resistance value $r$, round $i$, we have

[TABLE]

where the flow mapping is the tree-portal routing and its reverse.

Proof.

Fix any resistance value $r$ and iteration $i$, the set of remaining off-tree edges of resistance $r$ in iteration $i$ is $E^{r}_{i}$, and these edges have a total stretch at most $O(m\log n\log\log n)$ with $T$ by Lemma 4.3, and $E^{r}_{i}=m_{r,i}$. As we use FindPortal to get a set of $m/\kappa$ portal nodes ${\widehat{{V}}}$ in that iteration, by Lemma 4.4, for any edge $e^{\prime}$ on $T$, we have in $T+E^{r}_{i}$

[TABLE]

and

[TABLE]

We first look at the direction from from $T+{\overline{{E}}}^{r}_{i}$ to $T+{\widehat{{E}}}^{r}_{i}$. Let $\boldsymbol{\mathit{f}}$ be the flow in $T+{\overline{{E}}}^{r}_{i}$, and $\boldsymbol{\widehat{\mathit{f}}}$ be the tree-portal routing of $\boldsymbol{\mathit{f}}$. In the tree-portal routing, flow on tree edges is mapped to the same flow, while any flow along an off-tree edge ${\overline{{e}}}=(u,v)\in{\overline{{E}}}^{r}_{i}$ is rerouted along the tree-portal path $P_{T,{\widehat{{V}}}}(u,v)$. This rerouting clearly preserves the residue between $\boldsymbol{\mathit{f}},\boldsymbol{\widehat{\mathit{f}}}$, and if ${\widehat{{e}}}\in{\widehat{{E}}}^{r}_{i}$ is the image of $(u,v)$, its gradient $\boldsymbol{\mathit{g}}_{{\widehat{{e}}}}$ in ${\widehat{{E}}}^{r}_{i}$ is by construction set to be the value which preserves the linear term in the objective function for $\boldsymbol{\mathit{f}}$ and $\boldsymbol{\widehat{\mathit{f}}}$. The cost of $\ell^{2}_{2}$ and $\ell^{p}_{p}$ terms for $f_{{\overline{{e}}}}$ is the same as the corresponding costs for $\boldsymbol{\widehat{\mathit{f}}}_{{\widehat{{e}}}}$, since ${\widehat{{e}}}$ is only used for the rerouting of ${\overline{{e}}}$ (so $|\boldsymbol{\widehat{\mathit{f}}}_{{\widehat{{e}}}}|=|f_{{\overline{{e}}}}|$), and they have the same $\boldsymbol{\mathit{r}},\boldsymbol{\mathit{s}}$ values. Thus, the contribution to the $\ell^{2}_{2},\ell^{p}_{p}$ terms in objective function from the off-tree edges are the same for $\boldsymbol{\mathit{f}}$ and $\boldsymbol{\widehat{\mathit{f}}}$. The only extra cost comes from the $\ell^{2}_{2}$ and $\ell^{p}_{p}$ terms of tree edges for $\boldsymbol{\widehat{\mathit{f}}}$ since we put additional flow through them. First consider the sum of the $\ell^{p}_{p}$ terms over all tree edges for $\boldsymbol{\widehat{\mathit{f}}}$ in $T+{\widehat{{E}}}^{r}_{i}$. Recall we don’t scale the $s$ value for tree edges, so the scalar is still $s^{\mathcal{G}}$ on tree edges, while for off-tree edges in ${\overline{{E}}}^{r}_{i}$, the value $s$ is scaled to be \mathopen{}\mathclose{{}\left(\frac{\kappa m_{r,i}}{c_{1}m\log^{c_{2}}n}}\right)^{p}s^{\mathcal{G}}

[TABLE]

So the $\ell_{p}^{p}$ term goes up by at most a factor $(10c_{1}\log^{c_{2}}n)^{p-1}m$. Similar calculation shows that the $\ell_{2}^{2}$ term goes up by at most a constant factor by the tree-portal routing. Thus, we get $T+{\overline{{E}}}^{r}_{i}\preceq_{\widetilde{O}(m^{1/(p-1)})}T+{\widehat{{E}}}^{r}_{i}$. The other direction is symmetric using the reverse tree-portal routing, and the calculation stays the same since the tree-portal routing in reverse incurs the same load/congestion on tree edges. ∎

If we put the ${\overline{{G}}}^{r}_{i}$ over all resistance $r$’s and round $i$’s together, we get

Lemma 4.15**.**

[TABLE]

The sum is over all possible resistance value $r$’s, and over all iterations $i$ for $r$.

Proof.

By Lemma 4.14 and Lemma 2.12 we have

[TABLE]

Note the ${\overline{{G}}}^{r}_{i}$’s (and the ${\widehat{{E}}}^{r}_{i}$’s) are disjoint for different resistance values or different iterations, thus these edges contribution to the objective function value simply adds up. For the tree edges, since there are at most $\log^{2}n$ different pairs of resistance and iteration pairs, we have

[TABLE]

by considering the mapping that split flow on one tree edge to $\log^{2}n$ copies of it and the reverse mapping of merging. Note \mathopen{}\mathclose{{}\left|a_{1}}\right|^{x}+\ldots+\mathopen{}\mathclose{{}\left|a_{1}}\right|^{x}\leq(\mathopen{}\mathclose{{}\left|a_{1}}\right|+\ldots+\mathopen{}\mathclose{{}\left|a_{k}}\right|)^{x}\leq k(\mathopen{}\mathclose{{}\left|a_{1}}\right|^{x}+\ldots+\mathopen{}\mathclose{{}\left|a_{1}}\right|^{x}). This gives the final result we want. ∎

Note that ${\overline{{G}}}$ is the disjoint union of $T+\sum_{r,i}{\overline{{G}}}^{r}_{i}+{\overline{{G}}}_{rest}$, while $\mathcal{H}$ is the disjoint union of $T+\sum_{r,i}{\overline{{H}}}^{r}_{i}+{\overline{{G}}}_{rest}$. Thus, we can show the following

Lemma 4.16**.**

${\overline{{G}}}\preceq_{\widetilde{O}(m^{2/(p-1)})}\mathcal{H}\preceq_{\widetilde{O}(m^{2/(p-1)})}{\overline{{G}}}$.

Proof.

Recall for each resistance value $r$, in the $i$-th round, ${\widehat{{G}}}^{r}_{i,j}$ is the $j$-th uniform expander we find, and ${\overline{{H}}}^{r}_{i,j}$ is the sparsified graph of ${\widehat{{G}}}^{r}_{i,j}$.

[TABLE]

${\overline{{G}}}\preceq_{\widetilde{O}(m^{2/(p-1)})}\mathcal{H}$ follows by taking the composition of all the intermediate steps, and multiplying the approximation error by Lemma 2.10. The other direction is similar. ∎

Now we can prove the main ultra-sparsification theorem.

Proof of Theorem 3.6.

Other than the additive error terms and the self-loops, everything in the theorem statement follow directly from Lemma 4.12 (the number of off-tree edges), and composition of Lemma 4.13 with Lemma 4.16 (the approximation error). We explicitly spell out the flow mappings between $\mathcal{G}$ and $\mathcal{H}$. We start with the $\mathcal{G}$ to $\mathcal{H}$ direction. We break the flow in $\mathcal{G}$ as the sum of flow on disjoint edge subsets $T$,${\overline{{G}}}_{rest}$, and ${\widetilde{{E}}}^{r}_{i}$, specify the mapping from each piece to $\mathcal{H}$, and later take the sum of the mappings. The mapping from $T$ and ${\overline{{G}}}_{rest}$ to $\mathcal{H}$ is just the identity. For flow on ${\widetilde{{E}}}^{r}_{i}$, we get a flow on $T+{\widehat{{E}}}^{r}_{i}$ by tree-portal routing. As ${\widehat{{E}}}^{r}_{i}$ is the sum of ${\widehat{{G}}}^{r}_{i,j}$’s, for the flow mapped to ${\widehat{{G}}}^{r}_{i,j}$, we map it to a flow on ${\overline{{H}}}^{r}_{i,j}$ using the flow mapping in SampleAndFixGradient. We add these mapping over all $j$’s to get a mapping from the flow on $T+{\widehat{{E}}}^{r}_{i}$ to a flow on $\mathcal{H}$, and take the composition with the tree-portal routing to get a mapping from ${\widetilde{{E}}}^{r}_{i}$ to $\mathcal{H}$. Summing over all $r,i$ (together with the identity on $T$ and ${\overline{{G}}}_{rest}$ gives the mapping from $\mathcal{G}$ to $\mathcal{H}$. The mapping from $\mathcal{H}$ to $\mathcal{G}$ is symmetric, and in the part from ${\widehat{{E}}}^{r}_{i}$ to $T+{\widetilde{{E}}}^{r}_{i}$ we use the reverse of tree-portal routing.

All the subroutines take nearly linear time, and we have at most $\log n$ different $r$, and for each $r$ there are at most $\log m$ iterations, so the overall running time is $\widetilde{O}(m)$. The flow mappings can also be applied in $\widetilde{O}(m)$ time, and they are linear maps.

Now we look at the additive error terms. In particular, additive errors come in at two places. The first is when we round an original resistance to [math] when it is less than $\delta$, and we have Lemma 4.11 to bound the error (at that step). The second place is in Decompose (Algorithm 6), we may get an expander ${\widehat{{G}}}^{r}_{i,j}$ whose projected gradient is not $\alpha$-uniform but has tiny norm (i.e. case $3(b)$), and we zero out its projection to the cycle space before sampling to make it $1$-uniform. If we have a flow $f$ on such an ${\widehat{{G}}}^{r}_{i,j}$, the additive error is in the linear term, and is equal to the dot product of $f$ with the removed gradient. We let $\boldsymbol{\mathit{g}}^{r}_{i}$, $\boldsymbol{\mathit{g}}^{r}_{i,j}$ be the gradient on edges in ${\widehat{{G}}}^{r}_{i}$,${\widehat{{G}}}^{r}_{i,j}$ respectively, and $\boldsymbol{\widehat{\mathit{g}}}^{r}_{i}$, $\boldsymbol{\widehat{\mathit{g}}}^{r}_{i,j}$ as the projection of $\boldsymbol{\mathit{g}}^{r}_{i}$ (and $\boldsymbol{\mathit{g}}^{r}_{i,j}$) to the cycle space of ${\widehat{{G}}}^{r}_{i}$ (and ${\widehat{{G}}}^{r}_{i,j}$). We remove $\boldsymbol{\widehat{\mathit{g}}}^{r}_{i,j}$ from the gradient $\boldsymbol{\mathit{g}}^{r}_{i,j}$ when $\boldsymbol{\widehat{\mathit{g}}}^{r}_{i,j}\leq\delta^{\prime}\boldsymbol{\widehat{\mathit{g}}}^{r}_{i}$ for some parameter $\delta^{\prime}$, so the additive error we introduce is $f^{T}\boldsymbol{\widehat{\mathit{g}}}^{r}_{i,j}$, which is at most \mathopen{}\mathclose{{}\left\|f}\right\|_{2}\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{r}_{i,j}}\right\|_{2}, which is in turn at most \delta^{\prime}\mathopen{}\mathclose{{}\left\|f}\right\|_{2}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{g}}^{r}_{i}}\right\|_{2} as $\boldsymbol{\widehat{\mathit{g}}}^{r}_{i}$ is a projection of $\boldsymbol{\mathit{g}}^{r}_{i}$. Now we look at how this additive error propagates in terms of the overall approximation error between $\mathcal{G}$ and $\mathcal{H}$. We will get an additional factor $m$ when we combine the additive errors over all the individual expanders where we carry out this perturbation. Note we are not really introducing more error here, but simply because \sqrt{m}\mathopen{}\mathclose{{}\left\|\sum_{i}f_{i}}\right\|\geq\sum_{i}\mathopen{}\mathclose{{}\left\|f_{i}}\right\|\geq\mathopen{}\mathclose{{}\left\|\sum_{i}f_{i}}\right\| when $f_{i}$’s have disjoint support and total size $m$. The additive error is also amplified through the intermediate steps, but since the multiplicative approximation errors are $m^{O}{1/p}$, we lose at most another polynomial factor. Additional polynomial factor comes in because the norm of the gradient vector after tree-routing can be off by a polynomial factor comparing to the norm of the original gradient. However, overall the blowup is at most polynomial, and we use a polynomially smaller $\delta^{\prime}$ in Decompose to accommodate these factors to get the additive error in our final result. The same argument applies to the additive error introduced by resistance rounding (e.g. round to [math] when the gradient is at most $\delta/m^{c}$ for some large enough $c$). ∎

We brief go over the case when tree-portal routing gives self-loops. We treat self-loops the same way as the edges that are in the uniform expanders except they don’t go through the expander decomposition and sampling steps. Once we get a self-loop ${\widehat{{e}}}$ from tree-portal routing of some edge $e\in E^{\mathcal{G}}$, we add ${\widehat{{e}}}$ to $\mathcal{H}$, where the gradient on ${\widehat{{e}}}$ is set (the same way as non self-loops) to preserve the flow dot gradient term under tree-portal routing. We remove its pre-image $e$ from $E^{r}_{i}$, but if in some iteration, more than half of the edges in $E^{r}_{i}$ are mapped to self-loops by tree-portal routing, we skip the decomposition and sampling steps also for other edges, as we don’t have a dense enough graph between the portal nodes to sparsify. We still have the size of $E^{r}_{i}$ drop by at least $1/2$ across each iteration as before. The final caveat is that since self-loops don’t go through SampleAndFixGradient, and thus their $s$ values are not scaled to be the same as the rest of the edges in $\mathcal{H}$. This is not an issue because we will remove them from the instance and optimize them individually (see Lemma 3.5), so they won’t exist in the instance that we recursively solve, so uniform $s$ scalar is not required for them.

5 Decomposing into Uniform Expanders

In this section we prove our decomposition result necessary for finding large portions of edges that can be sampled. This and the subsequent sampling step in Appendix D are critical for reducing the number of edges between portal vertices, after they were routed there in Line 12 of UltraSparsify (Algorithm 4). The main algorithmic guarantees can be summarized as below in Theorem 4.9.

See 4.9

We will obtain the expansion properties via expander decompositions. Specifically we will invoke the following result from [SW18] as a black box.

Lemma 5.1**.**

There is a routine ExpanderDecompose that when given any graph $G$ and any degrees $\boldsymbol{\mathit{d}}$ such that $\boldsymbol{\mathit{d}}_{u}\geq deg_{G}(u)$ for all $u$, along with a parameter $0<\phi<1$, $\textsc{ExpanderDecompose}(G,\boldsymbol{\mathit{d}},\phi)$ returns a partition of the vertices of $G$ into $V_{1},V_{2},\ldots$ in $O(m\phi^{-1}\log^{4}n)$ time such that $G[V_{i}]$ has conductance at least $\phi$ w.r.t. $\boldsymbol{\mathit{d}}_{u}$, and the number of edges between the $V_{i}$s is at most $O(\sum_{u}\boldsymbol{\mathit{d}}_{u}\phi\log^{3}n)$.

Note that we explicitly introduce the $\boldsymbol{\mathit{d}}$ vector containing the degrees of the initial graph because we will repeatedly invoke this partition routine. This is due to our other half of the routine, which is to repeatedly project $\boldsymbol{\mathit{g}}$ among the remaining edges, and removing the ones that contribute to too much of its $\ell_{2}^{2}$-norm in order to ensure uniformity as given in Case 2a of Theorem 4.9. To see that this process makes progress, we need the key observation from Lemma 2.15 that projections can only decrease the $\ell_{2}^{2}$ norm of $\boldsymbol{\widehat{\mathit{g}}}$, the projection of the gradient.

This leads to an approach where we alternate between dropping the edges with high energy, and repartitioning the remaining edges into expanders. Pseudocode of this routine is in Algorithm 6, which calls a recursive routine, DecomposeRecursive shown in Algorithm 7 with a suitable value of $\phi$ and number of layers. Note that we also need to trim the initial graph so that we only work with large degree vertices.

We will also need the following result (Lemma 28 of [KLOS14], see also [KM11]).

Lemma 5.2**.**

Suppose $G$ is a unit weight graph with conductance $\phi$. Then the projection operations into cycle and potential flow spaces both have $\ell_{\infty}$ norms bounded by $O(\phi^{-2}\log{n})$:

[TABLE]

and

[TABLE]

Proof.

(of Theorem 4.9)

We start by bounding the qualities of the $\mathcal{G}$ pieces returned. As we only return pieces that are the outputs of ExpanderDecompose, all of them have conductance at least $\phi$ by Lemma 5.1. Also, since we only keep the non-trivial pieces containing edges, taking the singleton cuts gives that the degrees in these pieces are at least

[TABLE]

Now consider the quality of each $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}_{i}}$: if it was returned due to $i=L$, then the energy of the projected gradient must have been halved at least $L-\log{n}$ times, or by a factor of $2^{L-\log{n}}=2^{\log(1/\delta)}=1/\delta$. Thus we would have

[TABLE]

Otherwise, we must have terminated because both the energy and edge count did not decrease too much. An edge $e$ was kept in the trimmed set only if

[TABLE]

Combining this with the termination requirement of \mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}_{i}}}\right\|_{2}^{2}\geq\frac{1}{2}\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}}\right\|_{2}^{2} gives that the $\ell_{\infty}$ norm of the pre-projection gradient on $\mathcal{G}_{i}$, $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}_{E^{\mathcal{G}_{i}}}$ satisfies

[TABLE]

On the other hand, because $\mathcal{G}_{i}$ has expansion $\phi$, doing an orthogonal cycle projection on it can only increase the $\ell_{\infty}$-norm of a vector by a factor of $O(\phi^{-2}\log{n})$ by Lemma 5.2. Thus we have

[TABLE]

which is the desired (post-projection) uniformity bound.

We now bound the number of edges removed during all the recursive calls. The bound on $L$ means this recursion has at most $O(\log(n/\delta))$ levels. Lemma 5.1 gives that the number of edges between the expander pieces is

[TABLE]

so the setting of $\phi=\frac{1}{c_{partition}\log^{3}n\log(n/\delta))}$ gives at most $m/10$ edges between the pieces for an appropriate choice of $c_{partition}$.

Furthermore, as each edge’s contribution to $\boldsymbol{\widehat{\mathit{g}}}$ is non-negative, the number of edges whose relative contribution exceed $\frac{10L}{m}$ is at most $\frac{m}{10L}$. Summing this over all levels gives at most $m/10$ edges removed from the trimming step on Line 3 of DecomposeRecursive in Algorithm 7.

Finally, the running time is dominated by the expander decomposition calls. As there are $O(\log(n/\delta))$ levels of recursion and each level deals with edge-disjoint subsets, we obtain the total running time by substituting the value of $\phi$ into the runtime of expander decompositions as given in Lemma 5.1. ∎

Acknowledgements

This project would not have been possible without Dan Spielman’s optimism about the existence of analogs of numerical methods for $\ell_{p}$-norms, which he has expressed to us on multiple occasions over the past six years. We also thank Ainesh Bakshi, Jelani Nelson, Aaron Schild, and Junxing Wang for comments and suggestions on earlier drafts and presentations of these ideas.

As with many recent works in optimization algorithms on graphs, this project has its large share of influence by the late Michael B. Cohen. In fact, Michael’s first papers on recursive preconditioning [CKM*+*14] and $\ell_{p}$-norm preserving sampling of matrices [CP15] directly influenced the constructions of preconditioners (Section 4.4) and uniform expanders (Section 5 and Appendix D) respectively. While our overall algorithm falls short of what Michael would consider ‘snazzy’, it’s also striking how many aspects of it he predicted, including: the use of expander decompositions; the $p\rightarrow\infty$ case being different than the $p\rightarrow 1$ case; and the large initial dependence on $p$ that’s also eventually fixable (see Section 1.4).

Richard regrets not being able to convince Michael to systematically investigate preconditioning for $\ell_{p}$-norm flows. He is deeply grateful to Aleksander Ma̧dry, Jon Kelner, Ian Munro, Tom Cohen, Marie Cohen, Sebastian Bubeck, and Ilya Razenshteyn for many helpful conversations following Michael’s passing.

Appendix A Deferred Proofs from Prelims, Section 2.2

See 2.9

Proof.

Consider the map $\mathcal{M}_{\mathcal{G}\rightarrow\mathcal{G}}$ such that for every flow $\boldsymbol{\mathit{f}}^{\mathcal{G}}$ on $\mathcal{G},$ we have $\mathcal{M}_{\mathcal{G}\rightarrow\mathcal{G}}(\boldsymbol{\mathit{f}}^{\mathcal{G}})=\boldsymbol{\mathit{f}}^{\mathcal{G}}.$ Thus,

[TABLE]

Moreover $(\boldsymbol{\mathit{B}}^{\mathcal{G}})^{\top}\mathcal{M}_{\mathcal{G}\rightarrow\mathcal{G}}(\boldsymbol{\mathit{f}}^{\mathcal{G}})=\boldsymbol{\mathit{B}}^{\mathcal{G}}\boldsymbol{\mathit{f}}^{\mathcal{G}}.$ Thus, the claims follow. ∎

See 2.10

Proof.

It is easy to observe that the given mapping is linear. Given a flow $\boldsymbol{\mathit{f}}^{\mathcal{G}_{1}}$ on $\mathcal{G}_{1},$ we have,

[TABLE]

Moreover,

[TABLE]

The same proof works for $\preceq^{\text{cycle}}.$ ∎

See 2.12

Proof.

Let $\boldsymbol{\mathit{f}}^{\mathcal{H}}$ be a flow on $\mathcal{H}.$ We write $\boldsymbol{\mathit{f}}^{\mathcal{H}}=(\boldsymbol{\mathit{f}}^{\mathcal{H}_{1}},\boldsymbol{\mathit{f}}^{\mathcal{H}_{2}}).$ Let $\boldsymbol{\mathit{f}}^{\mathcal{G}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathcal{M}_{\mathcal{H}\rightarrow\mathcal{G}}(\boldsymbol{\mathit{f}}^{\mathcal{H}}).$ If $\boldsymbol{\mathit{f}}^{\mathcal{G}_{i}}$ denotes $\mathcal{M}_{\mathcal{H}_{i}\rightarrow\mathcal{G}_{i}}(\boldsymbol{\mathit{f}}^{\mathcal{H}_{i}})$ for $i=1,2,$ then we know that $\boldsymbol{\mathit{f}}^{\mathcal{G}}=(\boldsymbol{\mathit{f}}^{\mathcal{G}_{1}},\boldsymbol{\mathit{f}}^{\mathcal{G}_{2}}).$ Thus, the objectives satisfy

[TABLE]

For the residues, we have,

[TABLE]

Thus, $\mathcal{H}\preceq_{\kappa}\mathcal{G}.$ ∎

See 2.13

Proof.

Consider the map $\mathcal{M}_{\mathcal{H}\rightarrow\mathcal{G}}(\boldsymbol{\mathit{f}})=\boldsymbol{\mathit{f}}.$ Thus, since the underlying graphs are the same, we immediately have $(\boldsymbol{\mathit{B}}^{\mathcal{G}})^{\top}\boldsymbol{\mathit{f}}=(\boldsymbol{\mathit{B}}^{\mathcal{H}})^{\top}\boldsymbol{\mathit{f}}.$ For the objective, we have

[TABLE]

∎

Appendix B Deferred Proofs for Numerical Methods from

Section 3

The following simple lemma characterizes the change in smoothed $\ell_{p}$-norms under rescaling of the input vector.

Lemma B.1**.**

For all $\boldsymbol{\mathit{x}}\in\mathbb{R}^{m},\boldsymbol{\mathit{r}}\in\mathbb{R}^{m}_{\geq 0},s\in\mathbb{R}_{\geq 0},$ and $\lambda\in\mathbb{R}$, we have,

[TABLE]

Proof.

It suffices to prove the claim for $x\in\mathbb{R},r\in\mathbb{R}_{\geq 0},s\in\mathbb{R}_{\geq 0}$

[TABLE]

Since all terms are non-negative, we get,

[TABLE]

∎

See 3.1

Proof.

Note that all the terms are a sum over the coordinates. Thus, it suffices to prove the inequality for $x,\delta\in\mathbb{R},$ and $r,s\in\mathbb{R}_{\geq 0}.$ We have,

[TABLE]

where $\delta^{\prime}=\delta/x.$

Lemma B.2, proved later, proves that for all $\delta^{\prime},$ and $p\geq 2,$ we have,

[TABLE]

Thus,

[TABLE]

Lemma B.3, proved later, shows that for all $\delta^{\prime},$ and $p\geq 2,$ we have,

[TABLE]

[TABLE]

∎

Lemma B.2**.**

For all $\delta\in\mathbb{R},$ $p\geq 1,$ we have,

[TABLE]

Proof.

The proof has to consider several cases.

$\boldsymbol{\delta\geq 1}$.

Using mean-value theorem, we know there is some $z\in[0,\delta]$ such that

[TABLE]

$\boldsymbol{0\leq\delta\leq 1}$.

Using mean-value theorem, we know there is some $z\in[0,\delta]$ such that

[TABLE]

If $p\leq 2,$ we have $(1+\delta)^{p-1}$ is a concave function, and hence $(1+\delta)^{p-1}\leq 1+(p-1)\delta.$ If $p\geq 2,$ we have $(1+\delta)^{p-1}$ is a convex function, and hence for $\delta\in[0,1],$ we have $(1+\delta)^{p-1}\leq 1+(2^{p-1}-1)\delta.$ Thus,

[TABLE]

$\boldsymbol{-1\leq\delta\leq 0}$.

Using mean-value theorem, we know there is some z\in[-\mathopen{}\mathclose{{}\left|\delta}\right|,0] such that

[TABLE]

If $p\leq 2,$ we have $(1+\delta)^{p-1}$ is a concave function, and hence for $\delta\in[-1,1],$ we have $(1+\delta)^{p-1}\geq 1+\delta.$ If $p\geq 2,$ we have $(1+\delta)^{p-1}$ is a convex function, and hence $(1+\delta)^{p-1}\geq 1+(p-1)\delta.$ Thus,

[TABLE]

$\boldsymbol{\delta\leq-1}$.

We have,

[TABLE]

since \mathopen{}\mathclose{{}\left|\delta}\right|\geq 1. ∎

Lemma B.3**.**

For all $\delta\in\mathbb{R},p\geq 2,$ we have

[TABLE]

Proof.

Let $h(\delta)$ denote the function

[TABLE]

Thus, $h(0)=0.$ As for the previous proof, we consider several cases:

$\boldsymbol{\delta\geq 0}$.

We have $h(\delta)=(1+\delta)^{p}-1-p\delta-2^{-p}(\delta^{2}+\delta^{p}).$ Thus,

[TABLE]

Observe that since $p\geq 2,$ we have $(1+\delta)^{p-2}\geq\max\{1,\delta^{p-2}\}\geq 2^{-1}(1+\delta^{p-2}),$ and $p(p-1)\geq 2.$ Thus,

[TABLE]

Since $h(0)=h^{\prime}(0)=0,$ and $h^{\prime\prime}(\delta)\geq 0,$ for all $\delta\geq 0,$ we must have $h(\delta)\geq 0$ for all $\delta\geq 0.$

$\boldsymbol{-1\leq\delta\leq 0}$.

We have,

[TABLE]

Since $0\leq 1+\delta\leq 1,$ and $p-1\geq 1,$ we have $(1+\delta)^{p-1}\leq 1+\delta,$ and \mathopen{}\mathclose{{}\left|\delta}\right|^{p-1}\leq\mathopen{}\mathclose{{}\left|\delta}\right|^{1}=-\delta. Thus,

[TABLE]

$\boldsymbol{\delta\leq-1}$.

We have,

[TABLE]

Since \mathopen{}\mathclose{{}\left|\delta}\right|\geq 1, and $p-1\geq 1,$ we have -\delta=\mathopen{}\mathclose{{}\left|\delta}\right|\leq\mathopen{}\mathclose{{}\left|\delta}\right|^{p-1}. Thus,

[TABLE]

Now, observe that since $0\leq-1-\delta,$ and $p-1\geq 1,$

[TABLE]

Thus,

[TABLE]

and hence $h^{\prime}(\delta)\leq 0$ for $\delta\leq-1.$

For the last two cases, since $h(0)=0,$ and $h^{\prime}(\delta)\leq 0,$ for all $\delta\leq 0.$ Thus, we must have $h(\delta)\geq 0,$ for $\delta\leq 0.$ ∎

See 3.2

Proof.

Let $\boldsymbol{\mathit{x}}^{\star}$ to be an optimal solution to Problem 10. Consider $\widetilde{\boldsymbol{\delta}}=\boldsymbol{\mathit{x}}^{\star}-\boldsymbol{\mathit{x}}_{0}.$ Thus,

[TABLE]

Thus, $\widetilde{\boldsymbol{\delta}}$ is a feasible solution to Problem 13. Moreover, it satisfies,

[TABLE]

Now, given a feasible solution $\delta$ to Problem 13, we must have $\boldsymbol{\mathit{A}}\boldsymbol{\delta}=\boldsymbol{\mathbf{0}}.$ Thus, $\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}_{1}=\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}_{0}+2^{-3p}\boldsymbol{\mathit{A}}\boldsymbol{\delta}=b,$ and $\boldsymbol{\mathit{x}}_{1}$ is a feasible solution to Problem 10. Moreover,

[TABLE]

∎

Appendix C Elimination of Low-Degree Vertices, and Loops

In this section, we that the instance $\mathcal{H}$ returned by UltraSparsify can be reduced to a smaller graph by repeatedly eliminating vertices of degree at most 2. This step is analogous to the partial Cholesky factorization in the Laplacian solver of Spielman and Teng [ST14]. A slight technical issue is that if we run into a cycle where at most $1$ vertex on the cycle has edge(s) to the rest of the graph, the elimination of the degree $2$ nodes on the cycle essentially becomes an optimization problem on only the cycle edges that can be solved independently from the rest of the graph.

See 3.4

Proof.

We first observe that a self-loop $e\in E^{\mathcal{G}}$ on a vertex $v\in V^{\mathcal{G}}$ does not contribute to the residue at any vertex, including $v.$ Thus, the circulation constraint on a flow $\boldsymbol{\mathit{f}}^{\mathcal{G}}$ does not impose any constraint on $\boldsymbol{\mathit{f}}^{\mathcal{G}}_{e}.$ Moreover, since the objective $\alpha^{\mathcal{G}}$ can be written as a sum over the edges, for every self-loop $e,$ the variable $\boldsymbol{\mathit{f}}^{\mathcal{G}}_{e}$ is independent of all other variables. Thus, we can ignore the self-loops in remainder of the proof.

We first prove that we can repeatedly eliminate vertices of non-selfloop degree 1 in $G$ while preserving $\mathcal{E}$ exactly for a circulation. Consider one such vertex $v\in V^{\mathcal{G}},$ and let $e=({v,u})\in E^{\mathcal{G}}$ be the only non-selfloop edge incident on $v$ (the argument for the reverse direction is identical). Given any circulation $\boldsymbol{\mathit{f}}^{\mathcal{G}},$ since the only non-selfloop edge incident on $v$ is $e,$ we must have $\boldsymbol{\mathit{f}}^{\mathcal{G}}_{e}=0.$ Thus, we can drop $e$ entirely from the instance. Formally, we define

[TABLE]

We let the mapping $\mathcal{M}_{\mathcal{G}\rightarrow\mathcal{G}^{\prime}}$ to be just the projection on to $E^{\mathcal{G}^{\prime}}.$ Thus, $\mathcal{M}_{\mathcal{G}\rightarrow\mathcal{G}^{\prime}}(\boldsymbol{\mathit{f}}^{\mathcal{G}})=\boldsymbol{\mathit{f}}^{\mathcal{G}}|_{E^{\mathcal{G}^{\prime}}}.$ Since $\boldsymbol{\mathit{f}}^{\mathcal{G}}_{e}=0,$ we immediately get $\mathcal{E}^{\mathcal{G}^{\prime}}(\mathcal{M}_{\mathcal{G}\rightarrow\mathcal{G}^{\prime}}(\boldsymbol{\mathit{f}}^{\mathcal{G}}))=\mathcal{E}^{\mathcal{G}}(\boldsymbol{\mathit{f}}^{\mathcal{G}}).$ Thus, $\mathcal{G}\preceq^{\text{cycle}}_{1}\mathcal{G}^{\prime}.$

Now, consider the mapping $\mathcal{M}_{\mathcal{G}^{\prime}\rightarrow\mathcal{G}}$ that pads a circulation $\boldsymbol{\mathit{f}}^{\mathcal{G}^{\prime}}$ on $\mathcal{G}^{\prime}$ with 0 on $e,$ i.e.,

[TABLE]

Again, it is immediate that $\mathcal{E}^{\mathcal{G}}(\mathcal{M}_{\mathcal{G}^{\prime}\rightarrow\mathcal{G}}(\boldsymbol{\mathit{f}}^{\mathcal{G}^{\prime}}))=\mathcal{E}^{\mathcal{G}^{\prime}}(\boldsymbol{\mathit{f}}^{\mathcal{G}^{\prime}}).$ Thus, $\mathcal{G}^{\prime}\preceq^{\text{cycle}}_{1}\mathcal{G}.$

We can repeatedly apply the above transformation to eliminate all vertices of non-selfloop degree 1 in $G.$ For convenience, we let $\mathcal{G}^{\prime}$ denote the final instance obtained. Thus, we have, $\mathcal{G}^{\prime}\preceq^{\text{cycle}}_{1}\mathcal{G}\preceq^{\text{cycle}}_{1}\mathcal{G}^{\prime}.$

Now, we will replace maximal paths with all internal vertices of non-selfloop degree 2 with single edges. Consider such a path $P.$ Formally, $P$ is a path of length $l$ in $G,$ say $P=(v_{0},v_{1},\ldots,v_{l-1},v_{l}),$ with all of $v_{1},\ldots,v_{l-1}$ having degree exactly 2, and $v_{0},v_{l}$ have non-selfloop degree at least 3. For convenience, we assume that all edges $(v_{i-1},v_{i})$ are oriented in the same direction. Observe that for a circulation $\boldsymbol{\mathit{f}}^{\mathcal{G}^{\prime}},$ the flow on all the edges $(v_{i-1},v_{i})$ must be the same, i.e., $\boldsymbol{\mathit{f}}^{\mathcal{G}^{\prime}}_{(v_{i-1},v_{i})}$ must all be equal. Thus, we can replace $P$ with a single edge $e_{P}$ while preserving the amount of flow and the direction.

Formally, let $\mathcal{P}=\{P_{1},\ldots,P_{t}\}$ denote the set of all maximal paths in $G^{\prime}=(V^{\mathcal{G}^{\prime}},E^{\mathcal{G}^{\prime}})$ such that all their internal vertices have non-selfloop degree exactly 2 in $G^{\prime}.$ We replace each of these paths with a new edge connecting its endpoints. Let,

[TABLE]

We define the mapping $\mathcal{M}_{\mathcal{G}^{\prime}\rightarrow\mathcal{G}^{\prime\prime}}$ as follows

[TABLE]

We define $\mathcal{M}_{\mathcal{G}^{\prime\prime}\rightarrow\mathcal{G}^{\prime}}$ to be the inverse map of $\mathcal{M}_{\mathcal{G}^{\prime}\rightarrow\mathcal{G}^{\prime\prime}}.$

[TABLE]

It follows from the definitions that for every circulation $\boldsymbol{\mathit{f}}^{\mathcal{G}^{\prime}},$ letting $\boldsymbol{\mathit{f}}^{\mathcal{G}^{\prime\prime}}$ denote $\mathcal{M}_{\mathcal{G}^{\prime}\rightarrow\mathcal{G}^{\prime\prime}}(\boldsymbol{\mathit{f}}^{\mathcal{G}^{\prime}}),$ we have, $\mathcal{M}_{\mathcal{G}^{\prime\prime}\rightarrow\mathcal{G}^{\prime}}(\boldsymbol{\mathit{f}}^{\mathcal{G}^{\prime\prime}})=\boldsymbol{\mathit{f}}^{\mathcal{G}^{\prime}}.$ Moreover,

[TABLE]

where the last inequality follows since for every path of length $l$, the contribution to the $\ell_{p}^{p}$ changes by a factor of $l^{-1},$ and since the paths must be vertex-disjoint, $l\leq|V^{\mathcal{G}}|\leq n.$ The above inequalities imply,

[TABLE]

and hence $\mathcal{G}^{\prime}\preceq^{\text{cycle}}_{1}\mathcal{G}^{\prime\prime}.$ Combined with $\mathcal{G}\preceq^{\text{cycle}}_{1}\mathcal{G}^{\prime},$ we get $\mathcal{G}\preceq^{\text{cycle}}_{1}\mathcal{G}^{\prime\prime}.$ Moreover, we have, for $\kappa=n^{\frac{1}{p-1}},$

[TABLE]

Thus, $\mathcal{G}^{\prime\prime}\preceq^{\text{cycle}}_{\kappa}\mathcal{G}^{\prime}.$ Combining with $\mathcal{G}^{\prime}\preceq^{\text{cycle}}_{1}\mathcal{G},$ we obtain $\mathcal{G}^{\prime\prime}\preceq^{\text{cycle}}_{\kappa}\mathcal{G}.$ The final instance returned is $\mathcal{G}^{\prime\prime},$ giving us our theorem. ∎

See 3.5

Proof.

Let $E^{\prime}$ denote the set of all self-loops in $E^{\mathcal{G}}.$ Then, we define $\mathcal{G}_{1}$ to be the instance obtained by removing all edges in $E^{\prime}.$ Formally,

[TABLE]

We define $\mathcal{G}_{2}$ be the instance $\mathcal{G}$ restricted to $E^{\prime}.$ Thus,

[TABLE]

It is immediate that $\mathcal{G}=\mathcal{G}_{1}\cup\mathcal{G}_{2}.$ Since $\mathcal{G}_{2}$ only has self-loops, we have that for every $\boldsymbol{\mathit{f}}^{\mathcal{G}_{2}},$ we have $(\boldsymbol{\mathit{B}}^{\mathcal{G}_{2}})^{\top}\boldsymbol{\mathit{f}}^{\mathcal{G}_{2}}=\boldsymbol{\mathbf{0}}.$ Thus, the constraint $(\boldsymbol{\mathit{B}}^{\mathcal{G}_{2}})^{\top}\boldsymbol{\mathit{f}}^{\mathcal{G}_{2}}=\boldsymbol{\mathbf{0}}$ is vacuous.

Now, observe that in the absence of linear constraints on $\boldsymbol{\mathit{f}}^{\mathcal{G}_{2}},$ the variables $\boldsymbol{\mathit{f}}^{\mathcal{G}_{2}}_{e}$ are independent for all $e\in E^{\mathcal{G}_{2}}.$ Moreover, we have

[TABLE]

Thus, we can solve for each $\boldsymbol{\mathit{f}}^{\mathcal{G}_{2}}_{e}$ independently. Now, consider a fixed $e\in E^{\mathcal{G}_{2}}.$ We write $f_{e}$ for $\boldsymbol{\mathit{f}}^{\mathcal{G}_{2}}_{e}.$ We wish to solve

[TABLE]

Note that the objective function is concave. The gradient of $\mathcal{E}^{\mathcal{G}_{2}}_{e}(f_{e})$ with respect to $f_{e}$ is

[TABLE]

First observe that if $f^{\star}_{e}$ is the optimal solution, it must have the same sign as $\boldsymbol{\mathit{g}}_{e}^{\mathcal{G}_{2}}.$ Without loss of generality, we assume that $\boldsymbol{\mathit{g}}^{\mathcal{G}_{2}}_{e}\geq 0.$ Observe that for $f_{e}\geq\frac{\boldsymbol{\mathit{g}}^{\mathcal{G}_{2}}_{e}}{2\boldsymbol{\mathit{r}}_{e}^{\mathcal{G}_{2}}},$ we have $\frac{\mathsf{d}}{{\mathsf{d}}f_{e}}\mathcal{E}^{\mathcal{G}_{2}}_{e}\leq 0.$ Thus, $f^{\star}_{e}\leq\frac{\boldsymbol{\mathit{g}}^{\mathcal{G}_{2}}_{e}}{2\boldsymbol{\mathit{r}}_{e}^{\mathcal{G}_{2}}}.$ Similarly, we have, f^{\star}_{e}\leq\mathopen{}\mathclose{{}\left(\frac{\boldsymbol{\mathit{g}}^{\mathcal{G}_{2}}_{e}}{ps^{\mathcal{G}_{2}}}}\right)^{\frac{1}{p-1}}, where $f^{\star}_{e}$ is the optimal solution. Thus, if we define $z$ as

[TABLE]

then $f^{\star}_{e}\leq z.$

Moreover, for $f_{e}\leq\frac{z}{2},$ we have,

[TABLE]

Thus, $f^{\star}_{e}\geq\frac{z}{2},$ and hence $z$ gives a 2-approximation to $f^{\star}$ that can be computed in $O(1)$ time. Now, applying binary search allows us to find $f_{e}\in[(1-\nicefrac{{\delta}}{{p}})f_{e}^{\star},(1+\nicefrac{{\delta}}{{p}})f^{\star}_{e}]$ in $O(\log\nicefrac{{1}}{{\delta}})$ time. Now, we show that such an estimate is good enough. Consider the point $\frac{3}{4}z.$ We have

[TABLE]

Now,

[TABLE]

Rewriting, we get $\mathcal{E}^{\mathcal{G}_{2}}_{e}(f_{e})\geq(1-16\delta)\max_{f_{e}}\mathcal{E}^{\mathcal{G}_{2}}_{e}(f_{e}).$ Rescaling $\delta,$ we obtain our claim.

We can compute such an estimate for all the edges in O(\mathopen{}\mathclose{{}\left|E^{\mathcal{G}_{2}}}\right|\log\nicefrac{{1}}{{\delta}}) time. ∎

Appendix D Sparsifying Uniform Expanders

We now verify that sparsifying that sampling $\alpha$-uniform expanders preserve the objectives of the optimizations. Pseudocode of our routine and the flow maps constructed by it are in Algorithm 9. We remark that the maps are identical to the ones used for flow sparsifiers by Kelner et al. [KLOS14].

Note that the decomposition in Line (5) can be found by first computing \boldsymbol{\mathit{\psi}}=\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{B}}^{\mathcal{G}\top}\boldsymbol{\mathit{B}}^{\mathcal{G}}}\right)^{{\dagger}}\boldsymbol{\mathit{B}}^{\mathcal{G}\top}\boldsymbol{\mathit{g}}^{\mathcal{G}}. The only randomness in the Algorithm 9 is in the sampling in Line (3). In Lines (5), and (7)-(10), when the pseudo-inverse of a Laplacian is applied, we can rely deterministically on a high-accuracy approximation based on the fact that if the earlier sampling succeeded, both matrices are Laplacians of expanders, and hence well-conditioned. Alternatively, we can call a high-accuracy Laplacian solver. This encurs another small failure probability. In either case, we can ensure that an implicit representation of the operator is only computed once, and succeeds with high probability.

See 4.10

Remark D.1**.**

If in Algorithm 9, the input gradient $\boldsymbol{\mathit{g}}^{\mathcal{G}}$ has zero cycle-space projection, i.e. $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}=\boldsymbol{\mathbf{0}}$, then the cycle-space gradient terms in Lines (9) and (10) should be set to zero, so that $\mathcal{M}_{\mathcal{G}\rightarrow\mathcal{H}}$ is the map \boldsymbol{\mathit{f}}\to\boldsymbol{\mathit{B}}^{\mathcal{H}}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{B}}^{H\top}\boldsymbol{\mathit{B}}^{\mathcal{H}}}\right)^{{\dagger}}\boldsymbol{\mathit{B}}^{\mathcal{G}\top}\boldsymbol{\mathit{f}} and $\mathcal{M}_{\mathcal{H}\rightarrow\mathcal{G}}$ is the map \boldsymbol{\mathit{f}}\to\boldsymbol{\mathit{B}}^{\mathcal{G}}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{B}}^{\mathcal{G}\top}\boldsymbol{\mathit{B}}^{\mathcal{G}}}\right)^{{\dagger}}\boldsymbol{\mathit{B}}^{H\top}\boldsymbol{\mathit{f}}. The proof of this case is simpler, we omit all terms that deal with cycle-space projected gradients and everything else stays the same as in the proof given in this section.

To prove this theorem, we first collect a number of observations that will help us. The most basic of these is that Line (3) succeeds in producing a sparsifier in the spectral approximation sense, and with edge set $F$ satisfying $0.5\tau m\leq|F|\leq 2\tau m$. This is a direct consequence of matrix concentration bounds [Tro12].

Lemma D.2**.**

Consider the edge-vertex incidence matrices with gradients (projected via $\boldsymbol{\mathit{\psi}}^{\mathcal{G}}$) appended as an extra column for both $\mathcal{G}$ and $\mathcal{H}$, $[\boldsymbol{\mathit{B}}^{\mathcal{G}},\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}]$ and $[\boldsymbol{\mathit{B}}^{\mathcal{H}},\boldsymbol{\widetilde{\mathit{g}}}^{\mathcal{H}}]$. With high probability we have that for any vector $\boldsymbol{\mathit{x}}$

[TABLE]

and the edge set $F$ of $H$ satisfies $0.5\tau m\leq|F|\leq 2\tau m$.

Proof.

The bounds on \mathopen{}\mathclose{{}\left|F}\right| follow from a scalar Chernoff bound. For the matrix approximation bound, we will invoke matrix Chernoff bounds [Tro12], which give such a bound as long as the rows of $[\boldsymbol{\mathit{B}}^{\mathcal{G}},\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}]$ are sampled with probaiblity exceeding $c_{sample}\log{n}$ times their leverage scores.

So it suffices to bound the leverage scores of the rows of this matrix. As $\boldsymbol{\mathit{B}}^{\mathcal{G}}$ and $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}$ are orthogonal to each other, we can bound the leverage scores of the rows in these two matrices and add them.

The fact that the graph $(V^{\mathcal{G}},E^{\mathcal{G}})$ has expansion $\phi$ means that its normalized Laplacian has eigenvalue at least $\phi^{-2}$. So the leverage score of a row of $\boldsymbol{\mathit{B}}^{\mathcal{G}}$ is at least $\phi^{-2}d_{\min}$. The leverage score of $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}_{e}$ in $\boldsymbol{\widehat{\mathit{g}}}$ on the other hand is at most $\alpha/m$ due to the $\alpha$-uniform assumption. Thus, the sampling probablity $\tau$ meets the requirements of matrix Chernoff bounds, and we get the approximation with high probability.

∎

Corollary D.3**.**

Assuming Equation (18), the graphs underlying $\mathcal{G}$ and $\mathcal{H}$ (with resistances $r^{\mathcal{G}}$ and $r^{\mathcal{H}}$) are spectral approximations of each other:

[TABLE]

and the subset of gradient terms chosen after rescaling, $\boldsymbol{\widetilde{\mathit{g}}}^{\mathcal{H}}$, has $\ell_{2}^{2}$ norm that’s bigger by a factor of about $\tau$:

[TABLE]

Proof.

The approximation of graphs follows from considering vectors $\boldsymbol{\mathit{x}}$ with [math] in the last coordinate.

The approximation of $\ell_{2}^{2}$ norms of vectors follow from considering the indicator vector with $1$ in the last column and [math] everywhere else. ∎

From this spectral approximation, we can also conclude that $(V^{\mathcal{H}},E^{\mathcal{H}})$ must be an expander, as captured by the next corollary.

Corollary D.4**.**

Assuming Equation (18), $H$ has conductance at least $0.8\phi$.

Proof.

Let $C_{\mathcal{G}}(S)$ and $C_{\mathcal{H}}(S)$ denote the number of edges of $E^{\mathcal{G}}$ and $E^{\mathcal{H}}$ respectively crossing a cut $S\subseteq V^{\mathcal{G}}=V^{\mathcal{H}}$. Condition (19) implies that cuts are preserved between $(V^{\mathcal{G}},E^{\mathcal{G}})$ and $(V^{\mathcal{H}},E^{\mathcal{H}})$: For all $S\subseteq V$ $\tau C_{\mathcal{G}}(S)\approx_{0.1}C_{\mathcal{H}}(S)$, by computing the quadratic form in an indicator vector of $S$.

The degree of every vertex is also preserved, to up a scaling of $\tau$ and a multiplicative error $1\pm 0.1$, i.e. $\tau\deg_{\mathcal{H}}(v)\approx_{0.1}\deg_{\mathcal{G}}(v)$. This follows from considering the quadratic form of Condition (19) in the indicator vector of vertex $v$. This implies for any $S$,

[TABLE]

from which we conclude the conductance is preserved up to a factor of $0.8$. ∎

We can also conclude from this that $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}$ is well-spread.

Corollary D.5**.**

Assuming Equation (18), the projection of $\boldsymbol{\mathit{g}}^{\mathcal{H}}$ onto the cycle-space of $H$, $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}$ has $\ell_{1}$ and $\ell_{2}$ norms that are close to $\tau$ times the corresponding terms in $\mathcal{G}$:

[TABLE]

[TABLE]

And $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}$ is $O(\alpha\phi^{-6}\log^{2}{n})$-well spread, i.e. (as $|F|$ is the number of entries of $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}$)

[TABLE]

Proof.

We first show \mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}}\right\|_{2}^{2}\approx_{0.5}\tau\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}}\right\|_{2}^{2}. Note that \mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}}\right\|_{2}^{2}=\mathopen{}\mathclose{{}\left\|{\boldsymbol{\mathit{B}}^{\mathcal{G}}}\boldsymbol{\mathbf{0}}+\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}}\right\|_{2}^{2} Now, consider $\boldsymbol{\mathit{x}}=\begin{pmatrix}\boldsymbol{\mathbf{0}}\\ 1\end{pmatrix}$. By Equation (18),

[TABLE]

Thus \tau\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}}\right\|_{2}^{2}\geq 0.9\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}}\right\|_{2}^{2}.

The definition of $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}$ ensures \mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}}\right\|_{2}^{2}=\min_{\boldsymbol{\mathit{y}}}\mathopen{}\mathclose{{}\left\|{\boldsymbol{\mathit{B}}^{\mathcal{G}}}\boldsymbol{\mathit{y}}+\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}}\right\|_{2}^{2}. Letting \boldsymbol{\mathit{y}}^{\mathcal{H}}\in\arg\min_{\boldsymbol{\mathit{y}}}\mathopen{}\mathclose{{}\left\|{\boldsymbol{\mathit{B}}^{\mathcal{H}}}\boldsymbol{\mathit{y}}+\boldsymbol{\widetilde{\mathit{g}}}^{\mathcal{H}}}\right\|_{2}^{2}, we get from the definition of $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}$, that

[TABLE]

Thus we also have \tau\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}}\right\|_{2}^{2}\geq 0.9\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}}\right\|_{2}^{2}, allowing us to conclude that Equation (21) is satisfied.

Next, observe that by Lemma 5.2, since $H$ has conductance at least $0.8\phi$,

[TABLE]

where in the last step we also used \mathopen{}\mathclose{{}\left\|\boldsymbol{\widetilde{\mathit{g}}}^{\mathcal{H}}}\right\|_{\infty}\leq\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}}\right\|_{\infty}, since the former vector consists of a subset of the entries of the latter. Furthermore, by combining the above inequality with the assumption that $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}$ is $\alpha$-well-spread, and Equation (21) holds, we get

[TABLE]

As $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}$ has \mathopen{}\mathclose{{}\left|F}\right|\leq 2\tau m entries, this shows that it is $O(\alpha\phi^{-6}\log^{2}n)$-well-spread, which establishes Equation (23).

Next, to prove that Equation (22) holds, we first observe that for any $\gamma$-well-spread vector on $\boldsymbol{\mathit{x}}$ with $t$ coordinates, since \mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{x}}}\right\|_{1}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{x}}}\right\|_{\infty}\geq\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{x}}}\right\|_{2}^{2},

[TABLE]

So \frac{t^{1/2}}{\gamma}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{x}}}\right\|_{2}\leq\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{x}}}\right\|_{1}\leq t^{1/2}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{x}}}\right\|_{2}. As $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}$ and $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}$ are $\alpha$ and $O(\alpha\phi^{-6}\log^{2}n)$-well-spread respectively, we then get

[TABLE]

Combining this with \mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}}\right\|_{2}^{2}\approx_{0.2}\tau\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}}\right\|_{2}^{2}, and $0.5\tau\leq\frac{|F|}{m}\leq 2\tau$, we get

[TABLE]

From this we can directly conclude that Equation (22) holds. ∎

Also, we can show that for any $\boldsymbol{\mathit{b}}$ and $\theta$, the optimum $\ell_{2}$ as well as $\ell_{p}$ energies are close to the per degree lower bounds. We start with the lower bounds.

Lemma D.6**.**

Consider any graph with degrees $\boldsymbol{\mathit{D}}$, uniform $r$ and $s$, gradient $\boldsymbol{\mathit{g}}$ decomposable into $\boldsymbol{\mathit{g}}=\boldsymbol{\widehat{\mathit{g}}}+\boldsymbol{\mathit{B}}\boldsymbol{\mathit{\psi}}$ where $\boldsymbol{\widehat{\mathit{g}}}$ is the cycle-space projection of $\boldsymbol{\mathit{g}}$, and any $\theta$, any flow $\boldsymbol{\mathit{f}}$ such that:

$\boldsymbol{\mathit{f}}$* has residues $\boldsymbol{\mathit{b}}$: $\boldsymbol{\mathit{B}}^{\top}\boldsymbol{\mathit{f}}=\boldsymbol{\mathit{b}}$, and* 2. 2.

$\boldsymbol{\mathit{f}}$* has dot product $\theta+\boldsymbol{\mathit{b}}^{\top}\boldsymbol{\mathit{\psi}}$ with $\boldsymbol{\mathit{g}}$, i.e. $\boldsymbol{\mathit{g}}^{\top}\boldsymbol{\mathit{f}}=\theta+\boldsymbol{\mathit{b}}^{\top}\boldsymbol{\mathit{\psi}}$*

must satisfy

[TABLE]

Proof.

First, note that because $\boldsymbol{\mathit{B}}^{\top}\boldsymbol{\mathit{f}}=\boldsymbol{\mathit{b}}$, we have

[TABLE]

That is, the dot of $\boldsymbol{\mathit{f}}$ against $\boldsymbol{\widehat{\mathit{g}}}$ must be $\theta$.

The total energy is a sum of the $\ell_{2}^{2}$ and $\ell_{p}^{p}$ terms. First, we will give two different lower bounds on the $\ell_{2}^{2}$, and can hence also conclude that the average of the two lower bounds is another lower bound. We then do the same for the $\ell_{p}^{p}$ terms and add the lower bounds together for a lower bound on the overall objective. We will do so separately, by matching the $\ell_{2}^{2}$ terms to the electrical energy, and the $\ell_{p}^{p}$ terms to the minimum congestion.

\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|_{2}^{2}\geq\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{b}}}\right\|_{\boldsymbol{\mathit{D}}^{-1}}^{2}: here we use the fact that the minimum energy of the electrical flow is given by

[TABLE]

and that the graph Laplacian is dominated by twice its diagonal

[TABLE]

to get

[TABLE] 2. 2.

\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|_{2}^{2}\geq\frac{\theta^{2}}{\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}}\right\|_{2}^{2}} is by rearranging Cauchy-Schwarz inequality, which in its simplest form gives

[TABLE] 3. 3.

\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|_{p}^{p}\geq\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{D}}^{-1}\boldsymbol{\mathit{b}}}\right\|_{\infty}^{p} is because if we have a residue of $\boldsymbol{\mathit{b}}_{u}$ at some vertex, then some edge incident to $u$ must have flow at least

[TABLE]

on it. The $p$-th power of that lower bounds the overall $p$-norm energy. 4. 4.

\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|_{p}^{p}\geq(\frac{\mathopen{}\mathclose{{}\left|\theta}\right|}{\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}}\right\|_{1}})^{p} uses a similar lower bound on \mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|_{\infty}, except using Holder’s inequality on $\ell_{\infty}$ and $\ell_{1}$ norms to obtain

[TABLE]

which rearranges to give

[TABLE]

∎

Before proving upper bounds on the energy required to route a flow in the graph, we state a lemma that upper bounds the energy required to route the “electrical” component of the flow, i.e. the projection of the flow orthogonal to the cycle space.

Lemma D.7**.**

Consider a graph $G$ with degrees $\boldsymbol{\mathit{D}}$, conductance $\phi$, and edge-vertex incidence matrix $\boldsymbol{\mathit{B}}$, and any demand $\boldsymbol{\mathit{b}}\bot\boldsymbol{\mathbf{1}}$. Define the electrical flow \boldsymbol{\mathit{f}}=\boldsymbol{\mathit{B}}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{B}}^{\top}\boldsymbol{\mathit{B}}}\right)^{{\dagger}}\boldsymbol{\mathit{b}}. Then \sum_{e}\boldsymbol{\mathit{f}}_{e}^{2}\leq 2\phi^{-2}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{b}}}\right\|_{\boldsymbol{\mathit{D}}^{-1}}^{2}.

The proof relies on first Cheeger’s Inequality (e.g. see [Spi18]):

Theorem D.8**.**

*(Cheeger’s Inequality)

Consider a graph $G$ with degrees $\boldsymbol{\mathit{D}}$, conductance $\phi$, and adjancency matrix $\boldsymbol{\mathit{A}}$. Then*

[TABLE]

We will also need the following helpful fact.

Fact D.9**.**

Suppose $\boldsymbol{\mathit{M}}=\boldsymbol{\mathit{X}}\boldsymbol{\mathit{A}}\boldsymbol{\mathit{X}}^{\top}$ where $\boldsymbol{\mathit{A}}$ is symmetric and $\boldsymbol{\mathit{X}}$ is non-singular, and that $\boldsymbol{\mathit{P}}$ is the projection orthogonal to the kernel of $\boldsymbol{\mathit{M}}$, i.e. $\boldsymbol{\mathit{P}}=\boldsymbol{\mathit{M}}^{\dagger}\boldsymbol{\mathit{M}}=\boldsymbol{\mathit{M}}\boldsymbol{\mathit{M}}^{\dagger}.$ Then $\boldsymbol{\mathit{M}}^{\dagger}=\boldsymbol{\mathit{P}}\boldsymbol{\mathit{X}}^{-\top}\boldsymbol{\mathit{A}}^{\dagger}\boldsymbol{\mathit{X}}^{-1}\boldsymbol{\mathit{P}}$.

Proof of Lemma D.7.

Note $\boldsymbol{\mathit{B}}^{\top}\boldsymbol{\mathit{B}}=\boldsymbol{\mathit{D}}-\boldsymbol{\mathit{A}}=\boldsymbol{\mathit{L}}$, where $\boldsymbol{\mathit{A}}$ is the adjacency matrix of the graph, and $\boldsymbol{\mathit{L}}$ is its Laplacian. Also \sum_{e}\boldsymbol{\mathit{f}}_{e}^{2}=(\boldsymbol{\mathit{B}}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{B}}^{\top}\boldsymbol{\mathit{B}}}\right)^{{\dagger}}\boldsymbol{\mathit{b}})^{\top}\boldsymbol{\mathit{B}}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{B}}^{\top}\boldsymbol{\mathit{B}}}\right)^{{\dagger}}\boldsymbol{\mathit{b}}=\boldsymbol{\mathit{b}}^{\top}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{B}}^{\top}\boldsymbol{\mathit{B}}}\right)^{{\dagger}}\boldsymbol{\mathit{b}}=\boldsymbol{\mathit{b}}^{\top}\boldsymbol{\mathit{L}}^{{\dagger}}\boldsymbol{\mathit{b}}. By Theorem D.8, we get that for $\boldsymbol{\mathit{x}}\bot\boldsymbol{\mathit{D}}\boldsymbol{\mathbf{1}}$

[TABLE]

Substituting $\boldsymbol{\mathit{y}}=\boldsymbol{\mathit{D}}^{1/2}\boldsymbol{\mathit{x}}$ changes the constraint to $\boldsymbol{\mathit{D}}^{-1/2}\boldsymbol{\mathit{y}}\bot\boldsymbol{\mathit{D}}\boldsymbol{\mathbf{1}}$ i.e. $\boldsymbol{\mathit{y}}\bot\boldsymbol{\mathit{D}}^{1/2}\boldsymbol{\mathbf{1}}$. The inequality now states

[TABLE]

If we let $\boldsymbol{\mathit{Q}}$ denote the projection orthogonal to $\boldsymbol{\mathit{D}}^{1/2}\boldsymbol{\mathbf{1}}$, we can summarize the inequality and orthogonality constraint in one condition using the Loewner order as

[TABLE]

Note that the null space of $\boldsymbol{\mathit{D}}^{-1/2}\boldsymbol{\mathit{L}}\boldsymbol{\mathit{D}}^{-1/2}$ is spanned by $\boldsymbol{\mathit{D}}^{1/2}\boldsymbol{\mathbf{1}}$, as $\boldsymbol{\mathbf{1}}$ spans the null space of $\boldsymbol{\mathit{L}}$. So in fact $\boldsymbol{\mathit{Q}}\boldsymbol{\mathit{D}}^{-1/2}\boldsymbol{\mathit{L}}\boldsymbol{\mathit{D}}^{-1/2}\boldsymbol{\mathit{Q}}=\boldsymbol{\mathit{D}}^{-1/2}\boldsymbol{\mathit{L}}\boldsymbol{\mathit{D}}^{-1/2}$, and we can conclude

[TABLE]

From this we conclude that,

[TABLE]

as $\boldsymbol{\mathit{A}}\succeq\boldsymbol{\mathit{B}}$ implies $\boldsymbol{\mathit{A}}^{\dagger}\preceq\boldsymbol{\mathit{B}}^{\dagger}$ when $\boldsymbol{\mathit{A}}$ and $\boldsymbol{\mathit{B}}$ have the same null space. Hence by Fact D.9 and $\boldsymbol{\mathit{Q}}=\boldsymbol{\mathit{Q}}^{\dagger}$, we then get

[TABLE]

This we can rewrite as for all $\boldsymbol{\mathit{y}}\bot\boldsymbol{\mathit{D}}^{1/2}\boldsymbol{\mathbf{1}}$.

[TABLE]

Substituting $\boldsymbol{\mathit{z}}=\boldsymbol{\mathit{D}}^{1/2}\boldsymbol{\mathit{y}}$ changes the constraint to $\boldsymbol{\mathit{D}}^{-1/2}\boldsymbol{\mathit{z}}\bot\boldsymbol{\mathit{D}}^{1/2}\boldsymbol{\mathbf{1}}$ i.e. $\boldsymbol{\mathit{z}}\bot\boldsymbol{\mathbf{1}}$. Thus we have that for all $\boldsymbol{\mathit{z}}\bot\boldsymbol{\mathbf{1}}$.

[TABLE]

Taking $\boldsymbol{\mathit{z}}=\boldsymbol{\mathit{b}}$, we then get \boldsymbol{\mathit{b}}^{\top}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{B}}^{\top}\boldsymbol{\mathit{B}}}\right)^{{\dagger}}\boldsymbol{\mathit{b}}\leq 2\phi^{-2}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{b}}}\right\|_{\boldsymbol{\mathit{D}}^{-1}}^{2}. ∎

Lemma D.10**.**

Consider a graph $G$ on $n$ vertices with degrees $\boldsymbol{\mathit{D}}$, conductance $\phi$, and edge-vertex incidence matrix $\boldsymbol{\mathit{B}}$, and any demand $\boldsymbol{\mathit{b}}\bot\boldsymbol{\mathbf{1}}$. Define the electrical flow \boldsymbol{\mathit{f}}=\boldsymbol{\mathit{B}}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{B}}^{\top}\boldsymbol{\mathit{B}}}\right)^{{\dagger}}\boldsymbol{\mathit{b}}. Then \mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}}\right\|_{\infty}\leq O(\phi^{-3}\log(n))\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{D}}^{-1}\boldsymbol{\mathit{b}}}\right\|_{\infty}.

Proof.

We first note that if $\boldsymbol{\mathit{f}}^{*}$ is the optimal routing of $\boldsymbol{\mathit{b}}$ in $G$, then

[TABLE]

as per Example 1.4 of [She13]. Secondly, we note that by Lemma 5.2, the electrical flow \boldsymbol{\mathit{f}}^{\mathcal{E}}=\boldsymbol{\mathit{B}}\mathopen{}\mathclose{{}\left({\boldsymbol{\mathit{B}}}^{\top}\boldsymbol{\mathit{B}}}\right)^{{\dagger}}\boldsymbol{\mathit{b}}=\boldsymbol{\mathit{B}}\mathopen{}\mathclose{{}\left({\boldsymbol{\mathit{B}}}^{\top}\boldsymbol{\mathit{B}}}\right)^{{\dagger}}{\boldsymbol{\mathit{B}}}^{\top}\boldsymbol{\mathit{f}}^{*} satisfies

[TABLE]

∎

Lemma D.11**.**

On an expander $G$ with degrees $\boldsymbol{\mathit{D}}$, conductance $\phi$, and gradient $\boldsymbol{\mathit{g}}$ whose projection into the cycle space of $G$, $\boldsymbol{\widehat{\mathit{g}}}$ is $\alpha$-well-spread, for any demand $\boldsymbol{\mathit{b}}\bot\boldsymbol{\mathbf{1}}$ and dot $\theta$ with $\boldsymbol{\widehat{\mathit{g}}}$, the flow given by

[TABLE]

*satisfies *

[TABLE]

Proof.

We first bound the quadratic term $\sum_{e}r\boldsymbol{\mathit{f}}_{e}^{2}$. Let us write \boldsymbol{\mathit{f}}^{\mathcal{E}}=\boldsymbol{\mathit{B}}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{B}}^{\top}\boldsymbol{\mathit{B}}}\right)^{{\dagger}}\boldsymbol{\mathit{b}} and \boldsymbol{\mathit{f}}^{\mathcal{C}}=\frac{\theta}{\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}}\right\|_{2}^{2}}\boldsymbol{\widehat{\mathit{g}}}, and note (in fact, appealing to orthogonality would save an additional factor of 2)

[TABLE]

Then we observe by Lemma D.7 that \sum_{e}(\boldsymbol{\mathit{f}}^{\mathcal{E}}_{e})^{2}\leq\phi^{-2}\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{b}}}\right\|_{\boldsymbol{\mathit{D}}^{-1}}^{2}. Furthermore,

[TABLE]

Combining these equations gives

[TABLE]

We then bound the $p$-th power term,

[TABLE]

Now by Lemma D.10, we have \mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{f}}^{\mathcal{E}}}\right\|_{\infty}=\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{B}}\mathopen{}\mathclose{{}\left(\boldsymbol{\mathit{B}}^{\top}\boldsymbol{\mathit{B}}}\right)^{{\dagger}}\boldsymbol{\mathit{b}}}\right\|_{\infty}\leq O(\phi^{-3}\log n)\mathopen{}\mathclose{{}\left\|\boldsymbol{\mathit{D}}^{-1}\boldsymbol{\mathit{b}}}\right\|_{\infty} and by the $\alpha$-well-spreadness of $\boldsymbol{\widehat{\mathit{g}}}$

[TABLE]

where in the last step we used \mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}}\right\|_{1}\leq m^{1/2}\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}}\right\|_{2}.

∎

Proof of Theorem 4.10..

Refer to the pseudo-code in Algorithm 9. We first collect the facts that we have established about the sampling procedure. In Line 3, $E^{\mathcal{H}}$ is formed from $E^{\mathcal{G}}$ by sampling each edge independently with probability $\tau$. It follows that the expected number of edges in $E^{\mathcal{H}}$ is $\tau m$, and since $\tau>\log{n}/m$, a standard scalar Chernoff bound shows that \mathopen{}\mathclose{{}\left|E^{\mathcal{H}}}\right|\leq 2\tau m with high probability. The parameters $r^{\mathcal{H}}=\tau\cdot r^{\mathcal{G}}$ and $s^{\mathcal{H}}=\tau^{p}\cdot s^{\mathcal{G}}$ are set in Line 4. By Lemma D.2, with high probability the sampling in Line 3 guarantees Equation (18). Note also that

• •

Equation (19) implies $\tau\boldsymbol{\mathit{D}}^{\mathcal{H}}\approx_{0.1}\boldsymbol{\mathit{D}}^{\mathcal{G}}$, by considering the quadratic form in each of the standard basis vectors.

• •

By Corollary D.5, $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}$ is $O(\alpha\phi^{-6}\log^{2}{n})$-well-spread, and \mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}}\right\|_{2}^{2}\approx\tau\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}}\right\|_{2}^{2} and \mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}}\right\|_{1}\approx_{O(\alpha\phi^{-6}\log^{2}{n})}\tau\mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}}}\right\|_{1}.

• •

As $G$ has conductance at least $\phi$, by Corollary D.4, $H$ has conductance at least $0.8\phi$.

We can now establish $\mathcal{G}\preceq_{\kappa}\mathcal{H}$. Suppose $\boldsymbol{\mathit{f}}^{\mathcal{G}}$ is a flow in $G$ with $\boldsymbol{\mathit{B}}^{\mathcal{G}}\boldsymbol{\mathit{f}}^{\mathcal{G}}=\boldsymbol{\mathit{b}}$ and $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}\top}\boldsymbol{\mathit{f}}^{\mathcal{G}}=\theta$. Then $\boldsymbol{\mathit{g}}^{\mathcal{G}\top}\boldsymbol{\mathit{f}}^{\mathcal{G}}=\theta+\boldsymbol{\mathit{\psi}}^{\top}\boldsymbol{\mathit{b}}$. By Lemma D.6, we then get that

[TABLE]

Applying our flow map from $\mathcal{G}$ to $\mathcal{H}$

[TABLE]

We note that by construction, we can readily verify $\boldsymbol{\mathit{B}}^{\mathcal{H}}\boldsymbol{\mathit{f}}^{\mathcal{H}}=\boldsymbol{\mathit{b}}$, $\boldsymbol{\widehat{\mathit{g}}}^{H\top}\boldsymbol{\mathit{f}}^{\mathcal{H}}=\theta$, and $\boldsymbol{\mathit{g}}^{H\top}\boldsymbol{\mathit{f}}^{\mathcal{H}}=\theta+\boldsymbol{\mathit{\psi}}^{\top}\boldsymbol{\mathit{b}}$. So by applying Lemma D.11 to $\frac{1}{\kappa}\boldsymbol{\mathit{f}}^{\mathcal{H}}$, we get

[TABLE]

Our goal is to ensure

[TABLE]

Because the linear terms cancel out, we can use the upper and lower bounds established above to say that this inequality holds provided the following conditions are satisfied (for a $C_{p}$ which is a constant greater than 1 that depends on $p$):

• •

C_{p}\mathopen{}\mathclose{{}\left(\frac{\phi^{-1}}{\kappa}}\right)^{2}\leq 1/\kappa.

• •

$C_{p}\kappa^{-2}\leq 1/\kappa$.

• •

C_{p}\mathopen{}\mathclose{{}\left(\frac{m^{1/p}\phi^{-2}}{\kappa}}\right)^{p}\leq 1/\kappa.

• •

C_{p}\mathopen{}\mathclose{{}\left(\frac{m^{1/p}(\alpha\phi^{-6}\log^{2}{n})^{3/2}}{\kappa}}\right)^{p}\leq 1/\kappa.

Recalling that $\alpha$ is a constant, it follows that there exists a constant $C^{\prime}_{p}$ (depending on p), s.t. all of the above conditions are satisfied, provided

[TABLE]

And this in turn is implied by the stronger condition, $\kappa\geq C^{\prime}_{p}(m^{1/(p-1)}\phi^{-9}\log^{3}n)$, which is hence sufficient to ensure $\mathcal{G}\preceq_{\kappa}\mathcal{H}$.

We can then show $\mathcal{H}\preceq_{\kappa}\mathcal{G}$ with a very similar calculation. We include it for completeness. Suppose $\boldsymbol{\mathit{f}}^{\mathcal{H}}$ is a flow in $H$ with $\boldsymbol{\mathit{B}}^{\mathcal{H}}\boldsymbol{\mathit{f}}^{\mathcal{H}}=\boldsymbol{\mathit{b}}$ and $\boldsymbol{\widehat{\mathit{g}}}^{H\top}\boldsymbol{\mathit{f}}^{\mathcal{H}}=\theta$. Then $\boldsymbol{\mathit{g}}^{H\top}\boldsymbol{\mathit{f}}^{\mathcal{H}}=\theta+\boldsymbol{\mathit{\psi}}^{\top}\boldsymbol{\mathit{b}}$.

By Lemma D.6, we then get that

[TABLE]

Applying our flow map from $\mathcal{H}$ to $\mathcal{G}$

[TABLE]

Again, by construction, we have $\boldsymbol{\mathit{B}}^{\mathcal{G}}\boldsymbol{\mathit{f}}^{\mathcal{G}}=\boldsymbol{\mathit{b}}$, $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{G}\top}\boldsymbol{\mathit{f}}^{\mathcal{G}}=\theta$, and $\boldsymbol{\mathit{g}}^{\mathcal{G}\top}\boldsymbol{\mathit{f}}^{\mathcal{G}}=\theta+\boldsymbol{\mathit{\psi}}^{\top}\boldsymbol{\mathit{b}}$. So by applying Lemma D.11 to $\frac{1}{\kappa}\boldsymbol{\mathit{f}}^{\mathcal{G}}$, we get

[TABLE]

Now, we want to guarantee

[TABLE]

Again the linear terms agree, and termwise verification shows that $\kappa\geq C^{\prime}_{p}(m^{1/(p-1)}\log^{3}(n)\phi^{-9})$ is sufficient to give $\mathcal{H}\preceq_{\kappa}\mathcal{G}$.

∎

Appendix E Using Approximate Projections

Finally, we need to account for the errors in computing the cycle projections $\boldsymbol{\widehat{\mathit{g}}}$ of the gradients $\boldsymbol{\mathit{g}}$. This error arise due to the use of iterative methods in Laplacian solvers used to evaluate $(\boldsymbol{\mathit{B}}^{\top}\boldsymbol{\mathit{B}})^{{\dagger}}$. As we only perform such projections on expanders, we can in fact use iterative methods. However, a dependence of $\log(1/\epsilon)$ in the error $\epsilon$ still remain.

We first formalize the exact form of this error. Kelner et al. [KOSZ13] showed that a Laplacian solver can converge in error proportional to that of the electrical flow. That is, for a slightly higher overhead of $O(\log(n/\epsilon))$, we can obtain a vector $\boldsymbol{\widetilde{\mathit{g}}}$ such that

[TABLE]

This was also generalized to a black-box reduction between solvers for vertex solutions and flows subsequently [CKM*+*14]. As a result, we will work this guarantee with errors relative to $\boldsymbol{\widehat{\mathit{g}}}$.

For the partitioning stage, this error occurs in two places: for computing the norm of the projection, and for identifying edges with high contributions (aka. non-uniform) for removal.

For the former, a constant factor error in the norm of $\boldsymbol{\widehat{\mathit{g}}}$ will only lead to a constant factor increase in:

The uniformity of the true projected gradient, 2. 2.

The factor of decrease in the norm of the projected gradient from one step to next.

For both of these, such constant factor slowdowns can be absorbed by an increase in the thresholds, which in turn result in a higher uniformity parameter in decompositions returned. As this uniformity parameter only affects the number of edges sampled in Theorem 4.10, they only accumulate to a larger overhead in the $m^{O(\nicefrac{{1}}{{\sqrt{p}}})}$ term in the overall running time.

The other invocation of projections is in the sparsification of expanders in Algorithm 9. Here the decomposition of $\boldsymbol{\mathit{g}}^{\mathcal{G}}$ into a circulation and potential flows is necessary for the construction of the gradient of the sampled graph, $\mathcal{H}$.

While an approximate energy minimizing circulation $\boldsymbol{\widetilde{\mathit{g}}}$ will not have $\boldsymbol{\mathit{g}}-\boldsymbol{\widetilde{\mathit{g}}}$ being a potential flow, we can instead perturb $\boldsymbol{\mathit{g}}$ slightly in this instance. Specifically, we can also compute a set of approximate potentials $\boldsymbol{\widetilde{\mathit{\psi}}}$ so that

[TABLE]

That is, we can perturb the initial $\boldsymbol{\mathit{g}}$ based on the result of this solve so that we have an exact decomposition of it into a circulation and a potential flow. The error of this perturbation is then incorporated in the same manner as terminating when \mathopen{}\mathclose{{}\left\|\boldsymbol{\widehat{\mathit{g}}}}\right\|_{2} is too small in Case 2b of Theorem 4.9. Specifically, the additive error of this goes into the additive trailing terms of the guarantees of the ultra-sparsifier shown in Theorem 3.6.

Finally, the projection of the sampled gradient $\boldsymbol{\widetilde{\mathit{g}}}^{\mathcal{H}}$ into $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}$ also carries such an error term. By picking $\epsilon$ to be in the $1/poly(n)$ range, we ensure that both the $\ell_{2}$ and $\ell_{1}$ norms of $\boldsymbol{\widehat{\mathit{g}}}^{\mathcal{H}}$ is close to their true terms. This in turn leads to constant factor errors in the lower and upper bounds on objectives give in Lemmas D.6 and D.11, and thus a constant factor increase in the overall approximation factors.

Therefore, it suffices to set $\epsilon$ in these approximate projection algorithms to be within $poly(n)$ factors of the $\delta$ by which UltraSparsify is invoked by the overall recursive preconditioning scheme. The choice of parameters in Theorem 3.7 then gives that it suffices to have $\log(1/\epsilon)\leq\widetilde{O}(1)$ in all projection steps. In other words, all the projections can be performed in time nearly-liner in the sizes of the graphs.

Appendix F $\ell_{p}$-norm Semi-Supervised Learning on Graphs.

In this appendix, we briefly describe how to convert Problem (3), into a form that can be solved using our algorithm for smoothed $p$-norm flows as stated in Theorem 1.1.

Recall that formally, given a graph $G=(V,E)$ and a labelled subset of the nodes $T\subset V$ with labels $\boldsymbol{\mathit{s}}_{T}\in\mathbb{R}^{T}$, we can write the problem as

[TABLE]

Taking a Lagrange dual now results in the problem

[TABLE]

where $q=\frac{1}{1-1/p}$, and the gradient $\boldsymbol{\mathit{g}}$ is given by $\boldsymbol{\mathit{g}}=\boldsymbol{\mathit{B}}_{:,T}\boldsymbol{\mathit{s}}_{T}$. We cannot directly solve this formulation, since the net incoming flow at vertices in $T$ is unknown. However, notice that the flow is preserved at all other vertices, so summed across all of $T$, the net flow must be zero. Thus if we merge all the vertices in $T$ into one vertex, while turning edges in $T\times T$ into self-loops, the problem is now a circulation. Note that the optimal flow on each self-loop can be computed exactly. Now the resulting problem can be solved to high accuracy using Theorem 1.1. Meanwhile, mapping the flow back to the original flow, it can be shown that the optimal flow arises as a simple non-linear function of some voltages $\boldsymbol{\mathit{x}}$: $\boldsymbol{\mathit{f}}_{e}=(\boldsymbol{\mathit{B}}\boldsymbol{\mathit{x}})_{e}^{q}$. This means that if we have $\boldsymbol{\mathit{f}}$ to high enough accuracy, we can get an almost optimal set of voltages, e.g. by looking at flow along edges of a tree to compute a set of voltages $\boldsymbol{\mathit{x}}$ that are a $(1+1/{\textrm{poly}}(m))$ multiplicative accuracy solution to Problem (3). Since we call the algorithm of Theorem 1.1 using to solve a $\frac{p}{p-1}$-flow problem, where $p<2$, the running time will be on the order of $2^{O((\frac{p}{p-1})^{\nicefrac{{3}}{{2}}})}m^{1+O(\sqrt{\frac{p-1}{p}})}$. This in turn can be further simplified as $2^{O((\frac{1}{p-1})^{\nicefrac{{3}}{{2}}})}m^{1+O(\sqrt{p-1})},$ since $p<2$. For $p=1+\frac{1}{\sqrt{\log n}}$, this is time is bounded by $m^{1+o(1)}$.