# Solving Dense Linear Systems Faster than via Preconditioning

Michał Dereziński and Jiaming Yang University of Michigan
(derezin@umich.edu)University of Michigan<EMAIL_ADDRESS>
###### Abstract

We give a stochastic optimization algorithm that solves a dense $n\times n$
real-valued linear system $Ax=b$, returning $\tilde{x}$ such that
$\|A\tilde{x}-b\|\leq\epsilon\|b\|$ in time:

$\displaystyle\tilde{O}\big{(}(n^{2}+nk^{\omega-1})\log 1/\epsilon\big{)},$

where $k$ is the number of singular values of $A$ larger than $O(1)$ times its
smallest positive singular value, $\omega<2.372$ is the matrix multiplication
exponent, and $\tilde{O}$ hides a poly-logarithmic in $n$ factor. When
$k=O(n^{1-\theta})$ (namely, $A$ has a flat-tailed spectrum, e.g., due to
noisy data or regularization), this improves on both the cost of solving the
system directly, as well as on the cost of preconditioning an iterative method
such as conjugate gradient. In particular, our algorithm has an
$\tilde{O}(n^{2})$ runtime when $k=O(n^{0.729})$. We further adapt this result
to sparse positive semidefinite matrices and least squares regression.

Our main algorithm can be viewed as a randomized block coordinate descent
method, where the key challenge is simultaneously ensuring good convergence
and fast per-iteration time. In our analysis, we use theory of majorization
for elementary symmetric polynomials to establish a sharp convergence
guarantee when coordinate blocks are sampled using a determinantal point
process. We then use a Markov chain coupling argument to show that similar
convergence can be attained with a cheaper sampling scheme, and accelerate the
block coordinate descent update via matrix sketching.

###### Contents

  1. 1 Introduction
    1. 1.1 Main Results
    2. 1.2 Our Techniques
    3. 1.3 Background and Related Work
  2. 2 Preliminaries
  3. 3 Main Algorithms
    1. 3.1 Fast Kaczmarz-type Solver
    2. 3.2 Fast Coordinate Descent-type Solver
  4. 4 Overview of the Analysis
    1. 4.1 Expected Projection under DPP Sampling
    2. 4.2 Replacing DPP with Uniform Sampling
    3. 4.3 Approximate Sketch-and-Project
    4. 4.4 Proof of Theorem 3.1
  5. 5 Expected Projection under DPP Sampling
    1. 5.1 Elementary Symmetric Polynomials
    2. 5.2 Proof of Lemma 4.1
  6. 6 Replacing DPP with Uniform Sampling
    1. 6.1 Preprocessing
    2. 6.2 Successful DPP Sample
    3. 6.3 Proof of Lemma 4.2
    4. 6.4 Comparison with Direct $k$-DPP Sampling
  7. 7 Approximate Sketch-and-Project
    1. 7.1 Approximate Update by Solving Sub-Problem
    2. 7.2 Proof of Lemma 4.3
  8. 8 Positive Semidefinite Linear Systems
  9. 9 Implicit PSD Linear Systems
  10. 10 Least Squares Regression
  11. A Omitted Proofs

##  1 Introduction

Solving a linear system $\mathbf{A}\mathbf{x}=\mathbf{b}$ is one of the most
classical algorithmic tasks, with applications in scientific computing,
engineering, data science, optimization, and more. Many classical
deterministic algorithms are known for this task, from direct methods such as
Gaussian elimination to iterative methods such as Conjugate Gradient. More
recently, there has been much success in developing faster algorithms based on
randomized techniques such as sampling and sketching for solving linear
systems under specific structural conditions, e.g., graph-structured [ST14,
KMP12, CKK+18] and other classes of sparse linear systems [XXG12, Gra06], as
well as extremely over/under-determined systems [Sar06, Woo14]. However, the
benefits of using randomized algorithms for general dense linear systems
remain unclear.

Most of the work on improving the time complexity of direct methods for
solving linear systems has focused on the reduction to matrix multiplication
via divide-and-conquer approaches. This has led to the complexity of
$O(n^{\omega})$ for solving an $n\times n$ linear system
$\mathbf{A}\mathbf{x}=\mathbf{b}$, where $\omega<2.371552$ [WXXZ23] is the
current best known bound on the matrix multiplication exponent, developed over
a long line of works [Str69, Pan84, CW87, Wil12].

On the opposite end of the spectrum from direct methods are iterative methods,
which, through a sequence of cheap iterations, construct estimates that
gradually converge to the solution. These approaches are widely used in
practice for solving large linear systems. One of the most famous of those
methods is Conjugate Gradient (CG) / Lanczos algorithm [HS52]. This method
under exact arithmetic recovers the optimal solution in $n$ iterations, each
costing $O(n^{2})$ time for a dense system. However, CG is more commonly used
as an approximation algorithm for well-conditioned systems, taking
$O(\log(1/\epsilon))$ iterations to converge $\epsilon$-close when the
condition number $\kappa(\mathbf{A})$ (the ratio between largest and smallest
positive singular values) is bounded by a constant. These dramatically
opposing views of the complexity of CG can be unified through careful analysis
of its dependence on the singular value distribution of $\mathbf{A}$ [SW09].
In particular, CG will take $k+O(\log(1/\epsilon))$ iterations to converge,
when $\mathbf{A}$ has at most $k$ singular values larger than $O(1)$ times the
smallest positive singular value, or equivalently, when
$\kappa(\mathbf{A}-\mathbf{A}_{k})=O(1)$, where $\mathbf{A}_{k}$ is
$\mathbf{A}$’s best rank-$k$ approximation. This leads to the overall
complexity of $\tilde{O}(n^{2}k)$, where we use $\tilde{O}$ to hide the
logarithmic factors.

The above discussion indicates that solving linear systems with a matrix
having only $k=o(n)$ large singular values should be substantially easier than
in the worst-case. Understanding this phenomenon is motivated by the real-
world prevalence of such linear systems: either due to direct regularization
(e.g., $\mathbf{A}=\mathbf{B}+\lambda\mathbf{I}$), as in continuous
optimization [BV04], machine learning [ZDW13], and statistics [DW18], or
because the input matrix is distorted by isotropic noise (e.g.,
$\mathbf{A}=\mathbf{B}+\delta\mathbf{N}$, where $\mathbf{N}$ is the noise),
which can be caused by measurement error [LW11], data compression [LGW+21],
and privacy enforcement [DR14]. Yet, the $\tilde{O}(n^{2}k)$ complexity of CG
is highly unsatisfactory, as it only provides an improvement over direct
methods for small $k$, and does not lead to a near-linear time algorithm
unless $k=\tilde{O}(1)$. This leads to our central question:

> What is the complexity of solving an $n\times n$ linear system with $k$
> large singular values?

Arguably the most popular strategy of accelerating the complexity of linear
system solvers such as CG is via preconditioning. For instance, given access
to $\mathbf{A}$’s best rank $k$ approximation $\mathbf{A}_{k}$, we can easily
construct a preconditioner matrix $\mathbf{M}$ that allows CG to solve our
problem in $O(\log(1/\epsilon))$ iterations (since $\mathbf{A}_{k}$ describes
the entire ill-conditioned part of $\mathbf{A}$’s spectrum). While computing
$\mathbf{A}_{k}$ exactly is as expensive as directly solving the linear
system, much recent work has been dedicated to fast randomized algorithms for
low-rank approximation via matrix sketching. In particular, we can use power
iteration [HMT11] (or Krylov subspace iteration [MM15]) to construct a low-
rank approximation $\tilde{\mathbf{A}}_{k}$ with a spectral norm guarantee
$\|\tilde{\mathbf{A}}_{k}-\mathbf{A}\|_{2}=O(\sigma_{k+1}(\mathbf{A}))$,
sufficient to construct a good preconditioner for our problem. These
algorithms require (repeatedly) multiplying matrix $\mathbf{A}$ with an
$n\times k$ random matrix, leading to the runtime of
$\tilde{O}(n^{2}k^{\omega-2})$. This complexity can be further refined using
fast rectangular matrix multiplication [LG12] by replacing the exponent
$\omega$ with $\omega(1,1,\log_{n}(k))\leq\omega$. Nevertheless, for large
enough $k$, the resulting preconditioning algorithm will not run in near-
linear time.111Using current matrix multiplication algorithms, we have
$\omega(1,1,\alpha)>2$ for any $\alpha>0.321335$ [WXXZ23]. For simplicity,
from now on we will present all time complexities using just the square matrix
multiplication exponent $\omega$. Some randomized low-rank approximation
algorithms avoid power iteration and run in $\tilde{O}(n^{2}+nk^{\omega-1})$
time [CW13, CEM+15, CMM17], but they only guarantee a bound on the Frobenius
norm error $\|\tilde{\mathbf{A}}_{k}-\mathbf{A}\|_{F}$ which is worse by up to
a $\sqrt{n}$ factor than the spectral norm bound and thus does not in general
lead to a good preconditioner.

In summary, existing approaches suffer from an undesirable trade-off between
the cost of preconditioning and the cost of solving the system, where one of
those costs has to be larger than the input size by at least a
$\mathrm{poly}(k)$ factor.

Problem | Condition on $\mathbf{A}$ | Our results | Theorem
---|---|---|---
Dense linear system | $\mathbf{A}\in\mathbb{R}^{n\times n}$ | $\tilde{O}(n^{2}+nk^{\omega-1})$ | Thm 3.2
Least squares | $\mathbf{A}\in\mathbb{R}^{n\times d}$ | $\tilde{O}(\mathrm{nnz}(\mathbf{A})+d^{2}+dk^{\omega-1})$ | Thm 10.1
PSD linear system | $\mathbf{A}\in\mathcal{S}_{n}$ | $\tilde{O}(\mathrm{nnz}^{*}(\mathbf{A})+nk^{\omega-1})$ | Thm 8.1
Table 1: Summary of our complexity bounds for solving linear systems with at
most $k$ large singular values to high precision, hiding
$\mathrm{polylog}(n/\epsilon)$ factors, as a function of the square matrix
multiplication exponent $\omega$. In the least squares setting, we consider
$\mathbf{A}$ to be a tall matrix, i.e., $n\gg d$, and in the psd setting, we
define
$\mathrm{nnz}^{*}(\mathbf{A})=n\cdot\max_{i}\\{\mathrm{nnz}(\mathbf{A}_{i,:})\\}$.

###  1.1 Main Results

In this paper, we avoid the undesirable trade-off between the cost of
preconditioning and the cost of solving a dense linear system with $k$ large
singular values, by abandoning the preconditioning step in favor of stochastic
optimization.

######  Theorem 1.1 (Dense linear system, simplified Theorem 3.2).

Given $\mathbf{A}\in\mathbb{R}^{n\times n}$, $\mathbf{b}\in\mathbb{R}^{n}$,
$\epsilon>0$ and a constant $C=O(1)$, we can compute $\tilde{\mathbf{x}}$ such
that $\|\mathbf{A}\tilde{\mathbf{x}}-\mathbf{b}\|\leq\epsilon\|\mathbf{b}\|$
in time:

$\displaystyle\tilde{O}\Big{(}(n^{2}+nk^{\omega-1})\log 1/\epsilon\Big{)},$

where $k$ is the number of singular values of $\mathbf{A}$ larger than $C$
times its smallest positive singular value.

######  Remark 1.1.

We note that $\tilde{O}$ hides only $\mathrm{polylog}(n)$ factors, and the
algorithm can be implemented without the knowledge of $k$. Theorem 3.2 further
extends to rectangular and inconsistent systems.

The above theorem implies that we can solve an $n\times n$ linear system with
$k$ large singular values in near-linear time, $\tilde{O}(n^{2})$, as long as
$k=O(n^{\frac{1}{\omega-1}})=O(n^{0.729})$. Moreover, when we compare the time
complexity $\tilde{O}(n^{2}+nk^{\omega-1})$ against the cost of CG with low-
rank preconditioning, the proposed algorithm entirely avoids the
$\tilde{O}(n^{2}k^{\omega-2})$ cost of power iteration, without sacrificing
the complexity of the solver (up to logarithmic factors). See Table 1 for a
summary of our results.

#### Least squares.

Our approach naturally extends to both over- and under-determined systems. In
particular, Theorem 3.2 directly shows that we can solve a tall $n\times d$
least squares problem $\min_{\mathbf{x}}\|\mathbf{A}\mathbf{x}-\mathbf{b}\|$
with $k$ large singular values, in time $\tilde{O}(nd+dk^{\omega-1})$. This
can be further improved for sparse matrices, by combining our algorithm with
existing sketching-based preconditioning techniques for tall least squares.
Below, we use $\mathrm{nnz}(\mathbf{A})$ to denote the number of non-zeros in
$\mathbf{A}$.

######  Theorem 1.2 (Least squares, simplified Theorem 10.1).

Given $\mathbf{A}\in\mathbb{R}^{n\times d}$ with $k$ large singular values,
$\mathbf{b}\in\mathbb{R}^{n}$ and $\epsilon>0$, we can compute
$\tilde{\mathbf{x}}$ so that
$\|\mathbf{A}\tilde{\mathbf{x}}-\mathbf{b}\|^{2}\leq\min_{\mathbf{x}}\|\mathbf{A}\mathbf{x}-\mathbf{b}\|^{2}+\epsilon\|\mathbf{b}\|^{2}$
in time

$\tilde{O}\Big{(}(\mathrm{nnz}(\mathbf{A})+d^{2}+dk^{\omega-1})\log
1/\epsilon\Big{)}.$

#### Positive semidefinite systems.

An important application of our results is solving symmetric positive
semidefinite (psd) linear systems, which arise from data covariances, kernel
matrices, Hessian matrices, and more. While Theorem 1.1 already provides
improved time complexity guarantees for dense psd systems, a careful
adaptation of our algorithm allows it to exploit the sparsity structure of the
problem, as well as improving the $\mathrm{polylog}(n)$ factors in the bounds.
Similarly to prior stochastic psd linear system solvers [LS13], our complexity
bounds depend on the following notion of sparsity,
$\mathrm{nnz}^{*}(\mathbf{A})=n\cdot\max_{i}\\{\mathrm{nnz}(\mathbf{A}_{i,:})\\}\leq
n^{2}$, which is proportional to $\mathrm{nnz}(\mathbf{A})$ when the non-zeros
are distributed somewhat uniformly. For a psd matrix $\mathbf{A}$, we define
$\|x\|_{\mathbf{A}}=\sqrt{\mathbf{x}^{\top}\mathbf{A}\mathbf{x}}$.

######  Theorem 1.3 (PSD linear system, simplified Theorem 8.1).

Given a psd matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$ with $k$ large
eigenvalues, $\mathbf{b}\in\mathbb{R}^{n}$, and $\epsilon>0$, we can find
$\tilde{\mathbf{x}}$ such that
$\|\tilde{\mathbf{x}}-\mathbf{A}^{\dagger}\mathbf{b}\|_{\mathbf{A}}\leq\epsilon\|\mathbf{A}^{\dagger}\mathbf{b}\|_{\mathbf{A}}$
in time:

$\displaystyle\tilde{O}\Big{(}(\mathrm{nnz}^{*}(\mathbf{A})+nk^{\omega-1})\log
1/\epsilon\Big{)}.$

Our results naturally apply to regularized linear systems, such as kernel
ridge regression, i.e., $(\mathbf{K}+\lambda\mathbf{I})\mathbf{x}=\mathbf{b}$
for a psd matrix $\mathbf{K}$, and $l_{2}$-regularized least squares,
$\min_{\mathbf{x}}\|\mathbf{A}\mathbf{x}-\mathbf{b}\|^{2}+\lambda\|\mathbf{x}\|^{2}$
(as special cases of Theorem 1.3 and Theorem 1.2, respectively). In these
cases, we provide the first complexity bounds that scale with the number of
eigenvalues (of $\mathbf{K}$ or of $\mathbf{A}^{\top}\mathbf{A}$) larger than
the regularization parameter $\lambda$. This is in contrast to prior works
that directly exploit this explicit regularization structure [ACW17, FTU21],
where the time complexity scales with a notion of $\lambda$-effective
dimension which can be arbitrarily larger than our $k$. Thus, remarkably, we
achieve a direct runtime improvement over explicitly regularized linear system
solvers, even though our algorithms are entirely oblivious to the
regularization structure themselves (see Section 1.3 for discussion).

###  1.2 Our Techniques

A core inspiration for our algorithms is a broad family of stochastic
iterative methods for solving linear systems, which includes randomized
Kaczmarz [SV09], coordinate descent [LL10], sketch-and-project [GR15], and
many others [LS13]. Unlike CG-type methods which do a pass over the entire
matrix $\mathbf{A}$ in each iteration, these approaches randomly sub-sample a
small number of rows/columns of $\mathbf{A}$ (or construct a small sketch of
$\mathbf{A}$) and use this stochastic information to perform a sublinear
iterative update. There has been much success in characterizing the
convergence and complexity of stochastic methods in the ultra-small sample
regime (i.e., single row/column per iteration [SV09]), but the theoretical
analysis of the larger sample regime has proven much more elusive. While there
has been some progress in this direction over the last few years, it has come
with trade-offs: meaningful characterizations of the convergence behavior have
required relying on either constructing sophisticated and expensive sub-
sampling distributions (such as determinantal point processes [RK20, MDK20]),
or even more expensive dense sketching methods (such as the Gaussian sketch
[RN21, DR22]), rendering the entire procedure computationally impractical.

In light of this context, when designing and analyzing our stochastic
iterative solvers, we encountered three main bottlenecks which had to be
overcome to establish any improvement over CG-based complexity guarantees for
solving dense linear systems:

  1. 1.

Obtaining a sharp convergence guarantee that depends on the spectral structure
of $\mathbf{A}$;

  2. 2.

Relying on a fast sub-sampling/sketching scheme;

  3. 3.

Ensuring fast computation of the iterative updates.

1\. Convergence guarantee. As a starting point in our analysis we used recent
works on analyzing the convergence rate of block coordinate descent algorithms
using Determinantal Point Processes (DPP) as the sub-sampling method. DPPs
form a family of discrete distributions over subsets of a finite domain, e.g.,
$\\{1,...,n\\}$, where the probabilities of each subset are proportional to a
determinant of a certain sub-matrix (the uniform distribution over spanning
trees of a graph is one example). Exploiting certain formulas over sums of
determinants and adjugate matrices, recent works [RK20, MDK20] have
characterized the convergence rate of block coordinate descent algorithms with
DPP and k-DPP (DPP with size $k$ subsets) sub-sampling in terms of the
elementary symmetric polynomials (ESP) of the singular values of matrix
$\mathbf{A}$. Of particular relevance to our results are the guarantees
obtained by Rodomanov and Kropotov [RK20]. However, their convergence bounds
by themselves are not sufficiently sharp to obtain any improvement over CG-
type methods, even if we entirely ignored the remaining two bottlenecks. Thus,
as a first step, we improve the convergence guarantees of [RK20] for k-DPP-
based block coordinate descent by applying theory of majorization, via Schur-
concavity of the ratios of ESPs, leading to a factor $k$ improvement in the
convergence complexity. In particular, we show that for matrices with $k$
large singular values, using $\tilde{O}(k)$-sized blocks, we need
$\tilde{O}(n/k)$ iterations to solve a linear system.

2\. Fast sub-sampling. There has been much recent work on efficiently sampling
from discrete distributions, including strongly Rayleigh measures [AGR16] such
as k-DPPs. Most recently, Anari, Liu and Vuong [ALV22] provided nearly-tight
complexity analysis for a Markov chain Monte Carlo (MCMC) algorithm that
samples from a k-DPP (to within a small total variation distance). However,
the time complexity of their algorithm is still not sufficient to produce a
stochastic method that improves upon CG-type solvers for dense linear systems,
as it requires performing multiple singular value decompositions of $n\times
k$ matrices to produce a single k-DPP sample, at the cost of at least
$\tilde{O}(nk^{\omega-1})$ per iteration. To address this, in the next step of
our analysis, we show that the convergence guarantees of k-DPP sampling can be
retained (up to logarithmic factors) by a much simpler and cheaper sub-
sampling procedure, after preprocessing the linear system with a randomized
Hadamard transform [Tro11]. We show this by coupling a hypothetical run of the
MCMC algorithm from [ALV22] with our sub-sampling procedure, in such a way
that the k-DPP sample will (with high probability) be contained within our
sample.

3\. Fast updates. Standard block coordinate methods perform their iterations
by reducing the $n\times n$ linear system to a much smaller $k\times n$ linear
system, which is extremely under-determined, and then choose a solution with
the smallest distance from the current iterate. For $k=O(1)$, this is a very
simple and inexpensive update, but as $k$ grows, the cost of solving the
system exactly becomes prohibitively expensive. However, we can leverage the
fact that, unlike the original system, this new system is extremely
rectangular, which makes it much easier to construct an effective
preconditioner via sketching [RT08, MSM14, Woo14]. We prove that
$O(\log(n/k))$ (inner) iterations of preconditioned CG can be used to solve
this sub-problem to within sufficient accuracy so that the overall convergence
complexity of the stochastic method remains unchanged. This improves the cost
of the stochastic updates from $O(nk^{\omega-1})$ to
$\tilde{O}(nk+k^{\omega})$. Since we showed that $\tilde{O}(n/k)$ updates are
needed to converge, this leads to the final complexity result in Theorem 1.1.

###  1.3 Background and Related Work

Linear system solvers. Extensive literature has been dedicated to solving
linear systems both via direct and iterative methods. The direct approaches
are based on the strategy of applying a (possibly implicit) representation of
$\mathbf{A}^{-1}$ to the vector $\mathbf{b}$. Examples of this approach
include Gaussian elimination, Cholesky factorization and QR decomposition. A
key breakthrough in the time complexity of direct methods came from Strassen,
who showed that matrix inversion is equivalent to matrix multiplication
[Str69]. This kickstarted a long line of works on fast $O(n^{\omega})$ matrix
multiplication [Pan84], including the Coppersmith-Winograd method [CW87]. This
is still an active area of research [Wil12, LG12, AW21, WXXZ23].

There has also been much work on accelerating iterative linear system solvers,
particularly those based on Krylov subspace methods, such as conjugate
gradient / Lanczos. Here, most of the efforts have been concentrated on
constructing effective preconditioners, for example via specialized efficient
Cholesky factorizations [KLP+16, KS16], or on improving their numerical
stability [PV21, Nie22]. However, these approaches have focused on sparse
and/or structured systems, e.g., those decomposable into a sum of simple
components, particularly Laplacians of undirected graphs [Vai89, ST14, KMP12],
and their extensions [CKK+18]. Other examples of structured systems with
efficient solvers include Hankel/Toeplitz matrices [KKM79, XXG12] and
circulant matrices [Gra06].

Preconditioning regularized systems. Linear systems with an explicit
regularization term form another class of structured systems that can benefit
from effective preconditioning. The main examples here are: the
$l_{2}$-regularized least squares task,
$\min_{\mathbf{x}}\|\mathbf{A}\mathbf{x}-\mathbf{b}\|^{2}+\lambda\|\mathbf{x}\|^{2}$;
kernel ridge regression, $(\mathbf{K}+\lambda\mathbf{I})\mathbf{x}=\mathbf{b}$
where $\mathbf{K}$ is a psd kernel matrix; and the regularized Newton step,
$(\mathbf{H}+\lambda\mathbf{I})\mathbf{x}=\mathbf{g}$, where $\mathbf{H}$ is
the Hessian matrix, which is symmetric but may not be psd. A long line of
works have used randomized sketching, sampling, and the Nyström method [EAM14,
RCR17, MM17] to construct effective low-rank preconditioners for these systems
[ACW17, FTU21], relying on the explicit regularizer to avoid power iteration.
To construct a good preconditioner for, say, the matrix
$\mathbf{K}+\lambda\mathbf{I}$, those algorithms construct a decomposition of
$\mathbf{K}$ with rank at least proportional to the so-called
$\lambda$-effective dimension:
$d_{\lambda}={\operatorname{tr}}(\mathbf{K}(\mathbf{K}+\lambda\mathbf{I})^{-1})\leq
n$ [EAM14]. The best known cost of this preconditioning step, even exploiting
the structural properties of psd matrices, is at least
$\tilde{O}(nd_{\lambda}^{\omega-1})$ [MM17]. Similarly as our parameter $k$,
the $\lambda$-effective dimension exploits the flat tail of the spectrum of
$\mathbf{K}+\lambda\mathbf{I}$. In fact, it is easy to see that such a matrix
may have no more than $k=2d_{\lambda}$ eigenvalues larger than twice its
smallest eigenvalue. However the bound cannot be reversed: for any $n$, $k$
and $\lambda$, one can construct an $n\times n$ psd matrix $\mathbf{K}$ such
that $\mathbf{K}+\lambda\mathbf{I}$ has only $k$ eigenvalues larger than twice
its smallest eigenvalue, and yet its $\lambda$-effective dimension is at least
$n/2$, in which case, the preconditioning cost becomes $\Theta(n^{\omega})$,
as large as the cost of directly solving the system. Thus, our algorithms not
only recover but also improve upon the time complexity of these specialized
solvers, even without explicit access to the regularization term.

Stochastic solvers. The literature on stochastic iterative methods for solving
linear systems has its origins in the Kaczmarz method [Kac37], which performs
iterative updates while cycling through the rows of the system and projecting
the current iterate onto the hyperplane formed by an individual row. A
randomized version of this method was introduced and analyzed by Strohmer and
Vershynin [SV09], followed up by closely related randomized coordinate descent
methods for psd linear systems and least squares by Leventhal and Lewis
[LL10]. This triggered much research into stochastic iterative solvers,
leading to more general classes of methods that sample blocks of rows instead
of a single row [NW13, NT14, MDK20], or use matrix sketching [RN21, DR22] via
the Sketch-and-Project framework [GR15]. Despite these efforts, it is still
not fully understood when stochastic solvers are faster than deterministic
(e.g., Krylov-based) iterative methods. One such result was obtained by Lee
and Sidford [LS13], who showed that the convergence rate of an accelerated
stochastic coordinate descent solver exhibits better condition number
dependence than conjugate gradient for psd linear systems. However, this still
leads to poor performance for matrices that have an ill-conditioned top-$k$
part of the spectrum, whereas our algorithms avoid this by using blocks of
coordinates (rows or columns) as opposed to one coordinate at a time.

Matrix sketching and sub-sampling. Our techniques are closely related to
matrix sketching and sampling techniques commonly used in randomized linear
algebra [HMT11, Woo14, DM16, MT20, MDM+23], which is an area that has led to
improved randomized approximation algorithms for tasks such as least squares
regression [Sar06, RT08, CW13], low-rank approximation [CEM+15, CMM17, LW20],
$l_{p}$ regression [MM13, CP15, CD21, WW22], linear programming [CLS21] and
more [SWZ19, JKL+20]. Some of the key matrix sketching techniques we use in
this work are subspace embeddings [CW13, NN13, MM13, Coh16, CDDR23], which can
be used to construct effective preconditioners for highly over- or under-
determined linear systems (e.g., least squares). In particular, our fast least
squares solver is an accelerated version of the sketch-to-precondition
algorithm by Rokhlin and Tygert [RT08], which, using current subspace
embedding techniques, can solve this task in
$\tilde{O}(\mathrm{nnz}(\mathbf{A})+d^{\omega})$ time by sketching the tall
$n\times d$ matrix $\mathbf{A}$ down to an $\tilde{O}(d)\times d$ size, and
then performing a QR/SVD decomposition to produce a preconditioner. When the
matrix $\mathbf{A}$ has only $k$ large singular values, then this can be
improved to $\tilde{O}(\mathrm{nnz}(\mathbf{A})+d^{2}k^{\omega-2})$, by using
randomized low-rank SVD with power iteration [HMT11]. We instead use our
stochastic solver to avoid constructing an explicit preconditioner altogether
(see Table 1 for comparison). Finally, as mentioned earlier, many randomized
low-rank approximation methods mentioned above, such as [CEM+15, CMM17], avoid
power iteration and run in linear time. However, these approaches provide a
Frobenius norm error guarantee (as opposed to an approximation of all top $k$
singular values), which does not translate into a strong preconditioner for a
linear system.

One of the key matrix sub-sampling techniques we use (albeit mainly for our
theoretical analysis, rather than as an algorithmic routine) is called
determinantal point processes (DPP), also known as volume sampling. This is a
family of discrete distributions which has its origins in physics [Mac75] and
graph theory [Gue83]. More recently, DPPs have been used in randomized linear
algebra [DM21] to construct optimal coresets for the column subset selection
problem [DRVW06, GS12, DKM20] and least squares regression [DCMW19], as well
as for discrete and stochastic optimization [NST19, DLM20, RK20, MDK20]. Early
algorithmic approaches to DPP sampling required an exact SVD of the input
matrix [HKPV06, KT12], but a more recent line of works showed that this can be
avoided [Der19, DCV19], culminating in fast Markov chain Monte Carlo sampling
algorithms for DPPs [AGR16, AD20, ALV22]. Still, these methods lead to a
substantial overhead in time complexity for our algorithms, which is why we do
not use them explicitly except for certain special cases (namely, psd linear
systems).

##  2 Preliminaries

#### Problem Setup and Notation.

Throughout most of the paper, we will focus on a linear system
$\mathbf{A}\mathbf{x}=\mathbf{b}$ where $\mathbf{A}\in\mathbb{R}^{m\times n}$
and $\mathbf{b}\in\mathbb{R}^{m}$. Let $\mathrm{rank}(\mathbf{A})=r$ and
$\\{\sigma_{i}\\}_{i=1}^{r}$ be the decreasing singular values of
$\mathbf{A}$. There are two cases we consider. In the first case where
$\mathbf{b}\in\mathrm{range}(\mathbf{A})$, the linear system is consistent and
has one unique min-length solution
$\mathbf{x}^{*}=\mathbf{A}^{\dagger}\mathbf{b}$. In the second case where
$\mathbf{b}\notin\mathrm{range}(\mathbf{A})$, the linear system is
inconsistent and there is no exact solution, and we denote
$\mathbf{x}^{*}=\mathbf{A}^{\dagger}\mathbf{b}$ as the solution of the least
squares problem $\min_{\mathbf{x}}\|\mathbf{A}\mathbf{x}-\mathbf{b}\|$. As a
special case, we also consider solving a symmetric positive semidefinite
definite (psd) linear system $\mathbf{A}\mathbf{x}=\mathbf{b}$ where
$\mathbf{A}\in\mathcal{S}_{n}$, $\mathbf{b}\in\mathbb{R}^{n}$, and we use
$\mathcal{S}_{n}$ to denote the psd cone.

For vector $\mathbf{v}$, we use $\|\mathbf{v}\|$ to denote its Euclidean norm,
and for a psd matrix $\mathbf{M}$, we denote
$\|\mathbf{v}\|_{\mathbf{M}}:=\sqrt{\mathbf{v}^{\top}\mathbf{M}\mathbf{v}}$.
For matrix $\mathbf{A}$, we denote $\|\mathbf{A}\|_{2}$ and
$\|\mathbf{A}\|_{F}$ as its operator norm and Frobenius norm respectively, and
denote $\kappa(\mathbf{A})$ as its condition number. We also define the
averaged condition number of $\mathbf{A}$ as
$\bar{\kappa}(\mathbf{A}):=\frac{\|\mathbf{A}\|_{F}\|\mathbf{A}^{\dagger}\|_{2}}{\sqrt{\mathrm{rank}(\mathbf{A})}}\leq\kappa(\mathbf{A})$,
which arises in the analysis of our algorithms. For $1\leq
k\leq\mathrm{rank}(\mathbf{A})$, we denote the best rank-$k$ approximation of
$\mathbf{A}$ as
$\mathbf{A}_{k}:=\operatorname*{arg\,min}_{\mathrm{rank}(\hat{\mathbf{A}})=k}\|\mathbf{A}-\hat{\mathbf{A}}\|_{2}$,
and denote $\bar{\kappa}_{k}:=\bar{\kappa}(\mathbf{A}-\mathbf{A}_{k})$ (i.e.,
the averaged condition number of the tail of the spectrum). Let
$\\{\sigma_{i}\\}_{i=1}^{r}$ be the singular values of $\mathbf{A}$ in
decreasing order, if we assume that there exists $k$ such that
$\sigma_{i}/\sigma_{r}=O(1)$ holds for all $i\geq k$, then we have
$\bar{\kappa}_{k}=(\frac{1}{r-k}\sum_{j>k}^{r}\sigma_{j}^{2}/\sigma_{r}^{2})^{1/2}=O(1)$.
For simplicity of presentation, throughout the paper we assume that
$\bar{\kappa}_{k}<m^{4}$, which is a rather loose condition. For two $n\times
n$ psd matrices $\mathbf{A}$ and $\mathbf{B}$, we say
$\mathbf{A}\approx_{1+\epsilon}\mathbf{B}$ when
$\frac{1}{1+\epsilon}\mathbf{A}\preceq\mathbf{B}\preceq(1+\epsilon)\mathbf{A}$
where $\preceq$ denotes the matrix ordering.

For $1\leq\tau\leq m$, we use $\tbinom{[m]}{\tau}$ to denote the set of all
the $\tau$-element subsets of $[m]=\\{1,2,\ldots,m\\}$. For each
$S\in\tbinom{[m]}{\tau}$, we use $\mathbf{I}_{S}\in\mathbb{R}^{\tau\times m}$
to denote the matrix obtained from $m\times m$ identity matrix $\mathbf{I}$ by
only retaining the rows whose indices are in $S$. For matrix
$\mathbf{L}\in\mathbb{R}^{m\times m}$ and subset $S\in\tbinom{[m]}{\tau}$, we
use
$\mathbf{L}_{S,S}=\mathbf{I}_{S}\mathbf{L}\mathbf{I}_{S}^{\top}\in\mathbb{R}^{\tau\times\tau}$
to denote the principal submatrix located at the intersection of the rows and
columns with indices from $S$, we also denote
$\mathbf{L}_{:,S}=\mathbf{L}\mathbf{I}_{S}^{\top}$ and
$\mathbf{L}_{S,:}=\mathbf{I}_{S}\mathbf{L}$ respectively. For a square
$n\times n$ matrix $\mathbf{L}$, we use $\mathrm{adj}(\mathbf{L})$ to denote
its adjugate matrix, which is an $n\times n$ matrix such that
$\mathrm{adj}(\mathbf{L})_{ij}=(-1)^{i+j}\mathrm{det}(\mathbf{L}_{[n]\backslash
j,[n]\backslash i})$.

#### Elementary symmetric polynomials.

We next introduce some concepts that are needed for the convergence analysis
of our algorithms, starting with elementary symmetric polynomials (ESP) and
their properties.

######  Definition 2.1 (Elementary symmetric polynomial).

For a given vector
$\lambda=(\lambda_{1},\ldots,\lambda_{m})\in\mathbb{R}^{m}$, we define its
$\ell^{th}$ elementary symmetric polynomial as

$\displaystyle\sigma_{\ell}(\lambda):=\sum_{S\in\tbinom{[m]}{\ell}}\prod_{i\in
S}\lambda_{i}.$

######  Lemma 2.1 (Sum of principal minors).

Let $\mathbf{L}\in\mathbb{R}^{m\times m}$ be a symmetric matrix with
eigenvalues $\lambda=(\lambda_{1},\ldots,\lambda_{m})$ where
$\lambda_{1}\geq\cdots\geq\lambda_{m}$. Let $1\leq\ell\leq m$, then we have

$\displaystyle\sum_{S\in\tbinom{[m]}{\ell}}\mathrm{det}(\mathbf{L}_{S,S})=\sigma_{\ell}(\lambda).$

#### Determinantal point processes.

Given a PSD matrix with dimension $m$, a determinantal point process
($\mathrm{DPP}$) is a probability measure defined over $2^{[m]}$, i.e., all
the subsets of $[m]$. There is a special class of $\mathrm{DPP}$s where we
restrict the size of the subsets to be fixed, we call this a fixed-size
$\mathrm{DPP}$ or $k$-DPP, where $k$ is the size of the subsets. Formally we
have the following definitions.

######  Definition 2.2 ($k$-DPP).

Given a PSD matrix $\mathbf{L}\in\mathbb{R}^{m\times m}$, a $k$-$\mathrm{DPP}$
is a distribution on $\tbinom{[m]}{k}$ such that
$\Pr\\{S\\}\propto\mathrm{det}(\mathbf{L}_{S,S})$ where $S$ satisfies $|S|=k$.
We denote it as $k$-$\mathrm{DPP}(\mathbf{L})$.

An important special case of $k$-DPPs, used in some of the DPP sampling
algorithms, is known as Projection DPPs.

######  Definition 2.3 (Projection DPP).

Given a PSD matrix $\mathbf{L}\in\mathbb{R}^{m\times m}$, a $k$-$\mathrm{DPP}$
is called a Projection $\mathrm{DPP}$ iff $k=\mathrm{rank}(\mathbf{L})$. We
denote it as $\mathrm{P}$-$\mathrm{DPP}(\mathbf{L})$.

When the PSD matrix $\mathbf{L}$ is clear in the text, we omit $\mathbf{L}$
and denote the $k$-DPP and projection DPP as $k$-DPP and P-DPP respectively.
We have the following lemma which states that a $k$-DPP can be decomposed into
a mixture of projection DPPs. This result has been used to develop some of the
first DPP sampling algorithms, which require an eigendecomposition of
$\mathbf{L}$.

######  Lemma 2.2 (Algorithm 8 and Theorem 5.2 in [KT12]).

Given $\mathbf{L}\in\mathbb{R}^{m\times m}$ with
$\mathrm{rank}(\mathbf{L})=r$, let
$\mathbf{L}=\sum_{i}\lambda_{i}\mathbf{u}_{i}\mathbf{u}_{i}^{\top}$ be the
eigendecomposition of $\mathbf{L}$ where
$\lambda_{1}\geq\cdots\geq\lambda_{m}$. Suppose that we sample
$\gamma=(\gamma_{1},\ldots,\gamma_{r})$ where
$\gamma_{i}\sim\mathrm{Bernoulli}(\frac{\lambda_{i}}{\lambda_{i}+1})$ are
independent, and then sample
$S\sim\mathrm{P}$-$\mathrm{DPP}(\sum_{i}\gamma_{i}\mathbf{u}_{i}\mathbf{u}_{i}^{\top})$.
Then, conditioned on the event that $\sum_{i}\gamma_{i}=k$, we have
$S\overset{d}{=}S_{\mathrm{DPP}}$, where $S_{\mathrm{DPP}}\sim
k$-$\mathrm{DPP}(\mathbf{L})$.

#### Markov chain $k$-DPP sampling.

Recent works have shown that we can approximately sample from a $k$-DPP
without an eigendecomposition of the $\mathbf{L}$ matrix, by relying on Markov
chain Monte Carlo (MCMC) techniques, which apply more generally to all
strongly Rayleigh measures [ALV22].

######  Lemma 2.3 (Corollary 7 in [ALV22]).

Given an $m\times m$ PSD matrix $\mathbf{L}$, there is an algorithm that
outputs $s$ independent approximate samples from
$k$-$\mathrm{DPP}(\mathbf{L})$ in time $\tilde{O}(mk^{\omega-1}+sk^{\omega})$.

The above result is based on the following down-up walk sampling algorithm
(Algorithm 1). Given a discrete measure $\mu$ over $\tbinom{[m]}{\tau}$, this
procedure will converge to $\mu$ as its stationary distribution.

Algorithm 1 Down-up walk on the complement distribution for a discrete measure
$\mu$.

1:for $i=0,1,2,\ldots$ do

2: From all $t$-sized supersets of $S_{i}$, select one uniformly at random and
name it $T_{i}$;

3: Select among $k$-sized subsets of $T_{i}$ a random set $S_{i+1}$ with
$\Pr\\{S_{i+1}\\}\propto\mu(S_{i})$.

4:end for

#### Total variation distance.

In the analysis of sampling a $k$-DPP sample, we will need the following
definitions and properties of total variation distance.

######  Definition 2.4 (Total variation distance).

Let $\mu,\nu$ be probability measures defined on a measurable space
$(\Omega,\mathcal{F})$, the total variation distance between $\mu$ and $\nu$
is defined as

$\displaystyle\delta(\mu,\nu)=\sup_{A\in\mathcal{F}}|\mu(A)-\nu(A)|.$

######  Definition 2.5 (Coupling, see Definition 4.12 in [VH14]).

Let $\mu,\nu$ be two probability measures defined on $\Omega_{1},\Omega_{2}$
respectively. We say a joint distribution of random variable $Z=(X,Y)$ over
$\Omega_{1}\times\Omega_{2}$ is a coupling of $\mu$ and $\nu$, if it has
marginal distributions $X\sim\mu$ and $Y\sim\nu$.

The following standard lemma states that the total variation distance between
$\mu$ and $\nu$ is the minimum probability that random variables $X\sim\mu$
and $Y\sim\nu$ do not coincide.

######  Lemma 2.4 (Total variation, see Example 4.14 in [VH14]).

Let $\mu,\nu$ be probability measures defined on some $(\Omega,\mathcal{F})$.
Then, we have $\inf_{(X,Y)}\Pr\\{X\neq Y\\}=\delta(\mu,\nu)$, where the
infimum is taken over all couplings of $\mu$ and $\nu$.

#### Randomized Hadamard transform.

We use the randomized Hadamard transform (RHT) as a preprocessing step in our
algorithms, to ensure that a simple block coordinate sampling procedure works
as well as a determinantal point process.

######  Definition 2.6.

The Hadamard matrix $\mathbf{H}_{2^{k}}$ of dimension $2^{k}\times 2^{k}$ is
obtained using the following recurrence relation

$\displaystyle\mathbf{H}_{0}$
$\displaystyle=[1],\qquad\mathbf{H}_{2n}=\begin{bmatrix}\mathbf{H}_{n}&\mathbf{H}_{n}\\\
\mathbf{H}_{n}&-\mathbf{H}_{n}\end{bmatrix}.$

Throughout the paper, we will drop the subscript from $\mathbf{H}$, when the
dimensions are clear.

######  Definition 2.7.

The randomized Hadamard transform (RHT) of an $m\times n$ matrix $\mathbf{A}$
is the product $\frac{1}{\sqrt{m}}\mathbf{H}\mathbf{D}\mathbf{A}$, where
$\mathbf{D}$ is a random $m\times m$ diagonal matrix whose entries are
independent random $\pm 1$ signs. By padding $\mathbf{A}$ with zero rows if
necessary, we may assume that $n$ is a power of 2.

In the analysis of using randomized Hadamard transform, one key technique we
use is the tail bound for a Rademacher random vector, which is defined as
$\mathbf{d}=(d_{1},\ldots,d_{n})$ where each $d_{i}$ is an independent random
sign variable that equals $1$ with probability $1/2$ and $-1$ otherwise. We
have the following tail bound lemma.

######  Lemma 2.5 (Rademacher tail bound, e.g., [Tro11]).

Let $f$ be a convex function on vectors with Lipschitz bound
$|f(\mathbf{x})-f(\mathbf{y})|\leq L\|\mathbf{x}-\mathbf{y}\|$. Let
$\mathbf{d}$ be a Rademacher vector. Then, for all $t\geq 0$,

$\displaystyle\Pr\left\\{f(\mathbf{d})\geq{\mathbb{E}}[f(\mathbf{d})]+Lt\right\\}\leq
e^{-t^{2}/8}.$

#### Conjugate gradient.

We have the following algorithm and convergence guarantee of using conjugate
gradient method to solve positive semidefinite linear system
$\mathbf{A}\mathbf{A}^{\top}\mathbf{u}=\mathbf{b}$, where we are only given
$\mathbf{A}$ and $\mathbf{b}$ without knowing $\mathbf{A}\mathbf{A}^{\top}$.

######  Lemma 2.6 (Convergence of CG, [Hay18]).

Consider solving positive semidefinite linear system
$\mathbf{A}\mathbf{A}^{\top}\mathbf{u}=\mathbf{b}$ where
$\mathbf{A}\in\mathbb{R}^{\tau\times n},\mathbf{b}\in\mathbb{R}^{\tau}$ are
given. Let $\kappa:=\kappa(\mathbf{A})$ be the condition number and
$\mathbf{u}^{*}=(\mathbf{A}\mathbf{A}^{\top})^{\dagger}\mathbf{b}$ be the min-
length solution. Then, running conjugate gradient by starting from
$\mathbf{u}_{0}=\mathbf{0}$ will yield the following convergence result:

$\displaystyle\|\mathbf{u}_{s}-\mathbf{u}^{*}\|_{\mathbf{A}\mathbf{A}^{\top}}\leq
2\left(\frac{\kappa-1}{\kappa+1}\right)^{s}\|\mathbf{u}^{*}\|_{\mathbf{A}\mathbf{A}^{\top}}.$

We denote the output of conjugate gradient as
$\mathbf{u}_{s}=\mathrm{CG}(\mathbf{A},\mathbf{b},s)$.

#### Subspace embeddings.

To construct the preconditioner in our inner iteration for solving the linear
system, as well as for obtaining an improved input sparsity time algorithm for
least squares, we use the following result which guarantees a
$(1+\epsilon)$-subspace embedding with an optimal embedding dimension
[CDDR23]. We slightly adapt their result for our setting, so that it holds
with a high enough probability to apply a union bound across all iterations of
our algorithms.

######  Lemma 2.7 (Adapted from Theorem 1.4 in [CDDR23]).

Given $\mathbf{A}\in\mathbb{R}^{n\times d},\epsilon<1/2$, $\delta<1/2$, in
time
$O(\mathrm{nnz}(\mathbf{A})\log(d/\delta)/\epsilon+d^{2}\log^{4}(d/\delta)/\epsilon^{6})$
we can compute $\mathbf{\Phi}\mathbf{A}$ where
$\mathbf{\Phi}\in\mathbb{R}^{\phi\times n}$ is an embedding matrix with
$\phi=O((d+\log(1/\delta))/\epsilon^{2})$, such that with probability
$1-\delta$ we have

$\displaystyle\forall\mathbf{x}\in\mathbb{R}^{d},~{}~{}\frac{1}{1+\epsilon}\|\mathbf{A}\mathbf{x}\|\leq\|\mathbf{\Phi}\mathbf{A}\mathbf{x}\|\leq(1+\epsilon)\|\mathbf{A}\mathbf{x}\|.$

##  3 Main Algorithms

In this section we state our main algorithms and main results for solving a
dense linear system $\mathbf{A}\mathbf{x}=\mathbf{b}$ where
$\mathbf{A}\in\mathbb{R}^{m\times n}$. For the sake of completeness, and
because they provide somewhat different guarantees in somewhat different
regimes, we propose two algorithms.

The first algorithm (Algorithm 2) can be viewed as a Kaczmarz-type stochastic
solver, in that it sub-samples the rows of $\mathbf{A}$ and projects its
iterate onto the hyperplane defined by this smaller linear system. As other
Kaczmarz-type methods, this algorithm converges to the optimum solution of a
consistent linear system, and its convergence rate is naturally measured in
the Euclidean distance from the optimum, obtaining a guarantee of the form
$\|\tilde{\mathbf{x}}-\mathbf{x}^{*}\|\leq\epsilon\|\mathbf{x}^{*}\|$.

The second algorithm (Algorithm 3) can be viewed more as a coordinate descent-
type method, which sub-samples the columns of $\mathbf{A}$. This algorithm
actually solves an implicit psd linear system
$\mathbf{A}^{\top}\mathbf{A}\mathbf{x}=\mathbf{c}$, but by setting
$\mathbf{c}=\mathbf{A}^{\top}\mathbf{b}$, we can use it to solve the
$\mathbf{A}\mathbf{x}=\mathbf{b}$ system via the normal equations. This
approach has the advantage that when the original linear system is
inconsistent, then the procedure will converge to the least squares solution
$\mathbf{x}^{*}=\min_{\mathbf{x}}\|\mathbf{A}\mathbf{x}-\mathbf{b}\|=\mathbf{A}^{\dagger}\mathbf{b}$.
Another difference from the Kaczmarz-type algorithm is that the convergence
guarantees of Algorithm 3 are expressed as
$\|\mathbf{A}(\tilde{\mathbf{x}}-\mathbf{x}^{*})\|\leq\epsilon\|\mathbf{A}\mathbf{x}^{*}\|$.

In order to ensure fast convergence with a cheap sub-sampling scheme, we first
transform the linear system to an equivalent form by using the randomized
Hadamard transform (RHT). The transformed system enjoys the property that the
corresponding determinantal point process is nearly isotropic, which makes it
easier to sample, see Section 6. Note that our two algorithms use slightly
different preprocessing steps: the Kaczmarz-type method (Algorithm 2) applies
the RHT from the left,
$\mathbf{A}\leftarrow\frac{1}{\sqrt{m}}\mathbf{H}\mathbf{D}\mathbf{A}$,
because it uses row sampling; whereas the coordinate descent-type method
(Algorithm 3) applies the RHT from the right,
$\mathbf{A}\leftarrow\frac{1}{\sqrt{n}}\mathbf{A}\mathbf{D}\mathbf{H}$,
because it uses column sampling. In either case, this transformation can be
done in $\tilde{O}(mn)$ time, and since the RHT is an orthogonal matrix, we
can easily recover the solution to the original system.

Both algorithms use an inner solver to approximately solve the small linear
system obtained via sub-sampling (the projection step). We do this by
constructing a preconditioner for the smaller system and use CG as an inner
solver. The construction of the preconditioner is based on a subspace
embedding sketch, leveraging the fact that the small system is highly
rectangular, see Algorithm 6 and Lemma 2.7.

###  3.1 Fast Kaczmarz-type Solver

Next we state our first main algorithm and convergence result which solves a
consistent linear system. Recall that we define
$\bar{\kappa}(\mathbf{A})=\|\mathbf{A}\|_{F}\|\mathbf{A}^{\dagger}\|_{2}/\sqrt{r}\leq\kappa(\mathbf{A})$
as the averaged condition number for matrix $\mathbf{A}$, and define
$\bar{\kappa}_{k}(\mathbf{A})=\bar{\kappa}(\mathbf{A}-\mathbf{A}_{k})$ for
$k<r$. Denote $\\{\sigma_{i}\\}_{i=1}^{r}$ as the singular values of
$\mathbf{A}$ in decreasing order, then we have
$\bar{\kappa}_{k}(\mathbf{A})=(\frac{1}{r-k}\sum_{j>k}^{r}\sigma_{j}^{2}/\sigma_{r}^{2})^{1/2}\leq\sigma_{k+1}/\sigma_{r}$.
Moreover, notice that the statement “with probability $0.99$” depends only on
ensuring that the preprocessing construction
$\frac{1}{\sqrt{m}}\mathbf{H}\mathbf{D}\mathbf{A}$ is successful, and comes
from Lemma 6.4 by choosing $\delta=0.01$.

Algorithm 2 Fast Kaczmarz-type solver for consistent linear systems.

1:Input: matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$, vector
$\mathbf{b}\in\mathbb{R}^{m}$, sample size $\tau$, iterate $\mathbf{x}_{0}$,
inner iteration $s_{\max}$, outer iteration $t_{\max}$;

2:Compute
$\mathbf{A}\leftarrow\frac{1}{\sqrt{m}}\mathbf{H}\mathbf{D}\mathbf{A}$ and
$\mathbf{b}\leftarrow\frac{1}{\sqrt{m}}\mathbf{H}\mathbf{D}\mathbf{b}$;
$\triangleright$ Takes $O(mn\log m)$ time.

3:for $t=0,1,\ldots,t_{\max}-1$ do

4: Generate $\tilde{S}$ by calling Algorithm 5; $\triangleright$ Takes
$O(\tau)$ time.

5: Compute $\tilde{\mathbf{A}}\leftarrow\mathbf{I}_{\tilde{S}}\mathbf{A}$;
$\triangleright$ Subsampling: takes $O(n\tau)$ time.

6: Compute
$\tilde{\mathbf{b}}\leftarrow\mathbf{I}_{\tilde{S}}\mathbf{A}\mathbf{x}_{t}-\mathbf{I}_{\tilde{S}}\mathbf{b}$;
$\triangleright$ Takes $O(n\tau)$ time.

7: Generate preconditioner $\mathbf{M}$ by calling Algorithm 6 with
$\tilde{\mathbf{A}}$; $\triangleright$ Takes $\tilde{O}(n\tau+\tau^{\omega})$
time.

8: Compute
$\mathbf{u}_{s_{\max}}\leftarrow\mathrm{CG}(\mathbf{M}^{\top}\tilde{\mathbf{A}},\mathbf{M}^{\top}\tilde{\mathbf{b}},s_{\max})$
with Lemma 2.6; $\triangleright$ Takes $O(n\tau s_{\max})$ time.

9: Compute
$\mathbf{w}_{t}\leftarrow\tilde{\mathbf{A}}^{\top}\mathbf{M}\mathbf{u}_{s_{\max}}$;
$\triangleright$ Takes $O(n\tau)$ time.

10: Update $\mathbf{x}_{t+1}\leftarrow\mathbf{x}_{t}-\mathbf{w}_{t}$;
$\triangleright$ Takes $O(n)$ time.

11:end for

12:return $\tilde{\mathbf{x}}=\mathbf{x}_{t_{\max}}$; $\triangleright$ Solves
$\mathbf{A}\mathbf{x}=\mathbf{b}$.

######  Theorem 3.1 (Fast Kaczmarz-type solver).

Given matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$ with
$\mathrm{rank}(\mathbf{A})=r$ and $\mathbf{b}\in\mathbb{R}^{m}$, let
$\\{\sigma_{i}\\}_{i=1}^{r}$ be the singular values of $\mathbf{A}$ in
decreasing order. For $\log m\leq k<r$, define
$\bar{\kappa}_{k}:=(\frac{1}{r-k}\sum_{j>k}^{r}\sigma_{j}^{2}/\sigma_{r}^{2})^{1/2}$
and assume for simplicity that $\bar{\kappa}_{k}\leq m^{4}$. Assume the linear
system $\mathbf{A}\mathbf{x}=\mathbf{b}$ is consistent (i.e.,
$\mathbf{b}\in\mathrm{range}(\mathbf{A})$) and let
$\mathbf{x}^{*}=\mathbf{A}^{\dagger}\mathbf{b}$ be its min-length solution.
Assume $\mathbf{x}_{0}-\mathbf{x}^{*}$ is not in the nullspace of
$\mathbf{A}$. Then conditioned on an event that holds with probability $0.99$
(and only depends on $\frac{1}{\sqrt{m}}\mathbf{H}\mathbf{D}\mathbf{A}$), for
$C=O(1)$, running Algorithm 2 with choice $\tau\geq Ck\log^{3}m$ and
$s_{\max}=O(\log(r\bar{\kappa}_{k}/\tau))$ will yield the following
convergence result:

$\displaystyle{\mathbb{E}}\|\mathbf{x}_{t}-\mathbf{x}^{*}\|^{2}\leq\left(1-\frac{\tau/r}{10\bar{\kappa}_{k}^{2}\log^{3}m}\right)^{t}\cdot\|\mathbf{x}_{0}-\mathbf{x}^{*}\|^{2}.$

By further choosing $\tau=O(k\cdot\log^{3}m),\mathbf{x}_{0}=\mathbf{0}$ and
$t_{\max}=O(r\bar{\kappa}_{k}^{2}/k\cdot\log(1/\epsilon))$ for given
$\epsilon>0$, Algorithm 2 outputs $\tilde{\mathbf{x}}$ such that
$\|\tilde{\mathbf{x}}-\mathbf{x}^{*}\|^{2}\leq\epsilon\|\mathbf{x}^{*}\|^{2}$
holds with probability $0.98$ in time:

$\displaystyle O\Big{(}mn\log
m+(nr\log^{4}m+rk^{\omega-1}\log^{3\omega}m)\cdot\bar{\kappa}_{k}^{2}\log
1/\epsilon\Big{)}.$

###  3.2 Fast Coordinate Descent-type Solver

Next, we present our second main algorithm (fast coordinate descent-type
solver), and a corresponding convergence result, which solves a linear system
that does not necessarily need to be consistent. We note that while this
algorithm can be used directly on tall least squares problems, we provide an
improved variant for this setting, which can leverage the sparsity of
$\mathbf{A}$ to further improve the runtime (see Section 10).

######  Theorem 3.2 (Fast coordinate descent-type solver).

Given matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$ with
$\mathrm{rank}(\mathbf{A})=r$ and $\mathbf{b}\in\mathbb{R}^{m}$, let
$\\{\sigma_{i}\\}_{i=1}^{r}$ be the singular values of $\mathbf{A}$ in
decreasing order. For $\log n\leq k<r$, define
$\bar{\kappa}_{k}:=(\frac{1}{r-k}\sum_{j>k}^{r}\sigma_{j}^{2}/\sigma_{r}^{2})^{1/2}$
and assume for simplicity that $\bar{\kappa}_{k}\leq m^{4}$. Let
$\mathbf{x}^{*}=\mathbf{A}^{\dagger}\mathbf{b}$ be the least squares solution
of the (possibly inconsistent) linear system
$\mathbf{A}\mathbf{x}=\mathbf{b}$. Assume $\mathbf{x}_{0}-\mathbf{x}^{*}$ is
not in the nullspace of $\mathbf{A}$. Then conditioned on an event that holds
with probability $0.99$ (and only depends on
$\frac{1}{\sqrt{n}}\mathbf{A}\mathbf{D}\mathbf{H}$), for $C=O(1)$, running
Algorithm 3 with choice $\mathbf{c}=\mathbf{A}^{\top}\mathbf{b},\tau\geq
Ck\log^{3}n$ and $s_{\max}=O(\log(r\bar{\kappa}_{k}/\tau))$ will yield the
following convergence result:

$\displaystyle{\mathbb{E}}\|\mathbf{A}(\mathbf{x}_{t}-\mathbf{x}^{*})\|^{2}\leq\left(1-\frac{\tau/r}{10\bar{\kappa}_{k}^{2}\log^{3}n}\right)^{t}\cdot\|\mathbf{A}(\mathbf{x}_{0}-\mathbf{x}^{*})\|^{2}.$

By further choosing $\tau=O(k\cdot\log^{3}n),\mathbf{x}_{0}=\mathbf{0}$ and
$t_{\max}=O(r\bar{\kappa}_{k}^{2}/k\cdot\log(1/\epsilon))$ for given
$\epsilon>0$, Algorithm 3 outputs $\tilde{\mathbf{x}}$ such that
$\|\mathbf{A}(\tilde{\mathbf{x}}-\mathbf{x}^{*})\|^{2}\leq\epsilon\|\mathbf{A}\mathbf{x}^{*}\|^{2}$
holds with probability $0.98$ in time:

$\displaystyle O\Big{(}mn\log
n+(mr\log^{4}n+rk^{\omega-1}\log^{3\omega}n)\cdot\bar{\kappa}_{k}^{2}\log
1/\epsilon\Big{)}.$

######  Remark 3.1.

By further assuming that the linear system $\mathbf{A}\mathbf{x}=\mathbf{b}$
is consistent (i.e., $\mathbf{b}\in\mathrm{range}(\mathbf{A})$), we know that
$\mathbf{x}^{*}=\mathbf{A}^{\dagger}\mathbf{b}$ is indeed a solution of the
linear system and $\mathbf{A}\mathbf{x}^{*}=\mathbf{b}$. Suppose we choose
$k\geq\log n$ be such that $\frac{\sigma_{i}}{\sigma_{r}}=O(1)$ holds for all
$i\geq k$, then we have $\bar{\kappa}_{k}=O(1)$. Then, Algorithm 3 outputs
$\tilde{\mathbf{x}}$ such that
$\|\mathbf{A}\tilde{\mathbf{x}}-\mathbf{b}\|^{2}\leq\epsilon\|\mathbf{b}\|^{2}$
with probability $0.98$ in time:

$\displaystyle\tilde{O}\Big{(}mn+(mr+rk^{\omega-1})\log 1/\epsilon\Big{)}.$

If we do not know such $k$ in advance, we can use the following doubling trick
to find the right $k$: we start from a guess $\tilde{k}=1$ and run Algorithm 3
to obtain an output $\tilde{\mathbf{x}}$ and check if we have
$\|\mathbf{A}\tilde{\mathbf{x}}-\mathbf{b}\|^{2}\leq\epsilon\|\mathbf{b}\|^{2}$.
If this equation holds, then we can stop with the current value of $\tilde{k}$
and the algorithm solves the linear system; if not, we can simply double the
value of $\tilde{k}$ and run Algorithm 3 again. Notice that $\tilde{k}$
satisfies the definition of $k$ as long as $\tilde{k}\geq k$, since
$\sigma_{i}/\sigma_{r}=O(1)$ also holds for all $i\geq\tilde{k}$. In this way,
we need to run Algorithm 3 at most $\log r$ times (since $k\leq r$), where
each run costs $\tilde{O}((mn+rk^{\omega-1})\log 1/(\epsilon\delta))$ in order
to guarantee that the failure probability is at most $\delta$. Thus the final
cost is $\tilde{O}((mn+rk^{\omega-1})\log 1/(\epsilon\delta))$ with
$\tilde{O}(\cdot)$ hiding $\mathrm{polylog}(n)$ factors.

Algorithm 3 Fast coordinate descent-type solver for linear systems.

1:Input: matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$, vector
$\mathbf{c}\in\mathbb{R}^{n}$, uniform sample size $\tau$, iterate
$\mathbf{x}_{0}$, inner iteration $s_{\max}$, outer iteration $t_{\max}$;

2:Compute
$\mathbf{A}\leftarrow\frac{1}{\sqrt{n}}\mathbf{A}\mathbf{D}\mathbf{H}$ and
$\mathbf{c}\leftarrow\frac{1}{\sqrt{n}}\mathbf{H}\mathbf{D}\mathbf{c}$;
$\triangleright$ Takes $O(mn\log n)$ time.

3:Compute $\mathbf{y}_{0}\leftarrow\mathbf{A}\mathbf{x}_{0}$; $\triangleright$
Invariant: $\mathbf{y}_{t}=\mathbf{A}\mathbf{x}_{t}$.

4:for $t=0,1,\ldots,t_{\max}-1$ do

5: Generate $\tilde{S}$ by calling Algorithm 5; $\triangleright$ Takes
$O(\tau)$ time.

6: Compute
$\tilde{\mathbf{A}}\leftarrow\mathbf{A}\mathbf{I}_{\tilde{S}}^{\top}$;
$\triangleright$ Subsampling: takes $O(m\tau)$ time.

7: Compute
$\tilde{\mathbf{c}}\leftarrow\tilde{\mathbf{A}}^{\top}\mathbf{y}_{t}-\mathbf{I}_{\tilde{S}}\mathbf{c}$;
$\triangleright$ Takes $O(m\tau)$ time.

8: Generate preconditioner $\mathbf{M}$ by calling Algorithm 6 with
$\tilde{\mathbf{A}}^{\top}$; $\triangleright$ Takes
$\tilde{O}(m\tau+\tau^{\omega})$ time.

9: Compute
$\mathbf{w}_{s_{\max}}\leftarrow\mathrm{CG}(\mathbf{M}^{\top}\tilde{\mathbf{A}}^{\top},\mathbf{M}^{\top}\tilde{\mathbf{c}},s_{\max})$
with Lemma 2.6; $\triangleright$ Takes $O(m\tau s_{\max})$ time.

10: Compute
$\mathbf{x}_{t+1}\leftarrow\mathbf{x}_{t}-\mathbf{I}_{\tilde{S}}^{\top}\mathbf{w}_{s_{\max}}$;
$\triangleright$ Takes $O(\tau)$ time.

11: Compute
$\mathbf{y}_{t+1}\leftarrow\mathbf{y}_{t}-\tilde{\mathbf{A}}\mathbf{w}_{s_{\max}}$;
$\triangleright$ Takes $O(m\tau)$ time.

12:end for

13:return
$\tilde{\mathbf{x}}=\frac{1}{\sqrt{n}}\mathbf{D}\mathbf{H}\mathbf{x}_{t_{\max}}$;
$\triangleright$ Solves $\mathbf{A}^{\top}\mathbf{A}\mathbf{x}=\mathbf{c}$.

##  4 Overview of the Analysis

For simplicity, in this section we only provide analysis and proof sketch for
Algorithm 2 and Theorem 3.1, which is about solving a consistent linear system
$\mathbf{A}\mathbf{x}=\mathbf{b}$ where $\mathbf{A}\in\mathbb{R}^{m\times
n},\mathbf{b}\in\mathbb{R}^{m}$. For the analysis and proof of Algorithm 3 and
Theorem 3.2, which is similar, see Section 9.

###  4.1 Expected Projection under DPP Sampling

Let $\mathbf{I}_{S}\in\mathbb{R}^{\tau\times m}$ be a sampling matrix with
$\tau\leq n$. Our analysis is based on an iterative method to solve the linear
system which we will refer to as _sketch-and-project_ , following [GR15] (it
is also known as the block Kaczmarz method):

$\displaystyle\mathbf{x}_{t+1}=\mathbf{x}_{t}-\mathbf{A}^{\top}\mathbf{I}_{S}^{\top}(\mathbf{I}_{S}\mathbf{A}\mathbf{A}^{\top}\mathbf{I}_{S}^{\top})^{\dagger}(\mathbf{I}_{S}\mathbf{A}\mathbf{x}_{t}-\mathbf{I}_{S}\mathbf{b})=\mathbf{x}_{t}-\underbrace{(\mathbf{I}_{S}\mathbf{A})^{\dagger}(\mathbf{I}_{S}\mathbf{A}\mathbf{x}_{t}-\mathbf{I}_{S}\mathbf{b})}_{\mathbf{w}_{t}^{*}}.$
(4.1)

Denote $\mathbf{x}^{*}=\mathbf{A}^{\dagger}\mathbf{b}$ as the min-length
solution and denote
$\mathbf{P}_{S}=\mathbf{A}^{\top}\mathbf{I}_{S}^{\top}(\mathbf{I}_{S}\mathbf{A}\mathbf{A}^{\top}\mathbf{I}_{S}^{\top})^{\dagger}\mathbf{I}_{S}\mathbf{A}$.
We have

$\displaystyle\mathbf{x}_{t+1}-\mathbf{x}^{*}=(\mathbf{I}-\mathbf{P}_{S})\mathbf{x}_{t}+\mathbf{A}^{\top}\mathbf{I}_{S}^{\top}(\mathbf{I}_{S}\mathbf{A}\mathbf{A}^{\top}\mathbf{I}_{S}^{\top})^{\dagger}\mathbf{I}_{S}\mathbf{b}-\mathbf{x}^{*}=(\mathbf{I}-\mathbf{P}_{S})(\mathbf{x}_{t}-\mathbf{x}^{*}).$

Notice that $\mathbf{P}_{S}$ is a projection matrix and thus
$\mathbf{P}_{S}^{2}=\mathbf{P}_{S}$. Since $\mathbf{x}_{0}-\mathbf{x}^{*}$ is
not in the nullspace of ${\mathbb{E}}[\mathbf{P}_{S}]$, we know
$\mathbf{x}_{t}-\mathbf{x}^{*}$ will always not be in the nullspace of
${\mathbb{E}}[\mathbf{P}_{S}].$222For the detailed reason see proof of Lemma
4.3. Thus we have

$\displaystyle{\mathbb{E}}\|\mathbf{x}_{t+1}-\mathbf{x}^{*}\|^{2}=$
$\displaystyle~{}{\mathbb{E}}[(\mathbf{x}_{t}-\mathbf{x}^{*})^{\top}(\mathbf{I}-\mathbf{P}_{S})^{2}(\mathbf{x}_{t}-\mathbf{x}^{*})]$
$\displaystyle=$
$\displaystyle~{}(\mathbf{x}_{t}-\mathbf{x}^{*})^{\top}{\mathbb{E}}[(\mathbf{I}-\mathbf{P}_{S})^{2}](\mathbf{x}_{t}-\mathbf{x}^{*})$
$\displaystyle=$
$\displaystyle~{}(\mathbf{x}_{t}-\mathbf{x}^{*})^{\top}(\mathbf{I}-{\mathbb{E}}[\mathbf{P}_{S}])(\mathbf{x}_{t}-\mathbf{x}^{*})$
$\displaystyle\leq$
$\displaystyle~{}(1-\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{S}]))\cdot\|\mathbf{x}_{t}-\mathbf{x}^{*}\|^{2}.$

In order to obtain the convergence rate of sketch-and-project, we need to
lower bound the quantity $\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{S}])$,
which turns out to be quite challenging for most sampling schemes. To address
this, following recent prior work [RK20], we will consider $S$ sampled
according to $k$-$\mathrm{DPP}(\mathbf{A}\mathbf{A}^{\top})$. [RK20] provide a
lower bound of this quantity by relating the expectation under a determinantal
point process to a formula for the sum of adjugates of sub-matrices, expressed
in terms of the elementary symmetric polynomials of the eigenvalues of
$\mathbf{A}\mathbf{A}^{\top}$. Directly applying their results to our setting
with $k^{\prime}=O(k)$, we can obtain a lower bound of order
$\Omega(\frac{1}{r})$, where $r$ is the rank of $\mathbf{A}$ and we hide the
dependence on $\bar{\kappa}_{k}$. In the following lemma, we refine the proof
by providing a sharper bound on the expression derived by [RK20], using Schur
concavity of the ratios of elementary symmetric polynomials. Namely, we show
that if we use a slightly larger DPP sample size $k^{\prime}=2k$, then
$\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}])=\Omega(\frac{k}{r})$,
an improvement by a factor $k$ which ends up being crucial in the complexity
analysis (see Section 5). Going forward, to simplify the notation we will
often use $k$ instead of $k^{\prime}$, since they are within a factor 2 of one
another.

######  Lemma 4.1 (Improved bound).

Given matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$ with rank $r$. Let
$S_{\mathrm{DPP}}\sim k^{\prime}$-$\mathrm{DPP}(\mathbf{A}\mathbf{A}^{\top})$
where $k^{\prime}<r$, let
$\mathbf{I}_{S}=\mathbf{I}_{S_{\mathrm{DPP}}}\in\mathbb{R}^{k^{\prime}\times
m}$ be the corresponding sampling matrix and define the projection matrix
$\mathbf{P}_{S_{\mathrm{DPP}}}:=\mathbf{A}^{\top}\mathbf{I}_{S}^{\top}(\mathbf{I}_{S}\mathbf{A}\mathbf{A}^{\top}\mathbf{I}_{S}^{\top})^{\dagger}\mathbf{I}_{S}\mathbf{A}$.
For any $k<k^{\prime}$ we have

$\displaystyle\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}])\geq\frac{k^{\prime}-k}{k^{\prime}-k-1+(r-k)\bar{\kappa}_{k}^{2}}.$
(4.2)

###  4.2 Replacing DPP with Uniform Sampling

Despite many recent improvements, sampling a $k$-DPP is still too expensive
for our problem. Currently the best known algorithm for sampling $k$-DPPs is
obtained by [ALV22], which gives an $\tilde{O}(mk^{\omega-1}+sk^{\omega})$
complexity to approximately obtain $s$ samples according to a $k$-DPP based on
an $m\times m$ PSD matrix, see Lemma 2.3. Since in our case we are dealing
with a rectangular non-PSD matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$, after
adapting their algorithm to this setting, the cost of sampling becomes
$\tilde{O}(mnk^{\omega-2}+snk^{\omega-1})$, see Section 6.4. However this is
still too expensive. Can we get rid of $\tilde{O}(mnk^{\omega-2})$? Notice
that we do not precisely need a $k$-DPP sample, but rather, it suffices to
produce a sample that _contains_ a $k$-DPP, since the latter can actually
guarantee a faster convergence rate, see Lemma 6.3. We propose Algorithm 5
which only uses _uniform sampling_ with replacement, and show that, after
preprocessing the linear system with the randomized Hadamard transform, this
algorithm will output a sample $\tilde{S}$ that contains a $k$-DPP sample with
high probability. The size of such a uniform sample is only
$O((k+\log\frac{1}{\delta})\log^{3}m)$, which is almost linear in $k$, see
Lemma 6.2 and Section 6.

By combining Lemma 4.1, Lemma 6.2 and Lemma 6.3, we have the following main
result for our sampling algorithm. Notice that we first apply randomized
Hadamard transform $\mathbf{Q}=\frac{1}{\sqrt{m}}\mathbf{H}\mathbf{D}$ to
matrix $\mathbf{A}$ as a preprocessing step, in order to guarantee that with
high probability the marginal probabilities of the corresponding $k$-DPP are
approximately isotropic (see Lemma 6.4). For the full proof of Lemma 4.2 see
Section 6.3.

######  Lemma 4.2 (Sampling).

Given matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$ with rank $r$, denote
$\mathbf{Q}=\frac{1}{\sqrt{m}}\mathbf{H}\mathbf{D}$ as the RHT matrix and
$\bar{\mathbf{A}}=\mathbf{Q}\mathbf{A}$ as the preprocessed matrix.
Conditioned on an event that holds with probability $0.99$ (and only depends
on RHT), for $k\geq\log m$, if we draw $\tau$ uniform samples from $[m]$
(according to Algorithm 5) such that $\tau\geq Ck\log^{3}m$ for some constant
$C=O(1)$, then the resulting set $\tilde{S}\subseteq[m]$ satisfies

$\displaystyle\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{\tilde{S}}])\geq\frac{\tau/r}{2.5\bar{\kappa}_{k}^{2}\log^{3}m}$

where
$\mathbf{P}_{\tilde{S}}:=\bar{\mathbf{A}}^{\top}\mathbf{I}_{\tilde{S}}^{\top}(\mathbf{I}_{\tilde{S}}\bar{\mathbf{A}}\bar{\mathbf{A}}^{\top}\mathbf{I}_{\tilde{S}}^{\top})^{\dagger}\mathbf{I}_{\tilde{S}}\bar{\mathbf{A}}$
is the projection matrix and $\mathbf{I}_{\tilde{S}}\in\mathbb{R}^{\tau\times
m}$ is the row sampling matrix according to $\tilde{S}$.

###  4.3 Approximate Sketch-and-Project

In the update rule of sketch-and-project (see Eq.(4.1)), the remaining problem
is that we need to compute
$\mathbf{w}_{t}^{*}:=(\mathbf{I}_{\tilde{S}}\bar{\mathbf{A}})^{\dagger}(\mathbf{I}_{\tilde{S}}\bar{\mathbf{A}}\mathbf{x}_{t}-\mathbf{I}_{\tilde{S}}\bar{\mathbf{b}})$
in each iteration, where $\bar{\mathbf{A}}=\mathbf{Q}\mathbf{A}$ and
$\bar{\mathbf{b}}=\mathbf{Q}\mathbf{b}$. Directly computing this quantity
takes $O(n\tau^{\omega-1})$. Instead of directly computing
$\mathbf{w}_{t}^{*}$, we show that it is possible to _approximate_ this
quantity in $\tilde{O}(n\tau+\tau^{\omega})$ time (see Section 7). Denote
$\tilde{\mathbf{A}}=\mathbf{I}_{\tilde{S}}\bar{\mathbf{A}}\in\mathbb{R}^{\tau\times
n}$ and
$\tilde{\mathbf{b}}_{t}=\mathbf{I}_{\tilde{S}}\bar{\mathbf{A}}\mathbf{x}_{t}-\mathbf{I}_{\tilde{S}}\bar{\mathbf{b}}$,
let $\mathbf{M}\in\mathbb{R}^{\tau\times\tau}$ be a preconditioner such that
$\tilde{\kappa}:=\kappa(\tilde{\mathbf{A}}^{\top}\mathbf{M})=O(1)$ holds with
high probability. We first show that we can transform the problem of computing
$\mathbf{w}_{t}^{*}=\tilde{\mathbf{A}}^{\dagger}\tilde{\mathbf{b}}_{t}$ to
solving the following preconditioned normal equation of the second kind:

$\displaystyle(\mathbf{M}^{\top}\tilde{\mathbf{A}}\tilde{\mathbf{A}}^{\top}\mathbf{M})\mathbf{u}=\mathbf{M}^{\top}\tilde{\mathbf{b}}_{t},~{}~{}~{}\mathbf{w}=\tilde{\mathbf{A}}^{\top}\mathbf{M}\mathbf{u}.$

Then we use CG to solve this positive semidefinite linear system over
$\mathbf{u}$. Let
$\mathbf{u}^{*}=\mathbf{M}^{-1}(\tilde{\mathbf{A}}\tilde{\mathbf{A}}^{\top})^{\dagger}\tilde{\mathbf{b}}_{t}$
and $\tilde{\mathbf{w}}^{*}=\tilde{\mathbf{A}}^{\top}\mathbf{M}\mathbf{u}^{*}$
be the optimal solution, we show that
$\tilde{\mathbf{w}}^{*}=\mathbf{w}_{t}^{*}$. By using Lemma 2.6 and doing the
change of variable back from $\mathbf{u}$ to $\mathbf{w}$, we show in Lemma
7.1 that after $s$ iterations,

$\displaystyle\|\mathbf{w}_{s}-\mathbf{w}_{t}^{*}\|\leq
2\left(\frac{\tilde{\kappa}-1}{\tilde{\kappa}+1}\right)^{s}\|\mathbf{w}^{*}_{t}\|.$

The only part left is to construct a preconditioner
$\mathbf{M}\in\mathbb{R}^{\tau\times\tau}$ such that $\tilde{\kappa}=O(1)$,
and we construct it by using fast oblivious subspace embedding matrix proposed
by [CDDR23]. According to Lemma 2.7, by setting $\phi=O(\tau+\log(1/\delta))$
we have that $\tilde{\kappa}=O(1)$ holds with probability $1-\delta$, since we
can set $\epsilon$ to be a small constant. Next we analyze the time
complexity: given $\tilde{\mathbf{A}}$, computing
$\mathbf{\Theta}\tilde{\mathbf{A}}^{\top}$ takes
$O(n\tau\log(\tau/\delta)+\tau^{2}\log^{4}(\tau/\delta))$, computing the SVD
of $\mathbf{\Theta}\tilde{\mathbf{A}}^{\top}$ takes $O(\phi\tau^{\omega-1})$,
and computing $\mathbf{M}$ takes $O(\tau^{\omega})$. Thus we conclude that the
total time complexity is

$\displaystyle
O(n\tau\log(\tau/\delta)+\tau^{2}\log^{4}(\tau/\delta)+\tau^{\omega-1}\phi)=O(n\tau\log(\tau/\delta)+\tau^{\omega}+\tau^{\omega-1}\log(1/\delta))=O(n\tau\log(\tau/\delta)+\tau^{\omega}).$

With all the above analysis, we have the following convergence lemma for
approximate sketch-and-project. For the proof of Lemma 4.3 see Section 7.2.

######  Lemma 4.3 (Convergence of approximate sketch-and-project).

Given matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$ with rank $r$, denote
$\mathbf{Q}=\frac{1}{\sqrt{m}}\mathbf{H}\mathbf{D}$ as the RHT matrix and
$\bar{\mathbf{A}}=\mathbf{Q}\mathbf{A}$ as the preprocessed matrix. Let
$\mathbf{I}_{\tilde{S}}\in\mathbb{R}^{\tau\times m}$ be the sampling matrix
and
$\mathbf{P}_{\tilde{S}}:=\bar{\mathbf{A}}^{\top}\mathbf{I}_{\tilde{S}}^{\top}(\mathbf{I}_{\tilde{S}}\bar{\mathbf{A}}\bar{\mathbf{A}}^{\top}\mathbf{I}_{\tilde{S}}^{\top})^{\dagger}\mathbf{I}_{\tilde{S}}\bar{\mathbf{A}}$
be the corresponding projection matrix. Denote
$\rho=\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{\tilde{S}}])$, then, running
Algorithm 2 with choice $s_{\max}=O(\log\frac{1}{\rho})$ will yield the
following convergence result:

$\displaystyle{\mathbb{E}}\|\mathbf{x}_{t+1}-\mathbf{x}^{*}\|^{2}\leq\left(1-\frac{\rho}{4}\right)\cdot\|\mathbf{x}_{t}-\mathbf{x}^{*}\|^{2}.$

Since we already obtain a lower bound for the quantity
$\rho=\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{\tilde{S}}])$ in Lemma 4.2,
by combining it with Lemma 4.3 we are ready to prove our main result.

###  4.4 Proof of Theorem 3.1

###### Proof of Theorem 3.1.

Notice that in Algorithm 2 we only carry out randomized Hadamard transform
once. By combining Lemma 4.2 and Lemma 4.3 we know that conditioned on an
event that holds with probability $0.99$, if we set
$s_{\max}=O(\log\frac{1}{\rho})=O(\log(r\bar{\kappa}_{k}/\tau))$, then we have

$\displaystyle{\mathbb{E}}\|\mathbf{x}_{t}-\mathbf{x}^{*}\|^{2}\leq\left(1-\frac{\rho}{4}\right)^{t}\cdot\|\mathbf{x}_{0}-\mathbf{x}^{*}\|^{2}\leq\left(1-\frac{\tau/r}{10\bar{\kappa}_{k}^{2}\log^{3}m}\right)^{t}\cdot\|\mathbf{x}_{0}-\mathbf{x}^{*}\|^{2}.$

By further choosing $\mathbf{x}_{0}=\mathbf{0},\tau=O(k\cdot\log^{3}m)$ and
$t_{\max}=O(r\bar{\kappa}_{k}^{2}/k\cdot\log(100/\epsilon))$ and denoting
$\tilde{\mathbf{x}}=\mathbf{x}_{t_{\max}}$, we have
${\mathbb{E}}\|\tilde{\mathbf{x}}-\mathbf{x}^{*}\|^{2}\leq\frac{\epsilon}{100}\|\mathbf{x}_{0}-\mathbf{x}^{*}\|^{2}=\frac{\epsilon}{100}\|\mathbf{x}^{*}\|^{2}$.
By using Markov’s inequality and applying a union bound, we have that
$\|\tilde{\mathbf{x}}-\mathbf{x}^{*}\|^{2}\leq\epsilon\|\mathbf{x}^{*}\|^{2}$
holds with probability $0.98$, and the total time complexity is

$\displaystyle~{}O(mn\log m)+t_{\max}\cdot
O(n\tau\log(\tau/\delta_{1})+\tau^{\omega}+n\tau s_{\max})$ $\displaystyle=$
$\displaystyle~{}O(mn\log
m)+r\bar{\kappa}_{k}^{2}/k\cdot\log(100/\epsilon)\cdot
O(nk\log^{3}m\log(r\bar{\kappa}_{k})+k^{\omega}\log^{3\omega}m)$
$\displaystyle=$ $\displaystyle~{}O(mn\log
m+(nr\log^{3}m\log(r\bar{\kappa}_{k})+rk^{\omega-1}\log^{3\omega}m)\cdot\bar{\kappa}_{k}^{2}\log
1/\epsilon)$ $\displaystyle=$ $\displaystyle~{}O(mn\log
m+(nr\log^{4}m+rk^{\omega-1}\log^{3\omega}m)\cdot\bar{\kappa}_{k}^{2}\log
1/\epsilon)$

where we use the fact that $\delta_{1}=\rho/324=\Omega(k/r)$ as stated in the
proof of Lemma 4.3 in Section 7.2, and the assumption $\bar{\kappa}_{k}\leq
m^{4}$ which gives that $\log(r\bar{\kappa}_{k})=O(\log m)$. ∎

##  5 Expected Projection under DPP Sampling

In this section we give a lower bound for the quantity
$\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}])$ where
$S_{\mathrm{DPP}}\sim k^{\prime}$-$\mathrm{DPP}(\mathbf{A}\mathbf{A}^{\top})$.
This quantity is directly associated with the convergence rate of approximated
sketch-and-project (see Lemma 4.3), and later we will show that
$\rho=\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{\tilde{S}}])\geq\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}])$,
where $\tilde{S}$ is the sampling scheme used in our main algorithms.

###  5.1 Elementary Symmetric Polynomials

We start by presenting the following two lemmas which help us relate the
expectation formulas that arise in sampling from $\mathrm{DPP}(\mathbf{L})$
(for some psd matrix $\mathbf{L}$) to properties of elementary symmetric
polynomials defined by the eigenvalues of $\mathbf{L}$. These lemmas will then
be used in the proof of Lemma 4.1, the main result of this section.

First we introduce Lemma 5.1 which is later used to compute the numerator of
Eq.(5.1), in terms of elementary symmetric polynomials of eigenvalues. Recall
that we use $\mathrm{adj}(\mathbf{L})\in\mathbb{R}^{n\times n}$ to denote the
adjugate of an $n\times n$ matrix $\mathbf{L}$, and note that when
$\mathbf{L}$ is invertible, then we have
$\mathrm{adj}(\mathbf{L})=\mathrm{det}(\mathbf{L})\mathbf{L}^{-1}$, which is
how the lemma connects to computing expectations for DPP sampling.

######  Lemma 5.1 (Lemma 3.3 of [RK20]).

Let $\mathbf{L}\in\mathbb{R}^{m\times m}$ be a real symmetric matrix with
eigenvalues $\lambda=(\lambda_{1},\ldots,\lambda_{m})$ where
$\lambda_{1}\geq\cdots\geq\lambda_{m}\geq 0$, let
$\mathbf{L}=\mathbf{Q}{\operatorname{diag}}(\lambda)\mathbf{Q}^{\top}$ be its
eigendecomposition where $\mathbf{Q}\in\mathbb{R}^{m\times m}$ is orthogonal,
and let $1\leq\ell\leq m$. Then we have

$\displaystyle\sum_{S\in\tbinom{[m]}{\ell}}\mathbf{I}_{S}\cdot\mathrm{adj}(\mathbf{L}_{S,S})\cdot\mathbf{I}_{S}^{\top}=\mathbf{Q}{\operatorname{diag}}(\sigma_{\ell-1}(\lambda_{-1}),\ldots,\sigma_{\ell-1}(\lambda_{-m}))\mathbf{Q}^{\top}$

where $\lambda_{-i}$ denotes the vector $\lambda$ without the $i$-th element
for $1\leq i\leq n$.

Next we introduce Lemma 5.2 which provides an upper bound for the ratio of
elementary symmetric polynomials by using Schur-concavity property of this
quantity. Lemma 5.2 will be used in Eq.(5.4).

######  Lemma 5.2 (Lemma 3.1 of [GS12]).

Given vector $\lambda=(\lambda_{1},\ldots,\lambda_{n})$ with
$\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{n}\geq 0$, given positive
integers $\tau\geq j>0$, we have

$\displaystyle\frac{\sigma_{\tau+1}(\lambda)}{\sigma_{\tau}(\lambda)}\leq\frac{1}{\tau+1-j}\sum_{i=j+1}^{n}\lambda_{i}.$

###  5.2 Proof of Lemma 4.1

With the above two lemmas, we are ready to prove Lemma 4.1, which gives a
lower bound for
$\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}])$ when
$S_{\mathrm{DPP}}\sim k^{\prime}$-$\mathrm{DPP}(\mathbf{A}\mathbf{A}^{\top})$.

###### Proof of Lemma 4.1.

For simplicity we denote $S=S_{\mathrm{DPP}}\sim
k^{\prime}$-$\mathrm{DPP}(\mathbf{A}\mathbf{A}^{\top})$, then the quantity we
consider is
$\lambda_{\min}^{+}(\mathbf{A}^{\top}{\mathbb{E}}[\mathbf{I}_{S}^{\top}(\mathbf{I}_{S}\mathbf{A}\mathbf{A}^{\top}\mathbf{I}_{S}^{\top})^{\dagger}\mathbf{I}_{S}]\mathbf{A})$.
Denote $\mathbf{B}=\mathbf{A}\mathbf{A}^{\top}$, notice that
$\mathbf{B}\in\mathbb{R}^{m\times m}$ is a PSD matrix with
$\mathrm{rank}(\mathbf{B})=r$. Denote $\\{\lambda_{i}\\}_{i=1}^{m}$ as the
eigenvalues of $\mathbf{B}$, then we have $\lambda_{i}=\sigma_{i}^{2}$ for
$i\leq r$ and $\lambda_{i}=0$ for $r+1\leq i\leq m$. Denote the singular value
decomposition of $\mathbf{A}$ as
$\mathbf{A}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$ where
$\mathbf{U}\in\mathbb{R}^{m\times m},\mathbf{\Sigma}\in\mathbb{R}^{m\times n}$
and $\mathbf{V}\in\mathbb{R}^{n\times n}$, then we have eigendecomposition
$\mathbf{B}=\mathbf{U}(\mathbf{\Sigma}\mathbf{\Sigma}^{\top})\mathbf{U}^{\top}$.
According to [RK20], assuming for simplicity that all the $k^{\prime}\times
k^{\prime}$ submatrices of $\mathbf{B}$ are non-degenerate (the other case can
be reduced to this), by using Cramer’s rule
$\mathrm{det}(\mathbf{B}_{S,S})\cdot(\mathbf{B}_{S,S})^{-1}=\mathrm{adj}(\mathbf{B}_{S,S})$
we have

$\displaystyle{\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}]=\mathbf{A}^{\top}{\mathbb{E}}[\mathbf{I}_{S}^{\top}(\mathbf{I}_{S}\mathbf{B}\mathbf{I}_{S}^{\top})^{\dagger}\mathbf{I}_{S}]\mathbf{A}=\frac{\sum_{S\in\tbinom{[m]}{k^{\prime}}}\mathbf{A}_{S,:}^{\top}\cdot\mathrm{adj}\left(\mathbf{B}_{S,S}\right)\cdot\mathbf{A}_{S,:}}{\sum_{S\in\tbinom{[m]}{k^{\prime}}}\mathrm{det}\left(\mathbf{B}_{S,S}\right)}.$
(5.1)

If we do not assume that all the $k^{\prime}\times k^{\prime}$ submatrices are
non-degenerate, then the numerator becomes
$\sum_{S\in\tbinom{[m]}{k^{\prime}}:\mathrm{det}(\mathbf{B}_{S,S})\neq
0}\mathbf{I}_{S}\cdot\mathrm{adj}\left(\mathbf{B}_{S,S}\right)\cdot\mathbf{I}_{S}^{T}$.
Then, for any vector $\mathbf{w}\in\mathbb{R}^{n}$ we have

$\displaystyle\mathbf{w}^{\top}{\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}]\mathbf{w}=$
$\displaystyle~{}\frac{\sum_{S\in\tbinom{[m]}{k^{\prime}}}\mathrm{det}(\mathbf{B}_{S,S})\cdot\mathbf{w}^{\top}\mathbf{A}_{S,:}^{\top}(\mathbf{B}_{S,S})^{\dagger}\mathbf{A}_{S,:}\mathbf{w}}{\sum_{S\in\tbinom{[m]}{k^{\prime}}}\mathrm{det}\left(\mathbf{B}_{S,S}\right)}$
$\displaystyle=$
$\displaystyle~{}\frac{\sum_{S\in\tbinom{[m]}{k^{\prime}}:\mathrm{det}(\mathbf{B}_{S,S})\neq
0}\mathbf{w}^{\top}\mathbf{A}_{S,:}^{\top}\mathrm{adj}(\mathbf{B}_{S,S})\mathbf{A}_{S,:}\mathbf{w}}{\sum_{S\in\tbinom{[m]}{k^{\prime}}}\mathrm{det}\left(\mathbf{B}_{S,S}\right)}$
(5.2)

where the second step also follows from Cramer’s rule. Notice that for any
vector $\mathbf{u}$ and matrix $\mathbf{C}$ we have the following formula for
adjugate matrix (see Fact 2.14.2 in [Ber09]):

$\displaystyle\mathrm{det}(\mathbf{C}+\mathbf{u}\mathbf{u}^{\top})=\mathrm{det}(\mathbf{C})+\mathbf{u}^{\top}\mathrm{adj}(\mathbf{C})\mathbf{u}.$

By applying this formula, for any $S$ such that
$\mathrm{det}(\mathbf{B}_{S,S})=0$, we can compute each component of numerator
in Eq.(5.2) as

$\displaystyle\mathbf{w}^{\top}\mathbf{A}_{S,:}^{\top}\mathrm{adj}(\mathbf{B}_{S,S})\mathbf{A}_{S,:}\mathbf{w}=$
$\displaystyle~{}\mathrm{det}(\mathbf{B}_{S,S}+\mathbf{A}_{S,:}\mathbf{w}\mathbf{w}^{\top}\mathbf{A}_{S,:}^{\top})-\mathrm{det}(\mathbf{B}_{S,S})$
$\displaystyle=$
$\displaystyle~{}\mathrm{det}(\mathbf{A}_{S,:}\mathbf{A}_{S,:}^{\top}+\mathbf{A}_{S,:}\mathbf{w}\mathbf{w}^{\top}\mathbf{A}_{S,:}^{\top})$
$\displaystyle=$
$\displaystyle~{}\mathrm{det}(\mathbf{A}_{S,:}(\mathbf{I}+\mathbf{w}\mathbf{w}^{\top})\mathbf{A}_{S,:}^{\top})=0$

where the last step follows since
$\mathrm{rank}(\mathbf{A}_{S,:}(\mathbf{I}+\mathbf{w}\mathbf{w}^{\top})\mathbf{A}_{S,:}^{\top})\leq\mathrm{rank}(\mathbf{A}_{S,:})<k^{\prime}$,
according to our assumption that
$\mathrm{det}(\mathbf{B}_{S,S})=\mathrm{det}(\mathbf{A}_{S,:}\mathbf{A}_{S,:}^{\top})=0$.
By applying this result to Eq.(5.2) we have

$\displaystyle\mathbf{w}^{\top}{\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}]\mathbf{w}=$
$\displaystyle~{}\frac{\sum_{S\in\tbinom{[m]}{k^{\prime}}:\mathrm{det}(\mathbf{B}_{S,S})\neq
0}\mathbf{w}^{\top}\mathbf{A}_{S,:}^{\top}\mathrm{adj}(\mathbf{B}_{S,S})\mathbf{A}_{S,:}\mathbf{w}}{\sum_{S\in\tbinom{[m]}{k^{\prime}}}\mathrm{det}\left(\mathbf{B}_{S,S}\right)}$
$\displaystyle=$
$\displaystyle~{}\frac{\sum_{S\in\tbinom{[m]}{k^{\prime}}}\mathbf{w}^{\top}\mathbf{A}_{S,:}^{\top}\mathrm{adj}(\mathbf{B}_{S,S})\mathbf{A}_{S,:}\mathbf{w}}{\sum_{S\in\tbinom{[m]}{k^{\prime}}}\mathrm{det}\left(\mathbf{B}_{S,S}\right)}$

which shows that the degenerate case can indeed be reduced to the non-
degenerate case as in Eq.(5.1), since we only care about the smallest positive
eigenvalue of ${\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}]$. According to
Lemma 2.1 and Lemma 5.1, we can compute the numerator and denominator
separately and have

$\displaystyle{\mathbb{E}}[\mathbf{I}_{S}^{\top}(\mathbf{I}_{S}\mathbf{B}\mathbf{I}_{S}^{\top})^{\dagger}\mathbf{I}_{S}]=\frac{\mathbf{U}{\operatorname{diag}}\left(\sigma_{k^{\prime}-1}\left(\lambda_{-1}\right),\ldots,\sigma_{k^{\prime}-1}\left(\lambda_{-m}\right)\right)\mathbf{U}^{\top}}{\sigma_{k^{\prime}}(\lambda)}$
(5.3)

where $\sigma_{k^{\prime}}$ is the real elementary symmetric polynomial with
degree $k^{\prime}$, and $\lambda_{-i}\in\mathbb{R}^{m-1}$ denotes the vector
$\lambda$ without the $i$-th element. Since we want to lower bound the above
quantity, we need to lower bound each
$\frac{\sigma_{k^{\prime}-1}(\lambda_{-i})}{\sigma_{k^{\prime}}(\lambda)}$ for
$i\in[m]$. We apply Lemma 5.2 and know for all $k^{\prime}>k$ we have

$\displaystyle\frac{\sigma_{k^{\prime}}(\lambda)}{\sigma_{k^{\prime}-1}(\lambda)}\leq\frac{1}{k^{\prime}-k}\sum_{j=k+1}^{m}\lambda_{j}.$

Thus for all $i\in[m]$ and $k^{\prime}>k$, we have

$\displaystyle\frac{\sigma_{k^{\prime}-1}(\lambda_{-i})}{\sigma_{k^{\prime}}(\lambda)}=\frac{\sigma_{k^{\prime}-1}(\lambda_{-i})}{\sigma_{k^{\prime}}(\lambda_{-i})+\lambda_{i}\cdot\sigma_{k^{\prime}-1}(\lambda_{-i})}=\frac{1}{\lambda_{i}+\frac{\sigma_{k^{\prime}}(\lambda_{-i})}{\sigma_{k^{\prime}-1}(\lambda_{-i})}}\geq\frac{1}{\lambda_{i}+\frac{1}{k^{\prime}-k}\sum_{j>k,j\neq
i}\lambda_{j}}.$ (5.4)

We consider two different cases for the value of $i$ and $k$, and have

$\displaystyle\begin{cases}\frac{\sigma_{k^{\prime}-1}(\lambda_{-i})}{\sigma_{k^{\prime}}(\lambda)}\geq\frac{1}{\lambda_{i}+\frac{1}{k^{\prime}-k}\sum_{j>k}\lambda_{j}}~{}~{}~{}\text{for}~{}~{}~{}i\leq
k;\\\
\frac{\sigma_{k^{\prime}-1}(\lambda_{-i})}{\sigma_{k^{\prime}}(\lambda)}\geq\frac{1}{\frac{k^{\prime}-k-1}{k^{\prime}-k}\lambda_{i}+\frac{1}{k^{\prime}-k}\sum_{j>k}\lambda_{j}}~{}~{}~{}\text{for}~{}~{}~{}k+1\leq
i\leq m.\end{cases}$

Based on these two cases, we construct the following approximation matrix:

$\displaystyle\mathbf{B}_{k^{\prime}}:=\mathbf{U}{\operatorname{diag}}\left(\lambda_{1},\ldots,\lambda_{k},\frac{k^{\prime}-k-1}{k^{\prime}-k}\lambda_{k+1},\ldots,\frac{k^{\prime}-k-1}{k^{\prime}-k}\lambda_{m}\right)\mathbf{U}^{\top}+\frac{1}{k^{\prime}-k}\sum_{j>k}^{m}\lambda_{j}\mathbf{I}$

an have
${\mathbb{E}}[\mathbf{I}_{S}^{\top}(\mathbf{I}_{S}\mathbf{B}\mathbf{I}_{S}^{\top})^{\dagger}\mathbf{I}_{S}]\succeq\mathbf{B}_{k^{\prime}}^{-1}$.
By applying this result we can lower bound
$\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{S}])$.

$\displaystyle\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{S}])=$
$\displaystyle~{}\lambda_{\min}^{+}(\mathbf{A}^{\top}{\mathbb{E}}[\mathbf{I}_{S}^{\top}(\mathbf{I}_{S}\mathbf{B}\mathbf{I}_{S}^{\top})^{\dagger}\mathbf{I}_{S}]\mathbf{A})$
$\displaystyle\geq$
$\displaystyle~{}\lambda_{\min}^{+}(\mathbf{A}^{\top}\mathbf{B}_{k^{\prime}}^{-1}\mathbf{A})$
$\displaystyle=$
$\displaystyle~{}\lambda_{\min}^{+}(\mathbf{V}\mathbf{\Sigma}^{\top}\mathbf{U}^{\top}\mathbf{B}_{k^{\prime}}^{-1}\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top})$
$\displaystyle=$
$\displaystyle~{}\lambda_{\min}^{+}(\mathbf{V}{\operatorname{diag}}(\frac{\lambda_{1}}{\lambda_{1}+\frac{1}{k^{\prime}-k}\sum_{j>k}^{m}\lambda_{j}},\ldots,\frac{\lambda_{k}}{\lambda_{k}+\frac{1}{k^{\prime}-k}\sum_{j>k}^{m}\lambda_{j}},$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}\frac{\lambda_{k+1}}{\frac{k^{\prime}-k-1}{k^{\prime}-k}\lambda_{k+1}+\frac{1}{k^{\prime}-k}\sum_{j>k}^{m}\lambda_{j}},\ldots,\frac{\lambda_{n}}{\frac{k^{\prime}-k-1}{k^{\prime}-k}\lambda_{n}+\frac{1}{k^{\prime}-k}\sum_{j>k}^{m}\lambda_{j}})\mathbf{V}^{\top})$
$\displaystyle=$
$\displaystyle~{}\frac{\sigma_{r}^{2}}{\frac{k^{\prime}-k-1}{k^{\prime}-k}\sigma_{r}^{2}+\frac{1}{k^{\prime}-k}\sum_{j>k}^{r}\sigma_{j}^{2}}$

where the last step follows from $\lambda_{j}=\sigma_{j}^{2}$ for $j\in[r]$,
and that $\lambda_{r+1}=\cdots=\lambda_{m}=0$. Recall that we define
$\bar{\kappa}_{k}:=(\frac{1}{r-k}\sum_{j>k}^{r}\sigma_{j}^{2}/\sigma_{r}^{2})^{1/2}$,
thus we can simplify the above expression as

$\displaystyle\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{S}])\geq\frac{1}{\frac{k^{\prime}-k-1}{k^{\prime}-k}+\frac{1}{k^{\prime}-k}\sum_{j>k}^{r}\sigma_{j}^{2}/\sigma_{r}^{2}}=\frac{k^{\prime}-k}{k^{\prime}-k-1+(r-k)\bar{\kappa}_{k}^{2}}.$

∎

#### Discussion.

Let $\mathbf{B}=\mathbf{A}\mathbf{A}^{\top}\in\mathbb{R}^{m\times m}$ with
rank $r$, in comparison with our result, [RK20] constructed the
“$k^{\prime}$-coordinate approximation” of $\mathbf{B}$ as

$\displaystyle\hat{\mathbf{B}}_{k^{\prime}}:=\mathbf{U}{\operatorname{diag}}(\lambda_{1},\ldots,\lambda_{k^{\prime}},\lambda_{k^{\prime}},\ldots,\lambda_{k^{\prime}})\mathbf{U}^{\top}+\sum_{j=k^{\prime}+1}^{m}\lambda_{j}\mathbf{I}$

and showed that
$\mathbb{E}[\mathbf{I}_{S}^{\top}(\mathbf{I}_{S}\mathbf{B}\mathbf{I}_{S}^{\top})^{\dagger}\mathbf{I}_{S}]\succeq\hat{\mathbf{B}}_{k^{\prime}}^{-1}$.
By using their result, one could show the following weaker lower bound:

$\displaystyle\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{S}])\geq\lambda_{\min}^{+}(\mathbf{A}^{\top}\hat{\mathbf{B}}_{k^{\prime}}^{-1}\mathbf{A})=\frac{\lambda_{r}}{\sum_{j=k^{\prime}}^{m}\lambda_{j}}=\frac{\sigma_{r}^{2}}{\sum_{j=k^{\prime}}^{r}\sigma_{j}^{2}}=\frac{1}{(r-k^{\prime}+1)\bar{\kappa}_{k^{\prime}-1}^{2}}.$
(5.5)

By comparing this result with the lower bound we get in Eq.(4.2), we can see
that we have an extra parameter $k$ that only needs to satisfy $k<k^{\prime}$.
If we choose $k=k^{\prime}-1$ then we recover the same result as in [RK20].
However, if we consider the setting where $k$ is such that
$\frac{\sigma_{i}}{\sigma_{r}}=O(1)$ holds for all $i\geq k$ and assume that
$k\leq\frac{r}{2}$, then by choosing $k^{\prime}=2k-1<r$, Eq.(5.5) only gives
a lower bound of order $\Omega(\frac{1}{r-k})$. This is worse by a factor of
$k$ than our result which is of order $\Omega(\frac{k}{r})$.

##  6 Replacing DPP with Uniform Sampling

In Lemma 4.1 we give a lower bound for
$\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{S}])$ when $S$ is a $k$-DPP
sample.333Note that in this section we use $k$ to denote DPP size for
simplicity, which corresponds to $k^{\prime}$ in previous sections. However,
obtaining an exact DPP sample can be expensive. Currently, the best known
result of obtaining $s$ independent approximate samples takes
$\tilde{O}(mk^{\omega-1}+sk^{\omega})$, see Lemma 2.3, and adapting this
result to our rectangular matrix setting takes
$\tilde{O}(mnk^{\omega-2}+snk^{\omega})$, see Section 6.4. This sampling
method consists of two steps: a preprocessing step which takes
$\tilde{O}(mnk^{\omega-2})$ and a sampling step which takes
$\tilde{O}(nk^{\omega})$ for each sample. Notice that _both these two steps
are bottlenecks and unaffordable_ when we use the sketch-and-project method
with $k$-DPP sampling, since with the outer iteration
$t_{\max}=O(r\bar{\kappa}_{k}^{2}/k\cdot\log(1/\epsilon))$, the cost will be
$\tilde{O}(mnk^{\omega-2}+nr\bar{\kappa}_{k}^{2}k^{\omega-2}\log(1/\epsilon))$,
i.e., at least as much as standard preconditioning for CG via power iteration.

To conquer the first bottleneck, we use a different method of preprocessing.
The purpose of the preprocessing step is to make the marginal probabilities
$\Pr\\{j\in S_{\mathrm{DPP}}\\}$ almost the same. [ALV22] do this by
estimating all of the marginal probabilities, which itself requires repeatedly
sampling from the $k$-DPP distribution. We instead show in Lemma 6.4 that the
randomized Hadamard transform (RHT) can uniformize the marginals _in near-
linear time_ $O(mn\log m)$, which is much better than
$\tilde{O}(mnk^{\omega-2})$. With the approximately uniform marginals, we can
use the following result from [ALV22] to obtain an approximated $k$-DPP
sample.

######  Lemma 6.1 (Theorem 3 in [ALV22]).

Given a strongly Rayleigh distribution $\mu\in\mathbb{R}^{\tbinom{[m]}{k}}$
and marginal overestimates $q_{i}\geq\Pr\\{i\in T\\}$ for $i\in[m]$. Let
$K:=\sum_{i\in[m]}q_{i}$, there is an algorithm that outputs a sample from a
distribution with total variation distance $m^{-O(1)}$ from $\mu$ in time
bounded by $O(\log^{3}m)$ calls to $\mathcal{T}_{\mu}(O(K),k)$, where
$\mathcal{T}_{\mu}(t,k)$ is the time it takes to produce a sample from $\mu$
conditioned on all elements of the sample being a subset of a fixed set $T$ of
size $|T|=t$.

######  Remark 6.1.

In Algorithm 1, we give the sampling procedure claimed by Lemma 6.1 under the
additional “isotropic” assumption, i.e., that the marginal overestimates are
all uniform. Moreover, the $m^{-O(1)}$ term in Lemma 6.1 is specified as
$m^{-9}$ in the proof of [ALV22].

Algorithm 4 $k$-DPP sampling algorithm based on [ALV22] (Ghost Algorithm).

1:Input: $m,k,t,\mathrm{iter},\mathrm{call},\mathcal{O}$;

2:Initialize $S_{0}\in\tbinom{[m]}{k},T_{0}\leftarrow\varnothing$;
$\triangleright$ By calling $\tilde{O}(k)$ times the oracle.

3:for $i=0,1,2,\ldots,\mathrm{iter}-1$ do $\triangleright$
$\mathrm{iter}=O(\log^{3}m)$.

4: Initialize $T_{i,0}\leftarrow T_{i}$;

5: for $j=0,1,2,\ldots,\mathrm{call}-1$ do $\triangleright$ $\mathrm{call}\geq
Ck$ for some $C=O(1)$.

6: Sample $i_{j}$ by calling $\mathcal{O}$;

7: if $i_{j}\in S_{i}\cup T_{i,j}$ then:

8: $T_{i,j+1}=T_{i,j}$;

9: else:

10: $T_{i,j+1}=T_{i,j}\cup\\{i_{j}\\}$;

11: end if

12: end for

13: Update $T_{i+1}\leftarrow T_{i,\mathrm{call}}$; $\triangleright$ Up-walk
step by calling $\tilde{O}(t-k)$ times the oracle.

14: Update $S_{i+1}\subseteq S_{i}\cup T_{i+1}$; $\triangleright$ Down-walk
sampling from a $k$-DPP restricted to $S_{i}\cup T_{i+1}$.

15:end for

16:return $S_{\mathrm{iter}}$.

To conquer the second bottleneck, we observe that compared with a real $k$-DPP
sample, we are actually more satisfied with a larger set that _contains_ a
$k$-DPP in the sketch-and-project scheme, since it can guarantee an even
faster convergence rate, see Lemma 6.3. How do we obtain a sample that
contains a $k$-DPP as a subset? We show that we can achieve this by using
_only uniform sampling_. Suppose we have an oracle $\mathcal{O}$ that can
uniformly sample one element from $[m]$, and we would like to obtain
$\tilde{S}$ such that $\tilde{S}$ contains a $k$-DPP sample. Our idea is to
rewrite Algorithm 1 so that it can only access the elements via the oracle
$\mathcal{O}$: starting from $S_{0}$ with size $k$, we call the oracle several
times, if the sample is new then we add it to $S_{0}$, otherwise we discard
it. After we obtain $S_{0}$, we follow the same rule to obtain $T_{1}$ such
that (i) $|T_{1}|=t-k$ and (ii) $T_{1}\cap S_{0}=\varnothing$. Then, we down
sample $S_{1}$ such that (i) $|S_{1}|=k$ and (ii) $S_{1}\subseteq S_{0}\cup
T_{1}$. These two steps finish one iteration. Notice that only step 1 needs to
call the oracle, while step 2 does not.

These steps correspond to Algorithm 4, and we call it the “Ghost Algorithm”
since we do not implement it due to the high cost of the down-walk step (line
$14$ of Algorithm 4). Instead, we use the coupled Algorithm 5 and show that
the output of Algorithm 5 contains the output of Algorithm 4 as a subset. By
further showing that the output of Algorithm 4 is the same as a $k$-DPP sample
with high probability, we have the following theorem which shows that the
output of Algorithm 5 indeed contains a $k$-DPP sample. Notice that since we
only use uniform sampling in our algorithm, the total cost of sampling takes
only $O((k+\log m)\log^{3}m)=\tilde{O}(k)$ which is much better than
$\tilde{O}(nk^{\omega})$ and almost linear in $k$. For detailed discussion see
Section 6.4. In this whole section, we denote
$\mathbf{Q}=\frac{1}{\sqrt{m}}\mathbf{H}\mathbf{D}$ as the RHT matrix, where
$\mathbf{H}\in\mathbb{R}^{m\times m}$ is the Hadamard matrix and
$\mathbf{D}\in\mathbb{R}^{m\times m}$ is the diagonal matrix with Rademacher
entries. We denote $\bar{\mathbf{A}}=\mathbf{Q}\mathbf{A}$ as the preprocessed
matrix.

######  Lemma 6.2 (Coupling).

Given matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$, denote
$\mathbf{Q}=\frac{1}{\sqrt{m}}\mathbf{H}\mathbf{D}$ as the RHT matrix and
$\bar{\mathbf{A}}=\mathbf{Q}\mathbf{A}$ as the preprocessed matrix. Suppose we
have an oracle $\mathcal{O}$ that can uniformly sample one element from $[m]$.
Conditioned on an event that holds with probability $0.99$ (and only depends
on RHT), for $k\geq\log m$ and $\delta=m^{-O(1)}$, there is a coupling
$(\tilde{S},S_{\mathrm{DPP}})$ where $\tilde{S}$ is sampled according to
Algorithm 5 with choice $\tau\geq Ck\log^{3}m$ for some constant $C=O(1)$, and
$S_{\mathrm{DPP}}\sim
k$-$\mathrm{DPP}(\bar{\mathbf{A}}\bar{\mathbf{A}}^{\top})$, such that with
probability $1-\delta$, we have $S_{\mathrm{DPP}}\subseteq\tilde{S}$.

Algorithm 5 Sample by using uniform sampling oracle $\mathcal{O}$.

1:Input: $m,\tau,\mathcal{O}$;

2:Initialize $\tilde{S}_{0}\leftarrow\varnothing$;

3:for $\ell=0,1,2,\ldots,\tau-1$ do

4: Sample $i_{\ell}$ by calling $\mathcal{O}$;

5: if $i_{\ell}\in\tilde{S}_{\ell}$ then:

6: $\tilde{S}_{\ell+1}=\tilde{S}_{\ell}$;

7: else:

8: $\tilde{S}_{\ell+1}=\tilde{S}_{\ell}\cup\\{i_{\ell}\\}$;

9: end if

10:end for

11:return $\tilde{S}=\tilde{S}_{N}$.

To prove Lemma 6.2, we show the following in this section:

  1. 1.

After the preprocessing step, the marginals are approximately isotropic. See
Lemma 6.4.

  2. 2.

In Algorithm 4, with probability $1-\delta$, we can obtain each
$T_{i}\in\tbinom{[m]}{t-k}$ after $O(t-k+\log\frac{1}{\delta})$ calls to the
oracle. See Lemma 6.5.

  3. 3.

In Algorithm 4, if we set $\mathrm{call}\geq Ck$ for some constant $C=O(1)$
and $\mathrm{iter}=O(\log^{3}m)$, then with probability $1-\delta$,
$S_{\mathrm{iter}}$ is the same as $S_{\mathrm{DPP}}$. See Lemma 6.6.

#### Monotonicity of Projections.

We start our analysis by making the following simple observation, which
implies that in the sketch-and-project scheme, a larger set can only improve
the convergence guarantee, since it leads to a “larger” projection matrix in
the sense of matrix orders. Formally we have the following lemma.

######  Lemma 6.3.

Given matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$, for any subsets
$S_{1}\subseteq S_{2}\subseteq[m]$, define projection matrix
$\mathbf{P}_{S_{i}}:=\mathbf{A}^{\top}\mathbf{I}_{S_{i}}^{\top}(\mathbf{I}_{S_{i}}\mathbf{A}\mathbf{A}^{\top}\mathbf{I}_{S_{i}}^{\top})^{\dagger}\mathbf{I}_{S_{i}}\mathbf{A}$
where $\mathbf{I}_{S_{i}}\in\mathbb{R}^{|S_{i}|\times m}$ is the sampling
matrix. Then we have

$\displaystyle\mathbf{P}_{S_{1}}\preceq\mathbf{P}_{S_{2}}.$

###### Proof.

We only need to show that for any vector $\mathbf{x}\in\mathbb{R}^{n}$ such
that $\|\mathbf{x}\|=1$, we have
$\mathbf{x}^{\top}\mathbf{P}_{S_{1}}\mathbf{x}\leq\mathbf{x}^{\top}\mathbf{P}_{S_{2}}\mathbf{x}$.
Let $\mathrm{span}\langle\mathbf{A}_{S_{i}}\rangle$ be the subspace spanned by
the rows of $\mathbf{I}_{S_{i}}\mathbf{A}$. We have

$\displaystyle\mathbf{x}^{\top}\mathbf{P}_{S_{i}}\mathbf{x}=1-\mathbf{x}^{\top}(\mathbf{I}-\mathbf{P}_{S_{i}})\mathbf{x}=1-\|\mathbf{x}-\mathbf{P}_{S_{i}}\mathbf{x}\|^{2}=1-\min_{\mathbf{x}_{i}\in\mathrm{span}\langle\mathbf{A}_{S_{i}}\rangle}\|\mathbf{x}-\mathbf{x}_{i}\|^{2}.$

Note that since $S_{1}\subseteq S_{2}$, we know that
$\mathrm{span}\langle\mathbf{A}_{S_{2}}\rangle$ contains
$\mathrm{span}\langle\mathbf{A}_{S_{1}}\rangle$, thus we have:

$\displaystyle\min_{\mathbf{x}_{1}\in\mathrm{span}\langle\mathbf{A}_{S_{1}}\rangle}\|\mathbf{x}-\mathbf{x}_{1}\|^{2}\geq\min_{\mathbf{x}_{2}\in\mathrm{span}\langle\mathbf{A}_{S_{2}}\rangle}\|\mathbf{x}-\mathbf{x}_{2}\|^{2}.$

We conclude that
$\mathbf{x}^{\top}\mathbf{P}_{S_{1}}\mathbf{x}\leq\mathbf{x}^{\top}\mathbf{P}_{S_{2}}\mathbf{x}$
holds for all $\mathbf{x}\in\mathbb{R}^{n}$ such that $\|\mathbf{x}\|=1$, as
desired. ∎

According to Lemma 6.3, $\mathbf{P}_{S_{1}}\preceq\mathbf{P}_{S_{2}}$ holds
for all $S_{1}\subseteq S_{2}$. Thus we can take expectation over the
randomness of $S_{1}$ and $S_{2}$ and have
$\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{S_{1}}])\leq\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{S_{2}}])$
holds conditioned on the event $S_{1}\subseteq S_{2}$. Thus, if we can obtain
a sample $\tilde{S}$ that contains a $k$-DPP sample $S_{\mathrm{DPP}}$, then
we know that
$\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{\tilde{S}}])\geq\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}])$
which guarantees a faster convergence rate for our algorithm.

###  6.1 Preprocessing

In this section we show that the $k$-DPP marginals become approximately
isotropic (i.e., uniform) after we use the randomized Hadamard transform (RHT)
as a preprocessing step. We note that prior works have shown that the RHT can
uniformize the entries of a vector, or the row-norms or leverage scores of a
matrix [Tro11]. However, this analysis typically relies on the fact that all
of these quantities can be easily described as linear or quadratic functions
of the rows of the target matrix. On the other hand, the marginals of a
general fixed size DPP are not expressed with any simple formula (they can
only be described with expressions that involve elementary symmetric
polynomials of certain sub-matrices of the input matrix). We get around this
issue by relying on the fact that $k$-DPP marginals can be described as a
mixture of the marginals of simpler Projection DPPs (this fact has been
previously used in $k$-DPP sampling algorithms [KT12]).

First, we show that applying the RHT as in Algorithm 2 does not affect the
optimal solution of the linear system $\mathbf{A}\mathbf{x}=\mathbf{b}$.
Denote
$\bar{\mathbf{A}}=\mathbf{Q}\mathbf{A},\bar{\mathbf{b}}=\mathbf{Q}\mathbf{b}$
and consider the linear system $\bar{\mathbf{A}}\mathbf{x}=\bar{\mathbf{b}}$.
Since $\mathbf{Q}^{\top}\mathbf{Q}=\mathbf{I}$, we know the optimal solution
is

$\displaystyle\bar{\mathbf{A}}^{\dagger}\bar{\mathbf{b}}=(\mathbf{Q}\mathbf{A})^{\dagger}\mathbf{Q}\mathbf{b}=(\mathbf{A}^{\top}\mathbf{Q}^{\top}\mathbf{Q}\mathbf{A})^{\dagger}\mathbf{A}^{\top}\mathbf{Q}^{\top}\mathbf{Q}\mathbf{b}=(\mathbf{A}^{\top}\mathbf{A})^{\dagger}\mathbf{A}^{\top}\mathbf{b}=\mathbf{A}^{\dagger}\mathbf{b}=\mathbf{x}^{*},$

which is the same as the optimal solution of
$\mathbf{A}\mathbf{x}=\mathbf{b}$. To construct this system, we need to
compute $\mathbf{Q}\mathbf{A}$ and $\mathbf{Q}\mathbf{b}$ which takes
$O(mn\log m)=\tilde{O}(mn)$ time, this is much better than the
$\tilde{O}(mnk^{\omega-2})$ result in Section 6.4. Next we show that the
marginals of $k$-$\mathrm{DPP}(\bar{\mathbf{A}}\bar{\mathbf{A}}^{\top})$ are
approximately isotropic.

######  Lemma 6.4 (Approximately isotropic $k$-DPP).

Given matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$, denote
$\mathbf{Q}=\frac{1}{\sqrt{m}}\mathbf{H}\mathbf{D}$ as the RHT matrix and
$\bar{\mathbf{A}}=\mathbf{Q}\mathbf{A}$ as the preprocessed matrix. Let
$S_{\mathrm{DPP}}\sim
k$-$\mathrm{DPP}(\bar{\mathbf{A}}\bar{\mathbf{A}}^{\top})$, for any
$0<\delta<1$, with probability $1-\delta$ we have

$\displaystyle\Pr\\{j\in
S_{\mathrm{DPP}}\\}\leq\left(\sqrt{\frac{k}{m}}+\sqrt{\frac{8\log(m/\delta)}{m}}\right)^{2}~{}~{}~{}\text{for
all}~{}~{}~{}j\in[m].$

###### Proof.

Although there is no explicit expression for the marginals of general fixed
size DPP, there is one specific class of fixed size DPP called Projection DPP
which has simple expressions for the marginals. Moreover, sampling a
$k$-$\mathrm{DPP}$ is equivalent to sampling a mixture of Projection DPPs, see
Lemma 2.2. Let $\mathbf{A}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$ be the
SVD of $\mathbf{A}$ and denote $\bar{\mathbf{U}}=\mathbf{Q}\mathbf{U}$, then
we have eigendecomposition
$\bar{\mathbf{A}}\bar{\mathbf{A}}^{\top}=\bar{\mathbf{U}}(\mathbf{\Sigma}\mathbf{\Sigma}^{\top})\bar{\mathbf{U}}^{\top}$.
Denote
$r=\mathrm{rank}(\bar{\mathbf{A}}\bar{\mathbf{A}}^{\top})=\mathrm{rank}(\mathbf{A})$.
Lemma 2.2 shows that if we first sample
$\gamma=(\gamma_{1},\ldots,\gamma_{r})$ where
$\gamma_{i}\sim\mathrm{Bernoulli}(\frac{\lambda_{i}}{\lambda_{i}+1})$ are
independent, and then sample
$S_{\gamma}\sim\mathrm{P}$-$\mathrm{DPP}(\bar{\mathbf{U}}_{:,\gamma}\bar{\mathbf{U}}_{:,\gamma}^{\top})$,
where $\bar{\mathbf{U}}_{:,\gamma}$ is the submatrix of columns of
$\bar{\mathbf{U}}$ indexed by the non-zeros in $\gamma$, then conditioned on
the event that $\sum_{i}\gamma_{i}=k$, $S_{\gamma}$ is distributed the same as
$S_{\mathrm{DPP}}\sim
k$-$\mathrm{DPP}(\bar{\mathbf{A}}\bar{\mathbf{A}}^{\top})$.

The advantage of decomposing $k$-DPP into $\mathrm{P}$-DPPs is that we can
explicitly express marginals of $\mathrm{P}$-DPP. Since
$\bar{\mathbf{U}}_{:,\gamma}\bar{\mathbf{U}}_{:,\gamma}^{\top}=\bar{\mathbf{U}}_{:,\gamma}\bar{\mathbf{U}}_{:,\gamma}^{\dagger}$
is a projection matrix, we know that the marginals of
$\mathrm{P}$-$\mathrm{DPP}(\bar{\mathbf{U}}_{:,\gamma}\bar{\mathbf{U}}_{:,\gamma}^{\top})$
can be expressed as the diagonal entries of this kernel matrix as (e.g., see
[DM21])

$\displaystyle\Pr\\{j\in
S_{\gamma}\\}=(\bar{\mathbf{U}}_{:,\gamma}\bar{\mathbf{U}}_{:,\gamma}^{\top})_{j,j}=\|\mathbf{e}_{j}^{\top}\bar{\mathbf{U}}_{:,\gamma}\|^{2}~{}~{}~{}\text{where}~{}~{}~{}j\in[m].$

Based on this result, we can compute the marginals of $S_{\mathrm{DPP}}$ by
summing over the marginals of all possible P-DPPs. For $j\in[m]$ we have

$\displaystyle\Pr\\{j\in
S_{\mathrm{DPP}}\\}=\sum_{\gamma\in\\{0,1\\}^{r}}\left(\Pr\\{j\in
S_{\gamma}\\}\cdot\Pr\\{\gamma\mid\sum_{i}\gamma_{i}=k\\}\right)=\sum_{\gamma\in\\{0,1\\}^{r}}\left(\|\mathbf{e}_{j}^{\top}\bar{\mathbf{U}}_{:,\gamma}\|^{2}\cdot\Pr\\{\gamma\mid\sum_{i}\gamma_{i}=k\\}\right).$

Note that by definition we have

$\displaystyle\|\mathbf{e}_{j}^{\top}\bar{\mathbf{U}}_{:,\gamma}\|^{2}=\|\mathbf{e}_{j}^{\top}\mathbf{Q}\mathbf{U}_{:,\gamma}\|^{2}=\frac{1}{m}\|\mathbf{e}_{j}^{\top}\mathbf{H}\mathbf{D}\mathbf{U}_{:,\gamma}\|^{2}=\frac{1}{m}\|\mathbf{d}^{\top}{\operatorname{diag}}(\mathbf{e}_{j}^{\top}\mathbf{H})\mathbf{U}_{:,\gamma}\|^{2},$

where $\mathbf{d}$ is the vector of Rademacher variables along the diagonal of
$\mathbf{D}$. If we denote
$\mathbf{E}_{j}:={\operatorname{diag}}(\mathbf{e}_{j}^{\top}\mathbf{H})$ and
$p_{\gamma}:=\Pr\\{\gamma\mid\sum_{i}\gamma_{i}=k\\}$, then we have

$\displaystyle\Pr\\{j\in S_{\mathrm{DPP}}\\}=$
$\displaystyle~{}\sum_{\gamma\in\\{0,1\\}^{r}}p_{\gamma}\cdot\|\mathbf{e}_{j}^{\top}\bar{\mathbf{U}}_{:,\gamma}\|^{2}$
$\displaystyle=$
$\displaystyle~{}\frac{1}{m}\sum_{\gamma\in\\{0,1\\}^{r}}p_{\gamma}\cdot\|\mathbf{d}^{\top}\mathbf{E}_{j}\mathbf{U}_{:,\gamma}\|^{2}.$
$\displaystyle=$
$\displaystyle~{}\frac{1}{m}\sum_{\gamma\in\\{0,1\\}^{r}}p_{\gamma}\cdot\mathbf{d}^{\top}\mathbf{E}_{j}\mathbf{U}_{:,\gamma}\mathbf{U}_{:,\gamma}^{\top}\mathbf{E}_{j}\mathbf{d}$
$\displaystyle=$
$\displaystyle~{}\frac{1}{m}\mathbf{d}^{\top}\left(\sum_{\gamma\in\\{0,1\\}^{r}}p_{\gamma}\mathbf{E}_{j}\mathbf{U}_{:,\gamma}\mathbf{U}_{:,\gamma}^{\top}\mathbf{E}_{j}\right)\mathbf{d}.$

We define function
$f_{j}(\mathbf{x})=\sqrt{\frac{1}{m}\mathbf{x}^{\top}\mathbf{B}_{j}\mathbf{x}}$
where
$\mathbf{B}_{j}:=\sum_{\gamma\in\\{0,1\\}^{r}}p_{\gamma}\mathbf{E}_{j}\mathbf{U}_{:,\gamma}\mathbf{U}_{:,\gamma}^{\top}\mathbf{E}_{j}$.
We first analyze the properties of $\mathbf{B}_{j}$. Since
$\mathbf{E}_{j}^{2}=\mathbf{I}$ and
$\mathbf{U}_{:,\gamma}^{\top}\mathbf{U}_{:,\gamma}=\mathbf{I}$, we have

$\displaystyle\|\mathbf{B}_{j}\|\leq\sum_{\gamma\in\\{0,1\\}^{r}}p_{\gamma}\cdot\|\mathbf{E}_{j}\mathbf{U}_{:,\gamma}\mathbf{U}_{:,\gamma}^{\top}\mathbf{E}_{j}\|\leq\sum_{\gamma\in\\{0,1\\}^{r}}p_{\gamma}\cdot\|\mathbf{U}_{:,\gamma}\mathbf{U}_{:,\gamma}^{\top}\|=1$

and

$\displaystyle{\operatorname{tr}}(\mathbf{B}_{j})=\sum_{\gamma\in\\{0,1\\}^{r}}p_{\gamma}\cdot{\operatorname{tr}}(\mathbf{E}_{j}\mathbf{U}_{:,\gamma}\mathbf{U}_{:,\gamma}^{\top}\mathbf{E}_{j})=\sum_{\gamma\in\\{0,1\\}^{r}}p_{\gamma}\cdot{\operatorname{tr}}(\mathbf{U}_{:,\gamma}^{\top}\mathbf{U}_{:,\gamma})=k\cdot\left(\sum_{\gamma\in\\{0,1\\}^{r}}p_{\gamma}\right)=k.$

Next we analyze the properties of $f_{j}$. Note that $f_{j}$ is convex, with
Lipschitz constant as follows:

$\displaystyle|f_{j}(\mathbf{x})-f_{j}(\mathbf{y})|=\frac{1}{\sqrt{m}}|\|\mathbf{x}\|_{\mathbf{B}_{j}}-\|\mathbf{y}\|_{\mathbf{B}_{j}}|\leq\frac{1}{\sqrt{m}}\|\mathbf{x}-\mathbf{y}\|_{\mathbf{B}_{j}}\leq\frac{1}{\sqrt{m}}\|\mathbf{B}_{j}\|\cdot\|\mathbf{x}-\mathbf{y}\|\leq\frac{1}{\sqrt{m}}\|\mathbf{x}-\mathbf{y}\|.$

We consider the random variable $f_{j}(\mathbf{d})$ and bound the expectation
as follows:

$\displaystyle{\mathbb{E}}[f_{j}(\mathbf{d})]\leq({\mathbb{E}}[f_{j}(\mathbf{d})^{2}])^{1/2}=\left(\frac{1}{m}{\mathbb{E}}[\mathbf{d}^{\top}\mathbf{B}_{j}\mathbf{d}]\right)^{1/2}=\left(\frac{1}{m}{\operatorname{tr}}(\mathbf{B}_{j})\right)^{1/2}=\sqrt{\frac{k}{m}}.$

Then we use Lemma 2.5 on the concentration of a Lipschitz function of a
Rademacher random vector, with $t=\sqrt{8\log(m/\delta)}$, and obtain

$\displaystyle\Pr\left\\{\sqrt{\frac{1}{m}\mathbf{d}^{\top}\mathbf{B}_{j}\mathbf{d}}\geq\sqrt{\frac{k}{m}}+\sqrt{\frac{8\log(m/\delta)}{m}}\right\\}\leq\frac{\delta}{m}.$

Further notice that $\Pr\\{j\in
S_{\mathrm{DPP}}\\}=f_{j}^{2}(\mathbf{d})=\frac{1}{m}\mathbf{d}^{\top}\mathbf{B}_{j}\mathbf{d}$,
thus by applying union bound we have that with probability $1-\delta$,

$\displaystyle\Pr\\{j\in
S_{\mathrm{DPP}}\\}\leq\left(\sqrt{\frac{k}{m}}+\sqrt{\frac{8\log(m/\delta)}{m}}\right)^{2}~{}~{}~{}\text{for
all}~{}~{}~{}j\in[m].$

∎

######  Remark 6.2.

Since we assume that $k\geq\log m$, by setting $\delta=0.01$ we have

$\displaystyle\Pr\\{j\in
S_{\mathrm{DPP}}\\}\leq\left(\sqrt{\frac{k}{m}}+\sqrt{\frac{8\log(m/\delta)}{m}}\right)^{2}\leq\frac{2k}{m}+\frac{16\log(m/\delta)}{m}\leq\frac{c_{0}k}{m}$

for some constant $c_{0}=O(1)$. This shows that after the randomized Hadamard
transform, the marginals are approximately isotropic up to a constant factor
with probability $0.99$.

###  6.2 Successful DPP Sample

In this section we show that as long as the marginals are approximately
isotropic, then with high probability, the output of Algorithm 4 is the same
as a real $k$-$\mathrm{DPP}(\bar{\mathbf{A}}\bar{\mathbf{A}}^{\top})$ sample.
The main result is Lemma 6.6. Notice that Algorithm 4 generates new $T_{i}$’s
with given size by calling the uniform sampling oracle, thus we first compute
the number of calls to the oracle we need in Lemma 6.5. Since the proof of
Lemma 6.5 is a straightforward concentration argument, we defer the proof to
Appendix A.

######  Lemma 6.5.

Suppose we have an oracle that can uniformly sample one element from $[m]$.
Suppose we have $S\in\tbinom{[m]}{k}$ where $k\leq m$. Let $t=O(k)$ and
$0<\delta<1$, then with probability $1-\delta$, after at most
$O(t+\log\frac{1}{\delta})$ calls to the oracle, we can obtain
$T:=\cup_{j=1}^{t}\\{i_{j}\\}$ such that $T\subseteq[m]\backslash S$.

######  Lemma 6.6 (Ghost coupling).

Conditioned on an event that holds with probability $0.99$ (and only depends
on RHT). For $k\geq\log m$ and $\delta=m^{-O(1)}$, if we run Algorithm 4 with
choice $\mathrm{call}\geq Ck$ for some constant $C=O(1)$, and
$\mathrm{iter}=O(\log^{3}m)$, then with probability $1-\delta$,
$S_{\mathrm{iter}}$ is the same as $S_{\mathrm{DPP}}$, where
$S_{\mathrm{DPP}}\sim
k$-$\mathrm{DPP}(\bar{\mathbf{A}}\bar{\mathbf{A}}^{\top})$.

###### Proof.

According to Lemma 6.4 and Remark 6.2, with probability $0.99$ we have
$\Pr\\{j\in S_{\mathrm{DPP}}\\}\leq\frac{c_{0}k}{m}$ holds for all $j\in[m]$
for some constant $c_{0}=O(1)$, which gives

$\displaystyle K:=\sum_{j=1}^{m}\Pr\\{j\in S_{\mathrm{DPP}}\\}\leq c_{0}k.$

By combining this result with Lemma 6.1, we know that conditioned on an event
that holds with probability $0.99$, Algorithm 4 produces a sample from a
distribution with total variation distance $m^{-O(1)}$ from the $k$-DPP
distribution with $O(\log^{3}m)$ iterations of up-down walk, and size of the
super set needs to satisfies $|T_{i}|=O(K)$. Since we know that $K=O(k)$, the
constraints become $|T_{i}|=O(k)$ for all $i$. By Lemma 6.5, in order to
satisfy $|T_{i}|=t-k=O(k)$ for all $i$ holds with probability
$1-\frac{\delta}{2}$, we need to set

$\displaystyle\mathrm{call}=O(t-k+\log\frac{2\mathrm{iter}}{\delta})=O(k+\log\log
m+\log\frac{1}{\delta})=O(k+\log\frac{1}{\delta}).$

We set $\delta=m^{-O(1)}$ and have $\log\frac{1}{\delta}=O(\log m)$, thus
there exists a constant $C=O(1)$ such that it suffices to set
$\mathrm{call}\geq Ck$. Then we can apply a union bound and know that with
probability $1-\frac{\delta}{2}$, we can obtain all
$\\{T_{i}\\}_{i=1}^{\mathrm{iter}}$ such that $|T_{i}|=t-k$ and $T_{i}\cap
S_{i-1}=\varnothing$. Finally, by using Lemma 2.4 and applying another union
bound again, we know that
$\Pr\\{S_{\mathrm{iter}}=S_{\mathrm{DPP}}\\}=1-\delta$. ∎

#### Completing the coupling argument.

We now use the guarantee for the ghost algorithm from Lemma 6.6 to complete
the proof of Lemma 6.2 from the beginning of this section.

###### Proof of Lemma 6.2.

In Algorithm 5, in order to calculate $\tau$ which is the number of calls to
the oracle we need, we look at the coupled Algorithm 4. There are two parts
that need to be computed: (i) $\mathrm{call}$, the number of calls we need to
obtain each $T_{i}$ and (ii) $\mathrm{iter}$, the number of iterations we
need. Let $\delta=m^{-O(1)}$ be the failure probability. We set
$\mathrm{call}\geq Ck$ for some constant $C=O(1)$ and
$\mathrm{iter}=O(\log^{3}m)$, then we have
$\tau=\mathrm{call}\times\mathrm{iter}\geq Ck\log^{3}m$. According to Lemma
6.6, we know $S_{\mathrm{iter}}=S_{\mathrm{DPP}}$ holds with probability
$1-\delta$. Furthermore, since we almost surely have

$\displaystyle
S_{\mathrm{iter}}\subseteq\bigcup_{i=0}^{\mathrm{iter}-1}(S_{i}\cup
T_{i})=\tilde{S},$

thus we know $S_{\mathrm{DPP}}\subseteq\tilde{S}$ holds with probability
$1-\delta$. Thus we finish the proof. ∎

Now with Lemma 4.1, Lemma 6.2 and Lemma 6.3, we prove our main result as
follows.

###  6.3 Proof of Lemma 4.2

###### Proof of Lemma 4.2.

For $k\geq\log m$, we set $\tau\geq Ck\log^{3}m$ for some constant $C=O(1)$
which is specified in the proof of Lemma 6.6 and larger than $2$. We set
$k^{\prime}=\min\\{\tau/\log^{3}m,r\\}$ and consider $\mathrm{DPP}$ sample
$S_{\mathrm{DPP}}\sim
k^{\prime}$-$\mathrm{DPP}(\bar{\mathbf{A}}\bar{\mathbf{A}}^{\top})$, and let
$\mathbf{P}_{S_{\mathrm{DPP}}}$ be the corresponding projection matrix. Let
$\tilde{S}$ be the output of Algorithm 5 with given $\tau$, and let
$\mathbf{P}_{\tilde{S}}$ be the corresponding projection matrix. We define
$\mathcal{E}$ as the event $\\{S_{\mathrm{DPP}}\subseteq\tilde{S}\\}$, then
according to Lemma 6.2, event $\mathcal{E}$ holds with probability $1-\delta$.
Thus we have

$\displaystyle{\mathbb{E}}[\mathbf{P}_{\tilde{S}}]=$
$\displaystyle~{}{\mathbb{E}}[\mathbf{P}_{\tilde{S}}\mid\mathcal{E}]\Pr\\{\mathcal{E}\\}+{\mathbb{E}}[\mathbf{P}_{\tilde{S}}\mid\mathcal{E}^{c}]\Pr\\{\mathcal{E}^{c}\\}$
$\displaystyle\succeq$
$\displaystyle~{}{\mathbb{E}}[\mathbf{P}_{\tilde{S}}\mid\mathcal{E}]\Pr\\{\mathcal{E}\\}$
$\displaystyle\succeq$
$\displaystyle~{}{\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}\mid\mathcal{E}]\Pr\\{\mathcal{E}\\}$
$\displaystyle=$
$\displaystyle~{}{\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}]-{\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}\mid\mathcal{E}^{c}]\Pr\\{\mathcal{E}^{c}\\}$
$\displaystyle\succeq$
$\displaystyle~{}{\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}]-\Pr\\{\mathcal{E}^{c}\\}\cdot\mathbf{I}$
$\displaystyle=$
$\displaystyle~{}{\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}]-\delta\cdot\mathbf{I}$

where the first step follows from the law of total expectation, the second
step follows since the projection matrices are PSD, the third step follows
from the definition of $\mathcal{E}$ and Lemma 6.3, the fourth step follows
from the law of total expectation, the fifth step follows from Lemma 6.3 and
the last step follows since $\Pr\\{\mathcal{E}\\}\geq 1-\delta$. Notice that
since $\mathbf{Q}^{\top}\mathbf{Q}=\mathbf{I}$, we know that
$\bar{\mathbf{A}}=\mathbf{Q}\mathbf{A}$ has the same spectrum as $\mathbf{A}$
and thus
$\bar{\kappa}_{k}(\bar{\mathbf{A}})=\bar{\kappa}_{k}(\mathbf{A})=\bar{\kappa}_{k}$.
By applying Lemma 4.1 to $\bar{\mathbf{A}}$ we know that
$\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}])\geq\frac{k^{\prime}-k}{k^{\prime}-k-1+(r-k)\bar{\kappa}_{k}^{2}}$.
By combining this with the above result we have

$\displaystyle\lambda_{\min}({\mathbb{E}}[\mathbf{P}_{\tilde{S}}])\geq$
$\displaystyle~{}\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}])-\delta$
$\displaystyle\geq$
$\displaystyle~{}\frac{k^{\prime}-k}{k^{\prime}-k-1+(r-k)\bar{\kappa}_{k}^{2}}-\delta$
$\displaystyle>$
$\displaystyle~{}\frac{\tau-k\log^{3}m}{\tau-k\log^{3}m+(r-k)\bar{\kappa}_{k}^{2}\log^{3}m}-\delta$
$\displaystyle=$
$\displaystyle~{}\frac{\tau-k\log^{3}m}{\tau-(\bar{\kappa}_{k}^{2}+1)k\log^{3}m+r\bar{\kappa}_{k}^{2}\log^{3}m}-\delta$
(6.1)

where the third step follows from $k^{\prime}\leq\tau/\log^{3}m$. We can
verify that by treating $\bar{\kappa}_{k}$ as a constant and considering the
above expression (Eq.(6.3)) as a function of $k$, then this expression is
decreasing in $k$. By using $k\leq\tau/(C\log^{3}m)$, we have

$\displaystyle\lambda_{\min}({\mathbb{E}}[\mathbf{P}_{\tilde{S}}])\geq$
$\displaystyle~{}\frac{(C-1)\tau}{C\tau-(\bar{\kappa}_{k}^{2}+1)\tau+Cr\bar{\kappa}_{k}^{2}\log^{3}m}-\delta$
$\displaystyle=$
$\displaystyle~{}\frac{(C-1)\tau/r}{(C-1-\bar{\kappa}_{k}^{2})\tau/r+C\bar{\kappa}_{k}^{2}\log^{3}m}-\delta.$

We consider two cases here. If $\bar{\kappa}_{k}^{2}\geq C-1$, then by using
$C\geq 2$ we have

$\displaystyle\lambda_{\min}({\mathbb{E}}[\mathbf{P}_{\tilde{S}}])\geq\frac{(C-1)\tau/r}{C\bar{\kappa}_{k}^{2}\log^{3}m}-\delta\geq\frac{\tau/r}{2\bar{\kappa}_{k}^{2}\log^{3}m}-\delta.$

If $\bar{\kappa}_{k}^{2}<C-1$, then we use the assumption that
$\tau/\log^{3}m\leq r$ (otherwise we can simply choose $k^{\prime}=r$ which is
the trivial case) and have

$\displaystyle\lambda_{\min}({\mathbb{E}}[\mathbf{P}_{\tilde{S}}])\geq$
$\displaystyle~{}\frac{(C-1)\tau/r}{(C-1-\bar{\kappa}_{k}^{2})\log^{3}m+C\bar{\kappa}_{k}^{2}\log^{3}m}-\delta$
$\displaystyle=$
$\displaystyle~{}\frac{(C-1)\tau/r}{(C-1)\log^{3}m+(C-1)\bar{\kappa}_{k}^{2}\log^{3}m}-\delta$
$\displaystyle=$
$\displaystyle~{}\frac{\tau/r}{(\bar{\kappa}_{k}^{2}+1)\log^{3}m}-\delta\geq\frac{\tau/r}{2\bar{\kappa}_{k}^{2}\log^{3}m}-\delta$

where we use the fact that $\bar{\kappa}_{k}\geq 1$. Finally, according to
Remark 6.1 we can see that the failure probability is actually
$\delta=m^{-9}$. After some calculation we know that as long as
$\bar{\kappa}_{k}\leq m^{4}$, we have the following holds:

$\displaystyle\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{\tilde{S}}])\geq\frac{\tau/r}{2\bar{\kappa}_{k}^{2}\log^{3}m}-\frac{1}{m^{9}}\geq\frac{\tau/r}{2.5\bar{\kappa}_{k}^{2}\log^{3}m}.$

∎

###  6.4 Comparison with Direct $k$-DPP Sampling

In this section we compare our result with the results from [ALV22], and show
an improvement of at least $O(k^{\omega-2})$. The main result of approximately
sampling $k$-DPP in [ALV22] is Lemma 2.3, however our setting is different in
that we are not given the PSD matrix $\mathbf{A}\mathbf{A}^{\top}$, instead we
only know $\mathbf{A}\in\mathbb{R}^{m\times n}$. We can naively first compute
$\mathbf{A}\mathbf{A}^{\top}$ which takes $O(m^{\omega-1}n)$ then use Lemma
2.3, resulting in prohibitive $\tilde{O}(m^{\omega-1}n+mk^{\omega-1})$ cost of
sampling one $k$-DPP. As an alternative, we can improve this result by not
computing $\mathbf{A}\mathbf{A}^{\top}$ explicitly.

By slightly revising the proof in [ALV22], we show that given
$\mathbf{A}\in\mathbb{R}^{m\times n}$, we can obtain $s$ independent
approximate samples from $k$-$\mathrm{DPP}(\mathbf{A}\mathbf{A}^{\top})$ in
time $\tilde{O}(mnk^{\omega-2}+snk^{\omega-1})$. To be specific there are two
steps, here we give a brief proof sketch of each step without rigorous proof.

  * •

Preprocessing: according to [ALV22], we first need to obtain an exact $k$-DPP
sample defined on an $8k\times 8k$ PSD matrix. Constructing this matrix takes
$O(nk^{\omega-1})$, and sampling $k$-DPP defined by it takes
$\tilde{O}(k^{\omega})$, thus it takes
$\tilde{O}(nk^{\omega-1}+k^{\omega})=\tilde{O}(nk^{\omega-1})$. Then by using
Theorem 34 in [ALV22], Algorithm 1 produces marginal overestimates
$q_{i}\geq\Pr\\{i\in S_{\mathrm{DPP}}\\}$ with sum $\sum_{i\in[m]}q_{i}\leq
4k$ in time $\tilde{O}(\frac{m}{k}\cdot
nk^{\omega-1})=\tilde{O}(mnk^{\omega-2})$. After the preprocessing step, the
distribution is in approximately isotropic position.

  * •

Sampling: after the preprocessing step, we obtain an approximately isotropic
distribution, and the next step is to get independent approximate samples.
This step is relatively easy, we can directly use Lemma 6.1 and know that we
can obtain approximate samples from
$k$-$\mathrm{DPP}(\mathbf{A}\mathbf{A}^{\top})$ in time
$\tilde{O}(1)\cdot\tilde{O}(nk^{\omega-1})=\tilde{O}(nk^{\omega-1})$. Here we
use the fact that
$\mathcal{T}_{\mu}(4k,k)=\tilde{O}(n(4k)^{\omega-1}+(4k)^{\omega})=\tilde{O}(nk^{\omega-1})$
with $\mu=k$-$\mathrm{DPP}(\mathbf{A}\mathbf{A}^{\top})$ is the time of
sampling exact DPP.

By the above two steps, we conclude that we can obtain independent approximate
$k$-DPP samples in time $\tilde{O}(mnk^{\omega-2}+snk^{\omega-1})$. In
comparison, our Lemma 6.2 only needs an $\tilde{O}(mn)$ preprocessing time and
$O(k\log^{3}m)$ calls to an uniform sampling oracle to guarantee the same
failure probability, thus the total cost is only
$\tilde{O}(mn+sk\log^{3}m)=\tilde{O}(mn+sk)$ which is at least an
$O(k^{\omega-2})$ improvement.

##  7 Approximate Sketch-and-Project

In this section we detail the idea and discussion in Section 4.3, and show how
we approximate the update step of sketch-and-project by solving a
preconditioned normal equation as a sub-problem. For notation convenience we
use $\mathbf{A}\in\mathbb{R}^{m\times n}$ and $\mathbf{b}\in\mathbb{R}^{m}$ to
denote the ones after randomized Hadamard transform, which is the original
$\mathbf{Q}\mathbf{A}$ and $\mathbf{Q}\mathbf{b}$ respectively for
$\mathbf{Q}=\frac{1}{\sqrt{m}}\mathbf{H}\mathbf{D}$. Specifically,
$\mathbf{A}$ in this section corresponds to $\bar{\mathbf{A}}$ in Section 6.

###  7.1 Approximate Update by Solving Sub-Problem

In the update rule of sketch-and-project (see Eq.(4.1)), we need to compute
$\mathbf{w}_{t}^{*}:=(\mathbf{I}_{\tilde{S}}\mathbf{A})^{\dagger}(\mathbf{I}_{\tilde{S}}\mathbf{A}\mathbf{x}_{t}-\mathbf{I}_{\tilde{S}}\mathbf{b})$
in each iteration. We can directly compute this quantity: we first compute
$\mathbf{I}_{\tilde{S}}\mathbf{A}$ in $O(n\tau)$, then compute
$\mathbf{I}_{\tilde{S}}\mathbf{A}\cdot(\mathbf{I}_{\tilde{S}}\mathbf{A})^{\top}$
in $O(n\tau^{\omega-1})$, then compute
$(\mathbf{I}_{\tilde{S}}\mathbf{A}\mathbf{A}^{\top}\mathbf{I}_{\tilde{S}}^{\top})^{\dagger}$
in $O(\tau^{\omega})$, and finally apply $\mathbf{I}_{\tilde{S}}\mathbf{A}$
then
$(\mathbf{I}_{\tilde{S}}\mathbf{A}\mathbf{A}^{\top}\mathbf{I}_{\tilde{S}}^{\top})^{\dagger}$
to a vector, which takes $O(n\tau+\tau^{2})$ time. These steps take
$O(n\tau^{\omega-1}+\tau^{\omega})$ in total, which is prohibitive for our
per-iteration cost. Instead of directly computing $\mathbf{w}_{t}^{*}$, we
show that it is possible to _approximately_ compute this quantity in
$\tilde{O}(n\tau+\tau^{\omega})$ time, by relying on the fact that it can be
transformed into a highly under-determined least squares problem, and using a
variant of the sketch-to-precondition approach [RT08, MSM14].

Denote
$\tilde{\mathbf{A}}=\mathbf{I}_{\tilde{S}}\mathbf{A}\in\mathbb{R}^{\tau\times
n}$ and
$\tilde{\mathbf{b}}_{t}=\mathbf{I}_{\tilde{S}}\mathbf{A}\mathbf{x}_{t}-\mathbf{I}_{\tilde{S}}\mathbf{b}$,
notice that since $\tau\leq n$, we can transform the problem of computing
$(\mathbf{I}_{\tilde{S}}\mathbf{A})^{\dagger}(\mathbf{I}_{\tilde{S}}\mathbf{A}\mathbf{x}_{t}-\mathbf{I}_{\tilde{S}}\mathbf{b})$
to solving an under-determined least squares problem:

$\displaystyle\mathbf{w}_{t}^{*}=\arg\min\|\mathbf{w}\|_{2}~{}~{}~{}\text{subject
to}~{}~{}~{}\tilde{\mathbf{A}}\mathbf{w}=\tilde{\mathbf{b}}_{t}$ (7.1)

which has solution
$\mathbf{w}_{t}^{*}=\tilde{\mathbf{A}}^{\dagger}\tilde{\mathbf{b}}_{t}$. It is
equivalent to solving the second kind normal equation:

$\displaystyle\tilde{\mathbf{A}}\tilde{\mathbf{A}}^{\top}\mathbf{u}=\tilde{\mathbf{b}}_{t},~{}~{}~{}\mathbf{w}=\tilde{\mathbf{A}}^{\top}\mathbf{u}.$

Algorithm 6 Construct preconditioner for solving normal equation.

1:Input: matrix $\tilde{\mathbf{A}}\in\mathbb{R}^{\tau\times n}$;

2:Compute $\mathbf{\Phi}\tilde{\mathbf{A}}^{\top}$ by calling Lemma 2.7;
$\triangleright$ Sketch: takes $\tilde{O}(n\tau)$ time.

3:Compute compact size SVD as
$\mathbf{\Phi}\tilde{\mathbf{A}}^{\top}=\tilde{\mathbf{U}}\tilde{\mathbf{\Sigma}}\tilde{\mathbf{V}}^{\top}$;
$\triangleright$ Takes $O(\tau^{\omega})$ time.

4:Construct preconditioner
$\mathbf{M}\leftarrow\tilde{\mathbf{V}}\tilde{\mathbf{\Sigma}}^{-1}$;

5:return $\mathbf{M}$;

Notice that
$\tilde{\mathbf{A}}\tilde{\mathbf{A}}^{\top}\mathbf{u}=\tilde{\mathbf{b}}_{t}$
is a positive semidefinite linear system over $\mathbf{u}$, thus we can use CG
to solve it. However since we cannot bound $\kappa(\tilde{\mathbf{A}})$, we
need to precondition this linear system first. Notice that the subspace
embedding result of Lemma 2.7 directly implies a way to construct a
preconditioner. By calling Algorithm 6 and treating $\epsilon$ in Lemma 2.7 as
a fixed constant, we can obtain an invertible matrix
$\mathbf{M}\in\mathbb{R}^{\tau\times\tau}$ such that
$\tilde{\kappa}:=\kappa(\tilde{\mathbf{A}}^{\top}\mathbf{M})=O(1)$. Thus we
can precondition Problem(7.1) as

$\displaystyle\tilde{\mathbf{w}}^{*}=\arg\min\|\mathbf{w}\|_{2}~{}~{}~{}\text{subject
to}~{}~{}~{}\mathbf{M}^{\top}\tilde{\mathbf{A}}\mathbf{w}=\mathbf{M}^{\top}\tilde{\mathbf{b}}_{t}$
(7.2)

and similarly, we can reformulate Problem(7.2) to the following second kind
normal equation:

$\displaystyle(\mathbf{M}^{\top}\tilde{\mathbf{A}}\tilde{\mathbf{A}}^{\top}\mathbf{M})\mathbf{u}=\mathbf{M}^{\top}\tilde{\mathbf{b}}_{t},~{}~{}~{}\mathbf{w}=\tilde{\mathbf{A}}^{\top}\mathbf{M}\mathbf{u}.$

We use CG to solve this positive semidefinite linear system over $\mathbf{u}$,
and have the following lemma.

######  Lemma 7.1 (Convergence).

Given $\tilde{\mathbf{A}}\in\mathbb{R}^{\tau\times
n},\tilde{\mathbf{b}}\in\mathbb{R}^{\tau}$, let
$\mathbf{w}^{*}=\tilde{\mathbf{A}}^{\dagger}\tilde{\mathbf{b}}$ be the least
squares solution of $\tilde{\mathbf{A}}\mathbf{w}=\tilde{\mathbf{b}}$. Given
$\delta<1/2$, if the subspace embedding matrix
$\mathbf{\Phi}\in\mathbb{R}^{\phi\times n}$ from Lemma 2.7 satisfies
$\phi=O(\tau+\log(1/\delta))$, then with probability $1-\delta$, Algorithm 6
outputs $\mathbf{M}$ in time $O(n\tau\log(\tau/\delta)+\tau^{\omega})$ such
that $\tilde{\kappa}=\kappa(\tilde{\mathbf{A}}^{\top}\mathbf{M})=O(1)$, and
solving the following system

$\displaystyle(\mathbf{M}^{\top}\tilde{\mathbf{A}}\tilde{\mathbf{A}}^{\top}\mathbf{M})\mathbf{u}=\mathbf{M}^{\top}\tilde{\mathbf{b}},~{}~{}~{}\mathbf{w}=\tilde{\mathbf{A}}^{\top}\mathbf{M}\mathbf{u}$
(7.3)

by using CG and starting from $\mathbf{w}_{0}=\mathbf{0}$ will yield the
following convergence result:

$\displaystyle\|\mathbf{w}_{s}-\mathbf{w}^{*}\|\leq
2\left(\frac{\tilde{\kappa}-1}{\tilde{\kappa}+1}\right)^{s}\cdot\|\mathbf{w}^{*}\|.$

###### Proof.

Let $\mathbf{u}^{*}$ be the min-length solution of the PSD linear system
$(\mathbf{M}^{\top}\tilde{\mathbf{A}}\tilde{\mathbf{A}}^{\top}\mathbf{M})\mathbf{u}=\mathbf{M}^{\top}\tilde{\mathbf{b}}$.
Denote $\tilde{\mathbf{w}}^{*}$ as the solution of Problem(7.3), then we have

$\displaystyle\tilde{\mathbf{w}}^{*}=$
$\displaystyle~{}\tilde{\mathbf{A}}^{\top}\mathbf{M}\mathbf{u}^{*}$
$\displaystyle=$
$\displaystyle~{}\tilde{\mathbf{A}}^{\top}\mathbf{M}(\mathbf{M}^{\top}\tilde{\mathbf{A}}\tilde{\mathbf{A}}^{\top}\mathbf{M})^{\dagger}\mathbf{M}^{\top}\tilde{\mathbf{b}}$
$\displaystyle=$
$\displaystyle~{}\tilde{\mathbf{A}}^{\top}\mathbf{M}\mathbf{M}^{-1}(\tilde{\mathbf{A}}\tilde{\mathbf{A}}^{\top})^{\dagger}(\mathbf{M}^{\top})^{-1}\mathbf{M}^{\top}\tilde{\mathbf{b}}$
$\displaystyle=$
$\displaystyle~{}\tilde{\mathbf{A}}^{\top}(\tilde{\mathbf{A}}\tilde{\mathbf{A}}^{\top})^{\dagger}\tilde{\mathbf{b}}$
$\displaystyle=$
$\displaystyle~{}\tilde{\mathbf{A}}^{\dagger}\tilde{\mathbf{b}}=\mathbf{w}^{*}$

which shows that the optimal solution of Problem(7.3) is the same as that of
$\tilde{\mathbf{A}}\mathbf{w}=\tilde{\mathbf{b}}$. According to Lemma 2.7 we
know that with the choice of $\phi=O(\tau+\log(1/\delta))$,
$\tilde{\kappa}=O(1)$ holds with probability $1-\delta$. According to Lemma
2.6, the convergence guarantee is

$\displaystyle\|\tilde{\mathbf{A}}^{\top}\mathbf{M}(\mathbf{u}_{s}-\mathbf{u}^{*})\|\leq
2\left(\frac{\tilde{\kappa}-1}{\tilde{\kappa}+1}\right)^{s}\cdot\|\tilde{\mathbf{A}}^{\top}\mathbf{M}\mathbf{u}^{*}\|$

and by doing the change of variable back from $\mathbf{u}$ to $\mathbf{w}$ we
have

$\displaystyle\|\mathbf{w}_{s}-\mathbf{w}^{*}\|\leq
2\left(\frac{\tilde{\kappa}-1}{\tilde{\kappa}+1}\right)^{s}\cdot\|\mathbf{w}^{*}\|.$

Finally we compute the cost of constructing the preconditioner. Following
Lemma 2.7 and Algorithm 6, we first construct
$\mathbf{\Phi}\tilde{\mathbf{A}}^{\top}$ which takes
$O(n\tau\log(\tau/\delta)+\tau^{2}\log^{4}(\tau/\delta))$, then compute the
compact size SVD as
$\mathbf{\Phi}\tilde{\mathbf{A}}^{\top}=\tilde{\mathbf{U}}\tilde{\mathbf{\Sigma}}\tilde{\mathbf{V}}^{\top}$
which takes $O(\phi\tau^{\omega-1})$ where $\phi=O(\tau+\log(1/\delta))$, and
finally we set $\mathbf{M}=\mathbf{V}\tilde{\mathbf{\Sigma}}^{-1}$ which takes
$O(\tau^{2})$. Thus we conclude that the total time complexity is

$\displaystyle
O(n\tau\log(\tau/\delta)+\tau^{2}\log^{4}(\tau/\delta)+\tau^{\omega-1}\phi)=O(n\tau\log(\tau/\delta)+\tau^{\omega})=\tilde{O}(n\tau+\tau^{\omega}).$

∎

###  7.2 Proof of Lemma 4.3

We are now ready to provide the convergence analysis for sketch-and-project
with an approximate update, as it is implemented in Algorithm 2.

###### Proof of Lemma 4.3.

In outer iteration $t$, we analyze the update as follows.

$\displaystyle{\mathbb{E}}\|\mathbf{x}_{t+1}-\mathbf{x}^{*}\|^{2}=$
$\displaystyle~{}{\mathbb{E}}\|\mathbf{x}_{t}-\mathbf{w}_{s_{\max}}-\mathbf{x}^{*}\|^{2}$
$\displaystyle=$
$\displaystyle~{}{\mathbb{E}}\|(\mathbf{x}_{t}-\mathbf{x}^{*}-\mathbf{w}_{t}^{*})+(\mathbf{w}_{t}^{*}-\mathbf{w}_{s_{\max}})\|^{2}$
$\displaystyle=$
$\displaystyle~{}{\mathbb{E}}\|\mathbf{x}_{t}-\mathbf{x}^{*}-\mathbf{w}_{t}^{*}\|^{2}+\|\mathbf{w}_{t}^{*}-\mathbf{w}_{s_{\max}}\|^{2}+2{\mathbb{E}}[(\mathbf{x}_{t}-\mathbf{x}^{*}-\mathbf{w}_{t}^{*})^{\top}(\mathbf{w}_{t}^{*}-\mathbf{w}_{s_{\max}})]$
$\displaystyle\leq$
$\displaystyle~{}{\mathbb{E}}\|\mathbf{x}_{t}-\mathbf{x}^{*}-\mathbf{w}_{t}^{*}\|^{2}+\|\mathbf{w}_{t}^{*}-\mathbf{w}_{s_{\max}}\|^{2}+2\|\mathbf{w}_{t}^{*}-\mathbf{w}_{s_{\max}}\|\cdot({\mathbb{E}}\|\mathbf{x}_{t}-\mathbf{x}^{*}-\mathbf{w}_{t}^{*}\|^{2})^{1/2}.$

For the first part, as we already discussed in Section 4.1, we have

$\displaystyle{\mathbb{E}}\|\mathbf{x}_{t}-\mathbf{x}^{*}-\mathbf{w}_{t}^{*}\|^{2}=$
$\displaystyle~{}(\mathbf{x}_{t}-\mathbf{x}^{*})^{\top}{\mathbb{E}}[(\mathbf{I}-\mathbf{P}_{{\tilde{S}}})](\mathbf{x}_{t}-\mathbf{x}^{*})$
$\displaystyle\leq$
$\displaystyle~{}(1-\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{{\tilde{S}}}]))\cdot\|\mathbf{x}_{t}-\mathbf{x}^{*}\|^{2}.$
(7.4)

This is due to the update rule of sketch-and-project:

$\displaystyle\mathbf{x}_{t+1}=\mathbf{x}_{t}-\mathbf{w}_{t}^{*}=\mathbf{x}_{t}-\mathbf{A}^{\top}\mathbf{I}_{\tilde{S}}^{\top}(\mathbf{I}_{\tilde{S}}\mathbf{A}\mathbf{A}^{\top}\mathbf{I}_{\tilde{S}}^{\top})^{\dagger}(\mathbf{I}_{\tilde{S}}\mathbf{A}\mathbf{x}_{t}-\mathbf{I}_{\tilde{S}}\mathbf{b})$

and the fact that if $\mathbf{x}_{t}-\mathbf{x}^{*}$ is not in the null-space
of ${\mathbb{E}}[\mathbf{P}_{\tilde{S}}]$, then neither will be
$\mathbf{x}_{t+1}-\mathbf{x}^{*}$. This fact follows because

$\displaystyle{\mathbb{E}}[\mathbf{P}_{\tilde{S}}]=$
$\displaystyle~{}\mathbf{A}^{\top}{\mathbb{E}}[\mathbf{I}_{\tilde{S}}^{\top}(\mathbf{I}_{\tilde{S}}\mathbf{A}\mathbf{A}^{\top}\mathbf{I}_{\tilde{S}}^{\top})^{\dagger}\mathbf{I}_{\tilde{S}}]\mathbf{A}$

which implies that the nullspace of ${\mathbb{E}}[\mathbf{P}_{\tilde{S}}]$
contains the nullspace of $\mathbf{A}^{\top}\mathbf{A}$. At the same time, in
the proof of Lemma 4.2 we showed that
${\mathbb{E}}[\mathbf{P}_{\tilde{S}}]\succeq{\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}]-\delta\mathbf{I}\succeq\frac{1}{2}{\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}]$,
where $S_{\mathrm{DPP}}\sim k^{\prime}$-DPP$(\mathbf{A}\mathbf{A}^{\top})$, so
the nullspace of ${\mathbb{E}}[\mathbf{P}_{\tilde{S}}]$ is no larger than that
of ${\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}]$. But we know that the
nullspace of ${\mathbb{E}}[\mathbf{P}_{S_{\mathrm{DPP}}}]$ is exactly the same
as that of $\mathbf{A}^{\top}\mathbf{A}$ using Eq.(5.3). Finally, denoting the
singular value decomposition of $\mathbf{A}$ as
$\mathbf{A}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$ where
$\mathbf{U}\in\mathbb{R}^{m\times m},\mathbf{\Sigma}\in\mathbb{R}^{m\times n}$
and $\mathbf{V}\in\mathbb{R}^{n\times n}$, the update rule can be seen as
$\mathbf{x}_{t+1}=\mathbf{x}_{t}-\mathbf{V}\cdot\mathbf{v}_{t}$ for some
vector $\mathbf{v}_{t}$, thus if $\mathbf{x}_{t}-\mathbf{x}^{*}$ is not in the
null-space of ${\mathbb{E}}[\mathbf{P}_{\tilde{S}}]$, then
$\mathbf{x}_{t+1}-\mathbf{x}^{*}$ will never be in the null-space of
${\mathbb{E}}[\mathbf{P}_{\tilde{S}}]$. In particular, if we set
$\mathbf{x}_{0}=\mathbf{0}$, then notice that
$\mathbf{x}^{*}=\mathbf{A}^{\dagger}\mathbf{b}=\mathbf{V}\mathbf{\Sigma}^{\dagger}\mathbf{U}^{\top}\mathbf{b}$,
we know $\mathbf{x}_{0}-\mathbf{x}^{*}=\mathbf{x}^{*}$ is not in the nullspace
of ${\mathbb{E}}[\mathbf{P}_{\tilde{S}}]$. Thus we know Eq.(7.2) holds.

For the second and third part, denote
$\Delta_{t}:=\|\mathbf{w}_{t}^{*}-\mathbf{w}_{s_{\max}}\|$ and
$\rho:=\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{\tilde{S}}])$. Then, we
have

$\displaystyle{\mathbb{E}}\|\mathbf{x}_{t+1}-\mathbf{x}^{*}\|^{2}\leq$
$\displaystyle~{}(1-\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{\tilde{S}}]))\cdot\|\mathbf{x}_{t}-\mathbf{x}^{*}\|^{2}+2\Delta_{t}\sqrt{1-\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{\tilde{S}}])}\cdot\|\mathbf{x}_{t}-\mathbf{x}^{*}\|+\Delta_{t}^{2}$
$\displaystyle=$
$\displaystyle~{}(\sqrt{1-\lambda_{\min}^{+}({\mathbb{E}}[\mathbf{P}_{\tilde{S}}])}\cdot\|\mathbf{x}_{t}-\mathbf{x}^{*}\|+\Delta_{t})^{2}$
$\displaystyle=$
$\displaystyle~{}(\sqrt{1-\rho}\cdot\|\mathbf{x}_{t}-\mathbf{x}^{*}\|+\Delta_{t})^{2}.$
(7.5)

We define $\mathcal{E}$ as the event that
$\kappa(\tilde{\mathbf{A}}^{\top}\mathbf{M})=O(1)$, then according to Lemma
2.7, as long as we set $\phi=O(\tau+\log(1/\delta_{1}))$, then event