# Faster Dynamic Matrix Inverse for Faster LPs

Shunhua Jiang Zhao Song Omri Weinstein Hengjie Zhang<EMAIL_ADDRESS>Columbia University. Research supported by NSF CAREER award
<EMAIL_ADDRESS>Princeton University and Institute for Advanced
<EMAIL_ADDRESS>Columbia University. Research supported by NSF
CAREER award<EMAIL_ADDRESS>Columbia University. Research
supported by NSF CAREER award CCF-1844887.

Motivated by recent Linear Programming solvers, we design dynamic data
structures for maintaining the inverse of an $n\times n$ real matrix under
_low-rank_ updates, with polynomially faster amortized running time. Our data
structure is based on a recursive application of the Woodbury-Morrison
identity for implementing _cascading_ low-rank updates, combined with recent
sketching technology. Our techniques and amortized analysis of multi-level
partial updates, may be of broader interest to dynamic matrix problems.

This data structure leads to the fastest known LP solver for general (dense)
linear programs, improving the running time of the recent algorithms of (Cohen
et al.’19, Lee et al.’19, Brand’20) from
$O^{*}(n^{2+\max\\{\frac{1}{6},\omega-2,\frac{1-\alpha}{2}\\}})$ to
$O^{*}(n^{2+\max\\{\frac{1}{18},\omega-2,\frac{1-\alpha}{2}\\}})$, where
$\omega$ and $\alpha$ are the fast matrix multiplication exponent and its
dual. Hence, under the common belief that $\omega\approx 2$ and $\alpha\approx
1$, our LP solver runs in $O^{*}(n^{2.055})$ time instead of
$O^{*}(n^{2.16})$.

##  1 Introduction

_Dynamic matrix inverse problems_ ask to maintain the inverse of an $n\times
n$ matrix $M$ over some field (say $\mathbb{R}$), when $M$ undergoes a long
sequence of row/column updates. This data structure problem arises in many
important TCS applications, such as directed reachability in dynamic graphs
and maintaining the eigenvalues and rank of a matrix [San04], empirical risk
minimization [LSZ19], routing and electrical-flow computation (inverting
Laplacians) [ST04, Mad13], and in fastest known linear programming solvers
[LS15, CLS19, BLSS20]. While many variants of this problem have been
considered over the years (depending on the specific application), the most
common one requires the data structure to support _vector-queries_ , i.e.,
computing $\textsc{Query}(h)=M^{-1}h$ where $h$ comes from some
restricted111Indeed, supporting _arbitrary_ online queries
$h\in\mathbb{R}^{n}$ is conjectured to be impossible in truly sub-quadratic
time even in the _static_ case where $M$ remains fixed, see the “oMV
Conjecture” [HKNS15]. family of $\mathbb{R}^{n}$, under a long sequence of
_low-rank_ updates (e.g., row/column changes to $M$). Recomputing the inverse
from scratch upon each update in this setting incurs a daunting computational
overhead, and therefore the goal is to optimize the tradeoff between the
(amortized) update time $t_{u}$ and query time $t_{q}$ (or the total running
time). The first dynamic data structure for this problem, tracing back to the
1950’s [Woo49, Woo50] shows how to _explicitly_ maintain the inverse of $M$ in
$O(n^{2})$ time, which is already nontrivial. Substantial improvements to this
data structure were developed more recently, showing that _row-queries_ (which
are equivalent to answering $(M^{\top})^{-1}h$ for any _$1$ -sparse_ vector
$h$) under row/column updates, can be done in
$\max\\{t_{u},t_{q}\\}=O(n^{1.529})$ time [San04, SM10, BNS19]. Brand et al.
[BNS19] also showed a certain conditional $\Omega(n^{1.5})$ worst-case lower
bound for the _exact_ version of this problem, which is important for some of
the aforementioned applications (e.g., maintaining graph properties such as
directed reachability).

In this work we consider a more challenging dynamic inverse problem: The data
structure needs to support low-rank updates of the form
$\textsc{Update}(u,v)=M+uv^{\top}$ where $u,v$ come from some fixed set of
vectors of size $O(n)$, and needs to answer inverse queries $M^{-1}h^{(i)}$
under a _slowly changing_ vector sequence ($\|h^{(i)}-h^{(i-1)}\|_{0}\leq
O(1)$). In contrast to row-queries (a special case $\|h^{(i)}\|_{0}=1$), here
the online query vectors $h$ may be arbitrarily dense, but their _differences_
are sparse.222 Note that sparsity of $\Delta h$ is _not_ equivalent to
sparsity of queries $h$ themselves: If $M$ were _fixed_ , then by linearity,
computing $M^{-1}(\Delta h)$ would indeed suffice, but here $M$ is dynamically
changing so this standard trick doesn’t work. An important special case of
this dynamic problem is “projection maintenance" [CLS19], where updates to $M$
take the form $\textsc{Update}(D)=M+ADA^{\top}$ where $A$ is an arbitrary
fixed matrix and $D$ is a sparse diagonal matrix (hence
$\operatorname{rank}(ADA^{\top})\leq\operatorname{rank}(D)$ is small and can
be therefore written as $\sum_{i=1}^{rank(D)}A_{i_{1}}A_{i_{2}}^{\top}$).

Nevertheless, an important difference between this work and the aforementioned
ones on dynamic inverse maintenance (e.g., [San04, SM10, BNS19]), is allowing
_approximate_ answers, i.e., the data structure needs to compute $M^{-1}h$ up
to some small relative $\ell_{\infty}$ error. This enables the use of
randomized tools from _sketching and sparse recovery_ literature [GLPS10,
CW13, LNNT16, LSZ19]. Another point of departure is our use of very heavy
_amortization_ , augmenting the approach of [CLS19, LSZ19, Bra20] with a novel
algebraic technique and more sophisticated potential analysis (based on “high
order” martingales). A recurring theme in dynamic data structures (both upper
and lower bounds) is that amortized analysis is a different ballgame compared
to worst-case performance – Some classic examples are the amortized analysis
of fully-dynamic undirected connectivity [ST85] and its matching amortized
cell-probe lower bound [PD06], which required substantially new techniques
(and another decade) compared to the worst-case lower bound [FS89]. In
contrast to dynamic graph problems, whose amortized complexity has been
studied extensively, its counterpart in _dynamic matrix problems_ is far less
understood ([HKNS15, BNS19]), and we believe our work sheds further light on
the power of amortization.

##### Dynamic inverse in linear programming

The primary application and motivation of this work is the role of dynamic
matrix inverse data structures in speeding up _interior-point methods_ (IPMs)
for solving linear programs (LPs) in close to matrix-multiplication time.

Linear programming is one of the cornerstones of algorithm design and convex
optimization, dating back to as early as Fourier in 1827. LPs are the key
toolbox for (literally hundreds of) approximation algorithms, and a standard
subroutine in convex optimization problems. Dantzig’s 1947 _simplex algorithm_
[Dan47] was the first proposed solution for general LPs with $n$ variables and
$d$ constraints ($\min_{Ax=b,x\geq 0}c^{\top}x$). Despite its impressive
performance in practice, however, the simplex algorithm turned out to have
exponential worst-case running time (Klee and Minty [KM72]). The first
_polynomial time_ algorithm for general LPs was only developed in 1980, when
Khachiyan [Kha80] introduced the _Elliposid method_ , and showed that it runs
in $O(n^{6})$ time. Unfortunately, this algorithm is very slow in practice
compared to the simplex algorithm, raising a quest for LP solvers which are
efficient in both theory and practice.

This was the primary motivation behind the development of _interior point
methods_ (IPMs), which uses a primal-dual gradient descent approach to
iteratively converge to a feasible solution (Karmarkar, [Kar84]). An appealing
feature of IPM methods for solving LPs is that they are not only guaranteed to
run fast in theory, but also in practice [Str87]. In 1989, Vaidya proposed an
$O(n^{2.5})$ LP solver based on a specific implementation of IPMs, known as
the _Central Path algorithm_ [Vai87, Vai89]. Vaidya already observed that the
main bottleneck of this algorithm boils down to a dynamic data structure
problem of maintaining the _inverse matrix_ $M$ associated with the central
path equations (see Section 3), under a sequence of updates of the form
$(M+A\Delta A^{\top})^{-1}$, where $\Delta$ is a sparse diagonal matrix and
$A$ is the fixed LP constraint matrix. Since $\Delta$ is sparse, each
“gradient descent” iteration of this algorithm induces a _low-rank_ update to
$M$, hence it is conceivable to avoid recomputing the inverse matrix from
scratch—which would naively cost $n^{\omega}$ time per iteration—and gain
substantially from amortization. This data structure problem was the
centerpiece of the recent line of developments on IPM solvers [CLS19, LSZ19,
Bra20], which focused on designing faster dynamic inverse structures for
implementing the Central Path algorithm.

The fastest known algorithm for general (dense) LPs, based on this approach,
is due to Cohen, Lee and Song [CLS19], whose running time is

$O^{*}(n^{\omega}+n^{2.5-\alpha/2}+n^{2+1/6}),$

where $\omega<2.37$ is the fast matrix-multiplication exponent and
$\alpha>0.31$ is the dual matrix multiplication exponent.333The dual exponent
$\alpha$ is defined as the asymptotically maximum number $a\leq 1$ s.t
multiplying an $n\times n^{a}$ matrix by an $n^{a}\times n$ matrix can be done
is $n^{2+o(1)}$ time. The current best lower bound is $\alpha>0.31389$ [GU18].
Note that $n^{\omega}$ is the minimal time for merely inverting a matrix,
i.e., finding _any_ feasible solution ($Ax=b$) to the LP, hence it seems quite
remarkable that solving the full optimization problem ($\min_{Ax=b,x\geq
0}c^{\top}x$) may be done at virtually no extra cost. It is widely believed
that $\omega\approx 2$ [CKSU05, Wil12] and as such, $\alpha\approx 1$ (though
the only formal connection between these constants is
$\omega+(\omega/2)\alpha\leq 3$ [CGLZ20]). Assuming indeed that
$\omega<2+1/6\approx 2.166$ and $\alpha>0.66$, the runtime of the
aforementioned algorithms is $n^{2.166}$. Whether the additive $n^{2.166}$
term can be improved or completely removed was explicitly posed as an open
question in [CLS19] and [Son19].

Our main result is an affirmative answer to this open question, asserting that
LPs can be solved in matrix multiplication time for nearly any value of
$\omega$ (i.e., so long as $\omega>2.055$). We design an improved LP solver
which runs in time $O^{*}(n^{\omega}+n^{2.5-\alpha/2}+n^{2+1/18})$. In the
most notable (ideal) case that $\omega\approx 2$ and $\alpha\approx 1$, our
algorithm runs in $O^{*}(n^{2.055})$ time, instead of $O^{*}(n^{2.166})$ time
of previous IPM algorithms [CLS19, LSZ19, Bra20]. More precisely:

######  Theorem 1.1 (Main result, Informal statement of Theorem 4.1).

Let $\min_{Ax=b,x\geq 0}c^{\top}x$ be a linear program where
$A\in\mathbb{R}^{d\times n}$ and $d=\Omega(n)$. Then for any accuracy
parameter $\delta\in(0,1)$, there is a randomized algorithm that solves the LP
in expected time

$\displaystyle
O^{*}(n^{\omega}+n^{2.5-\alpha/2}+n^{2+1/18})\cdot\log(n/\delta).$

We achieve this result by designing a more efficient projection maintenance
data structure, speeding up both the update and query times of previous
algorithms. This is done via a new algebraic framework for bootstrapping lazy
updates (described next), combined with randomized compression techniques and
a sophisticated amortized analysis of the underlying dynamic process.

##### Organization

In Section 2 we provide a high-level description of our main technique, which
is the centerpiece of this paper. Section 3 contains some necessary background
and brief overview of previous related work. In Section 4, we provide a
detailed technical overview of the proof of Theorem 1.1. This 10-page
streamlined overview should be understood as an extended abstract of our
entire result, deferring technical proofs and calculations to the Appendix.

_Appendix Organization._ Section A contains preliminaries and notation. In
Section B we provide an analysis of the Stochastic Central Path algorithm,
postponing output-feasibility issues to Section I. We put preliminary part of
our data structure in Section C. We present our full “cascading data
structure” and prove its correctness in Section D, and we analyze its running
time in Sections E, F. Finally, Section G contains the final runtime analysis
of our LP solver using the full data structure. In Section H, we present more
details for Section 2.

##  2 Bootstrapping low-rank updates via cascading lazy updates

One of the main new ideas of this work is “bootstrapping”444This term refers
to a general approach for speeding up dynamic algorithms by repeating a
certain technique several times recursively [San04, BNS19]. We remark that
[BNS19] uses a different kind of bootstrapping to speed up exact inverse
maintenance under simpler row-updates and row-queries (via “one more leve” of
FMM). In particular, [BNS19] has nothing to do with _dynamic LU-
decompositions_ nor recursion on Woodburry’s identity. Nevertheless, it is
noteworthy that the bottleneck in both works is maintaining certain matrix
products for $>2$ “levels” of recursion.  the lazy updates technique of
[CLS19] by repeated application of the Woodburry-Morrison identity (Lemma
A.2), allowing faster low-rank operations via _cascading lazy updates_. For
ease of presentation, let us focus on the following simplified dynamic data
structure problem, which we henceforth call _low-rank inverse maintenance_
problem (Definition. H.5). The data structure is initially given a full rank
matrix $M^{(0)}=M\in\mathbb{R}^{n\times n}$, and a fixed vector
$h\in\mathbb{R}^{n}$. The data structure needs to support a sequence of rank-
one updates and vector-queries as follows. In the $t$-th Update operation, we
are given two vectors $u^{(t)},v^{(t)}\in\mathbb{R}^{n}$ and a real number
$c^{(t)}\in\mathbb{R}$, and need to perform a rank-1 update
$M^{(t)}=M^{(t-1)}+c^{(t)}\cdot u^{(t)}(v^{(t)})^{\top}$. A Query operation
asks to calculate $x=(M^{(T)})^{-1}\cdot h$, where $T$ is the number of
updates in the sequence so far. Note that this setting easily captures
rank-$k$ updates invoking $k$ consecutive rank-1 updates.

Achieving sub-quadratic update and query times for this problem requires some
restriction on the update vectors $u^{(t)},u^{(t)}$ (we show in Section H that
otherwise it would break the oMV Conjecture [HKNS15]). Our technique requires
one reasonable assumption, namely, that all updates $u^{(t)}$ and $v^{(t)}$
come from a fixed set $|S|=O(n)$. As noted in the introduction, LP projection
maintenance is a special case where $M=(ADA^{\top})$, $A$ is the fixed LP
matrix and $D$ is a diagonal matrix which is changing slowly under $\ell_{0}$
norm. Thus, a sparse update $\Delta_{i}$ to $D$ corresponds to $M\leftarrow
M+\Delta_{i}\cdot A_{i}A_{i}^{\top}$.

We also remark that the assumption that the query vector $h$ is fixed
throughout the sequence—while natural in many streaming applications—is only
for simplicity of exposition: Our data structure will actually support an
online sequence of slowly-changing queries $h$ (i.e.,
$\|h^{(t)}-h^{(t-1)}\|_{0}=o(n)$). Handling sparse updates to $h$ turns out to
be much easier than low-rank updates to $M$, hence we focus on the latter
task.

Our technique for solving the _low-rank inverse maintenance_ problem is based
on an algorithmic generalization of Woodbury’s identity to $K>1$ “levels” (to
be explained below), allowing for _recursive lazy updates_ of these $K$ levels
using different thresholds. The basic idea is to (dynamically) group updates
into $K$ “epochs” $0\leq t_{1}\leq\cdots\leq t_{K-1}\leq t_{K}=T$. The first
“level” maintains the first epoch $t_{1}$ and $M^{(t_{1})}$. Similarly, in
level $k\in\\{2,\cdots,K\\}$ we group all updates $u^{(t)},v^{(t)},c^{(t)}$
for which $t\in(t_{k-1},t_{k}]$, and partition the update sequence in terms of
these epochs. More formally, let
$U_{k},V_{k}\in\mathbb{R}^{n\times(t_{k}-t_{k-1})}$ and the diagonal matrix
$C_{k}\in\mathbb{R}^{(t_{k}-t_{k-1})\times(t_{k}-t_{k-1})}$ be, respectively,
the concatenation of all $u_{i},v_{i}$ and $c_{i}$ in the $k$th epoch, so that
$\sum_{i=t_{k-1}+1}^{t_{k}}c_{i}\cdot
u_{i}v_{i}^{\top}=U_{k}C_{k}V_{k}^{\top}$. Let $r_{k}$ be defined as the rank
of $U_{k}C_{k}V_{k}^{\top}$, the epochs are maintained under the invariant
that $r_{k}\leq n_{k}$, where $n=n_{1}\gg n_{2}\gg\cdots\gg n_{K}\geq 1$ are
predefined thresholds that decrease exponentially and can later be optimized.
In this terminology, for any $h\in\mathbb{R}^{n}$, the query answer

$\displaystyle
x:=\left(M^{(t_{1})}+\sum_{i=t_{1}+1}^{T}c_{i}u_{i}v_{i}^{\top}\right)^{-1}h$

can be equivalently re-written as the following linear system (by adding
‘dummy’ variables $\xi_{j}$):

$\displaystyle\begin{bmatrix}x\\\ \xi_{2}\\\ \xi_{3}\\\ \xi_{4}\\\ \vdots\\\
\xi_{K}\end{bmatrix}=\left(\underbrace{\begin{bmatrix}M^{(t_{1})}&U_{2}&U_{3}&U_{4}&\cdots&U_{K}\\\
V_{2}^{\top}&-C_{2}^{-1}&0&0&\cdots&0\\\
V_{3}^{\top}&0&-C_{3}^{-1}&0&\cdots&0\\\
V_{4}^{\top}&0&0&-C_{4}^{-1}&\cdots&0\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\\
V_{K}^{\top}&0&0&0&\cdots&-C_{K}^{-1}\end{bmatrix}}_{D}\right)^{-1}\cdot\begin{bmatrix}h\\\
0\\\ 0\\\ 0\\\ \vdots\\\ 0\end{bmatrix}$ (1)

Note that this equation is precisely Woodbury’s identity, written in a
$K$-block matrix form. Indeed, Woodbury’s identity is derived as the solution
$x$ to the linear system

$\begin{bmatrix}x\\\ \xi\end{bmatrix}=\begin{bmatrix}M^{(t_{1})}&U\\\
V^{\top}&-C^{-1}\end{bmatrix}^{-1}\cdot\begin{bmatrix}h\\\ 0\end{bmatrix},$

where here $\xi$, $U$, $V$, $C$ denote the concatenation of $\xi_{k}$,
$U_{k}$, $V_{k}$, $C_{k}$ respectively, over all $k\in\\{2,\cdots,K\\}$. We
now explain how the above generalization leads to an efficient data structure
for implementing low-rank updates.

##### LU-decomposition

At update time, the data structure will maintain an LU-decomposition (lower-
upper triangular factorization) of the matrix $D$ in (1).

$\displaystyle D=L\cdot U,$ (2)

where $L$ and $U$ are both $K$-block triangular matrices (which are uniquely
defined by the Gaussian-elimination algorithm and imposing the diagonals of
$L$ to be identity matrices, see Eq.(H.1) for an example of the case $K=3$).
As we will explain in the next paragraph, such decomposition is useful since
the inverse of triangular matrices ($L,U$) can be maintained efficiently. This
means that the (slowly changing) query answers
$U^{-1}L^{-1}[h,0]^{\top}=D^{-1}[h,0]^{\top}=[x,\xi]^{\top}$ can be maintained
efficiently with respect to the updated $D$.

##### Cascading lazy updates and query

The key part of our data structure is performing lazy updates recursively
w.r.t different thresholds, to speed up the amortized runtime. Recall that in
the latest $T$-th update, we are given $u_{T}$, $v_{T}$, $c_{T}$. In order to
solve the forthcoming queries using this framework, we need to update the
epoch in the bottom level $t_{K}$ to reflect the up-to-date $M^{(T)}$, by
updating its corresponding tuple $U_{K},V_{K},C_{K}$ and maintain the LU-
decomposition. If the update $t_{K}\leftarrow T$ violates the invariant
$r_{K}\leq n_{K}$ (Recall $r_{k}=t_{k}-t_{k-1}$ is the rank of the update
$U_{k}C_{k}V_{k}^{\top}$ in the $k$-th epoch), we “cascade” the update to the
next level by updating $t_{K-1}\leftarrow t_{K}\leftarrow T$ and also updating
$U_{K-1},V_{K-1},C_{K-1}$ to include $U_{k},V_{K},C_{K}$, and recurse on the
next level $t_{K-2}$ and so on, until the threshold invariant is restored. A
_crucial observation_ in the implementation of cascading lazy updates and
maintaining the LU-decomposition is that when level $k$ gets updated, the
upper-left $(k-1)\times(k-1)$ blocks of $L$ and $U$, which is the dominating
part compared to the entire matrix, _does not change_. Since $L$ and $U$ are
triangular matrices, the upper-left $(k-1)\times(k-1)$ blocks of $L^{-1}$ and
$U^{-1}$ _also_ remain intact. If we write the new $L^{-1},U^{-1}$ as
$(L^{-1})^{\operatorname{new}}=L^{-1}+\Delta L$ and
$(U^{-1})^{\operatorname{new}}=U^{-1}+\Delta U$, the non-zero part of $\Delta
L$ ($\Delta U$) is lower (upper) triangular of $n\times n_{k}$ submatrices
that reside on the bottom (right). Thus, the query answer can be maintained by
computing
$(U^{-1})^{\operatorname{new}}(L^{-1})^{\operatorname{new}}[h,0]^{\top}=(U^{-1}+\Delta
U)(L^{-1}+\Delta L)[h,0]^{\top}$. (The “heavy” part $U^{-1}L^{-1}[h,0]^{\top}$
can be reused, and the remaining components are relatively cheap to
calculate).

This argument implies that, at level $k$, we can afford time proportional to
its size $r_{k}\leq n_{k}$ to rebuild the lower-upper triangular matrix $L$
and $U$ on their changed part along with the vector
$U^{-1}L^{-1}[h,0]^{\top}$, to perform fast queries. We remark that
recomputing $\Delta L,\Delta U$ requires certain preprocessing which is where
we exploit the assumption that updates $u,v$ are from a fixed set.

##### Application to the LP setting.

We now explain how the above framework can be successfully applied to the LP
setting, i.e., to efficiently implement IPM algorithms. As explained in the
next section, the goal of each iteration in IPM solvers is to (re-)calculate
an approximate matrix-vector product $r$ of the following form, given new
disposition vectors $w^{\operatorname{appr}}$ and $h^{\operatorname{appr}}$
(see Algorithm 1):

$r:=P(w^{\operatorname{appr}})\cdot
h^{\operatorname{appr}}=\sqrt{W^{\operatorname{appr}}}A^{\top}(AW^{\operatorname{appr}}A^{\top})^{-1}A\sqrt{W^{\operatorname{appr}}}\cdot
h^{\operatorname{appr}}.$

We can use our cascading lazy updates technique to maintain the middle term
$(AW^{\operatorname{appr}}A^{\top})^{-1}\cdot h$, where for simplicity we
assume here that $h$ is some fixed vector. Indeed, letting
$w^{\operatorname{appr}}$ denote the $j$-th iteration update $w^{(j)}$,
recalculating $(AW^{\operatorname{appr}}A^{\top})^{-1}\cdot h$ corresponds to
maintaining a sequence of $T_{j}:=\|w^{(j)}-w^{(j-1)}\|_{0}$ of updates
$u_{i},v_{i},c_{i}$ followed by one query such that
$\sum_{i=1}^{T_{j}}c_{i}\cdot u_{i}v_{i}^{\top}=A(W^{(j)}-W^{(j-1)})A^{\top}$.

Pictorially, the cascading lazy updates process resembles the following
“chasing game”: Child number $t_{k-1}$ is chasing its friend $t_{k}$. Once the
distance between them is too large, $t_{k-1}$ updates his position to the
position of $t_{k}$ (Figure 1). This process generalizes [CLS19], in which
there’s only one child $t_{1}$ chasing its friend $t_{2}=T$. To analyze the
amortized cost of this process, we exploit the following special features of
the LP problem:

Figure 1: The cascading lazy updates process with per-level invariants
$\|w^{(t_{k})}-w^{(t_{k}-1)}\|_{0}\leq n_{k}$. Updates become more expensive
but less frequent as we move down the levels.

  1. 1.

The updates $u^{(t)},v^{(t)}$ is some row of the original (fixed) LP matrix
$A$.

  2. 2.

$w^{\operatorname{appr}}$ is slowly changing in each IPM iteration, which
imposes nontrivial _sparsity guarantees_ that can be used to determine the
cascading thresholds $\\{n_{k}\\}_{k=1}^{K}$: Informally speaking, since
$w^{\operatorname{appr}}$ is roughly a martingale, the rank $r_{k}$ of epoch
$k$ will be typically far less than its boundary condition ($r_{k}\ll
t_{k}-t_{k-1}$) – It takes about $\sqrt{n_{k}}$ LP-iterations for epoch $k$ to
exceed the threshold $n_{k}$, in which case we need to update epoch $k-1$ and
compute $\Delta L$, $\Delta U$ which are of size $n\times n_{k-1}$. (see
Sections 4.2, F for more details).

  3. 3.

The robustness of the Central Path algorithm implies that queries $M^{-1}\cdot
h$ can tolerate small relative $\ell_{\infty}$ error. This allows the use of
randomized compression (left-sketching $R^{\top}R\cdot
U^{-1}L^{-1}[h,0]^{\top}$) to further reduce the running time (see Section 4.1
for technique overview).

Therefore, assuming ideal matrix-multiplication constants ($\omega=2$,
$\alpha=1$) and that $\Delta L$ and $\Delta U$ can be recomputed in time
linear in their sparsity $O(n\cdot n_{k})$ for _every level_ $k\in[K]$, the
vector $(L^{-1})^{\operatorname{new}}[h,0]^{\top}=(L^{-1}+\Delta
L)[h,0]^{\top}$ can be maintained in $O(n\cdot n_{k})$ time. Furthermore, the
solution

$\displaystyle(U^{-1})^{\operatorname{new}}(L^{-1})^{\operatorname{new}}[h,0]^{\top}=$
$\displaystyle\leavevmode\nobreak\ (U^{-1}+\Delta U)(L^{-1}+\Delta
L)[h,0]^{\top}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\
U^{-1}L^{-1}[h,0]^{\top}+\Delta U\left((L^{-1}+\Delta
L)[h,0]^{\top}\right)+(U^{-1}\Delta L)[h,0]^{\top}$

can be maintained in the same $O(n\cdot n_{k})$ time since the non-zero part
of $\Delta L$ is an $(n_{k}\times n)$ block. The bottom $K$-th level is a
special one, as every level $k$ except the last one involves extra pre-
computation to support the next $(k+1)$-th level. As we show, with some extra
(non-trivial) effort, this fact enables maintaining the $K$-th level in only
$O(n_{K-1}\cdot n_{K})$ time instead of the brute-force $O(n\cdot n_{K})$
time. (Our algorithm implements $K=3$ levels of this framework. By convention,
in later sections we refer to the third level as the “query level”; The first
and second levels are referred as “two-level” updates).

In conclusion, the ideal time per LP-iteration according to our framework is
proportional to

$\displaystyle\sum_{k=1}^{K-1}\underbrace{(n\cdot
n_{k}/\sqrt{n_{k+1}})}_{\text{$K-1$ levels}}+\underbrace{n_{K-1}\cdot
n_{K}}_{K\text{-th\leavevmode\nobreak\ level}},$ (3)

where $n\cdot n_{k}/\sqrt{n_{k+1}}$ is the amortized update cost of level
$k\in[K-1]$, and $n_{K-1}\cdot n_{K}$ is the update cost of the last level
$K$. Carefully balancing the terms by setting geometrically decreasing
thresholds $n_{k}$’s (see Claim H.6), this runtime is optimized to be

$\displaystyle\widetilde{O}(K)\cdot
n^{2+\frac{1}{6\cdot(2^{K-1}-1)}}=\widetilde{O}(K)\cdot n^{2+O(2^{-K})},$ (4)

using the fact that the LP algorithm (see Algorithm 1) has
$\widetilde{O}(\sqrt{n})$ iterations. The formal calculation can be found in
Section H.3. Note that when $K=2$, this running time is precisely
$O^{*}(n^{2+1/6})$, matching [CLS19, LSZ19, Bra20]’s results. For $K=3$
levels, this matches our result $O^{*}(n^{2+1/18})$.

Achieving the runtime predicted by Eq. (4) for $K=\omega(1)$ levels would show
that LPs can be solved in the ideal time $O^{*}(n^{\omega})$. The main
challenge in implementing our technique beyond $K>3$ levels is maintaining the
LU-decomposition of (2) in the desired ($\sim n\cdot n_{k}$) time, and indeed
doing so even for $K=3$ is highly nontrivial, as this paper shows. While
completing this ambitious program is beyond the reach of this paper, we
believe the cascading lazy updates framework may have applications in future
developments in LP and SDP solvers, cutting plane methods, and more generally
for dynamic inverse problems, hence it is worthwhile presenting it at its full
generality.

##  3 Background

This section provides a brief overview of recent developments in LP solvers,
the optimization framework of IPMs and its relation to dynamic inverse
problems.

###  3.1 Recent developments in LP solvers

The recent work of [CLS19] improved the $O(n^{2.5})$ LP algorithm of [Vai89]
to555The running time is actually
$O^{*}(n^{\omega}+n^{2.5-a/2}+n^{1.5+a}+n^{\omega a+0.5}+n^{2a+0.5})$, and the
$n^{\omega a+0.5}$ and $n^{2a+0.5}$ terms also come from query time. But they
are dominated by other terms. In our paper we improved not only the $n^{1+a}$
term, but also these other two terms.

$\displaystyle O^{*}(n^{\omega}+n^{2.5-a/2}+n^{1.5+a}),$ (5)

where $a\leq\alpha$ is a tunable parameter, and $\omega,\alpha$ are the fast
matrix multiplication exponent and its dual, respectively. Note that for the
current values of $\omega\approx 2.38$ and $\alpha\approx 0.31$, this running
time is already $O^{*}(n^{\omega})$. However, under common belief that
$\omega=2$ and $\alpha=1$ [CKSU05, Wil12], the running time is
$O^{*}(n^{2+1/6})$, hence there is still a polynomial gap to the ideal running
time $O^{*}(n^{2})$.

The three main ingredients in [CLS19]’s algorithm are: (i) considering a
_stochastic_ version of the Central Path algorithm (see Algorithm 1 below),
and then leveraging the _robustness_ of this algorithm to design a more
efficient matrix maintenance data structure via _subsampling_ (sparsification
of the “gradient” vector) yielding $o(n^{2})$ query time per iteration. (ii)
_Lazy updates_ : Delaying updates to the projection matrix (associated with
the central path equation) via “soft thresholding” and analyzing their
amortized performance via martingale-based potential analysis. (iii) Using
fast rectangular matrix multiplication to gain extra speedup. Our data
structure will also take advantage of these building blocks.

###  3.2 Optimization: The stochastic central-path algorithm

We use a similar optimization framework as that of [CLS19] (see Section 2.1
for a more detailed explanation and context). Roughly speaking, the Central
Path (CP) algorithm maintains a primal-dual pair of vectors, $x^{(i)}$ and
$s^{(i)}$, and iteratively shrinks the _duality gap_
$\mu^{(i)}:=\sum_{j=1}^{n}x_{j}^{(i)}s_{j}^{(i)}$ by $\sim(1-1/\sqrt{n})$ in
each iteration, until converging to a feasible point ($\mu^{(i)}\approx 0$).
Hence, the Central Path algorithm has a total of $O(\sqrt{n})$ iterations. In
matrix notation, this algorithm essentially boils down to implementing the
following iterative algorithm [CLS19]:

Algorithm 1 Stochastic Central Path

1:$i\leftarrow 1$, initialize $x,s\in\mathbb{R}^{n}$

2:while $i<\sqrt{n}$ do $\triangleright$ In each iteration, we hope
$\mu\approx t$

3: $t\leftarrow t\cdot(1-1/\sqrt{n})$ $\triangleright$ target decrease of
duality gap

4: $\mu\leftarrow x\cdot s$ $\triangleright$ actual decrease in duality gap

5: Compute $\delta_{\mu}$ based on $-\frac{\mu}{\sqrt{n}}$ and the gradient
$-\nabla\Psi(\mu/t-1)$.

6:
$P\leftarrow\sqrt{\frac{X}{S}}A^{\top}(A\frac{X}{S}A^{\top})^{-1}A\sqrt{\frac{X}{S}}$
$\triangleright$ _matrix inverse_ , matrix-matrix mult.

7:
$\delta_{x}\leftarrow\frac{X}{\sqrt{XS}}(I-P)\frac{1}{\sqrt{XS}}\delta_{\mu}$,
$\delta_{s}\leftarrow\frac{S}{\sqrt{XS}}P\frac{1}{\sqrt{XS}}\delta_{\mu}$
$\triangleright$ matrix-vector mult.

8: $x\leftarrow x+\delta_{x}$, $s\leftarrow s+\delta_{s}$, $i\leftarrow i+1$

9:end while

Here, $X=\operatorname{diag}(x)$, $S=\operatorname{diag}(s)$ are the primal
and dual vectors, and $A\in\mathbb{R}^{n\times n}$ is the (fixed) LP
constraint matrix. $\Psi$ is a potential function measuring how close $\mu$ is
from $t$ (the “target” duality gap), and the vector $\delta_{\mu}$ has two
purposes: decreasing $\mu$ by a $(1-1/\sqrt{n})$ factor while keeping the
potential function bounded. The vectors $\delta_{x},\delta_{s}$ compute the
disposition of the primal and dual vectors in each iteration. $P$ is an
orthogonal _projection matrix_ ($P^{2}=P$ and $P=P^{\top}$), and the formulas
$\frac{X}{\sqrt{XS}}$, $\frac{S}{\sqrt{XS}}$, $\frac{1}{\sqrt{XS}}$, and
$\frac{X}{S}\in\mathbb{R}^{n\times n}$ are the diagonal matrices of the
corresponding vectors.

A key observation in [CLS19] is that this algorithm is _robust to small
perturbations_ along the central path: Denoting by $w$ the vector $x/s$, and
by $h$ the vector $\frac{\delta_{\mu}}{\sqrt{XS}}$, [CLS19] shows that in the
above algorithm, is enough to approximately maintain
$w^{\operatorname{appr}}\approx_{\epsilon_{\mathrm{mp}}}w,h^{\operatorname{appr}}\approx_{\epsilon_{\mathrm{mp}}}h,$
where $\epsilon_{\mathrm{mp}}<1/4$ and $\approx_{\epsilon_{\mathrm{mp}}}$
denotes coordinate-wise approximation.

###  3.3 Data structures: Projection maintenance

The main bottleneck of Algorithm 1 is to efficiently maintain the approximate
projection matrix

$\displaystyle
P(w^{\operatorname{appr}})=\sqrt{W^{\operatorname{appr}}}A^{\top}(AW^{\operatorname{appr}}A^{\top})^{-1}A\sqrt{W^{\operatorname{appr}}},$
(6)

recalculating the queries $r:=P(w^{\operatorname{appr}})h$ on line 7, where
$h:=\delta_{\mu}/\sqrt{XS}$. There are $O(\sqrt{n})$ iterations.

##### Lazy updates.

It was already observed in [Vai89] that since each iteration only changes the
$\ell_{2}$ mass of $w$ by a small amount (which can be turned into an
$\ell_{0}$ sparsity guarantee by “rounding” and absorbing a small error), most
of the time the queries can be answered efficiently by computing the _low-
rank_ incremental change in $P$, amortizing away the rare cases where too many
coordinates of $w$ have changed, which are handled using brute force fast
matrix multiplication. As noted above, [CLS19] further used the power of fast
rectangular matrix multiplication: By definition of $\alpha$, for any
threshold parameter $a\leq\alpha$, the complexity of multiplying an $n\times
n^{a}$ rectangular matrix by an $n^{a}\times n$ rectangular matrix is the same
as multiplying an $n\times 1$ vector with a $1\times n$ vector, so they only
update $P$ when _at least $n^{a}$ coordinates_ of $w$ have changed. [CLS19]
further make a clever use of _soft thresholding_ on $n^{a}$, which combined
with a potential function analysis, yields an amortized
$O(n^{\omega-1/2}+n^{2-a/2})$ update time per iteration. Note that the
$n^{2-a/2}$ term needs to be balanced with query time. [LSZ19] and [Bra20]
both follow the same update scheme.

##### Fast queries.

Computing the queries $r=P(w^{\operatorname{appr}})h$ in each iteration from
scratch takes $n^{2}$ time since the projection matrix may be very dense. The
three papers [CLS19, LSZ19, Bra20] proposed different techniques for speeding
up this matrix-vector multiplication to $o(n^{2})$. In [CLS19, LSZ19], the
authors use the idea of “iterating and sketching” [Son19], an adaptive version
of the classic “sketch and solve” paradigm [CW13]. In [Bra20], the author
maintains the query answer $r$ directly, exploiting the observation that the
vector $h$ is also slowly changing (and not just updates $w$). Both techniques
([CLS19, Bra20]) are essentially using _sparsification_ of the vector $h$,
hence involve a “right hand side” linear operation. In contrast, [LSZ19] uses
a “left-hand side” operation by _sketching the projection matrix_ itself,
effectively making it smaller.

##### Sampling on the right [CLS19].

Here the idea is to apply a $O(\sqrt{n})$-sparse _diagonal sampling matrix_
$D\in\mathbb{R}^{n\times n}$ on the right hand side of the maintained matrix:

$\displaystyle\sqrt{W^{\operatorname{appr}}}A^{\top}(AW^{\operatorname{appr}}A^{\top})^{-1}A\sqrt{W^{\operatorname{appr}}}\underbrace{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}D}}_{\text{sample\leavevmode\nobreak\
right}}h.$

After this sampling, the vector $Dh$ becomes $O(\sqrt{n})$-sparse, so
computing the multiplication of a matrix with $Dh$ takes $O(n^{1.5})$ time.
Also, since the rank of the change in $W$ is guaranteed to be at most $n^{a}$,
there is also a $O(n^{1+a})$ term in the query time.

##### Sketching on the left [LSZ19].

The idea here is to apply a (subsampled) Hadamard/Fourier transform matrix
[LDFU13, PSW17] $R\in\mathbb{R}^{\sqrt{n}\times n}$, by sketching on the left
hand side:

$\displaystyle\underbrace{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}R^{\top}R}}_{\text{sketch\leavevmode\nobreak\
left}}\sqrt{W^{\operatorname{appr}}}A^{\top}(AW^{\operatorname{appr}}A^{\top})^{-1}A\sqrt{W^{\operatorname{appr}}}h.$

In this way the matrix $RM$ has size $\sqrt{n}\times n$,666Intuitively,
$O(\sqrt{n})$ rows is the minimal sketch size one can hope for, as this
ensures that with $O(\sqrt{n})$ CP iterations, every coordinate $i\in[n]$ has
a constant chance to be sampled so multiplying this matrix with a vector takes
$O(n^{1.5})$ time. So the final query time is also $O(n^{1.5}+n^{1+a})$.
[LSZ19] algorithm is also maintaining some extra vectors, e.g., the
explicit/implicit version of $x,s$.

##### Maintaining query answers [Bra20].

Brand observed that is possible to maintain not only the inverse matrix
$P(w^{\operatorname{appr}})$ via lazy updates, but also the previously
computed query answers $r$, by observing that the vector $h$ it also changes
slowly. In each iteration, the new $r$ is computed as:

$\displaystyle
r+\sqrt{W^{\operatorname{appr}}}A^{\top}(AW^{\operatorname{appr}}A^{\top})^{-1}A\sqrt{W^{\operatorname{appr}}}\underbrace{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Delta
h}}_{\text{change in }h}.$

[Bra20] chooses a similar sparsity threshold $n^{a}$ for the vector and
updates the maintained $r$ when $\Delta h$ exceeds $n^{a}$, so $\Delta h$ is
guaranteed to be $n^{a}$-sparse at query time. As such, the query time is
$O(n^{1+a})$. It is noteworthy that this algorithm is deterministic, as it
avoids sketching/sampling altogether. Indeed, the main motivation of [Bra20]
was derandomizing [CLS19].

##### Our approach.

We show how to break the $O(n^{1+a})$ barrier for query time, by _combining_
both left-hand and right-hand linear transformations on $P$, together with the
_cascading lazy updates_ technique from Section 2. Such combination is needed
to further “compress” the matrix-vector multiplication when (re)calculating
$r$. It turns out that using two sources of randomness—sampling on the right
[CLS19] \+ sketching on the left [LSZ19]—obliterates the error analysis (which
needs to be controlled to ensure convergence), and this is intuitively where
[Bra20]’s deterministic alternative is useful for us as a substitute to right-
hand-sampling.

##  4 Detailed Technical Overview

This section provides a detailed, streamlined technical overview of the proof
of Theorem 1.1, which we restate formally below. This section should be
understood as a self-contained extended abstract of our entire algorithm.
Formal proofs of all technical claims can be found in the Appendix.

######  Theorem 4.1 (Main result).

Given a linear program $\min_{Ax=b,x\geq 0}c^{\top}x$ with no redundant
constraints. Assume that the polytope has diameter $R$ in $\ell_{1}$ norm,
namely, for any $x\geq 0$ with $Ax=b$, we have $\|x\|_{1}\leq R$.

Then, for any $\delta\in(0,1]$, Main$(A,b,c,\delta,a,\widetilde{a})$
(Algorithm 17) outputs $x\geq 0$ such that

$\displaystyle c^{\top}x\leq\min_{Ax=b,x\geq
0}c^{\top}x+\delta\|c\|_{\infty}R,\text{\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ and\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
}\|Ax-b\|_{1}\leq\delta\cdot(R\|A\|_{1}+\|b\|_{1})$

in expected time

$\displaystyle\leavevmode\nobreak\
O(n^{\omega+o(1)}+n^{2.5-a/2+o(1)}+n^{1.5+a-\widetilde{a}/2+o(1)}+n^{0.5+a+(\omega-1)\widetilde{a}})\cdot\log(n/\delta)$

for any $0<a\leq\alpha$ and $0<\widetilde{a}\leq\alpha a$. In particular, so
long as the constants of fast matrix multiplication satisfy $\omega>2.055$ and
$\alpha>5-2\omega$, general LPs can be solved in $O(n^{\omega+o(1)})$ time. In
the ideal case that $\omega=2$ and $\alpha=1$, the running time is
$n^{2+1/18}=n^{2.055}$.

The first two terms $n^{\omega}$ and $n^{2.5-a/2}$ of our running time are the
same as [CLS19], stemming from the amortized cost of lazy updates (for $K=1$
levels). The $n^{1.5+a-\widetilde{a}/2}$ term comes from the amortized cost of
our cascading lazy updates algorithm for the $K=2$nd level update. Finally,
the $n^{0.5+a+(\omega-1)\widetilde{a}}$ term is the worst-case cost of the
query algorithm. We note that a prerequisite for achieving the runtime of
Theorem 1.1 is removing the explicit $n^{1+a}$ term as well as the two
(implicit) $n^{a\omega}$, $n^{2a}$ terms in the query time of all previous
works. Below we elaborate on how this is achieved step by step, as shown in
Table 1.

Statement | Technique | Time | Ideal | Choice
---|---|---|---|---
[CLS19] |  Sec. 3 | $n^{2.5-a/2}+n^{0.5+2a}+n^{0.5+\omega a}+n^{1.5+a}$ | $n^{2+1/6}$ | $a=2/3$
Thm. 4.2 |  Sec. 4.1.1 | $n^{2.5-a/2}+n^{0.5+2a}+n^{0.5+\omega a}$ | $n^{2+1/10}$ | $a=4/5$
Thm. 4.3 |  Sec. 4.1.2 | $n^{2.5-a/2}+n^{0.5+2a}$ | $n^{2+1/10}$ | $a=4/5$
Thm. 4.1 |  Sec. 4.1.3 | $n^{2.5-a/2}+n^{1.5+a-\widetilde{a}/2}+n^{0.5+a+(\omega-1)\widetilde{a}}$ | $n^{2+1/18}$ | $a=8/9,\widetilde{a}=2/3$
Table 1: We ignore the $n^{\omega}$ term, and also ignore all the $n^{o(1)}$
terms. Ideal denotes the resulting running time when $\omega=2$ and
$\alpha=1$. (Note that the current values are $\omega\sim 2.38$ and
$\alpha\sim 0.31$). Figure 2: $(a)$: In each iteration, we are given
$w^{\operatorname{new}}$ which is changing slowly. The algorithm will find
$w^{\operatorname{appr}}$ such that $w^{\operatorname{appr}}$ is coordinate-
wise $\epsilon_{\mathrm{mp}}$ close to $w^{\operatorname{new}}$(pink wave
line), and it is also close to $v$ in $\ell_{0}$ norm(blue solid line). $(b)$:
Based on $(a)$, we add an intermediate level $\widetilde{v}$, such that in the
query, $\|v-\widetilde{v}\|_{0}\leq n^{a}$(blue solid line) and
$\|\widetilde{v}-w^{\operatorname{appr}}\|_{0}\leq n^{\widetilde{a}}$(green
dashed line).

We note that, although our better running time is achieved by breaking the
three bottlenecks of the _query_ time of previous works, it actually requires
non-trivial design of both the _query_ part and the _update_ part. Indeed, our
data structure involves an entirely new update subroutine for second-level of
updates.

In the next Section 4.1, we walk through the techniques for removing the three
Query bottlenecks one by one, while introducing the data structure members
that must be maintained in order to make queries faster. At a high level,
these members correspond to maintaining the components appearing in the LU-
decomposition part of our cascading framework for $K=2$ levels (described in
Section 2, see Figure 2 for illustration). In Section 4.2 we describe and
analyze the Update algorithm of our data strucutre. This part uses a
combination of “soft thresholding” for two levels of updates (Figure 6),
capturing the cascading lazy updates process. We remark that this type of
analysis may be of independent interest to “bootstrapping” techniques in
dynamic matrix problems.

###  4.1 Our Query Algorithm

Our dynamic data structure relies on the following three main ingredients, for
removing each of the three bottlenecks shown in Table 1 respectively:

  1. 1.

“Compressing” the projection matrix using linear operations on both sides
(Section 4.1.1): We show that by combining a sketching on the left to reduce
the matrix size from $n\times n$ to $\sqrt{n}\times n$, with the direct
maintenance of the previous query answer that makes the new query vector
sparse, already breaks the $n^{1+a}$ query time barrier of previous works.

  2. 2.

Faster queries through the “cascading lazy updates” framework (Section 4.1.2):
The core technique of previous works is exploiting the fact that the inverse
matrix and the query vector are both changing slowly over iterations, hence
the data structure can maintain all intermediate values from previous
iterations, and only re-compute the differences. As explained in the last part
of Section 2, we show that even the _“derivative”_ of updates and queries is
slowly changing (these are the “sparsity thresholds” we allude to in Section
2). It is therefore natural to maintain a second-level of values (“changes” of
the inverse matrix and the query vector), where each new iteration only
computes the changes in these second-level values. We update both levels
according to the _cascading lazy updates_ framework (recall Figure 1). We show
this technique removes the $n^{a\omega}$ term in the query time of previous
data structures (Table 1).

  3. 3.

Maintaining the components of the second-level structure efficiently (Section
4.1.3): This is done by essentially recursing the approach of maintaining the
first-level to the second-level, by carefully designing the maintained objects
(matrix-products, vectors, sets). This removes the $n^{2a}$ term in the query
time, which in fact appeared in _multiple_ places in previous algorithms.

We now turn to a detailed description of the query algorithm, using [CLS19] as
a baseline. Recall that each iteration of the CP algorithm generates two
vectors $w$ and $h$ from the previous outputs of the data structure. The data
structure needs to re-calculate

$\displaystyle\sqrt{W}A^{\top}(AWA^{\top})^{-1}A\sqrt{W}h.$

By robustness of the stochastic CP algorithm, it suffices to output an
approximated value

$\displaystyle
r=\underbrace{({\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}R})^{\top}}_{\rm
de-
sketch}\cdot\underbrace{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}R}\sqrt{W^{\operatorname{appr}}}A^{\top}(AW^{\operatorname{appr}}A^{\top})^{-1}A\sqrt{W^{\operatorname{appr}}}h^{\operatorname{appr}}}_{\rm
sketch},$

where for some fixed parameter $\epsilon_{\mathrm{mp}}<1/4$ the data structure
guarantees that $w^{\operatorname{appr}}\approx_{\epsilon_{\mathrm{mp}}}w$ and
$h^{\operatorname{appr}}\approx_{\epsilon_{\mathrm{mp}}}h$. We use the same
sketching matrix $R$ as that of [LSZ19] which satisfies
$\operatorname*{{\mathbb{E}}}[R^{\top}R]=I$ and a guarantee that the variance
and $\ell_{\infty}$ error of $(R^{\top}RPh-Ph)$ are both small. $R$ has size
$\sqrt{n}\times n$ – this value is natural, since we have $\sim\sqrt{n}$ CP
iterations, hence using $o(\sqrt{n})$ linear measurements would fail to even
detect changes in all $n$ coordinates. Also, the size of $R$ allows us to pre-
batch $O(\sqrt{n})$ copies of sketching matrices in the beginning, which only
takes $O(n^{\omega})$ time.

As in [CLS19, Bra20], our data structure maintains a member $v$ that serves as
a proxy for $w$, and a member $g$ that serves as a proxy for $h$. If the new
value $w^{\operatorname{appr}}$ is too far from $v$, then in addition to
computing the new displacement $r$, the data structure also updates $v$ along
with all of its members that depend on $v$. An analogous scheme is used for
$g$ w.r.t $h^{\operatorname{appr}}$. Since the new values
$w^{\operatorname{appr}}$ and $h^{\mathrm{sample}}$ show up in three places in
the output $r$, from now on we will refer to these three places as the left
part, the middle part, and the right part, and label them as
$a^{\operatorname{new}}\in\mathbb{R}^{\sqrt{n}\times n}$,
$b^{\operatorname{new}}\in\mathbb{R}^{n\times n}$,
$c^{\operatorname{new}}\in\mathbb{R}^{n}$ for ease of presentation.
Accordingly, we use $a\in\mathbb{R}^{\sqrt{n}\times n}$,
$b\in\mathbb{R}^{n\times n}$, $c\in\mathbb{R}^{n}$ to denote the values that
depend on $v$ and $g$:

$\displaystyle\underbrace{R\sqrt{W^{\operatorname{appr}}}}_{a^{\operatorname{new}}}\underbrace{A^{\top}(AW^{\operatorname{appr}}A^{\top})^{-1}A}_{b^{\operatorname{new}}}\underbrace{\sqrt{W^{\operatorname{appr}}}h^{\operatorname{appr}}}_{c^{\operatorname{new}}}\;\;\;,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\underbrace{R\sqrt{V}}_{a}\underbrace{A^{\top}(AVA^{\top})^{-1}A}_{b}\underbrace{\sqrt{V}g}_{c}.$

We also denote $\partial a:=a^{\operatorname{new}}-a$, $\partial
b:=b^{\operatorname{new}}-b$, $\partial c:=c^{\operatorname{new}}-c$.

For a tunable parameter $a\in(0,\alpha]$, the worst case query time per
iteration of [CLS19]’s data structure for implementing this process is
$t_{q}=n^{2a}+n^{\omega a}+n^{1+a}$, the three terms come from the cost of
recomputing the query $r$. In the remainder of this subsection, we describe
how this query time can be improved to

$\displaystyle t_{q}=n^{a+(\omega-1)\widetilde{a}}$ (7)

for any $a\in(0,\alpha]$ and $\widetilde{a}\in(0,\alpha a]$. We show this in
three steps, removing the terms $n^{1+a}$, $n^{\omega a}$ , and $n^{2a}$ one
by one.

####  4.1.1 Technique for removing $n^{1+a}$

We combine the “sketching on the left” technique of [LSZ19] and the “query
maintenance” technique of [Bra20] to remove this term. In [Bra20] the query is
computed as

$\displaystyle
r=\underbrace{\sqrt{W^{\operatorname{appr}}}}_{a^{\operatorname{new}}}\bigg{(}\underbrace{\beta_{2}}_{bc}+\underbrace{M}_{b}\underbrace{\big{(}\sqrt{W^{\operatorname{appr}}}h^{\operatorname{appr}}-\sqrt{V}g\big{)}}_{\partial
c}+\underbrace{\big{(}-M_{S}(\Delta_{S,S}^{-1}+M_{S,S})^{-1}(M_{S})^{\top}\big{)}}_{\partial
b}\underbrace{\sqrt{W^{\operatorname{appr}}}h^{\operatorname{appr}}}_{c^{\operatorname{new}}}\bigg{)},$

where $M:=A^{\top}(AVA^{\top})^{-1}A$,
$\beta_{2}:=M\sqrt{V}g\in\mathbb{R}^{n}$ are members that the data structure
maintains, $\Delta:=W^{\operatorname{appr}}-V$, and
$S:=\operatorname{supp}(w^{\operatorname{appr}}-v)$. The subscript $M_{S}$
means taking the sub matrix of columns of $M$ in the set $S$. Note that this
output is $a^{\operatorname{new}}(bc+b\cdot\partial c+\partial b\cdot
c^{\operatorname{new}})=a^{\operatorname{new}}b^{\operatorname{new}}c^{\operatorname{new}}$.

The $n^{1+a}$ term shows up in three places when computing different terms:

  * •

In $a^{\operatorname{new}}\cdot b\cdot\partial c$, multiplying a $n\times n$
matrix $M$ with a $n^{a}$-sparse vector
$\big{(}\sqrt{W^{\operatorname{appr}}}h^{\operatorname{appr}}-\sqrt{V}g\big{)}$.

  * •

In $a^{\operatorname{new}}\cdot\partial b\cdot c^{\operatorname{new}}$,
multiplying a $n\times n^{a}$ matrix $M_{S}$ with a $n^{a}\times 1$ vector
$(\Delta^{-1}_{S,S}+M_{S,S})^{-1}(M_{S})^{\top}\sqrt{W^{\operatorname{appr}}}h^{\operatorname{appr}}$.

  * •

In $a^{\operatorname{new}}\cdot\partial b\cdot c^{\operatorname{new}}$,
multiplying a $n^{a}\times n$ matrix $(M_{S})^{\top}$ with a $n\times 1$
vector $\sqrt{W^{\operatorname{appr}}}h^{\operatorname{appr}}$.

The last one is easy, and we deal it by splitting
$c^{\operatorname{new}}\in\mathbb{R}^{n}$ into $c+\partial c$ again:

$\displaystyle(M_{S})^{\top}\underbrace{\sqrt{W^{\operatorname{appr}}}h^{\operatorname{appr}}}_{c^{\operatorname{new}}}=(\beta_{2})_{S}+(M_{S})^{\top}\underbrace{(\sqrt{W^{\operatorname{appr}}}h^{\operatorname{appr}}-\sqrt{V}g)}_{\partial
c}.$

Since the $\partial c$ term is $n^{a}$-sparse, this computation only takes
$n^{2a}$ time now.

We deal with the first two $n^{1+a}$ terms by adding the sketching matrix $R$
on the left, and splitting $a^{\operatorname{new}}$ into $a+\partial a$. Aside
from maintaining $M=A^{\top}(AVA^{\top})^{-1}A$ and $\beta_{2}=M\sqrt{V}g$, we
further maintain their sketched versions:

$Q=R\sqrt{V}M\in\mathbb{R}^{\sqrt{n}\times n},\leavevmode\nobreak\
\leavevmode\nobreak\ \beta_{1}=R\sqrt{V}\beta_{2}\in\mathbb{R}^{\sqrt{n}}.$

In addition to the temporary variables
$S=\operatorname{supp}(w^{\operatorname{appr}}-v)$ and
$\Delta=W^{\operatorname{appr}}-V$, we also define
$\Gamma:=\sqrt{W^{\operatorname{appr}}}-\sqrt{V}$. We have the guarantee that
$\Delta$, $\Gamma$ and
$(\sqrt{W^{\operatorname{appr}}}h^{\operatorname{appr}}-\sqrt{V}g)$ are all
$n^{a}$-sparse and $|S|\leq n^{a}$, because otherwise the algorithm would
first update $V$ and $g$. The query part of our new data structure becomes

$\displaystyle r_{1}:=$ $\displaystyle\leavevmode\nobreak\
\underbrace{\beta_{1}}_{abc},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ r_{2}:=\leavevmode\nobreak\
(\underbrace{Q_{S}}_{ab}+\underbrace{R\Gamma M}_{\partial a\cdot
b})\underbrace{(\sqrt{W^{\operatorname{appr}}}h^{\operatorname{appr}}-\sqrt{V}g)}_{\partial
c},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
r_{3}:=\leavevmode\nobreak\ \underbrace{R\Gamma}_{\partial
a}\cdot\underbrace{\beta_{2}}_{bc}$ $\displaystyle r_{4}:=$
$\displaystyle\leavevmode\nobreak\ \hbox
to0.0pt{$\displaystyle\underbrace{\phantom{(Q_{S}+R\Gamma
M_{S})}}_{a^{\operatorname{new}}\cdot
M_{S}}$\hss}(Q_{S}+R\Gamma\cdot\overbrace{M_{S})(-(\Delta_{S,S}^{-1}+M_{S,S})^{-1})\hbox
to0.0pt{$\displaystyle\underbrace{\phantom{\Big{(}(M_{S})^{\top}\cdot\big{(}\sqrt{W^{\operatorname{appr}}}h^{\operatorname{appr}}-\sqrt{V}g\big{)}+\beta_{2,S}\Big{)}}}_{(M_{S})^{\top}\cdot
c^{\operatorname{new}}}$\hss}\Big{(}(M_{S})^{\top}}^{\partial
b}\cdot\big{(}\sqrt{W^{\operatorname{appr}}}h^{\operatorname{appr}}-\sqrt{V}g\big{)}+\beta_{2,S}\Big{)}$

Note that this output is

$\displaystyle r:=\leavevmode\nobreak\
\underbrace{abc}_{r_{1}}+\underbrace{(a+\partial a)b\cdot\partial
c}_{r_{2}}+\underbrace{\partial a\cdot
bc}_{r_{3}}+\underbrace{a^{\operatorname{new}}\cdot\partial b\cdot
c^{\operatorname{new}}}_{r_{4}}=\leavevmode\nobreak\
a^{\operatorname{new}}b^{\operatorname{new}}c^{\operatorname{new}}.$

Recall that the $n^{1+a}$ term stemmed from multiplying an $n\times n$ matrix
by an $n^{a}$-sparse vector when computing $a^{\operatorname{new}}\cdot
b\cdot\partial c$. We split this term as $a^{\operatorname{new}}\cdot
b\cdot\partial c:=a\cdot b\cdot\partial c+\partial a\cdot b\cdot\partial c$
(see the formula for $r_{2}$). The data structure will now maintain
$ab:=Q_{S}$, whose size is only $\sqrt{n}\times n^{a}$, hence computing
$ab\cdot\partial c$ only takes $n^{1/2+a}<n^{1.5}$ time now. Note that when
computing $\partial a\cdot b\cdot\partial c$, the $n\times n$ matrix $M$ is
“sandwiched” by a $n^{a}$-sparse diagonal matrix $\Gamma$ on the left and a
$n^{a}$-sparse vector
$(\sqrt{W^{\operatorname{appr}}}h^{\operatorname{appr}}-\sqrt{V}g)$ on the
right, thus computing the product $\Gamma
M(\sqrt{W^{\operatorname{appr}}}h^{\operatorname{appr}}-\sqrt{V}g)$ takes
$n^{2a}$ time.

The last $n^{1+a}$ bottleneck from [CLS19] is removed in a similar way,
yielding the following intermediate result:

######  Theorem 4.2 (Informal, first improvement).

For any $a\leq\alpha$, there is a randomized algorithm for solving general LPs
in expected time $O^{*}(n^{\omega}+n^{2.5-a/2}+n^{0.5+a\omega})$.

Note that for $\omega=2$ and $\alpha=1$, this approach already yields an
improved $n^{2+1/10}$ LP algorithm.

####  4.1.2 Technique for removing $n^{\omega a}$

The $n^{\omega a}$ term in previous data structures came from computing the
inverse of an $n^{a}\times n^{a}$ matrix $(\Delta_{S,S}^{-1}+M_{S,S})$, and
this inverse is still present in the intermediate algorithm described in
Section 4.1.1 (see the formula for $r_{4}$). This is where the _cascading lazy
updates_ technique comes useful – we shall remove the $n^{\omega a}$ by
showing how to implement it for $K=2$ “levels”. We now provide the details of
this data structure.

Once again, the main observation is that since we’ve already computed
$(\Delta_{S,S}^{-1}+M_{S,S})^{-1}$ in previous iterations, we do not need to
re-compute it from scratch. Instead, we only need to compute the difference
between the new inverse matrix and the old one. More concretely, we maintain a
second-level data structure member $\widetilde{v}$. $\widetilde{v}$ keeps a
_closer_ distance with $w^{\operatorname{appr}}$ than $v$, and is therefore
updated more frequently. We update $v$ whenever
$\|\widetilde{v}-v\|_{0}>n^{a}$ (for some $a\leq\alpha$) and update
$\widetilde{v}$ whenever
$\|w^{\operatorname{appr}}-\widetilde{v}\|>n^{\widetilde{a}}$ (for some
$\widetilde{a}\leq\alpha a$). By abuse of notation, we define

$\displaystyle S:=\operatorname{supp}(\widetilde{v}-v),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
S^{\operatorname{new}}:=\operatorname{supp}(w^{\operatorname{appr}}-v),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \partial
S:=\operatorname{supp}(w^{\operatorname{appr}}-\widetilde{v}).$ (8)

In this overview we can think of $S^{\operatorname{new}}=S\cup\partial S$.
Later we also handle the possibly non-empty set $S^{\prime}=(S\cup\partial
S)\backslash S^{\operatorname{new}}$, but the key ideas are the same. The
updates guarantee that in the query we always have that
$|S^{\operatorname{new}}|\leq n^{a}$ and $|\partial S|\leq n^{\widetilde{a}}$.
Our data structure maintains a second-level member

$\displaystyle B:=(\Delta_{S,S}^{-1}+M_{S,S})^{-1}\in\mathbb{R}^{n^{a}\times
n^{a}}.$

Observe that the new matrix
$((\Delta^{\operatorname{new}})_{S^{\operatorname{new}},S^{\operatorname{new}}}^{-1}+M_{S^{\operatorname{new}},S^{\operatorname{new}}})$
only differs from $B^{-1}=(\Delta_{S,S}^{-1}+M_{S,S})$ on entries in
$S^{\operatorname{new}}\times\partial S$ and $\partial S\times
S^{\operatorname{new}}$. (See the left part of Figure 3.)

Figure 3: An illustration of the construction of $U^{\prime},C,U^{\top}$. $N$
is defined as
$((\Delta^{\operatorname{new}})_{S^{\operatorname{new}},S^{\operatorname{new}}}^{-1}+M_{S^{\operatorname{new}},S^{\operatorname{new}}})-(\Delta_{S,S}^{-1}+M_{S,S})$.
$I$ denotes the identity matrix.

So we have the following decomposition that can be computed in
$O(n^{\widetilde{a}+a})$ time:

$\displaystyle
U^{\prime}CU^{\top}=((\Delta^{\operatorname{new}})_{S^{\operatorname{new}},S^{\operatorname{new}}}^{-1}+M_{S^{\operatorname{new}},S^{\operatorname{new}}})-(\Delta_{S,S}^{-1}+M_{S,S}),$

where $U^{\prime},U\in\mathbb{R}^{n^{a}\times n^{\widetilde{a}}}$, and
$C\in\mathbb{R}^{n^{\widetilde{a}}\times n^{\widetilde{a}}}$. In fact
$U^{\prime}$ and $U$ are both constructed by taking a submatrix from $M$ and
concatenate it with two identity matrices (see Figure 3), this will be useful
in the next section where we remove the $n^{2a}$ term. Now we can use Woodbury
identity to compute

$\displaystyle((\Delta^{\operatorname{new}})^{-1}_{S^{\operatorname{new}},S^{\operatorname{new}}}+M_{S^{\operatorname{new}},S^{\operatorname{new}}})^{-1}=\leavevmode\nobreak\
(B^{-1}+U^{\prime}CU^{\top})^{-1}=\leavevmode\nobreak\
B-BU^{\prime}(C^{-1}+U^{\top}BU^{\prime})^{-1}U^{\top}B.$

The most expensive part of this formula is to compute the multiplication
$U^{\top}B$ and $BU^{\prime}$, and the inverse
$(C^{-1}+U^{\top}BU^{\prime})^{-1}$. Since $\widetilde{a}\leq\alpha a$, using
fast rectangular matrix multiplication, multiplying a $n^{\widetilde{a}}\times
n^{a}$ matrix $U^{\top}$ with a $n^{a}\times n^{a}$ matrix $B$ takes
$O(n^{2a})$ time. Computing $BU^{\prime}$ takes the same time. Computing the
inverse of a $n^{\widetilde{a}}\times n^{\widetilde{a}}$ matrix
$(C^{-1}+U^{\top}BU^{\prime})$ takes $n^{\widetilde{a}\omega}={\cal
T}_{\mathrm{mat}}(n^{\widetilde{a}},n^{\widetilde{a}},n^{\widetilde{a}})\leq{\cal
T}_{\mathrm{mat}}(n^{a},n^{a},n^{\widetilde{a}})=n^{2a}$. All other parts of
the query algorithm remain the same as in Section 4.1.1. Hence, so far, the
running time is upper bounded by $O(n^{2a})$:

######  Theorem 4.3 (informal, second improvement).

For any $a\leq\alpha$, there is a randomized algorithm for solving general LPs
in expected time $O^{*}(n^{\omega}+n^{2.5-a/2}+n^{0.5+2a})$.

We remark that the second-level members $B$ and the local variables $U$, $C$
and $U^{\prime}$ that the data structure maintains, correspond to the
variables in the cascading lazy update framework of Section 2: The matrix $B$
here is exactly the same inverse $B$ as defined in Section 2 for $K=2$ (see
Eq.(H.1)). In the same vein, the block $C^{-1}$ here corresponds to the term
$-C_{2}^{-1}-V_{2}^{\top}M^{-1}U_{2}$, $U^{\top}$ corresponds to
$V_{2}^{\top}M^{-1}U_{1}$, and $U^{\prime}$ corresponds to
$V_{1}^{\top}M^{-1}U_{2}$ of Eq.(H.1). That said, since Section 2 deals with a
simplified version of our actual inverse problem, we will need to maintain
several other ad-hoc members to achieve the claimed running time. We turn to
describe those next.

####  4.1.3 Technique for removing $n^{2a}$

Now we show how to remove the $n^{2a}$ term from the current algorithm as
described in Section 4.1.1, with the inverse matrix of $r_{4}$ computed as
described in Section 4.1.2). The $n^{2a}$ term appears in the computation of
both $r_{2}$ and $r_{4}$. Since removing the $n^{2a}$ terms in $r_{4}$ is more
difficult, we mainly focus on the $r_{4}$ term in this section. In analogy to
the second-level member $\widetilde{v}$, we also maintain $\widetilde{g}$, and
update it whenever
$\|h^{\operatorname{appr}}-\widetilde{g}\|_{0}>n^{\widetilde{a}}$. Similar to
the definition of $S^{\operatorname{new}},S,\partial S$ (Eq. 8), we shall use
the following three variants of notations:

$\displaystyle\Delta^{\operatorname{new}}=$ $\displaystyle\leavevmode\nobreak\
W^{\operatorname{appr}}-V,$ $\displaystyle\leavevmode\nobreak\ \Delta=$
$\displaystyle\leavevmode\nobreak\ \widetilde{V}-V,$
$\displaystyle\leavevmode\nobreak\ \partial\Delta=$
$\displaystyle\leavevmode\nobreak\ W^{\operatorname{appr}}-\widetilde{V},$
$\displaystyle\Gamma^{\operatorname{new}}=$ $\displaystyle\leavevmode\nobreak\
\sqrt{W^{\operatorname{appr}}}-\sqrt{V},$ $\displaystyle\leavevmode\nobreak\
\Gamma=$ $\displaystyle\leavevmode\nobreak\ \sqrt{\widetilde{V}}-\sqrt{V},$
$\displaystyle\leavevmode\nobreak\ \partial\Gamma=$
$\displaystyle\leavevmode\nobreak\
\sqrt{W^{\operatorname{appr}}}-\sqrt{\widetilde{V}},$
$\displaystyle\xi^{\operatorname{new}}=$ $\displaystyle\leavevmode\nobreak\
\sqrt{W^{\operatorname{appr}}}h^{\operatorname{appr}}-\sqrt{V}g,$
$\displaystyle\leavevmode\nobreak\ \xi=$ $\displaystyle\leavevmode\nobreak\
\sqrt{\widetilde{V}}\widetilde{g}-\sqrt{V}g,$
$\displaystyle\leavevmode\nobreak\ \partial\xi=$
$\displaystyle\leavevmode\nobreak\
\sqrt{W^{\operatorname{appr}}}h^{\operatorname{appr}}-\sqrt{\widetilde{V}}\widetilde{g}.$

Figure 4: The three variants of notations. The left vector
$w^{\operatorname{appr}}-v$ corresponds to the notations
$X^{\operatorname{new}}$, and
$|\operatorname{supp}(w^{\operatorname{appr}}-v)|\leq n^{a}$. The middle
vector $\widetilde{v}-v$ corresponds to the notations $X$, and
$|\operatorname{supp}(\widetilde{v}-v)|\leq n^{a}$. The right vector
$w^{\operatorname{appr}}-\widetilde{v}$ corresponds to the notations $\partial
X$, and $|\operatorname{supp}(w^{\operatorname{appr}}-\widetilde{v})|\leq
n^{\widetilde{a}}$.

Intuitively, for $X\in\\{S,\Delta,\Gamma,\xi\\}$, $X^{\operatorname{new}}$
represents the difference between $w^{\operatorname{appr}}$ and the first-
level proxy $v$, this is the real “first derivative”, and it is what we need
to compute the output; $X$ represents the difference between the first-level
proxy $v$ and the second-level proxy $\widetilde{v}$, this is what we maintain
in the data structure, and we can think of it as an out-of-date “first
derivative”; $\partial X$ represents the difference between
$w^{\operatorname{appr}}$ and the second-level proxy $\widetilde{v}$. This
term is very sparse, and should be thought of as the “second derivative”. The
actual “first derivative” is the sum of the outdated “first derivative” plus
the “second derivative”. We can therefore split $X^{\operatorname{new}}$ as
$X+\partial X$ when computing the output. In this way we exploits both the
maintained members of the data structure and the _sparsity_ of $\partial X$,
compares to computing $X^{\operatorname{new}}$ from scratch.

We now turn to formalize the above intuition. Our data structure maintains the
second-level members $S\subseteq[n],\Delta\in\mathbb{R}^{n\times
n},\Gamma\in\mathbb{R}^{n\times n},\xi\in\mathbb{R}^{n}$. By design, our
update algorithm ensures that

$\displaystyle|S^{\operatorname{new}}|,\|\Delta^{\operatorname{new}}\|_{0},\|\Gamma^{\operatorname{new}}\|_{0},\|\xi^{\operatorname{new}}\|_{0},|S|,\|\Delta\|_{0},\|\Gamma\|_{0},\|\xi\|_{0}\leq
n^{a};\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ |\partial
S|,\|\partial\Delta\|_{0},\|\partial\Gamma\|_{0},\|\partial\xi\|_{0}\leq
n^{\widetilde{a}}.$

Now, recall that from Section 4.1.1 and Section 4.1.2, $r_{4}$ is computed as
follows

$\displaystyle
r_{4}=-\Big{(}Q_{S^{\operatorname{new}}}+R\underbrace{\Gamma^{\operatorname{new}}M_{S^{\operatorname{new}}}\Big{)}\cdot\Big{(}}_{3}\underbrace{BU^{\prime}}_{2}(C^{-1}+U^{\top}\underbrace{BU^{\prime}}_{2})^{-1}U^{\top}-I\Big{)}\cdot\underbrace{B\Big{(}(\beta_{2})_{S^{\operatorname{new}}}+(M_{S^{\operatorname{new}}})^{\top}\xi^{\operatorname{new}}\Big{)}}_{1},$

where $B=((\Delta_{S,S})^{-1}+M_{S,S})^{-1}$ and
$U^{\prime}CU^{\top}=((\Delta^{\operatorname{new}})^{-1}_{S^{\operatorname{new}},S^{\operatorname{new}}}+M_{S^{\operatorname{new}},S^{\operatorname{new}}})^{-1}-B$.
Observe that in the above formula the $n^{2a}$ term appears in the following
places:

  1. 1.

Multiplying a $n^{a}\times n^{a}$ matrix $B$ with a $n^{a}\times 1$ vector
$(\beta_{2})_{S^{\operatorname{new}}}+(M_{S^{\operatorname{new}}})^{\top}\xi^{\operatorname{new}}$.

  2. 2.

Multiplying a $n^{a}\times n^{a}$ matrix $B$ with a $n^{a}\times
n^{\widetilde{a}}$ matrix $U^{\prime}$.

  3. 3.

Multiplying a $n^{a}$-sparse diagonal matrix $\Gamma^{\operatorname{new}}$
with a $n\times n^{a}$ matrix $M_{S^{\operatorname{new}}}$ and then with a
$n^{a}\times 1$ vector that comes from later terms.

To remove these $n^{2a}$ terms, we maintain additional second-level data
structure members (precomputed matrix products):

$\displaystyle r_{4}=$ $\displaystyle\leavevmode\nobreak\
-\Big{(}Q_{S^{\operatorname{new}}}+\underbrace{R\Gamma M_{S}}_{F\text{:
member}}+R(\Gamma M_{\partial S\backslash S}+\partial\Gamma
M_{S^{\operatorname{new}}})\Big{)}\cdot\Big{(}(\underbrace{BU^{\prime}}_{U^{\mathrm{tmp}}}(C^{-1}+U^{\top}\underbrace{BU^{\prime}}_{U^{\mathrm{tmp}}})^{-1}U^{\top})-I\Big{)}\cdot$
$\displaystyle\leavevmode\nobreak\
\Big{(}\underbrace{B(\beta_{2})_{S}+B(M_{S})^{\top}\cdot\xi}_{\gamma_{1}\text{:
member}}+B(\beta_{2})_{\partial S\backslash S}+B(M_{\partial S\backslash
S})^{\top}\xi^{\operatorname{new}}+\underbrace{B(M_{S})^{\top}}_{E\text{:
member}}\partial\xi\Big{)}$

Each of the previous $n^{2a}$ terms are removed as follows.

  1. 1.

We maintain
$\gamma_{1}:=B(\beta_{2})_{S}+B(M_{S})^{\top}\xi\in\mathbb{R}^{n^{a}}$ so that
we only need to compute the difference

$(B(\beta_{2})_{S^{\operatorname{new}}}+B(M_{S^{\operatorname{new}}})^{\top}\xi^{\operatorname{new}})-\gamma_{1}=B(\beta_{2})_{\partial
S\backslash S}+B(M_{\partial S\backslash
S})^{\top}\xi^{\operatorname{new}}+B(M_{S})^{\top}\partial\xi.$

Since $(\beta_{2})_{\partial S\backslash S}$ is $n^{\widetilde{a}}$-sparse,
and $(M_{\partial S\backslash S})^{\top}$ only has $n^{\widetilde{a}}$ non-
zero rows, the first two terms $B(\beta_{2})_{\partial S\backslash S}$ and
$B(M_{\partial S\backslash S})^{\top}\xi^{\operatorname{new}}$ can both be
computed in $O(n^{a+\widetilde{a}})$ time. For the third term, we also
maintain $E:=B(M_{S})^{\top}\in\mathbb{R}^{n^{a}\times n^{a}}$. Since
$\partial\xi$ is also $n^{\widetilde{a}}$-sparse, it takes
$O(n^{a+\widetilde{a}})$ time to compute this term as well.

  2. 2.

The construction of $U^{\prime}$ has the following property:
$BU^{\prime}=[B_{(\partial S\backslash S)},B_{\partial S},BM_{S,{(\partial
S\backslash S)}}].$ Since we already maintain
$E:=B(M_{S})^{\top}\in\mathbb{R}^{n^{a}\times n}$, we define a local variable
$U^{\mathrm{tmp}}:=[B_{(\partial S\backslash S)},B_{\partial S},E_{(\partial
S\backslash S)}]\in\mathbb{R}^{n^{a}\times 3n^{\widetilde{a}}}$. Then
$U^{\mathrm{tmp}}=BU^{\prime}$, and it can be computed in the same time as its
size, which is $O(n^{a+\widetilde{a}})$.

  3. 3.

We maintain $F:=R\Gamma M_{S}\in\mathbb{R}^{\sqrt{n}\times n^{a}}$, and the
difference is
$R\Gamma^{\operatorname{new}}M_{S^{\operatorname{new}}}-F=R(\Gamma M_{\partial
S\backslash S}+\partial\Gamma M_{S^{\operatorname{new}}}).$ Multiplying $F$
with a $n^{a}\times 1$ vector that comes from later terms only takes
$n^{1/2+a}<n^{1.5}$ time. Also since $M_{\partial S\backslash S}$ only has
$n^{\widetilde{a}}$ non-empty columnns, and $\partial\Gamma$ is
$n^{\widetilde{a}}$-sparse, multiplying $(\Gamma M_{\partial S\backslash
S}+\partial\Gamma M_{S^{\operatorname{new}}})$ with the $n^{a}\times 1$ vector
that comes from later terms takes $n^{a+\widetilde{a}}$ time.

A full list of all the second-level members that we maintain is as follows:

$\displaystyle 1.S=\operatorname{supp}(\widetilde{v}-v),$ $\displaystyle
2.\Delta=\widetilde{V}-V,$ $\displaystyle
3.\Gamma=\sqrt{\widetilde{V}}-\sqrt{V},$ $\displaystyle
4.\xi=\sqrt{\widetilde{V}}\widetilde{g}-\sqrt{V}g,$ $\displaystyle
5.B=(\Delta_{S,S}^{-1}+M_{S,S})^{-1},$ $\displaystyle 6.E=B(M_{S})^{\top},$
$\displaystyle 7.F=R\Gamma M_{S},$ $\displaystyle
8.\gamma_{1}=B(\beta_{2})_{S}+B(M_{S})^{\top}\xi.$

Now whenever our data structure needs to multiply an $n^{a}\times n^{a}$
matrix with a $n^{a}\times 1$ vector, it is always the case that either the
vector is $n^{\widetilde{a}}$-sparse, or the matrix only has
$n^{\widetilde{a}}$ rows, so in both cases this operation takes
$O(n^{a+\widetilde{a}})$ time. We also avoid multiplying the $n^{a}\times
n^{a}$ matrix $B$ with the $n^{\widetilde{a}}\times n^{a}$ matrix $U^{\prime}$
directly by maintaining $E=B(M_{S})^{\top}$, and we can now extract
$U^{\mathrm{tmp}}=BU^{\prime}$ from $E$ efficiently. But still we need to
multiply a $n^{\widetilde{a}}\times n^{a}$ matrix $U^{\top}$ with a
$n^{a}\times n^{\widetilde{a}}$ matrix $U^{\mathrm{tmp}}$, which takes time
${\cal T}_{\mathrm{mat}}(n^{\widetilde{a}},n^{a},n^{\widetilde{a}})\leq
n^{a-\widetilde{a}}\cdot{\cal
T}_{\mathrm{mat}}(n^{\widetilde{a}},n^{\widetilde{a}},n^{\widetilde{a}})=n^{a+(\omega-1)\widetilde{a}}$.
This is our final running time for query, as presented in the last line of
Table 1.

Finally we remark that removing the last $n^{2a}$ term that stems from
computing $r_{2}$, can be done in an analogous way to the usage of
$\gamma_{1}$ for $r_{4}$ (i.e., by maintaining an additional member
$\gamma_{2}:=\Gamma M\xi$). We omit the formal details here.

Thus we finished the proof of how to get the
$O(n^{a+(\omega-1)\widetilde{a}})$ query time of Theorem 4.1.

###  4.2 Our Update Algorithm

Bounding the running time of our two-level update scheme requires both
algorithmic modifications and a more sophisticated amortized analysis than
that of [CLS19], in order to capture the cascading lazy updates process (as a
random process under the sketching of the CP). This section is organized as
follows. In Section 4.2.1 we describe the four update subroutines for
maintaining the two-level members for both $w$ and $h$, and present their
worst-case running time per call in Table 2. These subroutines correspond to
maintaining the LU-decomposition of the cascading updates algorithm for a
$K=3$ block matrix, described in Section 2 (see Figure 2). Section 4.2.2
describes the main algorithm deciding when to call each of these four
subroutines, using “two-level soft thresholding” (see Figure 6). In order to
synchronize two levels of soft thresholding, we introduce a new Adjust
function. Finally, Section 4.2.3 describes a potential-based amortized
analysis of our update algorithm, and the final amortized running time of the
four update subroutines is presented in Table 3.

####  4.2.1 Cascading updates subroutines

We now describe the four subroutines required to efficiently implement the
cascading lazy updates process for $K=3$ levels as described in the previous
section. Recall our data structure maintains two levels of proxies $v$ and
$\widetilde{v}$ for the input $w$: $v$ is the first-level proxy, and
$\widetilde{v}$ is the second-level proxy. $v$ keeps a larger distance of
$n^{a}$ with $w$ and is updated less frequently, while $\widetilde{v}$ keeps a
smaller distance of $n^{\widetilde{a}}$ with $w$ and is updated more
frequently. We define two subroutines to update $v$ and $\widetilde{v}$ (see
Figure 5 for illustration).

Level | Name | Time per call | Rank/Sparsity | Comment
---|---|---|---|---
1 | Matrix | ${\cal T}_{\mathrm{mat}}(n,n,k)$ | $k:=\|v^{\operatorname{new}}-v\|_{0}$ |  Update $v$ and $\widetilde{v}$ if $\|w^{\operatorname{appr}}-v\|_{0}\geq n^{a}$
2 | PartialMatrix | ${\cal T}_{\mathrm{mat}}(n,n^{a},\widetilde{k})$ | $\widetilde{k}:=\|\widetilde{v}^{\operatorname{new}}-\widetilde{v}\|_{0}$ |  Update $\widetilde{v}$ if $\|w^{\operatorname{appr}}-\widetilde{v}\|_{0}\geq n^{\widetilde{a}}$
1 | Vector | $pn+n^{2a}$ | $p:=\|g^{\operatorname{new}}-g\|_{0}$ |  Update $g$ and $\widetilde{g}$ if $\|h^{\operatorname{appr}}-g\|_{0}\geq n^{a}$
2 | PartialVector | $\widetilde{p}n+n^{2a}$ | $\widetilde{p}:=\|\widetilde{g}^{\operatorname{new}}-\widetilde{g}\|_{0}$ |  Update $\widetilde{g}$ if $\|h^{\operatorname{appr}}-\widetilde{g}\|_{0}\geq n^{\widetilde{a}}$
Table 2: Four update procedures

The first subroutine is PartialMatrixUpdate, which corresponds to second-level
updates we alluded to in Section 2: so long as the rank of the updates is
smaller than the second level threshold $n^{\widetilde{a}}$, it suffices to
only update the second level members as defined in Section 4.1.3. These
members are relatively cheap to update, and thus PartialMatrixUpdate has a
cost of ${\cal T}_{\mathrm{mat}}(n,n^{a},\widetilde{k})$ per call (third
column of Table 2). When the algorithm has exceeded the allowable changes in
$w$, we cascade to the first level update subroutine MatrixUpdate. This
subroutine must update all data structure members, and consequently is more
expensive: it has has a cost of ${\cal T}_{\mathrm{mat}}(n,n,k)$ per call
(third column of Table 2). The subroutines VectorUpdate and
PartialVectorUpdate play an analogous role for updating $h$. We proceed to
describe when to execute each of these subroutines.

Figure 5: An illustration of the two types of updates in our Update algorithm.
The blue leash connects $v$ with $\widetilde{v}$ and $w^{\operatorname{appr}}$
and is of length $n^{a}$. The pink leash connects $\widetilde{v}$ with
$w^{\operatorname{appr}}$ and is of length $n^{\widetilde{a}}$. (a) When all
leashes are loose, no update is performed. (b) The leash on $v$ is tight. In
this case we call MatrixUpdate. (c) The leash on $\widetilde{v}$ is tight. In
this case we call PartialMatrixUpdate.

####  4.2.2 Synchronizing two-level soft thresholding

For simplicity, we only focus on MatrixUpdate and PartialMatrixUpdate
(VectorUpdate and PartialVectorUpdate are analogous).

In order to bound the amortized cost of the two level updates, we use the
“soft thresholding” implicit777The soft thresholding is proposed by [CLS19],
and note that they embed it inside their update function since they only have
one update function. Since we use it four times, we give it an explicit name
SoftThreshold. in [CLS19] to determine when to invoke each update. The basic
idea is to use a smooth threshold (see Algorithm 7) to approximate the
discrete threshold shown in the last column of Table 2. Soft thresholding
ensures that there is a gap between errors of coordinates that are updated and
the errors of other coordinates, and is essential to guarantee a proper
decrease in potential. Our update algorithm invokes this subroutine
SoftThreshold twice – once for $\widetilde{v}$ and once for $v$ (see Figure
6). We now explain the 3 main components of combining the two SoftThreshold
subroutines (see Algorithm 2):

Algorithm 2 When to execute MatrixUpdate and PartialMatrixUpdate

1:$\widetilde{v}^{\operatorname{new}}\leftarrow\textsc{SoftThreshold}(y\leftarrow|w^{\operatorname{new}}/\widetilde{v}-1|,w^{\operatorname{new}},\widetilde{v},\epsilon_{\mathrm{mp}},n^{\widetilde{a}})$

2:$\textsc{Adjust}(\widetilde{v}^{\operatorname{new}},\epsilon_{\mathrm{mp}}/(100\log
n))$

3:if $\|\widetilde{v}^{\operatorname{new}}-v\|_{0}\geq n^{a}$ then

4:
$v^{\operatorname{new}}\leftarrow\textsc{SoftThreshold}(y\leftarrow|w^{\operatorname{new}}/v-1|+|w^{\operatorname{new}}/\widetilde{v}-1|,w^{\operatorname{new}},v,\frac{\epsilon_{\mathrm{mp}}}{100\log^{2}n},n^{a})$

5: $w^{\operatorname{appr}}\leftarrow v^{\operatorname{new}}$

6: MatrixUpdate()

7:else

8: $w^{\operatorname{appr}}\leftarrow\widetilde{v}^{\operatorname{new}}$

9: if $\|\widetilde{v}^{\operatorname{new}}-\widetilde{v}\|_{0}\geq
n^{\widetilde{a}}$ then

10: PartialMatrixUpdate()

11: end if

12:end if

##### Restoring threshold gaps via Adjust.

Our algorithm first computes $\widetilde{v}^{\operatorname{new}}$ to update
all the coordinates $i$ for which
$|w^{\operatorname{new}}_{i}/\widetilde{v}_{i}-1|$ exceeds the threshold
$\epsilon_{\mathrm{mp}}$. It then calls the Adjust function to restore all
updated coordinates $\widetilde{v}^{\operatorname{new}}_{i}$ whose new value
$\widetilde{v}^{\operatorname{new}}_{i}$ is within a distance of
$\epsilon_{\mathrm{mp}}/(100\log n)$ from $v_{i}$, back to the original value
$v_{i}$. Since $\widetilde{v}^{\operatorname{new}}$ is the new value of
$\widetilde{v}$, in this way we ensure that
$|\widetilde{v}_{i}-v_{i}|>\epsilon_{\mathrm{mp}}/(100\log n)$ for all
$i\in\operatorname{supp}(\widetilde{v}-v)$. Hence when
$\|\widetilde{v}-v\|_{0}$ exceeds its threshold, there is a large enough
decrease in our potential function
($\approx\|\widetilde{v}-v\|_{0}\cdot\epsilon_{\mathrm{mp}}/(100\log n)$) for
“charging” the update cost (of $v\leftarrow\widetilde{v}$). The purpose of
using a smaller threshold-error of $\epsilon_{\mathrm{mp}}/(100\log n)$ here
ensures that even when $\widetilde{v}^{\operatorname{new}}_{i}$ (which is the
new value of $\widetilde{v}$) is $v_{i}$ instead of
$w^{\operatorname{new}}_{i}$, the error still decreases by at least a
$(1-1/\log n)$ factor after updating
$\widetilde{v}\leftarrow\widetilde{v}^{\operatorname{new}}$.

##### Synchronizing the error function.

When updating $v$, we define the error as
$|w^{\operatorname{new}}/v-1|+|w^{\operatorname{new}}/\widetilde{v}-1|$ which
is a function of both $v$ and $\widetilde{v}$. This is because as long as
_one_ of $v_{i}$ and $\widetilde{v}_{i}$ is too far from
$w^{\operatorname{new}}_{i}$, we need to update both variables to be the same
as $w^{\operatorname{new}}_{i}$ .

##### Two error thresholds.

When updating $v$, we use a smaller error threshold of
$\epsilon_{\mathrm{mp}}/(100\log^{2}n)$ than that of the Adjust threshold.
This is because Adjust guarantees that
$v_{i}-\widetilde{v}_{i}\geq\epsilon_{\mathrm{mp}}/(100\log n)$ on all
coordinates $i$ for which $v_{i}\neq\widetilde{v}_{i}$, hence using an even
smaller threshold when updating $v$ ensures that all such coordinates are
counted as “error larger than threshold”. As such, after MatrixUpdate,
$\widetilde{v}$ is set back to be the same as $v$ on all coordinates.

Figure 6: Two-level vs. one-level soft thresholding: (a). [CLS19] has a single
level update scheme, implemented via one level of soft thresholding. (b). Our
Update algorithm has two cascading subroutines (MatrixUpdate,
PartialMatrixUpdate), each requires its own potential function and soft
threshold.

####  4.2.3 Amortized analysis based on high-order martingales

As noted in Section 2, the main source of amortization in the update algorithm
comes from the fact that the maintained variables in all levels are changing
slowly. More formally, the $j$-th iteration of the CP algorithm calls our data
structure with the following inputs: A vector $w^{(j+1)}\in\mathbb{R}^{n}$ and
a vector $h^{(j+1)}\in\mathbb{R}^{n}$ satisfying the following relative error
bounds (where the randomness is over the sketching matrix used to generate
$w^{(j+1)}$ and $h^{(j+1)}$):

$\displaystyle\|\operatorname*{{\mathbb{E}}}[w^{(j+1)}|w^{(j)}]/w^{(j)}-1\|_{2}\leavevmode\nobreak\
\leq O(1),\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\|\operatorname*{{\mathbb{E}}}[h^{(j+1)}|h^{(j)}]/h^{(j)}-1\|_{2}\leq O(1).$
(9)

The “evolution” of $w^{(j)}$ and $h^{(j)}$ is essentially a martingale process
with the above guarantee. Informally, this is true because we are using
_independent_ random sketches in each iteration. We analyze these random
processes by defining four different potential functions to capture our four
different update subroutines, as shown in Table 3. We next elaborate on the
careful design of potentials and what they aim to measure.

Level | Name | $\Phi\in\mathbb{R}$ |  $\phi\in\mathbb{R}^{n}$, $i>n^{a}$ (resp. $i>n^{\widetilde{a}}$) | Amortized
---|---|---|---|---
1 | Matrix | $\sum_{i=1}^{n}\phi_{i}\cdot(|w_{i}/v_{i}-1|+|w_{i}/\widetilde{v}_{i}-1|)$ | $\phi_{i}=i^{\frac{\omega-2}{1-a}-1}\cdot n^{-\frac{a(\omega-2)}{1-a}}$ | $n^{2-a/2}+n^{\omega-1/2}$
2 | PartialMatrix | $\sum_{i=1}^{n}\phi_{i}\cdot\left(|w_{i}/\widetilde{v}_{i}-1|\right)$ | $\phi_{i}=i^{\frac{a(\omega-2)}{a-\widetilde{a}}-1}\cdot n^{-\frac{a\widetilde{a}(\omega-2)}{a-\widetilde{a}}}$ | $n^{1+a-\widetilde{a}/2}$
1 | Vector | $\sum_{i=1}^{n}\phi_{i}\cdot\left(|h_{i}/g_{i}-1|+|h_{i}/\widetilde{g}_{i}-1|\right)$ | $\phi_{i}=1$ | $n^{1.5}$
2 | PartialVector | $\sum_{i=1}^{n}\phi_{i}\cdot(|h_{i}/\widetilde{g}_{i}-1|)$ | $\phi_{i}=i^{-1}$ | $n^{1.5}+n^{2a-\widetilde{a}/2}$
Table 3: The four different potential functions $\Phi$ are shown in the third
column. We always assume that the coordinates are sorted in decreasing order,
e.g., for PartialMatrix, we assume that
$|w_{i}/\widetilde{v}_{i}-1|\geq|w_{i+1}/\widetilde{v}_{i+1}-1|$. The vector
$\phi\in\mathbb{R}^{n}$ is non-increasing (in $i$): for level-1 subroutines,
$\phi_{i}:=n^{-a}$ when $i\leq n^{a}$, and for level-2 subroutines,
$\phi_{i}:=n^{-\widetilde{a}}$ for $i\leq n^{\widetilde{a}}$. The definition
of $\phi_{i}$ for $i>n^{a}$ (level 1) or $n^{\widetilde{a}}$ (level 2) are
shown in the fourth column. The vectors $\phi$ are designed to upper bound the
worst-case running time of the update procedures. The last column shows the
amortized cost of our four update subroutines.

1\. MatrixUpdate. The potential function for this subroutine is defined in row
1 of Table 3. Instead of splitting the change in the potential function into
“$w$ move” and “$v$ move” as [CLS19], we split the change into a “$w$ move”
part and a _“ $v$ and $\widetilde{v}$ move”_ part:

$\displaystyle\Phi^{\mathrm{mat}}_{j+1}-\Phi^{\mathrm{mat}}_{j}=\leavevmode\nobreak\
\text{($w$\leavevmode\nobreak\ move)}-\text{($v$\leavevmode\nobreak\
and\leavevmode\nobreak\ $\widetilde{v}$\leavevmode\nobreak\ move)}$

where the “$w$ move” part measures how much the potential can increase due to
the input changes from $w^{(j)}$ to $w^{(j+1)}$, and the “$v$ and
$\widetilde{v}$ move” part measures how much we can decrease the potential by
updating $v^{(j)},\widetilde{v}^{(j)}$ to $v^{(j+1)},\widetilde{v}^{(j+1)}$.
The formal details can be found in Section F.

Using Eq. (9) which upper bounds the expected relative error of $w^{(j+1)}$
with $w^{(j)}$, it is possible to upper bound the “$w$ move” term by
$O(\|\phi\|_{2})=O(n^{\omega-5/2}+n^{-a/2})$.

When entering MatrixUpdate, for some coordinate $i\in[k]$ (recall that
$k:=\|v^{\operatorname{new}}-v\|_{0}$, see Table 2), by design $v^{(j+1)}_{i}$
and $\widetilde{v}^{(j+1)}_{i}$ are both reset to $w^{(j+1)}_{i}$, hence the
potential
$\phi_{i}(|{w^{(j+1)}_{i}}/{v^{(j+1)}_{i}}-1|+|{w^{(j+1)}_{i}}/{\widetilde{v}^{(j+1)}_{i}}-1|)$
decreases to 0. This fact can be used to show that the “$v$ move” term in
$\Phi$ decreases by at least

$\displaystyle\Omega(k\cdot\phi_{k})\geq n^{-2}\cdot{\cal
T}_{\mathrm{mat}}(n,n,k).$

We also prove that PartialMatrixUpdate can only further decrease the “$v$
move” term. Since the “$v$ and $\widetilde{v}$ move” term is upper bounded by
the “$w$ move” term, and the cost per call of MatrixUpdate is ${\cal
T}_{\mathrm{mat}}(n,n,k)$, the amortized cost of MatrixUpdate per iteration is
bounded by $O(n^{\omega-1/2+o(1)}+n^{2-a/2+o(1)}).$

2\. PartialMatrixUpdate. The potential function for this subroutine is defined
in row 2 of Table 3. Once again, we split the change in the potential
function, this time into a “$w$ move” part and a “$\widetilde{v}$ move” part:

$\displaystyle\Phi_{j+1}-\Phi_{j}=\text{($w$\leavevmode\nobreak\
move)}-\text{($\widetilde{v}$\leavevmode\nobreak\ move)}$

The “$w$ move” term is upper bounded by
$O(\|\phi\|_{2})=O(n^{a\omega-5\widetilde{a}/2}+n^{-\widetilde{a}/2})$. To
lower bound the “$\widetilde{v}$ move” term, we observe that when entering
PartialMatrixUpdate, for any $i\in[\widetilde{k}]$ (recall that
$\widetilde{k}:=\|\widetilde{v}^{\operatorname{new}}-\widetilde{v}\|_{0}$, see
Table 2), the term
$|{w^{\operatorname{new}}_{i}}/{\widetilde{v}^{\operatorname{new}}_{i}}-1|$ is
decreased by at least a factor of $(1-1/\log n)$. This is where we use the
guarantees (and smaller threshold parameter) of Adjust and SoftThreshold for
$v$. Using this, we show that the “$\widetilde{v}$ move” term decreases by at
least

$\Omega(\widetilde{k}\cdot\phi_{\widetilde{k}})\geq n^{-1-a}\cdot{\cal
T}_{\mathrm{mat}}(n,n^{a},\widetilde{k}).$

Therefore, the amortized cost of PartialMatrixUpdate is
$O(n^{1+(\omega+1)a-5\widetilde{a}/2}+n^{1+a-\widetilde{a}/2})=O(n^{1+a-\widetilde{a}/2})$
since the cost per call of PartialMatrixUpdate is ${\cal
T}_{\mathrm{mat}}(n,n^{a},\widetilde{k})$.

3\. VectorUpdate. The potential function for this subroutine is defined in row
3 of Table 3. Note that the dominating term in the worst-case cost (per call)
of this subroutine is the $pn$ term (recall
$p:=\|g^{\operatorname{new}}-g\|_{0}$, see Table 2). Again, after amortization
in each iteration we have

$\sqrt{n}=\|\phi\|_{2}\geq\Omega(p\cdot\phi_{p})\geq n^{-1}(pn).$

Therefore, the amortized cost of VectorUpdate per iteration is $O(n^{1.5}).$

4\. PartialVectorUpdate. The potential function for this subroutine is defined
in row 4 of Table 3. Here, the dominating term in the worst-case cost (per
call) is the $n^{2a}$ term (Table 2). Since PartialVectorUpdate is invoked
only when $\widetilde{p}\geq n^{\widetilde{a}}$ (recall that
$\widetilde{p}:=\|\widetilde{g}^{\operatorname{new}}-\widetilde{g}\|_{0}$, see
Table 2), the amortized cost of the $j$-th iteration is
$n^{2a}\cdot\mathbf{1}_{\widetilde{p}_{j}>n^{\widetilde{a}}}$. Once again,
after amortization in each iteration we have

$n^{-\widetilde{a}}=\|\phi\|_{2}\geq\Omega(p\cdot\phi_{p})\geq\mathbf{1}_{\widetilde{p}_{j}>n^{\widetilde{a}}}.$

Thus the amortized cost of VectorUpdate per iteration is $n^{2a}\cdot
O(\|\phi\|_{2})=O(n^{2a-\widetilde{a}/2}).$

###  4.3 Putting it all together

Combining the query time $t_{q}=n^{a+(\omega-1)\widetilde{a}}$ (Eq. (7) in
Section 4.1) and the update time
$t_{u}=n^{\omega-1/2}+n^{2-a/2}+n^{1+a-\widetilde{a}/2}$ (see the last column
of Table 3), and since there are $O(\sqrt{n})$ iterations in total, we have
that the final running time of our LP algorithm is

$\displaystyle\sqrt{n}\cdot(t_{q}+t_{u})=n^{0.5+a+(\omega-1)\widetilde{a}}+n^{\omega}+n^{2.5-a/2}+n^{1.5+a-\widetilde{a}/2},$

matching the statement of Theorem 4.1.

Also note that in the ideal case where $\omega=2$ and $\alpha=1$, we can
choose $a=\frac{8}{9}$ and $\widetilde{a}=\frac{2}{3}$, and this leads to an
$O(n^{2+1/18})$ algorithm.

## References

  * [BCRL79] D. Bini, M. Capovani, F. Romani, and G. Lotti.  ${O}(n^{2.7799})$ complexity for $n\times n$ approximate matrix multiplication.  8(5):234–235, 1979.
  * [Ber24] Sergei Bernstein.  On a modification of chebyshev’s inequality and of the error formula of laplace.  Ann. Sci. Inst. Sav. Ukraine, Sect. Math, 1(4):38–49, 1924.
  * [BLSS20] Jan van den Brand, Yin Tat Lee, Aaron Sidford, and Zhao Song.  Solving tall dense linear programs in nearly linear time.  In STOC. https://arxiv.org/pdf/2002.02304.pdf, 2020.
  * [BNS19] Jan van den Brand, Danupon Nanongkai, and Thatchaphol Saranurak.  Dynamic matrix inverse: Improved algorithms and matching conditional lower bounds.  In 2019 IEEE 60th Annual Symposium on Foundations of Computer Science (FOCS), pages 456–480. IEEE, 2019.
  * [Bra20] Jan van den Brand.  A deterministic linear program solver in current matrix multiplication time.  In Proceedings of the Fourteenth Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 259–278. SIAM, 2020.
  * [CGLZ20] Matthias Christandl, François Le Gall, Vladimir Lysikov, and Jeroen Zuiddam.  Barriers for fast rectangular matrix multiplication.  In arXiv preprint. https://arxiv.org/pdf/2003.03019.pdf, 2020\.
  * [CKSU05] Henry Cohn, Robert Kleinberg, Balazs Szegedy, and Christopher Umans.  Group-theoretic algorithms for matrix multiplication.  In 46th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pages 379–388. IEEE, 2005.
  * [CLS19] Michael B Cohen, Yin Tat Lee, and Zhao Song.  Solving linear programs in the current matrix multiplication time.  In STOC. https://arxiv.org/pdf/1810.07896, 2019.
  * [Cop82] Don Coppersmith.  Rapid multiplication of rectangular matrices.  SIAM Journal on Computing, 11(3):467–471, 1982.
  * [Cop97] Don Coppersmith.  Rectangular matrix multiplication revisited.  Journal of Complexity, 13(1):42–49, 1997.
  * [CW82] Don Coppersmith and Shmuel Winograd.  On the asymptotic complexity of matrix multiplication.  SIAM Journal on Computing, 11(3):472–492, 1982.
  * [CW87] Don Coppersmith and Shmuel Winograd.  Matrix multiplication via arithmetic progressions.  In Proceedings of the nineteenth annual ACM symposium on Theory of computing (STOC), pages 1–6. ACM, 1987.
  * [CW13] Kenneth L. Clarkson and David P. Woodruff.  Low rank approximation and regression in input sparsity time.  In Symposium on Theory of Computing Conference (STOC), pages 81–90. https://arxiv.org/pdf/1207.6365, 2013.
  * [Dan47] George B Dantzig.  Maximization of a linear function of variables subject to linear inequalities.  Activity analysis of production and allocation, 13:339–347, 1947\.
  * [FS89] Michael Fredman and Michael Saks.  The cell probe complexity of dynamic data structures.  In Proceedings of the twenty-first annual ACM symposium on Theory of computing, pages 345–354, 1989.
  * [GLPS10] Anna C Gilbert, Yi Li, Ely Porat, and Martin J Strauss.  Approximate sparse recovery: optimizing time and measurements.  SIAM Journal on Computing 2012 (A preliminary version of this paper appears in STOC 2010), 41(2):436–453, 2010.
  * [GU18] Francois Le Gall and Florent Urrutia.  Improved rectangular matrix multiplication using powers of the coppersmith-winograd tensor.  In Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms(SODA), pages 1029–1046. https://arxiv.org/pdf/1708.05622.pdf, 2018.
  * [HKNS15] Monika Henzinger, Sebastian Krinninger, Danupon Nanongkai, and Thatchaphol Saranurak.  Unifying and strengthening hardness for dynamic problems via the online matrix-vector multiplication conjecture.  In Proceedings of the forty-seventh annual ACM symposium on Theory of computing (STOC), pages 21–30, 2015.
  * [Kar84] Narendra Karmarkar.  A new polynomial-time algorithm for linear programming.  In Proceedings of the sixteenth annual ACM symposium on Theory of computing (STOC), pages 302–311. ACM, 1984.
  * [Kha80] Leonid G Khachiyan.  Polynomial algorithms in linear programming.  USSR Computational Mathematics and Mathematical Physics, 20(1):53–72, 1980.
  * [KM72] Victor Klee and George J Minty.  How good is the simplex algorithm.  Inequalities, 3(3):159–175, 1972.
  * [KN12] Daniel M. Kane and Jelani Nelson.  Sparser johnson-lindenstrauss transforms.  In SODA, pages 1195–1206, 2012.
  * [LDFU13] Yichao Lu, Paramveer Dhillon, Dean P Foster, and Lyle Ungar.  Faster ridge regression via the subsampled randomized hadamard transform.  In Advances in neural information processing systems, pages 369–377, 2013.
  * [LG14] François Le Gall.  Powers of tensors and fast matrix multiplication.  In Proceedings of the 39th international symposium on symbolic and algebraic computation (ISSAC), pages 296–303. ACM, https://arxiv.org/pdf/1401.7714.pdf, 2014.
  * [LNNT16] Kasper Green Larsen, Jelani Nelson, Huy L Nguyen, and Mikkel Thorup.  Heavy hitters via cluster-preserving clustering.  In 57th Annual Symposium on Foundations of Computer Science (FOCS), pages 61–70. IEEE, https://arxiv.org/pdf/1604.01357, 2016.
  * [LS14] Yin Tat Lee and Aaron Sidford.  Path finding methods for linear programming: Solving linear programs in ${O}(\sqrt{rank})$ iterations and faster algorithms for maximum flow.  In 55th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pages 424–433. https://arxiv.org/pdf/1312.6677.pdf, https://arxiv.org/pdf/1312.6713.pdf, 2014.
  * [LS15] Yin Tat Lee and Aaron Sidford.  Efficient inverse maintenance and faster algorithms for linear programming.  In 56th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pages 230–249. https://arxiv.org/pdf/1503.01752.pdf, 2015.
  * [LSZ19] Yin Tat Lee, Zhao Song, and Qiuyi Zhang.  Solving empirical risk minimization in the current matrix multiplication time.  In COLT. https://arxiv.org/pdf/1905.04447, 2019.
  * [Mad13] Aleksander Madry.  Navigating central path with electrical flows: From flows to matchings, and back.  In 2013 IEEE 54th Annual Symposium on Foundations of Computer Science (FOCS), pages 253–262. IEEE, 2013.
  * [NN94] Yurii Nesterov and Arkadii Nemirovskii.  Interior-point polynomial algorithms in convex programming, volume 13.  Siam, 1994.
  * [Pan78] V Ya Pan.  Strassen’s algorithm is not optimal trilinear technique of aggregating, uniting and canceling for constructing fast algorithms for matrix operations.  In 19th Annual Symposium on Foundations of Computer Science (FOCS), pages 166–176. IEEE, 1978.
  * [PD06] Mihai Patrascu and Erik D Demaine.  Logarithmic lower bounds in the cell-probe model.  SIAM Journal on Computing, 35(4):932–963, 2006.
  * [PSW17] Eric Price, Zhao Song, and David P. Woodruff.  Fast regression with an ${\ell}_{\infty}$ guarantee.  In International Colloquium on Automata, Languages, and Programming (ICALP). https://arxiv.org/pdf/1705.10723.pdf, 2017.
  * [Ren88] James Renegar.  A polynomial-time algorithm, based on newton’s method, for linear programming.  Mathematical Programming, 40(1-3):59–93, 1988.
  * [Rom82] Francesco Romani.  Some properties of disjoint sums of tensors related to matrix multiplication.  SIAM Journal on Computing, 11(2):263–267, 1982.
  * [San04] Piotr Sankowski.  Dynamic transitive closure via dynamic matrix inverse.  In 45th Annual IEEE Symposium on Foundations of Computer Science, pages 509–517. IEEE, 2004.
  * [Sar06] Tamás Sarlós.  Improved approximation algorithms for large matrices via random projections.  In Proceedings of 47th Annual IEEE Symposium on Foundations of Computer Science (FOCS), 2006.
  * [Sch81] Arnold Schönhage.  Partial and total matrix multiplication.  SIAM Journal on Computing, 10(3):434–455, 1981.
  * [SM10] Piotr Sankowski and Marcin Mucha.  Fast dynamic transitive closure with lookahead.  Algorithmica, 56(2):180, 2010.
  * [Son19] Zhao Song.  Matrix Theory : Optimization, Concentration and Algorithms.  PhD thesis, The University of Texas at Austin, 2019.
  * [ST85] Daniel Dominic Sleator and Robert Endre Tarjan.  Self-adjusting binary search trees.  Journal of the ACM (JACM), 32(3):652–686, 1985.
  * [ST04] Daniel A Spielman and Shang-Hua Teng.  Nearly-linear time algorithms for graph partitioning, graph sparsification, and solving linear systems.  In Proceedings of the thirty-sixth annual ACM symposium on Theory of computing (STOC), pages 81–90, 2004.
  * [Str69] Volker Strassen.  Gaussian elimination is not optimal.  Numerische mathematik, 13(4):354–356, 1969.
  * [Str86] Volker Strassen.  The asymptotic spectrum of tensors and the exponent of matrix multiplication.  In 27th Annual Symposium on Foundations of Computer Science (FOCS), pages 49–54. IEEE, 1986.
  * [Str87] Gilbert Strang.  Karmarkar’s algorithm and its place in applied mathematics.  The Mathematical Intelligencer, 9(2):4–10, 1987.
  * [Str91] Volker Strassen.  Degeneration and complexity of bilinear maps: some asymptotic spectra.  J. reine angew. Math, 413:127–180, 1991.
  * [Vai87] Pravin M Vaidya.  An algorithm for linear programming which requires ${O}(((m+n)n^{2}+(m+n)^{1.5}n){L})$ arithmetic operations.  In 28th Annual IEEE Symposium on Foundations of Computer Science (FOCS), 1987.
  * [Vai89] Pravin M Vaidya.  Speeding-up linear programming using fast matrix multiplication.  In 30th Annual Symposium on Foundations of Computer Science (FOCS), pages 332–337. IEEE, 1989.
  * [Wil12] Virginia Vassilevska Williams.  Multiplying matrices faster than coppersmith-winograd.  In Proceedings of the forty-fourth annual ACM symposium on Theory of computing (STOC), pages 887–898. ACM, 2012.
  * [Woo49] Max A Woodbury.  The stability of out-input matrices.  Chicago, IL, 9, 1949.
  * [Woo50] Max A Woodbury.  Inverting modified matrices.  ., 1950.
  * [YTM94] Yinyu Ye, Michael J. Todd, and Shinji Mizuno.  An ${O}(\sqrt{n}{L})$-iteration homogeneous and self-dual linear programming algorithm.  Math. Oper. Res., 19(1):53–67, 1994.

###### Contents

  1. 1 Introduction
  2. 2 Bootstrapping low-rank updates via cascading lazy updates
  3. 3 Background
    1. 3.1 Recent developments in LP solvers
    2. 3.2 Optimization: The stochastic central-path algorithm
    3. 3.3 Data structures: Projection maintenance
  4. 4 Detailed Technical Overview
    1. 4.1 Our Query Algorithm
      1. 4.1.1 Technique for removing $n^{1+a}$
      2. 4.1.2 Technique for removing $n^{\omega a}$
      3. 4.1.3 Technique for removing $n^{2a}$
    2. 4.2 Our Update Algorithm
      1. 4.2.1 Cascading updates subroutines
      2. 4.2.2 Synchronizing two-level soft thresholding
      3. 4.2.3 Amortized analysis based on high-order martingales
    3. 4.3 Putting it all together
  5. A Preliminaries
  6. B Optimization
    1. B.1 Definitions
    2. B.2 Facts
    3. B.3 Bounding ${\delta}_{s}$, ${\delta}_{x}$, $\delta_{t}$, $\delta_{\Phi}$ and ${\delta}_{\mu}$
    4. B.4 Bounding $\overline{\mu}^{\operatorname{new}}-\overline{\mu}$
    5. B.5 Potential martingale
    6. B.6 Bounding the movement of $\overline{w}$
    7. B.7 Bounding the movement of $\overline{\mu}$
    8. B.8 One step of central path
  7. C Data structure : preliminary
    1. C.1 Preliminary and Definitions
    2. C.2 Facts
    3. C.3 Main result
  8. D Data structure : correctness
    1. D.1 Correctness of Query
    2. D.2 Correctness of UpdateV and UpdateG
    3. D.3 Correctness of MatrixUpdate
    4. D.4 Correctness of PartialMatrixUpdate
    5. D.5 Correctness of VectorUpdate
    6. D.6 Correctness of PartialVectorUpdate
    7. D.7 Correctness of Initialize
  9. E Data structure : time per call
    1. E.1 Sparsity guarantees
    2. E.2 Running time of Query
    3. E.3 Running time of MatrixUpdate
    4. E.4 Running time of PartialMatrixUpdate
    5. E.5 Running time of VectorUpdate
    6. E.6 Running time of PartialVectorUpdate
    7. E.7 Running time of Initialize
  10. F Data structure : amortized time
    1. F.1 Definitions and Preliminaries
    2. F.2 Facts based on Adjust and two level of SoftThreshold
    3. F.3 Amortized analysis for MatrixUpdate
      1. F.3.1 Definitions
      2. F.3.2 Main result
      3. F.3.3 $w$ move
      4. F.3.4 $v,\widetilde{v}$ move
      5. F.3.5 $\ell_{2}$-norm of $g$
    4. F.4 Amortized analysis for PartialMatrixUpdate
      1. F.4.1 Definitions
      2. F.4.2 Main result
      3. F.4.3 $w$ move
      4. F.4.4 $\widetilde{v}$ move
      5. F.4.5 $\ell_{2}$-norm of $g$
    5. F.5 Amortized analysis for VectorUpdate
    6. F.6 Amortized analysis for PartialVectorUpdate
    7. F.7 Potential function $\psi$
  11. G Combining data structure with optimization
  12. H Multi-level with more details
    1. H.1 LU-decomposition of Woodbury identity when $K=3$
    2. H.2 Online low-rank inverse and the oMV conjecture.
    3. H.3 Optimizing the parameters of Eq. (3)
  13. I A feasible algorithm
    1. I.1 Analysis
    2. I.2 Correctness of feasible algorithm
    3. I.3 Bounding $x$ and $\overline{x}$
    4. I.4 Running time of feasible data structure
  14. J History of Matrix Multiplication and LP

## Appendix

##  Appendix A Preliminaries

Throughout this paper when considering the linear program $\min_{Ax=b,x\geq
0}c^{\top}x$, we always assume that the input matrix $A\in\mathbb{R}^{d\times
n}$ is full-rank, and $d=\Omega(n)$, $d\leq n$. For convenience we also assume
that $n\geq 10$. In our algorithm we use the standard transformation of LP by
[YTM94] such that it is easy to obtain the initial $x$ and $s$ for the
transformed LP. We refer the readers to [YTM94] and [CLS19] for more details.

##### Standard notations.

For a positive integer $n$, we denote $[n]=\\{1,2,\cdots,n\\}$.

For a positive integer $n$, we use $I_{n}$ to denote the identity matrix of
size $n\times n$. We use standard definitions of hyperbolic functions
$\sinh(x)=\frac{e^{x}-e^{-x}}{2}$, $\cosh(x)=\frac{e^{x}+e^{-x}}{2}$.

For a vector $v\in\mathbb{R}^{n}$, we use the standard definition of
$\ell_{p}$ norms: $\forall p\geq 1$,
$\|v\|_{p}=(\sum_{i=1}^{n}|v_{i}|^{p})^{1/p}$. Specially,
$\|v\|_{\infty}=\max_{i\in[n]}|v_{i}|$.

A vector $v\in\mathbb{R}^{n}$ is called $k-$sparse if at most $k$ entries in
$v$ is non-zero.

For any random variable $x$, we use $\operatorname*{{\mathbb{E}}}[x]$ to
denote its expectation, and we use $\operatorname*{{\bf{Var}}}[x]$ to denote
its variance. We define $\operatorname*{{\bf{Sup}}}[x]$ to be the
deterministic maximum of $x$.

##### Approximations.

For any function $f$, we define
$\widetilde{O}(f)=f\cdot\operatorname{poly}\log(f)$, and $O^{*}(f)=f\cdot
f^{o(1)}$.

For vectors $a,b\in\mathbb{R}^{n}$ and accuracy parameter $\epsilon\in(0,1)$,
we use $a\approx_{\epsilon}b$ to denote that $(1-\epsilon)b_{i}\leq
a_{i}\leq(1+\epsilon)b_{i}$, for any $i\in[n]$. For constant $t$, we use
$a\approx_{\epsilon}t$ to denote that $(1-\epsilon)t\leq
a_{i}\leq(1+\epsilon)t$, for any $i\in[n]$.

##### Coordinate-wise operations.

For a vector $x\in\mathbb{R}^{n}$ and $s\in\mathbb{R}^{n}$, we denote
$xs\in\mathbb{R}^{n}$ as a length-$n$ vector whose $i$-th coordinate is
$(xs)_{i}=x_{i}s_{i}$, $\forall i\in[n]$. Similarly, we also define other
scalar operations on vectors as coordinate-wise operations.

For a scalar function $f:\mathbb{R}\rightarrow\mathbb{R}$ and a vector
$x\in\mathbb{R}^{n}$, we define
$f(x)=[f(x_{1}),f(x_{2}),\dots,f(x_{n})]^{\top}$.

##### Upper case as diagonal matrix.

Given vectors $x,s\in\mathbb{R}^{n}$, we use $X\in\mathbb{R}^{n\times n}$ and
$S\in\mathbb{R}^{n\times n}$ to denote the diagonal matrix of those two
vectors. We use $X/S$ to denote the diagonal matrix there the $i$-th entry on
the diagonal is $(X/S)_{i,i}=x_{i}/s_{i}$, $\forall i\in[n]$. Similarly, we
extend other scalar operations to diagonal matrix.

##### Matrix and vectors with subscripts.

For any matrix $M\in\mathbb{R}^{m\times n}$ where $m>1$ and $n>1$, and any
subset $S\subseteq[n]$, we define $M_{S}\in\mathbb{R}^{m\times|S|}$ to be the
submatrix of $M$ that only has columns in $S$.

For any subsets $S_{1}\subseteq[m],S_{2}\subseteq[n]$, we also define
$M_{S_{1},S_{2}}\in\mathbb{R}^{|S_{1}|\times|S_{2}|}$ to be the submatrix of
$M$ that only has rows in $S_{1}$ and columns in $S_{2}$.

For any vector $v\in\mathbb{R}^{n\times 1}$ where $n>1$, and any subset
$S\subseteq[n]$, we define $v_{S}\in\mathbb{R}^{|S|\times 1}$ to be the
subvector of $M$ that only has the entries in $S$.

##### Fast Matrix Multiplication.

We use $\omega$ to denote the exponent of matrix multiplication, which is
defined as the infimum number such that the multiplication of two $n\times n$
matrices can be done in $O(n^{\omega})$ time. We use $\alpha$ to denote the
the dual exponent of matrix multiplication, which is defined as the supremum
over all $a\geq 0$ such that the multiplication of a $n\times n^{a}$ matrix
with another $n^{a}\times n$ matrix can be done in $O(n^{2+o(1)})$ time.

We denote ${\cal T}_{\mathrm{mat}}(n,k,r)$ to be the time needed to multiply a
$n\times k$ matrix with a $k\times r$ matrix. ${\cal T}_{\mathrm{mat}}(n,k,r)$
has the following property.

######  Lemma A.1 ([Str91, GU18, CGLZ20]).

$\displaystyle{\cal T}_{\mathrm{mat}}(n,k,r)=O({\cal
T}_{\mathrm{mat}}(n,r,k))=O({\cal T}_{\mathrm{mat}}(k,n,r)).$

##### Woodbury identity.

We use the Woodbury matrix identity to calculate the inverse of a matrix $M$
under low-rank updates.

######  Fact A.2 (Woodbury matrix identity, [Woo49, Woo50]).

For matrices $M\in\mathbb{R}^{n\times n}$, $U\in\mathbb{R}^{n\times k}$,
$C\in\mathbb{R}^{k\times k}$, $V\in\mathbb{R}^{k\times n}$,

$\displaystyle(M+UCV)^{-1}=M^{-1}-M^{-1}U(C^{-1}+VM^{-1}U)^{-1}VM^{-1}.$

##### Probability tools.

We use the following well-known probability tools.

######  Lemma A.3 (Bernstein Inequality, [Ber24]).

Let $X_{1},\cdots,X_{n}$ be independent zero-mean random variables. Suppose
that $|X_{i}|\leq M$ almost surely, for all $i$. Then, for all $t>0$,

$\displaystyle\Pr\left[\sum_{i=1}^{n}X_{i}>t\right]\leq\exp\left(-\frac{t^{2}/2}{\sum_{j=1}^{n}\operatorname*{{\mathbb{E}}}[X_{j}^{2}]+Mt/3}\right).$

##### Central Path Method.

The stochastic central path method of [CLS19] consider the following LP with
$A\in\mathbb{R}^{d\times n}$: $\min_{Ax=b,x\geq 0}c^{\top}x$, and its dual:
$\max_{A^{\top}y\leq c}b^{\top}y$. By defining $s\in\mathbb{R}^{n}$ as the
slack variables, the optimal solution of the above LP must satisfy the
following constraints:

$\displaystyle x_{i}s_{i}=$ $\displaystyle\leavevmode\nobreak\
0,\leavevmode\nobreak\ \forall i,$ $\displaystyle Ax=$
$\displaystyle\leavevmode\nobreak\ b,$ $\displaystyle A^{\top}y+s=$
$\displaystyle\leavevmode\nobreak\ c,$ $\displaystyle x_{i},s_{i}\geq$
$\displaystyle\leavevmode\nobreak\ 0,\leavevmode\nobreak\ \forall i.$

where $\sum_{i=1}^{n}x_{i}s_{i}$ is the duality gap, and the other three
equations are the feasibility constraints. In the central path method, the
duality gap is parameterized by $t$ that decreases by a factor of
$1-O(\frac{1}{\sqrt{n}})$ in each iteration. Let $\mu:=xs$. [CLS19] designed a
potential function $\Psi(\mu/t-1)$ such that as long as
$\Psi(\mu/t-1)\leq\operatorname{poly}(n)$, $\mu\approx t$ is satisfied. The
change $\delta_{\mu}$ of $\mu$ has two parts: a term $-\frac{\mu}{\sqrt{n}}$
to decrease $\mu$, and a term $-\nabla\Psi(\mu/t-1)$ to bound the potential
function.

Given $\delta_{\mu}$, the changes $\delta_{x}$ and $\delta_{s}$ should satisfy
the following constraints (the second-order term $\delta_{x}\cdot\delta_{s}$
is ignored):

$\displaystyle X\delta_{s}+S\delta_{x}=$ $\displaystyle\leavevmode\nobreak\
\delta_{\mu},$ (10) $\displaystyle A\delta_{x}=$
$\displaystyle\leavevmode\nobreak\ 0,$ $\displaystyle
A^{\top}\delta_{y}+\delta_{s}=$ $\displaystyle\leavevmode\nobreak\ 0.$

The unique solution of the above linear equations is

$\displaystyle\delta_{x}=$ $\displaystyle\leavevmode\nobreak\
\frac{X}{\sqrt{XS}}(I-P)\frac{1}{\sqrt{XS}}\delta_{\mu},$
$\displaystyle\delta_{s}=$ $\displaystyle\leavevmode\nobreak\
\frac{S}{\sqrt{XS}}P\frac{1}{\sqrt{XS}}\delta_{\mu},$

where
$P=\sqrt{\frac{X}{S}}A^{\top}(A\frac{X}{S}A^{\top})^{-1}A\sqrt{\frac{X}{S}}$
is a projection matrix. Thus the stochastic central path method is summarized
as the Algorithm 1 presented in Section 3.2.

##  Appendix B Optimization

In this section we provide the error analysis of central path method. The
meaning of $x$, $s$, and the deviation of $\delta_{x}$, $\delta_{s}$ are
standard (see the Central Path Method paragraph of Section A).

The proof of this whole paper is induction based. The induction hypothesis
ensures that at the beginning of each iteration Assumption B.5 is satisfied.
In this section we use this induction hypothesis to prove the guarantees of
the central path method. Later in Section G we combine the central path
guarantees with the data structure guarantees to prove that Assumption B.5 is
still satisfied at the end of each iteration. More specifically in this
section we prove the following:

  1. 1.

Two guarantees Lemma B.28 and B.29 that are needed by the data structure given
in Section D. Then the properties of the data structure prove that Part 1 of
Assumption B.5 is still satisfied.

  2. 2.

An upper bound on the potential function (Lemma B.27) which proves that Part 2
of Assumption B.5 is still satisfied.

Lemma | Section | Comment
---|---|---
Lemma B.16 | Section B.3 | Bounding $\delta_{x},\delta_{s},\delta_{\mu},\delta_{t},\delta_{\Phi}$
Lemma B.21 | Section B.4 | Bounding $\mu^{\operatorname{new}}-\mu$
Lemma B.27 | Section B.5 | Bounding the expectation of potential function
Lemma B.28 | Section B.6 | Bounding the movement of ${w}$
Lemma B.29 | Section B.7 | Bounding the movement of ${\mu}$
Table 4: Summary of this section

###  B.1 Definitions

Some of the definitions are used as variables in the algorithm, some of the
definitions are used only for analysis, and some of the definitions are used
for both.

######  Assumption B.1.

We use the following two error parameters $\epsilon$ and
$\epsilon_{\mathrm{mp}}$:

$\displaystyle\epsilon\in(0,10^{-4}),\text{\leavevmode\nobreak\
and\leavevmode\nobreak\ }\epsilon_{\mathrm{mp}}\in(0,10^{-4}).$

We use the same potential function as [CLS19, LSZ19, Bra20],

######  Definition B.2 (Potential function).

For a parameter $\lambda>\log n$, we define function
$\Phi_{\lambda}:\mathbb{R}^{n}\rightarrow\mathbb{R}$:

$\displaystyle\Phi_{\lambda}(r):=\sum_{i=1}^{n}\cosh(\lambda r_{i}).$

######  Definition B.3 (Overline version of parameters).

At the beginning of each iteration, we have $t\in\mathbb{R}$, and
$\overline{x},\overline{s}\in\mathbb{R}^{n}$ from last iteration.
$t^{\operatorname{new}}\in\mathbb{R}$ is the new value of $t$. We define
$\overline{w}\in\mathbb{R}^{n}$ and $\overline{\mu}\in\mathbb{R}^{n}$:

$\displaystyle\overline{w}:=\overline{x}/\overline{s},\leavevmode\nobreak\
\overline{\mu}:=\overline{x}\cdot\overline{s}.$

We define overline version of projection matrix
$\overline{P}\in\mathbb{R}^{n\times n}$:

$\displaystyle\overline{P}:=\sqrt{\overline{W}}A^{\top}(A\overline{W}A^{\top})^{-1}A\sqrt{\overline{W}}.$

The change in $\overline{\mu}$ consists of two parts $\overline{\delta}_{t}$
and $\overline{\delta}_{\Phi}$:

$\displaystyle\overline{\delta}_{t}:=$ $\displaystyle\leavevmode\nobreak\
\Big{(}\frac{t^{\operatorname{new}}}{t}-1\Big{)}\overline{\mu},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\overline{\delta}_{\Phi}:=\leavevmode\nobreak\
-\frac{\epsilon}{2}t^{\operatorname{new}}\cdot\frac{\nabla\Phi_{\lambda}(\overline{\mu}/t-1)}{\|\nabla\Phi_{\lambda}(\overline{\mu}/t-1)\|_{2}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\overline{\delta}_{\mu}:=\leavevmode\nobreak\
\overline{\delta}_{t}+\overline{\delta}_{\Phi}.$

The changes to $\overline{x}$ and $\overline{s}$ is computed from
$\overline{\delta}_{\mu}$:

$\displaystyle\overline{\delta}_{x}:=\leavevmode\nobreak\
\frac{\overline{X}}{\sqrt{\overline{X}\overline{S}}}(I-\overline{P})\frac{1}{\sqrt{\overline{X}\overline{S}}}\overline{\delta}_{\mu},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\overline{\delta}_{s}:=\leavevmode\nobreak\
\frac{\overline{S}}{\sqrt{\overline{X}\overline{S}}}\overline{P}\frac{1}{\sqrt{\overline{X}\overline{S}}}\overline{\delta}_{\mu}.$

The data structure will take $\overline{w}\in\mathbb{R}^{n}$ and
$\overline{\mu}\in\mathbb{R}^{n}$ as inputs, and output some
$\widetilde{w}\in\mathbb{R}^{n}$ and $\widetilde{\mu}\in\mathbb{R}^{n}$ such
that
$\widetilde{w}\approx_{\epsilon_{\mathrm{mp}}}\overline{w},\leavevmode\nobreak\
\widetilde{\mu}\approx_{\epsilon_{\mathrm{mp}}}\overline{\mu}$.

######  Definition B.4 (Tilde version of parameters).

For a given $\widetilde{w}\in\mathbb{R}^{n}$ and
$\widetilde{\mu}\in\mathbb{R}^{n}$ that is returned by the data structure, we
define $\widetilde{x}\in\mathbb{R}^{n}$ and $\widetilde{s}\in\mathbb{R}^{n}$:

$\displaystyle\widetilde{x}:=\sqrt{\widetilde{w}\widetilde{\mu}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\widetilde{s}:=\sqrt{\widetilde{\mu}/\widetilde{w}}.$

Note that $\widetilde{x}/\widetilde{s}=\widetilde{w}$, and
$\widetilde{x}\widetilde{s}=\widetilde{\mu}$. We define tilde version of
projection matrix $\widetilde{P}\in\mathbb{R}^{n\times n}$:

$\displaystyle\widetilde{P}:=\sqrt{\widetilde{W}}A^{\top}(A\widetilde{W}A^{\top})^{-1}A\sqrt{\widetilde{W}}.$

Further, we define $\widetilde{\delta}_{t}$, $\widetilde{\delta}_{\Phi}$,
$\widetilde{\delta}_{\mu}\in\mathbb{R}^{n}$:

$\displaystyle\widetilde{\delta}_{t}:=\leavevmode\nobreak\
(\frac{t^{\operatorname{new}}}{t}-1)\widetilde{\mu},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\widetilde{\delta}_{\Phi}:=\leavevmode\nobreak\
-\frac{\epsilon}{2}t^{\operatorname{new}}\cdot\frac{\nabla\Phi_{\lambda}(\widetilde{\mu}/t-1)}{\|\nabla\Phi_{\lambda}(\widetilde{\mu}/t-1)\|_{2}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\widetilde{\delta}_{\mu}:=\leavevmode\nobreak\
\widetilde{\delta}_{t}+\widetilde{\delta}_{\Phi}.$

And we define $\widetilde{\delta}_{x}$,
$\widetilde{\delta}_{s}\in\mathbb{R}^{n}$:

$\displaystyle\widetilde{\delta}_{x}:=\leavevmode\nobreak\
\frac{\widetilde{X}}{\sqrt{\widetilde{X}\widetilde{S}}}(I-\widetilde{P})\frac{1}{\sqrt{\widetilde{X}\widetilde{S}}}\widetilde{\delta}_{\mu},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\widetilde{\delta}_{s}:=\leavevmode\nobreak\
\frac{\widetilde{S}}{\sqrt{\widetilde{X}\widetilde{S}}}\widetilde{P}\frac{1}{\sqrt{\widetilde{X}\widetilde{S}}}\widetilde{\delta}_{\mu}.$

Given these definitions, we state the following important assumptions. We
assume they are satisfied in the beginning of each iteration, and use them to
prove the correctness of the algorithm in this iteration. Later we will verify
that these assumptions are always satisfied by induction on iterations. Note
that the second assumption is true if the potential function is bounded by
$\operatorname{poly}(n)$.

######  Assumption B.5.

We make the following assumptions:

  1. 1.

$\widetilde{\mu}\approx_{\epsilon_{\mathrm{mp}}}\overline{\mu}$,
$\widetilde{w}\approx_{\epsilon_{\mathrm{mp}}}\overline{w}$, where
$\widetilde{\mu}$ and $\widetilde{w}$ are returned by the data structure, and
$\overline{\mu}$, $\overline{w}$ are the input to the data structure.

  2. 2.

$\overline{\mu}\approx_{0.1}t$.

We further make the following definitions that make use of sketching matrices:

######  Definition B.6 (Hat version of parameters).

Let $b=o(n)$. For any $\widetilde{x}\in\mathbb{R}^{n}$,
$\widetilde{s}\in\mathbb{R}^{n}$, $\widetilde{P}\in\mathbb{R}^{n\times n}$,
$\widetilde{\delta}_{\mu}\in\mathbb{R}^{n}$, and a sketching matrix
$R\in\mathbb{R}^{b\times n}$, we define $\widehat{\delta}_{x}$,
$\widehat{\delta}_{s}\in\mathbb{R}^{n}$:

$\displaystyle\widehat{\delta}_{x}:=\leavevmode\nobreak\
\frac{\widetilde{X}}{\sqrt{\widetilde{X}\widetilde{S}}}(I-R^{\top}R\widetilde{P})\frac{1}{\sqrt{\widetilde{X}\widetilde{S}}}\widetilde{\delta}_{\mu},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\widehat{\delta}_{s}:=\leavevmode\nobreak\
\frac{\widetilde{S}}{\sqrt{\widetilde{X}\widetilde{S}}}R^{\top}R\widetilde{P}\frac{1}{\sqrt{\widetilde{X}\widetilde{S}}}\widetilde{\delta}_{\mu}.$

######  Remark B.7.

In our case, $R\in\mathbb{R}^{b\times n}$ is a subsampled randomized Hadamard
matrix. See more details in Definition B.10.

######  Definition B.8 (New versions of definition).

We use a superscript “new” to denote the corresponding variables at the
beginning of the next iteration. Specifically, we define
$t^{\operatorname{new}}\in\mathbb{R}$ as follows:

$\displaystyle t^{\operatorname{new}}:=(1-\frac{\epsilon}{3\sqrt{n}})t.$

We define
$\overline{\mu}^{\operatorname{new}},\overline{w}^{\operatorname{new}}\in\mathbb{R}^{n}$
as follows:

$\displaystyle\overline{\mu}^{\operatorname{new}}:=(\overline{x}+\widehat{\delta}_{x})\cdot(\overline{s}+\widehat{\delta}_{s}),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\overline{w}^{\operatorname{new}}:=(\overline{x}+\widehat{\delta}_{x})/(\overline{s}+\widehat{\delta}_{s}).$

The following facts directly follow from the definitions and Assumption B.5.

######  Fact B.9.

We have the following properties:

  1. 1.

$\widetilde{X}\widehat{\delta}_{s}+\widetilde{S}\widehat{\delta}_{x}=\widetilde{\delta}_{\mu}=\widetilde{\delta}_{t}+\widetilde{\delta}_{\Phi}$,

  2. 2.

$\widetilde{x}\approx_{2\epsilon_{\mathrm{mp}}}\overline{x}$,
$\widetilde{s}\approx_{2\epsilon_{\mathrm{mp}}}\overline{s}$,

  3. 3.

$\|\widetilde{\delta}_{t}\|_{2}\leq 0.5\epsilon t$,
$\|\widetilde{\delta}_{\Phi}\|_{2}\leq 0.5\epsilon t$,
$\|\widetilde{\delta}_{\mu}\|_{2}\leq\epsilon t$.

###### Proof.

Part 1. From the definition of $\widehat{\delta}_{x}$ and
$\widehat{\delta}_{s}$ (Definition B.6) we have that

$\displaystyle\widetilde{X}\widehat{\delta}_{s}+\widetilde{S}\widehat{\delta}_{x}=$
$\displaystyle\leavevmode\nobreak\
\frac{\widetilde{X}\widetilde{S}}{\sqrt{\widetilde{X}\widetilde{S}}}R^{\top}R\widetilde{P}\frac{1}{\sqrt{\widetilde{X}\widetilde{S}}}\widetilde{\delta}_{\mu}+\frac{\widetilde{S}\widetilde{X}}{\sqrt{\widetilde{X}\widetilde{S}}}(I-R^{\top}R\widetilde{P})\frac{1}{\sqrt{\widetilde{X}\widetilde{S}}}\widetilde{\delta}_{\mu}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\frac{\widetilde{X}\widetilde{S}}{\sqrt{\widetilde{X}\widetilde{S}}}I\frac{1}{\sqrt{\widetilde{X}\widetilde{S}}}\widetilde{\delta}_{\mu}=\widetilde{\delta}_{\mu}=\leavevmode\nobreak\
\widetilde{\delta}_{t}+\widetilde{\delta}_{\Phi}.$

Part 2. By part 1 and part 2 of Assumption B.5, we have that
$\widetilde{\mu}\approx_{\epsilon_{\mathrm{mp}}}\overline{\mu}$ and
$\widetilde{w}\approx_{\epsilon_{\mathrm{mp}}}\overline{w}$, so we have

$\displaystyle\widetilde{x}=\sqrt{\widetilde{w}\cdot\widetilde{\mu}}\leavevmode\nobreak\
\approx_{2\epsilon_{\mathrm{mp}}}\leavevmode\nobreak\
\sqrt{\overline{w}\cdot\overline{\mu}}=\sqrt{(\overline{x}/\overline{s})\cdot(\overline{x}\cdot\overline{s})}=\overline{x},$

where the first step follows from Definition B.4, the second step follows from
the fact that if $a\approx_{\epsilon_{\mathrm{mp}}}a^{\prime}$ and
$b\approx_{\epsilon_{\mathrm{mp}}}b^{\prime}$, then
$ab\approx_{2\epsilon_{\mathrm{mp}}}a^{\prime}b^{\prime}$, and the third step
follows from Definition B.3.

Using a similar argument, we also have that
$\widetilde{s}\approx_{2\epsilon_{\mathrm{mp}}}\overline{s}$.

Part 3. We upper bound $\|\widetilde{\delta}_{t}\|_{2}$ as follows:

$\displaystyle\|\widetilde{\delta}_{t}\|_{2}=\Big{\|}\Big{(}\frac{t^{\operatorname{new}}}{t}-1\Big{)}\widetilde{\mu}\Big{\|}_{2}=\Big{\|}\frac{\epsilon}{3\sqrt{n}}\widetilde{\mu}\Big{\|}_{2}\leq(1+\epsilon_{\mathrm{mp}})\Big{\|}\frac{\epsilon}{3\sqrt{n}}\overline{\mu}\Big{\|}_{2}\leq
1.1(1+\epsilon_{\mathrm{mp}})\Big{\|}\frac{\epsilon t}{3\sqrt{n}}{\bf
1}\Big{\|}_{2}\leq 0.5\epsilon t.$

where the first step follows from the definition of $\widetilde{\delta}_{t}$
(Definition B.4), the second step follows from the definition of
$t^{\operatorname{new}}$ (Definition B.8), the third step follows from
$\widetilde{\mu}\approx_{\epsilon_{\mathrm{mp}}}\overline{\mu}$ (Part 1 of
Assumption B.5), the forth step follows from $\overline{\mu}\approx_{0.1}t$
(Part 2 of Assumption B.5), and the last step follows from
$\epsilon_{\mathrm{mp}}\leq 10^{-4}$ (Assumption B.1).

Then we upper bound $\|\widetilde{\delta}_{\Phi}\|_{2}$ as follows:

$\displaystyle\|\widetilde{\delta}_{\Phi}\|_{2}=\|\frac{\epsilon}{2}t^{\operatorname{new}}\cdot\frac{\nabla\Phi_{\lambda}(\widetilde{\mu}/t-1)}{\|\nabla\Phi_{\lambda}(\widetilde{\mu}/t-1)\|_{2}}\|_{2}=\frac{\epsilon}{2}t^{\operatorname{new}}=\frac{\epsilon}{2}(1-\frac{\epsilon}{3\sqrt{n}})t\leq
0.5\epsilon t,$

where the first step follows from the definition of
$\widetilde{\delta}_{\Phi}$ (Definition B.4), and the third step follows from
the definition of $t^{\operatorname{new}}$ (Definition B.8).

Finally, we can use triangle inequality to upper bound

$\|\widetilde{\delta}_{\mu}\|_{2}\leq\|\widetilde{\delta}_{t}\|_{2}+\|\widetilde{\delta}_{\Phi}\|_{2}\leq
0.5\epsilon t+0.5\epsilon t\leq\epsilon t.$

∎

###  B.2 Facts

Random sketching matrices are usually used to give subspace embedding and
approximate matrix product. Instead of using subspace embedding [Sar06] and
approximate matrix product [KN12], LP solver requires a different version of
embedding, it was from [PSW17] implicitly and defined in [LSZ19] explicitly.

######  Definition B.10 (Coordinate-wise embedding).

Let $\Pi$ denote a distribution on $b\times n$ matrices $R$. We say $\Pi$ is
an $(\alpha,\beta,\delta)$-coordinate-wise embedding if for any fixed vector
$h\in\mathbb{R}^{n}$, the following properties hold:

$\displaystyle 1.$
$\displaystyle\operatorname*{{\mathbb{E}}}_{R\sim\Pi}[R^{\top}Rh]=h,$
$\displaystyle 2.$
$\displaystyle\operatorname*{{\mathbb{E}}}_{R\sim\Pi}[(R^{\top}Rh)_{i}^{2}]\leq
h_{i}^{2}+\frac{\alpha}{b}\|h\|_{2}^{2},$ $\displaystyle 3.$
$\displaystyle\Pr_{R\sim\Pi}\Big{[}|(R^{\top}Rh)_{i}-h_{i}|>\|h\|_{2}\frac{\beta}{\sqrt{b}}\Big{]}\leq\delta.$

######  Lemma B.11 (Lemma A.1 in [CLS19]).

Let $x$ and $y$ be (possibly dependent) random variables such that $|x|\leq
c_{x}$ and $|y|\leq c_{y}$ almost surely. Then, we have

$\displaystyle\operatorname*{{\bf{Var}}}[xy]\leq
2c_{x}^{2}\cdot\operatorname*{{\bf{Var}}}[y]+2c_{y}^{2}\cdot\operatorname*{{\bf{Var}}}[x].$

######  Fact B.12 (Gradient and Hessian of potential function).

Let $\Phi_{\lambda}(r)=\sum_{i=1}^{n}\cosh(\lambda r_{i})$ for some
$\lambda>0$. The gradient and the Hessian of $\Phi_{\lambda}(r)$ are

$\displaystyle\nabla\Phi_{\lambda}(r)=$ $\displaystyle\leavevmode\nobreak\
\lambda\cdot\big{(}\sinh(\lambda r_{1}),\sinh(\lambda
r_{2}),\cdots,\sinh(\lambda r_{n})\big{)}^{\top}\in\mathbb{R}^{n},$
$\displaystyle\nabla^{2}\Phi_{\lambda}(r)=$ $\displaystyle\leavevmode\nobreak\
\operatorname{diag}\big{(}\lambda^{2}\cosh(\lambda
r_{1}),\lambda^{2}\cosh(\lambda r_{2}),\cdots,\lambda^{2}\cosh(\lambda
r_{n})\big{)}\in\mathbb{R}^{n\times n}.$

######  Lemma B.13 (Basic properties of potential function, Lemma 4.12 in
[CLS19]).

Let $\Phi_{\lambda}(r)=\sum_{i=1}^{n}\cosh(\lambda r_{i})$ for some
$\lambda>0$. For any vector $r\in\mathbb{R}^{n}$,
1\. For any vector $\|v\|_{\infty}\leq 1/\lambda$, we have that

$\displaystyle\Phi_{\lambda}(r+v)\leq\Phi_{\lambda}(r)+\langle\nabla\Phi_{\lambda}(r),v\rangle+2\|v\|_{\nabla^{2}\Phi_{\lambda}(r)}^{2}.$

2\.
$\|\nabla\Phi_{\lambda}(r)\|_{2}\geq\frac{\lambda}{\sqrt{n}}(\Phi_{\lambda}(r)-n).$
3\. $\left(\sum_{i=1}^{n}\lambda^{2}\cosh^{2}(\lambda
r_{i})\right)^{1/2}\leq\lambda\sqrt{n}+\|\nabla\Phi_{\lambda}(r)\|_{2}.$

######  Lemma B.14 (Appendix E in [LSZ19]).

Let $R\in\mathbb{R}^{b\times n}$ denote a subsample randomized Hadamard
transform, then it gives
$(\alpha=1,\beta=O(\log(n/\delta)),\delta)$-Coordinate-wise Embedding
(Definition B.10).

We state another property of potential function

######  Lemma B.15 (Basic properties of potential function, general version of
Lemma 5.14 of [Bra20]).

Let $\Phi_{\lambda}(r)=\sum_{i=1}^{n}\cosh(\lambda r_{i})$ for some
$\lambda>0$. If $\|v\|_{\infty}\leq 1/(30\lambda)$, then we have

$\displaystyle\Big{\langle}\nabla\Phi_{\lambda}(r),-\frac{\nabla\Phi_{\lambda}(r+v)}{\|\nabla\Phi_{\lambda}(r+v)\|_{2}}\Big{\rangle}\leq-0.9\|\nabla\Phi_{\lambda}(r)\|_{2}+0.1\lambda\sqrt{n}.$

The proof is very similar to that of [Bra20], we put it here for completeness.

###### Proof.

We have

$\displaystyle\langle\nabla\Phi_{\lambda}(r),-\nabla\Phi_{\lambda}(r+v)\rangle=$
$\displaystyle\>-\langle\nabla\Phi_{\lambda}(r+v),\nabla\Phi_{\lambda}(r+v)\rangle+\langle\nabla\Phi_{\lambda}(r+v)-\nabla\Phi_{\lambda}(r),\nabla\Phi_{\lambda}(r+v)\rangle$
$\displaystyle\leq$
$\displaystyle\>-\|\nabla\Phi_{\lambda}(r+v)\|_{2}^{2}+\|\nabla\Phi_{\lambda}(r)-\nabla\Phi_{\lambda}(r+v)\|_{2}\cdot\|\nabla\Phi_{\lambda}(r+v)\|_{2}$

where the last step follows from $\langle
a,b\rangle\leq\|a\|_{2}\cdot\|b\|_{2}$. Then we have that

$\displaystyle\leavevmode\nobreak\
\Big{\langle}\nabla\Phi_{\lambda}(r),-\frac{\nabla\Phi_{\lambda}(r+v)}{\|\nabla\Phi_{\lambda}(r+v)\|_{2}}\Big{\rangle}$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
-\|\nabla\Phi_{\lambda}(r+v)\|_{2}+\|\nabla\Phi_{\lambda}(r)-\nabla\Phi_{\lambda}(r+v)\|_{2}$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
-\Big{(}\|\nabla\Phi_{\lambda}(r)\|_{2}-\|\nabla\Phi_{\lambda}(r)-\nabla\Phi_{\lambda}(r+v)\|_{2}\Big{)}+\|\nabla\Phi_{\lambda}(r)-\nabla\Phi_{\lambda}(r+v)\|_{2}$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
-\|\nabla\Phi_{\lambda}(r)\|_{2}+2\|\nabla\Phi_{\lambda}(r)-\nabla\Phi_{\lambda}(r+v)\|_{2},$
(11)

where the second step follows from the fact that
$\|a\|_{2}\geq\|b\|_{2}-\|b-a\|_{2}$.

Now we need to upper bound the norm
$\|\nabla\Phi_{\lambda}(r)-\nabla\Phi_{\lambda}(r+v)\|_{2}$. Note that
$\nabla\Phi_{\lambda}(x)_{i}=\lambda\sinh(\lambda x_{i})$ and
$\sinh(x)=(e^{x}-e^{-x})/2$. We have the following property for
$|\sinh(x+y)-\sinh(x)|$:

$\displaystyle|\sinh(x+y)-\sinh(x)|=$ $\displaystyle\leavevmode\nobreak\
|e^{x}\cdot e^{y}-e^{-x}\cdot e^{-y}-(e^{x}-e^{-x})|/2$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\
|e^{x}\cdot(e^{y}-1)+e^{-x}\cdot(1-e^{-y})|/2$ $\displaystyle\leq$
$\displaystyle\leavevmode\nobreak\
(e^{x}\cdot|e^{y}-1|+e^{-x}\cdot|1-e^{-y}|)/2$ $\displaystyle\leq$
$\displaystyle\leavevmode\nobreak\
(e^{x}+e^{-x})/2\cdot\max\\{|e^{y}-1|,|1-e^{-y}|\\}$ $\displaystyle\leq$
$\displaystyle\leavevmode\nobreak\
(e^{x}+e^{-x})/2\cdot(e^{|y|}-1)=\leavevmode\nobreak\ \cosh(x)(e^{|y|}-1),$
(12)

where the first and the last step follows from the definitions of $\sinh$ and
$\cosh$, the third step follows from the triangle inequality of absolute
values, and the fifth step follows from $e^{y}+e^{-y}\geq 2$.

Thus we can upper bound the difference as follows

$\displaystyle\|\nabla\Phi_{\lambda}(r)-\nabla\Phi_{\lambda}(r+v)\|_{2}=$
$\displaystyle\leavevmode\nobreak\ \lambda\|\sinh(\lambda(r+v))-\sinh(\lambda
r)\|_{2}$ $\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\lambda\|\cosh(\lambda r)(e^{|\lambda v|}-1)\|_{2}$ $\displaystyle\leq$
$\displaystyle\leavevmode\nobreak\ \lambda\|\cosh(\lambda
r)\|_{2}(e^{\lambda\|v\|_{\infty}}-1)$ $\displaystyle\leq$
$\displaystyle\leavevmode\nobreak\
(\lambda\sqrt{n}+\|\nabla\Phi_{\lambda}(r)\|_{2})(e^{\lambda\|v\|_{\infty}}-1),$
(13)

where the second step follows from Eq. (B.2), the third step follows from the
fact that $\|a\cdot b\|_{2}\leq\|a\|_{2}\|b\|_{\infty}$, and the last step
follows from Part 3 of Lemma B.13.

We then have

$\displaystyle(e^{\lambda\|v\|_{\infty}}-1)<\leavevmode\nobreak\
e^{1/30}-1\leq\leavevmode\nobreak\ 0.05,$ (14)

where the first step follows from $\|v\|_{\infty}\leq 1/(30\lambda)$.

Finally, this allows us to obtain

$\displaystyle\Big{\langle}\nabla\Phi_{\lambda}(r),-\frac{\nabla\Phi_{\lambda}(r+v)}{\|\nabla\Phi_{\lambda}(r+v)\|_{2}}\Big{\rangle}\leq$
$\displaystyle\leavevmode\nobreak\
-\|\nabla\Phi_{\lambda}(r)\|_{2}+2\|\nabla\Phi_{\lambda}(r)-\nabla\Phi_{\lambda}(r+v)\|_{2}$
$\displaystyle<$ $\displaystyle\leavevmode\nobreak\
-\|\nabla\Phi_{\lambda}(r)\|_{2}+0.1(\lambda\sqrt{n}+\|\nabla\Phi_{\lambda}(r)\|_{2})$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
-0.9\|\nabla\Phi_{\lambda}(r)\|_{2}+0.1\lambda\sqrt{n}$

where the first step follows from Eq. (11), and the second step follows from
Eq. (B.2) and Eq. (14). ∎

###  B.3 Bounding ${\delta}_{s}$, ${\delta}_{x}$, $\delta_{t}$,
$\delta_{\Phi}$ and ${\delta}_{\mu}$

Quantity | Bound | Stated in Lem. B.16 | Used by Lem. B.21
---|---|---|---
$\|\overline{s}^{-1}\widetilde{\delta}_{s}\|_{2}$, $\|\overline{x}^{-1}\widetilde{\delta}_{x}\|_{2}$ | $\epsilon$ | Part 1 | Part 1
$\|\operatorname*{{\mathbb{E}}}[\overline{s}^{-1}\widehat{\delta}_{s}]\|_{2}$, $\|\operatorname*{{\mathbb{E}}}[\overline{x}^{-1}\widehat{\delta}_{x}]\|_{2}$ | $\epsilon$ | Part 1 | Part 1
$\|\overline{\mu}^{-1}(\widetilde{\delta}_{t}-\overline{\delta}_{t})\|_{2}$ | $\epsilon_{\mathrm{mp}}\cdot\epsilon$ | Part 1 | Part 1
$\|\overline{\mu}^{-1}(\overline{\delta}_{t}+\widetilde{\delta}_{\Phi})\|_{2}$ | $\epsilon$ | Part 1 | Part 4
$\operatorname*{{\bf{Var}}}[\overline{x}_{i}^{-1}\widehat{\delta}_{x,i}],\operatorname*{{\bf{Var}}}[\overline{s}_{i}^{-1}\widehat{\delta}_{s,i}]$ | $\epsilon^{2}/b$ | Part 2 | Part 1, 2