2\operatorname{Tr}(\mathbf{U}^{\top}\mathbf{U})-\operatorname{Tr}(\mathbf{U}^{\top}\mathbf{U})=\operatorname{Tr}(\mathbf{U}^{\top}\mathbf{U}).\qed$

######  Proposition 2.15.

Let $\mathbf{X}$ be a random $d\times d$ matrix which is a.s. PSD, and let
$\mathbf{A}$ be fixed matrix which is PSD and non-singular. Suppose that for
every $\mathbf{U}$ in $\mathbb{R}^{d\times m}$ it holds that

$\operatorname{Pr\/}\left(\operatorname{Tr}(\mathbf{U}^{\top}\mathbf{X}\mathbf{U})\leq\operatorname{Tr}(\mathbf{U}^{\top}\mathbf{A}\mathbf{U})\right)\geq
1-\delta.$

Then for every $\mathbf{V}$ in $\mathbb{R}^{d\times m}$ it also holds that

$\operatorname{Pr\/}\left(\operatorname{Tr}(\mathbf{V}^{\top}\mathbf{X}^{-1}\mathbf{V})\geq\operatorname{Tr}(\mathbf{V}^{\top}\mathbf{A}^{-1}\mathbf{V})\right)\geq
1-\delta.$

###### Proof.

Plugging $\mathbf{A}^{-1}\mathbf{V}$ into the hypothesis gives that for all
$\mathbf{V}$,

$\operatorname{Pr\/}(\operatorname{Tr}(\mathbf{V}^{\top}\mathbf{A}^{-1}\mathbf{X}\mathbf{A}^{-1}\mathbf{V})\leq\operatorname{Tr}(\mathbf{V}^{\top}\mathbf{A}^{-1}\mathbf{V}))\geq
1-\delta,$

or equivalently

$\operatorname{Pr\/}\left(\operatorname{Tr}((\mathbf{A}^{-1/2}\mathbf{V})^{\top}\mathbf{A}^{-1/2}\mathbf{X}\mathbf{A}^{-1/2}(\mathbf{A}^{-1/2}\mathbf{V}))\leq\operatorname{Tr}((\mathbf{A}^{-1/2}\mathbf{V})^{\top}(\mathbf{A}^{-1/2}\mathbf{V}))\right)\geq
1-\delta.$

By Proposition 2.14, this gives that for all $V$,

$\operatorname{Pr\/}\left(\operatorname{Tr}((\mathbf{A}^{-1/2}\mathbf{V})^{\top}\mathbf{A}^{1/2}\mathbf{X}^{-1}\mathbf{A}^{1/2}(\mathbf{A}^{-1/2}\mathbf{V}))\geq\operatorname{Tr}((\mathbf{A}^{-1/2}\mathbf{V})^{\top}(\mathbf{A}^{-1/2}\mathbf{V}))\right)\geq
1-\delta,$

which simplifies to the desired conclusion. ∎

######  Lemma 2.16.

Let $\mathbf{A}=[\mathbf{A}^{(1)};\ldots;\mathbf{A}^{(s)}]$ be non-singular
and let $\mathbf{S}^{(1)},\ldots\mathbf{S}^{(s)}$ be random sketching matrices
of appropriate dimension so that the products $S^{(i)}A^{(i)}$ are defined.
Assume that each $S^{(i)}$ satisfies the $(1,\delta,2)$-JL-moment property.
Let $\mathbf{V}\in\mathbb{R}^{d\times m}$. Then with probability at least
$1-\delta$,

$\operatorname{Tr}\left(\mathbf{V}^{\top}\left(\sum_{i=1}^{s}{\mathbf{A}^{(i)}}^{\top}{\mathbf{S}^{(i)}}^{\top}\mathbf{S}^{(i)}\mathbf{A}^{(i)}\right)^{-1}\mathbf{V}\right)\geq\frac{1}{2}\operatorname{Tr}\left(\mathbf{V}^{\top}(\mathbf{A}^{\top}\mathbf{A})^{-1}\mathbf{V}\right).$

###### Proof.

We apply Proposition 2.15. Let $\mathbf{U}\in\mathbb{R}^{d\times m}$ be an
arbitrary fixed matrix. Then we have

$\operatorname{Tr}\left(\mathbf{U}^{\top}\left(\sum_{i=1}^{s}{\mathbf{A}^{(i)}}^{\top}{\mathbf{S}^{(i)}}^{\top}\mathbf{S}^{(i)}\mathbf{A}^{(i)}\right)\mathbf{U}\right)=\sum_{i=1}^{s}\|\mathbf{S}^{(i)}\mathbf{A}^{(i)}\mathbf{U}\|_{F}^{2}.$

The block matrix $\mathbf{S}^{(1)}\oplus\cdots\oplus\mathbf{S}^{(s)}$ also has
the $(1,\delta,2)$-JL-moment property (see for example Lemma 13 of [Ahl+20]).
So with probability at least $1-\delta$,

$\sum_{i=1}^{s}\|\mathbf{S}^{(i)}\mathbf{A}^{(i)}\mathbf{U}\|_{F}^{2}\leq
2\|\mathbf{A}\mathbf{U}\|_{F}^{2}=2\operatorname{Tr}(\mathbf{U}^{\top}\mathbf{A}^{\top}\mathbf{A}\mathbf{U}).$

The claim now follows by Proposition 2.15. ∎

######  Lemma 2.17.

Let $\mathbf{A}=[\mathbf{A}^{(1)};\dots,\mathbf{A}^{(s)}]$, let
$\mathbf{S}^{(1)},\ldots,\mathbf{S}^{(s)}$ be sketching matrices satisfying
the $(1,\delta,2)$-JL-moment property, and let
$\widetilde{\mathbf{A}}=[\mathbf{S}^{(1)}\mathbf{A}^{(1)};\ldots;\mathbf{S}^{(s)}\mathbf{A}^{(s)}],$
where $\mathbf{S}_{k}=\mathbf{I}$ for a fixed $k$. Then with probability at
least $1-\delta$,
$\mathcal{L}_{k}(\widetilde{\mathbf{A}})\geq\frac{1}{2}\mathcal{L}_{k}(\mathbf{A}).$

###### Proof.

By Lemma 2.16 (and the hypothesis that $\mathbf{S}^{(k)}=\mathbf{I}$), we have

$\mathcal{L}_{k}(\widetilde{\mathbf{A}})=\operatorname{Tr}\left(\mathbf{A}^{(k)}\left(\sum_{i=1}^{s}{\mathbf{A}^{(i)}}^{\top}{\mathbf{S}^{(i)}}^{\top}{\mathbf{S}^{(i)}}^{\top}\mathbf{A}^{(i)}\right)^{-1}\mathbf{A}^{(k)}\right)\geq\frac{1}{2}\operatorname{Tr}\left(\mathbf{A}^{(k)}\left(\sum_{i=1}^{s}{\mathbf{A}^{(i)}}^{\top}\mathbf{A}^{(i)}\right)^{-1}\mathbf{A}^{(k)}\right)=\frac{1}{2}\mathcal{L}_{k}(\mathbf{A}).$

∎

Combining the two block sketching results Lemma 2.13 and Lemma 2.17 gives the
following.

######  Lemma 2.18.

Let $\mathbf{S}^{(1)},\ldots,\mathbf{S}^{(s)}$, each with $d$ columns, all be
(normalized) Rademacher with $O(k\log(s/\delta))$ rows. Then for each $i$,

$\mathcal{L}_{i}([\mathbf{S}^{(1)}\mathbf{A}^{(1)};\ldots;\mathbf{S}^{(s)}\mathbf{A}^{(s)}])\geq
C\min\left(k,\mathcal{L}_{i}([\mathbf{A}^{(1)};\ldots;\mathbf{A}^{(s)}])\right),$

with probability at least $1-\delta$.

######  Remark 2.19.

The sketches in the above result were taken to be Rademacher only for
convenience. The same argument applies to sparse sketches for example.

####  2.2.3 Estimating Block Leverage Scores

The sketch from the previous section shows that we can accurately
(over-)estimate a given block leverage score by sketching down to dimension
roughly $k.$ Unfortunately the block leverage scores can be as large as $d$
and we are unable to take a sketch of $d$ rows from all servers. Fortunately,
not many servers can have large block leverage score, so by iteratively
pruning off the ones that don’t, we can focus on the servers with the most
information.

Input.  For $i\in[s]$, each server $i$ has the block $\mathbf{A}^{(i)}$ of
$\mathbf{A}$

Output.  List $L$ where $L[i]$ estimates of the block leverage score of
$\mathbf{A}^{(i)}$

Initialize. $L\leftarrow[\perp,\ldots,\perp]$ of length $s$,
$\mathcal{S}_{0}\leftarrow\\{1,\ldots,s\\}$

For iterations $r=0,1,2,\dotsc,\lceil\log d\rceil$:

  1. 1

$k_{r}\leftarrow 2^{r}$

  2. 2

For  $i\in\mathcal{S}_{r}$

    1. (a)

Each server $i$ draws $\mathbf{S}_{r,i}$ a constant distortion
$k_{r}$-dimensional oblivious subspace embedding, and sends
$\mathbf{S}_{r,i}\mathbf{A}^{(i)}$ to the coordinator

  3. 3

Coordinator forms block matrix $\mathbf{A}_{\langle r\rangle}$ with blocks
given by $\mathbf{S}_{r,i}\mathbf{A}^{(i)}$ for $i$ in $\mathcal{S}_{r}$, and
where the blocks are indexed by $\mathcal{S}_{r}$

  4. 4

For all $i\in\mathcal{S}_{r}$, coordinator computes
$\widehat{\mathcal{L}}_{r,i}=\mathcal{L}_{i}(\mathbf{A}_{\langle r\rangle})$

  5. 5

$\mathcal{S}_{r+1}\leftarrow\\{i\in\mathcal{S}_{r}:\widehat{\mathcal{L}}_{r,i}\geq
Ck_{r}\\}$

  6. 6

For  $i$ in $\mathcal{S}_{r}\setminus\mathcal{S}_{r+1}$

    1. (a)

$L[i]\leftarrow\widehat{\mathcal{L}}_{r,i}$

For all $i\in[s]$, if $L[i]=\perp$, then set $L[i]\leftarrow d$

Return: $L$.

Algorithm 4: Block leverage score estimation.

alg]alg:block_lev_est

######  Theorem 2.20.

LABEL:alg:block_lev_est runs with $O(\log d)$ rounds of communication, and
returns a list $L$ satisfying

  1. (i)

$L[i]\geq C\mathcal{L}_{i}([\mathbf{A}^{(1)};\ldots;\mathbf{A}^{(s)}])$ for
all $i$

  2. (ii)

$\sum_{i=1}^{s}L[i]\leq O(d\log d)$

Moreover the servers collectively send at most $O(cs+cd\log d)$ vectors of
length $d$ to the coordinator.

###### Proof.

We start by bounding the number of servers which are active in a given round.
In round $0$, $|\mathcal{S}_{0}|=s.$ For $r\geq 1$, note that for every server
$i$ in $\mathcal{S}_{r}$, $\widehat{\mathcal{L}}_{r-1,i}\geq Ck_{r-1}.$ On the
other hand there cannot be many such servers, since by Proposition 2.11,

$\sum_{i\in\mathcal{S}_{r}}\widehat{\mathcal{L}}_{r-1,i}\leq d,$

which implies that $|\mathcal{S}_{r}|\leq\frac{d}{Ck_{r-1}}.$ This immediately
gives a bound on the communication cost. Summing the number of vectors
transmitted in each round gives a total of

$\sum_{i=0}^{\lceil\log d\rceil}ck_{r}|\mathcal{S}_{r}|\leq
c\left(sk_{0}+\sum_{i=1}^{\lceil\log
d\rceil}k_{r}\frac{d}{Ck_{r-1}}\right)=c\left(s+\frac{2d}{C}\lceil\log
d\rceil\right)$

vectors sent to the coordinator. Next we show that (i) holds. By the
algorithm, note that either $L[i]=d$ or on round $r$ we set
$L[i]=\widehat{\mathcal{L}}_{r,i}$. In the first case (i) is trivial since all
block leverage scores are at most $d.$ In the latter case,
$\widehat{\mathcal{L}}_{r,i}\leq Ck_{r}.$ But $\widehat{\mathcal{L}}_{r,i}\geq
C\min\left(k_{r},\mathcal{L}_{i}(\mathbf{A}_{\langle r\rangle})\right)$ by
Lemma 2.18, which is at least $\min(Ck_{r},C\mathcal{L}_{i}(\mathbf{A}))$ by
monotonicity. So we have $k_{r}\geq\mathcal{L}_{i}(\mathbf{A}),$ which implies
that

$\widehat{\mathcal{L}}_{r,i}\geq
C\min(k_{r},\mathcal{L}_{i}(A))=C\mathcal{L}_{i}(\mathbf{A}).$

Finally we show (ii). Since we have

$\sum_{i\in\mathcal{S}_{r}}\widehat{\mathcal{L}}_{r,i}=\sum_{i\in\mathcal{S}_{r}}\mathcal{L}_{i}(\mathbf{A}_{\langle
r\rangle})\leq d,$

it follows that the entries of $L$ which are set in round $r$ sum to at most
$d$. Hence the sum of the entries of $L$ set in the outer for-loop is at most
$(\lceil\log d\rceil+1)d.$ By the argument given above for the communication
cost, there are at most $\frac{d}{Ck_{r}}\leq\frac{1}{C}$ entries of $L$ which
are not set after the loop. These entries are set to $d$, which gives
$\sum_{i=1}^{s}L[i]\leq(\lceil\log d\rceil+1)d+\frac{d}{C}\leq O(d\log d).$ ∎

####  2.2.4 Sampling via Block Leverage Scores

Given the overestimates computed for the block leverage scores in the previous
section, a straightforward concentration bound allows to get a spectral
approximation via block leverage sampling. By combining with the algorithm for
estimating the block leverage scores, this immediately yields an algorithm
with $\tilde{O}(sd+d^{2}\varepsilon^{-2})$ communication in the coordinator
model, for computing an $\varepsilon$ distortion subspace embedding for the
columns of $\mathbf{A}.$

Input.  Sampling distribution $p$, number of samples $N$

Output.  A spectral approximation $\widehat{\mathbf{A}}$ of $\mathbf{A}$

Initialize. Coordinator sets $\widehat{\mathbf{A}}=\mathbf{0}$
$\in\mathbb{R}^{N\times d}$

For iterations $\ell=1,2,\dotsc,N$:

  1. 1

Sample server $j$ from the distribution $p$

  2. 2

Server $j$ generates a Rademacher random vector
$\mathbf{g}\in\mathbb{R}^{m_{j}}$ and sends
$\mathbf{g}^{\top}\mathbf{A}^{(j)}$ to coordinator

  3. 3

Coordinator appends $\mathbf{g}^{\top}\mathbf{A}^{(j)}$ to
$\widehat{\mathbf{A}}$

Return: $\widehat{\mathbf{A}}$.

Algorithm 5: Block leverage score sampling.

alg]alg:block_lev_sampling

As is standard for analyses of leverage score sampling, we rely on the Matrix
Chernoff bound (see [Woo+14] for example). We state a version here which
follows from [Tro12]. The version we use is slightly less general, but more
convenient for our purposes.

######  Theorem 2.21.

Let $\mathbf{X}_{1},\ldots,\mathbf{X}_{d}\in\mathbb{R}^{d\times d}$ be random
matrices which are independent and symmetric PSD, with
$\mu_{\min}I\leq\mathbb{E}\mathbf{X}_{i}\leq\mu_{\max}I,$ and
$\|\mathbf{X}_{i}\|\leq R$ a.s. Let
$\overline{\mathbf{X}}=\frac{1}{N}\sum_{i=1}^{N}\mathbf{X}_{i}.$ Then for all
$\delta\in[0,1),$

$\operatorname{Pr\/}\lparen\lambda_{\max}(\overline{\mathbf{X}})\geq(1+\delta)\mu_{\max}\rparen\leq
d\exp\left(-\delta^{2}\frac{N\mu_{\max}}{3R}\right).$

$\operatorname{Pr\/}\lparen\lambda_{\max}(\overline{\mathbf{X}})\leq(1-\delta)\mu_{\min}\rparen\leq
d\exp\left(-\delta^{2}\frac{N\mu_{\min}}{2R}\right).$

The quantities that we apply the matrix Chernoff bound to will have operator
norm given by a Hutchinson trace estimator [Hut89]. This standard application
of Matrix Chernoff is the core argument. The additional work simply fixes a
technical issue.

Hutchinson’s estimator may very occasionally be much larger than expected,
which would require $R$ in Theorem 2.21 to be undesirably large. Fortunately
Hutchinson’s estimator has exponential tail decay, and so these potential
large values may be safely ignored with high probability. To make this
precise, we will use the following technical fact, which is effectively a
restatement of results in [DM21].

######  Proposition 2.22.

Let $\mathbf{A}\in\mathbb{R}^{d\times d}$ be symmetric positive semidefinite,
and let $g\in\mathbb{R}^{d}$ be a Rademacher random vector. Let
$\mu=\mathbb{E}(\mathbf{g}^{\top}\mathbf{A}\mathbf{g})=\operatorname{Tr}(\mathbf{A}).$
Then

$\operatorname{Pr\/}(\mathbf{g}^{\top}\mathbf{A}\mathbf{g}\geq t)\leq
c_{1}e^{-c_{2}t/\mu},$

for all $t\geq c_{3}\mu.$ The $c_{i}$’s are positive absolute constants.

###### Proof.

We set $\ell=1$ in Claim A.3 of [DM21]. By bounding $\|\mathbf{A}\|_{2}$ and
$\|\mathbf{A}\|_{F}$ each by $\operatorname{Tr}(\mathbf{A})=\mu$, we may take
$\nu=c_{4}\mu$ and $\beta=c_{5}\mu$ in Claim A.3. By properties of
subexponential random variables given in [Wai15], it then follows that
$\operatorname{Pr\/}(\mathbf{g}^{\top}\mathbf{A}\mathbf{g}\geq\mu+t)\leq
2e^{-c_{6}t/\mu}$ for $t\geq\nu^{2}/\beta=c_{7}\mu,$ which by adjusting
constants rearranges to claim above. ∎

######  Theorem 2.23.

Suppose that the input to LABEL:alg:block_lev_sampling satisfies
$p_{i}\geq\beta\frac{\mathcal{L}_{i}(\mathbf{A})}{d}$ for some
$\beta\in(0,1]$, and with
$N\geq\Omega\left(\frac{d}{\beta\varepsilon^{2}}\log\lparen\frac{d}{\beta\varepsilon}\rparen\log
d\right),$ where $\varepsilon<1$. Then the output $\widehat{\mathbf{A}}$ of
LABEL:alg:block_lev_sampling satisfies

$(1-\varepsilon)\mathbf{A}^{\top}\mathbf{A}\leq\widehat{\mathbf{A}}^{\top}\widehat{\mathbf{A}}\leq(1+\varepsilon)\mathbf{A}^{\top}\mathbf{A}.$

###### Proof.

Let $X_{k}$ be distributed as
$\frac{1}{p_{j}}\mathbf{A}_{j}^{\top}\mathbf{g}\mathbf{g}^{\top}\mathbf{A}_{j}$
where the index $j$ is drawn from $p$, and $\mathbf{g}$ is independently drawn
as a Rademacher random vector. Then
$\widehat{\mathbf{A}}^{\top}\widehat{\mathbf{A}}$ is distributed as
$\frac{1}{N}\sum_{k=1}^{N}\mathbf{X}_{k}$, so we show concentration for this
average. As is standard in such arguments, we show that the following
equivalent statement holds with the desired probability:

$(1-\varepsilon)\mathbf{I}\leq(\mathbf{A}^{\top}\mathbf{A})^{-1/2}\mathbf{X}(\mathbf{A}^{\top}\mathbf{A})^{-1/2}\leq(1+\varepsilon)\mathbf{I}.$

First note that

$\mathbb{E}(\mathbf{X}_{k})=\sum_{i=1}^{s}p_{i}\mathbb{E}_{\mathbf{g}}(\frac{1}{p_{i}}{\mathbf{A}^{(i)}}^{\top}\mathbf{g}^{\top}\mathbf{g}{\mathbf{A}^{(i)}})=\sum_{i=1}^{s}{\mathbf{A}^{(i)}}^{\top}\mathbb{E}_{\mathbf{g}}(\mathbf{g}^{\top}\mathbf{g}){\mathbf{A}^{(i)}}=\sum_{i=1}^{s}{\mathbf{A}^{(i)}}^{\top}{\mathbf{A}^{(i)}}=\mathbf{A}^{\top}\mathbf{A},$

since $\mathbb{E}(\mathbf{g}^{\top}\mathbf{g})=I.$ Let
$\mathbf{Y}_{k}=(\mathbf{A}^{\top}\mathbf{A})^{-1/2}\mathbf{X}_{k}(\mathbf{A}^{\top}\mathbf{A})^{-1/2}$,
and note that by the above, $\mathbb{E}(\mathbf{Y}_{k})=\mathbf{I}.$ Next we
have

$\|\mathbf{Y}_{k}\|=\|(\mathbf{A}^{\top}\mathbf{A})^{-1/2}\lparen\frac{1}{p_{i}}{\mathbf{A}^{(i)}}^{\top}\mathbf{g}\mathbf{g}^{\top}{\mathbf{A}^{(i)}}\rparen(\mathbf{A}^{\top}\mathbf{A})^{-1/2}\|=\mathbf{g}^{\top}\lparen\frac{1}{p_{i}}{\mathbf{A}^{(i)}}(\mathbf{A}^{\top}\mathbf{A})^{-1}{\mathbf{A}^{(i)}}^{\top}\rparen\mathbf{g}.$

For fixed $i$, this latter expression is the classic Hutchinson’s trace
estimator for
$\frac{1}{p_{i}}{\mathbf{A}^{(i)}}(\mathbf{A}^{\top}\mathbf{A})^{-1}{\mathbf{A}^{(i)}}^{\top}$,
which has mean

$\operatorname{Tr}\left(\frac{1}{p_{i}}{\mathbf{A}^{(i)}}(\mathbf{A}^{\top}\mathbf{A})^{-1}{\mathbf{A}^{(i)}}^{\top}\right)=\frac{1}{p_{i}}\mathcal{L}_{i}(\mathbf{A})\leq\frac{d}{\beta}.$

So by Proposition 2.22,

$\operatorname{Pr\/}(\|\mathbf{Y}_{k}\|\geq t)\leq c_{1}e^{-c_{2}t/\mu},$
(2.2)

for $t\geq c_{3}\mu$, where we set $\mu=d/\beta.$ At this point we would like
to apply Matrix Chernoff to the $\mathbf{Y}_{k}$’s. Unfortunately we cannot
since the $\mathbf{Y}_{k}$’s are not bounded a.s. Therefore we let
$\widetilde{\mathbf{Y}}_{k}$ be the random variable obtained by conditioning
$\mathbf{Y}_{k}$ on the event that $\|\mathbf{Y}_{k}\|\leq M.$ We will show
below that taking $M=c_{5}\mu\log\frac{\mu}{\varepsilon}$ gives
$\|\mathbb{E}\mathbf{Y}_{k}-\mathbb{E}\widetilde{\mathbf{Y}}_{k}\|\leq\varepsilon.$
As a consequence this gives

$(1-\varepsilon)\mathbf{I}\leq\mathbb{E}(\widetilde{\mathbf{Y}}_{k})\leq(1+\varepsilon)\mathbf{I}.$

Given this choice of $M$, we apply the Matrix Chernoff bound to the
$\widetilde{\mathbf{Y}}_{k}$’s, which now satisfy
$\|\widetilde{\mathbf{Y}}_{k}\|\leq M$. Setting
$\overline{\widetilde{Y}}=\frac{1}{N}\sum_{i=1}^{N}\widetilde{Y}_{k}$, and
plugging into Theorem 2.21 gives

$\operatorname{Pr\/}(\|\overline{\widetilde{\mathbf{Y}}}-\mathbf{I}\|\geq
3\varepsilon)\leq d\exp\lparen-c\varepsilon^{2}\frac{N}{M}\rparen,$

which is bounded by $0.05$ for
$N\geq\Omega\lparen\frac{d}{\beta\varepsilon^{2}}\log\lparen\frac{d}{\beta\varepsilon}\rparen\log
d\rparen.$ (We replace $\varepsilon$ with $\varepsilon/3$ to recover the
statement in the theorem.) Possibly by adjusting the constant in the
definition of $M$, we can arrange so that with probability at least $0.95$,
all of the $N$ samples $\mathbf{Y}_{k}$ are such that $\|\mathbf{Y}_{k}\|\leq
M$ (this follows from the exponential tail bound on the $\mathbf{Y}_{k}$’s),
and hence indistinguishable from the $\widetilde{\mathbf{Y}}_{k}$’s. So with
probability at least $0.9$ the $\mathbf{Y}_{k}$’s enjoy the same concentration
bound as the $\widetilde{\mathbf{Y}}_{k}$’s above, which then implies the
conclusion of the theorem.

Finally we conclude the argument by showing that
$\mathbb{E}\widetilde{\mathbf{Y}_{k}}$ is approximately
$\mathbb{E}\mathbf{Y}_{k}.$ To simplify notation, let $\mathbf{Y}$ and
$\widetilde{\mathbf{Y}}$ be distributed as $\mathbf{Y}_{k}$ and
$\widetilde{\mathbf{Y}}_{k}$ respectively. We write

$\mathbb{E}Y=\mathbb{E}(\widetilde{\mathbf{Y}})\operatorname{Pr\/}(\|\mathbf{Y}\|\leq
M)+\mathbb{E}(\mathbf{Y}|\|\mathbf{Y}\|\geq
M)\operatorname{Pr\/}(\|\mathbf{Y}\|\geq M).$

Thus we have

$\displaystyle\|\mathbb{E}\mathbf{Y}-\mathbb{E}\widetilde{\mathbf{Y}}\|$
$\displaystyle\leq(1-\operatorname{Pr\/}(\|\mathbf{Y}\|\leq
M))\|\mathbb{E}\widetilde{\mathbf{Y}}\|+\operatorname{Pr\/}(\|\mathbf{Y}\|>M)\|\mathbb{E}(\mathbf{Y}|\|\mathbf{Y}\|\geq
M)\|$
$\displaystyle=\operatorname{Pr\/}(\|\mathbf{Y}\|>M)\|\mathbb{E}\widetilde{\mathbf{Y}}\|+\operatorname{Pr\/}(\|\mathbf{Y}\|>M)\|\mathbb{E}(\mathbf{Y}|\|\mathbf{Y}\|\geq
M)\|$
$\displaystyle\leq\operatorname{Pr\/}(\|\mathbf{Y}\|>M)\mathbb{E}\|\widetilde{\mathbf{Y}}\|+\operatorname{Pr\/}(\|\mathbf{Y}\|>M)\mathbb{E}(\|\mathbf{Y}\||\|\mathbf{Y}\|\geq
M)$
$\displaystyle\leq\operatorname{Pr\/}(\|\mathbf{Y}\|>M)\mathbb{E}\|\mathbf{Y}\|+\operatorname{Pr\/}(\|\mathbf{Y}\|>M)\mathbb{E}(\|\mathbf{Y}\||\|\mathbf{Y}\|\geq
M),$

where in the last step we observed that
$\mathbb{E}\|\widetilde{\mathbf{Y}}\|\leq\mathbb{E}\|\mathbf{Y}\|$. We bound
each of the relevant terms.

We will take $M\geq c_{3}\mu.$ As shown above,
$\mathbb{E}\|\mathbf{Y}\|\leq\mu.$ By Equation 2.2,
$\operatorname{Pr\/}(\|\mathbf{Y}\|>M)\leq c_{1}e^{-c_{2}M/\mu}$. To handle
the last term,

$\mathbb{E}\left(\|\mathbf{Y}\|\bigg{|}\|\mathbf{Y}\|\geq
M\right)=\int_{0}^{\infty}\operatorname{Pr\/}\left(\|\mathbf{Y}\|\geq
t|\|\mathbf{Y}\|\geq
M\right)\,dt=M+\int_{M}^{\infty}\operatorname{Pr\/}\left(\|\mathbf{Y}\|\geq
t|\|\mathbf{Y}\|\geq M\right)\,dt.$

So

$\displaystyle\operatorname{Pr\/}(\|\mathbf{Y}\|>M)\mathbb{E}(\|\mathbf{Y}\||\|\mathbf{Y}\|\geq
M)$
$\displaystyle=\operatorname{Pr\/}(\|\mathbf{Y}\|>M)M+\int_{M}^{\infty}\operatorname{Pr\/}\left(\|\mathbf{Y}\|\geq
t\right)\,dt$ $\displaystyle\leq
c_{1}e^{-c_{2}M/\mu}M+\int_{M}^{\infty}c_{1}e^{-c_{2}t/\mu}\,dt$
$\displaystyle=c_{1}e^{-c_{2}M/\mu}M+c_{4}\mu e^{-c_{2}M/\mu}.$

Putting the pieces together gives

$\|\mathbb{E}\mathbf{Y}-\mathbb{E}\widetilde{\mathbf{Y}}\|\leq c_{1}\mu
e^{-c_{2}M/\mu}+c_{1}Me^{-c_{2}M/\mu}+c_{4}\mu e^{-c_{2}M/\mu},$

which is bounded by $\varepsilon$ for $M\geq
c_{5}\mu\log\left(\frac{\mu}{\varepsilon}\right).$ ∎

###  2.3 Low-Rank Matrix Approximation

In the preceding section, we showed how the coordinator can learn a subspace
embedding matrix $\mathbf{S}$ for $\mathbf{A}$. We now show how to utilize
this embedding to learn with efficient communication a projection giving a
good rank-$k$ approximation to $\mathbf{A}$.

See 1.7

###### Proof.

We first reduce the number of columns of $\mathbf{A}$ by right-multiplying by
a Rademacher random matrix $\mathbf{R}$ with $O(k/\varepsilon)$ columns. Note
that this computation can be carried out locally on each server; each server
simply computes $\mathbf{A}^{(i)}\mathbf{R}$, where $\mathbf{R}$ is known
using shared randomness. Since $\mathbf{R}$ is Rademacher, the bit complexity
of each $\mathbf{A}^{(i)}\mathbf{R}$ is at most $L+O(\log d).$ By [CW09,
Theorem 4.2], we have that

$\min_{\textrm{rank}(\mathbf{X})=k}\|\mathbf{A}\mathbf{R}\mathbf{X}-\mathbf{A}\|_{F}\leq(1+\varepsilon/3)\|\mathbf{A}_{k}-\mathbf{A}\|_{F},$

where $\mathbf{A}_{k}$ is defined in Section 1.1.1 as the best rank-$k$
approximation of $\mathbf{A}.$ We would like to find an approximate minimizer
$\mathbf{X}$, since then $\mathbf{A}\mathbf{R}\mathbf{X}$ is a rank-$k$
approximation of $\mathbf{A}.$ To do this efficiently, we use our protocol
from Section 2.1666One could also apply our protocol from Section 2.2 (with
minor modifications when we apply the black-boxed Lemma 32 below as the
resulting matrix is not a true leverage score sampling matrix, since we sample
sketches of the rows from each block). to construct a row-sampling matrix
$\mathbf{S}$ from the leverage score distribution of $\mathbf{A}\mathbf{R}$
(up to constant factors on the probabilities), with
$\widetilde{O}(k/\varepsilon)$ rows. Running our protocol requires
$\widetilde{O}\left((\frac{sk}{\varepsilon}+\frac{k^{2}}{\varepsilon^{2}})\cdot
L\right)$ communication (note that the rows of $\mathbf{A}\mathbf{R}$ have
dimension $O(k/\varepsilon)$). Then Lemma 32 and Theorem 36 of [CW17] together
imply that

$\min_{\operatorname{\textrm{rank}}(\mathbf{X})=k}\|\mathbf{S}\mathbf{A}\mathbf{R}\mathbf{X}-\mathbf{S}\mathbf{A}\|_{F}=\min_{\operatorname{\textrm{rank}}(\mathbf{X})=k}\|\mathbf{S}(\mathbf{A}\mathbf{R}\mathbf{X}-\mathbf{A})\|_{F}\leq(1+\varepsilon/3)\min_{\operatorname{\textrm{rank}}(\mathbf{X})=k}\|\mathbf{A}\mathbf{R}\mathbf{X}-\mathbf{A}\|_{F}.$

The coordinator can learn $\mathbf{S}\mathbf{A}\mathbf{R}$ and
$\mathbf{S}\mathbf{A}$ using total communication
$\widetilde{O}(kdL/\varepsilon^{2}).$ This allows the coordinator to compute
$\widehat{\mathbf{X}}$ with

$\|\mathbf{A}\mathbf{R}\widehat{\mathbf{X}}-\mathbf{A}\|_{F}\leq(1+\varepsilon/3)\min_{\operatorname{\textrm{rank}}(\mathbf{X})=k}\|\mathbf{A}\mathbf{R}\mathbf{X}-\mathbf{A}\|_{F}\leq(1+\varepsilon)\|\mathbf{A}_{k}-\mathbf{A}\|_{F}.$

Let $\Pi$ be the orthogonal projection onto the row space of
$\mathbf{R}\widehat{\mathbf{X}}$ which has dimension at most $k$ by
construction. For a given row $\mathbf{a}_{i}$ of $\mathbf{A}$,
$\mathbf{a}_{i}^{\top}\mathbf{R}\widehat{\mathbf{X}}$ is in the row space of
$\mathbf{R}\widehat{\mathbf{X}}$. So
$\|\mathbf{a}_{i}^{\top}\Pi-\mathbf{a}_{i}^{\top}\|^{2}\leq\|\mathbf{a}_{i}^{\top}\mathbf{R}\widehat{\mathbf{X}}-\mathbf{a}_{i}^{\top}\|^{2}$,
and hence
$\|\mathbf{A}\Pi-\mathbf{A}\|_{F}\leq\|\mathbf{A}\mathbf{R}\widehat{\mathbf{X}}-\mathbf{A}\|_{F}\leq(1+\varepsilon)\|\mathbf{A}_{k}-\mathbf{A}\|_{F}.$

∎

##  3 High-Accuracy Linear Regression in the Coordinator Model

In this section, we present our result for the communication complexity of
computing high-accuracy solutions to linear regression problems as in the
setup of Problem 1.3. The primary result of this section is the following.

See 1.8

Our framework for achieving the results in Theorem 1.8 builds upon
Richardson’s iteration with preconditioning. This circumvents the need to send
to the coordinator the matrix $\mathbf{A}^{\top}\mathbf{A}$ — the approach in
the previous best result [VWW20] for this problem, which incurs a
communication cost $\Omega(sd^{2}L)$777To see this, we observe that each of
$s$ servers computes and sends to the coordinator
$\langle\mathbf{A}^{(i)},\mathbf{A}^{(i)}\rangle\in\mathbb{R}^{d\times d}$,
each of which takes $O(d^{2})$ bits., thus exceeding our targeted budget. To
ensure convergence with this Richardson-style iteration, the key conceptual
idea is to use a matrix $\mathbf{M}$ that spectrally approximates
$\mathbf{A}^{\top}\mathbf{A}$ as a preconditioner. We construct $\mathbf{M}$
via sampling with respect to overestimates of the leverage scores of
$\mathbf{A}$, which are computed using an iterative process we build upon the
refinement sampling framework of [CLMMPS15]. The main novelty of our algorithm
and analysis are in careful roundings and bit complexity analysis to guarantee
convergence while ensuring a small number of bits are communicated. In Section
3.1, we discuss the main components of our algorithm and analysis.

###  3.1 An Overview of Our Algorithm and Analysis

Our main algorithm for solving the linear regression problem to high-accuracy
is LABEL:alg:lin-reg-coordinator-poly-cond. It has three main components that
we explain next.

Input. A matrix $\mathbf{A}:=[\mathbf{A}^{(i)}]\in\mathbb{R}^{n\times d}$ and
vector $\mathbf{b}:=[\mathbf{b}^{(i)}]\in\mathbb{R}^{n}$, where the
$i^{\mathrm{th}}$ machine stores matrix
$\mathbf{A}^{(i)}\in\mathbb{R}^{n_{i}\times d}$ and vector
$\mathbf{b}^{(i)}\in\mathbb{R}^{n_{i}}$; accuracy parameter $0<\varepsilon<1$;
probability parameter $c$.

Output. Vector $\widehat{\mathbf{x}}\in\mathbb{R}^{d}$ such that

$\|\mathbf{A}\widehat{\mathbf{x}}-\mathbf{b}\|_{2}\leq\varepsilon\cdot\|\mathbf{A}(\mathbf{A}^{\top}\mathbf{A})^{-1}\mathbf{A}^{\top}\mathbf{b}\|_{2}+\min_{\mathbf{x}\in\mathbb{R}^{d}}\|\mathbf{A}\mathbf{x}-\mathbf{b}\|_{2}$

  1. 1

Run LABEL:alg:levscoresRefinementSampling with inputs $\mathbf{A}$ and $c$.
This returns $\widehat{\tau}\in\mathbb{R}^{n}_{\geq 0}$, which satisfies
$\widehat{\tau}\geq\tau(\mathbf{A})$ and $\|\widehat{\tau}\|_{1}\leq 9d$. We
store the coordinates of $\widehat{\tau}$ on the corresponding machines.

  2. 2

Using $\widehat{\tau}$ from Step 1, each machine $i\in[s]$ forms a diagonal
sampling matrix $\mathbb{R}^{n\times
n}\ni\mathbf{S}^{(i)}=\texttt{Sample}(\widehat{\tau},100,c)$ (cf. Definition
3.3) and sends to the coordinator the following objects: the nonzero entries
of $\mathbf{S}^{(i)}$ and the rows of $\mathbf{A}^{(i)}$ corresponding to the
nonzero entries of $\mathbf{S}^{(i)}$.

  3. 3

The coordinator forms the matrix
$\widetilde{\mathbf{A}}:=\frac{1}{\sqrt{1.1}}[\widetilde{\mathbf{S}}^{(i)}\mathbf{A}^{(i)}]$,
where $\widetilde{\mathbf{S}}^{(i)}$ is the matrix left after removing the
zero rows of $\mathbf{S}^{(i)}$. Set
$\mathbf{M}=\widetilde{\mathbf{A}}^{\top}\widetilde{\mathbf{A}}$.

  4. 4

Implement the protocol in LABEL:alg:richardson with $\mathbf{A}$,
$\mathbf{b}$, and $\mathbf{M}$ as input.

Algorithm 6: Protocol for linear regression in the coordinator setting.

alg]alg:lin-reg-coordinator-poly-cond

##### Leverage score computation.

The first step of LABEL:alg:lin-reg-coordinator-poly-cond, as shown in Step 1,
is to compute $\widehat{\tau}$, a vector of sufficiently accurate
overestimates of $\tau(\mathbf{A})$, the true leverage scores of $\mathbf{A}$.
We do this using LABEL:alg:levscoresRefinementSampling, which iteratively
refines our initial crude overestimates. In each iteration of
LABEL:alg:levscoresRefinementSampling, we sample a matrix
$\widetilde{\mathbf{A}}$, a spectral approximation of $\mathbf{A}$, and use
$(\widetilde{\mathbf{A}}^{\top}\widetilde{\mathbf{A}})^{\dagger}$ to compute a
new $\widehat{\tau}$ with improved accuracy. This reduces
$\|\widehat{\tau}\|_{1}$ by a constant factor in each iteration, and therefore
LABEL:alg:levscoresRefinementSampling runs for only $O(\log n)$ iterations
before $\|\widehat{\tau}\|_{1}\leq O(d)$, the required accuracy of
$\widehat{\tau}$. A detailed discussion is presented in Section 3.2.

Throughout LABEL:alg:levscoresRefinementSampling, the machines communicate to
the coordinator only the $\widetilde{O}(d)$ rows they sample locally according
to their leverage score overestimates. The coordinator forms
$\widetilde{\mathbf{A}}^{\top}\widetilde{\mathbf{A}}$, a constant-factor
spectral approximation to $\mathbf{A}^{\top}\mathbf{A}$ and sends to the
machines its sketch formed by a combination of rounding and a Johnson-
Lindenstrauss random projection. Then the machines use this sketch to update
their overestimates to $\widehat{\tau}^{\widetilde{\mathbf{A}}}$, an
approximation of $\tau^{\widetilde{\mathbf{A}}}(\mathbf{A})$, the generalized
leverage scores of $\mathbf{A}$ computed with respect to
$\widetilde{\mathbf{A}}$. While this approach is similar to that of
[CLMMPS15], our target communication complexity necessitates additional
approximations, e.g. via roundings, of the objects we communicate between the
coordinator and the machines (cf. Step 5 and Step 6 of
LABEL:alg:levscoresRefinementSampling), whereas [CLMMPS15] is analyzed in
exact arithmetic. We elaborate this in Proposition 3.6 and the overall
communication complexity bound for computing the leverage score overestimates
in Lemma 3.2.

##### Sampling a spectral approximation.

After computing a sufficiently accurate $\widehat{\tau}$, in the next step
(Step 2 of LABEL:alg:lin-reg-coordinator-poly-cond), each machine locally
samples a set of rows according to its leverage score overestimates and
communicates these to the coordinator. The coordinator uses these rows to form
$\mathbf{M}$, the final spectral approximation of
$\mathbf{A}^{\top}\mathbf{A}$, as stated in Step 3 of LABEL:alg:lin-reg-
coordinator-poly-cond.

##### Richardson-type iteration.

Finally, in Step 4 of LABEL:alg:lin-reg-coordinator-poly-cond,
$\mathbf{M}^{-1}$ is used as a preconditioner in a Richardson-type algorithm,
displayed in LABEL:alg:richardson. In each iteration of LABEL:alg:richardson,
a residual vector is computed on each machine and communicated to the
coordinator. In Step 2 and Step 3 of LABEL:alg:richardson, the coordinator
linearly combines these residual vectors from all the machines and
communicates back to the machines a vector carefully designed to optimize
communication complexity. The output of LABEL:alg:richardson is the solution
to the given regression problem at the specified accuracy, with a
communication cost of $\widetilde{O}(sd(L+\log\kappa)\log(\varepsilon^{-1}))$
(cf. Lemma 3.10). A detailed discussion is presented in Section 3.3

######  Remark 3.1.

The primary theme of all the components of LABEL:alg:lin-reg-coordinator-poly-
cond is the simple idea that we communicate only the bits that are necessary
for the convergence of our algorithms. This idea is perhaps most apparent in
the last component (LABEL:alg:richardson). To this end, we carefully round
intermediate vectors and matrices (i.e., discard bits with low place values)
and reuse some bits (e.g., by not communicating bits with high place values).
The latter is possible essentially because in an iterative algorithm, when the
solution is converging, the bits with high place values are the same from one
iteration to the next.

##### Discussion of our analysis.

We conclude this overview by reiterating that although our approach is
conceptually fairly simple, the overall algorithm and analysis are involved
due to careful bit modifications in intermediate steps such as the JL
projection, matrix inversions, and Richardson’s iteration to reduce the
communication complexity while ensuring convergence. The careful bit
modification and analysis in turn plays a crucial role in yielding
improvements compared to [VWW20], particularly for matrices with a condition
number of $e^{o(dL)}$, as summarized in Section 1.1. In subsequent sections,
we explain each of the components of our approach. We first discuss the
computation of leverage score overestimates in Section 3.2 and our iterative
preconditioning approach for solving the linear regression problem in Section
3.3. We finally prove the main result of this section in Section 3.4.

###  3.2 Leverage Score Overestimates

We now focus on the point-to-point communication complexity of computing a
vector of overestimates for leverage scores of the matrix. We emphasize that
although [VWW20] considers computing the leverage scores in the blackboard
model, this is significantly more complicated in the coordinator (point-to-
point) setting due to bit complexity issues involving inverse of matrices. The
main result of this section is Lemma 3.2.

######  Lemma 3.2.

Given the linear regression setting of Problem 1.3 with matrix
$\mathbf{A}=[\mathbf{A}^{(i)}]\in\mathbb{R}^{n\times d}$, and $n\geq 5$, there
is a randomized algorithm that, with

$\widetilde{O}(d^{2}L+sd(L+\log(\kappa)))\text{ bits of communication}$

and, with high probability, computes a vector
$\widehat{\tau}\in\mathbb{R}^{n}$ such that $\|\widehat{\tau}\|_{1}\leq 9d$
and $\widehat{\tau}_{i}\geq\tau_{i}(\mathbf{A})$, for all $i\in[n]$, where
$\tau_{i}(\mathbf{A})$ is the $i^{\mathrm{th}}$ leverage score of the matrix
$\mathbf{A}$. Each entry of $\widehat{\tau}$ is stored on the machine that
contains the corresponding row.

We achieve the results in Lemma 3.2 via LABEL:alg:levscoresRefinementSampling.
As described in Section 3.1, this algorithm is then used as a black-box in our
main algorithm (LABEL:alg:lin-reg-coordinator-poly-cond) as a pre-cursor to
construct a spectral approximation of $\mathbf{A}$.
LABEL:alg:levscoresRefinementSampling is based on the Refinement Sampling
algorithm of [CLMMPS15] to approximately compute leverage scores of a matrix,
with appropriate modifications for the coordinator setting. We provide an
overview of the distributed version of Refinement Sampling in Section 3.2.1
and prove Lemma 3.2 in Section 3.2.4. While the results of [CLMMPS15] are
provided under the exact arithmetic model, whereas we need to employ careful
bit complexity analysis under fixed-point arithmetic.

####  3.2.1 An Overview of Refinement Sampling

The Refinement Sampling algorithm returns a small-sized spectral approximation
to an input matrix $\mathbf{A}$ by iteratively refining (hence the name)
$\widehat{\tau}$, the leverage score overestimates of $\mathbf{A}$; when
$\|\widehat{\tau}\|_{1}\leq O(d)$, it is accurate enough for use (with
appropriate scaling) in sampling a matrix with the desired approximation
guarantee, leading to termination of the algorithm.

The algorithm starts by setting $\widehat{\tau}=\mathbf{1}$, the vector of all
ones. In each iteration, the algorithm uses $\widehat{\tau}$ to construct a
matrix $\widetilde{\mathbf{A}}\in\mathbb{R}^{\widetilde{n}\times d}$ composed
of a subset of $\widetilde{n}$ (rescaled) rows of $\mathbf{A}$, where each row
$i$ is sampled independently with probability
$\mathbf{p}_{i}\propto\alpha\widehat{\tau}_{i}$ for some sampling rate
$\alpha>0$ (cf. Definition 3.3 for the precise row sampling). The vector
$\widehat{\tau}$ is then updated to
$\tau^{\widetilde{\mathbf{A}}}(\mathbf{A})$, the generalized leverage scores
of $\mathbf{A}$ with respect to $\widetilde{\mathbf{A}}$ (cf. Definition
1.23), and this process continues iteratively.

We pick $\alpha$ such that the number of rows of $\widetilde{\mathbf{A}}$,
which is proportional to $\alpha\cdot\|\widehat{\tau}\|_{1}$, is
$\widetilde{O}(d)$ with high probability. This ensures (see Step 3) that the
number of rows communicated between the coordinator and machines is only
$\widetilde{O}(d)$, incurring a bit complexity of $\widetilde{O}(d^{2}L)$.
Therefore, iteratively reducing $\|\widehat{\tau}\|_{1}$ by a constant factor
enables a corresponding increase in $\alpha$, which eventually yields the
desired spectral approximation (Lemma 3.4 and Lemma 3.5, first stated in
[CLMMPS15]), all the while maintaining a row size of $\widetilde{O}(d)$ for
the matrix used to compute generalized leverage scores. We state these formal
guarantees next since they are used in our proof of correctness of
LABEL:alg:levscoresRefinementSampling.

######  Definition 3.3 (Sampling Function [CLMMPS15]).

Given a vector $\mathbf{u}\in\mathbb{R}^{n}_{\geq 0}$, a parameter $\alpha>0$,
and a positive constant $c$, we define vector
$\mathbf{p}\in\mathbb{R}^{n}_{\geq 0}$ as $\mathbf{p}_{i}=\min(1,\alpha c\log
d\cdot\mathbf{u}_{i})$. We define the function
$\texttt{Sample}(\mathbf{u},\alpha,c)$ to be one which returns a random
diagonal $n\times n$ matrix $\mathbf{S}$ with independently chosen entries:

$\mathbf{S}_{ii}=\left\\{\begin{array}[]{ll}\frac{1}{\sqrt{\mathbf{p}_{i}}}&\mbox{}\textrm{with
probability }\mathbf{p}_{i}\\\ 0&\mbox{}\textrm{otherwise}\end{array}\right..$

The spectral approximation guarantees of this sampling approach are formalized
in the following lemma.

######  Lemma 3.4 (Spectral Approximation via Leverage Score Sampling; Lemma 4
of [CLMMPS15]).

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$, a sampling rate
$\alpha>1$, and a fixed constant $c>0$. Let $\mathbf{u}\in\mathbb{R}^{n}_{\geq
0}$ be a vector of leverage score overestimates, that is,

$\mathbf{u}\geq\tau(\mathbf{A}),\,\,\text{ and
}\mathbf{S}:=\texttt{Sample}(\mathbf{u},\alpha,c)$

as in Definition 3.3. Then, with probability at least $1-d^{-c/3}-(3/4)^{d}$,
the following results hold:

$\textbf{nnz}({\mathbf{S}})=2c\alpha\|\mathbf{u}\|_{1}\log d\text{ and
}\frac{1}{\sqrt{1+\alpha^{-1/2}}}\mathbf{SA}\approx_{\left(\frac{1+\alpha^{-1/2}}{1-\alpha^{-1/2}}\right)}\mathbf{A}.$

The following lemma formalizes how $\|\mathbf{u}\|_{1}$ shrinks when we use a
spectral approximation $\widetilde{\mathbf{A}}$ obtained by the sampling
approach of Definition 3.3 (cf. Lemma 3.4) to update $\widehat{\tau}$ (the
leverage score estimates) to $\tau^{\widetilde{\mathbf{A}}}(\mathbf{A})$ (the
generalized leverage scores of $\mathbf{A}$ with respect to
$\widetilde{\mathbf{A}}$).

######  Lemma 3.5 (Leverage Score Estimate Update via Undersampling; Theorem 3
of [CLMMPS15]).

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ and an undersampling
parameter $\alpha\in(0,1]$, let $\mathbf{u}\in\mathbb{R}^{n}_{\geq 0}$ and
$\mathbf{u}^{(\textrm{new})}\in\mathbb{R}^{n}_{\geq 0}$ be vectors such that:

$\mathbf{u}\geq\tau(\mathbf{A}),\,\,\mathbf{S}:=\sqrt{\frac{3\alpha}{4}}\texttt{Sample}(\mathbf{u},9\alpha,c),\text{
and
}\mathbf{u}_{i}^{(\textrm{new})}:=\min\\{\tau_{i}^{\mathbf{S}\mathbf{A}}(\mathbf{A}),\mathbf{u}_{i}\\}\text{
for all }i\in[n].$

Then, with probability at least $1-d^{-c/3}-(3/4)^{d}$,
$\mathbf{u}_{i}^{(\textrm{new})}$ is a leverage score overestimate, i.e.,
$\mathbf{u}_{i}^{(\textrm{new})}\geq\tau_{i}(\mathbf{A})$. Furthermore,

$\|\mathbf{u}^{(\textrm{new})}\|_{1}\leq 3d/\alpha\text{ and
}\textbf{nnz}(\mathbf{S})=O(\alpha\cdot\|\mathbf{u}\|_{1}\cdot\log d).$

####  3.2.2 Invariance of Leverage Score Overestimates:
$\widehat{\tau}\geq\tau(\mathbf{A})$

The bulk of communication in LABEL:alg:levscoresRefinementSampling happens in
Step 3 and Step 5. In Step 3, we send only the rows of the original matrix and
the associated vector of probabilities. We construct the probabilities to be
powers of two so that they can be communicated with a small number of bits.
However, in Step 5, we need to send the product of the inverse of a matrix
with other matrices. To do this, we need to round the product and communicate
the rounded version. The following two technical lemmas (Proposition 3.6 and
Proposition 3.7) help us bound the error in computing the generalized leverage
scores arising from this rounding process.

Input. A matrix $\mathbf{A}:=[\mathbf{A}^{(i)}]\in\mathbb{R}^{n\times d}$,
where $n=\sum_{i=1}^{s}n_{i}$ and the $i^{\mathrm{th}}$ machine stores matrix
$\mathbf{A}^{(i)}\in\mathbb{R}^{n_{i}\times d}$; probability parameter $c$.

Output. Vector $\widehat{\tau}$ that satisfies
$\widehat{\tau}\geq\tau(\mathbf{A})$ and $\|\widehat{\tau}\|_{1}\leq 9d$. The
$i^{\mathrm{th}}$ machine stores $\widehat{\tau}^{(i)}$, the set of
coordinates of $\widehat{\tau}$ corresponding to the rows of $\mathbf{A}$
stored on that machine.

Initialize. Set the total number of iterations
$T=\left\lceil\log_{2}(n/d)\right\rceil$ and $r=10^{6}\cdot\log n$. Set the
leverage scores estimate vector $\widehat{\tau}=\mathbf{1}$. Set
$\lambda=\frac{1}{100\kappa^{2}}$, where $\kappa$ is the condition number of
$\mathbf{A}$.

For iterations $\ell=1,2,\dots,T$:

  1. 1

Each machine $i$ computes $t_{i}=\|\widehat{\tau}^{(i)}\|_{1}$, the sum of its
leverage score overestimates, and sends it to the coordinator.

  2. 2

The coordinator computes
$t\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sum_{i=1}^{s}t_{i}$ and sends it
to all the machines.

  3. 3

All machines set $\widehat{\alpha}$ to the smallest power of half that is at
least $\frac{25d\log d}{t}$. Set $\alpha=\min\\{1,\frac{\widehat{\alpha}}{\log
d}\\}$. Each machine $i$ forms an $n_{i}$-by-$n_{i}$ random diagonal matrix
$\mathbf{S}^{(i)}=\texttt{Sample}(1.01\widehat{\tau}^{(i)},9\alpha,c)$ as per
Definition 3.3 and sends to the coordinator those rows of $\mathbf{A}^{(i)}$
and coordinates of $\widehat{\tau}^{(i)}$ that correspond to the nonzero
entries of $\mathbf{S}^{(i)}$.

  4. 4

Using the received entries of $\widehat{\tau}^{(i)}$ for all $i\in[s]$, the
coordinator computes $\widetilde{\mathbf{S}}^{(i)}$, the non-zero rows of
$\mathbf{S}^{(i)}$ as per Definition 3.3. It then constructs the matrices
$\widetilde{\mathbf{A}}^{(i)}=\sqrt{\frac{3\alpha}{4}}\widetilde{\mathbf{S}}^{(i)}\mathbf{A}^{(i)}$
and sets
$\widetilde{\mathbf{A}}=[\widetilde{\mathbf{A}}^{(i)}]\in\mathbb{R}^{\widetilde{n}\times
d}$ and
$\mathbf{B}=[\widetilde{\mathbf{A}};\sqrt{\lambda}\cdot\mathbf{I}]\in\mathbb{R}^{(\widetilde{n}+d)\times
d}$.

  5. 5

The coordinator samples a JL sketching matrix
$\mathbf{G}\in\\{-1,+1\\}^{r\times(\widetilde{n}+d)}$ and uses $\mathbf{B}$
from the previous step to compute
$\widehat{\mathbf{G}}:=\frac{\sqrt{1.01}}{0.99\sqrt{r}}\mathbf{G}\mathbf{B}(\mathbf{B}^{\top}\mathbf{B})^{-1}\in\mathbb{R}^{r\times
d}$. It then generates $\mathbf{J}\in\mathbb{R}^{r\times d}$ by rounding the
entries of $\widehat{\mathbf{G}}$ so that for all $i_{1},i_{2}$,
$|\widehat{g}_{i_{1}i_{2}}-j_{i_{1}i_{2}}|<\frac{1}{10^{4}\cdot
n^{2}d\sqrt{r}\cdot 2^{L}}$ and sends $\mathbf{J}$ to all the machines.

  6. 6

The coordinator computes an integer basis for the kernel of
$[\widetilde{\mathbf{S}}^{(i)}\mathbf{A}^{(i)}]$ (e.g., by Gaussian
elimination). It computes a linear combination $\mathbf{v}$ of the kernel
basis by picking independent and uniformly random coefficients for in
$[-2^{dL},2^{dL}]$. It picks $z=\left\lceil 100\log(n)\right\rceil$ random
prime numbers $y_{1},\ldots,y_{z}$ less than $(dL)^{2}$. For each prime number
$y_{k}$, it sets $\mathbf{v}^{(k)}$ to be $\mathbf{v}$ modulo $y_{k}$. It then
sends all $\mathbf{v}^{(k)}$’s and $y_{k}$’s to all the machines.

  7. 7

If for all $k\in[z]$, $\mathbf{a}_{j}^{\top}\mathbf{v}^{(k)}=0$ modulo
$y_{k}$, then we set $\widehat{\tau}^{\widetilde{\mathbf{A}}}_{j}=\infty$ (on
the machine holding row $\mathbf{a}_{j}$). Otherwise, we set
$\widehat{\tau}^{\widetilde{\mathbf{A}}}_{j}=\|\mathbf{J}\mathbf{a}_{j}\|_{2}^{2}$.
Then we update $\widehat{\tau}_{j}^{\textrm{new}}$ to the smallest power of
two that is at least
$\max\left\\{\min\left\\{\widehat{\tau}_{j},\widehat{\tau}^{\widetilde{\mathbf{A}}}_{j}\right\\},\frac{1}{2n^{2}}\right\\}$.

  8. 8

Each machine updates its overestimates with
$\widehat{\tau}_{j}=1.01\cdot\widehat{\tau}_{j}^{\textrm{new}}$.

Algorithm 7: Protocol for computing leverage score overestimates in the
coordinator setting

alg]alg:levscoresRefinementSampling

######  Proposition 3.6.

Let $0<\widehat{\varepsilon}<1$ be an accuracy parameter, and let
$\mathbf{u},\widetilde{\mathbf{u}}\in\mathbb{R}^{d}$ be vectors satisfying,
for all $i\in[d]$, that $|u_{i}-\widetilde{u}_{i}|\leq\widehat{\varepsilon}$.
Then for any $\mathbf{v}\in\mathbb{R}^{d}$, we have
$|{\widetilde{\mathbf{u}}^{\top}\mathbf{v}-\mathbf{u}^{\top}\mathbf{v}}|\leq\sqrt{d}\cdot\widehat{\varepsilon}\|\mathbf{v}\|_{2}$.
Moreover if $\mathbf{B},\widetilde{\mathbf{B}}\in\mathbb{R}^{m\times d}$ such
that for all $i\in[m],j\in[d]$,
$|b_{ij}-\widetilde{b}_{ij}|\leq\widehat{\varepsilon}$, then

$|\|\mathbf{B}\mathbf{v}\|_{2}-\|\widetilde{\mathbf{B}}\mathbf{v}\|_{2}|\leq\sqrt{md}\cdot\widehat{\varepsilon}\|\mathbf{v}\|_{2}.$

######  Proposition 3.7 (Approximating Ridge Leverage Scores with Leverage
Scores).

Let $\tau(\mathbf{A})$ and $\tau^{\lambda}(\mathbf{A})$ be the vector of
leverage scores and $\lambda$-ridge leverage scores (see Definition 1.25),
respectively, of $\mathbf{A}\in\mathbb{R}^{n\times d}$, $n\geq d$. Let
$\lambda\geq 0$ and $\sigma_{\min}$ be the smallest _nonzero_ singular value
of $\mathbf{A}$. Then

$\frac{\sigma_{\min}^{2}+\lambda}{\sigma_{\min}^{2}}\cdot\tau^{\lambda}_{j}(\mathbf{A})\geq\tau_{j}(\mathbf{A})\geq\tau^{\lambda}_{j}(\mathbf{A})\text{
for all $j\in[n]$. }$

Equipped with these technical results, we prove the different components of
Lemma 3.2, starting with the invariant of
LABEL:alg:levscoresRefinementSampling that the leverage score overestimates
vector $\widehat{\tau}$ indeed always remains larger than
$\tau^{\lambda}(\mathbf{A})$, the $\lambda$-ridge leverage scores vector of
$\mathbf{A}$. As we see in Lemma 3.2, this implies
$\widehat{\tau}\geq\tau(\mathbf{A})$.

######  Lemma 3.8 (Invariance of Leverage Score Overestimates).

Assume that at the start of each iteration $\ell\in[T]$ of
LABEL:alg:levscoresRefinementSampling, we have
$\widehat{\tau}\geq\tau^{\lambda}(\mathbf{A})$. Then, at the end of this
iteration, we have
$\widehat{\tau}^{\textrm{(new)}}\geq\tau^{\lambda}(\mathbf{A})$.

###### Proof of Lemma 3.8.

The random vector $\mathbf{v}$ computed in Step 6 of
LABEL:alg:levscoresRefinementSampling satisfies
$\mathbf{v}\in\mathcal{N}(\widetilde{\mathbf{A}})$ (we remove the
$\sqrt{\frac{3\alpha}{4}}$ coefficient only to be sure that the matrix is an
integer matrix and we can find an integer basis for the kernel). Therefore if
a vector $\mathbf{w}\in\mathcal{N}(\widetilde{\mathbf{A}})^{\perp}$, then
$\mathbf{w}^{\top}\mathbf{v}=0$ as well as
$\mathbf{w}^{\top}\mathbf{v}^{(k)}=0$. Otherwise, if
$\mathbf{w}\notin\mathcal{N}(\widetilde{\mathbf{A}})^{\perp}$, then with high
probability, $\mathbf{w}^{\top}\mathbf{v}\neq 0$. Since both $\mathbf{v}$ and
$\mathbf{a}_{j}$ are integer vectors, we can look at
$\mathbf{a}_{j}^{\top}\mathbf{v}$ modulo $y_{k}$. Since we are using Gaussian
elimination to compute the integer basis for the kernel and since the
$\mathbf{A}$ is an integer matrix with bit complexity $L$,
$|\mathbf{a}_{j}^{\top}\mathbf{v}|\leq\operatorname{poly}(n)\cdot 2^{dL}$.
Therefore the number of prime factors of $|\mathbf{a}_{j}^{\top}\mathbf{v}|$
is $O(dL\log n)$. Therefore if we select a random prime number less than
$(dL)^{2}$, then with a large probability, we pick a prime $y_{k}$ that is not
a factor of $|\mathbf{a}_{j}^{\top}\mathbf{v}|\leq\operatorname{poly}(n)\cdot
2^{dL}$. Therefore if $|\mathbf{a}_{j}^{\top}\mathbf{v}|\neq 0$, then with
high probability, $|\mathbf{a}_{j}^{\top}\mathbf{v}|$ modulo $y_{k}$ is also
not zero. In Step 6 of LABEL:alg:levscoresRefinementSampling, we select
multiple random primes independently to boost this probability even further.
Then by taking union bound over all the rows of $\mathbf{A}$, the algorithm,
with high probability, can detect which ones are orthogonal to the kernel of
$\widetilde{\mathbf{A}}$. Thus if
$\mathbf{a}_{j}\notin\mathcal{N}(\widetilde{\mathbf{A}})^{\perp}$, then with
high probability, we set $\widehat{\tau}^{\widetilde{\mathbf{A}}}_{j}=\infty$,
and therefore (from Step 7), we have the desired inequality

$\widehat{\tau}_{j}^{(\textrm{new})}\geq\max\left\\{\min\left\\{\widehat{\tau}_{j},\widehat{\tau}^{\widetilde{\mathbf{A}}}_{j}\right\\},\frac{1}{2n^{2}}\right\\}\geq\widehat{\tau}_{j}\geq\tau^{\lambda}_{j}(\mathbf{A}),$

where the final inequality follows from the lemma’s assumption.

Next, consider the case in which
$\mathbf{a}_{j}\in\mathcal{N}(\widetilde{\mathbf{A}})^{\perp}$. In this case,
Step 7 of LABEL:alg:levscoresRefinementSampling sets
$\widehat{\tau}^{\widetilde{\mathbf{A}}}_{j}=\|\mathbf{J}\mathbf{a}_{j}\|_{2}^{2}$
for $\mathbf{J}$ as defined in Step 5 of
LABEL:alg:levscoresRefinementSampling. We now proceed to show
$\|\mathbf{J}\mathbf{a}_{j}\|_{2}^{2}\geq\tau^{\lambda}_{j}(\mathbf{A})$;
since the lemma assumes
$\widehat{\tau}_{j}\geq\tau^{\lambda}_{j}(\mathbf{A})$, we then have that

$\|\mathbf{J}\mathbf{a}_{j}\|_{2}^{2}\geq\tau^{\lambda}_{j}(\mathbf{A})\implies\widehat{\tau}_{j}^{\textrm{new}}\geq\max\left\\{\min\left\\{\widehat{\tau}_{j},\widehat{\tau}^{\widetilde{\mathbf{A}}}_{j}\right\\},\frac{1}{2n^{2}}\right\\}\geq\tau^{\lambda}_{j}(\mathbf{A}),$

which is the desired inequality. In the rest of this proof, we show
$\|\mathbf{J}\mathbf{a}_{j}\|_{2}^{2}\geq\tau^{\lambda}_{j}(\mathbf{A})$.

First note that
$\sqrt{\frac{3\alpha}{4}}\mathbf{S}^{(i)}:=\sqrt{\frac{3\alpha}{4}}\texttt{Sample}(1.01\cdot\widehat{\tau}^{(i)},9\alpha,c)$
(as set in Step 3 of LABEL:alg:levscoresRefinementSampling) has a distribution
equivalent to
$\widehat{\mathbf{S}}^{(i)}:=\sqrt{\frac{3}{4}}\texttt{Sample}(1.01\cdot\widehat{\tau}^{(i)},9,c)$
with some nonzero entries set to zero. This is because in Step 3 of
LABEL:alg:levscoresRefinementSampling, we set $\alpha\leq 1$, which implies
that $\texttt{Sample}(1.01\cdot\widehat{\tau}^{(i)},9\alpha,c)$ has a smaller
sampling rate compared to $\texttt{Sample}(1.01\cdot\widehat{\tau}^{(i)},9,c)$
for all entries, but the value of any entry selected in both matrices (i.e.,
nonzero entry) is the same in $\sqrt{\frac{3\alpha}{4}}\mathbf{S}^{(i)}$ and
$\widehat{\mathbf{S}}^{(i)}$. Therefore by Lemma 3.4 and since
$1.01\cdot\widehat{\tau}$ is a vector of overestimates for leverage scores of
$\mathbf{A}$, with high probability,

$\widetilde{\mathbf{A}}^{\top}\widetilde{\mathbf{A}}\preceq\mathbf{A}^{\top}\widehat{\mathbf{S}}^{\top}\widehat{\mathbf{S}}\mathbf{A}\preceq\mathbf{A}^{\top}\mathbf{A},$

where $\widehat{\mathbf{S}}\in\mathbb{R}^{n\times n}$ is the diagonal matrix
obtained by putting $\widehat{\mathbf{S}}^{(i)}$ together (as block-diagonals
of $\mathbf{S}$). Therefore, for
$\mathbf{B}=[\mathbf{A};\sqrt{\lambda}\mathbf{I}]$ as defined in Step 4, we
have

$\mathbf{B}^{\top}\mathbf{B}=\widetilde{\mathbf{A}}^{\top}\widetilde{\mathbf{A}}+\lambda\cdot\mathbf{I}\preceq\mathbf{A}^{\top}\mathbf{A}+\lambda\cdot\mathbf{I}.$

Thus
$(\mathbf{B}^{\top}\mathbf{B})^{-1}\succeq(\mathbf{A}^{\top}\mathbf{A}+\lambda\cdot\mathbf{I})^{-1}$,
and $\tau^{\mathbf{B}}_{j}(\mathbf{A})\geq\tau^{\lambda}_{j}(\mathbf{A})$, for
all $j\in[n]$. Moreover,

$\tau^{\mathbf{B}}_{j}(\mathbf{A})=\mathbf{a}_{j}^{\top}(\mathbf{B}^{\top}\mathbf{B})^{\dagger}\mathbf{a}_{j}=\mathbf{a}_{j}^{\top}(\mathbf{B}^{\top}\mathbf{B})^{\dagger}\mathbf{B}^{\top}\mathbf{B}(\mathbf{B}^{\top}\mathbf{B})^{\dagger}\mathbf{a}_{j}=\|\mathbf{B}(\mathbf{B}^{\top}\mathbf{B})^{\dagger}\mathbf{a}_{j}\|_{2}^{2}.$

Consequently, by Lemma 1.20, for $\widehat{\mathbf{G}}$, as defined in Step 5
of LABEL:alg:levscoresRefinementSampling, and any $j\in[n]$, since
$\frac{(0.01^{2}-0.01^{3})\cdot r}{4}>20\log n$,

$\frac{101}{99}\tau^{\lambda}_{j}(\mathbf{A})\leq\frac{101}{99}\tau^{\mathbf{B}}_{j}(\mathbf{A})\leq\|\widehat{\mathbf{G}}\mathbf{a}_{j}\|_{2}^{2}\leq\frac{101^{2}}{99^{2}}\cdot\tau^{\mathbf{B}}_{j}(\mathbf{A}),$

with probability at least $1-2n^{-20}$. Taking the union bound, with
probability of at least $1-2n^{-19}$, Section 3.2.2 holds for all $j\in[n]$,
where $\mathbf{a}_{j}\in\mathcal{N}(\widetilde{\mathbf{A}})^{\perp}$. Now by
Proposition 3.6, we have

$|\|\widehat{\mathbf{G}}\mathbf{a}_{j}\|_{2}-\|\mathbf{J}\mathbf{a}_{j}\|_{2}|\leq\sqrt{rd}\cdot\varepsilon\|\mathbf{a}_{j}\|_{2}\leq
d\sqrt{r}\cdot\frac{1}{10^{4}\cdot n^{2}d\sqrt{r}2^{L}}\cdot
2^{L}=\frac{1}{10^{4}\cdot n^{2}}.$

We now consider two cases.

Case 1: $\frac{1}{2n^{2}}\geq\tau_{j}^{\lambda}(\mathbf{A})$.
Then by the construction in Step 7 of LABEL:alg:levscoresRefinementSampling we
have

$\widehat{\tau}^{(\textrm{new})}_{j}\geq\frac{1}{2n^{2}}\geq\tau^{\lambda}_{j}(\mathbf{A}),$

which is the claim of the lemma.

Case 2: $\tau^{\lambda}_{j}(\mathbf{A})>\frac{1}{2n^{2}}$.
Recall that our goal is to show
$\|\mathbf{J}\mathbf{a}_{j}\|_{2}^{2}\geq\tau^{\lambda}_{j}(\mathbf{A}).$ For
the sake of contradiction, suppose instead that
$\|\mathbf{J}\mathbf{a}_{j}\|_{2}^{2}<\tau^{\lambda}_{j}(\mathbf{A})$; then we
have

$\left(\|\mathbf{J}\mathbf{a}_{j}\|_{2}+\frac{1}{10^{4}\cdot
n^{2}}\right)^{2}<\left(\tau^{\lambda}_{j}(\mathbf{A})+\frac{1}{10^{4}\cdot
n^{2}}\right)^{2}\leq\tau^{\lambda}_{j}(\mathbf{A})+\frac{1}{5000n^{2}}+\frac{1}{10^{8}n^{4}}\leq\tau^{\lambda}_{j}(\mathbf{A})+\frac{1}{2500n^{2}},$

where the penultimate inequality follows from the fact that
$\tau^{\lambda}_{j}(\mathbf{A})\leq 1$. Further, observe that we may combine
Section 3.2.2 and Section 3.2.2 to conclude

$\frac{101}{99}\tau^{\lambda}_{j}(\mathbf{A})\leq\|\widehat{\mathbf{G}}\mathbf{a}_{j}\|_{2}^{2}\leq\left(\|\mathbf{J}\mathbf{a}_{j}\|_{2}+\frac{1}{10^{4}\cdot
n^{2}}\right)^{2}.$

Combining Section 3.2.2 with Section 3.2.2, we have
$\frac{2}{99}\tau^{\lambda}_{j}(\mathbf{A})<\frac{1}{2500n^{2}},$ which is a
contradiction to our assumption that
$\tau^{\lambda}_{j}(\mathbf{A})>\frac{1}{2n^{2}}$. This shows
$\|\mathbf{J}\mathbf{a}_{j}\|_{2}^{2}\geq\tau^{\lambda}_{j}(\mathbf{A})$,
which then finishes the proof because of Section 3.2.2. ∎

####  3.2.3 $\ell_{1}$ Bound on Leverage Score Overestimates:
$\|\widehat{\tau}\|_{1}\leq O(d)$

######  Lemma 3.9.

The vector $\widehat{\tau}$ returned by LABEL:alg:levscoresRefinementSampling
satisfies $\|\widehat{\tau}\|_{1}\leq O(d)$.

###### Proof.

Define the vector

$\mathbf{u}^{(\textrm{new})}_{j}:=\min\\{\tau_{j}^{\widetilde{\mathbf{A}}},\widehat{\tau}_{j}\\}.$

Note that the vector $\mathbf{u}^{\textrm{(new)}}_{j}$ is not the same as
$\min\\{\widehat{\tau}^{\widetilde{\mathbf{A}}}_{j},\widehat{\tau}_{j}\\}$ as
seen in Step 7 of LABEL:alg:levscoresRefinementSampling. In Step 7 of
LABEL:alg:levscoresRefinementSampling, if
$\mathbf{a}_{j}\notin\mathcal{N}(\widetilde{\mathbf{A}})^{\perp}$, then we set
$\widehat{\tau}^{\widetilde{\mathbf{A}}}_{j}=\tau^{\widetilde{\mathbf{A}}}_{j}=\infty$,
with high probability. Otherwise, we have
$\mathbf{a}_{j}\in\mathcal{N}(\widetilde{\mathbf{A}})^{\perp}$, and we set
$\widehat{\tau}_{j}=\|\mathbf{J}\mathbf{a}_{j}\|_{2}^{2}$. Then, by Section
3.2.2 and the fact that
$\widetilde{\mathbf{A}}^{\top}\widetilde{\mathbf{A}}\preceq\widetilde{\mathbf{A}}^{\top}\widetilde{\mathbf{A}}+\lambda\cdot\mathbf{I}:=\mathbf{B}^{\top}\mathbf{B}$
and $\mathbf{a}_{j}\in\mathcal{N}(\widetilde{\mathbf{A}})^{\perp}$, we have

$\|\widehat{\mathbf{G}}\mathbf{a}_{j}\|_{2}^{2}\leq\frac{101^{2}}{99^{2}}\cdot\tau^{\mathbf{B}}_{j}\leq\frac{101^{2}}{99^{2}}\cdot\tau^{\widetilde{\mathbf{A}}}_{j}.$

Case 1:
$\|\mathbf{J}\mathbf{a}_{j}\|_{2}<\|\widehat{\mathbf{G}}\mathbf{a}_{j}\|_{2}$.
By combining the assumed inequality with Section 3.2.3, we get

$\|\mathbf{J}\mathbf{a}_{j}\|_{2}^{2}\leq\|\widehat{\mathbf{G}}\mathbf{a}_{j}\|_{2}^{2}\leq\frac{101^{2}}{99^{2}}\cdot\tau^{\widetilde{\mathbf{A}}}_{j}.$

Therefore, we may combine Step 7 of LABEL:alg:levscoresRefinementSampling with
the value of $\widehat{\tau}^{\widetilde{\mathbf{A}}}_{j}$, definition of
$\mathbf{u}_{j}^{\textrm{new}}$ and Section 3.2.3 to get

$\widehat{\tau}^{(\textrm{new})}_{j}\leq
2\cdot\min\\{\widehat{\tau}_{j},\widehat{\tau}^{\widetilde{\mathbf{A}}}_{j}\\}+2\cdot\frac{1}{2n^{2}}=2\cdot\min\\{\widehat{\tau}_{j},\|\mathbf{J}\mathbf{a}_{j}\|_{2}^{2}\\}+\frac{1}{n^{2}}\leq
2\cdot\frac{101^{2}}{99^{2}}\mathbf{u}^{(\textrm{new})}_{j}+\frac{1}{n^{2}}\leq
2.1\cdot\mathbf{u}^{(\textrm{new})}_{j}+\frac{1}{n^{2}}.$

Case 2:
$\|\mathbf{J}\mathbf{a}_{j}\|_{2}\geq\|\widehat{\mathbf{G}}\mathbf{a}_{j}\|_{2}$.
In this case, we may combine Section 3.2.2 and Section 3.2.3 to obtain

$\|\mathbf{J}\mathbf{a}_{j}\|_{2}^{2}\leq\left(\|\widehat{\mathbf{G}}\mathbf{a}_{j}\|_{2}+\frac{1}{10^{4}\cdot
n^{2}}\right)^{2}\leq\left(\frac{101}{99}\cdot\sqrt{\tau^{\widetilde{\mathbf{A}}}_{j}}+\frac{1}{10^{4}\cdot
n^{2}}\right)^{2}\leq
1.05\cdot\tau^{\widetilde{\mathbf{A}}}_{j}+\frac{1}{2500n^{2}},$

where the last inequality follows from $\sqrt{x}\leq\max(x,1)$. Therefore by
the construction of Step 7 of LABEL:alg:levscoresRefinementSampling, the value
of $\widehat{\tau}^{\widetilde{\mathbf{A}}}_{j}$, definition of
$\mathbf{u}_{j}^{\textrm{new}}$, and Section 3.2.3, we get

$\widehat{\tau}^{(\textrm{new})}_{j}\leq
2\cdot\min\\{\widehat{\tau}_{j},\widehat{\tau}^{\widetilde{\mathbf{A}}}_{j}\\}+2\cdot\frac{1}{2n^{2}}=2\cdot\min\\{\widehat{\tau}_{j},\|\mathbf{J}\mathbf{a}_{j}\|_{2}^{2}\\}+\frac{1}{n^{2}}\leq
2\cdot\frac{101^{2}}{99^{2}}\mathbf{u}^{(\textrm{new})}_{j}+\frac{1}{n^{2}}\leq
2.1\cdot\mathbf{u}^{(\textrm{new})}_{j}+\frac{2502}{2500}\cdot\frac{1}{n^{2}}.$

Combining Section 3.2.3 and Section 3.2.3 from the two cases and using Lemma
3.5, we have

$\|\widehat{\tau}^{(\textrm{new})}\|_{1}\leq
2.1\cdot\|\mathbf{u}^{(\textrm{new})}\|_{+}\frac{2502}{2500}\cdot\frac{1}{n}\leq
2.1\cdot\frac{3d}{\alpha}+\frac{2502}{2500}\cdot\frac{1}{n}.$

If $\alpha=1$, since $n,d\geq 1$,

$\|\widehat{\tau}^{(\textrm{new})}\|_{1}\leq 8d.$

Otherwise, $\alpha\geq\frac{25d}{\|\widehat{\tau}\|_{1}}$. Therefore since
$\|\widehat{\tau}\|_{1}\geq 1$ and $n\geq 5$,

$\|\widehat{\tau}^{(\textrm{new})}\|_{1}\leq\frac{6.3}{25}\cdot\|\widehat{\tau}\|_{1}+\frac{2502}{2500}\cdot\frac{1}{n}\leq\frac{1}{2}\cdot\|\widehat{\tau}\|_{1}.$

Therefore after $\left\lceil\log_{2}(n/d)\right\rceil$ iterations, either
$\alpha$ becomes one, which means $\|\widehat{\tau}^{(\textrm{new})}\|_{1}\leq
8d$, or $\|\widehat{\tau}\|_{1}$ is cut by half in each iteration which means
since at the beginning $\|\widehat{\tau}\|_{1}=n$, at the end
$\|\widehat{\tau}\|_{1}\leq d$. Finally in Step 8, we multiply the vector
$\widehat{\tau}$ by $1.01$ which by Proposition 3.7 and choice of $\lambda$ is
guaranteed to be a vector of leverage score overestimates. Moreover by the
above discussion for the final $\widehat{\tau}$, we have
$\|\widehat{\tau}\|_{1}\leq 1.01\cdot 8d\leq 9d$.

∎

####  3.2.4 Correctness and Communication Complexity of Leverage Score
Overestimates

We are now ready to prove the result stated at the start of Section 3.2, which
we first restate below. See 3.2

###### Proof of Lemma 3.2.

We first prove the correctness of LABEL:alg:levscoresRefinementSampling, and
then we bound its communication complexity.

##### Correctness.

The algorithm initializes $\widehat{\tau}=\mathbf{1}$, which makes it a vector
of leverage score overestimates. In Lemma 3.8, we show that throughout the
algorithm, we maintain the invariant $\widehat{\tau}\geq\tau^{\lambda}$. This
implies for each $j\in[n]$:

$1.01\cdot\widehat{\tau}_{j}\geq\frac{\sigma_{\min}^{2}+\lambda}{\sigma_{\min}^{2}}\cdot\widehat{\tau}_{j}\geq\frac{\sigma_{\min}^{2}+\lambda}{\sigma_{\min}^{2}}\cdot\tau^{\lambda}_{j}\geq\tau_{j},$

where the first inequality follows from $\sigma_{\min}\geq\frac{1}{\kappa}$
and $\lambda\leq\frac{1}{100\kappa^{2}}$, the second inequality is by the
invariant $\widehat{\tau}_{j}\geq\tau^{\lambda}_{j}$, and the third inequality
is by Proposition 3.7. This implies that $1.01\cdot\widehat{\tau}$, which is
used in Step 3 of LABEL:alg:levscoresRefinementSampling for sampling
$\widetilde{\mathbf{A}}$, satisfies $1.01\cdot\widehat{\tau}\geq\tau$.
Finally, Lemma 3.9 implies the claimed bound $\|\widehat{\tau}\|_{1}\leq 9d$.

##### Communication complexity.

Having proved the correctness of the returned output, we next bound the
communication complexity of LABEL:alg:levscoresRefinementSampling. First, note
that by construction, the entries of $\widehat{\tau}$ are always powers of two
and in the interval $[\frac{1}{2n^{2}},1]$. Therefore communicating their
summations in Step 1 and Step 2 of LABEL:alg:levscoresRefinementSampling can
be done with $O(s\log n)$ bits in each iteration. Therefore the cost of these
steps over the course of the algorithm is $O(s\log^{2}n)$. The part of Step 3
that forms $\widehat{\alpha}$ and $\alpha$ does not pose any communication
cost since the machines can compute them given $t$.

In Step 3, by choice of $\alpha$, with high probability, we only need to
communicate $O(d\log d)$ rows of $\mathbf{A}$ to the coordinator — note that
the extra factor of $1.01$ for the sampling only increases the number of
sampled rows by a constant factor. Since each row is a vector of size $d$ with
bit complexity $L$, the cost of sending the rows to the coordinator is
$O(d^{2}L\log d)$ in each iteration. Moreover, as we mentioned, since the
entries of $\widehat{\tau}$ are always powers of two and in the interval
$[\frac{1}{2n^{2}},1]$, we can send the entries of leverage score
overestimates to the coordinator with $O(d\log(d)\log(n))$ bits in each
iteration. Therefore the total cost of this step over the course of the
algorithm is $O(d^{2}L\log(d)\log^{2}(n))$. Step 4 does not incur any
communication cost.

Note that $\mathbf{B}^{\top}\mathbf{B}\succeq\lambda\cdot\mathbf{I}$.
Therefore
$(\mathbf{B}^{\top}\mathbf{B})^{-1}\preceq\frac{1}{\lambda}\cdot\mathbf{I}$.
Moreover since $\widehat{\tau}\geq\frac{1}{2n^{2}}$ and by choice of
$\lambda$, we have

$\|\widehat{\mathbf{G}}\|_{F}\leq\operatorname{poly}(n)\cdot
2^{L}\cdot\frac{1}{\lambda}=\operatorname{poly}(n)\cdot\kappa^{2}2^{L}.$

Therefore due to error tolerance for $\mathbf{J}$ described in Step 5, the bit
complexity of $\mathbf{J}$ is $O((L+\log\kappa)\cdot\log n)$. Therefore since
the number of rows and columns of $\mathbf{J}$ is $O(\log n)$ and $d$,
respectively, it can be sent to all the machines with
$O(sd\cdot(L+\log\kappa)\cdot\log n)$ bits of communication, which concludes
the cost of Step 5 in each iteration.

In Step 6, since the primes are less than $(dL)^{2}$, the primes and the
entries of the vectors $\mathbf{v}^{(k)}$ only need $O(\log(dL))$ bits of
communication. Therefore the total cost of Step 6 in each iteration is
$O(d\log(dL)\log n)$.

Finally, note that Step 7 and Step 8 do not pose any communication cost. The
total communication complexity of LABEL:alg:levscoresRefinementSampling is
then bounded by the above discussion.

∎

###  3.3 Richardson-Type Iteration with Preconditioning

The main export of this section is Lemma 3.10, wherein we use an approximate
preconditioner (obtained by sampling according to leverage score overestimate
computed using LABEL:alg:levscoresRefinementSampling and Lemma 3.2) within a
Richardson-type iteration to solve our linear regression problem iteratively.
As alluded to earlier, in LABEL:alg:richardson, in addition to appropriate
roundings, we communicate only the functions of difference of consecutive
vectors in the computation, i.e., $\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}$. This
prevents redundant communication and permits the claimed communication bounds
in Lemma 3.10.

LABEL:alg:richardson differs from the classic Richardson’s iteration in a few
key aspects, which allow for smaller communication complexity. Richardson’s
method requires $\log(\varepsilon^{-1})$ iterations, in each of which, we
multiply the current vector by the Gram matrix of the original matrix and
inverse of the preconditioner. Such a multiplication increases the bit
complexity of the vector by $\log(\kappa)+L$ even if we use a rounded version
of the inverse of the preconditioner. Therefore at iteration $k$, the bit
complexity of the vector is $k(L+\log\kappa)$. Communicating such vectors
between the coordinator and the machines would result in a communication
complexity of

$\Omega(sd\cdot\sum_{k=1}^{\left\lceil\log\tfrac{1}{\varepsilon}\right\rceil}k(L+\log\kappa))=\Omega(sd\cdot(L+\log\kappa)\log^{2}(\varepsilon^{-1})).$

In this section, we show that one of the $\log(\varepsilon^{-1})$ factors in
the above expression can be avoided. Step 1 of LABEL:alg:richardson is a
classic Richardson’s step which gives $\widetilde{\mathbf{x}}^{(k+1)}_{i}$ —
see the example in Equation 3.14, which we display here again. After this
computation, we modify the bits of the $i^{\mathrm{th}}$ coordinate
$\widetilde{\mathbf{x}}^{(k+1)}_{i}$ to obtain $\mathbf{x}^{(k+1)}_{i}$ by
setting the bits with small place values to their values in
$\mathbf{x}^{(k+1)}_{i}$ (cf. Step 2 of LABEL:alg:richardson). The underlying
principle is that since the place value of these bits is small, this
modification does not affect the convergence of Richardson. Moreover, since
our algorithm is converging (proven in Lemma 3.10), the bits with high place
values do not change from one iteration to the next. These bits are
represented by blue in the following equation. Thus, since the machines store
$\mathbf{x}^{(k)}_{i}$, upon receiving the bits of $\mathbf{x}^{(k+1)}_{i}$
represented by red, they can construct $\mathbf{x}^{(k+1)}_{i}$. Our proof
essentially boils down to showing the number of bits in this red middle part
of $\mathbf{x}^{(k+1)}_{i}$ is $O(L+\log\kappa)$ and that the modification of
bits with small place value does not affect the convergence. This then implies
both the communication complexity bounds and correctness as stated in Lemma
3.10. We note that (3.14) does not completely reflect all cases. For example,
it is possible that $\mathbf{x}^{(k)}_{i}=1.0000$ and
$\widetilde{\mathbf{x}}^{(k+1)}_{i}=0.1111$. In this case, it might appear
that the bits with high place values are not equal. However, in the difference
$\widetilde{\mathbf{x}}^{(k+1)}_{i}-\mathbf{x}^{(k)}_{i}=0.0001$, the bits
with high place value are equal to zero. This is carefully written and
analyzed in our algorithms and proofs.

$\begin{split}\widetilde{\mathbf{x}}^{(k+1)}_{i}&=\underbrace{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}101.0110101001}}_{\text{is
same as
$\mathbf{x}^{(k)}_{i}$}}\underbrace{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}011011010010101}}_{\text{is
different than
$\mathbf{x}^{(k)}_{i}$}}\underbrace{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0010000010101\cdots}}_{\text{is
different than $\mathbf{x}^{(k)}_{i}$}}\\\
\mathbf{x}^{(k+1)}_{i}&=\underbrace{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}101.0110101001}}_{\text{is
same as
$\mathbf{x}^{(k)}_{i}$}}\underbrace{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}011011010010101}}_{\text{is
same as
$\widetilde{\mathbf{x}}^{(k+1)}_{i}$}}\underbrace{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}1010100001111\cdots}}_{\text{is
same as $\mathbf{x}^{(k)}_{i}$}}\end{split}$ (3.14)

Input. A matrix $\mathbf{A}:=[\mathbf{A}^{(i)}]\in\mathbb{R}^{n\times d}$,
where $n=\sum_{i=1}^{s}n_{i}$, and vector
$\mathbf{b}:=[\mathbf{b}^{(i)}]\in\mathbb{R}^{n}$, where the $i^{\mathrm{th}}$
machine stores matrix $\mathbf{A}^{(i)}\in\mathbb{R}^{n_{i}\times d}$ and
vector $\mathbf{b}^{(i)}\in\mathbb{R}^{n_{i}}$. The matrix $\mathbf{M}$ with
$\mathbf{A}^{\top}\mathbf{A}\preceq\mathbf{M}\preceq\lambda\mathbf{A}^{\top}\mathbf{A}$,
with $\lambda\geq 1$, stored on the coordinator machine; accuracy parameter
$0<\varepsilon<1$.

Output. A vector $\widehat{\mathbf{x}}\in\mathbb{R}^{d}$ on all machines such
that
$\|\mathbf{A}\widehat{\mathbf{x}}-\mathbf{b}\|_{2}\leq\varepsilon\cdot\|\mathbf{A}(\mathbf{A}^{\top}\mathbf{A})^{-1}\mathbf{A}\mathbf{b}\|_{2}+\min_{\mathbf{x}}\|\mathbf{A}\mathbf{x}-\mathbf{b}\|_{2}$.

Initialize. Set $\mathbf{x}^{(0)}=\vec{0}$,
$\widehat{\varepsilon}=\frac{1}{2\lambda}$, and $T\in\mathbb{N}$ be the
smallest integer such that $(1-\frac{1}{2\lambda})^{T}\leq\varepsilon$.
Moreover each machine $i\in[s]$, computes
$(\mathbf{A}^{(i)})^{\top}\mathbf{b}^{(i)}$ and sends it to the coordinator.

For iterations $k=0,1,\dots,T-1$:

  1. 1

The coordinator computes
$\widetilde{\mathbf{x}}^{(k+1)}=\mathbf{x}^{(k)}-\mathbf{M}^{-1}(\mathbf{A}^{\top}\mathbf{A}\mathbf{x}^{(k)}-\mathbf{A}^{\top}\mathbf{b})$.

  2. 2

The coordinator sets each coordinate $j$ of $\mathbf{x}^{(k+1)}$ equal to the
corresponding coordinate $j$ of $\widetilde{\mathbf{x}}^{(k+1)}$ at all bits
except those with a place value less than
$\widehat{\varepsilon}\cdot\frac{\|\widetilde{\mathbf{x}}^{(k+1)}-\mathbf{x}^{(k)}\|_{\mathbf{M}}}{\lambda\cdot
nd^{2}\cdot 2^{2L+1}}$, at which bits the coordinator sets the bit to the
corresponding bit of coordinate $j$ of $\mathbf{x}^{(k)}$.

  3. 3

The coordinator sends the vector
$\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}\in\mathbb{R}^{d}$ to the machines.

  4. 4

Each machine $i\in[s]$ computes the vector
$(\mathbf{A}^{(i)})^{\top}\mathbf{A}^{(i)}(\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)})\in\mathbb{R}^{d}$
and sends it to the coordinator.

Algorithm 8: Protocol for computing the solution of a linear regression
problem given a preconditioner.

alg]alg:richardson

In the next result, we assume the matrix is full-rank. However, this is not a
limitation of our approach since we can concatenate the matrix with a small-
scaled identity matrix — see Appendix A in [GPV23].

######  Lemma 3.10 (Richardson-Type Iteration).

Given $\varepsilon>0$ and the linear regression setting of Problem 1.3 with
input matrix $\mathbf{A}=[\mathbf{A}_{i}]\in\mathbb{R}^{n\times d}$ and vector
$\mathbf{b}=[\mathbf{b}_{i}]\in\mathbb{R}^{n}$, let
$\mathbf{M}\in\mathbb{R}^{d\times d}$ be a matrix stored in the coordinator
such that
$\mathbf{A}^{\top}\mathbf{A}\preceq\mathbf{M}\preceq\lambda\cdot\mathbf{A}^{\top}\mathbf{A}$,
for constant $\lambda\geq 1$. Then there is an algorithm that outputs a vector
$\widehat{\mathbf{x}}$ such that

$\|\mathbf{A}\widehat{\mathbf{x}}-\mathbf{b}\|_{2}^{2}\leq\varepsilon\cdot\|\mathbf{A}(\mathbf{A}^{\top}\mathbf{A})^{-1}\mathbf{A}^{\top}\mathbf{b}\|_{2}^{2}+\min_{\mathbf{x}}\|\mathbf{A}\mathbf{x}-\mathbf{b}\|_{2}^{2}$

using $\widetilde{O}(sd(L+\log\kappa)\log(\varepsilon^{-1}))$ bits of
communication. Moreover, the vector $\widehat{\mathbf{x}}$ is available on all
the machines at the end of the algorithm.

###### Proof.

We show that LABEL:alg:richardson returns the specified solution with the
specified communication complexity bounds. As a first step, we argue that the
vectors that are required for computation at the coordinator can indeed be
computed from the information that is communicated to the coordinator.

First, note that
$\mathbf{A}^{\top}\mathbf{b}=\sum_{i=1}^{s}(\mathbf{A}^{(i)})^{\top}\mathbf{b}^{(i)}.$
The quantity $\mathbf{A}^{\top}\mathbf{b}$ needs to be communicated to the
coordinator only once (since it does not change through the algorithm). Since
the bit complexities of $\mathbf{A}$ and $\mathbf{b}$ are $L$ and each machine
sends to the coordinator $(\mathbf{A}^{(i)})^{\top}\mathbf{b}^{(i)}$, it takes
a total of $O(sdL\log n)$ to communicate $\mathbf{A}^{\top}\mathbf{b}$.

We now show by induction that the coordinator can compute
$\mathbf{A}^{\top}\mathbf{A}\mathbf{x}^{(k)}$. Initially $\mathbf{x}^{(0)}=0$.
Therefore $\mathbf{A}^{\top}\mathbf{A}\mathbf{x}^{(0)}$, and the coordinator
has this information. Then, assume as the base case of the induction that the
coordinator has the value of $\mathbf{A}^{\top}\mathbf{A}\mathbf{x}^{(k)}$ and
the machines are sending the vectors
$(\mathbf{A}^{(i)})^{\top}\mathbf{A}^{(i)}(\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)})$
to the coordinator (see Step 4 of LABEL:alg:richardson). We observe that

$\displaystyle\mathbf{A}^{\top}\mathbf{A}\mathbf{x}^{(k+1)}$
$\displaystyle=\mathbf{A}^{\top}\mathbf{A}\mathbf{x}^{(k)}+\mathbf{A}^{\top}\mathbf{A}(\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)})=\mathbf{A}^{\top}\mathbf{A}\mathbf{x}^{(k)}+\sum_{i=1}^{s}(\mathbf{A}^{(i)})^{\top}\mathbf{A}^{(i)}(\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}).$

Therefore, if the coordinator has
$\mathbf{A}^{\top}\mathbf{A}\mathbf{x}^{(k)}$ then Step 4 enables it to
compute $\mathbf{A}^{\top}\mathbf{A}\mathbf{x}^{(k+1)}$. Moreover, since
$\mathbf{M}$ and $\mathbf{x}^{(k)}$ are also stored on the coordinator, it can
perform the computation in Step 1.

We now prove that LABEL:alg:richardson converges. The optimal vector
$\mathbf{x}^{\ast}=\arg\min_{\mathbf{x}}\|\mathbf{A}\mathbf{x}-\mathbf{b}\|_{2}$
satisfies
$\mathbf{x}^{\ast}=(\mathbf{A}^{\top}\mathbf{A})^{-1}\mathbf{A}^{\top}\mathbf{b}$.
Subtracting it from the iterate $\widetilde{\mathbf{x}}^{(k+1)}$ and applying
the update from Step 1 of LABEL:alg:richardson in each iteration, we have

$\displaystyle\widetilde{\mathbf{x}}^{(k+1)}-\mathbf{x}^{\ast}$
$\displaystyle=\mathbf{x}^{(k)}-\mathbf{M}^{-1}(\mathbf{A}^{\top}\mathbf{A}\mathbf{x}^{(k)}-\mathbf{A}^{\top}\mathbf{b})-\mathbf{x}^{\ast}=(\mathbf{I}-\mathbf{M}^{-1}\mathbf{A}^{\top}\mathbf{A})(\mathbf{x}^{(k)}-\mathbf{x}^{\ast}).$
(3.15)

To obtain a bound on the rate of shrinkage of the distance (in
$\mathbf{M}$-norm) of the algorithm’s iterate from the true optimizer, we
study the squared $\mathbf{M}$-norm of the right-hand side in Equation 3.15.

$\displaystyle\|(\mathbf{I}-\mathbf{M}^{-1}\mathbf{A}^{\top}\mathbf{A})(\mathbf{x}^{(k)}-\mathbf{x}^{\ast})\|_{\mathbf{M}}^{2}$
$\displaystyle=(\mathbf{x}^{(k)}-\mathbf{x}^{\ast})^{\top}(\mathbf{I}-\mathbf{A}^{\top}\mathbf{A}\mathbf{M}^{-1})\cdot\mathbf{M}\cdot(\mathbf{I}-\mathbf{M}^{-1}\mathbf{A}^{\top}\mathbf{A})(\mathbf{x}^{(k)}-\mathbf{x}^{\ast}).$
(3.16)

To simplify notion, we define
$\mathbf{H}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathbf{M}^{-1/2}\cdot\mathbf{A}^{\top}\mathbf{A}\cdot\mathbf{M}^{-1/2}$.
Therefore, we can check the computation

$(\mathbf{I}-\mathbf{A}^{\top}\mathbf{A}\mathbf{M}^{-1})\cdot\mathbf{M}\cdot(\mathbf{I}-\mathbf{M}^{-1}\mathbf{A}^{\top}\mathbf{A})=\mathbf{M}^{1/2}(\mathbf{I}-\mathbf{H})^{2}\mathbf{M}^{1/2}.$

Moreover, by the assumption
$0\prec\mathbf{A}^{\top}\mathbf{A}\preceq\mathbf{M}\preceq\lambda\cdot\mathbf{A}^{\top}\mathbf{A}$,
the definition of $\mathbf{H}$, and the properties stated in Fact A.1, we have

$\frac{1}{\lambda}\mathbf{I}=\frac{1}{\lambda}\mathbf{M}^{-1/2}\mathbf{M}\mathbf{M}^{-1/2}\preceq\mathbf{H}\preceq\mathbf{M}^{-1/2}\mathbf{M}\mathbf{M}^{-1/2}=\mathbf{I}.$

Therefore
$0\preceq\mathbf{I}-\mathbf{H}\preceq(1-\tfrac{1}{\lambda})\mathbf{I}$;
applying to this Fact A.1, we have

$0\preceq(\mathbf{I}-\mathbf{A}^{\top}\mathbf{A}\mathbf{M}^{-1})\cdot\mathbf{M}\cdot(\mathbf{I}-\mathbf{M}^{-1}\mathbf{A}^{\top}\mathbf{A})\preceq(1-\lambda^{-1})^{2}\cdot\mathbf{M}.$

Plugging Section 3.3 into Equation 3.15 and Equation 3.16 gives

$\|\widetilde{\mathbf{x}}^{(k+1)}-\mathbf{x}^{\ast}\|_{\mathbf{M}}^{2}=\|(\mathbf{I}-\mathbf{M}^{-1}\mathbf{A}^{\top}\mathbf{A})(\mathbf{x}^{(k)}-\mathbf{x}^{\ast})\|_{\mathbf{M}}^{2}\leq(1-\lambda^{-1})^{2}\|\mathbf{x}^{(k)}-\mathbf{x}^{\ast}\|_{\mathbf{M}}^{2}.$

By the update rule in Step 1 of LABEL:alg:richardson and the expression for
$\mathbf{x}^{\ast}$, we have

$\displaystyle\widetilde{\mathbf{x}}^{(k+1)}-\mathbf{x}^{(k)}=-\mathbf{M}^{-1}(\mathbf{A}^{\top}\mathbf{A}\mathbf{x}^{(k)}-\mathbf{A}^{\top}\mathbf{b})=-\mathbf{M}^{-1}(\mathbf{A}^{\top}\mathbf{A}\mathbf{x}^{(k)}-\mathbf{A}^{\top}\mathbf{A}\mathbf{x}^{\ast})=-\mathbf{M}^{-1}\mathbf{A}^{\top}\mathbf{A}(\mathbf{x}^{(k)}-\mathbf{x}^{\ast}).$
(3.19)

By applying the assumption $\mathbf{A}^{\top}\mathbf{A}\preceq\mathbf{M}$ to
Equation 3.19 and using Fact A.1, we get

$\displaystyle\|\widetilde{\mathbf{x}}^{(k+1)}-\mathbf{x}^{(k)}\|_{\mathbf{M}}^{2}$
$\displaystyle=(\mathbf{x}^{(k)}-\mathbf{x}^{\ast})^{\top}\mathbf{A}^{\top}\mathbf{A}\mathbf{M}^{-1}\mathbf{A}^{\top}\mathbf{A}(\mathbf{x}^{(k)}-\mathbf{x}^{\ast})$
$\displaystyle\leq(\mathbf{x}^{(k)}-\mathbf{x}^{\ast})^{\top}\mathbf{A}^{\top}\mathbf{A}(\mathbf{A}^{\top}\mathbf{A})^{-1}\mathbf{A}^{\top}\mathbf{A}(\mathbf{x}^{(k)}-\mathbf{x}^{\ast})$
$\displaystyle=(\mathbf{x}^{(k)}-\mathbf{x}^{\ast})^{\top}\mathbf{A}^{\top}\mathbf{A}(\mathbf{x}^{(k)}-\mathbf{x}^{\ast})$
$\displaystyle\leq\|\mathbf{x}^{(k)}-\mathbf{x}^{\ast}\|_{\mathbf{M}}^{2}.$
(3.20)

We then have

$\displaystyle\|\mathbf{x}^{(k+1)}-\widetilde{\mathbf{x}}^{(k+1)}\|_{\mathbf{M}}\leq\lambda\cdot
nd\cdot
2^{2L}\|\mathbf{x}^{(k+1)}-\widetilde{\mathbf{x}}^{(k+1)}\|_{2}\leq\widehat{\varepsilon}\|\widetilde{\mathbf{x}}^{(k+1)}-\mathbf{x}^{(k)}\|_{\mathbf{M}}\leq\widehat{\varepsilon}\|\mathbf{x}^{(k)}-\mathbf{x}^{\ast}\|_{\mathbf{M}},$
(3.21)

where the first step uses the assumptions that $\mathbf{A}$ has a bit
complexity of $L$ and associated Fact 1.26 and also that
$\mathbf{M}\preceq\lambda\mathbf{A}^{\top}\mathbf{A}$, the second step uses
Step 2 from LABEL:alg:richardson (specifically, the example visually depicted
in Equation 3.14), and the third step uses Equation 3.20. Combining Section
3.3 and Equation 3.21 along with the triangle inequality, we have

$\displaystyle\|\mathbf{x}^{(k+1)}-\mathbf{x}^{\ast}\|_{\mathbf{M}}\leq(1-\lambda^{-1}+\widehat{\varepsilon})\|\mathbf{x}^{(k)}-\mathbf{x}^{\ast}\|_{\mathbf{M}}=(1-\tfrac{1}{2\lambda})\|\mathbf{x}^{(k)}-\mathbf{x}^{\ast}\|_{\mathbf{M}}.$
(3.22)

Since $\lambda$ is a constant, this implies that after
$O(\log(\frac{1}{\varepsilon}))$ iterations, we achieve the required accuracy.
Starting from the zero vector and since
$\mathbf{A}^{\top}\mathbf{A}\preceq\mathbf{M}\preceq\lambda\mathbf{A}^{\top}\mathbf{A}$,
we have

$\|\mathbf{x}^{(T)}-\mathbf{x}^{\ast}\|_{\mathbf{A}^{\top}\mathbf{A}}\leq\varepsilon\|\mathbf{x}^{\ast}\|_{\mathbf{A}^{\top}\mathbf{A}}.$

We then have

$\displaystyle\|\mathbf{A}\mathbf{x}^{(T)}-\mathbf{b}\|_{2}^{2}$
$\displaystyle=\|\mathbf{A}\mathbf{x}^{(T)}-\mathbf{A}(\mathbf{A}^{\top}\mathbf{A})^{-1}\mathbf{A}^{\top}\mathbf{b}-(\mathbf{I}-\mathbf{A}(\mathbf{A}^{\top}\mathbf{A})^{-1}\mathbf{A}^{\top})\mathbf{b}\|_{2}^{2}$
$\displaystyle=\|\mathbf{A}\mathbf{x}^{(T)}-\mathbf{A}(\mathbf{A}^{\top}\mathbf{A})^{-1}\mathbf{A}^{\top}\mathbf{b}\|_{2}^{2}+\|(\mathbf{I}-\mathbf{A}(\mathbf{A}^{\top}\mathbf{A})^{-1}\mathbf{A}^{\top})\mathbf{b}\|_{2}^{2}$
$\displaystyle=\|\mathbf{A}(\mathbf{x}^{(T)}-\mathbf{x}^{\ast})\|_{2}^{2}+\min_{\mathbf{x}}\|\mathbf{A}\mathbf{x}-\mathbf{b}\|_{2}^{2}.$

Then applying Section 3.3 to the above inequality, we have the claimed
convergence bound:

$\displaystyle\|\mathbf{A}\mathbf{x}^{(T)}-\mathbf{b}\|_{2}^{2}\leq\varepsilon^{2}\cdot\|\mathbf{A}(\mathbf{A}^{\top}\mathbf{A})^{-1}\mathbf{A}^{\top}\mathbf{b}\|_{2}^{2}+\min_{\mathbf{x}}\|\mathbf{A}\mathbf{x}-\mathbf{b}\|_{2}^{2}.$

We now bound the bit complexity of each iteration of LABEL:alg:richardson.
Denoting $\sigma_{\textrm{min}}$ to be the minimum singular value of
$\mathbf{M}$, we have

$\displaystyle\|\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}\|_{2}$
$\displaystyle\leq\frac{1}{\sigma_{\textrm{min}}}\|\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}\|_{\mathbf{M}}$
$\displaystyle\leq\frac{1}{\sigma_{\textrm{min}}}\left(\|\mathbf{x}^{(k+1)}-\widetilde{\mathbf{x}}^{(k+1)}\|_{\mathbf{M}}+\|\widetilde{\mathbf{x}}^{(k+1)}-\mathbf{x}^{(k)}\|_{\mathbf{M}}\right)$
$\displaystyle\leq\frac{1+\widehat{\varepsilon}}{\sigma_{\textrm{min}}}\cdot\|\widetilde{\mathbf{x}}^{(k+1)}-\mathbf{x}^{(k)}\|_{\mathbf{M}},$

where the first step uses
$\|\mathbf{u}\|_{\mathbf{M}}\geq\sigma_{\textrm{min}}\|\mathbf{u}\|_{2}$ for
any vector $\mathbf{u}$, the second step is by triangle inequality, and the
third step uses the rounding we performed in Step 2 of LABEL:alg:richardson.
Moreover note that if $\|\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}\|_{2}\neq 0$, by
construction, it is more than
$\widehat{\varepsilon}\cdot\frac{\|\widetilde{\mathbf{x}}^{(k+1)}-\mathbf{x}^{(k)}\|_{\mathbf{M}}}{\lambda\cdot
nd^{2}2^{2L+2}}$. Let $p$ be an integer such that
$2^{p}\leq\widehat{\varepsilon}\cdot\frac{\|\widetilde{\mathbf{x}}^{(k+1)}-\mathbf{x}^{(k)}\|_{\mathbf{M}}}{\lambda\cdot
nd^{2}2^{2L+2}}<2^{p+1}$. Then we have for any $j$ corresponding to a nonzero
entry of $\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}$,

$\displaystyle 1\leq 2^{-p}|(\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)})_{j}|\leq
2^{-p}\cdot\frac{1+\widehat{\varepsilon}}{\sigma_{\min}}\cdot\|\widetilde{\mathbf{x}}^{(k+1)}-\mathbf{x}^{(k)}\|_{\mathbf{M}}<\frac{(1+\widehat{\varepsilon})}{\widehat{\varepsilon}\cdot\sigma_{\min}}\cdot\lambda\cdot
nd^{2}2^{2L+3}.$

Therefore we only need to communicate
$\log(p)+\log(\frac{(1+\widehat{\varepsilon})}{\widehat{\varepsilon}\cdot\sigma_{\min}}\cdot\lambda\cdot
nd^{2}2^{2L+3})$ bits to each machine in Step 4. Finally, in Step 4, when the
machines multiply $\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}$ by
$\mathbf{A}^{\top}\mathbf{A}$ and send it back to the coordinator this just
adds $O(L\log n)$ bits. ∎

###  3.4 Proof of Main Result on High-Accuracy Linear Regression

In this section, we prove our main result for the communication complexity of
linear regression in the point-to-point model of communication.

###### Proof of Theorem 1.8.

We show that LABEL:alg:lin-reg-coordinator-poly-cond returns the correct
output and satisfies the communication complexity bounds. We start with the
proof of correctness.

##### Correctness.

By Lemma 3.2, Step 1 of LABEL:alg:lin-reg-coordinator-poly-cond returns a
vector of overestimates $\widehat{\tau}\geq\tau$ with
$\|\widehat{\tau}\|_{1}\leq 9d$. Then by Lemma 3.4, sampling according to
these leverage score overestimates with $\alpha=100$ produces the matrix
$\widetilde{\mathbf{A}}$ in Step 3 which is a $\frac{1.1}{0.9}$-spectral
approximation of the matrix $\mathbf{A}$. Finally, in Step 4, we use this
spectral approximation to return a high-accuracy solution as proved in Lemma
3.10.

##### Communication complexity.

First, observe that Step 3 incurs no communication. Next, the communication
complexities of Step 1 and Step 4 directly follow from Lemma 3.2 and Lemma
3.10 and are $\widetilde{O}(d^{2}L+sd(L+\log\kappa))$ and
$\widetilde{O}(sd(L+\log\kappa)\log(\varepsilon^{-1}))$, respectively.
Finally, Step 2 requires the communication of the selected rows and their
probabilities, which we now compute. Since Step 1 provides vector
$\widehat{\tau}$ of leverage score overestimates that satisfies
$\|\widehat{\tau}\|_{1}\leq 9d$, we are guaranteed by Lemma 3.4 that, with
high probability, there are only $O(d\log d)$ non-zero rows in Step 2. These
can be communicated with $\widetilde{O}(d^{2}L)$ bits since the bit complexity
of the matrix $\mathbf{A}$ is $L$ and each vector has $d$ coordinates. As per
Definition 3.3, the sampling probabilities can be computed at the coordinator
by knowing the leverage score overestimates for the selected rows at the
coordinator. Therefore since before the last step of
LABEL:alg:levscoresRefinementSampling, the leverage score overestimates are
powers of two, we can communicate these powers of two (that are in the
interval $[\frac{1}{2n^{2}},1]$) with $O(d\log(d)\cdot\log(n))$, bits and then
the coordinator can obtain the leverage score overestimates by multiplying
these by $1.01$ and computing the probabilities accordingly. ∎

##  4 Linear Programming in the Coordinator Model

In this section, we bound the communication complexity of solving the
following linear program in the setup of Problem 1.3:

$\displaystyle\min_{\mathbf{x}\in\mathbb{R}^{n}_{\geq
0}:\mathbf{A}^{\top}\mathbf{x}=\mathbf{b}}\mathbf{c}^{\top}\mathbf{x},$ (4.1)

with $\mathbf{A}=[\mathbf{A}^{(i)}]_{i\in[s]}$,
$\mathbf{c}=[\mathbf{c}^{(i)}]_{i\in[s]}$, the $i^{\mathrm{th}}$ machine
holding $\mathbf{A}^{(i)}\in\mathbb{R}^{n_{i}\times d}$ and
$\mathbf{c}^{(i)}\in\mathbb{R}^{n_{i}}$. We assume that the bit complexities
of $\mathbf{A},\mathbf{b},\mathbf{c}$ are bounded by $L$, and the vector
$\mathbf{b}$ is available to all machines. The main result of this section is
the following.

See 1.10

Our framework for achieving the results in Theorem 1.10 is an adaptation of
the algorithm of [BLSS20] to the coordinator setting, which in turn builds
upon the techniques of [LS14]. While we describe our procedure in more detail
in Section 4.1, essentially each iteration of our interior-point method
involves solving a linear system of the form
$\mathbf{A}^{\top}\mathbf{W}\mathbf{A}\mathbf{v}=\mathbf{A}^{\top}\mathbf{W}^{1/2}\mathbf{g}$,
where $\mathbf{v}$ is the variable vector and $\mathbf{W}$ is a nonnegative
diagonal matrix — see Section 4.1. This is equivalent to solving a linear
regression of the form
$\min_{\mathbf{v}}\|\mathbf{W}^{1/2}\mathbf{A}\mathbf{v}-\mathbf{g}\|_{2}^{2}$.
To this end, we use a randomized approach based on leverage score sampling.
Since the randomness of previous iterations determines the solutions of next
iterations, this might cause adaptive adversary issues — see [LS15]. To
prevent such an issue, we solve the linear regression instance in each
iteration to high-accuracy by Theorem 1.8.

###  4.1 An Overview of Our Algorithm

Set in the framework of primal-dual path-following interior-point method
[Ren88, Ren01], our algorithm maintains a primal feasible point
$\mathbf{x}\in\mathbb{R}^{n}_{\geq 0}$ satisfying (the primal feasibility
condition) $\mathbf{A}^{\top}\mathbf{x}=\mathbf{b}$ and a dual feasible point
$\mathbf{y}\in\mathbb{R}^{d}_{\geq 0}$ satisfying (the dual feasibility
condition) $\mathbf{A}\mathbf{y}-\mathbf{c}\geq 0$ (note that all vector
inequalities in this section are coordinate-wise.) The goal is to steadily
decrease the primal-dual gap

$\mathbf{c}^{\top}\mathbf{x}-\mathbf{b}^{\top}\mathbf{y}=(\mathbf{A}\mathbf{y}+\mathbf{s})^{\top}\mathbf{x}-\mathbf{b}^{\top}\mathbf{y}=\mathbf{s}^{\top}\mathbf{x},$

where $\mathbf{s}:=\mathbf{c}-\mathbf{A}\mathbf{y}$ is the _slack_ variable,
and the other terms cancel out due to the aforementioned feasibility
conditions. At each iteration, a primal-dual path-following algorithm trades
off decreasing the gap $\mathbf{x}^{\top}\mathbf{s}$ against maintaining
feasibility. Collecting all these requirements yields the following optimality
conditions we want the algorithm to satisfy in each iteration:

$\mathbf{x}\odot\mathbf{s}=\mu\cdot\tau(\mathbf{x},\mathbf{s}),\,\,\mathbf{A}^{\top}\mathbf{x}=\mathbf{b},\,\,\mathbf{A}\mathbf{y}+\mathbf{s}=\mathbf{c},\,\,\mathbf{x},\mathbf{s}\geq
0,$

where we use $\odot$ to denote coordinate-wise product of two vectors,
$\tau(\mathbf{x},\mathbf{s})$ is a carefully designed function that tracks the
primal-dual gap, and $\mu$ is a parameter the algorithm gradually decreases.
The (unique) set $(\mathbf{x},\mathbf{y},\mathbf{s})$ satisfying this system
of equations is said to follow the weighted central path. Decreasing $\mu$
puts more emphasis on decreasing the primal-dual gap (i.e., making progress on
the objective value), and the algorithm alternates between decreasing $\mu$
and taking a Newton-like step from the current approximate solution of Section
4.1 towards that with the updated $\mu$.

The choice of weight function $\tau$ plays an immensely critical role in the
convergence rate of the overall algorithm. For example, setting $\tau$ to be
the all-ones vector yields the standard log-barrier-based central path [Ren88,
Meh92, Gon92, NN94, YTM94] which takes $\widetilde{O}(\sqrt{n})$ iterations
for convergence [Kar84, Vai89a]. Drawing on geometric connections between
Lewis weights and ellipsoidal approximations of polytopes, the breakthrough
work of [LS14] developed an algorithm with
$\widetilde{O}(\sqrt{\mathrm{rank}(\mathbf{A})})$ iterations by choosing
$\tau$ to be $\ell_{p}$ Lewis weights of a certain matrix for $p=\Omega(\log
n)$. The work of [BLSS20] simplified this to the leverage scores of a certain
matrix, with an added regularizer for simplicity of analysis, and this is the
weight function we also use:

$\tau(\mathbf{x},\mathbf{s}):=\sigma(\mathbf{S}^{-1/2-\alpha}\mathbf{X}^{1/2-\alpha}\mathbf{A})+\frac{d}{n}\cdot\bm{1},$

where $\sigma(\mathbf{B})$ denotes the vector of leverage scores of matrix
$\mathbf{B}$, the uppercase $\mathbf{X}$ and $\mathbf{S}$ denote,
respectively, the diagonal matrices formed using the vectors $\mathbf{x}$ and
$\mathbf{s}$, and the parameter $\alpha=O(\tfrac{1}{\log(n/d)})$.

In order to improve the proximity of the next point to the true central path
described by Section 4.1, the algorithm attempts to ensure, for the updated
$\mu$, that $\mathbf{x}\odot\mathbf{s}\approx\mu\tau(\mathbf{x},\mathbf{s})$
defined in Section 4.1. This proximity is called _centrality_ [LS14], and to
measure and track it, prior works [CLS21, LSZ19, Bra20, BLSS20] have
successfully used the following “soft-max”-like potential

$\Phi(\mathbf{v})=\sum_{i=1}^{n}\exp(\lambda(\mathbf{v}_{i}-1))+\exp(-\lambda(\mathbf{v}_{i}-1)),$

with $v=\tfrac{\mu\tau(\mathbf{x},\mathbf{s})}{\mathbf{x}\odot\mathbf{s}}$ and
$\lambda=O(\varepsilon^{-1})$; this is what we also therefore use. Making fast
progress along the central path therefore entails updating $\mathbf{x}$ and
$\mathbf{s}$ so that Section 4.1 decreases sufficiently fast.

To achieve this goal, one observes that the updated points
$\mathbf{x}+\delta_{\mathbf{x}}$ and $\mathbf{s}+\delta_{\mathbf{s}}$ must
also satisfy Section 4.1; considering the fact that $\mathbf{x}$ and
$\mathbf{s}$ also satisfy this equation, one infers that the updates
$\delta_{\mathbf{x}}$ and $\delta_{\mathbf{s}}$ must satisfy

$\mathbf{X}\delta_{\mathbf{s}}+\mathbf{S}\delta_{\mathbf{s}}=\delta_{\widetilde{\mu}},\,\,\mathbf{A}^{\top}\delta_{\mathbf{x}}=0,\,\,\mathbf{A}\delta_{\mathbf{y}}+\delta_{\mathbf{s}}=0,$

where we have introduced the notation
$\widetilde{\bm{\mu}}:=\mu\cdot\tau(\mathbf{x},\mathbf{s})=\mu\cdot(\sigma(\mathbf{S}^{-1/2-1/p}\mathbf{X}^{1/2-1/p}\mathbf{A})+\frac{d}{n}1)$,
and
$\delta_{\widetilde{\mu}}\approx\nabla\Phi(\tfrac{\widetilde{\bm{\mu}}}{\mathbf{x}\odot\mathbf{s}})$
is the update step we choose to decrease the potential, thereby improving the
centrality of the next point. Plugging this value of
$\delta_{\widetilde{\bm{\mu}}}$ back into Section 4.1 yields the solution

$\displaystyle\mathbf{\bm{\delta}_{x}}=\sqrt{\tfrac{\mathbf{X}}{\mathbf{S}}}(\mathbf{I}-\mathbf{P})\tfrac{\mathbf{I}}{\sqrt{\mathbf{X}\mathbf{S}}}\cdot\delta_{\widetilde{\mu}}\text{
and
}\mathbf{\bm{\delta}_{s}}=\sqrt{\tfrac{\mathbf{S}}{\mathbf{X}}}\mathbf{P}\tfrac{\mathbf{I}}{\sqrt{\mathbf{X}\mathbf{S}}}\cdot\delta_{\widetilde{\mu}},$
(4.6)

where $\mathbf{P}$ is the orthogonal projection matrix given by the following
closed-form expression

$\displaystyle\mathbf{P}$
$\displaystyle=\sqrt{\tfrac{\mathbf{X}}{\mathbf{S}}}\mathbf{A}\left(\mathbf{A}^{\top}\tfrac{\mathbf{X}}{\mathbf{S}}\mathbf{A}\right)^{-1}\mathbf{A}^{\top}\sqrt{\tfrac{\mathbf{X}}{\mathbf{S}}}.$
(4.7)

Our high-level strategy then is to first compute an initial feasible solution
(Section 4.2) and then iteratively compute the updates described in Equation
4.6. To perform these updates efficiently, we use inverse maintenance, whose
iteration count and communication complexity we discuss in Section 4.3.
Finally in Section 4.4, we put these together to prove Theorem 1.10.

###  4.2 Finding an Initial Feasible Point

We construct our initial iterates following the approach of [BLSS20]. We first
obtain an initial set of primal and dual points $\mathbf{x}\approx\mathbf{1}$
and $\mathbf{s}\approx\mathbf{1}$, which are feasible for the modified linear
program in Definition 4.1, and as stated in Lemma 4.2 (proved in [BLSS20]),
these vectors may be easily transformed to be a set of feasible points for the
original linear program.

######  Definition 4.1 (Theorem 12 of [BLSS20]).

For a linear program
$\min_{\mathbf{A}^{\top}\mathbf{x}=\mathbf{b},\mathbf{x}\geq
0}\mathbf{c}^{\top}\mathbf{x}$ with outer radius $R$, and any
$\varepsilon\in(0,1]$, we define the modified linear program to
$\min_{\overline{\mathbf{A}}^{\top}\overline{\mathbf{x}}=\overline{\mathbf{b}},\overline{\mathbf{x}}\geq
0}\overline{\mathbf{c}}^{\top}\overline{\mathbf{x}},$ where

$\overline{\mathbf{A}}=\begin{bmatrix}\mathbf{A}&\|\mathbf{A}\|_{\mathrm{F}}\cdot\bm{1}_{n}\\\
0&\|\mathbf{A}\|_{\mathrm{F}}\\\
\frac{1}{R}\mathbf{b}^{\top}-\bm{1}^{\top}\mathbf{A}&0\end{bmatrix}\in\mathbb{R}^{(n+2)\times(d+1)},\overline{\mathbf{b}}=\begin{bmatrix}\frac{1}{R}\mathbf{b}\\\
(n+1)\|\mathbf{A}\|_{\mathrm{F}}\end{bmatrix}\in\mathbb{R}^{d+1},\overline{\mathbf{c}}=\begin{bmatrix}\frac{\varepsilon}{\|\mathbf{c}\|_{2}}\cdot\mathbf{c}\\\
0\\\ 1\end{bmatrix}\in\mathbb{R}^{n+2}.$

######  Lemma 4.2 (Theorem 12 of [BLSS20]).

The following are feasible primal and dual vectors for the modified linear
program stated in Definition 4.1:

$\overline{\mathbf{x}}=\bm{1}\in\mathbb{R}^{n+2},\,\,\overline{\mathbf{y}}=\begin{bmatrix}\bm{0}_{d}\\\
-1\end{bmatrix}\in\mathbb{R}^{d+1},\,\,\overline{\mathbf{s}}=\begin{bmatrix}\bm{1}_{n}+\frac{\varepsilon}{\|c\|_{2}}\cdot\mathbf{c}\\\
1\\\ 1\end{bmatrix}\in\mathbb{R}^{n+2}.$

Let $(\widehat{\mathbf{x}},\widehat{\mathbf{y}},\widehat{\mathbf{s}})$ be an
arbitrary set of primal dual vectors of the modified linear program, and let
$\rho:=\frac{1}{\mu}\cdot\|\overline{\mathbf{A}}^{\top}\widehat{\mathbf{x}}-\overline{\mathbf{b}}\|^{2}_{(\overline{\mathbf{A}}^{\top}\widehat{\mathbf{X}}\widehat{\mathbf{S}}^{-1}\overline{\mathbf{A}})^{-1}}=O(1)$.
Then,

$\|\widehat{\mathbf{x}}\|_{\infty}=O(n)\cdot(1+\rho/\mu).$

Moreover for $\mu<\frac{\varepsilon^{2}}{8d}$ and
$\widehat{\mathbf{x}},\widehat{\mathbf{s}}$ such that
$\widehat{\mathbf{x}}\odot\widehat{\mathbf{s}}\approx_{0.5}\mu\cdot\tau(\widehat{\mathbf{x}},\widehat{\mathbf{s}})$
(where $\tau$ is defined in Section 4.1), if we set
$\mathbf{x}=R\cdot\widehat{\mathbf{x}}_{1:n}$ (where
$\widehat{\mathbf{x}}_{1:n}$ is the vector of the first $n$ coordinates of
$\widehat{\mathbf{x}}$), then $\mathbf{x}\geq 0$ is an approximate solution to
the original linear program in the following sense:

$\mathbf{c}^{\top}\mathbf{x}\leq\min_{\mathbf{A}^{\top}\mathbf{x}=\mathbf{b},\mathbf{x}\geq
0}\mathbf{c}^{\top}\mathbf{x}+O(nR\cdot\|\mathbf{c}\|_{2})\cdot(\sqrt{\rho}+\varepsilon),\,\,\|\mathbf{A}^{\top}\mathbf{x}-\mathbf{b}\|_{2}\leq
O(n^{2})\cdot(\|\mathbf{A}\|_{\mathrm{F}}\cdot
R+\|\mathbf{b}\|_{2})\cdot(\sqrt{\rho}+\varepsilon).$

It is critical to note that although
$(\overline{\mathbf{x}},\overline{\mathbf{y}},\overline{\mathbf{s}})$ in Lemma
4.2 is feasible, it is not necessarily on the central path (or even near it)
since $\overline{\mathbf{x}}\odot\overline{\mathbf{s}}$ is not necessarily
near $\mu\cdot\tau(\overline{\mathbf{x}},\overline{\mathbf{s}})$. To enforce
this centrality with $\mathbf{x}=\mathbf{1}_{n+2}$ and $\mu=1$, we must have
$\mathbf{s}\approx\tau(\mathbf{1},\mathbf{s}),$ which corresponds roughly to
the definition of regularized $\ell_{p}$ Lewis weights of $\mathbf{A}$ for
$p=\tfrac{1}{1+\alpha}.$ The fact that $\alpha>0$ corresponds to $p<1$, which
is a range of $p$ for which $\ell_{p}$ Lewis weights may be efficiently
computed via a simple fixed-point iteration [CP15]. We note that for the
interval complementary to that in [CP15], [FLPS22] provides an efficient (but
different) algorithm for computing $\ell_{p}$ Lewis weights. Since in our
setting, $p<1$, we use the algorithm by [CP15], whose main guarantee we
restate below.

######  Lemma 4.3 (Lemma 3.2 of [CP15]).

Given $\mathbf{A}\in\mathbb{R}^{n\times d}$, $0<p<4$, and vectors
$\mathbf{v},\mathbf{w}\in\mathbb{R}^{n}$ such that
$\mathbf{v}\approx_{\zeta}\mathbf{w}$, define
$\widehat{\mathbf{v}}_{i}=(\mathbf{a}_{i}^{\top}(\mathbf{A}^{\top}\mathbf{V}^{1-2/p}\mathbf{A})^{-1}\mathbf{a}_{i})^{p/2}$
and
$\widehat{\mathbf{w}}_{i}=(\mathbf{a}_{i}^{\top}(\mathbf{A}^{\top}\mathbf{W}^{1-2/p}\mathbf{A})^{-1}\mathbf{a}_{i})^{p/2}.$
Then $\widehat{\mathbf{v}}$ and $\widehat{\mathbf{w}}$ satisfy the
approximation

$\widehat{\mathbf{v}}\approx_{\zeta\cdot{|p/2-1|}}\widehat{\mathbf{w}},$

where, for some given $\alpha>0$, we use the notation $x\approx_{\alpha}y$ to
denote $y\cdot e^{-\alpha}\leq x\leq y\cdot e^{\alpha}$.

The proof of [CP15] for convergence to the Lewis weights, for $p\in(0,4)$,
through the iterative process described above is based on the fact that
$|p/2-1|<1$ and hence results in a contraction. We also use this observation
to bound the communication complexity of computing the (approximate)
regularized Lewis weights in Lemma 4.7, adapted from Corollary 3.4 of [CP15]
and Theorem 13 of [BLSS20]. Before proving this lemma, we first establish some
notation for easier readability, followed by two technical results, which we
invoke in order to prove Lemma 4.7.

######  Definition 4.4.

Given a vector $\mathbf{x}$ and a diagonal matrix $\mathbf{X}$ formed using
$\mathbf{x}$, we define the functions

$\mathbf{q}_{i}(\mathbf{x}):=\mathbf{a}_{i}^{\top}(\mathbf{A}^{\top}\mathbf{X}^{1-2/p}\mathbf{A})^{-1}\mathbf{a}_{i}\text{
and
}\mathcal{T}_{i}(\mathbf{x}):=(\mathbf{q}_{i}(\mathbf{x})+\eta\mathbf{x}_{i}^{2/p-1})^{p/2},$

and the following approximation to $\mathcal{T}$, with both additive and
multiplicative error,

$\widetilde{\mathcal{T}}_{i}(\mathbf{x}):=(\widetilde{\mathbf{q}}_{i}(\mathbf{x})+\eta\mathbf{x}_{i}^{2/p-1})^{p/2},\text{
where
}\mathbf{q}_{i}(\mathbf{x})+e^{-\varepsilon_{1}}\cdot\frac{\eta^{2/p}}{n^{2}}\leq\widetilde{\mathbf{q}}_{i}(\mathbf{x})\leq
e^{\varepsilon_{1}}\cdot\mathbf{q}_{i}(\mathbf{x})+\frac{\eta^{2/p}}{n^{2}}.$

We remark that our $\mathcal{T}$ is identical to $T$ in Theorem $13$ of
[BLSS20].

######  Claim 4.5 ([BLSS20]).

For some $\alpha>0$, let $x$ and $y$ be positive numbers satisfying
$x\approx_{\alpha}y$. Then, for $\lambda,u,k>0$, we have $\lambda
x\approx_{\alpha}\lambda y$, $(x+u)\approx_{\alpha}(y+u)$, and
$x^{k}\approx_{\alpha k}y^{k}$.

Next, we extend Claim 4.5 to the following lemma, which essentially says that
if two vectors $\mathbf{x}$ and $\mathbf{y}$ are $\alpha$-approximations of
each other (in the sense stated in Lemma 4.3), then
$\widetilde{\mathcal{T}}(\mathbf{x})$ and
$\widetilde{\mathcal{T}}(\mathbf{y})$ are also close. We use this lemma to
inductively prove Lemma 4.7.

######  Lemma 4.6.

Consider the notation from Definition 4.4. Suppose for some $\alpha>0$, the
vectors $\mathbf{x}$ and $\mathbf{y}$ satisfy
$\mathbf{x}\approx_{\alpha}\mathbf{y}$. Then we have
$\widetilde{\mathcal{T}}_{i}(\mathbf{x})\approx_{p/2(\varepsilon_{1}+\alpha|1-2/p|)}\widetilde{\mathcal{T}}_{i}(\mathbf{y})$,
where $\varepsilon_{1}=\begin{cases}p\cdot\varepsilon/4&\text{ if }p\leq 2\\\
\frac{\varepsilon\cdot(4-p)}{4}&\text{ if }p>2\end{cases}$.

###### Proof.

Since $\mathbf{x}\approx_{\alpha}\mathbf{y}$, by Claim 4.5, we have

$\mathbf{q}_{i}(\mathbf{y})\approx_{\alpha|1-2/p|}\mathbf{q}_{i}(\mathbf{x})\text{
and }\mathbf{y}_{i}^{2/p-1}\approx_{\alpha|1-2/p|}\mathbf{x}_{i}^{2/p-1}.$

Let $\gamma=\varepsilon_{1}+\alpha|1-2/p|$ for some $\varepsilon_{1}>0.$ Then
the above approximations may equivalently be expressed as

$\mathbf{q}_{i}(\mathbf{x})\approx_{\gamma-\varepsilon_{1}}\mathbf{q}_{i}(\mathbf{y})\text{
and
}\mathbf{x}_{i}^{2/p-1}\approx_{\gamma-\varepsilon_{1}}\mathbf{y}_{i}^{2/p-1}.$

This implies, for the above choice of $\gamma,$ the following inequalities
hold:

$e^{\varepsilon_{1}}\mathbf{q}_{i}(\mathbf{y})+\frac{\eta^{2/p}}{n^{2}}+\eta\mathbf{y}_{i}^{2/p-1}\leq
e^{\gamma}\mathbf{q}_{i}(\mathbf{x})+e^{\gamma-\varepsilon_{1}}\frac{\eta^{2/p}}{n^{2}}+\eta
e^{\gamma}\mathbf{x}_{i}^{2/p-1},$
$e^{\varepsilon_{1}}\mathbf{q}_{i}(\mathbf{x})+\frac{\eta^{2/p}}{n^{2}}+\eta\mathbf{x}_{i}^{2/p-1}\leq
e^{\gamma}\mathbf{q}_{i}(\mathbf{y})+e^{\gamma-\varepsilon_{1}}\cdot\frac{\eta^{2/p}}{n^{2}}+\eta
e^{\gamma}\mathbf{y}_{i}^{2/p-1}.$

Chaining these inequalities with the definition of $\widetilde{\mathcal{T}}$
in Definition 4.4 finishes the proof. ∎

######  Lemma 4.7.

Given $1>\eta>0$, $0<\varepsilon<0.5$, $p\in(0,4)$ and the setting of Problem
1.3 with input matrix $\mathbf{A}=[\mathbf{A}^{(i)}]\in\mathbb{R}^{n\times
d}$, there is a randomized algorithm that, with high probability, outputs a
vector $\widehat{\mathbf{w}}$ such that

$\widehat{\mathbf{w}}\approx_{\varepsilon}\sigma(\widehat{\mathbf{W}}^{1/2-1/p}\mathbf{A})+\eta\cdot\bm{1}.$

The number of bits of communication this algorithm uses is

$\widetilde{O}\left((d^{2}L+sd(L+\log\kappa+p^{-1}\log(\eta^{-1})))\cdot\frac{\log(\eta^{-1}\varepsilon^{-1}p^{-1})}{1-|p/2-1|}\right).$

Each machine has access to $\widehat{\mathbf{w}}_{i}\in\mathbb{R}^{n_{i}}$,
which is the part of $\widehat{\mathbf{w}}$ corresponding to the rows of
$\mathbf{A}^{(i)}$.

###### Proof.

For conciseness inside this proof, we use the notation from Definition 4.4.
The algorithm that achieves this lemma’s stated guarantee constructs iterates
$\widehat{\mathbf{w}}^{(t)}$ as follows:

$\widehat{\mathbf{w}}^{(0)}=\eta\cdot\bm{1}\,\,\text{ and
}\,\,\widehat{\mathbf{w}}^{(t)}=\widetilde{\mathcal{T}}_{i}(\widehat{\mathbf{w}}^{(t-1)}).$

Similarly, the algorithm constructs iterates $\mathbf{w}^{(t)}$ as follows:

$\mathbf{w}^{(0)}=\eta\cdot\bm{1}\,\,\text{ and
}\,\,\mathbf{w}^{(t)}_{i}=\mathcal{T}_{i}(\widehat{\mathbf{w}}^{(t-1)}),$

noting that the difference from Section 4.2 is in the use of $\mathcal{T}$
instead of $\widetilde{\mathcal{T}}$. Starting with the definition of
$\widehat{\mathbf{w}}^{(0)}$ and update rule in Section 4.2, we obtain the
following lower bound on $\widehat{\mathbf{w}}_{i}^{(1)}$:

$\displaystyle\widehat{\mathbf{w}}_{i}^{(1)}$
$\displaystyle=(\widetilde{\mathbf{q}}_{i}(\widehat{\mathbf{w}}_{i}^{(0)})+\eta^{2/p})^{p/2}\geq(\eta^{2/p-1}\cdot\mathbf{a}_{i}^{\top}(\mathbf{A}^{\top}\mathbf{A})^{-1}\mathbf{a}_{i}+\eta^{2/p})^{p/2}\geq\eta=\widehat{\mathbf{w}}_{i}^{(0)}.$
(4.12)

To obtain an upper bound on $\widehat{\mathbf{w}}_{i}^{(1)}$, we start with
the definition of $\widehat{\mathbf{w}}^{(1)}$, followed by the upper bound
from Definition 4.4 with $\varepsilon_{1}<0.04$:

$\displaystyle\widehat{\mathbf{w}}_{i}^{(1)}$
$\displaystyle=(\widetilde{\mathbf{q}}_{i}(\widehat{\mathbf{w}}_{i}^{(0)})+\eta^{2/p})^{p/2}\leq(e^{\varepsilon_{1}}\cdot\eta^{2/p-1}\mathbf{a}_{i}^{\top}(\mathbf{A}^{\top}\mathbf{A})^{-1}\mathbf{a}_{i}+\frac{\eta^{2/p}}{n^{2}}+\eta^{2/p})^{p/2}$
$\displaystyle\leq\eta\cdot(1+1.1\cdot\frac{1}{\eta}+\frac{1}{n^{2}})^{p/2}\leq\eta\cdot\exp((p/2)\cdot(1.1\cdot\frac{1}{\eta}+\frac{1}{n^{2}}))$
$\displaystyle=\widehat{\mathbf{w}}^{(0)}_{i}\cdot\exp((p/2)\cdot(1.1\cdot\frac{1}{\eta}+\frac{1}{n^{2}})),$
(4.13)

where the second step uses Definition 4.4, the third step uses the upper bound
$\mathbf{a}_{i}(\mathbf{A}^{\top}\mathbf{A})^{-1}\mathbf{a}_{i}\leq 1$, and
the fourth step uses the inequality $1+x\leq e^{x}$. Therefore, from Equation
4.12 and Equation 4.13, we may conclude that
$\widehat{\mathbf{w}}_{i}^{(1)}\approx_{\beta}\widehat{\mathbf{w}}_{i}^{(0)}$
for $\beta=(p/2)\cdot(1.1\cdot\frac{1}{\eta}+\frac{1}{n^{2}})$. We may now
invoke Lemma 4.6 on $\widehat{\mathbf{w}}^{(1)}$ and
$\widehat{\mathbf{w}}^{(0)}$ inductively to claim that

$\widehat{\mathbf{w}}^{(t+1)}\approx_{\beta\cdot|p/2-1|^{t}+\theta_{t}}\widehat{\mathbf{w}}^{(t)},$

where $\theta_{t}=\varepsilon_{1}\cdot(p/2)\cdot\sum_{j=1}^{t}|p/2-1|^{j-1}$.
If $p\leq 2$, then this formula implies that $\theta_{t}\leq\varepsilon_{1}$;
if $p>2$, then we have $\theta_{t}\leq\varepsilon_{1}\cdot\frac{p}{4-p}$. For
$p\leq 2$, we pick $\varepsilon_{1}=p\cdot\varepsilon/4$, and if $p>2$, we
pick $\varepsilon_{1}=\frac{\varepsilon\cdot(4-p)}{4}$. Since $0<p<4$,
$0<\eta<1$, and $|\frac{p}{2}-1|\leq\exp(|\frac{p}{2}-1|-1)$, after
$O(\frac{\log(\eta^{-1}\varepsilon^{-1}p^{-1})}{1-|p/2-1|})$ iterations, we
have $k$ such that

$\widehat{\mathbf{w}}^{(k)}\approx_{p\cdot\varepsilon/2}\widehat{\mathbf{w}}^{(k-1)}.$

Next, based on the update rule for $\mathbf{w}_{i}^{(t)}$ in Section 4.2, we
infer
$\mathbf{w}^{(k+1)}_{i}=(\widehat{\mathbf{w}}^{(k)}_{i})^{1-p/2}\cdot(\sigma_{i}((\widehat{\mathbf{W}}^{(k)})^{1/2-1/p}\mathbf{A})+\eta)^{p/2}$.
By Lemma 4.6 and arguments similar to above,

$\mathbf{w}^{(k+1)}\approx_{p\cdot\varepsilon/2}\widehat{\mathbf{w}}^{(k)}.$

Therefore

$(\widehat{\mathbf{w}}^{(k)}_{i})^{1-p/2}\cdot(\sigma_{i}((\widehat{\mathbf{W}}^{(k)})^{1/2-1/p}\mathbf{A})+\eta)^{p/2}\approx_{p\cdot\varepsilon/2}\widehat{\mathbf{w}}^{(k)}_{i}.$

Thus

$(\sigma_{i}((\widehat{\mathbf{W}}^{(k)})^{1/2-1/p}\mathbf{A})+\eta)^{p/2}\approx_{p\cdot\varepsilon/2}(\widehat{\mathbf{w}}^{(k)}_{i})^{p/2}.$

Raising above to the power of $2/p$, we have

$\sigma_{i}((\widehat{\mathbf{W}}^{(k)})^{1/2-1/p}\mathbf{A})+\eta\approx_{\varepsilon}\widehat{\mathbf{w}}^{(k)}_{i}.$

The communication complexity follows from the number of iterations and Lemma
3.2. Note that the guarantee of Definition 4.4 follows from the proof of Lemma
3.2 and the term $sd\cdot\frac{\log(\eta^{-1})}{p})$ appears because we need a
$\frac{\eta^{2/p}}{n^{2}}$ additive error for the leverage score computation.
∎

Input.  A matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ and vectors
$\mathbf{b}\in\mathbb{R}^{d},\mathbf{c}\in\mathbb{R}^{n}$ with parameters in
Theorem 1.10; Error parameters $0<\varepsilon<1$.

Output.  A vector $\widehat{\mathbf{x}}\in\mathbb{R}^{n}_{\geq 0}$ satisfying
Theorem 1.10.

  1. 1

Let $\alpha:=\tfrac{1}{4\log(4n/d)}$,
$\lambda:=\tfrac{\alpha}{32000}\cdot\log(2^{16}n\tfrac{\sqrt{d}}{\alpha^{2}})$,
and $\gamma:=\min\\{\tfrac{\alpha}{64000},\tfrac{\alpha}{50\lambda}\\}$,
$\widehat{\varepsilon}=0.1$

  2. 2

Let
$\overline{\mathbf{A}},\overline{\mathbf{b}},\overline{\mathbf{c}},\overline{\mathbf{x}},\overline{\mathbf{y}},\overline{\mathbf{s}}$
be as defined in Definition 4.1 for the modified linear program.

  3. 3

Set $\widehat{\mathbf{s}}\in\mathbb{R}^{n+2}$ to a vector with
$\widehat{\mathbf{s}}\approx_{\widehat{\varepsilon}}\sigma(\widehat{\mathbf{S}}^{1/2-\alpha}\mathbf{A})+\frac{d}{n}\cdot\bm{1}$
(see Lemma 4.7).

  4. 4

Set $\widehat{\mathbf{c}}=\widehat{\mathbf{s}}$,
$\widehat{\bm{\tau}}=\widehat{\mathbf{s}}$, $\widehat{\mathbf{y}}=\bm{0}$,
$\widehat{\mu}=1$

  5. 5

Let
$(\widehat{\mathbf{x}}^{(\text{end})},\widehat{\mathbf{y}}^{(\text{end})},\widehat{\mathbf{s}}^{(\text{end})},\widehat{\bm{\tau}}^{(\text{end})},\widehat{\mu}^{(\text{end})})=\textsc{IPM}(\overline{\mathbf{A}},\overline{\mathbf{b}},\widehat{\mathbf{c}},\overline{\mathbf{x}},\widehat{\mathbf{y}},\widehat{\mathbf{s}},\widehat{\bm{\tau}}^{(\text{end})},\widehat{\mu},\theta(n^{2}\sqrt{d}/(\gamma\alpha^{2})))$

  6. 6

Set $\mathbf{x}^{(0)}=\overline{\mathbf{x}}^{(\text{end})}$,
$\mathbf{y}^{(0)}=\overline{\mathbf{y}}^{(\text{end})}$ and
$\mathbf{s}^{(0)}=\overline{\mathbf{s}}^{(\text{end})}+\overline{\mathbf{c}}-\widehat{\mathbf{c}}$,
$\bm{\tau}^{(0)}=\widehat{\bm{\tau}}^{(\text{end})}$,
$\mu^{(0)}=\widehat{\mu}^{(\text{end})}$

  7. 7

Let
$(\mathbf{x}^{(\text{end})},\mathbf{y}^{(\text{end})},\mathbf{s}^{(\text{end})},\bm{\tau}^{(\text{end})},\mu^{(\text{end})})=\textsc{IPM}(\overline{\mathbf{A}},\overline{\mathbf{b}},\overline{\mathbf{c}},\mathbf{x}^{(0)},\mathbf{y}^{(0)},\mathbf{s}^{(0)},\bm{\tau}^{(0)},\mu^{(0)},\varepsilon^{2}/(512\cdot
n^{4}d))$

  8. 8

return $R\cdot\mathbf{x}^{(\text{end})}_{1:n}$

Algorithm 9: LP Solver

alg]alg:ipm-init

Input. A matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ and vectors
$\mathbf{d}^{(0)},\mathbf{d}^{(1)},\ldots,\mathbf{d}^{(r)}\in\mathbb{R}^{n}$
and their corresponding diagonal matrices
$\mathbf{D}^{(0)},\mathbf{D}^{(1)},\ldots,\mathbf{D}^{(r)}$. For $i\geq 1$,
each $\mathbf{d}^{(i)}$ is received after we returned the output for
$\mathbf{d}^{(i-1)}$.

Output. For $i\geq 0$, the output is a spectral approximation to
$(\mathbf{A}^{\top}\mathbf{D}^{(i)}\mathbf{A})^{-1}$.

  1. 1

Let $\gamma=1000C\cdot\log d$ (where $C$ is some absolute constant)

  2. 2

Compute (with high probability using LABEL:alg:levscoresRefinementSampling)
$\tau^{(\text{apr})}\in\mathbb{R}^{n}$ such that
$0.99\cdot\tau^{(\text{apr})}_{i}\leq\tau_{\mathbf{A}}(\mathbf{d}^{(0)})_{i}\leq
1.01\cdot\tau^{(\text{apr})}_{i}$

  3. 3

Set $\mathbf{d}^{(\text{old})}=\mathbf{d}^{(0)}$ and
$\tau^{(\text{old})}=\tau^{(\text{apr})}$

  4. 4

For each $i\in[n]$, let

$\mathbf{h}_{i}^{(0)}=\begin{cases}\mathbf{d}_{i}^{(0)}/\min\\{1,\gamma\cdot\tau_{i}^{\text{apr}}\\}&,\text{
with probability }\min\\{1,\gamma\cdot\tau_{i}^{\text{apr}}\\}\\\ 0&,\text{
otherwise}\end{cases}$

  5. 5

Set $\mathbf{K}^{(0)}=(\mathbf{A}^{\top}\mathbf{H}^{(0)}\mathbf{A})^{-1}$

  6. 6

For $j=1,\ldots,r$

    1. (a)

Use Richardson’s iteration and $\mathbf{K}^{(j-1)}$ to compute JL sketchings
of the form
$(\mathbf{A}^{\top}\mathbf{D}^{(j)}\mathbf{A})^{-1}\mathbf{A}^{\top}(\mathbf{D}^{(j)})^{1/2}\mathbf{G}$
(Each machine computes its own part of
$\mathbf{A}^{\top}(\mathbf{D}^{(j)})^{1/2}\mathbf{G}$ and sends it to the
coordinator and then there is a back and forth to solve the linear system
using Richardson). Then compute
$(\mathbf{D}^{(j)})^{1/2}\mathbf{A}(\mathbf{A}^{\top}\mathbf{D}^{(j)}\mathbf{A})^{-1}\mathbf{A}^{\top}(\mathbf{D}^{(j)})^{1/2}\mathbf{G}$
and use it to compute $\tau^{(\text{apr})}\in\mathbb{R}^{n}$ such that
$0.99\cdot\tau^{(\text{apr})}_{i}\leq\tau_{\mathbf{A}}(\mathbf{d}^{(j)})_{i}\leq
1.01\cdot\tau^{(\text{apr})}_{i}$

    2. (b)

For each $i\in[n]$, if
$|\tau^{(\text{apr})}_{i}-\tau^{(\text{old})}_{i}|/\tau^{(\text{old})}_{i}>0.1$
or
$|\mathbf{d}^{(j)}_{i}-\mathbf{d}^{(\text{old})}_{i}|/\mathbf{d}^{(\text{old})}_{i}>0.1$,
then

      1. i.

Set $\mathbf{d}^{(\text{old})}_{i}=\mathbf{d}^{(j)}_{i}$,
$\tau^{(\text{old})}_{i}=\tau^{(\text{apr})}_{i}$, and

$\mathbf{h}_{i}^{(j)}=\begin{cases}\mathbf{d}_{i}^{(j)}/\min\\{1,\gamma\cdot\tau_{i}^{\text{apr}}\\}&,\text{
with probability }\min\\{1,\gamma\cdot\tau_{i}^{\text{apr}}\\}\\\ 0&,\text{
otherwise}\end{cases}$

    3. (c)

else

      1. i.

$\mathbf{h}_{i}^{(j)}=\mathbf{h}_{i}^{(j-1)}$

    4. (d)

Set $\mathbf{K}^{(j)}=(\mathbf{A}^{\top}\mathbf{H}^{(j)}\mathbf{A})^{-1}$

Algorithm 10: Inverse Maintenance for LP

alg]alg:inv-maintenance

As we discussed, after computing $\widehat{\mathbf{s}}$ such that
$\widehat{\mathbf{s}}\approx_{\widehat{\varepsilon}}\sigma(\widehat{\mathbf{S}}^{-1/2-\alpha}\overline{\mathbf{A}})+\frac{d}{n}\bm{1}$,
we have a point near the central path for the modified linear program for the
modified objective vector $\widehat{\mathbf{c}}=\widehat{\mathbf{s}}$.
Therefore, we can run the IPM so that $\mu$ is small enough. This part is
illustrated as Step 7 of LABEL:alg:ipm-init. By notation of LABEL:alg:ipm-
init, and the guarantees of our IPM, we have

$\widehat{\mathbf{x}}^{(\text{end})}\odot\widehat{\mathbf{s}}{(\text{end})}=\widehat{\mu}^{(\text{end})}\cdot\widehat{\bm{\tau}}^{(\text{end})}.$

Since we have
$\mathbf{A}\widehat{\mathbf{y}}^{(\text{end})}+\widehat{\mathbf{s}}^{(\text{end})}=\widehat{\mathbf{c}}$
(this equality is exact since we do not explicitly update
$\widehat{\mathbf{y}}$) and $\widehat{\mathbf{c}}=\widehat{\mathbf{s}}$,
$\mathbf{A}\widehat{\mathbf{y}}^{(\text{end})}+\mathbf{s}^{(0)}-\overline{\mathbf{c}}+\widehat{\mathbf{c}}=\widehat{\mathbf{c}}$.
Thus
$\mathbf{A}\widehat{\mathbf{y}}^{(\text{end})}+\mathbf{s}^{(0)}=\overline{\mathbf{c}}$.
Moreover, by construction
$\mathbf{s}^{(0)}=\overline{\mathbf{s}}^{(\text{end})}+\overline{\mathbf{c}}-\widehat{\mathbf{c}}$.
Therefore
$\frac{\mathbf{s}^{(0)}-\mathbf{s}^{(\text{end})}}{\mathbf{s}^{(0)}}=\frac{\overline{\mathbf{c}}-\widehat{\mathbf{c}}}{\mathbf{s}^{(0)}}$.
Since by the guarantees of the IPM
$\mathbf{s}^{(0)}\odot\mathbf{x}^{(0)}\approx_{\widehat{\varepsilon}}\mu^{(0)}\cdot\bm{\tau}^{(0)}$,
and $\widehat{\varepsilon}\leq 0.5$, for each $i\in[n]$, we have

$\mathbf{s}^{(0)}_{i}\geq\frac{\mu^{(0)}}{2\mathbf{x}^{(0)}_{i}}\cdot\frac{d}{n}\geq\frac{\mu}{\Omega(n^{2})},$

where the second inequality follows from Lemma 4.2 by
$\|\mathbf{x}\|_{\infty}=O(n)$. Moreover, note that
$\|\overline{\mathbf{c}}\|_{\infty}\leq 1$ and
$\|\widehat{\mathbf{c}}\|_{\infty}=\|\widehat{\mathbf{s}}\|_{\infty}\leq 3$.
Therefore by picking the appropriate constant in Step 7 of LABEL:alg:ipm-init,
for any constant $\beta$, we have

$\|\frac{\mathbf{s}^{(0)}-\mathbf{s}^{(\text{end})}}{\mathbf{s}^{(0)}}\|_{\infty}=\|\frac{\overline{\mathbf{c}}-\widehat{\mathbf{c}}}{\mathbf{s}^{(0)}}\|_{\infty}\leq\frac{16n^{2}}{\mu^{(0)}}\leq\frac{\gamma\cdot\alpha^{2}}{\beta\cdot\sqrt{d}}.$

Therefore, we can pick $\beta$ small enough so that
$\mathbf{x}^{(0)}\odot\mathbf{s}^{(0)}\approx_{2\widehat{\varepsilon}}\mu^{(0)}\bm{\tau}^{(0)}$,
where

$\bm{\tau}^{(0)}\approx_{\gamma/4}\sigma((\mathbf{S}^{(0)})^{-1/2-\alpha}(\mathbf{X}^{(0)})^{1/2-\alpha}).$

Thus after Step 9 of LABEL:alg:ipm-init, we have
$\mathbf{x}^{(\text{end})}\odot\mathbf{s}^{(\text{end})}\approx_{\widehat{\varepsilon}}\mu^{(\text{end})}\cdot\bm{\tau}^{(\text{end})}$.

We finally need to bound the condition number of the matrix
$\overline{\mathbf{A}}$. To do this, we use the following result.

######  Lemma 4.8 (Lemma 5.15 of [GPV23]).

Let $\mathbf{A}\in\mathbb{R}^{n\times d}$, $n>d$, be a matrix with full column
rank. Moreover let $\mathbf{g}\in\mathbb{R}^{n}$. Suppose $\kappa>1$, and

$\|\mathbf{A}^{\top}\mathbf{A}\|_{\mathrm{F}},\|(\mathbf{A}^{\top}\mathbf{A})^{-1}\|_{\mathrm{F}},\|\mathbf{g}\|_{2},1/\|(\mathbf{I}-\mathbf{A}(\mathbf{A}^{\top}\mathbf{A})^{-1}\mathbf{A}^{\top})\mathbf{g}\|_{2}\leq\kappa.$

Then
$\|\overline{\mathbf{A}}^{\top}\overline{\mathbf{A}}\|_{\mathrm{F}},\|(\overline{\mathbf{A}}^{\top}\overline{\mathbf{A}})^{-1}\|_{\mathrm{F}}\leq
8\kappa^{7}$, where
$\overline{\mathbf{A}}=\begin{bmatrix}\mathbf{A}|\mathbf{g}\end{bmatrix}.$

Note that setting

$\widetilde{\mathbf{A}}=\begin{bmatrix}\mathbf{A}\\\ 0\\\
\frac{1}{R}\mathbf{b}^{\top}-\bm{1}^{\top}\mathbf{A}\end{bmatrix},$

since
$\widetilde{\mathbf{A}}^{\top}\widetilde{\mathbf{A}}=\mathbf{A}^{\top}\mathbf{A}+(\frac{1}{R}\mathbf{b}^{\top}-\bm{1}^{\top}\mathbf{A})^{\top}(\frac{1}{R}\mathbf{b}^{\top}-\bm{1}^{\top}\mathbf{A})$,
the condition number of $\widetilde{\mathbf{A}}$ is smaller than the condition
number of $\mathbf{A}$. Now by Lemma 4.8, we have

$\|\overline{\mathbf{A}}^{\top}\overline{\mathbf{A}}\|_{\mathrm{F}},\|(\overline{\mathbf{A}}^{\top}\overline{\mathbf{A}})^{-1}\|_{\mathrm{F}}\leq
8(\kappa(\mathbf{A})+2\|\mathbf{A}\|_{\mathrm{F}})^{7}=O((\kappa(\mathbf{A})+2^{L}\sqrt{nd})^{7}).$

###  4.3 Inverse Maintenance for IPM

In this section, we present the subprocedure IPM of LABEL:alg:ipm-init. This
is presented in LABEL:alg:lin-prog-coordinator-poly-cond. Essentially, in each
iteration, in Step 4b, each machine computes its own part of the gradient of
the potential at a specific point and sends it to the coordinator. The
coordinator then in Step 4c sums the vectors obtained from the machines and
pre-multiplies this sum by a matrix $\mathbf{H}^{-1}$ that spectrally
approximates $(\mathbf{A}^{\top}\mathbf{X}\mathbf{S}^{-1}\mathbf{A})^{-1}$ and
sends the result to all the machines. In Step 4d, each machine uses this
vector to update its own part of the primal and slack vectors. As discussed in
Lemma 4.2, we do not need exact feasibility, and therefore, the vectors sent
to the coordinator or to the machine can be rounded down. Essentially, for
each entry, we only need to send $L+\log(\kappa R/r)$ of its bits to guarantee
convergence (see [BLSS20, Theorem 32]).

The only remaining part for bounding the communication complexity of the
algorithm is to show that Step 4f of LABEL:alg:lin-reg-coordinator-poly-cond
does not resample too many rows over the course of the algorithm.

To improve the running time of the IPMs based on the Lee-Sidford barrier,
[LS15] introduced the following $\sigma$-stability property. As discussed in
Lemma 21 of [BLSS20], the primal and slack vectors in LABEL:alg:lin-reg-
coordinator-poly-cond satisfy this property. This then allows us to bound the
number of resampled rows in Step 4f of LABEL:alg:lin-reg-coordinator-poly-
cond.

######  Definition 4.9 ($\sigma$-stability assumption).

We say that the inverse maintenance problem satisfies the $\sigma$-stability
assumption if for each $k\in[r]$ (where $r$ is the number of rounds/iterations
of the algorithm), we have

$\|\log(\mathbf{d}^{(k)})-\log(\mathbf{d}^{(k-1)})\|_{\sigma(\mathbf{D}^{(k)}\mathbf{A})}\leq
0.1,$ $\|\log(\mathbf{d}^{(k)})-\log(\mathbf{d}^{(k-1)})\|_{\infty}\leq 0.1,$

and

$\beta^{-1}\mathbf{A}^{\top}\mathbf{D}^{(0)}\mathbf{A}\preceq\mathbf{A}^{\top}\mathbf{D}^{(k)}\mathbf{A}\preceq\beta\mathbf{A}^{\top}\mathbf{D}^{(0)}\mathbf{A},$

for $\beta=\operatorname{poly}(n)$.

Input. A matrix $\mathbf{A}:=[\mathbf{A}^{(i)}]\in\mathbb{R}^{n\times d}$,
vector $\mathbf{b}\in\mathbb{R}^{d}$, and vector
$\mathbf{c}:=[\mathbf{c}^{(i)}]\in\mathbb{R}^{n}$, where the $i^{\mathrm{th}}$
machine stores matrix $\mathbf{A}^{(i)}\in\mathbb{R}^{n_{i}\times d}$; initial
feasible primal, dual, slack and weight vectors $\mathbf{x}^{(0)}$,
$\mathbf{y}^{(0)}$, $\mathbf{s}^{(0)}$, $\bm{\tau}^{(0)}$, respectively.
Initial and final centrality parameter $\mu^{(0)}$ and $\mu^{(\text{end})}$.

Output. Vector $\widehat{\mathbf{x}}\in\mathbb{R}^{d}$.

procedure
$\textsc{IPM}(\mathbf{A},\mathbf{b},\mathbf{c},\mathbf{x}^{(0)},\mathbf{y}^{(0)},\mathbf{s}^{(0)},\bm{\tau}^{(0)},\mu^{(0)},\mu^{(\text{end})})$:

  1. 1

Coordinator sets $\alpha=\tfrac{1}{4\log(4n/d)}$,
$\lambda=\tfrac{2}{\varepsilon}\log(\frac{2^{16}n\sqrt{d}}{\alpha^{2}})$,
$\gamma=\min(\frac{\varepsilon}{4},\frac{\alpha}{50\lambda})$,
$\mu=\mu^{(0)}$.

  2. 2

Each machine sets its components of the vectors as follows:
$\overline{\mathbf{x}}=\mathbf{x}^{(0)}$,
$\overline{\mathbf{s}}=\mathbf{s}^{(0)}$,
$\overline{\bm{\tau}}=\bm{\tau}^{(0)}$,
$\mathbf{x}^{(\text{tmp})}=\mathbf{x}^{(0)}$,
$\mathbf{s}^{(\text{tmp})}=\mathbf{s}^{(0)}$,
$\bm{\tau}^{(\text{tmp})}=\bm{\tau}^{(0)}$,
$\overline{\mathbf{w}}=\overline{\mathbf{X}}\overline{\mathbf{s}}$.

  3. 3

The coordinator computes $\mathbf{H}\in\mathbb{R}^{d\times d}$ with
$\mathbf{H}\approx_{\tfrac{\beta\widehat{\varepsilon}}{d^{1/4}\log^{3}n}}\mathbf{A}^{\top}\tfrac{\overline{\mathbf{X}}}{\overline{\mathbf{S}}}\mathbf{A}$
and $\mathbf{H}^{-1}$.

  4. 4

while true do

    1. (a)

Let
$\Phi(\mathbf{v}):=\exp(\lambda(\mathbf{v}-1))+\exp(-\lambda(\mathbf{v}-1))$
for $\mathbf{v}\in\mathbb{R}^{n}$.

    2. (b)

Each machine $i$ computes
${(\mathbf{S}^{(i)})}^{-1}\delta_{\widetilde{\mu}}^{(i)}\in\mathbb{R}^{n_{i}}$
and
${(\mathbf{A}^{(i)})}^{\top}{(\mathbf{S}^{(i)})}^{-1}\delta_{\widetilde{\mu}}^{(i)}\in\mathbb{R}^{d}$,
where
$\delta_{\widetilde{\mu}}\approx\nabla\Phi(\overline{\mathbf{W}}^{-1}\widetilde{\bm{\mu}})$,
and sends it to the coordinator.

    3. (c)

The coordinator computes
$\mathbf{H}^{-1}\sum_{i=1}^{s}{\mathbf{A}^{(i)}}^{\top}{(\mathbf{S}^{(i)})}^{-1}\delta_{\widetilde{\mu}}^{(i)}\in\mathbb{R}^{d}$
and sends it to all of the machines.

    4. (d)

Each machine $i$ computes
$\frac{\overline{\mathbf{X}}^{(i)}}{\overline{\mathbf{S}}^{(i)}}\mathbf{A}^{(i)}\mathbf{H}^{-1}{\mathbf{A}}^{\top}{\overline{\mathbf{S}}}^{-1}\delta_{\widetilde{\mu}}\in\mathbb{R}^{n_{i}}$
and uses it to compute $\mathbf{\bm{\delta}_{x}}^{(i)}$ and
$\mathbf{\bm{\delta}_{s}}^{(i)}$ to update, respectively,
$\overline{\mathbf{x}}^{(i)}$ and $\overline{\mathbf{s}}^{(i)}$.

    5. (e)

Each machine sets its portion of the vectors as the following:
$\mathbf{x}^{(\text{tmp})}_{i}=\overline{\mathbf{x}}_{i}$ if
$\overline{\mathbf{x}}_{i}\approx_{\gamma/8}\mathbf{x}^{(\text{tmp})}_{i}$;
$\mathbf{s}^{(\text{tmp})}_{i}=\overline{\mathbf{s}}_{i}$ if
$\overline{\mathbf{s}}_{i}\approx_{\gamma/8}\mathbf{s}^{(\text{tmp})}_{i}$;
$\bm{\tau}^{(\text{tmp})}_{i}=\overline{\bm{\tau}}_{i}$ if
$\overline{\bm{\tau}}_{i}\approx_{\gamma/8}\bm{\tau}^{(\text{tmp})}_{i}$. Set
$\overline{\mathbf{w}}=\overline{\mathbf{X}}\overline{\mathbf{s}}$ and
$\overline{\mathbf{v}}=\mu\overline{\mathbf{W}}^{-1}\overline{\bm{\tau}}$.

    6. (f)

For any $i$ where either of
$\mathbf{x}^{(\text{tmp})}_{i},\mathbf{s}^{(\text{tmp})}_{i}$ or
$\bm{\tau}^{(\text{tmp})}_{i}$ has changed, we resample the $i$’th row
according to its leverage scores and send it to the coordinator with the
corresponding probability. The coordinator updates $\mathbf{H}^{-1}$
accordingly and sends a sketch of it to the machines.

    7. (g)

The machines use the sketch to update their leverage scores.

    8. (h)

If $\mu>\mu^{(\text{end})}$, set
$\mu=\max\\{\mu^{(\text{end})},(1-\frac{\gamma\alpha}{2^{15}\sqrt{d}})\mu\\}$.
Otherwise, set
$\mu=\max\\{\mu^{(\text{end})},(1+\frac{\gamma\alpha}{2^{15}\sqrt{d}})\mu\\}$.

    9. (i)

If $\mu=\mu^{(\text{end})}$ and
$\Phi(\overline{\mathbf{v}})\leq\frac{2^{16}n\sqrt{d}}{\alpha^{2}}$, break.

  5. 5

Return
$(\overline{\mathbf{x}},\overline{\mathbf{y}},\overline{\mathbf{s}},\overline{\bm{\tau}},\mu)$

Algorithm 11: Protocol for linear programming in the coordinator setting

alg]alg:lin-prog-coordinator-poly-cond

We use LABEL:alg:inv-maintenance for our inverse maintenance (i.e., to
maintain a spectral approximation of
$(\mathbf{A}^{\top}\mathbf{X}\mathbf{S}^{-1}\mathbf{A})^{-1}$ in the
coordinator). This algorithm is inspired by Algorithm 3 of [LS15] and based on
the following which is [LS15, Lemma 15], the number of changes in
LABEL:alg:inv-maintenance is bounded by
$\widetilde{O}(d\log(\varepsilon^{-1}))$ when used with LABEL:alg:lin-reg-
coordinator-poly-cond since the number of iterations of LABEL:alg:lin-reg-
coordinator-poly-cond is
$\widetilde{O}(\sqrt{d}\cdot\log^{2}(\varepsilon^{-1}))$.

######  Lemma 4.10 ([LS15]).

Suppose changes of $\mathbf{d}$ and the error occurred in computing leverage
scores is independent of the sampled matrix. Under the $\sigma$ stability
guarantee, during the first $r$ iterations of LABEL:alg:inv-maintenance, the
expected number of coordinate changes in $\mathbf{H}^{(k)}$ over all
iterations $k\in[r]$ is $O(r^{2}\log(d))$.

###  4.4 Proof of Main Result on High-Accuracy Linear Programming

We are now equipped to prove the main theorem for the communication complexity
of linear programming in the point-to-point model of communication. The
correctness and number of iterations of the algorithm follow from [BLSS20].
Therefore, we focus on discussing only the communication complexity bounds.

###### Proof of Theorem 1.10.

First note that, the only parts of LABEL:alg:ipm-init that requires
communication are Step 4, Step 5, Step 7. By Lemma 4.7, Step 4 only requires
$\widetilde{O}(d^{2}L+sd(L+\log\kappa))$ bits of communication.

Step 5 and Step 7 of LABEL:alg:ipm-init both call the IPM procedure of
LABEL:alg:lin-prog-coordinator-poly-cond. The only parts of this algorithm
with communication are Step 3, Step 4b, Step 4c, Step 4f, Step 4i. By Lemma
3.2, the communication cost of Step 3 is
$\widetilde{O}(d^{2}L+sd(L+\log(R\kappa/r)))$. For each iteration, the
communication cost of Step 4b and Step 4c is
$\widetilde{O}(sd(L+\log(R\kappa/(r\varepsilon))))$. Therefore since there are
$\widetilde{O}(\sqrt{d}\log(\varepsilon^{-1}))$ iterations, the total cost of
these steps is
$\widetilde{O}(sd^{1.5}(L+\log(R\kappa/(r\varepsilon)))\log(\varepsilon^{-1}))$.
By Lemma 4.10 and [BLSS20, Lemma 21], the total communication cost of Step 4f
is $\widetilde{O}(d^{2}L\log^{2}(\varepsilon^{-1}))$ because there are at most
a total of $\widetilde{O}(d\log^{2}(\varepsilon^{-1}))$ rows that need to be
sent to the coordinator over the course of the algorithm. For Step 4i, note
that each machine needs to compute a number and send it to the coordinator, so
the coordinator be able to compute $\Phi(\overline{\mathbf{v}})$. Since there
are $\widetilde{O}(\sqrt{d}\log(\varepsilon^{-1}))$ iterations, the total cost
of this over the whole course of the algorithm is
$\widetilde{O}(s\sqrt{d}(L+\log(\kappa
Rr^{-1}\varepsilon^{-1}))\cdot\log(\varepsilon^{-1}))$. ∎

##  5 Finite-Sum Minimization in the Blackboard Model

In this section, we study finite-sum minimization (in the distributed
setting), i.e.,

$\text{minimize}_{\mathbf{x}\in\mathbb{R}^{d}}\sum_{i=1}^{s}f_{i}(\mathbf{x}),$

where each $f_{i}$ is convex, Lipschitz, and supported on only $d_{i}$
coordinates of $\mathbf{x}$. While finite-sum minimization itself is a general
problem class encompassing, for example, empirical risk minimization, with
each $f_{i}$ measuring the loss incurred by the $i^{\mathrm{th}}$ data point
from the training set, the additional structural assumption (of dependence on
$d_{i}$ coordinates) is also seen in prominent problem classes like
decomposable submodular function minimization [AKMSV21]. There exist numerous
fast variants of stochastic gradient descent [RM51, BC03, Zha04, Bot12] for
Section 5 such as [RSB12, SZ13, JZ13, MZJ13, DBL14, Mai15, AY16, HL16, SLB17]
but most of these algorithms depend on the problem’s condition number, which
could be quite large (and hence undesirable) for non-smooth $f_{i}$.
Conversely, both the cutting-plane methods [LSW15] and robust interior-point
methods [LV21] exchange their reliance on condition number for worse
dependencies on the problem dimension.

In this section, our goal is to solve this problem with efficient
communication complexity. We now formally state our problem setup and all the
main results of this section.

See 1.13

In order to obtain our communication bound, we first derive the following
fine-grained cost (in terms of certain weights), which also yields improved
rates for submodular function minimization.

See 1.14

See 1.15

###  5.1 An Overview of Our Algorithm

The goal of this section is to prove Theorem 1.13. We obtain this result via
LABEL:alg:min-sum-convex-blackboard obtained by generalizing a technique
introduced in [DJLPY22] and then setting it in the distributed framework.

Every machine holds a copy of all the data (i.e., the current variable
$\mathbf{x}$ , and each machine $i$ holds the corresponding function’s
subgradient oracle $\mathcal{O}_{i}$. Per the technique of [DJLPY22], we first
use the standard epigraph trick to reduce Section 5 to the following
formulation,

$\begin{array}[]{ll}\mbox{minimize}&\mathbf{c}^{\top}\mathbf{x},\\\
\mbox{subject
to}&\mathbf{x}_{\textrm{i}}\in\mathcal{K}_{i}\subseteq\mathbb{R}^{d_{i}+1}\;\forall
i\in[s]\\\ &\mathbf{A}\mathbf{x}=\mathbf{b}.\end{array}$

where $\mathbf{x}=[\mathbf{x}_{i}]$ concatenates the $s$ vectors
$\mathbf{x}_{i}\in\mathbb{R}^{d_{i}}$, and all the sets $\mathcal{K}_{i}$ are
disjoint and convex. Section 5.1 transfers the overlap in support between the
original $f_{i}$’s into $\mathbf{A}\mathbf{x}=\mathbf{b}$. This reduction
requires only the knowledge of support of each $f_{i}$, so this reduction can
be done using $O(\sum_{i=1}^{s}d_{i}L)$ bits of communications. After the
reduction, all the machines hold all the data — vectors
$\mathbf{c},\mathbf{x},\mathbf{b}$, matrix $\mathbf{A}$ — and the
$i^{\mathrm{th}}$ machine holds the separation oracle $\mathcal{O}_{i}$ for
the $i^{\mathrm{th}}$ set $\mathcal{K}_{i}$. This oracle is the only means to
access the sets $\mathcal{K}_{i}$ and is obtained via a reduction from the
corresponding subgradient oracle for $f_{i}$, see [LSV18]. Specifically, at
any queried point $\mathbf{z}_{i}$, the oracle either asserts
$\mathbf{z}_{i}\in\mathcal{K}_{i}$, or returns a separating hyperplane that
separates $\mathbf{z}_{i}$ from $\mathcal{K}_{i}$. Formally, we prove the
theorem below.

######  Theorem 5.1 (Main theorem for Section 5.1).

Consider the convex program described in Section 5.1, with every machine
holding all the data and the $i^{\mathrm{th}}$ machine having a separation
oracle access to $\mathcal{K}_{i}$. Denote
$\mathcal{K}=\mathcal{K}_{1}\times\mathcal{K}_{2}\times\dotsc\times\mathcal{K}_{s}$.
Assuming we have

  * outer radius $R$: For any $\mathbf{x}_{i}\in\mathcal{K}_{i}$, we have $\|\mathbf{x}_{i}\|_{2}\leq R$, and

  * inner radius $r$: There exists a $\mathbf{z}\in\mathbb{R}^{d}$ such that $\mathbf{A}\mathbf{z}=\mathbf{b}$ and $\mathcal{B}(\mathbf{z},r)\subset\mathcal{K}$,

then, for any $0<\varepsilon<\frac{1}{2}$, we can find a point
$\mathbf{x}\in\mathcal{K}$ satisfying $\mathbf{A}\mathbf{x}=\mathbf{b}$ and

$\mathbf{c}^{\top}\mathbf{x}\leq\min_{\begin{subarray}{c}\mathbf{x}_{\textrm{i}}\in\mathcal{K}_{i}\subseteq\mathbb{R}^{d_{i}+1}\forall
i\in[s],\\\
\mathbf{A}\mathbf{x}=\mathbf{b}\end{subarray}}\mathbf{c}^{\top}\mathbf{x}+\varepsilon\cdot\|\mathbf{c}\|_{2}\cdot
R,$