A Deterministic Linear Program Solver in Current Matrix Multiplication Time

Introduction

Fast algorithms for solving linear programs have a long history in computer science. Solving linear programs was first proven to be in $P$ in 1979 by Khachiyan [Kha79]; and later Karmarkar [Kar84] found the first polynomial time algorithm that was feasible in practice. This initiated the long line of work of solving linear programs using interior point algorithms, motivated by the fact that many problems can be stated as linear programs and solved using efficient solvers. [Ren88, Vai87, Vai89b, Vai89a, Meg89, NN89, NN91, VA93, Ans96, NT97, Ans99, LS14, LS15]

For linear programs of the form $\min_{Ax=b,x\geq 0}c^{\top}x$ with $n$ variables, $d$ constraints and $nnz(A)$ non-zero entries, the current fastest algorithms are $\widetilde{O}(\sqrt{d}(nnz(A)+d^{2}))$ [LS14, LS15] Here $\widetilde{O}(\cdot)$ hides $\text{polylog}(n)$ and $\text{polylog}(1/\delta)$ terms. and $\widetilde{O}(n^{\omega})$-time [CLS18] The algorithm of [CLS18] runs in $O((n^{\omega}+n^{2.5-\alpha/2+o(1)}+n^{2+1/6})\log(n)\log(n/\delta))$ time, where $\delta$ is the relative accuracy and $\alpha$ is the dual matrix exponent. The dual exponent $\alpha$ is the largest $a$ such that an $n\times n$ matrix can be multiplied with an $n\times n^{a}$ matrix in $n^{2+o(1)}$ arithmetic operations. For current $\omega\approx 2.38$ and $\alpha\approx 0.31$ this time complexity is just $O(n^{\omega}\log(n)\log(n/\delta))$. , where the bound $O(n^{\omega})$ is the number of arithmetic operations required to multiply two $n\times n$ matrices. The parameter $\omega$ is also called matrix exponent. For the generic case $d=\Omega(n)$, the latter complexity is essentially optimal as all known linear system solvers require up to $O(n^{\omega})$ time for solving $Ax=b$. As the complexity is essentially optimal, but the algorithm is randomized, a typical next step (e.g. [KT19, Cha00, PR02, MRSV17]) is to attempt to derandomize this algorithm. Derandomizing algorithms has the benefit that the required analysis can lead to further understanding of the studied problem. There have been precedences where derandomizing algorithms required developing new techniques, which then allowed for improvements in other settings. For example, in order to derandomize Karger’s edge connectivity algorithm [Kar00], Kawarabayashi and Thorup [KT19] had to develop new techniques, which then lead to new results in the distributed setting [DHNS19].

In the related area of linear program solvers in the real RAM model (i.e. when analyzing the complexity only in terms of the dimension, but not the bit-complexity of the input), a lot of effort has been put in derandomization and finding fast deterministic algorithms (see e.g. [CM96, BCM99, Cha16]). Yet, there is still a wide gap between the best randomized and deterministic complexity bounds. For an overview see [Cha16]. The fastest deterministic algorithm requires $O(nd^{d(1/2+o(1))})$ time [Cha16], while with randomized techniques an $O(nd^{2}+\exp(O(\sqrt{d\log d})))$ time bound is possible (a combination of [Cla95, Kal92, MSW96]). The same observation can be made in our setting, when analyzing the complexity with respect to the bit-complexity of the input, where the best deterministic bounds are $\widetilde{O}(\sqrt{n}\cdot nnz(A)+nd^{1.38})$ [Kar84] When using the $\widetilde{O}(\sqrt{n})$-iterations short step method., $\widetilde{O}(\sqrt{n}\cdot nnz(A)+n^{1.34}d^{1.15})$ [Vai89b] and $\widetilde{O}(d\cdot nnz(A)+d^{\omega+1})$ [Vai89a]. For curious readers we recommend [LS15]. They give a brief overview of these algorithms and offer a helpful graph that shows which algorithm is fastest for which range of $n,d,nnz(A)$. For $d=\Omega(n)$, all deterministic algorithms are stuck at $\Omega(n^{2.5})$ time. Further, these bounds are at least 30 years old and all new algorithms, that have been able to improve upon these bounds, crucially use randomized techniques. This raises the question: Is there a deterministic algorithm that can close the gap between deterministic and randomized complexity bounds, or at least break the 30 years old $\Omega(n^{2.5})$ barrier?

We are able to answer this question affirmatively by derandomizing the algorithm of Cohen et al. [CLS18]. Our deterministic algorithm is not just able to break the 30 years old barrier, it even matches one of the fastest randomized bounds of $\widetilde{O}(n^{\omega})$. This closes the complexity gap between randomized and deterministic algorithms for large $d$. More formally, we prove the following result:

Let $\min_{Ax=b,x\geq 0}c^{\top}x$ be a linear program without redundant constraints. Let $R$ be a bound on $\|x\|_{1}$ for all $x\geq 0$ with $Ax=b$. Then for any $0<\delta\leq 1$ we can compute $x\geq 0$ such that

in time $O(n^{\omega}\log^{2}(n)\log(n/\delta)),$ for the current matrix multiplication time with $\omega\approx 2.38$ [Wil12, Gal14].

which can be simplified to $O(n^{\omega}\log^{2}(n)\log(n/\delta))$ for current values of $\omega\approx 2.38$ [Wil12, Gal14], $\alpha\approx 0.31$ [GU18]. For integral $A,b,c$ the parameter $\delta=2^{-O(L)}$ is enough to round the approximate solution of Theorem 1.1 to an exact solution. Here $L=\log(1+\det_{\max}+\|c\|_{\infty}+\|b\|_{\infty})$ is the bit-complexity, where $\det_{\max}$ is the largest determinant of any square submatrix of $A$. [Ren88, LS13]

Derandomizing the $\widetilde{O}(n^{\omega})$ algorithm of [CLS18] was stated by Song as an open question in [Son19]. In addition to answering this open question, our techniques also allow us to simplify the analysis of the central path method used in [CLS18], reducing the length by roughly half.

Interior point algorithms must typically repeatedly compute the projection of a certain vector $v$, i.e. they must compute $Pv$ where $P$ is a projection matrix. It suffices to use an approximation $\widetilde{P}$ of $P$, and in each iteration the matrix $P$ changes only a bit, which allowed previous results to maintain the approximation $\widetilde{P}$ quickly (See for example [Kar84, NN89, Vai89b, LS15]). A natural barrier for improving linear program solvers is the fact that computing $\widetilde{P}v$ requires $\Omega(n^{2})$ for dense $v$. This leads to the $\Omega(n^{2.5})$ barrier for linear program solvers, because a total of $\Omega(\sqrt{n})$ projections must be computed. Cohen et al. were able to break this barrier in [CLS18] by sparsifying $v$ to some approximate $\widetilde{v}$ via random sampling and computing $\widetilde{P}\widetilde{v}$ instead of $\widetilde{P}v$. Our new approach for breaking this barrier deterministically is to maintain the product $\widetilde{P}\widetilde{v}$ directly for some approximation $\widetilde{v}$ of $v$. This is the key difference of our linear program solver compared to previous results, which only maintained $\widetilde{P}$.

We will now outline the difference between our deterministic $\widetilde{O}(n^{\omega})$ solver for linear programs and the randomized result of Cohen et al. [CLS18]. They managed to obtain a fast solver for linear programs by computing the projection $\widetilde{P}\widetilde{v}$ in subquadratic time using two clever tools:

They created a data-structure to maintain $\widetilde{P}$ in sub-quadratic time, amortized over $\sqrt{n}$ iterations.

They created a novel stochastic central path method which can sparsify the vector $v$ to some approximate $\widetilde{v}$ via random sampling. Thus the projection $\widetilde{P}\widetilde{v}$ could be computed in sub-quadratic worst-case time.

Derandomizing this algorithm seems like a difficult task as it is not clear how to obtain a deterministic sparsification of $v$. Recently [LSZ19] derandomized the central path method (2), so they could extend their linear program solver to the problem of Empirical Risk Minimization. However, in order to achieve $\widetilde{O}(n^{\omega})$ total time, they had to reduce some dimension in the representation of $\widetilde{P}$ via random sketching, which resulted in randomizing the data-structure (1).

In this paper we show how to completely derandomize the algorithm of [CLS18] via a data-structure that can maintain the projection $\widetilde{P}\widetilde{v}$ directly for some dense approximate $\widetilde{v}\approx v$, instead of just maintaining the matrix $\widetilde{P}$ as in [CLS18]. This result can be obtained in two different ways. One option is to use a dynamic linear system algorithm (e.g. [San04, vdBNS19]) via a black-box reduction, or alternatively one can interpret the resulting data-structure as a surprisingly simple extension of the data-structure used in [CLS18]. Indeed the algorithmic description of the data-structure (1) of [CLS18] grows only by a few lines (see Algorithm 1 in Section 4).

The high level idea of our new data-structure is that the vector $v$ can be written as some function $v_{i}=f(w_{i})$, where the argument vector $w$ does not change much between two iterations of the central path method. By approximating $w$ by some $\widetilde{w}$ we can re-use information of the previous iteration when computing the projection of $\widetilde{v}=f(\widetilde{w})$. One difficulty, that we must overcome, is that $\widetilde{v}:=f(\widetilde{w})$ is not an approximation of $v=f(w)$ in the classical sense (i.e. we can not satisfy $\|v-\widetilde{v}\|\leq\varepsilon\|v\|$ or even $\widetilde{v}_{i}\approx v_{i}$), even if $\widetilde{w}$ is an approximation of $w$, because for non-monotonous $f$, the vectors $v$ and $\widetilde{v}$ could point in opposite directions. This is a problem for the short step central path method, because these algorithm can be interpreted as some gradient descent and here $v$ depends on the gradient of some potential function. So if $\widetilde{v}$ points in the opposite direction, then the algorithm will actually increase the potential function instead of decreasing it.

Outline

In this section we outline how our algorithm works and how we adapt existing ideas such as the short step central path method and the projection maintenance. Readers only interested in verifying our algorithm can skip this overview, but reading it can help provide some intuition for how the different parts of our algorithm interact and what difficulties must be solved in order for our algorithm to work.

We start the outline with a brief summary of the short step central path method, which motivates why we must maintain a certain projection. Readers already familiar with the short step central path method can skip ahead to the next subsection 2.2.

In Section 2.2 we describe the task of the projection maintenance, and how we are able to perform this task quickly by using a certain notion of approximation (details in Section 4). The next Section 2.3 of the outline explains the difficulties that we encounter by using this type of approximation, and how we are able to solve these problems (details in Section 5).

We first give a brief summary of the short step central path method. Readers familiar with these types of algorithms can skip ahead to the next subsection 2.2.

Consider the linear program $\displaystyle\min_{Ax=b,x\geq 0}c^{\top}x$ and its dual program $\displaystyle\max_{A^{\top}y\leq c}b^{\top}y$. Given a feasible dual solution $y$ (a vector $y$ s.t. $A^{\top}y\leq c$), we can define the slack vector $s:=c-A^{\top}y$. Based on the complementary slackness condition (see e.g. [PS82]) we know a triple $(x,y,s)$ is optimal, if and only if

If only the last three conditions are satisfied, then we call the triple $(x,y,s)$ feasible. Given such a feasible triple, we define the vector $\mu$ such that $\mu_{i}:=x_{i}s_{i}$ and the complementary slackness theorem motivates why we should try to minimize the entries of $\mu$.

It is known how to transform the LP in such a way, that we can easily construct a feasible solution triple $(x,y,s)$ with $x_{i}s_{i}\approx 1$ for all $i=1,...,n$ (e.g. Lemma A.3 [YTM94]). Thus for $t:=1$ we have $\mu_{i}\approx t$. The idea is to repeatedly decrease $t$ and to modify the solution $x\leftarrow x+\delta_{x},y\leftarrow y+\delta_{y},s\leftarrow s+\delta_{s}$ in such a way, that the entries of $\mu$ stay close to $t$. The change of $\mu$ is given by $\mu^{\operatorname{new}}_{i}=(x+\delta_{x})_{i}(s+\delta_{s})_{i}=\mu_{i}+x_{i}\delta_{s,i}+s_{i}\delta_{x,i}+\delta_{x,i}\delta_{s,i}$ and if $\delta_{x},\delta_{s}$ are small enough, this can be approximated via $\mu^{\operatorname{new}}_{i}\approx\mu_{i}+x_{i}\delta_{s,i}+s_{i}\delta_{x,i}$. Thus to change $\mu$ by (approximately) $\delta_{\mu}$, we can solve the following linear system

where $X=\operatorname{diag}(x)$ and $S=\operatorname{diag}(s)$ are diagonal matrices with the entries of $x$ and $s$ on the diagonal respectively. The solution to this system is given by the following lemma:

The solution for $\delta_{x},\delta_{s}$ in (1) is given by

A typical choice for the decrement of $t$ is to multiply it by $1-O(\frac{1}{\sqrt{n}})$, which means it takes about $O(\sqrt{n}/\delta)$ iterations until $t$ reaches some desired accuracy parameter $\delta>0$ [Ren88, Vai87].

For the short step central path method the distance between $\mu$ and $t$ is typically measured by $\sum_{i=1}^{n}(\mu_{i}-t)^{2}=\|\mu-t\|_{2}^{2}$ and one tries to maintain $x,s$ in such a way that $\|\mu-t\|_{2}^{2}\leq O(t^{2})$. This can be modelled via the potential function $\Phi(x)=\|x\|_{2}^{2}$, and then one tries to maintain $\mu$ such that $\Phi(\mu/t-1)=O(1)$, which is equivalent to $\|\mu-t\|_{2}^{2}\leq O(t^{2})$. Thus a good choice for $\delta_{\mu}$ would be a vector with the same direction as $-\nabla\Phi(\mu/t-1)$, as this allows us to decrease the potential, which then means the distance between $\mu$ and $t$ is reduced.

2 Projection Maintenance (Details in Section 4)

In this subsection we outline one of the main results of this paper and sketch its proof. As described in the previous section, we must repeatedly compute $Pv$ for $P:=\sqrt{\frac{X}{S}}A^{\top}\left(A\frac{X}{S}A^{\top}\right)^{-1}A\sqrt{\frac{X}{S}}$ and $v:=\frac{\delta_{\mu}}{\sqrt{XS}}$, where the matrix $A$ describes the constraints of the linear program, $X$ and $S$ are diagonal matrices that depend on some current feasible solution and $\delta_{\mu}$ is some vector.

Our main result is to maintain an approximation of $Pv$ deterministically in $\widetilde{O}(n^{\omega-0.5}+n^{2.5-\alpha})$ amortized time, where $\omega$ is the current matrix multiplication exponent and $\alpha$ is the dual exponent. This new data-structure is a simple extension of the data-structure presented in [CLS18], which was able to maintain an approximation of $P$ within the same time bound, but their data-structure required up to $O(n^{2})$ time for computing $Pv$ for dense $v$.

The exact statement of our result involves various details, for example how $P$ and $v$ change over time. So we first want to describe the task of maintaining $Pv$ in more detail.

The matrix $P=\sqrt{\frac{X}{S}}A^{\top}\left(A\frac{X}{S}A^{\top}\right)^{-1}A\sqrt{\frac{X}{S}}$ shares a lot of structure between two iterations. Indeed only the diagonal matrices $X$ and $S$ change, while the matrix $A$ stays fixed. Thus for simplicity we define $U:=X/S$, in which case $P:=\sqrt{U}A^{\top}\left(AUA^{\top}\right)^{-1}A\sqrt{U}$ and only the diagonal matrix $U$ changes from one iteration to the next one.

For this task we would wish for a data-structure that can compute $Pv$ for any vector $v$ in $O(n^{\omega-0.5})$ time, which with $O(\sqrt{n})$ iterations would then result in an $O(n^{\omega})$-time solver for linear programs. More accurately, we hope for an algorithm that solves the following task:

Initialize$(A,u,v)$: Given matrix $A$ and two $n$ dimensional vectors $u,v$ we preprocess the matrix and return

where $U=diag(u)$ is the diagonal matrix with $u$ on the diagonal.

Update$(u,v)$: Given two $n$ dimensional vectors $u,v$, we must compute

It is not clear whether a data-structure exists for this task with $O(n^{\omega-0.5})$ update time, but due to the very first short step linear program solver by Karmarkar [Kar84] it is known that one can relax the requirements. Indeed it is enough to use an approximation of $P$.

Relaxation and result

Due to [Kar84] it is known, that it is enough to use an approximation $\widetilde{P}:=\sqrt{\widetilde{U}}A^{\top}(A\widetilde{U}A^{\top})^{-1}A\sqrt{\widetilde{U}}$ for $(1-\varepsilon)U\leq\widetilde{U}\leq(1+\varepsilon)U$, instead of the exact matrix $P$. We show later in Section 5, that it is also enough to approximate the vector $v$ via some $\widetilde{v}$. The type of approximation for $v$ is a bit different: We show in Section 5 that we can write $v=\delta_{\mu}/\sqrt{XS}$ as a function of $\mu/t$, so $v=f(\mu/t)$. We then “approximate” $v$ via some $\widetilde{v}:=f(\widetilde{\mu}/t)$, where $(1-\varepsilon)\mu\leq\widetilde{\mu}\leq(1+\varepsilon)\mu$. Note that thus $\widetilde{v}$ itself is not necessarily an approximation of $v$ in the classical sense (i.e. $\|v-\widetilde{v}\|_{2}\gg\varepsilon\|v\|_{2}$) and the two vectors might even point in opposite directions.

Motivated by these observations we want to maintain $\widetilde{P}\widetilde{v}$ instead of $Pv$. This idea allows for a speed-up, because in each iteration we only need to change the entries of $\widetilde{u}$ and $\widetilde{\mu}$ for which the $(1+\varepsilon)$-approximation condition is broken. Thus if the vectors $u$ and $\mu$ do not change much per iteration, then we only need to change few entries of $\widetilde{u}$ and $\widetilde{\mu}$. We prove in Section 5.2 that $u$ and $\mu$ satisfy the following condition:

Let $(u^{k})_{k\geq 1}$ be the sequence of vectors $u$, generated by the central path method. Then $\|(u^{k+1}-u^{k})/u^{k}\|_{2}\leq C$ for all $k$ and some constant $C$. (A similar statement can be made for $\mu$)

Thus, while we are not able to solve 2.2 exactly, we do obtain a data-structure that (i) maintains the solution approximately, and (ii) is fast if $\|(u^{k+1}-u^{k})/u^{k}\|_{2}$ and $\|(\mu^{k+1}-\mu^{k})/\mu^{k}\|_{2}$ are small.

Initialize$(A,u,f,v,\varepsilon_{mp})$: The data-structure preprocesses the given two $n$ dimensional vectors $u,v$, the $d\times n$ matrix $A$ the function $f$ in $O(n^{2}d^{\omega-2})$ time. The given parameter $\varepsilon_{mp}>0$ specifies the accuracy of the approximation.

Update$(u,v)$: Given two $n$ dimensional vectors $u,v$. Then the data-structure returns four vectors

Here $\widetilde{U}$ is the diagonal matrix $\operatorname{diag}(\widetilde{u})$ and $\widetilde{v}$ a vector such that

If the update sequence $u^{(1)},...,u^{(T)}$ (and likewise $v^{(1)},...,v^{(T)}$) satisfies

for all $k=1,...,T$ then the total time for the first $T$ updates is

There are two equivalent ways to prove Theorem 2.4: One could use the data-structures of [San04, vdBNS19] which maintain $M^{-1}b$ for some non-singular matrix $M$ and some vector $b$. Via a black-box reduction these data-structures would then be able to maintain $\widetilde{P}\widetilde{v}$ and applying the tools of [CLS18] for optimizing the amortized complexity would then result in Theorem 2.4.

If one tries to write down a pseudo-code description of the resulting data-structure, then the code is very similar to the data-structure from [CLS18]. This is because all these data-structures are based on the Sherman-Morrison-Woodburry identity. Hence an alternative way to prove Theorem 2.4 is to take the data-structure from [CLS18], which already maintains $\widetilde{P}$, and extend such that it also maintains $\widetilde{P}\widetilde{v}$.

In this paper we present the second option, where we modify the existing data-structure of [CLS18]. This is because we want to highlight that our derandomization result can be obtained from a simple modification of the existing randomized algorithm. Though for the curious reader we also give a sketch of the first variant in Appendix B.

Proof sketch (Details in Section 4)

We now outline how to obtain Theorem 2.4 by extending the data-structure of [CLS18] to also maintain $\widetilde{P}\widetilde{v}$, instead of just $\widetilde{P}$. Their data-structure internally has three matrices $M,L,R$ with the property

where $M$ is some $n\times n$ matrix and $L,R$ are rectangular matrices with some $m\ll n$ columns. With each update, the matrices $L,R$ change and the number of their columns may increase. This way the $n^{2}$ entries of the matrix $\widetilde{P}$ are not explicitly computed and a sub-quadratic update time can be achieved.

As the number of columns $m$ of $L,R$ grows, the data-structure will become slower and slower. Once these matrices have too many columns, the data-structure performs a “reset”. This means we set

and the matrices $L,R$ are set to be empty (so zero columns). Thus after the reset we have $\widetilde{P}=M+LR^{\top}=M$, so (3) is still satisfied. Such a reset requires $\Omega(n^{2})$ time, but it does not happen too often so the cost is small on average. Section 5 of [CLS18] is about bounding this amortized cost.

One can now easily maintain $\widetilde{P}f(\widetilde{v})$ as follows: Assume we already know $Mf(\widetilde{v})$, then a new solution $\widetilde{P}f(\widetilde{v})$ is given by

because of (3). Here the term $LR^{\top}f(\widetilde{v})$ can be computed in $O(nm)\ll O(n^{2})$ time, because $L,R$ have $m\ll n$ columns. The assumption, that $Mf(\widetilde{v})$ is known, can be satisfied easily: During the initialization of the algorithm we compute this value, and whenever $M$ changes (i.e. during the reset (4)) we can compute the new $Mf(\widetilde{v})$ in $O(n^{2})$ time. This does not affect the complexity of the data-structure, because a reset does already require $\Omega(n^{2})$ time to compute the new $M$.

At last, we must handle the case where entries of $\widetilde{v}$ are changed. Let’s say $\widetilde{v}^{\operatorname{new}}\leftarrow\widetilde{v}+\delta_{v}$, then

where the last equality comes from (3). The complexity can be bounded as follows: The term $\widetilde{P}f(\widetilde{v})$ is computed as described in (5). The second term $M(f(\widetilde{v}^{\operatorname{new}})-f(\widetilde{v}))$ can be computed quickly because on average $\widetilde{v}^{\operatorname{new}}$ and $\widetilde{v}$ differ in only few entries, because of the small change to $v$ per iteration (as given by (2) of Theorem 2.4). The last term $LR^{\top}(f(\widetilde{v}^{\operatorname{new}})-f(\widetilde{v}))$ is again computed quickly because the matrices $L,R$ have very few columns.

3 Adapting the Central Path Method for Approximate Projection Maintenance (Details in Section 5)

We now outline difficulties that occur, if one tries to use the projection maintenance algorithm (Theorem 2.4, outlined in Section 2.2) in the classical central path method (outlined in Section 2.1), and how we are able to solve these issues in Section 5.

The central path method can be interpreted as some gradient descent, where we try to minimize some potential. When we use the data-structure of Theorem 2.4, then we are essentially performing this gradient descent while using some approximate gradient. This approximation is of such low quality, that the approximate gradient occasionally points in a completely wrong direction, effectively increasing the potential instead of decreasing it. By adapting the potential function, we are able to prove that the approximate gradient only points in the wrong direction when the potential is small. Whenever the potential is large, the approximate gradient points in the correct direction. (A formal proof of this will be in Section 5.3, Lemma 5.14.) This adaption to the short step central path method allows us to handle these faulty approximate gradients. Before we can outline why this is true, we must first explain why we obtain these faulty approximate gradients in the first place.

The central path method tries to maintain some vector $\mu$ close to a scalar $t$, where the relative distance is measured via some potential function $\Phi(\mu/t-1)$. The central path method tries to minimize this potential function by solving some linear system that depends on the gradient $\nabla\Phi(\mu/t-1)$.

Adapting the central path method

The majority of the proof that this choice for $\Phi$ works, is adapted from [CLS18]. For their stochastic central path method, Cohen et al. sparsify the gradient $\nabla\Phi(\mu/t-1)$ via randomly sampling its entries. This sparsification could be interpreted as some type of approximation of the gradient, which allows us to adapt their proof to our new notion of “approximating” the gradient via $\nabla\Phi(\widetilde{\mu}/t-1)$ for $\widetilde{\mu}\approx\mu$. The main difference is that in [CLS18], the exact and approximate gradient always point in the same direction (i.e. their inner product is positive), so in [CLS18] it was a bit easier to show that the potential $\Phi(\mu/t-1)$ decreases in each iteration. For comparison, when using our approximation, the inner product of exact gradient $\nabla\Phi(\mu/t-1)$ and the approximation $\nabla\Phi(\widetilde{\mu}/t-1)$ may become negative. So we must spend some extra effort in Section 5.3 to show that the approximate gradient points in the correct direction, whenever $\Phi(\mu/t-1)$ is large (this will be proven in Lemma 5.14). Intuitively, this is true because when $\Phi(\mu/t-1)$ is large, then there are many indices $i$ such that $\mu_{i}$ is further from $t$ than some $(1\pm\varepsilon_{mp})$-factor. As outlined before, this means the approximate gradient tries to change the $i$th coordinate of $\mu$ in the correct direction, i.e. $i$th entry of the exact and approximate gradient have the same sign.

We also want to point out, that our approach of using a gradient w.r.t the approximate $\widetilde{\mu}\approx\mu$ means it is enough to maintain an approximate $\widetilde{x}\approx x$, $\widetilde{s}\approx s$, so $\widetilde{x}\widetilde{s}=:\widetilde{\mu}\approx\mu$. The same observation was made independently in [LSZ19], where that property was exploited to compute the steps $\delta_{x}$, $\delta_{s}$ approximately via random sketching. The analysis of their central path method is based on modifying the standard newton steps to be a variant of gradient descent in some hessian norm space. In comparison our proof is arguably simpler, as we perform a typical gradient descent w.r.t $\Phi(\widetilde{\mu}/t)$.

Preliminaries

For the linear program $\min_{Ax=b,x\geq 0}c^{\top}x$ we assume there are no redundant constraints, i.e. the matrix $A$ is of rank $d$ and $n\geq d$.

For two $n$ dimensional vectors $v,w$ their inner products is written as $v^{\top}w$ or alternatively $\langle v,w\rangle$. We write $vw$ for the entry-wise product, so $(vw)_{i}:=v_{i}w_{i}$. The same is true for all other arithmetic operations, for example $(v/w)_{i}:=v_{i}/w_{i}$ and $(\sqrt{v})_{i}:=\sqrt{v_{i}}$. For a scalar $s$ the product $sv$ is the typical entry-wise product and analogously we define $v-s$ as the entry-wise difference, so $(v-s)_{i}:=v_{i}-s$.

Inequalities

We write $v\leq w$ if $v_{i}\leq w_{i}$ for all $i=1,...,n$ and we use the notation $v\approx_{\varepsilon}w$ to express a $(1\pm\varepsilon)$ approximation, defined as $(1-\varepsilon)w\leq v\leq(1+\varepsilon)w$. Note that $v\approx_{\varepsilon}w$ is not symmetric, but it implies $w\approx_{2\varepsilon}v$ for $\varepsilon\leq 1/2$.

Relative error and multiplicative change

Let $v,w,\delta_{v},\delta_{w}$ be vectors, such that $v^{\operatorname{new}}=v+\delta_{v}$, $w^{\operatorname{new}}=w+\delta_{w}$ then

Here the last term can be bounded via $\|\frac{\delta_{v}}{v}\frac{\delta_{w}}{w}\|_{2}\leq\|\frac{\delta_{v}}{v}\|_{\infty}\|\frac{\delta_{w}}{w}\|_{2}\leq\|\frac{\delta_{v}}{v}\|_{2}\|\frac{\delta_{w}}{w}\|_{2}$ ∎

Further, if $v^{\operatorname{new}}$ has small multiplicative change compared to $v$, then the same is true for $1/v^{\operatorname{new}}$ and $1/v$.

Fast Matrix Multiplication

We write $O(n^{\omega})$ for the arithmetic complexity of multiplying two $n\times n$ matrices. Computing the inverse has the same complexity. The exponent $\omega$ is also called the matrix exponent. We call $\alpha$ the dual matrix exponent, which is the largest value such that multiplying a $n\times n$ matrix with an $n\times n^{\alpha}$ requires $O(n^{2+o(1)})$ time. The current best bounds are $\omega\approx 2.38$ [Wil12, Gal14] and $\alpha\approx 0.31$ [GU18].

Projection Maintenance

Initialize$(A,u,f,v,\varepsilon_{mp})$: The data-structure preprocesses the given two $n$ dimensional vectors $u,v$, the $d\times n$ matrix $A$ the function $f$ in $O(n^{2}d^{\omega-2})$ time. The given parameter $\varepsilon_{mp}>0$ specifies the accuracy of the approximation.

Update$(u,v)$: Given two $n$ dimensional vectors $u,v$. Then the data-structure returns four vectors

Here $\widetilde{U}$ is the diagonal matrix $\operatorname{diag}(\widetilde{u})$ and $\widetilde{v}$ a vector such that

If the update sequence $u^{(1)},...,u^{(T)}$ (and likewise $v^{(1)},...,v^{(T)}$) satisfies

for all $k=1,...,T$ then the total time for the first $T$ updates is

This section is split into three parts: We first present the algorithm and give a high-level description in Section 4.1. The next subsection (Section 4.2) proves that the algorithm returns the correct result, and at last in Section 4.3 we bound the complexity of the algorithm.

Algorithm 1 describes a data-structure, so we have variables that persist between calls to its function Update. What these variables represent might be a bit hard to deduce from just reading the pseudo-code, so we want to give a brief outline of Algorithm 1 here. This outline is not required for verifying the proofs, but it might help for understanding how the algorithm works.

The internal variables are $n$-dimensional vectors $\widetilde{u},\widetilde{v},w$ and an $n\times n$ matrix $M$. The relationship between them is

where $\widetilde{U}=\operatorname{diag}(\widetilde{u})$.

These internal variables are useful because of the following reason: In each call to Update, the data-structure receives two new vectors $u^{\operatorname{new}},v^{\operatorname{new}}$ and for $U^{\operatorname{new}}=\operatorname{diag}(u^{\operatorname{new}})$ the task is to return an approximation of $\sqrt{U^{\operatorname{new}}}A^{\top}(AU^{\operatorname{new}}A^{\top})^{-1}A\sqrt{U^{\operatorname{new}}}f(v^{\operatorname{new}})$ by $(1\pm\varepsilon_{mp})$-approximating $U^{\operatorname{new}}$ and $v^{\operatorname{new}}$. Thus if

then $\sqrt{\widetilde{U}}w$ would be the desired approximate result. If this $(1+\varepsilon_{mp})$-approximation condition (8) is not satisfied, then we can define two new valid approximations

for all $i=1,...,n$. If $\widetilde{u}^{\operatorname{new}}$ and $\widetilde{u}$ (and respectively $\widetilde{v}^{\operatorname{new}}$, $\widetilde{v}^{\operatorname{new}}$) differ in at most $k$ many entries, then it is known (via Sherman-Morrison-Woodbury identity Lemma 4.3) that one can quickly construct two $n\times k$ matrices $R,L$ such that

This in turn means that we can get the desired approximate result as follows:

Here each term can be computed in at most $O(nk)$ time, because the first term is the already known vector $w$, the vector of the second term is sparse, and $R,L$ are $n\times k$ matrices.

Thus, if $k$ is small, then we can maintain the solution quickly. In [CLS18], Cohen et al. have developed a strategy with low amortized cost, that specifies when to recompute $M$ for some new $\widetilde{u}$, such that $k$ stays small. In their algorithm they do not maintain the matrix-vector product, so their data-structure does not have the internal variables $w$ and $\widetilde{v}$. We extend their strategy to also recompute $w$ for some new $\widetilde{v}$, such that the above outlined procedure has low amortized cost.

2 Correctness

The task of this subsection is to prove the following lemma, which says that the vectors returned by Algorithm 1 are as specified in Lemma 4.1.

After every update to Algorithm 1 with input $(u^{\operatorname{new}},v^{\operatorname{new}})$ the returned vectors $\widetilde{u}^{\operatorname{new}},\widetilde{v}^{\operatorname{new}},f(\widetilde{v}^{\operatorname{new}}),r$ satisfy $u^{\operatorname{new}}\approx_{\varepsilon_{mp}}\widetilde{u}^{\operatorname{new}}$, $v^{\operatorname{new}}\approx_{\varepsilon_{mp}}\widetilde{v}^{\operatorname{new}}$ and

Before we can prove this lemma, we must first prove that the internal variables of the data-structure save the correct values, i.e. we want to prove that equation (7) is correct. For this we must first state the following lemma from [CLS18], based on Sherman-Morison-Woodbury identity.

If $M=A^{\top}(A\widetilde{U}A^{\top})^{-1}A$ at the start of the update of Algorithm 1 and $M_{S},M_{S,S},\Delta_{S,S}$ are chosen as described in Algorithm 1, then we have

We can now prove that the internal variables store the correct values.

At the start of every update to Algorithm 1 we have

If it is the first update after the initialization, then the claim is true by definition of the procedure Initialize. Next, we prove that at the end of every call to Update we satisfy (9), if (9) was satisfied at the start of Update. This then implies Lemma 4.4. If $k\geq n^{a}$, then 33 makes sure that $M=A^{\top}(A\widetilde{U}^{\operatorname{new}}A^{\top})^{-1}A$ (see Lemma 4.3). The next lines set $w\leftarrow M\sqrt{\widetilde{U}^{\operatorname{new}}}f(v^{\operatorname{new}})$, $\widetilde{v}\leftarrow v^{\operatorname{new}}$ and $\widetilde{u}\leftarrow\widetilde{u}^{\operatorname{new}}$. Thus (9) is satisfied for the case $k\geq n^{a}$. If $|T|\geq n^{a}$, then we compute $w\leftarrow M\sqrt{\widetilde{U}}f(v^{\operatorname{new}})$ and set $\widetilde{v}\leftarrow v^{\operatorname{new}}$. The matrices $M$ and $\widetilde{U}$ are not modified, so (9) is satisfied. If $|T|<n^{a}$, then we do not change $M,\widetilde{u},\widetilde{v}$ or $r$, so (9) is satisfied.

We now prove the correctness of Algorithm 1 by proving Lemma 4.2.

Note that we always have $u^{\operatorname{new}}\approx_{\varepsilon_{mp}}\widetilde{u}^{\operatorname{new}}$ by 27.

In 33 we have set $M=A^{\top}(A\widetilde{U}^{\operatorname{new}}A^{\top})^{-1}A$ (see Lemmas 4.3 and 4.4). Hence by setting $r\leftarrow\sqrt{\widetilde{U}^{\operatorname{new}}}w=\sqrt{\widetilde{U}^{\operatorname{new}}}M\sqrt{\widetilde{U}^{\operatorname{new}}}f(v^{\operatorname{new}})$, and $\widetilde{v}^{\operatorname{new}}\leftarrow v^{\operatorname{new}}$, all claims of Lemma 4.2 are satisfied.

Case |T|≥na|T|\geq n^{a}:

In this case we set $r$ to the following value:

Here the equality comes from Lemmas 4.3 and 4.4. Further, we set $\widetilde{v}^{\operatorname{new}}\leftarrow v^{\operatorname{new}}$, so Lemma 4.2 is correct for the case $|T|\geq n^{a}$.

Case |T|<na|T|<n^{a}:

Here $v^{\operatorname{new}}\approx_{\varepsilon_{mp}}\widetilde{v}^{\operatorname{new}}$ by 43, so we are left with verifying $r$. First note that $w=M\sqrt{\widetilde{U}}f(\widetilde{v})$ by Lemma 4.4, so $w+M(\sqrt{\widetilde{U}^{\operatorname{new}}}f(\widetilde{v}^{\operatorname{new}})-\sqrt{\widetilde{U}}f(\widetilde{v}))=M\sqrt{\widetilde{U}^{\operatorname{new}}}f(\widetilde{v}^{\operatorname{new}})$. Thus $r$ is set to the following term:

Where the last equality comes from Lemmas 4.3 and 4.4.

3 Complexity

In this section we will bound the complexity of Algorithm 1, proving the stated complexity bound in Lemma 4.1:

If the updates to Algorithm 1 satisfy the condition (6) as stated in Lemma 4.1, then after $T$ updates the total update time of Algorithm 1 is

The preprocessing requires $O(n^{2}d^{\omega-2})$ time.

As our data-structure is a modification of the data-structure presented in [CLS18], we must only analyze the complexity of the modified part. To bound the complexity of the unmodified sections of our algorithm, we will here refer to [CLS18]. The complexity analysis in [CLS18] requires an entire section (about 7 pages) via analysis of some complicated potential function. In the Appendix (Lemma A.2) we present an alternative simpler proof.

The preprocessing requires $O(n^{2}d^{\omega-2})$ time. After $T$ updates the total time of all updates of Algorithm 1, when ignoring the branch for $k<n^{a}$ (so we assume that branch of 36 has cost $0$), is

When ignoring the branch of 36, then our algorithm performs the same operations as [CLS18][Algorithm 3] and we both maintain $M$ in the exact same way. The only difference is that we also compute the vector $r$ in 34, but this requires only $O(n^{2})$ time and is subsumed by the complexity of 33. Thus our time complexity (when ignoring the branch of 36) can be bounded by the update complexity of [CLS18][Algorithm 3], which is the complexity stated in Lemma 4.6. In the same fashion we can bound the complexity of the preprocessing. The preprocessing of [CLS18][Algorithm 3] takes $O(n^{2}d^{\omega-2})$ time, where their algorithm computes only the matrix $M$. The only difference in our algorithm is that we also compute the vector $w$ in 13. The required $O(n^{2})$ time to compute $w$ is subsumed by computing $M$. ∎

In order to prove Lemma 4.5 we only need to bound the complexity of the branch for the case $k<n^{a}$. The time required by all other steps of Algorithm 1 is already bounded by Lemma 4.6.

In every update we must compute $(\Delta_{S,S}^{-1}+M_{S,S})^{-1}$, which takes $O(n^{a\cdot\omega})$ time via the assumption $k<n^{a}$. Additionally, if $|T|<n^{a}$, then one update requires additional $O(n^{1+a})$ operations to compute $r$ and $w$, because $(f(\widetilde{v}^{\operatorname{new}})-f(\widetilde{v}))$ and $(\sqrt{\widetilde{U}}-\sqrt{\widetilde{U}^{\operatorname{new}}})$ both have at most $n^{a}$ non-zero entries and $M_{S}$ is a $n\times n^{a}$ matrix.

If $T\geq n^{a}$, then computing $r$ and $w$ can take up to $O(n^{2})$ operations. This can happen at most every $O(n^{a/2}\varepsilon_{mp}/C)$ updates by Lemma A.1, because $\sum_{i=1}^{n}((v_{i}^{\operatorname{new}}-v_{i})/v_{i})^{2}\leq C^{2}$, Hence the amortized time per update is $O(n^{2-a/2}C/\varepsilon_{mp})$.

Note that by assuming $a\leq\alpha\leq 1$ the term $O(n^{a\cdot\omega})$ is subsumed by $O(n^{1+a})$, because $\omega\leq 3-\alpha$, so $a\cdot\omega\leq a(3-\alpha)\leq a(3-a)\leq 1+a$. ∎

Central Path Method

In this section we prove the main result Theorem 1.1, by showing how to use the projection maintenance algorithm of Section 4 to obtain a fast deterministic algorithm for solving linear programs.

The algorithm for Theorem 1.1 is based on the short step central path method, outlined in Section 2.1: We construct some feasible solution triple $(x,y,s)$ with $xs=:\mu\approx 1$ and then repeatedly decrease $t$ while maintaining $x,s$ such that $\mu$ stays close to $t$. Once $t$ is small enough, we have a good approximate solution. This is a high-level summary of Algorithm 2, which first constructs a solution, and then runs a while-loop until $t$ is small enough. The actual hard part, maintaining the solution pair $x,s$ with $\mu\approx t$, is done in Algorithm 3. For this task, Algorithm 3 solves a linear system (similar to (1) in Section 2.1) via the data-structure of Lemma 4.1. The majority of this section is dedicated to proving that Algorithm 3 does not require too much time and does indeed maintain the solution pairs $(x,s)$ with $\mu\approx t$. For this we must verify the following three properties:

Algorithm 3 does solve an approximate variant of the linear system (1).

We do not change the linear system too much between two calls to Algorithm 3. Otherwise the data-structure of Lemma 4.1 would become too slow.

The approximate result obtained in Algorithm 3 is good enough to maintain $x,s$ such that $\mu$ is close to $t$.

The proof for this is based on the stochastic central path method by Cohen et al. [CLS18]. In [CLS18], they randomly sampled a certain vector, while in our algorithm this vector will be approximated deterministically via the data-structure of Algorithm 1. This derandomization has the nice side-effect, that we can skip many steps of Cohen et al.’s proof. For example they had to bound the variance of random vectors, which is no longer necessary for our algorithm.

The outline of this section is as follows. We first explain in more detail how Algorithm 3 works in Section 5.1, where we also verify the first requirement, that Algorithm 3 does indeed solve the system (1) approximately. In the next Section 5.2, we check that the input parameters for the data-structure of Lemma 4.1, used by Algorithm 3, do not change too much per iterations. The last Section 5.3 verifies, that we indeed always have $\mu\approx t$. We also consolidate all results in the last subsection by proving the main result Theorem 1.1.

In this section we outline how Algorithm 3 works and we prove that it does indeed solve the linear system (1) (outlined in Section 2.1) in some approximate way. The high-level idea of Algorithm 3 is as follows: In order to maintain $\mu$ close to $t$, we want to measure the distance via some potential function $\Phi(\mu/t-1)$. As we want to minimize the distance, it makes sense to change $\mu$ by some $\delta_{\mu}$, which points in the same direction as $-\nabla\Phi(\mu/t-1)$. We can find out how to change $x$ and $s$, in order to change $\mu$ by approximately $\delta_{\mu}$, by solving the linear system (1) via the data-structure of Lemma 4.1.

In reality, we choose $\delta_{\mu}$ to be slightly different:

where $t^{\operatorname{new}}:=(1-\frac{\varepsilon}{3\sqrt{n}})$ is the new smaller value that we want to set $t$ to, and $\varepsilon$ is a parameter for how large our step size should be for decreasing $t$.

This choice of $\delta_{\mu}$ is motivated by the fact, that the first term $(\frac{t^{{\operatorname{new}}}}{t}-1)\mu$ leads to some helpful cancellations in later proofs. The second term is the one pointing in the direction of $-\nabla\Phi(\mu/t-1)$, which is motivated by decreasing $\Phi(\mu/t-1)$.

One can split $\delta_{\mu}$ into the two terms $\delta_{t}=(\frac{t^{{\operatorname{new}}}}{t}-1)\mu$ and $\delta_{\Phi}=\frac{\varepsilon}{2}\cdot t^{{\operatorname{new}}}\cdot\frac{\nabla\Phi(\mu/t-1)}{\|\nabla\Phi(\mu/t-1)\|_{2}}$.

Algorithm 3 approximates both vectors in a different way. Specifically, given $x,s,\mu,\delta_{t},\delta_{\Phi},\delta_{\mu}$, Algorithm 3 internally maintains approximations $\widetilde{x},\widetilde{s},\widetilde{\mu},\widetilde{\delta}_{t},\widetilde{\delta}_{\Phi},\widetilde{\delta}_{\mu}$ with the following properties (here $\varepsilon_{mp}>0$ is the accuracy parameter for Lemma 4.1)

and for these approximate values, we solve the following system (which is the same as (1), but using the approximate values):

We prove in two steps that Algorithm 3 does indeed solve (11) for approximate values as in (10): First we prove in Lemma 5.2 that the approximations are as stated in (10), then we show in Lemma 5.3 that we indeed solve the linear system (11).

Note that $\widetilde{\delta}_{\Phi}$ is not an approximation of $\delta_{\Phi}$ in the classical sense (likewise $\widetilde{\delta}_{\mu}$ and $\delta_{\mu}$) and the vectors could point in completely different directions. They are only “approximate” in the sense that their definition is the same, but for $\widetilde{\delta}_{\Phi}$ we replace $\mu$ by the approximate $\widetilde{\mu}$.

As outlined in the overview Section 2.3, this results in our algorithm not always decreasing the difference between $\mu$ and $t$. We prove in Section 5.3 that this is not a problem, if we use the following potential function $\Phi$, accuracy parameter $\varepsilon_{mp}$ (for Lemma 4.1) and step size $\varepsilon$.

where $\cosh(x):=(e^{x}+e^{-x})/2$, $\lambda=40\ln n$. For the step size $\varepsilon$ and the accuracy parameter $\varepsilon_{mp}$ for Lemma 4.1, assume $0<\varepsilon_{mp}\leq\varepsilon\leq 1/(1500\ln n)$.

The computed vectors $\widetilde{x},\widetilde{s},\widetilde{\mu},\widetilde{\delta}_{t},\widetilde{\delta}_{\Phi},\widetilde{\delta}_{\mu}$ in Algorithm 3 satisfy the following properties: Let $\widetilde{\mu}/t$ be the approximation of $\mu/t$ maintained internally by $\text{mp}_{\Phi}$, then $\mu\approx_{\varepsilon_{mp}}\widetilde{\mu}$ and $\widetilde{\delta}_{\Phi}=-\frac{\varepsilon}{2}\cdot t^{\operatorname{new}}\cdot\frac{\nabla\Phi_{\lambda}(\widetilde{\mu}/t-1)}{\|\nabla\Phi_{\lambda}(\widetilde{\mu}/t-1)\|}$. Further $\delta_{t}\approx_{\varepsilon_{mp}}\widetilde{\delta}_{t}$, $x\approx_{\varepsilon_{mp}}\widetilde{x}$, $s\approx_{2\varepsilon_{mp}}\widetilde{s}$.

The returned vector $m$ in line 16 is an approximation in the sense that $\mu/t\approx_{\varepsilon_{mp}}m$, which means $\mu\approx_{\varepsilon_{mp}}mt=:\widetilde{\mu}$. We have $\frac{x}{s}=:u\approx_{\varepsilon_{mp}}\widetilde{u}$, hence $\frac{u}{\widetilde{u}}\approx_{\varepsilon_{mp}}1$ and $1\approx_{\varepsilon_{mp}}\frac{\widetilde{u}}{u}$. Thus $x\approx_{\varepsilon_{mp}}x\sqrt{\frac{\widetilde{\mu}}{\mu}\frac{\widetilde{u}}{u}}=\widetilde{x}$ and $s\approx_{2\varepsilon_{mp}}s\sqrt{\frac{\widetilde{\mu}}{\mu}\frac{u}{\widetilde{u}}}=\widetilde{s}$.

As potential function we have chosen $\Phi_{\lambda}(x)=\sum_{i=1}^{n}\cosh(x_{i})$, so $(\nabla\Phi_{\lambda}(x-1)/\sqrt{x})_{i}=\lambda\sinh(\lambda(x_{i}-1))/\sqrt{x_{i}}.$ This means $\lambda\sinh(\lambda(x-1))/\sqrt{x}$ for $x=\mu/t$ is $\nabla\Phi_{\lambda}(\mu/t-1)/\sqrt{\mu/t}$ and $w=\lambda\sinh(\lambda(m-1))/\sqrt{m}=\nabla\Phi_{\lambda}(\widetilde{\mu}/t-1)/\sqrt{\widetilde{\mu}/t}.$ Hence we have that $\widetilde{\delta}_{\Phi}=-\frac{\varepsilon}{2}\cdot t^{\operatorname{new}}\cdot\frac{\sqrt{\widetilde{\mu}/t}w}{\|\nabla\Phi_{\lambda}(\widetilde{\mu}/t-1)\|_{2}}=-\frac{\varepsilon}{2}\cdot t^{\operatorname{new}}\cdot\frac{\nabla\Phi_{\lambda}(\widetilde{\mu}/t-1)}{\|\nabla\Phi_{\lambda}(\widetilde{\mu}/t-1)\|_{2}}$. We also have $\delta_{t}=(\frac{t^{{\operatorname{new}}}}{t}-1)\mu$, $\mu\approx_{\varepsilon_{mp}}v^{2}$ and $\mu\approx_{\varepsilon_{mp}}\widetilde{\mu}$, so $\mu\approx_{\varepsilon_{mp}}v\sqrt{\widetilde{\mu}}$ which implies $\delta_{t}\approx_{\varepsilon_{mp}}(\frac{t^{{\operatorname{new}}}}{t}-1)v\sqrt{\widetilde{\mu}}=:\widetilde{\delta}_{t}$.

The computed vectors in Algorithm 3 satisfy the following linear system:

We define the following projection matrix:

Hence the change to $x$ and $s$ is given by

2 Bounding the change per iteration

Algorithm 3 uses the data-structure of Lemma 4.1. The complexity of this data-structure depends on how much the input parameters (in our case $u:=x/s$, $\mu$ and $\mu/t$) change per iteration. In this section we prove:

Assume $\mu\approx_{0.1}t$. Let $\mu^{\operatorname{new}}:=(x+\widetilde{\delta}_{x})(s+\widetilde{\delta}_{s})$, the value of $\mu$ in the upcoming iteration, and let $u:=\frac{x}{s}$, $u^{\operatorname{new}}:=\frac{x+\widetilde{\delta}_{x}}{s+\widetilde{\delta}_{s}}$, then

In order to prove this lemma, we must assume that $\mu$ is currently a good approximation of $t$. We assume the following proposition, which is proven in the next subsection.

For the input to Algorithm 3 we have $\mu\approx_{0.1}t$

How much we change $x,s$ depends on how long the vector $\delta_{\mu}$ is, so we start by bounding that length.

$\|\delta_{t}\|_{2}\leq 1.1\frac{\varepsilon}{3}t$, $\|\delta_{\Phi}\|_{2}\leq\frac{\varepsilon}{2}t$, $\|\delta_{\mu}\|_{2}\leq\varepsilon t$ $\|\widetilde{\delta}_{t}\|_{2}\leq 1.2\frac{\varepsilon}{3}t$, $\|\widetilde{\delta}_{\Phi}\|_{2}\leq\frac{\varepsilon}{2}t$, $\|\widetilde{\delta}_{\mu}\|_{2}\leq\varepsilon t$