Faster Algorithms for Structured Linear and Kernel Support Vector Machines
Authors: Yuzhou Gu
Zhao Song
Lichen Zhang
Venue: ICLR 2025 Poster (2025)
Keywords: interior point method
support vector machine
data structure
Primary area: optimization
TL;DR: We design the first nearly-linear time, interior point algorithm for approximately solving support vector machine programs.

Abstract: Quadratic programming is a ubiquitous prototype in convex programming. Many machine learning problems can be formulated as quadratic programming, including the famous 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 $\Theta(n^2)$, where $n$ is the number of variables. Assuming the Strong Exponential Time Hypothesis ($\textsf{SETH}$), it is known that no $O(n^{2-o(1)})$ time algorithm exists when the quadratic objective matrix is positive semidefinite (Backurs, Indyk, and Schmidt, NeurIPS'17). However, problems such as SVMs usually admit much smaller input sizes: one is given $n$ data points, each of dimension $d$, and $d$ is oftentimes much smaller than $n$. Furthermore, the SVM program has only $O(1)$ equality linear constraints. This suggests that faster algorithms are feasible, provided the program exhibits certain structures.

In this work, we design the first nearly-linear time algorithm for solving quadratic programs whenever the quadratic objective admits a low-rank factorization, and the number of linear constraints is small. Consequently, we obtain results for SVMs:

* For linear SVM when the input data is $d$-dimensional, our algorithm runs in time $\widetilde O(nd^{(\omega+1)/2}\log(1/\epsilon))$ where $\omega\approx 2.37$ is the fast matrix multiplication exponent;
   
* For Gaussian kernel SVM, when the data dimension $d = O(\log n)$ and the squared dataset radius is sub-logarithmic in $n$, our algorithm runs in time $O(n^{1+o(1)}\log(1/\epsilon))$. We also prove that when the squared dataset radius is at least $\Omega(\log^2 n)$, then $\Omega(n^{2-o(1)})$ time is required. This improves upon the prior best lower bound in both the dimension $d$ and the squared dataset radius.

Published as a conference paper at ICLR 2025
FASTER ALGORITHMS FOR STRUCTURED LINEAR AND KERNEL SUPPORT VECTOR MACHINES
Yuzhou Gu
Zhao Song
New York University yuzhouguênyu. nyu.
Simons Institute for the Theory of Computing, UC Berkeley magic. iinukdda@gma com
Lichen Zhang
MIT CSAIL 1 Lichenzecsail.oo wit.aiiiiii....................m mit.edu
ABSTRACT
Quadratic programming is ubiquitous prototype in convex programming. Many machine learning problems can be formulated as quadratic programming, including the famous Support Vector Machines (SVMs). Linear and kernel SVMs have been among the most popular models in machine learning over the past three decades, prior the deep learning era.
Generally, quadratic program has an input size of ((n)), where n is the number variables. Assuming the Strong Exponential Time Hypothesis (SETH), it is known that no (nn2a(1)) time algorithm exists when the quadratic objective matrix is positive semidefinite (Backurs, Indyk, and Schmidt, NeurlPS' 17). However, problems such as SVMs usually admit much smaller input sizes: one is given data points, each of dimension d, and is oftentimes much smaller than n. Furthermore, the SVM program has only (11) equality linear constraints. This suggests that faster algorithms are feasible, provided the program exhibits certain structures.
In this work, we design the first nearly-linear time algorithm for solving quadratic programs whenever the quadratic objective admits low-rank factorization, and the number of linear constraints is small. Consequently, we obtain results for SVMs:
For linear SVM when the input data is d'dimensional, our algorithm runs time (ndd 1)/2 log(1/c)) where w z 2.37 is the fast matrix multiplication exponent;
For Gaussian kernel SVM, when the data dimension d O n) and the squared dataset radius is sub-logarithmic in n, our algorithm runs in time O log(1/c)). We also prove that when the squared dataset radius is at least 2log n)., then ((n==o(1)) time is required. This improves upon the prior best lower bound in both the dimension d and the squared dataset radius.
INTRODUCTION
Quadratic programming (QP) represents class of convex optimization problems that optimize quadratic objective over the intersection of an affine subspace and the non- nog-niettive orthant'. QPs naturally extend linear programming by incorporating quadratic objective, and they find extensive applications in operational research, theoretical computer science, and machine learning (Kozlov al., 1979; Wright, 1999; Gould & Toint, 2000; Gould et al., 2001; Propato & Uber, 2004; Cor- nuejols & Tütüncü, 2006). The quadratic objective introduces challenges: QPs with general (not necessarily positive semidefinite) symmetric quadratic objective matrix are NP-hard to solve (Sahni, 1974; Pardalos & Vavasis, 1991). When the quadratic objective matrix is positive semidefinite, the problem becomes weakly polynomial-tin solvable, as can be reduced to convex empirical risk minimization (Lee al., 2019) (refer Section C for further discussion).
There classes of QPs with quadratic constraints well. However, in this paper, we focus cases where the constraints linear.

---

Published as a conference paper at ICLR 2025
Formally, the QP problem is defined follows:
Definition 1.1 (Quadratic Programming). Given an x n symmetric, positive semidefinite objective matrix Q. a vector e E RR", and a polytope described by a pair (A R=x", b E R."), the linearly constrained quadratic programming (LCQP) or simply quadratic programming (QP) problem seeks solve the following optimization problem:
min =ER"
r Qx + c.
2
8.1. Ax
= b
x 2 0.
A classic application of QP is the Support Vector Machine (SVM) problem (Boser et al., 1992; Cortes & Vapnik, 1995). In SVMs, dataset .I...., IRd is provided, along with corresponding labels .1....., 3/ E {=1]. The objective is to identify 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 (Muller et al., 2001). 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 linear programming (Cohen et al., 2019; Jiang et al., 2021; Brand, 2020; Song & Yu, 2021), convex empirical risk minimization (Lee et al., 2019; Qin al., 2023), and semidefinite programming (Jiang et al., 2020; Huang et al., 2022; Gu & Song, 2022). In this work, we take pioneering step in developing fast and robust interior point-type algorithms for solving QPs. 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 the hard-margin linear SVMs:
Definition 1.2 (Linear SVM), Given a dataset X E RRad and a collection of labels y1,..., Un each in +1, the linear SVM problem requires solving the following quadratic program:
1
max aER"
1, a 2
(yy" T x
(2)
s.t. a T y 0
a a20. 0.
where 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 Rd such that
min - 11 EER 2
S.t. y: x, b,
Vi E [n].
where b E I is the bias term. Given the solution solution: w" = Dim of yer, At first glance, one might be inclined to solve the primal problem directly, especially cases where d & n, as it presents lower-dimensional optimization problem compared to the dual. The dual formulation becomes particularly advantageous when solving the kernel SVM, which maps features a high or potentially infiite-dimmensional space.
R2, can conveniently convert it to primal
Definition 1.3 (Kernel SVM). Given dataset X E Raxd and positive definite kernel function
K Rd x Rd + R, let R.x= denote the kernel matrix, where KJ,, = K(x,, 2y). With a collection of labels .i......8 Yn each in {=1}. the kernel SVM problem requires solving the following quadratic program:
max
1 a
1
a
T
(wy"
a K)a,
(3)
aER"
2
s.t. a y = 0,
a > 0.
2

---

Published as a conference paper at ICLR 2025
The positive definite kernel function K corresponds to feature mapping, implying that K(z;; 25) ((2,) ((z)))(())) for some Φ R + R². 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 '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 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 single equality constraint of the form a 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 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, oftentimes lack tight theoretical analysis and the estimation of iteration count is usually pessimistic (Chang & Lin, 2011). Theoretically, Joachims (2006) systematically analyzed this class of algorithms. demonstrating that to achieve an c-approximation solution, Ö(c-²B B. nnz(4)) time is sufficient, where B is the squared-radiue of the dataset and nnz(.) denotes the number of nonzero entries. This is subsequently improved in Shalev-Shwartz et al. (2011) 1) with a subgradient-baa method that runs in Ö d) time. Unfortunately, the polynomial dependence on the precision c- makes them hard to be adapted for even moderately small c. For example, when is set to be 10-3 to account for the usual machine precision errors, these algorithms would require at least 10 iterations to converge.
To develop a nigh-precision algorithm with poly log(e-1) dependence instead of poly(e-1), we focus on second-order methods for QPs. A variety of approaches have been explored in previous works, including the interior point method (Karmarkar, 1984), active set methods (Murty, 1988), augmented Lagrangian techniques (Delbos & Gilbert, 2003), conjugate gradient, gradient projection, and extensions of the simplex algorithm (Dantzig, 1955; Wolfe, 1959; Murty, 1988). Our interest is particularly piqued by the interior point method (IPM). Recent advances in the robust IPM framework have led significant successes for convex programming problems (Cohen et al.. 2019; Lee et al., 2019; Brand, 2020; Jiang et al., 2020; Brand et al., 2020; Jiang et al., 2021; Song & Yu, 2021; Jiang et al., 2022; Huang et al., 2022: Gu & Song. 2022: Qin et al., 2023; Song al., 2023). These successes are a result of combining robust analysis of IPM with dedicated data structure design.
Applying IPM to solve QPs with a constant number of constraints is not entirely novel; existing work (Ferris& & Munson, 2002) has already adapted IPM to solve the linear SVM problem. However, the runtime of their algorithm is sub-optimal. Each iteration of their algorithm requires multiplying a d d' n matrix with an n x d matrix in ((d time, where w ~ 2.37 is the fast matrix multiplication exponent (Duan et al., 2023; Williams et al., 2024; Gall, 2024; Alman et al., 2025). Moreover, the IPM requires O(v log((/))) iterations to converge. This ends up with an overall runtime - log(1/) (1/e)). which is super-linear in the dataset size even when the dimension d is small. In practical scenarios where n is usually large, the 1.5 dependence becomes prohibitive. Therefore, it is crucial to develop an algorithm with almost- or nearly-lineal dependence on n and logarithmic dependence on c-!
For linear SVM, we propose a nearly-lineer time algorithm with nigh-precision guarantees, applicable when the dimension of the dataset is smaller than the number of points:
Theorem 1.4 (Low-rank QP and Linear SVM, informal version of Theorem E. Given quadratic program as defined in Definition 1.1, and assuming a low-rank factorization of the quadratic objective matrix Q = UV where U, V E BB there exists an algorithm that can solve the program (1) up to e-error2 in Öm(k + m) mn)lw1)//2 log(n/e)) time.
Specificully, for linear SVM (as per Definition 1.2) with d S n, one can solve program (2) up to e-error in ((nd(++1)/2 log(n/e)) time.
We say algorithm that solves the program up to e-error if it returns approximate solution vector ā whose objective value most more than the optimal objective value.
3

---

Published as a conference paper at ICLR 2025
While learly time algorithm for linear SVMs is appealing, most applications look kernel SVMs as they provide more expressive power to the linear classifier. This poses significant challenge in algorithm design, forming the kernel matrix exactly would require 22(2²) time. Moreover, the kernel matrix could be full-rank without any structural assumptions, rendering our low-rank QP solver inapplicable. In fact, it has been shown that for data dimension d = (log n), algorithm can approximately solve kernel SVM within an error exp( (log"" n)) in time ( assuming the famous Strong Exponential Time Hypothesis (SETH)² (Backurs et al., 2017).
Conversely, long line of works aim to speed up computation with the kernel matrix faster than quadratic, especially when the kernel has certain smooth and Lipschitz properties (Alman et al., 2020; Aggarwal & Alman, 2022; Bakshi et al., 2023; Charikar et al., 2024). 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 Aggarwal & Alman (2022) showed that for dimension d Ollo n) and squared dataset radius (defined as maXi, ||=; - xy ( B = o(log n)), there exists low-rank matrices U, VVE R=xn"()) such that for any vector E R", [I(K UV )all
UV]lal S <||ær||1. They subsequently develop an algorithm to solve the Batch Gaussian KDE problem in +o()) time.
Based on this dichotomy in fast kernel matrix algebra, we establish two results: 1) Solving Gaussian kernel SVM in (n2 log(1/e)) time is feasible when B = ol log
Tog ). and 2) Assuming SETH, no sub-quadratic time algorithm exists for B = 2ilog" n) in SVMs without bias and = 2 n) in SVMs with bias. This improves the lower bound established by Backurs et al. (2017) in terms of dimension d.
Theorem 1.5 (Gaussian Kernel SVM, informal version of Theorem G.7 and G.12). Given a dataset
X E Rn with dimension d and squared radius denoted by B, let K(z;, 2zy) exp(-Ilæ; ||a, 2,11B) be the Gaussian kernel function. Then, for the kernel SVM problem defined in Definition 1.3,
* O(log n), B = ol ( log
log log n. ). there exists an algorithm that solves the Gaussian kernel SVM up to c-error in time O +o(1) log(1/<));
If d= Silog n), B Qlogg n). then assuming SETH, any algorithm that solves the Gaussian kernel SVM without bias term up to exp(-w(log ni) error would require ( time;
Ifd= 22ogg n), B = 2(og n), then assuming SETH, any algorithm that solves the Gaus- sian kernel SVM with a bias term up to exp(-w(log" n)) error would require 2(n² o(1) time.
To our knowledge, this the first Aimost-linear time algorithm for Gaussian kernel SVM even when d log n and the radius is small. Our algorithm effectively utilizes the rank-n factorization of the Gaussian kernel matrix alongside our low-rank QP solver.
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 its algorithmic speedup. For linear SVM, Joachims (2006) offers first-order algorithm that solves its QP in nearly-lincar time, but with a runtime dependence of 2 limiting use in high precision settings. This runtime is later significantly improved by Shalev- descent algorithm. For SVM classification, existing algorithms such SVM-Light (Joachims, 1999), SMO (Platt, 1998). LIBSVM (Chang & Lin, 201 1), and SVM-Torch (Collobert & Bengio, 2001) perform well in high-dimensional data settings. However, their runtime scales super-lined linearly with n, making them less viable for large datasets. Previous investigations into solving linear SVM via interior point methods (Ferris & Munson, 2002) have been somewhat basic, leading to an overall
et al. (2011) to Ō(c-ld) d) via a stochastic subgradient
SETH is a standard complexity theoretical assumption (Impagliazzo et al., 1998; Impagliazzo & Paturi, 2001). Informally, it states that for Conjunctive Normal Form (CNF) formula with m clauses and n variables, there is no algorithm for checking its feasibility in time less than O(e" poly(m)) for 2.
4

---

Published as a conference paper at ICLR 2025
runtime of Of For comprehensive survey on efficient algorithms for SVM, refer to Cervantes et al. (2020). On the hardness side, Backurs et al. (2017) provides an efficient reduction from the Bichromatic Closest Pair problem to Gaussian kereel SVM, establishing an almost-quadratic lower bound assuming SETH.
Interior Point Method. The interior point method, well-established approach for solving convex programs under constraints, was first proposed by Karmarkar (1984) as (weakly) polynomial-ting algorithm for linear programs, later improved by Vaidya (1989) in terms of runtime. Recent work by Cohen et al. (2019) has shown how to solve linear programs with interior point methods in the current matrix multiplication time, utilizing robust IPM framework. Subsequent studies (Lee et al., 2019; Brand, 2020; Jiang et al., 2021; Song & Yu, 2021; Huang et al., 2022; Jiang al., 2022; Qin et al., 2023) have further refined their algorithm or applied to different optimization problems.
Kernel Matrix Algebra. Kernel methods, fundamental in machine learning, enable feature map- pings to potentially infinite dimensions for data points in d dimensions. The kernel matrix, crucial component of kernel methods, often has 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-vechs -vettor multiplication, spectral sparsification, and Laplacian system solving. The study by Alman et al. (2020) introduces comprehensive toolkit for solving these problems in almost-linear time for small dimensions, leveraging techniques like polynomial methods and ultra Johnson- Linionn--nnnnssssuss transforms. Alternatively, Backurs et al. (2021); Bakshi et al. (2023) reduce various kernel matrix algebra tasks to kernel density estimation (KDE), which recent advancements in KDE data structures (Charikar & Siminelakis, 2017: Backurs al., 2018; Charikar et al., 2020) have made more efficient. A recent contribution by Aggarwal & Alman (2022) provides tighter characterization of the low-degree polynomial approximation for the e function, leading to more efficient algorithms for the Batch Gaussian KDE problem. Another line of works aim apply oblivious sketching the kernel matrix, so that spectral approximation can be computed in time nearly- -ineer in the dataset size and a low-rank factorization with rank being the statistical dimension of the problem can computed (Avron et al., 2014: Ahle al., 2020; Zandich al., 2021; Woodruff & Zandieh, 2022). Unfortunately, for the statistical dimension to be smaller than n, one has to add a regularization parameter X, which usually appears when one wants to solve a kerne ridge regression. The work (Song et al., 2021) does not require the parameter X, but obtains runtime improvements for d 2 n, and the sketching dimension is super-linear in n.
2 TECHNIQUE OVERVIEW
In this section, provide an overview of the techniques employed in development of two nearly-linear time algorithms for structured QPs. In Section 2.1, we detail the robust IPM framework, which forms the foundation of our algorithms. Subsequent section, namely Section 2.2 and 2.3, delves into dedicated data structures designed for efficiently solving low-treewidth and low-rank QPs, respectively. Finally, in Section 2.4, we discuss the adaptation of these advanced QP solvers for both linear and kernel SVMs.
Due to the heavill-technical nature, we recommend that in the first read, the audience can skip Section 2.2 and 2.3.
2.1 GENERAL STRATEGY
Our algorithm is built upon the robust IPM framework, an efficient variant of the primal-dual central path method (Renegar, 1988). This framework maintains a primal-dual solution pair (I, E R2" x To understand the central path for QPs, we first consider the central path equations for linear programming (see Cohen et al. (2019); Lee et al. (2019) for reference):
s/t VΦ(x) μr
Ax b,
A y + C,
5

---

Published as a conference paper at ICLR 2025
where is the primal variable, is the slack variable, is the dual variable, o(x)) is self-concordant barrier function, and jl denotes the error. The central path is defined by the trajectory (æ, s) as approaches 0.
In quadratic programming, we modify these equations:
s/t + VΦ(z) = p4
Ar
Qu + A + c,
where Q is the positive semidefinite objective matrix. The key difference in the central path equations for LP and QP the inclusion of the Qr term in the third equation, significantly affecting algorithm design.
Fundamentally, IPM is Newton's 's method in which we update the variables x, y and s through the second-order information from the self-concordan barrier function. We derive the update rules for QP (detailed derivation in Section F.2):
8, tM- -
8y 1(AM-IAT) AM-
6, = t6 -HHM 1PMM
where
H = V2d(T),
M M
P
-AM--AT--AA
where 6g, Ky. 8s. and 6, are the incremental steps for æ, y. 8, and μu, respectively.
The robust IPM approximates these updates rather than computing them exactly. It maintains an approximate primal-dual solution pair (I, 3) R" R" and computes the steps using this approximation. Provided the approximation is sufficiently accurate, it can be shown (see Section F for more details) that the algorithm converges efficiently to the optimal solution along the robust central path.
Therefore, the critical challenge lies efficiently maintaining (7. 3), an approximation to (ar, s). when (x, s) evolves following the robust central path steps. The primary difficulty is that explicitly managing the primal-dual solution pair (xr s) 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 (z, s) are only required at the algorithm's 's conclusion, not during its execution. Instead, we can identify entries with significant changes and update the approximation (x, 8) correspondingly, With IPM's lazy updates, only a nearly-linear number of entries are adjusted throughout the algorithm:
Σ |z(00) ((1-1))lo + ||s(0) - ((--1) = log(//e)
where T T= Öōvá log(1/c)) is the number of iterations for IPM convergence. This indicates that, average, each entry of E and ø is updated log(1/e) times, facilitating rapid updates to these quantities and, consequently, H.
In the special case where Q = 0, the path reverts to the I.P case, with M = tH being a diagonal ma- trix, allowing for efficient computation and updates of - This simplifies maintaining AM-IAT, as updates to -1 correspond to row and column scaling of A. However, in the QP scenario, where M is symmetric positive semidefinite, maintaining the term AM-1A 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 M-1/2 in the robust central path steps. While the matrix Woodbury identity could be considered, falls short when maintaining square root term. Despite these hurdles, we construct efficient data structures for - 1/2 maintenance when Q possesses succinct representations, such as low-rank.
Before diving into the particular techniques for low-rank QPs, we start by exploring the low-treewidth QPs, which could be viewed as a structured sparsity condition. provides valuable insights for the low-rank scenario.

---

Published as a conference paper at ICLR 2025
2.2
SETTING: How TO LEVERAGE SPARSITY
Treewidth is parameter for graphs that captures the sparsity pattern. Given a graph G - (V, E) with vertices and m edges, free decomposition of G arranges its vertices into bags, which collectively form a tree structure. For any two bags X,, Xj, if a vertex is present in both, it must also be included in all bags along the path between X; and Xj. Additionally, each pair of adjacent vertices in the graph must be present together in at least one bag. The treewidth is defined as the maximum size of bag minus one. Intuitively, a graph G with small treewidth implies structure akin to tree. For formal definition, see Definition A. When relating this combinatorial structure to linear algebra, we could treat the quadratic objective matrix Q as a generalized adjacency matrix, where we put vertex nonzero. The low-treewidth structure of the graph corresponds to a sparsity pattern that allows one to compute c column-spant Cholesky factorization of Q. Since M = Q tH and H is diagonal, can decompose M = LL T into sparse Cholesky factors*.
i-th row of Q, and put an edge {us, v2} whenever the entry Qi, is
Under any coordinate update I, M is updated on only one diagonal entry, enabling efficient updates to L. The remaining task to use this Cholesky decomposition to maintain the central path step.
By expanding the central path equations and substituting M = LL", we derive
6² = tM 1/2 P)M-1 8/28,
tM 8, - 1M A (AM (AM--ATT-1 AM AM
tL -L- L- - tL L L---LLLAAALLLLL
6,-
1-L--AAL-T 1011
8. = z6 - H2
- -1 1/25
t6 -L-LL--LL- L- + A (ALT(ALLTT) LAL- -AL---AA--AA L L
Updates to the diagonal of M do not change L's nonzero pattern, allowing for efficient utiliza- tion of the sparse factor and maintenance of L-1 E R"xm and L L-15,, , RR. Terms like (AL-TL-1AA - 1 E R" can also be explicitly maintained.
With this approach, we propose the following implicit representation for maintaining s):
r= == ã + H W/ (h2z h,3 + +zz),
(4)
s s= = + 1/2 c.B., Be, - = (h3, - h13, + ts),
(5)
where 2, W=
Ruxn R. h L E Rax Bz, Bg E IR" Ez, E R.. All quantities except for W can be explicitly maintained. For linear programming. the implicit representation is as follows:
L p R". c, H- 1/28, R", Br. B., Sc, -- 18,
x = 2 + H 1/2 Boco H 1/2 (Bah ez)
$ = + /
c.),
with L.----AH-- maintained implicitly and the other terms explicitly.
The representation in (4) and (5) enables us to maintain the central path step using a combination of "coefficients" h+ h+h and "basis" WT. We need to detect entries of I that deviate significantly from and capture these changes with |/x 221 We maintain this vector using zo + WT (h- + h.2)). Here, W acts as wavelet basis and the vector h + h as its multiscale coefficients. While computing and maintaining W h seems challenging, leveraging column-sparsin of L-1 is possible through contraction with a vector v:
T =(Wu)
= -LI-'
By applying the Johnson-Lindenstraun transform (JL) in place of v, we can quickly approximate |[WT hil2 by maintaining T for JL matrix 4. Similarly, we handle WTh3, h3, by explicitly computing A B. and using the sparsity of L for h32.
'Note that adding non-negative diagonal matrix does not change its sparsity pattern, hence M also retains the treewidth 7.
7

---

Published as a conference paper at ICLR 2025
We focus on entries significantly deviating from zu, the heavy entries of WT (h + hB,). Here, the treewidth-7 decomposition enables quick computation of an elimination tree based on L L- sparsity, facilitating efficient estimation of II(W' (h (h+))))))xc ll2 for any subtree x(v)². With an elimination tree of height (T)) we can employ heavy-light decomposition (Sleator & Tarjan, 1981) for an O(log nhehig tree.
Using these data structures, convergence is established using the robust IPM framework (Ye, 2020; Lee & Vempala, 2021). While the framework is generally applicable to QPs, computing initial point remains a challenge. We propose a simpler objective zi- = arg millzeRr >i=1 @x (2;) with Φ; as the log-barrier function, resembling the initial point reduction in Lee et al. (2019). 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 form of sparsity, allowing for a sparse factorization
M = LL Another significant structure arises when the matrix Q admits a low-rank factorization. Let Q UV where U, V Ra*k and k < n, then M Q tH UV + tH. Although Q has a low-rank structure, M may not be low-rank due to the diagonal matrix being dense. However, in the central path equations, we need only handle M-1, which can be efficiently maintained using the matrix Woodbury identity:
M = -
Given that H is diagonal, the complex term (1 -- can be quickly updated under sparse changes to H 11 by simply scaling rows of U and V. With only nearly-linear number of updates to H H_ the total update time across Ölvā log(1/e)) iterations is bounded by Ö(nk--1- + h")) We modify the (Ir, s) implicit representation follows:
++ H 1/2 + H- + 1-1
(6)
8 s= = + h,, + H 3 + H
where x. RR 1-122 R. : H-1/2U E Roak and H- R with 3z. 3. Rm and Bz, 3, E R. The nontrivial terms to maintain are and h, but both can be managed straightforwardly: updates to 1-1 correspond to scaling rows of U and A TT, and can he performed in total (nkk) and (mm) time, respectively. The key observation is that we never explicitly form M 1/2, hence matrix Woodbury identity suffices for fast updates.
The remaining task 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 's leaves. While a heavy-hitter data structure could lead to improvements in poly-logarithmic and sub-logarithmie 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 when the data dimension d is much smaller than n. However, for kernel SVM, forming the kernel matrix exactly would take O(n²) time. Fortunately, advancements in kernel matrix algebra (Alman et al., 2020; Backurs al.. 2021: Aggarwal Alman, 2022: Bakshi et al., 2023) have enabled sub-quadratic algorithms when the data dimension d is small or the kernel matrix has relatively large minimum entry. Both Alman et al. (2020) and Bakshi et (2023) introduce algorithms for spectral sparsification, generating an approximate matrix K E R= such that (1 11-6. c) K K N (1 + c) K, with K having only O 2n log n) nonzero entries. Alman et al. (2020) achieves this in O(n +o(1) time for multiplicatively Lipschitz kernels when d Olog n). while Bakshi et al. (2023) overcomes limitations for Gaussian kernels by basing their algorithm on KDE and the magnitude of the minimum entry of the kernel matrix, parameterized by Their algorithm for Gaussian kernlls runs in time
Given any node v, we use x(v) to denote the subtree rooted at
8

---

Published as a conference paper at ICLR 2025
O(nd/r3.173t0)) Unfortunately, spectral sparsifiers do not aid our primitives since sparsifier only reduces the number of nonzero entries, but not the rank of the kernel matrix.
Besides spectral sparsification, Alman al. (2020); Aggarwal & Alman (2022) also demonstrate that with d = d = (log r2) and suitable kernels, there exists an O(n O(n!=(1)) time algorithm to multiply the kernel matrix with an arbitrary vector e E IR.. This operation is crucial in Batch KDE as shown in Aggarwal & Alman (2022). Moreover, Aggarwal & Alman (2022) establishes an almost-quadratic lower bound for this operation when the squared dataset radius B = w(log n), assuming SETH. These results rely on computing a rank-n factorization for the Gaussian kernel matrix. The function e can be approximated by low-degree polynomial of degree
log (1/c) g O(max v B log(1/e), log(log (olgg(11e)))) L/e)/B)
for x E [0, B]. Using this polynomial, one can create matrices U, V with rank (2 = no(1) in time On Given this factorization, multiplying with a vector U U(V b) takes On
time. Let K = UVT where K. = f(11=; - 1,)13), we have for any (t, E [n]] x [n].
|f( - xy1l²) exp(- =,, Se
and for any row i E [n],
(Kv), - (Kt) = 1 uy (f(1lz= -
-x, =,) - exp(-11x3 x,112))
j=1
$ (max [[[1 x,,113) exp(-1l || xy 111111111111
jE(n)
S <|lu||a,
using Hölder's inequality. This provides an -ugaranne of the error vector (K K)v, useful for Batch Gaussian KDE, Transforming this -uaraante into spectral approximator yields
(1 - cn) K k K < (1 cn) K.
Setting e 1/n², the low-rank factorization offers an adequate spectral approximation to the exact kernel matrix K.
Given K UVT for U.V E Raxn(()) we can solve program (3) with K using our low-rank QP algorithm in time Oon1+ log(1/e)). This is the first almost-linear time algorithm for Gaussian kernel SVM, even in low-precision settings, as prior works either lack machinery to approximately form the kernel matrix efficiently, or do not possess faster convex optimization solvers for solving structured quadratic program associated with a kernel SVM.
The requirements d = O(log n) and B = ol log Tog n may seem restrictive, but they are necessary, as no sub-quadratic time algorithm exists for Gaussian kernel SVM without bias when d 2(log n) and B = 2logg" rt), and with bias when B = 2(ogg n). assuming SETH, This is based on reduction from Bichromatic Closet Pair to Gaussian kernel SVM, as established by Backurs et al. (2017). Our assumptions on d and B are therefore justified for seeking almost-lineat time algorithms.
We note that in other variants of definitions for Gaussian kernels, one requires an additional parameter called the kernel width, and the kernel function is defined as exp( [ 1 In commonly used heuristics (Ramdas et al., 2015), O(vd), hence we could without loss of generality assuming by requiring the squared radius to be B/d.
CONCLUSION
On the algorithmic front, we introduce the first nearly- linear time algorithms for low-rank convex quadratic programming, leading to nearly-linead time algorithms for linear SVMs. For Gaussian kernel
Additiional requirement: ol log
See Section G for further discussion.

---

Published as a conference paper at ICLR 2025
SVMs, we utilize low-rank approximation from Aggarwal & Alman (2022) when d Ollog n) and the squared dataset radius is small, enabling an almost-liness time algorithm. On the hardness aspect, we establish that when d = (log r2), if the squared dataset radius is sufficiently large (1o n) without bias and 2loog 6 n) with bias), then assuming SETH, no sub-quadratic algorithm exists. As work is theoretical in nature, we do not foresee any potential negative societal impact. Several open problems arise from our work:
Better dependence on k for low-rank QPs. Our low-rank QP solver exhibits dependence of k(=+1)/2 on the rank k. Given the precomputed factorization, can we improve the exponents on k? Ideally, an algorithm with nearly-lined dependence on would align more closely with 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 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 achieve result akin to that of Lee & Sidford (2019), 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 2 n) without bias and B = 12(lo n) with bias. This contrasts with our algorithm, which requires B have sub-logarithmic dependence on n. For Batch Gaussian KDE, Aggarwal & Alman (2022) demonstrated that fast algorithms are feasible for B nl) with no sub-quadratic time algorithms for = (log n2) assuming SETH. Can a stronger lower bound be shown for SVM programs with bias term, reflecting more natural setting?
for Gaussian kernel SVM when B = 22log
ACKNOWLEDGEMENT
Part of the work was done while Yuzhou Gu was supported by the National Science Foundation under Grant No. 1MS-1926686. Lichen Zhang was supported in part by NSF CCF-195521 and NSF DMS-2022448.
REFERENCES
Amol Aggarwal and Josh Alman. Optimal-degree polynomial approximations for exponentials and gaussian kernel density estimation. In Proceedings of the 37th Computational Complexity Conference, CCC 22, Dagstuhl, DEU, 2022. Schloss Dagstuhl-Leibniz- Zenituu fuer Informatik.
Thomas D. Ahle, Michael Kapralov, Jakob Bæk Tejs Knudsen, Rasmus Pagh, Ameya Velingker, David P. Woodruff, and Amir Zandieh. Oblivious sketching of high-degree polynomial kernels. In Proceedings of the 2020 ACM-SIAM Symposium on Discrete Algorithms (SODA), pp. 141-160, 2020.
Josh Alman, Timothy Chu, Aaron Schild, and Zhao Song. Algorithms and hardness for linear algebra geometric graphs. In 2020 IEEE 61 Annual Symposium on Foundations of Computer Science (FOCS), 2020.
Josh Alman, Ran Duan, Virginia Vassilevska Williams, Yinzhan Xu, Zixuan Xu, and Renfei Zhou. More asymmetry yields faster matrix multiplication. In Proceedings of the 2025 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), pp. 2005-2039, 2025.
Haim Avron, Huy L. Nguyen, and David P. Woodruff. Subspace embeddings for the polynomial kemel. In NeurlPS, 2014.
Arturs Backurs, Piotr Indyk, and Ludwig Schmidt. On the fine- tine-grained complexity of empirical risk minimization: Kernel methods and neural networks. Advances in Neural Information Processing Systems, 30, 2017.
Arturs Backurs, Moses Charikar, Piotr Indyk, and Paris Siminelakis. Efficient density evaluation for smooth kernels. In 2018 IEEE 59th Anmual Symposium on Foundations of Computer Science (FOCS), 2018.
10

---

Published as a conference paper at ICLR 2025
Arturs Backurs, Piotr Indyk, Cameron Musco, and Tal Wagner. Faster kernel matrix algebra via density estimation. In Marina Meila and Tong Zhang (eds.), Proceedings of the 38th International Conference Machine Learning, Proceedings of Machine Learning Research. PMLR, 2021.
Ainesh Bakshi, Piotr Indyk, Praneeth Kacham, Sandeep Silwal, and Samson Zhou. Sub-quadratic algorithms for kernel matrices via kernel density estimation. In International Conference on Learning Representation, ICLR 23, 2023.
Aaron Bernstein, Maximilian Probst Gutenberg, and Thatchaphol Saranurak. Deterministic decre- mental sssp and approximate min-cost flow in almost-linear time. In 2021 JEEE 62nd Annual Symposium Foundations of Computer Science (FOCS), pp. 1000-1008. IEEE, 2022.
Bernhard E. Boser, Isabelle M. Guyon, and Vladimir N. Vapnik. A training algorithm for optimal margin classifiers. In Proceedings of the Fifth Annual Workshop on Computational Learning Theory, COLT '92, 1992.
Jan van den Brand, A deterministic linear program solver in current matrix multiplication time. In Proceedings of the Thirty-First Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 20, 2020.
Jan van den Brand, Yin Tat Lee, Aaron Sidford, and Zhao Song. Solving tall dense linear programs in nearly linear time. In Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing, pp. 775-788, 2020.
Jair Cervantes, Farid Garcia Lamont, Lisbeth Rodriguez Mazahua, and Asdrubal Lopez. A com- prehensive survey on support vector machine classification: Applications, challenges and trends. Neurocomputing, 108:189-215: 2020.
Chih-Chung Chang and Chih-Jen Lin. Training v-support vector classifiers: Theory and algorithms. Neural Comput., sep 2001.
Chih-Chung Chang and Chih-Jen Lin. Libsvm: A library for support vector machines. ACM Trans. Intell. Syst. Technol., may 2011.
Moses Charikar and Paris Siminelakis. basshing-base--sstiaaee for kernel density in high di- mensions. In 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS), 2017.
Moses Charikar, Michael Kapralov, Navid Nouri, and Paris Siminelakis. Kernel density estimation through density constrained near neighbor search. In 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS). IEEE Computer Society, 2020.
Moses Charikar, Michael Kapralov, and Erik Waingarten. A quasi-monte carlo data structure for smooth kernel evaluations. In Proceedings of the 2024 Anmual ACM-SIAM Symposium on Discrete Algorithms (SODA), pp. 5118-5144, 2024.
Michael Cohen, Yin Tat Lee, and Zhao Song. Solving linear programs in the current matrix multiplication time. STOC, 2019.
Ronan Collobert and Samy Bengio. Symtorch: Support vector machines for large-scale regression problems. Journal of machine learning research, 1(Feb):143-160, 2001.
Gerard Cornuejols and Reha Tütüncü. Optimization methods in finance, volume Cambridge University Press, 2006.
Corinna Cortes and Vladimir Vapnik. Support-vector networks. Mach. Learn., 1995.
George B Dantzig. Linear programming under uncertainty. Management science, 1(3-4):197-206, 1955.
Timothy A Davis and William W Hager. Modifying sparse cholesky factorization. SIAM Journal on Matrix Analysis and Applications, 20(3):606-627, 1999.
Frédéric Delbos and Jean Charles Gilbert. Global linear convergence of an augmented Lagrangian algorithm for solving convex quadratic optimization problems. PhD thesis, INRIA, 2003.
11

---

Published as a conference paper at ICLR 2025
James R Driscoll, Neil Sarnak, Daniel D Sleator, and Robert E Tarjan. Making data structures persistent. Journal of computer and system sciences, 38(1):86-145) 124, 1989.
Ran Duan, Hongxun Wu, and Renfei Zhou. Faster matrix multiplication via asymmetric hashing, In FOCS, 2023.
Michael C Ferris and Todd S Munson. Interior-point methods for massive support vector machines. SIAM Journal on Optimization, 13(3):883-804, 2002.
Francois Le Gall. Faster rectangular matrix multiplication by combination loss analysis. In Pro- ceedings of the Thirty-Fifth Annual ACM-SIAM Symposium on Discrete Algorithms, S0DA'24, 2024.
Nicholas IM Gould and Philippe L Toint. A quadratic programming bibliography. Numerical NumericalAnallstt Group Internal Report, 1:32, 2000.
Nicholas IM Gould, Mary E Hribar, and Jorge Nocedal. On the solution of equality constrained quadratic programming problems arising in optimization. SIAM Journal on Scientific Computing. 1374)11376-1355, 2001.
Yuzhou Gu and Zhao Song. A faster small treewidth sdp solver. arXiv preprint arXiv:2211.06033, 2022.
Baihe Huang, Shunhua Jiang. Zhao Song, Runzhou Tao, and Ruizhe Zhang. Solving sdp faster: A robust ipm framework and efficient implementation. In FOCS, 2022.
R. Impagliazzo, R. Paturi, and F. Zane. Which problems have strongly exponential complexity? In Proceedings 39th Annual Symposium on Foundations of Computer Science, 1998.
Russell Impagliazzo and Ramamohan Paturi. On the complexity of k-sat. J. Comput. Syst. Sci., 62: 367-375, mar 2001.
Haotian Jiang, Tarun Kathuria, Yin Tat Lee, Swati Padmanabhan, and Zhao Song. A faster interior point method for semidefinite programming. In 2020 IEEE 6lst annual symposium on foundations of computer science (FOCS), pp. 910-918. IEEE, 2020.
Shunhua Jiang, Zhao Song, Omri Weinstein, and Hengjie Zhang. A faster algorithm for solving general 1ps. In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing. pp. 823-832, 2021.
Shunhua Jiang, Bento Natura, and Omri Weinstein. A faster interior-point method for sum-of- squares optimization. In 49th EATCS International Conference on Automata, Languages, and Programming, LIPIcs. Leibniz Int. Proc. Inform., 2022.
Thorsten Joachims. Making large-scale support vector machine learning practical. In Advances in kernel methods: support vector learning, pp. 169. MIT press, 1999.
Thorsten Joachims. Training linear svms in linear time. In Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 217-226, 2006.
William B Johnson and Joram Lindenstrauss. Extensions of lipschitz mappings into a hilbert space. Contemporary mathematics, 26(189-206):.. 1984.
Narendra Karmarkar. A new polynomial- -ime algorithm for linear programming. In Proceedings of the sixteenth annual ACM symposium Theory of computing, pp. 302-311, 1984.
Mikhail K Kozlov, Sergei Pavlovich Tarasov, and Leonid Genrikhovich Khachiyan. Polynomial solvability of convex quadratic programming. In Doklady Akademii Nauk, pp. 1049-1051. Russian Academy of Sciences, 1979.
Yin Tat Lee and Aaron Sidford. Solving linear programs with sqrt (rank) linear system solves. arXiv preprint arXiv:1910.08033. 2019.
Yin Tat Lee and Santosh S. Vempala. Tutorial on the robust interior point method, 2021.
12

---

Published as a conference paper at ICLR 2025
Yin Tat Lee, Zhao Song, and Qiuyi Zhang. Solving empirical risk minimization in the current matrix multiplication time. In Conference on Learning Theory. pp. 2140-2157. PMLR, 2019.
K.-R. Muller, S. Mika, G. Ratsch, K. Tsuda, and B. Scholkopf. An introduction to kernel- -ased learning algorithms. IEEE Transactions Neural Networks, 12, 2001.
Katta G Murty. Linear complementarity, linear and nonlinear programming, volume 3. Citeseer, 1988.
Yurii Nesterov. Introductory lectures on convex programming volume i: Basic course. Lecture noles, 3(4):5, 1998.
Panos M Pardalos and Stephen A Vavasis. Quadratic programming with one negative eigenvalue is np-hard. Journal of Global optimization, 1()):15-2,, 1991.
John Platt. Sequential minimal optimization: A fast algorithm for training support vector machines. MSR-TR-98-14 1998.
Marco Propato and James Uber. Booster system design using mixed-intege quadratic programming. Journal of Water Resources Planning and Management, 130(4):348- -3248 2004.
Lianke Qin, Zhao Song, Lichen Zhang, and Danyang Zhuo An online and unified algorithm for projection matrix vector multiplication with application to empirical risk minimization. In AISTATS, 2023.
Aaditya Ramdas, Sashank J. Reddi, Barnabás Póczos, Aarti Singh, and Larry Wasserman. On the decreasing power of kernel and distance based nonparametric hypothesis tests in high dimensions. In Proceedings of the Twenty-Ninth AAAI Conference Artificial Intelligence, AAAI'15, 15, pp. 3571-3577 AAAI Press, 2015. ISBN 026251 1290.
James Renegar. A polynomial-tine time algorithm, based on newton's method, for linear programming. Math. Program., 40(1-3):55 93, jan 1988.
Aviad Rubinstein. Hardness of approximate nearest neighbor search. In Proceedings of the 50th Anmual ACM SIGACT Symposium on Theory of Computing, STOC 2018, pp. 1260-1268, New York, NY, USA, 2018. Association for Computing Machinery.
Sartaj Sahni. Computationally related problems. SIAM Journal computing, 3(4):262-279, 1974.
Robert Schreiber. new implementation of sparse gaussian elimination. ACM Transactions on Mathematical Software (TOMS), 83)):266-27,, 1982.
Shai Shalev-Shwartz, Yoram Singer, Nathan Srebro, and Andrew Cotter. Pegasos: primal estimated sub-gradient solver for svm. Math. Program., 127, 2011.
Daniel D Sleator and Robert Endre Tarjan. A data structure for dynamic trees. In Proceedings of the thirteenth anmual ACM symposium on Theory of computing, pp. 114-122, 1981.
Zhao Song and Zheng Yu. Oblivious sketching-base central path method for linear programming. In International Conference on Machine Learning, pp.
PMLR, 2021.
Zhao Song, David Woodruff, Zheng Yu, and Lichen Zhang. Fast sketching of polynomial kernels of polynomial degree. In International Conference on Machine Learning, pp. 9812-9823. PMLR, 2021.
Zhao Song, Xin Yang, Yuanyuan Yang, and Lichen Zhang. Sketching meets differential privacy: fast algorithm for dynamic kronecker projection maintenance. In International Conference on Machine Learning (ICML), pp. 32418-32462. PMLR, 2023.
Pravin M Vaidya. Speeding-up up linear programming using fast matrix multiplication. 30th Annual Symposium Foundations of Computer Science, pp. 332-337. IEEE, 1989.
Virginia Vassilevska Williams, Yinzhan Xu, Zixuan Xu. and Renfei Zhou. New bounds for matrix multiplication: from alpha to omega. In Proceedings of the Thirty-Fith Fifth Annual ACM-SIAM Symposium Discrete Algorithms, 20DA'24, 2024.
13

---

Published as a conference paper at ICLR 2025
Philip Wolfe. The simplex method for quadratic programming. Econometrica: Journal of the Econometric Sociery, pp. 382-398, 1959.
David Woodruff and Amir Zandieh. Leverage score sampling for tensor product matrices in input sparsity time. In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato (eds.), Proceedings of the 39th International Conference on Machine Learning, Proceedings of Machine Learning Research. PMLR, 2022.
Stephen Wright. Continuous optimization (nonlinear and linear programming). Foundations of Computer-Aided Process Design, 1999.
Guanghao Ye. Fast algorithm for solving structured convex programs. The University of Washington, Master Thesis, 2020.
Amir Zandich, Insu Han, Haim Avron, Neta Shoham, Chaewon Kim, and Jinwoo Shin. Scaling neural langent kernels via sketching and random features. A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan (eds.), Advances in Neural Information Processing Systems, 2021.
14

---

Published as a conference paper at ICLR 2025
Roadmap. In Section A, present some basic definitions and tools that will be used in the reminder the paper. In Section B, we introduce a few more SVM formulations, including classi- fication and distribution estimation. In Section C, we show convex quadratic programming reduced to convex empirical risk minimization, and therefore can be solved the current matrix multiplication time owing Lee et al. (2019). In Section D and E, prove results low-treewidth and low-rank QPs, respectively. In Section F, we present robust IPM framework for QPs, generalize beyond LPs and convex ERMs with linear objective. In Section G, we present our algorithms for Gaussian kernel SVMs, with complementary lower bound.
be
A PRELIMINARY
A.1 NOTATIONS
For positive integer we use [n] to denote the set f1,2,...,n)) For matrix A, we use A to denote transpose. For matrix A, we define ||Al|p-> == sup, || / When p = =9- = we recover the spectral norm.
We define the entrywise n--ormm of a matrix A as ||Al|p :- ( |A
For any function f + N and n E N, we use (n)))) to denote (f poly log f(n)). We use 1{E) to denote the indicator for event E, i... if E happens, 1(E) and otherwise it's 0.
A.2 TREEWIDTH
Treewidth captures the sparsity and tree-like structures of graph.
Definition A.1 (Tree Decomposition and Treewidth). Let G = (V, E) be graph, tree decomposi-
tion G is tree with vertices, and sets JI...., Jo S V (called bags), satisfying the following properties:
* For every edge (u, v) E E, there exists [b] such that u, Jz:
For every vertex
V. (j [M] : Jj} is non-empty subtree of T.
The treewidth of G is defined as the minimum value of max{|J5] ; [b]] 1 over all tree decompositions.
near-optimal tree decomposition graph can computed in almost linear time.
Theorem A.2 (Bernstein et al. (2022)). Given a graph G, there is an O(m2+((1) time algorithm
that produces a free decomposition of G of maximum bag size O(T log" n). where T is the actual (unknown) treewidih of G.
Therefore, when = 911) we can compute an ((r)-size tree decomposition in time (lm)())) which is negligible in the final running time of Theorem D.I.
A.3 SPARSE CHOLESKY DECOMPOSITION
In this section we state few results on sparse Cholesky decomposition. Fast sparse Cholesky decomposition algorithms are based on the concept of elimination tree, introduced in Schreiber (1982).
Definition A.3 (Elimination tree). Let G be an undirected graph on n vertices. An elimination tree T is a rooted tree on VG)) together with an ordering æ of V (G) such that for any vertex v, its parent is the smallest (under π)) element such that there exists path P from v u, such that ((w) S T(v) for all w E P u.
The following lemma relates the elimination tree and the structure of Cholesky factors.
Lemma A.4 (Schreiber (1982)). Let M be a PSD matrix and T be an elimination tree of the adjacency graph of M (i.e., (i.j) F(G)) iff Majj 0) together with elimination ordering
15

---

Published as a conference paper at ICLR 2025
Let P be the permutation matrix Ps,. 1{v T()).. Then the Cholesky factor L of PMP (i.e.,
PMP LL") satisfies Lij * 0 only if (i) is ancestor of (()).
The following result the current best result for computing a sparse Cholesky decomposition. Lemma A.5 ((Gu & Song, 2022, Lemma 8,4)). Let M E R be a PSD matrix whose adjacency
graph has treewidth т. Then we can compute the Cholesky factorization M = LLL in O(nT*11) time.
The following result is the current best result for updating a sparse Cholesky decomposition.
Lemma A.6 (Davis & Hager (1999)). Let M E Rax= be a PSD matrix whose adjacency graph has treewidth т. Assume that we are given the Cholseky factorization M LL Let w E IR" vector such that M + www has the same adjacency graph M. Then we compute Ar, E gaxx such that L + Ag the Cholesky factor of M + www in ((T)) time.
Throughout our algorithm, we need to compute matrix-vector multiplications involving Cholesky factors. We use the following results from Gu & Song (2022).
Lemma A.7 (GGu & Song, 2022, Lemma 4.7)). Let M E Rexx be a PSD matrix whose adjacency graph has treewidth Assume that we given the Cholseky factorization M LL Then have the following running time for matrix-vector multiplications.
(i) For o R", computing Lv, LTd, L L, T takes O(nT) time.
(ii For o E RR", computing Lv takes O(1vo|OOT) time.
(iii) For o E R", computing - takes O(IL= time.
(iv) For v E IR", ifu is supported on path in the elimination tree, then computing L-1y takes O(T²) time.
(v) For v E RR", computing WT v takes (nT) time, where W = 1--1 with H E RRxx is a non-negative diagonal matrix.
Lemma A.8 ((Gu & Song, 2022, Lemma 4.8)). Let M E Rax* be a PSD matrix whose adjacency graph has treewidth т. Assume that we are given the Cholseky factorization M LLT, Then we have the following running time for mutrix-vector multiplications, when only need result for subset of coordinates.
(i) Let S be a path in the elimination tree whose one endpoint is the root. For v E IR", computing (L-Tv)s takes O(72) time.
(ii) For U E IR", for E [n], computing (WTx) v); takes O(7²) time, where W = 1-2H1/2 with H. E Rnxn be a non-negative diagonal matrix.
A.4 JOHNSON-LINDENTR..N
LEMMA
We recall the Lindess--indenssrruuus lemma, a powerful algorithmic primitive that reduces dimension while preserving E2 norms.
Lemma A.9 (Johnson & Lindenstrauss (1984)). Let e E (0,1) be the precision parameter. Let
5 E (0, 1) be the failure probability. Let A Rmxn n be real matrix. Let r = c loggmm//)). For
R E R whose entries are i.i.d N0, }),, the following holds with probability at least 5:
(1 - 1-aalla llz I Ra,ll2 ll S (1 + c)lla,12, Vi E [m],
where for matrix A, a. denotes the i-th row of matrix A E mx
A.S HEAVY-LIGHT DECOMPOSITION
re-balance binary tree with height Ollog n). Lemma A.10 (Sleator & Tarjan (1981)). Given a rooted tree T with n vertices, we can construct in ((n) time an ordering of the vertices such that (1) every path in T can be decomposed into Olog n) contiguous subseqeuences under n, and (2) every subtree in is a single contiguous subsequence under
Heavy-light decomposition is useful when
16

---

Published as a conference paper at ICLR 2025
B SVM FORMULATIONS
In this section, we review a list of formulations of SVM. These formulations have been implemented in the LIBSVM library Chang & Lin (2011).
Throughout this section, we use Φ : IR4 + RR to denote the feature mapping, to denote the associated kernel function and K R denote the kernel matrix. For linear SVM, is just the identity mapping. We will focus on the dual quadratic program formulation as usual. We will also assume for each problem, dataset X IRaxd is given together with binary labels RR Let Q = (yy ) K.
B.1 C-SUPPORT VECTOR CLASSIFICATION
This formulation also referred as the seft-margin SVM. It can be viewed as imposing regularization on the primal program to allow mis-classification. Definition B.1 (C-Support Vector Classification). Given a parameter > 0, the C-support vector classification (C-SVC) is defined
1 T QQ
1 a
aER
2
s.t. a T = O 1,
B.2 P-SUPPORT VECTOR CLASSIFICATION
The C-SVC (Definition B.1) penalizes large values of a by limiting the magnitude of it. The v-SVC (Definition B.2) turns 1 a from an objective into a constraint on ε norm.
Definition B.2 (v-Support Vector Classification). Given parameter > 0, the v-support vector
classification (v-SVC) is defined as
1
min OEE"
Qu 2a
s.t. a y =0
1 a = v,
0 a 1 1
One can interpret this formulation as to find a vector that lives in the orthogonal complement that is non-negative, each entry is at most and its l norm is Clearly, must have (0, 1]. More specifically, let k4 be the number of positive labels and k_ be the number of negative labels. shown by Chang & Lin (2001) that the above problem is feasible if and only if
2 min (k_ k+}
V 5
n
B.3 DISTRIBUTION ESTIMATION
SVM widely-used for predicting binary labels. can also be used estimate the support of high-dimensional distribution. The formulation is similar to v-SVC, except the PSD matrix Q is label-less.
Definition B.3 (Distribution Estimation). Given parameter 0, the y-distribution estimation problem is defined as
min oER"
2 a T Ka
s.t. 0 a < 1
a
17

---

Published as a conference paper at ICLR 2025
B.4 E-SUPPORT VECTOR REGRESSION
In addition to classification, support vector framework can also be adapted for regression.
Definition B.4 (c-Support Vector Regression), Given parameters c, c,C 0, the c-support vector regression (e-VVR) is defined as
- min 3a--a') (la-aa)"(aaa'))-
a²) + u gT (a= - -a))
a,a Elt
s.t.
(a a²) =0,
0 03a a C 1n
0 0 SC 1n.
B.5 U-SUPPORT VECTOR REGRESSION
One can similar adapt the parameter to control the & norm of the regression. Definition B,5 (v-Support Vector Regression), Given parameters v, C > O, the v-support vector regression (v-VVR) defined as
2 (a a²) aaTTa(aaa a("r-a)) + (a(aa- a')
1
min 0,a"
s.t. 1 (a a²) =
1 (a a²) CV,
C
0 <
1
C
0 0 <
1p
B.6 ONE EQUALITY CONSTRAINT
We classify C-SVC (Definition B.1), e-SVR (Definition B.4) and p-distribution estimation (Defini- tion B.3) into the following generic form:
1 T Qa +p
min oER
2 a
s.t. a a =
0 0SaSC. 1,.
Note that C-SVC (Definition B.1) and distribution estimation (Definition B.3) are readily in this form. For c-SV (Definition B.4), we need to perform a simple transformation:
Set â a 2 then it can be written
a*
1
Q
-Q
el, + T
min 6e
-a
T
Q
Q
El,
y
â
1,
0
-1
0 à <
C
12n
B.7 TWO EQUALITY CONSTRAINTS
Both VVC (Definition B.2) and v-SV (Definition B.5) require one extra constraint. They can be formulated as follows:
1 -a
min agR
QQ + p
18

---

Published as a conference paper at ICLR 2025
s.t. 1 a A y gaa a A2.
0 0 a a 1,
For v-SVR (Definition B.5), one can leverage a similar transformation as c-SVR (Definition B.4). Remark B.6. All variants of SVM-related quadratic programs can all be solved using our QP solvers for three cases:
Linear SVM with n >> d, we can solve it in ((nd(w+1)/2 log/:///////////)))// log(1/e)) time;
Linear SVM with small treewidth decomposition with width on XX"", we can solve in O(r log(1/<)) time;
Gaussian kernel SVM with = (loggn) and = o( (nn++(()) log(1/<)) time.
a log log we can solve in
Even though solvers have relatively bad dependence the number of equality constraints, for all these SVM formulations, at most 2 equality constraints are presented and thus can be solved very fast.
C ALGORITHMS FOR GENERAL QP
In this section, 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 Lee et al. (2019).
C.I LCQP IN THE CURRENT MATRIX MULTIPLICATION TIME
Proposition C.1. There an algorithum that solves LCQP (Definition 1.1) up to e error in Ö((r n 2.5- + n2+1/0) log(1/e)) time, where w S 2.373 the matrix multiplication constant and 0.32 is the dual matrix multiplication constant.
Proof. Let Q = PDP be an eigen-decomposition of Q where diagonal and P is orthogonal. Let = p-1g. Then suffices to solve
1 TT DE c T PT
min
2
s.t. APт =
PF 2
By adding non-negative variables and constraints PT 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 into two. 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) where Q diagonal. Let y Qi. be the i-th element on the diagonal. Then the QP is equivalent to the following program
min x +x+q"
s.t. Ar b
4 2 प o4
* 0
Note that the set ((1,, 4:) E IR:
4
| is convex set. So apply Lee et al. (2019) here
with variables, each in the convex set { b) T22
0,b }
19

---

Published as a conference paper at ICLR 2025
C.2 ALGORITHM FOR QCQP
Our algorithm for LCQP in the previous section can be generalized to quadratically constrained quadratic programs (QCQP). QCQP is defined follows.
Definition C.2 (QCQP), Let Qo. .... Qm R be PSD matrices. Let go..... 4m E R Let
r E R" Let E [ b E R The quadratically constrained quadratic programming (QCQP)
problem asks the solve the following program.
1
min ZER"
x
Qux + x
2
1
s.t.
Qix + q x
E [m]
2
Ax
b
x 2 0
Proposition C.3. There is an algorithm that solves QCQP (Definition C.2) up to error in Ó((mn)) + (mm))2.5-a? + (mm)2+1107 log(1/e)) time, where w I 2.373 is the matrix multi- plication constant and 2 0.32 is the dual matrix multiplication constant.
Proof. We first rewrite the program as following.
min
To
s.t. Q + q. g
Y4 S i
Ax x 20
Write Qi = P.D,P be an cigen-decoopositiin of Q; where D; is diagonal and Pa is orthogonal.
Let T; E R" be defined as _; == P.'aa x. Then we can rewrite the program as
min
-r
8.1.
2 r
D + 4, P.u. +550
v4 S < m
Ax =
T. Pi- 1
For 0 < i m and j [m]. define variable tijj E R. Then we can rewrite the program as
min
ro
s.t.
Σ
Disi + P. +r 50
44 i $
jE(m
Ar = b
G; P,-'x
t. 2 o
2 0
We can consider (Zi.j, ligs (oa.)oossmmee as (m 1)n variables in the convex set {(a b): :bbbbb> b> f Then we can apply Lee et al. (2019).
D ALGORITHM FOR COW-TREEWIDTH QP
In this section we present nearly-linea time algorithm for solving low-treewidth QP with small number of linear constraints. We briefty describe the outline of this section.
20

---

Published as a conference paper at ICLR 2025
In Section D.I, present the main statement of Section D.
In Section D.2, we present the main data structure CENTRALIAAEEMALTEENNANCC
In Section D.3, we present several data structures used in CENTRALPATHHMANTENAAAC including EXACTDS (Section D.3.1), APPROXDS (Section D.3.2), BATCHSKETCH (Sec- tion D.3.3), VECTORSKETCH (Section D.3.4), BALANCEDSKETCI (Section D.3.5).
In Section D.4, we prove correctness and running time PENTRALPATIMAAINENNNCCC data structure.
In Section D.5, we prove the main result (Theorem D.1).
D.1 MAIN STATEMENT
consider programs of the form (16), i.c.,
1
min
Qu + cTx
sER*
2
s.t. Ax = b
3; E Kr
E [n]]
where Q E Sheet C C I Mtos, A E Muneer b be Rm, m Kr C I is convex set. For simplicity, we assume that n, 01) for all [n].
Theorem D.1. Consider the convex program (16). Let di KG + IR be u,-seff-ooncordmm barrier for all E [n]. Suppose the program satisfies the following properties:
Inner radius v: There exists E Ti such that Az and B(2,r) E K.
Outer radius R: K. S B(0, R) where T
Lipschitz constant L:: [Q[|2-<2 L, ||e||2 S T..
Treewidth т: Treewidth (Definition A.1) of the adjacency graph of Q is at most T.
Let (wi)se(n) IR2, and Lie(n) wiys. Given any 0 S 1. we can find an approximate
solution IC satisfiying
1 T Qx cc
{ Tu +ELR(R+ 1),
min Ar=b,zeR
2
|Ax - bll1 S 3e(R]|AAll + (|b|12).
expected time
Ō((Vkn-17 11/2 + (og(RR((re)) =n(r2 + Tm? 2/2 + rm + m m22////
When max v4 O(1, wi m ((7 the running time simplifies to ((nrlw+1)) 11/2 222 1/2 log( (((rr)))
D.2 ALGORITHM STRUCTURE AND CENTRAL PATH MAINTENANCE
Our algorithm based on the robust Interior Point Method (robust IPM). Details of the robust IPM will be given in Section F. During the algorithm, we maintain primal-dual solution pair (x, 8) E Reot x Reat on the robust central path. In addition, we maintain a sparsely-changing approximation (Ir, Meor Reat to (I, 8). In each iteration, we implicitly perform update
I I x x 7B
1/2 [I -P 1/2 w.zl
w.z,t
+ 16, - Hm. B-1f2
P Pw...E B 1/12
W.F.
where
H n (2
(see Eq. (24))
21

---

Published as a conference paper at ICLR 2025
B W, Q wH
(see Eq. (25))
B-1/2 A P n,,,, w. w.z.! w. W.F.
(AB
A
(see Eq. (26)
and explicitly maintain (2, such that they remain close to (x, 8) in lx -isttacce
This task is handled by the PENTRAI.PATMMAITTEN structure. The robust IPM algorithm (Algorithm 19, 20) directly calls it in every iteration.
data structure, which our main data
The
data structure (Algorithm 1) has two main sub data structures, EXACTDS (Algorithm 2, 3) and APPROXDS (Algorithm 4, 5). EXACTDS is used to maintain (a, s), and APPROXDS used maintain (7, 8).
Algorithm Our main data structure for low treewidth QP solver.
1: data structure CENTRALPATHMAINTENAN
Theorem D.2
2: private: member
3: EXACTDS exact
Algorithm 2, 3 > Algorithm
4: APPROXDS approx
5: E N
6: end members
7:
procedure INITIALIZE(Z, E RR R4, (0,
8: exact. S,
D Algorithm
9:
e 0
10:
W t- Vmaxi N t- Vílog log asR
or
11:
+ (7 fm²)- ( + mm 11/1
12:
Eapx,* f- E, Co f- 20x, dapx f
13:
Eapx, t- E t, So t- 3at
14:
approx. INTTIALIITTILIIEE((, 8, s, h Co. H.? H H
1/2 Bo, Bear
3z. Ba, q, &eexact, Capx,z, fapx.,; Sapx)
Algorithm 4. P.Parmmeters from x to 3g come from exact. &exxact is pointer to exact
16: end procedure
17: procedure MULTIPLYANDMOVEIC E R4)
18:
e 1+1+1
19:
if 7 t| > & or > then
20:
Jr f- exact. OuuTCO)
Algorithm
21:
INITIALIZE(E,
end if
23:
Be, 3, Bz, Bo, Be,,B,,,
+ exVEE,
A Algorithm
24:
Se, 8z appox.MOVEANDQUERYY Boo Bc, , Be, Bo)
Algorithm
25:
òn 57. Se,, S,,, /11/220 8 11/22g 8c, & exact. UPDATE( 5z)
A Algorithm
26:
approx. CEBTOX.UPDATE(555 (5r, 5n of de de,, 8 H2/32 0H 2/2333 Sc,)
& Algorithm
27:
end procedure
28: procedure OUTPUT()
29: return exact. OUTTUT)
Algorithm
30: end procedure
31: end data structure
Theorem D.2. Data structure CENTRAL WWWAAAAAAAMMMAANNAAAAM (Algorithm 1) implicitly maintains central path primal-dual solution pair (x, 8) E Rurto x Burtet and explicitly maintains its approximation (Ie 8) E Rear x Teoss using the following functions:
INITIALIZE(Z E RRae E R. to E Ryo.c E (0, 1)): Initializes the data structure with initial primal-dual solution pair (x, 8) E TRos x IR Moos, initial central path timestamp to Roo in 0 + mw- time.
22

---

Published as a conference paper at ICLR 2025
MULTIPLYANDMOVE(t Rzo): It implicitly maintains
I - I -1B11/2 (1
1) W,,,t (
w.z,t
8t Ed p. H B 1/22
w.z.
where Hw,.. B B P.ra. are defined in Eq. (24)(25)(66) respectively, and some timestamp satisfying [ 5- ct.
It also explicitly maintains (I, mo X such that | - 24IT S and 3 silla. S tēw; for all i [m]] with probability least 0.9.
Assuming the function is called at most N times and decreases from mmx to bminn the total running time
Ö((Nn- + log(tmax//mmnn) n(r2m rm?) (1"0((/w1111 Tm + m= 1 122).
OUTPUT: Computes (Ir, 8) E IR x Bres exactly and oulputs them in (nTm) time.
D.3 DATA STRUCTURES USED IN CENTALPATTIMANNTEAAN
In this section we present several data structures used in CENTRALPATEMAANTTNNANCCCC including:
EXACTDS (Section D.3.I): ): This data structure maintains implicit representation of the
primal-dual solution pair (x, s). This directly used by
APPROXDS (Section D.3.2): This data structure explicitly maintains an approximation [z, 3) of (I, s). This data structure is directly used by CENTRALPPTHMMINTENAAN
BATCHSKETCH (Section D.3.3): This data structure maintains sketch (r, s). This data structure is used by APPROXDS.
VECTORSKETCH (Section D.3.4): This data structure maintains sketch of sparsely- changing vectors. This data structure used by BATCHSKETCHI
BALANCEDSKETCH (Section D.3.5): This data structure maintains sketch of vectors of form W u, where is sparsely-changing. This data structure is used by BATCHŠKETCH.
Notation: In this section, for simplicity, we write By for B and L for the Cholesky factor of Bz, i.c. B = LIL
D.3.1 EXACTDS
In this section we present the data structure EXACTDS. It maintains an implicit representation of the primal-dual solution pair (x. 8) by maintaining several sparsely-changing vectors (see Eq. (8)))). This data structure has similar spirit as EXACTDS in Gu & Song (2022), but we have different representation from the previous works because are working with quadratic programming rather than linear programming.
Theorem D.3. Data structure EXACTDS (Algorithm 2, 3) implicitly maintains the primal-dual pair (z, 8) E I Teos x RRoe computable via the expression
x = + H 1/2 (1,3, -13, +
(8)
n.
= 3 + 2122 CaBe, - H2 (h,, - hh3, Ex);;
(9)
w.
where 3,3 Rsoe, W - 112 a. mo Xnsos , h 1 E Ramt, Co w.z
1122
E Boot
Bz. B2y Bc, RR h= A= AT
T xm, Bz. Bg Tmmm Cr, E Dtot
The data structure supports the following functions:
INITIALIZE(I, s. 2.3 R 7 Ryo): Initializes the data structure in Ö(n +nom- nm="11) 1) time, with initial value of the primal-duu ddual pair (x, s), its initial approximation (xr, s), and initial approximate timestamp T.
23

---

Published as a conference paper at ICLR 2025
MooE()): Performs robust central path step
2 マ + 5/ - A AAA,
6,s,
(10)
s s + f6 - B5 + A - 8 in ((m)) time by updating its implicit representation.
(11)
UPDATE(ŌT, 8z R: Updates the approximation pair (I,3) to (Inew I + 85 R Fuew = = 8z E RRtoot) ) in (((=m + Tm² ([] + (102[10)) time, and output the changes in variables 6 81/2 0i Ög_. 0g. 2/.0 deg 1- 0z,, 0z, 6e,-
Furthermore, h, Eg change in ( + (|021l0)) coordinates, changes in ((Tm(1]ōzllo ((o=1lo)) coordinates, and 2 H 1/1 c, change in O + (| coordinates.
OUTPUT(): Output x and s in Ó(nmm) time.
QUERYa(i E [n]): Output x4 in m(7²m) time. This function is used by APPROXDS.
QUERYS(i E [n]): Output 8; in m(72m) time. This function is used by APPROXDS.
Proof of Theurem D.3. By combining Lemma D.4 and D.5.
Lemma D.4. EXACTDS correctly maintains an implicit representation of (x,, 8), iee, invariant
x= E H W.
1 TT (h13g - hB + cz),
+ 1/2
waz CoBe. - H1/22 (h33 - hB3 + (2),
w.
= 1
C. = 1/23
c-
h= I A
w. pm
h hTh,
u u= T
Σ 1E|n]
Wi cosli" (oshii(~77ll 7, (,00," 8,)), Wi
8 a /2 8 8,
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 INITIAL- IZE.
MOVE:
By comparing the implicit representation (8)()) and the robust central path step (10)(11), we that MOVE updates (z. s) correctly.
UPDATE: We would like prove that UPDATE correctly updates the values of h, cas h, u, oz, 8, while preserving the values of (z,s). 8).
First note that HT,,, Bz, La are updated correctly. The remaining updates are separated into two steps: UPDATEh and UPDATEh.
Step UPDATE4; The first few lines of UPDATE/ updates `े and 6, correctly.
We define H Hw, + AHm... Bnew Bx + AB_- Lew
Lo + AL. Op new
=z 08g
and Immediately after Algorithm 3, Line 26, we have
8N + L 11 d3 (Lg - 10 L h 6n -,
= =(1=² (LT A Dr. L= 1= new
L new nnow
C. + 80 1122 8, + 1/122 + 68, + H 1122 8b m.z.
24

---

Published as a conference paper at ICLR 2025
Algorithm 2 The EXACTDS data structure used in Algorithm 1.
data structure EXACTDS
Theorem D.3
2: members
3:
2, E R 7 E H..Z, Bz, L R Poot X Airot
4: 2, Ex, Ez: E Ruos h E R Be, So, Bc, E IR, Be, 3, E R
5:
TH E Rmxmm u E Im O &ER,O, E
6: EEN
7: end members
8: procedure INITIALIZE(Z, 8, S, P 3 EE R+,)
9:
It 2+2,E-8,1+1
10:
2, 0, 0. + 0, t- 0, Br t- 0, B, t 0, Sc, t 0
B2 Hw, B. Q w.
Compute lower Cholesky factor Lg where Li B3 INITIALIZEh(Z, H Hu,7, Lz) end procedure
13:
15:
procedure INITIALIZEN(T, 3 Ro HW, LE
16:
for E [n] do
(ōw); f- sinh Tc(3,3.7)
17:
(I,
7(3,,,11)
18:
+ ā a w cosh 2 ( (,
19: end for
20: h Lg Ton. h t Lg AT C. + H 1/2g
ü+ TT, u Th
22: end procedure
23: procedure MOVE()
24:
B₂ + 32 + (a) 112
B + (α) 1/2 ū u
BG, 8, + -1/2
-12
27:
B. - 3, +7 (a) (a)-1/2
28:
B. B. +7. (a)- l1
return 8z, 3,, Bo,, 3n.3. 3s
30: end procedure
31: procedure OUTPUT()
32: return ~ + H w
1/2 WT (hgr -h3, +x),8++ Ho.... - 2/ (h3, - -h3. e.)
33: end procedure
34: procedure QUERY E [n])
35:
return T4 + H
H-1/2 (...2. (W (hB, -1,,, + ex)));
36: end procedure
37: procedure QUERYs(i [n])
38: return 3; + 1/2
(i,,( Ca,s + 112 (1,) (WT333 113,-33,
-h3, + e,));
39: end procedure
end data structure
(Hnew) )
- 1/2 1/2 new
h 53 AT (Iz + ( AT
= - (Tog + A AL, T-11) AAT
Lew A
So h, C., h are updated correctly. Also
ū 0 = TT + 87 (h + = ( T
(h+ 8z). u + 8, = 7v h + o. (h 6n) + on = (h + 87,) (h+ dx).
25

---

Published as a conference paper at ICLR 2025
Algorithm 3 Algorithm 2 continued.
1: data structure EXACTDS
Theorem D.3
2: procedure UPDATE(ŌZ dz E 1
3: AH... < V2 (I + 85) Hw, D AHH non-zero only for diagonal blocks (i, for
which dz,á # 0
4: Compute ALr where (Lz + (LT) + AL=) = Bz + H.H
5:
UPDATEh( 8z, AHo.: ALT)
6:
UPDATEW( Hu..: AL,)
7:
-+ 8z, 3 1--34 15.33430 3
8:
HE. t Hw, + AH BJ + B7 + tD Ha, La t- Lo + A
9:
return on, 6x dea Se, H1 1= dc,
end procedure
11: procedure UPDATEh( 0g E R A A. Aer E IR Reos Xfott)
12:
S - E [n] δzti + Or 8,,i + 0}
13:
0g +0
14:
for
i S do
15:
Let 7/- (T 3, t), n = 7/ + 8z 3 8z, f).
- 6g, 3 + 8,,
16:
a ā wi coshi" 2( wy 1;) + w; cosh 2 2 Aneww
83 t a sinh 7/ ) ylow /. 8 p.
18:
end for
L - 6/P (L + A A Lr(h (h - ),,))
20:
dc, t A 11 1/2 (5p- + 8 + H w.z 83
1/2
21:
8
T + A...)) -
22:
dg -- Se,
23:
de, 6n,32 + 0g
24:
d + - 5,,3, + 87,30
25:
8n 8 (h + 8%) + 0ū
h
26:
8n - h
(h + 8n) + on
27:
8, 4 8, + 6, h t- + Sn h +- Te + 8 6z t ex + &c, cs + to + 8,, + ū + 8z.
u. - u uuuu 6n 28: end procedure
29: procedure UPDATEW(A Hwzz Aiz E R
30:
dc, T I (h03 h3, + (2)
31: de, A (h33 h3g + Ex)
32: Ex 6 + de, Cs + C, + dc.
33: end procedure
34: end data structure
So ū and maintained correctly. Furthermore, immediately after Algorithm 3, Line 26, we have
(i + L TT new 3, - new Bg + +emew) - (ā + LE T (h3, - - + ex))
I T==T(bne (6n,22 5g,3, + +
Therefore, after UPDATEh finishes, we have
x x= LE (h32 - h32 + cx).
For s, we have
neew + (H 1/2 new Bo, 1 T [h 3g
keew 3. + new ))
(3 H1122 Cs, T (113 T. + ex))
26

---

Published as a conference paper at ICLR 2025
83 8g0e, 0.
-8,33 16 de,)
(ôh33
Therefore, after UPDATEh finishes, we have
$ = (Hw 1/2 ) coße. L (h.3. ( -h3, + c²).
So and are both updated correctly. This proves the correctness of UPDATE/. Step UPDATEW: Define c == Ez + 0 cee = c + de,- Immediately after Algorithm 3, Line 31, we have
( + (
1e (1.32 - h3g + (zw)) - (8 L (h3 - h3. + ex)
((T4)))
1 (11 h13z + 6z) + (I.2 ) 1
=0,
(8 (Hw,))
1/2 Bc, (I.2 new -Th,, h,3, + 6 new)
(5+ (H
1/2 1/2cc/////// CoBc, I (h3,
T - h3, + (x))
(-(L
Lg T (h, h.33 + (x) (Lz new TT 8e,
Therefore, after UPDATEW finishes, we have
= 2 + ( - (h32
- 1 h3, + En),
+ " 1/2 cg.3. - (L - (hh3, - -Ti3g + ea)
So and are both updated correctly. This proves the correctness of UPDATEW.
Lemma D.5. We bound the running time of EXACTDS as following.
(i) EXACTDS. INITIALIE (Algorithm runs in Ö(nr"~1 + nrm nmm-11 time.
(i)) EXACTDS.Move (Algorithm 2) runs in Ö(m²) time.
(iii) EXACTDS.OUTPUT (Algorithum 2) runs in Ö(TT)) time and correctly outputs (x, s).
(iv) EXACTDS. QUEERY and EXACTDS.QUERT (Algorithm 2) runs in m(72m) time and returns the correct answer.
(v) SEACTDS.UPDATES (Algorithm 2) runs in Ö((rmm Fm3) ((=(1))) time. Further- more, h, Ez. Ca change in ( flo + (())l))) coordinates, h changes in ((mm(1)5zll ()ol)))) coordinates, and H=|²3
1222 H 1/2g cs change in O(1|ž=lo (16511o) coordinates.
Proof.
Computing Ly takes Ö( time by Lemma A.5. Computing h and takes ( by Lemma A.7(i)? Computing ū and u takes Tmat (m, m) O(nm-=1) time. All other operations are cheap.
(ii) Computing 7 takes Ö(m²) time. All other operations take (mm) time.
(iii) Running time by Lemma A.7(v). Correctness by Lemma D.4.
(iv) Running time by Lemma A.8(i). Correctness by Lemma D.4.
(v) Computing AL, takes Ō(7 (z)1)) |lo) time by Lemma A.6. It easy that nnz( (A..)) (0l)))) and unz(Ar, ) = ( 2| (x It remains to analyze UPDATEh and UPDATEW. For simplicity, write 8x11o l + 6 in this proof only.
Herre we compute h by computing hi = (At.)) T for i E [m] independently. Using fast rectangular matrix multiplication is possible to improve this running time and other similar places. We keep the current hounds for simplicity.
27

---

Published as a conference paper at ICLR 2025
UPDATEh: Updating a and computing 65, takes Okk) time. Also, 1165, O(k).
Computing 6д takes k(7*k) time by Lemma A.7(). Also, òn is supported on (
paths the elimination tree, so ||x1|lo ((Tk). Similarly we see that computing 0%
take mk) time and nnz((;,) Ö(Tmk).
Computing de, and 63 takes ()7²k) time and || ||ó3|lo = O(k).
Computing 6z, and de, takes ((Tmk) time after computing 6, and og. Furthermore. || ll,, | (
Computing 6g takes maa (m, Tk, m) ((m??k) time. Computing 6, takes Ö(rmk) time.
UPDATEW: To compute de, and de,, we first compute Lh h32 + ex) and Ti= (,,3--h,,+ee +c,), where S S [Mtot] is the row support of Ar. which can be decom- posed into at most (110)))) paths. This takes O(72m ( time by Lemma A.8(i) (the extra factor is due to h).
Combining everything finish the proof of running time of EXACTDS.
D.3.2 APPROXDS
In this section we present the data structure APPROXDS. Given BATCHSKETCH, a data structure maintaining a sketch the primal-dual pair (z, E neor x R APPROXDS maintains a sparsely- changing foo -approximation of (I, This data structure is slight variation of APPROXDS in Gu & Song (2022).
Algorithm 4 The APPROXDS data structure used in Algorithm 1.
1: data structure APPROXDS
Theorem D.6
2: private: members
3: Capx.a; fapx.s EK
4: EEN
5: BATCHSKETCH bs
This maintains sketch of H w. and H- s. See Algorithm 6, 7, 8. This a pointer the EXACTDS (Algorithm 2, maintain in parallel
6: EXACTDS* exact
APPROXDS.
7: T. IR & (T. 3) is a sparsely- approximation of (x, s). They have the same value as
(I. s). but for these local variables we use (z, 3) to avoid confusion.
8: end members
9: procedure INITIALIZE(,,
Ra
Ri h
Thoo Xm E2S E,, H1/27 H 1/2 8,
P N,,,,,,, Be, E ERR Ge, E R",,"" EXACTDS* exact, Capa,as Capa, Sape E (R))
10:
E (+ 2+0, 1
11:
Eapx,* Cupx,ar Saps.. Capx,a
12:
INIIIIIAL, h, Es, H H1/32, H=1/3, c. Bo, Be, Bo, 8apx/9)
Algorithm 6
13:
I T,
14:
exact + exact
15:
end procedure
procedure UPDATE( E noot Teot 8ú R Beer Ben, //1/221 1-122/ de, E R
17:
bS.UPDATE(&2, dn, 0R deg de, 0 u 1/2g1 6cn)
Algorithm 7
18:
é (+(+1
19:
end procedure
20: procedure MOVEANDQUERY(C.. Bo, Bo, E ER 3z, G. 1R))
21:
bs.MovE(BB,... 3r, Bc. 3r, 3e)
Algorithm 7. Do not update yet
22:
83 < QUERYE((APx..C log 1))
Algorithm 5
23:
8, QUERYS (Eapx. (2 log
Algorithm
F T+
25:
return (Ōz, 5z)
26:
end procedure
27:
end data structure
Theorem D.6. Given parameters Capx,r, Capx.s E (0, 1), Supx E (0, 1), So So E R. such that
1/2 x
222 20
(++1))
Ce, Cas H
ur 8
1/2
H H-122 1/2 8 (1+1)
||2 - Co
28

---

Published as a conference paper at ICLR 2025
Algorithm 5 APPROXDS Algorithm 4 continued.
data structure APPROXDS
Theorem D.6
2: private:
3: procedure QUERY E IR)
4: I
5: for jorj= 0 - llog2 do
6: if l mod 2 = 0 then
7:
I TU QUERYQ(E 2 1,e)
b Algorithm
8:
end if
9:
end for
10:
òz
11:
for
all i E I do
34 <- exact. QUERY
Algorithm
13:
if [ xz zi > then
14:
òz Fy 3z.
15: end if
end for
17:
return &
end procedure
19: procedure (UERYs(E IR)
20:
Same as QUERY., except for replacing x, with 8, s, and replacing "I[ă; x4 lalla," in Line 13 with "!
21: end procedure
22: end data structure
for all E 4 1,, data structure APPROXDS (Algorithm 4 and Algorithm 5) supports the following operations:
INITIALIZE(I,S E R h Ri h Blltot Xm C2, Esr H H. H T Bs.8,. E &, E R" . N, EXACTDS* exact, Capa, Capx, Sapx E R): Initialize the data structure in Ö(n nTm) time.
MOVEANDQUERY((... Be, Be, R. Be, R'"): Update values of Bz, Bo, Bee Ber 3 by calling BATCHSKETCH.MOVE. This effectively moves (z))) to (x(1+1) (el(t)),((tt1) s((+l) while keeping ( unchanged.
Then return two sets
(E) I1(0) T.2
C [m] where
_ 2 {i [n]
|
maree 2
1/2
H 1/2
w.2(0) 2
(041) 2 Eapx
20
- 2 fie
[n]
w.5
1/2 (c)
si
HH -HHH1 1/2 (1+1)
2(1)
S
1|2 2 capx.s)
satisfying
Σ
[
05859-1
Σ
[T(00
Ō(E 22.00
05259-1
For every query, with probability at least 1 Sapx /g. the return values are correct. Furthermore, total time cost over all queries is at most
δ ((E + + 55p 22)
UPDATE(Ō E T òn E Riat 67 Ruitt dey de,, 6 0/122 6 11 de, E IR
Updute sketches of 222 H2/2 ((+1) and H 1122 ((+1)) by calling BATCHŠKETCH. UPPAATE "
29

---

Published as a conference paper at ICLR 2025
This effectively moves z!0) to x((-1) while keeping (zltt1),s(tt)) s(++!)) unchanged. Then ad- vance timestamp
Each update costs
(1 +
lónllo + 11o + 116,, llo + llo) + 116 +1 #/////l + 116 ||12=// + 116c, lo)
time.
Proof. The proof is essentially the same as proof of (Gu & Song, 2022, Theorem 4.18). For the running time claims, we plug in Theorem D.8 when necessary.
D.3.3 BATCHSKETCH
In this section we present the data structure BATCHSKETCH, It maintains sketch of H x and HF s. It variation of BATCUSKETCII in Gu & Song (2022).
We recall the following definition from Gu & Song (2022).
Definition D.7 (Partition tree). A partition tree (S, x) of IR is a constant degree rooted tree = (V, E) and labeling of the vertices X + 2|m], such that
(root)
[n]:
if a leaf S, then Ix(v) 1;
for any non-leaf node E V, the set {x(c) is child of o} is partition of (vv).
Algorithm 6 The BATCHŠKETCH data structure used by Algorithm 4 and 5.
1:
data structure BATCHSKETCI
Theorem D.8
2:
members
3: Prxmno
All sketches need to share the same sketching matrix
4: S. partition tree
5: EEN
Current timestamp sketch cz, sketchWT c, Algorithm 10
6:
BALANCEDSKETCH " h, sketchW
7: VECTORSKETCH sketch /2, H1/2,, sketch 1/2g sketchc,
b Algorithm
8: Be,3e,& Sc, E R, Fe. 3, E R"
9: (history [t)
Snapshot of data at timestamp t. See Remark D.9.
10:
end members
11: procedure INITIALIZETT E Rio h E 1R h
Rh Xm Ex: Esr 1122 Hm11/2 3, Cs
1 Be. Be, So,, IR, Bz, B, E 12mh Öapx R) Construct partition tree (S, x) x) as in Definition D.11 + (log" (Meot) (1(0pp))))))
12:
13:
14:
Initialize & E 1Rex with iid (O,
B2 +- Bz, Bg +- BS. 9c, + Be Br &- Brg Bg x- 3.
16:
sketch INTTIALLIIZ 1.
Algorithm 10
X, 4,T,h)
Algorithm 10
18:
sketchW Eg. INITIALLIE (S, 4,T,6))
Algorithm 10
sketchWT IN.IIIIAIII.IIIIIIIIIIIIIIE 2,T,6,
Algorithm 10
20:
sketchH 1/2
w, LNITIALII
H. T
A Algorithm
21:
sketchH 1/2 INNIIIIIIII X, $. H-122/
A Algorithm
w.F
22:
INIEIAC,, x. +,C,)
Algorithm
23:
L 1+0. 0. Make snapshot history
b Remark D.9
24: end procedure
25: end data structure
Theorem D.8. Data structure BATCHSKETCH (Algorithm 6, 8) supports the following operations:
30

---

Published as a conference paper at ICLR 2025
Algorithm 7 BATCUSKETCCI Algorithm continued.
data structure BATCHSKETCH
Theorem D.8
2: procedure MOVE(3,, S,, Pe, E IR, 32. 3. E )
3: B2 Bz. B. + Ba, Sc. t Be,, Be Bus Bg Bg
Do not update & yet
4: end procedure
5: procedure UPDATE(Ō TRMoot Sn E Rooor .6 R de, de, 8 11/32 a. 6 H 1/2g dc,
6:
sketchW/ h. UPDATE(ÓZ 5n))
Algorithm 11
7:
sketchW UPDATE( %7)
Algorithm
8:
eketchW" UPDATE ÖZ de,)
Algorithm I1
9:
AketchW/T.... C..D:///T.UUUUTUTTET UPDATE Öz. 6e,)
Algorithm 11
10:
3 UPDATE. 1123
2.
1/2g.
Algorithm
11:
sketchH 1/2
MMetchH,"!/PUUPDAATEE S.UDATE (1111)))))
Algorithm
12:
sketche,.UPDATEE06.))
Algorithm
13:
++1+1
Make snapshot history[4]
Remark D.9
15:
end procedure
end data structure
INITIALIZE(T
Rm R h
IR Rio Xmh Ces Ess H1123, H_1ss,,,
c, R): Initialize the data structure in
IR Pinre Bar Be, E R. Bn R Saps Ön nTm) time.
MOVE( Bar Be, E ER, Bo ERT):: Update values of Bu, Ber Be, B2, in O(m) time. This effectively moves (xl0), ((0) to (x(e+1), (((tt), while keeping Z(() unchanged.
UPDATE( E Re 6b E Root .6 E R xm Se,, Se, 8 11 11 3. de, E R
Update sketches of 1/2 anee (f+1) and 1-1 while keeping (xlf+1), ((++)) unchanged. Then advance timestamp
w.ale) H s(et1), This effectively moves x to x(++1)
Each update costs
(( + 18allo + l6g 116e, llo + (10. + 116 11 116 112221l 116c, ll))
time.
QUERYa(!' E N, E R): Given timestamp E', return set S C [n]] where
$ $2 fi [n]
1/2
x
H/2
I
(E+1)
l| 2 c)
and
|S] 0 (l +1)
Σ
H
Σ
|it(t)
(1--1) (120)
E'SISE
Esist-1
where is the current timestamp.
For every query, with probability least most
8, the return values are correct, and costs
Ō (c- +1)
Σ
H |H
1/2 x (141)
|| (1+11) (12,0))
201
E'SiSt
running time.
QUERYs(E' E N,e E R): Given timestamp E, return set S C [n]] where
S2 {i
[n]
||
(1)) 8
1/2
H - 112 1/2 (8+
||2 2 c)
w.
31

---

Published as a conference paper at ICLR 2025
Algorithm 8 BATCHISKETCI Algorithm 6, 7 continued.
Algoritn
data structure BATCHSKETCH
Theorem D.8
2: private:
3: procedure QUERY,SKETCH(( S)
Return the value of & x((=) (H w. x(v)
4: return sketch 2
M.rcchWT.QUERRYY sketchW h. QUERY())
B2 + sketchW (MeechhW"T..QUERRR QUERY(
> Algorithm 9, 10
5: end procedure
6: procedure QUERYSSKETCH(: S)
> Return the value of 4b x()) H ( s)x(1)
1/2
7: return sketch S. (U) sketchc,.QUERE QUERY(U) Se, sketchW h. QUERY())-
B. + tssetchh "...RV/WW.QUUU((QO B. sketchW" QU.QUEY(())
& Algorithm 9, 10
8: end procedure
9: public:
procedure QUERY E E, R)
11:
Lo (root(S))
12:
S
13:
for d 0 oo do
14:
La = 0 then
15:
return S
16:
end
La+1 <- 0
18:
for v E La do
v is leaf node then
20:
S Su
else
22:
for u child of v do
23:
if [QUERYYSKETTCN()) - histor/eleuuuxxsseerrtollll > 0.9e then
24:
L Lat1 U [u}
25:
end if end for end end for
26:
27:
28:
29:
end for
end procedure
31:
procedure QUERYs(L' E N, R)
Same as QUERYr, except for replacing QUERY2SSKETCH in Line with QUERTSSKETCH.
33: end procedure
end structure
and
|S| (1-
Σ
||
1/2 50
-112/2 1/2 w.
HI
(t-1)
|z(0 ((+ (2,00
EESIEE
e'Sit-1
where is the current timestamp.
For every query, with probability least most Ō(72 (e-
8. the return values are correct, and costs
Σ 8'31SE
1/2 s()
H 1/
(+11)
Σ
|z(1) x (1+1) (l,,)))
210)
running time.
Proof. The proof is essentially the same as proof of (Gu & Song, 2022, Theorem 4.21). For the running time claims, we plug in Lemma D.10 and D.12 when necessary.
Remark D.9 (Snapshot). As in previous works, we use persistent data structures (e.g., Driscoll et al. (1989)) to keep snapshot of the data structure after every update. This allows support query
32

---

Published as a conference paper at ICLR 2025
historical data. This incurs Ollog Noot) (1) multiplicative factor in all running times, which we ignore in our analysis.
D.3.4 VECTORSKETCH
VECTORSKETCH a data structure used to maintain sketches sparsely-changing vectors. It is a direct application of segment trees. For completeness, include code (Algorithm 9) from (Gu & Song, 2022, Algorithm 10).
Algorithm (Gu & Song, 2022, Algorithm 10). Used in Algorithm 6, 7, 8.
1: data structure SKETORSKETCH
Di Lemma D.10
2: private: members
3: & E R
4: Partition (S,x)
5: E Rf
6: Segment tree T on [n]] with values in R"
7: end members
8: procedure INITIALIZE(S,) X partition tree, grxtt 2 E ghott)
9:
(S, + (8,x).
10:
I +- x
11:
Order leaves of S (variable blocks) such that every node x(e) corresponds to contiguous
interval C [x]]
12:
Build segment tree T on [xi] such that each segment tree interval / C [n]] maintains
1Ι E R².
end procedure
14:
procedure UPDATE(O_ E (
for all [N4or] such that 6z,4 + do
16:
Let E [n] be such that is in j-th block
17:
Update T at j-th coordinate Djæg + 4,0, + 0; 8z,5-
18:
Ti < + oz
19:
end for
20: end procedure
21: procedure QUERY(U E V(S)
22: Find interval I corresponding to (1)
23: return range sum of T on interval
24: end procedure
25: end data structure
Lemma D.10 (Guu & Song, 2022, Lemma 4.23)). Given a partition tree (S. x) of R." and
JL sketching matrix & Rr the data structure VECTORSKETCH (Algorithm 9) maintains "x(e) xx(o) for all nodes in the partition tree implicitly through the following functions:
INITIALIZE(S, X, 4): Initializes the data structure in O(rDeat) time.
UPDATE(Ó__ E I Maintains the data structure for & & I + dz in O(rlló logn) time.
QUERY((U E (S)): Outputs 2x(v)Tx(v) in O(r logn) time.
D.3.5 BALANCEDSKETCH
a data structure for maintaining sketch of vector of form Th,, h, where W = I.= H w.= and h E Rse is sparsely-changing vector.
this section, we present data structure BALANCEDSKETCIL.
This is variation of BLOCKBALANCEDSKFTOO
in Gu & Song (2022).
We the following construction of partition tree.
Definition D.11 (Construction of Partition Tree). We fix ordering of [m] using the heavy-light decomposition (Lemma A. 10). Let S complete binary tree with leaf set [n] and ordering л. Let x map node to the set of leaves in its subtree. Then (S, x) is valid partition tree.
33

---

Published as a conference paper at ICLR 2025
Algorithm 10 The BALANCEDSKETCII data structure is used in Algorithm 6, 7,
data structure BALANCEDSKETCH
Lemma D.12
2:
private: members
3: DER*XO Rex
4: Partition tree (S,x) with balanced binary tree B
5: I JEN
6: HE Ruhttt, E E Rhoor, Hw E RRot Riot
7:
{L14 E T Xfroot )120
8:
(Jo E R FUES
9:
{Zg E R JUEB
10: {UT TR")
11: {tr E NNEEB
end members
13: procedure INITIALIZE(S, x: partition tree, E Rexmtt , Rutt h E Ru Xk)
(S, + (S,x). (
15:
+ t
H + Vzdz)) B3 + QH+HH n.5
17:
Compute lower Cholesky factor [z((] of Bz do
for all
J f- 4 x(v) 1/2
19:
20:
end for
21:
for all
E do
22:
Zp f- J -T]
23:
yu + Zo (I- LA(w))h
24:
to +- t
25:
end for
26:
end procedure
27:
procedure QUERY(U E S)
28:
if S then
29:
return J. L 11|1-T/
end if
31:
ALT + (Lalt] - Lx(tu]) IA(e)
32:
dz. t- -[=]((]- Aiz Z.
33:
Z. - Z, + 82.
34:
& 8z. (I- (I-Ix hh
35:
D + 30 + +dgg
36:
to t
37:
De <- Zu IN( h
38:
return vf + +
39:
end procedure
40:
end data structure
Lemma D.12. Given an elimination T with height n, a matrix & RRrxn and partition tree (S, x) constructed in Definition D.11 with height Ó(1), the data structure BALANCEDSKETCH (Algorithm 10, 11, 12), maintains & xx() (WTh)x(o operations
h)x(o) for each E V(S) through the following
INITIALIZE((S,)) x) parriiiio tree, & g Teae data structure in Ö(r(n7w-1 +nrk)) time.
R(moox)):: Initializes the
UPDATE(55 E Mo 5n E IR.
Pior xk); Updates all sketches in implicitly to reflect (W, h)
updating to (Wpuew, hewww ) in (vr72k) time.
QUERY(( E S): Outputs & x(() (WT h)x(m) in Ō(r72k) time.
Proof. The proof is almost same as the proof of (Gu & Song, 2022, Lemma 4.24). (In fact, our W simpler than the used in Gu & Song (2022).)
34

---

Published as a conference paper at ICLR 2025
Algorithm 11 BALANCEDSKETOH Algorithm 10 continued. This is used in Algorithm 6, 7, 8.
data structure BALANCEDSKETCH
2:
procedure UPDATE( E R 5b E Tio Xk)
3: for [n] where ò + do
4: UPDATET((2.))
5: end for
6: for all dn.; + do
7:
< (1)
8:
for all PB (0) do
9:
y - 4 + Zu [() 8n.
10: end for
11: end for
h + h 6n
13:
end procedure procedure UPDATET(OZ, (R.))
15:
++ I 0z,
(H.,z.(1.)) (T,) HW, (k.,))
18:
Compute A fx such that [x +- Lgst 1] Ar br is the lower Cholesky factor of A w.
A A
19:
S + PB (A (1)))
(1)))
20:
UPDATEL ALT)
21: UPDATE. (t,)..
22: end procedure
23: end data structure
For INITIALIZE running time, we note that computing Z, for all e B takes Ö(rn7"-1) time by (Gu & Song, 2022, Lemma 8.3). Because Zu is supported on the path from v to the root in 'T, we know that nnz Z) = (rnT)). Therefore computing v for all E B takes (rnrk)) time.
Remaining claims follow from combining proof of (Gu & Song. 2022, Lemma 4.24) and (Gu Song, 2022, Lemma 8.3).
D.4 ANALYSIS OF CENTRALPATIMMANNTEAANCEE
Lemma D.13 (Correctness of CENTRALLATHHAANEENNNCE,,, Algorithm implicitly maintains the primal-ducl dual solution pair (z, s) via representation Eq. (8)().. It also explicitly maintains (τ,z) RRono x R teot such that [ Iz Lilla, < and || - sillz, S tew; for all i E [m] with probability at least 0.9.
Proof. We correctly maintain the implicit representation because of correctness of exact. UPDATE (Theorem D.3).
We show that | x ||,, < E and |ل - sill < LEw; for all i E [m]] (c.f. Algorithm 20, Line 16).
approx maintains an foe approximation of H 1.2 For E < q, we have
1113 x (0) ||2 1|6 S a Sx
where the first step from definition of || I|,, the second step follows from Lemma F.11, the third step follows from definition of Su-
By Theorem D.6, with probability at least - Sapx. approx correctly maintains v such that || H 11/22 5 capx,a E. Then
21 S w. 1/2 1/11/1111//2
11/2201000
H
S <w;" -1/2 <
35

---

Published as a conference paper at ICLR 2025
Algorithm 12 BALANCEDSKETCII Algorithm 10, 11 continued. This is used in Algorithm 6, 7, 8.
data structure BALANCEDSKETCH
Lemma D.12
2: private:
3: procedure UPDATEL(S B, ALT Bioto XRoot)
4: for all E S do
5:
8z. s - (La - (L= -1 --s((e)) IN(O) ZT)
6:
8z t- - [I [1]1]-1 L (Zx + )
7:
Z. + Zu + 8z. + +8%
8:
bg + (8z. + 8'2,) (I I((D)
9:
V. + yd + d5g
10: to 6 t
11: end for
12: end procedure
13: private:
14: procedure UPDATEH( E [n], A (((i.) <R.xx)))
15: Find such that x(1)= = {i}
16: A (1.))) + (H.,z,(1.) + An (t.... 11
(3..)
17:
5J.-
f 4,- - (22
18:
for all E (u) do
Jn f- Ju + 6J.
20:
if E B then
21:
8z. + ò5, Ltlto] -T
22:
Z. Z. + dz.
23:
5g + 8z. (I Ja(D)) h
Yo y døg
25:
end if
26:
end for
27:
Hw, H A A (w....
28:
end procedure
end data structure
Note that the last step is loose by a factor w 1/2 When w;s are large, we could improve running time by using tighter choice of capx,z. did in Gu & Song (2022). Here we use loose bound for simplicity presentation. Same remark applies to s.
The proof for is similar. We have
17
HHH1
1/2 W.F.
5,112
5g
11
S a S
8
and
3 - - S
1 11112 8l1oc S w w/// aaxx, ST. wy-
1/2
u.z
Lemma 1.14. We bound the running time of CENTAILPATIMAITTENNCCC as following.
TENTRALPTTIMANNTENNNN..NTITT takes Ö(1 nTm + nm==1)) time.
If total running time
is called N times, then it has
Ō((Nn-1/2 log(tmax/ttmin) n + 1 + 1m1-111/2)
CENTRALLPAHHMMINTEEANNEEE.O takes Ö(TT)) time.
Proof. INITIALIZE part: By Theorem D.3 and D.6.
OUTPUT part: By Theorem D.3.
36

---

Published as a conference paper at ICLR 2025
MULTIPLYANDMONE part: Between two restarts, the total size of Ly returned by approx. QUERY bounded by 2 by Theorem D.6. By plugging in Çe 2a, Capx,a E, have Eiety LE Ö(q²). Similarly, for s have Eie(s) - Ögg²).
Update time: By Theorem D.3 and D.6, in a sequence of q updates, total cost for update is ( + rm?)). So the amortized update cost per iteration is Ö(q(7 + rm?)). The total update cost is
number iterations time per iteration O(N((T(Ng))) rm²)).
Init/restart time:
We restart the data structure whenever & > q or | t > te,, so there are O(N/q + (og((max//min)ds restarts in total. By Theorem D.3 and D.6, time cost per restart O(n(7"11 + tm + m)) So the total initialization time is
number of restarts time per restart Ō((N/q log(fmax/tminn n( + rm + m²-
Combine everything: Overall running time is Ö(Ng(r2m + rm² + (N/q + mes(tmax/tminler -1( Tm + m m
Taking & c, the optimal choice for is
n + Tm?) -1/ 1/2 1/2 (7W-1 Tm + m 1w-131/2
achieving overall running time
Ö((Nn 1/2 - Proof of Theorem D.2. Combining Lemma D.13 and D.14.
)) nn(r2 + Tm² rm?) 1/2 + Tm + m 1112)).
D.5 PROOF OF MAIN STATEMENT
Proof of Theorem D.I. Use CENTRALPATHMAINTENAAC (Algorithm 1) as the maintenance data
structure in Algorithm 20. Combining Theorem D.2 and Theorem F.1 finishes the proof.
E ALGORITHM FOR LOW-RANK QP
In this section we present 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 E. we present the main statement of Section E.
In Section E.2, we present the main data structure MANTRALPATHMAINTEEMT
In Section E.3, we present several data structures used in CENTRALLPAHHMAINTEEAAA including EXACTDS (Section E.3.1), APPROXDS (Section E.3.2), BATCHSKETCH (Sec- tion E.3.3).
In Section E.4, we prove correctness and running time of LENTRALPPTHMAIITTEERR data structure.
In Section E.5, we prove the main result (Theorem F.
E.1 MAIN STATEMENT
We consider programs of the form (16), i.e.,
1 T Qx c
min BER
37

---

Published as a conference paper at ICLR 2025
s.t. Ar b Ij E
Vi [x]
where Q E Sthot C Rsoo, A E R. b E R", Kr C I is convex set. For simplicity, assume that n, = 0(1) for all [n].
Theorem E.1. Consider the convex program (16). Let à; K2 - IR be u.-seff-concordamt barrier for all E [n]. Suppose the program satisfies the following properties:
Inner radius T: There exists E mR such that Az b and B(z,r) K.
Outer radius R: K. _ B00, R) where 0 E Rus.
Lipschitz constant L: [[Ql2=>2 L, ||d|2 S L.
Low rank: We are given factorization UV where U.V V E R
Let (wi))e(n) E RR and n. = Eie(n) Wi4y. Given any 0 < 6 I 1 2. we can find an approximate
solution & E K. satisfiying
Qz tc min b.zkk
{
1 === +
2 x
|Ax - -111 s (c||R|| + l[b|[1),
in expected time
Ö((Vkn 11/2 log(R/(re)) n(k m)w++))/2 When max 14 O()), wi the running time simplifies to
Ö(n(k + +m)(+1)/22) (og(R((re)))
E.2 ALGORITHM STRUCTURE AND CENTRAL PATH MAINTENANCE
Similar to the low-treewidth case, our algorithm is based on the robust IPM. Details of the robust IPM will be given Section F, During the algorithm, we maintain a primal-dual solution pair (2r, Rrot R Root on the robust central path. In addition, maintain pparsely-changee approximation (z, E Reot x Reot to (x, s). In each iteration, we implicitly perform update
I I+ x x+
iB
1/2 Pu... P....,)
BB B 1/22 W...7
- s 16 - H
1/2 8n B B-LBLLB 1/211 P F. w.r,1
w
w.z
where
H = (iow((
(see Eq. (24))
B B = + iH
(see Eq. (25))
P Pu.z,i B-1/2 (AB-1 1B221/2 - W.F.F
(see Eq. (26))
",I,
and explicitly maintain (z, 3) such that they remain close to (x, 8) in -distanacc
This task is handled by the SETTRAA.PATMMAINTENNN structure. The robust IPM algorithm (Algorithm 19, 20) directly calls it in every iteration.
data structure, which main data
The CENTRAL SAAHPAAAAAIINNNNNCCCC data structure (Algorithm 13) has two main sub data structures, EXACTDS (Algorithm 14, 15) and APPROXDS (Algorithm 16). EXACTDS is used to maintain (I, s). and APPROXDS used maintain (I, 3).
(Algorithm 13) implicitly maintains the central path primal-dual solution pair (x. 8) E R x R and explicitly maintains its approximation (I. 3) E T yoor x TR Teoe using the following functions:
Theorem E.2. Data structure
INITIALIZE(Z E Reat E R to E RR.. E (0, 1)): Initializes the data structure with initial primal-dual solution pair (2, 8) E To x IR Tiss, initial central path timestamp to Rh in ( m=-1)) time.
38

---

Published as a conference paper at ICLR 2025
Algorithm 13 Main algorithm for low-rank QP.
data structure CENTRALPATHMAANTEAAAGC
Theorem E.2
2:
private: : members
3: EXACTDS exact
Algorithm 14, 15 b Algorithm 16
4: APPROXDS approx
5: EEN
6: end members
7: procedure INITIALIZE(Z, 8 E RRot E & (0,1))
8:
Exact.INITIALIZEI 3, 8, ++00
Algorithm 14
9:
10:
maxi N Vk log log nsl Pr q 5 (k2 + (d 14-1)112 6apx,2 +- 5x 20x, dapx
12:
13:
Eapx,s E So 3at
14:
approx. ApprOX.NIIIIALIZEEE,.... h, ".,.H', 1222 w,z H 3. 3o Bor Brr B B. q. &cxxact,
Eapx,a Sapx,s: 0apx)
> Algorithm 16. Parameeess from æ to 3, come from exact. &rexact is pointer exact
16:
end procedure procedure MULTIPLYANDMOVE( e IR+))
18:
e+ ++1
19:
if [ - 1] > 1. e, or & then
20:
2, 4- exact.OuTPUTI INITIALIZE(:, 8, 1,
Algorithm 15
21:
22:
end if
23:
Bz, Bo, 3e. 3., B. + MOVE()
> Algorithm 14
24:
Sz, t approx. aepraa. WWW.MMOVVANDQERY(C.. Bo, Bz. 3zr Bz. B.)
> Algorithm 16
25:
Ōn. oz. 8H!22 8 1H +- exact. UPDATE Oay 6z)
H
Algorithm 15
26:
approx. approx.UPDATEEN UPDATE( 5n ox d7 62//2 (H1112)
Algorithm 16
27:
end procedure
28: procedure OUTPUT()
29: return exact. OUTPUO)
P. Algorithm 15
30: end procedure
31: end data structure
MMLLIPLYANDMOVEE Ryo): It implicitly maintains
8, 1, I x -TB -12 Pw.,,.)) B-1122 w
1/2
3.i
1/2 8 s E6 12 Hw
u.. 1/2 P
W.J.F
where Ha,z. s. B P P.... are defined in Eq. (24)(25)(26) respectively, and some timestamp satisfying [ t| S E, T.
It also explicitly maintains (I, 8) E Toot X such that -14-=zll7, S and |3 -i Lew; for all i (n) with probability at least 0.9.
Assuming the function called most N times and decreases from mmax to bmin the total ranning time is
Ō((Nn-112 112 llog(mmxx//min)) n(k(w+1)/2 + m (w+1/)2)
OUTPUT: Computes (x, 8) E R Bror x Bittot exactly and outputs them in ( m)) time.
39

---

Published as a conference paper at ICLR 2025
E.3 DATA STRUCTURES USED IN CENTRALPATHMAINEENNNE
In this section we present several data structures used in CENTRALPATHMAINTENACCC including:
EXACTDS (Section E.3.1): This data structure maintains an implicit representation of the primal-dua solution pair (z, s). This is directly used by CENTRAL,AAHHMAITENAA
APPRoXDS (Section E.3.2): This data structure explicitly maintains an approximation (z, 8) of (z, s). This data structure is directly used by CENTRALPATHHAAINEE.NNCC
BATCHSKETCH (Section E.3.3): This data structure maintains sketch of (I, s). This data structure is used by APPROXDS.
E.3.1 EXACTDS
In this section we present the data structure EXACTDS. It maintains implicit representation of the primal-dual solution pair (z, s) by maintaining several
vectors (see Eq. (12)(13)). Theorem E.3. Data structure EXACTDS (Algorithm 14, 15) implicitly maintains the primal-dual pair (x, 8) E Boot x Rroe, computable via the expression
2 H 1+H11/2
hB 1 + 211hhh, + H H 1/2 h3,
1/2
(12)
6.
3
11 w
+ 11/23 + 11/22
(13)
w.
where 2,5 Roos h 3z 3a R. Bu, E Rt Be, 3. E R.
2
8
R Miss, h 1 T e Re xk h H 1-112
AT
E R Poos Xm
2.
w3
W.
The data structure supports the following functions:
INITIALIZE(1, 8. x. E BAo , E Rso): Initializes the data structure in Ómn(k" + +m))) time, with initial value of the primal-dual pair (x. 8), initial approximation (x, 8), and initial approximate timestamp t.
MOVE(): Performs robust central path step
I - -B -B AT (AB,) A 8,s,
(14)
$ st s + i5 BB + + A - 8
(15)
in O(K" + m²) time by updating its implicit representation.
UPDATE(Ō 8z :: Updates the approximation pair (I,3) to (w E + 85 TRae gnow + 8g E R. in O(k2 + m²) + [|x&(1o) time, and output the 1/2
changes in variables h, h, h, H 1 H w.F
Furthermore, h, H H. H= 1/2 changes in O(1 + (dsllo) coordinates, changes in O1:(1 lo + ((51l)o) coordinates, h changes in (m + [|اIllo)) coordinates.
OUTPUT(): Output x and s in Ö + m)) time.
QUERYa(i e [n]): Output x; in Ö(k m) time, This function used by APPROXDS.
QUERYs(i [n]): Output 8; in Ō(k m) time. This function is used by APPROXDS.
Proof of Theorem E.3. By combining Lemma E.4 and E.5.
Lemma E.4. EXACTDS correctly maintains an implicit representation of (x, s), i.e. invariant
H 1/11/2 1/2,,,8 HH-1/22/ + -H 1/TT
$ H1/2 h/2 2/2 1/2
h H-1122 8p E R Ptot hn
H
H-1121
Thot Xd To H-112 AT R Thot Xm
40

---

Published as a conference paper at ICLR 2025
Algorithm 14 This is used in Algorithm 13.
data structure EXACTDS
Theorem E.3
2: members
3: E, E R EE R+, Hu, E R Rtot Kfoor
4:
2. Ru h E R xk h E Xm
Bz. B. E R. Be. E. Bu, B, R
5:
W1: U2 RRxx, Ug 1mxx u4 Rm Ug RR ug 1Rick
6: a E R. 8p E ER"
7: KEN
8: end members
9: procedure INITIALIZE(E, 8, Root I Re)
10:
I+ T. 1
11:
8r f- 0, Bg < 0, 3. - 0, + 0, Ar t- 0, + 0
12: Hw. - V 2. (T)
13:
INITIALIZE/<Z, 3, Hw,7)
14: end procedure
15: procedure INITIALIZEN(E, & Rho Hw,x E RR X
16:
for E [n] do
(op)s f- asinn( my (((,,,
17:
(I,3,,)
+ ur +u, cosh 2 8,
18:
19:
end for
H 1/27 Nn t- W.z
H-112U h + H 1/2
H-1122 AT w.z
20:
W. 21: u1 UH A U2 VH AT, U33
14 - AH Jus us &- V H-1 W.
AH 7. AT
22:
6ns 46 x V UT
w.r
23: end procedure
24: procedure MOVE(
25:
Up
+ t E Rext
26:
U1
U3 1 11 T
110
12
R mx
27:
U2 t t 14 - -2,, "u 1
Ug E R
28:
B + 8r + -1/2
29:
3 3, (a) 1/2
E
U5
(1)-1/2.5
111.t-110 -t' U1 U2 02
3 3 (x)
1/2 01 1/
31:
8, 3,
32:
8, 3, + (
1/2
1/
U5 (a)
1/2
10
12 02
33:
B, 3, + 10 1/2 10
12
34:
return 8., 3r, Br. B,, Brn B
35: end procedure
36: end data structure
W1 UHELA Rdxm : =VH adm Ug AH-LA E mm
44 = AH E R" Us = =VH E Rd, up = VH Rdxd
w cosh? A
(x,,,,7))
Wh
iE[m]
8, 1/2 55 (T,
always holds after every external call, and return values of the queries are correct.
Proof. INITIALIZE: By checking the definitions see that all invariants are satisfied after INITIAL- IZE.
MOVE:
By the invariants, have
Un 1 + +
41

---

Published as a conference paper at ICLR 2025
Algorithm 15 Algorithm 14 continued.
data structure EXACTDS
Theorem E.3
2: procedure OUTPUT()
3: return 2 + H 1/2 h33 + H 13 H= $ H + H + H Hh
4: end procedure
5: procedure QUERY [n])
6: return T. +
H
11-11211 11,,82 Br + H 1121 1/2, hr + Hw 4-122 T,,..33
10,2
7: end procedure
8: procedure QUERY E [n])
9: return + H 1 3g + H 1 Bg + H H
10: end procedure
procedure UPDATE(OZ, dz E R
12:
AH.... <- 62 (I: 0z) Hw,
& A Aaa is non-zero only for diagonal blocks (7, i) for
which Sz,i #
13:
8 + EiE [n] δari + or 82, 0}
14:
6,
+0
15:
for
E 5 do
16:
Let 7 YA (T, 1), 71 =/(T 8z 3 8z. 1),
ps new
14(T ōg, 8g,
wi cosh 2 ( 7;) w. cosh?
+
mz
18:
3 a sinh 2 7/ now) 7: 14 new 8 p.
end for
20:
Ōn +
H 1/2 (5 - - + H
w. 83
og A H wiR
1/2UT
22:
6 - A W 1/2
23:
8 f (5n,3, B. + 6r (1
24:
0g r (6 6 3, + 0 Bz)
25:
+ + 8n, h t h + 6q. t h + 6
2 8 8
26:
41 + + U A H A
27:
+ VA H AT
28:
u3 u3 + AA AT
29:
44 - 44 + A(A H 1 (5 + 55, + H 1
30:
45 t 4s + V(A H W.
(3 + 83,) ) + H.I
31:
26 - u6 + V 11 T
32:
- & 8
83
33:
HEr + H + A...
return 8n, 87 57 n///2g 8 1/23
35: end procedure
36: end data structure
01 = 7 --AHH AT -T ( -V
10.2
= ABB AT
0/ 7-
-
-1 3
w.
= =AB;1_ 2,
By implicit representation (12),
8 -H HH11/2/0 h5 + H 112
11ō + H H= 1210
-H 8p (a) -a)-1/2 1/2
H- (ax) 1/22
1/ 1 (-us us u2 1 U2)
42

---

Published as a conference paper at ICLR 2025
A (a) (1)-11201]) 02
+ 1. ( w0.5
( 10.2
(AB='AT) 1 (n))
-H A ATAAATAAA (AB= A -
= 8. (T ( H - H
(I+T VH U 66
-i(E - w. - 2 w.
U
(I
1
V (AB AB 6p
VH
=i iR= A (ABE
A 5p.
Comparing with the robust central path step (14), we see that x updated correctly.
For from implicit representation 13 we have
6, = H1/2h0g + H h..
w.. 8. + H 1/12//5 hd
= --U ( (aa) (-us + u2 U2) + A (a) 1/2 Fej
15p - + A A")
Comparing with robust central path step (15), we see that s is updated correctly.
UPDATE: We would like to prove that UPDATE correctly updates the values &, h, h h, U1 U2, ug, un us, u6s a, ons while preserving the values of (z, In fact, by checking the defi- nitions, it is easy to see that h h, h, u1, U2, u3, U4, Us, ug, 6 are updated correctly. Furthermore
8,- 8z + 8r + 8q B, + 8 =0, 0, B,
8, 53 + 6, + 8 + 8 =
So values of (x. s) are preserved.
Lemma E.S. We bound the running time of EXACTDS as following.
(i) SRACTDDS.NNTIAAE (Algorithm 14) runs in Ö(n(k"-1 + m)-1)) time.
(i)) EXACTDS. MOVE (Algorithm 14) runs in Ö(K" + +m²) time.
(ii) EXACTDS. (Algorithm 15) runs in Ö(n(k + m)) time and correctly outputs (x, s).
(iv) EXACTDS. EXACTDS.QUERY& returns the correct answer.
QUERY= and EXACTDS, QUERRY (Algorithm 15) runs in Ö(k + m) time and
(v) EXACTDS. EXACTDS.UPDATT UPDATE (Algorithm 15) runs in Ō((k + m + (|5zllo) time. Fur- thermore, ||0,|lo. &a o; [ lo O(1]&zloo + [|бæ]lo), mnz(h) Odd(1)=1lo + [[0allo)).
mnz(h) (mm[[1621 + (lo.1o)))
Proof.
EXACTDS..NNTIALIE Computing un and U2 takes Tmatt (k, n, m) = Ō(((
m=-1) time. Computing U3 takes Tmai (m,, m) = 0(um*11) time. Computing takes ((mm) time. Computing u5 takes ((n)) time. Computing ue takes Tinat (ki, n, n,,k).
(nk time. All other computations are cheaper.
(ii) EXACTDS. Computing -1 takes Ö(K*)) time. Computing Vi takes ((m)) time. All other computations cheaper.
(iii) EXACTDS.OUTPUT Takes Ö(n(k m)) time.
(iv) EXACTDS. QUERYY: and EXACTDS.QUERYS. Takes Ō(k m) time.
43

---

Published as a conference paper at ICLR 2025
(v) EXACTDS. UPDATE: For simplicity, write t = ||ōzllo + ||6z|lo- Computing on takes Ö(t) time. Computing o takes O(tk) time. Computing 6x takes ((tm) time. Computing óą and 03 takes ([((k + m)) time. The sparsity statements follow directly. Computing in and u2 takes (tkmm) time. Computing uy takes Ö(tm²) time. Computing u4 takes ((m)) time. Computing us takes (tk) time. Computing us takes (th22) time.
E.3.2 APPROXDS
In this section we present the data structure APPROXDS, Given BATCHSKETCH, a data structure maintaining sketch of the primal-dual pair (x, 8) E R x Rioe, APPROXDS maintains cparsely-changing & approximation of (z.
Algorithm 16 This is used in Algorithm 13.
1: data structure APPROXDS
Theorem
2: private: : members
3: Eapx,at; Eague, R
4: LEN
5: BATCHSKETCH bs
1/22 See Algorithm 18. 6: EXACTDS" exact This is a pointer to the EXACTDS (Algorithm 14, 15) we maintain in parallel to
This maintains sketch of H and H
APPROXDS.
7: H. Te
(3. 3) is sparsely-changing approximation of (z, They have the same value (), but for these local variables we (3. 8) avoid confusion.
8: end members
9: procedure INITIALIZE(1, E Ru h E R h TR xk E R. Xm H 1/2g, H H= 1/2
Re Bz, 3, ER, RR 8, RR, ER" , EXACTDS exact, Capa,ay Capx.sy Sape R.)
L e 2+0, 1 Eapx,= f Eapx,a, Capx f Eapx,
11:
IN.TIIIIALIEEEE. R H 1/2 1=113,, Br, Bn, Bas Ba, Br, Ba, 8mpx/9)
Algorithm 17
13:
J 2,
14: exact
end procedure
16:
procedure UPDATE(OS E R .ō "eot 8g E R R x 88 E R xm
41/2
1/2g E &
u.
17:
bs.UPDATE(SR, 6n, og 5g ///22'0 n-y)))
Algorithm 17
18:
1+(+1
end procedure
20:
procedure MOVEANDQUERY(S3 3, R, 3,,3, RR, 1,, B. Rm)
21:
bs. Ds.MovE(Bz, 3o B. B. 3z 3.)
b Algorithm 17. Do not update yet
22:
QUERY (Fapx.. log 1))
Algorithm 16
23:
5z t (EE
log 1))
b Algorithm 16
24:
I I+ 8z. & 8 83
25:
return (6z, 8z)
26: end procedure
27: procedure QUERYZ(E R)
28: Same as Algorithm 5, QUERYZ.
29: end procedure
30: procedure (UERYs(E R
31: Same Algorithm 5, QUERYS.
32: end procedure
33: end data structure
Theorem E.6. Given parameters Capx, Capor, E (0, 1), Sapx E (0. 1), Sz, E R such that
1/2 w. (1)
1/2 w 2 ((+1) S | 1/2 1/2 w
111 s (1+1) ||2 S.
for all & E {o...., 1}, data structure APPROXDS (Algorithm 16) supports the following operations:
INITIALIZE
E Root h
E
Phot h
E
Rea xk T E
B Reus 11/2 H 1/2 1-12 E Me Br,
E
R, 3
3
E
R
Rmm
.
44

---

Published as a conference paper at ICLR 2025
N, EXACTDS* exact, Capx, Capx,s: Oapx R): Initialize the data structure in O(((k+m))
time.
MOVEANDQUERY($2..) 3. E R2, IR 32,3,.33 3, Bu, B2, Br, Br, Br, Bg by calling BATCHSKETCH.MOL This effectively moves (xl(), s(0) to (xle+1), s((+1)) while keeping I(() unchanged.
RT"): Update values of
Then return two sets L6 L.a [n] where
- 2 (i E [n] |1 1/2
m 1 - 11122 1et11112 2 Eapx, )
11,310 2i
4 2 1i E |n]
1|
w Si
1/2 (8) 11 1/2
w,2(1) 51
(e+2)
1|2 Eapx,
w
satisfying
Σ
1 = Ö Öc S29")
0S45q-1
Σ
[
05PE9-1
For every query, with probability at least 1 Sapx/0 / the return values are correct.
Furthermore, total time cost over all queries is at most
Ô ((E + <x (k: + m))
UPDATE(05 "oos 6n E RT 0g E R 0g E R 1///2 8 H- R
Update sketches of H1/2 x(1 and H 1/2 s(++1) by calling BATCHSKETCH. UPDATE.
w. w
This effectively moves x(f) to x(e-1) while keeping (x(t+1), s(t+!)) unchanged. Then ad- vance timestamp
Each update costs
Ó(118aloo nnz(5;) + mnz(6r)) +
1//3 + w.z
||H-1/28
time.
Proof. The proof is essentially the same as proof of (Gu & Song, 2022, Theorem 4.18). For the running time claims, plug in Theorem E.7 when necessary.
E.3.3 BATCHSKETCH
In this section we present the data strueture BATCHSKETCH. It maintains sketch of H x and 1 1/2 s. a variation of BATCUSKETCII in Gu Song (2022). Theorem E.7. Data structure BATCHSKETCH (Algorithm 17, 18) supports the following operations:
. INITIALIZE(F E Riot h Biot h E B xk h
xm -H H 1/2 E 8apx E R): Initialize the data structure E TRm R
TR Be, Gg E RR Ba, 3. E RR Re, B, in Ö(n(k + m)) time.
MOVE(32.. 3o E R. Br, E RR, R B. E R" Update values of B2, 3z. B., 32.32..2 Be in Ol: + m) time. This effectively moves (z((), s(0) to (x(e+1), g((+1) while keeping (1)) unchanged.
UPDATE(52 R.. Peos, .6p E 1Rthot,8p E X 5g E Re xm 8 51//22 8 1/22 (e
Update sketches of 1/2 (4+11 and H 1-112 ((++)) This effectively moves (1) to (f+1)
w. o while keeping (z((+1) s(++1)) unchanged. Then advance timestamp Each update costs
w, o
Ō(16allo + nnz(fr) mnz(dr) + 12/27 ||T
45

---

Published as a conference paper at ICLR 2025
Algorithm 17 This is used by Algorithm 16.
data structure BATCHSKETCH
Theorem E.7
2:
members
3: Rexaot
All sketches need to share the same sketching matrix
4: S.X partition tree
5: E N
Current timestamp
6:
VECTORSKETCH sketch 1/2 sketch/f "
u,, sketchh, sketchh, sketchh > Algorithm
7: 8r.,3, 33 E ER., Ba, , B. E Rd R 33 E Rm
8: (history[t] [t])220
Snapshot of data at timestamp t. See Remark D.9.
9: end members
procedure INITIALIZETT R Ru T Ru xk h R Xm 12 1/12
1 B2, R. 32, 3, Rd Be, 3g Rm Supx E R)
11:
Construct any valid partition tree (S.x)
12:
<- (log' (Manc) log(1/ (1/6pp))
13:
Initialize + E TRee with iid (O,
B₂ B2, Bg Br, Br. 3, + 33, Bg t B.. 3.
15:
sketch 2/2 INITIALIZE X, 1 1/22g)
b Algorithm
w.z -1 sketchH, wiF INITIALIZE(S, X, ④ H=)3
16:
D Algorithm
17:
sketchi.. x, 1,
D. Algorithm
18:
sketch/.. D,, 4,1)
Algorithm 9. Here we construct one sketch for h. for
every E [k]]
19:
INeTchh.IIITIALIIEE(,,, h,) every E [m].
Algorithm 9. Here we construct one sketch for h. for
20:
21: Make snapshot history [4]
b Remark D.9
22: end procedure
23: procedure MOVE(3,, 8, R. Bz, E 3r, R)
B₂ Bz, Bg
B., 9r + Br, 3r + Ba, Be t Br.Ba + B.
Do not update & yet
25: end procedure
procedure UPDATE($ E TR Sb. RR 8g E R. xk 8g Br xm 8
#1/2
8
H 1/2g w
27:
sketch//
1/2g. UPDATE(6 13/30) 11.7
Algorithm
28:
- SUPDATEE (0
sketch w.z
A Algorithm
29:
Sketchl..UPDATTE(M))
D. Algorithm
30:
sketchh..PPAATE (dz)
D Algorithm
31:
sketchh.UPDATT((
A Algorithm
32:
{t ++1++
33:
Make snapshot history[[]
b Remark D.9
34:
end procedure
35:
end data structure
QUERYa((' N E R); Given timestamp E', return a set S C [n] where
S fie
1/2
(
H 1/2 (f+1)
(0+1)
1|2 2 2c.
and
|S| 0 2 - +1) Σ
1/2
112
(++1)
Σ
(i(t)
=(c-1)||20)
PSiSt-1
where & is the current timestamp.
46

---

Published as a conference paper at ICLR 2025
Algorithm 18 BATCUSKETCHI Algorithm 17 continued. This is used by Algorithm 16.
data structure BATCHSKETCH
> Theorem E.7
2: private:
3: procedure QUERY,SKETCH(C S)
Return the value of @ xx(v) (H 1/2 x(v)
4: return
sketchh..UERY/CE sketchh. (UERY())))) 3,
sketchh. (UERCh).Q
Algorithm
5: end procedure
6: procedure QUERYSSKETCH(: E S)
> Return the value of t x(0) H ( x)x(1)
7: return 2 S.QUER((U) + sketchh..UERR((t) B, + sketchh Qurrchh),QuE(())) 8
1/2
sketchh. (u)
> Algorithm
8: end procedure
9: public:
procedure QUERY E N,e E IR)
11: Same as Algorithm 7, QUERYZ, using QUERYZSKETCH defined here instead of the one in
Algorithm 7.
12: end procedure
13:
procedure QUERY E N,,c E R)
14:
Same Algorithm 7, QUERYS, using QKETTH defined here instead of the in
Algorithm 7.
end procedure
16: end structure
For every query, with probability at least 6, the return values are correct, and costs at most
Ō((k + mm) 3 - 1) Σ
1/2 1/2 H
F0 20)
20
(1+1
Σ
|
x((+1) (12,))
F(o)
rumning time.
QUERYS((' E N,e E R): Given timestamp E, return set S C [m]] where
$2 [n)
||-1// (8))
w.x(e") 8 ex $
1-122 1/2 (f+1)
||2 2 c)
and
JS| -0 +1)
Σ
11-1/2
H- (1-1)
Σ
|
(++)) 2200)
w
B'SESE
e'SIt-1
where is the current timestamp.
For every query, with probability least most
8. the return values are correct, and costs at
Õ(kk + m) (e- +1)
Σ
1/2
H (011)
Σ
||x(0)
x x(++1)
(1.00)
201
P'ESS
PS13e-1
running time.
Proof. The proof is essentially the same proof of (Gu & Song, 2022, Theorem 4.21).
E.4 ANALYSIS OF
Algorithm 13 implicitly maintains the primal-dual solution pair (x., 8) via representation Eq. (12)(13). It also explicitly maintains (Z. 8) E Tooe TR such that [ Fillz, I 7 and || 8 sil < tow; for all [r]] with probability at least 0.9.
Lemma E.8 (Correctness of CENTRALPATHMANTENNAN
Proof. Same proof of Lemma D.13.
47

---

Published as a conference paper at ICLR 2025
Lemma E.9. We bound the running time of CENTRALPATHNAANTENNNNEE PANNAAAPPPTMMANNNNNNNEENNTTTTTAT takes ((m(k"-1 + m. mm time.
as following.
If SATTAAAAEEEELLEEENNNNNAETTAAAAAA is called N times, then it has
total running time
Ó((Nn 11/2 (og(4max/mmin)) nl + m) (w+t))/2)..
CANTTAAPAATTMMNNNNNNNAAACCOOOOTT takes Ō(n(k: + m)) time.
Proof. INITIALIZE part: By Theorem E.3 and E.6. OUTPUT part: By Theorem E.3.
MULTIPLYANDMOVE part: Between two restarts, the total size of Lal returned by approx. QUERY is bounded by O(q2 (px,.)) by Theorem E.6. By plugging in 6x = 2a, Capx,2 E, we have (g Similarly, for 8 we have Eee(gl [ ( = (g)²)
ExEly T
Update time: By Theorem E.3 and E.6, in sequence of updates, total cost for update is (k2 m²)). So the amortized update cost per iteration is Ō(g(k² + m²)). The total update cost is
number of iterations time per iteration O(Nq(k2 + +²)).
We restart the data structure whenever K > or [ - > tcg, so there are O(N/q + log(tmax/tunin ') restarts in total. By Theorem E.3 and E.6, time cost per restart is O(n(k=11 + m)- So the total initialization time is
Init/restart time:
number restarts time per restart = ((NN/q + log(1 maltm/mm/mmmm//ee men (et n(hw + m m)=
Combine everything: Overall running time is
Ō(Ng(k² + m²) + (N/q + log((max/ttmmg tmin) 1) n(kww1 mw-1))
Taking & = 2/2, the optimal choice for is
(02))-12(()))) q = n/ (k2 + m²) 1/2 (k + m2-12112
achieving overall running time
Ō((Nn-1/3 1/2 + (og(fmx//mmm)) nkk + 1/2 + m 121111
(w+ =((Nn-11 11/ log(4max/ (aiix//maa)) nik: ls mn( (w+))
Proof of Theorem E.2. Combining Lemma E.8 and E.9.
E.5 PROOF OF MAIN STATEMENT
Proof of Theorem E.I. Use CENTRALPATHMANTENNNCCE (Algorithm 13) the maintenance data structure in Algorithm 20. Combining Theorem E.2 and Theorem F.1 finishes the proof.
F ROBUST IPM ANALYSIS
In this section present robust IPM algorithm for quadratic programming. The algorithm is modification of previous robust IPM algorithms for linear programming Lee et al. (2019); Lee Vempala (2021).
48

---

Published as a conference paper at ICLR 2025
Convention: Variables are in blocks of dimension n, (2 [n]).. Total dimension is Neot EeE ny We write (2) (Z1.... ..l.....00 Zn) " where 2; 12 We consider programs of the following form:
min #ER"
1
+ cc+
(16)
2
s.t. Ax =b=
E EK;
vi E [n]]
where Q Strot C RRottt, R.x b Rm, Kr
C R is convex set. Let = Iliemn Ki.
Theorem F.1. Consider the convex program (16). Let 0; K2 + R be a s:-seff-concortant barrier for all i E [n]. Suppose the program satisfies the following properties:
Inner radius r: There exists E Riot such that A= and B(2,r) K.
Outer radius R: K _ B(0, R) where E Bo
Lipschitz constant L: 11Q122-2 & L. ||cl2 SL.
Let (w)))e(m) 12 and Eie(n way, For any < < 2.2 Algorithm 19 outputs approximate steps, satisfying
solution x in ( log n log nilt er
1 T Qx c. x Ar=b,zE
r
ELR(R+ 1), AAx bll2 S (e| 11b|12).
min
(2" FxTQo+cTa
2
x E K.
Algorithm 19 Our main algorithm
1: procedure ROBUSTQPIPMIQ S C E Retot RmxRR E m (òi K.
R R)eii)) ml) W E RTt)
2:
7 Initial point reduction %/
3: p LR R + 1), a (0)
<- arg min, Eve(n) midh(x;),, (00) <- ep(c =p((c+Qx(0)
(0)
(0)
Q Epp 0
4:
1 +-
0
A - [A] 6 Az)
w
5:
0. divi E [n]], 2 (x) = - -log x - log(2 - x:)
1
6: (I, + CENTERING (o.)e W.E. tstart = 1. tend 2 1x
7:
return (Zim, 31m)
8:
end procedure
F.1 PRELIMINARIES
Previous works on linear programming (e.g. Lee al. (2019), Lee & Vempala (2021)) use the following path:
s/t Vou (x) μ, Ar b,
AT y + =
where (x) = Ei= (xr))
For quadratic programming, we modify the above central path as following:
8/t + V th (x) 14
Ax b, Qu A y +
We make the following definitions.
49

---

Published as a conference paper at ICLR 2025
Algorithm 20 Subroutine used by Algorithm 19
1: procedure CENTERING(Q E Stho A E M (di K; - R)ie(n), w E R E
R 8 E R Estart E -Bsotten Lend R2o)
2:
1 Parameters #/
3: À 64 1og(266n EiE(n) w:), E 1440 X,
4:
G (minig]n) with: ) h = 64VE
5:
1s Definitions
6:
)u(x) = Eie(n) (x;)
7:
11 = =x/l + (o)) Vi E [n]
Eq. (17)
8:
7i(,,8,1) - 11, 8) . Vi E [m]
Eq. (18)
sinh( A m.(2,s,))
9:
a(I,8,t)= 70(2,5,t) Eyetm "5"" cosh" (x,,)))))
Vi E [n]
Eq. (22)
10:
H = (z)
> Eq. (24)
Bua = Q. t1l tlw
Eq. (25)
12:
Pu.e,:== = B AT (AB 1, AT) - w.ztt
Eq. (26)
1E Main loop
14:
7+ 1 5- fstart. I t 8
15:
while > fend do
16:
Maintain T, l such that [o| -il1T, S, ||: -sili tew; and It - 8 S Eet
17:
6,,5 +- c(x, 8, 1) p4(x, 8, 1), vi [n]
Eq. (21)
18:
Pick 82 and 8, such that Ad = 0, 6, Qd Rango(AT) and
| - iB 1/2 I P m,a, (., B 1/2
w,,,7 5p w| S Ca,
w
|| ||E'd, 8, - (6p -H,,, 1/2 (1 P 1B, w f l| S Ea.
B-12
1/2
19:
t < max h)), tend), + 8mg +
20:
end while
21: return (2)
22:
end procedure
Definition F.2. For each i E [n]. we define the i-th coordinate error
(((x,,,) = 8i + w. (x;
(17)
t
We define 14's norm as
(12 = (2 s, I)E
(18)
We define the soft-max function by
V Vx(r) Σ
cosh (A
(19)
Wi
for some and the potential function the soft-max of the norm of the error of each coordinate
(1,, Vx(7(2,8,1)
(20)
We choose the step direction 6,p
8 -a -a-(=s)c((ssstt)) s,(xs,t)
(21)
where
sinh( 2 (x,
C1 8,1)
(22)
melo, s, tt) Eye(m wj wj² cosh? A 72(2,, 8, t))
50

---

Published as a conference paper at ICLR 2025
We define induced norms following. Note that we include the weight vector in the subscript to avoid confusion.
Definition F.3. For each block K, we define
||l|a, = ||
|ell = (+))))))))
for R": For the whole domain K. 11:=1 Krs we define
|| == (ollva
w/ w/||uall2 (1, 32/2
im1
||lli = [1]]][(33))))))))
w;
"(Ivallz, (2))
for 1hoo
The Hessian matrices of the barrier functions appear lot in the computation. Definition F.4. We define matrices Hz, xn; and H Ha, Reor Rtot as
H2 == V² (xx), H = (x)
(23)
(24)
From the definition, we that
H (w,., W. W.HH,i.
The following equations are immediate from definition. Claim F.5. Let H E Ro Koss be defined as Definition F.4. For R" we have
||ol| = 1//2 022 0l12,
| = 2/
Claim F.6. For each i E [n], let Hz, be defined as Definition F44. For E RP', E [n], we have
|ull, |1112,/211
|elle, ||, 0012--
We define matrices B and P used in the algorithm.
Definition F.7. Let A, Q denote two fixed matrices. Let Harz E Piot Khoot be defined as Defini- tion F.4. We define matrix B I Plus Moot as
Bu == Q Hwv
(25)
We define projection matrix P Purzet E Biltot Xfannt as
Pu.z, + B A (AB A) AB-1/2
(26)
F.2 DERIVING THE CENTRAL PATH STEP
In this section explain how derive the central path step.
We follow the central path
8/t ( P
51

---

Published as a conference paper at ICLR 2025
Ar
Qx + A y+s=c
We perform gradient descent on with step dμ. Then Newton step gives
8 V2 (z) = 6,
(27)
Ad 0
(28)
Q8 A by + 8, 0
(29)
where 6g (resp. 6y. 5g) is the step taken by (resp. y, s). For simplicity, we define E IR to represent V dw (x))
From Eq. (27) we get
5, = =ō - tHH
(30)
Plug the above equation into Eq. (29) we get
Q6 + A 8y + tō tH = O.
(31)
Let B Q + tH, multiply by AB- get
A6 + AB A 6, +
+AB 16, O.
Using Eq. (28) we get
AB-lATo 5y LAB 15, = 0.
Solve for & (assuming that AB-1A is invertible), we get
6, = - AB-18 8
Plug into Eq. (31) we get
Bo, - -A -BB-A 8, + 16, 0.
AB
Solve for 82 we get
(AB-1A 5, = =tB - -B- AAB-IAAAAAAAAAAAAA AB- 8
-12-1/2 - P)
where P B-1/2A (AB-(AA-AA))) (A AT)
AB- 1/2 -A
the projection matrix. Solve for 8, in Eq. (30)
get
8, = 16, - 1 (1-B-
In summary, we have
8, 1B-1/2 P)BB//1
8, = -AA AB- 8,,,
6, = t6, -B B-
P B 1/2 (AB -
These equations will guide the design of our actual algorithm.
IIn this section, and in this section only, we omit the subscript in H, B, P for simplicity.
52

---

Published as a conference paper at ICLR 2025
F.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, Line 18).
Assumption F.8. We make the following assumptions on O, E RRot and o₂ E RRos
100 iB-1/2 (1 P x,,,,))BB 1/2 B/2 1/220y wl J S Ea, w,z,t
177'5, 5, - -6p TH B-1/2 Par B S Eo.
The following lemma bounds the movement potential function assuming bound on 8y. Lemma F.9 ((Ye, 2020, Lemma A.5)), For any E R and R2 Let and X denote the parameters that are satisfying S a S 1 Let Ep E R denote a vector satiisfying
8X'
m 222 S a/8.
Suppose that vector T E Root is satisfying the following property
Ir r.
Wi
Vi
[n]
82
We define vector 8, E Rutt as follows:
-a- sinhi wy
5r,3
Cra
Viin wj² cosh? 2 m (,,,,)
Then, we have that
V (1(+
ai 2
w
cosh? (2 "i))" ))2/ + aA()
E i=1
W 1/1//
Wl
i=1
The following lemma bounds the norm of 6
Lemma F.10 (Bounding norm of 6,²).
( aa.
Proof.
(3,3
Σ w
a Σ wi (7. 3,7) (1)
sE(m]
= a
Σ
w (z. 8,7) 1 1/2 (1(T, 8, 8)1
xE(n]
o Σ
m. o (2r 8, (x,
iE(m)
w sinh 2 wy (,,,
= a² Σ
(x,
(7, Σ Σje(n) ue, coshi" ( 25 (2,3,,1))
53

---

Published as a conference paper at ICLR 2025
E Eje(mi W/ sinh" ( 1T, 3, ()) ΣjE(N) Wj coshi my 7s (,, 3,1)
a
S aaa
where the first step follows from Definition F.3, the second step follows from (T, (x, 8, 1) p1(2, 8, t). the third step follows from norm of F; (see Definition F.3), the forth step follows
from (x 8, 1) = 7 (1, 8, 1)12 (see Eq. (18)), the fifth step follows from ca(2, 8, 1)2
sinh (
(2,,3,1) ]M "ij cosh?" br ()( (see Eq. (22)), the sixth step follows from canceling the term
7: (T,,, 3, E), and the last step follows from cosh (z) 2 (x) for all x.
The following lemma bounds the norm of 82 and 5s.
Lemma F.11. For euch E [n], we define ax; == |60,s|1T7 Then, we have
9 Σ wra s 8
||5
(32)
iE(m))
In particular, we have ag < a.. Similarly, for da we have
17 (18m))))))))
Σ w;
a.
(33)
iE(m)
Proof. For Su, we have
16 3 ||tH| 1/2
10.5 B 1/2
B-12 P R, ) 1/2 1 8ul12 + Ea
1.z,
$ |1e /2 - P PB...)BB11 1,5,7 8y + Ea
11 1/2 B W
B-122
5 2112 + Ca
| 1123
1/ 6ul12 + Eα
0 + 7a
VA
0
First step follows from Assumption F.8. Second step is because i. 3 B Third step is
because P a. is a projection matrix. Fourth step is because T. J B B,,,77 Fifth step is by Lemma F.10. Sixth step is because -
For 8g. we have
|| S |1ō + |1 Hw..ZBB12(11
1/2 w.7,7 (w. "B, BB2 1/2
W.F.F. " 11 Eail
$ a7 + at + Eat 17 at
First step is by triangle inequality and the assumption that
8, ~J t8, - t² Hw, W.E.t
1/2 1-2 P P,. B-BBB11 1/2 w,
1/28
Second step is by same analysis as the analysis for d... Third step is by t <
33 and ج
<
32
32
The following lemma shows that new is close to μ + 6,, under an approximate step. Lemma F.12 (Variation of (Ye, 2020, Lemma A.9)). For each i [n], we define
3 ||
54

---

Published as a conference paper at ICLR 2025
For each i E [nl,, let
Ms (ewz s 14((x,stt) + 5n,s + Gp..-
Then, we have
wi 12) <
156a.
Proof.
The proof similar as (Ye, 2020, Lemma A.9), except for changing the definitions of 6 and
2
61 == H H1/2
w 5 a E. H 1 1-12
1/2
B-12 N.... _T, W p.:
W.
w
62 t H n
1/2 8, 2
t.. B-2
1/2 u.5..
P u. 1BB33 10.5.
One key step in the proof of Ye (2020) is the following property:
5μti 7ー 5,,i + Hw.z
+
Under our new definition E and e2, the above property still holds. Remaining parts of the proof are similar and we omit the details here.
The following lemma shows that error ya(ar, 8, 1) on the robust central path close to error p( 8, I) on the ideal central path. Furthermore, norms of errors 7xx, and 7 (I, 3, are also close to each other.
Lemma F.13 ((Ye, 2020, Lemma A.)))). Assume that r(x, 8, 1) S for all i. For all i E [n]. we have
(1, p,( 3, 3Ew
Furthermore, we have that
h/(2, s,t) -Va (T 1,7,3,7 7. < 5ew,
Proof. Same as proof of (Ye, 2020, Lemma 10).
The following lemma bounds the change of y under robust IPM step. Lemma F.14 ((Ye, 2020, Lemma A.12)). Assume ((I, s, S cosh(A). A). For E [n], we define Cry == (ew geew) 7/(x, 8,t) + a c, ((,,,,7)
Then, we have
w. 2/12 < - da + 4. max(w, (xx s,t) a.
(E(n)
Proof. The proof is similar the proof of (Ye, 2020, Lemma A.12). By replacing corresponding references in Ye (2020) by our versions (Lemma F.11, F.12, F.13) we get proof of this lemma.
Finally, the following theorem bounds the movement of potential function & under one robust IPM step.
Theorem F.15 (Variation of (Ye, 2020, Theorem A.15)). Assume 4(x, s, s cosh(A/64). Then for any 0 < h < 64 /1 urses we have
ad
new new new S (1
4 + ad
Σ w
Vieei) wi
kE(n)
In particular, for any cosh(A/128) < ( s, S co6)()/64), have that Φ( new (uew) < &t),s
Proof. Similar to the proof of (Ye, 2020, Theorem A.15), but replacing lemmas with the correspond- ing QP versions.
55

---

Published as a conference paper at ICLR 2025
F.4 INITIAL POINT REDUCTION
In this section, we propose an initial point reduction scheme for quadratic programming. Our scheme closer Lee al. (2019) rather than Ye (2020); Lee Vempala (2021). The reason is that Lee & Vempala (2021)'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 F.16 (Nesterov, 1998, Theorem 4.1.7 and Lemma 4.2.4)). Let o be u-seof-concordant barrier. Then for any I, E dom(@), we have
(Vo(z), 2) Su
(Vo(y) - Vo(z), - x) 2 |1 1
+ 11 11
Let x= arg min, (()). For any x E R such that || x² 1|=- I 1, we have that x E dom(().
Lemma F.17 (QP version of (Lee al., 2019, Lemma D.2)). Work under the setting of Theorem F.I.
Let (0)) arg min, Eigjn] (z..( Let LRiR+Ij. For any 0 < < 21 the modified program
min Ar=boEKRRgo
x Q
with
cpp
0]
A= [A|b Ax:
b= b,
EpC
0
0
satisfies the following:
2
T = 0 E RD and 3 = Teple (0)
Σ
are feasible primal dual vectors with
$ + (
where o (x) Ei= w.d, (T.) log(Tnn+).
For any Te E K. x R satisfying AT = and
0 Q + +
x Qr +
min
(34)
2
the vector Imm Zion is the first n coordinates of T is an approximate solution to the original convex program in the following sense:
1
(1 Qx + xTxz)) + tep
T lam Qui.n + Fim S min
2
A==b,aEK
||T
bll S 3c (R|| [[b]]]),
Fin E K.
Proof. First bullet point: Direct computation shows that (x, y. 8) is feasible.
Let us compute || + V v E- We have
3 (T) = |lep(c+ (17)d())))))))
Lemma F.16 says that for x E R" with ||x - (0)) (1 I 1. we have x K. because
(0))) = arg ming on (2). Therefore for any u such that v. V2 (x. )D < 1, we have (0))) + u E K.
and hence (|z + ul[2 < R. This implies ||l||2 I R for any v (z(0)) s 1. Hence
(z = R² I. So we have
|| (.) llep(c
s EpRllc+ + Qz (0)22
<pR[[[c[|2 ||Q| 2- 11:20
56

---

Published as a conference paper at ICLR 2025
<pRRLL+ LR)
$ E,
Second bullet point: We define
N fe(QotoTo)
OPT: = min
(35)
Az=
OPT
min Az=b,TeXXRgo
Qr Qx+c=x)
(36)
For any feasible the original problem (35), x Therefore we have
is feasible in the modified problem (36).
OPT < 5p(5" 1 Qu + cJx) =p OPT.
Dimm
Given a feasible Z satisfying (34), we write I
for some T 2 0. Then we have
- 5 Ii Quin + c (mm) + S OPT + e S Ep OPT
Therefore
1
20 Im Qunn + t Tim s OPT +ep
We have
r $ -epl 1:n Quin + Eim) + ep- OPT +e² s 3e
because 1 Qx + c a| LR(R + 1) for all x K.
Note that V satisfies 13 + (b Ax Yr So
|A b11 S |bb Az 7.
This finishes the proof.
The following lemma is generalization of (Lee et al., 2019, Lemma D.3) to quadratic program, and with weight vector w.
Lemma F.18 (QP version of (Lee al., 2019, Lemma D.3)). Work under the setting of Theorem F.I. Suppose we have &. + wy (z;) Mi for all i E [n]. -Qr + A A + = and Ar = b. If || I S wi for all i [n], then we have
z x Qu + T x 5 Qx" + tT 4
Tx
where x² arg Ae (BoTQr QTTQrre x.).
Let = (1-a)x+ar* for some a to be chosen. By Lemma F. 16, we have (Vow x" Ia) k. (Note that &w is k-self-comcorddat barrier for K.) Therefore we have
Proof.
k0
(Vow z²)
a
(Von (Maa) Vohno (2), Sa - x/) + (p , x²)
|xa, xz 122
1
+ (1, - x)
- (c- A ATy y
Qz, x))
|a,i 2z 2iller
VE [x]
57

---

Published as a conference paper at ICLR 2025
a² a² |0 34 11l
I 1| xilla. & (c Qx, x² - -x).
+ 2 =illa
RE(n]
First step is because vv (Fa), -Ra) I Second step because y == + V Third step by Lemma F.16 and Qu A + Fourth step is by Cauchy-Schwae and Ar Ar.
So we get
(x cx+c(x)
a||= x;ll2,
$ (x Qx" xxxa)) + 1-a + Σ
xilla.
Σ wi
1 alla;
- -2illz;
iE(n]
iEin]
1 1
a ||
2g 11:
$
2
Qr
x Qr + +c c.T. x²) +
+ Σ
w/ - 24 Ila:
W
- a
a |
2 18
iE(n)
(E
1
1
203
2 11
Qr
Qx" r +c'x") +
Σ
W0
a
a
|
iE(n)
1.1 1 $ 2
Qr
I 2
TQ= +c x²) +
T Qr
K
+
W/
1 0
iE(n)
1 1 Qx 7"0 x Qx +c=x") + a) $ 2/22"
First step by rearranging terms in the previous inequality. Second step is by AM-GM inequality and || ll. I w;. Third step is by merging the last two terms. Fourth step is by bounding the last term. Fifth step is by Eie[m] w; Lie(n] Wiy4 =
Finally,
2 r Qr t c x5 xx Txx" TTx kt
a(1 a)1
1 n50000 Qx" +c + 4xt.
First step is by rearranging terms in the previous inequality, Second step is by taking 2. This finishes the proof.
F,5 PROOF OF THEOREM F.
In this section combine everything and prove Theorem F.1.
Proof of Theorem F.l. Lemma F.17 shows that the initial and satisfies
|| S
This implies w. Il,, <
because w; 2 1.
Because e S 1. we have
(xr, cosh(Xu; 1144 coshs(1 - cosh(A/64)
sE(n
for the initial x and by the choice of X.
Using Theorem F.15, we see that
(xx, s.t) S cosh(A/64)
during the entire algorithm.
So at the end the algorithm, we have w/ 1 lullli;; < for all E [n).
for all E [n]. In particular, || S w;
64
58

---

Published as a conference paper at ICLR 2025
Therefore, applying Lemma F.18 get
2 x Qz Qx+cc
x
x Qu" +c=x= 4tk
x TT Qu" +c +
where we used the stop condition for at the end.
So Lemma F.17 shows how to get an approximate solution for the original quadratic program with error eLR(R+ + 1).
The number of iterations is because we decrease é by a factor of 1 - every iteration, and the choice GIVK'
G GAUSSIAN KERNEL SVM: ALMOST-LINEAR TIME ALGORRTHM AND HARDNESS
In this section, we provide both algorithm and hardness for Gaussian kernel SVM problem. For the algorithm, we utilize result due to Aggarwal & Alman (2022) in conjunction with our low-rank QP solver to obtain an O(n!~o(1) l-(()) (og(1 time algorithm. For the hardness, we build upon the framework outlined in Backurs al. (2017) and improve their results terms dependence on dimension d.
We start by proving a simple lemma that shows that if UVT for low-rank U, V, then the quadratic objective K (yv") also admits such factorization via simple scaling.
Lemma G.I. Let U, V E gank and y E R2 Then, there exists pair of matrices U, V gnxk such that
UV (UVT) o (yy")
moreover, Ü, V can be computed in time (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
Ao (yx) = D,
where Dy, Dr E B are diagonal matrices that put y, x on their diagonals. Thus, we can simply compute U, V as follows:
Ú = D
V D.V.
consequently,
ÚVT DrUVT Dj
(yy ) (UVT) = (UVT)= (vv
as desired. Moreover, the diagonal scaling of U, can be indeed performed in Onnk) time, as advertised.
Throughout this section, we will let denote the squared radius of the dataset.
G. ALMOST-LINEAR TIME ALGORITHM FOR GAUSSIAN KERNEL SVM
We state result due to Aggarwal & Alman (2022), in which they present optimal-degree polyno- mial approximation to the function e e" 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=:
59

---

Published as a conference paper at ICLR 2025
Definition G.2. Let [0, B] - I, we let (B:e(f) N denote the minimum degree of
constant polynomial p(z) such that
p(z) - f(x)] s =E|O,
sup
Utilizing the Chebyshev polynomial machinery together with the orthgonal polynomial families, Ag- garwal & Alman (2022) provides the following characterization on qB (f)):
Theorem G.3 (Theorem 1.2 of Aggarwal & Alman (2022)). Let B 2 and 6 E (0, 1). Then log(L/e) gB: (e: O(max v B (og(e/)), log( B log(1/e)) 1)
Theorem G.4 (Corollary 1.7 of Aggarwal & Alman (2022)), Let , E Rd be dataser with squared radius B and e E (0, 1). Let (e: Let K E Ruxa be the Gaussian kernel matrix formed by . Tn- Finally, let k = (2d+20). Then, there exists deterministic algorithm that computes a pair of matrices U, V R such that for any vector V E R",
|Ko UV u/lso S cclool1..
Moreover, matrices U, V can be computed in time O(nkd).
Even though &oo error in terms of & norm 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 G.5. Let K R PSD kernel matrix and (0,1) be a parameter. Let K gaxn be an approximation to K with the guarantee that for any E IR",
|| Kullx S Sellu|11,
then
10 Ke Ke- - Kel < e||o|| e
Proof. The proof a simple application of Hölder's 's inequality:
|u (Kr Ke) (0,, Kr Ke)
||1|/2|1 [ullallKo Kulloo
$ <||0|]?
S en|lu|la
where the second step is by Hölder's inequality, and the last step is by Cauchy-Schwarz. 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 G.6 (Restatement of Definition 1.3). Given data matrix X guxd and labels e R.. Let Q E x denote a matrix where Qi,; = K(x,, xz)) y:3y for i,j E [m]. The hard-mragin kernel SVM problem with bias asks to solve the following program.
max wER"
1, a TQQ
s.t. y=
a 220
Theorem G.7. Let Gaussian kernel SVM training problem be defined above with kernel function K(æ;, 2z)) = exp( -||2; x2,12). Suppose the dataset has squared radius B 2 1, and let e E (0, be the precision parameter. Suppose the program satisfies the following:
60

---

Published as a conference paper at ICLR 2025
There exists point E R such that there is Euclidean ball with radius centered at that is contained 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 E RR" such that 0, moreover,
1 TTQa OPT
1 1 a
2
l 111 S
where OPT denote the optimal cost of the objective function. Let (0:0(c/nR3)(e and k (2d+29). Then, the vector â can be computed in expected time
Proof. Throughout the proof, é1 O(e/(nn²)), We will craft algorithm that first computes an approximate Gaussian kernel together with a proper low-rank factorization, then use this proxy kernel matrix to solve the quadratie 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 The- orem G.4 using data matrix with accuracy parameter ei. We let K UVT denote this approximate kernel matrix, and we let Q D2UVT D, denote the approximate quadratic matrix. Owing to Lemma G.5, we know that for any vector x E R"
x T (Q - ))11 1(D T(K (D)
S 1 D
where we use the fact that E {=1)". This also implies that
| - Q11 S 617
(37)
this simple bound will come in handy later.
Solving the approximate program to high precision. Given Q, we solve the following program:
1
Qa
xER
s.t. y= 0 2 0
by invoking Theorem E.1. To do so, we need a bound on the Lipschitz constant of the program, i.e.. the spectral norm of Q and & norm of 1. The latter is clearly vn, we shall show the first term is at most (1 (1) n.
Note that
|Q = ID, S tr(D,KD,)
tr[K]
<
where we use K is PSD. Combining with Eq. (37) and triangle inequality, we have
| s ||인|| + Q ~1
(1 4)
61

---

Published as a conference paper at ICLR 2025
With these Lipschitz constants, examine the error guarantee provided by Theorem E.1: it produces a vector IR" such that
T a Q6
(1, 2" T Qa) - O(emnR))
max at a=
| ull S O(1mR),
we mainly focus on the first error bound, as we need to understand the quality of T when plug into the program with Q.
We will follow a chain of triangle inequalities, so we first bound
Ja (Ô -Q)ai S en en|lãll
S RnR
Next, let
a' == arg max
1 a
Qa,
y=0,220
2
1 2
max a" = arg
1
Qa,
a 0=0,220
then we have the following
1 a T Qa' 2 a, [aa)) Qa
2
2
1 a- 2 (aat)
Qa ((anR²)
OPT ((EInR))
where the second step by applying Lemma G.5 to a Now we are ready to bound the final error:
1 a
1 a Qa O(anR?3
a
1,
2 a a Qa' O(emni²))
TTQQaa
OPT O(Ein(²)).
The final error guarantec follows by the choice of E and indeed design an algorithm that outputs a vector with
a TTQQ Qa 2 OPT
2
la 1ll E.
Runtime analysis. It remains to analyze the runtime of our proposed algorithm. We first compute an approximate kernel K with parameter c,, owing to Theorem G.4, we have
log(nR/e)
[B:c; () (mmax BB Blg((R//c) log B log(nR/c))
=
then by setting h = (2d+20), the matrix K can be computed in time O(nkd). Given this rank-k factorization, the program can then be solved with precision Ej in time
2(nkl++1/1lognnll
desired.
Remark G.8. To understand the value range of k, let consider the set of parameters:
Ollog n., 1/ poly n, R poly n, B O(1),
62

---

Published as a conference paper at ICLR 2025
under this setting, O(log(nR/e)) Ollog n)) and the degree is the rank is then
q O(Vlogn)
2d + 2q
k =
29
((log n)3 vlogn)
0 log nvlogn)
consequentially, our algorithm runs in almost-linear time in n:
Ö(n ( log
is worth noting to achieve the almost-linear runtime, the data radius can be further relaxed. In fact, as long as
log
B =
log log
we can ensure that & = (1) and subsequently the almost-linear runtime.
The runtime we obtain can be viewed as consequence of the "phase transition" phenomenon observed Aggarwal Alman (2022), in which they prove that if B w(lg n). then quadratic time in is essentially needed approximate the Gaussian kernel assuming SETH.
G.2 HARDNESS OF GAUSSIAN KERNEL SVM WITH LARGE RADIUS
In this section, we show that for (log n). any algorithm that solves the program associated to hard-margin Gaussian kernel SVM would require 2(n2 time for B = wlog n). This justifies the choice of B in Remark G.8. To prove the hardness result, need introduce the approximate Hamming nearest neighbor problem.
Definition G.9. For of vectors, and let & E {0, I,..., d} be threshold. The (1 o)-Aprroximate Hamming Nearest Neighbor problem asks distinguish the following two cases:
and d let fai..... an) (bI..... b,, S 10, be two sets
If there exists some O; and by such that ||a; byll I t, output "Yes";
If for any i,j E [m]] we have ||a; b (1 + 8) 4, output "No".
Note that the algorithm can output any value if the datasets fall in neither of these two cases. will utilize the following hardness result due to Rubinstein.
Theorem G.10 (Rubinstein (2018)). Assuming SETH, for every > 0, there exists 5 > 0 and C > such that (1 Apprrxxmmaee Hamming Nearest Neighbor in dimension d = C'log requires time 2(n2-9)
final ingredient is rewriting of the dual SVM into its primal form, without resorting to optimize over an infiniteddimensional hyperplane.
Lemma G.11. Consider the dual hard-margin kernel SVM defined as
1 aER
max
Σ
a.a.u,,K(w..wj)
[.jEm]]([b]
s.t. a y 00
a a20. O.
The primal program can be writlen as
min aER 2
1
Σ
K(w,,, w)))
63

---

Published as a conference paper at ICLR 2025
s.t. yef(w;) 1,
a 2 0,
where f Rd + R is defined as
f(w)= Σ aouj w) - b.
1=1
Moreover, the primal and dual program has no duality gap and the optimal solution to both programs are the same.
Proof. Recall that the primal hard-margin SVM is the following program:
min 2 11ul22
s.t. y: ((w,) b) 2 1,
where ò E R is the bias term and é Rd - RKKK is the feature mapping corresponding to the kernel in sense that K(w,, = ((w;) ((w)). The optimal weight e Σi=1 (a))(())))) where or IR is the optimal solution to the dual program. Consequently,
||u12 (w,))
Σ
a;a,)1y;((u;) oww))
Σ
a.a,,,,,K wz)
[jE(m|x[m]
Qa,
where the matrix Q is the usual
Q (yyT) K,
the constraint can be rewritten
(i(u ((w;) 9:(( a:y:b(w:)) (eo((w))T<(w)))) b)
(5b))(w3)) ((ww))) (T)()) - yeb
ya
ayy, K(w,, wy)) 4:b
where / Rdd - R is defined as
f(w) = Σ
Koy,K(ws, b.
1=1
Thus, alternatively write the primal
min R 2
a Qa.
T
s.t. vif((w:) 2
For the strong duality and optimal solution, see, e.g., Muller et al. (2001).
64

---

Published as a conference paper at ICLR 2025
We will prove the almost-quadratic lower bound for Guassian kernel SVM. Our proof strategy is similar to that of Backurs et al. (2017) with different set of parameters. It is also worth noting that the Backurs et al. (2017) construction
Requires the dimension d = (llgg n2); Requires the squared dataset radius B = Olog"* n).
We will improve both of these results.
Theorem G.12. Assuming SETH, for every > 0, there exists hard-margin Gaussian kernel SVM without the bias term as defined in Definition 1.3 with d = @(logn) and error ε = exp( (sp(=O(log' n)) for inputs whose squared radius is at most B = (log n) requiring time (2(2²99) to solve.
Proof. Let = -(c'))))o n. We will provide a reduction from (1 + Approxxxaee Hamming Nearest Neighbor to Gaussian kernel SVM. Let A := {a),..., an J. B == (bi..... bb) S {0, 1j be the datasets, we assign label to vectors a, and label -1 to all vectors by. moreover, we scale both and by this results in two datasets with points in (0.1)4 The squared radius of this dataset then
B max max lla la 112. max ||b, 1by max |la, 16,112) 1.j
< 12d
=66 log" n).
To simplify the notation, we will implictly assume A and B are scaled by / without explicitly writing out la,, lby. Now consider the following three programs:
Classifying A:
min aERg
1
Σ
asa,K(ai,az)
2
s.t. Σ ==1
a()
Vi [n]
(38)
. Classifying B:
Σ b,),
min SERe
s.t. B,K(b,,))) -1, [i
(39)
Classifying both A and B:
min of
Σ
((j,)) + 1
Σ B,8,K(b,,b))- Σ
a,3,K((s,bb))
2
20
xjjE/n|([n]
[JE[n]][N]
a K(aK,, - Σ 3yK(as, by) 21,
We [n],
j=1
a
Σ
LayK(b, ag)
Σ (,K(b,,,))
Vie in]
(40)
We will prove that the optimal solution at's 's and 3t's are both lower and upper bounded. Use Va)(A), Val(B) and Val(A, B) to denote the value of program (38), (39) and (40) respectively, then note that
Val(A) < 2
2
65

---

Published as a conference paper at ICLR 2025
by plugging in 1 and setting all vectors to be the same. On the other hand,
1
Val(AA) 2 (a))
(aa))
2 1 (a)))
(a))
2 i=1
Combining these two, we can conclude that for any aç, it must be ag S n. For the lower bound, consider the inequality constraint for the i-th point:
of Σ a)) 21,
jti
to estimate K(a,, ay), note that ||a a, 18 |a; 0,1l2 2 1 for j * i,9 and
K(a;, az) exp( -121a,; a l12)
exp(- -((c)) log logn/la, a,ll1)
I exp(-22 (xp(-2((c6) -1 logg
S n 10 10-10/1100 /100.
combining with 0z S n, we have
of 2 1 Σ o(;k()))))) jti
10 1000
n
2 1/2.
This lower bound on of is helpful when we attempt to lower bound Val(A, B) with (Al(A)+ +aa(()) Following the outline of Backurs et al. (2017), we consider the three dual programs:
Dual of classifying A:
max EERA
Cj
Σ in
ajayK(an,
(41)
2
20
Dual classifying B:
max Σ B
Σ
3.B,K(b,,bb)
(42)
BERR
2
1=1
b2
Dual of classifying A and B:
max B.
3
Kaa,Klai,
1
Σ b
B,B,K(b.,bb) + Σ 3 Ko,8,K(a;,40 bj)
2
2
2.3
(43)
as the SVM program exhibits strong duality, we know that the optimal value of the primal equals the dual, so can alternatively bound Val(A, B) using the dual program. Plug in a", 8 to program (43), we have
Val(A, B) > Σ
Σ
B 2
Σ
KajkK(a,,,, - Σ 3, by) Σ K0j
2=1
1=1
Val(A) + Val(B) Σ Ka
i.j
We without loss of generality that during prepracess, we have remove duplicates in and B.
66

---

Published as a conference paper at ICLR 2025
to bound the third term, consider the pair b such that ||as 630 l12 s t 1, and note that
Eat by) 2 ag 87 12,K((t ((asos 630)
Σ 4
exp((2(cc)) -log (t 1)).
To wrap up, have
Va)(A. 2 Val(A) + Val(B) exp( -log (t 1))
We now prove the "No" case, where for any ai, b3. ||; byll L. We have
K(a; by) exp(-12 by1
exp(-2(c' log t),
t), set a' a* + 10m2m and 3' 8" + 10n2m, we let V to denote the value when evaluating program (40) with a', f3. We will essentially show that
let m
Val(A, SV
and
V Val(A) + Val( + 400m m
chaining these two gives us a certificate for the "No" case. To prove the first assertion, we show that a', ,² are feasible solution to program (40) since
Σ a))((an)))=
Σ ((zikaanaa) 10n'm0(aa, a,)) j=1
= a + 10n2m
Σ [a} 10n 10n²m))(an,0))
2 a + Σ
a
a K(a;; a,) + 10n2m
Ji
10n2mm m +
Σ
a a3(((a,
j=1
2 10n2m.
where we use of satisfy the inequality constraint of program (38). We compute upper bound on Σj=1 B,K(a,,b;).
Σ 8/K(a, by) S Σ
2nm
j=1
= =nn
where we use the fact that = -cxp(_((c)) n f) S n n-10/10 therefore 3 + 10n2m S 2n. Thus, must be the case that
12
ac) Σ ByK(a,,b)) (,),,,,)) 2 8m2m +1
j=1
21,
desired. The other linear constraint follows by symmetric argument. This indeed shows that a' Bt are feasible solutions to program (40) and Val B) < V.
To prove upper bound on V. note that
V Σ aaa
a.) + 1 Σ
K, b,) - Cd'B,K(ar.b)
67

---

Published as a conference paper at ICLR 2025
E afa) + 1 + B,8gK(bb..b))
we bound the first quantity, as the second follows similarly:
= Σ
(aja} + m0n3m + aa) + 110000 (00mmm²)K(())))) az)
I Val(A) + Σ In0n az) + Σ
100n m K(a,, az)
s Val(A) 10n'm + m0
S Val(A) 200n
we can thus conclude
S Val(A) + Val(B) + 400n
Chaining these two, we obtain the following threshold for the "No" case:
Val(A, < Val(A) + Val(B) + 400n m.
Finally, we observe that
400n exp(-2(c(5) -ogn logn 1) < exp(-2(c66) llogn n (t 1)),
we can therefore distinguish these two cases.
Note that when considers solving the program with additive error, need to make that the error is smaller than the distinguishing threshold, i.e.,
c < 4
exp( exp(-2(c6))- log -(t 1))
exp( ((0(-2(c)))
= exp( exp(-0(log n)))
where we use < d and d = (log n). This concludes the proof.
Remark G.13. Our proof can be interpreted as using stronger complexity theoretical tool in place of the one used by Backurs et al. (2017), to obtain a better dependence on dimension and B. We also note that the construction due to Backurs al. (2017) has the relation that B O(d logn), this is because in order to lower bound Val(A, B). one has to lower bound the optimal values of a 's and 13 s.. To so, one needs to further scale up a, and by's so that within datasets and B, the radius least (log n). This is in contrast the Batch Gaussian KDE studied in Aggarwal Alman (2022), where they show the almost-quadratic lower bound can be achieved for both d, B = Olog n).
Similar to Backurs et al. (2017), we obtain hardness results for hard-margin kernel SVM with bias, and soft-margin kernel SVM with bias.
Corollary G.14. Assuming SETH. for every q > O. there exists a hard-margin Gaussian kernel SVM with the bias term with d @(log n) and error exp(-x n)) for inputs whose squared radius most B Olog" 6 n) requiring time 2(n=299 to solve.
Proof. The proof is similar to Backurs et al. (2017). Given a hard instance of Theorem 12, except we append O(log n)) entries with magnitude O(log² n) instead of distributing the values across (log n) entries. Rest the proof follows exactly the same as Backurs et al. (2017).
Corollary G.15. Assuming SETH, for every > there exists soft-margin Gaussian kernel SVM with the bias term with d @(log and error xpoo= n)) for inputs whose squared radius at most B O(log" n) requiring time ((n2-9) to solve.
68

---

Published as a conference paper at ICLR 2025
Remark G.16. Compared to the construction of Backurs et al. (2017) in which they distribute total mass of (log" n) across (logg" n) entries so that they ensure after the reduction, the vectors take values in {-1, 0, 1}. we instead distribute the mass across (log n) entries so that each entry has magnitude (log* 2 n.). To make the reduction work, the total mass of (log" n)) is needed, and for Backurs et al. (2017), it is fine to append an extra @log" n) entries as their hardness result for hard-margin SVM without bias does require at the price each entry has a larger magnitude of ( nt). This blows up the squared radius from log² n to log n. We note that the Backurs et al. (2017) construction has squared radius log" n.
( n). For us, we need to restrict d = (log n)
69


Reviews:
[{"review_id": "nl1nwQOQtm", "reviewer": "Reviewer_Vy1p", "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.", "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.\n2. 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. \n3. 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.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 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.", "weaknesses": "1. No experiment to demonstrate the superiority to validate the new claimed contributions.\n2. 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. \n3. What if the kernel is a combination of kernels, say K_new=K1+K2 where K1 and K2 are Guassian kernel?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "PKuOh8s6Rc", "reviewer": "Reviewer_WSha", "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.", "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.\n\nThe 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?\n\nIs there a reference for the result on $T$ on line 291?\n\nWhat is the graph structure on the bags on line 310?\n\nCan the authors give references for the sparse Cholesky factory result on Line 318?\n\nWhat 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?\n\n---\n\n[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.\n\n[2] Alman, Josh, and Ryan Williams. \"Probabilistic polynomials and hamming nearest neighbors.\" 2015 IEEE 56th Annual Symposium on Foundations of Computer Science. IEEE, 2015.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 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.", "weaknesses": "While this paper addresses a theoretical problem in complexity of SVMs, it'd still be interesting to see some actual experimental results.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Umv0k8MOSs", "reviewer": "Reviewer_6s2J", "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.", "questions": "1.\tIn 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.\n\n2.\tThough 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.\n\n3.\tIt 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. \n\n4.\tOverall, 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.\n\n5.\tSmall typos: for example, in line 1540, is it “memebers” or “members”? Please double check.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "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.\n\nRelevance: 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.\n\nComplexity 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.", "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.\n\nAssumptions 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.\n\nComparison 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "HK4v2lXj16", "reviewer": "Reviewer_CgdP", "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.", "questions": "see strengths and weaknesses.\n\n**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.", "limitations": "", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- The authors provide intensive mathematical analysis of robust interior point methods for quadratic optimization problems.", "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.\n- 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.\n- 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? \n- 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)