Determinantal Point Processes in Randomized Numerical Linear Algebra

Introduction

Randomized Numerical Linear Algebra (RandNLA), is an area which uses randomness, most notably random sampling and random projection methods, to develop improved algorithms for ubiquitous matrix problems. It began as a niche area in theoretical computer science about fifteen years ago , and since then the area has exploded. Matrix problems are central to much of applied mathematics, from traditional scientific computing and partial differential equations to statistics, machine learning, and artificial intelligence. Generalizations and variants of matrix problems are central to many other areas of mathematics, via more general transformations and algebraic structures, nonlinear optimization, infinite-dimensional operators, etc. Much of the work in RandNLA has been propelled by recent developments in machine learning, artificial intelligence, and large-scale data science, and RandNLA both draws upon and contributes back to both pure and applied mathematics.

A seemingly different topic, but one which has a long history in pure and applied mathematics, is that of Determinantal Point Processes (DPPs). A DPP is a stochastic point process, the probability distribution of which is characterized by sub-determinants of some matrix. Such processes were first studied to model the distribution of fermions at thermal equilibrium . In the context of random matrix theory, DPPs emerged as the eigenvalue distribution for standard random matrix ensembles , and they are of interest in other areas of mathematics such as graph theory, combinatorics and quantum mechanics . More recently, DPPs have also attracted significant attention within machine learning and statistics as a tractable probabilistic model that is able to capture a balance between quality and diversity within data sets and that admits efficient algorithms for sampling, marginalization, conditioning, etc. . This resulted in practical application of DPPs in experimental design , recommendation systems , stochastic optimization and more.

Until very recently, DPPs have had little if any presence within RandNLA. However, recent work has uncovered deep connections between these two topics. The purpose of this article is to provide an overview of RandNLA, with an emphasis on discussing and highlighting these connections with DPPs. In particular, we will show how random sampling with a DPP leads to new kinds of unbiased estimators for the classical RandNLA task of least squares regression, enabling a more refined statistical and inferential understanding of RandNLA algorithms. We will also demonstrate that a DPP is, in some sense, an optimal randomized method for low-rank approximation, another ubiquitous matrix problem. Finally, we also discuss how a standard RandNLA technique, called leverage score sampling, can be derived as the marginal distribution of a DPP, as well as the algorithmic consequences this has for efficient DPP sampling.

We start (in Section 2) with a brief review of a prototypical RandNLA algorithm, focusing on the ubiquitous least squares problem and highlighting key aspects that will put in context the recent work we will review. In particular, we discuss the trade-offs between standard sampling methods from RandNLA, including uniform sampling, norm-squared sampling, and leverage score sampling. Next (in Section 3), we introduce the family of DPPs, highlighting some important sub-classes and the basic properties that make them appealing for RandNLA. Then (in Section 4), we describe the fundamental connections between certain classes of DPPs and the classical RandNLA tasks of least squares regression and low-rank approximation, as well as the relationship between DPPs and the RandNLA method of leverage score sampling. Finally (in Section 5), we discuss the algorithmic aspects of both leverage scores and DPPs. We conclude (in Section 6) by briefly mentioning several other connections between DPPs and RandNLA, as well as a recently introduced class of random matrices, called determinant preserving, which has proven useful in this line of research.

RandNLA: Randomized Numerical Linear Algebra

Many different approaches have been established for randomly down-sizing data matrices $\mathbf{X}$ (see for a detailed survey). While some methods randomly zero out most of the entries of the matrix, most randomly keep only a small random subset of rows and/or columns. In either case, however, the choice of randomness is crucial in preserving the structure of the data. For example, if the data matrix $\mathbf{X}$ contains a few dominant entries/rows/columns (e.g., as measured by their absolute value or norm or some other “importance” score), then we should make sure that our sketch is likely to retain the information they carry. This leads to data-dependent sampling distributions that will be the focus of our discussion. However, data-independent sketching techniques, which involve applying a random transformation $\mathbf{S}$ to the matrix $\mathbf{X}$ (typically called a “random rotation” or a “random projection,” even if it is not precisely a rotation or projection in the linear algebraic sense), have also proven very successful . Among the most common examples of such random transformations are i.i.d. Gaussian matrices, fast Johnson-Lindenstraus transforms and count sketches , all of which provide different trade-offs between efficiency and accuracy. These “data-oblivious random projections” can be interpreted either in terms of the Johnson-Lindenstraus lemma or as a preconditioner for the “data aware random sampling” methods we discuss .

Most RandNLA techniques can be divided into one of two settings, depending on the dimensionality or aspect ratio of $\mathbf{X}$, and on the desired size of the sketch (see Figure 1):

Rank-preserving sketch. When $\mathbf{X}$ is a tall full-rank matrix (i.e., $n\gg d$), then we can reduce the larger dimension while preserving the rank.

Low-rank approximation. When $\mathbf{X}$ has comparably large dimensions (i.e., $n\sim d$), then the sketch typically has a much lower rank than $\mathbf{X}$.

The least squares solution can be computed exactly, using the Moore-Penrose pseudoinverse, $\mathbf{w}^{*}=\mathbf{X}^{\dagger}\mathbf{y}$. In a traditional RandNLA setup, in order to avoid solving the full problem, our goal is to use the sketch-and-solve paradigm to obtain an $(\epsilon,\delta)$-approximation of $\mathbf{w}^{*}$, i.e., $\widehat{\mathbf{w}}$ such that:

Imposing statistical modeling assumptions on the vector $\mathbf{y}$ leads to different objectives, such as the mean squared error (MSE):

where $\boldsymbol{\xi}$ is a noise vector with a known distribution. There has been work on statistical aspects of RandNLA methods , and these statistical objectives pose different challenges than the standard RandNLA guarantees (some of which can be addressed by DPPs, see Section 4).

To illustrate the types of guarantees achieved by RandNLA methods on the least squares task, we will focus on row sampling, i.e., sketches consisting of a small random subset of the rows of $\mathbf{X}$, in the case that $n\gg d$. Concretely, the considered meta-strategy is to draw random i.i.d. row indices $j_{1},...,j_{k}$ from $\{1,...,n\}$, with each index distributed according to $(p_{1},...,p_{n})$, and then solve the subproblem formed from those indices:

where \widetilde{\mathbf{x}}_{i}^{\scriptscriptstyle{\top}}=\sqrt{{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}|}}\!\!\!\smash{\text{\tfrac{1}{kp_{j_{i}}}}}}\,\mathbf{x}_{j_{i}}^{\scriptscriptstyle{\top}} and \widetilde{y}_{i}=\sqrt{{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}|}}\!\!\!\smash{\text{\tfrac{1}{kp_{j_{i}}}}}}\,y_{j_{i}}, for $i=1,...,k$, denote the $i^{th}$ row of $\widetilde{\mathbf{X}}$ and entry of $\widetilde{\mathbf{y}}$, respectively. The rescaling is introduced to account for the biases caused by non-uniform sampling. We consider the following standard sampling distributions:

Squared norms: $p_{i}=\|\mathbf{x}_{i}\|^{2}/\|\mathbf{X}\|_{F}^{2}$, where $\|\cdot\|_{F}$ is the Frobenius (Hilbert-Schmidt) norm.

Leverage scores: $p_{i}=l_{i}/d$, where we let $l_{i}=\|\mathbf{x}_{i}\|^{2}_{(\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X})^{-1}}$ denote the $i^{th}$ leverage score of $\mathbf{X}$ and $\|\mathbf{v}\|_{\mathbf{A}}=\sqrt{\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{v}}$ is the Mahalanobis norm.

Both squared norm and leverage score sampling are standard RandNLA techniques used in a variety of applications . The following theorem (which, for convenience, we state with a failure probability of $\delta=0.1$) puts together the results that allow us to compare each row sampling distribution in the context of least squares.

Estimator $\widehat{\mathbf{w}}$ constructed as in (2) is an $(\epsilon,0.1)$-approximation, as in (1), if:

$k=\Omega(\mu d\log d+\mu d/\epsilon)$ for Uniform, where $\mu$ denotes the matrix coherence of $\mathbf{X}$.

$k=\Omega(\kappa d\log d+\kappa d/\epsilon)$ for Squared norms, where $\kappa$ is the condition number of $\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X}$.

$k=\Omega(d\log d+d/\epsilon)$ for Leverage scores.

Recall that (when considering least squares) we typically assume that $n\gg d$, so any of the three sample sizes $k$ may be much smaller than $n$. Thus, each sampling method offers a potentially useful guarantee for the number of rows needed to achieve a $(1+\epsilon)$-approximation. However, in the case of both uniform and squared norm sampling, the sample size depends not only on the dimension $d$, but also on other data-dependent quantities. For uniform sampling, that quantity is matrix coherence $\mu$, which measures the degree of non-uniformity among the data points in terms of their leverage scores: $\mu=\max_{i}\frac{n}{d}l_{i}$. For squared norm sampling, that quantity is the condition number $\kappa$, which is the ratio between the largest and the smallest eigenvalue of the $d\times d$ data covariance $\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X}$. To address the dependence on the condition number, one can replace the standard Euclidean norm with a Mahalanobis norm that takes into account the data covariance, obtaining the leverage scores.

We now briefly discuss a key structural property, called the subspace embedding, which is needed to show the guarantees of Theorem 1. This important property—first introduced into RandNLA for data-aware random sampling by and then for data-oblivious random projection by —is ubiquitous in the analysis of many RandNLA techniques. Remarkably, most DPP results do not rely on subspace embedding techniques, which is an important differentiating factor for this class of sampling distributions.

A sketching matrix $\mathbf{S}$ is a $(1\pm\epsilon)$ subspace embedding for the column space of $\mathbf{X}$ if:

The matrix $\widetilde{\mathbf{X}}$ used in (2) for constructing $\widehat{\mathbf{w}}$ can be written as $\widetilde{\mathbf{X}}=\mathbf{S}\mathbf{X}$, by letting the $i^{th}$ row of $\mathbf{S}$ be the scaled standard basis vector \sqrt{{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}|}}\!\!\!\smash{\text{\tfrac{1}{kp_{j_{i}}}}}}\,\mathbf{e}_{j_{i}}. Relying on established measure concentration results for random matrices (see, e.g., ), we can show that $\mathbf{S}$ is a subspace embedding (up to some failure probability) for each of the i.i.d sampling methods from Theorem 1. However, only leverage score sampling (or random projections, which precondition the input to have approximately uniform leverage scores ) achieves this for $O(d\log d)$ samples, independent of any input-specific quantities such as the coherence of condition number.

The least squares task formulated in (1), as well as the subspace embedding condition, require using rank-preserving sketches. However, these ideas can be naturally extended to the task of low-rank approximation. In particular, as discussed in more detail in Section 4.3, the standard error metric for low-rank sketches can be reduced to the least squares error defined on a subspace with the same rank as the sketch . Similarly, leverage score sampling has been extended to adapt to the low-rank setting. In Section 4.4, we discuss one of these extensions, called ridge leverage scores, and its connections to DPPs.

In the following sections, we show how non-i.i.d. sampling via DPPs goes beyond the standard RandNLA analysis. Among other things, this will permit us to obtain approximation guarantees with fewer than $d\log d$ samples and without a failure probability.

DPPs: Determinantal Point Processes

In this section, we define DPPs and related families of distributions (see Figure 2 for a diagram), including some basic properties and intuitions. A detailed introduction to DPPs can be found in . For a more general treatment that includes sampling from continuous domains, see .

Negative correlation: if $i\neq j$ and $\mathbf{K}_{ij}\neq 0$, then $\Pr(i\in S\mid j\in S)<\Pr(i\in S)$.

In the context of linear algebra, a slightly more restrictive definition of DPPs has proven useful.

Namely, this sampling probability is proportional to the squared $|S|$-dimensional volume of the parallelepiped spanned by the rows of $\mathbf{X}$ indexed by $S$. This immediately implies that the size of $S$ will never exceed the rank of $\mathbf{X}$ (which is bounded by $d$). Furthermore, such a distribution ensures that the set of sampled rows will be non-degenerate: no row can be obtained as a linear combination of the others. Intuitively, this property is desirable for RandNLA sampling as it avoids redundancies. Also, all else being equal, rows with larger norms are generally preferred as they contribute more to the volume.

While the subset size of a DPP is in most cases a random variable, it is easy to constrain the cardinality to some fixed value $k$. The resulting distribution is not a DPP, in the sense of Definition 2, but it retains many useful properties of proper DPPs.

While a $k$-DPP is not, in general, a DPP in the sense of Definition 2, both families belong to a broader class of negatively correlated distributions called Strongly Rayleigh (SR) measures, which were introduced by . See Figure 2.

At the intersection of DPPs and $k$-DPPs lies a family of distributions called Projection DPPs. This family is of particular importance to RandNLA.

DPPs in RandNLA

In this section, we demonstrate the fundamental connections between DPPs and standard RandNLA tasks, as well as the new kinds of RandNLA guarantees that can be achieved via these connections. Our discussion focuses on two types of DPP-based sketches (that were illustrated in Figure 1):

Projection DPPs as a rank-preserving sketch;

L-ensemble DPPs as a low-rank approximation.

These sketches can be efficiently constructed using DPP sampling algorithms which we discuss later (in Section 5). We also discuss the close relationship between DPPs and the RandNLA method of leverage score sampling, shedding light on why these two different randomized techniques have proven effective in RandNLA. For the summary, see Table 1.

We now define the least squares estimators that naturally arise from row sampling with Projection DPPs and L-ensembles. The definitions are motivated by the fact that the estimators are unbiased, relative to the solutions of the full least squares problems. Importantly, this property is not shared by i.i.d. row sampling methods used in RandNLA.

This seemingly simple identity relies on the negative correlations between the samples in a DPP, and thus cannot hold for any i.i.d. row sampling method.

Theorem 2 has an analogue in the context of low-rank approximation, where both the dimensions of $\mathbf{X}$ are comparably large (i.e., $n\sim d$), and so the desired sample size is typically much smaller than $d$. When the selected subproblem $(\mathbf{X}_{S},\mathbf{y}_{S})$ has fewer than $d$ rows (i.e., it is under-determined) then it has multiple exact solutions. A standard way to address this is picking the solution with smallest Euclidean norm, defined via the Moore-Penrose inverse: $\widehat{\mathbf{w}}=\mathbf{X}_{S}^{\dagger}\mathbf{y}_{S}$. To sample the under-determined subproblem, we use a scaled L-ensemble DPP with the expected sample size controlled by a parameter $\lambda>0$ .

Thus, the minimum norm solution of the under-determined subproblem is an unbiased estimator of the Tikhonov-regularized least squares problem, i.e., ridge regression. Ridge regression is a natural extension of the standard least squares task, particularly useful when $n\sim d$ or $n\ll d$.

2 Exact error analysis

The error analysis for DPP sampling differs significantly from the standard RandNLA techniques discussed in Section 2. In particular, approximation guarantees are formulated in terms of the expected error, without relying on measure concentration results. This means that we avoid failure probabilities such as the one present in Theorem 1, and the analysis is often much more precise, sometimes even exact. Furthermore, because of the non-i.i.d. nature of DPPs, these guarantees can be achieved with smaller sample sizes than for RandNLA sampling methods. Of course, this comes with computational trade-offs, which we discuss in Section 5.

and factor $d+1$ cannot in general be improved.

This exact error analysis is particularly useful in statistical modeling, where under additional assumptions about the vector $\mathbf{y}$, we wish to estimate accurately the generalization error of our estimator. Specifically, consider the following noisy linear model of the vector $\mathbf{y}$:

A number of extensions to Theorems 4 and 5 have been proposed, covering larger sample sizes as well as different statistical models .

3 Optimal approximation guarantees

As we have seen above, the non-i.i.d. nature of DPP sampling can lead to improved approximation guarantees, compared to standard RandNLA methods, when we wish to minimize the size of the down-sampled problem. We next discuss this in the low-rank approximation setting, i.e., when $n\sim d$. Here, cardinality constrained L-ensembles are known to achieve optimal $(1+\epsilon)$-approximation guarantees.

In Section 4.1, we used low-rank sketches to construct unbiased estimates for regularized least squares, given matrix $\mathbf{X}$ and a vector $\mathbf{y}$. However, even without introducing $\mathbf{y}$, a natural low-rank approximation objective for sketching $\mathbf{X}$ can be defined via a reduction to least squares . Namely, we can measure the error in reconstructing the $i^{th}$ row of $\mathbf{X}$ by finding the best fit among all linear combinations of the rows of the sketch $\mathbf{X}_{S}$. Repeating this over all rows of $\mathbf{X}$, we get:

and the size $r+r/\epsilon-1$ is worst-case optimal.

The task of low-rank approximation is often formulated in the context of symmetric positive semi-definite (p.s.d.) matrices. Let $\mathbf{L}$ be an $n\times n$ p.s.d. matrix. We briefly discuss how DPPs can be applied in this setting via the Nyström method, which constructs a rank $k$ approximation of $\mathbf{L}$ by using the eigendecomposition of a small $k\times k$ submatrix $\mathbf{L}_{S,S}$ for some index subset $S$.

We define the Nyström approximation of $\mathbf{L}$ based on a subset $S$ as the $n\times n$ matrix $\widetilde{\mathbf{L}}(S)=\mathbf{L}_{\cdot,S}\mathbf{L}_{S,S}^{\dagger}\mathbf{L}_{S,\cdot}$.

Originally developed in the context of obtaining numerical solutions to integral equations , this method has found applications in a number of areas such as machine learning , Gaussian Process regression and Indenpendent Component Analysis . The following result follows similarly as Theorem 6 and also provides the optimal sample size. We use $\|\cdot\|_{*}$ to denote the nuclear (trace) norm and $\mathbf{L}_{(r)}$ as the best rank $r$ approximation.

and the size $r+r/\epsilon-1$ is worst-case optimal.

4 Connections to RandNLA methods

The natural applicability of DPPs in the RandNLA tasks of least squares regression and low-rank approximation discussed above raises the question of how DPPs relate to traditional RandNLA sampling methods used for this task. As discussed in Section 2, one of the main RandNLA techniques for constructing rank-preserving sketches (i.e., relative to a tall matrix $\mathbf{X}$ with $n\gg d$) is i.i.d. leverage score sampling. (From this perspective, random projections can be seen as preprocessing or preconditioning the input so that leverage scores are approximately uniform, thereby enabling uniform sampling—in the randomly-transformed space—to be successfully used.) Even though leverage score sampling was developed independently of DPPs, this method can be viewed as an i.i.d. counterpart of the Projection DPP from Theorem 2.

This fact can be easily derived from the DPP properties discussed in Section 3. Recall that the marginal kernel of the Projection DPP is the projection matrix

The marginal probabilities of this distribution lie on the diagonal of the marginal kernel, which also contains the leverage scores of $\mathbf{X}$.

Thus, leverage score sampling can be obtained as a distribution constructed from the marginals of the Projection DPP. Naturally, when going from i.i.d. to non-i.i.d. sampling, we lose all the negative correlations between the points in a DPP sample, and therefore the expectation formulas and inequalities from the preceeding sections no longer hold for leverage score sampling. Furthermore, recall that to achieve a rank-preserving sketch (e.g., for least squares) with leverage score sampling for a full rank matrix $\mathbf{X}$ we require at least $d\log d$ rows (Theorem 1), whereas a Projection DPP generates only $d$ samples and also provides a rank-preserving sketch (Theorem 4). This shows that losing the negative correlations costs us a factor of $\log d$ in the sample size.

Another connection between leverage scores and Projection DPPs emerges in the reverse direction, i.e., going from i.i.d. to non-i.i.d. samples. Namely, a leverage score sample of size at least $2d\log d$ contains a Projection DPP with probability at least $1/2$ .

Let $j_{1},j_{2},...$ be a sequence of i.i.d. leverage score samples from matrix $\mathbf{X}$. There is a random set $T\subseteq\{1,2,...\}$ of size $d$ s.t. $\max\{i\in T\}\leq 2d\log d$ with probability at least $1/2$, and:

Many extensions of leverage scores have been proposed for use in the low-rank approximation setting (i.e., when $n\sim d$). Arguably the most popular one is called ridge leverage scores . Ridge leverage scores can be recovered as the marginals of an L-ensemble.

The typical sample size required for low-rank approximation with ridge leverage scores is at least $d_{\lambda}\log d_{\lambda}$, where $d_{\lambda}$ is the ridge effective dimension and also the expected size of the L-ensemble. Once again, the logarithmic factor appears as a trade-off coming from i.i.d. sampling.

A reverse connection analogous to Theorem 9, i.e., going from i.i.d. to non-i.i.d. sampling, can also be obtainined for ridge leverage scores , although only a weaker version, with $O(d_{\lambda}^{2})$ instead of $O(d_{\lambda}\log d_{\lambda})$, is currently known in this setting.

Let $j_{1},j_{2},...$ be a sequence of i.i.d. $\lambda$-ridge leverage score samples from matrix $\mathbf{X}$. There is a random set $T\subseteq\{1,2,...\}$ such that $\max\{i\in T\}\leq 2d_{\lambda}^{2}$ with probability at least $1/2$, and:

Sampling algorithms

One of the key considerations in RandNLA is computational efficiency of constructing random sketches. For example, the i.i.d. leverage score sampling sketch defined in Section 2 requires pre-computing all of the leverage scores. If done naïvely, this costs as much as performing the singular value decomposition (SVD) of the data. However, by employing fast RandNLA projection methods, one obtains efficient near-linear time complexity approximation algorithms for leverage score sampling .

In the case of DPPs, the challenge may seem even more daunting, since the naïve algorithm has exponential time complexity relative to the data size. However, the connections between leverage scores and DPPs (summarized in Table 1) have recently played a crucial role in the algorithmic improvements for DPP sampling. In particular, recent advances in DPP sampling have resulted in several algorithmic techniques which are faster than SVD, and in some regimes even approach the time complexity of fast leverage score sampling algorithms. See Table 2 for an overview.

We start by discussing fast sketching methods for approximating leverage scores more rapidly than by naïvely computing them via the SVD or a QR decomposition. This is a good illustration of RandNLA algorithmic techniques, and it is also relevant in our later discussion of DPP sampling.

For simplicity, we focus on constructing leverage scores for rank-preserving sketches (i.e., for a tall $n\times d$ matrix $\mathbf{X}$), but similar ideas apply to the low-rank approximation setup . Recall that the $i^{th}$ leverage score of $\mathbf{X}$ is given by $l_{i}=\mathbf{x}_{i}^{\scriptscriptstyle{\top}}(\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X})^{-1}\mathbf{x}_{i}$, which can be expressed as the squared norm of the $i^{th}$ row of the matrix $\mathbf{X}(\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X})^{-\frac{1}{2}}$. Assuming that $n\gg d$, the primary computational cost of obtaining this matrix involves two expensive matrix multiplications: first, computing $\mathbf{R}=\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X}$ (or a similarly expensive operation such as a QR decomposition or the SVD); and second, computing $\mathbf{X}\mathbf{R}^{-\frac{1}{2}}$. While each of these steps costs $O(nd^{2})$ arithmetic operations, showed that both of them can be approximated using efficient randomized sketching techniques.

The first step, which involves producing a matrix $\mathbf{\widetilde{R}}=\widetilde{\mathbf{X}}^{\scriptscriptstyle{\top}}\widetilde{\mathbf{X}}\approx\mathbf{R}$, requires using a rank-preserving sketch $\widetilde{\mathbf{X}}$, since we must ensure the invertibility of $\mathbf{\widetilde{R}}$. This can be achieved, e.g., using the Subsampled Randomized Hadamard Transform sketch , which is a random sketching matrix $\mathbf{S}$ that, with high probability, satisfies the subspace embedding property (Definition 1) on matrix $\mathbf{X}$. Crucially, the matrix multiplication $\widetilde{\mathbf{X}}=\mathbf{S}\mathbf{X}$ can be performed in time $O(nd\log d)$. The resulting sketch $\widetilde{\mathbf{X}}$ has $O(d\log n)$ rows, so computing $\mathbf{\widetilde{R}}=\widetilde{\mathbf{X}}^{\scriptscriptstyle{\top}}\widetilde{\mathbf{X}}$ takes $O(d^{3}\log n)$ time.

2 DPPs: Exact sampling

In this and subsequent sections, we discuss several algorithmic techniques for sampling from DPPs and $k$-DPPs. We focus on the general parameterization of a DPP via an $n\times n$ kernel matrix (either the correlation kernel $\mathbf{K}$ or the L-ensemble kernel $\mathbf{L}$), but we also discuss how these techniques can be applied to sampling from DPPs defined on a tall $n\times d$ matrix $\mathbf{X}$, which we used in Section 4.

We start with an important result of , which shows that any DPP can be decomposed into a mixture of Projection DPPs.

where $\mathbf{v}_{i}^{\scriptscriptstyle{\top}}$ denotes the $i^{th}$ row of the matrix $\mathbf{U}$. (Note that since $\mathbf{U}$ is an orthogonal matrix, its squared row-norms are also its leverage scores.) This gives us an easy way to sample the first point to be included in $S$. The key idea of how to continue the procedure is similar to computing the volume of a parallelogram spanning a pair of vectors $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$: first, pick one of them, say $\mathbf{v}_{1}$, and compute its length; then, multiply that by the length of $\mathbf{v}_{2}$ computed along the direction orthogonal to $\mathbf{v}_{1}$. An algorithm proposed by implements this probabilistically:

Sample one point $\mathbf{v}_{i}$ with $\Pr(i)\propto\|\mathbf{v}_{i}\|^{2}$;

Project all points $\mathbf{v}_{j}$ onto the subspace orthogonal to the sampled point $\mathbf{v}_{i}$;

Update the probabilities and go to step 1.

The overall sampling procedure suggested by is very efficient, if we are given the eigendecomposition of kernel $\mathbf{K}$ or of the L-ensemble kernel $\mathbf{L}$. It can also be easily adapted to sampling cardinality constrained DPPs. However, obtaining the eigendecomposition itself can be a significant bottleneck: it costs $O(n^{3})$ time for a general $n\times n$ kernel. If we are given a tall $n\times d$ matrix $\mathbf{X}$ such that $\mathbf{L}=\mathbf{X}\mathbf{X}^{\scriptscriptstyle{\top}}$, as was the case in Section 4, then the sampling cost can be reduced to $O(nd^{2})$ . There have been a number of attempts at avoiding the eigendecomposition in this procedure, leading to several approximate algorithms. Finally, approaches using other factorizations of the kernel matrix have been proposed , and these offer computational advantages in certain settings.

3 DPPs: Intermediate sampling

The DPP sampling algorithms from Section 5.2 can be accelerated with a recently introduced technique , which uses leverage score sampling to reduce the size of the $n\times n$ kernel matrix, without distorting the underlying DPP distribution. Recall from Section 4.4 that i.i.d. leverage score sampling can be viewed as an approximation of a DPP in which we ignore the negative correlations between the sampled points. Naturally, in most cases such a sample as a whole will be a very poor approximation of a DPP. However, with high probability, it contains a DPP of a smaller size (Theorems 9 and 11).

Sample $i_{1},...,i_{t}$ i.i.d. $\propto(\|\mathbf{v}_{1}\|^{2},...,\|\mathbf{v}_{n}\|^{2})$;

Construct $\widetilde{\mathbf{U}}$ from the $t$ rows $\frac{1}{\|\mathbf{v}_{i_{j}}\|}\mathbf{v}_{i_{j}}^{\scriptscriptstyle{\top}}$;

With prob. 1-\det\!\big{(}\frac{1}{t}\widetilde{\mathbf{U}}^{\scriptscriptstyle{\top}}\widetilde{\mathbf{U}}\big{)} go back to 1;

4 DPPs: Monte Carlo sampling

The advantages of this sampling procedure over the algorithm of are that we are not required to perform the eigendecomposition and that the computational cost of the MCMC chain scales linearly with $n$. The disadvantages are that the sampling is approximate and that we have to run the entire chain every time we wish to produce a new sample $S$.

In addition to the approach of for $k$-DPPs, proposed an MCMC designed specifically for Projection DPPs, and another Markov chain was developed for unconstrained L-ensembles by . Some of these MCMC approaches also apply beyond DPPs, to Strongly Rayleigh measures (see Figure 2).

Conclusions

We have briefly surveyed two established research areas which exhibit deep connections that have only recently began to emerge:

In particular, we discussed recent developments in applying DPPs to classical tasks in RandNLA, such as least squares regression and low-rank approximation; and we surveyed recent results on sampling algorithms for DPPs, comparing and contrasting several different approaches.

We expect that these connections will be fruitful more generally. As an example of this, we briefly mention a recently proposed mathematical framework for studying determinants, which played a key role in obtaining some of these results.

Determinant preserving random matrices. A square random matrix $\mathbf{A}$ is determinant preserving (d.p.) if all of its sub-determinants commute with taking the expectation, i.e., if:

for all index subsets $S$, $T$ of the same size. Not all random matrices satisfy this property (e.g., take $X\,\mathbf{I}_{2}$ for $X$ standard Gaussian), however there are many non-trivial examples. Moreover, this class of random matrices possesses a useful algebraic structure: if $\mathbf{A}$ and $\mathbf{B}$ are independent and d.p., then both $\mathbf{A}+\mathbf{B}$ and $\mathbf{A}\mathbf{B}$ are also determinant preserving. The first examples of d.p. matrices where given by (used in the analysis of a fast DPP sampling algorithm) and (used for eliminating bias in distributed optimization), and further discussion can be found in .

Of course, our survey of the applications of DPPs necessarily excluded many areas where this family of distributions appears. Here, we briefly discuss some other applications of DPPs which are relevant in the context of numerical linear algebra and RandNLA but did not fit in the scope of this work.

Implicit regularization. In many optimization tasks (e.g., in machine learning), the true minimizer of a desired objective is not unique or not computable exactly, so that the choice of the optimization procedure affects the output. Implicit regularization occurs when these algorithmic choices provide an effect similar to explicitly introducing a regularization penalty into the objective. This has been observed for approximate solutions returned by stochastic and combinatorial optimization algorithms , but a precise characterization of this phenomenon for RandNLA sampling methods has proven challenging. Recently, DPPs have been used to derive exact expressions for implicit regularization in RandNLA algorithms , connecting it to a phase transition called the double descent curve .

Experimental design. In statistics, the task of selecting a subset of data points for a down-stream regression task is referred to as experimental design . In this context, it is often assumed that the coefficients $y_{i}$ (or responses) are random variables obtained as a linear transformation of the vector $\mathbf{x}_{i}$ distorted by some mean zero noise. A number of optimality criteria have been considered for selecting the best subsets in experimental design. DPP subset selection has been shown to provide useful guarantees for some of the most popular criteria (such as for A-optimality and D-optimality), leading to new approximation algorithms .

Stochastic optimization. Randomized selection of mini-batches of data or subsets of parameters has been very successful in speeding up many iterative optimization algorithms. Here, non-uniform sampling can be used to reduce the variance in the iteration steps. In particular, showed that using a DPP for sampling mini-batches in stochastic gradient descent improves the convergence rate of the optimizer, whereas used a DPP sampler for the Stochastic Dual Newton Ascent method, showing an improved convergence analysis.

Monte Carlo integration. DPPs have been shown to achieve theoretically improved guarantees for numerical integration, i.e., using a weighted sum of function evaluations to approximate an integral. In particular, constructed a DPP for which the root mean squared errors of Monte Carlo integration decrease as $n^{-(1+1/d)/2}$, where $n$ is the number of function evaluations and $d$ is the dimension. This is faster than the typical $n^{-1/2}$ rate. See for other results on Monte Carlo integration with DPPs.

In conclusion, despite having been studied for at least forty five years, DPPs are enjoying an explosion of renewed interest, with novel applications emerging on a regular basis. Their rich connections to RandNLA, which we have only briefly summarized and which offer a nice example of how deep mathematics informs practical problems and vice versa, provide a particularly fertile ground for future work.

We would like to acknowledge DARPA, NSF (via the TRIPODS program), and ONR (via the BRC on RandNLA) for providing partial support for this work.

References