Computing Flows in Subquadratic Space

We present the following result for computing min-cost flows in the multi-pass streaming model.

While the streaming area often focuses on small number of passes (e.g., constant or polylog), it was proven that polynomial number of passes may be necessary for flows in the weighted case. In the presence of capacities as large as sub-exponential in $n$, [6] showed that a streaming algorithm for computing the size of the max-flow in $p$ passes needs $\mathrm{\Omega }(n^2/p^5)$ space. This lower bound was later improved to $\mathrm{\Omega }(n^2/p^3)$ space [8]. For smaller polynomially-sized capacities, Theorem 1 is the first upper bound for computing the exact size of the max-flow in the weighted case (when interested in the size, we can simply round the result to nearest integer). It is also the first result to return the flow $f\in {ℝ}^E$ with high accuracy rather than the cut or cost/size of the flow.

Previous work was either highly inaccurate with polynomial error when returning the flow [22], or when exact, previous work was only for the unweighted undirected case [58] and returned only the cut but not the flow. Observe that for unweighted graphs, any max flow uses at most $O(n^{1.5})$ edges (see full version for details), but for weighted graphs the number of used edges can be as large as $O(n^2)$. So minimizing space complexity is generally harder for the weighted case.

Our algorithm is based on solving the linear program representation of min-cost flow. The linear program representation of related problems such as bipartite matching or transshipment to obtain streaming algorithms has been used before, e.g., in [9, 51, 2]. A natural approach for matching is to maintain the smaller $n$-dimensional dual iterate (i.e., fractional vertex cover). However, this does not extend to min-cost/maximum flows, because the capacity constraints increase the dimension of the dual solution from $n$ to $m+n$.

Our algorithm works by implementing the interior point method/central path method by [19], while maintaining the primal solution (i.e., flow) in low space, rather than the dual. One can show that for flows on the central path, an $\tilde{O}(n)$ space representation is possible³³3Technical details: Centered flows are minimizers of some potential function $\mathrm{\Phi }(x)$ subject to ${𝐀}^{\top }x=b$ for edge-vertex incidence matrix $𝐀$ and demand $b\in {ℝ}^V$. Minimizers satisfy $\nabla \mathrm{\Phi }(x)=𝐀y$ for some $y\in {ℝ}^V$, thus by storing $y$ we can reconstruct $x$.. Unfortunately, this observation does not directly imply a streaming algorithm because in each iteration, the central path method constructs intermediate flows that are not sufficiently centered to be stored in $\tilde{O}(n)$ space. This issue is exacerbated when using the fast $\tilde{O}(\sqrt{n})$ iteration central path of the Lee-Sidford-barrier [50] as used in [19]. Their definition of centrality uses Lewis-weights whose efficient computation uses sketching. Due to the resulting approximation error, one cannot compute sufficiently centered flows to store them in $\tilde{O}(n)$ space. Thus we develop a method that can store non-centered flows in small space.

We do this by storing the steps of each iteration rather than the flow $x$. In the context of flows, each step of the central path method is equivalent to augmenting the current flow $x$ by an electric circulation. Rather than storing flow $x$ directly, we store the sequence of $\tilde{O}(\sqrt{n})$ augmenting electric flows. Each circulation is an $m$-dimensional vector (i.e., every edge carries some flow), but we show that each circulation can be stored in only $\tilde{O}(n)$ space. Thus we obtain a $\tilde{O}(n^{1.5})$ space representation of the min-cost flow.
The main technical ideas/ingredients for storing electric circulations are as follows.

1. (i)

   The $\tilde{O}(\sqrt{n})$ iterations framework [50, 19] requires computation of Lewis-weights, a generalization of leverage scores and effective resistances. Via sketching, we can store the Lewis-weights in only $\tilde{O}(n)$ space. Given any of the $O(n^2)$ edges, we can then query/compute the (approximate) Lewis-weight from only the $\tilde{O}(n)$ space representation.

2. (ii)

   While electric flows can be stored via vertex potentials (i.e., in $\tilde{O}(n)$ space), electric circulations generally need $\mathrm{\Omega }(m)$ space. That is because electric circulations are a combination/difference $g-f$ of an arbitrary flow $g\in {ℝ}^E$ and an electric flow $f\in {ℝ}^E$, where $f$ and $g$ satisfy the same vertex demands. To store the electric circulation $g-f$ in low space, we must store $g$ in low space. Via the robust interior point framework [19], one can reduce the dimension of $g$ from $m$ to $\tilde{O}(1)$, and can thus be stored in small space. The main barrier is showing that this projection $E\to \{1,\mathrm{\dots }.,\tilde{O}(1)\}$ can be computed from the stored information.

3. (iii)

   Electric flows can only be stored via vertex potentials if we know the edges’ resistances. Retrieving the flow on edge $e=(u,v)$ from vertex potentials $y\in {ℝ}^V$ (i.e., $f_e=(y_v-y_u)/r_e$) requires knowing the edge’s resistance $r_e$. In [50, 19], the resistance $r_e$ is defined w.r.t. the edge’s Lewis-weights, and the current congestion of that edge. This allows for an inductive argument: Given an edge $e$ and its capacity $c_e$, assuming we know the amount of flow on edge $e$ after $t$ electric circulations, then we can compute the resistance $r_e$ of that edge from the current congestion and the stored Lewis-weight. From the resistance, we can then obtain how much flow the next electric circulation will route through $e$. This then yields the flow on edge $e$ after $t+1$ electric circulations. After $\tilde{O}(\sqrt{n})$ iterations of this argument, we know the final amount of flow on that edge, i.e., how much flow is routed on $e$ by the min-cost flow.

Before the semi-streaming lower bounds of $\mathrm{\Omega }(n^2/p^5)$ [6] and $\mathrm{\Omega }(n^2/p^3)$ [8] for $p$ passes, there have been $n^{1+\mathrm{\Omega }(1/p)}$ space lower bounds [21, 39]. For small number $o(\sqrt{\log n})$ passes, there is an $n^{2-o(1)}$ space lower bound for $st$-reachability, which extends to directed flows [25].

While our focus is on $st$-flows and thus $st$-min cuts, there has been work on global min cuts in the streaming model. There are $\tilde{O}(1)$ pass $\tilde{O}(n)$ space algorithms for global min cuts [56, 7, 58]. These also transfer to the 2-party communication setting, where $\mathrm{\Theta }(n^{1.5})$ upper and lower bound are known for minimum vertex cut [17].

Our algorithm is based on solving linear programs in the streaming model. Other work using this approach (though not for flows, but instead for other linear programs such as bipartite matching) include [51, 2, 10, 23]. From the communication perspective, linear programs have been studied in [59, 36]. The crucial aspect is that previous work studies linear programs of form $\min c^{\top }$ subject to ${𝐀}^{\top }x=b,x\ge 0$ for tall matrices $𝐀$, but max-flow requires additional capacity constraints $x\le u$. Encoding those within the matrix $𝐀$ turns the matrix from tall into an almost square matrix, thus losing any space/communication gains from the small width of $𝐀$.

Here we outline how to store a flow $x\in {ℝ}^E$ in low-space, such that when given any edge $e\in E$, its cost $c_e$, and capacity $u_e$, then we can query/compute $x_e$.

Consider the linear programming definition of minimum cost flow. For an $E\times V$ incidence matrix $𝐀\in {\{-1,0,1\}}^{E\times V}$, edge-cost vector $c\in {ℝ}^E$, demand vector $b\in {ℝ}^V$, and edge capacities $u\in {ℝ}^E$, the following linear program models minimum-cost flows.

A common idea for computing minimum cost flows is to start with some initial flow, then repeatedly augment the flow via negative cost cycles or circulations. When solving minimum cost flows via interior point methods such as [52, 50, 53, 13], these circulations are so called “electric circulations.” That is, for some feasible flow $x\in {ℝ}^E$, we augment it by a circulation ${\delta }_x\in {ℝ}^E$ where

Here $g_x\in {ℝ}^E$ is some flow and ${𝐑}_x$ is an $E\times E$ diagonal matrix with ${({𝐑}_x)}_{e,e}$ being resistance of edge $e$. (Both $g_x$ and ${𝐑}_x$ will be defined later to depend on $x$). Flow ${\delta }_x$ is a circulation, because it satisfies 0-demand ${𝐀}^{\top }{\delta }_x=0$. It is an electric circulation because

is an electric flow for demand $d=-{𝐀}^{\top }{g}_x\in {ℝ}^V$ and resistances ${𝐑}_x$ (i.e., among all demand $d$ flows, $f$ minimizes the energy ${\sum }_{e\in E}{({𝐑}_x)}_{e,e}{f}_e^2$).

The idea of our low space representation is that we only store some initial flow $x^0\in {ℝ}_{\ge 0}^E$, and then store all the augmenting circulations ${\delta }_x^1,{\delta }_x^2,\mathrm{\dots }$. We will argue that each ${\delta }_x^t\in {ℝ}^E$ can be stored in only $\tilde{O}(n)$ space. It was proven in [50] that it takes $\tilde{O}(\sqrt{n}\log (W/ϵ))$ iterations/augmentations to find a min-cost flow whose cost is at most an additive $ϵ$ off from optimal. Thus $\tilde{O}(n^{1.5}\log (W/ϵ))$ total space suffices⁴⁴4For the rest of the overview, we will hide the $\log (W/ϵ)$ term in $\tilde{O}$-notation.. We can then query the amount of flow on any edge $e\in E$ by computing $x_e^0+{\sum }_{t=1}^{\tilde{O}(\sqrt{n})}{({\delta }_x^t)}_e$.

Crucial for our low space representation is that the $(t+1)$-th augmenting circulation ${\delta }_x^{t+1}$ depends on the $t$-th flow $x^t=x^0+{\sum }_{k=1}^t{\delta }_x^k$. The fact that we can store circulation ${\delta }_x^{t+1}$ in only $\tilde{O}(n)$ space relies heavily on the following property: when querying/computing ${({\delta }_x^{t+1})}_e$ for some edge $e$, we already know the flow $x_e^t$ on that edge. This assumption is given by induction: if we know $x_e^{t-1}$, then we can compute ${({\delta }_x^t)}_e$, thus we know $x_e^t=x_e^{t-1}+{({\delta }_x^t)}_e$, and so forth. This way we iteratively reconstruct $x_e^t$ and ${({\delta }_x^t)}_e$ from $t=1,2,3,\mathrm{\dots }$ all the way to the final/optimal flow on that edge $e$.

To outline these low-space representations, let us start with the initial flow $x^0$.

As context and motivation for our method, we start by discussing the method of [50]. Readers familiar with the $\tilde{O}(\sqrt{n})$ iteration interior point methods from [50, 19] can skip to Section 3.3.

For matrices $𝐀\in {ℝ}^{m\times n}$, vectors $b\in {ℝ}^n$, $c\in {ℝ}^m$, and upper bounds $u\in {ℝ}^m$, this IPM solves linear programs of the form

The central path method of [50] repeatedly solves the following minimization problem. In each iteration, parameter $\mu \in {ℝ}_{>0}$ is decreased, and then the new minimizer $x_{\mu }$ is found by performing one Newton-step from the old minimizer $x_{{\mu }^{\prime }}$.

where $\tau \in {ℝ}_{>0}^m$ and ${\varphi }_e:{ℝ}_{>0}\to ℝ$ is defined as ${\varphi }_e(z)=-\log (u_e-z)-\log (z)$. (For simplicity, we also write $\varphi (x)\in {ℝ}^m$ for the vector with entries ${\varphi }_e(x_e)$ for $e=1,..,m$.) Note that as $x_e$ approaches 0 or capacity $u_e$, the ${\varphi }_e(x_e)$ term goes to $\mathrm{\infty }$, which keeps the minimizer $x_{\mu }$ within the feasible region. As $\mu \to 0$, the term $c^{\top }x$ starts to dominate; thus $x_{\mu }$ converges towards the optimal solution of the linear program for $\mu \to \mathrm{\infty }$. Here $\tau \in {ℝ}^m$ is a weight measuring the importance of each row of $𝐀$. [50] chooses $\tau$ to be the regularized $l_p$ Lewis-weight for $p=1-\frac{1}{4\log (4m/n)}$. This is defined to be the solution to

where ${\mathrm{\Phi }}^{\prime \prime }(x)$ is a diagonal matrix with the 2nd derivatives ${\varphi }_e^{\prime \prime }(x_e)=1/x_e^2-1/{(u_e-x_e)}^2$ on the diagonal for $e=1,\mathrm{\dots },m$.

The curve $x_{\mu }$ for $\mu \in (\mathrm{\infty },0)$ is referred to as the central path. It turns out that it is easy to find a point $x$ close to $x_{\mu }$ on this path for large $\mu$ (see full version for details). Starting from this point, we iteratively take steps to decrease $\mu$ and move $x$ closer to the new minimizer $x_{\mu }$.

To do this, we must define a measure for how far $x$ is from $x_{\mu }$. To satisfy $x={argmin}_{{𝐀}^{\top }x=b}{f}_{\mu }(x)$ (Equation 5), optimality conditions tell us that we need ${𝐀}^{\top }x=b$ and the gradient $\nabla f_{\mu }(x)=c+\mu 𝐓{\varphi }^{\prime }(x)$ must be orthogonal to the feasible hyperplane $\{x\mid {𝐀}^{\top }x=b\}$; i.e., $c+𝐀y+\mu 𝐓{\varphi }^{\prime }(x)=0$ for some $y\in {ℝ}^n$. We can rewrite this as $s=\mu 𝐓{\varphi }^{\prime }(x)=0$, where $s:=c+𝐀y$. Now, to define what it means to follow the central path, we say that a triple $(x,s,\mu )$ is central if

As we cannot expect to be exactly centered, we measure centrality via the following potential function:

where $\lambda =O(\log n)$. Note that the numerator is 0 (which is equivalent to minimizing (8)) if and only if (7) is satisfied.

Our goal is now to take steps where we alternate between decreasing $\mu$ (which increases (8)) and updating $(x,y,s)$ such that (8) decreases again. This way we approximately trace the curve $x_{\mu }$.

In the central path method by [19], the improvements in each iterations are of the following form. Given $(x^{t-1},s^{t-1},{\mu }^{t-1})$ at iteration number $t-1$, the next values for iteration $t$ are given by

To compute the movements ${\delta }_x^t,{\delta }_s^t$ we first calculate the weights ${\tau }^t\in {ℝ}^m$ of our current point $x^{t-1}$, and set the vector $g^t\in {ℝ}^m$ to be the gradient of (8), which intuitively moves $x$ as close as possible to $x_{\mu }$ on the central path. [19] have shown that instead of exact calculation of ${\tau }^t$, it suffices to compute weight ${\overline{\tau }}^t$ that satisfy recurrence (6) only with some $(1\pm 1/O(\log n))$-approximation. Then, the movements are given by

where ${𝐇}^t$ is a spectral approximation, allowing us to solve the linear system approximately (i.e., we can use fast Laplacian system solvers). Because of this approximation, we can no longer guarantee ${𝐀}^{\top }x=b$, but inclusion of vector ${\delta }_c^t$ guarantees ${𝐀}^{\top }x$ stays close to $b$. Here $g^t$ is essentially the gradient of our potential function $\mathrm{\Psi }$ (since we want to reduce that potential in each iteration), but we are taking a step that is bounded in a certain $(\tau +\mathrm{\infty })$-norm⁵⁵5${\Vert v\Vert }_{\tau +\mathrm{\infty }}:={\Vert v\Vert }_{\mathrm{\infty }}+C\log (4m/n)\cdot {\Vert v\Vert }_{\tau }$ where ${\Vert v\Vert }_{\tau }={({\sum }_i{\tau }_i\cdot v_i^2)}^{1/2}$ and $C$ is some large constant. (see full version for details). Calculation of this $g^t$ in low space will later be one of the main issues discussed in Section 3.3.

In [50] it was shown that $\tilde{O}(\sqrt{n})$ iterations of the central path method (9) suffice, i.e. one must compute (9) for $t=1,2,\mathrm{\dots },\tilde{O}(\sqrt{n})$. The naive implementation of (9) takes $\mathrm{\Omega }(m)$ space which is $\mathrm{\Omega }(n^2)$ in worst-case, since the vectors $x^t,s^t,g^t\in {ℝ}^m$ are $m$-dimensional. However, $O(n^2)$ space for min-cost flow is trivial, because then we could solve the problem in a single pass by storing the entire input graph. Therefore, our goal is to implement (9) using less than $O(n^2)$ space.

In the following, we outline how we can implement all calculations in (9) in the semi-streaming model. In particular, we show that for any one $t$, (9) can be calculated in just $\tilde{O}(1)$ passes over the input. Since we have $t=1,2,\mathrm{\dots },\tilde{O}(\sqrt{n})$, our algorithm takes $\tilde{O}(\sqrt{n})$ passes in total.

Our goal is to calculate (9) a total of $\tilde{O}(\sqrt{n})$ times in a semi-streaming setting with $\tilde{O}(1)$ passes over the input per iteration, and using only $\tilde{O}(n^{1.5})$ total space.

We already outlined in Section 3.1 that we can store the values of the $t$-th iteration in total space $\tilde{O}(nt)$ by storing each iteration in $\tilde{O}(n)$ space. The precise values being stored are:

• $\mathrm{■}$

  $y^{\mathrm{\ell }},{\delta }_y^{\mathrm{\ell }},{\delta }_c^{\mathrm{\ell }}\in {ℝ}^n$ for $\mathrm{\ell }=1,\mathrm{\dots },t.$ This uses $O(nt)$ space.
  This implicitly represents $s^{\mathrm{\ell }}\in {ℝ}^m$ because for any edge $e=(w,z)$ and its cost $c_e$, we can compute $s_e^{\mathrm{\ell }}=c_e-{(𝐀y^{\mathrm{\ell }})}_e=c_e-(y_z^{\mathrm{\ell }}-y_w^{\mathrm{\ell }})$.

• $\mathrm{■}$

  We store ${\tau }^{\mathrm{\ell }}$ for $\mathrm{\ell }=1,\mathrm{\dots },t$ implicitly in $\tilde{O}(nt)$ space via a certain recursive JL-sketch.

• $\mathrm{■}$

  We store $g^t$ for $\mathrm{\ell }=1,\mathrm{\dots },t$ implicitly in $\tilde{O}(t)$ space. (Details to be discussed further below).
  This implicitly represents $x^{\mathrm{\ell }}\in {ℝ}^m$ because for any edge $e=(w,z)$ and its cost/capacity $c_e,u_e$, we can recursively compute

  $$x_e^{\mathrm{\ell }}=x_e^{\mathrm{\ell }-1}+{\varphi }_e^{\prime \prime }{(x_e^{\mathrm{\ell }-1})}^{-1/2}{g}_e^{\mathrm{\ell }}-{({\overline{\tau }}_e^{\mathrm{\ell }})}^{-1}{\varphi }_e^{\prime \prime }(x_e^{\mathrm{\ell }-1}){(𝐀({\delta }_y^{\mathrm{\ell }}+{\delta }_c^{\mathrm{\ell }}))}_e.$$