# Derandomizing Directed Random Walks in Almost-Linear Time

Rasmus Kyng,
ETH Zurich
<EMAIL_ADDRESS>
The research leading to these results has received funding from grant no.
200021 204787 of the Swiss National Science Foundation. Simon
Meierhans††footnotemark: ,
ETH Zurich
<EMAIL_ADDRESS>Maximilian Probst Gutenberg††footnotemark: ,
ETH Zurich
<EMAIL_ADDRESS>
###### Abstract

In this article, we present the first deterministic directed Laplacian
$\boldsymbol{\mathit{L}}$ systems solver that runs in time almost-linear in
the number of non-zero entries of $\boldsymbol{\mathit{L}}$. Previous
reductions imply the first deterministic almost-linear time algorithms for
computing various fundamental quantities on directed graphs including
stationary distributions, personalized PageRank, hitting times and escape
probabilities.

We obtain these results by introducing _partial symmetrization_ , a new
technique that makes the Laplacian of an Eulerian directed graph “less
directed” in a useful sense, which may be of independent interest. The
usefulness of this technique comes from two key observations: Firstly, the
partially symmetrized Laplacian preconditions the original Eulerian Laplacian
well in Richardson iteration, enabling us to construct a solver for the
original matrix from a solver for the partially symmetrized one. Secondly, the
undirected structure in the partially symmetrized Laplacian makes it possible
to sparsify the matrix _very crudely_ , i.e. with large spectral error, and
still show that Richardson iterations convergence when using the sparsified
matrix as a preconditioner. This allows us to develop deterministic
sparsification tools for the partially symmetrized Laplacian.

Together with previous reductions from directed Laplacians to Eulerian
Laplacians, our technique results in the first deterministic almost-linear
time algorithm for solving linear equations in directed Laplacians. To
emphasize the generality of our new technique, we show that two prominent
existing (randomized) frameworks for solving linear equations in Eulerian
Laplacians can be derandomized in this way: the squaring-based framework of
[CKP+17] and the sparsified Cholesky-based framework of [PS22].

##  1 Introduction

The development of spectral graph sparsification and nearly-linear time
solvers for Laplacian linear equations initiated by the seminal article of
Spielman and Teng [ST04], that has since been split into three parts [ST13,
ST11, ST14], is foundational for algorithmic spectral graph theory and one of
the success stories in the design of graph algorithms. Until recently,
spectral techniques were mainly used for analysing undirected graphs. While
many of the techniques developed in algorithmic spectral graph theory heavily
use that Laplacians of undirected graphs are symmetric, positive semi-definite
matrices, Cohen-Kelner-Peebles-Peng-Sidford-Vladu [CKP+16a] first demonstrated
that the structure of directed Laplacians can be used to accelerate linear
equation solvers. Together with Rao they later introduced the first notion of
spectral approximation for directed Laplacians [CKP+17], which enabled them to
develop sparsification tools in the directed setting. These tools were used to
build the first almost-linear time directed Laplacian solver, operating in the
framework of an undirected Laplacian solver by Peng and Spielman [PS14]. The
runtime was improved to nearly-linear time by Cohen-Kelner-Kyng-Peebles-Peng-
Rao-Sidford [CKK+18] and further improved by Peng and Song [PS22] very
recently. Both operate within sparsified-Cholesky frameworks for solving
undirected Laplacian linear equations: [CKK+18] uses the framework of Kyng and
Sachdeva [KS16] and [PS22] uses the framework of Kyng-Lee-Peng-Sachdeva-
Spielman [KLP+16].

All previous fast algorithms for solving directed Laplacian linear equations
rely on sampling for globally sparsifying directed graphs while retaining a
spectral $(1\pm\epsilon)$-approximation guarantee for some $\epsilon<1$,
according to the notion of approximation by [CKP+17]111When we refer to
$(1\pm\epsilon)$-approximations for Eulerian Laplacians in the introduction
and overview we mean $\epsilon$-approximations as in Definition 3.1. We use
this naming for convenience when comparing to undirected approximations.. This
suggests two main approaches to derandomizing these solvers: (1) developing a
fast deterministic $(1\pm\epsilon)$-approximate spectral sparsification
routine or (2) introducing cruder forms of approximation and adapting the
algorithms to cope with weaker guarantees. However, even for undirected
graphs, no deterministic almost-linear time algorithms that achieve
$(1\pm\epsilon)$-approximate spectral sparsification are known, and this is a
major obstacle to approach (1). We circumvent this issue by instead taking the
route (2): we develop a new way of measuring approximation via the suitability
as a preconditoner in Richardson, which allows cruder guarantees. We achieve
this by introducing a “robustification” step before sparsification, which we
call $\beta$-partial symmetrization. Partial symmetrization counteracts the
fragility of directed approximations and is at the core of our new crude
deterministic sparsification procedure for Eulerian Laplacians. It lets us
derandomize the frameworks of [CKP+17] and [PS22] with an almost-linear
runtime.

###  1.1 Prior Work

#### Undirected Laplacian Solvers.

The first nearly-linear time Laplacian solver [ST13, ST11, ST14] sparked the
development of a field producing a whole host of distinct algorithms for
solving Laplacian linear equations [KMP14, KMP11, KOSZ13, CKP+14, PS14,
KLP+16, KS16, JS21]. Of these, our algorithm is most similar to [PS14] and its
directed counterpart [CKP+17], which works by repeatedly squaring the
adjacency matrix. In Appendix C we consider a Cholesky-factorisation based
framework, which is motivated by [KLP+16] in the undirected setting and was
recently translated to the directed setting by [PS22].

#### Spectral Sparsification.

A crucial building block for spectral graph algorithms, including linear
equation solvers, is the ability to sparsify undirected graphs while
preserving their spectral properties. This is a stronger notion of
sparsification than cut-sparsifiers, which only preserve the approximate size
of cuts. Such spectral sparsifiers were first introduced in [ST11], and later
strengthened and simplified by Spielman and Srivastava [SS11]. It is known
that for every $n$-vertex graph there exists a spectral sparsifier with
$O(n)$-edges, and such sparsifiers can be constructed deterministically, as
shown by Batson, Spielman, and Srivastava [BSS12]. However, no deterministic
almost-linear time algorithms are known for $(1\pm\epsilon)$-spectral
sparsification. Recently, Chuzhoy-Gao-Li-Nanongkai-Peng-Saranurak [CGL+20a]
presented an almost-linear time algorithm achieving $n^{o(1)}$-spectral
sparsifiers for undirected graphs via deterministic expander decompositions.
Our algorithm relies on both their sparsification and their expander
decomposition results.

#### Directed Laplacian Solvers.

Before [CKP+16a, CKP+17] it was unclear that directed Laplacians, which are a
natural generalization of undirected Laplacians to directed graphs, also
exhibit properties that allow for useful notions of sparsification and/or
accelerated solving of linear equations. The reduction to strongly connected
Eulerian Laplacians recovered some of the spectral properties of undirected
Laplacians, and allowed for the development of a useful notion of
sparsification. Current fast algorithms for solving Eulerian Laplacian linear
equations either follow a squaring [CKP+17] or a sparsified-Cholesky approach
[CKK+18, PS22]. Both rely on spectral sparsification techniques developed in
[CKP+17].

#### Low Space Algorithms.

Recently, the first deterministic $\tilde{O}(\log N)$-space solver for
Eulerian Laplacians was introduced by Ahmadinejad-Kelner-Murtagh-Peebles-
Sidford-Vadhan [AKM+20]. They show that the Rozeman and Vadhan [RV05]
deterministic squaring conserves approximation under squaring for a new,
stronger measure of approximation. The directed to Eulerian reduction remains
a major obstacle to solving directed Laplacian linear equations in small
space.

#### Applications of Directed Laplacian Solvers.

There are numerous applications of directed Laplacian solvers given in Section
7 of [CKP+16c], most of which are deterministic reductions to directed
Laplacian system solving. The deterministic ones include

  * •

solving large classes of linear systems,

  * •

computing personalised PageRank vectors,

  * •

estimating the stationary distribution,

  * •

and simulating random walks.

Since the reductions are deterministic, we obtain deterministic almost-linear
time algorithms for all these problems222Note that we require an upper bound
$\kappa$ on the condition number, and hence mixing time..

###  1.2 Our Contributions

We introduce the first notion of crude approximation for Eulerian Laplacians
$\boldsymbol{\mathit{L}}$. It is defined via the suitability as a
preconditioner in the Richardson iteration and sparse approximations are
constructed using _partial symmetrization_ to increase robustness. This
technique allows us to trade off additional Richardson iterations for a
behaviour more akin to the undirected setting. The obtained _deterministic_
crude global sparsification routine ultimately allows us to derandomize two
prominent frameworks for solving directed Laplacian linear equations. We
summarize our main result assuming polynomially bounded edge weights and
condition number.

######  Theorem 1.1 (Informal Version of Theorem 7.1).

Given an Eulerian Laplacian $\boldsymbol{\mathit{L}}_{\underrightarrow{G}}$
associated with a strongly connected Eulerian Graph $\underrightarrow{G}$ with
$n$ vertices and $m$ edges, a vector
$\boldsymbol{\mathit{b}}\in\operatorname*{im}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}})$,
and a parameter $\epsilon\in(0,1)$ the algorithm
$\textsc{SolveEulerian}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}},\epsilon)$
in time $m^{1+o(1)}\log\epsilon^{-1}$ computes a vector
$\boldsymbol{\mathit{x}}$ satisfying

$\displaystyle\left\|\boldsymbol{\mathit{x}}-\boldsymbol{\mathit{L}}_{\underrightarrow{G}}^{+}\boldsymbol{\mathit{b}}\right\|_{\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}}\leq\epsilon\left\|\boldsymbol{\mathit{L}}_{\underrightarrow{G}}^{+}\boldsymbol{\mathit{b}}\right\|_{\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}}$

where
$\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{A}}}=(\boldsymbol{\mathit{A}}+\boldsymbol{\mathit{A}}^{T})/2$
for any square matrix $\boldsymbol{\mathit{A}}$.

Previous reductions allow us to reduce general directed Laplacian solvers to
$\log^{O(1)}(n\kappa^{-1}\epsilon^{-1})$ Eulerian solves with polynomially
bounded condition number and edge weights. We state the main theorem as a
corollary.

######  Theorem 1.2.

Given a directed Laplacian
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}}=\boldsymbol{\mathit{D}}_{\underrightarrow{G}}-\boldsymbol{\mathit{A}}_{\underrightarrow{G}}^{T}$
associated with a directed Graph $\underrightarrow{G}$ with $n$ vertices and
$m$ edges, a vector
$\boldsymbol{\mathit{b}}\in\operatorname*{im}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}})$,
a parameter $\epsilon\in(0,1)$ and an upper bound
$\kappa\geq\max(\kappa(\boldsymbol{\mathit{D}}_{\underrightarrow{G}}),\kappa(\boldsymbol{\mathit{L}}_{\underrightarrow{G}}))$
there is an algorithm
$\textsc{SolveFull}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}},\epsilon)$
that in time $m^{1+o(1)}\log^{O(1)}(\kappa\epsilon^{-1})$ computes a vector
$\boldsymbol{\mathit{x}}$ satisfying

$\displaystyle\left\|\boldsymbol{\mathit{x}}-\boldsymbol{\mathit{L}}_{\underrightarrow{G}}^{+}\boldsymbol{\mathit{b}}\right\|_{2}\leq\epsilon\left\|\boldsymbol{\mathit{L}}_{\underrightarrow{G}}^{+}\boldsymbol{\mathit{b}}\right\|_{2}$

where
$\kappa(\boldsymbol{\mathit{A}})=\left\|\boldsymbol{\mathit{A}}\right\|_{2}\cdot\left\|\boldsymbol{\mathit{A}}^{+}\right\|_{2}$
denotes the condition number of a given matrix $\boldsymbol{\mathit{A}}$.

A key application of these results is the deterministic simulation of random
walks in almost-linear time as presented in Sections 7.4 and 7.5 of [CKP+16c].

#### Crude Global Sparsification.

Our global sparsification routine is based on first increasing the robustness
via $\beta$-partial symmetrization. Then we bucket by edge weight and layer
partitions of undirected and unweighted graphs into expanders. First, the
directed part in such an expander can be sparsified in a crude way that
conserves the degrees. Then, the undirected part can be replaced with a sparse
expander. Since the undirected part is scaled with a factor $\beta$ by
$\beta$-partial symmetrization it dominates the directed part which allows us
to bound the error of the first step.

#### Open Questions.

Our techniques, in particular $\beta$-partial symmetrization, might be of
interest for deterministically simulating random walks in small space.
Further, the development of a composable notion of crude approximation for
Eulerian Laplacians is an interesting task. However, there are Eulerian
Laplcians that do not precondtion each other regardless of stepsize, which is
a major obstruction (See Appendix D).

##  2 Overview

Existing randomized algorithms for solving linear equations in directed
Eulerian Laplacians can be classified as belonging to one of two general
frameworks: squaring-solvers and sparsified-cholesky-solvers. Both heavily
rely on $(1\pm\epsilon)$-approximate graph sparsification techniques for two
separate tasks: 1) _globally sparsifying_ a directed graph
$\underrightarrow{G}$ with $m$-edges in time close to $m$ and 2) _sparsifying
the square graph_ $\underrightarrow{G}^{2}$ in time close to $m$. The square
graph $\underrightarrow{G}^{2}$ is the graph with adjacency matrix
$\boldsymbol{\mathit{A}}_{\underrightarrow{G}}\boldsymbol{\mathit{D}}^{-1}_{\underrightarrow{G}}\boldsymbol{\mathit{A}}_{\underrightarrow{G}}$
where $\boldsymbol{\mathit{A}}_{\underrightarrow{G}}$ denotes the adjacency
matrix of $\underrightarrow{G}$ and
$\boldsymbol{\mathit{D}}_{\underrightarrow{G}}$ denotes the diagonal matrix
containing the degrees of the Eulerian graph $\underrightarrow{G}$. For some
graphs the running time of 2) is sublinear in the number of edges of the
sparsified graph $\underrightarrow{G}^{2}$ since squaring can drastically
increase the density. We provide derandomized algorithms for both tasks,
allowing us to derandomize two prominent frameworks. We focus on the squaring-
solver [CKP+17], and sketch the sparsified-Cholesky solver [PS22] in Appendix
C333[PS22] also relies on randomness for finding a linear sized almost
independent subset. We provide a simple procedure to derandomize this step in
Section C.3..

While we match the $(1\pm\epsilon)$-approximation and roughly the runtime of
randomized routines for 2), our algorithm for globally sparsifying directed
graphs only achieves a crude $n^{o(1)}$-approximation guarantee. This is not
surprising, since deterministic $(1\pm\epsilon)$-approximate sparsification in
almost linear time is an open problem even for undirected graphs. However, the
current notions of spectral approximation for directed graphs due to [CKP+17]
break down for approximation factor $\epsilon$ larger than $1$. To address
this problem, we directly define approximation via the suitability as a
preconditioner in the Richardson iteration. In particular, suppose we want to
solve a linear equation in $\boldsymbol{\mathit{M}}$, by preconditioning with
a matrix $\boldsymbol{\mathit{N}}$. Roughly speaking, if, in an appropriate
norm $\left\|\cdot\right\|$, we have
$\left\|\boldsymbol{\mathit{I}}_{\operatorname*{im}(\boldsymbol{\mathit{M}})}-\boldsymbol{\mathit{N}}^{+}\boldsymbol{\mathit{M}}\right\|\leq
1-\alpha$, then we can solve the linear equation in $\boldsymbol{\mathit{M}}$
to high accuracy using $\widetilde{O}(1/\alpha)$ preconditioned Richardson
iterations where we apply $\boldsymbol{\mathit{M}}$ and
$\boldsymbol{\mathit{N}}^{+}$ once per iteration. A priori, it may be
surprising that our approach relies on distinguishing the case
${1-\left\|\boldsymbol{\mathit{I}}_{\operatorname*{im}(\boldsymbol{\mathit{M}})}-\boldsymbol{\mathit{N}}^{+}\boldsymbol{\mathit{M}}\right\|\geq
n^{o(-1)}}$ from the case
${1-\left\|\boldsymbol{\mathit{I}}_{\operatorname*{im}(\boldsymbol{\mathit{M}})}-\boldsymbol{\mathit{N}}^{+}\boldsymbol{\mathit{M}}\right\|\approx
0}$, however, our partial symmetrization technique and a careful choice of
norm makes this possible.

Our notion of high-error approximation is not directly composable, but it does
algorithmically compose: If $\boldsymbol{\mathit{A}}$ is precondioned by
$\boldsymbol{\mathit{B}}$ and $\boldsymbol{\mathit{B}}$ is preconditioned by
$\boldsymbol{\mathit{C}}$, then a solver for $\boldsymbol{\mathit{A}}$ can be
constructed with access to $\boldsymbol{\mathit{B}}$ and
$\boldsymbol{\mathit{C}}^{+}$ using two layers of preconditioned Richardson.
The development of a directly composable notion of spectral approximation in
directed graphs for approximation factors larger than $1$ remains an
interesting open problem.

###  2.1 Global Sparsification

Given an Eulerian graph $\underrightarrow{G}=\underrightarrow{G}_{0}$ our
sparsification routine does not directly produce a sparse graph
$\tilde{\underrightarrow{G}}=\underrightarrow{G}_{3}$, but does so via a
detour involving two other graphs $\underrightarrow{G}_{1}$ and
$\underrightarrow{G}_{2}$:

  * •

$\underrightarrow{G}_{0}$ is the initial graph $\underrightarrow{G}$.

  * •

$\underrightarrow{G}_{1}$ is obtained by $\beta$-partially-symmetrizing
$\underrightarrow{G}_{0}$. It has a directed and a undirected part.

  * •

$\underrightarrow{G}_{2}$ is obtained from $\underrightarrow{G}_{1}$ by
sparsifying its directed part without touching the undirected part.

  * •

$\underrightarrow{G}_{3}$ is obtained by also sparsifying the undirected part.

We call such a bundle of four directed graphs a quadruple. Our construction
ensures that $\boldsymbol{\mathit{L}}^{+}_{\underrightarrow{G}_{i+1}}$ is a
suitable preconditioner for
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{i}}$ with respect to the norm
$\left\|\cdot\right\|_{\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{i+1}}}\rightarrow\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{i+1}}}}$
for $i=0,1,2$, meaning that
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{i}}^{+}$ can be (approximately)
applied by applying $\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{i+1}}^{+}$
for $N=Exp(O((\log n)^{1/10}))$ times via the preconditioned Richardson
iteration444Because not all approximations refer to the same norm, this does
unfortunately not ensure that
$\boldsymbol{\mathit{L}}^{+}_{\underrightarrow{G}_{3}}$ is a suitable
preconditioner for $\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{0}}$.. Given
an oracle for applying $\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{3}}^{+}$
we can use it $N$ times to apply
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{2}}^{+}$, and then in turn
apply $\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{1}}^{+}$ via $N$
applications of $\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{2}}^{+}$, until
we can finally apply $\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{0}}^{+}$.
This yields a procedure for applying
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{0}}^{+}$ that relies on
$N^{3}=Exp(O((\log n)^{1/10}))$ applications of
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{3}}^{+}$. Next we describe the
way the graphs $\underrightarrow{G}_{1},\underrightarrow{G}_{2}$ and
$\underrightarrow{G}_{3}$ are constructed in more detail. Our full global
sparsification results are presented in Section 4.

Figure 1: Obtaining $\underrightarrow{G}_{1}$ from $\underrightarrow{G}_{0}$
via $\beta$-partial symmetrization. Given the graph $\underrightarrow{G}_{0}$
as depicted on the left hand side, we add $\beta$ times its undirectification.
This could double the amount of edges, but the graph $\underrightarrow{G}_{1}$
has much less directed structure, which we will use to sparsify it in the next
steps.

#### From $\underrightarrow{G}_{0}$ to $\underrightarrow{G}_{1}$: Robustness
through Partial Symmetrization.

Let $\mathcal{U}(\underrightarrow{G}_{0})$ denote the undirectification of the
graph $\underrightarrow{G}_{0}$, i.e. the graph obtained by replacing each
directed edge with an undirected edge of half the weight. Then we simply
obtain
$\underrightarrow{G}_{1}=\beta\cdot\mathcal{U}(\underrightarrow{G}_{0})+\underrightarrow{G}_{0}$
where $\beta=n^{o(1)}$ is a sub-polynomial factor. Surprisingly, even though
$\underrightarrow{G}_{1}$ removes a lot of directed structure,
$\boldsymbol{\mathit{L}}^{+}_{\underrightarrow{G}_{0}}$ can be applied via
$O(\beta)$ applications of
$\boldsymbol{\mathit{L}}^{+}_{\underrightarrow{G}_{1}}$. Although
$\underrightarrow{G}_{1}$ might seem similar to
$\beta\cdot\mathcal{U}(\underrightarrow{G}_{0})$, in fact, they behave very
differently. The key technical observation is that

$\displaystyle\left\|\boldsymbol{\mathit{L}}_{\mathcal{U}(\underrightarrow{G}_{1})}^{+/2}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{0}}-\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{1}})\boldsymbol{\mathit{L}}_{\mathcal{U}(\underrightarrow{G}_{1})}^{+/2}\right\|_{2}=\left\|\boldsymbol{\mathit{L}}_{(1+\beta)\cdot\mathcal{U}(\underrightarrow{G}_{0})}^{+/2}\boldsymbol{\mathit{L}}_{\beta\cdot\mathcal{U}(\underrightarrow{G}_{0})}\boldsymbol{\mathit{L}}_{(1+\beta)\cdot\mathcal{U}(\underrightarrow{G}_{0})}^{+/2}\right\|_{2}=\frac{1}{\beta+1}$

which relies on the fact that the directed graphs _cancel out_ additively,
only leaving us with symmetric Laplacians. See Figure 1 for an illustration.

Figure 2: Obtaining $\underrightarrow{G}_{2}$ from $\underrightarrow{G}_{1}$
via patching. The dashed box contains the complete bipartite graph, which is a
good expander. We use this information to drastically sparsify the directed
part on the same vertex set $V^{\prime}$, only ensuring that the degrees match
up by increasing the weight of the edges. Since the induced subgraph on
$V^{\prime}$ remains a good expander when summing up, this does not alter the
spectral structure of the graph by much.

#### From $\underrightarrow{G}_{1}$ to $\underrightarrow{G}_{2}$: Patching
Expander Parts.

The graph
$\underrightarrow{G}_{1}=\beta\cdot\mathcal{U}(\underrightarrow{G}_{0})+\underrightarrow{G}_{0}$
is made up of two parts: an undirected graph
$\beta\cdot\mathcal{U}(\underrightarrow{G}_{0})$ and a directed graph
$\underrightarrow{G}_{0}$. In this step we aim to sparsify the directed part
by constructing a sparse graph $\underrightarrow{R}$ with the same in- and
out-degrees as $\underrightarrow{G}_{0}$. Then
$\underrightarrow{G}_{2}=\beta\cdot\mathcal{U}(\underrightarrow{G}_{0})+\underrightarrow{R}$.
Our strategy for obtaining $\underrightarrow{R}$ relies on the structure of
the undirected graph $\mathcal{U}(\underrightarrow{G}_{0})$ and deterministic
expander decompositions as presented in [CGL+20a]. Our main technical
observation here is that given a vertex set $V^{\prime}$ so that
$\mathcal{U}(\underrightarrow{G}_{0})[V^{\prime}]$ is a good enough expander,
we can replace the directed graph $\underrightarrow{G}_{0}[V^{\prime}]$ in
$\underrightarrow{G}_{1}$ with _any other_ directed graph
$\underrightarrow{R}^{\prime}$, as long as
$\underrightarrow{G}_{0}[V^{\prime}]$ and $\underrightarrow{R}^{\prime}$ have
exactly the same in- and out-degrees, retaining a close approximation between
$\underrightarrow{G}_{1}-\underrightarrow{G}_{0}[V^{\prime}]+\underrightarrow{R}^{\prime}$
and $\underrightarrow{G}_{1}$. Our algorithm for constructing
$\underrightarrow{R}$ first removes the weighted structure from
$\mathcal{U}(\underrightarrow{G}_{0})$ by bucketing by edge weight, and then
layers deterministic undirected and unweighted expander decompositions. These
are partitions of undirected and unweighted graphs into expanding parts and a
remainder, and the layering is achieved by recursing on the remainder. A
single expanding part can be greedily sparsified using any sparse greedily
constructed graph $\underrightarrow{R}^{\prime}$ as introduced above. Summing
up all of these yields $\underrightarrow{R}$. Crucially, while the error sums
up over layers and buckets, the disjointness of the expander parts ensures
this is not the case within a single decomposition. See Figure 2 for an
illustration.

Figure 3: Obtaining $\underrightarrow{G}_{3}$ from $\underrightarrow{G}_{2}$
via undirected sparsification. The undirected part
$\beta\cdot\mathcal{U}(\underrightarrow{G})$ of $\underrightarrow{G}_{2}$ is
sparsified using previously known spectral sparsification routines for
undirected graphs. Since both remaining parts are sparse, their sum is a
sparse graph.

#### From $\underrightarrow{G}_{2}$ to $\underrightarrow{G}_{3}$: Scaling an
Undirected Sparsifier.

Finally we sparsify the undirected part
$\beta\cdot\mathcal{U}(\underrightarrow{G}_{0})$ obtaining a sparse undirected
graph $H$ and let $\underrightarrow{G}_{3}=H+\underrightarrow{R}$. This can be
achieved via known theorems for deterministic low-accuracy undirected graph
sparsification provided by [CGL+20b]. To obtain a good preconditioner we have
to scale $H$, but not $\underrightarrow{R}$, with the inverse of an
appropriate rate $\eta$. This relies crucially on our observation that when
the approximation error is only on the undirected part, learning rates can be
leveraged more effectively than for general Eulerian approximation. We are
left with a sum of two sparse graphs $H$ and $\underrightarrow{R}$, which is a
sparse directed graph. See Figure 3 for an illustration.

###  2.2 Sparsified Squaring

Given a directed graph $\underrightarrow{G}$ with Laplacian
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}}=\boldsymbol{\mathit{D}}_{\underrightarrow{G}}-\boldsymbol{\mathit{A}}^{T}_{\underrightarrow{G}}\in\mathbb{R}^{n\times
n}$, the Laplacian of its square $\underrightarrow{G}^{2}$ is given by

$\displaystyle\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{2}}=\boldsymbol{\mathit{D}}_{\underrightarrow{G}}-\underbrace{\boldsymbol{\mathit{A}}^{T}_{\underrightarrow{G}}\boldsymbol{\mathit{D}}_{\underrightarrow{G}}^{-1}\boldsymbol{\mathit{A}}^{T}_{\underrightarrow{G}}}_{\text{adjacency
matrix}}=\boldsymbol{\mathit{D}}_{\underrightarrow{G}}-\sum_{i=1}^{n}\underbrace{\frac{1}{\boldsymbol{\mathit{D}}_{\underrightarrow{G}}(i,i)}(\boldsymbol{\mathit{A}}(i,:))^{T}\cdot(\boldsymbol{\mathit{A}}(:,i))^{T}}_{\boldsymbol{\mathit{A}}_{i}^{T}}.$

We consider the directed product graphs
$\boldsymbol{\mathit{L}}_{i}=\boldsymbol{\mathit{D}}_{i}-\boldsymbol{\mathit{A}}_{i}^{T}$
with adjacency matrix $\boldsymbol{\mathit{A}}_{i}$. Consider the matrix

$\displaystyle\boldsymbol{\mathit{L}}=\begin{pmatrix}\operatorname*{diag}(\boldsymbol{\mathit{A}}_{i}\boldsymbol{1})&\boldsymbol{0}\\\
\boldsymbol{0}&\operatorname*{diag}(\boldsymbol{\mathit{A}}_{i}^{T}\boldsymbol{1})\end{pmatrix}-\begin{pmatrix}\boldsymbol{0}&\boldsymbol{\mathit{A}}_{i}\\\
\boldsymbol{\mathit{A}}_{i}^{T}&\boldsymbol{0}\end{pmatrix}$

which is the Laplacian of a bipartite product graph $G$. Such graphs are
constant expanders, which allows for a simple trick. First, we sparsify the
undirected bipartite product graph $\boldsymbol{\mathit{L}}$ to high accuracy
$\epsilon$, by adapting a procedure from [KLP+15] to be degree preserving. We
call the sparse bipartite graph we obtain $\tilde{\boldsymbol{\mathit{L}}}$.
Then we obtain a sparsified version $\tilde{\boldsymbol{\mathit{A}}}_{i}^{T}$
of $\boldsymbol{\mathit{A}}_{i}^{T}$ by simply taking the bottom left block of
$\tilde{\boldsymbol{\mathit{L}}}$. We show that
$\tilde{\boldsymbol{\mathit{L}}}_{i}=\boldsymbol{\mathit{D}}_{i}-\tilde{\boldsymbol{\mathit{A}}}_{i}^{T}$
is a $\epsilon/\Phi^{2}$-approximation of $\boldsymbol{\mathit{L}}_{i}$, where
$\Phi$ is the expansion of $G$. Using the fact that bipartite product graphs
are constant expanders lets us directly translate the approximation guarantee,
up to a constant overhead in runtime. We obtain an
$(1\pm\epsilon)$-approximation
$\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{2}}=\sum_{i=1}^{n}\tilde{\boldsymbol{\mathit{L}}}_{i}$
of $\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{2}}$ with
$\operatorname{nnz}(\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{2}})\leq
O(\operatorname{nnz}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}})\epsilon^{-4})$
in almost linear time. Similar deterministic squaring techniques were
developed for small space algorithms [RV05, AKM+20], albeit the slightly
different guarantees. We believe its likely that their approach to be adapted
to our setting, but for convenience we adapt a more direct approach.

###  2.3 The Squaring Framework

Consider an Eulerian Laplacian
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}}=\boldsymbol{\mathit{D}}_{\underrightarrow{G}}-\boldsymbol{\mathit{A}}_{\underrightarrow{G}}^{T}$
where $\boldsymbol{\mathit{A}}_{\underrightarrow{G}}$ is the adjacency matrix
of an Eulerian graph $\underrightarrow{G}$ and
$\boldsymbol{\mathit{D}}_{\underrightarrow{G}}$ is the diagonal matrix
containing the degrees. Then the normalised Laplacian is given by
$\boldsymbol{\mathit{N}}=\boldsymbol{\mathit{I}}-\boldsymbol{\mathcal{A}}^{T}$
where
$\boldsymbol{\mathcal{A}}=\boldsymbol{\mathit{D}}_{\underrightarrow{G}}^{+/2}\boldsymbol{\mathit{A}}_{\underrightarrow{G}}^{T}\boldsymbol{\mathit{D}}_{\underrightarrow{G}}^{+/2}$,
having the property that $\left\|\boldsymbol{\mathcal{A}}\right\|_{2}\leq 1$.
The Neumann series expansion yields

$\displaystyle(\boldsymbol{\mathit{I}}-\boldsymbol{\mathcal{A}})^{+}\boldsymbol{\mathit{b}}=\sum_{i=0}^{\infty}\boldsymbol{\mathcal{A}}^{i}\boldsymbol{\mathit{b}}=\prod_{i=0}^{\infty}\left(\boldsymbol{\mathit{I}}+\boldsymbol{\mathcal{A}}^{2^{k}}\right)\boldsymbol{\mathit{b}}$

for $\boldsymbol{\mathit{b}}$ orthogonal to the kernel of
$(\boldsymbol{\mathit{I}}-\boldsymbol{\mathcal{A}})^{+}$. Given a
$1/\operatorname{poly}(n)$ lower bound on the smallest eigenvalue
$\lambda_{*}$ of $\boldsymbol{\mathit{I}}-\boldsymbol{\mathcal{A}}$,
truncating the product expansion after $\Theta(\log n)$ factors (and hence
squarings) yields a constant relative error. A very convenient equality in the
same spirit is given by

$\displaystyle(\boldsymbol{\mathit{I}}-\boldsymbol{\mathcal{A}})^{+}\boldsymbol{\mathit{b}}=(\boldsymbol{\mathit{I}}-\boldsymbol{\mathcal{A}}^{2})^{+}(\boldsymbol{\mathit{I}}+\boldsymbol{\mathcal{A}})\boldsymbol{\mathit{b}}$
(1)

for $\boldsymbol{\mathit{b}}$ orthogonal to the kernel of
$(\boldsymbol{\mathit{I}}-\boldsymbol{\mathcal{A}})^{+}$ which is at the
center of the squaring mechanism of [CKP+16b], the algorithm our squaring
solver resembles most. Their squaring scheme is in turn inspired by the
squaring solver for symmetric Laplacians presented in [PS14].

This leaves us with the task of solving linear equations in
$\boldsymbol{\mathit{I}}-\boldsymbol{\mathcal{A}}^{2}$, which is another
normalised Laplacian. However, this is the normalised Laplacian of the square
graph $\underrightarrow{G}^{2}$, and it can be shown that squaring drastically
improves the condition number of the problem, such that after $k=\Theta(\log
n)$ squaring steps linear equations can be solved to high accuracy quickly via
a simple iterative scheme. To avoid periodic behaviours we consider the
normalised adjacency matrix
$\boldsymbol{\mathcal{A}}_{l}^{(\alpha)}:=\alpha\boldsymbol{\mathit{I}}+(1-\alpha)\boldsymbol{\mathcal{A}}_{j}$,
which can be interpreted as adding self loops proportional to the out-degrees
of $\underrightarrow{G}$.

#### Sparsified Squaring.

Since squaring not only improves the condition number, but may also quickly
increase the density of the graph, we let
$\boldsymbol{\mathcal{A}}_{0}=\boldsymbol{\mathcal{A}}$ and iteratively obtain
$\boldsymbol{\mathcal{A}}_{j+1}$ by implicitly sparsifying
$(\boldsymbol{\mathcal{A}}_{j}^{(\alpha)})^{2}$ using our sparsified squaring
technique with accuracy parameter $\epsilon$. We can conclude from (1) that

$\displaystyle(\boldsymbol{\mathit{I}}-\boldsymbol{\mathcal{A}})^{+}\boldsymbol{\mathit{b}}\approx\underbrace{(1-\alpha)^{d-1}\left(\boldsymbol{\mathit{I}}-\boldsymbol{\mathcal{A}}_{d}\right)^{+}\left(\boldsymbol{\mathit{I}}+(\boldsymbol{\mathit{A}}_{d-1}^{(\alpha)})^{2}\right)\cdots\left(\boldsymbol{\mathit{I}}-(\boldsymbol{\mathit{A}}_{0}^{(\alpha)})^{2}\right)}_{:=\boldsymbol{\mathit{Z}}}\boldsymbol{\mathit{b}}$

as we can ensure that $\boldsymbol{\mathit{b}}$ is orthogonal to the known
kernel of $(\boldsymbol{\mathit{I}}-\boldsymbol{\mathcal{A}})^{+}$. However,
the repeated sparsification accumulates an error proportional to $\epsilon
e^{d}$, and it is imperative that it stays below $1$ such that
$\boldsymbol{\mathit{Z}}$ is an approximate pseudoinverse of
$\boldsymbol{\mathit{I}}-\boldsymbol{\mathcal{A}}$. Therefore, we have to
choose $\epsilon$ proportional to $e^{-d}$. We conclude from the previous
subsection that
$\operatorname{nnz}(\boldsymbol{\mathit{A}}_{d})=O(\operatorname{nnz}(\boldsymbol{\mathcal{A}})e^{4d^{2}})$.
Unfortunately, we cannot set $d=\Theta(\log n)$ without ending up with
potentially dense matrices. Therefore, we set $d=\Theta((\log n)^{1/3})$ and
have $e^{4d^{2}}=n^{o(1)}$.

#### Global Sparsification and Chains of Sparse Matrices.

Since $d=\Theta((\log n)^{1/3})$ squarings do not sufficiently decrease the
condition number, we globally sparsify after $d$ sparsified squarings and
repeat. Given
$\boldsymbol{\mathcal{A}}_{0}^{(0)}=\boldsymbol{\mathit{D}}_{\underrightarrow{G}}^{+/2}\boldsymbol{\mathit{A}}_{\underrightarrow{G}}\boldsymbol{\mathit{D}}_{\underrightarrow{G}}^{+/2}$
for $i=0,...,\Theta((\log n)^{2/3})$ we iteratively construct:

  * •

Given $\boldsymbol{\mathcal{A}}_{0}^{(i)}$, construct
$\boldsymbol{\mathcal{A}}_{0}^{(i)},...,\boldsymbol{\mathcal{A}}_{d}^{(i)}$ by
sparsified squaring as described in the previous paragraph.

  * •

Let $\underrightarrow{H}$ be the graph with adjacency matrix
$\boldsymbol{\mathit{D}}_{\underrightarrow{G}}^{1/2}\boldsymbol{\mathcal{A}}_{d}^{(i)}\boldsymbol{\mathit{D}}_{\underrightarrow{G}}^{1/2}$.
Globally sparsify $\underrightarrow{H}$ obtaining
$\tilde{\underrightarrow{H}}$. Then let
$\boldsymbol{\mathcal{A}}_{0}^{(i+1)}=\boldsymbol{\mathit{D}}_{\underrightarrow{G}}^{+/2}\boldsymbol{\mathit{A}}_{\tilde{\underrightarrow{H}}}\boldsymbol{\mathit{D}}_{\underrightarrow{G}}^{+/2}$.

We discuss these collections of squaring chains linked by global
sparsification in Section 6. For our algorithm to run in almost-linear time,
it is imperative that these chains are constructed once, and then our
algorithm operates recursively on them.

#### The Recursive Algorithm.

Our global sparsification routine allows us to solve linear equations in
$\boldsymbol{\mathit{I}}-\boldsymbol{\mathcal{A}}_{0}^{(i)}$ by solving
$\operatorname{Exp}(O((\log n)^{1/10}))$ linear equations in
$\boldsymbol{\mathit{I}}-\boldsymbol{\mathcal{A}}_{0}^{(i+1)}$. Since linear
equations in
$\boldsymbol{\mathit{I}}-\boldsymbol{\mathcal{A}}_{0}^{(\Theta(\log
n)^{2/3})}$ are easy to solve using standard iterative procedures, this is the
depth of our recursion. Therefore, the total amount of branches is
$\operatorname{Exp}(O((\log n)^{2/3+1/10}))=n^{o(1)}$. Since all involved
matrices contain an almost linear amount of entries, this gives an almost
linear time deterministic algorithm for solving linear equations involving
Eulerian Laplacians. That is, because preconditioned Richardson only does
matrix vector multiplications with the matrix and the preconditioner.

###  2.4 The Sparsified-Cholesky Framework

Very recently [PS22] showed that the framework of [KLP+16] directly works for
Eulerian Laplacians by developing new tools for analysing the accumulation of
error. Our sparsification tools can also be used to derandomize this
algorithm. Unlike the squaring framework, which makes progress by improving
the condition number, sparsified-Cholesky frameworks operate by eliminating
rows and columns like Gaussian elimination. Such elimination steps can be
directly interpreted as deleting a vertex from the graph and adding a weighted
clique.

Some algorithms eliminate one vertex at a time [KS16, CKK+18], but [KLP+16]
and [PS22] eliminate a large set of $\Omega(n)$ vertices together. To do so, a
linear sized $\rho$-_row-column-diagonally-dominant_ ($\rho$-RCDD) subset
$V^{\prime}$ of the vertices is chosen for some constant $\rho$. A set of
vertices is $\rho$-RCDD, if for each vertex a $\rho$-fraction of the weighted
in-edges come from $V\setminus V^{\prime}$ and a $\rho$-fraction of the
weighted out-edges go to $V\setminus V^{\prime}$555Notice that [KLP+15] is
concerned with the undirected case.. Then the vertices belonging to this set
can be eliminated using $O(\log\log n)$ sparsified squaring operations. While
randomized algorithms using this paradigm can afford to globally sparsify
after each squaring step, we have to allow for some build up of the edge
count. Namely, we wait for a subpolynomial number of elimination rounds, and
then globally sparsify and recurse. As previously, global sparsifications
correspond to branching points when applying the inverse. We give a more
detailed description in Appendix C.

###  2.5 Reduction to the Eulerian Setting with bounded Condition Number

Previous work [CKP+16c, CKP+16b] reduced solving linear equations in directed
Laplacians
$\boldsymbol{\mathit{L}}=\boldsymbol{\mathit{D}}-\boldsymbol{\mathit{A}}^{T}$
to solving $\log^{O(1)}(n\kappa^{-1}\epsilon^{-1})$ systems involving Eulerian
Laplacians with polynomially bounded condition number and edge weights to
constant accuracy, where $\kappa$ is an upper bound on the maximum of
$\kappa(\boldsymbol{\mathit{D}})$ and $\kappa(\boldsymbol{\mathit{L}})$ (See
Appendix D and F of [CKP+16b] and Sections 5, 7.1 and 7.3 of [CKP+16c]). They
use that edge weights are polynomially bounded in the proof of Lemma C.3 in
Appendix C of [CKP+16b]. Different reductions to the Eulerian case were
presented by Ahmadinejad-Jambulapati-Saberi-Sidford [AJSS19] and in the thesis
of Peebles [Pee19].

##  3 Preliminaries

###  3.1 Linear Algebra

#### Matrices.

We denote matrices as bold upper case letters $\boldsymbol{\mathit{A}}$. For a
matrix $\boldsymbol{\mathit{A}}\in\mathbb{R}^{n\times n}$, we let
$\operatorname{nnz}(\boldsymbol{\mathit{A}})$ denote its number of non-zero
entries and for $X\subseteq[n]$, $Y\subseteq[n]$ we let
$\boldsymbol{\mathit{A}}(X,Y)=\boldsymbol{\mathit{A}}_{XY}$ denote the
$|X|\times|Y|$ submatrix containing the entries with index in $X\times Y$. If
$X=Y$, we write $\boldsymbol{\mathit{A}}[X]$ as a shorthand for
$\boldsymbol{\mathit{A}}(X,X)$. When selecting submatrices, we let $l:u$
denote the set $\\{l,l+1,\dots,u\\}$ for $u\geq l$ and $:$ the set of all
columns/rows, e.g. $\boldsymbol{\mathit{A}}(1:3,:)$ denotes the submatrix of
$\boldsymbol{\mathit{A}}$ consisting of the first $3$ rows of
$\boldsymbol{\mathit{A}}$. Further, we let $\boldsymbol{\mathit{A}}^{+}$
denote the Moore-Penrose-Pseudoinverse of matrix $\boldsymbol{\mathit{A}}$.
Finally, $\boldsymbol{\mathit{I}}$ denotes the identity matrix. Sometimes we
denote its dimension with $\boldsymbol{\mathit{I}}_{n}$.

#### Vectors.

We denote vectors as bold lower case letters $\boldsymbol{\mathit{v}}$.
Further, we let $\boldsymbol{1}$ denote the all ones vector, $\boldsymbol{0}$
the zero vector and $\boldsymbol{\mathit{e}}_{i}$ denotes the $i$-th vector of
the standard basis. Sometimes we indicate the dimension with
$\boldsymbol{1}_{n}$ and $\boldsymbol{0}_{n}$.

#### The Loewner-order.

For symmetric matrices $\boldsymbol{\mathit{A}}$ and
$\boldsymbol{\mathit{B}}$, we let
$\boldsymbol{\mathit{A}}\preceq\boldsymbol{\mathit{B}}$ iff for all vectors
$\boldsymbol{\mathit{x}}$:
$\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}\leq\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{B}}\boldsymbol{\mathit{x}}$.
We define $\prec,\succeq$ and $\succ$ analogously. If a symmetric matrix
$\boldsymbol{\mathit{A}}$ satisfies
$\boldsymbol{0}\preceq\boldsymbol{\mathit{A}}$ we call it _PSD_. For a _PSD_
matrix $\boldsymbol{\mathit{A}}$, we let $\boldsymbol{\mathit{A}}^{1/2}$
denote the unique matrix so that
$\boldsymbol{\mathit{A}}^{1/2}\boldsymbol{\mathit{A}}^{1/2}=\boldsymbol{\mathit{A}}$.

#### Norms.

For every vector $\boldsymbol{\mathit{x}}$, we let
$\left\|\boldsymbol{\mathit{x}}\right\|_{\boldsymbol{{H}}}:=\sqrt{\boldsymbol{\mathit{x}}^{T}\boldsymbol{{H}}\boldsymbol{\mathit{x}}}$
for a _PSD_ matrix $\boldsymbol{{H}}$. We let
$\left\|\boldsymbol{\mathit{M}}\right\|_{\boldsymbol{{H}}\rightarrow\boldsymbol{{H}}}:=\max_{\boldsymbol{\mathit{x}}\neq
0}\frac{\left\|\boldsymbol{\mathit{M}}\boldsymbol{\mathit{x}}\right\|_{\boldsymbol{{H}}}}{\left\|\boldsymbol{\mathit{x}}\right\|_{\boldsymbol{{H}}}}$
for a PSD matrix $\boldsymbol{{H}}$. Notice that
$\left\|\boldsymbol{\mathit{M}}\right\|_{\boldsymbol{{H}}\rightarrow\boldsymbol{{H}}}=\left\|\boldsymbol{{H}}^{1/2}\boldsymbol{\mathit{M}}\boldsymbol{{H}}^{+/2}\right\|_{2}$.
Further, we let $\left\|\boldsymbol{\mathit{M}}\right\|_{1}$ and
$\left\|\boldsymbol{\mathit{M}}\right\|_{\infty}$ denote the maximum
$\ell_{1}$ norm of a column and row of $\boldsymbol{\mathit{M}}$ respectively.

#### Condition Number.

For a matrix $\boldsymbol{\mathit{A}}\in\mathbb{R}^{n\times n}$ we let
$\kappa(\boldsymbol{\mathit{A}}):=\left\|\boldsymbol{\mathit{A}}\right\|_{2}\left\|\boldsymbol{\mathit{A}}^{+}\right\|_{2}$
denote its condition number. Further, for PSD matrices
$\boldsymbol{\mathit{A}}$ and $\boldsymbol{\mathit{B}}$ with the same kernel
we let
$\kappa(\boldsymbol{\mathit{A}},\boldsymbol{\mathit{B}}):=\kappa(\boldsymbol{\mathit{A}}^{+/2}\boldsymbol{\mathit{B}}\boldsymbol{\mathit{A}}^{+/2})$.

#### Misc.

We let $\operatorname{Exp}(x):=e^{x}$. In this paper $\tilde{O}(\cdot)$
suppresses poly-logarithmic factors in $n$. For a PSD matrix
$\boldsymbol{\mathit{A}}$, we let $\lambda_{*}(\boldsymbol{\mathit{A}})$
denote its smallest nonzero eigenvalue.

###  3.2 Graphs

#### General Notation.

We let $G=(V,E,\omega)$ denote an undirected graph where
$\omega(e)=\omega(u,v)$ denotes the weight of edge $e=(u,v)$. Further, we let
$\underrightarrow{G}=(V,E,\omega)$, where the edge weight
$\omega(e)=\omega(u,v)$ of edge $e$ now depends on its direction. Sometimes we
omit $\omega$ for unit weight (aka unweighted) graphs. When sometimes also
write $V(\underrightarrow{G}),E(\underrightarrow{G})$ and
$\omega_{\underrightarrow{G}}$ to avoid ambiguity. We let $\omega^{\max}$ and
$\omega^{\min}$ denote the maximum and minimum edge weight respectively.

#### Undirectification.

For a directed graph $\underrightarrow{G}=(V,E,\omega)$, we let
$\mathcal{U}(\underrightarrow{G})=(V,E^{\prime},\omega)$, with $\\{u,v\\}\in
E^{\prime}$ iff $(u,v)\in E$ or $(v,u)\in E$ and
$\omega{\\{u,v\\}}=\frac{1}{2}(\omega(u,v)+\omega(u,v))$, denote its
undirectification (where we use the convention $\omega(u,v)=0$ for
$(u,v)\notin E$).

#### Induced Subgraphs.

For $G=(V,E,\omega)$ and $X\subseteq V$ we let $G[X]$ denote the induced
subgraph on $X$. For a directed graph $\underrightarrow{G}$ we define
$\underrightarrow{G}[X]$ analogously.

###  3.3 Graph Laplacians

#### General Notation.

For a undirected graph $G=(V,E,\omega)$ we denote its (graph) Laplacian as
$\boldsymbol{\mathit{L}}_{G}=\boldsymbol{\mathit{D}}_{G}-\boldsymbol{\mathit{A}}_{G}$
where $\boldsymbol{\mathit{D}}_{G}$ is the diagonal matrices containing the
degrees and
$\boldsymbol{\mathit{A}}_{G}(i,j)=\boldsymbol{\mathit{A}}_{G}(j,i)=\omega((i,j))$.
Generalizing this notation to directed graphs
$\underrightarrow{G}=(V,E,\omega)$, we let
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}}=\boldsymbol{\mathit{D}}_{\underrightarrow{G}}-\boldsymbol{\mathit{A}}^{T}_{\underrightarrow{G}}$
where the adjacency matrix is given by
$\boldsymbol{\mathit{A}}_{\underrightarrow{G}}(i,j)=\omega((i,j))$ and the
diagonal matrix
$\boldsymbol{\mathit{D}}_{\underrightarrow{G}}(i,i)=\sum_{j}\boldsymbol{\mathit{A}}_{\underrightarrow{G}}(i,j)$
contains the out degrees. Naturally
$\boldsymbol{1}^{T}\boldsymbol{\mathit{L}}_{\underrightarrow{G}}=\boldsymbol{0}$.
For an undirected graph $G$, we let $\deg_{G}(v)$ be the (weighted) degree of
vertex $v$. For directed graphs we let $\deg^{-}_{G}(v)$ and $\deg^{+}_{G}(v)$
denote in- and out-degree respectively.

#### Eulerian Laplacians.

If
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}}\boldsymbol{1}=\boldsymbol{0}$
we call a Laplacian Eulerian. This correspond to the underlying directed graph
$\underrightarrow{G}$ being Eulerian, i.e. having equal in- and out-degree for
each vertex.

#### Symmetrization.

For any matrix $\boldsymbol{\mathit{A}}$ we denote
$\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{A}}}=\boldsymbol{\mathit{U}}(\boldsymbol{\mathit{A}}):=\frac{1}{2}(\boldsymbol{\mathit{A}}+\boldsymbol{\mathit{A}}^{T})$.
Notice that for an Eulerian Laplacian
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}}$ we have
$\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}=\boldsymbol{\mathit{L}}_{\mathcal{U}(\underrightarrow{G})}$.
This is a crucial fact exploited by all algorithms for directed Laplacians
including ours. Symmetric Laplacians are _PSD_.

#### Induced subgraphs and submatrices.

The reader should note that for a Laplacian
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}}$ the matrices
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}[X]}$ and
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}}[X]$ are not equivalent unless
there is no edge from $X$ to $V\setminus X$ or vice versa. Specifically, while
the off-diagonal entries are equal, we have
$\boldsymbol{\mathit{D}}_{\underrightarrow{G}[X]}(i,i)\leq\boldsymbol{\mathit{D}}_{\underrightarrow{G}}[X](i,i)$.

###  3.4 Directed Graph Approximation

We will use the notions of approximation for directed graphs introduced in
[CKP+16b].

######  Definition 3.1 (Asymmetric Matrix Approximation, Definition 3.1 in
[CKP+16b]).

A (possibly asymmetric) matrix $\boldsymbol{\mathit{\tilde{A}}}$ is said to be
an $\epsilon$-matrix-approximation of $\boldsymbol{\mathit{A}}$ if

  1. 1.

$\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{A}}}$ is a symmetric PSD matrix,
with
$\ker(\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{A}}})\subseteq\ker(\boldsymbol{\mathit{\tilde{A}}}-\boldsymbol{\mathit{A}})\cap\ker((\boldsymbol{\mathit{\tilde{A}}}-\boldsymbol{\mathit{A}})^{T})$.

  2. 2.

$\left\|\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{A}}}^{+/2}(\boldsymbol{\mathit{\tilde{A}}}-\boldsymbol{\mathit{A}})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{A}}}^{+/2}\right\|_{2}\leq\epsilon$

######  Remark 3.2.

Notice that Definition 3.1 is not symmetric.

###  3.5 Expanders

Expander graphs, or expanders for short, play a central role in both the
workings and analysis of our algorithms. Notice that while our algorithm
operates on directed graphs, we only use expanders in the context of
undirected graphs. Whenever we use the term expander, we mean expanders with
respect to _conductance_. We follow the notational conventions of [CGL+20b].

######  Definition 3.3 (Conductance).

For a weighted but undirected graph $G=(V,E,\omega)$, given a set
$\emptyset\subset S\subset V$ we let $\delta_{G}(S):=\sum_{(u,v)\in E:u\in
S,v\notin S}\omega(u,v)$ and $\operatorname{vol}_{G}(S):=\sum_{v\in
S}\sum_{u\in V}\omega(u,v)$. Then we define the conductance

$\displaystyle\Phi_{G}(S)=\frac{\delta_{G}(S)}{\min\\{\operatorname{vol}_{G}(S),\operatorname{vol}_{G}(V\setminus
S)\\}}.$

######  Definition 3.4 (Expander).

We call a graph $G=(V,E,\omega)$ a $\Phi$-expander (with respect to
conductance) if $\min_{S:\emptyset\subset S\subset V}\Phi_{G}(S)\geq\Phi$. We
further let $\Phi_{G}=\min_{S:\emptyset\subset S\subset V}\Phi_{G}(S)$.

Given a vector $\boldsymbol{\mathit{d}}\in\mathbb{R}^{n}_{>0}$ we let
$G(\boldsymbol{\mathit{d}})$ denote the weighted and undirected graph on $n$
vertices with
$\omega(i,j)=\frac{\boldsymbol{\mathit{d}}(i)\cdot\boldsymbol{\mathit{d}}(j)}{\left\|\boldsymbol{\mathit{d}}\right\|_{1}}$.
Note that
$\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{d}})}=\operatorname*{diag}(\boldsymbol{\mathit{d}})-\frac{\boldsymbol{\mathit{d}}\boldsymbol{\mathit{d}}^{T}}{\left\|d\right\|_{1}}$.

######  Fact 3.5 (Observation 6.6 in [CGL+20b]).

$G(\boldsymbol{\mathit{d}})$ is a $1/2$-expander.

We will frequently use that Laplacians
$\boldsymbol{\mathit{L}}_{H}=\boldsymbol{\mathit{D}}_{H}-\boldsymbol{\mathit{A}}_{H}^{T}$
of good expanders $H$ are well approximated by
$\boldsymbol{\mathit{L}}_{G(\operatorname*{diag}(\boldsymbol{\mathit{D}}_{H}))}$.
To establish this, we need the following lemma.

######  Lemma 3.6 (Lemma 6.7 in [CGL+20b]).

Let $G$ and $H$ be two undirected weighted n-vertex graphs on the same vertex
set, such that $\boldsymbol{\mathit{D}}_{G}=\boldsymbol{\mathit{D}}_{H}$.
Assume further that $\Phi(G),\Phi(H)\geq\Phi$ for some threshold $\Phi$. Then
for any real vector $\boldsymbol{\mathit{x}}\in\mathbb{R}^{n}$:
$\frac{\Phi^{2}}{4}\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{L}}_{G}\boldsymbol{\mathit{x}}\leq\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{L}}_{H}\boldsymbol{\mathit{x}}^{T}\leq\frac{4}{\Phi^{2}}\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{L}}_{G}\boldsymbol{\mathit{x}}$.

We state a direct corollary in a form that turns out to be convenient for our
arguments.

######  Corollary 3.7 (Non Uniform Degree Bound).

Given a $\Phi$-expander $G$ with Laplacian
$\boldsymbol{\mathit{L}}_{G}=\boldsymbol{\mathit{D}}_{G}-\boldsymbol{\mathit{A}}^{T}_{G}$
we have

$\displaystyle\frac{\Phi^{2}}{4}(\boldsymbol{\mathit{D}}_{G}-\frac{\boldsymbol{\mathit{d}}_{G}\boldsymbol{\mathit{d}}_{G}^{T}}{\left\|\boldsymbol{\mathit{d}}_{G}\right\|_{1}})\preceq\boldsymbol{\mathit{L}}_{G}\preceq\frac{4}{\Phi^{2}}(\boldsymbol{\mathit{D}}_{G}-\frac{\boldsymbol{\mathit{d}}_{G}\boldsymbol{\mathit{d}}_{G}^{T}}{\left\|\boldsymbol{\mathit{d}}_{G}\right\|_{1}})$

for
$\boldsymbol{\mathit{d}}_{G}=\operatorname*{diag}(\boldsymbol{\mathit{D}}_{G})$
and $\Phi\leq 1/2$.

###### Proof.

Directly follows from 3.5 and Lemma 3.6. ∎

Next we state the definition of the expander decomposition for unweighted and
undirected graphs in the notation of [CGL+20b].

######  Definition 3.8 (See Section 6 of [CGL+20b], Proposed by [KVV04,
GR99]).

A $(\epsilon,\Phi)$-expander decomposition of a undirected and unweighted
graph $G=(V,E)$ is a partition $\mathcal{P}=\\{V_{1},...,V_{k}\\}$ of the
vertex set $V$ such that for all $i\in[k]$ the conductance of $G[V_{i}]$ is at
least $\Phi$ and
$\sum_{i=1}^{k}\delta_{G}(V_{i})\leq\epsilon\operatorname{vol}_{G}(V)$.

The next theorem shows that expander decompositions can be computed in almost
linear time.

######  Theorem 3.9 (Corollary 7.7 in [CGL+20b]).

There is a deterministic algorithm that, given an undirected and unweighted
graph $G=(V,E)$ with $m$ edges and parameters $\epsilon\in(0,1]$ and $1\leq
r\leq O(\log m)$, computes a $(\epsilon,\Phi)$-expander decomposition of $G$
with $\Phi\geq\Omega(\epsilon/(\log m)^{O(r^{2})})$ in time
$O\left(m^{1+O(1/r)+o(1)}\cdot(\log m)^{O(r^{2})}\right)$.

Finally, we adapt the previous theorem to interface more conveniently with our
statements.

######  Corollary 3.10 (Expander Decomposition).

There is a deterministic algorithm ExpDecomp($G,\gamma$) that, given an
undirected and unweighted graph $G=(V,E)$ with $m$ edges and a constant
$\gamma\in(0,1)$, computes a $(1/2,\frac{1}{\operatorname{Exp}((\log
n)^{\gamma})})$-expander decomposition in time $m^{1+o(1)}$.

###### Proof.

Follows directly from Theorem 3.9. ∎

###  3.6 Preconditioned Richardson

The next lemma analyses preconditioned Richardson (Algorithm 1) for asymmetric
matrices.

######  Lemma 3.11 (Preconditioned Richardson, Lemma 4.2 in [CKP+16b]).

Let $\boldsymbol{\mathit{b}}\in\mathbb{R}^{n}$ and
$\boldsymbol{\mathit{M}},\boldsymbol{\mathit{Z}},\boldsymbol{\mathit{U}}\in\mathbb{R}^{n\times
n}$ such that $\boldsymbol{\mathit{U}}$ is symmetric positive definite,
$\ker(\boldsymbol{\mathit{U}})\subseteq\ker(\boldsymbol{\mathit{M}})=\ker(\boldsymbol{\mathit{M}}^{T})=\ker(\boldsymbol{\mathit{Z}})=\ker(\boldsymbol{\mathit{Z}}^{T})$,
and $\boldsymbol{\mathit{b}}\in\operatorname*{im}(\boldsymbol{\mathit{M}})$.
Then $N$ iterations of preconditioned Richardson with step size $\eta>0$,
result in a vector
$\boldsymbol{\mathit{x}}_{N}=\textsc{PreconRichardson}(\boldsymbol{\mathit{M}},\boldsymbol{\mathit{Z}},\boldsymbol{\mathit{b}},\eta,N)$
so that

$\displaystyle\left\|\boldsymbol{\mathit{x}}_{N}-\boldsymbol{\mathit{M}}^{+}\boldsymbol{\mathit{b}}\right\|_{\boldsymbol{\mathit{U}}}\leq\left\|\boldsymbol{\mathit{I}}_{\operatorname*{im}(\boldsymbol{\mathit{M}})}-\eta\boldsymbol{\mathit{Z}}\boldsymbol{\mathit{M}}\right\|^{N}_{\boldsymbol{\mathit{U}}\rightarrow\boldsymbol{\mathit{U}}}\left\|\boldsymbol{\mathit{M}}^{+}\boldsymbol{\mathit{b}}\right\|_{\boldsymbol{\mathit{U}}}.$

Furthermore preconditioned Richardson implements a linear operator, in the
sense that
$\boldsymbol{\mathit{x}}_{N}=\boldsymbol{\mathit{Z}}_{N}\boldsymbol{\mathit{b}}$
for some matrix $\boldsymbol{\mathit{Z}}_{N}$ only depending on
$\boldsymbol{\mathit{Z}},\boldsymbol{\mathit{M}},\eta$ and $N$.

The previous lemma leads to the notion of an approximate pseudoinverse by
measuring the suitability of a matrix as a preconditioner.

######  Definition 3.12 (Approximate Pseudoinverse, Definition 4.3 in
[CKP+16b]).

Matrix $\boldsymbol{\mathit{Z}}$ is an $\epsilon$-approximate-pseudoinverse of
matrix $\boldsymbol{\mathit{M}}$ with respect to a _PSD_ matrix
$\boldsymbol{\mathit{U}}$, if
$\ker(\boldsymbol{\mathit{U}})\subseteq\ker(\boldsymbol{\mathit{M}})=\ker(\boldsymbol{\mathit{M}})^{T}=\ker(\boldsymbol{\mathit{Z}})=\ker(\boldsymbol{\mathit{Z}}^{T})$,
and

$\displaystyle\left\|\boldsymbol{\mathit{I}}_{\operatorname*{im}(\boldsymbol{\mathit{M}})}-\boldsymbol{\mathit{Z}}\boldsymbol{\mathit{M}}\right\|_{\boldsymbol{\mathit{U}}\rightarrow\boldsymbol{\mathit{U}}}\leq\epsilon.$

1 $\boldsymbol{\mathit{x}}_{0}=\boldsymbol{0}$

2 for  _i = 0, …, N - 1_ do

3
$\boldsymbol{\mathit{x}}_{i+1}=\boldsymbol{\mathit{x}}_{i}+\eta\boldsymbol{\mathit{Z}}(\boldsymbol{\mathit{b}}-\boldsymbol{\mathit{M}}\boldsymbol{\mathit{x}}_{i})$

4

5 end for

return _$\boldsymbol{\mathit{x}}_{N}$_

Algorithm 1
PreconRichardson($\boldsymbol{\mathit{M}},\boldsymbol{\mathit{Z}},\boldsymbol{\mathit{b}},\eta,N$)

Finally, we state three more lemmas that are crucial for arguing about
approximate pseudoinverses. The first two are often applied consecutively to
upper bound
$\left\|\boldsymbol{\mathit{I}}_{\operatorname*{im}(\boldsymbol{\mathit{M}})}-\boldsymbol{\mathit{Z}}\boldsymbol{\mathit{M}}\right\|_{\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{Z}}}\rightarrow\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{Z}}}}$
with a variational form.

######  Lemma 3.13 (Part of Lemma B.9 in [CKP+16b]).

If $\boldsymbol{\mathit{L}}$ is a matrix with
$\ker(\boldsymbol{\mathit{L}})=\ker(\boldsymbol{\mathit{L}}^{T})=\ker(\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}})$,
and $\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}}$ is positive
semidefinite, then for any matrix $\boldsymbol{\mathit{A}}$ with the same left
and right kernels as $\boldsymbol{\mathit{L}}$ we have

$\displaystyle\left\|\boldsymbol{\mathit{A}}\right\|_{\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}}\rightarrow\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}}}\leq\left\|\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}}^{+/2}\boldsymbol{\mathit{L}}\boldsymbol{\mathit{A}}\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}}^{+/2}\right\|_{2}$

######  Lemma 3.14 (Part of Lemma B.2 in [CKP+16b]).

For all $\boldsymbol{\mathit{A}}\in\mathbb{R}^{n\times n}$ and symmetric _PSD_
$\boldsymbol{\mathit{M}},\boldsymbol{\mathit{N}}\in\mathbb{R}^{n\times n}$
such that
$\ker(\boldsymbol{\mathit{M}})\subseteq\ker(\boldsymbol{\mathit{A}}^{T})$ and
$\ker(\boldsymbol{\mathit{N}})\subseteq\ker(\boldsymbol{\mathit{A}})$ we have

$\displaystyle\left\|\boldsymbol{\mathit{M}}^{+/2}\boldsymbol{\mathit{A}}\boldsymbol{\mathit{N}}^{+/2}\right\|_{2}=2\max_{\boldsymbol{\mathit{x}},\boldsymbol{\mathit{y}}\neq
0}\frac{\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{A}}\boldsymbol{\mathit{y}}}{\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{M}}\boldsymbol{\mathit{x}}+\boldsymbol{\mathit{y}}^{T}\boldsymbol{\mathit{N}}\boldsymbol{\mathit{y}}}$

where we define $0/0$ to be $0$.

######  Lemma 3.15 (Part of Lemma B.4 in [CKP+16b]).

For a _PSD_ diagonal matrix $\boldsymbol{\mathit{D}}$ and any matrix
$\boldsymbol{\mathit{M}}\in\mathbb{R}^{n\times n}$

$\displaystyle\left\|\boldsymbol{\mathit{D}}^{-1/2}\boldsymbol{\mathit{M}}\boldsymbol{\mathit{D}}^{-1/2}\right\|_{2}\leq\max\\{\left\|\boldsymbol{\mathit{D}}^{-1}\boldsymbol{\mathit{M}}\right\|_{\infty},\left\|\boldsymbol{\mathit{D}}^{-1}\boldsymbol{\mathit{M}}^{T}\right\|_{\infty}\\}=\max\\{\left\|\boldsymbol{\mathit{M}}^{T}\boldsymbol{\mathit{D}}^{-1}\right\|_{1},\left\|\boldsymbol{\mathit{M}}\boldsymbol{\mathit{D}}^{-1}\right\|_{1}\\}.$

##  4 Global Sparsification for Directed Laplacians

In this section we describe our low accuracy global sparsification routine.
This constitutes the backbone of our algorithm. We first formally define the
$\beta$-partial-symmetrization of an Eulerian graph $\underrightarrow{G}$. See
Figure 4 for an illustration.

######  Definition 4.1.

For an Eulerian directed graph $\underrightarrow{G}$ we call
$\mathcal{U}^{(\beta)}(\underrightarrow{G}):=\beta\cdot\mathcal{U}(\underrightarrow{G})+\underrightarrow{G}$
its $\beta$-partial-symmetrization.

######  Remark 4.2.

$\boldsymbol{\mathit{L}}_{\mathcal{U}^{(\beta)}(\underrightarrow{G})}=\beta\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}+\boldsymbol{\mathit{L}}_{\underrightarrow{G}}$.

There are three steps in our sparsification procedure.

  1. 1.

The first step relies on what may be the most crucial observation. Given an
Eulerian directed graph $\underrightarrow{G}=\underrightarrow{G}_{0}$, we let
$\underrightarrow{G}_{1}=\mathcal{U}^{(\beta)}(\underrightarrow{G})=\beta\cdot\mathcal{U}(\underrightarrow{G})+\underrightarrow{G}$
be the graph obtained from $\underrightarrow{G}_{0}$ by $\beta$-partial
symmetrization. Then, surprisingly,
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{1}}$ can be used as a
preconditioner for solving linear equations in
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{0}}$ in time $O(\beta)$ using
Richardson. We call $\mathcal{U}^{(\beta)}(\underrightarrow{G})$ the
$\beta$-partial-symmetrization of $\underrightarrow{G}$.

  2. 2.

A partial-symmetrization is naturally interpreted as the sum of an undirected
graph $\beta\cdot\mathcal{U}(\underrightarrow{G})$ and a directed graph
$\underrightarrow{G}$. We expander decompose the undirected graph
$\beta\cdot\mathcal{U}(\underrightarrow{G})$ into parts
$V_{1},V_{2},...,V_{k}$. Then a simple greedy patching scheme can be used to
sparsify the induced subgraphs $\underrightarrow{G}[V_{i}]$ leveraging the
expander structure for error control. In our actual algorithm we additionally
bucket by edge weight and layer expander decompositions. The former allows us
to treat the graph as unweighted and the latter ensures that every edge is in
an expander after $O(\log n)$ layers. The Laplacian of the obtained graph
$\underrightarrow{G}_{2}=\beta\cdot\mathcal{U}(\underrightarrow{G})+\underrightarrow{R}$
can then be used as a preconditioner for
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{1}}$

  3. 3.

Lastly, the undirected graph $\beta\cdot\mathcal{U}(\underrightarrow{G})$ can
be sparsified via previously known deterministic algorithms presented in
[CGL+20b]. We obtain
$\underrightarrow{G}_{3}=\frac{\beta}{\eta}\tilde{G}+\underrightarrow{R}$
which in turn is a preconditioner for $\underrightarrow{G}_{2}$.

Figure 4: This drawing illustrates the concept of an $\beta$-partial-
symmetrization. We use boldness to roughly indicate edge weights. On the left
hand side an Eulerian graph $\underrightarrow{G}$ is depicted. On the right
hand side its $\beta$-partial-symmetrization is drawn, alongside its natural
split into an undirected and directed part.

We define pseudoinverse sparsification quadruples. Constructing these is at
the core of our global sparsification routine.

######  Definition 4.3.

(Sparsification quadruple) We call strongly connected $n$-vertex Eulerian
graphs
$\underrightarrow{G}_{0},\underrightarrow{G}_{1},\underrightarrow{G}_{2},\underrightarrow{G}_{3}$
a $(\gamma,\beta,\eta)$-quadruple for some constant $\gamma\in(0,1)$ if

  1. 1.

$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{i}}^{+}$ is a
$\left(1-\frac{1}{\operatorname{Exp}(O((\log
n)^{\gamma}))}\right)$-approximate pseudoinverse of
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{i-1}}$ with respect to
$\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{i}}}$
for $i=1,2,3$.

  2. 2.

$\frac{1}{\operatorname{Exp}(O((\log
n)^{\gamma}))}\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{i-1}}}\preceq\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{i}}}\preceq\operatorname{Exp}(O((\log
n)^{\gamma}))\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{i-1}}}$
for $i\in\\{1,2,3\\}$.

  3. 3.

$|E(\underrightarrow{G}_{3})|=\tilde{O}(n)$. $|E(\underrightarrow{G}_{i})|\leq
2|E(\underrightarrow{G}_{0})|+\tilde{O}(n)$ for $i=1,2$.

  4. 4.

For all vertices $v$:
$\deg_{\underrightarrow{G}_{0}}^{+}(v)=\deg_{\underrightarrow{G}_{0}}^{-}(v)=(1+\beta)\deg_{\underrightarrow{G}_{i}}^{+}(v)=(1+\beta)\deg_{\underrightarrow{G}_{i}}^{-}(v)$
for $i=1,2$ and
$\deg_{\underrightarrow{G}_{0}}^{+}(v)=(1+\frac{\beta}{\eta})\deg_{\underrightarrow{G}_{3}}^{+}(v)=(1+\frac{\beta}{\eta})\deg_{\underrightarrow{G}_{3}}^{-}(v)$.

We then state the main lemma of this section. It shows that it is possible to
construct a ($\gamma,\beta,\eta$)-quadruple in almost linear time.

######  Lemma 4.4 (Global Sparsification).

For every $m$-edge, strongly connected Eulerian graph
$\underrightarrow{G}=\underrightarrow{G}_{0}$ and a constant $\gamma\in(0,1)$
the routine
$\underrightarrow{G}_{1},\underrightarrow{G}_{2},\underrightarrow{G}_{3}=\textsc{GlobalSparsification}(\underrightarrow{G},\gamma)$
yields a $(\gamma,\beta=\operatorname{Exp}(O((\log
n)^{\gamma})),\eta=\operatorname{Exp}(-3(\log n)^{\gamma}))$-quadruple
$\underrightarrow{G}_{0},\underrightarrow{G}_{1},\underrightarrow{G}_{2},\underrightarrow{G}_{3}$.
The runtime is $m^{1+o(1)}$.

######  Remark 4.5.

While in the description of our algorithm, $\beta$ scales linearly in
$O\left(\log\left(\frac{\omega^{max}_{\underrightarrow{G}}}{\omega^{min}_{\underrightarrow{G}}}\right)\right)$,
we assume throughout the paper that
$\beta=\tilde{O}(1)\cdot\operatorname{Exp}(2(\log n)^{\gamma})$ is fixed to a
global upper bound as
$\log\left(\frac{\omega^{max}_{\underrightarrow{G}}}{\omega^{min}_{\underrightarrow{G}}}\right)=\tilde{O}(1)$
for all graphs ${\underrightarrow{G}}$ we work with. This avoids clutter in
the analysis.

1 $\beta=L\cdot\operatorname{Exp}(2\cdot(\log n)^{\gamma})$ for $L=128\cdot
20\cdot P\cdot\log n$ and
$P=\left\lceil\log\left(\frac{\omega^{\max}_{\underrightarrow{G}}}{\omega^{\min}_{\underrightarrow{G}}}\right)\right\rceil$.

2 $\eta=\operatorname{Exp}(-3\cdot(\log n)^{\gamma})$

$\underrightarrow{G}_{1}=\beta\mathcal{U}(\underrightarrow{G})+\underrightarrow{G}$
;

// Note that
$\underrightarrow{G}_{1}=\mathcal{U}^{(\beta)}(\underrightarrow{G})$.

3 $\underrightarrow{R}=\textsc{SparsifyDirected}(\underrightarrow{G},\gamma)$

4
$\underrightarrow{G}_{2}=\beta\mathcal{U}(\underrightarrow{G})+\underrightarrow{R}$

5
$\tilde{G}=\textsc{SpectralSparsifyDeg}(\mathcal{U}(\underrightarrow{G}),\gamma)$

6 $\underrightarrow{G}_{3}=\frac{\beta}{\eta}\tilde{G}+\underrightarrow{R}$

return
_$\underrightarrow{G}_{1},\underrightarrow{G}_{2},\underrightarrow{G}_{3}$_

Algorithm 2 GlobalSparsification($\underrightarrow{G},\gamma$)

###  4.1 Preconditioning with the Partial-Symmetrization

Our next lemma shows that
$\boldsymbol{\mathit{L}}_{\mathcal{U}^{(\beta)}(\underrightarrow{G})}$ is a
good preconditioner in terms of Lemma 3.11.

######  Lemma 4.6.

For every Eulerian Laplacian
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{0}}$ the matrix
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{1}}^{+}$ is an
$(1-\frac{1}{1+\beta})$-approximate pseudoinverse of
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{0}}$ with respect to
$\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{1}}}$
if $\underrightarrow{G}_{1}=\mathcal{U}^{(\beta)}(\underrightarrow{G}_{0})$.

###### Proof.

First recall that
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{1}}^{+}=\boldsymbol{\mathit{L}}_{\mathcal{U}^{(\beta)}(\underrightarrow{G})}^{+}$.
We have

$\displaystyle\left\|\boldsymbol{\mathit{I}}_{\operatorname*{im}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}})}-\boldsymbol{\mathit{L}}_{\mathcal{U}^{(\beta)}(\underrightarrow{G})}^{+}\boldsymbol{\mathit{L}}_{\underrightarrow{G}}\right\|_{\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\mathcal{U}^{(\beta)}(\underrightarrow{G})}}\rightarrow\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\mathcal{U}^{(\beta)}(\underrightarrow{G})}}}\leq\left\|\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\mathcal{U}^{(\beta)}(\underrightarrow{G})}}^{+/2}(\boldsymbol{\mathit{L}}_{\mathcal{U}^{(\beta)}(\underrightarrow{G})}-\boldsymbol{\mathit{L}}_{\underrightarrow{G}})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\mathcal{U}^{(\beta)}(\underrightarrow{G})}}^{+/2}\right\|_{2}$

by Lemma 3.13. With Remark 4.2 we conclude

$\displaystyle\left\|\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\mathcal{U}^{(\beta)}(\underrightarrow{G})}}^{+/2}(\boldsymbol{\mathit{L}}_{\mathcal{U}^{(\beta)}(\underrightarrow{G})}-\boldsymbol{\mathit{L}}_{\underrightarrow{G}})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\mathcal{U}^{(\beta)}(\underrightarrow{G})}}^{+/2}\right\|_{2}$
$\displaystyle=\beta\left\|\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\mathcal{U}^{(\beta)}(\underrightarrow{G})}}^{+/2}\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\mathcal{U}^{(\beta)}(\underrightarrow{G})}}^{+/2}\right\|_{2}$
$\displaystyle=\frac{\beta}{1+\beta}\left\|\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}\right\|_{2}$
$\displaystyle=\frac{\beta}{1+\beta}=1-\frac{1}{1+\beta}$

where we use that
$\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\mathcal{U}^{(\beta)}(\underrightarrow{G})}}=\boldsymbol{\mathit{L}}_{\mathcal{U}(\mathcal{U}^{(\beta)}(\underrightarrow{G}))}=(1+\beta)\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}$.
The lemma follows from chaining the calculations. ∎

###  4.2 Sparsifying the Directed Part

Given
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{1}}=\beta\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}+\boldsymbol{\mathit{L}}_{\underrightarrow{G}}$,
we aim to obtain a sparse directed graph $\underrightarrow{R}$ with the same
in- and out-degrees as $\underrightarrow{G}$ so that the directed Laplacian
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{2}}=\beta\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}+\boldsymbol{\mathit{L}}_{\underrightarrow{R}}$
preconditions $\boldsymbol{\mathit{L}}_{\underrightarrow{G}}$. Our strategy
closely follows common strategies for sparsifying undirected graphs via
expander decompositions. First we get rid of most of the weighted structure by
bucketing by edge weight. We obtain $\tilde{O}(1)$ graphs
$\underrightarrow{G}^{(i)}$ with close to uniform edge weight such that
$\sum_{i}\underrightarrow{G}^{(i)}=\underrightarrow{G}$.

We let $H^{(i)}$ denote the unweighted and undirected graph with the same
edges as $\underrightarrow{G}^{(i)}$. Then we layer $j=1,...,O(\log n)$
undirected and unweighted expander decompositions on this graph, where each of
them peels of at least $1/2$ of the remaining edges $E_{r}^{(i,j)}$. This
procedure computes $O(\log n)$ partitions
$V^{(i,j)}_{1},...,V^{(i,j)}_{k(i,j)}$ of the vertex set, such that for each
component $V^{(i,j)}_{p}$ the graph $H^{(i)}[V^{(i,j)}_{p}]$ is an expander.
In the $j$-th layer, we put the remaining edges of the directed graph
$\underrightarrow{G}^{(i)}$ that do not go across sets in the partition
$V^{(i,j)}_{1},...,V^{(i,j)}_{k(i,j)}$ into the graph
$\underrightarrow{G}^{(i,j)}$ and remove them from the set of remaining edges.

The expander structure allows us to sparsify $\underrightarrow{G}^{(i,j)}$ via
a greedy patching scheme obtaining $\tilde{\underrightarrow{G}}^{(i,j)}$.
Finally, we sum up across layers and buckets and obtain
$\underrightarrow{R}=\sum_{i,j}\tilde{\underrightarrow{G}}^{(i,j)}$. We
leverage the robustness introduced by partial symmetrization to bound the
error. See Algorithm 3 for detailed pseudocode. We first state the main lemma
of this subsection, which analyses this algorithm.

######  Lemma 4.7.

Let
$\underrightarrow{R}=\textsc{SparsifyDirected}(\underrightarrow{G},\gamma)$
for $\gamma\in(0,1)$ constant. Then
$\boldsymbol{\mathit{L}}^{+}_{\underrightarrow{G}_{2}}=(\beta\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}+\boldsymbol{\mathit{L}}_{\underrightarrow{R}})^{+}$
is a $1/2$-approximate pseudoinverse of
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{1}}=\boldsymbol{\mathit{L}}_{\mathcal{U}^{\beta}(\underrightarrow{G})}$
with respect to
$\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{2}}}$
for $\beta=\tilde{O}(1)\cdot\operatorname{Exp}(2\cdot(\log n)^{\gamma})$.
Further, the graph $\underrightarrow{R}$ has $\tilde{O}(n)$ edges and the same
in- and out-degrees as $\underrightarrow{G}$.

1 Algorithm  _SparsifyDirected(_$\underrightarrow{G},\gamma$_)_

2
$P=\left\lceil\log\left(\frac{\omega^{\max}_{\underrightarrow{G}}}{\omega^{\min}_{\underrightarrow{G}}}\right)\right\rceil$

3 for  _$i=1,...,P$_ do

4 $E^{(i)}=\\{e\in
E(\underrightarrow{G}):\omega_{\underrightarrow{G}}^{min}\cdot
2^{i-1}\leq\omega_{\underrightarrow{G}}(e)<\omega_{\underrightarrow{G}}^{min}\cdot
2^{i}\\}$

5
$\tilde{\underrightarrow{G}}^{(i)}=\textsc{Sparsify}(\underrightarrow{G}^{(i)}=(V(\underrightarrow{G}),E^{(i)},\omega_{\underrightarrow{G}}))$

6

7 end for

8 return
_$\underrightarrow{R}=\sum_{i=1}^{P}\tilde{\underrightarrow{G}}^{(i)}$_

9

10

11Procedure  _Sparsify(_$\underrightarrow{G}^{(i)}$_)_

12 Let $H^{(i)}$ denote the unweighted and undirected graph with the same
edges as $\mathcal{U}(\underrightarrow{G}^{(i)})$

13 for  _$j=1,...,10\log n$_ do

14 $E_{r}^{(i,j)}=E({H}^{(i)})-\bigcup_{l=1}^{j-1}E({H}^{(i,j)})$;
$\underrightarrow{E}_{r}^{(i,j)}=E(\underrightarrow{G}^{(i)})-\bigcup_{l=1}^{j-1}E(\underrightarrow{G}^{(i,j)})$

15
$V^{(i,j)}_{1},...,V_{k(i,j)}^{(i,j)}=\textsc{ExpDecomp}((\boldsymbol{\mathit{V}}({H}^{(i)}),E_{r}^{(i,j)}),\gamma)$

16 $E({H}^{(i,j)})=\bigcup_{p=1}^{k(i,j)}\\{(u,v)\in E_{r}^{(i,j)}:u\in
V^{(i,j)}_{p}\land v\in V^{(i,j)}_{p}\\}$

17
$E(\underrightarrow{G}^{(i,j)})=\bigcup_{p=1}^{k(i,j)}\\{(u,v)\in\underrightarrow{E}_{r}^{(i,j)}:u\in
V^{(i,j)}_{p}\land v\in V^{(i,j)}_{p}\\}$

18 $H^{(i,j)}=(V({H}^{(i)}),E({H}^{(i,j)}))$;
$\underrightarrow{G}^{(i,j)}=(V(\underrightarrow{G}^{(i)}),E(\underrightarrow{G}^{(i,j)}),\omega_{\underrightarrow{G}^{(i)}})$

19
$\tilde{\underrightarrow{G}}^{(i,j)}=\sum_{p=1}^{k(i,j)}\textsc{Patch}(\underrightarrow{G}^{(i,j)}[V^{(i,j)}_{p}])$

20 end for

21 return _$\tilde{\underrightarrow{G}}^{(i)}=\sum_{j=1}^{10\log
n}\tilde{\underrightarrow{G}}^{(i,j)}$_

22

23

24Procedure  _Patch(_$\underrightarrow{H}$_)_

Let
$\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}}\in\boldsymbol{\mathit{R}}^{n}_{\geq
0}$ so that $\boldsymbol{\mathit{a}}(v)=\deg^{+}_{\underrightarrow{H}}(v)$ and
$\boldsymbol{\mathit{b}}(v)=\deg^{-}_{\underrightarrow{H}}(v)$. // Note
$\left\|\boldsymbol{\mathit{a}}\right\|_{1}=\left\|\boldsymbol{\mathit{b}}\right\|_{1}$.

25 $E(\tilde{\underrightarrow{H}})=\emptyset$;
$\omega_{\tilde{\underrightarrow{H}}}(e)=0$ for all $e$.

26 while  _$\left\|\boldsymbol{\mathit{a}}\right\|\neq 0$_ do

27 Let $i$ and $j$ be arbitrary such that ${\boldsymbol{\mathit{a}}}(i)>0$ and
${\boldsymbol{\mathit{b}}}(j)>0$.

28 $w=\min\\{\boldsymbol{\mathit{a}}(i),\boldsymbol{\mathit{b}}(j)\\}$;
$\boldsymbol{\mathit{a}}(i)=\boldsymbol{\mathit{a}}(i)-w$;
$\boldsymbol{\mathit{b}}(j)=\boldsymbol{\mathit{b}}(j)-w$

29
$E(\tilde{\underrightarrow{H}})=E(\tilde{\underrightarrow{H}})\cup\\{(i,j)\\}$;
$\omega(i,j)=w$

30 end while

31 return _$\tilde{\underrightarrow{H}}$_

32

Algorithm 3 SparsifyDirected($\underrightarrow{G},\gamma$) and subroutines
Sparsify() and Patch()

The proof of Lemma 4.7 relies on analysing the error incurred by sparsifying
each individual expander decomposition, i.e. the cost of replacing
$\underrightarrow{G}^{(i,j)}$ with $\tilde{\underrightarrow{G}}^{(i,j)}$. Then
we conclude by just summing up the error. The next lemma carefully analyses
the amount of error sparsifying such an expander decomposition creates. It
crucially relies on the fact that the expander parts form a disjoint
partition, and therefore the error does not scale in the number of expanders.

######  Lemma 4.8.

In the context of Sparsify() in Algorithm 3 we have

$\displaystyle\left\|\boldsymbol{\mathit{U}}_{L_{\underrightarrow{G}}}^{+/2}(\underrightarrow{G}^{(i,j)}-\tilde{\underrightarrow{G}}^{(i,j)})\boldsymbol{\mathit{U}}_{L_{\underrightarrow{G}}}^{+/2}\right\|_{2}\leq
128\cdot\operatorname{Exp}(2\cdot(\log n)^{\gamma})$

for every edge weight bucket $i$ and expander decomposition layer $j$.

###### Proof.

$\underrightarrow{G}^{(i,j)}$ and $\tilde{\underrightarrow{G}}^{(i,j)}$ are
directed graphs with the same in and out degrees since Patch() in Algorithm 3
preserves degrees exactly. Therefore $\boldsymbol{1}$ is in both the left and
right kernel of
$\underrightarrow{G}^{(i,j)}-\tilde{\underrightarrow{G}}^{(i,j)}$. We apply
Lemma 3.14 twice and obtain

$\displaystyle\left\|\boldsymbol{\mathit{U}}_{L_{\underrightarrow{G}}}^{+/2}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}}-\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{(i,j)}})\boldsymbol{\mathit{U}}_{L_{\underrightarrow{G}}}^{+/2}\right\|_{2}$
$\displaystyle=2\max_{\boldsymbol{\mathit{x}},\boldsymbol{\mathit{y}}\neq\boldsymbol{0}}\frac{\boldsymbol{\mathit{x}}^{T}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}}-\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{(i,j)}})\boldsymbol{\mathit{y}}}{\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{U}}_{L_{\underrightarrow{G}}}\boldsymbol{\mathit{x}}+\boldsymbol{\mathit{y}}\boldsymbol{\mathit{U}}_{L_{\underrightarrow{G}}}\boldsymbol{\mathit{y}}^{T}}$
$\displaystyle\stackrel{{\scriptstyle
i)}}{{\leq}}2\max_{\boldsymbol{\mathit{x}},\boldsymbol{\mathit{y}}\neq\boldsymbol{0}}\frac{\boldsymbol{\mathit{x}}^{T}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}}-\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{(i,j)}})\boldsymbol{\mathit{y}}}{\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{L}}_{\mathcal{U}(\underrightarrow{G}^{(i,j)})}\boldsymbol{\mathit{x}}+\boldsymbol{\mathit{y}}\boldsymbol{\mathit{L}}_{\mathcal{U}(\underrightarrow{G}^{(i,j)})}\boldsymbol{\mathit{y}}^{T}}$
$\displaystyle\stackrel{{\scriptstyle
ii)}}{{\leq}}\frac{2}{\omega^{\min}_{\underrightarrow{G}}\cdot
2^{i-2}}\max_{\boldsymbol{\mathit{x}},\boldsymbol{\mathit{y}}\neq\boldsymbol{0}}\frac{\boldsymbol{\mathit{x}}^{T}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}}-\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{(i,j)}})\boldsymbol{\mathit{y}}}{\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{L}}_{H^{(i,j)}}\boldsymbol{\mathit{x}}+\boldsymbol{\mathit{y}}\boldsymbol{\mathit{L}}_{H^{(i,j)}}\boldsymbol{\mathit{y}}^{T}}$
$\displaystyle=\frac{1}{\omega^{\min}_{\underrightarrow{G}}\cdot
2^{i-2}}\left\|\boldsymbol{\mathit{L}}_{H^{(i,j)}}^{+/2}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}}-\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{(i,j)}})\boldsymbol{\mathit{L}}_{H^{(i,j)}}^{+/2}\right\|_{2}$
(2)

where i) follows since
$\boldsymbol{\mathit{L}}_{\mathcal{U}(\underrightarrow{G}^{(i,j)})}\preceq\boldsymbol{\mathit{U}}_{L_{\underrightarrow{G}}}$
because $\mathcal{U}(\underrightarrow{G}^{(i,j)})$ is a subgraph of
$\mathcal{U}(\underrightarrow{G})$ and ii) uses that all edge weights in
$\mathcal{U}(\underrightarrow{G}^{i,j})$ are in
$[\omega^{min}_{\underrightarrow{G}}\cdot
2^{i-2},\omega^{min}_{\underrightarrow{G}}\cdot 2^{i}]$. Next we can use that
the expander parts $V^{(i,j)}_{p}$ are disjoint in both
$\underrightarrow{G}^{(i,j)}$ and $\tilde{\underrightarrow{G}}^{(i,j)}$ to
bound

$\displaystyle\left\|\boldsymbol{\mathit{L}}_{H^{(i,j)}}^{+/2}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}}-\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{(i,j)}})\boldsymbol{\mathit{L}}_{H^{(i,j)}}^{+/2}\right\|_{2}$
$\displaystyle\leq\max_{p}\left\|\boldsymbol{\mathit{L}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+/2}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}[V^{(i,j)}_{p}]}-\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{(i,j)}[V^{(i,j)}_{p}]})\boldsymbol{\mathit{L}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+/2}\right\|_{2}$
(3)

since the spectral norm of a block diagonal matrix is upper bounded by the
maximum spectral norm of a block (See A.1). Then, for every
$p\in\\{1,...,k(i,j)\\}$ we have

$\displaystyle\left\|\boldsymbol{\mathit{L}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+/2}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}[V^{(i,j)}_{p}]}-\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{(i,j)}[V^{(i,j)}_{p}]})\boldsymbol{\mathit{L}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+/2}\right\|_{2}\stackrel{{\scriptstyle
i)}}{{=}}2\max_{\boldsymbol{\mathit{x}},\boldsymbol{\mathit{y}}\neq\boldsymbol{0}}\frac{\boldsymbol{\mathit{x}}^{T}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}[V^{(i,j)}_{p}]}-\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{(i,j)}[V^{(i,j)}_{p}]})\boldsymbol{\mathit{y}}}{\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{L}}_{H^{(i,j)}[V^{(i,j)}_{p}]}\boldsymbol{\mathit{x}}+\boldsymbol{\mathit{y}}\boldsymbol{\mathit{L}}_{H^{(i,j)}[V^{(i,j)}_{p}]}\boldsymbol{\mathit{y}}^{T}}$
$\displaystyle\stackrel{{\scriptstyle ii)}}{{\leq}}8\cdot 2^{2(\log
n)^{\gamma}}\max_{\boldsymbol{\mathit{x}},\boldsymbol{\mathit{y}}\neq\boldsymbol{0}}\frac{\boldsymbol{\mathit{x}}^{T}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}[V^{(i,j)}_{p}]}-\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{(i,j)}[V^{(i,j)}_{p}]})\boldsymbol{\mathit{y}}}{\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{d}}^{(i,j)}_{p})}\boldsymbol{\mathit{x}}+\boldsymbol{\mathit{y}}\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{d}}^{(i,j)}_{p})}\boldsymbol{\mathit{y}}^{T}}$
(4)

for $\boldsymbol{\mathit{d}}^{(i,j)}_{p}$ being the degree vector of
$H^{(i,j)}[V^{(i,j)}_{p}]$, where i) is by Lemma 3.14 and ii) is by the
expansion of $H^{(i,j)}[V^{(i,j)}_{p}]$ and Corollary 3.7. Let
$\boldsymbol{\mathit{x}},\boldsymbol{\mathit{y}}\perp\boldsymbol{1}$ be
maximising the right hand side of the previous inequality. Then, also
$\boldsymbol{\mathit{x}}^{\prime}=\boldsymbol{\mathit{x}}-\frac{\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{d}}^{(i,j)}_{p}}{\left\|\boldsymbol{\mathit{d}}^{(i,j)}_{p}\right\|_{2}}\boldsymbol{1}$
and
$\boldsymbol{\mathit{y}}^{\prime}=\boldsymbol{\mathit{y}}-\frac{\boldsymbol{\mathit{y}}^{T}\boldsymbol{\mathit{d}}^{(i,j)}_{p}}{\left\|\boldsymbol{\mathit{d}}^{(i,j)}_{p}\right\|_{2}}\boldsymbol{1}$
are maximizing. We obtain

$\displaystyle
2\frac{\boldsymbol{\mathit{x}}^{T}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}[V^{(i,j)}_{p}]}-\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{(i,j)}[V^{(i,j)}_{p}]})\boldsymbol{\mathit{y}}}{\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{d}}^{(i,j)}_{p})}\boldsymbol{\mathit{x}}+\boldsymbol{\mathit{y}}\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{d}}^{(i,j)}_{p})}\boldsymbol{\mathit{y}}^{T}}$
$\displaystyle=2\frac{\boldsymbol{\mathit{x}}^{\prime
T}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}[V^{(i,j)}_{p}]}-\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{(i,j)}[V^{(i,j)}_{p}]})\boldsymbol{\mathit{y}}^{\prime}}{\boldsymbol{\mathit{x}}^{\prime
T}\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{d}}^{(i,j)}_{p})}\boldsymbol{\mathit{x}}^{\prime}+\boldsymbol{\mathit{y}}^{\prime}\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{d}}^{(i,j)}_{p})}\boldsymbol{\mathit{y}}^{\prime
T}}$ $\displaystyle=2\frac{\boldsymbol{\mathit{x}}^{\prime
T}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}[V^{(i,j)}_{p}]}-\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{(i,j)}[V^{(i,j)}_{p}]})\boldsymbol{\mathit{y}}^{\prime}}{\boldsymbol{\mathit{x}}^{\prime
T}\boldsymbol{\mathit{D}}_{H^{(i,j)}[V^{(i,j)}_{p}]}\boldsymbol{\mathit{x}}^{\prime}+\boldsymbol{\mathit{y}}^{\prime}\boldsymbol{\mathit{D}}_{H^{(i,j)}[V^{(i,j)}_{p}]}\boldsymbol{\mathit{y}}^{\prime
T}}$
$\displaystyle=\left\|\boldsymbol{\mathit{D}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+/2}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}[V^{(i,j)}_{p}]}-\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{(i,j)}[V^{(i,j)}_{p}]})\boldsymbol{\mathit{D}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+/2}\right\|$
(5)

where the last equality is by Lemma 3.14. Next we upper bound

$\displaystyle\left\|\boldsymbol{\mathit{D}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+/2}\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}[V^{(i,j)}_{p}]}\boldsymbol{\mathit{D}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+/2}\right\|$

and

$\displaystyle\left\|\boldsymbol{\mathit{D}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+/2}\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{(i,j)}[V^{(i,j)}_{p}]}\boldsymbol{\mathit{D}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+/2}\right\|.$

By Lemma 3.15 we have

$\displaystyle\left\|\boldsymbol{\mathit{D}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+/2}\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}[V^{(i,j)}_{p}]}\boldsymbol{\mathit{D}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+/2}\right\|\leq\max\left\\{\left\|\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}[V^{(i,j)}_{p}]}\boldsymbol{\mathit{D}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+}\right\|_{1},\left\|\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}[V^{(i,j)}_{p}]}^{T}\boldsymbol{\mathit{D}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+}\right\|_{1}\right\\}.$

Since for every $v$, the undirected graph
$\omega^{\min}_{\underrightarrow{G}}\cdot 2^{i+1}\cdot
H^{(i,j)}[V^{(i,j)}_{p}]$ satisfies

$\displaystyle\deg_{\omega^{\min}_{\underrightarrow{G}}\cdot 2^{i+1}\cdot
H^{(i,j)}[V^{(i,j)}_{p}]}(v)\geq\max\left\\{\deg^{+}_{\underrightarrow{G}^{(i,j)}[V^{(i,j)}_{p}]}(v),\deg^{-}_{\underrightarrow{G}^{(i,j)}[V^{(i,j)}_{p}]}(v)\right\\}$

we have

$\displaystyle\left\|\boldsymbol{\mathit{D}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+/2}\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}[V^{(i,j)}_{p}]}\boldsymbol{\mathit{D}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+/2}\right\|\leq\omega^{\min}_{\underrightarrow{G}}\cdot
2^{i+2}.$ (6)

Since we only used the in and out degrees of
$\underrightarrow{G}^{(i,j)}[V^{(i,j)}_{p}]$ in the above, and
$\tilde{\underrightarrow{G}}^{(i,j)}[V^{(i,j)}_{p}]$ has exactly the same
degrees, we analogously conclude

$\displaystyle\left\|\boldsymbol{\mathit{D}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+/2}\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{(i,j)}[V^{(i,j)}_{p}]}\boldsymbol{\mathit{D}}_{H^{(i,j)}[V^{(i,j)}_{p}]}^{+/2}\right\|\leq\omega^{\min}_{\underrightarrow{G}}\cdot
2^{i+2}.$ (7)

Chaining inequalities (2), (3), (4), (5), (6) and (7) yields

$\displaystyle\left\|(\boldsymbol{\mathit{U}}_{L_{\underrightarrow{G}}})^{+/2}(\underrightarrow{G}^{(i,j)}-\tilde{\underrightarrow{G}}^{(i,j)})(\boldsymbol{\mathit{U}}_{L_{\underrightarrow{G}}})^{+/2}\right\|_{2}\leq
128\cdot\operatorname{Exp}(2\cdot(\log n)^{\gamma})$

which concludes our proof. ∎

Next we analyse the number of edges of $\tilde{\underrightarrow{G}}^{(i,j)}$

######  Lemma 4.9.

$|E(\tilde{\underrightarrow{G}}^{(i,j)})|=O(n)$.

###### Proof.

It is easy to see that the patching routine adds at most $2|V^{(i,j)}_{p}|$
edges to graph $\tilde{\underrightarrow{G}}^{(i,j)}[V^{(i,j)}_{p}]$, since
each added edge repairs either the desired in-degree or the desired out-degree
of a vertex. The result follows since $\sum_{p}|V^{(i,j)}_{p}|=n$. ∎

We show the main lemma of this subsection by summing up the parts.

###### Proof of Lemma 4.7.

Notice that each expander decomposition peels of half of the edges, and thus
every edge is part of an unique expander part by the end of the procedure
Sparsify() in Algorithm 3. Thus our algorithm exactly preserves the in- and
out-degrees of
$\underrightarrow{G}_{1}=\beta\cdot\mathcal{U}(\underrightarrow{G})+\underrightarrow{G}$,
since each individual patching exactly preserves degrees. Since the graph
$\underrightarrow{G}_{2}=\beta\cdot\mathcal{U}(\underrightarrow{G})+\underrightarrow{R}$
remains connected the null-spaces are unaltered. Further, since
$\underrightarrow{R}$ is the sum of $\tilde{O(1)}$ graphs with $O(n)$ edges
the total amount of edges of $\underrightarrow{R}$ is bounded by
$\tilde{O}(n)$.

Finally, we show the approximation bound. By Lemma 3.13 we have

$\displaystyle\left\|\boldsymbol{\mathit{I}}_{\operatorname*{im}(\mathcal{U}^{(\beta)}(\underrightarrow{G}})-\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{2}}^{+}\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{1}}\right\|_{\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{2}}}\rightarrow\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{2}}}}\leq\left\|\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{2}}}^{+/2}(\boldsymbol{\mathit{L}}_{\underrightarrow{R}}-\boldsymbol{\mathit{L}}_{\underrightarrow{G}})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{2}}}^{+/2}\right\|.$

We use Lemma 3.14 to obtain

$\displaystyle\left\|\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{2}}}^{+/2}(\boldsymbol{\mathit{L}}_{\underrightarrow{R}}-\boldsymbol{\mathit{L}}_{\underrightarrow{G}})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{2}}}^{+/2}\right\|$
$\displaystyle=2\max_{\boldsymbol{\mathit{x}},\boldsymbol{\mathit{y}}\neq
0}\frac{\boldsymbol{\mathit{x}}^{T}(\boldsymbol{\mathit{L}}_{\underrightarrow{R}}-\boldsymbol{\mathit{L}}_{\underrightarrow{G}})\boldsymbol{\mathit{y}}}{\boldsymbol{\mathit{x}}^{T}(\beta\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}+\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}})\boldsymbol{\mathit{x}}+\boldsymbol{\mathit{y}}^{T}(\beta\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}+\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}})\boldsymbol{\mathit{y}}}$
$\displaystyle\leq 2\max_{\boldsymbol{\mathit{x}},\boldsymbol{\mathit{y}}\neq
0}\frac{\boldsymbol{\mathit{x}}^{T}(\boldsymbol{\mathit{L}}_{\underrightarrow{R}}-\boldsymbol{\mathit{L}}_{\underrightarrow{G}})\boldsymbol{\mathit{y}}}{\beta\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}\boldsymbol{\mathit{x}}+\beta\boldsymbol{\mathit{y}}^{T}\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}\boldsymbol{\mathit{y}}}.$

Using Lemma 3.14 again we have

$\displaystyle 2\max_{\boldsymbol{\mathit{x}},\boldsymbol{\mathit{y}}\neq
0}\frac{\boldsymbol{\mathit{x}}^{T}(\boldsymbol{\mathit{L}}_{\underrightarrow{R}}-\boldsymbol{\mathit{L}}_{\underrightarrow{G}})\boldsymbol{\mathit{y}}}{\beta\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}\boldsymbol{\mathit{x}}+\beta\boldsymbol{\mathit{y}}^{T}\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}\boldsymbol{\mathit{y}}}$
$\displaystyle=\frac{1}{\beta}\left\|\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}(\boldsymbol{\mathit{L}}_{\underrightarrow{R}}-\boldsymbol{\mathit{L}}_{\underrightarrow{G}})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}\right\|_{2}$
$\displaystyle=\frac{1}{\beta}\left\|\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}\sum_{i,j}(\boldsymbol{\mathit{L}}_{\underrightarrow{\tilde{G}}^{(i,j)}}-\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}\right\|_{2}$
$\displaystyle\leq\frac{1}{\beta}\sum_{i,j}\left\|\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}(\boldsymbol{\mathit{L}}_{\underrightarrow{\tilde{G}}^{(i,j)}}-\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}\right\|_{2}$
$\displaystyle\stackrel{{\scriptstyle i)}}{{\leq}}\frac{10\cdot P\cdot\log
n\cdot 128\cdot\operatorname{Exp}(2\cdot(\log
n)^{\gamma})}{\beta}\stackrel{{\scriptstyle ii)}}{{\leq}}\frac{1}{2}$

where i) follows from Lemma 4.8 and ii) is by $\beta\geq 128\cdot 20\cdot
P\cdot\log n\cdot\operatorname{Exp}(2\cdot(\log n)^{\gamma})$. Chaining the
inequalities shows the approximation statement and concludes our proof. ∎

###  4.3 Sparsifying the Undirected Part

The final task we have left is to sparsify the undirected graph
$\mathcal{U}(\underrightarrow{G})$. This can be more or less directly achieved
by employing a sparsification theorem presented in [CGL+20b]. Mainly for
notational convenience in our algorithm, we adapt this sparsifcation technique
to be degree preserving in Appendix B and state the resulting lemma here. See
Lemma B.2 for the proof.

######  Lemma 4.10.

(Degree Preserving Sparsification) There is a deterministic algorithm
$\textsc{SpectralSparsifyDeg}(G,\gamma)$ (Algorithm 10) that given a parameter
$\gamma\in(0,1)$ and an undirected graph $G=(V,E,\omega)$ with $n$ vertices
and $m$ edges such that that $P:=\frac{\max_{e\in E}\omega(e)}{\min_{e\in
E}\omega(e)}=poly(n)$ computes $\tilde{G}$ satsifying

  1. 1.

$\operatorname{Exp}(-(\log
n)^{\gamma})\boldsymbol{\mathit{L}}_{G}\preceq\boldsymbol{\mathit{L}}_{\tilde{G}}\preceq\operatorname{Exp}((\log
n)^{\gamma})\boldsymbol{\mathit{L}}_{G}$

  2. 2.

$\operatorname{nnz}(\boldsymbol{\mathit{A}})=\tilde{O}(n)$

in time

$\displaystyle\tilde{O}(m^{1+O(1/(\log n)^{\gamma/2})}\cdot(\log m)^{O((\log
n)^{\gamma})})=m^{1+o(1)}.$

The graph $\tilde{G}$ has self loops and exactly the same degrees as $G$.

Next we apply degree preserving sparsification together with an appropriate
scaling to sparsify the undirected part.

######  Lemma 4.11.

There exists a routine
$\tilde{G}=\textsc{SpectralSparsifyDeg}(\mathcal{U}(\underrightarrow{G}),\gamma)$
(Algorithm 10) that given the undirected graph
$\mathcal{U}(\underrightarrow{G})$ computes an undirected graph $\tilde{G}$
with $\tilde{O}(n)$ edges so that
$\boldsymbol{\mathit{L}}^{+}_{\underrightarrow{G}_{3}}=(\frac{\beta}{\eta}\boldsymbol{\mathit{L}}_{\tilde{G}}+\boldsymbol{\mathit{L}}_{\underrightarrow{R}})^{+}$
is an $\left(1-\frac{1}{2\cdot\operatorname{Exp}(4\cdot(\log
n)^{\gamma})}\right)$-approximate pseudoinverse of
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{2}}=\beta\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}+\boldsymbol{\mathit{L}}_{\underrightarrow{R}}$
with respect to
$\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{3}}}$.
The degrees of $\tilde{G}$ and $\mathcal{U}(\underrightarrow{G})$ are the
same.

###### Proof.

The sparsity and degree preservation follows directly from Lemma 4.10. To show
the approximation property, we use Lemma 3.13 and obtain

$\displaystyle\left\|\boldsymbol{\mathit{I}}_{\operatorname*{im}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}})}-\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{3}}^{+}\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{2}}\right\|_{\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{3}}}\rightarrow\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{3}}}}\leq\left\|\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{3}}}^{+/2}\left(\frac{\beta}{\eta}\boldsymbol{\mathit{L}}_{\tilde{G}}-\beta\boldsymbol{\mathit{L}}_{\mathcal{U}(\underrightarrow{G})}\right)\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{3}}}^{+/2}\right\|_{2}.$
(8)

Then we use Lemma 3.14 twice with
$\frac{\beta}{\eta}\boldsymbol{\mathit{L}}_{\tilde{G}}\preceq\frac{\beta}{\eta}\boldsymbol{\mathit{L}}_{\tilde{G}}+\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}}$
to obtain

$\displaystyle\left\|\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{3}}}^{+/2}\left(\frac{\beta}{\eta}\boldsymbol{\mathit{L}}_{\tilde{G}}-\beta\boldsymbol{\mathit{L}}_{\mathcal{U}(\underrightarrow{G})}\right)\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{3}}}^{+/2}\right\|_{2}$
$\displaystyle=2\max_{\boldsymbol{\mathit{x}},\boldsymbol{\mathit{y}}\neq
0}\frac{\boldsymbol{\mathit{x}}^{T}\left(\frac{\beta}{\eta}\boldsymbol{\mathit{L}}_{\tilde{G}}-\beta\boldsymbol{\mathit{L}}_{\mathcal{U}(\underrightarrow{G})}\right)\boldsymbol{\mathit{y}}}{\boldsymbol{\mathit{x}}^{T}(\frac{\beta}{\eta}\boldsymbol{\mathit{L}}_{\tilde{G}}+\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}})\boldsymbol{\mathit{x}}+\boldsymbol{\mathit{y}}^{T}(\frac{\beta}{\eta}\boldsymbol{\mathit{L}}_{\tilde{G}}+\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}})\boldsymbol{\mathit{y}}}$
$\displaystyle\leq 2\max_{\boldsymbol{\mathit{x}},\boldsymbol{\mathit{y}}\neq
0}\frac{\boldsymbol{\mathit{x}}^{T}\left(\frac{\beta}{\eta}\boldsymbol{\mathit{L}}_{\tilde{G}}-\beta\boldsymbol{\mathit{L}}_{\mathcal{U}(\underrightarrow{G})}\right)\boldsymbol{\mathit{y}}}{\boldsymbol{\mathit{x}}^{T}(\frac{\beta}{\eta}\boldsymbol{\mathit{L}}_{\tilde{G}})\boldsymbol{\mathit{x}}+\boldsymbol{\mathit{y}}^{T}(\frac{\beta}{\eta}\boldsymbol{\mathit{L}}_{\tilde{G}})\boldsymbol{\mathit{y}}}$
$\displaystyle=\left\|\boldsymbol{\mathit{L}}_{\tilde{G}}^{+/2}(\boldsymbol{\mathit{L}}_{\tilde{G}}-\eta\boldsymbol{\mathit{L}}_{\mathcal{U}(\underrightarrow{G})})\boldsymbol{\mathit{L}}_{\tilde{G}}^{+/2}\right\|.$
(9)

Next we compute

$\displaystyle\left\|\boldsymbol{\mathit{L}}_{\tilde{G}}^{+/2}(\boldsymbol{\mathit{L}}_{\tilde{G}}-\eta\boldsymbol{\mathit{L}}_{\mathcal{U}(\underrightarrow{G})})\boldsymbol{\mathit{L}}_{\tilde{G}}^{+/2}\right\|=\left\|\boldsymbol{\mathit{I}}_{\operatorname*{im}(\boldsymbol{\mathit{L}}_{\tilde{G}})}-\eta\boldsymbol{\mathit{L}}_{\tilde{G}}^{+/2}\boldsymbol{\mathit{L}}_{\mathcal{U}(\underrightarrow{G})}\boldsymbol{\mathit{L}}_{\tilde{G}}^{+/2}\right\|_{2}.$
(10)

We use the standard strategy of bounding the square. Let
$\boldsymbol{\mathit{M}}:=\boldsymbol{\mathit{L}}_{\tilde{G}}^{+/2}\boldsymbol{\mathit{L}}_{\mathcal{U}(\underrightarrow{G})}\boldsymbol{\mathit{L}}_{\tilde{G}}^{+/2}$
be a shorthand. Then we have

$\displaystyle\left\|\boldsymbol{\mathit{I}}_{\operatorname*{im}(\boldsymbol{\mathit{L}}_{\tilde{G}})}-\eta\boldsymbol{\mathit{M}}\right\|_{2}$
$\displaystyle=\max_{\boldsymbol{\mathit{x}}\in\operatorname*{im}(\boldsymbol{\mathit{L}}_{\tilde{G})}:\left\|\boldsymbol{\mathit{x}}\right\|_{2}=1}\boldsymbol{\mathit{x}}^{T}(\boldsymbol{\mathit{I}}_{\operatorname*{im}{\mathcal{U}(\underrightarrow{G})}}-\eta\boldsymbol{\mathit{M}})^{T}(\boldsymbol{\mathit{I}}_{\operatorname*{im}{\mathcal{U}(\underrightarrow{G})}}-\eta\boldsymbol{\mathit{M}})\boldsymbol{\mathit{x}}$
$\displaystyle=1+\max_{\boldsymbol{\mathit{x}}\in\operatorname*{im}{\boldsymbol{\mathit{L}}_{\tilde{G}}}:\left\|\boldsymbol{\mathit{x}}\right\|_{2}=1}-2\eta\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{M}}\boldsymbol{\mathit{x}}+\eta^{2}\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{M}}^{2}\boldsymbol{\mathit{x}}$
$\displaystyle\leq
1-2\eta\lambda_{*}(\boldsymbol{\mathit{M}})+\eta^{2}\left\|\boldsymbol{\mathit{M}}\right\|_{2}^{2}$

where $\lambda_{*}(\boldsymbol{\mathit{M}})$ denotes the smallest non-zero
eigenvalue of $\boldsymbol{\mathit{M}}$. We first lower bound
$\lambda_{*}(\boldsymbol{\mathit{M}})$. By Lemma 4.10 we have

$\displaystyle\frac{1}{2^{(\log
n)^{\gamma}}}\boldsymbol{\mathit{L}}_{\tilde{G}}$
$\displaystyle\preceq\boldsymbol{\mathit{L}}_{\mathcal{U}(\underrightarrow{G})}$
$\displaystyle\frac{1}{2^{(\log
n)^{\gamma}}}\boldsymbol{\mathit{I}}_{\operatorname*{im}(\tilde{G})}$
$\displaystyle\preceq\boldsymbol{\mathit{L}}_{\tilde{G}}^{+/2}\boldsymbol{\mathit{L}}_{\mathcal{U}(\underrightarrow{G})}\boldsymbol{\mathit{L}}_{\tilde{G}}^{+/2}$

and thus $\operatorname{Exp}(-(\log
n)^{\gamma})\leq\lambda_{*}(\boldsymbol{\mathit{M}})$. We obtain
$\left\|\boldsymbol{\mathit{M}}\right\|_{2}\leq\operatorname{Exp}(2\cdot(\log
n)^{\gamma})$ in an analogous way from Lemma 4.10. Then, we use
$\eta=\operatorname{Exp}(-3(\log
n)^{\gamma})\leq\frac{\lambda_{*}(\boldsymbol{\mathit{M}})}{\left\|M\right\|_{2}^{2}}$
to conclude.

$\displaystyle\left\|\boldsymbol{\mathit{I}}_{\operatorname*{im}(\boldsymbol{\mathit{L}}_{\tilde{G}})}-\eta\boldsymbol{\mathit{M}}\right\|_{2}\leq
1-\frac{1}{2\cdot\operatorname{Exp}(4\cdot(\log n)^{\gamma})}$ (11)

where we use $\sqrt{1-\epsilon}\leq 1-\epsilon/2$ for $\epsilon\in(0,1)$.
Chaining inequalities (8), (9), (10) and (11) shows the desired approximate
pseudoinverse property and finishes the proof. ∎

###  4.4 Loss in Condition Number

The overarching goal of the recursion of a squaring solver is to reduce the
condition number of the undirectification of the graph. When a constant
condition number is reached, the recursion terminates. However, since our
sparsification has relatively low accuracy, we have to ensure that our losses
do not outweigh our gains when we call it. In this subsection we bound the
loss in condition number incurred by our global sparsification routine. We
address the loss in condition number incurred by all of the three steps. To do
so we consider spectral bounds.

#### By Partial Symmetrization.

The first part is the easiest. Since
$\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\mathcal{U}^{(\beta)}(\underrightarrow{G})}}=(1+\beta)\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}$
$\beta$-partial-symmetrization does not change the condition number. Further
we observe

######  Observation 4.12.

$\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}\preceq\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\mathcal{U}^{(\beta)}(\underrightarrow{G})}}\preceq(1+\beta)\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}$

#### By Sparsifying the Directed Part.

Next we analyse the loss incurred by patching undirected expanders. This part
crucially relies on a spectral upper bound on
$\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}}$. We
show the following lemma

######  Lemma 4.13.

$\beta\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}\preceq\beta\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}+\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}}\preceq(\beta+\tilde{O}(1)\cdot\operatorname{Exp}(2(\log
n)^{\gamma}))\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}$

###### Proof.

The first inequality,
$\beta\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}\preceq\beta\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}+\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}}$,
is clear since
$\boldsymbol{0}\preceq\boldsymbol{\mathit{L}}_{\mathcal{U}(\underrightarrow{R})}=\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}}$.
To upper bound
$\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}}$ we
compute

$\displaystyle\left\|\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}(\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}}-\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}\right\|_{2}$
$\displaystyle\leq
2\left\|\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}}-\boldsymbol{\mathit{L}}_{\underrightarrow{R}})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}\right\|_{2}$
$\displaystyle=2\left\|\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}\sum_{i,j}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}}-\boldsymbol{\mathit{L}}_{\underrightarrow{\tilde{G}}^{(i,j)}})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}\right\|_{2}$
$\displaystyle\leq
2\sum_{i,j}\left\|\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{(i,j)}}-\boldsymbol{\mathit{L}}_{\underrightarrow{\tilde{G}}^{(i,j)}})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}\right\|_{2}$
$\displaystyle\leq\tilde{O}(1)\cdot\operatorname{Exp}(2(\log n)^{\gamma})$

where we use that the sum is over $\tilde{O}(1)$ edge weight buckets $i$ and
$\tilde{O}(1)$ expander decomposition layers $j$. We conclude

$\displaystyle\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}(\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}}-\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}$
$\displaystyle\preceq\tilde{O}(1)\cdot\operatorname{Exp}(2(\log
n)^{\gamma})\boldsymbol{\mathit{I}}_{\operatorname*{im}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}})}$
$\displaystyle\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}}\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}^{+/2}-\boldsymbol{\mathit{I}}_{\operatorname*{im}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}})}$
$\displaystyle\preceq\tilde{O}(1)\cdot\operatorname{Exp}(2(\log
n)^{\gamma})\boldsymbol{\mathit{I}}_{\operatorname*{im}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}})}$
$\displaystyle\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}}$
$\displaystyle\preceq\tilde{O}(1)\cdot\operatorname{Exp}(2(\log
n)^{\gamma})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}.$

This concludes our proof. ∎

#### By Sparsifying the Undirected Part.

Finally, we analyse the loss incurred by undirected sparsification. This part
directly follows from the approximation guarantee of Lemma 4.10.

######  Lemma 4.14.

For
$\tilde{G}=\textsc{SpectralSparsifyDeg}(\mathcal{U}(\underrightarrow{G}),\gamma)$
we have

$\displaystyle\operatorname{Exp}(-4(\log
n)^{\gamma})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{2}}}\preceq\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{3}}}\preceq\operatorname{Exp}(4(\log
n)^{\gamma})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{2}}}$

###### Proof.

Recall that
$\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{2}}}=\beta\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}+\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}}$
and
$\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{3}}}=\frac{\beta}{\eta}\boldsymbol{\mathit{L}}_{\tilde{G}}+\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}}$.
We show the following stronger statement

$\displaystyle\operatorname{Exp}(-4(\log
n)^{\gamma})\beta\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}+\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}}\preceq\frac{\beta}{\eta}\boldsymbol{\mathit{L}}_{\tilde{G}}+\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}}\preceq\operatorname{Exp}(4(\log
n)^{\gamma})\beta\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}+\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}}.$

After substracting
$\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{R}}}$ and
dividing by $\beta$ we have

$\displaystyle\operatorname{Exp}(-4(\log
n)^{\gamma})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}\preceq\frac{1}{\eta}\boldsymbol{\mathit{L}}_{\tilde{G}}\preceq\operatorname{Exp}(4(\log
n)^{\gamma})\boldsymbol{\mathit{U}}_{\boldsymbol{\mathit{L}}_{\underrightarrow{G}}}$

which directly follows from Lemma 4.10. ∎

###  4.5 Proof of Lemma 4.4

Now that we have assembled the pieces we are ready to prove Lemma 4.4, the
main lemma of this section.

###### Proof of Lemma 4.4.

We show each point in the enumeration seperately.

  1. 1.

The statement for $i=1$ is by Lemma 4.6 and
$\beta=\tilde{O}(1)\cdot\operatorname{Exp}(2(\log n)^{\gamma})$, for $i=2$ it
is by Lemma 4.7 and for $i=3$ it is by Lemma 4.11.

  2. 2.

The statement for $i=1$ is by 4.12, for $i=2$ it is by Lemma 4.13 and for
$i=3$ it is by Lemma 4.14.

  3. 3.

Clearly, $\underrightarrow{G}_{1}$ has at most twice as many edges as
$\underrightarrow{G}=\underrightarrow{G}_{0}$. Further,
$|E(\underrightarrow{G}_{2})|\leq|E(\underrightarrow{G}_{1})|+|E(\underrightarrow{R})|$
and $|E(\underrightarrow{R})|=\tilde{O}(n)$ by Lemma 4.7. Finally
$|E(\underrightarrow{G}_{3})|\leq|E(\tilde{G})|+|E(\underrightarrow{R})|=\tilde{O}(n)$
by Lemma 4.11 and Lemma 4.7. .

  4. 4.

Adding $\beta$ times the undirectifaction to an Eulerian graph increases the
in- and out- degrees by exactly a factor of $(1+\beta)$. Then the statement
follows from the degree preservation of our sparsification routines shown in
Lemma 4.7 and Lemma 4.11 and the extra scaling by $\frac{1}{\eta}$ of the
undirected part of $\underrightarrow{G}_{3}$

Finally, the runtime follows from Corollary 3.10. ∎

##  5 Approximately Squaring Directed Laplacians

Given a directed Eulerian Laplacian
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}}=\boldsymbol{\mathit{D}}_{\underrightarrow{G}}-\boldsymbol{\mathit{A}}^{T}_{\underrightarrow{G}}$
we call
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{2}}=\boldsymbol{\mathit{D}}_{\underrightarrow{G}}-\boldsymbol{\mathit{A}}^{T}_{\underrightarrow{G}}\boldsymbol{\mathit{D}}_{\underrightarrow{G}}^{-1}\boldsymbol{\mathit{A}}^{T}_{\underrightarrow{G}}$
its square. The squaring operation not only preserves degrees, but also
Eulerianness. In this section we introduce a deterministic algorithm for high
accuracy sparsification of the square of an Eulerian Laplacian in roughly
$\operatorname{nnz}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}})$ time.
Since squaring can drastically increase the density, the run-time of this
routine sometimes has to be sub-linear in the size of the graph it is
approximating. Thus it is clear that we can only afford work with implicit
representations of $\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{2}}$. Our
routines are based on and inspired by [KLP+15] and [PS21]. A similar
sparisification routine is developed in [AKM+20]. We first state the main
lemma of this section.

######  Lemma 5.1.

For every strongly connected Eulerian Laplacian
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}}=\boldsymbol{\mathit{D}}_{\underrightarrow{G}}-\boldsymbol{\mathit{A}}_{\underrightarrow{G}}^{T}$
the algorithm SparseSquare($\underrightarrow{G},\epsilon$) computes a strongly
connected Eulerian graph $\tilde{\underrightarrow{G}}^{2}$ in time
$O(\operatorname{nnz}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}})^{1+o(1)}/\epsilon^{4})$
so that
$\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}^{2}}=\boldsymbol{\mathit{D}}_{\tilde{\underrightarrow{G}}^{2}}-\boldsymbol{\mathit{A}}^{T}_{\tilde{\underrightarrow{G}}^{2}}$
is an $\epsilon$-graph-approximation of
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{2}}=\boldsymbol{\mathit{D}}_{\underrightarrow{G}}-\boldsymbol{\mathit{A}}^{T}_{\underrightarrow{G}}\boldsymbol{\mathit{D}}_{\underrightarrow{G}}^{-1}\boldsymbol{\mathit{A}}^{T}_{\underrightarrow{G}}$,
$\operatorname{nnz}(\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{2}})=O(\operatorname{nnz}(\boldsymbol{\mathit{L}}_{\tilde{\underrightarrow{G}}})/\epsilon^{4})$
and degrees are preserved.

As evident from the pseudocode of SparseSquare() (Algorithm 4) our approach to
implicitly sparsifying $\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{2}}$
breaks down its adjacency matrix into rank one contributions

$\displaystyle\boldsymbol{\mathit{A}}^{T}_{\underrightarrow{G}}\boldsymbol{\mathit{D}}_{\underrightarrow{G}}^{-1}\boldsymbol{\mathit{A}}^{T}_{\underrightarrow{G}}=\sum_{i=1}^{n}\underbrace{\boldsymbol{\mathit{D}}_{\underrightarrow{G}}^{-1}(i,i)(\boldsymbol{\mathit{A}}_{\underrightarrow{G}}(i,:))^{T}\cdot(\boldsymbol{\mathit{A}}_{\underrightarrow{G}}(:,i))^{T}}_{=:\boldsymbol{\mathit{A}}_{\underrightarrow{G}_{i}}}$

as commonplace in randomized squaring sparsification algorithms [PS21]. Then
SparseProduct() (Algorithm 6) computes
$\boldsymbol{\mathit{A}}_{\underrightarrow{H}_{i}}$, a suitable sparse
approximation of $\boldsymbol{\mathit{A}}_{\underrightarrow{G}_{i}}$. The full
adjacency matrix is obtained by simply summing up these contributions. Note
that SparseProduct() will preserve in- and out-degrees exactly, ensuring that
the resulting approximation preserves degrees, which in turn ensures
$\boldsymbol{\mathit{L}}_{\underrightarrow{G}^{2}}$ being Eulerian.

1 for  _$i=1,...,n$_ do

2
$\boldsymbol{\mathit{A}}_{\underrightarrow{H}_{i}}=\textsc{SparseProduct}(\boldsymbol{\mathit{A}}_{\underrightarrow{G}}(i,:),\boldsymbol{\mathit{A}}_{\underrightarrow{G}}(:,i),\epsilon)$

3

4 end for

5$\boldsymbol{\mathit{A}}_{\tilde{\underrightarrow{G}}^{2}}=\sum_{i=1}^{n}\boldsymbol{\mathit{A}}_{\underrightarrow{H}_{i}}$

6 return _$\tilde{\underrightarrow{G}}^{2}$_

Algorithm 4 SparseSquare($\underrightarrow{G},\epsilon$)

###  5.1 Deterministically Sparsifying Bipartite Product Graphs

High accuracy spectral sparsifiers for undirected bipartite product graphs are
known to the literature [KLP+15]. We will first adapt these to our setting by
showing that a simple greedy patching scheme can be employed to ensure the
conservation of degrees. Then we will directly use sparsified bipartite
product graphs to build SparseProduct() (Algorithm 6), our routine for
sparsifying $\boldsymbol{\mathit{L}}_{\underrightarrow{G}_{i}}$. We start by
defining the central object of this subsection.

######  Definition 5.2.

For $\boldsymbol{\mathit{a}}\in\mathbb{R}^{n}_{\geq 0}$ and
$\boldsymbol{\mathit{b}}\in\mathbb{R}^{n}_{\geq 0}$ positive vectors so that
$\left\|\boldsymbol{\mathit{a}}\right\|_{1}=\left\|\boldsymbol{\mathit{b}}\right\|_{1}=d$,
let

$\displaystyle\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}:=\begin{pmatrix}\operatorname*{diag}(\boldsymbol{\mathit{a}})&\boldsymbol{0}\\\
\boldsymbol{0}&\operatorname*{diag}(\boldsymbol{\mathit{b}})\end{pmatrix}-\begin{pmatrix}\boldsymbol{0}&\frac{\boldsymbol{\mathit{b}}\boldsymbol{\mathit{a}}^{T}}{d}\\\
\frac{\boldsymbol{\mathit{a}}\boldsymbol{\mathit{b}}^{T}}{d}&\boldsymbol{0}\end{pmatrix}$

be the Laplacian of the undirected bipartite product graph
$G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})$ of
$\boldsymbol{\mathit{a}}$ and $\boldsymbol{\mathit{b}}$.

######  Definition 5.3.

We call an undirected graph Laplacian

$\displaystyle\boldsymbol{\mathit{L}}_{\tilde{G}(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}=\begin{pmatrix}\operatorname*{diag}(\boldsymbol{\mathit{a}})&\boldsymbol{0}\\\
\boldsymbol{0}&\operatorname*{diag}(\boldsymbol{\mathit{b}})\end{pmatrix}-\begin{pmatrix}\boldsymbol{0}&\boldsymbol{\mathit{C}}\\\
\boldsymbol{\mathit{C}}^{T}&\boldsymbol{0}\end{pmatrix}$

with $\boldsymbol{\mathit{C}}\in\mathbb{R}^{n\times n}_{\geq 0}$ an
$\epsilon$-bipartite-sparsifier of
$\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}$
if

$\displaystyle(1-\epsilon)\boldsymbol{\mathit{L}}_{\tilde{G}(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}\preceq\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}\preceq(1+\epsilon)\boldsymbol{\mathit{L}}_{\tilde{G}(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}$

and
$\operatorname{nnz}(\boldsymbol{\mathit{L}}_{\tilde{G}(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})})=O((\operatorname{nnz}(\boldsymbol{\mathit{a}})+\operatorname{nnz}(\boldsymbol{\mathit{b}}))/\epsilon^{4})$.

Next we introduce the previously known spectral sparsification lemma we will
use.

######  Lemma 5.4 (Lemma G.16 in [KLP+15]).

There is a routine
$\textsc{WeightedBipartiteExpander}(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}},\epsilon)$
such that for any vectors
$\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}}\in\mathbb{R}^{n}_{\geq 0}$
and a parameter $\epsilon\in(0,1/2]$, it returns in time
$O((\operatorname{nnz}(\boldsymbol{\mathit{a}})+\operatorname{nnz}(\boldsymbol{\mathit{b}}))\epsilon^{-4})$
a bipartite graph $\tilde{G}(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})$
on the same bipartition as
$G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})$ with
$O((\operatorname{nnz}(\boldsymbol{\mathit{a}})+\operatorname{nnz}(\boldsymbol{\mathit{b}}))\epsilon^{-4})$
edges such that

$\displaystyle(1-\epsilon)\boldsymbol{\mathit{L}}_{\tilde{G}(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}\preceq\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}\preceq(1+\epsilon)\boldsymbol{\mathit{L}}_{\tilde{G}(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}$

The routine $\textsc{SparseBipartite}()$ (Algorithm 5) combines
WeightedBipartiteExpander() from Lemma 5.4 with the greedy patching procedure
BipartitePatching() (Algorithm 5) to construct an $\epsilon$-bipartite-
sparsifier. We first show the desired runtime and preservation of degrees.

1 Algorithm
_SparseBipartite(_$\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}},\epsilon$_)_

2 $\epsilon^{\prime}=\epsilon/128$

3
$H_{1}=\textsc{WeightedBipartiteExpander}(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}},\epsilon^{\prime})$

4 $H_{2}=H_{1}/(1+\epsilon^{\prime})$

5
$\boldsymbol{\mathit{d}}=\operatorname*{diag}(\boldsymbol{\mathit{D}}_{H_{2}})$

6
$\boldsymbol{\mathit{L}}_{R}=\textsc{BipartitePatching}(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}},\boldsymbol{\mathit{d}}(1:n),\boldsymbol{\mathit{d}}((n+1):2n))$

7 return _$\boldsymbol{\mathit{L}}_{H_{2}}+\boldsymbol{\mathit{L}}_{R}$_

8

9

10Procedure
_BipartitePatching(_$\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}},\boldsymbol{\mathit{d}}^{A},\boldsymbol{\mathit{d}}^{B}$_)_

11
$\tilde{\boldsymbol{\mathit{a}}}=\boldsymbol{\mathit{a}}-\boldsymbol{\mathit{d}}^{A}$

12
$\tilde{\boldsymbol{\mathit{b}}}=\boldsymbol{\mathit{b}}-\boldsymbol{\mathit{d}}^{B}$

13 $\boldsymbol{\mathit{L}}_{R}=\boldsymbol{0}^{2n\times 2n}$

14 while  _$\tilde{\boldsymbol{\mathit{a}}}\neq\boldsymbol{0}$_ do

15 Let $i$ and $j$ be arbitrary such that
$\tilde{\boldsymbol{\mathit{a}}}(i)>0$ and
$\tilde{\boldsymbol{\mathit{b}}}(j)>0$.

16
$w=\min\\{\tilde{\boldsymbol{\mathit{a}}}(i),\tilde{\boldsymbol{\mathit{b}}}(j)\\}$

17 $\tilde{\boldsymbol{\mathit{a}}}(i)=\tilde{\boldsymbol{\mathit{a}}}(i)-w$

18 $\tilde{\boldsymbol{\mathit{b}}}(j)=\tilde{\boldsymbol{\mathit{b}}}(j)-w$

19
$\boldsymbol{\mathit{L}}_{R}=\boldsymbol{\mathit{L}}_{R}+w(\boldsymbol{\mathit{e}}_{i}-\boldsymbol{\mathit{e}}_{j+n})(\boldsymbol{\mathit{e}}_{i}-\boldsymbol{\mathit{e}}_{j+n})^{T}$

20 end while

21 return _$\boldsymbol{\mathit{L}}_{R}$_

22

Algorithm 5
SparseBipartite($\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}},\epsilon$)
and
BipartitePatching($\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}},\boldsymbol{\mathit{d}}^{A},\boldsymbol{\mathit{d}}^{B}$)

######  Lemma 5.5.

The routine
$\boldsymbol{\mathit{L}}_{\tilde{G}(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}=\textsc{SparseBipartite}(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}},\epsilon)$
with $\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}}\in R^{n}_{\geq 0}$ so
that
$\left\|\boldsymbol{\mathit{a}}\right\|_{1}=\left\|\boldsymbol{\mathit{b}}\right\|_{1}$,
and $\epsilon\leq 1/2$ terminates in time
$O((nnz(\boldsymbol{\mathit{a}})+\operatorname{nnz}(\boldsymbol{\mathit{b}}))/\epsilon^{4})$
and ensures
$\boldsymbol{\mathit{D}}_{\tilde{G}(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}=\boldsymbol{\mathit{D}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}$.

###### Proof.

Since for any graph $G$ we have
$\boldsymbol{\mathit{e}}_{v}^{T}\boldsymbol{\mathit{L}}_{G}\boldsymbol{\mathit{e}}_{v}=\deg_{G}(v)$
for all $v$ it is easy to see that the spectral approximation of
WeightedBipartiteExpander() from Lemma 5.4 preserves degrees in the sense that
for all $v$

$\displaystyle(1-\epsilon^{\prime})\deg_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}(v)\leq\deg_{H_{1}}\leq(1+\epsilon^{\prime})\deg_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}.$

From this it is immediate that
$\deg_{H_{2}}(v)\leq\deg_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}(v)$.
Then the routine BipartitePatching() iteratively computes the Laplacian of a
bipartite patching graph so that
$\deg_{R}(v)=\deg_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}(v)-\deg_{H_{2}}(v)$
for all $v$. It is clear that both $H_{2}$ and
$G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})$ being bipartite ensures
the invariant that
$\left\|\tilde{\boldsymbol{\mathit{a}}}\right\|_{1}=\left\|\tilde{\boldsymbol{\mathit{b}}}\right\|_{1}$
throughout the execution of BipartitePatching() and
$\tilde{\boldsymbol{\mathit{a}}}\neq\boldsymbol{0}\implies\tilde{\boldsymbol{\mathit{b}}}\neq\boldsymbol{0}$.
Since each edge $(u,v)$ added to $R$ either ensures
$\tilde{\boldsymbol{\mathit{a}}}(u)=0$ or
$\tilde{\boldsymbol{\mathit{b}}}(v)=0$ the routine terminates after at most
$\operatorname{nnz}(\boldsymbol{\mathit{a}})+\operatorname{nnz}(\boldsymbol{\mathit{b}})$
edges are added since
$\operatorname{nnz}(\tilde{\boldsymbol{\mathit{a}}})+\operatorname{nnz}(\tilde{\boldsymbol{\mathit{b}}})\leq\operatorname{nnz}(\boldsymbol{\mathit{a}})+\operatorname{nnz}(\boldsymbol{\mathit{b}})$
initially. The result follows from these observations and Lemma 5.4. ∎

To conclude this subsection we would like to use Lemma 5.4 to show that
$\textsc{SparseBipartite}()$ (Algorithm 5) computes an $\epsilon$-bipartite-
sparsifier of
$\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}$.
The remaining obstruction is the quantification of the error introduced by the
BipartitePatching() (Algorithm 5). Our analysis of this error relies on the
high expansion of $G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})$.

Figure 5: This figure depicts a bipartite graph
$G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})$, where the set $S$ as
defined in Lemma 5.6 consits of the vertices highlighted with black dots. Then
the set is split into parts $S_{A}$ and $S_{B}$ as indicated.

######  Lemma 5.6.

$G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})$ is a $1/2$-expander.

###### Proof.

Consider any cut $\emptyset\subset S\subset V$. Without loss of generality,
let
$\operatorname{vol}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}(S)\leq\operatorname{vol}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}(V\setminus
S)$. We distinguish between the two parts $A$ and $B$ of the bipartite graph
and let $S_{A}\subseteq A$ and $S_{B}\subseteq B$ so that $S_{A}\cup S_{B}=S$.
We first introduce some terminology that will prove useful:

$\displaystyle
d:=\left\|\boldsymbol{\mathit{a}}\right\|_{1}(=\left\|\boldsymbol{\mathit{b}}\right\|_{1})\hskip
20.0ptt^{A}_{S_{A}}:=\sum_{v\in S_{A}}\boldsymbol{\mathit{a}}(v)\hskip
20.0ptt^{B}_{S_{B}}:=\sum_{v\in S_{B}}\boldsymbol{\mathit{b}}(v).$

See Figure 5 for a drawing of the situation. Then we have

$\displaystyle\delta_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}(S)$
$\displaystyle=\sum_{u\in S_{A}}\sum_{v\in B\setminus
S_{B}}\frac{\boldsymbol{\mathit{a}}(u)\boldsymbol{\mathit{b}}(v)}{\left\|\boldsymbol{\mathit{a}}\right\|_{1}}+\sum_{v\in
S_{B}}\sum_{u\in A\setminus
S_{a}}\frac{\boldsymbol{\mathit{a}}(u)\boldsymbol{\mathit{b}}(v)}{\left\|\boldsymbol{\mathit{a}}\right\|_{1}}$
$\displaystyle=\frac{1}{d}(t^{A}_{S_{A}}(d-t^{B}_{S_{B}})+t^{B}_{S_{B}}(d-t^{A}_{S_{A}}))$

and

$\displaystyle\operatorname{vol}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}(S)$
$\displaystyle=\sum_{u\in S_{A}}\sum_{v\in
B}\frac{\boldsymbol{\mathit{a}}(u)\boldsymbol{\mathit{b}}(v)}{\left\|\boldsymbol{\mathit{a}}\right\|_{1}}+\sum_{v\in
S_{B}}\sum_{u\in
A}\frac{\boldsymbol{\mathit{a}}(u)\boldsymbol{\mathit{b}}(v)}{\left\|\boldsymbol{\mathit{a}}\right\|_{1}}$
$\displaystyle=t^{A}_{S_{A}}+t^{B}_{S_{B}}.$

Finally we show

$\displaystyle\frac{\delta_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}(S)}{\operatorname{vol}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}(S)}\geq
1/2\iff(t^{A}_{S_{A}}(d-t^{B}_{S_{B}})+t^{B}_{S_{B}}(d-t^{A}_{S_{A}}))$
$\displaystyle\geq\frac{d}{2}(t^{A}_{S_{A}}+t^{B}_{S_{B}})$
$\displaystyle\iff\frac{d}{2}(t^{A}_{S_{A}}+t^{B}_{S_{B}})$ $\displaystyle\geq
2t^{B}_{S_{B}}t^{A}_{S_{A}}.$

Using
$\operatorname{vol}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}(S)\leq\operatorname{vol}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}(V\setminus
S)$ and
$d=\operatorname{vol}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}(V)$
we have $d\geq t^{A}_{S_{A}}+t^{B}_{S_{B}}$ and thus

$\displaystyle
d(t^{A}_{S_{A}}+t^{B}_{S_{B}})\geq(t^{A}_{S_{A}})^{2}+(t^{A}_{S_{A}})^{2}+2t^{A}_{S_{A}}t^{B}_{S_{B}}\geq
4t^{A}_{S_{A}}t^{B}_{S_{B}}$

by since for any real numbers $a$ and $b$ we have $a^{2}+b^{2}\geq 2ab$. ∎

Using the expansion of
$\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}$,
we give a spectral upper bound on the error by patching next.

######  Lemma 5.7.

In the context of
SparseBipartite($\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}},\epsilon$)
with $\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}}\in\mathbb{R}^{n}_{\geq
0}$,
$\left\|\boldsymbol{\mathit{a}}\right\|_{1}=\left\|\boldsymbol{\mathit{b}}\right\|_{1}$,
we have $\boldsymbol{\mathit{L}}_{R}\preceq
128\epsilon^{\prime}\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}$.

###### Proof.

As observed in the proof of Lemma 5.5 we have
$\deg_{R}(v)=\deg_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}(v)-\deg_{H_{2}}(v)$
for all $v$. Since
$\deg_{H_{2}}(v)\geq\frac{1-\epsilon^{\prime}}{1+\epsilon^{\prime}}\deg_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}(v)$
by Lemma 5.4 we have

$\displaystyle\deg_{R}(v)\leq\deg_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}-\frac{1-\epsilon^{\prime}}{1+\epsilon^{\prime}}\deg_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}(v)\leq
4\epsilon^{\prime}\deg_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}(v).$
(12)

In the following, let
$\boldsymbol{\mathit{d}}=\operatorname*{diag}(\boldsymbol{\mathit{D}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})})$
be the degree vector of $G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})$.
We use Lemma 3.14 and arrive at

$\displaystyle\left\|\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}^{+/2}\boldsymbol{\mathit{L}}_{R}\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}^{+/2}\right\|_{2}=2\max_{\boldsymbol{\mathit{x}},\boldsymbol{\mathit{y}}\neq\boldsymbol{0}}\frac{\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{L}}_{R}\boldsymbol{\mathit{y}}}{\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}\boldsymbol{\mathit{x}}+\boldsymbol{\mathit{y}}^{T}\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}\boldsymbol{\mathit{y}}}$
(13)

Let $\boldsymbol{\mathit{x}}$ and $\boldsymbol{\mathit{y}}$ be maximizing the
right hand side of (13) and
$\boldsymbol{\mathit{x}},\boldsymbol{\mathit{y}}\perp\boldsymbol{1}$. Then,
$\boldsymbol{\mathit{x}}^{\prime}=\boldsymbol{\mathit{x}}-\frac{\boldsymbol{\mathit{x}}^{T}\boldsymbol{\mathit{d}}}{\left\|\boldsymbol{\mathit{d}}\right\|_{2}}\boldsymbol{1}$
and
$\boldsymbol{\mathit{y}}^{\prime}=\boldsymbol{\mathit{y}}-\frac{\boldsymbol{\mathit{y}}^{T}\boldsymbol{\mathit{d}}}{\left\|\boldsymbol{\mathit{d}}\right\|_{2}}\boldsymbol{1}$
are also maximizing since the all-ones vector is in the kernel of Laplacian
matrices. By Corollary 3.7 and Lemma 5.6 we have
$\frac{1}{16}(\operatorname*{diag}(\boldsymbol{\mathit{d}})-\frac{\boldsymbol{\mathit{d}}\boldsymbol{\mathit{d}}^{T}}{\left\|\boldsymbol{\mathit{d}}\right\|_{1}})\preceq\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}$.
We calculate

$\displaystyle\left\|\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}^{+/2}\boldsymbol{\mathit{L}}_{R}\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}^{+/2}\right\|_{2}$
$\displaystyle=2\frac{\boldsymbol{\mathit{x}}^{\prime
T}\boldsymbol{\mathit{L}}_{R}\boldsymbol{\mathit{y}}^{\prime}}{\boldsymbol{\mathit{x}}^{\prime
T}\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}\boldsymbol{\mathit{x}}^{\prime}+\boldsymbol{\mathit{y}}^{\prime
T}\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}\boldsymbol{\mathit{y}}^{\prime}}$
$\displaystyle\leq 32\frac{\boldsymbol{\mathit{x}}^{\prime
T}\boldsymbol{\mathit{L}}_{R}\boldsymbol{\mathit{y}}^{\prime}}{\boldsymbol{\mathit{x}}^{\prime
T}(\operatorname*{diag}(\boldsymbol{\mathit{d}})-\frac{\boldsymbol{\mathit{d}}\boldsymbol{\mathit{d}}^{T}}{\left\|\boldsymbol{\mathit{d}}\right\|_{1}})\boldsymbol{\mathit{x}}^{\prime}+\boldsymbol{\mathit{y}}^{\prime
T}(\operatorname*{diag}(\boldsymbol{\mathit{d}})-\frac{\boldsymbol{\mathit{d}}\boldsymbol{\mathit{d}}^{T}}{\left\|\boldsymbol{\mathit{d}}\right\|_{1}})\boldsymbol{\mathit{y}}^{\prime}}$
$\displaystyle=32\frac{\boldsymbol{\mathit{x}}^{\prime
T}\boldsymbol{\mathit{L}}_{R}\boldsymbol{\mathit{y}}^{\prime}}{\boldsymbol{\mathit{x}}^{\prime
T}\operatorname*{diag}(\boldsymbol{\mathit{d}})\boldsymbol{\mathit{x}}^{\prime}+\boldsymbol{\mathit{y}}^{\prime
T}\operatorname*{diag}(\boldsymbol{\mathit{d}})\boldsymbol{\mathit{y}}^{\prime}}$
$\displaystyle\leq
16\left\|\operatorname*{diag}(\boldsymbol{\mathit{d}})^{+/2}\boldsymbol{\mathit{L}}_{R}\operatorname*{diag}(\boldsymbol{\mathit{d}})^{+/2}\right\|$

where the last line follows from another application of Lemma 3.14. By Lemma
3.15 and (12) we have

$\displaystyle\left\|\operatorname*{diag}(\boldsymbol{\mathit{d}})^{+/2}\boldsymbol{\mathit{L}}_{R}\operatorname*{diag}(\boldsymbol{\mathit{d}})^{+/2}\right\|\leq\left\|\boldsymbol{\mathit{L}}_{R}\boldsymbol{\mathit{D}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}\right\|_{1}\leq
8\epsilon^{\prime}.$

We conclude

$\displaystyle\boldsymbol{\mathit{L}}_{R}\preceq
128\epsilon^{\prime}\boldsymbol{\mathit{L}}_{G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})}.$

∎

Finally, we assemble the pieces and show that SparseBipartite() sparsifies
$G(\boldsymbol{\mathit{a}},\boldsymbol{\mathit{b}})$ while preserving degrees.

######  Lemma 5.8.

There is a routine