Patrick Flandrin

Arthur Gretton

Judith Rousseau

Pierre-Olivier Amblard

Peter Grün- Wald

Joaquim Ortega

Mylène, Subhro, Meixia, Raphaël (×2 Adam

Pascal

Benjamin

Chris Grégory

Julien Adrien Pierre

Ayoub, Diala, Yoann, Martin Hugo Barbara

Alexandre, Michaël, Xiaoyi Juan Manuel

my one past research engineer Guillaume again, who gets to be thanked twice in one sentence. Thanks for all your hard work and the stimulating discussions, on science, on problems of the world, or on the wacky random topics on which we like to debate in SigMA. Where else would I have learnt that you can drown in space, or that 18th century Moroccan monarch Mulay Ismail is allegedly second only to Genghis Khan in the number of children he fathered (868+). I would like to see the face of the future historian who will look into Arnaud's browsing history.

My gratitude extends to the whole SigMA team, among whom Pierre and Patrick stand out as mentors. Working among such friendly and knowledgeable people is a luxury. Beyond SigMA, the CRIStAL lab has been a great work environment, and many interesting discussions have been had with Guille, Charles, Pierre, Rayan, Laurentiu, Odalric, Émilie, Alessandro, Aurélien, and others. Our Esprit building only misses a music room and a rooftop with a barista, but hopefully a suitable commission will make sure we get these at some point. Thanks also to École Centrale, who hosted our team until we had our new building. They did not have a rooftop coffee bar either, but they actually had a music room, even if getting the keys was tricky and the admin department right above were complaining. A big thank you also goes to all the people running structures such as CRIStAL, École Centrale, Université de Lille, and DR18, in particular to all the staff who have been helping with my grants, Fatima, Christine, Myriam, Sonia, and Dominique; I know it is not always easy to administer researchers. I also have kind thoughts for my colleagues from the maths lab Painlevé; long live our joint point process workgroups! At the national level, I would like to thank INS2I for their support throughout my early career, and all the scientific communities to which I belong, in particular people in stats, signal, stochastic geometry, or physics, around GdRs like ISIS and GeoSto. Interactions with so many smart and friendly people is also what contributes to a researcher's happiness.

Many people above were already or have become friends, some even godfathers to my kids. This brings me to the source of my happiness, namely my family, with a special mention for the three people who share my life. First my two wonderful kids: Aliénor, who regularly offers to add some color whenever she sees a point configuration drawn on a draft at home, and Basile, who believes that my job consists in drinking coffee with friends. And, last but not least, my beloved Amélie, who has been continuously supporting me. I hope that she can forgive me for accidentally putting (another) one of her fragile sweaters in the washing machine yesterday evening.

Finally, to keep things relatively short and because I present in the manuscript my activity while in Lille, I have focused here on thanking people related to that period of time. Yet, I certainly do not forget all the wonderful people met in Oxford, and before that in Orsay, to whom I send my warm thoughts if they come across this page.

Chapter 1

Introduction " Million-to-one chances... crop up nine times out of ten. "

Terry Pratchett, Equal Rites, 1987 This document is the manuscript of my habilitation à diriger des recherches (HDR), described in the official texts1 as validating the high scientific level of the candidate, the originality of his scientific programme, his ability to master a research strategy in a sufficiently vast scientific or technological domain, and his ability to supervise young researchers. The typical form of the manuscript is a high-level summary of one's work after the PhD, emphasizing the research programme and directions, and I will follow that template.

My thesis at the departments of computer science (LRI) and particle physics (LAL) of Univ. Paris-Sud2 dates from December 2012. After that, I spent two years as a postdoc at the department of statistics of the University of Oxford, working on Monte Carlo methodology for big data (Bardenet, Doucet, andHolmes, 2014, 2017a), as well as Bayesian inference for cardiac electrophysiology (Johnstone et al., 2016a,b;[START_REF] Beattie | Sinusoidal voltage protocols for rapid characterization of ion channel kinetics[END_REF][START_REF] Ridder | A systematic strategy for estimating herg block potency and its implications in a new cardiac safety paradigm block potency and its implications in a new cardiac safety paradigm[END_REF]. In early 2015, I crossed back the Channel, and started as a CNRS junior permanent researcher (CR) at CRIStAL (short for research centre on computer science, signal processing, and automatic control of Lille). The major part of my work since 2015 has dealt with using repulsive point processes as computational tools, on which I now focus.

Many tasks in statistics, machine learning, or signal processing, can be cast as sampling, i.e., summarizing a large (possibly infinite) collection of items by a finite subset of these items. Variable selection in regression, for instance, aims at selecting a few representative features to reduce the dimensionality of a dataset. Monte Carlo integration aims at approximating integrals over a measure space by weighted sums of evaluations of the integrand at a handful of points in the support of the target measure. Finally, digital signal processing requires first to encode a square-integrable function by its orthogonal projections onto a finite number of reference functions. In all these examples, random sampling, i.e., drawing at random either features, quadrature nodes, or the reference functions on which to project, has a role to play. Loosely said, one motivation for random instead of deterministic sampling is to obtain simpler guarantees on the performance of the algorithm of interest under weaker assumptions. The output of standard Markov chain Monte Carlo algorithms, for instance, often satisfies a central limit theorem, with the classical Monte Carlo rate, as soon as the integrand is square-integrable. In contrast, deterministic quadrature guarantees in large dimension typically involve very smooth integrands and dimension-dependent rates; see e.g. quasi-Monte Carlo methods [START_REF] Dick | Digital nets and sequences: discrepancy theory and quasi-Monte Carlo integration[END_REF]. It is also common that, for performance guarantees of the same order, the best known randomized algorithms are computationally cheaper than their deterministic counterparts. For instance, we shall see in Chapter 4 that for feature selection with both Frobenius and spectral guarantees, we do not know of a better deterministic algorithm yet than enumerating all subsets of columns of a (fat) matrix, when randomized feature selection algorithms can simultaneously provide both types of guarantees in polynomial time. That is not to say that randomized algorithms are necessarily the final answer to every, and actually any, problem. Whether randomness actually helps is actually a deep and open question in theoretical computer science [START_REF] Impagliazzo | Hardness as randomness: a survey of universal derandomization[END_REF]. Moreover, it is sometimes possible to de-randomize a randomized algorithm, while transferring the bounds obtained with probabilistic tools to the obtained deterministic algorithm; see [START_REF] Deshpande | Efficient volume sampling for row/column subset selection[END_REF] for an example in column subset selection. To put it briefly, randomized sampling algorithms are a useful object of study, whether you accept them as a solution to your problem or only as an intermediate step.

If you agree with me that random sampling is an interesting research area,3 then my objective in this manuscript is to convince you that negatively dependent random sampling is also a promising one. In essence, drawing items from a large set independently leaves holes, and the corresponding sample looks far from anything regular. One illustrative example is shown in Figure 1.1. The left panel shows four draws from a homogenous Poisson process on the real line, intersected with [0, 1]. Equivalently, each row is independently obtained by first drawing a Poisson random variable N with some large enough mean, and conditionally on N , drawing N locations independently from the uniform distribution on [0, 1]. In each row, it is striking how independent locations fail to fill the unit interval evenly. If we evaluated a function f : [0, 1] → R at the samples shown in any of the rows of this left panel, the set of evaluations would poorly represent f , in that any oscillation of f in the observed holes would go unnoticed. We argue that by introducing dependence between the sampled items, and actually negative dependence or repulsiveness, we can obtain random samples that represent the initial set of items much more accurately. For instance, the right panel of Figure 1.1 shows samples of a translation-invariant determinantal point process. 4 We have been able to make more precise points on how negative dependence can solve sampling tasks at polynomial expense.

In Monte Carlo integration, we have shown that a specific class of negatively dependent yet computationally tractable distributions, called determinantal point processes, can provide central limit theorems with faster-than-square-root rates [START_REF] Bardenet | Monte Carlo with Determinantal Point Processes[END_REF]. We have even shown arbitrarily fast rates of decay for the mean square error of a related estimator under additional assumptions on the integrand [START_REF] Belhadji | Kernel quadrature with DPPs[END_REF](Belhadji et al., , 2020a)). Moreover, determinantal point processes, while providing tunable repulsiveness, can be sampled in polynomial time, which has made them a workhorse in my research programme. In feature selection for linear regression, we have shown that, again with a carefully designed determinantal point process, one could obtain a tight bound on the reconstruction error of a data matrix, as well as guarantees for linear regression based on the obtained column sample (Belhadji, Bardenet, and Chainais, 2020b). Finally, in signal processing, we have studied the repulsive point processes that appear as the zeros of specific random analytic functions, which result from applying time-frequency transforms to white noise [START_REF] Bardenet | On the zeros of the spectrogram of white noise[END_REF][START_REF] Bardenet | Monte Carlo with Determinantal Point Processes[END_REF]. One of my objectives in this manuscript, is to present these selected works under a common umbrella, namely random sampling with negative dependence.

In Chapter 2, I present some background material on repulsive point processes, i.e., random configurations of items in which similar or close-by items have little chance of simultaneously appearing. My flagship examples will be determinantal point processes and zeros of Gaussian analytic functions. My first own contributions appear in Chapter 3, where I discuss our results on Monte Carlo integration with determinantal point processes. In Chapter 4, I develop our results on column-subset selection and variable selection in linear regression. In Chapter 5, I explain our results on the zeros of time-frequency transforms of white noise, and discuss how they fit in my narrative of negatively dependent sampling. Each of Chapters 3 to 5 ends with chapter-specific open questions. I conclude in Chapter 6 by taking a step back and presenting higher-level open questions, as well as shortly discussing the part of my work since 2013 that I chose to omit in the manuscript for conciseness. Finally, Appendix A contains a history of my teaching, supervision, and grant management since 2013.5

Notation

This section is here for reference.

We shall constantly work with a complete metric space X, equipped with µ a Borel measure on the Borel sets of X. Our point processes will be random subsets of that space X. For our purpose, it is enough to think of either X ⊂ R d equipped with the Lebesgue measure, or X = [n] := {1, . . . , n} a finite set, equipped with its counting measure.

We shall also often consider submatrices. We write matrices in bold font (K, L, etc.). For a matrix K, a subset A of |A| = a indices of rows of K, and a subset B of |B| = b indices of columns of K, we write K A,B for the a × b submatrix of K formed by the rows indexed by A and the columns indexed by B. We also write K A instead of K A,A , and K :,A when we keep all rows but only the columns indexed by A.

Chapter-specific notation will be introduced on-the-fly. Note that because Chapters 3, 4, and 5 correspond to work published in different communities, some conventional notations are bound to conflict. For instance, k is the kernel of an RKHS in Chapter 3, while in Chapter 4, k ∈ N is the dimension of the reference principal subspace. These notations are so classical in the respective literatures that I refrained from changing, hoping that the context will keep it clear what k is in which chapter.

Chapter 2

Repulsive point processes " 'Y'know,' he said, 'it's very hard to talk quantum using a language originally designed to tell other monkeys where the ripe fruit is.' " Let us fix a complete metric space X, with µ a Borel measure on the Borel sets of X. For this manuscript, it is enough to bear in mind the two cases of X ⊂ R d equipped with the Lebesgue measure, and a finite set X equipped with its counting measure. Define the space of configurations as Conf(X) := S ⊂ X : |S ∩ A| < ∞ for all compact A ⊂ X , where |B| stands for the cardinality of set B. In other words, configurations are locally finite subsets of X.

Definition 1 A point process S on X is a random configuration 1 of points in X.

1 Stricto sensu, we are defining here simple point processes, i.e., such that with probability 1, S does not contain a given point of X more than once. Since all point processes in this work are simple, we take this shortcut and identify a (simple) point process with a random configuration.

Correlation functions

Correlation functions are a common way of defining and studying point processes. For any k 1, the k-point correlation function

ρ k : X k → [0, ∞] satisfies, when it exists, E     x 1 ,...,x k ∈S x i =x j if i =j f (x 1 , . . . , x k )     = X k f (x 1 , . . . , x k ) ρ k (x 1 , . . . , x k ) k i=1 dµ(x i ) (2.1)
for any bounded (or positive) and measurable test function f . Here

X k = k X × • • • × X
and the symbol E in (2.1) stands for the expectation under the law P of the random variable S. Thus, the k-point correlation function ρ k encodes the distribution of k-tuples of points from S. Indeed, an informal rewriting of (2.1) reads

ρ k (x 1 , . . . , x k )dµ(x 1 ) . . . dµ(x k ) = P There are at least k points in S, one around each of x 1 , . . . , x k . . (2.2)
It is common to actually define a point process S by a sequence (ρ k ) of compatible correlation functions; see [START_REF] Daley | An introduction to the theory of point processes[END_REF]. By compatible, we mean that not every sequence (ρ k ) actually defines a point process, and that usually a mathematical argument for existence is necessary.

Of particular practical significance in the description of point processes are the one-and two-point correlation functions. The one-point correlation function ρ 1 describes marginal information about the particles, and is called the intensity of the point process. In particular, (2.2) applied to f = 1 A implies that the average number of points falling in A is A ρ 1 dµ. The second correlation function ρ 2 describes pairwise interactions, and is often discussed in its normalized form

g(x, y) = ρ 2 (x, y) ρ 1 (x)ρ 1 (y) . (2.3)
Finally, when X = R d and the distribution of S is translation-and rotation-invariant, g only depends on r = xy , so that we write g(x, y) = g 0 (r). The function g 0 is called the pair correlation function of S.

Example 1 Let λ : X → R + be locally integrable, that is, B λ(x)dµ(x) < ∞ for every bounded B ⊂ X. Further assume for simplicity that the measure λdµ has no atom. 2 The point process with correlation functions

ρ k (x 1 , . . . , x k ) = λ(x 1 ) . . . λ(x k ), k 1, (2.4)
always exists and is called the Poisson point process with parameter function λ.

For a Poisson process with parameter function λ, the one-point correlation function is by definition ρ 1 = λ. The pair correlation function (2.3) is constant, g(x, y) = 1, which indicates that the relative positions of x and y do not alter the marginal probability of seeing them co-occur. In a sense, the factorized form of (2.4) indicates the lowest level of correlation among points. When g(x, y) > 1, Equation (2.2) indicates that pairs are more likely to occur around (x, y) than under a Poisson process with the same intensity. Similarly g(x, y) < 1 indicates that pairs are less likely to occur. Consequently, one way to define a repulsive point process3 is to require that g(x, y) < 1 for all x, y.

Figure 2.1 shows samples of three translation-invariant point processes in R, along with their pair correlation functions. The three point processes share the same (constant) one-point correlation function, and the left panel corresponds to a Poisson point process. Note how the pair correlation function in Figure 2.1(f) keep below 1, and even close to zero around zero, indicating fewer small pairwise distances than the reference Poisson point process in Figure 2.1(d). This is a sign of a very regular, more grid-like distribution of the points, as seen in Figure 2.1(c). On the contrary, the pair correlation function of Figure 2.1(e) shows more small pairwise distances than the reference Poisson process: points are lumped together, as confirmed by Figure 2.1(b). The latter process is actually a permanental point process, and was introduced by Macchi ( 2017) is her 1972 thesis. 4It is a Poisson process with a random (smooth) parameter function, which is why we see lumps. We refer the reader interested in permanental point processes to our recent survey on point processes in quantum optics [START_REF] Bardenet | From point processes to quantum optics and back[END_REF].

We now define a flexible class of repulsive point process, of which the right-hand panel of Figure 2.1 shows an instance: determinantal point processes.

Determinantal point processes (DPPs)

DPPs were also introduced by [START_REF] Macchi | Point processes and coincidences -Contributions to the theory, with applications to statistical optics and optical communication, augmented with a scholion by[END_REF] in her 1972 thesis, and first featured in [START_REF] Bénard | Detection and "emission" processes of quantum particles in a "chaotic" state[END_REF][START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF], this time to model detection times in quantum beams of fermions.

Definition 2 A point process is determinantal if there exists K : X × X → R or C such that the n-point correlation function (2.2) exists for every n and reads

ρ n (x 1 , . . . , x n ) = det K(x k , x ) n k, =1 , x 1 , . . . , x n ∈ X. (2.5)
We say that K is the kernel of the DPP, µ its reference measure, and we denote the point process by DPP(K, µ).

Just like a covariance kernel is enough to describe a zero-mean Gaussian process, the kernel K of a DPP thus encodes how the points in the random configurations interact. However, the existence of a point process with (2.5) as its correlation functions is, in general, a difficult question. It is easy to see that the kernel has to be positive definite, so that the right-hand side of (2.5) is always non-negative. Henceforth, all DPP kernels are thus assumed to be positive definite. But non-negativity is not sufficient for (2.5) to consistently define a point process.

Projection DPPs

A canonical way to ensure existence of a DPP is to restrict oneself to so-called projection DPPs, which generate configurations of a constant number N < ∞ of points Palmost surely, i.e. S = {x 1 , . . . , x N }. More precisely, consider N orthonormal functions ϕ 0 , . . . , ϕ N -1 in L 2 (µ), that is, ϕ k (x)ϕ (x)µ(dx) = δ k , and take for kernel

K(x, y) = N -1 k=0 ϕ k (x)ϕ k (y).
(2.6)

In this setting, it turns out that the (permutation invariant) random variables x 1 , . . . , x N with joint probability distribution

1 N ! det K(x i , x ) N i, =1 × N i=1 dµ(x i ) (2.7)
generate a DPP with kernel K(x, y). (Hough et al., 2006, Section 2) gives an elementary proof that (2.7) satisfies (2.5). The resulting DPP is called a projection DPP, in reference to the fact that its kernel defines a projection operator.

Example 2 (CUE) Take X = [0, 1], µ to be the uniform distribution, and the orthonormal functions ϕ to be the Fourier basis ϕ (x) = e 2iπ x , 0. The projection DPP with kernel

x, y → N -1 =0 e 2iπ (x-y) = e iπ(N -1)(x-y) sin πN (x -y) sin π(x -y)
given by (2.6) is called the circular unitary ensemble (CUE).

Note that from the definition (2.5) of the correlation functions of a DPP, for any valid kernel and any non-vanishing function f , the kernel x, y → f (x)K(x, y)/f (y) defines the same DPP as K. The kernel of the CUE in Example 2 can thus be simplified to

x, y → sin πN (xy) sin π(xy) .

(2.8)

A CUE sample is shown in the right panel of Figure 2.2. Note how regularly spread the points of each realization of the CUE are on the circle. In comparison, i.i.d. draws from µ = dx in the left panel seem less well organized, creating clusters and holes. Historically, the CUE was introduced as the law of the eigenvalues of a Haar-distributed random unitary matrix, as we shall recall in Section 3.1. Since the Haar distribution on unitary matrices is actually easy to sample, this connection with random matrices can actually be seen as a O(N 3 ) exact sampling algorithm. Thanks to a scaling argument, this also provides an approximate sampling algorithm for another well-known DPP, the sine process. Indeed, if (θ i ) N i=1 ∼ CUE, then the average spacing between the points of (θ i ) is of order 1/N . Scaling the point process to increase this mean distance to 1, the random measure δ N θ i converges in distribution (a.k.a. weakly) to the DPP on R with kernel

-1 0 1 -1 0 1 (a) Uniform -1 0 1 -1 0 1 (b) Haar
K sine (x, y) = sin(x -y) x -y , (2.9)
a DPP called the sine process. Approximate samples obtained by scaling CUE samples are shown in Figure 2.1(c). The sine process is known for its universality: it appears as the scaling limit not only of the CUE, but of many DPPs. Another reason for its fame is its link with the zeros of the Riemann zeta function; see e.g. [START_REF] Rudnick | Zeros of principal L-functions and random matrix theory[END_REF] and references therein. We refer to (Lambert, 2016, Section 3.4) for a survey of the properties of the sine process and precise statements on scaling limits of DPPs. We now switch to another reference DPP from random matrix theory.

Example 3 (Ginibre ensemble) Take X = C, the reference measure µ to be the Gaussian measure π -1 e -|z| 2 , and the orthonormal functions to be ϕ (z) := (k!) -1/2 z k . The projection DPP with kernel (2.6) is called the Ginibre ensemble of size N .

Like the CUE, the Ginibre ensemble actually rose to attention in random matrix theory, this time as the distribution of the eigenvalues of an N ×N random matrix filled with i.i.d. unit complex Gaussians; see e.g. [START_REF] Anderson | An introduction to random matrices[END_REF]. One can actually let N → ∞ and define a limiting DPP, called the Ginibre ensemble, see (Hough et al., 2006, Example 15).

A geometric chain rule for projection DPPs

Besides their constant cardinality, one key property for computational applications of projection DPPs is that (2.7) can easily be rewritten as a chain rule, and this chain rule brings insight on the link between DPPs and geometry in Hilbert spaces.

We first illustrate our point by doing explicit computations when the rank of the kernel (2.6) is N = 2, and then discuss the chain rule for arbitrary N . Let thus ϕ 0 , ϕ 1 ∈ L 2 (µ) be two orthonormal functions. Consider the projection kernel (2.6), which simply reads

K(x, y) = ϕ 0 (x)ϕ 0 (y) + ϕ 1 (x)ϕ 1 (y).
Samples from the corresponding DPP will thus have N = 2 points almost surely, drawn from (2.14), which becomes

1 2 det K(x, x) K(x, y) K(y, x) K(y, y) dµ(x)dµ(y) = 1 2 K(x, x)dµ(x) × K(y, y) - K(x, y) 2 K(x, x) dµ(y)
(2.10) Now, it is routine to check that both terms in the product in the RHS of (2.10) integrate to 1. This implies that (2.10) is a chain rule. In other words, to obtain {x, y} ∼ DPP(K, µ), it is enough to first draw x ∼ 1 2 K(x, x)dµ(x), and then

y|x ∼ K(y, y) - K(x, y) 2 K(x, x) dµ(y).
As long as the two sampling steps can be implemented, e.g. using rejection sampling, this is already an interesting remark. But an important remark is that this chain rule is actually a simple area computation using the base times height formula. Indeed, consider the parallelotope formed in Span(ϕ 0 , ϕ 1 ) by the so-called feature functions K(x, •) and K(y, •). The squared length of K(x, •) (the base) is

K(x, •)K(x, •)dµ(x) = |ϕ 0 (x)| 2 |ϕ 0 | 2 dµ + |ϕ 1 (x)| 2 |ϕ 1 | 2 dµ (2.11) = K(x, x),
(2.12) while the squared length of the projection of K(y, •) onto the orthogonal of Span(K(x, •))

(the height) is K(y, y) - K(y, •) K(x, •) K(x, x) dµ(x) 2 = K(y, y) - K(x, y) 2 K(x, x) .
(2.13) From (2.12) and (2.13), it should now be clear that the chain rule (2.10) computes the squared volume of a parallelotope. For a projection kernel K of rank N 1, similar computations lead to an Ndimensional volume computation through a generalized "base times height" formula, namely

(2.7) = N i=1 1 N -i + 1 P H i-1 K(x i , •) 2 L 2 (µ)
dµ(x i ).

(2.14)

In (2.14), P H is the orthogonal projection onto a subspace H of L 2 (µ),

H 0 = Span(ϕ 0 , . . . , ϕ N -1 ),
and for every i > 1, H i-1 is the orthocomplement in H 0 of Span (K(x , •), 1 i -1) .

Like (2.10), (2.14) is a chain rule since each term in the product is a probability measure (Hough et al., 2006, Proposition 19).

Beyond sampling algorithms, we believe that the geometric interpretation of the chain rule (2.14) is paramount to understanding DPPs: the favoured point configurations in a projection DPP are those that correspond to feature vectors that span a large volume. Several contributions described in this manuscript come from understanding what feature vectors to choose in a particular problem involving linear structures, so that large volumes correspond to a small value of the problem's objective.

DPPs with Hermitian kernels

The projection kernels of Section 2.2.1 are a special case of Hermitian kernels, that is K(x, y) = K(y, x). For Hermitian kernels, the Macchi-Soshnikov theorem gives a necessary and sufficient condition for the existence of the associated DPP. The Hermiticity of the kernel guarantees that the integration operator Σ : L 2 (dµ) → L 2 (dµ) defined by

Σ : f → K(•, y)f (y)dµ(y) is self-adjoint.
Theorem 1 (Macchi-Soshnikov;[START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF][START_REF] Soshnikov | Determinantal random point fields[END_REF] If Σ is selfadjoint and trace-class, then existence of the associated DPP is equivalent to all the eigenvalues of Σ lying is [0, 1].

Being trace-class is a technical, but unrestrictive assumption, which guarantees in particular that Σ is compact, so that we can talk about its eigenvalues. Actually, a slight modification of Theorem 1 allows treating locally trace-class operators, i.e., such that P C ΣP C is trace-class for any compact C ⊂ X, with P C the multiplication operator by the indicator function 1 C . Checking that a given operator is locally trace-class can be hard work, but there are useful criteria in the literature; see e.g. [START_REF] Gohberg | Classes of linear operators[END_REF]. For instance, since our kernels in this section are positive definite and Hermitian, a sufficient condition for being locally trace-class is that K is continuous and C K(x, x)dµ(x) < ∞ for any compact C ⊂ X. In particular, the sine kernel defined in (2.9) is locally trace-class. This is also the criterion used in spatial statistics by [START_REF] Lavancier | Determinantal point process models and statistical inference[END_REF], who show how covariance functions for Gaussian processes give rise to DPPs in X = R d using the Macchi-Soshnikov theorem. The constraint in Theorem 1 on the spectrum of the integration operator then further constrains the parameters of the covariance function. A simple such example is the Gaussian DPP.

Example 4 (Gaussian kernel; [START_REF] Lavancier | Determinantal point process models and statistical inference[END_REF]) Take µ to be the Lebesgue measure on X = R d and K(x, y) = α exp(-x/ 2 ).

Provided the parameters are restricted to α ( √ π ) -d , DPP(K, µ) exists (Lavancier et al., 2014, Theorem 2.5).

It is interesting to note that the locally trace-class condition guarantees that the average number of points that fall in a bounded A ⊂ X is finite, since the definition of the one-point correlation function of a DPP leads to

E x∈S 1 A (x) = A K(x, x)dµ(x) < ∞.
(2.15)

In particular, for the Gaussian DPP of Example 4, the expected number of points in a given compact C is α times the Lebesgue measure of C. The parameter α thus controls the intensity of the point process. Meanwhile, the pair correlation function reads

g(x, y) = 1 -exp -2 x -y 2 / 2 ,
and thus only depends on , which controls repulsiveness. The fact that parameters of K can be easily interpreted is part of the appeal of DPPs in statistical modeling.

For the rest of this section, we turn to examining Hermitian DPPs in the special case when the cardinality |X| = n is finite. The kernel is now automatically trace-class, and the Macchi-Soshnikov theorem thus only requires Spec(K) ∈ [0, 1], where the kernel matrix is defined as K = (( K(i, j) )) 1 i,j n .

Example 5 (L-ensemble) Take X = {1, . . . , n} and µ to be the uniform distribution over X. Let L ∈ R n×n be symmetric and positive definite, then

K = (I + L) -1 L (2.16)
is the kernel matrix of a DPP.

The construction in Example 5 satisfies the Macchi-Soshnikov conditions by design, and it is convenient to define finite DPPs through any positive definite matrix L. Indeed, by construction, (2.16) ensures that K is symmetric, with eigenvalues in [0, 1]. Furthermore, in that case, L plays the role of a likelihood kernel, in the sense that if S is drawn from the corresponding L-ensemble,

P(S = A) = det L A det I + L , ∀A ⊂ {1, . . . , n}, (2.17)
where L A is the |A| × |A| submatrix of L indexed by the rows and columns in A; see e.g. [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF]. This is to be compared to K being the kernel that encodes the correlation functions, since in the finite setting of Example 5, (2.5) is equivalent to

P(A ⊂ S) = det K A , ∀A ⊂ {1, . . . , n}.
L-ensembles are popular DPPs in machine learning, where a typical construction is as follows (Kulesza and Taskar, 2012, Section 4.2.1). Consider a large ground set of n objects, indexed by X = {1, . . . , n}. Say, n is the number of sentences in a large text, and you want to extract a few (random) sentences from the text that are representative of the whole text.

Example 6 (Summary extraction; [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF] In the finite setting of Example 5, let L ij = q i Φ T i Φ j q j , 1 i, j n, where

• q i ∈ R measures the marginal relevance of sentence i. For instance, it can be chosen as an increasing function of the number of nouns and a decreasing function of the number of pronouns in the sentence. A sentence with many pronouns is unlikely to be useful in a summary.

• Φ i ∈ R p is a normed vector of features of sentence i, such that two sentences can be judged dissimilar if and only if Φ T i Φ j is close to zero. For instance, take Φ i to have one entry per word w in the whole text, and Φ i (w) to be an increasing function of how frequent word w is in sentence i and how rare it is in the text. Pairs of features with a large inner product will likely indicate sentences with similar content.

A draw from the corresponding DPP will be a set of indices of sentences that are marginally relevant, but jointly dissimilar; in other words, a potentially good summary.

To help with choosing the model, it is natural to actually parametrize (q i ) and (Φ i ) and learn them from a set of human-made summaries; see (Kulesza and Taskar, 2012, Section 4.2.1). Learning can be achieved by likelihood-based methods, since L gives the likelihood (2.17), although one has to care for the cost of evaluating the normalization constant in (2.17). However, in my opinion, a more severe downside of modeling L directly is that the parameters of L are usually hard to interpret. In particular, the number of items in a realization of S is hard to control. In that respect, working with K is more natural, as we noticed when interpreting the parameters of the Gaussian DPP in Example 4. One way to recover interpretability and a control over the number of points in a realization is to work with projection DPPs as in [START_REF] Gautier | Zonotope hit-and-run for efficient sampling from projection DPPs[END_REF], taking the number p of features in Example 6 to be the number of desired items in a summary. If the features are linearly independent, the resulting projection kernel has rank p. This construction using projection DPPs has however not been favoured in machine learning so far. MLers typically prefer using k-DPPs [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF].

For a given k ∈ N, a k-DPP is a DPP conditioned to have cardinality k. While in general a k-DPP is not a DPP, a projection DPP is a k-DPP, since it has constant cardinality almost surely. Usually, k-DPPs are built by conditioning an L-ensemble, i.e., replacing (2.17) by

P(S = A) ∝ det L A 1 |S|=k .
(2.18)

In summary extraction applications such as (Kulesza and Taskar, 2012, Section 4.2.1), which inspired Example 6, the L-ensemble is usually replaced by a k-DPP with the same L matrix. Interestingly, k-DPPs are statistical mixtures of projection DPPs; see e.g. [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF]. Combined to the chain rule for projection DPPs in (2.14), this observation yields gives a generic sampling algorithm for k-DPPs. We end this section with a fundamental example of k-DPP.

Example 7 (Volume sampling; [START_REF] Deshpande | Efficient volume sampling for row/column subset selection[END_REF] Let X be an n × d data matrix, with d n, take L = X T X, and let k n. The k-DPP (2.18) on {1, . . . , d} is called volume sampling.

Volume sampling was introduced by Deshpande and Vempala (2006) as a distribution for randomized column subset selection. Indeed L ij is the inner produce of columns i and j of X, and intuitively, sampling S proportionally to size-k subdeterminants of L should choose a diverse set of columns; see Chapter 4 for more details.

Sampling DPPs

To conclude on Hermitian kernels, still assuming that K is Hermitian and that Σ is traceclass, one can write DPP(K, µ) as a statistical mixture of projection DPPs (Hough et al., 2006, Theorem 7), with the mixture weights being simple functions of the eigenvalues of the kernel. Together with the chain rule for projection DPPs (2.14), this provides a generic abstract algorithm to sample DPPs with Hermitian kernels. It is generic because it applies widely. Is it abstract because sampling from the conditionals in the chain rule is not straightforward in general, and rejection sampling routines have to be carefully designed for each kernel; see e.g. (Gautier, Bardenet, and Valko, 2019b).

DPPs are (usually) repulsive

As statistical models, DPPs are often thought of as repulsive point processes. By definition (2.7), the pair correlation function (2.3) of a DPP reads g(x, y) = 1 -K(x, y)K(y, x) K(x, x)K(y, y) .

If K is Hermitian and does not vanish, this indeed implies that g < 1, so that the DPP is repulsive in the sense of Section 2.1. Note however that many DPP kernels in random matrix theory are not Hermitian (Anderson et al., 2010, Section 4.2.9), and thus should not be thought of as repulsive. Moreover, non-Hermitian DPPs are starting to be used as statistical models that capture attractiveness as well as repulsive behaviour, e.g. in recommendation tasks such as basket completion in online retail [START_REF] Gartrell | Learning Nonsymmetric Determinantal Point Processes[END_REF].

To explain the idea, if you just put a coffee machine in your basket, then you should be recommended items such as coffee grains (attractiveness), but the seller should avoid offering a second coffee machine (repulsiveness). On a historical note, DPPs were first formalized in the 1972 thesis of [START_REF] Macchi | Point processes and coincidences -Contributions to the theory, with applications to statistical optics and optical communication, augmented with a scholion by[END_REF]. Her initial objective was a point process description of a beam of fermions [START_REF] Bénard | Detection and "emission" processes of quantum particles in a "chaotic" state[END_REF], using correlation functions (2.2) and the by-then recent quantum coherence formalism of [START_REF] Glauber | Quantum theory of coherence[END_REF]. Importantly, the repulsive behaviour of the resulting (Hermitian) DPPs corresponds to a physical phenomenon for fermions known as anti-bunching.

We refer to our survey [START_REF] Bardenet | From point processes to quantum optics and back[END_REF] for more on the physical roots of DPPs. The connections between DPPs and quantum (electronic) optics have been more or less dormant since the 70s, and DPPs had to wait for the 2000s to be developed and popularized by random matrix theory [START_REF] Anderson | An introduction to random matrices[END_REF]. After that, they have been in particular investigated as repulsive statistical models in spatial statistics [START_REF] Lavancier | Determinantal point process models and statistical inference: Extended version[END_REF]Bardenet et al., 2017b) and machine learning [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF]. For more information on determinantal point processes, we refer the reader to [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF][START_REF] Hough | Determinantal Processes and Independence[END_REF][START_REF] Johansson | Random matrices and determinantal processes[END_REF][START_REF] Soshnikov | Determinantal random point fields[END_REF][START_REF] Lyons | Determinantal probability measures[END_REF][START_REF] Kulesza | Determinantal point processes for machine learning[END_REF][START_REF] Lavancier | Determinantal point process models and statistical inference: Extended version[END_REF], as well as to the theses of our students Gautier (2020), [START_REF] Belhadji | Subspace sampling using determinantal point processes[END_REF], and the HDR manuscript of [START_REF] Hardy | Attractive repulsion[END_REF]. In the rest of this section, we introduce some examples of DPPs that will play a role in this manuscript.

Zeros of Gaussian analytic functions

Besides DPPs, a natural way to define a repulsive point process on X = C ≈ R 2 equipped with the Lebesgue measure is by looking at the zeros of a random analytic function. Intuitively, an analytic function is very smooth, so that its zeros should be far from each other. Restricting oneself to analytic-valued Gaussian processes further allows to describe the correlation functions (2.2) of the zeros using the covariance kernel of the Gaussian process.

In the following, N C (0, σ 2 ) refers to the law of a centered Gaussian variable on C with variance σ 2 , namely the law of Z := σ(X + iY )/ √ 2 with X, Y identically distributed and independent (i.i.d.) real standard N (0, 1) variables; note that the √ 2 is here to have

Var(Z) := E|Z| 2 -|E(Z)| 2 = σ 2 .
More generally, given m 1, a complex Gaussian vector of mean µ ∈ C m and covariance matrix Σ := AA * , with A ∈ M n (C), is a random variable with the same distribution as AX + µ, where Z := (Z 1 , . . . , Z d ) with (Z i ) m i=1 i.i.d. N C (0, 1) random variables. Given an open subset Λ ⊂ C, we let A(Λ) be the space of analytic functions on Λ.

A constructive definition of GAFs

A Gaussian analytic function (GAF) on an open subset Λ ⊂ C is a random variable GAF taking its values in the space A(Λ) of analytic functions such that, for any m 1 and z 1 , . . . , z m ∈ Λ, the vector (GAF(z 1 ), . . . , GAF(z m )) is a complex Gaussian vector. In this work, we actually slightly restrict the definition.

Definition 3 All GAFs in this work are defined by

GAF(z) = ∞ k=0 ξ k Ψ k (z), (2.19) where (ξ k ) k∈N is a sequence of i.i.d. N C (0, 1) random variables and (Ψ k ) k∈N is a sequence in A(Λ) satisfying sup z∈K ∞ k=0 |Ψ k (z)| 2 < ∞ (2.20) for any compact K ⊂ Λ.
Condition (2.20) ensures that GAF converges almost surely on every compact subset of Λ, and thus that GAF is a well-defined analytic function, see (Hough et al., 2009, Lemma 2.2.3). As a side note, the proof of the latter lemma, while elementary, is a good example of the beautiful interplay of analytic and probabilistic notions that is behind GAFs.

As a centered Gaussian process on the space of analytic functions A(Λ), GAF is completely characterized by its covariance structure, namely

C(z, w) := E GAF(z)GAF(w) = ∞ k=0 Ψ k (z)Ψ k (w), z, w ∈ Λ.
(2.21)

Our particular interest lies in the random set GAF -1 (0) of zeros of GAF. More precisely, our interest lies in the random zeros of GAF, that is, zeros that are not present with probability one. Indeed, z → (z -1)GAF(z) is still a GAF, but the (deterministic) zero in 1 is not of interest to us. (Hough et al., 2009, Lemma 2.4.1) states that the random zeros of any nonzero GAF are pairwise distinct with probability one. Thus, the random zeros form a (simple) point process in the sense of Definition 1. Like the DPPs of Section 2.2, the correlation functions (2.2) of the zeros of a GAF are known in closed form (Hough et al., 2009, Corollary 3.4.2). GAFs are very constrained models. For instance, it is known that the law of the point process GAF -1 (0) is characterized by its first intensity ρ 1 solely and moreover that, if two GAFs, say GAF and GAF, share the same zero set distribution, namely GAF -1 (0)

law = GAF -1
(0), then there exists a non-vanishing, deterministic ϕ ∈ A(Λ) such that GAF = ϕ × GAF, see (Hough et al., 2009, Theorem 2.5.2).

Three canonical GAFs

There are three prototypical spaces in complex geometry: the complex plane C, the hyperbolic plane H, and the sphere S. For instance, the so-called uniformization theorem says that any simply connected Riemann surface is conformally equivalent to C, H or S. The complex plane is a Riemannian manifold with null curvature when equipped with the Euclidean metric g flat . The hyperbolic plane H can be modeled as the unit disc

D := {z ∈ C : |z| < 1} with the hyperbolic metric (1 -|z| 2 ) -2 g flat ,
giving a Riemannian manifold of constant negative curvature. It can alternatively be modeled as the upper-half plane C + := {z ∈ C : Im(z) > 0} with the metric Im(z) -1 g flat . The sphere S can be modeled as the extended complex plane C ∪ {∞} with the spherical metric (1 + |z| 2 ) -2 g flat , yielding a Riemannian manifold of constant positive curvature. This metric is the image by the stereographic projection of the natural metric on the unit sphere of R 3 induced by the Euclidean structure of R 3 . The group of isometries for each of these three spaces is well-known, and it is natural to ask which GAFs have their zeros invariant under these isometries. It turns out that, up to multiplication by a non-vanishing and deterministic analytic function, there is only a one-parameter family of GAFs satisfying this condition for each space; see [START_REF] Hough | Zeros of Gaussian analytic functions and determinantal point processes[END_REF]. We shall meet these so-called planar, hyperbolic, and spherical GAFs in Chapter 5. Meanwhile, a realization of the zero sets of each of the three invariant GAFs is shown in Figure 2.3 as white dots. For readers unfamiliar with complex geometry, Figure 2.3 can also serve as an illustration of the previous paragraph. Roughly speaking, on all three subfigures, the white dots are uniformly spread, with roughly equal distance between each pair of nearest neighbors. Figure 5.2 represents a portion of the complex plane, and the flat metric corresponds to the intuitive idea of a uniform distribution of zeros on this portion. Figure 5.3 represents a portion of C + with the hyperbolic metric Im(z) -1 g flat . Again, the white dots are uniformly spread for this metric, with distances between nearest neighbors roughly equal. Figure 5.4 similarly represents the sphere with the spherical metric.

Chapter 3

Monte Carlo integration with determinantal point processes " God does not play dice with the universe; He plays an ineffable game of His own devising, which might be compared, from the perspective of any of the other players [i.e. everybody], to being involved in an obscure and complex variant of poker in a pitch-dark room, with blank cards, for infinite stakes, with a Dealer who won't tell you the rules, and who smiles all the time. My initial motivation to look into DPPs came from discussions around a swimming pool with my long-time friend Adrien Hardy in the Summer holidays of 2013, while we were both postdocs, respectively in Oxford and Stockholm. Adrien was trying to give me an overview of the big problems of random matrix theory, while I was trying to do the same with Monte Carlo methodology. What caught our attention is that fundamental theorems about the asymptotic behaviour of the spectrum of certain random matrices could be interpreted as yielding fast Monte Carlo estimators in R and C. By fast, I mean estimators with mean square error decreasing faster than the inverse of the number of evaluations of the integrand, that is, faster than the usual Monte Carlo error. One natural question was then to understand what variance cancellation mechanisms made this fast decay possible, and try to construct quadrature nodes in arbitrary dimension that featured the same limiting behaviour. This turned out to be possible, though quite technically involved, using specific DPPs. It gave [START_REF] Bardenet | Monte Carlo with Determinantal Point Processes[END_REF], which is the starting point of my work on repulsive point processes, although the publication date does not reflect it. Section 3.1 describes a classical theorem by [START_REF] Johansson | On random matrices from the compact classical groups[END_REF] that provided motivation for our work, and our main results are given in Section 3.2.

After I presented the results of [START_REF] Bardenet | Monte Carlo with Determinantal Point Processes[END_REF] at a seminar in Bristol, UK, Mathieu Gerber, who was in the audience, thought that he had seen something similar somewhere. In the afternoon following the seminar, we engaged in serious scientific archaeology to find out the source of the remark he remembered. Following a chain of progressively more obscure references, we ended up with an English translation of a paper by [START_REF] Ermakov | Polynomial Approximations and the Monte-Carlo Method[END_REF], two Soviet mathematicians that actually had proposed using a projection DPP for Monte Carlo, almost 60 years before us, and some 15 years before DPPs where even formally defined by [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF]. Of course, [START_REF] Ermakov | Polynomial Approximations and the Monte-Carlo Method[END_REF] did not have the mathematical toolbox back then that we had at our disposal in [START_REF] Bardenet | Monte Carlo with Determinantal Point Processes[END_REF], and did not prove a CLT, for instance, but they already had the intuition to sample quadrature nodes proportionally to the square of some determinant, motivated by the fact that it gave good interpolation nodes. For a discussion of their results in modern parlance, see (Gautier, Bardenet, and Valko, 2019c) or Section 3.3 below. Revisiting this connection with interpolation with Pierre Chainais and our PhD student Ayoub Belhadji, we came up with a dual approach to [START_REF] Bardenet | Monte Carlo with Determinantal Point Processes[END_REF]. Instead of trying to find a DPP yielding a good estimator for as wide a function class as possible, we asked the following question: if I give you a DPP with kernel K, is that a good interpolating set of nodes for some function space, and what functions is this DPP supposed to integrate well? The approach naturally connects with recent work on quadrature in reproducing kernel Hilbert spaces [START_REF] Bach | On the equivalence between kernel quadrature rules and random feature expansions[END_REF], and gave [START_REF] Belhadji | Kernel quadrature with DPPs[END_REF](Belhadji et al., , 2020a)); see Section 3.4.

To sum up, this chapter is based on our works [START_REF] Bardenet | Monte Carlo with Determinantal Point Processes[END_REF][START_REF] Belhadji | Kernel quadrature with DPPs[END_REF]Gautier et al., 2019b;Belhadji et al., 2020a).

A central limit theorem from random matrix theory

Much of random matrix theory deals with the behaviour of the eigenvalues of certain random matrices, as can be gathered from the reference book [START_REF] Anderson | An introduction to random matrices[END_REF]. These eigenvalues are highly correlated random variables, and they can have surprising properties to the scientist more used to independence. The following example is due to [START_REF] Johansson | On random matrices from the compact classical groups[END_REF].

Let U N ⊂ M N (C) be the group of N ×N unitary matrices, that is,

A ∈ U N ⇔ A ĀT = ĀT A = I N . There is a unique probability measure µ H on U N such that µ H (A) = µ H ({ga; a ∈ A}) (3.1)
for any g ∈ U N and any Borel set A ⊂ U N . This measure is called the Haar measure on U N , see for instance (Anderson et al., 2010, Theorem F.13). In other words, the Haar measure µ H plays on U N the role of the Lebesgue measure on the additive group R, so that µ H is often thought of as the natural equivalent of the uniform probability distribution. Now, let U ∼ µ H , and since the eigenvalues of a unitary matrix are complex numbers of unit modulus, let us further denote the N eigenvalues of U by their arguments [START_REF] Johansson | On random matrices from the compact classical groups[END_REF] used a result called the strong Szegő theorem to prove that if g : [0, 2π] → R is in L 1 , and its Fourier coefficients satisfy

(θ i ) ∈ [0, 2π] N .
σ 2 2 ∞ k=1 k|ĝ k | 2 < ∞, (3.2) then N i=1 g(θ i ) -N 2π 0 g(θ) dθ 2π → N (0, σ 2 ), (3.3)
where the convergence is in distribution when N → ∞. The central limit theorem (CLT) in Equation (3.3) is surprising by its lack of normalization: the standard CLT for i.i.d. variables would have a factor √ N in the left-hand side. In other words, if you think of

I N N -1 N i=1 g(θ i
) as a Monte Carlo estimator of (2π) -1 g(θ)dθ, then (3.3) implies that the mean square error of I N is asymptotically equivalent to 1/N 2 . This contrasts with traditional Monte Carlo variance of order 1/N .

Efficient cancellation happens in the variance of the left-hand side of (3.3). Actually, the arguments (θ i ) of U form a DPP on [0, 2π], namely the circular unitary ensemble defined in Example 2. One way to get some intuition for the fast convergence is as follows. The eigenvalues of U tend to be very regularly spaced on the unit circle compared to uniform samples, as if they "repelled" each other; see Figure 2.2 If we think of the eigenvalues as a random grid on which the signal g is measured, repulsiveness is akin to the grid being rigid. Rigidity has two positive features. First, the number of points that fall in any fixed Borel subset of [0, 2π] does not vary much from one realization of U to the other, hence the small variance of g(θ i ). Second, g essentially contains low frequencies by assumption (3.2), so that repulsiveness ensures both that g(θ i ) has small bias and that we make the most of the smoothness of g by spreading out samples (θ i ).

So how do we take a result like Johansson's in (3.3) and make it a general numerical integration result? In other words, we want to reverse-engineer (3.3) and assess whether repulsive random variables as seen in random matrix theory can be a general device for variance reduction in Monte Carlo computation. There are two issues to consider. First, eigenvalues of random matrices are bound to be in R or C, while we would like to integrate over R d for arbitrary d 1. Second, we would like to integrate against an arbitrary measure, whereas random matrix results such as (3.3) implicitly choose one measure (here, uniform over [0, 2π]). We add to these requirements that we need to preserve both the fast convergence compared to traditional Monte Carlo methods, and the interpretability of the asymptotic variance (3.2) as a measure of the decay of the Fourier coefficients of the integrand.

From random matrices to numerical integration

Let µ be a positive Borel measure supported in [-1, 1] d . We want a quadrature

f dµ ≈ N i=1 w i f (x i ) (3.4)
for a wide class of functions f .

Multivariate orthogonal polynomial ensembles

When d = 1, letting ϕ k in (2.6) be the orthonormal polynomials w.r.t. µ results in a DPP called an orthogonal polynomial ensemble (OPE, [START_REF] König | Orthogonal polynomial ensembles in probability theory[END_REF]). When d > 1, orthonormal polynomials can still be uniquely defined by applying the Gram-Schmidt procedure to a set of monomials, provided the base measure is not pathological. However, unlike for d = 1, there is no natural order on multivariate monomials, so we pick an order b : N → N d before we apply Gram-Schmidt to the first N monomials in that order. In [START_REF] Bardenet | Monte Carlo with Determinantal Point Processes[END_REF] we use the graded lexicographic order, defined by ordering multi-indices (α 1 , . . . , α d ) by their maximum degree max α i , and for constant maximum degree, by the usual lexicographic order. There is some limited freedom to change the order in our proofs.

Denoting by Φ k the corresponding multivariate orthonormal polynomials, then by multivariate OPE we mean the DPP with base measure µ(dx) = ω(x)dx and kernel

K OPE (x, y) = Φ k (x)Φ k (y)
, where the sum runs over multi-indices k such that 0 b(k) N -1. This is a projection kernel, and thus a valid DPP kernel; see Section 2.2.1. Note that for a separable base measure, i.e. if

ω(x) = ω 1 (x 1 ) . . . ω d (x d ), (3.5)
where x = (x 1 , . . . , x d ) ∈ R d , multivariate orthonormal polynomials are products of univariate ones. In that case,

K OPE (x, y) = N -1 b(k)=0 d j=1 ϕ j k j (x j )ϕ j k j (y j ), (3.6)
where (ϕ j k ) k 0 are the orthonormal polynomials w.r.t. ω j .

A stochastic Gaussian quadrature

Let µ be separable as in (3.5), and {x 1 , . . . , x N } be drawn from the multivariate OPE given by (3.6). If we integrate (2.7) N -1 times, we obtain that any point in our DPP has marginal pdf K OPE (x, x)ω(x). For f ∈ L 1 (dµ), we deduce that

ÎN N i=1 f (x i ) K OPE (x i , x i ) (3.7)
is an unbiased estimator of f (x)ω(x)dx. Note that (3.3) does not contain nodedependent weights, since the kernel of the CUE has constant diagonal, K CUE (x, x) ∝ 1.

Note also that the estimator (3.7) should be familiar to users of Gaussian quadrature. [START_REF] Gauss | Methodus nova integralium valores per approximationem inveniendi[END_REF] indeed showed that if you restrict to d = 1, and replace the DPP samples in (3.7) by the zeros of ϕ N , then ÎN is deterministic and has zero error when f is a polynomial of degree less than 2N -1. This is remarkable as ϕ N has degree N and thus only N zeros, while integrating polynomials of degree twice that. In a sense, our estimator (3.7) is a stochastic relaxation of Gaussian quadrature.

To further visualize the similarity with Gaussian quadrature, we plot in Figure 3.1 three samples. The left plot depicts an i.i.d. sample from a product measure with Jacobi marginals. The marginals are plotted in green on each axis. The middle plot contains the Cartesian product of two Gaussian quadratures, one on each axis with the corresponding green marginal as target measure. Each node x i is depicted with a circle, the area of which is proportional to the node's weight 1/K OPE (x i , x i ) in (3.7) and in Gaussian quadrature. The right plot shows a sample from the corresponding multivariate OPE, with the same weights. Well-spread points and large weights accumulate in the bulk of ω in both the middle and right plots.

Our generalization of Johansson's CLT

Theorem 2 [START_REF] Bardenet | Monte Carlo with Determinantal Point Processes[END_REF] Let µ(dx) = ω(x)dx with ω separable, C 1 , positive on the open set (-1, 1) d , and satisfying a technical regularity assumption. If x 1 , . . . , x N stands for the associated multivariate OPE, then for every g that is C 1 and vanishing outside [-1 + ε, 1ε] d for some ε > 0,

N 1+1/d ÎN -g(x)µ(dx) law ----→ N →∞ N 0, Ω 2 g,ω ,
where

Ω 2 g,ω = 1 2 ∞ k 1 ,...,k d =0 (k 1 + • • • + k d ) gω ω ⊗d eq (k 1 , . . . , k d ) 2 , (3.8)
and

ω eq (x) = π -1 (1 -x 2 ) -1/2 .
Looking back at our requirements in Section 3.1, we have extended (3.3) to an arbitrary dimension d, and an arbitrary choice of target measure. We have preserved the fast convergence of (3.3), although the gain w.r.t. the usual Monte Carlo rate decreases with dimension, and we have preserved the interpretability of the asymptotic variance as a measure of decay of the Fourier coefficients of the integrand, including the target pdf. The additional ω eq term is related to the projection of the circle onto the x-axis: ω eq is the "marginal" distribution of the uniform distribution on the circle. As a last comment, we have a CLT that essentially only assumes that the integrand is C 1 , while typical fasterthan-Monte-Carlo algorithms like quasi-Monte Carlo [START_REF] Dick | Digital nets and sequences: discrepancy theory and quasi-Monte Carlo integration[END_REF] require the smoothness of the integrand to grow with d. To sum up, the fact that the repulsiveness in our DPP is tailored to the integration problem at hand allows for weaker assumptions, and gives a fast CLT with interpretable asymptotic variance.

An importance sampling version

One downside of Theorem 2 is that in practice, one rarely can evaluate the orthogonal polynomials w.r.t. µ. Indeed, this would imply that we know all moments of µ, which is unrealistic in, say, Bayesian inference, where access to µ is limited to pointwise evaluation up to a multiplicative constant [START_REF] Robert | Monte Carlo statistical methods[END_REF]. In [START_REF] Bardenet | Monte Carlo with Determinantal Point Processes[END_REF], Theorem 2.9), we show that if x 1 , . . . , x N is the multivariate OPE associated to a pdf q that satisfies the assumption of Theorem 2, then

ĨN N i=1 g(x i ) K OPE (x i , x i ) ω(x i ) q(x i )
satisfies the same CLT as in Theorem 2, with the same asymptotic variance. In other words, you do not incur any cost, asymptotically, for having replaced ω by another reference distribution. This is highly unusual from a Monte Carlo perspective, where you would expect the asymptotic variance to contain a term that measures the mismatch between q and ω. From an orthogonal point of view, this is a consequence of the universality of Chebyshev polynomials: whenever the underlying measure satisfies some weak assumptions (technically, is Nevai-class; see [START_REF] Bardenet | Monte Carlo with Determinantal Point Processes[END_REF], the corresponding orthogonal polynomials.

Soviets already did it

In a classical twist, we later realized that Soviet mathematicians were actually the first to consider drawing quadrature nodes from projection DPPs. Fifteen years before [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF] even defined DPPs, [START_REF] Ermakov | Polynomial Approximations and the Monte-Carlo Method[END_REF] had already proven the following result, which we restate in modern jargon.

Theorem 3 [START_REF] Ermakov | Polynomial Approximations and the Monte-Carlo Method[END_REF] Consider f ∈ L 2 (µ), and N functions ϕ 0 , . . . , ϕ N -1 ∈ L 2 (µ) orthonormal w.r.t. µ. Let {x 1 , . . . , x N } be drawn from the DPP with reference measure µ and projection kernel (2.6). Consider the linear system

   ϕ 0 (x 1 ) . . . ϕ N -1 (x 1 ) . . . . . . ϕ 0 (x N ) . . . ϕ N -1 (x N )       y 1 . . . y N    =    f (x 1 ) . . . f (x N )    . (3.9)
Then the solution of (3.9) is unique, µ-almost surely, with coordinates

y k = det Φ ϕ k-1 ,f (x 1:N ) det Φ(x 1:N ) ,
where Φ ϕ k-1 ,f (x 1:N ) is the matrix obtained by replacing the k-th column of Φ(x 1:N ) by f (x 1:N ). Moreover, for all 1 k = j N , the solution vector satisfies Of course, back then, [START_REF] Ermakov | Polynomial Approximations and the Monte-Carlo Method[END_REF] did not have the DPP machinery that we now have in our toolbox to study, e.g., the asymptotic fluctuations of their estimator, but Theorem 3 is still good food for thought. We investigated some the consequences of Theorem 3 in (Gautier, Bardenet, and Valko, 2019c).

E[y k ] = f, ϕ k-1 , Var[y k ] = f 2 - N -1 =0 f, ϕ 2 ,
For starters, it provides unbiased estimators of the N Fourier-like coefficients f, ϕ k , k N -1, and these N estimators are uncorrelated, with the same variance. Moreover, the faster the decay of the Fourier-like coefficients, the smaller the variance. In particular, if f is in the span of the eigenfunctions of the kernel, the variance is zero. In the setting of multivariate OPEs, the first (in the order prescribed by b) orthonormal polynomial ϕ 0 is constant, equal to µ [-1, 1] d -1/2 . Hence, a direct application of Theorem (3) yields that

I N (f ) µ [-1, 1] d 1/2 det Φ ϕ 0 ,f (x 1:N ) det Φ(x 1:N ) (3.11)
is an unbiased estimator of [-1,1] d f dµ. Developing the numerator w.r.t. the column that contains evaluations of f , I N (f ) is seen to be a quadrature estimator

I N (f ) = N i=1 w EZ i (x 1 , . . . , x N )f (x i ),
which we term the EZ estimator after [START_REF] Ermakov | Polynomial Approximations and the Monte-Carlo Method[END_REF]. Unlike the variance of our Gaussian quadrature estimator ÎN (f ) in (3.7), the variance of I N (f ) clearly reflects the accuracy of the approximation of f by its projection onto Span(ϕ 1 , . . . , ϕ N ).

In particular, it allows to interpolate and perfectly integrate polynomials up to degree b -1 (N -1).

Nonetheless, the limiting theoretical properties of I N (f ), like a CLT, look hard to establish. In particular, unlike our estimator ÎN (f ) in (3.7), the EZ estimator is not a linear statistic of the underlying DPP, since each weight w EZ i (x 1 , . . . , x N ) depends on all nodes. Thus, it does not fit the assumptions of traditional CLTs for DPPs [START_REF] Soshnikov | Determinantal random point fields[END_REF]. More general CLTs are available [START_REF] Błaszczyszyn | Limit theory for geometric statistics of point processes having fast decay of correlations[END_REF], but only for DPPs with sufficiently short-range correlation.

Experiments in (Gautier et al., 2019c) suggest that the EZ estimator is an excellent replacement for I N (f ) when the residual of the projection of f onto Span(ϕ 1 , . . . , ϕ N ) is small. In other words, if you can interpolate the integrand with a few basis functions, then you should form a projection kernel using these basis functions, draw from the corresponding DPP, and use the EZ estimator. However, the performance of EZ rapidly decreases once the integrand spreads even a small part of its norm on all functions ϕ k , k > N . In particular, the variance of the EZ estimator on a non-polynomial but still C ∞ function is seen in (Gautier et al., 2019c) to quickly fall down to the classical O(1/N ) as d grows, while the variance of our estimator I N (f ) steadily decreases as the predicted rate O(1/N 1+1/d ) from Theorem 2.

Interpolation and quadrature in RKHSs

To take the best of both our stochastic Gaussian quadrature estimator [START_REF] Bardenet | Monte Carlo with Determinantal Point Processes[END_REF] and the interpolative quadrature of [START_REF] Ermakov | Polynomial Approximations and the Monte-Carlo Method[END_REF], it is natural to ask what functions can be well interpolated using their evaluation on a DPP sample. Intuitively, since repulsiveness in DPPs is measured by a kernel, one would like to consider functions, the smoothness of which is measured by the same kernel. This is naturally embodied by requiring the function to belong to a reproducing kernel Hilbert space (RKHS), for which interpolation and integration are tightly linked.

Like the EZ estimator of Theorem 3, our quadrature scheme shall come from interpolating a certain function on a DPP sample. But the interpolated function will not directly be the integrand, rather a representative of the target measure in the space of functions. Finally, relating the RKHS kernel (governing the smoothness of the integrand) and the DPP kernel (governing the repulsiveness of the points) will lead us to an algorithm that realizes tight interpolation and quadrature rates in RKHSs (Belhadji, Bardenet, andChainais, 2019, 2020b).

Reproducing kernel Hilbert spaces

Let k : X × X → R + be a positive definite kernel, and consider the set of real linear combinations of the so-called feature functions

F 0 = Span k(x, •), x ∈ X .
Now let F be the completion of F 0 for the norm derived from k(x, •), k(y, •) := k(x, y). The resulting F is a Hilbert space of functions such that

f, k(x, •) = f (x),
(3.12) for all f ∈ F. F is called the reproducing kernel Hilbert space of kernel k (RKHS; [START_REF] Berlinet | Reproducing kernel Hilbert spaces in probability and statistics[END_REF]. Two famous examples are the RKHS with Gaussian kernel, and Sobolev spaces of sufficiently high smoothness; see again [START_REF] Berlinet | Reproducing kernel Hilbert spaces in probability and statistics[END_REF].

The RKHS kernel k itself might not be a valid DPP kernel. One natural, "PCAflavoured" strategy to build a valid DPP kernel from k is to write k in diagonal form (also called a Mercer decomposition), and build a projection kernel onto the span of some of its eigenfunctions, as in Section 2.2.1. To do that, we need a few assumptions. Still assuming that X is equipped with a Borel measure µ, we further assume that it is supported on all of X, and that X k(x, x)dµ(x) < +∞.

(3.13)

As a consequence, we get that the embedding operator I F : F → L 2 (µ) is compact and injective [START_REF] Steinwart | Support Vector Machines[END_REF], and that the operator

Σf (•) = X k(•, y)f (y)dµ(y), f ∈ L 2 (µ) (3.14)
is self-adjoint, positive semi-definite, and trace-class [START_REF] Simon | Trace Ideals and Their Applications[END_REF]. Thus, there exists an orthonormal family (e m ) of L 2 (µ) formed by eigenfunctions of Σ. The basis is ordered by requiring that the corresponding eigenvalues (σ m ) are non-decreasing. Our assumptions imply a Mercer-type decomposition of k,

k(x, y) = m∈N * σ m e m (x)e m (y), (3.15)
where N * = N \ {0} and the convergence is pointwise [START_REF] Steinwart | Mercer's theorem on general domains: on the interaction between measures, kernels, and RKHSs[END_REF]. Moreover, for m ∈ N * , we write e F m = √ σ m e m . Since I F is injective, (e F m ) m∈N * is an orthonormal basis of F [START_REF] Steinwart | Mercer's theorem on general domains: on the interaction between measures, kernels, and RKHSs[END_REF]. Unless explicitly stated, we assume that F is dense in L 2 (dµ), so that (e m ) m∈N * is an orthonormal basis of L 2 (dµ). For more intuition, under our assumptions,

f ∈ F ⇔ m σ -1 m f, e m 2 L 2 (dω) < ∞.
(3.16)

Again, this shows that the smoothness in an RKHS is governed by the kernel, in the sense that the modulus | f, e m | of the Fourier-like coefficients of f ∈ F has to decrease fast enough compared to the eigenvalues of the kernel to make (3.16) converge.

3.4.2 Quadrature in a RKHS with a projection DPP Now let

K(x, y) = N m=1
e m (x)e m (y) (3.17) be the projection kernel onto the first N eigenfunctions of the RKHS kernel k. Unlike k, the projection kernel K is a valid DPP kernel, so that we can consider quadrature nodes x = {x 1 , . . . , x N } drawn from the projection DPP with reference measure µ and kernel (3.17). To get a quadrature scheme, it remains to define the weights in (3.4).

To this end, we note that, given f ∈ F and g ∈ L 2 (µ), Cauchy-Schwarz in F yields

X f g dµ - i∈[N ] w i f (x i ) f F × ν g - N i=1 w i k(x i , •) F , (3.18)
where ν g = X g(x)k(x, .)dµ(x)= Σg is the so-called embedding of g in the RKHS F. From now on, we fix an arbitrary g ∈ L 2 (µ).

By (3.18), an upper bound on the approximation error of ν g in F implies an upper bound on the integration error of f g that is uniform over any bounded subset of F. This observation has sparked intense research on the kernel approximation of embeddings ν g . On our end, we note that thanks to the reproducing property (3.12), the approximation error in the RHS of (3.18) rewrites

ν g - N i=1 w i k(x i , •) F = ν g 2 -2 N i=1 w i ν g (x i ) + N i,j=1
w i w j k(x i , x j ).

(3.19)

If the nodes x 1 , . . . , x N are fixed and in general position, i.e. det((k(x i , x j ))) > 0, then (3.19) entails that there is a unique set of weights that minimizes the upper bound (3.18), namely

   ŵ1 . . . ŵN    = k(x i , x j ) -1    ν g (x 1 ) . . . ν g (x N )    . (3.20)
Note that these optimal weights also guarantee that N i=1 ŵi k(x i , •) interpolates ν g at the nodes x i .

Theorem 4 [START_REF] Belhadji | Kernel quadrature with DPPs[END_REF] Under the assumptions of Section 3.4, let (x 1 , . . . x N ) be drawn from the projection DPP with reference measure µ and kernel K in (3.17). Then det((k(x i , x j )) > 0 almost surely, so that we can consider the optimal weights (3.20).

Let

r N = m N +1 σ m , then E sup g µ 1 ν g - N i=1 w i k(x i , •) 2 F 2σ N +1 + 2N 2 r N . (3.21)
This version of the theorem is actually a slight improvement over our the original result in [START_REF] Belhadji | Kernel quadrature with DPPs[END_REF]; see the PhD manuscript of [START_REF] Belhadji | Subspace sampling using determinantal point processes[END_REF]. Note that I present a uniform bound in g for simplicity, but we have more precise g-dependent bounds; see (Belhadji, 2020, Chapter 4) for the sharpest version. The important thing to notice here is that if σ n decreases fast enough, say geometrically as for the Gaussian kernel, then both terms in the bound are of the same order. But if the spectrum decays polynomially, like for Sobolev RKHSs, then N 2 r N dominates. In experiments, however, we have actually observed that the expected approximation error decays as σ N +1 for any classical kernel, so that it should be possible to prove a bound in O(σ N +1 ). While we haven't been able to prove it under a DPP, we have obtained this faster rate for a mixture of projection DPPs called volume sampling.

Quadrature in an RKHS with continuous volume sampling

In Section 3.4.2, we built a DPP kernel K from the RKHS kernel k. If we momentarily forget about DPPs, it is also natural to define a (permutation invariant) distribution over nodes that simply measures the squared volume in F of the parallelotope built onto the feature functions k(x i , •), i.e. For N ∈ N * , we have

x 1 , . . . , x N ∼ 1 Z det k(x i , x j dµ(x 1 ) . . . dµ(x N ). ( 3
det k(x i , x j ) ∝ |U |=N u∈U σ u 1 N ! det K U (x i , x j ) . (3.24)
Each component in the mixture (3.24) is a projection DPP onto a choice of N eigenfunctions of k, with weight proportional to the product of the corresponding eigenvalues. The component with the largest weight is thus precisely the DPP that we use in Section 3.4.2. We can thus expect similar performances for kernel quadrature, at least when the spectrum of k has a sharp cutoff. Actually, we can even improve on the bound in Theorem 4.

Theorem 5 (Belhadji, Bardenet, and Chainais, 2020b) Under the assumptions of Section 3.4, let (x 1 , . . . x N ) be drawn from the volume sampling distribution (3.22) with reference measure µ and kernel k. Then det k(x i , x j > 0 almost surely, so that we can consider the optimal weights (3.20). Moreover,

sup g µ 1 E ν g - N i=1 ŵi k(x i , •) 2 F σ N (1 + β N ) , (3.25)
where

β N = min M ∈[2:N ] [(N -M + 1)σ N ] -1 m M σ m .
We note that β N is bounded as soon as the spectrum (σ n ) decays polynomially, so that the bound in Theorem 5 is quite generic. It is actually the first such generic bound on interpolation and quadrature in RKHSs: give it any kernel, the theorem gives you a sharp bound. Indeed, a result by (Pinkus, 2012, Theorem 2.2, Chapter 4) shows that the worst-case (in g) approximation error is lower-bounded by a constant times σ N +1 ; in that sense, our bound is sharp. Note also that, unlike Theorem 4, the supremum is outside of the expectation in (3.25); this comes from the fact that a direct computation of the expected squared error is possible under volume sampling, so that bounds appear later in the proof. We refer to the discussion in (Belhadji et al., 2020b) for details, and numerical experiments that illustrate Theorem 5.

Closing remarks

I think that the bound in Theorem 5 is a solid milestone. It demonstrates what DPPs and mixtures thereof can reach in terms of quadrature and interpolation error in a favourable setting. By favourable, I mean that the integrands live in a Hilbert space, that their smoothness is controlled by the spectral properties of a known kernel, and that the performance measure is the residual of an orthogonal projection.

Indeed, in Theorems 4 and 5, the weighted sum 1 i N ŵi k(x i , •) is the orthogonal projection of ν g onto the subspace Span(k(x i , •), 1 i N ) formed by linear combinations of feature functions anchored at our quadrature design. In words, the setting is favourable to DPPs because sampling quadrature nodes proportionally to the determinant of det K(x i , x j ) in Theorem 4 guarantees both that the functions k(x i , •) are almost orthogonal in F, and that their span is close to that of the first N eigenfunctions of the integration operator Σ, thus capturing a large part of the norm of any f ∈ F. More technically, the alignment between the two spans is controlled by the so-called principal angles between the two subspaces. Taking expectations under our DPP yields tractable upper bounds on symmetric quantities of these angles; see the detailed discussion in the thesis of (Belhadji, 2020, Section 4.4). In particular, there is room for improvement in controlling non-symmetric quantities of the principal angles, which would likely lead to a similar sharp bound as for volume sampling.

Rather than controlling principal angles, volume sampling allows for a more direct approach. The proof of Theorem 5 works by decomposing the interpolation error onto the eigenfunctions of the integration operator Σ, with the coefficients giving a prominent role to the continuous equivalent of the leverage scores from regression analysis. Volume sampling controls these leverage scores in expectation, thus yielding our Theorem 5. I don't think there is much room for improvement here.

My current belief is that DPPs and volume sampling are the right way to perform quadrature as long as the target functions are in an RKHS. The randomization helps in two ways. Algorithmically, on the one hand, one can sample the nodes in polynomial time as long as one knows the Mercer decomposition of the kernel, while maximizing the kind of Gramian determinants appearing in the densities of our nodes is a hard non-convex optimization problem. Mathematically, on the other hand, as often in Monte Carlo versus, say, quasi-Monte Carlo [START_REF] Dick | Digital nets and sequences: discrepancy theory and quasi-Monte Carlo integration[END_REF], it is easier to derive and interpret bounds on the mean square error of an estimator than on the error of a deterministic quadrature. I believe that Theorem 5 is a solid translation of these beliefs, and that we should now turn to methodology. In terms of Bayesian inference applications, for instance, it is not clear what role the function g ∈ L 2 (µ) in Theorem 5 should play: should it be part of the loss function, part of the posterior pdf? How can we use that degree of freedom to further optimize the constants in the (g-dependent) version of our bound?

Similarly, in Bayesian practice, we rarely know the kernel of an RKHS to which our integrands belong. Can we learn one as we sample nodes using some MCMC sampler, and progressively turn on repulsiveness in sampling? What is a parametric class of kernels that is large enough to capture the integrands in a particular Bayesian application? How robust are our bounds to misspecifying the kernel? Reversing the viewpoint, can we force the integrand to belong to a pre-specified RKHS by limiting the choices of the modeler when designing her likelihood, so that she can trade off modeling and integration accuracy?

On the other hand, how do we actually sample from the point processes involved in our theorems, if we don't have access to the Mercer decomposition of the kernel? The advantage of volume sampling over any projection DPP like (3.17) is that evaluating the density (3.22) of volume sampling does not require a spectral decomposition, but only pointwise evaluations of the kernel. MCMC approaches thus seem natural [START_REF] Rezaei | A polynomial time MCMC method for sampling from continuous determinantal point processes[END_REF], but one would ideally need a sampler that reaches an error in total variation of the order of σ N in a reasonable number of iterations, say quadratic in N . I believe that this requires a dedicated sampler that uses the kernel in a smart way.

Chapter 4

Feature selection with determinantal point processes In Chapter 3, we saw how projection DPPs can select a finite number of feature functions among {k(x, •), x ∈ X}, so as to yield a small average reconstruction error of some target function belonging to a Hilbert space. We shall see here the same general principle in action. The role of the feature functions is played by the columns of a large matrix X ∈ R N ×d , d N , of which we want to select s d columns. This is the general problem of feature selection, or variable selection, in statistics and learning.

One common motivation is to reduce dimensionality while preserving the interpretability of the selected features, unlike principal component analysis (PCA), which defines new mixed features. If the dimensionality reduction is meant as a preprocessing for, say, a linear regression only on the p selected columns, then one potential performance metric is the excess risk of regressing onto the p columns, compared to the original regression vector. However, if one wants to have a performance metric for selecting a subset S of column indices that is independent of the regression method, another natural choice is some matrix norm X -Π ν S X ν (4.1) of the residual of the projection of X onto the span of the selected columns. The index ν in (4.1) denotes the matrix norm that is used in measuring the size of the residual and in defining the projection Π ν S ; see Section 4.1. If the residual is measured in terms of Frobenius or spectral norm, i.e. ν ∈ {Fr, 2}, then the rank-k projection Π that minimizes the norm of the residual X -ΠX ν is the projection Π k onto the subspace spanned by the first k principal components of X.

After quickly surveying the state-of-the-art in column subset selection in Sections 4.2 and 4.3, we shall first show in Section 4.4 that a suitably chosen projection DPP yields a reconstruction error (4.1) in both Frobenius and spectral norm that is as close as possible to the PCA gold standard. In that sense, DPPs yield a variable selection that is as efficient as PCA, while preserving the semantics of the original features, since -unlike PCA-the new features (columns) are a subset of the original ones. Then, we will show that the same DPP allows to control the excess risk in linear regression; see Section 4.5. As in Chapter 3, we shall finish by underlining the limitations of the approach and current work.

Finally, this chapter is based on part of Ayoub Belhadji's PhD thesis, in particular the paper (Belhadji, Bardenet, and Chainais, 2020a), to which I refer the reader for details. On top of all proofs and illustrative experiments, the paper contains a detailed discussion of past work on column-subset selection that is of independent interest.

Apples and oranges: column subsets and PCA

The result of column subset selection will usually be compared to the result of projecting onto the first k principal components of X, and in this chapter, we reserve the notation k ∈ N for that reference size. We denote by Π k X the best rank-k approximation to X. The sense of best can be understood in either Frobenius or spectral norm, as both give the same result. But this coincidence is lost when projecting onto the span of a subset of columns of X, and some care must be taken in defining the quantities to be compared in a fair way.

For a given subset S ⊂ {1, . . . , d} of size |S| = s and ν ∈ {2, Fr}, let

Π ν S,k X = arg min A X -A ν ,
where the minimum is taken over all matrices A = X :,S B such that B ∈ R s×d and rkB k; in words, the minimum is taken over matrices of rank at most k that lie in the column space of C = X :,S . When |S| = k, we simply write Π ν S X = Π ν S,k X. In practice, the Frobenius projection can be computed as Π Fr S X = CC + X, yet there is no simple expression for Π 2 S X. However, Π Fr S X can be used as a proxy for Π 2 S X since

X -Π 2 S X 2 X -Π Fr S X 2 √ 2 X -Π 2 S X 2 , (4.2)
see (Boutsidis et al., 2011, Lemma 2.3).

Multinomial sampling with and without SVD

Note that because my overall narrative is about explaining when and where DPPs can help, I take it for granted that we want to apply randomized algorithms. Randomized approaches to numerical linear algebra tend to be more scalable, and come with more interpretable error bounds. For instance, naively optimizing a criterion like X -Π S,k ν X over subsets of {1, . . . , d} requires exhaustive enumeration. This is prohibitive when d 1, while randomized approaches typically have polynomial costs. Less naive deterministic column subset selection methods have been investigated, in particular rank-revealing QR factorizations in the vein of [START_REF] Golub | Numerical methods for solving linear least squares problems[END_REF], but even there, randomized approaches [START_REF] Boutsidis | An improved approximation algorithm for the column subset selection problem[END_REF] can match the performance and cost of state-of-the-art factorizations [START_REF] Gu | Efficient algorithms for computing a strong rank-revealing QR factorization[END_REF] in spectral norm, while also providing sharp guarantees in Frobenius norm.

The first randomized algorithms to have been investigated are naturally independent draws from the set of columns (Deshpande and Vempala, 2006;[START_REF] Drineas | Relative-error cur matrix decompositions[END_REF][START_REF] Boutsidis | Near optimal column-based matrix reconstruction[END_REF]. More specifically, a prototypical example corresponds to sampling s columns of X i.i.d. from a multinomial distribution of parameter p ∈ R d . This parameter p can be the squared norms of each column [START_REF] Drineas | Clustering large graphs via the singular value decomposition[END_REF], for instance, or the more subtle k-leverage scores [START_REF] Drineas | Relative-error cur matrix decompositions[END_REF], defined as follows.

Definition 4 (k-leverage scores) Let X = U ΣV T be the singular value decomposition of the original matrix X. The k-leverage scores, denoted by k 1 , . . . , k d , are the squared Euclidean norms of the d columns of the first k rows of V T . This definition is further illustrated in Figure 4.2. Geometrically, k i is the squared cosine of the angle between the i-th vector of the canonical basis of R d and the k-principal subspace, i.e., the subspace spanned by the vectors with canonical coordinates the first k columns of V . We refer the reader to (Belhadji et al., 2020a, Section 3.4) for a detailed discussion of k-leverage scores and principal angles. It thus makes intuitive sense to select columns marginally according to their k-leverage score if one wants to approximate the k-principal subspace. Note however in passing that there may be redundancy here, and that it would be even more natural to select k columns jointly, according to how close their span is to the k-principal subspace. This is precisely what the DPP in Section 4.4 shall do.

The k-leverage scores are expensive to evaluate, since they call for a truncated SVD of order k, but they come with the best known bound for multinomial sampling on the ratio of their expected approximation error over that of PCA.

Theorem 6 (Drineas et al., 2007, Theorem 3)

Let ε ∈ (0, 1]. If the number s of sampled columns satisfies s 3200k 2 ε 2 , (4.3)
then, under i.i.d. sampling from the k-leverage scores distribution, 2007) also considered replacing multinomial with Bernoulli sampling, still using the k-leverage scores. The expected number of columns needed for (4.4) to hold is then lowered to O( k log k ε 2 ). A natural question is then to understand how low the number of columns can be, while still guaranteeing a multiplicative bound like (4.4). A partial answer has been given by Deshpande and Vempala (2006).

P X -Π Fr S,k X 2 Fr (1 + ε) X -Π k X 2
Proposition 2 (Deshpande and Vempala, 2006, Proposition 4) d+1) such that for any S ⊂ {1, . . . , d},

Given ε > 0, k, d ∈ N such that dε 2k, there exists a matrix X ε ∈ R kd×k(
X ε -Π Fr S,k X ε 2 Fr (1 + ε) X ε -X ε k 2 Fr . (4.5)
This suggests that a lower bound for the number of columns is 2k/ε, at least in the worst case sense of Proposition 2. Interestingly, the k-leverage scores distribution of the matrix X ε in the proof of Proposition 2 is uniform, so that k-leverage score sampling boils down to simple uniform sampling.

Repulsive sampling needs fewer columns

A multiplicative bound on the expected approximation error using only k columns was obtained by [START_REF] Deshpande | Matrix approximation and projective clustering via volume sampling[END_REF] using so-called volume sampling, which we introduced in Example 7.

Theorem 7 (Deshpande et al., 2006) Let S be a random subset of {1, . . . , d}, chosen with probability

P(S) ∝ det(X T :,S X :,S ) 1 {|S|=k} . (4.6) Then E X -Π Fr S X 2 Fr (k + 1) X -Π k X 2 Fr (4.7)
and

E X -Π 2 S X 2 2 (d -k)(k + 1) X -Π k X 2 Fr . (4.8)
Sampling according to (4.6) can be done in polynomial time; see [START_REF] Deshpande | Efficient volume sampling for row/column subset selection[END_REF] or the mixture argument for k-DPPs in Chapter 2. Using a worst case example, Deshpande et al. (2006) proved that the k + 1 factor in (4.7) cannot be improved.

Proposition 3 (Deshpande et al., 2006) Let ε > 0. There exists a (k + 1) × (k + 1) matrix X ε such that for every subset S of k columns of X ε ,

X ε -Π Fr S X ε 2 Fr > (1 -ε)(k + 1) X ε -Π k X ε 2
Fr .

(4.9)

4.4 DPPs take the best of leverage scores and repulsiveness

In (Belhadji et al., 2020a), our contribution is to introduce the DPP on {1, . . . , d} with kernel

K = V :,1:k V T :,1:k , (4.10)
and reference measure the counting measure on [d] := {1, . . . , d}. Note that this kernel matrix is a projection (K 2 = K) and has rank k, so that the corresponding DPP is a projection DPP, and samples S thus satisfy |S| = k almost surely; see Chapter 2. We are thus on par with volume sampling in Theorem 7, in the sense that we sample only k columns. We shall see that our DPP reaches the (k + 1) factor of volume sampling in Theorem 7, and even improves over volume sampling under a certain sparsity assumption.

To see why such a DPP is natural, note that the diagonal of K is formed by the k-leverage scores of Definition 4. Since, by definition (2.5) of the one-point correlation function of a DPP, the diagonal of a DPP kernel is proportional to the marginal probability of inclusion of each element, the DPP with kernel K is thus an extension of k-leverage score sampling. Moreover, the DPP enforces diversity among the selected columns, in the sense of favouring sets of k columns S ⊂ {1, . . . , d} such that the k feature vectors V i,1:k , i ∈ S, span a large volume in R k . Note that the square of this volume is also the product of the squared cosines of the principal angles between the span of the selected columns and the k-principal subspace; see (Belhadji et al., 2020a, Appendix C). Methodologically, this makes sure that we have a set of columns that jointly approximate the k-principal subspace. Mathematically, this is all the more interesting that performance metrics that compare residuals of projections to that of PCA involve the same set of principal angles. This is precisely the same argument as for the DPP we used in Section 3.4.2 to perform quadrature in an RKHS. For this reason, the proofs of our results follow the same abstract layout. First, we decompose the error as the contribution of the principal subspace and a mismatch term that measures the alignment between the selected columns and that principal subspace. Second, we show that the mismatch term is bounded by a quantity that is symmetric in the principal angles between the two subspaces; taking expectations, we get a simple closed-form bound.

To present our results, we first denote the number of nonzero k-leverage scores by p d, and we quantify the decay of the singular values of X by the flatness parameter

β = σ 2 k+1   1 d -k j k+1 σ 2 j   -1 . (4.11)
In words, β ∈ [1, dk] measures the flatness of the spectrum of X above the cut-off at k. Indeed, (4.11) is the ratio of the largest term in a mean to that mean. The closer β is to 1, the more similar the terms in the sum in the denominator of (4.11). At the extreme,

β = d -k when σ 2 k+1 > 0 while σ 2 j = 0, ∀j k + 2.
Proposition 4 (Belhadji et al., 2020a) Under the DPP of kernel (4.10),

E X -Π 2 S X 2 2 (1 + k(p -k)) X -Π k X 2 2 (4.12)
and

E X -Π Fr S X 2 Fr 1 + β p -k d -k k X -Π k X 2 Fr . (4.13)
The bound in (4.12) compares favorably with volume sampling (4.8) since the dimension d has been replaced by the sparsity level p. For β close to 1, the bound in (4.13) is better than the bound (4.7) of volume sampling since (pk)/(dk) 1. Again, the sparser the k-leverage scores, the smaller the bounds. Finally, if needed, bounds in high probability easily follow from Proposition 4 using Markov's inequality. Now, one could argue that, in practice, sparsity is never exact: it can well be that p = d while there still are a lot of small k-leverage scores. We introduce two ideas to further tighten the bounds of Proposition 4. First, we define an effective sparsity level in the vein of [START_REF] Papailiopoulos | Provable deterministic leverage score sampling[END_REF] as follows. Let π be a permutation of {1, . . . , d} such that leverage scores are reordered

k π 1 k π 2 ... k π d . (4.14)
For δ ∈ {1, . . . , d}, let T δ = {π δ , . . . , π d }. Let θ 1 and

p eff (θ) = min    q ∈ {1, . . . , d} i q k π i k -1 + 1 θ    . (4.15)
Note that by definition, the k-leverage scores sum to k, so that the partial sum defining p eff (θ) should be close to the total sum. Second, we condition the DPP on a favourable event with controlled probability.

Theorem 8 (Belhadji et al., 2020a) Let A θ be the event {S ∩ T p eff (θ) = ∅}. Then, under the DPP of kernel (4.10), on the one hand,

P (A θ ) 1 θ , (4.16)
and, on the other hand,

E X -Π 2 S X 2 2 A θ (1 + (p eff (θ) -k + 1)(k -1 + θ)) X -Π k X 2 2 (4.17) and E X -Π Fr S X 2 Fr A θ 1 + β (p eff (θ) + 1 -k) d -k (k -1 + θ) X -Π k X 2 Fr . (4.18)
In Theorem 8, the effective sparsity level p eff (θ) replaces the sparsity level p of Proposition 4. The key is to condition on S not containing any index of column with too small a k-leverage score, that is, the event A θ . In practice, this is achieved by rejection sampling: we repeatedly and independently sample S ∼ DPP(K) until S ∩ T p eff (θ) = ∅. The caveat of any rejection sampling procedure is a potentially large number of samples required before acceptance. But in the present case, Equation (4.16) guarantees that the expectation of that number of samples is less than θ. The free parameter θ thus interestingly controls both the "energy" threshold in (4.15), and the complexity of rejection sampling. The approximation bounds suggest picking θ close to 1, which implies a compromise with the value of p eff (θ) that should not be too large either. We have empirically observed that the performance of the DPP is relatively insensitive to the choice of θ. Finally, we have observed that the bound in Theorem 8 is representative of the actual behaviour of the mean error, and that many real datasets from natural language processing had sparse leverage scores, so that our DPP is a practical improvement over volume sampling in that setting; see the numerical experiments in (Belhadji et al., 2020a, Section 6).

A guarantee on the excess risk in sketched linear regression

So far, we have focused on approximation bounds in spectral or Frobenius norm for the residual X -Π ν S,k X. This is a reasonable generic measure of error as long as it is not known what the practitioner wants to do with the submatrix X :,S . In this section, we assume that the ultimate goal is to perform linear regression of some y ∈ R N onto X. Other measures of performance then become of interest, such as the excess risk incurred by regressing onto X :,S rather than X. We use here the framework of [START_REF] Slawski | On principal components regression, random projections, and column subsampling[END_REF], further assuming well-specification for simplicity.

For every i ∈ {1, . . . , N }, assume y i = X i,: w * + ξ i , where w * is some ground truth vector of regression coefficients, and the noise variables ξ i are i.i.d. real variables with mean 0 and variance v. For a given estimator w = w(X, y), its excess risk is defined as

E(w) = E ξ Xw * -Xw 2 2 N . (4.19)
In particular, it is easy to show that the ordinary least squares (OLS) estimator ŵ = X + y has excess risk

E( ŵ) = v × rk(X) N . (4.20)
For convenience, for S ⊂ {1, . . . , d} of cardinality k, we introduce the matrix S ∈ {0, 1} d×k such that AS = A :,S for any A ∈ R N ×d . Selecting k columns indexed by S in X prior to performing linear regression yields w S = (XS) + y ∈ R k . We are interested in the excess risk of the corresponding sparse vector

ŵS := Sw S = S(XS) + y ∈ R d ,
which has all coordinates zero, except those indexed by S.

Proposition 5 (Theorem 9, Mor-Yosef and Avron, 2019

) Let S ⊂ [d], such that |S| = k. Let (θ i (S)) i∈{1,.
..,k} be the principal angles between SpanS and SpanV k , see (Belhadji et al., 2020a, Appendix C). Then

E( ŵS ) 1 N 1 + max i∈{1,...,k} tan 2 θ i (S) w * 2 σ 2 k+1 + vk N . (4.21)
Compared to the excess risk (4.20) of the OLS estimator, the second term of the righthand side of (4.21) replaces rankX by k. But the price is the first term of the right-hand side of (4.21), which we loosely term bias. To interpret this bias term, we first look at the excess risk of the principal component regressor

w * k ∈ arg min w∈SpanV :,k E ξ y -Xw 2 /N . (4.22)
Proposition 6 (Corollary 11, Mor-Yosef and Avron, 2019)

E(w * k ) w * 2 σ 2 k+1 N + vk N . (4.23)
The right-hand side of (4.23) is almost that of (4.21), except that the bias term in the CSS risk (4.21) is larger by a factor that measures how well the subspace spanned by S is aligned with the principal eigenspace V k . This makes intuitive sense: the performance of CSS will match PCR if selecting columns yields almost the same eigenspace. Again, everything is prepared for DPPs to yield a small bound in expectation, since the right-hand side of (4.21) is a symmetric function of the principal angles between the span of the selected columns and the k-principal subspace. Jumping directly to the effective sparsity setting of Theorem 8, we use the same rejection sampling trick, and arrive at the following.

Theorem 9 (Belhadji et al., 2020a) Using the notation of Theorem 8 for leverage scores and their indices, it holds

E E( ŵS ) A θ 1 + k -1 + θ p eff (θ) -k + 1 w * 2 σ 2 k+1 N + vk N . (4.24)
The price of interpretable dimensionality reduction w.r.t. principal component regression is thus the same O(1 + k(pk)) prefactor as in (4.17) in front of the PCA bias term from (4.23). The price goes down as the k-leverage scores get sparser. Finally, to the best of our knowledge, bounding the excess risk in linear regression has not been investigated under volume sampling.

Closing remarks

Seeing a projection DPP as a random subspace with an alignment constraint is maybe the most intuitive viewpoint to understand how to build a kernel for a given subsampling application, in my opinion, as long as the target performance measure is somehow a function of that alignment. By expressing alignment in terms of principal angles, it becomes obvious that the determinant of the Gram matrix somehow controls how close a data-aligned subspace (the span of selected columns in this chapter, the span of the feature functions k(x, •) anchored at x a quadrature node in Chapter 3) is to a predefined subspace (the k-principal subspace in this chapter, the subspace of polynomials of bounded degree in Chapter 3, the span of the first eigenfunctions of the RKHS kernel in Section 3.4). We refer to [START_REF] Kassel | Determinantal probability measures on Grassmannians[END_REF] for a more abstract take on this geometric view of DPPs.

In particular, we have been using the same alignment technique to provide a supervised version of the DPP-based column selection presented in this chapter; see [START_REF] Mai | Un processus ponctuel déterminantal pour la sélection de variables supervisée[END_REF] for preliminary results. The idea is to sample columns from a DPP which aligns the span of the selected columns with a Krylov subspace, which can be interpreted as a supervised counterpart to PCA.

At a lower level, sampling our PCA-based DPP requires computing cross-leverage scores, i.e. our kernel, which requires running a PCA. If we could lower that cost, our algorithms would become more than a replacement for PCA. There are workarounds to computing leverage scores, like estimating them by random projections [START_REF] Drineas | Fast approximation of matrix coherence and statistical leverage[END_REF][START_REF] Boutsidis | Near optimal column-based matrix reconstruction[END_REF]. There are also double-phase techniques, where one first preselects marginally relevant columns before finding a diverse subset of the selected ones. The latter family of methods has given algorithms that are extremely efficient in practice [START_REF] Boutsidis | An improved approximation algorithm for the column subset selection problem[END_REF]), yet come with loose bounds. It would be interesting to find the right way to pre-screen items that does not impact too much the theoretical performance of a DPP; see e.g. [START_REF] Dereziński | Exact sampling of determinantal point processes with sublinear time preprocessing[END_REF] for results on DPP sampling that relate to this idea.

Chapter 5

Time-frequency filtering with zeros of GAFs " Silence fell like a hammer made of feathers. It left holes in the shape of the sound of the sea. When I arrived in Lille at the beginning of 2015, I joined a signal processing team. From coffee break to lunch time at the restaurant universitaire with my excellent colleagues, I traded long monologues on DPPs against a crash course in signal processing. In particular, I developed an interest in time-frequency and time-scale analysis listening to Pierre Chainais and Julien Flamant. When we all went to GRETSI'15 in Lyon, the main French-speaking event in signal and image processing, I saw a plenary talk by Patrick Flandrin, who presented his recent interest in the zeros of spectrograms [START_REF] Flandrin | Time-frequency filtering based on spectrogram zeros[END_REF], as opposed to the traditional focus on spectrogram maxima.

As an aside, spectrograms are a cornerstone of time-frequency analysis [START_REF] Flandrin | Time-frequency/time-scale analysis[END_REF]. They are quadratic time-frequency representations of a signal f : R → C (Gröchenig, 2001, Chapter 4), associating to each time and frequency a real number that measures the energy content of a signal at that time and frequency, unlike global-in-time tools such as the Fourier transform. Since it is natural to expect that there is more energy where there is more information or signal, most methodologies have focused on detecting and processing the local maxima of the spectrogram; see e.g. [START_REF] Flandrin | Time-frequency/time-scale analysis[END_REF]. Usual techniques include ridge extraction, e.g., to identify chirps, reassignment and synchrosqueezing, to better localize the maxima before further quantitative analysis. Using the zeros of a spectrogram instead of higher-level sets thus may sound counter-intuitive.

Yet, at GRETSI'15, Patrick Flandrin explained that the locations of the zeros of a spectrogram in the time-frequency plane almost completely characterize the spectrogram, and he proposed to use the point pattern formed by the zeros in tasks such as filtering and reconstruction of signals in noise. This proposition stems from the empirical observation that the zeros of the spectrogram of white noise are uniformly spread over the timefrequency plane, and tend not to cluster, as if they repelled each other. In the presence of a signal, zeros are absent in the time-frequency support of the signal, thus creating large holes that appear to be very rare when observing pure white noise. This leads to testing the presence of signal by looking at statistics of the point pattern of zeros, and trying to identify holes. Patrick's talk connected to all my pet subjects of the time: time-frequency analysis, repulsive point processes, and even reproducing kernel Hilbert spaces. This was the trigger for what has become a large part of my research programme.

The similarity of the mathematical backbone of DPPs with the analytical tools behind time-frequency analysis [START_REF] Gröchenig | Foundations of time-frequency analysis[END_REF] made me conjecture during GRETSI'15 that the point process formed by the zeros of the spectrogram of white Gaussian noise was a DPP. As we shall see in Section 5.1 and detailed in [START_REF] Bardenet | On the zeros of the spectrogram of white noise[END_REF], this is not the case; the zeros rather follow the distribution of the zeros of the planar Gaussian analytic function introduced in Section 2.3. By the time we finished working on [START_REF] Bardenet | On the zeros of the spectrogram of white noise[END_REF], I was invited to give a talk at GRETSI'19 on DPPs and signal processing. How frustrating that my flagship time-frequency example was not a DPP! However, once all the maths were laid out to characterize the zeros of random spectrograms, I realized that I could have obtained a well-known DPP if I had replaced the Hermite functions by the inverse Fourier transform of the Laguerre functions in the proofs. A crucial discussion with Patrick Flandrin helped me realize that my modified proof corresponded to characterizing the zeros not of the classical spectrogram of white noise, but of the analytic wavelet transform [START_REF] Daubechies | Time-frequency localization operators -a geometric phase space approach: Ii. the use of dilations[END_REF]) of a particular white noise. This is sketched in Section 5.2. A few days later, Adrien Hardy, passing though my office for some other reason (most likely for one of his usual pranks), listened to my musings of the moment and immediately realized the key role implicitly played by generating functions of orthogonal polynomials. This allowed us to further generalize our results on the identification of zeros to more time-frequency transforms, and cleanly define all sorts of white noises [START_REF] Bardenet | Time-frequency transforms of white noises and Gaussian analytic functions[END_REF]. Section 5.3 gives this construction, of which Sections 5.1 and 5.2 can be seen as special cases. In Section 5.4, I show how this generic construction led us to define a new time-frequency transform, acting on finitedimensional vectors [START_REF] Bardenet | Time-frequency transforms of white noises and Gaussian analytic functions[END_REF][START_REF] Pascal | A covariant, discrete time-frequency representation tailored for zero-based signal detection[END_REF]. The topic of zeros of spectrograms has now gained momentum, and there are many open questions, some of which I list in Section 5.5.

Finally, while this chapter may seem out of my initial narrative for this manuscript (i.e., repulsive point processes provide small residuals of orthogonal projections in expectation), I show in Section 5.6 that this is not the case, and that we are somehow pursuing the same goal as in Chapters 3 and 4, only in a different form.

For reference, this chapter is based on our works [START_REF] Bardenet | On the zeros of the spectrogram of white noise[END_REF][START_REF] Bardenet | Time-frequency transforms of white noises and Gaussian analytic functions[END_REF], and [START_REF] Pascal | A covariant, discrete time-frequency representation tailored for zero-based signal detection[END_REF].

From the STFT to the planar GAF

Let f ∈ L 2 (R) and g be the pdf of a centered Gaussian, with variance normalized so that its L 2 norm is g 2 = 1. The short-time Fourier transform (STFT) V g f of f with Gaussian window g, is defined by

V g f (u, v) = f (t)g(t -u)e -2iπtv dt = f, M v T u g , (u, v) ∈ R 2 , (5.1) with •, • denoting the inner product in L 2 (R), M v f = e 2iπv• f (•) and T u f = f (• -u).
The spectrogram of f is simply the squared modulus of its STFT, i.e.

u, v → |V g f (u, v)| 2 .
From (5.1), it is intuitive that the spectrogram is large in (u, v) if frequency v is contributing to the signal around time u. An example spectrogram is shown in Figure 5.1(a). The signal is a so-called linear chirp, as can be seen guessed from the linear relation between time and frequency. The (numerical) zeros of the spectrogram look well-spread, with a slight accumulation around the signal support, and a complete absence in the support, which the spectrogram takes large values.

To investigate the zeros of the STFT of white noise, our first step was to rigorously define the STFT of complex white Gaussian noise (WGN), in a manner that both satisfied the intuition we may have of white noise, and allowed us to speak of the zeros of the spectrogram. In (Bardenet et al., 2018, Section 3.1), we used a classical definition of real WGN as a random tempered distribution with Gaussian characteristic function. A complex WGN ξ is then defined as having independent white Gaussian noises ξ 1 , ξ 2 as its real and imaginary parts.

This allowed us to pointwise compute ξ, M v T u g since M v T u g is smooth. By pointwise, we mean for a given time u and frequency v. Then in (Bardenet et al., 2018, Proposition 3), we stitched all these pointwise evaluations together and concluded that with probability 1 on the white noise, the zeros of the function where a k = ( ξ 1 , h k + i ξ 2 , h k )/ √ 2 are i.i.d. unit complex Gaussians, and h k are the Hermite functions. To read this chapter, it is enough to know that Hermite functions are a basis of L 2 (R) built using orthogonal polynomials with respect to the Gaussian g.

z = u + iv → ξ, M v T u g are those of so-called planar GAF, GAF C (z) := ∞ k=0 a k π k/2 z k √ k! , (5.2)
All the above was the formalisation of the following heuristic. A complex WGN should have unit complex Gaussian coefficients in any basis of L 2 (R). Among such bases, we pick the Hermite basis, because Hermite functions essentially map to complex monomials through the STFT [START_REF] Gröchenig | Foundations of time-frequency analysis[END_REF]. We thus expect an analytic function of the form (5.2) as the STFT of ξ.

Identifying the zeros of the STFT of white noise with the zeros of a random entire function on the time-frequency plane is satisfactory, since entire functions have isolated zeros. It thus makes sense to speak of the point process of the zeros. Furthermore, the zeros of (5.2) are actually a well-studied stationary point process [START_REF] Hough | Zeros of Gaussian analytic functions and determinantal point processes[END_REF], called the zeros of the planar Gaussian analytic function (GAF), and which we introduced in Section 2.3. We show a sample of these zeros as white dots in Figure 5.2(a), with the heatmap showing the corresponding spectrogram realization. 1 Combining known probabilistic results on this point process and standard statistical methodology for spatial point processes, we obtained in (Bardenet et al., 2018, Section 5) powerful signal detection procedures that strengthen the seminal work of [START_REF] Flandrin | Time-frequency filtering based on spectrogram zeros[END_REF]. Without entering into details, spatial statistics provides a sophisticated toolbox for one-sample testing with stationary point processes [START_REF] Baddeley | Spatial point patterns : methodology and applications with R[END_REF]. In Figure 5.1(b), the beige pixels correspond to the centers of unusually large Euclidean balls that do not contain any zero of the spectrogram. By unusual, I mean that under the hypothesis that we are observing the zeros of the spectrogram of pure white noise, the corresponding empty balls appear with very small probability. These empty balls have caused to reject the hypothesis of pure white noise. One can then apply any heuristic to extract an envelope of these beige points; see Figure 5.1(c), and invert the STFT after multiplication with a smooth indicator function of the support. This yields a denoised, or "filtered", signal; see Figure 5.1(d).

To conclude, one can show (Hough et al., 2009, Section 5.1) that the zeros of the planar GAF are not a DPP with Hermitian kernel. Indeed, the pair correlation function of the zeros of the planar GAF is known in closed form, and it takes values above 1. For a DPP with Hermitian kernel, this is forbidden; see Section 2.2.5. Strictly speaking, they could still be a DPP with non-Hermitian kernel.

From a DPP to analytic wavelets

We were slightly frustrated of not finding a well-known DPP behind the STFT of white noise in [START_REF] Bardenet | On the zeros of the spectrogram of white noise[END_REF]. In particular, for α > -1, consider the random analytic function

GAF (α) D (z) := k∈N a k Γ(k + α + 1) k! z k , |z| < 1, (5.3)
1 Actually, this specific figure was obtained through a different time-frequency transform giving the same distribution; see [START_REF] Bardenet | Time-frequency transforms of white noises and Gaussian analytic functions[END_REF]. where the a i are again i.i.d. unit complex Gaussians. The series (5.3) defines the socalled hyperbolic GAF [START_REF] Hough | Zeros of Gaussian analytic functions and determinantal point processes[END_REF]. When α = 0, it simply reads a k z k and its zeros are a DPP. A natural 2 question is then to know what signal processing problem we should have tried to solve in order to end up with the zeros of the hyperbolic GAF instead of the planar GAF. The gamma coefficients are a hint that the Hermite functions should somehow be replaced by Laguerre functions. This backward reasoning is the starting point of [START_REF] Bardenet | Time-frequency transforms of white noises and Gaussian analytic functions[END_REF].

The same heuristic as in Section 5.1 holds: we need a basis of some Hilbert space of signals, a white noise that can be decomposed onto that basis, and a time-frequency transform that sends elements of the basis onto the right monomials.

The space of signals

Laguerre functions are a basis of L 2 (R + ) built using the orthogonal polynomials with respect to the pdf of a gamma random variable. Thus, we are looking for a space of signals where the half-line R + plays a role. A natural3 example is the space of "analytic" signals

H 2 (R) := {f ∈ L 2 (R, C) : Supp( f ) ⊂ R + }.
The inverse Fourier transforms (f k ) of the Laguerre functions form a basis of H 2 (R).

The time-scale transform

Now, we need a transform that maps the basis (f k ) to monomials. Let β > -1/2 and set

ψ β (t) := 1 (t + i) β , t ∈ R.
The Fourier transform of ψ β is essentially a gamma pdf

1 √ 2π R + x β e -x e itx dx = - 1 √ 2π (-i) β+1 (t + i) β+1 Γ(β + 1),
so that ψ β is intuitively a fine window if we want to use Laguerre calculus in the Fourier domain. The Daubechies-Paul wavelet transform of f ∈ H 2 (R) is defined as (5.4) where u ∈ R and s ∈ R * + are thought of as time and scale, T u f (t) := f (tu) and D s f (t) := s -1/2 f (t/s) are the usual translation and dilation operators.

W β f (u, s) := f, T u D s ψ β ,
We show in (Bardenet and Hardy, 2019, Theorem 2.3) that for f ∈ H 2 (R) and up to a non-vanishing term,

W β f (-u, s) ∝ L (β) D f (ϕ(u + is)), (5.5)
where u + is ∈ C + is in the time-scale half-plane, ϕ(w) := (2wi)/(2w + i) is a conformal mapping from C + to D, and L (β) D is a transform that precisely sends the inverse Fourier transforms of the Laguerre functions onto the right monomials

L (β) D f k (z) = Γ(k + 2β + 1) k! z k . (5.6)
Everything is in place but the white noise. Heuristically, if we were able to somehow decompose a white noise ξ on H 2 (R + ) as k 0 a k f k , where a k are i.i.d. unit complex Gaussians, then (5.5) and (5.6) should help us conclude that

W β ξ(z) ∝ k 0 a k Γ(k + 2β + 1) k! z k = GAF (2β) D (z)
up to non-vanishing terms. This would characterize the zeros of the analytic wavelet transform of white noise.

The white noise

The series k 0 a k f k almost surely diverges in H 2 (R). To remedy this, we follow [START_REF] Gross | Abstract Wiener Spaces[END_REF] in his formalisation of abstract Wiener spaces. First, we define a slightly weaker norm on H 2 (R) that ensures the series converges

f 2 Θ := k∈N 1 1 + k 2 f, f k 2 .
(5.7) Indeed, it follows that (5.8) since the right-hand side has finite expectation. Now define Θ as the completion of H 2 (R) for the norm (5.7). Leonard Gross' construction [START_REF] Gross | Abstract Wiener Spaces[END_REF] then gives a unique probability measure on Θ, called white Gaussian noise, such that Θ e iλ(θ) µ(dθ) = e -h λ 2 /2 , λ ∈ Θ * , (5.9)

k∈N a k f k 2 Θ = k∈N |ξ k | 2 1 + k 2 < ∞ a.s.
where Θ * is the topological dual of Θ and

h λ ∈ H 2 (R) is associated to λ ∈ Θ * ⊂ H 2 (R) by
Riesz' representation theorem. Equation (5.9) is to be interpreted as specifying the Fourier transform of the white noise on H 2 (R). In particular, it leads to a k f k converging almost surely in Θ to µ, which is our goal. We can even give rates of convergence (Bardenet and Hardy, 2019, Section 5). We also note in passing that µ is a rigorous definition of an analytic white noise.

With the white noise now defined, we prove in (Bardenet and Hardy, 2019, Theorem 2.3) that L Finally, we note that [START_REF] Abreu | Filtering with wavelet zeros and Gaussian analytic functions[END_REF] have also independently proved that the zeros of W β ξ are the zeros of the hyperbolic GAF, with a slightly different definition of white noise.

Mapping signals to GAFs: a generic template

Having extended the STFT case to the analytic wavelet transform, we then realized that the mathematical construction in Section 5.2 was much more general than these two cases. The whole idea is to first fix a Hilbert space H of signals, fix an orthonormal basis (f k ) of H, and complete H into Θ as in Section 5.2.3, so as to be able to define a white Gaussian noise on Θ. Then pick a sequence (Ψ k ) k∈N of holomorphic functions on a domain Λ satisfying sup

z∈K k∈N |Ψ k (z)| 2 < ∞ (5.10)
for any compact K ⊂ Λ. Now define a "time-frequency" transform L acting on H by

L f (z) := ∞ k=0 f k , f Ψ k (z), z ∈ Λ.
(5.11)

In signal processing terms, the signal f ∈ H is analyzed by L in the basis (f k ), and reconstructed as an analytic function on Λ using the dictionary Ψ k . Note that dominated (5.12)

The construction of white noise in Section 5.2.3 guarantees that L extends uniquely to Θ and that L ξ has the same law as the GAF a k Ψ k . This is the main result of [START_REF] Bardenet | Time-frequency transforms of white noises and Gaussian analytic functions[END_REF].

The STFT with Gaussian window and the analytic wavelet transform of [START_REF] Daubechies | Time-frequency localization operators -a geometric phase space approach: Ii. the use of dilations[END_REF] are particular cases of L , where the representation (5.12) is provided by identities for the generating functions T of classical orthogonal polynomials. This prompted us to plug other generating functions for orthogonal polynomials as kernels, and check whether the resulting transform had a time-frequency/scale interpretation.

From the spherical GAF to the Kravchuk transform

Armed with the general pattern in Section 5.3, it is easy to find, given a GAF, a space H of signals, an analysis basis (f k ) and a reconstruction family (ψ k ) to map white Gaussian noise in H to any GAF. To the signal processer, the interesting question is whether the resulting transform has any time-frequency interpretation.

Take for instance the spherical GAF

GAF (N ) S (z) := N k=0 ξ k N k z k , z ∈ C.
(5.13)

We plot in Figure 5.4 a sample of the spherical GAF and its zeros, mapped onto a sphere through stereographic projection. This rendering highlights the invariance of the set of zeros to isometries of the sphere, just like the zeros of the planar and hyperbolic GAFs are invariant to isometries of the plane and the Poincaré half-plane, respectively. Following Section 5.3, we propose to take H = C N the space of finite signals, along with the basis of Kravchuk functions. In words, Kravchuk functions are to the binomial distribution what Hermite functions are to the Gaussian distribution, and what Laguerre functions are to the gamma distribution. Taking

ψ k (z) = N k z k
then leads to a transform L in (5.11) that maps the white noise on C N to the spherical GAF. Moreover, the kernel T in (5.12) of this transform has a closed form, due again to a generating identity for Kravchuk polynomials; see the discussion in [START_REF] Bardenet | Time-frequency transforms of white noises and Gaussian analytic functions[END_REF]). Yet, it is not easy to find a time-frequency interpretation of such a transform. To that end, in [START_REF] Pascal | A covariant, discrete time-frequency representation tailored for zero-based signal detection[END_REF], we further modified the Kravchuk transform L to T y(z) = (1 + |z| 2 ) -N × L y(z).

(5.14)

Note that T y is not analytic anymore, but that it has the same zeros as L y. Moreover, T is an isometry in the sense that 5.15) where (ϑ, ϕ) ∈ [0, π] × [0, 2π] are the spherical coordinates parametrizing the sphere S 2 , dµ(ϑ, ϕ) = sin(ϑ)dϑdϕ is the uniform measure on S 2 , and we write T y(ϑ, ϕ) for T y evaluated at the stereographic projection z = cot(ϑ/2)e iϕ . More importantly, using the formalism of coherent states, we show in [START_REF] Pascal | A covariant, discrete time-frequency representation tailored for zero-based signal detection[END_REF] that T is covariant under the action of the group of rotations SO(3), meaning that 5.16) where R u (resp. R u ) denotes a particular action of a rotation parametrized by the unitary vector u ∈ R 3 on vectors of size N + 1 (resp. on points of the unit sphere). We refer the reader to the supplementary material to [START_REF] Pascal | A covariant, discrete time-frequency representation tailored for zero-based signal detection[END_REF] for a concise introduction to the representation theory of SO(3). This covariance property (5.16) is meaningful to the signal processer. The classical STFT is covariant with respect to translations in the time-frequency plane [START_REF] Gröchenig | Foundations of time-frequency analysis[END_REF]. In words, translating a signal in time and frequency before computing the STFT yields the same result as evaluating the STFT of the original signal at a translated time and frequency. Even if the new variables (ϑ, ϕ) do not have as intuitive a meaning as time and frequency, they are still variables in a phase space (the sphere, instead of the plane), and there is a group of isometries (rotations, instead of time-frequency translations, a.k.a.

y 2 2 = (4π) -1 S 2 |T y(ϑ, ϕ)| 2 dµ(ϑ, ϕ), (
T [R u y](ϑ, ϕ) = T y(R u (ϑ, ϕ)), (
the Weyl-Heisenberg group), with respect to which the Kravchuk transform is covariant. That makes the Kravchuk transform close in spirit to the classical short-time Fourier transform.

Opening remarks

We close this section with a few thoughts, and remarks on ongoing work.

Removing the analyticity constraint

There are some fundamental limitations to the correspondence between GAFs and timefrequency transforms. The three GAFs we introduced -planar, hyperbolic, and sphericalare canonical in the sense that they are the only GAFs that are invariant to isometries of the complex plane, the Poincaré half-plane, and the sphere. From a statistical point of view, such invariances are crucial, as they make the point process of the zeros stationary and allow estimating functional statistics, a key step in designing detection/reconstruction algorithms (Bardenet et al., 2018, Section 5). Extending the correspondence will require either more exotic geometries, or dropping the Gaussian/analytic aspects on the GAF side.

Another point in favour of removing analyticity is that the Gaussian window is the only window that makes the STFT map L 2 (R) to analytic functions [START_REF] Ascensi | Model space results for the Gabor and Wavelet transforms[END_REF]. When the window is not Gaussian but still Schwartz, the STFT maps to polyanalytic functions, and [START_REF] Haimi | Zeros of Gaussian Weyl-Heisenberg functions and hyperuniformity of charge[END_REF] have studied the zeros of the STFT of white noise in that case. Interestingly, the analytic wavelet transform of Section 5.2 is also the only continuous wavelet transform that maps to analytic functions [START_REF] Holighaus | Characterization of analytic wavelet transforms and a new phaseless reconstruction algorithm[END_REF]. Besides working with polyanalytic functions, discrete transforms such as the Kravchuk transform [START_REF] Bardenet | Time-frequency transforms of white noises and Gaussian analytic functions[END_REF][START_REF] Pascal | A covariant, discrete time-frequency representation tailored for zero-based signal detection[END_REF] do bypass constraints on the window while still mapping to analytic functions.

Making the Kravchuk transform a practical tool

We believe that the Kravchuk transform of Section 5.4 has the potential of becoming a standard item in the practitioner's toolbox. When describing the classical STFT, we have swept under the rug the fact that any practical computation is done by approximating the involved integrals by Riemann sums, and then applying fast Fourier transforms. When Riemann sums are a poor approximation to the corresponding integrals, say when the sampling frequency is lower than typical frequencies in the signal, the numerical zeros of the spectrograms we compute in practice are likely to be far from the zeros of the underlying GAFs. In contrast, the Kravchuk transform takes as inputs finite signals, i.e. (finite-length) vectors of complex numbers. If we accept to work with the sphere as a phase space, the Kravchuk transform avoids any sophisticated approximation procedure of functions by finite sequences of numbers. We refer to [START_REF] Pascal | A covariant, discrete time-frequency representation tailored for zero-based signal detection[END_REF] for numerical experiments that demonstrate that signal detection tests for small signals are more powerful when using the Kravchuk transform than the STFT. This situation is reminiscent of what the finite Fourier transform is to the continuous Fourier transform, where replacing the latter by the former requires some intellectual gymnastics (e.g., consider only band-limited functions, use Shannon's sampling theorem, and truncate the coefficients) and edge effects that are not easy to quantify. We are also currently working with Barbara Pascal and Julien Flamant on fast algorithms to compute and invert the Kravchuk transform, leveraging analogies with spin spherical harmonics.

The point process of maxima

Maxima of spectrograms are also the zeros of a Gaussian process, but with non-analytic samples. Abreu (2022) shows that the maxima still form a point process, are stationary, and that their intensity if 1/3 that of zeros, in line with a conjecture of [START_REF] Flandrin | Time-frequency filtering based on spectrogram zeros[END_REF]. A question of practical importance is if, when, and how there is an interest in using maxima, zeros, or both for a given signal processing task. Preliminary experiments suggest that maxima alone are more efficient than zeros alone at, say, estimating the instantaneous frequency of a simple signal with large signal-to-noise ratio and little interference with other signal components. On the contrary, the more rigid structure of zeros (due to being the zeros of an analytic function) makes them more robust than maxima at low signal-tonoise ratio, or when there is strong interference caused by multiple signals with overlapping time-frequency supports. This is the topic of ongoing work with my former postdoc Arnaud Poinas, and my current postdoc Juan Manuel Miramont, along with co-supervisors Pierre Chainais and François Auger; see [START_REF] Miramont | A public benchmark for denoising and detection methods[END_REF] for preliminary results on a numerical benchmark where we compare maxima-based and zero-based detection and reconstruction.

Ridges and level sets

For estimation of instantaneous frequencies, it is actually customary to work with the maxima not of the spectrogram, but of each restriction of the spectrogram to a constanttime line in the time-frequency plane; sort of synchrosqueezed maxima. This gives curves of maxima, so-called ridges, and the random set of ridges of a (random) spectrogram would be an interesting object to study. Similarly, many geometric properties of the spectrogram could be used for filtering or instantaneous frequency estimation. [START_REF] Ghosh | Signal analysis via the stochastic geometry of spectrogram level sets[END_REF], for instance, showed that level sets of spectrograms can be used to detect signals made of single Hermite functions with a quantifiable testing power.

Despite appearances, this chapter fits in my narrative

My narrative is that functions can be efficiently reconstructed if sampled at carefully chosen repulsive point processes. In particular, functions that belong to a reproducing kernel Hilbert space (RKHS), sampled at a DPP or a mixture of DPPs, were shown to yield small residuals in expectation under a suitable DPP and a mixture of DPPs in Chapter 3. We are in a very similar setting here, and we could even make the settings coincide.

Let me consider the STFT V g from (5.1), with Gaussian window g for concreteness. The range V g (L 2 ) of the STFT with Gaussian window is an RKHS -the so-called Bargmann-Fock space of quantum physics. Second, time-frequency filtering is usually carried out by projection onto a subset of time-frequency atoms, just like we approximated mean elements by sums of feature functions in Chapter 3. Indeed, as we saw in Figure 5.1, filtering works by localizing a signal, i.e. (i) isolating a subset Ω ⊂ R 2 of the time-frequency plane (the assumed "support" of the signal, delimited in green in Figure 5.1(c)), and (ii) reconstructing a signal f through the restriction of the STFT of the (noisy) signal f to that subset,

f = V -1 g [1 Ω (u, v)V g f (u, v)] .
(5.17)

In words, we select the (typically uncountable) subset of time-frequency atoms indexed by Ω, and approximate the original signal by a linear combination of these atoms. The proposition of [START_REF] Flandrin | Time-frequency filtering based on spectrogram zeros[END_REF] can be paraphrased as defining the subset Ω as a function of the pattern of zeros of V g (f ). Just like in Chapters 3 and 4, the span of the timefrequency atoms indexed by Ω is indeed a random subset that is defined through a repulsive point process. Three differences with Chapter 3 are that the selected random subset is uncountable, that the localization operator in (5.17) is not an orthogonal projection, and that the involved point processes are not DPPs.

There are a number of steps that we could take to prove a result like those in Chapter 3 and 4 here. First, we could use the zeros (u i , v i ) of the spectrogram of white noise to design a (countable) random dictionary (M v i T u i g) of time-frequency atoms, and prove that it gives a (random) frame with good properties [START_REF] Gröchenig | Foundations of time-frequency analysis[END_REF]. Second, we could study the frame given by a DPP in the time-frequency plane with kernel that of the Bargmann-Fock space V g (L 2 ). It turns out that this DPP is precisely the Ginibre ensemble introduced in Example 3. Note that [START_REF] Ghosh | Determinantal processes and completeness of random exponentials: the critical case[END_REF] gave a clever indirect proof that the set of timefrequency atoms indexed by the Ginibre ensemble is a complete set; a natural question is thus to know whether it is a frame with tighter constants than the independent draws of the seminal [START_REF] Bass | Relevant sampling of band-limited functions[END_REF]. Playing with the spectrum of the Ginibre kernel, say using the truncated Ginibre ensemble of Example 3, we should also end up with a finite set of time-frequency atoms with properties resembling that of the infinite Ginibre frame. Such a result would be of practical importance, and is the subject of ongoing collaboration with Pierre Chainais and our former PhD student Ayoub Belhadji.

On algorithmic aspects of repulsive point processes, we have focused on sampling algorithms. Indeed, for a projection DPP of rank N , the cost of the standard algorithm is at least cubic in N ; see Chapter 2. Like for other kernel machines, this cubic cost typically becomes prohibitive as soon as the required N goes beyond 10 4 . There are several avenues around this bottleneck. On our side, we have worked on approximate sampling through Markov chain Monte Carlo for specific DPPs in the thesis of Guillaume Gautier [START_REF] Gautier | On sampling determinantal point processes[END_REF]. In particular, we have proposed an MCMC algorithm to sample from finite projection DPPs [START_REF] Gautier | Zonotope hit-and-run for efficient sampling from projection DPPs[END_REF], which we have experimentally demonstrated to mix remarkably fast. We have proposed a similar MCMC kernel for βensembles [START_REF] Gautier | Fast sampling of β-ensembles[END_REF]. Finally, with my postdoc Michaël Fanuel, we have studied a physics-inspired DPP on the edges of a graph [START_REF] Fanuel | Sparsification of the regularized magnetic Laplacian with multi-type spanning forests[END_REF], for which exact sampling is possible without relying on the "base times height" formula of Chapter 2. This DPP is a natural generalization of the edges of a uniform spanning tree, a well-studied DPP [START_REF] Lyons | Probabilities on trees and networks[END_REF]. Our generalization replaces the Laplacian of the graph by a so-called connection Laplacian, and has applications in, e.g., ranking. The reader can also check our recent preprint [START_REF] Jaquard | Smoothing complex-valued signals on graphs with Monte-Carlo[END_REF] for an application of that connection-graph DPP to graph signal processing.

On statistical aspects, we have worked first on inference for DPPs, where the latter are used as models for repulsive behaviour. Repulsive point processes are natural models for, e.g., the location of trees in a forest, which compete for nutrients, or sentences in a summary, whose diversity guarantees coverage. Among repulsive point processes, DPPs with Hermitian kernels stand out as exceptionally tractable: sampling is polynomial, and all correlation functions are available in closed form, which allows using them in inference procedures. In particular, we have proposed MCMC and variational algorithms for Bayesian inference in parametric DPP models [START_REF] Bardenet | Inference for determinantal point processes without spectral knowledge[END_REF]. More recently, building on a novel representation for non-negative functions in RKHSs [START_REF] Marteau-Ferey | Non-parametric models for non-negative functions[END_REF], we have proposed the first non-parametric algorithm to learn a continuous kernel from DPP samples [START_REF] Fanuel | Nonparametric inference of continuous DPPs[END_REF]. Still within statistics, but closer to our motto of using DPPs as subsampling tools, we have investigated DPPs for continuous experimental design [START_REF] Poinas | On proportional volume sampling for experimental design in general spaces[END_REF]. Finally, we have also investigated statistical diagnostics of hyperuniformity [START_REF] Hawat | On estimating the structure factor of a point process, with applications to hyperuniformity[END_REF]. Informally, a hyperuniform point process is one for which the variance of Monte Carlo estimators decreases faster than the classical Monte Carlo rate. Being able to infer that property can help us focus our mathematical efforts on promising point processes.

On experimental and software aspects, we developed a software library1 intended to collect all known sampling algorithms for DPPs, exact and approximate (Gautier, Bardenet, Polito, and Valko, 2019a). We did the same for estimators of the structure factor of a point process, 2 a key quantity to test hyperuniformity, which is especially useful to material scientists. Finally, we are developing a benchmarking suite for signal detection and reconstruction with time-frequency tools. 3 One of our objectives there is to numeri-cally assess in what regime algorithms based on the zeros of the spectrogram shine, and how to combine zero-and maxima-based algorithms.

Other work and a few perspectives

The Monte Carlo algorithms implicitly described in Chapter 3 are, in a sense, importance sampling algorithms with non-i.i.d. proposals. The resulting mean square errors, for volume sampling, or asymptotic variance, in the case of orthogonal polynomial ensembles, decreases fast with the size of the sample, but actually scales exponentially in the dimension of the space over which we wish to integrate. It is thus highly desirable to investigate whether a dynamic version of our results, involving e.g. the trajectory of a Markov process leaving a repulsive point process invariant, could lead to better scaling with the dimension. Indeed, having asymptotic variances that grow polynomially with the dimension is one of the theoretical reasons for the success of MCMC algorithms [START_REF] Chernozhukov | An MCMC approach to classical estimation[END_REF]. With that in mind, Mylène Maida (Laboratoire de mathématiques Paul Painlevé, Univ. Lille) and I are investigating piecewise deterministic Markov processes that leave repulsive point processes invariant. Relatedly, with our new jointly-supervised doctoral student Martin Rouault (2022-), Mylène and I are investigating other types of Gibbs point processes than DPPs for Monte Carlo integration. These other point processes are partly chosen for their amenability to MCMC sampling.

With Subhro Ghosh (Dept. of mathematics and dept. of data science, NUS Singapore) and aspiring PhD student Hugo Simon, we are investigating variance reduction for discrete DPPs in machine learning, in the vein of our work on stochastic gradient [START_REF] Bardenet | Determinantal point processes based on orthogonal polynomials for sampling minibatches in SGD[END_REF]. With Subhro again, we are also investigating methodological alternatives to sampling DPPs [START_REF] Bardenet | Learning from DPPs via sampling: Beyond HKPV and symmetry[END_REF], and topological data analysis applied to the zeros of time-frequency transforms.

With PhD student Diala Hawat (2020-) and her co-supervisor Raphaël Lachièze-Rey (Laboratoire de mathématiques appliquées MAP5, Univ. Paris), we are currently investigating how to obtain variance reduction with repulsiveness applied a posteriori. Indeed, to bypass the cost of sampling a DPP with large cardinality, one option is to sample points independently, and then push them apart. Diala is currently proving (and is almost there) that a well-chosen gradient step applied to all particles, pushing them away from each other, reduces the variance of Monte Carlo estimators based on these particles.

After reading Odile Macchi's thesis, where she introduced DPPs as models to explain the anti-bunching properties of beams of fermionic particles, I grew very enthusiastic about the connections between DPPs and quantum electronic optics. We organized a two-day event 4 in 2019 in Lille, France, opened by Odile Macchi and featuring both theoretical and experimental physicists exposing their view of fermionic coherence to an audience of DPP users across mathematics, computer science, and signal processing. The workshop having been met with cross-disciplinary enthusiasm, we made two decisions. The first was to organize an ambitious, two-week follow-up to the workshop,5 which took place in 2022 in Lyon, France. The second decision was to unite forces with physicists and write a survey on the links between point processes and optics, to pin down a common ground for discussions. With my postdoc Alexandre Feller, I led a team of dedicated speakers and organizers of the aforementioned workshops in writing -actually, Part I of-a survey on DPPs and fermions [START_REF] Bardenet | From point processes to quantum optics and back[END_REF], fifty years after Macchi's 1972 thesis. Part II of the survey is work in progress, and will feature chosen topics, like fermionic techniques in combinatorics, experimental measurements of (anti-)bunching particles, or quantum electronic signal processing. This work has already generated many puzzling questions that could easily become research lines at the intersection of physics, probability, and time-frequency signal processing, see e.g. (Bardenet et al., 2022, Section 6). For instance, physical fermions are described in quantum field theory by both a quantum state and a measurement. From the same state, a lot of point processes can be obtained by varying the measurement process: how do all these point processes relate to each other? How is the amount of repulsiveness in a point process related to the corresponding quantum measurement and state? What are the consequences of observables not commuting, and of Heisenberg's uncertainty bound, on the resulting point processes?

The last research line that I want to mention here has me take a step back. My initial motivation to study repulsive point processes is Monte Carlo integration. In turn, my motivation for Monte Carlo integration is Bayesian statistics, where every question asked by a decision-maker, be they biologists, physicists, or government, is turned into the computation of integrals. But after years of teaching Bayesian learning and reading about the foundations of Bayesian statistics, I have finally grown suspicious of part of the Bayesian thought process. For instance, there is no fully convincing case, in my opinion, for conditioning (i.e., using conditional expectations) to represent the acquisition of knowledge with new observations. In particular in machine learning, when the focus is on prediction accuracy at least as much as on internal coherence, and when models are necessarily misspecified (i.e., incomplete descriptions of the data generating process), I believe that blindly conditioning a large model as one acquires more data is hardly justified, and is not what practitioners do, in any case. I would like to spend time thinking of an axiomatic derivation of how to simultaneously update knowledge and make decisions online, which trades off internal coherence and prediction accuracy. No doubt that my upcoming semester at Univ. Pompeu Fabra Barcelona, amid online learning experts, will fuel my thoughts on this topic. 2015-2019 At École centrale de Lille, France, I taught advanced machine learning topics, with applications to bankruptcy prediction, to master-level engineering students specializing in data science (by then called the "DAD curriculum"). This was a joint course with Éric Séverin, a professor in management at Univ. Lille. I organized this course around a machine learning competition, see e.g. the 2019 edition. 2 The topics of the lectures were partly decided together with students as they progressed through the competition.

2013-2014 At Univ. Oxford, UK, I taught Advanced simulation on Monte Carlo methods to 4th year statistics students, jointly with Pierre Jacob (by then Univ. Oxford, now ESSEC Paris, France).

A.2 Supervision

Since my PhD, I have been involved in the supervision of

• 5 theses (3 completed, 1 discontinued, 1 running). Two more PhD theses are starting these days, as I am writing this document, and do not appear here.

• 6 postdocs (3 completed, 3 running). One more is also being signed these days, and does not appear here.

• 1 research engineer.

This amount of supervision is maybe unusual in France before an HDR, but it is a consequence of me having access to substantial funding, and the doctoral school in Lille consequently giving me pre-HDR license to supervise PhD students.

A.2.1 Postdocs

Figure 1 . 1 :

 11 Figure 1.1: The intersection with [0, 1] of independent samples of (a) a homogenous Poisson point process, and (b) a determinantal point process, more precisely a scaled version of the circular unitary ensemble to approximate the sine process; see Chapter 2 for details.

Figure 2

 2 Figure 2.1: (a, b, c) The intersection with W = [0, 1] of four independent samples of each of the following types of point processes: Poisson, permanental, and determinantal. All three point processes are translation-invariant and are scaled to have the same expected number of points falling in W . (d, e, f) The pair correlation function (2.3) of the three point processes, and histograms of pairwise distances from 100 independent samples.

Figure 2

 2 Figure 2.2: (a) An i.i.d. sample of the uniform distribution on the circle, compared to (b) a sample of the eigenvalues of a Haar-distributed unitary matrix. The arguments of the complex variables in (b) form a draw of the CUE.

  Figure 2.3: White dots are sample zero sets of (a) the planar GAF, (b) the spherical GAF, and (c) the hyperbolic GAF.

"

  Terry Pratchett and Neil Gaiman, Good Omens, 1990 Contents 3.1 A central limit theorem from random matrix theory . . . . . . . 20 3.2 From random matrices to numerical integration . . . . . . . . . . 22 3.2.1 Multivariate orthogonal polynomial ensembles . . . . . . . . . . . . . 22 3.2.2 A stochastic Gaussian quadrature . . . . . . . . . . . . . . . . . . . 22 3.2.3 Our generalization of Johansson's CLT . . . . . . . . . . . . . . . . . 23 3.2.4 An importance sampling version . . . . . . . . . . . . . . . . . . . . 24 3.3 Soviets already did it . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.4 Interpolation and quadrature in RKHSs . . . . . . . . . . . . . . 26 3.4.1 Reproducing kernel Hilbert spaces . . . . . . . . . . . . . . . . . . . 27 3.4.2 Quadrature in a RKHS with a projection DPP . . . . . . . . . . . . 28 3.4.3 Quadrature in an RKHS with continuous volume sampling . . . . . 29 3.5 Closing remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

Figure 3

 3 Figure 3.1: (a) I.i.d. sample of a product of beta marginals, (b) Product Gaussian quadrature, (c) bivariate OPE sample.

  and Cov[y j , y k ] = 0. (3.10)

  .22) Hadamard's inequality and the finite trace condition (3.13) guarantee that the normalization constant Z is finite, so that (3.22) is a bona fide probability distribution. By the Cauchy-Binet identity, it turns out that (3.22) actually defines a statistical mixture of projection DPPs. Proposition 1 Recall the Mercer decomposition k(x, y) = n 1 σ m e m (x)e m (y). For U ⊂ N * , let us further define the projection kernel K U (x, y) = u∈U e u (x)e u (y).(3.23)

Fr

  Figure 4.1: Two multinomial distributions that have been investigated for column subset selection: the squared norms of the columns of X, here denoted as length square distribution, and the k-leverage scores. Figure taken from the PhD thesis of Belhadji (2020).

  the STFT to the planar GAF . . . . . . . . . . . . . . . . . 45 5.2 From a DPP to analytic wavelets . . . . . . . . . . . . . . . . . . . 47 5.2.1 The space of signals . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 5.2.2 The time-scale transform . . . . . . . . . . . . . . . . . . . . . . . . 49 5.2.3 The white noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.3 Mapping signals to GAFs: a generic template . . . . . . . . . . . 50 5.4 From the spherical GAF to the Kravchuk transform . . . . . . . 52 5.5 Opening remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.5.1 Removing the analyticity constraint . . . . . . . . . . . . . . . . . . 54 5.5.2 Making the Kravchuk transform a practical tool . . . . . . . . . . . 54 5.5.3 The point process of maxima . . . . . . . . . . . . . . . . . . . . . . 55 5.5.4 Ridges and level sets . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.6 Despite appearances, this chapter fits in my narrative . . . . . . 55

Figure 5 . 1 :

 51 Figure 5.1: The different steps of a zero-based reconstruction algorithm, applied to a noisy signal. The darker the blue, the larger the spectrogram value, while zeros are shown as white dots.

  Figure 5.2: White dots are the zeros of the spectrogram of (a) a realization of white noise, and (b) a realization of white noise plus an Hermite function.

D

  ξ can be unambiguously defined as a random analytic function in the unit open disk D. Furthermore, it is equal in distribution to GAF (2β) D . In particular, the case β = 0 yields zeros that are a DPP! A sample of this DPP obtained as the numerical zeros of a Paul-Daubechies scalogram is shown in Figure 5.3.

  Figure 5.3: White dots are a sample zero set of the hyperbolic GAF, superimposed on an analytic scalogram of an analytic white noise.

  Figure 5.4: White dots are a sample zero set of the spherical GAF, superimposed on a Kravchuk spectrogram, see (Bardenet and Hardy, 2019; Pascal and Bardenet, 2022) for definitions.

  2021-I am supervising Michaël Fanuel's postdoc on Repulsive point processes for inverse problems.2021-I am supervising Alexandre Feller's postdoc on Determinantal point processes and fermions in quantum optics.2021-I am co-supervising Juan Manuel Miramont's postdoc with François Auger (Univ. Nantes), Pierre Chainais (Univ. Lille), and Patrick Flandrin (ENS Lyon), on Zerobased detection and reconstruction of signals.2021-2022 I supervised Xiaoyi Mai's postdoc on Repulsive point processes in supervised learning. Xiaoyi then got a permanent associate professorship at the department of mathematics (IMT) of Univ. Toulouse, France.

  Correlation functions . . . . . . . . . . . . . . . . . . . . . . . . . .

	Terry Pratchett, Night Watch, 2002
	Contents
	2.1

6 2.2 Determinantal point processes (DPPs) . . . . . . . . . . . . . . . 7 2.2.1 Projection DPPs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2.2 A geometric chain rule for projection DPPs . . . . . . . . . . . . . . 11 2.2.3 DPPs with Hermitian kernels . . . . . . . . . . . . . . . . . . . . . . 12 2.2.4 Sampling DPPs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2.5 DPPs are (usually) repulsive . . . . . . . . . . . . . . . . . . . . . . 15 2.3 Zeros of Gaussian analytic functions . . . . . . . . . . . . . . . . . 16 2.3.1 A constructive definition of GAFs . . . . . . . . . . . . . . . . . . . 16 2.3.2 Three canonical GAFs . . . . . . . . . . . . . . . . . . . . . . . . . . 17

  Johansson's result (3.3) does this for d = 1 and µ the uniform measure on [0, 2π], using nodes (x i ) drawn from a DPP with kernel the projection kernel onto the first N Fourier exponentials, which are orthonormal w.r.t. the uniform measure on S 1 . Betting that this structure is what yields (3.3), and focussing now on a generic measure µ(dx) = ω(x)dx supported in the cube [-1, 1] d , we look for a projection kernel onto N functions that are orthonormal w.r.t. the target measure µ in (3.4).

https://www.legifrance.gouv.fr/loda/id/JORFTEXT000000298904/

Now Paris-Saclay.

And if you don't, I hope the Pratchett quote opening this chapter tips the balance in my favour.

Actually approximate samples of the sine process; see Example 2

For a hierarchical list of publications, a more detailed track record, etc. please refer to my CV at http://rbardenet.github.io/pdf/cvFull_bardenet.pdf.

Without this assumption, a Poisson point process would not necessarily be a simple point process.

While we simply use g < 1 as a token of repulsiveness in this manuscript, much more can be said on quantifying repulsiveness in point processes[START_REF] Biscio | Quantifying repulsiveness of determinantal point processes[END_REF][START_REF] Møller | Couplings for determinantal point processes and their reduced Palm distributions with a view to quantifying repulsiveness[END_REF].

The date mismatch with the reference[START_REF] Macchi | Point processes and coincidences -Contributions to the theory, with applications to statistical optics and optical communication, augmented with a scholion by[END_REF] comes from the fact that the latter is a translated reprint.

This time, to signal processers!

https://github.com/guilgautier/DPPy

https://github.com/For-a-few-DPPs-more/structure-factor

https://github.com/jmiramont/benchmark-test

https://dpp-fermions.sciencesconf.org/

https://indico.in2p3.fr/event/25182/

http://rbardenet.github.com/pdf/cvFull_bardenet.pdf

https://www.kaggle.com/c/dad-bankruptcy-prediction-challenge-2019/

Acknowledgments

" 'Rain don't fall on a witch if she doesn't want it to, although personally

Chapter 6

Other work and perspectives " Nanny's philosophy of life was to do what seemed like a good idea at the time, and do it as hard as possible. It had never let her down. I start by reviewing work on point processes that I chose not to include in the main body of this manuscript for conciseness, but that still relates to the overall investigation of the computational and statistical properties of repulsive point processes. Then I discuss current and future work, with a few perspectives that progressively go beyond my focus on point processes.

Past work on repulsive point processes

I have voluntarily focused on ideas and mathematical results in the main body of this manuscript, to sketch the common umbrella above my work. That being said, I have also worked on methodological, algorithmic, statistical, and experimental aspects of repulsive point processes.

On methodology, recent work that could easily have made it to the main body includes [START_REF] Bardenet | Determinantal point processes based on orthogonal polynomials for sampling minibatches in SGD[END_REF]. We show there how to use the variance reduction that we obtained for multivariate orthogonal polynomial ensembles in [START_REF] Bardenet | Monte Carlo with Determinantal Point Processes[END_REF], see also Chapter 3, but this time for subsampling a finite-but-large set. Our target application is stochastic gradient learning, where random subsets of a large training set (so-called minibatches) are used to estimate a gradient, at each iteration of a gradient descent applied to an empirical loss function. Since the variance of the gradient estimator controls the mean square error of interest [START_REF] Moulines | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF], variance reduction for gradient estimators has been a rich topic of research. We propose in [START_REF] Bardenet | Determinantal point processes based on orthogonal polynomials for sampling minibatches in SGD[END_REF] a DPP-based estimator of the gradient of the empirical loss, with a mean square error that decreases faster than the inverse of the square root of the cardinality of the minibatch. This paper was selected as a spotlight at NeurIPS'21.

Appendices

Appendix A Teaching, supervision, grant management

" 'If only we had laboratories to produce self-replicating scientists, to ex- plore all the worlds. Ah, but we do! They're called university campuses.' "

Terry Pratchett and Stephen Baxter, The Long War, 2013

In the main body of the manuscript, I have focused on scientific content. For an HDR manuscript, I believe that it is also relevant to mention some of my parallel activities since my thesis, say teaching, supervision, and grant management. For more complete information, a detailed CV with track record, hierarchical publication list, awards, editorial duties, etc. can be found on my webpage. 1

A.1 Teaching

CNRS positions come with no teaching duty. Out of personal inclination, I maintain a small teaching activity at the master level. I find the live interaction with students extremely stimulating, and a great complement to research. On top of this, having to teach a topic is probably the best way to make sure you understand it deeply. More prosaically, teaching at the master-level is also the best way to recruit PhD students. Since the end of my PhD, I have been in charge of the following master-level courses.

2019-At ENS Paris-Saclay, France, I am teaching Bayesian machine learning, jointly with Julyan Arbel (Inria Grenoble, France). The master programme is called MVA ("Mathematics for vision and learning"), and it has been so far the most demanded machine learning master nationwide.