Faster Algorithms for Structured Linear and Kernel Support Vector Machines

## Abstract

Quadratic programming is a ubiquitous prototype in convex programming. Many combinatorial optimizations on graphs and machine learning problems can be formulated as quadratic programming; for example, Support Vector Machines (SVMs). Linear and kernel SVMs have been among the most popular models in machine learning over the past three decades, prior to the deep learning era. Generally, a quadratic program has an input size of Θ ​ ( n 2 ) Θ superscript 𝑛 2 \Theta(n^{2}) , where n 𝑛 n is the number of variables. Assuming the Strong Exponential Time Hypothesis ( SETH ), it is known that no O ​ ( n 2 − o ​ ( 1 ) ) 𝑂 superscript 𝑛 2 𝑜 1 O(n^{2-o(1)}) algorithm exists (Backurs, Indyk, and Schmidt, NIPS’17). However, problems such as SVMs usually feature much smaller input sizes: one is given n 𝑛 n data points, each of dimension d 𝑑 d , with d ≪ n much-less-than 𝑑 𝑛 d\ll n . Furthermore, SVMs are variants with only O ​ ( 1 ) 𝑂 1 O(1) linear constraints. This suggests that faster algorithms are feasible, provided the program exhibits certain underlying structures. In this work, we design the first nearly-linear time algorithm for solving quadratic programs whenever the quadratic objective has small treewidth or admits a low-rank factorization, and the number of linear constraints is small. Consequently, we obtain a variety of results for SVMs: • For linear SVM, where the quadratic constraint matrix has treewidth τ 𝜏 \tau , we can solve the corresponding program in time O ~ ​ ( n ​ τ ( ω + 1 ) / 2 ​ log ⁡ ( 1 / ϵ ) ) ~ 𝑂 𝑛 superscript 𝜏 𝜔 1 2 1 italic-ϵ \widetilde{O}(n\tau^{(\omega+1)/2}\log(1/\epsilon)) ; • For linear SVM, where the quadratic constraint matrix admits a low-rank factorization of rank- k 𝑘 k , we can solve the corresponding program in time O ~ ​ ( n ​ k ( ω + 1 ) / 2 ​ log ⁡ ( 1 / ϵ ) ) ~ 𝑂 𝑛 superscript 𝑘 𝜔 1 2 1 italic-ϵ \widetilde{O}(nk^{(\omega+1)/2}\log(1/\epsilon)) ; • For Gaussian kernel SVM, where the data dimension d = Θ ​ ( log ⁡ n ) 𝑑 Θ 𝑛 d=\Theta(\log n) and the squared dataset radius is sub-logarithmic in n 𝑛 n , we can solve it in time O ​ ( n 1 + o ​ ( 1 ) ​ log ⁡ ( 1 / ϵ ) ) 𝑂 superscript 𝑛 1 𝑜 1 1 italic-ϵ O(n^{1+o(1)}\log(1/\epsilon)) . We also prove that when the squared dataset radius is at least Ω ​ ( log 2 ⁡ n ) Ω superscript 2 𝑛 \Omega(\log^{2}n) , then Ω ​ ( n 2 − o ​ ( 1 ) ) Ω superscript 𝑛 2 𝑜 1 \Omega(n^{2-o(1)}) time is required. This improves upon the prior best lower bound in both the dimension d 𝑑 d and the squared dataset radius.

## 1 Introduction

Quadratic programming (QP) represents a class of convex optimization problems that optimize a quadratic objective over the intersection of an affine subspace and the non-negative orthant(^1^11There are classes of QPs with quadratic constraints as well. However, in this paper, we focus on cases where the constraints are linear.). QPs naturally extend linear programming by incorporating a quadratic objective, and they find extensive applications in operational research, theoretical computer science, and machine learning [, , , , , ]. The quadratic objective introduces challenges: QPs with a general (not necessarily positive semidefinite) symmetric quadratic objective matrix are NP-hard to solve [, ]. When the quadratic objective matrix is positive semidefinite, the problem becomes weakly polynomial-time solvable, as it can be reduced to convex empirical risk minimization [] (refer to Section [6](#S6) for further discussion).

Formally, the QP problem is defined as follows:

###### Definition 1.1 (Quadratic Programming) .

Given an $n\times n$ symmetric, positive semidefinite objective matrix $Q$, a vector $c\in\mathbb{R}^{n}$, and a polytope described by a pair $(A\in\mathbb{R}^{m\times n},b\in\mathbb{R}^{m})$, the linearly constrained quadratic programming (LCQP) or quadratic programming (QP) problem seeks to solve the following optimization problem:

$$ $\displaystyle\min_{x\in\mathbb{R}^{n}}~{}$ $\displaystyle~{}\frac{1}{2}x^{\top}Qx+c^{\top}x$ (1) $\displaystyle\mathrm{s.t.~{}}$ $\displaystyle~{}Ax=b$ $\displaystyle~{}x\geq 0.$ $$

A classic application of QP is the support vector machine (SVM) problem [, ]. In SVMs, a dataset ${x_{1},\ldots,x_{n}}\subseteq\mathbb{R}^{d}$ is provided, along with corresponding labels ${y_{1},\ldots,y_{n}}\in{\pm 1}$. The objective is to identify a hyperplane that separates the two groups of points with opposite labels, while maintaining a large margin from both. Remarkably, this popular machine learning problem can be formulated as a QP and subsequently solved using specialized QP solvers []. Thus, advancements in QP algorithms could potentially lead to runtime improvements for SVMs.

Despite its practical and theoretical significance, algorithmic quadratic programming has garnered relatively less attention compared to its close relatives in convex programming, such as linear programming [, , , ], convex empirical risk minimization [, ], and semidefinite programming [, , ]. In this work, we aim to take a pioneering step in developing a fast and robust interior point-type algorithm for solving QPs to high precision. We particularly focus on improving the runtime for high-precision hard- and soft-margin SVMs. For the purposes of this discussion, we will concentrate on hard-margin SVMs, with the understanding that our results naturally extend to soft-margin variants. We begin by introducing hard-margin linear SVMs:

###### Definition 1.2 (Linear SVM) .

Given a dataset $X\in\mathbb{R}^{n\times d}$ and a collection of labels $y_{1},\ldots,y_{n}$ each in ${\pm 1}$, the *linear SVM problem* requires solving the following quadratic program:

$$ $\displaystyle\max_{\alpha\in\mathbb{R}^{n}}$ $\displaystyle~{}{\bf 1}_{n}^{\top}\alpha-\frac{1}{2}\alpha^{\top}(yy^{\top}\circ XX^{\top})\alpha,$ (2) $\displaystyle{\rm s.t.~{}}$ $\displaystyle~{}\alpha^{\top}y=0,$ $\displaystyle~{}\alpha\geq 0.$ $$

where $\circ$ denotes the Hadamard product.

It should be noted that this formulation is actually the *dual* of the SVM optimization problem. The primal program seeks a vector $w\in\mathbb{R}^{d}$ such that

$$ $\displaystyle\min_{w\in\mathbb{R}^{d}}$ $\displaystyle~{}\frac{1}{2}\|w\|_{2}^{2},$ s.t. $\displaystyle~{}y_{i}(w^{\top}x_{i}-b)\geq 1,\qquad\forall i\in[n],$ $$

where $b\in\mathbb{R}$ is the bias term. Given the solution $\alpha\in\mathbb{R}^{n}$, one can conveniently convert it to a primal solution: $w^{*}=\sum_{i=1}^{n}\alpha_{i}^{*}y_{i}x_{i}$. At first glance, one might be inclined to solve the primal problem directly, especially in cases where $d\ll n$, as it presents a lower-dimensional optimization challenge compared to the dual. However, the dual formulation becomes particularly advantageous when solving the kernel SVM, which maps features to a high or potentially infinite-dimensional space.

###### Definition 1.3 (Kernel SVM) .

Given a dataset $X\in\mathbb{R}^{n\times d}$ and a positive definite kernel function ${\sf K}:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$, let $K\in\mathbb{R}^{n\times n}$ denote the kernel matrix, where $K_{i,j}={\sf K}(x_{i},x_{j})$. With a collection of labels $y_{1},\ldots,y_{n}$ each in $\{\pm 1\}$, the *kernel SVM problem* requires solving the following quadratic program:

$$ $\displaystyle\max_{\alpha\in\mathbb{R}^{n}}~{}$ $\displaystyle~{}{\bf 1}_{n}^{\top}\alpha-\frac{1}{2}\alpha^{\top}(yy^{\top}\circ K)\alpha,$ (3) $\displaystyle{\rm s.t.~{}}$ $\displaystyle~{}\alpha^{\top}y=0,$ $\displaystyle~{}\alpha\geq 0.$ $$

The positive definite kernel function ${\sf K}$ corresponds to a feature mapping, implying that ${\sf K}(x_{i},x_{j})=\phi(x_{i})^{\top}\phi(x_{j})$ for some $\phi:\mathbb{R}^{d}\rightarrow\mathbb{R}^{s}$. Thus, solving the primal SVM can be viewed as solving the optimization problem on the transformed dataset. However, the primal program’s dimension depends on the (transformed) data’s dimension $s$, which can be infinite. Conversely, the dual program, with dimension $n$, is typically easier to solve. Throughout this paper, when discussing the SVM program, we implicitly refer to the dual quadratic program, not the primal.

One key aspect of the SVM program is its minimal equality constraints. Specifically, for both linear and kernel SVMs, there is only a single equality constraint of the form $\alpha^{\top}y=0$. This constraint arises naturally from the bias term in the primal SVM formulation and its Lagrangian. The limited number of constraints enables us to design QP solvers with favorable dependence on the number of data points $n$, albeit with a higher dependence on the number of constraints $m$, thus offering effective end-to-end guarantees for SVMs.

Previous efforts to solve the SVM program efficiently typically involve breaking down the large QP into smaller, constant-sized QPs. These algorithms, while easy to implement and well-suited to modern hardware architectures, are widely adopted in popular SVM libraries, such as LibSVM []. Theoretically, [] systematically analyzed this class of algorithms, demonstrating that to achieve an $\epsilon$-approximation solution, $\widetilde{O}(\epsilon^{-2}B\cdot\operatorname{nnz}(A))$ time is sufficient, where $B$ is the squared-radius of the dataset and $\operatorname{nnz}(\cdot)$ denotes the number of nonzero entries. Unfortunately, these algorithms are inherently first-order and require $\operatorname{poly}(\epsilon^{-1})$ iterations to converge, rendering them impractical for applications requiring high precision.

To develop a high-precision algorithm with $\operatorname{poly}\log(\epsilon^{-1})$ dependence instead of $\operatorname{poly}(\epsilon^{-1})$, we focus on second-order methods for QPs. A variety of approaches have been explored in previous works, including the interior point method [], active set methods [], augmented Lagrangian techniques [], conjugate gradient, gradient projection, and extensions of the simplex algorithm [, , ]. Our interest is particularly piqued by the interior point method (IPM). Recent advances in the robust IPM framework have led to significant successes for convex programming problems [, , , , , , , , , , ]. These successes are a result of combining robust analysis of IPM with dedicated data structure design.

The application of IPM to solve QPs with a constant number of constraints is not entirely novel; existing work [] has already demonstrated the application of IPM in solving the linear SVM problem. However, the runtime efficiency of their algorithm is sub-optimal. Each iteration requires multiplying a $d\times n$ matrix with an $n\times d$ matrix, taking $O(nd^{\omega-1})$ time, where $\omega\approx 2.37$ is the matrix multiplication exponent [, , ]. Moreover, the IPM necessitates $O(\sqrt{n}\log(1/\epsilon))$ iterations for convergence. This results in an overall runtime of $O(n^{1.5}d^{\omega-1}\log(1/\epsilon))$, which is super-linear in the dataset size, even when the data dimension $d$ is small. In practical scenarios, where $n$ is typically large, the $n^{1.5}$ runtime dependence becomes prohibitive. Therefore, it is crucial to develop an algorithm with almost- or nearly-linear dependence on $n$ and (poly-)logarithmic dependence on $1/\epsilon$. This will be our primary objective.

For linear SVM, we propose two nearly-linear time algorithms with high-precision guarantees, applicable when the dataset exhibits low-rank or low-treewidth properties. Our first focus is on the low-rank scenario.

###### Theorem 1.4 (Low-rank QP and Linear SVM) .

Given a quadratic program as defined in Definition [1.1](#S1.Thmtheorem1), and assuming a low-rank factorization of the quadratic objective matrix $Q=UV^{\top}$, where $U,V\in\mathbb{R}^{n\times k}$, there exists an algorithm that can solve the program ([1](#S1.E1)) up to $\epsilon$ error in $\widetilde{O}(n(k+m)^{(\omega+1)/2}\log(1/\epsilon))$ time.

Specifically, for linear SVM (as per Definition [1.2](#S1.Thmtheorem2)) with $d\ll n$, one can solve program ([2](#S1.E2)) up to $\epsilon$ error in $\widetilde{O}(nd^{(\omega+1)/2}\log(1/\epsilon))$ time.

Considering linear SVM where the dataset $X\in\mathbb{R}^{n\times d}$ is a tall, skinny matrix with $n\gg d$, this naturally leads to a low-rank factorization of $Q$ as $Q=(yy^{\top})\circ(XX^{\top})$. Another significant structure is *sparsity*, reflected in the *treewidth* of $XX^{\top}$.

###### Theorem 1.5 (Low-treewidth QP and Linear SVM) .

For a quadratic program as outlined in Definition [1.1](#S1.Thmtheorem1), suppose we have a tree decomposition(^2^22If the treewidth $\tau$ is known but an effective tree decomposition is not given, [] can be used to obtain a suitable decomposition in $O(m^{1+o(1)})$ time. Further details are in Section [4.2](#S4.SS2).) of the adjacency graph $Q$ with size $\tau$. An algorithm then exists that can solve program ([1](#S1.E1)) up to $\epsilon$ error in $\widetilde{O}(n(\tau+m)^{(\omega+1)/2}\log(1/\epsilon))$ time.

Particularly, for linear SVM (as in Definition [1.2](#S1.Thmtheorem2)) with a provided tree decomposition of size $\tau$, the program ([2](#S1.E2)) can be solved up to $\epsilon$ error in $\widetilde{O}(n\tau^{(\omega+1)/2}\log(1/\epsilon))$ time.

It is noteworthy that the runtime of our solver aligns with the state-of-the-art small treewidth linear program solver []. However, Theorem [1.5](#S1.Thmtheorem5) is not directly comparable to [] as we impose the treewidth assumption on the quadratic objective matrix $Q$, while the latter applies it to the adjacency graph induced by the constraint matrix $A$. In the context of linear SVM, where only a constant number of equality constraints exist, applying the assumption to the constraints is less meaningful. Conversely, the quadratic objective matrix $Q$ effectively captures the correlation between data and labels. Hence, it is more reasonable to assume $Q$ has low-treewidth.

While both of the above results imply nearly-linear time algorithms for structured linear SVMs, the situation changes significantly when considering kernel SVMs. Firstly, while the matrix $XX^{\top}$ might be sparse and have low-treewidth, the kernel matrix $K$ is inherently dense, comprising $\Theta(n^{2})$ nonzero entries, thereby imposing an $\Omega(n)$ lower bound on the treewidth. Additionally, without structural assumptions, the kernel matrix can be full-rank, rendering the low-rank QP solver inapplicable. In fact, it has been demonstrated that for data dimensions $d=\omega(\log n)$, no algorithm can approximate kernel SVM within error $\exp(-\omega(\log^{2}n))$ in time $O(n^{2-o(1)})$ [].

Conversely, recent efforts aim to expedite computation with respect to the kernel matrix faster than quadratic [, , ], especially when the kernel possesses certain smooth and Lipschitz properties. For instance, when kernel functions are sufficiently smooth, efficient approximation using low-degree polynomials is feasible, leading to an approximate low-rank factorization. A prime example is the Gaussian RBF kernel, where [] showed that for dimension $d=\Theta(\log n)$ and squared dataset radius $B=o(\log n)$, there exist low-rank matrices $U,V\in\mathbb{R}^{n\times n^{o(1)}}$ such that for any vector $x\in\mathbb{R}^{n}$, $\|(K-UV^{\top})x\|_{\infty}\leq\epsilon\|x\|_{1}$. Utilizing this factorization, they developed an algorithm to solve the Batch Gaussian KDE problem in $O(n^{1+o(1)})$ time.

Based on this dichotomy in fast kernel matrix algebra, we establish two results: 1) Solving Gaussian kernel SVM in $O(n^{1+o(1)}\log(1/\epsilon))$ time is feasible when $B=o(\frac{\log n}{\log\log n})$, and 2) Assuming SETH, no sub-quadratic time algorithm exists for $B=\Omega(\log^{2}n)$ in SVMs without bias and $B=\Omega(\log^{6}n)$ in SVMs with bias. This improves the lower bound established by [] in terms of dimension $d$.

###### Theorem 1.6 (Gaussian Kernel SVM) .

Given a dataset $X\in\mathbb{R}^{n\times d}$ with $d=\Theta(\log n)$ and a squared radius denoted by $B$. Let ${\sf K}(x_{i},x_{j})=\exp(-\|x_{i}-x_{j}\|_{2}^{2})$ be the Gaussian kernel function. Then, for the kernel SVM problem defined in Definition [1.3](#S1.Thmtheorem3),

- •
If $B=o(\frac{\log n}{\log\log n})$, there exists an algorithm that solves the Gaussian kernel SVM up to $\epsilon$ error in time $O(n^{1+o(1)}\log(1/\epsilon))$;
- •
If $B=\Omega(\log^{2}n)$, then assuming SETH, any algorithm that solves the Gaussian kernel SVM without a bias term up to $\exp(-\omega(\log^{2}n))$ error would require $\Omega(n^{2-o(1)})$ time;
- •
If $B=\Omega(\log^{6}n)$, then assuming SETH, any algorithm that solves the Gaussian kernel SVM with a bias term up to $\exp(-\omega(\log^{2}n))$ error would require $\Omega(n^{2-o(1)})$ time.

To our knowledge, this is the first almost-linear time algorithm for Gaussian kernel SVM even when $d=\log n$ and the radius is small. Our algorithm effectively utilizes the rank-$n^{o(1)}$ factorization of the Gaussian kernel matrix alongside our low-rank QP solver. The design template of our algorithm is powerful, as it leverages almost-linear time algorithms for computing spectral sparsifiers of the kernel matrix with only $O(n\log n)$ nonzero entries. We leave developing algorithms for sparse kernel matrices as a future direction.

###### Definition 1.1 (Quadratic Programming) .

###### Definition 1.2 (Linear SVM) .

###### Definition 1.3 (Kernel SVM) .

###### Theorem 1.4 (Low-rank QP and Linear SVM) .

###### Theorem 1.5 (Low-treewidth QP and Linear SVM) .

###### Theorem 1.6 (Gaussian Kernel SVM) .

### 1.1 Related Work

##### Support Vector Machines

SVM, one of the most prominent machine learning models before the rise of deep learning, has a rich literature dedicated to its algorithmic speedup. For linear SVM, [] offers a first-order algorithm that solves its QP in nearly linear time, but with a runtime dependence of $\epsilon^{-2}$, limiting its use in high precision settings. For SVM classification, existing algorithms such as SVM-Light [], SMO [], LIBSVM [], and SVM-Torch [] perform well in high-dimensional data settings. However, their runtime scales super-linearly with $n$, making them less viable for large datasets. Previous investigations into solving linear SVM QP via interior point methods [] have been somewhat basic, leading to an overall runtime of $O(n^{1.5}d^{\omega-1}\log(1/\epsilon))$. For a more comprehensive survey on efficient algorithms for SVM, refer to []. On the hardness side, [] provides an efficient reduction from the Hamming nearest neighbor problem to Gaussian kernel SVM, establishing an almost-quadratic lower bound assuming SETH.

##### Interior Point Method

The interior point method, a well-established approach for solving convex programs under constraints, was first proposed by [] as a (weakly) polynomial-time algorithm for linear programs, later improved by [] in terms of runtime. Recent work by Cohen, Lee, and Song [] has shown how to solve linear programs with interior point methods in the current matrix multiplication time, utilizing a robust IPM framework. Subsequent studies [, , , , , , ] have further refined their algorithm or applied it to different optimization problems.

##### Kernel Matrix Algebra

Kernel methods, fundamental in machine learning, enable feature mappings to potentially infinite dimensions for $n$ data points in $d$ dimensions. The kernel matrix, a crucial component of kernel methods, often has a prohibitive quadratic size for explicit formation. Recent active research focuses on computing and approximating kernel matrices and related tasks in sub-quadratic time, such as kernel matrix-vector multiplication, spectral sparsification, and Laplacian system solving. The study by [] introduces a comprehensive toolkit for solving these problems in almost-linear time for small dimensions, leveraging techniques like polynomial methods and ultra Johnson-Lindenstrauss transforms. Alternatively, [, ] reduce various kernel matrix algebra tasks to kernel density estimation (KDE), which recent advancements in KDE data structures [, , ] have made more efficient. A recent contribution by [] provides a tighter characterization of the low-degree polynomial approximation for the $e^{-x}$ function, leading to more efficient algorithms for the Batch Gaussian KDE problem.

## 2 Technique Overview

In this section, we provide an overview of the techniques employed in our development of two nearly-linear time algorithms for structured QPs. In Section [2.1](#S2.SS1), we detail the robust IPM framework, which forms the foundation of our algorithms. Subsequent sections, namely Section [2.2](#S2.SS2) and Section [2.3](#S2.SS3), delve into dedicated data structures designed for efficiently solving low-treewidth and low-rank QPs, respectively. Finally, in Section [2.4](#S2.SS4), we discuss the adaptation of these advanced QP solvers for both linear and kernel SVMs.

### 2.1 General Strategy

Our algorithm is built upon the robust IPM framework, an efficient variant of the primal-dual central path method []. This framework maintains a primal-dual solution pair $(x,s)\in\mathbb{R}^{n}\times\mathbb{R}^{n}$. To understand the central path for QPs, we first consider the central path equations for linear programming (see [, ] for reference):

$$ $\displaystyle s/t+\nabla\phi(x)=$ $\displaystyle~{}\mu,$ $\displaystyle Ax=$ $\displaystyle~{}b,$ $\displaystyle A^{\top}y+s=$ $\displaystyle~{}c,$ $$

where $x$ is the primal variable, $s$ is the slack variable, $y$ is the dual variable, $\phi(x)$ is a self-concordant barrier function, and $\mu$ denotes the error. The central path is defined by the trajectory of $(x,s)$ as $t$ approaches 0 0.

In quadratic programming, we modify these equations:

$$ $\displaystyle s/t+\nabla\phi(x)=$ $\displaystyle~{}\mu,$ $\displaystyle Ax=$ $\displaystyle~{}b,$ $\displaystyle-Qx+A^{\top}y+s=$ $\displaystyle~{}c,$ $$

where $Q$ is the positive semidefinite objective matrix. The key difference in the central path equations for LP and QP is the inclusion of the $-Qx$ term in the third equation, significantly affecting algorithm design.

The Newton step yields the following update rules (detailed derivation in Section [9.2](#S9.SS2)):

$$ $\displaystyle\delta_{x}$ $\displaystyle=tB^{-1/2}(I-P)B^{-1/2}\delta_{\mu},$ $\displaystyle\delta_{y}$ $\displaystyle=-t(AB^{-1}A^{\top})^{-1}AB^{-1}\delta_{\mu},$ $\displaystyle\delta_{s}$ $\displaystyle=t\delta_{\mu}-t^{2}HB^{-1/2}(I-P)B^{-1/2}\delta_{\mu},$ $\displaystyle\text{where}\qquad H$ $\displaystyle=\nabla^{2}\phi(x),\qquad B=Q+tH,$ $\displaystyle P$ $\displaystyle=B^{-1/2}A^{\top}(AB^{-1}A^{\top})^{-1}AB^{-1/2},$ $$

where $\delta_{x}$, $\delta_{y}$, $\delta_{s}$, and $\delta_{\mu}$ are the incremental steps for $x$, $y$, $s$, and $\mu$, respectively.

The robust IPM approximates these updates rather than applying them exactly. It maintains an approximate primal-dual solution pair $(\overline{x},\overline{s})\in\mathbb{R}^{n}\times\mathbb{R}^{n}$ and computes the steps using this approximation. Provided the approximation is sufficiently accurate, it can be shown (see Section [9](#S9) for more details) that the algorithm converges efficiently to the optimal solution along the robust central path.

Therefore, the critical challenge lies in efficiently maintaining $(\overline{x},\overline{s})$, an approximation to $(x,s)$, when $(x,s)$ evolves following the robust central path steps. The primary difficulty is that explicitly managing the primal-dual solution pair $(x,s)$ is inefficient due to potential dense changes. Such changes can lead to dense updates in $H$, slowing down the computation of steps. The innovative aspect of robust IPM is recognizing that $(x,s)$ are only required at the algorithm’s conclusion, not during its execution. Instead, we can identify entries with significant changes and update the approximation $(\overline{x},\overline{s})$ correspondingly. With IPM’s lazy updates, only a nearly-linear number of entries are adjusted throughout the algorithm:

$$ $\displaystyle\sum_{t=1}^{T}\|\overline{x}^{(t)}-\overline{x}^{(t-1)}\|_{0}+\|\overline{s}^{(t)}-\overline{s}^{(t-1)}\|_{0}=$ $\displaystyle~{}\widetilde{O}(n\log(1/\epsilon))$ $$

where $T=\widetilde{O}(\sqrt{n}\log(1/\epsilon))$ is the number of iterations for IPM convergence. This indicates that, on average, each entry of $\overline{x}$ and $\overline{s}$ is updated $\log(1/\epsilon)$ times, facilitating rapid updates to these quantities and, consequently, to $H$.

In the special case where $Q=0$, the path reverts to the LP case, with $B=tH$ being a diagonal matrix, allowing for efficient computation and updates of $B^{-1}$. This simplifies maintaining $AB^{-1}A^{\top}$, as updates to $B^{-1}$ correspond to row and column scaling of $A$. However, in the QP scenario, where $B$ is symmetric positive semidefinite, maintaining the term $AB^{-1}A^{\top}$ becomes more complex. Nevertheless, when the number of constraints is small, as in SVMs, this issue is less problematic. Yet, even with this simplification, the challenge is far from trivial, given the presence of terms like $B^{-1/2}$ in the robust central path steps. While the matrix Woodbury identity could be considered, it falls short when maintaining a square root term. Despite these hurdles, we construct efficient data structures for $B^{-1/2}$ maintenance when $Q$ possesses succinct representations, such as low-rank and low-treewidth. In subsequent discussions on structured QPs, we first explore the low-treewidth settings, which also provide valuable insights for the low-rank scenario.

### 2.2 Low-Treewidth Setting: How to Leverage Sparsity

In scenarios with low treewidth, we assume that the adjacency graph of $Q\in\mathbb{R}^{n\times n}$, denoted as $G_{Q}$, has a treewidth of $\tau$. The adjacency graph of $Q$ is constructed by placing $n$ vertices $v_{1},\ldots,v_{n}$ and forming an edge $\{v_{i},v_{j}\}$ whenever the entry $Q_{i,j}$ is non-zero. A *tree decomposition* of $G_{Q}$ arranges its vertices into bags, which collectively form a tree structure. For any two bags $X_{i}$ and $X_{j}$, if a vertex $v$ is present in both, it must also be included in all bags along the path between $X_{i}$ and $X_{j}$. Additionally, each pair of adjacent vertices in the graph must be present together in at least one bag. The treewidth $\tau$ is defined as the maximum size of a bag minus one. Intuitively, a graph $G_{Q}$ with a small treewidth $\tau$ implies a structure akin to a tree. For a formal definition, refer to Definition [4.1](#S4.Thmtheorem1). When relating this combinatorial structure back to linear algebra, a low-treewidth graph corresponds to a (generalized) adjacency matrix whose sparsity pattern allows for a column-sparse Cholesky factorization. Symbolically, we represent this as $B=LL^{\top}$ to denote the sparse Cholesky factors(^3^33Note that adding a non-negative diagonal matrix to $Q$ does not change its sparsity pattern, hence $B$ also retains the treewidth $\tau$.).

Under any coordinate update to $\overline{x}$, $B$ is updated on only one diagonal entry, enabling efficient updates to $L$. The remaining task is to use this Cholesky decomposition to maintain the central path step.

By expanding the central path equations and substituting $B=LL^{\top}$, we derive

$$ $\displaystyle\delta_{x}$ $\displaystyle=tB^{-1/2}(I-P)B^{-1/2}\delta_{\mu}$ $\displaystyle=tB^{-1}\delta_{\mu}-tB^{-1}A^{\top}(AB^{-1}A^{\top})^{-1}AB^{-1}\delta_{\mu}$ $\displaystyle=tL^{-\top}L^{-1}\delta_{\mu}-tL^{-\top}L^{-1}A^{\top}(AL^{-\top}L^{-1}A^{\top})^{-1}AL^{-\top}L^{-1}\delta_{\mu},$ $\displaystyle\delta_{s}$ $\displaystyle=t\delta_{\mu}-t^{2}HB^{-1/2}(I-P)B^{-1/2}\delta_{\mu}$ $\displaystyle=t\delta_{\mu}-t^{2}L^{-\top}L^{-1}\delta_{\mu}+t^{2}L^{-\top}L^{-1}A^{\top}(AL^{-\top}L^{-1}A^{\top})^{-1}AL^{-\top}L^{-1}\delta_{\mu}.$ $$

Updates to the diagonal of $B$ do not change $L$’s nonzero pattern, allowing for efficient utilization of the sparse factor and maintenance of $L^{-1}A^{\top}\in\mathbb{R}^{n\times m}$ and $L^{-1}\delta_{\mu}\in\mathbb{R}^{n}$. Terms like $(AL^{-\top}L^{-1}A^{\top})^{-1}AL^{-\top}L^{-1}\delta_{\mu}\in\mathbb{R}^{m}$ can also be explicitly maintained.

With this approach, we propose the following implicit representation for maintaining $(x,s)$:

$$ $\displaystyle x$ $\displaystyle=\widehat{x}+H^{-1/2}\mathcal{W}^{\top}(h\beta_{x}-\widetilde{h}\widetilde{\beta}_{x}+\epsilon_{x}),$ (4) $\displaystyle s$ $\displaystyle=\widehat{s}+H^{1/2}c_{s}\beta_{c_{s}}-H^{1/2}\mathcal{W}^{\top}(h\beta_{s}-\widetilde{h}\widetilde{\beta}_{s}+\epsilon_{s}),$ (5) $$

where $\widehat{x},\widehat{s}\in\mathbb{R}^{n}$, $\mathcal{W}=L^{-1}H^{1/2}\in\mathbb{R}^{n\times n}$, $h=L^{-1}\overline{\delta}_{\mu}\in\mathbb{R}^{n}$, $c_{s}=H^{-1/2}\overline{\delta}_{\mu}\in\mathbb{R}^{n}$, $\beta_{x},\beta_{s},\beta_{c_{s}}\in\mathbb{R}$, $\widetilde{h}=L^{-1}A^{\top}\in\mathbb{R}^{n\times m}$, $\widetilde{\beta}_{x},\widetilde{\beta}_{s}\in\mathbb{R}^{m}$, $\epsilon_{x},\epsilon_{s}\in\mathbb{R}^{n}$. All quantities except for $\mathcal{W}$ can be explicitly maintained. For linear programming, the implicit representation is as follows:

$$ $\displaystyle x$ $\displaystyle=\widehat{x}+H^{-1/2}\beta_{x}c_{x}-H^{-1/2}\mathcal{W}^{\top}(\beta_{x}h+\epsilon_{x})$ $\displaystyle s$ $\displaystyle=\widehat{s}+H^{1/2}\mathcal{W}^{\top}(\beta_{s}h+\epsilon_{s}),$ $$

with $\mathcal{W}=L^{-1}AH^{-1/2}$ maintained implicitly and the other terms explicitly.

The representation in ([4](#S2.E4)) and ([5](#S2.E5)) enables us to maintain the central path step using a combination of “coefficients” $h+\widetilde{h}\widetilde{\beta}_{x}$ and “basis” $\mathcal{W}^{\top}$. We need to detect entries of $\overline{x}$ that deviate significantly from $x$ and capture these changes with $\|H^{1/2}(\overline{x}-x)\|_{2}$. We maintain this vector using $x_{0}+\mathcal{W}^{\top}(h+\widetilde{h}\widetilde{\beta}_{x})$. Here, $\mathcal{W}^{\top}$ acts as a wavelet basis and the vector $h+\widetilde{h}\widetilde{\beta}_{x}$ as its multiscale coefficients. While computing and maintaining $\mathcal{W}^{\top}h$ seems challenging, leveraging column-sparsity of $L^{-1}$ is possible through contraction with a vector $v$:

$$ $\displaystyle v^{\top}\mathcal{W}^{\top}=$ $\displaystyle~{}(\mathcal{W}v)^{\top}$ $\displaystyle=$ $\displaystyle~{}(L^{-1}H^{1/2}v)^{\top}.$ $$

By applying the Johnson-Lindenstrauss transform (JL) in place of $v$, we can quickly approximate $\|\mathcal{W}^{\top}h\|_{2}$ by maintaining $\Phi\mathcal{W}^{\top}$ for a JL matrix $\Phi$. Similarly, we handle $\mathcal{W}^{\top}\widetilde{h}\widetilde{\beta}_{x}$ by explicitly computing $A^{\top}\widetilde{\beta}_{x}$ and using the sparsity of $L^{-1}$ for $\widetilde{h}\widetilde{\beta}_{x}$.

We focus on entries significantly diverging from $x_{0}$, the heavy entries of $\mathcal{W}^{\top}(h+\widetilde{h}\widetilde{\beta}_{x})$. Here, the treewidth-$\tau$ decomposition enables quick computation of an elimination tree based on $L^{-1}$’s sparsity, facilitating efficient estimation of $\|(\mathcal{W}^{\top}(h+\widetilde{h}\widetilde{\beta}_{x}))_{\chi(v)}\|_{2}$ for any subtree $\chi(v)$(^4^44Given any tree node $v$, we use $\chi(v)$ to denote the subtree rooted at $v$.). With an elimination tree of height $\widetilde{O}(\tau)$, we can employ heavy-light decomposition [] for an $O(\log n)$-height tree.

Using these data structures, convergence is established using the robust IPM framework [, ]. While the framework is generally applicable to QPs, initializing the optimization remains a challenge, as prior literature lacks a method for QPs. We propose a simpler objective $x_{0}=\arg\min_{x\in\mathbb{R}^{n}}\sum_{i=1}^{n}\phi_{i}(x_{i})$ with $\phi_{i}$ as the log-barrier function, resembling the initialization in []. This initial point enables us to solve an augmented quadratic program that increases dimension by 1.

### 2.3 Low-Rank Setting: How to Utilize Small Factorization

The low-treewidth structure can be considered a form of sparsity, allowing for a sparse factorization $B=LL^{\top}$. Another significant structure arises when the matrix $Q$ admits a low-rank factorization. For instance, let $Q=UV^{\top}$ where $U,V\in\mathbb{R}^{n\times k}$ and $k\ll n$, then $B=Q+tH=UV^{\top}+tH$. Although $Q$ has a low-rank structure, $B$ may not be low-rank due to the diagonal matrix being dense. However, in the central path equations, we need only handle $B^{-1}$, which can be efficiently maintained using the matrix Woodbury identity:

$$ $\displaystyle B^{-1}=$ $\displaystyle~{}t^{-1}H^{-1}-t^{-2}H^{-1}U(I+t^{-1}V^{\top}H^{-1}U)^{-1}V^{\top}H^{-1},$ $$

Given that $H$ is diagonal, the complex term $(I+t^{-1}V^{\top}H^{-1}U)^{-1}$ can be quickly updated under sparse changes to $H^{-1}$ by simply scaling rows of $U$ and $V$. With only a nearly-linear number of updates to $H^{-1}$, the total update time across $\widetilde{O}(\sqrt{n}\log(1/\epsilon))$ iterations is bounded by $\widetilde{O}(nk^{\omega-1}+k^{\omega})$. We modify the $(x,s)$ implicit representation as follows:

$$ $\displaystyle x$ $\displaystyle=\widehat{x}+H^{-1/2}h\beta_{x}+H^{-1/2}\widehat{h}\widehat{\beta}_{x}+H^{-1/2}\widetilde{h}\widetilde{\beta}_{x},$ (6) $\displaystyle s$ $\displaystyle=\widehat{s}+H^{1/2}h\beta_{s}+H^{1/2}\widehat{h}\widehat{\beta}_{s}+H^{1/2}\widetilde{h}\widetilde{\beta}_{s},$ (7) $$

where $\overline{x},\overline{s}\in\mathbb{R}^{n}$, $h=H^{-1/2}\overline{\delta}_{\mu}\in\mathbb{R}^{n}$, $\widehat{h}=H^{-1/2}U\in\mathbb{R}^{n\times k}$, and $\widetilde{h}=H^{-1/2}A^{\top}\in\mathbb{R}^{n\times m}$, with $\widetilde{\beta}_{x},\widetilde{\beta}_{s}\in\mathbb{R}^{m}$. The nontrivial terms to maintain are $\widehat{h}$ and $\widetilde{h}$, but both can be managed straightforwardly: updates to $H^{-1/2}$ correspond to scaling rows of $U$ and $A^{\top}$, and can be performed in total $\widetilde{O}(nk)$ and $\widetilde{O}(nm)$ time, respectively.

The remaining task is to design a data structure for detecting heavy entries. Instead of starting with an elimination tree and re-balancing it through heavy-light decomposition, we construct a balanced tree on $n$ nodes, hierarchically dividing length-$n$ vectors by their indices. Sampling is then performed by traversing down to the tree’s leaves. While a heavy-hitter data structure could lead to improvements in poly-logarithmic and sub-logarithmic factors, we primarily focus on polynomial dependencies on various parameters and leave this enhancement for future exploration.

### 2.4 Gaussian Kernel SVM: Algorithm and Hardness

Our specialized QP solvers provide fast implementations for linear SVMs where the data dimension $d$ is much smaller than $n$ (low-rank QP) or where the matrix $XX^{\top}$ has small treewidth (low-treewidth QP). However, for kernel SVM, forming the kernel matrix exactly would take $\Theta(n^{2})$ time. Fortunately, advancements in kernel matrix algebra [, , , ] have enabled sub-quadratic algorithms when the data dimension $d$ is small or the kernel matrix has a relatively large minimum entry. Both [] and [] introduce algorithms for spectral sparsification, generating an approximate matrix $\widetilde{K}\in\mathbb{R}^{n\times n}$ such that $(1-\epsilon)\cdot K\preceq\widetilde{K}\preceq(1+\epsilon)\cdot K$, with $\widetilde{K}$ having only $O(\epsilon^{-2}n\log n)$ nonzero entries. [] achieves this in $O(n^{1+o(1)})$ time for multiplicatively Lipschitz kernels when $d=\Theta(\log n)$, while [] overcomes limitations for Gaussian kernels by basing their algorithm on KDE and the magnitude of the minimum entry of the kernel matrix, parameterized by $\tau$. Their algorithm for Gaussian kernels runs in time $\widetilde{O}(nd/\tau^{3.173+o(1)})$. Unfortunately, spectral sparsifiers do not aid our primitives since a sparsifier for a complete graph is an expander with large treewidth, and thus, even though $\widetilde{K}$ is sparse, it does not enable the deployment of our fast solvers.

Besides spectral sparsification, [, ] also demonstrate that with $d=\Theta(\log n)$ and suitable kernels, there exists an $O(n^{1+o(1)})$ time algorithm to multiply the kernel matrix with an arbitrary vector $v\in\mathbb{R}^{n}$. This operation is crucial in Batch KDE as shown in []. Moreover, [] establishes an almost-quadratic lower bound for this operation when the squared dataset radius $B=\omega(\log n)$, assuming SETH. These results rely on computing a rank-$n^{o(1)}$ factorization for the Gaussian kernel matrix. The function $e^{-x}$ can be approximated by a low-degree polynomial of degree

$$ $\displaystyle q:=$ $\displaystyle~{}\Theta(\max\{\sqrt{B\log(1/\epsilon)},\frac{\log(1/\epsilon)}{\log(\log(1/\epsilon)/B)}\})$ $$

for $x\in[0,B]$. Using this polynomial, one can create matrices $U,V$ with rank $\binom{2d+2q}{2q}=n^{o(1)}$ in time $O(n^{1+o(1)})$. Given this factorization, multiplying it with a vector $v$ as $U(V^{\top}v)$ takes $O(n^{1+o(1)})$ time. Let $\widetilde{K}=UV^{\top}$ where $\widetilde{K}_{i,j}=f(\|x_{i}-x_{j}\|_{2}^{2})$, we have for any $(i,j)\in[n]\times[n]$,

$$ $\displaystyle|f(\|x_{i}-x_{j}\|_{2}^{2})-\exp(-\|x_{i}-x_{j}\|_{2}^{2})|\leq$ $\displaystyle~{}\epsilon,$ $$

and for any row $i\in[n]$,

$$ $\displaystyle|(\widetilde{K}v)_{i}-(Kv)_{i}|=$ $\displaystyle~{}|\sum_{j=1}^{n}v_{j}(f(\|x_{i}-x_{j}\|_{2}^{2})-\exp(-\|x_{i}-x_{j}\|_{2}^{2}))|$ $\displaystyle\leq$ $\displaystyle~{}(\max_{j\in[n]}|f(\|x_{i}-x_{j}\|_{2}^{2})-\exp(-\|x_{i}-x_{j}\|_{2}^{2})|)\|v\|_{1}$ $\displaystyle\leq$ $\displaystyle~{}\epsilon\|v\|_{1},$ $$

using Hölder’s inequality. This provides an $\ell_{\infty}$-guarantee of $\|(\widetilde{K}-K)v\|_{\infty}$, useful for Batch Gaussian KDE. Transforming this $\ell_{\infty}$-guarantee into a spectral approximator yields

$$ $\displaystyle(1-\epsilon n)\cdot K\preceq\widetilde{K}\preceq(1+\epsilon n)\cdot K.$ $$

Setting $\epsilon=1/n^{2}$, the low-rank factorization offers an adequate spectral approximation to the exact kernel matrix $K$.

Given $\widetilde{K}=UV^{\top}$ for $U,V\in\mathbb{R}^{n\times n^{o(1)}}$, we can solve program ([3](#S1.E3)) with $\widetilde{K}$ using our low-rank QP algorithm in time $O(n^{1+o(1)}\log(1/\epsilon))$.(^5^55Additional requirement: $B=o(\frac{\log n}{\log\log n})$. See Section [10](#S10) for further discussion.) This represents the first almost-linear time algorithm for Gaussian kernel SVM, even in low-precision settings. The main hurdle has been the lack of advancements in fast kernel matrix algebra primitives and optimization frameworks capable of exploiting these primitives.

The requirements $d=\Theta(\log n)$ and $B=o(\frac{\log n}{\log\log n})$ may seem restrictive, but they are necessary, as no sub-quadratic time algorithm exists for Gaussian kernel SVM without bias when $d=\Theta(\log n)$ and $B=\Omega(\log^{2}n)$, and with bias when $B=\Omega(\log^{6}n)$, assuming SETH. This is based on a reduction from Hamming nearest neighbor to Gaussian kernel SVM, as established by []. Our assumptions on $d$ and $B$ are therefore justified for seeking sub-quadratic or almost-linear time algorithms.

## 3 Discussion

On the algorithmic front, we introduce the first nearly-linear time algorithms for low-rank and low-treewidth convex quadratic programming, leading to nearly-linear time algorithms for structured linear SVMs. For Gaussian kernel SVMs, we utilize a low-rank approximation from [] when $d=\Theta(\log n)$ and the squared dataset radius is small, enabling an almost-linear time algorithm. On the hardness aspect, we establish that when $d=\Theta(\log n)$, if the squared dataset radius is sufficiently large ($\Omega(\log^{2}n)$ without bias and $\Omega(\log^{6}n)$ with bias), then assuming SETH, no sub-quadratic algorithm exists.

Several open problems arise from our work:

Nearly-linear time algorithm for more general sparse QPs. While we assume $Q$ has low-treewidth, an extension of standard sparsity, we note that with a well-behaved kernel function and additional parameters, a spectral sparsifier can be computed in $O(n^{1+o(1)})$ time with $O(n\log n)$ nonzero entries. Can algorithms be designed for sparse kernel matrices without the treewidth assumption?

Better dependence on $k$ for low-rank QPs. Our low-rank QP solver exhibits a dependence of $k^{(\omega+1)/2}$ on the rank $k$. Given the precomputed factorization, can we improve the exponents on $k$? Ideally, an algorithm with nearly-linear dependence on $k$ would align more closely with the input size.

Better dependence on $m$ for general QPs. Focusing on SVMs with a few equality constraints, our QP solvers do not exhibit strong dependence on the number of equality constraints $m$. Without structural assumptions on the constraint matrix $A$, this is expected. However, many QPs, particularly in graph contexts, involve large $m$. Is there a pathway to an algorithm with better dependence on $m$? More broadly, can we achieve a result akin to that of Lee and Sidford for LP [], where the number of iterations depends on the square root of the rank of $A$, with minimal per-iteration cost?

Stronger lower bound in terms of $B$ for Gaussian kernel SVMs. We establish hardness results for Gaussian kernel SVM when $B=\Omega(\log^{2}n)$ without bias and $B=\Omega(\log^{6}n)$ with bias. This contrasts with our algorithm, which requires $B$ to have sub-logarithmic dependence on $n$. For Batch Gaussian KDE, [] demonstrated that fast algorithms are feasible for $B=o(\log n)$, with no sub-quadratic time algorithms for $B=\omega(\log n)$ assuming SETH. Can a stronger lower bound be shown for SVM programs with a bias term, reflecting a more natural setting?

##### Roadmap.

In Section [4](#S4), we present some basic definitions and tools that will be used in the reminder of the paper. In Section [5](#S5), we introduce a few more SVM formulations, including classification and distribution estimation. In Section [6](#S6), we show convex quadratic programming can be reduced to convex empirical risk minimization, and therefore can be solved in the current matrix multiplication time owing to []. In Section [7](#S7) and [8](#S8), we prove results on low-treewidth and low-rank QPs, respectively. In Section [9](#S9), we present a robust IPM framework for QPs, generalize beyond LPs and convex ERMs with linear objective. In Section [10](#S10), we present our algorithms for Gaussian kernel SVMs, with complementary lower bound.

## 4 Preliminary

### 4.1 Notations

For a positive integer $n$, we use $[n]$ to denote the set $\{1,2,\cdots,n\}$. For a matrix $A$, we use $A^{\top}$ to denote its transpose.
For a matrix $A$, we define $\|A\|_{p\rightarrow q}:=\sup_{x}\|Ax\|_{q}/\|x\|_{p}$. When $p=q=2$, we recover the spectral norm.

We define the entrywise $\ell_{p}$-norm of a matrix $A$ as $\|A\|_{p}:=(\sum_{i,j}|A_{i,j}|^{p})^{1/p}$.

For any function $f:\mathbb{N}\rightarrow\mathbb{N}$ and $n\in\mathbb{N}$, we use $\widetilde{O}(f(n))$ to denote $O(f(n)\operatorname{poly}\log f(n))$. We use $\mathbbm{1}\{E\}$ to denote the indicator for event $E$, i.e., if $E$ happens, $\mathbbm{1}\{E\}=1$ and otherwise it’s 0.

### 4.2 Treewidth

Treewidth captures the sparsity and tree-like structures of graphs.

###### Definition 4.1 (Tree Decomposition and Treewidth) .

Let $G=(V,E)$ be a graph, a tree decomposition of $G$ is a tree $T$ with $b$ vertices, and $b$ sets $J_{1},\ldots,J_{b}\subseteq V$ (called bags), satisfying the following properties:

- •
For every edge $(u,v)\in E$, there exists $j\in[b]$ such that $u,v\in J_{j}$;
- •
For every vertex $v\in V$, $\{j\in[b]:v\in J_{j}\}$ is a non-empty subtree of $T$.

The treewidth of $G$ is defined as the minimum value of $\max\{|J_{j}|:j\in[b]\}-1$ over all tree decompositions.

A near-optimal tree decomposition of a graph can be computed in almost linear time.

###### Theorem 4.2 ( [ 4 ] ) .

Given a graph $G$, there is an $O(m^{1+o(1)})$ time algorithm that produces a tree decomposition of $G$ of maximum bag size $O(\tau\log^{3}n)$, where $\tau$ is the actual (unknown) treewidth of $G$.

Therefore, when $\tau=m^{\Theta(1)}$, we can compute an $\widetilde{O}(\tau)$-size tree decomposition in time $O(m\tau^{o(1)})$, which is negligible in the final running time of Theorem [1.5](#S1.Thmtheorem5).

###### Definition 4.1 (Tree Decomposition and Treewidth) .

###### Theorem 4.2 ( [ 4 ] ) .

### 4.3 Sparse Cholesky Decomposition

In this section we state a few results on sparse Cholesky decomposition.
Fast sparsed Cholesky decomposition algorithms are based on the concept of elimination tree, introduced in [].

###### Definition 4.3 (Elimination tree) .

Let $G$ be an undirected graph on $n$ vertices. An elimination tree $\mathcal{T}$ is a rooted tree on $V(G)$ together with an ordering $\pi$ of $V(G)$ such that for any vertex $v$, its parent is the smallest (under $\pi$) element $u$ such that there exists a path $P$ from $v$ to $u$, such that $\pi(w)\leq\pi(v)$ for all $w\in P-u$.

The following lemma relates the elimination tree and the structure of Cholesky factors.

###### Lemma 4.4 ( [ 52 ] ) .

Let $M$ be a PSD matrix and $\mathcal{T}$ be an elimination tree of the adjacency graph of $M$ (i.e., $(i,j)\in E(G)$ iff $M_{i,j}\neq 0$) together with an elimination ordering $\pi$. Let $P$ be the permutation matrix $P_{i,v}=\mathbbm{1}\{v=\pi(i)\}$.
Then the Cholesky factor $L$ of $PMP^{\top}$ (i.e., $PMP^{\top}=LL^{\top}$) satisfies $L_{i,j}\neq 0$ only if $\pi(i)$ is an ancestor of $\pi(j)$.

The following result is the current best result for computing a sparse Cholesky decomposition.

###### Lemma 4.5 ( [ 28 , Lemma 8.4] ) .

Let $M\in\mathbb{R}^{n\times n}$ be a PSD matrix whose adjacency graph has treewidth $\tau$.
Then we can compute the Cholesky factorization $M=LL^{\top}$ in $\widetilde{O}(n\tau^{\omega-1})$ time.

The following result is the current best result for updating a sparse Cholesky decomposition.

###### Lemma 4.6 ( [ 22 ] ) .

Let $M\in\mathbb{R}^{n\times n}$ be a PSD matrix whose adjacency graph has treewidth $\tau$.
Assume that we are given the Cholseky factorization $M=LL^{\top}$.
Let $w\in\mathbb{R}^{n}$ be a vector such that $M+ww^{\top}$ has the same adjacency graph as $M$.
Then we can compute $\Delta_{L}\in\mathbb{R}^{n\times n}$ such that $L+\Delta_{L}$ is the Cholesky factor of $M+ww^{\top}$ in $O(\tau^{2})$ time.

Throughout our algorithm, we need to compute matrix-vector multiplications involving Cholesky factors. We use the following results from [].

###### Lemma 4.7 ( [ 28 , Lemma 4.7] ) .

Let $M\in\mathbb{R}^{n\times n}$ be a PSD matrix whose adjacency graph has treewidth $\tau$.
Assume that we are given the Cholseky factorization $M=LL^{\top}$.
Then we have the following running time for matrix-vector multiplications.

- (i)
For $v\in\mathbb{R}^{n}$, computing $Lv$, $L^{\top}v$, $L^{-1}v$, $L^{-\top}v$ takes $O(n\tau)$ time.
- (ii)
For $v\in\mathbb{R}^{n}$, computing $Lv$ takes $O(\|v\|_{0}\tau)$ time.
- (iii)
For $v\in\mathbb{R}^{n}$, computing $L^{-1}v$ takes $O(\|L^{-1}v\|_{0}\tau)$ time.
- (iv)
For $v\in\mathbb{R}^{n}$, if $v$ is supported on a path in the elimination tree, then computing $L^{-1}v$ takes $O(\tau^{2})$ time.
- (v)
For $v\in\mathbb{R}^{n}$, computing $\mathcal{W}^{\top}v$ takes $O(n\tau)$ time, where $\mathcal{W}=L^{-1}H^{1/2}$ with $H\in\mathbb{R}^{n\times n}$ be a non-negative diagonal matrix.

###### Lemma 4.8 ( [ 28 , Lemma 4.8] ) .

Let $M\in\mathbb{R}^{n\times n}$ be a PSD matrix whose adjacency graph has treewidth $\tau$.
Assume that we are given the Cholseky factorization $M=LL^{\top}$.
Then we have the following running time for matrix-vector multiplications, when we only need result for a subset of coordinates.

- (i)
Let $S$ be a path in the elimination tree whose one endpoint is the root. For $v\in\mathbb{R}^{n}$, computing $(L^{-\top}v)_{S}$ takes $O(\tau^{2})$ time.
- (ii)
For $v\in\mathbb{R}^{n}$, for $i\in[n]$, computing $(\mathcal{W}^{\top}v)_{i}$ takes $O(\tau^{2})$ time, where $\mathcal{W}=L^{-1}H^{1/2}$ with $H\in\mathbb{R}^{n\times n}$ be a non-negative diagonal matrix.

###### Definition 4.3 (Elimination tree) .

###### Lemma 4.4 ( [ 52 ] ) .

###### Lemma 4.5 ( [ 28 , Lemma 8.4] ) .

###### Lemma 4.6 ( [ 22 ] ) .

###### Lemma 4.7 ( [ 28 , Lemma 4.7] ) .

###### Lemma 4.8 ( [ 28 , Lemma 4.8] ) .

### 4.4 Johnson-Lindentrauss Lemma

We recall the Johnson-Lindenstrauss lemma, a powerful algorithmic primitive that reduces dimension while preserving $\ell_{2}$ norms.

###### Lemma 4.9 ( [ 32 ] ) .

Let $\epsilon\in(0,1)$ be the precision parameter. Let $\delta\in(0,1)$ be the failure probability. Let $A\in\mathbb{R}^{m\times n}$ be a real matrix.
Let $r=\epsilon^{-2}\log(mn/\delta)$.
For $R\in\mathbb{R}^{r\times n}$ whose entries are i.i.d $\mathcal{N}(0,\frac{1}{r})$, the following holds with probability at least $1-\delta$:

$$ $\displaystyle(1-\epsilon)\|a_{i}\|_{2}\leq\|Ra_{i}\|_{2}\leq(1+\epsilon)\|a_{i}\|_{2},~{}\forall i\in[m],$ $$

where for a matrix $A$, $a_{i}^{\top}$ denotes the $i$-th row of matrix $A\in\mathbb{R}^{m\times n}$.

###### Lemma 4.9 ( [ 32 ] ) .

### 4.5 Heavy-Light Decomposition

Heavy-light decomposition is useful when one wants to re-balance a binary tree with height $O(\log n)$.

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

Given a rooted tree $\mathcal{T}$ with $n$ vertices, we can construct in $O(n)$ time an ordering $\pi$ of the vertices such that (1) every path in $\mathcal{T}$ can be decomposed into $O(\log n)$ contiguous subseqeuences under $\pi$, and (2) every subtree in $\mathcal{T}$ is a single contiguous subsequence under $\pi$.

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

## 5 SVM Formulations

In this section, we review a list of formulations of SVM. These formulations have been implemented in the LIBSVM library [].

Throughout this section, we use $\phi:\mathbb{R}^{d}\rightarrow\mathbb{R}^{s}$ to denote the feature mapping, ${\sf K}$ to denote the associated kernel function and $K\in\mathbb{R}^{n\times n}$ to denote the kernel matrix. For linear SVM, $\phi$ is just the identity mapping. We will focus on the dual quadratic program formulation as usual. We will also assume for each problem, a dataset $X\in\mathbb{R}^{n\times d}$ is given together with binary labels $y\in\mathbb{R}^{n}$. Let $Q:=(yy^{\top})\circ K$.

### 5.1 C 𝐶 C -Support Vector Classification

This formulation is also referred as the *soft-margin SVM*. It can be viewed as imposing a regularization on the primal program to allow mis-classification.

###### Definition 5.1 ( $C$ -Support Vector Classification) .

Given a parameter $C>0$, the $C$-support vector classification ($C$-SVC) is defined as

$$ $\displaystyle\max_{\alpha\in\mathbb{R}^{n}}~{}$ $\displaystyle~{}{\bf 1}_{n}^{\top}\alpha-\frac{1}{2}\alpha^{\top}Q\alpha$ s.t. $\displaystyle~{}\alpha^{\top}y=0,$ $\displaystyle~{}0\leq\alpha\leq C\cdot{\bf 1}_{n}.$ $$

###### Definition 5.1 ( C 𝐶 C -Support Vector Classification) .

### 5.2 ν 𝜈 \nu -Support Vector Classification

The $C$-SVC (Definition [5.1](#S5.Thmtheorem1)) penalizes large values of $\alpha$ by limiting the magnitude of it. The $\nu$-SVC (Definition [5.2](#S5.Thmtheorem2)) turns ${\bf 1}_{n}^{\top}\alpha$ from an objective into a constraint on $\ell_{1}$ norm.

###### Definition 5.2 ( $\nu$ -Support Vector Classification) .

Given a parameter $\nu>0$, the $\nu$-support vector classification ($\nu$-SVC) is defined as

$$ $\displaystyle\min_{\alpha\in\mathbb{R}^{n}}~{}$ $\displaystyle~{}\frac{1}{2}\alpha^{\top}Q\alpha$ s.t. $\displaystyle~{}\alpha^{\top}y=0,$ $\displaystyle~{}{\bf 1}_{n}^{\top}\alpha=\nu,$ $\displaystyle~{}0\leq\alpha\leq\frac{1}{n}\cdot{\bf 1}_{n}.$ $$

One can interpret this formulation as to find a vector that lives in the orthogonal complement of $y$ that is non-negative, each entry is at most $\frac{1}{n}$ and its $\ell_{1}$ norm is $\nu$. Clearly, we must have $\nu\in(0,1]$. More specifically, let $k_{+}$ be the number of positive labels and $k_{-}$ be the number of negative labels. It is shown by [] that the above problem is feasible if and only if

$$ $\displaystyle\nu\leq$ $\displaystyle~{}\frac{2\min\{k_{-},k_{+}\}}{n}.$ $$

###### Definition 5.2 ( ν 𝜈 \nu -Support Vector Classification) .

### 5.3 Distribution Estimation

SVM is widely-used for predicting binary labels. It can also be used to estimate the support of a high-dimensional distribution. The formulation is similar to $\nu$-SVC, except the PSD matrix $Q$ is *label-less*.

###### Definition 5.3 (Distribution Estimation) .

Given a parameter $\nu>0$, the $\nu$-distribution estimation problem is defined as

$$ $\displaystyle\min_{\alpha\in\mathbb{R}^{n}}~{}$ $\displaystyle~{}\frac{1}{2}\alpha^{\top}K\alpha$ s.t. $\displaystyle~{}0\leq\alpha\leq\frac{1}{n}\cdot{\bf 1}_{n},$ $\displaystyle~{}{\bf 1}_{n}^{\top}\alpha=\nu.$ $$

###### Definition 5.3 (Distribution Estimation) .

### 5.4 ϵ italic-ϵ \epsilon -Support Vector Regression

In addition to classification, support vector framework can also be adapted for regression.

###### Definition 5.4 ( $\epsilon$ -Support Vector Regression) .

Given parameters $\epsilon,C>0$, the $\epsilon$-support vector regression ($\epsilon$-SVR) is defined as

$$ $\displaystyle\min_{\alpha,\alpha^{*}\in\mathbb{R}^{n}}~{}$ $\displaystyle~{}\frac{1}{2}(\alpha-\alpha^{*})^{\top}K(\alpha-\alpha^{*})+\epsilon{\bf 1}_{n}^{\top}(\alpha+\alpha^{*})+y^{\top}(\alpha-\alpha^{*})$ s.t. $\displaystyle~{}{\bf 1}_{n}^{\top}(\alpha-\alpha^{*})=0,$ $\displaystyle~{}0\leq\alpha\leq C\cdot{\bf 1}_{n},$ $\displaystyle~{}0\leq\alpha^{*}\leq C\cdot{\bf 1}_{n}.$ $$

###### Definition 5.4 ( ϵ italic-ϵ \epsilon -Support Vector Regression) .

### 5.5 ν 𝜈 \nu -Support Vector Regression

One can similar adapt the parameter $\nu$ to control the $\ell_{1}$ norm of the regression.

###### Definition 5.5 ( $\nu$ -Support Vector Regression) .

Given parameters $\nu,C>0$, the $\nu$-support vector regression ($\nu$-SVR) is defined as

$$ $\displaystyle\min_{\alpha,\alpha^{*}\in\mathbb{R}^{n}}~{}$ $\displaystyle~{}\frac{1}{2}(\alpha-\alpha^{*})^{\top}K(\alpha-\alpha^{*})+y^{\top}(\alpha-\alpha^{*})$ s.t. $\displaystyle~{}{\bf 1}_{n}^{\top}(\alpha-\alpha^{*})=0,$ $\displaystyle~{}{\bf 1}_{n}^{\top}(\alpha+\alpha^{*})=C\nu,$ $\displaystyle~{}0\leq\alpha\leq\frac{C}{n}\cdot{\bf 1}_{n},$ $\displaystyle~{}0\leq\alpha^{*}\leq\frac{C}{n}\cdot{\bf 1}_{n}.$ $$

###### Definition 5.5 ( ν 𝜈 \nu -Support Vector Regression) .

### 5.6 One Linear Constraint

We classify $C$-SVC (Definition [5.1](#S5.Thmtheorem1)), $\epsilon$-SVR (Definition [5.4](#S5.Thmtheorem4)) and $\nu$-distribution estimation (Definition [5.3](#S5.Thmtheorem3)) into the following generic form:

$$ $\displaystyle\min_{\alpha\in\mathbb{R}^{n}}~{}$ $\displaystyle~{}\frac{1}{2}\alpha^{\top}Q\alpha+p^{\top}\alpha$ s.t. $\displaystyle~{}\alpha^{\top}y=\Delta$ $\displaystyle~{}0\leq\alpha\leq C\cdot{\bf 1}_{n}.$ $$

Note that $C$-SVC (Definition [5.1](#S5.Thmtheorem1)) and distribution estimation (Definition [5.3](#S5.Thmtheorem3)) are readily in this form. For $\epsilon$-SVR (Definition [5.4](#S5.Thmtheorem4)), we need to perform a simple transformation:

Set $\widehat{\alpha}=\begin{bmatrix}\alpha\\
\alpha^{*}\end{bmatrix}\in\mathbb{R}^{2n}$, then it can be written as

$$ $\displaystyle\min_{\widehat{\alpha}\in\mathbb{R}^{2n}}~{}$ $\displaystyle~{}\frac{1}{2}\widehat{\alpha}^{\top}\begin{bmatrix}Q&-Q\\ -Q&Q\end{bmatrix}\widehat{\alpha}+\begin{bmatrix}\epsilon{\bf 1}_{n}+y\\ \epsilon{\bf 1}_{n}-y\end{bmatrix}^{\top}\widehat{\alpha}$ s.t. $\displaystyle~{}\begin{bmatrix}{\bf 1}_{n}\\ -{\bf 1}_{n}\end{bmatrix}^{\top}\widehat{\alpha}=0$ $\displaystyle~{}0\leq\widehat{\alpha}\leq C\cdot{\bf 1}_{2n}.$ $$

### 5.7 Two Linear Constraints

Both $\nu$-SVC (Definition [5.2](#S5.Thmtheorem2)) and $\nu$-SVR (Definition [5.5](#S5.Thmtheorem5)) require one extra constraint. They can be formulated as follows:

$$ $\displaystyle\min_{\alpha\in\mathbb{R}^{n}}~{}$ $\displaystyle~{}\frac{1}{2}\alpha^{\top}Q\alpha+p^{\top}\alpha$ s.t. $\displaystyle~{}{\bf 1}_{n}^{\top}\alpha=\Delta_{1},$ $\displaystyle~{}y^{\top}\alpha=\Delta_{2},$ $\displaystyle~{}0\leq\alpha\leq C\cdot{\bf 1}_{n}.$ $$

For $\nu$-SVR (Definition [5.5](#S5.Thmtheorem5)), one can leverage a similar transformation as $\epsilon$-SVR (Definition [5.4](#S5.Thmtheorem4)).

###### Remark 5.6 .

All variants of SVM-related quadratic programs can all be solved using our QP solvers for three cases:

- •
Linear SVM with $n\gg d$, we can solve it in $\widetilde{O}(nd^{(\omega+1)/2}\log(1/\epsilon))$ time;
- •
Linear SVM with a small treewidth decomposition with width $\tau$ on $XX^{\top}$, we can solve it in $\widetilde{O}(n\tau^{(\omega+1)/2}\log(1/\epsilon))$ time;
- •
Gaussian kernel SVM with $d=\Theta(\log n)$ and $B=o(\frac{\log n}{\log\log n})$, we can solve it in $O(n^{1+o(1)}\log(1/\epsilon))$ time.

Even though our solvers have relatively bad dependence on the number of equality constraints, for all these SVM formulations, at most 2 equality constraints are presented and thus can be solved very fast.

###### Remark 5.6 .

## 6 Algorithms for General QP

In this section, we discuss algorithms for general (convex) quadratic programming.
We show that they can be solved in the current matrix multiplication time via reduction to linear programming with convex constraints [].

### 6.1 LCQP in the Current Matrix Multiplication Time

###### Proposition 6.1 .

There is an algorithm that solves LCQP (Definition [1.1](#S1.Thmtheorem1)) up to $\epsilon$ error in $\widetilde{O}((n^{\omega}+n^{2.5-\alpha/2}+n^{2+1/6})\log(1/\epsilon))$ time, where $\omega\leq 2.373$ is the matrix multiplication constant and $\alpha\geq 0.32$ is the dual matrix multiplication constant.

###### Proof.

Let $Q=PDP^{\top}$ be an eigen-decomposition of $Q$ where $D$ is diagonal and $P$ is orthogonal.
Let $\widetilde{x}:=P^{-1}x$. Then it suffices to solve

$$ $\displaystyle\min~{}$ $\displaystyle~{}\frac{1}{2}\widetilde{x}^{\top}D\widetilde{x}+c^{\top}P\widetilde{x}$ $\displaystyle\mathrm{s.t.~{}}$ $\displaystyle~{}AP\widetilde{x}=b$ $\displaystyle~{}P\widetilde{x}\geq 0.$ $$

By adding $n$ non-negative variables and $n$ constraints $x=P\widetilde{x}$ we can make all constraints equality constraints. There are $n$ non-negative variables and $n$ unconstrained variables.
If we want to ensure all variables are non-negative, we need to split every coordinate of $\widetilde{x}$ into two.
In this way the coefficient matrix $Q$ will be block diagonal with block size $2$.

We perform the above reduction, and assume that we have a program of form ([1](#S1.E1)) where $Q$ is diagonal.
Let $q_{i}:=Q_{i,i}$ be the $i$-th element on the diagonal.
Then the QP is equivalent to the following program

$$ $\displaystyle\min~{}$ $\displaystyle~{}c^{\top}x+q^{\top}t$ $\displaystyle\mathrm{s.t.~{}}$ $\displaystyle~{}Ax=b$ $\displaystyle~{}t_{i}\geq\frac{1}{2}x_{i}^{2}\qquad\forall i\in[n]$ $\displaystyle~{}x\geq 0$ $$

Note that the set $\{(x_{i},t_{i})\in\mathbb{R}^{2}:t_{i}\geq\frac{1}{2}x_{i}^{2}\}$ is a convex set. So we can apply [] here with $n$ variables, each in the convex set $\{(a,b)\in\mathbb{R}^{2}:a\geq 0,b\geq\frac{1}{2}a^{2}\}$.
∎

###### Proposition 6.1 .

###### Proof.

### 6.2 Algorithm for QCQP

Our algorithm for LCQP in the previous section can be generalized to quadratically constrained quadratic programs (QCQP).
QCQP is defined as follows.

###### Definition 6.2 (QCQP) .

Let $Q_{0},\ldots,Q_{m}\in\mathbb{R}^{n\times n}$ be PSD matrices.
Let $q_{0},\ldots,q_{m}\in\mathbb{R}^{n}$.
Let $r\in\mathbb{R}^{m}$.
Let $A\in\mathbb{R}^{d\times n}$, $b\in\mathbb{R}^{d}$.
The quadratically constrained quadratic programming (QCQP) problem asks the solve the following program.

$$ $\displaystyle\min_{x\in\mathbb{R}^{n}}~{}$ $\displaystyle~{}\frac{1}{2}x^{\top}Q_{0}x+q_{0}^{\top}x$ $\displaystyle\mathrm{s.t.~{}}$ $\displaystyle~{}\frac{1}{2}x^{\top}Q_{i}x+q_{i}^{\top}x+r_{i}\leq 0\qquad\forall i\in[m]$ $\displaystyle~{}Ax=b$ $\displaystyle~{}x\geq 0$ $$

###### Proposition 6.3 .

There is an algorithm that solves QCQP (Definition [6.2](#S6.Thmtheorem2)) up to $\epsilon$ error in $\widetilde{O}(((mn)^{\omega}+(mn)^{2.5-\alpha/2}+(mn)^{2+1/6})\log(1/\epsilon))$ time, where $\omega\leq 2.373$ is the matrix multiplication constant and $\alpha\geq 0.32$ is the dual matrix multiplication constant.

###### Proof.

We first rewrite the program as following.

$$ $\displaystyle\min~{}$ $\displaystyle~{}-r_{0}$ $\displaystyle\mathrm{s.t.~{}}$ $\displaystyle~{}\frac{1}{2}x^{\top}Q_{i}x+q_{i}^{\top}x+r_{i}\leq 0\qquad\forall 0\leq i\leq m$ $\displaystyle~{}Ax=b$ $\displaystyle~{}x\geq 0$ $$

Write $Q_{i}=P_{i}D_{i}P_{i}^{\top}$ be an eigen-decomposition of $Q_{i}$ where $D_{i}$ is diagonal and $P_{i}$ is orthogonal.
Let $\widetilde{x}_{i}\in\mathbb{R}^{n}$ be defined as $\widetilde{x}_{i}:=P_{i}^{-1}x$.
Then we can rewrite the program as

$$ $\displaystyle\min~{}$ $\displaystyle~{}-r_{0}$ $\displaystyle\mathrm{s.t.~{}}$ $\displaystyle~{}\frac{1}{2}\widetilde{x}_{i}^{\top}D_{i}\widetilde{x}_{i}+q_{i}^{\top}P_{i}\widetilde{x}_{i}+r_{i}\leq 0\qquad\forall 0\leq i\leq m$ $\displaystyle~{}Ax=b$ $\displaystyle~{}\widetilde{x}_{i}=P_{i}^{-1}x$ $\displaystyle~{}x\geq 0$ $$

For $0\leq i\leq m$ and $j\in[n]$, define variable $t_{i,j}\in\mathbb{R}$.
Then we can rewrite the program as

$$ $\displaystyle\min~{}$ $\displaystyle~{}-r_{0}$ $\displaystyle\mathrm{s.t.~{}}$ $\displaystyle~{}\sum_{j\in[n]}D_{i,(j,j)}t_{i,j}+q_{i}^{\top}P_{i}\widetilde{x}_{i}+r_{i}\leq 0\qquad\forall 0\leq i\leq m$ $\displaystyle~{}Ax=b$ $\displaystyle~{}\widetilde{x}_{i}=P_{i}^{-1}x$ $\displaystyle~{}t_{i,j}\geq\widetilde{x}_{i,j}^{2}$ $\displaystyle~{}x\geq 0$ $$

We can consider $(\widetilde{x}_{i,j},t_{i,j})_{0\leq i\leq m,j\in[n]}$ as $(m+1)n$ variables in the convex set $\{(a,b):b\geq\frac{1}{2}a^{2}\}$.
Then we can apply [].
∎

###### Definition 6.2 (QCQP) .

###### Proposition 6.3 .

###### Proof.

## 7 Algorithm for Low-Treewidth QP

In this section we present a nearly-linear time algorithm for solving low-treewidth QP with low number of linear constraints. We briefly describe the outline of this section.

- •
In Section [7.1](#S7.SS1), we present the main statement of Section [7](#S7).
- •
In Section [7.2](#S7.SS2), we present the main data structure CentralPathMaintenance.
- •
In Section [7.3](#S7.SS3), we present several data structures used in CentralPathMaintenance, including ExactDS (Section [7.3.1](#S7.SS3.SSS1)), ApproxDS (Section [7.3.2](#S7.SS3.SSS2)), BatchSketch (Section [7.3.3](#S7.SS3.SSS3)), VectorSketch (Section [7.3.4](#S7.SS3.SSS4)), BalancedSketch (Section [7.3.5](#S7.SS3.SSS5)).
- •
In Section [7.4](#S7.SS4), we prove correctness and running time of CentralPathMaintenance data structure.
- •
In Section [7.5](#S7.SS5), we prove the main result (Theorem [7.1](#S7.Thmtheorem1)).

### 7.1 Main Statement

We consider programs of the form ([16](#S9.E16)), i.e.,

$$ $\displaystyle\min_{x\in\mathbb{R}^{n}}~{}$ $\displaystyle~{}\frac{1}{2}x^{\top}Qx+c^{\top}x$ $\displaystyle\mathrm{s.t.~{}}$ $\displaystyle~{}Ax=b$ $\displaystyle~{}x_{i}\in\mathcal{K}_{i}\qquad\forall i\in[n]$ $$

where $Q\in\mathcal{S}^{n_{\mathrm{tot}}}$, $c\in\mathbb{R}^{n_{\mathrm{tot}}}$,
$A\in\mathbb{R}^{m\times n_{\mathrm{tot}}}$, $b\in\mathbb{R}^{m}$, $\mathcal{K}_{i}\subset\mathbb{R}^{n_{i}}$ is a convex set.
For simplicity, we assume that $n_{i}=O(1)$ for all $i\in[n]$.

###### Theorem 7.1 .

Consider the convex program ([16](#S9.E16)).
Let $\phi_{i}:\mathcal{K}_{i}\to\mathbb{R}$ be a $\nu_{i}$-self-concordant barrier for all $i\in[n]$.
Suppose the program satisfies the following properties:

- •
Inner radius $r$: There exists $z\in\mathbb{R}^{n_{\mathrm{tot}}}$ such that $Az=b$ and $B(z,r)\in\mathcal{K}$.
- •
Outer radius $R$: $\mathcal{K}\subseteq B(0,R)$ where $0\in\mathbb{R}^{n_{\mathrm{tot}}}$.
- •
Lipschitz constant $L$: $\|Q\|_{2\to 2}\leq L$, $\|c\|_{2}\leq L$.
- •
Treewidth $\tau$: Treewidth (Definition [4.1](#S4.Thmtheorem1)) of the adjacency graph of $Q$ is at most $\tau$.

Let $(w_{i})_{i\in[n]}\in\mathbb{R}_{\geq 1}^{n}$ and $\kappa=\sum_{i\in[n]}w_{i}\nu_{i}$.
Given any $0<\epsilon\leq\frac{1}{2}$, we can find an approximate solution $x\in\mathcal{K}$ satisfiying

$$ $\displaystyle\frac{1}{2}x^{\top}Qx+c^{\top}x$ $\displaystyle\leq\min_{Ax=b,x\in\mathcal{K}}\left(\frac{1}{2}x^{\top}Qx+c^{\top}x\right)+\epsilon LR(R+1),$ $\displaystyle\|Ax-b\|_{1}$ $\displaystyle\leq 3\epsilon(R\|A\|_{1}+\|b\|_{1}),$ $$

in expected time

$$ $\displaystyle\widetilde{O}((\sqrt{\kappa}n^{-1/2}+\log(R/(r\epsilon)))\cdot n(\tau^{2}m+\tau m^{2})^{1/2}(\tau^{\omega-1}+\tau m+m^{\omega-1})^{1/2}).$ $$

When $\max_{i\in[n]}\nu_{i}=\widetilde{O}(1)$, $w_{i}=1$, $m=\widetilde{O}(\tau^{\omega-2})$, the running time simplifies to

$$ $\displaystyle\widetilde{O}(n\tau^{(\omega+1)/2}m^{1/2}\log(R/(r\epsilon))).$ $$

###### Theorem 7.1 .

### 7.2 Algorithm Structure and Central Path Maintenance

Our algorithm is based on the robust Interior Point Method (robust IPM). Details of the robust IPM will be given in Section [9](#S9).
During the algorithm, we maintain a primal-dual solution pair $(x,s)\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$ on the robust central path. In addition, we maintain a sparsely-changing approximation $(\overline{x},\overline{s})\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$ to $(x,s)$.
In each iteration, we implicitly perform update

$$ $\displaystyle x$ $\displaystyle\leftarrow x+\overline{t}B_{w,\overline{x},\overline{t}}^{-1/2}(I-P_{w,\overline{x},\overline{t}})B_{w,\overline{x},\overline{t}}^{-1/2}\delta_{\mu}$ $\displaystyle s$ $\displaystyle\leftarrow s+\overline{t}\delta_{\mu}-\overline{t}^{2}H_{w,\overline{x}}B_{w,\overline{x},\overline{t}}^{-1/2}(I-P_{w,\overline{x},\overline{t}})B_{w,\overline{x},\overline{t}}^{-1/2}\delta_{\mu}$ $$

where

$$ $\displaystyle H_{w,\overline{x}}$ $\displaystyle=\nabla^{2}\phi_{w}(\overline{x})$ (see Eq. ( 24 )) $\displaystyle B_{w,\overline{x},\overline{t}}$ $\displaystyle=Q+\overline{t}H_{w,\overline{x}}$ (see Eq. ( 25 )) $\displaystyle P_{w,\overline{x},\overline{t}}$ $\displaystyle=B_{w,\overline{x},\overline{t}}^{-1/2}A^{\top}(AB_{w,\overline{x},\overline{t}}^{-1}A^{\top})^{-1}AB_{w,\overline{x},\overline{t}}^{-1/2}$ (see Eq. ( 26 )) $$

and explicitly maintain $(\overline{x},\overline{s})$ such that they remain close to $(x,s)$ in $\ell_{\infty}$-distance.

This task is handled by the CentralPathMaintenance data structure, which is our main data structure. The robust IPM algorithm (Algorithm [19](#alg19), [20](#alg20)) directly calls it in every iteration.

The CentralPathMaintenance data structure (Algorithm [1](#alg1)) has two main sub data structures, ExactDS (Algorithm [2](#alg2), [3](#alg3)) and ApproxDS (Algorithm [4](#alg4), [5](#alg5)).
ExactDS is used to maintain $(x,s)$, and ApproxDS is used to maintain $(\overline{x},\overline{s})$.

Figure: Algorithm 1 Our main data structure for low-treewidth QP solver.

###### Theorem 7.2 .

Data structure CentralPathMaintenance (Algorithm [1](#alg1)) implicitly maintains the central path primal-dual solution pair $(x,s)\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$ and explicitly maintains its approximation $(\overline{x},\overline{s})\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$ using the following functions:

- •
Initialize$(x\in\mathbb{R}^{n_{\mathrm{tot}}},s\in\mathbb{R}^{n_{\mathrm{tot}}},t_{0}\in\mathbb{R}_{>0},\epsilon\in(0,1))$: Initializes the data structure with initial primal-dual solution pair $(x,s)\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$, initial central path timestamp $t_{0}\in\mathbb{R}_{>0}$ in $\widetilde{O}(n(\tau^{\omega-1}+\tau m+m^{\omega-1}))$ time.
- •
MultiplyAndMove$(t\in\mathbb{R}_{>0})$: It implicitly maintains
$\displaystyle x$
$\displaystyle\leftarrow x+\overline{t}B_{w,\overline{x},\overline{t}}^{-1/2}(I-P_{w,\overline{x},\overline{t}})B_{w,\overline{x},\overline{t}}^{-1/2}\delta_{\mu}(\overline{x},\overline{s},\overline{t})$
$\displaystyle s$
$\displaystyle\leftarrow s+\overline{t}\delta_{\mu}-\overline{t}^{2}H_{w,\overline{x}}B_{w,\overline{x},\overline{t}}^{-1/2}(I-P_{w,\overline{x},\overline{t}})B_{w,\overline{x},\overline{t}}^{-1/2}\delta_{\mu}(\overline{x},\overline{s},\overline{t})$
where $H_{w,\overline{x}}$, $B_{w,\overline{x},\overline{t}}$, $P_{w,\overline{x},\overline{t}}$ are defined in Eq. ([24](#S9.E24))([25](#S9.E25))([26](#S9.E26)) respectively, and $\overline{t}$ is some timestamp satisfying $|\overline{t}-t|\leq\epsilon_{t}\cdot\overline{t}$.
It also explicitly maintains $(\overline{x},\overline{s})\in\mathbb{R}^{n_{\mathrm{tot}}\times n_{\mathrm{tot}}}$ such that $\|\overline{x}_{i}-x_{i}\|_{\overline{x}_{i}}\leq\overline{\epsilon}$ and $\|\overline{s}_{i}-s_{i}\|^{*}_{\overline{x}_{i}}\leq t\overline{\epsilon}w_{i}$ for all $i\in[n]$ with probability at least $0.9$.
Assuming the function is called at most $N$ times and $t$ decreases from $t_{\max}$ to $t_{\min}$, the total running time is
$\displaystyle\widetilde{O}((Nn^{-1/2}+\log(t_{\max}/t_{\min}))\cdot n(\tau^{2}m+\tau m^{2})^{1/2}(\tau^{\omega-1}+\tau m+m^{\omega-1})^{1/2}).$
- •
Output: Computes $(x,s)\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$ exactly and outputs them in
$\widetilde{O}(n\tau m)$
time.

###### Theorem 7.2 .

### 7.3 Data Structures Used in CentralPathMaintenance

In this section we present several data structures used in CentralPathMaintenance, including:

- •
ExactDS (Section [7.3.1](#S7.SS3.SSS1)): This data structure maintains an implicit representation of the primal-dual solution pair $(x,s)$. This is directly used by CentralPathMaintenance.
- •
ApproxDS (Section [7.3.2](#S7.SS3.SSS2)): This data structure explicitly maintains an approximation $(\overline{x},\overline{s})$ of $(x,s)$. This data structure is directly used by CentralPathMaintenance.
- •
BatchSketch (Section [7.3.3](#S7.SS3.SSS3)): This data structure maintains a sketch of $(x,s)$. This data structure is used by ApproxDS.
- •
VectorSketch (Section [7.3.4](#S7.SS3.SSS4)): This data structure maintains a sketch of sparsely-changing vectors. This data structure is used by BatchSketch.
- •
BalancedSketch (Section [7.3.5](#S7.SS3.SSS5)): This data structure maintains a sketch of vectors of form $\mathcal{W}^{\top}v$, where $v$ is sparsely-changing. This data structure is used by BatchSketch.

Notation: In this section, for simplicity, we write $B_{\overline{x}}$ for $B_{w,\overline{x},\overline{t}}$, and $L_{\overline{x}}$ for the Cholesky factor of $B_{\overline{x}}$, i.e., $B_{\overline{x}}=L_{\overline{x}}L_{\overline{x}}^{\top}$.

#### 7.3.1 ExactDS

In this section we present the data structure ExactDS. It maintains an implicit representation of the primal-dual solution pair $(x,s)$ by maintaining several sparsely-changing vectors (see Eq. ([8](#S7.E8))([9](#S7.E9))).
This data structure has a similar spirit as ExactDS in [], but we have a different representation from the previous works because we are working with quadratic programming rather than linear programming.

###### Theorem 7.3 .

Data structure ExactDS (Algorithm [2](#alg2), [3](#alg3)) implicitly maintains the primal-dual pair $(x,s)\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$, computable via the expression

$$ $\displaystyle x$ $\displaystyle=\widehat{x}+H_{w,\overline{x}}^{-1/2}\mathcal{W}^{\top}(h\beta_{x}-\widetilde{h}\widetilde{\beta}_{x}+\epsilon_{x}),$ (8) $\displaystyle s$ $\displaystyle=\widehat{s}+H_{w,\overline{x}}^{1/2}c_{s}\beta_{c_{s}}-H_{w,\overline{x}}^{1/2}\mathcal{W}^{\top}(h\beta_{s}-\widetilde{h}\widetilde{\beta}_{s}+\epsilon_{s}),$ (9) $$

where
$\widehat{x},\widehat{s}\in\mathbb{R}^{n_{\mathrm{tot}}}$, $\mathcal{W}=L_{\overline{x}}^{-1}H_{w,\overline{x}}^{1/2}\in\mathbb{R}^{n_{\mathrm{tot}}\times n_{\mathrm{tot}}}$, $h=L_{\overline{x}}^{-1}\overline{\delta}_{\mu}\in\mathbb{R}^{n_{\mathrm{tot}}}$, $c_{s}=H_{w,\overline{x}}^{-1/2}\overline{\delta}_{\mu}\in\mathbb{R}^{n_{\mathrm{tot}}}$ $\beta_{x},\beta_{s},\beta_{c_{s}}\in\mathbb{R}$, $\widetilde{h}=L_{\overline{x}}^{-1}A^{\top}\in\mathbb{R}^{n_{\mathrm{tot}}\times m}$, $\widetilde{\beta}_{x},\widetilde{\beta}_{s}\in\mathbb{R}^{m}$, $\epsilon_{x},\epsilon_{s}\in\mathbb{R}^{n_{\mathrm{tot}}}$.

The data structure supports the following functions:

- •
Initialize$(x,s,\overline{x},\overline{s}\in\mathbb{R}^{n_{\mathrm{tot}}},\overline{t}\in\mathbb{R}_{>0})$: Initializes the data structure in
$\widetilde{O}(n\tau^{\omega-1}+n\tau m+nm^{\omega-1})$
time, with initial value of the primal-dual pair $(x,s)$, its initial approximation $(\overline{x},\overline{s})$, and initial approximate timestamp $\overline{t}$.
- •
Move$()$: Performs robust central path step
$\displaystyle x$
$\displaystyle\leftarrow x+\overline{t}B_{\overline{x}}^{-1}\delta_{\mu}-\overline{t}B_{\overline{x}}^{-1}A^{\top}(AB_{\overline{x}}^{-1}A^{\top})^{-1}AB_{\overline{x}}^{-1}\delta_{\mu},$
(10)
$\displaystyle s$
$\displaystyle\leftarrow s+\overline{t}\delta_{\mu}-\overline{t}^{2}B_{\overline{x}}^{-1}\delta_{\mu}+\overline{t}^{2}B_{\overline{x}}^{-1}A^{\top}(AB_{\overline{x}}^{-1}A^{\top})^{-1}AB_{\overline{x}}^{-1}\delta_{\mu}$
(11)
in $O(m^{\omega})$ time by updating its implicit representation.
- •
Update$(\delta_{\overline{x}},\delta_{\overline{s}}\in\mathbb{R}^{n_{\mathrm{tot}}})$: Updates the approximation pair $(\overline{x},\overline{s})$ to $(\overline{x}^{\mathrm{new}}=\overline{x}+\delta_{\overline{x}}\in\mathbb{R}^{n_{\mathrm{tot}}},\overline{s}^{\mathrm{new}}=\overline{s}+\delta_{\overline{s}}\in\mathbb{R}^{n_{\mathrm{tot}}})$ in
$\widetilde{O}((\tau^{2}m+\tau m^{2})(\|\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{s}}\|_{0}))$
time, and output the changes in variables $\delta_{H_{w,\overline{x}}^{1/2}\widehat{x}}$, $\delta_{h}$, $\delta_{\beta_{x}}$, $\delta_{\widetilde{h}}$, $\delta_{\widetilde{\beta}_{x}}$, $\delta_{\epsilon_{x}}$, $\delta_{H_{w,\overline{x}}^{-1/2}\widehat{s}}$, $\delta_{\beta_{s}}$, $\delta_{\widetilde{\beta}_{s}}$, $\delta_{\epsilon_{s}}$.
Furthermore, $h,\epsilon_{x},\epsilon_{s}$ change in $O(\tau(\|\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{s}}\|_{0}))$ coordinates,
$\widetilde{h}$ changes in $\widetilde{O}(\tau m(\|\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{s}}\|_{0}))$ coordinates,
and $H_{\overline{x}}^{1/2}\widehat{x},H_{\overline{x}}^{-1/2}\widehat{s},c_{s}$ change in $O(\|\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{s}}\|_{0})$ coordinates.
- •
Output$()$: Output $x$ and $s$ in $\widetilde{O}(n\tau m)$ time.
- •
Query$x(i\in[n])$: Output $x_{i}$ in $\widetilde{O}(\tau^{2}m)$ time.
This function is used by ApproxDS.
- •
Query$s(i\in[n])$: Output $s_{i}$ in $\widetilde{O}(\tau^{2}m)$ time.
This function is used by ApproxDS.

Figure: Algorithm 2 The ExactDS data structure used in Algorithm [1](#alg1).

Figure: Algorithm 3 Algorithm [2](#alg2) continued.

###### Proof of Theorem 7.3 .

By combining Lemma [7.4](#S7.Thmtheorem4) and [7.5](#S7.Thmtheorem5).
∎

###### Lemma 7.4 .

ExactDS correctly maintains an implicit representation of $(x,s)$, i.e., invariant

$$ $\displaystyle x=\widehat{x}+H_{w,\overline{x}}^{-1/2}\mathcal{W}^{\top}(h\beta_{x}-\widetilde{h}\widetilde{\beta}_{x}+\epsilon_{x}),$ $\displaystyle s=\widehat{s}+H_{w,\overline{x}}^{1/2}c_{s}\beta_{c_{s}}-H_{w,\overline{x}}^{1/2}\mathcal{W}^{\top}(h\beta_{s}-\widetilde{h}\widetilde{\beta}_{s}+\epsilon_{s}),$ $\displaystyle h=L_{\overline{x}}^{-1}\overline{\delta}_{\mu},\qquad c_{s}=H_{w,\overline{x}}^{-1/2}\overline{\delta}_{\mu},\qquad\widetilde{h}=L_{\overline{x}}^{-1}A^{\top},$ $\displaystyle\widetilde{u}=\widetilde{h}^{\top}\widetilde{h},\qquad u=\widetilde{h}^{\top}h,$ $\displaystyle\overline{\alpha}=\sum_{i\in[n]}w_{i}^{-1}\cosh^{2}(\frac{\lambda}{w_{i}}\gamma_{i}(\overline{x},\overline{s},\overline{t})),$ $\displaystyle\overline{\delta}_{\mu}=\overline{\alpha}^{1/2}\delta_{\mu}(\overline{x},\overline{s},\overline{t})$ $$

always holds after every external call, and return values of the queries are correct.

###### Proof.

Initialize:
By checking the definitions we see that all invariants are satisfied after Initialize.

Move:
By comparing the implicit representation ([8](#S7.E8))([9](#S7.E9)) and the robust central path step ([10](#S7.E10))([11](#S7.E11)), we see that Move updates $(x,s)$ correctly.

Update:
We would like to prove that Update correctly updates the values of $h$, $c_{s}$, $\widetilde{h}$, $\widetilde{u}$, $u$, $\overline{\alpha}$, $\overline{\delta}_{\mu}$, while preserving the values of $(x,s)$.

First note that $H_{w,\overline{x}}$, $B_{\overline{x}}$, $L_{\overline{x}}$ are updated correctly.
The remaining updates are separated into two steps: Update$h$ and Update$h$.

Step Update$h$:
The first few lines of Update$h$ updates $\overline{\alpha}$ and $\overline{\delta}_{\mu}$ correctly.

We define $H_{w,\overline{x}}^{\mathrm{new}}:=H_{w,\overline{x}}+\Delta_{H_{w,\overline{x}}}$, $B_{\overline{x}}^{\mathrm{new}}:=B_{\overline{x}}+\Delta_{B_{\overline{x}}}$, $L_{\overline{x}}^{\mathrm{new}}:=L_{\overline{x}}+\Delta_{L_{\overline{x}}}$, $\overline{\delta}_{\mu}^{\mathrm{new}}:=\overline{\delta}_{\mu}+\delta_{\overline{\delta}_{\mu}}$, and so on.
Immediately after Algorithm [3](#alg3), Line [26](#alg3.l26), we have

$$ $\displaystyle h+\delta_{h}$ $\displaystyle=L_{\overline{x}}^{-1}\overline{\delta}_{\mu}+L_{\overline{x}}^{-1}\delta_{\overline{\delta}_{\mu}}-(L_{\overline{x}}+\Delta_{L_{\overline{x}}})^{-1}\Delta_{L_{\overline{x}}}(L_{\overline{x}}^{-1}\overline{\delta}_{\mu}+L_{\overline{x}}^{-1}\delta_{\overline{\delta}_{\mu}})$ $\displaystyle=(L_{\overline{x}}^{-1}-(L_{\overline{x}}+\Delta_{L_{\overline{x}}})^{-1}\Delta_{L_{\overline{x}}}L_{\overline{x}}^{-1})\overline{\delta}_{\mu}^{\mathrm{new}}$ $\displaystyle=L_{\overline{x}}^{\mathrm{new}}\overline{\delta}_{\mu}^{\mathrm{new}},$ $\displaystyle c_{s}+\delta_{c_{s}}$ $\displaystyle=H_{w,\overline{x}}^{-1/2}\overline{\delta}_{\mu}+\Delta_{H_{w,\overline{x}}^{-1/2}}(\overline{\delta}_{\mu}+\delta_{\overline{\delta}_{\mu}}+H_{w,\overline{x}}^{-1/2}\delta_{\overline{\delta}_{\mu}}$ $\displaystyle=(H_{w,\overline{x}}^{\mathrm{new}})^{-1/2}\overline{\delta}_{\mu}^{\mathrm{new}},$ $\displaystyle\widetilde{h}+\delta_{\widetilde{h}}$ $\displaystyle=L_{\overline{x}}^{-1}A^{\top}-(L_{\overline{x}}+\Delta_{L_{\overline{x}}})^{-1}\Delta_{L_{\overline{x}}}A^{\top}$ $\displaystyle=(L_{\overline{x}}^{-1}-(L_{\overline{x}}+\Delta_{L_{\overline{x}}})^{-1}\Delta_{L_{\overline{x}}}L_{\overline{x}}^{-1})A^{\top}$ $\displaystyle=L_{\overline{x}}^{\mathrm{new}}A^{\top}.$ $$

So $h,c_{s},\widetilde{h}$ are updated correctly.
Also

$$ $\displaystyle\widetilde{u}+\delta_{\widetilde{u}}$ $\displaystyle=\widetilde{h}^{\top}\widetilde{h}+\delta_{\widetilde{h}}^{\top}(\widetilde{h}+\delta_{\widetilde{h}})+\widetilde{h}^{\top}\delta_{\widetilde{h}}=(\widetilde{h}+\delta_{\widetilde{h}})^{\top}(\widetilde{h}+\delta_{\widetilde{h}}),$ $\displaystyle u+\delta_{u}$ $\displaystyle=\widetilde{h}^{\top}h+\delta_{\widetilde{h}}^{\top}(h+\delta_{h})+\widetilde{h}^{\top}\delta_{h}=(\widetilde{h}+\delta_{\widetilde{h}})^{\top}(h+\delta_{h}).$ $$

So $\widetilde{u}$ and $u$ are maintained correctly.
Furthermore, immediately after Algorithm [3](#alg3), Line [26](#alg3.l26), we have

$$ $\displaystyle~{}(\widehat{x}+L_{\overline{x}}^{-\top}(h^{\mathrm{new}}\beta_{x}-\widetilde{h}^{\mathrm{new}}\widetilde{\beta}_{x}+\epsilon_{x}^{\mathrm{new}}))-(\widehat{x}+L_{\overline{x}}^{-\top}(h\beta_{x}-\widetilde{h}\widetilde{\beta}_{x}+\epsilon_{x}))$ $\displaystyle=$ $\displaystyle~{}L_{\overline{x}}^{-\top}(\delta_{h}\beta_{x}-\delta_{\widetilde{h}}\widetilde{\beta}_{s}+\delta_{\epsilon_{x}})$ $\displaystyle=$ $\displaystyle~{}0.$ $$

Therefore, after Update$h$ finishes, we have

$$ $\displaystyle x=\widehat{x}+L_{\overline{x}}^{-\top}(h\beta_{x}-\widetilde{h}\widetilde{\beta}_{x}+\epsilon_{x}).$ $$

For $s$, we have

$$ $\displaystyle~{}(\widehat{s}^{\mathrm{new}}+(H_{w,\overline{x}}^{\mathrm{new}})^{1/2}c_{s}^{\mathrm{new}}\beta_{c_{s}}-L_{\overline{x}}^{-\top}(h^{\mathrm{new}}\beta_{s}-\widetilde{h}^{\mathrm{new}}\widetilde{\beta}_{s}+\epsilon_{s}^{\mathrm{new}}))$ $\displaystyle-(\widehat{s}+H_{w,\overline{x}}^{1/2}c_{s}\beta_{c_{s}}-L_{\overline{x}}^{-\top}(h\beta_{s}-\widetilde{h}\widetilde{\beta}_{s}+\epsilon_{s}))$ $\displaystyle=$ $\displaystyle~{}\delta_{\widehat{s}}+\delta_{\overline{\delta}}\beta_{c_{s}}-L_{\overline{x}}^{-\top}(\delta_{h}\beta_{s}-\delta_{\widetilde{h}}\widetilde{\beta}_{s}+\delta_{\epsilon_{s}})$ $\displaystyle=$ $\displaystyle~{}0.$ $$

Therefore, after Update$h$ finishes, we have

$$ $\displaystyle s=\widehat{s}+(H_{w,\overline{x}}^{\mathrm{new}})^{1/2}c_{s}\beta_{c_{s}}-L_{\overline{x}}^{-\top}(h\beta_{s}-\widetilde{h}\widetilde{\beta}_{s}+\epsilon_{s}).$ $$

So $x$ and $s$ are both updated correctly. This proves the correctness of Update$h$.

Step $\textsc{Update}\mathcal{W}$:
Define $\epsilon_{x}^{\mathrm{new}}:=\epsilon_{x}+\delta_{\epsilon_{x}}$, $\epsilon_{s}^{\mathrm{new}}:=\epsilon_{s}+\delta_{\epsilon_{s}}$.
Immediately after Algorithm [3](#alg3), Line [31](#alg3.l31), we have

$$ $\displaystyle~{}(\widehat{x}+(L_{\overline{x}}^{\mathrm{new}})^{-\top}(h\beta_{x}-\widetilde{h}\widetilde{\beta}_{x}+\epsilon_{x}^{\mathrm{new}}))-(\widehat{x}+L_{\overline{x}}^{-\top}(h\beta_{x}-\widetilde{h}\widetilde{\beta}_{x}+\epsilon_{x}))$ $\displaystyle=$ $\displaystyle~{}((L_{\overline{x}}^{\mathrm{new}})^{-\top}-L_{\overline{x}}^{-\top})(h\beta_{x}-\widetilde{h}\widetilde{\beta}_{x}+\epsilon_{x})+(L_{\overline{x}}^{\mathrm{new}})^{-\top}\delta_{\epsilon_{x}}$ $\displaystyle=$ $\displaystyle~{}0,$ $\displaystyle~{}(\widehat{s}+(H_{w,\overline{x}}^{\mathrm{new}})^{1/2}c_{s}\beta_{c_{s}}-(L_{\overline{x}}^{\mathrm{new}})^{-\top}(h\beta_{s}-\widetilde{h}\widetilde{\beta}_{s}+\epsilon_{s}^{\mathrm{new}}))$ $\displaystyle~{}-(\widehat{s}+(H_{w,\overline{x}}^{\mathrm{new}})^{1/2}c_{s}\beta_{c_{s}}-L_{\overline{x}}^{-\top}(h\beta_{s}-\widetilde{h}\widetilde{\beta}_{s}+\epsilon_{s}))$ $\displaystyle=$ $\displaystyle~{}(-(L_{\overline{x}}^{\mathrm{new}})^{-\top}+L_{\overline{x}}^{-\top})(h\beta_{s}-\widetilde{h}\widetilde{\beta}_{s}+\epsilon_{s})-(L_{\overline{x}}^{\mathrm{new}})^{-\top}\delta_{\epsilon_{s}}$ $\displaystyle=$ $\displaystyle~{}0.$ $$

Therefore, after Update$\mathcal{W}$ finishes, we have

$$ $\displaystyle x$ $\displaystyle=\widehat{x}+(L_{\overline{x}}^{\mathrm{new}})^{-\top}(h\beta_{x}-\widetilde{h}\widetilde{\beta}_{x}+\epsilon_{x}),$ $\displaystyle s$ $\displaystyle=\widehat{s}+(H_{w,\overline{x}}^{\mathrm{new}})^{1/2}c_{s}\beta_{c_{s}}-(L_{\overline{x}}^{\mathrm{new}})^{-\top}(h\beta_{s}-\widetilde{h}\widetilde{\beta}_{s}+\epsilon_{s}).$ $$

So $x$ and $s$ are both updated correctly. This proves the correctness of Update$\mathcal{W}$.
∎

###### Lemma 7.5 .

We bound the running time of ExactDS as following.

- (i)
ExactDS.Initialize (Algorithm [2](#alg2)) runs in $\widetilde{O}(n\tau^{\omega-1}+n\tau m+nm^{\omega-1})$ time.
- (ii)
ExactDS.Move (Algorithm [2](#alg2)) runs in $\widetilde{O}(m^{\omega})$ time.
- (iii)
ExactDS.Output (Algorithm [2](#alg2)) runs in $\widetilde{O}(n\tau m)$ time and correctly outputs $(x,s)$.
- (iv)
ExactDS.Query$x$ and ExactDS.Query$s$ (Algorithm [2](#alg2)) runs in $\widetilde{O}(\tau^{2}m)$ time and returns the correct answer.
- (v)
ExactDS.Update (Algorithm [2](#alg2)) runs in $\widetilde{O}((\tau^{2}m+\tau m^{2})(\|\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{s}}\|_{0}))$ time.
Furthermore, $h,\epsilon_{x},\epsilon_{s}$ change in $O(\tau(\|\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{s}}\|_{0}))$ coordinates,
$\widetilde{h}$ changes in $\widetilde{O}(\tau m(\|\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{s}}\|_{0}))$ coordinates,
and $H_{\overline{x}}^{1/2}\widehat{x},H_{\overline{x}}^{-1/2}\widehat{s},c_{s}$ change in $O(\|\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{s}}\|_{0})$ coordinates.

###### Proof.

- (i)
Computing $L_{\overline{x}}$ takes $\widetilde{O}(n\tau^{\omega-1})$ time by Lemma [4.5](#S4.Thmtheorem5).
Computing $h$ and $\widetilde{h}$ takes $\widetilde{O}(n\tau m)$ by Lemma [4.7](#S4.Thmtheorem7)[(i)](#S4.I2.i1).(^6^66Here we compute $\widetilde{h}$ by computing $\widetilde{h}_{*,i}=L_{\overline{x}}^{-1}(A_{i,*})^{\top}$ for $i\in[m]$ independently. Using fast rectangular matrix multiplication is possible to improve this running time and other similar places. We keep the current bounds for simplicity.)
Computing $\widetilde{u}$ and $u$ takes ${\cal T}_{\mathrm{mat}}(m,n,m)=\widetilde{O}(nm^{\omega-1})$ time.
All other operations are cheap.
- (ii)
Computing $\widetilde{u}^{-1}$ takes $\widetilde{O}(m^{\omega})$ time. All other operations take $O(m^{2})$ time.
- (iii)
Running time is by Lemma [4.7](#S4.Thmtheorem7)[(v)](#S4.I2.i5). Correctness is by Lemma [7.4](#S7.Thmtheorem4).
- (iv)
Running time is by Lemma [4.8](#S4.Thmtheorem8)[(ii)](#S4.I3.i2). Correctness is by Lemma [7.4](#S7.Thmtheorem4).
- (v)
Computing $\Delta_{L_{\overline{x}}}$ takes $\widetilde{O}(\tau^{2}\|\delta_{\overline{x}}\|_{0})$ time by Lemma [4.6](#S4.Thmtheorem6). It is easy to see that $\operatorname{nnz}(\Delta_{H_{w,\overline{x}}})=O(\|\delta_{\overline{x}}\|_{0})$ and $\operatorname{nnz}(\Delta_{L_{\overline{x}}})=\widetilde{O}(\tau^{2}\|\delta_{\overline{x}}\|_{0})$. It remains to analyze Update$h$ and Update$\mathcal{W}$.
For simplicity, we write $k=\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{s}}\|_{0}$ in this proof only.
•
Update$h$: Updating $\overline{\alpha}$ and computing $\delta_{\overline{\delta}_{\mu}}$ takes $O(k)$ time.
Also, $\|\delta_{\overline{\delta}_{\mu}}\|_{0}=O(k)$.
Computing $\delta_{h}$ takes $\widetilde{O}(\tau^{2}k)$ time by Lemma [4.7](#S4.Thmtheorem7)[(i)](#S4.I2.i1).
Also, $\delta_{h}$ is supported on $O(k)$ paths in the elimination tree, so $\|\delta_{h}\|_{0}=\widetilde{O}(\tau k)$.
Similarly we see that computing $\delta_{\widetilde{h}}$ take $\widetilde{O}(\tau^{2}mk)$ time and $\operatorname{nnz}(\delta_{\widetilde{h}})=\widetilde{O}(\tau mk)$.
Computing $\delta_{c_{s}}$ and $\delta_{\widehat{s}}$ takes $O(\tau^{2}k)$ time and $\|\delta_{c_{s}}\|_{0},\|\delta_{\widehat{s}}\|_{0}=O(k)$.
Computing $\delta_{\epsilon_{x}}$ and $\delta_{\epsilon_{s}}$ takes $O(\tau mk)$ time after computing $\delta_{h}$ and $\delta_{\widetilde{h}}$. Furthermore, $\|\delta_{\epsilon_{x}}\|_{0},\|\delta_{\epsilon_{s}}\|_{0}=O(\tau k)$.
Computing $\delta_{\widetilde{u}}$ takes ${\cal T}_{\mathrm{mat}}(m,\tau k,m)=\widetilde{O}(\tau m^{2}k)$ time.
Computing $\delta_{u}$ takes $\widetilde{O}(\tau mk)$ time.
•
Update$\mathcal{W}$:
To compute $\delta_{\epsilon_{x}}$ and $\delta_{\epsilon_{s}}$, we first compute $L_{\overline{x}}^{-\top}(h\beta_{x}-\widetilde{h}\widetilde{\beta}_{x}+\epsilon_{x})$ and $L_{\overline{x}}^{-\top}(h\beta_{s}-\widetilde{h}\widetilde{\beta}_{s}+\epsilon_{s})$, where $S\subseteq[n_{\mathrm{tot}}]$ is the row support of $\Delta_{L_{\overline{x}}}$, which can be decomposed into at most $O(\|\delta_{\overline{x}}\|_{0})$ paths.
This takes $\widetilde{O}(\tau^{2}m\|\delta_{\overline{x}}\|_{0})$ time by Lemma [4.8](#S4.Thmtheorem8)[(i)](#S4.I3.i1) (the extra $m$ factor is due to $\widetilde{h}$).
Combining everything we finish the proof of running time of ExactDS.Update.

∎

###### Theorem 7.3 .

###### Proof of Theorem 7.3 .

###### Lemma 7.4 .

###### Proof.

###### Lemma 7.5 .

###### Proof.

#### 7.3.2 ApproxDS

In this section we present the data structure ApproxDS.
Given BatchSketch, a data structure maintaining a sketch of the primal-dual pair $(x,s)\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$, ApproxDS maintains a sparsely-changing $\ell_{\infty}$-approximation of $(x,s)$. This data structure is a slight variation of ApproxDS in [].

Figure: Algorithm 4 The ApproxDS data structure used in Algorithm [1](#alg1).

Figure: Algorithm 5 ApproxDS Algorithm [4](#alg4) continued.

###### Theorem 7.6 .

Given parameters $\epsilon_{\mathrm{apx},x},\epsilon_{\mathrm{apx},s}\in(0,1),\delta_{\mathrm{apx}}\in(0,1)$, $\zeta_{x},\zeta_{s}\in\mathbb{R}$ such that

$$ $\displaystyle\|H_{w,\overline{x}^{(\ell)}}^{1/2}x^{(\ell)}-H_{w,\overline{x}^{(\ell)}}^{1/2}x^{(\ell+1)}\|_{2}\leq\zeta_{x},~{}~{}~{}\|H_{w,\overline{x}^{(\ell)}}^{-1/2}s^{(\ell)}-H_{w,\overline{x}^{(\ell)}}^{-1/2}s^{(\ell+1)}\|_{2}\leq\zeta_{s}$ $$

for all $\ell\in\{0,\ldots,q-1\}$,
data structure ApproxDS (Algorithm [4](#alg4) and Algorithm [5](#alg5)) supports the following operations:

- •
$\textsc{Initialize}(x,s\in\mathbb{R}^{n_{\mathrm{tot}}},h\in\mathbb{R}^{n_{\mathrm{tot}}},\widetilde{h}\in\mathbb{R}^{n_{\mathrm{tot}}\times m},\epsilon_{x},\epsilon_{s},H_{w,\overline{x}}^{1/2}\widehat{x},H_{w,\overline{x}}^{-1/2}\widehat{s},c_{s}\in\mathbb{R}^{n_{\mathrm{tot}}},\beta_{x},\beta_{s},\beta_{c_{s}}\in\mathbb{R},\widetilde{\beta}_{x},\widetilde{\beta}_{s}\in\mathbb{R}^{m},q\in\mathbb{N},\textsc{ExactDS*}~{}\mathsf{exact},\epsilon_{\mathrm{apx},x},\epsilon_{\mathrm{apx},s},\delta_{\mathrm{apx}}\in\mathbb{R})$: Initialize the data structure in $\widetilde{O}(n\tau^{\omega-1}+n\tau m)$ time.
- •
$\textsc{MoveAndQuery}(\beta_{x},\beta_{s},\beta_{c_{s}}\in\mathbb{R}$, $\widetilde{\beta}_{x},\widetilde{\beta}_{s}\in\mathbb{R}^{m})$:
Update values of $\beta_{x},\beta_{s},\beta_{c_{s}},\widetilde{\beta}_{x},\widetilde{\beta}_{s}$ by calling $\textsc{BatchSketch}.\textsc{Move}$.
This effectively moves $(x^{(\ell)},s^{(\ell)})$ to $(x^{(\ell+1)},s^{(\ell+1)})$ while keeping $\overline{x}^{(\ell)}$ unchanged.
Then return two sets $L_{x}^{(\ell)},L_{s}^{(\ell)}\subset[n]$ where
$\displaystyle L_{x}^{(\ell)}$
$\displaystyle\supseteq\{i\in[n]:\|H_{w,\overline{x}^{(\ell)}}^{1/2}x^{(\ell)}_{i}-H_{w,\overline{x}^{(\ell)}}^{1/2}x^{(\ell+1)}_{i}\|_{2}\geq\epsilon_{\mathrm{apx},x}\},$
$\displaystyle L_{s}^{(\ell)}$
$\displaystyle\supseteq\{i\in[n]:\|H_{w,\overline{x}^{(\ell)}}^{-1/2}s^{(\ell)}_{i}-H_{w,\overline{x}^{(\ell)}}^{-1/2}s^{(\ell+1)}_{i}\|_{2}\geq\epsilon_{\mathrm{apx},s}\},$
satisfying
$\displaystyle\sum_{0\leq\ell\leq q-1}|L_{x}^{(\ell)}|=\widetilde{O}(\epsilon_{\mathrm{apx},x}^{-2}\zeta_{x}^{2}q^{2}),$
$\displaystyle\sum_{0\leq\ell\leq q-1}|L_{s}^{(\ell)}|=\widetilde{O}(\epsilon_{\mathrm{apx},s}^{-2}\zeta_{s}^{2}q^{2}).$
For every query, with probability at least $1-\delta_{\mathrm{apx}}/q$, the return values are correct.
Furthermore, total time cost over all queries is at most
$\displaystyle\widetilde{O}\left((\epsilon_{\mathrm{apx},x}^{-2}\zeta_{x}^{2}+\epsilon_{\mathrm{apx},s}^{-2}\zeta_{s}^{2})q^{2}\tau^{2}m\right).$
- •
$\textsc{Update}(\delta_{\overline{x}}\in\mathbb{R}^{n_{\mathrm{tot}}},\delta_{h}\in\mathbb{R}^{n_{\mathrm{tot}}},\delta_{\widetilde{h}}\in\mathbb{R}^{n_{\mathrm{tot}}\times m},\delta_{\epsilon_{x}},\delta_{\epsilon_{s}},\delta_{H_{w,\overline{x}}^{1/2}\widehat{x}},\delta_{H_{w,\overline{x}}^{-1/2}\widehat{s}},\delta_{c_{s}}\in\mathbb{R}^{n_{\mathrm{tot}}})$: Update sketches of $H_{w,\overline{x}^{(\ell)}}^{1/2}x^{(\ell+1)}$ and $H_{w,\overline{x}^{(\ell)}}^{-1/2}s^{(\ell+1)}$ by calling BatchSketch.Update.
This effectively moves $\overline{x}^{(\ell)}$ to $\overline{x}^{(\ell+1)}$ while keeping $(x^{(\ell+1)},s^{(\ell+1)})$ unchanged. Then advance timestamp $\ell$.
Each update costs
$\displaystyle\widetilde{O}(\tau^{2}(\|\delta_{\overline{x}}\|_{0}+\|\delta_{h}\|_{0}+\|\delta_{\widetilde{h}}\|_{0}+\|\delta_{\epsilon_{x}}\|_{0}+\|\delta_{\epsilon_{s}}\|_{0})+\|\delta_{H_{w,\overline{x}}^{1/2}\widehat{x}}\|_{0}+\|\delta_{H_{w,\overline{x}}^{-1/2}\widehat{s}}\|_{0}+\|\delta_{c_{s}}\|_{0})$
time.

###### Proof.

The proof is essentially the same as proof of [, Theorem 4.18].
For the running time claims, we plug in Theorem [7.8](#S7.Thmtheorem8) when necessary.
∎

###### Theorem 7.6 .

###### Proof.

#### 7.3.3 BatchSketch

In this section we present the data structure BatchSketch. It maintains a sketch of $H_{\overline{x}}^{1/2}x$ and $H_{\overline{x}}^{-1/2}s$. It is a variation of BatchSketch in [].

We recall the following definition from [].

###### Definition 7.7 (Partition tree) .

A partition tree $(\mathcal{S},\chi)$ of $\mathbb{R}^{n}$ is a constant degree rooted tree $\mathcal{S}=(V,E)$ and a labeling of the vertices $\chi:V\to 2^{[n]}$, such that

- •
$\chi(\text{root})=[n]$;
- •
if $v$ is a leaf of $\mathcal{S}$, then $|\chi(v)|=1$;
- •
for any non-leaf node $v\in V$, the set $\{\chi(c):\text{$c$ is a child of $v$}\}$ is a partition of $\chi(v)$.

Figure: Algorithm 6 The BatchSketch data structure used by Algorithm [4](#alg4) and [5](#alg5).

Figure: Algorithm 7 BatchSketch Algorithm [6](#alg6) continued.

Figure: Algorithm 8 BatchSketch Algorithm [6](#alg6), [7](#alg7) continued.

###### Theorem 7.8 .

Data structure BatchSketch (Algorithm [6](#alg6), [8](#alg8)) supports the following operations:

- •
$\textsc{Initialize}(\overline{x}\in\mathbb{R}^{n_{\mathrm{tot}}},h\in\mathbb{R}^{n_{\mathrm{tot}}},\widetilde{h}\in\mathbb{R}^{n_{\mathrm{tot}}\times m},\epsilon_{x},\epsilon_{s},H_{w,\overline{x}}^{1/2}\widehat{x},H_{w,\overline{x}}^{-1/2}\widehat{s},c_{s}\in\mathbb{R}^{n_{\mathrm{tot}}},\beta_{x},\beta_{s},\beta_{c_{s}}\in\mathbb{R},\widetilde{\beta}_{x},\widetilde{\beta}_{s}\in\mathbb{R}^{m},\delta_{\mathrm{apx}}\in\mathbb{R})$: Initialize the data structure in $\widetilde{O}(n\tau^{\omega-1}+n\tau m)$ time.
- •
$\textsc{Move}(\beta_{x},\beta_{s},\beta_{c_{s}}\in\mathbb{R},\widetilde{\beta}_{x},\widetilde{\beta}_{s}\in\mathbb{R}^{m})$: Update values of $\beta_{x},\beta_{s},\beta_{c_{s}},\widetilde{\beta}_{x},\widetilde{\beta}_{s}$ in $O(m)$ time. This effectively moves $(x^{(\ell)},s^{(\ell)})$ to $(x^{(\ell+1)},s^{(\ell+1)})$ while keeping $\overline{x}^{(\ell)}$ unchanged.
- •
$\textsc{Update}(\delta_{\overline{x}}\in\mathbb{R}^{n_{\mathrm{tot}}},\delta_{h}\in\mathbb{R}^{n_{\mathrm{tot}}},\delta_{\widetilde{h}}\in\mathbb{R}^{n_{\mathrm{tot}}\times m},\delta_{\epsilon_{x}},\delta_{\epsilon_{s}},\delta_{H_{w,\overline{x}}^{1/2}\widehat{x}},\delta_{H_{w,\overline{x}}^{-1/2}\widehat{s}},\delta_{c_{s}}\in\mathbb{R}^{n_{\mathrm{tot}}})$:
Update sketches of $H_{w,\overline{x}^{(\ell)}}^{1/2}x^{(\ell+1)}$ and $H_{w,\overline{x}^{(\ell)}}^{-1/2}s^{(\ell+1)}$.
This effectively moves $\overline{x}^{(\ell)}$ to $\overline{x}^{(\ell+1)}$ while keeping $(x^{(\ell+1)},s^{(\ell+1)})$ unchanged.
Then advance timestamp $\ell$.
Each update costs
$\displaystyle\widetilde{O}(\tau^{2}(\|\delta_{\overline{x}}\|_{0}+\|\delta_{h}\|_{0}+\|\delta_{\widetilde{h}}\|_{0}+\|\delta_{\epsilon_{x}}\|_{0}+\|\delta_{\epsilon_{s}}\|_{0})+\|\delta_{H_{w,\overline{x}}^{1/2}\widehat{x}}\|_{0}+\|\delta_{H_{w,\overline{x}}^{-1/2}\widehat{s}}\|_{0}+\|\delta_{c_{s}}\|_{0})$
time.
- •
$\textsc{Query}x(\ell^{\prime}\in\mathbb{N},\epsilon\in\mathbb{R})$:
Given timestamp $\ell^{\prime}$, return a set $S\subseteq[n]$ where
$\displaystyle S$
$\displaystyle\supseteq\{i\in[n]:\|H_{w,\overline{x}^{(\ell^{\prime})}}^{1/2}x^{(\ell^{\prime})}_{i}-H_{w,\overline{x}^{(\ell)}}^{1/2}x^{(\ell+1)}_{i}\|_{2}\geq\epsilon\},$
and
$\displaystyle|S|$
$\displaystyle=O(\epsilon^{-2}(\ell-\ell^{\prime}+1)\sum_{\ell^{\prime}\leq t\leq\ell}\|H_{w,\overline{x}^{(t)}}^{1/2}x^{(t)}-H_{w,\overline{x}^{(t)}}^{1/2}x^{(t+1)}\|_{2}^{2}+\sum_{\ell^{\prime}\leq t\leq\ell-1}\|\overline{x}^{(t)}-\overline{x}^{(t+1)}\|_{2,0})$
where $\ell$ is the current timestamp.
For every query, with probability at least $1-\delta$, the return values are correct, and costs at most
$\displaystyle\widetilde{O}(\tau^{2}\cdot(\epsilon^{-2}(\ell-\ell^{\prime}+1)\sum_{\ell^{\prime}\leq t\leq\ell}\|H_{\overline{x}^{(t)}}^{1/2}x^{(t)}-H_{\overline{x}^{(t)}}^{1/2}x^{(t+1)}\|_{2}^{2}+\sum_{\ell^{\prime}\leq t\leq\ell-1}\|\overline{x}^{(t)}-\overline{x}^{(t+1)}\|_{2,0}))$
running time.
- •
$\textsc{Query}s(\ell^{\prime}\in\mathbb{N},\epsilon\in\mathbb{R})$:
Given timestamp $\ell^{\prime}$, return a set $S\subseteq[n]$ where
$\displaystyle S$
$\displaystyle\supseteq\{i\in[n]:\|H_{w,\overline{x}^{(\ell^{\prime})}}^{-1/2}s^{(\ell^{\prime})}_{i}-H_{w,\overline{x}^{(\ell)}}^{-1/2}s^{(\ell+1)}_{i}\|_{2}\geq\epsilon\}$
and
$\displaystyle|S|$
$\displaystyle=O(\epsilon^{-2}(\ell-\ell^{\prime}+1)\sum_{\ell^{\prime}\leq t\leq\ell}\|H_{w,\overline{x}^{(t)}}^{-1/2}s^{(t)}-H_{w,\overline{x}^{(t)}}^{-1/2}s^{(t+1)}\|_{2}^{2}+\sum_{\ell^{\prime}\leq t\leq\ell-1}\|\overline{x}^{(t)}-\overline{x}^{(t+1)}\|_{2,0})$
where $\ell$ is the current timestamp.
For every query, with probability at least $1-\delta$, the return values are correct, and costs at most
$\displaystyle\widetilde{O}(\tau^{2}\cdot(\epsilon^{-2}(\ell-\ell^{\prime}+1)\sum_{\ell^{\prime}\leq t\leq\ell}\|H_{\overline{x}^{(t)}}^{1/2}s^{(t)}-H_{\overline{x}^{(t)}}^{1/2}x^{(t+1)}\|_{2}^{2}+\sum_{\ell^{\prime}\leq t\leq\ell-1}\|\overline{x}^{(t)}-\overline{x}^{(t+1)}\|_{2,0}))$
running time.

###### Proof.

The proof is essentially the same as proof of [, Theorem 4.21].
For the running time claims, we plug in Lemma [7.10](#S7.Thmtheorem10) and [7.12](#S7.Thmtheorem12) when necessary.
∎

###### Remark 7.9 (Snapshot) .

As in previous works, we use persistent data structures (e.g., []) to keep a snapshot of the data structure after every update. This allows us to support query to historical data.
This incurs an $O(\log n_{\mathrm{tot}})=\widetilde{O}(1)$ multiplicative factor in all running times, which we ignore in our analysis.

###### Definition 7.7 (Partition tree) .

###### Theorem 7.8 .

###### Proof.

###### Remark 7.9 (Snapshot) .

#### 7.3.4 VectorSketch

VectorSketch is a data structure used to maintain sketches of sparsely-changing vectors.
It is a direct application of segment trees. For completeness, we include code (Algorithm [9](#alg9)) from [, Algorithm 10].

Figure: Algorithm 9 [, Algorithm 10]. Used in Algorithm [6](#alg6), [7](#alg7), [8](#alg8).

###### Lemma 7.10 ( [ 28 , Lemma 4.23] ) .

Given a partition tree $(\mathcal{S},\chi)$ of $\mathbb{R}^{n}$, and a JL sketching matrix $\Phi\in\mathbb{R}^{r\times n_{\mathrm{tot}}}$, the data structure VectorSketch (Algorithm [9](#alg9)) maintains $\Phi_{\chi(v)}x_{\chi(v)}$ for all nodes $v$ in the partition tree implicitly through the following functions:

- •
Initialize$(\mathcal{S},\chi,\Phi)$: Initializes the data structure in $O(rn_{\mathrm{tot}})$ time.
- •
Update$(\delta_{x}\in\mathbb{R}^{n_{\mathrm{tot}}})$: Maintains the data structure for $x\leftarrow x+\delta_{x}$ in $O(r\|\delta_{x}\|_{0}\log n)$ time.
- •
Query$(v\in V(\mathcal{S}))$: Outputs $\Phi_{\chi(v)}x_{\chi(v)}$ in $O(r\log n)$ time.

###### Lemma 7.10 ( [ 28 , Lemma 4.23] ) .

#### 7.3.5 BalancedSketch

In this section, we present data structure BalancedSketch. It is a data structure for maintaining a sketch of a vector of form $\mathcal{W}^{\top}h$, where $\mathcal{W}=L_{\overline{x}}^{-1}H_{w,\overline{x}}^{1/2}$ and $h\in\mathbb{R}^{n_{\mathrm{tot}}}$ is a sparsely-changing vector. This is a variation of BlockBalancedSketch in [].

We use the following construction of a partition tree.

###### Definition 7.11 (Construction of Partition Tree) .

We fix an ordering $\pi$ of $[n]$ using the heavy-light decomposition (Lemma [4.10](#S4.Thmtheorem10)).
Let $\mathcal{S}$ be a complete binary tree with leaf set $[n]$ and ordering $\pi$.
Let $\chi$ map a node to the set of leaves in its subtree.
Then $(\mathcal{S},\chi)$ is a valid partition tree.

Figure: Algorithm 10 The BalancedSketch data structure is used in Algorithm [6](#alg6), [7](#alg7), [8](#alg8).

Figure: Algorithm 11 BalancedSketch Algorithm [10](#alg10) continued. This is used in Algorithm [6](#alg6), [7](#alg7), [8](#alg8).

Figure: Algorithm 12 BalancedSketch Algorithm [10](#alg10), [11](#alg11) continued. This is used in Algorithm [6](#alg6), [7](#alg7), [8](#alg8).

###### Lemma 7.12 .

Given an elimination tree $\mathcal{T}$ with height $\eta$, a JL matrix $\Phi\in\mathbb{R}^{r\times n_{\mathrm{tot}}}$, and a partition tree $(\mathcal{S},\chi)$ constructed as in Definition [7.11](#S7.Thmtheorem11)
with height $\widetilde{O}(1)$, the data structure BalancedSketch (Algorithm [10](#alg10), [11](#alg11), [12](#alg12)), maintains $\Phi_{\chi(v)}(\mathcal{W}^{\top}h)_{\chi(v)}$ for each $v\in V(\mathcal{S})$ through the following operations

- •
Initialize$((\mathcal{S},\chi):\text{partition tree},\Phi\in\mathbb{R}^{n_{\mathrm{tot}}},\overline{x}\in\mathbb{R}^{n_{\mathrm{tot}}},h\in\mathbb{R}^{n_{\mathrm{tot}}\times k})$: Initializes the data structure in $\widetilde{O}(r(n\tau^{\omega-1}+n\tau k))$ time.
- •
Update$(\delta_{\overline{x}}\in\mathbb{R}^{n_{\mathrm{tot}}},\delta_{h}\in\mathbb{R}^{n_{\mathrm{tot}}\times k})$: Updates all sketches in $\mathcal{S}$ implicitly to reflect $(\mathcal{W},h)$ updating to $(\mathcal{W}^{\mathrm{new}},h^{\mathrm{new}})$ in $\widetilde{O}(r\tau^{2}k)$ time.
- •
Query$(v\in\mathcal{S})$: Outputs $\Phi_{\chi(v)}(\mathcal{W}^{\top}h)_{\chi(v)}$ in $\widetilde{O}(r\tau^{2}k)$ time.

###### Proof.

The proof is almost same as the proof of [, Lemma 4.24].
(In fact, our $\mathcal{W}$ is simpler than the one used in [].)

For Initialize running time, we note that computing $Z_{v}$ for all $v\in\mathcal{B}$ takes $\widetilde{O}(rn\tau^{\omega-1})$ time by [, Lemma 8.3].
Because $Z_{v}$ is supported on the path from $v$ to the root in $\mathcal{T}$, we know that $\operatorname{nnz}(Z)=O(rn\tau)$.
Therefore computing $y_{v}^{\triangledown}$ for all $v\in\mathcal{B}$ takes $\widetilde{O}(rn\tau k)$ time.

Remaining claims follow from combining proof of [, Lemma 4.24] and [, Lemma 8.3].
∎

###### Definition 7.11 (Construction of Partition Tree) .

###### Lemma 7.12 .

###### Proof.

### 7.4 Analysis of CentralPathMaintenance

###### Lemma 7.13 (Correctness of CentralPathMaintenance ) .

Algorithm [1](#alg1) implicitly maintains the primal-dual solution pair $(x,s)$ via representation Eq. ([8](#S7.E8))([9](#S7.E9)).
It also explicitly maintains $(\overline{x},\overline{s})\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$ such that $\|\overline{x}_{i}-x_{i}\|_{\overline{x}_{i}}\leq\overline{\epsilon}$ and $\|\overline{s}_{i}-s_{i}\|_{\overline{x}_{i}}^{*}\leq t\overline{\epsilon}w_{i}$ for all $i\in[n]$ with probability at least $0.9$.

###### Proof.

We correctly maintain the implicit representation because of correctness of $\mathsf{exact}$.Update (Theorem [7.3](#S7.Thmtheorem3)).

We show that $\|\overline{x}_{i}-x_{i}\|_{\overline{x}_{i}}\leq\overline{\epsilon}$ and $\|\overline{s}_{i}-s_{i}\|_{\overline{x}_{i}}^{*}\leq t\overline{\epsilon}w_{i}$ for all $i\in[n]$ (c.f. Algorithm [20](#alg20), Line [16](#alg20.l16)).
$\mathsf{approx}$ maintains an $\ell_{\infty}$ approximation of $H_{w,\overline{x}}^{1/2}x$.
For $\ell\leq q$, we have

$$ $\displaystyle\|H_{w,\overline{x}}^{1/2}x^{(\ell+1)}-H_{w,\overline{x}}^{1/2}x^{(\ell)}\|_{2}=\|\delta_{x}\|_{w,\overline{x}}\leq\frac{9}{8}\alpha\leq\zeta_{x}$ $$

where the first step from definition of $\|\cdot\|_{w,\overline{x}}$, the second step follows from Lemma [9.11](#S9.Thmtheorem11), the third step follows from definition of $\zeta_{x}$.

By Theorem [7.6](#S7.Thmtheorem6), with probability at least $1-\delta_{\mathrm{apx}}$, $\mathsf{approx}$ correctly maintains $\overline{x}$ such that $\|H_{w,\overline{x}}^{1/2}\overline{x}-H_{w,\overline{x}}^{1/2}x\|_{\infty}\leq\epsilon_{\mathrm{apx},x}\leq\overline{\epsilon}$.
Then

$$ $\displaystyle\|\overline{x}_{i}-x_{i}\|_{\overline{x}_{i}}\leq w_{i}^{-1/2}\|H_{w,\overline{x}}^{1/2}\overline{x}-H_{w,\overline{x}}^{1/2}x\|_{\infty}\leq w_{i}^{-1/2}\overline{\epsilon}\leq\overline{\epsilon}.$ $$

Note that the last step is loose by a factor of $w_{i}^{1/2}$.
When $w_{i}$s are large, we could improve running time by using a tighter choice of $\epsilon_{\mathrm{apx},x}$, as did in []. Here we use a loose bound for simplicity of presentation. Same remark applies to $s$.

The proof for $s$ is similar.
We have

$$ $\displaystyle\|H_{w,\overline{x}}^{-1/2}\delta_{s}\|_{2}=\|\delta_{s}\|_{w,\overline{x}}^{*}\leq\frac{17}{8}\alpha\cdot t\leq\zeta_{s}$ $$

and

$$ $\displaystyle\|\overline{s}_{i}-s_{i}\|_{\overline{x}_{i}}^{*}\leq w_{i}^{1/2}\|H_{w,\overline{x}}^{-1/2}\overline{s}-H_{w,\overline{x}}^{-1/2}s\|_{\infty}\leq w_{i}^{1/2}\epsilon_{\mathrm{apx},s}\leq\overline{\epsilon}\cdot\overline{t}\cdot w_{i}.$ $$

∎

###### Lemma 7.14 .

We bound the running time of CentralPathMaintenance as following.

- •
CentralPathMaintenance.Initialize takes $\widetilde{O}(n\tau^{\omega-1}+n\tau m+nm^{\omega-1})$ time.
- •
If CentralPathMaintenance.MultiplyAndMove is called $N$ times, then it has total running time
$\displaystyle\widetilde{O}((Nn^{-1/2}+\log(t_{\max}/t_{\min}))\cdot n(\tau^{2}m+\tau m^{2})^{1/2}(\tau^{\omega-1}+\tau m+m^{\omega-1})^{1/2}).$
- •
CentralPathMaintenance.Output takes $\widetilde{O}(n\tau m)$ time.

###### Proof.

Initialize part:
By Theorem [7.3](#S7.Thmtheorem3) and [7.6](#S7.Thmtheorem6).

Output part:
By Theorem [7.3](#S7.Thmtheorem3).

MultiplyAndMove part:
Between two restarts, the total size of $|L_{x}|$ returned by $\mathsf{approx}$.Query is bounded by $\widetilde{O}(q^{2}\zeta_{x}^{2}/\epsilon_{\mathrm{apx},x}^{2})$ by Theorem [7.6](#S7.Thmtheorem6).
By plugging in $\zeta_{x}=2\alpha$, $\epsilon_{\mathrm{apx},x}=\overline{\epsilon}$, we have
$\sum_{\ell\in[q]}|L_{x}^{(\ell)}|=\widetilde{O}(q^{2}).$
Similarly, for $s$ we have
$\sum_{\ell\in[q]}|L_{s}^{(\ell)}|=\widetilde{O}(q^{2}).$

Update time:
By Theorem [7.3](#S7.Thmtheorem3) and [7.6](#S7.Thmtheorem6), in a sequence of $q$ updates, total cost for update is $\widetilde{O}(q^{2}(\tau^{2}m+\tau m^{2}))$.
So the amortized update cost per iteration is $\widetilde{O}(q(\tau^{2}m+\tau m^{2}))$.
The total update cost is

$$ $\displaystyle\text{number of iterations}\cdot\text{time per iteration}=\widetilde{O}(Nq(\tau^{2}m+\tau m^{2})).$ $$

Init/restart time:
We restart the data structure whenever $k>q$ or $|\overline{t}-t|>\overline{t}\epsilon_{t}$, so there are $O(N/q+\log(t_{\max}/t_{\min})\epsilon_{t}^{-1})$ restarts in total.
By Theorem [7.3](#S7.Thmtheorem3) and [7.6](#S7.Thmtheorem6), time cost per restart is $\widetilde{O}(n(\tau^{\omega-1}+\tau m+m^{\omega-1}))$.
So the total initialization time is

$$ $\displaystyle\text{number of restarts}\cdot\text{time per restart}=\widetilde{O}((N/q+\log(t_{\max}/t_{\min})\epsilon_{t}^{-1})\cdot n(\tau^{\omega-1}+\tau m+m^{\omega-1})).$ $$

Combine everything:
Overall running time is

$$ $\displaystyle\widetilde{O}(Nq(\tau^{2}m+\tau m^{2})+(N/q+\log(t_{\max}/t_{\min})\epsilon_{t}^{-1})\cdot n(\tau^{\omega-1}+\tau m+m^{\omega-1})).$ $$

Taking $\epsilon_{t}=\frac{1}{2}\overline{\epsilon}$, the optimal choice for $q$ is

$$ $\displaystyle q=n^{1/2}(\tau^{2}m+\tau m^{2})^{-1/2}(\tau^{\omega-1}+\tau m+m^{\omega-1})^{1/2},$ $$

achieving overall running time

$$ $\displaystyle\widetilde{O}((Nn^{-1/2}+\log(t_{\max}/t_{\min}))\cdot n(\tau^{2}m+\tau m^{2})^{1/2}(\tau^{\omega-1}+\tau m+m^{\omega-1})^{1/2}).$ $$

∎

###### Proof of Theorem 7.2 .

Combining Lemma [7.13](#S7.Thmtheorem13) and [7.14](#S7.Thmtheorem14).
∎

###### Lemma 7.13 (Correctness of CentralPathMaintenance ) .

###### Proof.

###### Lemma 7.14 .

###### Proof.

###### Proof of Theorem 7.2 .

### 7.5 Proof of Main Statement

###### Proof of Theorem 7.1 .

Use CentralPathMaintenance (Algorithm [1](#alg1)) as the maintenance data structure in Algorithm [20](#alg20).
Combining Theorem [7.2](#S7.Thmtheorem2) and Theorem [9.1](#S9.Thmtheorem1) finishes the proof.
∎

###### Proof of Theorem 7.1 .

## 8 Algorithm for Low-Rank QP

In this section we present a nearly-linear time algorithm for solving low-rank QP with small number of linear constraints. We briefly describe the outline of this section.

- •
In Section [8.1](#S8.SS1), we present the main statement of Section [8](#S8).
- •
In Section [8.2](#S8.SS2), we present the main data structure CentralPathMaintenance.
- •
In Section [8.3](#S8.SS3), we present several data structures used in CentralPathMaintenance, including ExactDS (Section [8.3.1](#S8.SS3.SSS1)), ApproxDS (Section [8.3.2](#S8.SS3.SSS2)), BatchSketch (Section [8.3.3](#S8.SS3.SSS3)).
- •
In Section [8.4](#S8.SS4), we prove correctness and running time of CentralPathMaintenance data structure.
- •
In Section [8.5](#S8.SS5), we prove the main result (Theorem [8.1](#S8.Thmtheorem1)).

### 8.1 Main Statement

We consider programs of the form ([16](#S9.E16)), i.e.,

$$ $\displaystyle\min_{x\in\mathbb{R}^{n}}~{}$ $\displaystyle~{}\frac{1}{2}x^{\top}Qx+c^{\top}x$ $\displaystyle\mathrm{s.t.~{}}$ $\displaystyle~{}Ax=b$ $\displaystyle~{}x_{i}\in\mathcal{K}_{i}\qquad\forall i\in[n]$ $$

where $Q\in\mathcal{S}^{n_{\mathrm{tot}}}$, $c\in\mathbb{R}^{n_{\mathrm{tot}}}$,
$A\in\mathbb{R}^{m\times n_{\mathrm{tot}}}$, $b\in\mathbb{R}^{m}$, $\mathcal{K}_{i}\subset\mathbb{R}^{n_{i}}$ is a convex set.
For simplicity, we assume that $n_{i}=O(1)$ for all $i\in[n]$.

###### Theorem 8.1 .

Consider the convex program ([16](#S9.E16)).
Let $\phi_{i}:\mathcal{K}_{i}\to\mathbb{R}$ be a $\nu_{i}$-self-concordant barrier for all $i\in[n]$.
Suppose the program satisfies the following properties:

- •
Inner radius $r$: There exists $z\in\mathbb{R}^{n_{\mathrm{tot}}}$ such that $Az=b$ and $B(z,r)\in\mathcal{K}$.
- •
Outer radius $R$: $\mathcal{K}\subseteq B(0,R)$ where $0\in\mathbb{R}^{n_{\mathrm{tot}}}$.
- •
Lipschitz constant $L$: $\|Q\|_{2\to 2}\leq L$, $\|c\|_{2}\leq L$.
- •
Low rank: We are given a factorization $Q=UV^{\top}$ where $U,V\in\mathbb{R}^{n_{\mathrm{tot}}\times k}$.

Let $(w_{i})_{i\in[n]}\in\mathbb{R}_{\geq 1}^{n}$ and $\kappa=\sum_{i\in[n]}w_{i}\nu_{i}$.
Given any $0<\epsilon\leq\frac{1}{2}$, we can find an approximate solution $x\in\mathcal{K}$ satisfiying

$$ $\displaystyle\frac{1}{2}x^{\top}Qx+c^{\top}x$ $\displaystyle\leq\min_{Ax=b,x\in\mathcal{K}}\left(\frac{1}{2}x^{\top}Qx+c^{\top}x\right)+\epsilon LR(R+1),$ $\displaystyle\|Ax-b\|_{1}$ $\displaystyle\leq 3\epsilon(R\|A\|_{1}+\|b\|_{1}),$ $$

in expected time

$$ $\displaystyle\widetilde{O}((\sqrt{\kappa}n^{-1/2}+\log(R/(r\epsilon)))\cdot n(k+m)^{(\omega+1)/2}).$ $$

When $\max_{i\in[n]}\nu_{i}=\widetilde{O}(1)$, $w_{i}=1$, the running time simplifies to

$$ $\displaystyle\widetilde{O}(n(k+m)^{(\omega+1)/2})\log(R/(r\epsilon))).$ $$

###### Theorem 8.1 .

### 8.2 Algorithm Structure and Central Path Maintenance

Similar to the low-rank case, our algorithm is based on the robust Interior Point Method. Details of the robust IPM will be given in Section [9](#S9).
During the algorithm, we maintain a primal-dual solution pair $(x,s)\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$ on the robust central path. In addition, we maintain a sparsely-changing approximation $(\overline{x},\overline{s})\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$ to $(x,s)$.
In each iteration, we implicitly perform update

$$ $\displaystyle x$ $\displaystyle\leftarrow x+\overline{t}B_{w,\overline{x},\overline{t}}^{-1/2}(I-P_{w,\overline{x},\overline{t}})B_{w,\overline{x},\overline{t}}^{-1/2}\delta_{\mu}$ $\displaystyle s$ $\displaystyle\leftarrow s+\overline{t}\delta_{\mu}-\overline{t}^{2}H_{w,\overline{x}}B_{w,\overline{x},\overline{t}}^{-1/2}(I-P_{w,\overline{x},\overline{t}})B_{w,\overline{x},\overline{t}}^{-1/2}\delta_{\mu}$ $$

where

$$ $\displaystyle H_{w,\overline{x}}$ $\displaystyle=\nabla^{2}\phi_{w}(\overline{x})$ (see Eq. ( 24 )) $\displaystyle B_{w,\overline{x},\overline{t}}$ $\displaystyle=Q+\overline{t}H_{w,\overline{x}}$ (see Eq. ( 25 )) $\displaystyle P_{w,\overline{x},\overline{t}}$ $\displaystyle=B_{w,\overline{x},\overline{t}}^{-1/2}A^{\top}(AB_{w,\overline{x},\overline{t}}^{-1}A^{\top})^{-1}AB_{w,\overline{x},\overline{t}}^{-1/2}$ (see Eq. ( 26 )) $$

and explicitly maintain $(\overline{x},\overline{s})$ such that they remain close to $(x,s)$ in $\ell_{\infty}$-distance.

This task is handled by the CentralPathMaintenance data structure, which is our main data structure. The robust IPM algorithm (Algorithm [19](#alg19), [20](#alg20)) directly calls it in every iteration.

The CentralPathMaintenance data structure (Algorithm [13](#alg13)) has two main sub data structures, ExactDS (Algorithm [14](#alg14), [15](#alg15)) and ApproxDS (Algorithm [16](#alg16)).
ExactDS is used to maintain $(x,s)$, and ApproxDS is used to maintain $(\overline{x},\overline{s})$.

Figure: Algorithm 13

###### Theorem 8.2 .

Data structure CentralPathMaintenance (Algorithm [13](#alg13)) implicitly maintains the central path primal-dual solution pair $(x,s)\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$ and explicitly maintains its approximation $(\overline{x},\overline{s})\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$ using the following functions:

- •
Initialize$(x\in\mathbb{R}^{n_{\mathrm{tot}}},s\in\mathbb{R}^{n_{\mathrm{tot}}},t_{0}\in\mathbb{R}_{>0},\epsilon\in(0,1))$: Initializes the data structure with initial primal-dual solution pair $(x,s)\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$, initial central path timestamp $t_{0}\in\mathbb{R}_{>0}$ in $\widetilde{O}(n(k^{\omega-1}+m^{\omega-1}))$ time.
- •
MultiplyAndMove$(t\in\mathbb{R}_{>0})$: It implicitly maintains
$\displaystyle x$
$\displaystyle\leftarrow x+\overline{t}B_{w,\overline{x},\overline{t}}^{-1/2}(I-P_{w,\overline{x},\overline{t}})B_{w,\overline{x},\overline{t}}^{-1/2}\delta_{\mu}(\overline{x},\overline{s},\overline{t})$
$\displaystyle s$
$\displaystyle\leftarrow s+\overline{t}\delta_{\mu}-\overline{t}^{2}H_{w,\overline{x}}B_{w,\overline{x},\overline{t}}^{-1/2}(I-P_{w,\overline{x},\overline{t}})B_{w,\overline{x},\overline{t}}^{-1/2}\delta_{\mu}(\overline{x},\overline{s},\overline{t})$
where $H_{w,\overline{x}}$, $B_{w,\overline{x},\overline{t}}$, $P_{w,\overline{x},\overline{t}}$ are defined in Eq. ([24](#S9.E24))([25](#S9.E25))([26](#S9.E26)) respectively, and $\overline{t}$ is some timestamp satisfying $|\overline{t}-t|\leq\epsilon_{t}\cdot\overline{t}$.
It also explicitly maintains $(\overline{x},\overline{s})\in\mathbb{R}^{n_{\mathrm{tot}}\times n_{\mathrm{tot}}}$ such that $\|\overline{x}_{i}-x_{i}\|_{\overline{x}_{i}}\leq\overline{\epsilon}$ and $\|\overline{s}_{i}-s_{i}\|^{*}_{\overline{x}_{i}}\leq t\overline{\epsilon}w_{i}$ for all $i\in[n]$ with probability at least $0.9$.
Assuming the function is called at most $N$ times and $t$ decreases from $t_{\max}$ to $t_{\min}$, the total running time is
$\displaystyle\widetilde{O}((Nn^{-1/2}+\log(t_{\max}/t_{\min}))\cdot n(k^{(\omega+1)/2}+m^{(\omega+1)/2})).$
- •
Output: Computes $(x,s)\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$ exactly and outputs them in
$\widetilde{O}(n(k+m))$
time.

###### Theorem 8.2 .

### 8.3 Data Structures Used in CentralPathMaintenance

In this section we present several data structures used in CentralPathMaintenance, including:

- •
ExactDS (Section [8.3.1](#S8.SS3.SSS1)): This data structure maintains an implicit representation of the primal-dual solution pair $(x,s)$. This is directly used by CentralPathMaintenance.
- •
ApproxDS (Section [8.3.2](#S8.SS3.SSS2)): This data structure explicitly maintains an approximation $(\overline{x},\overline{s})$ of $(x,s)$. This data structure is directly used by CentralPathMaintenance.
- •
BatchSketch (Section [8.3.3](#S8.SS3.SSS3)): This data structure maintains a sketch of $(x,s)$. This data structure is used by ApproxDS.

Notation: In this section, for simplicity, we write $B_{\overline{x}}$ for $B_{w,\overline{x},\overline{t}}$, and $L_{\overline{x}}$ for the Cholesky factor of $B_{\overline{x}}$, i.e., $B_{\overline{x}}=L_{\overline{x}}L_{\overline{x}}^{\top}$.

#### 8.3.1 ExactDS

In this section we present the data structure ExactDS. It maintains an implicit representation of the primal-dual solution pair $(x,s)$ by maintaining several sparsely-changing vectors (see Eq. ([12](#S8.E12))([13](#S8.E13))).

###### Theorem 8.3 .

Data structure ExactDS (Algorithm [14](#alg14), [15](#alg15)) implicitly maintains the primal-dual pair $(x,s)\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$, computable via the expression

$$ $\displaystyle x$ $\displaystyle=\widehat{x}+H_{w,\overline{x}}^{-1/2}h\beta_{x}+H_{w,\overline{x}}^{-1/2}\widehat{h}\widehat{\beta}_{x}+H_{w,\overline{x}}^{-1/2}\widetilde{h}\widetilde{\beta}_{x},$ (12) $\displaystyle s$ $\displaystyle=\widehat{s}+H_{w,\overline{x}}^{1/2}h\beta_{s}+H_{w,\overline{x}}^{1/2}\widehat{h}\widehat{\beta}_{s}+H_{w,\overline{x}}^{1/2}\widetilde{h}\widetilde{\beta}_{s},$ (13) $$

where
$\widehat{x},\widehat{s}\in\mathbb{R}^{n_{\mathrm{tot}}}$, $h=H_{w,\overline{x}}^{-1/2}\overline{\delta}_{\mu}\in\mathbb{R}^{n_{\mathrm{tot}}}$, $\widehat{h}=H_{w,\overline{x}}^{-1/2}U^{\top}\in\mathbb{R}^{n_{\mathrm{tot}}\times k}$, $\widetilde{h}=H_{w,\overline{x}}^{-1/2}A^{\top}\in\mathbb{R}^{n_{\mathrm{tot}}\times m}$, $\beta_{x},\beta_{s}\in\mathbb{R}$, $\widehat{\beta}_{x},\widehat{\beta}_{s}\in\mathbb{R}^{k}$, $\widetilde{\beta}_{x},\widetilde{\beta}_{s}\in\mathbb{R}^{m}$.

The data structure supports the following functions:

- •
Initialize$(x,s,\overline{x},\overline{s}\in\mathbb{R}^{n_{\mathrm{tot}}},\overline{t}\in\mathbb{R}_{>0})$: Initializes the data structure in
$\widetilde{O}(n(k^{\omega}+m^{\omega}))$
time, with initial value of the primal-dual pair $(x,s)$, its initial approximation $(\overline{x},\overline{s})$, and initial approximate timestamp $\overline{t}$.
- •
Move$()$: Performs robust central path step
$\displaystyle x$
$\displaystyle\leftarrow x+\overline{t}B_{\overline{x}}^{-1}\delta_{\mu}-\overline{t}B_{\overline{x}}^{-1}A^{\top}(AB_{\overline{x}}^{-1}A^{\top})^{-1}AB_{\overline{x}}^{-1}\delta_{\mu},$
(14)
$\displaystyle s$
$\displaystyle\leftarrow s+\overline{t}\delta_{\mu}-\overline{t}^{2}B_{\overline{x}}^{-1}\delta_{\mu}+\overline{t}^{2}B_{\overline{x}}^{-1}A^{\top}(AB_{\overline{x}}^{-1}A^{\top})^{-1}AB_{\overline{x}}^{-1}\delta_{\mu}$
(15)
in $O(k^{\omega}+m^{\omega})$ time by updating its implicit representation.
- •
Update$(\delta_{\overline{x}},\delta_{\overline{s}}\in\mathbb{R}^{n_{\mathrm{tot}}})$: Updates the approximation pair $(\overline{x},\overline{s})$ to $(\overline{x}^{\mathrm{new}}=\overline{x}+\delta_{\overline{x}}\in\mathbb{R}^{n_{\mathrm{tot}}},\overline{s}^{\mathrm{new}}=\overline{s}+\delta_{\overline{s}}\in\mathbb{R}^{n_{\mathrm{tot}}})$ in
$\widetilde{O}((k^{2}+m^{2})(\|\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{s}}\|_{0}))$
time, and output the changes in variables $h,\widehat{h},\widetilde{h},H_{w,\overline{x}}^{1/2}\widehat{x},H_{w,\overline{x}}^{-1/2}\widehat{s}$.
Furthermore, $h,H_{w,\overline{x}}^{1/2}\widehat{x},H_{w,\overline{x}}^{-1/2}\widehat{s}$ changes in $O(\|\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{s}}\|_{0})$ coordinates,
$\widehat{h}$ changes in $O(k(\|\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{s}}\|_{0}))$ coordinates,
$\widetilde{h}$ changes in $O(m(\|\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{s}}\|_{0}))$ coordinates.
- •
Output$()$: Output $x$ and $s$ in $\widetilde{O}(n(k+m))$ time.
- •
Query$x(i\in[n])$: Output $x_{i}$ in $\widetilde{O}(k+m)$ time.
This function is used by ApproxDS.
- •
Query$s(i\in[n])$: Output $s_{i}$ in $\widetilde{O}(k+m)$ time.
This function is used by ApproxDS.

Figure: Algorithm 14 This is used in Algorithm [13](#alg13).

Figure: Algorithm 15 Algorithm [14](#alg14) continued.

###### Proof of Theorem 8.3 .

By combining Lemma [8.4](#S8.Thmtheorem4) and [8.5](#S8.Thmtheorem5).
∎

###### Lemma 8.4 .

ExactDS correctly maintains an implicit representation of $(x,s)$, i.e., invariant

$$ $\displaystyle x=\widehat{x}+H_{w,\overline{x}}^{-1/2}h\beta_{x}+H_{w,\overline{x}}^{-1/2}\widehat{h}\widehat{\beta}_{x}+H_{w,\overline{x}}^{-1/2}\widetilde{h}\widetilde{\beta}_{x},$ $\displaystyle s=\widehat{s}+H_{w,\overline{x}}^{1/2}h\beta_{s}+H_{w,\overline{x}}^{1/2}\widehat{h}\widehat{\beta}_{s}+H_{w,\overline{x}}^{1/2}\widetilde{h}\widetilde{\beta}_{s},$ $\displaystyle h=H_{w,\overline{x}}^{-1/2}\overline{\delta}_{\mu}\in\mathbb{R}^{n_{\mathrm{tot}}},\widehat{h}=H_{w,\overline{x}}^{-1/2}U^{\top}\in\mathbb{R}^{n_{\mathrm{tot}}\times d},\widetilde{h}=H_{w,\overline{x}}^{-1/2}A^{\top}\in\mathbb{R}^{n_{\mathrm{tot}}\times m},$ $\displaystyle u_{1}=UH_{w,\overline{x}}^{-1}A^{\top}\in\mathbb{R}^{d\times m},u_{2}=VH_{w,\overline{x}}^{-1}A^{\top}\in\mathbb{R}^{d\times m},u_{3}=AH_{w,\overline{x}}^{-1}A^{\top}\in\mathbb{R}^{m\times m},$ $\displaystyle u_{4}=AH_{w,\overline{x}}^{-1}\overline{\delta}_{\mu}\in\mathbb{R}^{m},u_{5}=VH_{w,\overline{x}}^{-1}\overline{\delta}_{\mu}\in\mathbb{R}^{d},u_{6}=VH_{w,\overline{x}}^{-1}U^{\top}\in\mathbb{R}^{d\times d},$ $\displaystyle\overline{\alpha}=\sum_{i\in[n]}w_{i}^{-1}\cosh^{2}(\frac{\lambda}{w_{i}}\gamma_{i}(\overline{x},\overline{s},\overline{t})),$ $\displaystyle\overline{\delta}_{\mu}=\overline{\alpha}^{1/2}\delta_{\mu}(\overline{x},\overline{s},\overline{t})$ $$

always holds after every external call, and return values of the queries are correct.

###### Proof.

Initialize:
By checking the definitions we see that all invariants are satisfied after Initialize.

Move:
By the invariants, we have

$$ $\displaystyle v_{0}$ $\displaystyle=I+\overline{t}^{-1}VH_{w,\overline{x}}^{-1}U^{\top},$ $\displaystyle v_{1}$ $\displaystyle=\overline{t}^{-1}AH_{w,\overline{x}}^{-1}A^{\top}-\overline{t}^{-1}AH_{w,\overline{x}}^{-1}U^{\top}(I+\overline{t}^{-1}VH_{w,\overline{x}}^{-1}U^{\top})^{-1}VH_{w,\overline{x}}A^{\top}$ $\displaystyle=AB_{\overline{x}}^{-1}A^{\top}$ $\displaystyle v_{2}$ $\displaystyle=\overline{t}^{-1}AH_{w,\overline{x}}^{-1}\overline{\delta}_{\mu}-\overline{t}^{-1}AH_{w,\overline{x}}^{-1}U^{\top}(I+\overline{t}^{-1}VH_{w,\overline{x}}^{-1}U^{\top})^{-1}VH_{w,\overline{x}}\overline{\delta}_{\mu}$ $\displaystyle=AB_{\overline{x}}^{-1}\overline{\delta}_{\mu}.$ $$

By implicit representation ([12](#S8.E12)),

$$ $\displaystyle\delta_{x}=$ $\displaystyle~{}H_{w,\overline{x}}^{-1/2}h\delta_{\beta_{x}}+H_{w,\overline{x}}^{-1/2}\widehat{h}\delta_{\widehat{\beta}_{x}}+H_{w,\overline{x}}^{-1/2}\widetilde{h}\delta_{\widetilde{\beta}_{x}}$ $\displaystyle=$ $\displaystyle~{}H_{w,\overline{x}}^{-1}\overline{\delta}_{\mu}\cdot(\overline{\alpha})^{-1/2}$ $\displaystyle~{}+H_{w,\overline{x}}^{-1}U^{\top}\cdot(\overline{\alpha})^{-1/2}\overline{t}^{-1}v_{0}^{-1}(-u_{5}+u_{2}v_{1}^{-1}v_{2})$ $\displaystyle~{}-H_{w,\overline{x}}^{-1}A^{\top}\cdot(\overline{\alpha})^{-1/2}v_{1}^{-1}v_{2}$ $\displaystyle=$ $\displaystyle~{}H_{w,\overline{x}}^{-1}\delta_{\mu}$ $\displaystyle~{}+H_{w,\overline{x}}^{-1}U^{\top}\overline{t}^{-1}(I+\overline{t}^{-1}VH_{w,\overline{x}}^{-1}U^{\top})^{-1}(-VH_{w,\overline{x}}^{-1}\delta_{\mu}+VH_{w,\overline{x}}^{-1}A^{\top}(AB_{\overline{x}}^{-1}A^{\top})^{-1}AB_{\overline{x}}^{-1}\delta_{\mu})$ $\displaystyle~{}-H_{w,\overline{x}}^{-1}A^{\top}(AB_{\overline{x}}^{-1}A^{\top})^{-1}AB_{\overline{x}}^{-1}\delta_{\mu}$ $\displaystyle=$ $\displaystyle~{}\overline{t}\cdot(\overline{t}^{-1}H_{w,\overline{x}}^{-1}-\overline{t}^{-2}H_{w,\overline{x}}^{-1}U^{\top}(I+\overline{t}^{-1}VH_{w,\overline{x}}^{-1}U^{\top})^{-1}VH_{w,\overline{x}}^{-1})\delta_{\mu}$ $\displaystyle~{}-\overline{t}(\overline{t}^{-1}H_{w,\overline{x}}^{-1}-\overline{t}^{-2}\overline{t}^{-2}H_{w,\overline{x}}^{-1}U^{\top}(I+\overline{t}^{-1}VH_{w,\overline{x}}^{-1}U^{\top})^{-1}VH_{w,\overline{x}}^{-1})A^{\top}(AB_{\overline{x}}^{-1}A^{\top})^{-1}AB_{\overline{x}}^{-1}\delta_{\mu}$ $\displaystyle=$ $\displaystyle~{}\overline{t}B_{\overline{x}}^{-1}\delta_{\mu}-\overline{t}B_{\overline{x}}^{-1}A^{\top}(AB_{\overline{x}}^{-1}A^{\top})^{-1}AB_{\overline{x}}^{-1}\delta_{\mu}.$ $$

Comparing with the robust central path step ([14](#S8.E14)), we see that $x$ is updated correctly.

For $s$, from implicit representation [13](#S8.E13) we have

$$ $\displaystyle\delta_{s}=$ $\displaystyle~{}H_{w,\overline{x}}^{1/2}h\delta_{\beta_{x}}+H_{w,\overline{x}}^{1/2}\widehat{h}\delta_{\widehat{\beta}_{x}}+H_{w,\overline{x}}^{1/2}\widetilde{h}\delta_{\widetilde{\beta}_{x}}$ $\displaystyle=$ $\displaystyle~{}-U^{\top}\cdot(\overline{\alpha})^{-1/2}\cdot v_{0}^{-1}(-u_{5}+u_{2}v_{1}^{-1}v_{2})+A^{\top}\cdot(\overline{\alpha})^{-1/2}\cdot\overline{t}v_{1}^{-1}v_{2}$ $\displaystyle=$ $\displaystyle~{}\overline{t}\delta_{\mu}-\overline{t}^{2}B_{\overline{x}}^{-1}\delta_{\mu}+\overline{t}^{2}B_{\overline{x}}^{-1}A^{\top}(AB_{\overline{x}}^{-1}A^{\top})^{-1}AB_{\overline{x}}^{-1}\delta_{\mu}.$ $$

Comparing with robust central path step ([15](#S8.E15)), we see that $s$ is updated correctly.

Update:
We would like to prove that Update correctly updates the values of $\widehat{x},\widehat{s},h,\widehat{h},\widetilde{h}$, $u_{1},u_{2},u_{3},u_{4},u_{5},u_{6}$, $\overline{\alpha}$, $\overline{\delta}_{\mu}$, while preserving the values of $(x,s)$.
In fact, by checking the definitions, it is easy to see that $h,\widehat{h},\widetilde{h}$, $u_{1},u_{2},u_{3},u_{4},u_{5},u_{6}$, $\overline{\alpha}$, $\overline{\delta}_{\mu}$ are updated correctly.
Furthermore

$$ $\displaystyle\delta_{x}$ $\displaystyle=\delta_{\widehat{x}}+\delta_{h}\beta_{x}+\delta_{\widehat{h}}\widehat{\beta}_{x}+\delta_{\widetilde{h}}\widetilde{\beta}_{x}=0,$ $\displaystyle\delta_{s}$ $\displaystyle=\delta_{\widehat{s}}+\delta_{h}\beta_{s}+\delta_{\widehat{h}}\widehat{\beta}_{s}+\delta_{\widetilde{h}}\widetilde{\beta}_{s}=0.$ $$

So values of $(x,s)$ are preserved.

∎

###### Lemma 8.5 .

We bound the running time of ExactDS as following.

- (i)
ExactDS.Initialize (Algorithm [14](#alg14)) runs in $\widetilde{O}(n(k^{\omega-1}+m^{\omega-1}))$ time.
- (ii)
ExactDS.Move (Algorithm [14](#alg14)) runs in $\widetilde{O}(k^{\omega}+m^{\omega})$ time.
- (iii)
ExactDS.Output (Algorithm [15](#alg15)) runs in $\widetilde{O}(n(k+m))$ time and correctly outputs $(x,s)$.
- (iv)
ExactDS.Query$x$ and ExactDS.Query$s$ (Algorithm [15](#alg15)) runs in $\widetilde{O}(k+m)$ time and returns the correct answer.
- (v)
ExactDS.Update (Algorithm [15](#alg15)) runs in $\widetilde{O}((k^{2}+m^{2})(\|\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{s}}\|_{0}))$ time.
Furthermore,
$\|\delta_{h}\|_{0},\|\delta_{\widehat{x}}\|_{0},\|\delta_{\widehat{s}}\|_{0}=O(\|\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{s}}\|_{0})$,
$\operatorname{nnz}(\widehat{h})=O(d(\|\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{s}}\|_{0}))$,
$\operatorname{nnz}(\widetilde{h})=O(m(\|\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{s}}\|_{0}))$.

###### Proof.

- (i)
ExactDS.Initialize:
Computing $u_{1}$ and $u_{2}$ takes ${\cal T}_{\mathrm{mat}}(k,n,m)=\widetilde{O}(n(k^{\omega-1}+m^{\omega-1}))$ time.
Computing $u_{3}$ takes ${\cal T}_{\mathrm{mat}}(m,n,m)=\widetilde{O}(nm^{\omega-1})$ time.
Computing $u_{4}$ takes $O(nm)$ time.
Computing $u_{5}$ takes $O(nk)$ time.
Computing $u_{6}$ takes ${\cal T}_{\mathrm{mat}}(k,n,k)=\widetilde{O}(nk^{\omega-1})$ time.
All other computations are cheaper.
- (ii)
ExactDS.Move: Computing $v_{0}^{-1}$ takes $\widetilde{O}(k^{\omega})$ time.
Computing $v_{1}^{-1}$ takes $\widetilde{O}(m^{\omega})$ time.
All other computations are cheaper.
- (iii)
ExactDS.Output: Takes $\widetilde{O}(n(k+m))$ time.
- (iv)
ExactDS.Query$x$ and ExactDS.Query$s$: Takes $\widetilde{O}(k+m)$ time.
- (v)
ExactDS.Update:
For simplicity, write $t=\|\delta_{\overline{x}}\|_{0}+\|\delta_{\overline{x}}\|_{0}$.
Computing $\delta_{h}$ takes $\widetilde{O}(t)$ time.
Computing $\delta_{\widehat{h}}$ takes $\widetilde{O}(tk)$ time.
Computing $\delta_{\widetilde{h}}$ takes $\widetilde{O}(tm)$ time.
Computing $\delta_{\widehat{x}}$ and $\delta_{\widehat{s}}$ takes $\widetilde{O}(t(k+m))$ time.
The sparsity statements follow directly.
Computing $u_{1}$ and $u_{2}$ takes $\widetilde{O}(tkm)$ time.
Computing $u_{3}$ takes $\widetilde{O}(tm^{2})$ time.
Computing $u_{4}$ takes $\widetilde{O}(tm)$ time.
Computing $u_{5}$ takes $\widetilde{O}(tk)$ time.
Computing $u_{6}$ takes $\widetilde{O}(tk^{2})$ time.

∎

###### Theorem 8.3 .

###### Proof of Theorem 8.3 .

###### Lemma 8.4 .

###### Proof.

###### Lemma 8.5 .

###### Proof.

#### 8.3.2 ApproxDS

In this section we present the data structure ApproxDS.
Given BatchSketch, a data structure maintaining a sketch of the primal-dual pair $(x,s)\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$, ApproxDS maintains a sparsely-changing $\ell_{\infty}$-approximation of $(x,s)$.

Figure: Algorithm 16 This is used in Algorithm [13](#alg13).

###### Theorem 8.6 .

Given parameters $\epsilon_{\mathrm{apx},x},\epsilon_{\mathrm{apx},s}\in(0,1),\delta_{\mathrm{apx}}\in(0,1)$, $\zeta_{x},\zeta_{s}\in\mathbb{R}$ such that

$$ $\displaystyle\|H_{w,\overline{x}^{(\ell)}}^{1/2}x^{(\ell)}-H_{w,\overline{x}^{(\ell)}}^{1/2}x^{(\ell+1)}\|_{2}\leq\zeta_{x},~{}~{}~{}\|H_{w,\overline{x}^{(\ell)}}^{-1/2}s^{(\ell)}-H_{w,\overline{x}^{(\ell)}}^{-1/2}s^{(\ell+1)}\|_{2}\leq\zeta_{s}$ $$

for all $\ell\in\{0,\ldots,q-1\}$,
data structure ApproxDS (Algorithm [16](#alg16)) supports the following operations:

- •
$\textsc{Initialize}(x,s\in\mathbb{R}^{n_{\mathrm{tot}}},h\in\mathbb{R}^{n_{\mathrm{tot}}},\widehat{h}\in\mathbb{R}^{n_{\mathrm{tot}}\times k},\widetilde{h}\in\mathbb{R}^{n_{\mathrm{tot}}\times m},H_{w,\overline{x}}^{1/2}\widehat{x},H_{w,\overline{x}}^{-1/2}\widehat{s}\in\mathbb{R}^{n_{\mathrm{tot}}},\beta_{x},\beta_{s}\in\mathbb{R},\widehat{\beta}_{x},\widehat{\beta}_{s}\in\mathbb{R}^{k},\widetilde{\beta}_{x},\widetilde{\beta}_{s}\in\mathbb{R}^{m},q\in\mathbb{N},\textsc{ExactDS*}~{}\mathsf{exact},\epsilon_{\mathrm{apx},x},\epsilon_{\mathrm{apx},s},\delta_{\mathrm{apx}}\in\mathbb{R})$: Initialize the data structure in $\widetilde{O}(n(k+m))$ time.
- •
$\textsc{MoveAndQuery}(\beta_{x},\beta_{s}\in\mathbb{R},\widehat{\beta}_{x},\widehat{\beta}_{s}\in\mathbb{R}^{d},\widetilde{\beta}_{x},\widetilde{\beta}_{s}\in\mathbb{R}^{m})$:
Update values of $\beta_{x},\beta_{s},\widehat{\beta}_{x},\widehat{\beta}_{s},\widetilde{\beta}_{x},\widetilde{\beta}_{s}$ by calling $\textsc{BatchSketch}.\textsc{Move}$.
This effectively moves $(x^{(\ell)},s^{(\ell)})$ to $(x^{(\ell+1)},s^{(\ell+1)})$ while keeping $\overline{x}^{(\ell)}$ unchanged.
Then return two sets $L_{x}^{(\ell)},L_{s}^{(\ell)}\subset[n]$ where
$\displaystyle L_{x}^{(\ell)}$
$\displaystyle\supseteq\{i\in[n]:\|H_{w,\overline{x}^{(\ell)}}^{1/2}x^{(\ell)}_{i}-H_{w,\overline{x}^{(\ell)}}^{1/2}x^{(\ell+1)}_{i}\|_{2}\geq\epsilon_{\mathrm{apx},x}\},$
$\displaystyle L_{s}^{(\ell)}$
$\displaystyle\supseteq\{i\in[n]:\|H_{w,\overline{x}^{(\ell)}}^{-1/2}s^{(\ell)}_{i}-H_{w,\overline{x}^{(\ell)}}^{-1/2}s^{(\ell+1)}_{i}\|_{2}\geq\epsilon_{\mathrm{apx},s}\},$
satisfying
$\displaystyle\sum_{0\leq\ell\leq q-1}|L_{x}^{(\ell)}|=\widetilde{O}(\epsilon_{\mathrm{apx},x}^{-2}\zeta_{x}^{2}q^{2}),$
$\displaystyle\sum_{0\leq\ell\leq q-1}|L_{s}^{(\ell)}|=\widetilde{O}(\epsilon_{\mathrm{apx},s}^{-2}\zeta_{s}^{2}q^{2}).$
For every query, with probability at least $1-\delta_{\mathrm{apx}}/q$, the return values are correct.
Furthermore, total time cost over all queries is at most
$\displaystyle\widetilde{O}\left((\epsilon_{\mathrm{apx},x}^{-2}\zeta_{x}^{2}+\epsilon_{\mathrm{apx},s}^{-2}\zeta_{s}^{2})q^{2}(k+m)\right).$
- •
$\textsc{Update}(\delta_{\overline{x}}\in\mathbb{R}^{n_{\mathrm{tot}}},\delta_{h}\in\mathbb{R}^{n_{\mathrm{tot}}},\delta_{\widehat{h}}\in\mathbb{R}^{n_{\mathrm{tot}}\times d},\delta_{\widetilde{h}}\in\mathbb{R}^{n_{\mathrm{tot}}\times m},\delta_{H_{w,\overline{x}}^{1/2}\widehat{x}},\delta_{H_{w,\overline{x}}^{-1/2}\widehat{s}}\in\mathbb{R}^{n_{\mathrm{tot}}})$: Update sketches of $H_{w,\overline{x}^{(\ell)}}^{1/2}x^{(\ell+1)}$ and $H_{w,\overline{x}^{(\ell)}}^{-1/2}s^{(\ell+1)}$ by calling BatchSketch.Update.
This effectively moves $\overline{x}^{(\ell)}$ to $\overline{x}^{(\ell+1)}$ while keeping $(x^{(\ell+1)},s^{(\ell+1)})$ unchanged. Then advance timestamp $\ell$.
Each update costs
$\displaystyle\widetilde{O}(\|\delta_{h}\|_{0}+\operatorname{nnz}(\delta_{\widehat{h}})+\operatorname{nnz}(\delta_{\widetilde{h}})+\|H_{w,\overline{x}}^{1/2}\widehat{x}\|_{0}+\|H_{w,\overline{x}}^{-1/2}\widehat{s}\|_{0})$
time.

###### Proof.

The proof is essentially the same as proof of [, Theorem 4.18].
For the running time claims, we plug in Theorem [8.7](#S8.Thmtheorem7) when necessary.
∎

###### Theorem 8.6 .

###### Proof.

#### 8.3.3 BatchSketch

In this section we present the data structure BatchSketch. It maintains a sketch of $H_{\overline{x}}^{1/2}x$ and $H_{\overline{x}}^{-1/2}s$. It is a variation of BatchSketch in [].

Figure: Algorithm 17 This is used by Algorithm [16](#alg16).

Figure: Algorithm 18 BatchSketch Algorithm [17](#alg17) continued. This is used by Algorithm [16](#alg16).

###### Theorem 8.7 .

Data structure BatchSketch (Algorithm [17](#alg17), [18](#alg18)) supports the following operations:

- •
$\textsc{Initialize}(\overline{x}\in\mathbb{R}^{n_{\mathrm{tot}}},h\in\mathbb{R}^{n_{\mathrm{tot}}},\widehat{h}\in\mathbb{R}^{n_{\mathrm{tot}}\times k},\widetilde{h}\in\mathbb{R}^{n_{\mathrm{tot}}\times m},H_{w,\overline{x}}^{1/2}\widehat{x},H_{w,\overline{x}}^{-1/2}\widehat{s}\in\mathbb{R}^{n_{\mathrm{tot}}},\beta_{x},\beta_{s}\in\mathbb{R},\widehat{\beta}_{x},\widehat{\beta}_{s}\in\mathbb{R}^{k},\widetilde{\beta}_{x},\widetilde{\beta}_{s}\in\mathbb{R}^{m},\delta_{\mathrm{apx}}\in\mathbb{R})$: Initialize the data structure in $\widetilde{O}(n(k+m))$ time.
- •
$\textsc{Move}(\beta_{x},\beta_{s}\in\mathbb{R},\widehat{\beta}_{x},\widehat{\beta}_{s}\in\mathbb{R}^{k},\widetilde{\beta}_{x},\widetilde{\beta}_{s}\in\mathbb{R}^{m})$: Update values of $\beta_{x},\beta_{s},\widehat{\beta}_{x},\widehat{\beta}_{s},\widetilde{\beta}_{x},\widetilde{\beta}_{s}$ in $O(k+m)$ time. This effectively moves $(x^{(\ell)},s^{(\ell)})$ to $(x^{(\ell+1)},s^{(\ell+1)})$ while keeping $\overline{x}^{(\ell)}$ unchanged.
- •
$\textsc{Update}(\delta_{\overline{x}}\in\mathbb{R}^{n_{\mathrm{tot}}},\delta_{h}\in\mathbb{R}^{n_{\mathrm{tot}}},\delta_{\widehat{h}}\in\mathbb{R}^{n_{\mathrm{tot}}\times k},\delta_{\widetilde{h}}\in\mathbb{R}^{n_{\mathrm{tot}}\times m},\delta_{H_{w,\overline{x}}^{1/2}\widehat{x}},\delta_{H_{w,\overline{x}}^{-1/2}\widehat{s}}\in\mathbb{R}^{n_{\mathrm{tot}}})$:
Update sketches of $H_{w,\overline{x}^{(\ell)}}^{1/2}x^{(\ell+1)}$ and $H_{w,\overline{x}^{(\ell)}}^{-1/2}s^{(\ell+1)}$.
This effectively moves $\overline{x}^{(\ell)}$ to $\overline{x}^{(\ell+1)}$ while keeping $(x^{(\ell+1)},s^{(\ell+1)})$ unchanged.
Then advance timestamp $\ell$.
Each update costs
$\displaystyle\widetilde{O}(\|\delta_{h}\|_{0}+\operatorname{nnz}(\delta_{\widehat{h}})+\operatorname{nnz}(\delta_{\widetilde{h}})+\|H_{w,\overline{x}}^{1/2}\widehat{x}\|_{0}+\|H_{w,\overline{x}}^{-1/2}\widehat{s}\|_{0})$
time.
- •
$\textsc{Query}x(\ell^{\prime}\in\mathbb{N},\epsilon\in\mathbb{R})$:
Given timestamp $\ell^{\prime}$, return a set $S\subseteq[n]$ where
$\displaystyle S$
$\displaystyle\supseteq\{i\in[n]:\|H_{w,\overline{x}^{(\ell^{\prime})}}^{1/2}x^{(\ell^{\prime})}_{i}-H_{w,\overline{x}^{(\ell)}}^{1/2}x^{(\ell+1)}_{i}\|_{2}\geq\epsilon\},$
and
$\displaystyle|S|$
$\displaystyle=O(\epsilon^{-2}(\ell-\ell^{\prime}+1)\sum_{\ell^{\prime}\leq t\leq\ell}\|H_{w,\overline{x}^{(t)}}^{1/2}x^{(t)}-H_{w,\overline{x}^{(t)}}^{1/2}x^{(t+1)}\|_{2}^{2}+\sum_{\ell^{\prime}\leq t\leq\ell-1}\|\overline{x}^{(t)}-\overline{x}^{(t+1)}\|_{2,0})$
where $\ell$ is the current timestamp.
For every query, with probability at least $1-\delta$, the return values are correct, and costs at most
$\displaystyle\widetilde{O}((k+m)\cdot(\epsilon^{-2}(\ell-\ell^{\prime}+1)\sum_{\ell^{\prime}\leq t\leq\ell}\|H_{\overline{x}^{(t)}}^{1/2}x^{(t)}-H_{\overline{x}^{(t)}}^{1/2}x^{(t+1)}\|_{2}^{2}+\sum_{\ell^{\prime}\leq t\leq\ell-1}\|\overline{x}^{(t)}-\overline{x}^{(t+1)}\|_{2,0}))$
running time.
- •
$\textsc{Query}s(\ell^{\prime}\in\mathbb{N},\epsilon\in\mathbb{R})$:
Given timestamp $\ell^{\prime}$, return a set $S\subseteq[n]$ where
$\displaystyle S$
$\displaystyle\supseteq\{i\in[n]:\|H_{w,\overline{x}^{(\ell^{\prime})}}^{-1/2}s^{(\ell^{\prime})}_{i}-H_{w,\overline{x}^{(\ell)}}^{-1/2}s^{(\ell+1)}_{i}\|_{2}\geq\epsilon\}$
and
$\displaystyle|S|$
$\displaystyle=O(\epsilon^{-2}(\ell-\ell^{\prime}+1)\sum_{\ell^{\prime}\leq t\leq\ell}\|H_{w,\overline{x}^{(t)}}^{-1/2}s^{(t)}-H_{w,\overline{x}^{(t)}}^{-1/2}s^{(t+1)}\|_{2}^{2}+\sum_{\ell^{\prime}\leq t\leq\ell-1}\|\overline{x}^{(t)}-\overline{x}^{(t+1)}\|_{2,0})$
where $\ell$ is the current timestamp.
For every query, with probability at least $1-\delta$, the return values are correct, and costs at most
$\displaystyle\widetilde{O}((k+m)\cdot(\epsilon^{-2}(\ell-\ell^{\prime}+1)\sum_{\ell^{\prime}\leq t\leq\ell}\|H_{\overline{x}^{(t)}}^{1/2}s^{(t)}-H_{\overline{x}^{(t)}}^{1/2}x^{(t+1)}\|_{2}^{2}+\sum_{\ell^{\prime}\leq t\leq\ell-1}\|\overline{x}^{(t)}-\overline{x}^{(t+1)}\|_{2,0}))$
running time.

###### Proof.

The proof is essentially the same as proof of [, Theorem 4.21].
∎

###### Theorem 8.7 .

###### Proof.

### 8.4 Analysis of CentralPathMaintenance

###### Lemma 8.8 (Correctness of CentralPathMaintenance ) .

Algorithm [13](#alg13) implicitly maintains the primal-dual solution pair $(x,s)$ via representation Eq. ([12](#S8.E12))([13](#S8.E13)).
It also explicitly maintains $(\overline{x},\overline{s})\in\mathbb{R}^{n_{\mathrm{tot}}}\times\mathbb{R}^{n_{\mathrm{tot}}}$ such that $\|\overline{x}_{i}-x_{i}\|_{\overline{x}_{i}}\leq\overline{\epsilon}$ and $\|\overline{s}_{i}-s_{i}\|_{\overline{x}_{i}}^{*}\leq t\overline{\epsilon}w_{i}$ for all $i\in[n]$ with probability at least $0.9$.

###### Proof.

Same as proof of Lemma [7.13](#S7.Thmtheorem13).
∎

###### Lemma 8.9 .

We bound the running time of CentralPathMaintenance as following.

- •
CentralPathMaintenance.Initialize takes $\widetilde{O}(n(k^{\omega-1}+m^{\omega-1}))$ time.
- •
If CentralPathMaintenance.MultiplyAndMove is called $N$ times, then it has total running time
$\displaystyle\widetilde{O}((Nn^{-1/2}+\log(t_{\max}/t_{\min}))\cdot n(k+m)^{(\omega+1)/2}).$
- •
CentralPathMaintenance.Output takes $\widetilde{O}(n(k+m))$ time.

###### Proof.

Initialize part:
By Theorem [8.3](#S8.Thmtheorem3) and [8.6](#S8.Thmtheorem6).

Output part:
By Theorem [8.3](#S8.Thmtheorem3).

MultiplyAndMove part:
Between two restarts, the total size of $|L_{x}|$ returned by $\mathsf{approx}$.Query is bounded by $\widetilde{O}(q^{2}\zeta_{x}^{2}/\epsilon_{\mathrm{apx},x}^{2})$ by Theorem [8.6](#S8.Thmtheorem6).
By plugging in $\zeta_{x}=2\alpha$, $\epsilon_{\mathrm{apx},x}=\overline{\epsilon}$, we have
$\sum_{\ell\in[q]}|L_{x}^{(\ell)}|=\widetilde{O}(q^{2}).$
Similarly, for $s$ we have
$\sum_{\ell\in[q]}|L_{s}^{(\ell)}|=\widetilde{O}(q^{2}).$

Update time:
By Theorem [8.3](#S8.Thmtheorem3) and [8.6](#S8.Thmtheorem6), in a sequence of $q$ updates, total cost for update is $\widetilde{O}(q^{2}(k^{2}+m^{2}))$.
So the amortized update cost per iteration is $\widetilde{O}(q(k^{2}+m^{2}))$.
The total update cost is

$$ $\displaystyle\text{number of iterations}\cdot\text{time per iteration}=\widetilde{O}(Nq(k^{2}+m^{2})).$ $$

Init/restart time:
We restart the data structure whenever $K>q$ or $|\overline{t}-t|>\overline{t}\epsilon_{t}$, so there are $O(N/q+\log(t_{\max}/t_{\min})\epsilon_{t}^{-1})$ restarts in total.
By Theorem [8.3](#S8.Thmtheorem3) and [8.6](#S8.Thmtheorem6), time cost per restart is $\widetilde{O}(n(k^{\omega-1}+m^{\omega-1}))$.
So the total initialization time is

$$ $\displaystyle\text{number of restarts}\cdot\text{time per restart}=\widetilde{O}((N/q+\log(t_{\max}/t_{\min})\epsilon_{t}^{-1})\cdot n(k^{\omega-1}+m^{\omega-1})).$ $$

Combine everything:
Overall running time is

$$ $\displaystyle\widetilde{O}(Nq(k^{2}+m^{2})+(N/q+\log(t_{\max}/t_{\min})\epsilon_{t}^{-1})\cdot n(k^{\omega-1}+m^{\omega-1})).$ $$

Taking $\epsilon_{t}=\frac{1}{2}\overline{\epsilon}$, the optimal choice for $q$ is

$$ $\displaystyle q=n^{1/2}(k^{2}+m^{2})^{-1/2}(k^{\omega-1}+m^{\omega-1})^{1/2},$ $$

achieving overall running time

$$ $\displaystyle~{}\widetilde{O}((Nn^{-1/2}+\log(t_{\max}/t_{\min}))\cdot n(k^{2}+m^{2})^{1/2}(k^{\omega-1}+m^{\omega-1})^{1/2})$ $\displaystyle=$ $\displaystyle~{}\widetilde{O}((Nn^{-1/2}+\log(t_{\max}/t_{\min}))\cdot n(k+m)^{(\omega+1)/2}).$ $$

∎

###### Proof of Theorem 8.2 .

Combining Lemma [8.8](#S8.Thmtheorem8) and [8.9](#S8.Thmtheorem9).
∎

###### Lemma 8.8 (Correctness of CentralPathMaintenance ) .

###### Proof.

###### Lemma 8.9 .

###### Proof.

###### Proof of Theorem 8.2 .

### 8.5 Proof of Main Statement

###### Proof of Theorem 8.1 .

Use CentralPathMaintenance (Algorithm [13](#alg13)) as the maintenance data structure in Algorithm [20](#alg20).
Combining Theorem [8.2](#S8.Thmtheorem2) and Theorem [9.1](#S9.Thmtheorem1) finishes the proof.
∎

###### Proof of Theorem 8.1 .

## 9 Robust IPM Analysis

In this section we present a robust IPM algorithm for quadratic programming.
The algorithm is a modification of previous robust IPM algorithms for linear programming [, ].

Convention: Variables are in $n$ blocks of dimension $n_{i}$ ($i\in[n]$).
Total dimension is $n_{\mathrm{tot}}=\sum_{i\in[n]}n_{i}$.
We write $x=(x_{1},\ldots,x_{n})\in\mathbb{R}^{n_{\mathrm{tot}}}$ where $x_{i}\in\mathbb{R}^{n_{i}}$.
We consider programs of the following form:

$$ $\displaystyle\min_{x\in\mathbb{R}^{n}}~{}$ $\displaystyle~{}\frac{1}{2}x^{\top}Qx+c^{\top}x$ (16) $\displaystyle\mathrm{s.t.~{}}$ $\displaystyle~{}Ax=b$ $\displaystyle~{}x_{i}\in\mathcal{K}_{i}\qquad\forall i\in[n]$ $$

where $Q\in\mathcal{S}^{n_{\mathrm{tot}}}$, $c\in\mathbb{R}^{n_{\mathrm{tot}}}$,
$A\in\mathbb{R}^{m\times n_{\mathrm{tot}}}$, $b\in\mathbb{R}^{m}$, $\mathcal{K}_{i}\subset\mathbb{R}^{n_{i}}$ is a convex set.
Let $\mathcal{K}=\prod_{i\in[n]}\mathcal{K}_{i}$.

###### Theorem 9.1 .

Consider the convex program ([16](#S9.E16)).
Let $\phi_{i}:\mathcal{K}_{i}\to\mathbb{R}$ be a $\nu_{i}$-self-concordant barrier for all $i\in[n]$.
Suppose the program satisfies the following properties:

- •
Inner radius $r$: There exists $z\in\mathbb{R}^{n_{\mathrm{tot}}}$ such that $Az=b$ and $B(z,r)\in\mathcal{K}$.
- •
Outer radius $R$: $\mathcal{K}\subseteq B(0,R)$ where $0\in\mathbb{R}^{n_{\mathrm{tot}}}$.
- •
Lipschitz constant $L$: $\|Q\|_{2\to 2}\leq L$, $\|c\|_{2}\leq L$.

Let $(w_{i})_{i\in[n]}\in\mathbb{R}^{n}_{\geq 1}$ and $\kappa=\sum_{i\in[n]}w_{i}\nu_{i}$.
For any $0<\epsilon\leq\frac{1}{2}$, Algorithm [19](#alg19) outputs an approximate solution $x$ in $O(\sqrt{\kappa}\log n\log\frac{n\kappa R}{\epsilon r})$ steps, satisfying

$$ $\displaystyle\frac{1}{2}x^{\top}Qx+c^{\top}x$ $\displaystyle\leq\min_{Ax=b,x\in\mathcal{K}}\left(\frac{1}{2}x^{\top}Qx+c^{\top}x\right)+\epsilon LR(R+1),$ $\displaystyle\|Ax-b\|_{1}$ $\displaystyle\leq 3\epsilon(R\|A\|_{1}+\|b\|_{1}),$ $\displaystyle x$ $\displaystyle\in\mathcal{K}.$ $$

Figure: Algorithm 19 Our main algorithm

Figure: Algorithm 20 Subroutine used by Algorithm [19](#alg19)

###### Theorem 9.1 .

### 9.1 Preliminaries

Previous works on linear programming (e.g. [], []) use the following path:

$$ $\displaystyle s/t+\nabla\phi_{w}(x)=$ $\displaystyle~{}\mu,$ $\displaystyle Ax=$ $\displaystyle~{}b,$ $\displaystyle A^{\top}y+s=$ $\displaystyle~{}c$ $$

where
$\phi_{w}(x):=\sum_{i=1}^{n}w_{i}\phi_{i}(x_{i}).$

For quadratic programming, we modify the above central path as following:

$$ $\displaystyle s/t+\nabla\phi_{w}(x)=$ $\displaystyle~{}\mu,$ $\displaystyle Ax=$ $\displaystyle~{}b,$ $\displaystyle-Qx+A^{\top}y+s=$ $\displaystyle~{}c.$ $$

We make the following definitions.

###### Definition 9.2 .

For each $i\in[n]$, we define the $i$-th coordinate error

$$ $\displaystyle\mu_{i}(x,s,t):=\frac{s_{i}}{t}+w_{i}\nabla\phi_{i}(x_{i})$ (17) $$

We define $\mu_{i}$’s norm as

$$ $\displaystyle\gamma_{i}(x,s,t):=\|\mu_{i}(x,s,t)\|_{x_{i}}^{*}.$ (18) $$

We define the soft-max function by

$$ $\displaystyle\Psi_{\lambda}(r):=\sum_{i=1}^{m}\cosh(\lambda\frac{r_{i}}{w_{i}})$ (19) $$

for some $\lambda>0$ and the potential function is the soft-max of the norm of the error of each coordinate

$$ $\displaystyle\Phi(x,s,t)=\Psi_{\lambda}(\gamma(x,s,t))$ (20) $$

We choose the step direction $\delta_{\mu}$ as

$$ $\displaystyle\delta_{\mu,i}:=-\alpha\cdot c_{i}(x,s,t)\cdot\mu_{i}(x,s,t)$ (21) $$

where

$$ $\displaystyle c_{i}(x,s,t):=\frac{\sinh(\frac{\lambda}{w_{i}}\gamma_{i}(x,s,t))}{\gamma_{i}(x,s,t)\sqrt{\sum_{j\in[n]}w_{j}^{-1}\cosh^{2}(\frac{\lambda}{w_{j}}\gamma_{j}(x,s,t))}}$ (22) $$

We define induced norms as following. Note that we include the weight vector $w$ in the subscript to avoid confusion.

###### Definition 9.3 .

For each block $\mathcal{K}_{i}$, we define

$$ $\displaystyle\|v\|_{x_{i}}:=$ $\displaystyle~{}\|v\|_{\nabla^{2}\phi_{i}(x_{i})},$ $\displaystyle\|v\|_{x_{i}}^{*}:=$ $\displaystyle~{}\|v\|_{(\nabla^{2}\phi_{i}(x_{i}))^{-1}}$ $$

for $v\in\mathbb{R}^{n_{i}}$.

For the whole domain $\mathcal{K}=\prod_{i=1}^{n}\mathcal{K}_{i}$, we define

$$ $\displaystyle\|v\|_{w,x}:=$ $\displaystyle~{}\|v\|_{\nabla^{2}\phi_{w}(x)}=(\sum_{i=1}^{n}w_{i}\|v_{i}\|_{x_{i}}^{2})^{1/2},$ $\displaystyle\|v\|_{w,x}^{*}:=$ $\displaystyle~{}\|v\|_{(\nabla^{2}\phi_{w}(x))^{-1}}=(\sum_{i=1}^{n}w_{i}^{-1}(\|v_{i}\|_{x_{i}}^{*})^{2})^{1/2}$ $$

for $v\in\mathbb{R}^{n_{\mathrm{tot}}}$.

The Hessian matrices of the barrier functions appear a lot in the computation.

###### Definition 9.4 .

We define matrices $H_{x,i}\in\mathbb{R}^{n_{i}\times n_{i}}$ and $H_{w,x}\in\mathbb{R}^{n_{\mathrm{tot}}\times n_{\mathrm{tot}}}$ as

$$ $\displaystyle H_{x,i}$ $\displaystyle:=\nabla^{2}\phi_{i}(x_{i}),$ (23) $\displaystyle H_{w,x}$ $\displaystyle:=\nabla^{2}\phi_{w}(x).$ (24) $$

From the definition, we see that

$$ $\displaystyle H_{w,x,(i,i)}$ $\displaystyle=w_{i}H_{x,i}.$ $$

The following equations are immediate from definition.

###### Claim 9.5 .

Let $H_{w,x}\in\mathbb{R}^{n_{\mathrm{tot}}\times n_{\mathrm{tot}}}$ be defined as Definition [9.4](#S9.Thmtheorem4).
For $v\in\mathbb{R}^{n_{\mathrm{tot}}}$, we have

$$ $\displaystyle\|v\|_{w,x}$ $\displaystyle=\|H_{w,x}^{1/2}v\|_{2},$ $\displaystyle\|v\|_{w,x}^{*}$ $\displaystyle=\|H_{w,x}^{-1/2}v\|_{2}.$ $$

###### Claim 9.6 .

For each $i\in[n]$, let $H_{x,i}$ be defined as Definition [9.4](#S9.Thmtheorem4).
For $v\in\mathbb{R}^{n_{i}}$, $i\in[n]$, we have

$$ $\displaystyle\|v\|_{x_{i}}$ $\displaystyle=\|H_{x,i}^{1/2}v\|_{2},$ $\displaystyle\|v\|_{x_{i}}^{*}$ $\displaystyle=\|H_{x,i}^{-1/2}v\|_{2}.$ $$

We define matrices $B$ and $P$ used in the algorithm.

###### Definition 9.7 .

Let $A,Q$ denote two fixed matrices.
Let $H_{w,x}\in\mathbb{R}^{n_{\mathrm{tot}}\times n_{\mathrm{tot}}}$ be defined as Definition [9.4](#S9.Thmtheorem4).
We define matrix $B_{w,x,t}\in\mathbb{R}^{n_{\mathrm{tot}}\times n_{\mathrm{tot}}}$ as

$$ $\displaystyle B_{w,x,t}:=Q+t\cdot H_{w,x}$ (25) $$

We define projection matrix $P_{w,x,t}\in\mathbb{R}^{n_{\mathrm{tot}}\times n_{\mathrm{tot}}}$ as

$$ $\displaystyle P_{w,x,t}\leftarrow B_{w,x,t}^{-1/2}A^{\top}(AB_{w,x,t}^{-1}A^{\top})^{-1}AB_{w,x,t}^{-1/2}.$ (26) $$

###### Definition 9.2 .

###### Definition 9.3 .

###### Definition 9.4 .

###### Claim 9.5 .

###### Claim 9.6 .

###### Definition 9.7 .

### 9.2 Deriving the Central Path Step

In this section we explain how to derive the central path step.

We follow the central path

$$ $\displaystyle s/t+\nabla\phi_{w}(x)$ $\displaystyle=\mu$ $\displaystyle Ax$ $\displaystyle=b$ $\displaystyle-Qx+A^{\top}y+s$ $\displaystyle=c$ $$

We perform gradient descent on $\mu$ with step $\delta_{\mu}$. Then Newton step gives

$$ $\displaystyle\frac{1}{t}\delta_{s}+\nabla^{2}\phi_{w}(x)\delta_{x}=$ $\displaystyle~{}\delta_{\mu}$ (27) $\displaystyle A\delta_{x}=$ $\displaystyle~{}0$ (28) $\displaystyle-Q\delta_{x}+A^{\top}\delta_{y}+\delta_{s}=$ $\displaystyle~{}0$ (29) $$

where $\delta_{x}$ (resp. $\delta_{y}$, $\delta_{s}$) is the step taken by $x$ (resp. $y$, $s$).

For simplicity, we define $H\in\mathbb{R}^{n_{\mathrm{tot}}\times n_{\mathrm{tot}}}$ to represent $\nabla^{2}\phi_{w}(x)$.(^7^77In this section, and in this section only, we omit the subscript in $H$, $B$, $P$ for simplicity.)

From Eq. ([27](#S9.E27)) we get

$$ $\displaystyle\delta_{s}=t\delta_{\mu}-tH\delta_{x}.$ (30) $$

Plug the above equation into Eq. ([29](#S9.E29)) we get

$$ $\displaystyle-Q\delta_{x}+A^{\top}\delta_{y}+t\delta_{\mu}-tH\delta_{x}=0.$ (31) $$

Let $B=Q+tH$,
multiply by $AB^{-1}$ we get

$$ $\displaystyle-A\delta_{x}+AB^{-1}A^{\top}\delta_{y}+tAB^{-1}\delta_{\mu}=0.$ $$

Using Eq. ([28](#S9.E28)) we get

$$ $\displaystyle AB^{-1}A^{\top}\delta_{y}+tAB^{-1}\delta_{\mu}=0.$ $$

Solve for $\delta_{y}$ (assuming that $AB^{-1}A$ is invertible), we get

$$ $\displaystyle\delta_{y}=-t(AB^{-1}A^{\top})^{-1}AB^{-1}\delta_{\mu}.$ $$

Plug into Eq. ([31](#S9.E31)) we get

$$ $\displaystyle-B\delta_{x}-tA^{\top}(AB^{-1}A^{\top})^{-1}AB^{-1}\delta_{\mu}+t\delta_{\mu}=0.$ $$

Solve for $\delta_{x}$ we get

$$ $\displaystyle\delta_{x}$ $\displaystyle=tB^{-1}\delta_{\mu}-tB^{-1}A^{\top}(AB^{-1}A^{\top})^{-1}AB^{-1}\delta_{\mu}$ $\displaystyle=tB^{-1/2}(I-P)B^{-1/2}\delta_{\mu}$ $$

where $P=B^{-1/2}A^{\top}(AB^{-1}A^{\top})^{-1}AB^{-1/2}$ is the projection matrix.
Solve for $\delta_{s}$ in Eq. ([30](#S9.E30)) we get

$$ $\displaystyle\delta_{s}=t\delta_{\mu}-t^{2}HB^{-1/2}(I-P)B^{-1/2}\delta_{\mu}.$ $$

In summary, we have

$$ $\displaystyle\delta_{x}$ $\displaystyle=tB^{-1/2}(I-P)B^{-1/2}\delta_{\mu},$ $\displaystyle\delta_{y}$ $\displaystyle=-t(AB^{-1}A^{\top})^{-1}AB^{-1}\delta_{\mu},$ $\displaystyle\delta_{s}$ $\displaystyle=t\delta_{\mu}-t^{2}HB^{-1/2}(I-P)B^{-1/2}\delta_{\mu},$ $\displaystyle P$ $\displaystyle=B^{-1/2}A^{\top}(AB^{-1}A^{\top})^{-1}AB^{-1/2}.$ $$

These equations will guide the design of our actual algorithm.

### 9.3 Bounding Movement of Potential Function

The goal of this section is to bound the movement of potential function during the robust IPM algorithm.

In robust IPM, we do not need to follow the ideal central path exactly over the entire algorithm. Instead, we only use an approximate version.
For convenience of analysis we state two assumptions (see Algorithm [20](#alg20), Line [18](#alg20.l18)).

###### Assumption 9.8 .

We make the following assumptions on $\delta_{x}\in\mathbb{R}^{n_{\mathrm{tot}}}$ and $\delta_{s}\in\mathbb{R}^{n_{\mathrm{tot}}}$.

$$ $\displaystyle\|\delta_{x}-\overline{t}B_{w,\overline{x},\overline{t}}^{-1/2}(I-P_{w,\overline{x},\overline{t}})B_{w,\overline{x},\overline{t}}^{-1/2}\delta_{\mu}\|_{w,\overline{x}}$ $\displaystyle\leq\overline{\epsilon}\alpha,$ $\displaystyle\|\overline{t}^{-1}\delta_{s}-(\delta_{\mu}-\overline{t}H_{w,\overline{x}}B_{w,\overline{x},\overline{t}}^{-1/2}(I-P_{w,\overline{x},\overline{t}})B_{w,\overline{x},\overline{t}}^{-1/2}\delta_{\mu})\|_{w,\overline{x}}^{*}$ $\displaystyle\leq\overline{\epsilon}\alpha.$ $$

The following lemma bounds the movement of potential function $\Psi$ assuming bound on $\delta_{\gamma}$.

###### Lemma 9.9 ( [ 59 , Lemma A.5] ) .

For any $r\in\mathbb{R}^{n_{\mathrm{tot}}}$, and $w\in\mathbb{R}^{n_{\mathrm{tot}}}_{\geq 1}$. Let $\alpha$ and $\lambda$ denote the parameters that are satisfying $0\leq\alpha\leq\frac{1}{8\lambda}$.

Let $\epsilon_{r}\in\mathbb{R}^{n_{\mathrm{tot}}}$ denote a vector satisfying

$$ $\displaystyle(\sum_{i=1}^{n}w_{i}^{-1}\epsilon_{r,i}^{2})^{1/2}\leq\alpha/8.$ $$

Suppose that vector $\overline{r}\in\mathbb{R}^{n_{\mathrm{tot}}}$ is satisfying the following property

$$ $\displaystyle|r_{i}-\overline{r}_{i}|\leq\frac{w_{i}}{8\lambda},~{}~{}~{}\forall i\in[n]$ $$

We define vector $\delta_{r}\in\mathbb{R}^{n_{\mathrm{tot}}}$ as follows:

$$ $\displaystyle\delta_{r,i}:=\frac{-\alpha\cdot\sinh(\frac{\lambda}{w_{i}}\overline{r}_{i})}{\sqrt{\sum_{j=1}^{n}w_{j}^{-1}\cosh^{2}(\frac{\lambda}{w_{j}}\overline{r}_{j})}}+\epsilon_{r,i}.$ $$

Then, we have that

$$ $\displaystyle\Psi_{\lambda}(r+\delta_{r})\leq\Psi_{\lambda}(r)-\frac{\alpha\lambda}{2}(\sum_{i=1}^{n}w_{i}^{-1}\cosh^{2}(\lambda\frac{r_{i}}{w_{i}}))^{1/2}+\alpha\lambda(\sum_{i=1}^{n}w_{i}^{-1})^{1/2}$ $$

The following lemma bounds the norm of $\delta_{\mu}.$

###### Lemma 9.10 (Bounding norm of $\delta_{\mu}$ ) .

$$ $\displaystyle\|\delta_{\mu}(\overline{x},\overline{s},\overline{t})\|_{w,\overline{x}}^{*}\leq\alpha.$ $$

###### Proof.

$$ $\displaystyle(\|\delta_{\mu}(\overline{x},\overline{s},\overline{t})\|_{w,\overline{x}}^{*})^{2}=$ $\displaystyle~{}\sum_{i=1}^{n}w_{i}^{-1}(\|\delta_{\mu,i}(\overline{x},\overline{s},\overline{t})\|_{\overline{x}_{i}}^{*})^{2}$ $\displaystyle=$ $\displaystyle~{}\alpha^{2}\sum_{i\in[n]}w_{i}^{-1}c_{i}^{2}(\overline{x},\overline{s},\overline{t})\cdot\|\mu_{i}(\overline{x},\overline{s},\overline{t})\|_{\overline{x}_{i}}^{2}$ $\displaystyle=$ $\displaystyle~{}\alpha^{2}\sum_{i\in[n]}w_{i}^{-1}c_{i}^{2}(\overline{x},\overline{s},\overline{t})\cdot\|H_{\overline{x},i}^{-1/2}\mu_{i}(\overline{x},\overline{s},\overline{t})\|_{2}^{2}$ $\displaystyle=$ $\displaystyle~{}\alpha^{2}\sum_{i\in[n]}w_{i}^{-1}c_{i}^{2}(\overline{x},\overline{s},\overline{t})\cdot\gamma_{i}^{2}(\overline{x},\overline{s},\overline{t})$ $\displaystyle=$ $\displaystyle~{}\alpha^{2}\sum_{i\in[n]}\frac{w_{i}^{-1}\sinh^{2}(\frac{\lambda}{w_{i}}\gamma_{i}(\overline{x},\overline{s},\overline{t}))}{\gamma_{i}^{2}(\overline{x},\overline{s},\overline{t})\cdot{\sum_{j\in[n]}w_{j}^{-1}\cosh^{2}(\frac{\lambda}{w_{j}}\gamma_{j}(\overline{x},\overline{s},\overline{t}))}}\cdot\gamma_{i}^{2}(\overline{x},\overline{s},\overline{t})$ $\displaystyle=$ $\displaystyle~{}\alpha^{2}\frac{\sum_{j\in[n]}w_{j}^{-1}\sinh^{2}(\frac{\lambda}{w_{j}}\gamma_{j}(\overline{x},\overline{s},\overline{t}))}{\sum_{j\in[n]}w_{j}^{-1}\cosh^{2}(\frac{\lambda}{w_{j}}\gamma_{j}(\overline{x},\overline{s},\overline{t}))}$ $\displaystyle\leq$ $\displaystyle~{}\alpha^{2}.$ $$

where the first step follows from Definition [9.3](#S9.Thmtheorem3), the second step follows from $\delta_{\mu,i}(\overline{x},\overline{s},\overline{t})=-\alpha\cdot c_{i}(\overline{x},\overline{s},\overline{t})\cdot\mu_{i}(\overline{x},\overline{s},\overline{t})$, the third step follows from norm of $\overline{x}_{i}$ (see Definition [9.3](#S9.Thmtheorem3)), the forth step follows from $\gamma_{i}(\overline{x},\overline{s},\overline{t})=\|H_{\overline{x},i}^{-1/2}\mu_{i}(\overline{x},\overline{s},\overline{t})\|_{2}$ (see Eq. ([18](#S9.E18))), the fifth step follows from $c_{i}(\overline{x},\overline{s},\overline{t})^{2}=\frac{\sinh^{2}(\frac{\lambda}{w_{i}}\gamma_{i}(\overline{x},\overline{s},\overline{t}))}{\gamma_{i}^{2}(\overline{x},\overline{s},\overline{t})\sum_{j\in[n]}w_{j}^{-1}\cosh^{2}(\frac{\lambda}{w_{j}}\gamma_{j}(\overline{x},\overline{s},\overline{t}))}$ (see Eq. ([22](#S9.E22))), the sixth step follows from canceling the term $\gamma_{i}^{2}(\overline{x},\overline{s},\overline{t})$, and the last step follows from $\cosh^{2}(x)\geq\sinh^{2}(x)$ for all $x$.

∎

The following lemma bounds the norm of $\delta_{x}$ and $\delta_{s}$.

###### Lemma 9.11 .

For each $i\in[n]$, we define
$\alpha_{i}:=\|\delta_{x,i}\|_{\overline{x}_{i}}.$
Then, we have

$$ $\displaystyle\|\delta_{x}\|_{w,\overline{x}}=(\sum_{i\in[n]}w_{i}\alpha_{i}^{2})^{1/2}\leq\frac{9}{8}\alpha.$ (32) $$

In particular, we have
$\alpha_{i}\leq\frac{9}{8}\alpha.$
Similarly, for $\delta_{s}$, we have

$$ $\displaystyle\|\delta_{s}\|_{w,\overline{x}}^{*}=\sqrt{\sum_{i\in[n]}w_{i}^{-1}(\|\delta_{s,i}\|_{\overline{x}_{i}}^{*})^{2}}\leq\frac{17}{8}\alpha\cdot t.$ (33) $$

###### Proof.

For $\delta_{x}$, we have

$$ $\displaystyle\|\delta_{x}\|_{w,\overline{x}}$ $\displaystyle\leq\|\overline{t}H_{w,\overline{x}}^{1/2}B_{w,\overline{x},\overline{t}}^{-1/2}(I-P_{w,\overline{x},\overline{t}})B_{w,\overline{x},\overline{t}}^{-1/2}\delta_{\mu}\|_{2}+\overline{\epsilon}\alpha$ $\displaystyle\leq\|\overline{t}^{1/2}(I-P_{w,\overline{x},\overline{t}})B_{w,\overline{x},\overline{t}}^{-1/2}\delta_{\mu}\|_{2}+\overline{\epsilon}\alpha$ $\displaystyle\leq\|\overline{t}^{1/2}B_{w,\overline{x},\overline{t}}^{-1/2}\delta_{\mu}\|_{2}+\overline{\epsilon}\alpha$ $\displaystyle\leq\|H_{w,\overline{x}}^{-1/2}\delta_{\mu}\|_{2}+\overline{\epsilon}\alpha$ $\displaystyle\leq\alpha+\overline{\epsilon}\alpha$ $\displaystyle\leq\frac{9}{8}\alpha.$ $$

First step follows from Assumption [9.8](#S9.Thmtheorem8).
Second step is because $\overline{t}H_{w,\overline{x}}\preceq B_{w,\overline{x},\overline{t}}$.
Third step is because $P_{w,\overline{x},\overline{t}}$ is a projection matrix.
Fourth step is because $\overline{t}H_{w,\overline{x}}\preceq B_{w,\overline{x},\overline{t}}$.
Fifth step is by Lemma [9.10](#S9.Thmtheorem10).
Sixth step is because $\overline{\epsilon}\leq\frac{1}{8}$.

For $\delta_{s}$, we have

$$ $\displaystyle\|\delta_{s}\|_{w,\overline{x}}^{*}$ $\displaystyle\leq\|\overline{t}\delta_{\mu}\|_{w,\overline{x}}^{*}+\|\overline{t}^{2}H_{w,\overline{x}}B_{w,\overline{x},\overline{t}}^{-1/2}(I-P_{w,\overline{x},\overline{t}})B_{w,\overline{x},\overline{t}}^{-1/2}\delta_{\mu}\|_{w,\overline{x}}^{*}+\overline{\epsilon}\alpha\overline{t}$ $\displaystyle\leq\alpha\overline{t}+\alpha\overline{t}+\overline{\epsilon}\alpha\overline{t}$ $\displaystyle\leq\frac{17}{8}\alpha\cdot t.$ $$

First step is by triangle inequality and the assumption that

$$ $\displaystyle\delta_{s}\approx\overline{t}\delta_{\mu}-\overline{t}^{2}H_{w,\overline{x}}B_{w,\overline{x},\overline{t}}^{-1/2}(I-P_{w,\overline{x},\overline{t}})B_{w,\overline{x},\overline{t}}^{-1/2}\delta_{\mu}.$ $$

Second step is by same analysis as the analysis for $\delta_{x}$.
Third step is by $\overline{t}\leq\frac{33}{32}t$ and $\overline{\epsilon}\leq\frac{1}{32}$.
∎

The following lemma shows that $\mu^{\mathrm{new}}$ is close to $\mu+\delta_{\mu}$ under an approximate step.

###### Lemma 9.12 (Variation of [ 59 , Lemma A.9] ) .

For each $i\in[n]$, we define

$$ $\displaystyle\beta_{i}:=\|\epsilon_{\mu,i}\|_{x_{i}}^{*}$ $$

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

$$ $\displaystyle\mu_{i}(x^{\mathrm{new}},s^{\mathrm{new}},t)=\mu_{i}(x,s,t)+\delta_{\mu,i}+\epsilon_{\mu,i}.$ $$

Then, we have

$$ $\displaystyle(\sum_{i=1}^{n}w_{i}^{-1}\beta_{i}^{2})^{1/2}\leq 15\overline{\epsilon}\alpha.$ $$

###### Proof.

The proof is similar as [, Lemma A.9], except for changing the definitions of $\epsilon_{1}$ and $\epsilon_{2}$:

$$ $\displaystyle\epsilon_{1}$ $\displaystyle:=H_{w,\overline{x}}^{1/2}\delta_{x}-\overline{t}\cdot H_{w,\overline{x}}^{1/2}B_{w,\overline{x},\overline{t}}^{-1/2}(I-P_{w,\overline{x},\overline{t}})B_{w,\overline{x},\overline{t}}^{-1/2}\delta_{\mu},$ $\displaystyle\epsilon_{2}$ $\displaystyle:=\overline{t}^{-1}H_{w,\overline{x}}^{-1/2}\delta_{s}-H_{w,\overline{x}}^{-1/2}(\delta_{\mu}-\overline{t}H_{w,\overline{x}}B_{w,\overline{x},\overline{t}}^{-1/2}(I-P_{w,\overline{x},\overline{t}})B_{w,\overline{x},\overline{t}}^{-1/2}\delta_{\mu}).$ $$

One key step in the proof of [] is the following property:

$$ $\displaystyle\delta_{\mu,i}=\overline{t}^{-1}\cdot\delta_{s,i}+H_{w,\overline{x}}\delta_{x,i}-H_{w,\overline{x}}^{1/2}(\epsilon_{1}+\epsilon_{2}).$ $$

Under our new definition of $\epsilon_{1}$ and $\epsilon_{2}$, the above property still holds.
Remaining parts of the proof are similar and we omit the details here.
∎

The following lemma shows that error $\mu(\overline{x},\overline{s},\overline{t})$ on the robust central path is close to error $\mu(x,s,t)$ on the ideal central path. Furthermore, norms of errors $\gamma_{i}(x,s,t)$ and $\gamma_{i}(\overline{x},\overline{s},\overline{t})$ are also close to each other.

###### Lemma 9.13 ( [ 59 , Lemma A.10] ) .

Assume that $\gamma_{i}(x,s,t)\leq w_{i}$ for all $i$. For all $i\in[n]$, we have

$$ $\displaystyle\|\mu_{i}(x,s,t)-\mu_{i}(\overline{x},\overline{s},\overline{t})\|_{x_{i}}^{*}\leq 3\overline{\epsilon}w_{i}.$ $$

Furthermore, we have that

$$ $\displaystyle|\gamma_{i}(x,s,t)-\gamma_{i}(\overline{x},\overline{s},\overline{t})|\leq 5\overline{\epsilon}w_{i}.$ $$

###### Proof.

Same as proof of [, Lemma A.10].
∎

The following lemma bounds the change of $\gamma$ under one robust IPM step.

###### Lemma 9.14 ( [ 59 , Lemma A.12] ) .

Assume $\Phi(x,s,t)\leq\cosh(\lambda)$. For all $i\in[n]$, we define

$$ $\displaystyle\epsilon_{r,i}:=\gamma_{i}(x^{\mathrm{new}},s^{\mathrm{new}})-\gamma_{i}(x,s,t)+\alpha\cdot c_{i}(\overline{x},\overline{s},\overline{t})\cdot\gamma_{i}(\overline{x},\overline{s},\overline{t}).$ $$

Then, we have

$$ $\displaystyle(\sum_{i=1}^{n}w_{i}^{-1}\epsilon_{r,i}^{2})^{1/2}\leq 90\cdot\overline{\epsilon}\cdot\lambda\alpha+4\cdot\max_{i\in[n]}(w_{i}^{-1}\gamma_{i}(x,s,t))\cdot\alpha.$ $$

###### Proof.

The proof is similar to the proof of [, Lemma A.12]. By replacing corresponding references in [] by our versions (Lemma [9.11](#S9.Thmtheorem11), [9.12](#S9.Thmtheorem12), [9.13](#S9.Thmtheorem13)) we get proof of this lemma.
∎

Finally, the following theorem bounds the movement of potential function $\Phi$ under one robust IPM step.

###### Theorem 9.15 (Variation of [ 59 , Theorem A.15] ) .

Assume $\Phi(x,s,t)\leq\cosh(\lambda/64)$. Then for any $0\leq h\leq\frac{\alpha}{64\sqrt{\sum_{i\in[n]}w_{i}\nu_{i}}}$, we have

$$ $\displaystyle\Phi(x^{\mathrm{new}},s^{\mathrm{new}},t^{\mathrm{new}})\leq(1-\frac{\alpha\lambda}{\sqrt{\sum_{i\in[n]}w_{i}}})\cdot\Phi(x,s,t)+\alpha\lambda\sqrt{\sum_{i\in[n]}w_{i}^{-1}}.$ $$

In particular, for any $\cosh(\lambda/128)\leq\Phi(x,s,t)\leq\cosh(\lambda/64)$, we have that

$$ $\displaystyle\Phi(x^{\mathrm{new}},s^{\mathrm{new}},t^{\mathrm{new}})\leq\Phi(x,s,t).$ $$

###### Proof.

Similar to the proof of [, Theorem A.15], but replacing lemmas with the corresponding QP versions.
∎

###### Assumption 9.8 .

###### Lemma 9.9 ( [ 59 , Lemma A.5] ) .

###### Lemma 9.10 (Bounding norm of δ μ subscript 𝛿 𝜇 \delta_{\mu} ) .

###### Proof.

###### Lemma 9.11 .

###### Proof.

###### Lemma 9.12 (Variation of [ 59 , Lemma A.9] ) .

###### Proof.

###### Lemma 9.13 ( [ 59 , Lemma A.10] ) .

###### Proof.

###### Lemma 9.14 ( [ 59 , Lemma A.12] ) .

###### Proof.

###### Theorem 9.15 (Variation of [ 59 , Theorem A.15] ) .

###### Proof.

### 9.4 Initial Point Reduction

In this section, we propose an initial point reduction scheme for quadratic programming.
Our scheme is closer to [] rather than [, ]. The reason is that []’s initial point reduction requires an efficient algorithm for finding the optimal solution to an unconstrained program, which may be difficult in quadratic programming.

###### Lemma 9.16 ( [ 44 , Theorem 4.1.7 and Lemma 4.2.4] ) .

Let $\phi$ be a $\nu$-self-concordant barrier. Then for any $x,y\in\mathrm{dom}(\phi)$, we have

$$ $\displaystyle\langle\nabla\phi(x),y-x\rangle\leq$ $\displaystyle~{}\nu,$ $\displaystyle\langle\nabla\phi(y)-\nabla\phi(x),y-x\rangle\geq$ $\displaystyle~{}\frac{\|y-x\|_{x}^{2}}{1+\|y-x\|_{x}}.$ $$

Let $x^{*}=\arg\min_{x}\phi(x)$. For any $x\in\mathbb{R}^{n}$ such that $\|x-x^{*}\|_{x^{*}}\leq 1$, we have that $x\in\mathrm{dom}(\phi)$.

###### Lemma 9.17 (QP version of [ 40 , Lemma D.2] ) .

Work under the setting of Theorem [9.1](#S9.Thmtheorem1).
Let $x^{(0)}=\arg\min_{x}\sum_{i\in[n]}w_{i}\phi_{i}(x_{i})$.
Let $\rho=\frac{1}{LR(R+1)}$.
For any $0<\epsilon\leq\frac{1}{2}$, the modified program

$$ $\displaystyle\min_{\overline{A}\overline{x}=\overline{b},\overline{x}\in\mathcal{K}\times\mathbb{R}_{\geq 0}}\left(\frac{1}{2}\overline{x}^{\top}\overline{Q}\overline{x}+\overline{c}^{\top}\overline{x}\right)$ $$

with

$$ $\displaystyle\overline{Q}=\begin{bmatrix}\epsilon\rho Q&0\\ 0&0\end{bmatrix},\qquad\overline{A}=[A~{}|~{}b-Ax^{(0)}],\qquad\overline{b}=b,\qquad\overline{c}=\begin{bmatrix}\epsilon\rho c\\ 1\end{bmatrix}$ $$

satisfies the following:

- •
$\overline{x}=\begin{bmatrix}x^{(0)}\\
1\end{bmatrix}$, $\overline{y}=0\in\mathbb{R}^{m}$ and $\overline{s}=\begin{bmatrix}\epsilon\rho(c+Qx^{(0)})\\
1\end{bmatrix}$ are feasible primal dual vectors with $\|\overline{s}+\nabla\overline{\phi}_{w}(\overline{x})\|_{\overline{x}}^{*}\leq\epsilon$ where $\overline{\phi}_{w}(\overline{x})=\sum_{i=1}^{n}w_{i}\phi_{i}(\overline{x}_{i})-\log(\overline{x}_{n+1})$.
- •
For any $\overline{x}\in\mathcal{K}\times\mathbb{R}_{\geq 0}$ satisfying $\overline{A}\overline{x}=\overline{b}$ and
$\displaystyle\frac{1}{2}\overline{x}^{\top}\overline{Q}\overline{x}+\overline{c}^{\top}\overline{x}\leq\min_{\overline{A}\overline{x}=\overline{b},\overline{x}\in\mathcal{K}\times\mathbb{R}_{\geq 0}}\left(\frac{1}{2}\overline{x}^{\top}\overline{Q}\overline{x}+\overline{c}^{\top}\overline{x}\right)+\epsilon^{2},$
(34)
the vector $\overline{x}_{1:n}$ ( $\overline{x}_{1:n}$ is the first $n$ coordinates of $\overline{x}$ ) is an approximate solution to the original convex program in the following sense:
$\displaystyle\frac{1}{2}\overline{x}_{1:n}^{\top}Q\overline{x}_{1:n}+c^{\top}\overline{x}_{1:n}$
$\displaystyle\leq\min_{Ax=b,x\in\mathcal{K}}\left(\frac{1}{2}x^{\top}Qx+c^{\top}x\right)+\epsilon\rho^{-1},$
$\displaystyle\|A\overline{x}_{1:n}-b\|_{1}$
$\displaystyle\leq 3\epsilon\cdot(R\|A\|_{1}+\|b\|_{1}),$
$\displaystyle\overline{x}_{1:n}$
$\displaystyle\in\mathcal{K}.$

###### Proof.

First bullet point:
Direct computation shows that $(\overline{x},\overline{y},\overline{s})$ is feasible.

Let us compute $\|\overline{s}+\nabla\overline{\phi}_{w}(\overline{x})\|_{\overline{x}}^{*}$.
We have

$$ $\displaystyle\|\overline{s}+\nabla\overline{\phi}_{w}(\overline{x})\|_{\overline{x}}^{*}=\|\epsilon\rho(c+Qx^{(0)})\|_{\nabla^{2}\phi_{w}(x^{(0)})^{-1}}$ $$

Lemma [9.16](#S9.Thmtheorem16) says that for all $x\in\mathbb{R}^{n}$ with $\|x-x^{(0)}\|_{w,x^{(0)}}\leq 1$, we have $x\in\mathcal{K}$, because $x^{(0)}=\arg\min_{x}\phi_{w}(x)$.
Therefore for any $v$ such that $v^{\top}\nabla^{2}\phi_{w}(x^{(0)})v\leq 1$, we have $x^{(0)}\pm v\in\mathcal{K}$ and hence $\|x^{(0)}\pm v\|_{2}\leq R$. This implies $\|v\|_{2}\leq R$ for any $v^{\top}\nabla^{2}\phi_{w}(x^{(0)})v\leq 1$. Hence $(\nabla^{2}\phi_{w}(x^{(0)}))^{-1}\preceq R^{2}\cdot I$. So we have

$$ $\displaystyle\|\overline{s}+\nabla\overline{\phi}_{w}(\overline{x})\|_{\overline{x}}^{*}$ $\displaystyle=\|\epsilon\rho(c+Qx^{(0)})\|_{\nabla^{2}\phi_{w}(x^{(0)})^{-1}}$ $\displaystyle\leq\epsilon\rho R\|c+Qx^{(0)}\|_{2}$ $\displaystyle\leq\epsilon\rho R(\|c\|_{2}+\|Q\|_{2\to 2}\|x^{(0)}\|_{2})$ $\displaystyle\leq\epsilon\rho R(L+LR)$ $\displaystyle\leq\epsilon.$ $$

Second bullet point:
We define

$$ $\displaystyle\operatorname{OPT}$ $\displaystyle:=\min_{Ax=b,x\in\mathcal{K}}\left(\frac{1}{2}x^{\top}Qx+c^{\top}x\right),$ (35) $\displaystyle\overline{\operatorname{OPT}}$ $\displaystyle:=\min_{\overline{A}\overline{x}=\overline{b},\overline{x}\in\mathcal{K}\times\mathbb{R}_{\geq 0}}\left(\frac{1}{2}\overline{x}^{\top}\overline{Q}\overline{x}+\overline{c}^{\top}\overline{x}\right).$ (36) $$

For any feasible $x$ in the original problem ([35](#S9.E35)), $\overline{x}=\begin{bmatrix}x\\
0\end{bmatrix}$ is feasible in the modified problem ([36](#S9.E36)).
Therefore we have

$$ $\displaystyle\overline{\operatorname{OPT}}\leq\epsilon\rho(\frac{1}{2}x^{\top}Qx+c^{\top}x)=\epsilon\rho\cdot\operatorname{OPT}.$ $$

Given a feasible $\overline{x}$ satisfying ([34](#S9.E34)), we write $\overline{x}=\begin{bmatrix}\overline{x}_{1:n}\\
\tau\end{bmatrix}$ for some $\tau\geq 0$.
Then we have

$$ $\displaystyle\epsilon\rho(\frac{1}{2}\overline{x}_{1:n}^{\top}Q\overline{x}_{1:n}+c^{\top}\overline{x}_{1:n})+\tau\leq\overline{\operatorname{OPT}}+\epsilon^{2}\leq\epsilon\rho\cdot\operatorname{OPT}+\epsilon^{2}.$ $$

Therefore

$$ $\displaystyle\frac{1}{2}\overline{x}_{1:n}^{\top}Q\overline{x}_{1:n}+c^{\top}\overline{x}_{1:n}\leq\operatorname{OPT}+\epsilon\rho^{-1}.$ $$

We have

$$ $\displaystyle\tau\leq-\epsilon\rho(\frac{1}{2}\overline{x}_{1:n}^{\top}Q\overline{x}_{1:n}+c^{\top}\overline{x}_{1:n})+\epsilon\rho\cdot\operatorname{OPT}+\epsilon^{2}\leq 3\epsilon$ $$

because
$\left|\frac{1}{2}x^{\top}Qx+c^{\top}x\right|\leq LR(R+1)$
for all $x\in\mathcal{K}$.

Note that $\overline{x}$ satisfies $A\overline{x}_{1:n}+(b-Ax^{(0)})\tau=b$.
So

$$ $\displaystyle\|A\overline{x}_{1:n}-b\|_{1}\leq\|b-Ax^{(0)}\|_{1}\cdot\tau.$ $$

This finishes the proof.
∎

The following lemma is a generalization of [, Lemma D.3] to quadratic program, and with weight vector $w$.

###### Lemma 9.18 (QP version of [ 40 , Lemma D.3] ) .

Work under the setting of Theorem [9.1](#S9.Thmtheorem1).
Suppose we have $\frac{s_{i}}{t}+w_{i}\nabla\phi_{i}(x_{i})=\mu_{i}$ for all $i\in[n]$, $-Qx+A^{\top}y+s=c$ and $Ax=b$.
If $\|\mu_{i}\|_{x_{i}}^{*}\leq w_{i}$ for all $i\in[n]$, then we have

$$ $\displaystyle\frac{1}{2}x^{\top}Qx+c^{\top}x\leq\frac{1}{2}x^{*\top}Qx^{*}+c^{\top}x^{*}+4t\kappa$ $$

where $x^{*}=\arg\min_{Ax=b,x\in\mathcal{K}}\left(\frac{1}{2}x^{\top}Qx+c^{\top}x\right)$.

###### Proof.

Let $x_{\alpha}=(1-\alpha)x+\alpha x^{*}$ for some $\alpha$ to be chosen.
By Lemma [9.16](#S9.Thmtheorem16), we have $\langle\nabla\phi_{w}(x_{\alpha}),x^{*}-x_{\alpha}\rangle\leq\kappa$.
(Note that $\phi_{w}$ is a $\kappa$-self-concordant barrier for $\mathcal{K}$.)
Therefore we have

$$ $\displaystyle\frac{\kappa\alpha}{1-\alpha}$ $\displaystyle\geq\langle\nabla\phi_{w}(x_{\alpha}),x_{\alpha}-x\rangle$ $\displaystyle=\langle\nabla\phi_{w}(x_{\alpha})-\nabla\phi_{w}(x),x_{\alpha}-x\rangle+\langle\mu-\frac{s}{t},x_{\alpha}-x\rangle$ $\displaystyle\geq\sum_{i\in[n]}w_{i}\frac{\|x_{\alpha,i}-x_{i}\|_{x_{i}}^{2}}{1+\|x_{\alpha,i}-x_{i}\|_{x_{i}}}+\langle\mu,x_{\alpha}-x\rangle-\frac{1}{t}\langle c-A^{\top}y+Qx,x_{\alpha}-x\rangle$ $\displaystyle\geq\sum_{i\in[n]}w_{i}\frac{\alpha^{2}\|x_{i}^{*}-x_{i}\|_{x_{i}}^{2}}{1+\alpha\|x_{i}^{*}-x_{i}\|_{x_{i}}}-\alpha\sum_{i\in[n]}\|\mu_{i}\|_{x_{i}}^{*}\|x_{i}^{*}-x_{i}\|_{x_{i}}-\frac{\alpha}{t}\langle c+Qx,x^{*}-x\rangle.$ $$

First step is because $\langle\nabla\phi_{w}(x_{\alpha}),x^{*}-x_{\alpha}\rangle\leq\nu$.
Second step is because $\mu=\frac{s}{t}+\nabla\phi_{w}(x)$.
Third step is by Lemma [9.16](#S9.Thmtheorem16) and $c=-Qx+A^{\top}y+s$.
Fourth step is by Cauchy-Schwartz and $Ax_{\alpha}=Ax$.

So we get

$$ $\displaystyle~{}\frac{1}{t}(x^{\top}Qx+c^{\top}x)$ $\displaystyle\leq$ $\displaystyle~{}\frac{1}{t}(x^{\top}Qx^{*}+c^{\top}x^{*})+\frac{\kappa}{1-\alpha}+\sum_{i\in[n]}\|\mu_{i}\|_{x_{i}}^{*}\|x_{i}^{*}-x_{i}\|_{x_{i}}-\sum_{i\in[n]}w_{i}\frac{\alpha\|x_{i}^{*}-x_{i}\|_{x_{i}}^{2}}{1+\alpha\|x_{i}^{*}-x_{i}\|_{x_{i}}}$ $\displaystyle\leq$ $\displaystyle~{}\frac{1}{t}(\frac{1}{2}x^{\top}Qx+\frac{1}{2}x^{*\top}Qx^{*}+c^{\top}x^{*})+\frac{\kappa}{1-\alpha}+\sum_{i\in[n]}w_{i}\|x_{i}^{*}-x_{i}\|_{x_{i}}-\sum_{i\in[n]}w_{i}\frac{\alpha\|x_{i}^{*}-x_{i}\|_{x_{i}}^{2}}{1+\alpha\|x_{i}^{*}-x_{i}\|_{x_{i}}}$ $\displaystyle=$ $\displaystyle~{}\frac{1}{t}(\frac{1}{2}x^{\top}Qx+\frac{1}{2}x^{*\top}Qx^{*}+c^{\top}x^{*})+\frac{\kappa}{1-\alpha}+\sum_{i\in[n]}w_{i}\frac{\|x_{i}^{*}-x_{i}\|_{x_{i}}}{1+\alpha\|x_{i}^{*}-x_{i}\|_{x_{i}}}$ $\displaystyle\leq$ $\displaystyle~{}\frac{1}{t}(\frac{1}{2}x^{\top}Qx+\frac{1}{2}x^{*\top}Qx^{*}+c^{\top}x^{*})+\frac{\kappa}{1-\alpha}+\sum_{i\in[n]}\frac{w_{i}}{\alpha}$ $\displaystyle\leq$ $\displaystyle~{}\frac{1}{t}(\frac{1}{2}x^{\top}Qx+\frac{1}{2}x^{*\top}Qx^{*}+c^{\top}x^{*})+\frac{\kappa}{\alpha(1-\alpha)}.$ $$

First step is by rearranging terms in the previous inequality.
Second step is by AM-GM inequality and $\|\mu_{i}\|_{x_{i}}^{*}\leq w_{i}$.
Third step is by merging the last two terms.
Fourth step is by bounding the last term.
Fifth step is by $\sum_{i\in[n]}w_{i}\leq\sum_{i\in[n]}w_{i}\nu_{i}=\kappa$.

Finally,

$$ $\displaystyle\frac{1}{2}x^{\top}Qx+c^{\top}x$ $\displaystyle\leq\frac{1}{2}x^{*\top}Qx^{*}+c^{\top}x^{*}+\frac{\kappa t}{\alpha(1-\alpha)}$ $\displaystyle\leq\frac{1}{2}x^{*\top}Qx^{*}+c^{\top}x^{*}+4\kappa t.$ $$

First step is by rearranging terms in the previous inequality.
Second step is by taking $\alpha=\frac{1}{2}$.
This finishes the proof.
∎

###### Lemma 9.16 ( [ 44 , Theorem 4.1.7 and Lemma 4.2.4] ) .

###### Lemma 9.17 (QP version of [ 40 , Lemma D.2] ) .

###### Proof.

###### Lemma 9.18 (QP version of [ 40 , Lemma D.3] ) .

###### Proof.

### 9.5 Proof of Theorem 9.1

In this section we combine everything and prove Theorem [9.1](#S9.Thmtheorem1).

###### Proof of Theorem 9.1 .

Lemma [9.17](#S9.Thmtheorem17) shows that the initial $x$ and $s$ satisfies

$$ $\displaystyle\|\mu\|_{w,x}^{*}\leq\epsilon.$ $$

This implies $w_{i}^{-1}\|\mu_{i}\|_{x_{i}}^{*}\leq\epsilon$ because $w_{i}\geq 1$.

Because $\epsilon\leq\frac{1}{\lambda}$, we have

$$ $\displaystyle\Phi(x,s,t)=\sum_{i\in[n]}\cosh(\lambda w_{i}^{-1}\|\mu_{i}\|_{x_{i}}^{*})\leq n\cosh(1)\leq\cosh(\lambda/64)$ $$

for the initial $x$ and $s$, by the choice of $\lambda$.

Using Theorem [9.15](#S9.Thmtheorem15), we see that

$$ $\displaystyle\Phi(x,s,t)\leq\cosh(\lambda/64)$ $$

during the entire algorithm.

So at the end of the algorithm, we have
$w_{i}^{-1}\|\mu_{i}\|_{x_{i}}^{*}\leq\frac{1}{64}$
for all $i\in[n]$.
In particular, $\|\mu_{i}\|_{x_{i}}^{*}\leq w_{i}$ for all $i\in[n]$.

Therefore, applying Lemma [9.18](#S9.Thmtheorem18) we get

$$ $\displaystyle\frac{1}{2}x^{\top}Qx+c^{\top}x$ $\displaystyle\leq\frac{1}{2}x^{*\top}Qx^{*}+c^{\top}x^{*}+4t\kappa$ $\displaystyle\leq\frac{1}{2}x^{*\top}Qx^{*}+c^{\top}x^{*}+\epsilon^{2}$ $$

where we used the stop condition for $t$ at the end.

So Lemma [9.17](#S9.Thmtheorem17) shows how to get an approximate solution for the original quadratic program with error $\epsilon LR(R+1)$.

The number of iterations is because we decrease $t$ by a factor of $1-h$ every iteration, and the choice $h=\frac{\alpha}{64\sqrt{\kappa}}$.
∎

###### Proof of Theorem 9.1 .

## 10 Gaussian Kernel SVM: Almost-Linear Time Algorithm and Hardness

In this section, we provide both algorithm and hardness for Gaussian kernel SVM problem. For the algorithm, we utilize a result due to [] in conjunction with our low-rank QP solver to obtain an $O(n^{1+o(1)}\log(1/\epsilon))$ time algorithm. For the hardness, we build upon the framework outlined in [] and improve their results in terms of dependence on dimension $d$.

We start by proving a simple lemma that shows that if $K=UV^{\top}$ for low-rank $U,V$, then the quadratic objective $K\circ(yy^{\top})$ also admits such a factorization via a simple scaling.

###### Lemma 10.1 .

Let $U,V\in\mathbb{R}^{n\times k}$ and $y\in\mathbb{R}^{n}$. Then, there exists a pair of matrices $\widetilde{U},\widetilde{V}\in\mathbb{R}^{n\times k}$ such that

$$ $\displaystyle\widetilde{U}\widetilde{V}^{\top}=$ $\displaystyle~{}(UV^{\top})\circ(yy^{\top})$ $$

moreover, $\widetilde{U},\widetilde{V}$ can be computed in time $O(nk)$.

###### Proof.

The proof relies on the following identity for Hadamard product: for any matrix $A$ and conforming vectors $x,y$ (all real), one has

$$ $\displaystyle A\circ(yx^{\top})=$ $\displaystyle~{}D_{y}AD_{x}$ $$

where $D_{y},D_{x}\in\mathbb{R}^{n\times n}$ are diagonal matrices that put $y,x$ on their diagonals. Thus, we can simply compute $\widetilde{U},\widetilde{V}$ as follows:

$$ $\displaystyle\widetilde{U}=$ $\displaystyle~{}D_{y}U,$ $\displaystyle\widetilde{V}=$ $\displaystyle~{}D_{y}V,$ $$

consequently,

$$ $\displaystyle\widetilde{U}\widetilde{V}^{\top}=$ $\displaystyle~{}D_{y}UV^{\top}D_{y}$ $\displaystyle=$ $\displaystyle~{}(yy^{\top})\circ(UV^{\top})$ $\displaystyle=$ $\displaystyle~{}(UV^{\top})\circ(yy^{\top}),$ $$

as desired. Moreover, the diagonal scaling of $U,V$ can be indeed performed in $O(nk)$ time, as advertised.
∎

Throughout this section, we will let $B$ denote the squared radius of the dataset.

###### Lemma 10.1 .

###### Proof.

### 10.1 Almost-Linear Time Algorithm for Gaussian Kernel SVM

We state a result due to [], in which they present an optimal-degree polynomial approximation to the function $e^{-x}$ and consequentially, this produces an efficient approximate scheme to the Batch Gaussian Kernel Density Estimation problem.

We start by introducing a notion that captures the minimum degree polynomial that well-approximates $e^{-x}$:

###### Definition 10.2 .

Let $f:[0,B]\rightarrow\mathbb{R}$, we let $q_{B;\epsilon}(f)\in\mathbb{N}$ denote the minimum degree of a non-constant polynomial $p(x)$ such that

$$ $\displaystyle\sup_{x\in[0,B]}|p(x)-f(x)|\leq\epsilon$ $$

Utilizing the Chebyshev polynomial machinery together with the orthgonal polynomial families, [] provides the following characterization on $q_{B;\epsilon}(f)$:

###### Theorem 10.3 (Theorem 1.2 of [ 1 ] ) .

Let $B\geq 1$ and $\epsilon\in(0,1)$. Then

$$ $\displaystyle q_{B;\epsilon}(e^{-x})=$ $\displaystyle~{}\Theta(\max\{\sqrt{B\log(1/\epsilon)},\frac{\log(1/\epsilon)}{\log(B^{-1}\log(1/\epsilon))}\})$ $$

###### Theorem 10.4 (Corollary 1.7 of [ 1 ] ) .

Let $x_{1},\ldots,x_{n}\in\mathbb{R}^{d}$ be a dataset with squared radius $B$ and $\epsilon\in(0,1)$. Let $q=q_{B;\epsilon}(e^{-x})$. Let $K\in\mathbb{R}^{n\times n}$ be the Gaussian kernel matrix formed by $x_{1},\ldots,x_{n}$. Finally, let $k=\binom{2d+2q}{2q}$. Then, there exists a deterministic algorithm that computes a pair of matrices $U,V\in\mathbb{R}^{n\times k}$ such that for any vector $v\in\mathbb{R}^{n}$,

$$ $\displaystyle\|Kv-UV^{\top}v\|_{\infty}\leq$ $\displaystyle~{}\epsilon\|v\|_{1}.$ $$

Moreover, matrices $U,V$ can be computed in time $O(nkd)$.

Even though $\ell_{\infty}$ error in terms of $\ell_{1}$ norm of vector $v$ seems quite weak, it can be conveniently translated into more standard guarantees, e.g., spectral norm error. The following lemma provides a conversion of errors that come in handy later when integrating the kernel approximation to our low-rank QP solver.

###### Lemma 10.5 .

Let $K\in\mathbb{R}^{n\times n}$ be a PSD kernel matrix and $\epsilon\in(0,1)$ be a parameter. Let $\widetilde{K}\in\mathbb{R}^{n\times n}$ be an approximation to $K$ with the guarantee that for any $v\in\mathbb{R}^{n}$,

$$ $\displaystyle\|Kv-\widetilde{K}v\|_{\infty}\leq$ $\displaystyle~{}\epsilon\|v\|_{1},$ $$

then

$$ $\displaystyle|v^{\top}Kv-v^{\top}\widetilde{K}v|\leq\epsilon\|v\|_{1}^{2}\leq\epsilon n\|v\|_{2}^{2}.$ $$

###### Proof.

The proof is a simple application of Hölder’s inequality:

$$ $\displaystyle|v^{\top}(Kv-\widetilde{K}v)|=$ $\displaystyle~{}|\langle v,Kv-\widetilde{K}v\rangle|$ $\displaystyle\leq$ $\displaystyle~{}\|v\|_{1}\|Kv-\widetilde{K}v\|_{\infty}$ $\displaystyle\leq$ $\displaystyle~{}\epsilon\|v\|_{1}^{2}$ $\displaystyle\leq$ $\displaystyle~{}\epsilon n\|v\|_{2}^{2},$ $$

where the second step is by Hölder’s inequality, and the last step is by Cauchy-Schwartz. This completes the proof.
∎

We can now combine the Gaussian kernel low-rank decomposition with our low-rank QP solver to provide an almost-linear time algorithm for Gaussian kernel SVM. We restate the kernel SVM formulation here.

###### Definition 10.6 (Restatement of Definition 1.3 ) .

Given a data matrix $X\in\mathbb{R}^{n\times d}$ and labels $y\in\mathbb{R}^{n}$.
Let $Q\in\mathbb{R}^{n\times n}$ denote a matrix where $Q_{i,j}=\mathsf{K}(x_{i},x_{j})\cdot y_{i}y_{j}$ for $i,j\in[n]$. The hard-mragin kernel SVM problem with bias asks to solve the following program.

$$ $\displaystyle\max_{\alpha\in\mathbb{R}^{n}}~{}$ $\displaystyle~{}{\bf 1}_{n}^{\top}\alpha-\frac{1}{2}\alpha^{\top}Q\alpha$ $\displaystyle{\rm s.t.~{}}$ $\displaystyle~{}\alpha^{\top}y=0$ $\displaystyle~{}\alpha\geq 0.$ $$

###### Theorem 10.7 .

Let Gaussian kernel SVM training problem be defined as above with kernel function ${\sf K}(x_{i},x_{j})=\exp(-\|x_{i}-x_{j}\|_{2}^{2})$. Suppose the dataset has squared radius $B\geq 1$, and let $\epsilon\in(0,1)$ be the precision parameter. Suppose the program satisfies the following:

- •
There exists a point $z\in\mathbb{R}^{n}$ such that there is an Euclidean ball with radius $r$ centered at $z$ that is contained in the constraint set.
- •
The constraint set is enclosed by an Euclidean ball of radius $R$, centered at the origin.

Then, there exists a randomized algorithm that outputs an approximate solution $\widehat{\alpha}\in\mathbb{R}^{n}$ such that $\widehat{\alpha}\geq 0$, moreover,

$$ $\displaystyle{\bf 1}_{n}^{\top}\widehat{\alpha}-\frac{1}{2}\widehat{\alpha}^{\top}Q\widehat{\alpha}\geq$ $\displaystyle~{}{\rm OPT}-\epsilon,$ $\displaystyle\|\widehat{\alpha}^{\top}y\|_{1}\leq$ $\displaystyle~{}3\epsilon,$ $$

where ${\rm OPT}$ denote the optimal cost of the objective function. Let $q=q_{B;\Theta(\epsilon/nR^{2})}(e^{-x})$ and $k=\binom{2d+2q}{2q}$. Then, the vector $\widehat{\alpha}$ can be computed in expected time

$$ $\displaystyle\widetilde{O}(nk^{(\omega+1)/2}\log(nR/(\epsilon r))).$ $$

###### Proof.

Throughout the proof, we set $\epsilon_{1}=O(\epsilon/(nR^{2}))$. We will craft an algorithm that first computes an approximate Gaussian kernel together with a proper low-rank factorization, then use this proxy kernel matrix to solve the quadratic program. We will use $K$ to denote the exact Gaussian kernel matrix, $Q$ to denote the exact quadratic matrix.

Approximate the Gaussian kernel matrix with finer granularity. We start by invoking Theorem [10.4](#S10.Thmtheorem4) using data matrix $X$ with accuracy parameter $\epsilon_{1}$. We let $\widetilde{K}=UV^{\top}$ to denote this approximate kernel matrix, and we let $\widetilde{Q}=D_{y}UV^{\top}D_{y}$ to denote the approximate quadratic matrix. Owing to Lemma [10.5](#S10.Thmtheorem5), we know that for any vector $x\in\mathbb{R}^{n}$,

$$ $\displaystyle|x^{\top}(Q-\widetilde{Q})x|=$ $\displaystyle~{}|(D_{y}x)^{\top}(K-\widetilde{K})(D_{y}x)|$ $\displaystyle\leq$ $\displaystyle~{}\epsilon_{1}n\|D_{y}x\|_{2}^{2}$ $\displaystyle=$ $\displaystyle~{}\epsilon_{1}n\|x\|_{2}^{2},$ $$

where we use the fact that $y\in\{\pm 1\}^{n}$. This also implies that

$$ $\displaystyle\|Q-\widetilde{Q}\|\leq$ $\displaystyle~{}\epsilon_{1}n$ (37) $$

this simple bound will come in handy later.

Solving the approximate program to high precision. Given $\widetilde{Q}$, we solve the following program:

$$ $\displaystyle\max_{\alpha\in\mathbb{R}^{n}}$ $\displaystyle~{}{\bf 1}_{n}^{\top}\alpha-\frac{1}{2}\alpha^{\top}\widetilde{Q}\alpha$ s.t. $\displaystyle~{}\alpha^{\top}y=0$ $\displaystyle~{}\alpha\geq 0$ $$

by invoking Theorem [8.1](#S8.Thmtheorem1). To do so, we need a bound on the Lipschitz constant of the program, i.e., the spectral norm of $\widetilde{Q}$ and $\ell_{2}$ norm of ${\bf 1}$. The latter is clearly $\sqrt{n}$, we shall show the first term is at most $(1+\epsilon_{1})\cdot n$.

Note that

$$ $\displaystyle\|Q\|=$ $\displaystyle~{}\|D_{y}KD_{y}\|$ $\displaystyle\leq$ $\displaystyle~{}\operatorname{tr}[D_{y}KD_{y}]$ $\displaystyle=$ $\displaystyle~{}\operatorname{tr}[K]$ $\displaystyle\leq$ $\displaystyle~{}n,$ $$

where we use $K$ is PSD. Combining with Eq. ([37](#S10.E37)) and triangle inequality, we have

$$ $\displaystyle\|\widetilde{Q}\|\leq$ $\displaystyle~{}\|Q\|+\|Q-\widetilde{Q}\|$ $\displaystyle\leq$ $\displaystyle~{}(1+\epsilon_{1})\cdot n.$ $$

With these Lipschitz constants, we examine the error guarantee provided by Theorem [8.1](#S8.Thmtheorem1): it produces a vector $\widehat{\alpha}\in\mathbb{R}^{n}$ such that

$$ $\displaystyle{\bf 1}_{n}^{\top}\widehat{\alpha}-\frac{1}{2}\widehat{\alpha}^{\top}\widetilde{Q}\widehat{\alpha}\geq$ $\displaystyle~{}\max_{\alpha^{\top}y=0,x\geq 0}({\bf 1}_{n}^{\top}\alpha-\frac{1}{2}\alpha^{\top}\widetilde{Q}\alpha)-O(\epsilon_{1}nR^{2}),$ $\displaystyle\|\widehat{\alpha}^{\top}y\|_{1}\leq$ $\displaystyle~{}O(\epsilon_{1}nR),$ $$

we mainly focus on the first error bound, as we need to understand the quality of $\widehat{x}$ when plug into the program with $Q$.

We will follow a chain of triangle inequalities, so we first bound

$$ $\displaystyle|\widehat{\alpha}^{\top}(\widetilde{Q}-Q)\widehat{\alpha}|\leq$ $\displaystyle~{}\epsilon n\|\widehat{\alpha}\|_{2}^{2}$ $\displaystyle\leq$ $\displaystyle~{}\epsilon nR^{2}.$ $$

Next, let

$$ $\displaystyle\alpha^{\prime}:=$ $\displaystyle~{}\arg\max_{\alpha^{\top}y=0,\alpha\geq 0}{\bf 1}_{n}^{\top}\alpha-\frac{1}{2}\alpha^{\top}\widetilde{Q}\alpha,$ $\displaystyle\alpha^{*}:=$ $\displaystyle~{}\arg\max_{\alpha^{\top}y=0,\alpha\geq 0}{\bf 1}_{n}^{\top}\alpha-\frac{1}{2}\alpha^{\top}Q\alpha,$ $$

then we have the following

$$ $\displaystyle{\bf 1}_{n}^{\top}\alpha^{\prime}-\frac{1}{2}\alpha^{\prime\top}\widetilde{Q}\alpha^{\prime}\geq$ $\displaystyle~{}{\bf 1}_{n}^{\top}\alpha^{*}-\frac{1}{2}(\alpha^{*})^{\top}\widetilde{Q}\alpha^{*}$ $\displaystyle\geq$ $\displaystyle~{}{\bf 1}_{n}^{\top}\alpha^{*}-\frac{1}{2}(\alpha^{*})^{\top}Q\alpha^{*}-O(\epsilon_{1}nR^{2})$ $\displaystyle=$ $\displaystyle~{}{\rm OPT}-O(\epsilon_{1}nR^{2}),$ $$

where the second step is by applying Lemma [10.5](#S10.Thmtheorem5) to $\alpha^{*}$. Now we are ready to bound the final error:

$$ $\displaystyle{\bf 1}_{n}^{\top}\widehat{\alpha}-\frac{1}{2}\widehat{\alpha}^{\top}Q\widehat{\alpha}\geq$ $\displaystyle~{}{\bf 1}_{n}^{\top}\widehat{\alpha}-\frac{1}{2}\widehat{\alpha}^{\top}\widetilde{Q}\widehat{\alpha}-O(\epsilon_{1}nR^{2})$ $\displaystyle\geq$ $\displaystyle~{}{\bf 1}_{n}^{\top}\alpha^{\prime}-\frac{1}{2}\alpha^{\prime\top}\widetilde{Q}\alpha^{\prime}-O(\epsilon_{1}nR^{2})$ $\displaystyle\geq$ $\displaystyle~{}{\rm OPT}-O(\epsilon_{1}nR^{2}).$ $$

The final error guarantee follows by the choice of $\epsilon_{1}$, and we indeed design an algorithm that outputs a vector $\widehat{\alpha}$ with

$$ $\displaystyle{\bf 1}^{\top}\widehat{\alpha}-\frac{1}{2}\widehat{\alpha}^{\top}Q\widehat{\alpha}\geq$ $\displaystyle~{}{\rm OPT}-\epsilon,$ $\displaystyle\|\widehat{\alpha}^{\top}y\|_{1}\leq$ $\displaystyle~{}\epsilon.$ $$

Runtime analysis. It remains to analyze the runtime of our proposed algorithm. We first compute an approximate kernel $\widetilde{K}$ with parameter $\epsilon_{1}$, owing to Theorem [10.4](#S10.Thmtheorem4), we have

$$ $\displaystyle q_{B;\epsilon_{1}}(e^{-x})=$ $\displaystyle~{}\Theta(\max\{\sqrt{B\log(nR/\epsilon)},\frac{\log(nR/\epsilon)}{\log(B^{-1}\log(nR/\epsilon))}\})$ $$

then by setting $k=\binom{2d+2q}{2q}$, the matrix $\widetilde{K}$ can be computed in time $O(nkd)$. Given this rank-$k$ factorization, the program can then be solved with precision $\epsilon_{1}$ in time

$$ $\displaystyle\widetilde{O}(nk^{(\omega+1)/2}\log(nR/(\epsilon r))),$ $$

as desired.
∎

###### Remark 10.8 .

To understand the value range of $k$, let us consider the set of parameters:

$$ $\displaystyle d=\Theta(\log n),\epsilon=1/\operatorname{poly}n,R=\operatorname{poly}n,B=\Theta(1),$ $$

under this setting, $O(\log(nR/\epsilon))=O(\log n)$ and the degree $q$ is

$$ $\displaystyle q=$ $\displaystyle~{}\Theta(\sqrt{\log n})$ $$

the rank $k$ is then

$$ $\displaystyle k=$ $\displaystyle~{}\binom{2d+2q}{2q}$ $\displaystyle\leq$ $\displaystyle~{}\Theta((\log n)^{\frac{1}{2}\sqrt{\log n}})$ $\displaystyle=$ $\displaystyle~{}\Theta(2^{\Theta(\log\log n\sqrt{\log n})})$ $\displaystyle=$ $\displaystyle~{}n^{o(1)},$ $$

consequentially, our algorithm runs in almost-linear time in $n$:

$$ $\displaystyle\widetilde{O}(n^{1+o(1)}\log n).$ $$

It is worth noting to achieve the almost-linear runtime, the data radius $B$ can be further relaxed. In fact, as long as

$$ $\displaystyle B=$ $\displaystyle~{}o\left(\frac{\log n}{\log\log n}\right),$ $$

we can ensure that $k=n^{o(1)}$ and subsequently the almost-linear runtime.

The runtime we obtain can be viewed as a consequence of the “phase transition” phenomenon observed in [], in which they prove that if $B=\omega(\log n)$, then quadratic time in $n$ is essentially needed to approximate the Gaussian kernel assuming SETH.

###### Definition 10.2 .

###### Theorem 10.3 (Theorem 1.2 of [ 1 ] ) .

###### Theorem 10.4 (Corollary 1.7 of [ 1 ] ) .

###### Lemma 10.5 .

###### Proof.

###### Definition 10.6 (Restatement of Definition 1.3 ) .

###### Theorem 10.7 .

###### Proof.

###### Remark 10.8 .

### 10.2 Hardness of Gaussian Kernel SVM with Large Radius

In this section, we show that for $d=\Theta(\log n)$, any algorithm that solves the program associated to hard-margin Gaussian kernel SVM would require $\Omega(n^{2-o(1)})$ time for $B=\omega(\log n)$. This justifies the choice of $B$ in Remark [10.8](#S10.Thmtheorem8). To prove the hardness result, we need to introduce the approximate Hamming nearest neighbor problem.

###### Definition 10.9 .

For $\delta>0$ and $n,d\in\mathbb{N}$, let $\{a_{1},\ldots,a_{n}\},\{b_{1},\ldots,b_{n}\}\subseteq\{0,1\}^{d}$ be two sets of vectors, and let $t\in\{0,1,\ldots,d\}$ be a threshold. The *$(1+\delta)$-Approximate Hamming Nearest Neighbor* problem asks to distinguish the following two cases:

- •
If there exists some $a_{i}$ and $b_{j}$ such that $\|a_{i}-b_{j}\|_{1}\leq t$, output “Yes”;
- •
If for any $i,j\in[n]$, we have $\|a_{i}-b_{j}\|_{1}>(1+\delta)\cdot t$, output “No”.

Note that the algorithm can output any value if the datasets fall in neither of these two cases. We will utilize the following hardness result due to Rubinstein.

###### Theorem 10.10 ( [ 50 ] ) .

Assuming SETH, for every $q>0$, there exists $\delta>0$ and $C>0$ such that $(1+\delta)$-Approximate Hamming Nearest Neighbor in dimension $d=C\log n$ requires time $\Omega(n^{2-q})$.

A final ingredient is a rewriting of the dual SVM into its primal form, without resorting to optimize over an infinite-dimensional hyperplane.

###### Lemma 10.11 .

Consider the dual hard-margin kernel SVM defined as

$$ $\displaystyle\max_{\alpha\in\mathbb{R}^{n}}$ $\displaystyle~{}{\bf 1}^{\top}\alpha-\frac{1}{2}\sum_{i,j\in[n]\times[n]}\alpha_{i}\alpha_{j}y_{i}y_{j}{\sf K}(w_{i},w_{j}),$ s.t. $\displaystyle~{}\alpha^{\top}y=0,$ $\displaystyle~{}\alpha\geq 0.$ $$

The primal program can be written as

$$ $\displaystyle\min_{\alpha\in\mathbb{R}^{n}}$ $\displaystyle~{}\frac{1}{2}\sum_{i,j\in[n]\times[n]}\alpha_{i}\alpha_{j}y_{i}y_{j}{\sf K}(w_{i},w_{j}),$ s.t. $\displaystyle~{}y_{i}f(w_{i})\geq 1,$ $\displaystyle~{}\alpha\geq 0,$ $$

where $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is defined as

$$ $\displaystyle f(w)=$ $\displaystyle~{}\sum_{j=1}^{n}\alpha_{j}y_{j}{\sf K}(w_{j},w)-b.$ $$

Moreover, the primal and dual program has no duality gap and the optimal solution $\alpha$ to both programs are the same.

###### Proof.

Recall that the primal hard-margin SVM is the following program:

$$ $\displaystyle\min_{v}$ $\displaystyle~{}\frac{1}{2}\|v\|_{2}^{2},$ s.t. $\displaystyle~{}y_{i}(v^{\top}\phi(w_{i})-b)\geq 1,$ $$

where $b\in\mathbb{R}$ is the bias term and $\phi:\mathbb{R}^{d}\rightarrow\mathbb{R}^{K}$ is the feature mapping corresponding to the kernel in the sense that ${\sf K}(w_{i},w_{j})=\phi(w_{i})^{\top}\phi(w_{j})$. The optimal weight $v=\sum_{i=1}^{n}\alpha_{i}y_{i}\phi(w_{i})$ where $\alpha\in\mathbb{R}^{n}$ is the optimal solution to the dual program. Consequently,

$$ $\displaystyle\|v\|_{2}^{2}=$ $\displaystyle~{}(\sum_{i=1}^{n}\alpha_{i}y_{i}\phi(w_{i}))^{2}$ $\displaystyle=$ $\displaystyle~{}\sum_{i,j\in[n]\times[n]}\alpha_{i}\alpha_{j}y_{i}y_{j}\phi(w_{i})^{\top}\phi(w_{j})$ $\displaystyle=$ $\displaystyle~{}\sum_{i,j\in[n]\times[n]}\alpha_{i}\alpha_{j}y_{i}y_{j}{\sf K}(w_{i},w_{j})$ $\displaystyle=$ $\displaystyle~{}\alpha^{\top}Q\alpha,$ $$

where the matrix $Q$ is the usual

$$ $\displaystyle Q=$ $\displaystyle~{}(yy^{\top})\circ K,$ $$

the constraint can be rewritten as

$$ $\displaystyle y_{i}(v^{\top}\phi(w_{i})-b)=$ $\displaystyle~{}y_{i}((\sum_{i=1}^{n}\alpha_{i}y_{i}\phi(w_{i}))^{\top}\phi(w_{i})-b)$ $\displaystyle=$ $\displaystyle~{}y_{i}(\sum_{j=1}^{n}\alpha_{j}y_{j}\phi(w_{j})^{\top}\phi(w_{i}))-y_{i}b$ $\displaystyle=$ $\displaystyle~{}y_{i}(\sum_{j=1}^{n}\alpha_{j}y_{j}{\sf K}(w_{i},w_{j}))-y_{i}b$ $\displaystyle=$ $\displaystyle~{}y_{i}f(w_{i}),$ $$

where $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is defined as

$$ $\displaystyle f(w)=$ $\displaystyle~{}\sum_{j=1}^{n}\alpha_{j}y_{j}{\sf K}(w_{j},w)-b.$ $$

Thus, we can alternatively write the primal as

$$ $\displaystyle\min_{\alpha\in\mathbb{R}^{n}}$ $\displaystyle~{}\frac{1}{2}\alpha^{\top}Q\alpha,$ s.t. $\displaystyle~{}y_{i}f(w_{i})\geq 1.$ $$

For the strong duality and optimal solution, see, e.g., [].
∎

We will now prove the almost-quadratic lower bound for Guassian kernel SVM. Our proof strategy is similar to that of [] with different set of parameters. It is also worth noting that the [] construction

- •
Requires the dimension $d=\Theta(\log^{3}n)$;
- •
Requires the squared dataset radius $B=\Theta(\log^{4}n)$.

We will improve both of these results.

###### Theorem 10.12 .

Assuming SETH, for every $q>0$, there exists a hard-margin Gaussian kernel SVM without the bias term as defined in Definition [1.3](#S1.Thmtheorem3) with $d=\Theta(\log n)$ and error $\epsilon=\exp(-\Theta(\log^{2}n))$ for inputs whose squared radius is at most $B=\Theta(\log^{2}n)$ requiring time $\Omega(n^{2-q})$ to solve.

###### Proof.

Let $l=\sqrt{2(c^{\prime}\delta)^{-1}\log n}$. We will provide a reduction from $(1+\delta)$-Approximate Hamming Nearest Neighbor to Gaussian kernel SVM. Let $A:=\{a_{1},\ldots,a_{n}\},B:=\{b_{1},\ldots,b_{n}\}\subseteq\{0,1\}^{d}$ be the datasets, we assign label 1 to all vectors $a_{i}$ and label $-1$ to all vectors $b_{j}$, moreover, we scale both $A$ and $B$ by $l$, this results in two datasets with points in $\{0,l\}^{d}$. The squared radius of this dataset is then

$$ $\displaystyle B=$ $\displaystyle~{}\max\{\max_{i,j}\|la_{i}-la_{j}\|_{2}^{2},\max_{i,j}\|lb_{i}-lb_{j}\|_{2}^{2},\max_{i,j}\|la_{i}-lb_{j}\|_{2}^{2}\}$ $\displaystyle\leq$ $\displaystyle~{}l^{2}d$ $\displaystyle=$ $\displaystyle~{}\Theta(\delta^{-1}\log^{2}n).$ $$

To simplify the notation, we will implictly assume $A$ and $B$ are scaled by $l$ without explicitly writing out $la_{i}$, $lb_{j}$. Now consider the following three programs:

- •
Classifying $A$:
$\displaystyle\min_{\alpha\in\mathbb{R}^{n}_{\geq 0}}$
$\displaystyle~{}\frac{1}{2}\sum_{i,j\in[n]\times[n]}\alpha_{i}\alpha_{j}{\sf K}(a_{i},a_{j}),$
s.t.
$\displaystyle~{}\sum_{j=1}^{n}\alpha_{j}{\sf K}(a_{i},a_{j})\geq 1,\qquad\forall i\in[n]$
(38)
- •
Classifying $B$:
$\displaystyle\min_{\beta\in\mathbb{R}^{n}_{\geq 0}}$
$\displaystyle~{}\frac{1}{2}\sum_{i,j\in[n]\times[n]}\beta_{i}\beta_{j}{\sf K}(b_{i},b_{j}),$
s.t.
$\displaystyle~{}-\sum_{j=1}^{n}\beta_{j}{\sf K}(b_{i},b_{j})\leq-1,\qquad\forall i\in[n]$
(39)
- •
Classifying both $A$ and $B$:
$\displaystyle\min_{\alpha,\beta\in\mathbb{R}_{\geq 0}^{n}}$
$\displaystyle~{}\frac{1}{2}\sum_{i,j\in[n]\times[n]}\alpha_{i}\alpha_{j}{\sf K}(a_{i},a_{j})+\frac{1}{2}\sum_{i,j\in[n]\times[n]}\beta_{i}\beta_{j}{\sf K}(b_{i},b_{j})-\sum_{i,j\in[n]\times[n]}\alpha_{i}\beta_{j}{\sf K}(a_{i},b_{j}),$
s.t.
$\displaystyle~{}\sum_{j=1}^{n}\alpha_{j}{\sf K}(a_{i},a_{j})-\sum_{j=1}^{n}\beta_{j}{\sf K}(a_{i},b_{j})\geq 1,\qquad\forall i\in[n],$
$\displaystyle~{}\sum_{j=1}^{n}\alpha_{j}{\sf K}(b_{i},a_{j})-\sum_{j=1}^{n}\beta_{j}{\sf K}(b_{i},b_{j})\leq-1,\qquad\forall i\in[n]$
(40)

We will prove that the optimal solution $\alpha^{*}_{i}$’s and $\beta^{*}_{i}$’s are both lower and upper bounded. Use ${\rm Val}(A),{\rm Val}(B)$ and ${\rm Val}(A,B)$ to denote the value of program ([• ‣ 10.2](#S10.Ex423)), ([• ‣ 10.2](#S10.Ex424)) and ([• ‣ 10.2](#S10.Ex425)) respectively, then note that

$$ $\displaystyle{\rm Val}(A)\leq$ $\displaystyle~{}\frac{n^{2}}{2}$ $$

by plugging in $\alpha={\bf 1}$ and setting all vectors to be the same. On the other hand,

$$ $\displaystyle{\rm Val}(A)\geq$ $\displaystyle~{}\frac{1}{2}\sum_{i=1}^{n}(\alpha_{i}^{*})^{2}{\sf K}(a_{i},a_{i})$ $\displaystyle=$ $\displaystyle~{}\frac{1}{2}\sum_{i=1}^{n}(\alpha_{i}^{*})^{2}.$ $$

Combining these two, we can conclude that for any $\alpha_{i}^{*}$, it must be $\alpha_{i}^{*}\leq n$. For the lower bound, consider the inequality constraint for the $i$-th point:

$$ $\displaystyle\alpha_{i}^{*}+\sum_{j\neq i}\alpha_{j}^{*}{\sf K}(a_{i},a_{j})\geq$ $\displaystyle~{}1,$ $$

to estimate ${\sf K}(a_{i},a_{j})$, note that $\|a_{i}-a_{j}\|_{2}^{2}=\|a_{i}-a_{j}\|_{1}\geq 1$ for $j\neq i$,(^8^88We without loss of generality that during preprocess, we have remove duplicates in $A$ and $B$.) and

$$ $\displaystyle{\sf K}(a_{i},a_{j})=$ $\displaystyle~{}\exp(-l^{2}\|a_{i}-a_{j}\|_{2}^{2})$ $\displaystyle=$ $\displaystyle~{}\exp(-2(c^{\prime}\delta)^{-1}\log n\|a_{i}-a_{j}\|_{1})$ $\displaystyle\leq$ $\displaystyle~{}\exp(-2(c^{\prime}\delta)^{-1}\log n)$ $\displaystyle\leq$ $\displaystyle~{}n^{-10}/100,$ $$

combining with $\alpha_{j}^{*}\leq n$, we have

$$ $\displaystyle\alpha_{i}^{*}\geq$ $\displaystyle~{}1-\sum_{j\neq i}\alpha_{j}^{*}{\sf K}(a_{i},a_{j})$ $\displaystyle\geq$ $\displaystyle~{}1-n\cdot n\cdot n^{-10}/100$ $\displaystyle\geq$ $\displaystyle~{}1/2.$ $$

This lower bound on $\alpha_{i}^{*}$ is helpful when we attempt to lower bound ${\rm Val}(A,B)$ with ${\rm Val}(A)+{\rm Val}(B)$. Following the outline of [], we consider the three dual programs:

- •
Dual of classifying $A$:
$\displaystyle\max_{\alpha\in\mathbb{R}^{n}_{\geq 0}}$
$\displaystyle~{}\sum_{i=1}^{n}\alpha_{i}-\frac{1}{2}\sum_{i,j}\alpha_{i}\alpha_{j}{\sf K}(a_{i},a_{j})$
(41)
- •
Dual of classifying $B$:
$\displaystyle\max_{\beta\in\mathbb{R}^{n}_{\geq 0}}$
$\displaystyle~{}\sum_{i=1}^{n}\beta_{i}-\frac{1}{2}\sum_{i,j}\beta_{i}\beta_{j}{\sf K}(b_{i},b_{j})$
(42)
- •
Dual of classifying $A$ and $B$:
$\displaystyle\max_{\alpha,\beta\in\mathbb{R}^{n}_{\geq 0}}$
$\displaystyle~{}\sum_{i=1}^{n}\alpha_{i}+\sum_{i=1}^{n}\beta_{i}-\frac{1}{2}\sum_{i,j}\alpha_{i}\alpha_{j}{\sf K}(a_{i},a_{j})-\frac{1}{2}\sum_{i,j}\beta_{i}\beta_{j}{\sf K}(b_{i},b_{j})+\sum_{i,j}\alpha_{i}\beta_{j}{\sf K}(a_{i},b_{j})$
(43)

as the SVM program exhibits strong duality, we know that the optimal value of the primal equals to the dual, so we can alternatively bound ${\rm Val}(A,B)$ using the dual program. Plug in $\alpha^{*},\beta^{*}$ to program ([43](#S10.E43)), we have

$$ $\displaystyle{\rm Val}(A,B)\geq$ $\displaystyle~{}\sum_{i=1}^{n}\alpha_{i}^{*}+\sum_{i=1}^{n}\beta_{i}^{*}-\frac{1}{2}\sum_{i,j}\alpha_{i}^{*}\alpha_{j}^{*}{\sf K}(a_{i},a_{j})-\frac{1}{2}\sum_{i,j}\beta_{i}^{*}\beta_{j}^{*}{\sf K}(b_{i},b_{j})+\sum_{i,j}\alpha_{i}^{*}\beta_{j}^{*}{\sf K}(a_{i},b_{j})$ $\displaystyle=$ $\displaystyle~{}{\rm Val}(A)+{\rm Val}(B)+\sum_{i,j}\alpha_{i}^{*}\beta_{j}^{*}{\sf K}(a_{i},b_{j}),$ $$

to bound the third term, we consider the pair $(a_{i_{0}},b_{j_{0}})$ such that $\|a_{i_{0}}-b_{j_{0}}\|_{1}\leq t-1$, and note that

$$ $\displaystyle\sum_{i,j}\alpha_{i}^{*}\beta_{j}^{*}{\sf K}(a_{i},b_{j})\geq$ $\displaystyle~{}\alpha_{i_{0}}^{*}\beta_{j_{0}}^{*}{\sf K}(a_{i_{0}},b_{j_{0}})$ $\displaystyle\geq$ $\displaystyle~{}\frac{1}{4}\exp(-2(c^{\prime}\delta)^{-1}\log n\cdot(t-1)).$ $$

To wrap up, we have

$$ $\displaystyle{\rm Val}(A,B)\geq$ $\displaystyle~{}{\rm Val}(A)+{\rm Val}(B)+\frac{1}{4}\exp(-2(c^{\prime}\delta)^{-1}\log n\cdot(t-1))$ $$

We now prove the “No” case, where for any $a_{i},b_{j}$, $\|a_{i}-b_{j}\|_{1}\geq t$. We have

$$ $\displaystyle{\sf K}(a_{i},b_{j})=$ $\displaystyle~{}\exp(-l^{2}\|a_{i}-b_{j}\|_{2}^{2})$ $\displaystyle\leq$ $\displaystyle~{}\exp(-2(c^{\prime}\delta)^{-1}\log n\cdot t),$ $$

we let $m:=\exp(-2(c^{\prime}\delta)^{-1}\log n\cdot t)$, set $\alpha^{\prime}=\alpha^{*}+10n^{2}m$ and $\beta^{\prime}=\beta^{*}+10n^{2}m$, we let $V$ to denote the value when evaluating program ([• ‣ 10.2](#S10.Ex425)) with $\alpha^{\prime},\beta^{\prime}$. We will essentially show that

$$ $\displaystyle{\rm Val}(A,B)\leq$ $\displaystyle~{}V$ $$

and

$$ $\displaystyle V\leq$ $\displaystyle~{}{\rm Val}(A)+{\rm Val}(B)+400n^{6}m$ $$

chaining these two gives us a certificate for the “No” case. To prove the first assertion, we show that $\alpha^{\prime},\beta^{\prime}$ are feasible solution to program ([• ‣ 10.2](#S10.Ex425)) since

$$ $\displaystyle\sum_{j=1}^{n}\alpha^{\prime}_{j}{\sf K}(a_{i},a_{j})=$ $\displaystyle~{}\sum_{j=1}^{n}(\alpha^{*}_{j}{\sf K}(a_{i},a_{j})+10n^{2}m{\sf K}(a_{i},a_{j}))$ $\displaystyle=$ $\displaystyle~{}\alpha^{*}_{i}+10n^{2}m+\sum_{j\neq i}(\alpha^{*}_{j}+10n^{2}m){\sf K}(a_{i},a_{j})$ $\displaystyle\geq$ $\displaystyle~{}\alpha_{i}^{*}+\sum_{j\neq i}\alpha_{j}^{*}{\sf K}(a_{i},a_{j})+10n^{2}m$ $\displaystyle=$ $\displaystyle~{}10n^{2}m+\sum_{j=1}^{n}\alpha^{*}_{j}{\sf K}(a_{i},a_{j})$ $\displaystyle\geq$ $\displaystyle~{}10n^{2}m+1$ $$

where we use $\alpha_{i}^{*}$ satisfy the inequality constraint of program ([• ‣ 10.2](#S10.Ex423)). We compute an upper bound on $\sum_{j=1}^{n}\beta_{j}^{\prime}{\sf K}(a_{i},b_{j})$:

$$ $\displaystyle\sum_{j=1}^{n}\beta_{j}^{\prime}{\sf K}(a_{i},b_{j})\leq$ $\displaystyle~{}\sum_{j=1}^{n}2nm$ $\displaystyle=$ $\displaystyle~{}2n^{2}m,$ $$

where we use the fact that $m=\exp(-2(c^{\prime}\delta)^{-1}\log n\cdot t)\leq n^{-10}/10$ therefore $\beta^{*}+10n^{2}m\leq 2n$. Thus, it must be the case that

$$ $\displaystyle\sum_{j=1}\alpha_{j}^{\prime}{\sf K}(a_{i},a_{j})-\sum_{j=1}^{n}\beta^{\prime}_{j}{\sf K}(a_{i},b_{j})\geq$ $\displaystyle~{}8n^{2}m+1$ $\displaystyle\geq$ $\displaystyle~{}1,$ $$

as desired. The other linear constraint follows by a symmetric argument. This indeed shows that $\alpha^{\prime},\beta^{\prime}$ are feasible solutions to program ([• ‣ 10.2](#S10.Ex425)) and ${\rm Val}(A,B)\leq V$.

To prove an upper bound on $V$, we note that

$$ $\displaystyle V=$ $\displaystyle~{}\frac{1}{2}\sum_{i,j}\alpha_{i}^{\prime}\alpha_{j}^{\prime}{\sf K}(a_{i},a_{j})+\frac{1}{2}\sum_{i,j}\beta_{i}^{\prime}\beta_{j}^{\prime}{\sf K}(b_{i},b_{j})-\sum_{i,j}\alpha^{\prime}\beta_{j}^{\prime}{\sf K}(a_{i},b_{j})$ $\displaystyle\leq$ $\displaystyle~{}\frac{1}{2}\sum_{i,j}\alpha_{i}^{\prime}\alpha_{j}^{\prime}{\sf K}(a_{i},a_{j})+\frac{1}{2}\sum_{i,j}\beta_{i}^{\prime}\beta_{j}^{\prime}{\sf K}(b_{i},b_{j}),$ $$

we bound the first quantity, as the second follows similarly:

$$ $\displaystyle\frac{1}{2}\sum_{i,j}\alpha_{i}^{\prime}\alpha_{j}^{\prime}{\sf K}(a_{i},a_{j})=$ $\displaystyle~{}\frac{1}{2}\sum_{i,j}(\alpha_{i}^{*}\alpha_{j}^{*}+10n^{2}m(\alpha_{i}^{*}+\alpha_{j}^{*})+100n^{4}m^{2}){\sf K}(a_{i},a_{j})$ $\displaystyle\leq$ $\displaystyle~{}{\rm Val}(A)+\sum_{i,j}10n^{3}m{\sf K}(a_{i},a_{j})+\sum_{i,j}100n^{4}m^{2}{\sf K}(a_{i},a_{j})$ $\displaystyle\leq$ $\displaystyle~{}{\rm Val}(A)+10n^{5}m+100n^{6}m^{2}$ $\displaystyle\leq$ $\displaystyle~{}{\rm Val}(A)+200n^{6}m,$ $$

we can thus conclude

$$ $\displaystyle V\leq$ $\displaystyle~{}{\rm Val}(A)+{\rm Val}(B)+400n^{6}m.$ $$

Chaining these two, we obtain the following threshold for the “No” case:

$$ $\displaystyle{\rm Val}(A,B)\leq$ $\displaystyle~{}{\rm Val}(A)+{\rm Val}(B)+400n^{6}m.$ $$

Finally, we observe that

$$ $\displaystyle 400n^{6}\exp(-2(c^{\prime}\delta)^{-1}\log n\cdot t)\ll$ $\displaystyle~{}\frac{1}{4}\exp(-2(c^{\prime}\delta)^{-1}\log n\cdot(t-1)),$ $$

we can therefore distinguish these two cases.

Note that when one considers solving the program with additive error, we need to make sure that the error is smaller than the distinguishing threshold, i.e.,

$$ $\displaystyle\epsilon\leq$ $\displaystyle~{}\frac{1}{4}\exp(-2(c^{\prime}\delta)^{-1}\log n\cdot(t-1))$ $\displaystyle\leq$ $\displaystyle~{}\frac{1}{4}\exp(-2(c^{\prime}\delta)^{-1}d\log n)$ $\displaystyle=$ $\displaystyle~{}\exp(-\Theta(\log^{2}n)),$ $$

where we use $t\leq d$ and $d=\Theta(\log n)$. This concludes the proof.
∎

###### Remark 10.13 .

Our proof can be interpreted as using a stronger complexity theoretical tool in place of the one used by [], to obtain a better dependence on dimension $d$ and $B$. We also note that the construction due to [] has the relation that $B=\Theta(d\log n)$, this is because in order to lower bound ${\rm Val}(A,B)$, one has to lower bound the optimal values of $\alpha_{i}^{*}$’s and $\beta_{j}^{*}$’s. To do so, one needs to further scale up $a_{i}$’s and $b_{j}$’s so that within datasets $A$ and $B$, the radius is at least $\Theta(\log n)$. This is in contrast to the Batch Gaussian KDE studied in [], where they show the almost-quadratic lower bound can be achieved for both $d,B=\Theta(\log n)$.

Similar to [], we obtain hardness results for hard-margin kernel SVM with bias, and soft-margin kernel SVM with bias.

###### Corollary 10.14 .

Assuming SETH, for every $q>0$, there exists a hard-margin Gaussian kernel SVM with the bias term with $d=\Theta(\log n)$ and error $\epsilon=\exp(-\Theta(\log^{2}n))$ for inputs whose squared radius is at most $B=\Theta(\log^{6}n)$ requiring time $\Omega(n^{2-q})$ to solve.

###### Proof.

The proof is similar to []. Given a hard instance of Theorem [10.12](#S10.Thmtheorem12), except we append $\Theta(\log n)$ entries with magnitude $\Theta(\log^{2}n)$ instead of distributing the values across $\Theta(\log^{3}n)$ entries. Rest of the proof follows exactly the same as [].
∎

###### Corollary 10.15 .

Assuming SETH, for every $q>0$, there exists a soft-margin Gaussian kernel SVM with the bias term with $d=\Theta(\log n)$ and error $\epsilon=\exp(-\Theta(\log^{2}n))$ for inputs whose squared radius is at most $B=\Theta(\log^{6}n)$ requiring time $\Omega(n^{2-q})$ to solve.

###### Remark 10.16 .

Compared to the construction of [] in which they distribute a total mass of $\Theta(\log^{3}n)$ across $\Theta(\log^{3}n)$ entries so that they ensure after the reduction, the vectors take values in $\{-1,0,1\}$, we instead distribute the mass across $\Theta(\log n)$ entries so that each entry has magnitude $\Theta(\log^{2}n)$. To make the reduction work, the total mass of $\Theta(\log^{3}n)$ is needed, and for [], it is fine to append an extra $\Theta(\log^{3}n)$ entries as their hardness result for hard-margin SVM without bias does require $d=\Theta(\log^{3}n)$. For us, we need to restrict $d=\Theta(\log n)$ at the price of each entry has a larger magnitude of $\Theta(\log^{2}n)$. This blows up the squared radius from $\log^{2}n$ to $\log^{6}n$. We note that the [] construction has squared radius $\log^{4}n$.

###### Definition 10.9 .

###### Theorem 10.10 ( [ 50 ] ) .

###### Lemma 10.11 .

###### Proof.

###### Theorem 10.12 .

###### Proof.

###### Remark 10.13 .

###### Corollary 10.14 .

###### Proof.

###### Corollary 10.15 .

###### Remark 10.16 .

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. I know Eq. (1) is quadratic programming, but I believe the inequality constraint should be generalized into $Ax \le 0$, instead of simple $x\le 0$, though it is a special case.
2. For example, when ε is set to be 10^{−3} to account for the usual machine precision errors, these algorithms would require at least 10^6 iterations to converge. I don't understand this, as the complexity is O(d/ε), I thought it should be 1000*d. 
3. I know this is a theoretic paper, but I still hope the authors can compare the convergence with others which will be more convincing to reviewers and audiences. My experience told me that though libsvm is old, the speed is amazingly fast. I definitely hope the proposed algorithm can outperform libsvm.
soundness: 3
strengths: The paper clearly states the contribution, the story is easy to follow. This is a solid work on solving SVM dual problem based on quadratic programming.
summary: The authors design the first nearly-linear time algorithm for solving SVM problems via quadratic programs whenever the quadratic objective admits a low-rank factorization, and the number of linear constraints is small. They have derived the algorithm complexity analysis with different cases, which improves the existing result.
weaknesses: 1. No experiment to demonstrate the superiority to validate the new claimed contributions.
2. The result works for vanilla SVM and Gaussian Kernel, however, the analysis on other kernel options such as polynomial is blank, thus, whether it can be generalized into broader applications remains unknown. 
3. What if the kernel is a combination of kernels, say K_new=K1+K2 where K1 and K2 are Guassian kernel?

---

confidence: 4
contribution: 4
initial_rating: 8
presentation: 4
questions: The $O(\epsilon^{-2})$ time copmlexity mentioned on Line 121 was sharpened to $O(\epsilon^{-1})$ in [1]. Of course the point that the author makes about the logarithmic dependency still stands, but I think the improvement in [1] should be mentioned since the issue raised on line 124 is alleviated to some degree.

The complexity lower bound on the SVMs considered also implies complexity lower bound on QPs, correct? How does this compare to existing complexity lower bounds for QPs, if any?

Is there a reference for the result on $T$ on line 291?

What is the graph structure on the bags on line 310?

Can the authors give references for the sparse Cholesky factory result on Line 318?

What happens in the kernel SVM case when $d$ is constant (both for the upper bound and lower bound)? Why must $d$ grow with $n$ in order for the result to hold? For the lower bound, I presume it has to do with leveraging hardness results in [2]. Can the authors explain more about the $d = \Theta (\log n)$ assumption?

---

[1] Shalev-Shwartz, Shai, Yoram Singer, and Nathan Srebro. "Pegasos: Primal estimated sub-gradient solver for svm." Proceedings of the 24th international conference on Machine learning. 2007.

[2] Alman, Josh, and Ryan Williams. "Probabilistic polynomials and hamming nearest neighbors." 2015 IEEE 56th Annual Symposium on Foundations of Computer Science. IEEE, 2015.
soundness: 4
strengths: SVMs are classical subjects in machine learning, and thoroughly understanding the time complexity of optimizing SVMs is important. This paper takes an important step toward this. The technique for solving QP faster using IPM with the efficient Hessian updates could be beneficial for other problems in machine learning.
summary: The authors develop faster quadratic program solvers using interior point methods for problems where the quadratic form matrix has low rank. The technique is then applied to yield SVM solvers with complexity linearly/nearly-linearly in the number of samples for linear/kernel SVMs.
weaknesses: While this paper addresses a theoretical problem in complexity of SVMs, it'd still be interesting to see some actual experimental results.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 3
questions: 1.	In section 1, the authors denote $B$ as the squared radius, then in section 2, the authors redefine $B = Q + t H$. Please make sure the notation is constant.

2.	Though the theoretical results of this paper are very strong, it is highly recommended that the authors could do some experiments to verify the theoretical results.

3.	It seems the assumption that square radius $B = O(log(n)/loglog(n))$ is very strong. It would be great if the authors could share some kernel SVM problems which satisfy this assumption. 

4.	Overall, in section 2, The technical details are very heavy but not clear. It would great if the authors could give more intuitions about these approaches. For example, consider a more detailed discussion on how the low-rank factorization impacts specific types of SVM problems or datasets.

5.	Small typos: for example, in line 1540, is it “memebers” or “members”? Please double check.
soundness: 4
strengths: Novel Approach: The paper proposes an innovative nearly-linear time algorithm for solving quadratic programming problems with applications to kernel SVMs. The use of robust interior point method, low-rank factorization and low-rank matrix approximation make the approach computationally appealing.

Relevance: Given the broad usage of SVMs in machine learning, the proposed methodology could provide substantial computational improvements if the assumptions of the method are satisfied.

Complexity Analysis: The paper includes a rigorous complexity analysis and explores different scenarios for kernel SVMs under various dataset radius conditions. These insights provide a clear understanding of the expected performance under different cases. It also adheres to previous literature.
summary: The paper proposed a nearly-linear time algorithm for solving quadratic programming problems, especially designed for kernel SVMs. It combines the robust interior point method, low-treewidth setting, low-rank factorization and low-rank kernel matrix approximation techniques to enable the almost-linear time convergence. Though the assumptions to achieve this convergence rate is relatively strong, it is the first almost-linear convergence algorithm for quadratic programming problems or kernel SVMs.
weaknesses: Empirical Validation: While the theoretical contribution is significant, empirical results demonstrating the practical runtime improvements and accuracy on real datasets would be beneficial. This would help validate the theoretical findings and highlight the algorithm’s scalability and precision.

Assumptions and Constraints: The assumptions, especially regarding the squared radius, may limit applicability in some scenarios. The paper could benefit from discussing the practical implications of these assumptions on dataset characteristics.

Comparison with Existing Methods: Though the theoretical basis is sound, a more detailed comparison with existing SVM solvers (e.g., LIBSVM, SVM-Light) in terms of computational complexity and memory usage could strengthen the contribution by illustrating the specific benefits and trade-offs of the proposed algorithm.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: see strengths and weaknesses.

**Comment after the rebuttal period**: the authors' sincere responses resolved some of my concerns. I would like to express my gratitude for their efforts. Therefore, I decided to raise my score. However, I still believe that adding numerical experiments is necessary to ensure acceptance.
soundness: 3
strengths: - The authors provide intensive mathematical analysis of robust interior point methods for quadratic optimization problems.
summary: This paper studies fast algorithms for structured linear and kernel support vector machines (SVMs). SVMs are a type of machine learning model used for tasks such as classification and regression. The main contribution of the paper is to provide an algorithm with almost linear time for solving SVM problems with low tree width or low-rank structure. This is related to sparse patterns in the Cholesky decomposition of a matrix, etc. The paper also considers how to apply the proposed method to SVMs using Gaussian kernels, and provides an estimate of the computational complexity. No numerical investigation is given into the effectiveness of the method for real-world problems, etc.
weaknesses: - The paper states that the dimensionality constraint is d=O(log(n)), but this seems like a very strict assumption for real-world data analysis. It would be good to consider how the computational complexity changes when this constraint is relaxed.
- The effectiveness of the proposed method is unclear because there are no numerical experiments. It would be useful for practitioners to see things like implementation ideas and estimates of the data size that can actually be handled.
- As only the Gaussian kernel is used as the kernel, the range of applicability is narrow. Is it possible to apply the algorithm to other kernels, while maintaining the computational efficiency? How does the computation time of the algorithm relate to the kernel width of the Gaussian kernel? 
- For kernel methods, low-rank approximation methods such as the Nystrom approximation have been developed. In addition, it is known that when there are a large number of data points, it is more computationally efficient to use deep neural networks, which have developed significantly in recent years. Given this situation, it is necessary to explain the situations in which the proposed method in this paper is effective based on real-world problems.

Meta-review: This paper investigates solving SVMs via quadratic programming and proposes a novel, efficient, nearly-linear time algorithm for cases where the quadratic objective admits a low-rank factorization and the number of linear constraints is small—settings that are typical for SVMs. The results presented are of broad interest to the ML and optimization communities.

The reviewers reached a consensus that the contribution of this work is significant, and I therefore recommend accepting the paper. However, as suggested by **all** reviewers, it would be beneficial if the authors could provide numerical results, in some form, to demonstrate the practical numerical advantages of the proposed algorithm. 
Including such results would likely enhance the paper’s impact within the community.

Decision: Accept (Poster)