Claire Launay
email: claire.launay@parisdescartes.fr

Determinantal point processes for image processing Determinantal point processes for image processing

Keywords: Determinantal point processes, repulsion, subsampling, image, pixels, patches, stationarity, shot noise, inference

Determinantal point processes (DPPs) are probabilistic models of configurations that favour diversity or repulsion. They have recently gained influence in the machine learning community, mainly because of their ability to elegantly and efficiently subsample large sets of data. In this paper, we consider DPPs from an image processing perspective, meaning that the data we want to subsample are pixels or patches of a given image. To this end, our framework is discrete and finite. First, we adapt their basic definition and properties to DPPs defined on the pixels of an image, that we call determinantal pixel processes (DPixPs). We are mainly interested in the repulsion properties of such a process and we apply DPixPs to texture synthesis using shot noise models. Finally, we study DPPs on the set of patches of an image. Because of their repulsive property, DPPs provide a strong tool to subsample discrete distributions such as that of image patches.

Introduction

Determinantal point processes (DPPs) are models of random sets of points that are currently of great interest in several mathematical fields. They were initially studied in probability, in particular for the modeling of fermions in quantum mechanics [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF] or as a process describing the spectrum of random matrices [START_REF] Soshnikov | Determinantal random point fields[END_REF], and recently gained influence in the machine learning community [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF], mainly due to the necessity to handle high amounts of information [START_REF] Tremblay | Determinantal point processes for coresets[END_REF], [START_REF] Celis | Fair and diverse DPP-based data summarization[END_REF] and often high-dimensional data.

The main feature of DPPs is that they provide a family of models of random configurations that favour diversity or repulsion, in the sense that the probability of observing two points close or similar to each other is lower than in the case of the Poisson process whose points are independent [START_REF] Biscio | Quantifying repulsiveness of determinantal point processes[END_REF]. DPPs are completely determined by their correlation function. Unlike Gibbs' point processes, the moments of the DPPs are all known from their correlation function. That is the reason why these processes are an elegant model to reduce the dimensionality of data or to select representative samples from huge sets of points.

The amount of images and video content available is overwhelming. To be handled, to be processed, it needs to be sorted and summarized. That is the purpose of recommendation systems. Some methods using DPPs have been developed to cope with this issue and to enforce diverse subsets, for images selection [START_REF] Kulesza | k-DPPs: Fixed-size determinantal point processes[END_REF], [START_REF] Celis | Fair and diverse DPP-based data summarization[END_REF] or video recommendation [START_REF] Wilhelm | Practical diversified recommendations on Youtube with determinantal point processes[END_REF]. Moreover, images and videos are now in very high resolution, but remain intrinsically redundant. The strategies for video summarization intend to extract meaningful and representative frames using sequential DPPs. This is a type of DPP taking into account the temporal dependencies of video frames [START_REF] Gong | Diverse sequential subset selection for supervised video summarization[END_REF], [START_REF] Mahasseni | Unsupervised video summarization with adversarial LSTM networks[END_REF]. Besides, Chen et al. [START_REF] Chen | Dimensionality reduction based on determinantal point process and singular spectrum analysis for hyperspectral images[END_REF] prove that DPPs can be an appropriate tool to reduce the dimensionality of hyperspectral images, to select representative pixels from these images and be able to process such large-scale data.

In this paper, we will consider DPPs from an image processing perspective, meaning that the data we want to subsample are pixels or patches of a given image. To this end, our framework is discrete and finite. In this setting, the correlation function associated to a DPP is a matrix K that will be called its kernel. First, we need to introduce some notation. The initial data set is denoted by Y and the data items are assimilated with their index, meaning that for a set of N elements, we consider Y = {1, . . . , N }. The cardinality of a set A is denoted by |A|. When M is a N × N matrix, we denote by M A×B , for all subsets A, B ⊂ Y, the matrix (M (x, y)) (x,y)∈A×B and we use the short notation M A = M A×A . When focusing on a specific couple of points, for instance x, y ∈ Y, we sometimes identify M (x, y) and M xy for clarity purpose. If A and B are subsets of Y such that |A| = |B|, the determinant det(M A×B ) is called a minor of M and in case B = A, det(M A ) is called a principal minor of M .

Let us define DPPs in the general discrete case. Consider the set Y = {1, . . . , N } and a Hermitian matrix K of size N × N such that 0 K I, meaning that the eigenvalues of K are in [0, 1]. Then a random set X ⊂ Y is called a determinantal point process with kernel K if it is defined by

P(A ⊂ X) = det(K A ), ∀A ⊂ Y. (1)
We will denote X ∼ DPP(K). For a detailed presentation of discrete DPPs, their properties and some applications to machine learning, we recommend the article of Kulesza and Taskar [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF].

Notice that a DPP is simple: two points of the point process can't coincide. The diagonal coefficients of K define the marginal probabilities of any singleton: ∀x ∈ Y, P(x ∈ X) = K(x, x), [START_REF] Agarwal | Notes on using determinantal point processes for clustering with applications to text clustering[END_REF] and the off-diagonal coefficients of K give the similarity between points. The repulsion property becomes clear when observing the marginal probabilities of pairs of points. The more similar two points are, the less likely they belong to the DPP simultaneously:

∀{x, y} ⊂ Y, P({x, y} ⊂ X) = K(x, x)K(y, y) -|K(x, y)| 2 .

(3)

Moreover, let λ 1 , . . . , λ N be the eigenvalues of K. In this paper, we will adapt these definitions and properties to particular DPPs, defined on the pixels of an image. One of the key contributions of this paper is the definition and the study of this new kind of DPP over pixels, that we present in the following three sections and that we call determinantal pixel process (DPixP). The second section presents the setting adapted to the pixel of an image for the use of DPPs. In the third section, we apply DPixPs to texture synthesis through shot noise models and we present how to adapt a DPixP kernel to a given spot function. The fourth section investigates the estimation of DPixP kernels from one or several samples. We also gather answers to the issue of equivalence classes of DPPs. As two different matrices may generate the same DPP, we give conditions for DPP kernels to be equivalent, in particular the DPixP kernels. In these sections, we provide several algorithms to use DPixPs: sampling, inference, adapted shot noise sampling. Finally in the last section, we study DPPs over the set of patches of an image. We detail how to use DPPs in this setting and how to define an appropriate kernel, depending on the purpose of the subsampling.

Determinantal pixel processes: DPPs on pixels

Notations and definitions

In the following sections, we will consider DPPs defined on the pixels of an image. Let us first define any image as a function u : Ω → R d (d = 1 for gray-scale images and d = 3 for color images), where Ω = {0, ..., N 1 -1} × {0, ..., N 2 -1} ⊂ Z 2 is a finite grid representing the image domain. We use a common matrix convention for the notation of the indices, meaning that u(0, 0) is in the top left corner and N 1 and N 2 are respectively the height and the width of the image. The cardinality of Ω, that is the number of pixels in the image, is denoted by N = |Ω| = N 1 N 2 . Note that, if necessary, the pixels of an image are ordered and they are considered column by column. For any image u : Ω → R d , and y ∈ Z 2 , the translation τ y u of u by the vector y is defined using periodic boundary conditions by

∀ x = (x 1 , x 2 ) ∈ Ω, τ y u(x 1 , x 2 ) := u(x 1 -y 1 mod N 1 , x 2 -y 2 mod N 2 ).
In the following, we consider the Fourier domain Ω = -N 1 2 , . . . , N 1 2 -1 × -N 2 2 , . . . , N 2 2 -1 if N 1 and N 2 are even (otherwise, we consider -N i -1 2 , . . . , N i -1

2
if N i is odd), so that the frequency 0 is centered. We define the discrete Fourier transform of a function f : Ω → C by, for all ξ ∈ Ω,

f (ξ) = F(f )(ξ) = x∈Ω f (x)e -2iπ x,ξ , with x, ξ = x 1 ξ 1 N 1 + x 2 ξ 2 N 2 .
This transform is inverted using the inverse discrete Fourier transform:

∀x ∈ Ω, f (x) = F -1 f (x) = 1 N ξ∈ Ω f (ξ)e 2iπ x,ξ .
Besides, the Parseval formula asserts that for any function f : Ω → C,

f 2 2 = x∈Ω |f (x)| 2 = 1 N ξ∈ Ω | f (ξ)| 2 = 1 N f 2 2 .
Following these conventions, note that given a function f defined on Ω, we consider it is extended by periodicity to Z 2 . For any function f defined on Ω, we set f -(x) := f (-x). The convolution of two functions f and g defined on Ω is given by

∀x ∈ Ω, f * g(x) = y∈Ω f (x -y)g(y),
where again the boundary conditions are considered periodic. Then, f * g can be computed in the Fourier domain, since

∀ξ ∈ Ω, f * g(ξ) = f (ξ) g(ξ).
The autocorrelation of a function f is denoted by R f . It is defined for all x ∈ Ω by

R f (x) = f * f -(x) = y∈Ω f (x -y)f (-y).
Let us consider a DPP with kernel K defined on Ω × Ω. In this work, we will focus on the modeling of textures, which are often characterized by the repetition of a pattern, or small objects which may be indistinguishable individually. Their homogeneous aspect can be naturally modeled by a stationary random field. Thus we will suppose that the point processes under study are stationary and periodic. This hypothesis amounts to consider that the correlation between two pixels x and y only depends on the difference x -y: the distribution is invariant by translation, while assuming periodic boundary conditions. The pixels are ordered column by column so that the ordered index of a pixel x = (x 1 , x 2 ) ∈ Ω is x 1 + 1 + x 2 N 1 . Thus the kernel matrix K is a block-circulant matrix with circulant blocks, entirely characterized by its first column.

Definition 2 A block-circulant matrix with circulant blocks K verifies for all x = (x 1 , x 2 ), y = (y 1 , y 2 ) ∈ Ω, for all τ = (τ 1 , τ 2 ) ∈ Ω,

K(x + τ, y + τ ) = K(x, y),
where we still consider periodic boundary conditions.

Let us define a correlation function C : Ω → C, extended to Z 2 by periodicity, such that K(x, y) = C(x -y), ∀ x, y ∈ Ω. As it entirely characterizes K, it also characterizes the associated DPP. As circulant matrices, block-circulant matrices with circulant blocks are diagonalized in the Fourier basis [START_REF] Philip | Circulant matrices[END_REF][START_REF] Combescure | Block-circulant matrices with circulant blocks, Weil sums, and mutually unbiased bases. II. The prime power case[END_REF], and the eigenvalues of K are the Fourier coefficients of C. See Appendix A for an illustration and some details on the diagonalization of DPixP kernels.

In this new framework, we can define DPPs from their correlation function C, they are now called determinantal pixel processes (DPixP). A DPixP kernel has two representations: C defined on Ω or the initial matrix K defined on Ω × Ω which corresponds to the blockcirculant matrix with circulant blocks whose first column is C. Definition 3 (Stationary DPixP) Let C : Ω → C be a function defined on Ω, extended by periodicity to Z 2 , such that

∀ξ ∈ Ω, C(ξ) is real and 0 ≤ C(ξ) ≤ 1. ( 4
)
Such a function is called an admissible kernel. Any random subset X ⊂ Ω is called a (stationary) DPixP with kernel C and denoted

X ∼ DPixP(C) if ∀A ⊂ Ω, P(A ⊂ X) = det(K A ),
where

K A = (C(x -y)) x,y∈A is a |A| × |A| matrix.

Properties

The next proposition is directly deduced from properties of general DPPs (Proposition 1).

Proposition 4 (Distribution of the cardinality) The cardinality |X| of a DPixP is distributed as the sum

ξ∈ Ω

Ber( C(ξ))

, where for all ξ ∈ Ω, Ber( C(ξ)) are independent Bernoulli random variables with parameters C(ξ). In particular,

E(|X|) = ξ∈ Ω C(ξ) = N C(0) and Var(|X|) = ξ∈ Ω C(ξ)(1 -C(ξ)).
One can notice that it is easy to know and control the expected number of points in the point process. In the following, when comparing different DPixP kernels, we will consider a fixed expected cardinality n, meaning that we will fix C(0) = n N .

Proposition 5 (Separable kernel) Let C 1 and C 2 be two discrete kernels, of dimension 1, defined respectively on {0, ..., N 1 -1} and {0, ..., N 2 -1}, both verifying Equation (4) (for the 1D discrete Fourier transform). Then the kernel C given by ∀x

= (x 1 , x 2 ) ∈ Ω, C(x) = C(x 1 )C(x 2
) is an admissible DPixP kernel. Such a kernel C will be called separable. Notice that in the general discrete case, the first example corresponds to the case where K is diagonal and the second one corresponds to the case where the eigenvalues of K are either equal to 0 or to 1. It is also called a projection DPP and the cardinality of the point process is equal to the number of non-zero eigenvalues, i.e. the rank of K.

Proof Notice that for all ξ = (ξ 1 , ξ 2 ) ∈ Ω, C(ξ) = N 1 -1 x 1 =0 N 2 -1 x 2 =0 C 1 (x 1 )C 2 (x 2 )e -2iπ x 1 ξ 1 N 1 + x 2 ξ 2 N 2 = C 1 (ξ 1 ) C 2 (ξ 2 ). Thus,
Figure 1 presents two samples of these particular cases. Clearly, the projection DPixP enables a more "regular" distribution of the points in the square, tends to avoid both holes and clusters.

The common algorithm to sample exactly general determinantal processes is the spectral algorithm [START_REF] Hough | Determinantal processes and independence[END_REF]. This is a two steps strategy which relies on an eigendecomposition {(λ x , v x )} 1≤x≤N of the matrix K. Indeed, define (B x ) 1≤x≤N , N independent random variables such that

B x ∼ Ber(λ x ) and K B = x∈Ω B x v *
x v x , where v * x denotes the conjugate transpose of the vector v x . Such a matrix K B is a random version of K and Hough and al. [START_REF] Hough | Zeros of Gaussian Analytic Functions and Determinantal Point Processes[END_REF] proved that DPP(K) = DPP(K B ). Hence, the spectral algorithm consists in first drawing N independent Bernoulli random variables of parameters λ x : these variables select n active eigenvalues and eigenvectors, where n is distributed as 1≤x≤N B x . Then, it samples the n points from a projection DPP, obtained from the selected eigenvectors, thanks to a Gram-Schmidt procedure.

Recall that in our discrete stationary periodic framework, the eigenvalues of the matrix K are the Fourier coefficients of C and its eigenvectors are the elements of the Fourier basis (Appendix A). Then an eigendecomposition of a DPixP of kernel C is computed using the 2D Fast Fourier Transform (FFT2) algorithm. Algorithm 1 presents the classic spectral algorithm [START_REF] Hough | Determinantal processes and independence[END_REF][START_REF] Lavancier | Determinantal point process models and statistical inference[END_REF], adapted to sample a DPixP. In this algorithm, (ϕ ξ ) ξ∈ Ω denotes the columns of the unitary discrete Fourier transform matrix, the eigenvectors of K:

∀ ξ ∈ Ω, ∀ x ∈ Ω, ϕ ξ (x) = 1 √ N e -2iπ x,ξ . (5) Algorithm 1 Spectral simulation of X ∼ DPixP(C) • Sample a random field U = (U ξ ) ξ∈ Ω where the U ξ are i.i.d. uniform on [0, 1].
• Define the "active frequencies" {ξ 1 , . . . , ξ n } = {ξ ∈ Ω; U (ξ) ≤ C(ξ)}, and denote

∀x ∈ Ω, v(x) = (ϕ ξ 1 (x), . . . , ϕ ξn (x)) ∈ C n .
• Sample X 1 uniformly on Ω, and define e

1 = v(X 1 )/ v(X 1 ) = N n v(X 1 ). • For k = 2 to n do:
-Sample X k from the probability density p k on Ω, defined by

∀x ∈ Ω, p k (x) = 1 n -k + 1   n N - k-1 j=1 |e * j v(x)| 2   -Define e k = w k / w k where w k = v(X k ) - k-1 j=1 e * j v(X k )e j .
• Return X = (X 1 , . . . , X n ).

Note that the first point is supposed to be sampled from the distribution with density

p(x) = v(x) 2 n , ∀x ∈ Ω. Yet, in this framework, for all x ∈ Ω, v(x) 2 = n
N , thus the first point of the realization is chosen uniformly. This is in accordance with the stationarity of DPixP(C). Because of the eigendecomposition of a matrix of size |Ω| × |Ω| the initial spectral algorithm runs in O(|Ω| 3 ), yet thanks to the use of the Fast Fourier Transform algorithm to compute the Fourier coefficients of C, sampling DPixPs costs O(|Ω| log |Ω|). Whereas in general the spectral algorithm is heavy when dealing with a huge data set, in this setting, it is very efficient. This allows us to handle large images. Thus, in addition to the explicit computation of marginals and of moments of a DPixP from its kernel, this exact sampler is one more asset of this family of point processes with respect to Gibbs processes, which are another common repulsive model.

In spatial statistics, the pair-correlation function (p.c.f.) g X associated to a point process X is used to describe interactions between pairs of points. It characterizes the local repulsiveness of X [START_REF] Biscio | Quantifying repulsiveness of determinantal point processes[END_REF]. For any discrete stationary point process on Ω, it is defined for all x ∈ Ω by

g X (x) = P ({0, x} ⊂ X) ρ 2 ,
where ρ is the intensity of the point process, ρ = E(|X|)

|Ω|

= P(0 ∈ X). It quantifies the degree of interaction between two points separated by a gap x: the closest g is to 1, the less correlated they are. If g(x) > 1, the points are considered to attract each other, whereas if g(x) < 1 the points are considered to repel each other. Notice that if X ∼ DPixP(C),

g X (x) = C(0) 2 -|C(x)| 2 C(0) 2 = 1 - |C(x)| 2 |C(0)| 2 .
Thus, if X is a Bernoulli point process, for all x = 0, g X (x) = 1: there is no interaction between the points. Note also that for any DPixP, g X ≤ 1. During the sequential step of the sampling, each time a pixel is selected, a "repulsion zone" appears around it, where the probability for a pixel to be selected is low and whose shape depends on the kernel function C. This local "repulsion zone" is clearly retrieved in the pair correlation function computation.

Lavancier et al. [START_REF] Lavancier | Determinantal point process models and statistical inference[END_REF] have also studied determinantal point processes in a stationary framework. Their work is based on a continuous setting, with kernels defined on R d , and due to the stationarity, it is also strongly related to the Fourier transform. They study diverse statistics to compare several stationary isotropic kernels. Similarly to our study in the next subsections, they are interested in the quantification of the repulsiveness of different kernels. In particular, they obtain that the most repulsive DPP kernels are projection kernels.

Hardcore repulsion

Gibbs processes are often used as their definition enables to precisely characterize the repulsion. Besides, they can allow us for hard core repulsion, meaning that the points are prohibited from being closer than a certain distance. To compare with this family of point processes, we investigate the possibility of hard core repulsion in the case of DPixPs. First, we study a hard core repulsion for pairs of points. Specifically, if x ∈ Ω and e ∈ Ω (for instance e = (1, 0) or (0, 1)), is there a DPixP kernel such that x and x + e can't belong simultaneously to the sample? Proposition 6 Let us consider X ∼ DPixP(C) on Ω and e ∈ Ω. Then the following propositions are equivalent:

1. For all x ∈ Ω, the probability that x and x + e belong simultaneously to X is zero.

2. For all x ∈ Ω, the probability that x and x + λe belong simultaneously to X is zero, for λ ∈ Q such that λe ∈ Ω.

3. There exists θ ∈ R such that the only frequencies ξ ∈ Ω such that C(ξ) is non-zero are located on the discrete line defined by e, ξ = θ.

X contains almost surely at most one point on every discrete line of direction e.

This is called directional repulsion. Proof Let X be a DPixP defined on Ω with kernel C. First, let us prove that 1 ⇔ 3.

Recall that for all x ∈ Ω, P {x, x + e} ⊂ X = C(0) 2 -|C(e)| 2 . We deduce from the triangle inequality that

|C(e)| = 1 |Ω| ξ∈ Ω C(ξ)e 2iπ e,ξ ≤ 1 |Ω| ξ∈ Ω C(ξ) = C(0),
and the equality holds if and only if all non-zero elements of the left-hand side sum have equal argument. Thus, P {x, x + e} ⊂ X = 0 if and only if there exists θ ∈ R such that for all ξ ∈ Ω, either C(ξ) = 0, or e, ξ = θ. Hence, for all x ∈ Ω, the probability that x and

x + e belong simultaneously to X is zero if and only if the only non-zeros Fourier coefficients of C are aligned in the orthogonal direction of e. Second, let us prove that 2 ⇔ 3. Consider λ ∈ Q such that λe ∈ Ω. Similarly, P {x, x + λe} ⊂ X = 0 if and only if there exists θ ∈ R such that for all ξ ∈ Ω, either C(ξ) = 0, or λe, ξ = θ, meaning that e, ξ = θ λ , which also is the equation of a discrete line orthogonal to e. Finally, suppose that X contains almost surely at most one point on every discrete line of direction e. Then, for all x ∈ Ω, the probability that x and x + e belong to X is zero so 4 ⇒ 1 ⇔ 3. Now assume that the only non-zero Fourier coefficients of C are aligned on a discrete line that is orthogonal to e. As 2 ⇔ 3 for all λ ∈ Q such that λe ∈ Ω, P {x, x + λe} ⊂ X = 0. Hence, X contains at most one point on any line of direction e, which can be described as a hard core repulsion of direction e. The third figure presents a capture of the conditional density while the simulation is in progress, after 15 pixels already sampled. In each pixel, the probability that it is the next point selected is represented by the gray scale: the lighter a pixel is, the greater its probability of being the next point sampled. One can see that as soon as a pixel x is sampled, all the pixels belonging to the horizontal line passing through x have a zero probability of being sampled next. Note that this proposition is not limited to the horizontal or vertical cases.

Proposition 7 Let X ∼ DPixP(C) verifying the properties of Proposition 6, with e = (1, 0), meaning that X contains at most one point on any horizontal line and all non-zero Fourier coefficients of C are aligned on a vertical line. Then C is separable in the sense of Proposition 5. Besides, the associated vertical point process is a DPixP of dimension 1 and conditionally to the drawn ordinates, the associated horizontal point process consists of a single point chosen uniformly and independently from the other horizontal point processes. The same proposition holds for e = (0, 1) and vertical hardcore repulsion (inverting the terms horizontal and vertical).

Proof Consider an admissible DPixP kernel C such that all its Fourier coefficients are either zero either aligned on a vertical line, positioned in c ∈ -N 1 2 , . . . , N 1 2 -1 (here we assume that N 1 is even, the proof is similar if N 1 is odd). Thus we can define two functions

C 1 = 1 c and C 2 = C(c, .) such that for all ξ = (ξ 1 , ξ 2 ) ∈ Ω, C(ξ) = C 1 (ξ 1 ) C 2 (ξ 2 ) = C 2 (ξ 2 )1 c (ξ 1 ). Notice that C = F -1 ( C 1 )F -1 ( C 2 ) = C 1 C 2 .
Such a function C 1 corresponds to an admissible DPixP projection kernel defined in one dimension, drawing one point and remember that the first point of a DPixP is drawn uniformly. Furthermore, C is a separable kernel.

Note that as soon as a pair of points configuration is prohibited, the whole direction is prohibited. As imposing a minimum distance between points is equivalent to prohibiting pair of points configurations in all directions, we deduce that the only DPixP imposing a minimum distance between the points is the degenerate DPixP, consisting of a single pixel. Hence, there is no DPixP with hard core repulsion in the broad sense, as it exists for Gibbs processes.

Shot noise models based on DPPs

Shot noise models and micro-textures

In the following section, we study discrete shot noise models driven by a DPixP. Shot noise models naturally appear to model phenomena such as the superposition of impulses occurring at independent and random times or positions. These models have been introduced in the computer graphics domain with the work of van Wijk [START_REF] Van Wijk | Spot noise texture synthesis for data visualization[END_REF]. Notice that van Wijk uses the expression spot noise texture as the spatial counterpart of 1D shot noise models yet the term shot noise is commonly employed for general models. Thus, in the rest of the section, we use this more general expression. Shot noise models are frequently used to approximate Gaussian textures as they are well-defined and simple mathematical models that allows us for fast synthesis [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF], [START_REF] Galerne | Gabor noise by example[END_REF], [START_REF] Galerne | Texton noise[END_REF]. Here, we are interested in the discrete version of these models on the finite grid Ω = {0, . . . , N 1 -1} × {0, . . . , N 2 -1} ⊂ Z 2 . Definition 8 (Shot noise models based on a discrete point process) Consider X a discrete point process with intensity ρ and g a (deterministic) function defined on Ω, periodically extended to Z 2 . Then, the shot noise random field S based on the points X and the spot g is defined by

∀x ∈ Ω, S(x) = x i ∈X g(x -x i ).
In general, discrete shot noise models are based on a set of n i.i.d. random variables: it amounts to summing n randomly shifted versions of the spot. These models are particularly interesting for Gaussian texture synthesis as they have a Gaussian limit [START_REF] Galerne | Random phase textures: Theory and synthesis[END_REF]. Indeed, in that case, the shot noise is the sum of n i.i.d. random images so that thanks to the Central Limit Theorem, we obtain a Gaussian limit. We study here shot noise models based on DPixPs. At the end of the section, we prove that there is a similar central limit theorem for shot noise models based on DPixPs that needs a modified framework but that ensures a Gaussian limit.

From now on, we consider an admissible kernel C and we suppose that X is the DPixP of kernel C. We study the interactions between the kernel C and the spot function g. To compute the moments of a shot noise model S based on X and a given spot, we need a moment formula ( [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Process[END_REF], [START_REF] Baccelli | Stochastic Geometry and Wireless Networks, Volume I -Theory[END_REF]), also known as the Campbell or Slivnyak-Mecke formula, adapted to our discrete setting in the following proposition.

Proposition 9 (Moments formula for DPixPs) Let X be a DPixP of kernel C defined on Ω, let us consider k ≥ 1 an integer and f a function defined on Ω k . We have

E   = x 1 ,...,x k ∈X f (x 1 , . . . , x k )   = y 1 ,...,y k ∈Ω f (y 1 , . . . , y k ) det((C(y i -y j ) 1≤i,j≤k ), (6) where
= x 1 ,...,x k ∈X
means that the (x i ) are all different. In particular, for k = 1, we have

E x∈X f (x) = C(0) y∈Ω f (y).
Proof By definition of the DPixP of kernel C, for any y 1 , . . . , y k in Ω, we have

P({y 1 , . . . , y k } ⊂ X) = det((C(y i -y j ) 1≤i,j≤k ).
Therefore, by the Slivnyak-Mecke formula [START_REF] Baccelli | Stochastic Geometry and Wireless Networks, Volume I -Theory[END_REF], as we have

E   = x i 1 ,...,x i k ∈X f (x i 1 , . . . , x i k )   = y 1 ,...,y k ∈Ω f (y 1 , . . . , y k )P({y 1 , . . . , y k } ⊂ X),
we obtain the formula of the proposition.

Since X ∼ DPixP(C) is stationary, S as defined in 8 is also stationary, so that E(S(x) k ) = E(S(0) k ) for all x ∈ Ω and for all k ≥ 1.

Proposition 10 (First and second order moments) Let S be a shot noise model based on X ∼ DPixP(C) and the spot g. Recall that R g is the autocorrelation of g. We have

E(S(0)) = C(0) y∈Ω g(y), and for all x ∈ Ω, Γ S (x) := Cov(S(0), S(x)) = C(0)R g (x) -(R g * |C| 2 )(x).
In particular,

Var(S(0)) = C(0) y∈Ω g(y) 2 -(R g * |C| 2 )(0) and ∀ξ ∈ Ω, Γ S (ξ) = | g(ξ)| 2 (C(0) -|C| 2 (ξ)).
Proof First, let us compute the mean value of such a shot noise model S. Using the periodicity of g,

E(S(0)) = E x∈X g(-x) = y∈Ω g(-y)C(0) = C(0) y∈Ω g(y).
Second, let us compute the covariance function of S for all x ∈ Ω,

Γ S (x) = Cov (S(0), S(x)) = E ((S(0)S(x)) -E (S(0)) 2 = E   x 1 ∈X g(-x 1 ) x 2 ∈X g(x -x 2 )   -E (S(0)) 2 = E   = x 1 ,x 2 ∈X g(-x 1 )g(x -x 2 )   + E   x 1 ∈X g(-x 1 )g(x -x 1 )   -E (S(0)) 2 = y 1 ,y 2 ∈Ω g(-y 1 )g(x -y 2 ) C(0) 2 -|C(y 2 -y 1 )| 2 + y∈Ω g(-y)g(x -y)C(0) -E (S(0)) 2 = C(0)g * g -(x) -(g * g - * |C| 2 )(x)

Extreme cases of variance

We set

N = |Ω| = N 1 N 2 ∈ N, the integer n ≤ N and C n the set of admissible kernels such that C(0) = n N . If X ∼ DPixP(C), with C ∈ C n , notice that E(|X|) = |Ω|C(0) = n.
Given a spot function g, we are looking for admissible kernels C ∈ C n that generate shot noise models S of maximal and minimal variance. Indeed, the value Var(S(0)) quantifies a repulsion "in the sense of g" or the regularity of the shot noise. It is related to the superposition of the spot and thus to a particular spatial arrangement of the points which can be adapted to the spot g. The case of a shot noise S based on a spot function g defined as an indicator function gives one some intuition into this idea. If Var(S(0)) is low, the values taken by S are close to its mean value: there are few zones with no spot and few zones with many superpositions of the spot. This means that the points sampled from DPixP(C) tend to be far from one another, according to the shape of the function g and S appears more homogeneous. The repulsion is maximal. On the contrary, when Var(S(0)) is high, S may take high values, so there can be many points in the same area. In that case, the repulsion is minimal.

Proposition 11 (Extreme cases of variance) Fix g : Ω → R + and an integer n ≤ N . The variance of the shot noise model S is maximal if it is based on the Bernoulli DPixP that belongs to C n , meaning that its kernel C is such that C(0) = n N and for all x = 0, C(x) = 0. The variance of the shot noise model S is minimal when it is based on the projection DPixP of n points, such that the n frequencies {ξ 1 , ..., ξ n } associated with the non-zero Fourier coefficients of its kernel maximize

ξ,ξ ∈{ξ 1 ,...,ξn} | g(ξ -ξ )| 2 . Proof Given a fixed integer n ≤ N , let us consider C ∈ C n that maximizes or minimizes Var(S(0)) = C(0)g * g -(0) -(g * g - * |C| 2 )(0) = n |Ω| 2 ξ | g(ξ)| 2 - 1 |Ω| 2 ξ,ξ | g(ξ -ξ )| 2 C(ξ) C(ξ ).
If we identify the function C to a vector of R N , the question becomes finding C ∈ C n that maximizes or minimizes F : R N → R, where

F ( C) = ξ,ξ | g(ξ -ξ )| 2 C(ξ) C(ξ ).

Maximal variance:

We define a scalar product associated to

g for all v, w ∈ R N , by v, w g = ξ,ξ ∈Ω | g(ξ -ξ )| 2 v ξ w ξ = v T Gw where G is the N × N matrix such that G = ( | g(ξ -ξ )| 2 ) ξ,ξ ∈ Ω .
This scalar product is well defined as it is bilinear, symmetric and for all

v ∈ R N , N ξ,ξ =1 | g(ξ -ξ )| 2 v ξ v ξ = R g * | v| 2 (0) ≥ 0
, where R g is the autocorrelation of g, and

v, v g = 0 ⇔ v = 0.
Notice that since G is symmetric positive definite then F : C → C, C g is strictly convex. The case of maximal variance is achieved for the vector C that minimizes this strictly convex function on the convex set C n : the problem has at most one solution [START_REF] Boyd | Convex Optimization[END_REF].

According to the Cauchy-Schwarz inequality, we have for all

v, w ∈ R N , | v, w g | ≤ v g w g .
Let us pick v = C, the vector whose components are the Fourier coefficients of a kernel C ∈ C n and w = 1 (= (1, 1, . . . , 1) the constant vector of size N ). We have

v 2 g = F ( C) and w 2 g = ξ,ξ | g(ξ -ξ )| 2 = ξ,ξ R g (ξ -ξ ) = N 2 R g (0). Hence v g w g = N 2 F ( C)(g * g -)(0) and | v, w g | = ξ,ξ | g(ξ -ξ )| 2 C(ξ) = ξ C(ξ) ξ | g(ξ -ξ )| 2 = n N R g (0). Thus, F ( C) ≥ n 2 R g (0) and F ( C) is minimal if and only if C is proportional to w: necessarily, for all ξ ∈ Ω, C(ξ) = n N .
Hence, C is a Bernoulli process. This kernel maximizes the variance of any shot noise S, independently of the spot g. It it the least repulsive DPixP.

Minimal variance:

Let us characterize the kernel C that maximizes the function F on the convex set C n . F is quadratic so that solutions are on the boundaries of C n , meaning that for all kernel

C * ∈ C * F := {argmax C (F ( C))}, ξ C * (ξ) = n and ∀ξ ∈ Ω, C * (ξ)(1 -C * (ξ)) = 0.
Thus, the solutions are the projection DPixP kernels C * with exactly n frequencies

{ξ 1 , ..., ξ n } ⊂ Ω such that C * (ξ i ) = 1 chosen so that ξ,ξ ∈{ξ 1 ,...,ξn} | g(ξ -ξ )| 2 is maximal.
In the end, to determine the maximal repulsion kernel, one needs to maximize a quadratic function, which is NP-hard in general. In practice, it amounts to solve a combinatorial problem. It is possible to approximate the solution thanks to a greedy algorithm: first, one chooses two frequencies Figure 3 presents some results of this algorithm. This figure shows that a projection DPixP adapted to g generates shot noise models with very few spot superpositions. Recall that in section 2, we proved that it was impossible to completely prevent superpositions. Yet, it is possible to characterize the least and the most repulsive DPixPs according to a specific desired repulsion. These extreme cases are coherent with the results of Biscio and Lavancier [START_REF] Biscio | Quantifying repulsiveness of determinantal point processes[END_REF] who quantified the repulsion of stationary DPPs defined on R d and stated that the least repulsive DPP is the Poisson point process whereas the most repulsive family of DPP contains the kernels C such that their Fourier transform F(C) is the indicator function of a Borel set, an analog to the projection DPixPs defined here.

ξ 1 , ξ 2 maximizing | g(ξ 1 -ξ 2 )| 2 then, recursively, one chooses the kth frequency ξ k , 2 < k ≤ N, such that it maximizes ξ∈{ξ 1 ,...,ξ k-1 } | g(ξ -ξ k )| 2 .

Convergence to Gaussian processes

Shot noise models driven by DPixP enable more diverse types of textures than the usual shot noise models, based on points drawn uniformly and independently. As with usual shot noise models based on discrete Poisson processes, it is appealing to study the behavior of the model when the density of the point process increases and tends to infinity. The model presented here takes into account the case where the points are sampled independently, when the shot noise is based on a Bernoulli process. Yet, usual discrete shot noise models, as defined in [START_REF] Galerne | Random phase textures: Theory and synthesis[END_REF], are based on a point process that is not simple, meaning that the points can coincide. Hence it is possible to study the asymptotic behavior of the shot noise model when the intensity of the point process grows. On the contrary, in this work, the DPixPs we study are by definition simple: the points can't coincide. Thus, the framework needs to be adapted so that Ω expands to Z 2 . In spatial statistics, this study is called increasing-domain asymptotics [START_REF] Cressie | Statistics for spatial data[END_REF] and is commonly used when the data is defined on a lattice. In the following subsection, we will consider a fixed DPP kernel and spot function expanding with the image domain. A related configuration is to consider Ω as a grid in [0, 1] 2 = T 2 , the torus of dimension 2, that is refined. In that setting, called infill asymptotics [START_REF] Cressie | Statistics for spatial data[END_REF], the kernel of the DPP is rescaled so that the points are allowed to be increasingly close and the number of points inside [0, 1] 2 tends to infinity. In these configurations, it is possible to characterize asymptotic behaviors and to derive limit theorems such as a Law of Large Numbers or a Central Limit Theorem. To this end, let us consider stationary determinantal point processes on Z 2 [START_REF] Shirai | Fermion Process and Fredholm Determinant[END_REF][START_REF] Lyons | Stationary determinantal processes: Phase multiplicity, Bernoullicity, entropy, and domination[END_REF], that we will also call determinantal pixel processes. This point process is defined by a discrete bounded operator K on 2 (Z 2 ). That means that K :

2 (Z 2 ) → 2 (Z 2 ), f → Kf such that ∀ t ∈ Z 2 , Kf (t) = s∈Z 2 K(t, s)f (s). We suppose that this DPP is stationary: we define a kernel function C : Z 2 → C, such that K(t, s) = C(s -t) and C ∈ 2 (Z 2 ). Then for all t ∈ Z 2 , Kf (t) = s∈Z 2 C(s -t)f (s): such a K is a convolution operator.
As C belongs to 2 (Z 2 ), there exists a function

C ∈ L 2 (T 2 ) such that C : T 2 → [0, 1], ∀t ∈ Z 2 , C(t) = T 2 C(x)e 2iπt•x dx and C is equal to x → t∈Z 2 C(t)e -2iπt•x in the sense of L 2 (T 2
). Finally, the point process X ∼ DPixP(C) is defined by ∀A ⊂ Z 2 , a finite subset,

P(A ⊂ X) = det(C A ), where C A = (C(x i -x j )) x i ,x j ∈A .
This new definition of DPixPs on Z 2 is an extension of the point process defined on Ω. The main properties of DPixPs are preserved and it allows us to study the asymptotic behavior of shot noise models driven by DPixPs, when the grid is refined or equivalently when the support of the spot is spread out. To do so, we need to consider spot functions defined on R 2 .

Shirai and Takahashi [START_REF] Shirai | Random point fields associated with certain Fredholm determinants II: Fermion shifts and their ergodic and Gibbs properties[END_REF] state the following limit theorems. Some guidelines for the proofs can be found in [START_REF] Shirai | Random point fields associated with certain Fredholm determinants II: Fermion shifts and their ergodic and Gibbs properties[END_REF] for the Z 2 case and in [START_REF] Shirai | Fermion Process and Fredholm Determinant[END_REF] and [START_REF] Shirai | Random point fields associated with certain Fredholm determinants. I. Fermion, Poisson and boson point processes[END_REF] for its continuous counterpart.

Proposition 12 (Limit theorems for shot noise models and DPixPs [START_REF] Shirai | Random point fields associated with certain Fredholm determinants II: Fermion shifts and their ergodic and Gibbs properties[END_REF]) Let f be a bounded measurable function on R 2 with compact support, and X ∼ DPixP(C) with C some admissible kernel on Z 2 . Then, we have the following Law of Large Numbers

1 M 2 x∈X f x M ----→ M →∞ C(0) R 2 f (x)dx, a.e and in L 1 .
Moreover, assume that f is continuous and

R 2 f (x)dx = 0. Then, ∀t ∈ R, lim M →∞ E exp i √ M 2 x∈X f x M = exp - 1 2 σ(C) 2 f 2 2 where σ(C) 2 = C(0) - x∈Z 2
|C(x)| 2 , and consequently,

1 √ M 2 x∈X f x M D ----→ M →∞ N (0, σ(C) 2 f 2 2 ).
In the following, let g be a spot function, that we assume continuous, with compact support, and M > 0. Denote the M -normalized shot noise S M associated to g defined for

all y ∈ Z 2 by S M (y) = 1 M 2 x∈X g y - x M
. We obtain a Law of Large Numbers for the shot noise driven by DPixPs:

S M (0) = 1 M 2 x∈X g - x M ----→ M →∞ C(0) R 2 g(x)dx, a.e and in L 1 . (7)
Finally, it is also possible to obtain a multidimensional central limit theorem thanks to the previous formulations.

Proposition 13 (Central limit theorem for shot noise models) Let g be a continuous function on R 2 with zero mean and compact support, X ∼ DPixP(C) and the related shot

noise S M : S M (y) = 1 M 2 x∈X g y - x M , ∀y ∈ Z 2 . Then, ∀x 1 , ..., x m ∈ Z 2 , √ M 2 (S M (x 1 ), • • • , S M (x m )) D ----→ M →∞ N (0, Σ(C))
where for all k, l ∈ {1, • • • , m}, Proof Consider the M -normalized shot noise S M associated to g:

Σ(C)(k, l) = C(0) -C 2 2 R 2 g(x k -t)g(x l -t)dt = C(0) -C 2 2 R g (x l -x k ).
∀y ∈ Z 2 , S M (y) = 1 M 2 x∈X g y - x M . By setting ∀u ∈ R m , ∀x 1 , ..., x m ∈ Z 2 , ∀x ∈ R 2 , f (x) = u 1 g(x 1 -x) + u 2 g(x 2 -x) + • • • + u m g(x m -x),
f is continuous on R 2 , with compact support such that R 2 f (x)dx = 0 so it is possible to apply the limit theorem 12 and the Levy's continuity theorem.

Thus, shot noise models driven by a DPixP also converge to a Gaussian limit whose covariance is related to the spot and to the kernel C of the point process. Note that, in the previous proposition, the limit variance Σ(C) is equal to the product of a constant depending on the kernel C and the autocorrelation of the spot g. Similarly, a normalized Poisson shot noise associated to the spot g converges towards the distribution N (0, R g ), where R g is the autocorrelation of g [START_REF] Galerne | Random phase textures: Theory and synthesis[END_REF]. As the Bernoulli case corresponds to the kernel function C = δ 0 , we retrieve the same result here.

Notice the similar work in a continuous framework of Poinas et al. on the limit distribution of sums of functionals of DPPs defined on R d [START_REF] Poinas | Mixing properties and central limit theorem for associated point processes[END_REF]. Figure 4 presents the asymptotic behavior of shot noise models driven by the continuous spot function shown in (a), whose integral is zero, and a DPixP on Z 2 with a kernel whose Fourier coefficients are given by an isotropic 2D Gaussian function. When the grid is refined, the shot noise as defined in this section tends to a Gaussian texture associated to the spot and the kernel of the DPixP.

There are several perspectives to extend these results. Note that there is no more interaction between the spot and the kernel in the limit and the higher the repulsion is, in the sense of high kernel coefficients, the lower the variance is. In particular, for projection kernels C, the asymptotic variance Σ(K) vanishes, as 0). This suggests that the convergence rate must be faster for projection DPixPs. The convergence of determinantal shot noise models towards a Gaussian process has also been proven by Soshnikov [START_REF] Soshnikov | Gaussian limit for determinantal random point fields[END_REF], with some relaxed assumptions. In his work, the spot function g does not need to be continuous or to have a zero mean but it needs to be bounded with precompact support and the expectation and variance of the shot noise are constrained. We plan to adapt these results to our framework to lighten the assumptions on the spot function.

C 2 2 = T 2 | C(x)| 2 dx = T 2 C(x)dx = C(

Identifiability and Inference for Determinantal Pixel processes

One of the purposes of statistical inference is to fit a predetermined model to data that can be represented by points, using information on their global or local behaviour. When the data are assumed independent and well represented by a homogeneous point process, one can use Poisson point processes. Yet, some data may present attraction or repulsion, they may also have an anisotropic structure. DPixP models can be suitable for representing 2-dimensional discrete data points with repulsion. For instance, the positions of plant seeds [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Process[END_REF] or trees in a forest [START_REF] Lavancier | Determinantal point process models and statistical inference[END_REF] often exhibit repulsion because of limited shared supply, but also anisotropy due to environmental factors as wind orientation or ground steepness. In [START_REF] Perrin | Point processes in forestry : an application to tree crown detection[END_REF], Perrin et al. study statistical inference on repulsive point processes, Strauss processes, to detect tree crowns on aerial images of forests. Similarly, Descamps et al. [START_REF] Descamps | Automatic flamingo detection using a multiple birth and death process[END_REF] use Gibbs point processes to automatically detect flamingos on images of colonies. Once one has inferred the parameters of an appropriate model, it is possible to reproduce similar data, to detect objects or distinguish different regions by statistical testing.

Learning the parameters of a determinantal point process, either the whole underlying kernel K as in [START_REF] Kulesza | Learning with Determinantal Point Processes[END_REF][START_REF] Affandi | Learning the parameters of determinantal point process kernels[END_REF] or a few parameters encoding the kernel as in [START_REF] Bardenet | Inference for determinantal point processes without spectral knowledge[END_REF][START_REF] Biscio | Contrast estimation for parametric stationary determinantal point processes[END_REF], is still considered as a difficult task, first because the likelihood is often non-convex, and most of all because it is complex to compute as it uses the determinant of a huge matrix. Most papers studying inference for DPPs overcome this difficult computation by using restrictive hypothesis on the kernel such as in the papers [START_REF] Kulesza | Learning determinantal point processes[END_REF] or [START_REF] Affandi | Learning the parameters of determinantal point process kernels[END_REF]. Bardenet and Titsias [START_REF] Bardenet | Inference for determinantal point processes without spectral knowledge[END_REF] develop bounds on the likelihood and use Markov Chain Monte Carlo methods to infer the parameters of the kernel. On the other hand, using descriptive statistics to fit the models to the data enables to cope with this difficult computation and to obtain more efficient inference algorithm. It is the approach that we choose in this paper. Some authors try to infer first order characteristics such as the intensity of the point process [START_REF] Biscio | Standard and robust intensity parameter estimation for stationary determinantal point processes[END_REF], which provides the average number of points in a given area. In our finite and discrete setting, we can obtain a direct estimation of the intensity, as the ratio between the number of points and the size of the domain. Several second order characteristics are used to describe a sample, for example the empty space distance, the cumulative nearest-neighbor function, the pair correlation function (p.c.f. in short), presented above, or the Ripley's K function, closely related to the p.c.f (see [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Process[END_REF] for a detailed presentation). These statistics provide information on the interactions between points. Møller and Waagepetersen [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Process[END_REF] present these different statistics and state that higher order characteristics may be less stable if the number of points is low. In the following, we choose to focus on a quantity related to the p.c.f. It has several advantages: it is easy to

C 1 C 2 C 3

Equivalence classes of DPP and DPixP kernels

A model is not identifiable if two different parametrizations are equivalent. Here, it would correspond to several different kernel functions generating the same DPixP. Lavancier et al. [START_REF] Lavancier | Adaptive estimating function inference for nonstationary determinantal point processes[END_REF] proved that in a continuous setting, assuming that its intensity function is positive, a DPP kernel is uniquely defined. Yet, in the discrete case, DPPs and DPixP are not identifiable, as illustrates Figure 5. It is crucial, in particular for estimation purposes, to characterize these equivalence classes of kernels. Of course this question is also decisive in more general cases, when the kernel matrix K is Hermitian, with real or complex coefficients. We propose here a brief synthesis of what is known on this question, and we add a study on DPixP kernels.

The distribution of a DPP is entirely defined by all its principal minors (see Equation 1), thus characterizing DPP kernel equivalences classes is equivalent to understanding the consequences of equal principal minors on matrices, in the symmetric or Hermitian cases, and in the DPixP framework where the matrix is Hermitian circulant.

Notice that the characteristic polynomial of a matrix can be written as a function of its principal minors: det(tI

+ K) = N k=0 (-1) k   A⊆Y,|A|=k det K A   t n-k .
Hence, two matrices with equal principal minors have equal characteristic polynomial so they have the same eigenvalues, with the same algebraic multiplicity. Thus, two kernel matrices generating the same DPP have the same spectrum. A key notion here is the diagonal similarity between two matrices: two square matrices M 1 , M 2 are called diagonally similar if there exists a diagonal matrix D such that M 2 = D -1 M 1 D. In the following, we also need the notion of the directed graph associated to a matrix ( [START_REF] Engel | Matrices diagonally similar to a symmetric matrix[END_REF], [START_REF] Hartfiel | On matrices having equal corresponding principal minors[END_REF], [START_REF] Kulesza | Learning with Determinantal Point Processes[END_REF]). Consider a matrix M of size N × N . Its associated directed graph G M contains the N vertices Y = {1, . . . , N } and an edge between the vertices x and y if and only if M (x, y) = 0. The matrix M is called irreducible if G M is strongly connected, meaning that there exists a sequence a path from any vertex to any other one. In the opposite case, the matrix is called reducible, which is equivalent to being permutation-similar to a block upper triangular matrix. Besides, it is called completely reducible if it is permutation-similar to a block diagonal matrix with irreducible blocks, meaning that there exists a permutation matrix P such that

P T M P =   M 1 0 . . . 0 Mr   , M 1 , .
. . , M r irreducible. Notice that a Hermitian matrix is either irreducible or completely reducible.

Let us consider two general admissible DPP kernels K 1 and K 2 , admissible meaning that they are Hermitian and their eigenvalues are in [0, 1]. Thanks to basic determinant properties, notice that if there exists a diagonal matrix D such that

K 2 = D -1 K 1 D or K T 2 = D -1 K 1 D
, then K 1 and K 2 have same principal minors, that is, the equivalence class of a DPP kernel contains all the admissible matrices of which the kernel matrix itself or its transpose is diagonally similar. In the case where the DPP kernel is real, Kulesza [START_REF] Kulesza | Learning with Determinantal Point Processes[END_REF] proved the following proposition. The proof of this proposition is in two parts. First, the author demonstrates the relation when all coefficients of the matrices are non-zero. Then, using graph theory, Kulesza extend this proof to matrices associated to a connected graph and finally to a disconnected graph, when the matrix is reducible. This equivalence property for real DPP kernels has impacted serveral learning strategies as in [START_REF] Rising | An efficient algorithm for the symmetric principal minor assignment problem[END_REF], [START_REF] Brunel | Learning signed determinantal point processes through the principal minor assignment problem[END_REF], [START_REF] Urschel | Learning determinantal point processes with moments and cycles[END_REF] or [START_REF] Brunel | Rates of estimation for determinantal point processes[END_REF] which try to estimate real DPP kernels from several i.i.d. samples. In particular, the first two papers intend to solve the so-called principal minor assignment problem for symmetric matrices, and Brunel et al. [START_REF] Brunel | Rates of estimation for determinantal point processes[END_REF] maximize a log-likelihood depending on the equivalence class of DPP kernels. Urschel et al. [START_REF] Urschel | Learning determinantal point processes with moments and cycles[END_REF] obtain a bound on a distance between the estimated kernels L * and the equivalence class of the original kernel: min D L * -D -1 LD F , on diagonal matrices D with coefficients only equal to 1 or -1.

In the paper [START_REF] Stevens | Equivalent symmetric kernels of determinantal point processes[END_REF], Stevens characterizes equivalence classes of real or complex symmetric DPP kernels. We would like to characterize DPP equivalence classes in a more general setting, where the DPP kernels are no longer real or symmetric but complex and Hermitian. Schneider, Saunders and Engel ([51], [START_REF] Engel | Matrices diagonally similar to a symmetric matrix[END_REF]) worked on the relation between equal principal minors and diagonal similarity through graph theory: see for instance [START_REF] Saunders | Flows on graphs applied to diagonal similarity and diagonal equivalence for matrices[END_REF] for links between equality of cyclic products and diagonal similarity, or [START_REF] Engel | Matrices diagonally similar to a symmetric matrix[END_REF] where they deal with real symmetric matrices. In 1986, Loewy [START_REF] Loewy | Principal and diagonal similarity of matrices[END_REF] gives several sufficient conditions ensuring that if two square matrices have equal principal minors, one is diagonally similar to the other one or to the conjugate of the other one. We adapt these conditions to Hermitian DPP kernels in Theorem 16. In the following, we define D N ⊂ M N (C) as the set of diagonal matrices of size N × N such that its coefficients are of modulus one.

Lemma 15

Let K 1 and K 2 be two irreducible Hermitian matrices and assume that there exists an invertible diagonal matrix D such that

K 2 = D -1 K 1 D or K T 2 = D -1 K 1 D.
Then all the coefficients of D have the same modulus so one can choose D in D N .

Proof Assume that K 1 and K 2 are two irreducible Hermitian matrices and there exists a diagonal matrix D such that

K 2 = D -1 K 1 D or K T 2 = D -1 K 1 D. First, let us suppose that K 2 = D -1 K 1 D.
For all x, y ∈ Y such that K 1 (x, y) = 0, we have also K 2 (x, y) = 0 and

K 2 (x, y) = 1 d x K 1 (x, y)d y .
As K 2 is Hermitian,

K 2 (x, y) = K 2 (y, x) = 1 d y K 1 (y, x)d x = d x d y K 1 (x, y). Then d y d x = d x d y , hence for all x, y ∈ Y such that K 1 (x, y) = 0, |d x | = |d y |. Now recall that K 1 is irreducible.
Its associated graph is connected and every node is reachable from any other node so it is possible to propagate this equality so that for all x, y ∈ Y, |d x | = |d y | = λ. Then without loss of generality, changing if necessary to 1 λ D, we can choose D as the matrix such that

K 2 = D -1 K 1 D with diagonal coefficients of modulus equal to 1. The proof is similar if K T 2 = D -1 K 1 D.
Now we can prove the following theorem on the equivalence classes of Hermitian DPP kernels.

Theorem 16 (Identifiability for Hermitian DPP kernels) Let N be a positive integer and let Y = {1, . . . , N }. Suppose that K 1 , K 2 ∈ M N (C) are two Hermitian admissible DPP kernels and that K 1 is irreducible. If N ≥ 4, suppose furthermore that, for every partition of Y into subsets α, β such that |α| ≥ 2, |β| ≥ 2, rank (K 1 ) α×β ≥ 2. Then, the following propositions are equivalent:

(i) DPP(K 1 ) = DPP(K 2 ), (ii) There exists a diagonal matrix D such that K 2 = D -1 K 1 D or K T 2 = D -1 K 1 D, (iii) There exists a diagonal matrix D ∈ D N such that K 2 = D -1 K 1 D or K T 2 = D -1 K 1 D.
= D -1 K 1 D or K T 2 = D -1 K 1 D.
Notice that these two theorems, making the distinction between rank(K 1 ) α×β and rank(K 1 ) β×α , are equivalent in this Hermitian setting. Then (i) is equivalent to (ii). Besides, clearly (iii) implies (ii) and under these assumptions, by Lemma 15, (ii) implies (iii).

In this general setting, assuming that K 1 is irreducible is crucial. Indeed, Hartfiel and Loewy [START_REF] Hartfiel | On matrices having equal corresponding principal minors[END_REF] provide counterexamples of two admissible kernels generating the same DPP distribution without being diagonally similar.

We now turn to the special case of DPixP defined on Ω, the image domain of size N 1 × N 2 . Their kernel matrices are Hermitian block-circulant with circulant blocks. Recall that matrices generating DPixPs have all the same eigenvectors, the vectors of the Fourier basis. We also know that two matrices generating the same DPixP distribution have the same eigenvalues, so there is at most (N 1 N 2 )! different kernels associated to one DPixP model. In the following proposition and remark, we prove that in most cases, the class of equivalence is much more constrained.

Proposition 17 (Identifiability for DPixP)

Let Ω be a finite grid of size N 1 × N 2 , and C 1 , C 2 be two admissible DPixP kernels on Ω in the sense of Definition 3, generating the block-circulant matrices K 1 and K 2 that satisfy the hypothesis of Theorem 16. Then, DPixP(C 1 ) = DPixP(C 2 ) if and only if there exists a translation mapping the Fourier coefficients of C 2 to the Fourier coefficient of C 1 or to their symmetry with respect to (0, 0), meaning that

DPixP(C 1 ) = DPixP(C 2 ) ⇐⇒ ∃ τ ∈ Ω s.t. either ∀ξ ∈ Ω, C 2 (ξ) = C 1 (ξ -τ ) or ∀ξ ∈ Ω, C 2 (ξ) = C 1 (-ξ -τ ). ( 8
)

Proof

As K 1 and K 2 satisfy the hypothesis of Therorem 16, there exists an invertible diagonal matrix D such that

K 2 = D -1 K 1 D or K T 2 = D -1 K 1 D
, where D ∈ D N , meaning that D is a diagonal matrix with coefficients of modulus equal to one. First, assume that

K 2 = D -1 K 1 D. Define for all x ∈ Ω, θ x ∈ [0, 2π[ such that D(x, x) = e iθx .
The goal is to prove that there exists τ such that θ x = 2π x, τ , for all x ∈ Ω. Notice that, by changing D into 1 D(0,0) D, we can assume that θ 0 = 0, that is D(0, 0) = 1. By assumption, we obtain ∀x, y ∈ Ω, K 2 (x, y) = C 2 (y -x) = e -iθx K 1 (x, y)e iθy = e i(θy-θx) C 1 (y -x),

and C 2 (x) = C 2 (x -0) = e iθx C 1 (x).
Recall, thanks to Equations ( 2) and (3), that C 1 (0) = C 2 (0) and that, for all x ∈ Ω,

|C 1 (x)| = |C 2 (x)|. As C 2 (x) = 0 if and only if C 1 (x) = 0, for such x ∈ Ω, any value θ x is valid. Consider the set Ω * C = {x ∈ Ω; C 1 (x) = 0}
. For all z ∈ Ω, and all x ∈ Ω, we have

C 2 (z) = e iθz C 1 (z) = C 2 (z + x -x) = e i(θ z+x -θx) C 1 (z + x -x) = e i(θ z+x -θx) C 1 (z).
Denote for all x ∈ Ω, α(x) = e iθx . Thus, for all z ∈ Ω * C , for all x ∈ Ω, α(z

) = α(z + x)α(x), meaning that α(x) = α(z + x)α(z). For all ξ ∈ Ω, for all z ∈ Ω * C , we have α(ξ) = x∈Ω α(x)e -2iπ x,ξ = x∈Ω α(z)α(z + x)e -2iπ
x,ξ = α(z)e 2iπ<z,ξ> α(ξ).

As α is not the zero function, consider τ ∈ Ω such that α(τ ) is non-zero. Then, for all z ∈ Ω * C , α(z) = e 2iπ<z,τ > . Thus, for all z ∈ Ω * C , C 2 (z) = e 2iπ<z,τ > C 1 (z), which is also true for z such that C 1 (z) = 0. To conclude, for all z ∈ Ω, C 2 (z) = e 2iπ<z,τ > C 1 (z). In the second case when

K T 2 = D -1 K 1 D,
< 2. If rank(K 1 ) α×β = 0, that is (K 1 ) α×β = 0.
There exists a permutation matrix such that K 1 is permutation similar to a block diagonal matrix, which is in contradiction with the irreducible hypothesis. Hence, rank(K 1 ) α×β = 1. This means that there exist two vectors u ∈ C |α| \ {0} and v ∈ C |β| \ {0} such that (K 1 ) α×β = u T v. In practice, as K 1 is Hermitian and the Fourier coefficients of C are real, the coefficients of the matrix K 1 are tightly constrained. The matrix is determined by a small number of modulus and arguments. Then, when assuming that K 1 and K 2 are equivalent, as DPixP kernels, the matrices are even more constrained. See Appendix B.1 for a simple example of this configuration. Notice that in the 1D case of dimension 5, two equivalent DPixP kernels K 1 and K 2 in this configuration still verify that there exists a diagonal matrix

D ∈ D N such that K 2 = D -1 K 1 D or K T 2 = D -1 K 1 D.
Our conjecture is that this is always the case, whatever the dimension of Ω. Thus, this assumption on the rank of the submatrix (K 1 ) α×β leads to degenerate kernels that are numerically "rare".

Case 3: K 1 is not irreducible. Then, as a Hermitian or circulant matrix, K 1 is necessarily completely reducible, meaning that there exists a permutation matrix P such that K 1 is permutation similar to a block diagonal matrix with irreducible blocks. We prove in Appendix B.2 that these blocks are copies of one Hermitian block-circulant sub-matrix, that we can call the canonical block: they all have equal size and the coefficients are identical. Note that restricting DPP to a subset A define also a DPP on this subset A [34, Section 2.3]. Furthermore, as each block matrix is still circulant, each one defines a sub-DPixP defined on the associated subset of pixels. By assumption, these blocks are irreducible so they are either in the first or in the second configuration. Let us consider K 2 a DPixP kernel equivalent to K 1 . Thanks to the modulus equality, K 2 is similar to a block diagonal matrix with blocks of same size, using the same permutation matrix. If the canonical block is in the first configuration, verifying the rank hypothesis of Theorem 16, the final diagonal matrix D is simply the concatenation and rearrangement of all the diagonal sub-matrices D i associated to its respective i-th block. Notice that as the block submatrices are identical to the canonical block and each one concerns a different set of pixels, all submatrices are in the same configuration, meaning that either for all submatrices

K 1i of K 1 , K 1i = D i K 2i D i or for all submatrices K 1i , K 1i = D i K 2i D i .
On the other side, if the canonical block is in the second configuration, we can't conclude on the similarity of both matrices K 1 and K 2 in the general case yet. Notice that this completely reducible hypothesis is quite degenerate. It corresponds to a DPixP defined on an image domain that can be partitioned in groups of pixels evenly spaced with independence from one group to the other: that means that the pixels are independent to their immediate neighbors. A typical example of this model would be image domain partitioned following a grid. As DPixPs deals with spatial repulsion, there seems to be few applications of such models.

It is important to notice that the size of the equivalence classes we characterized in Proposition 17 is small and known: given a DPixP kernel verifying the appropriate hypothesis, it admits at most 2|Ω| equivalent kernels, generating the same DPixP distribution. Moreover, we have shown previously how a kernel that does not verify the hypothesis of the proposition is quite degenerate: in practice, when dealing with kernels adapted to a given problem, these hypothesis are always verified. In fact, we were not able to find an example of DPixP kernels C 1 and C 2 such that DPixP(C 1 ) = DPixP(C 2 ) and which does not verify the right hand side of Equivalence [START_REF] Biscio | Quantifying repulsiveness of determinantal point processes[END_REF]. We conjecture that Equivalence (8) holds for all DPixP kernels, regardless of the hypotheses of Theorem 16. Characterizing equivalence classes of DPPs and DPixPs is crucial for the estimation of DPixP kernels from point process samples. This is what we investigate in the next subsection.

Learning a DPixP kernel from one realization

First, we address the question of inference from one single realization. Consider one set of points Y on Ω, the finite and discrete grid of size N 1 × N 2 = N and assume that Y has been sampled from a certain DPixP of kernel C 0 . Note that in general, one realization does not provide enough information to characterize a model. Yet, due to the stationarity of the kernels we consider, all the translations of Y can also be seen as samples drawn by the same DPixP kernel C 0 .

Let n = |Y | denote the cardinality of Y . The problem is to find C e an admissible DPixP kernel that estimates C 0 , the original one. Equivalently, we want to find the Fourier coefficients C e ∈ [0, 1] N the closest to C 0 , in a sense defined below. In the following, we will work in the Fourier domain.

Let C be any admissible kernel on Ω and X ∼ DPixP(C). As before, we will consider C either as a function from Ω to [0, 1], or as a vector in [0, 1] N . Recall that the intensity of the point process is given by

E(|X|) Ω = 1 Ω ξ∈ Ω C(ξ) = C(0).
In case of a kernel estimation from one sample, it is natural to consider that the expected cardinality of the point process to be estimated is the cardinality of this unique sample. Thus, a straightforward estimation of the intensity of the point process is

C e (0) = n N (9)
or equivalently

ξ∈ Ω C e (ξ) = n. Now, we want to determine the estimator C e (x), for all

x ∈ Ω \ {0} denoted Ω * . Let us consider p C (x) =      P(x ∈ X| 0 ∈ X) = P({0, x} ⊂ X) P(0 ∈ X) = C(0) - |C(x)| 2 C(0) if x = 0, 0 if x = 0. ( 10
)
Now, from the realization Y , we can obtain θ(x) the empirical estimator of p C (x) by

θ(x) =      1 n y∈Ω 1 Y (y)1 Y (y + x) if x = 0 0 if x = 0. ( 11
)
For optimization purposes, we express all the quantities in function of C e . In the following computations, we consider that the vectors are column vectors. Let us denote the set of admissible functions by

C n = { C ∈ R N such that ξ∈ Ω C(ξ) = n and ∀ ξ ∈ Ω, 0 ≤ C(ξ) ≤ 1}.
We are looking for C e such that

C e = argmin C∈ Cn p C -θ 2 2 = argmin C∈ Cn x∈Ω *   n N - N n |F -1 ( C)(x)| 2 - 1 n y∈Y 1 Y (y)1 Y (y + x)   2 = argmin C∈ Cn x∈Ω *   n 2 N 2 - 1 N y∈Y 1 Y (y)1 Y (y + x) -|F -1 ( C)(x)| 2   2 = argmin C∈ Cn x∈Ω * b(x) -g( C)(x) 2 = argmin C∈ Cn E( C),
where, for all C ∈ R N , and for all x ∈ Ω * ,

g( C)(x) = |F -1 ( C)(x)| 2 and b(x) = n 2 N 2 - 1 N y∈Ω 1 Y (y)1 Y (y + x). ( 12
)
We want to minimize E on C n a non empty closed convex set so we can use the projected gradient algorithm. To project on the set of constraints C n , we use a classic adapted version of the algorithm to project onto the simplex [28, p. 77-78], integrating a maximum bound constraint, denoted "proj". From a vector u ∈ R N , this mapping returns the vector proj(u) ∈ R N such that proj(u) x = min(max(u x -a, 0), 1), where a is chosen such that N x=1 proj(u) x = n. Let us compute the gradient of the energy E we want to minimize.

As

g : R N → R N -1 , C → |F -1 ( C)(x)| 2 x∈Ω * , we have ∀x ∈ Ω * , ∀ξ ∈ Ω, ∂g( C)(x) ∂ C(ξ) = 1 N F -1 ( C)(x)e 2iπ x,ξ + 1 N F -1 ( C)(x)e -2iπ x,ξ = 2 N Re F -1 ( C)(x)e -2iπ x,ξ ,
and moreover

∇E( C) = Dg( C) T 2 g( C) -b .
Notice that given a vector u = (u 0 , . . . , u N -1 ) T ∈ R Ω , we let u * be (u 1 , . . . , u N -1 ) T the restriction of u to Ω * . For all ξ ∈ Ω, Then

Dg( C) T u * ξ = 2 N x∈Ω * u x Re F -1 ( C)(x)e -2iπ x,ξ = 2 N Re   x∈Ω u x F -1 ( C)(x) e -2iπ x,ξ -u 0 C(0)   .
Dg( C) T u * = 2 N Re F u F -1 ( C) - 2n N 2 u 0
, where refers to the componentwise product of vectors. Finally we obtain

∇E( C) = 4 N Re F |F -1 ( C)| 2 -b F -1 ( C) , by setting b(0) = n 2 N 2 . ( 13
)
In particular, computing ∇E( C) only requires two FFT calls. The projected gradient descent algorithm is recalled and adapted to this problem in Algorithm 2. Note that the energy we want to minimize is not convex and it has several local minima: the initialization of the algorithm is crucial. Indeed, if the algorithm is initialized with a random matrix C init , the results can be far from the original target. We propose to initialize the algorithm with

C init = proj F √ b , ( 14
)
which is believed to be quite close to a solution of the optimization and provides good results, as observed in the experiments. Note that b can be negative, so applying a square root to b may produce complex coefficients to which we apply the Fourier transform. This enables the initialization kernel C init to be asymmetric. The step of the gradient descent is chosen as a constant depending on the order of magnitude of the gradient of the energy E. Figure 6 presents two initialization kernels computed from a given realization. Figures 7 and8 (column 3) provide some results of this algorithm, from one realization generated by different DPixP kernels. One realization seems enough to retrieve the Fourier coefficients of a simple symmetric projection kernel (see Figure 7, a, b whose non-zero Fourier coefficients form a convex set). Even though for most projection kernels a predominant shape appears in the estimation, as soon as the kernel is more complex, one sample does not provide enough information.

Learning the kernel of a DPixP from several realizations

A unique realization does not provide enough information for our proposed algorithm to estimate the Fourier coefficients of a DPixP kernels but if several realizations are available, combining them provides better results. Assume that we have J realizations, J ∈ N * , each of cardinality n j , that we suppose generated by the same DPixP kernel. • Compute for all

x ∈ Ω * , b(x) = n 2 N 2 -1 N y∈Y 1 Y (y)1 Y (y + x), b(0) = 0 (12).
• Set C 0 = C init [START_REF] Celis | Fair and diverse DPP-based data summarization[END_REF].

• for k = 1,. . . , k max -Compute ∇E( C k-1 ) ( 13).

-

Set C k = proj C k-1 -t∇E( C k-1 ) . Output: C K .
We propose to combine all the realizations to produce a better empirical estimator θ J of p C . First, the expected number of points is approximated by the mean number of points in the realizations, n =

n 1 + • • • + n J J .
If we have J realizations (Y i ) i∈{1,...,J} , Equation ( 11) is replaced by:

∀x ∈ Ω, θ J (x) =        1 nJ J i=1 y∈Ω 1 Y i (y)1 Y i (y + x) if x = 0, 0 if x = 0.
The rest of the procedure remains similar as we want to minimize the function p C -θ J 2 2 , in particular, the initialization kernel is

C init = proj   F    n 2 N 2 - 1 N J J i=1 y∈Ω 1 Y i (y)1 Y i (y + x)       .
Figures 7 and8 present some experiments on several DPixP kernels, using the second strategy presented here and combining all the samples in one estimation process. First, Figure 7 presents the results of this estimation procedure with projection kernels, meaning that the Fourier coefficients of these kernels are zero or one. It shows how 10 realizations provide enough information to retrieve a kernel close to the original one. Using 100 realizations enables to obtain satisfying results. We have seen in the previous subsection that any translation of the estimated Fourier coefficients or a symmetry with respect to (0, 0) of the estimated Fourier coefficients generate the same DPixP. Thus, in Figures 7 and8, we display a centered version of the estimation.

Figure 8 presents some results of this algorithm for non-projection DPixP kernels. Kernel a) is a Bernoulli kernel: all the Fourier coefficients are equal to n N . As expected, no specific structure appears from the estimation, regardless of the number of samples used. The estimations b) and c) are much noisier than their projection equivalent (Figure 7(a,e)) even if the shape formed by the Fourier coefficients (which directly impacts the local repulsion of the point process) seems retrieved. To conclude, the algorithm presented in this subsection provides satisfying estimations if the original kernel is a projection DPixP kernel, in particular when we have more than 10 samples. Indeed, as we have seen in Section 3.2 and as the authors of [START_REF] Biscio | Quantifying repulsiveness of determinantal point processes[END_REF] noted, projection determinantal processes can be seen as the most repulsive DPPs. Thus, within a sample, the characteristics of the repulsion, and of the kernel, are more accessible. Nevertheless, if we deal with a general complex kernel, the algorithm fails to retrieve meaningful information.This question of DPixP kernel estimation is complex and could be the subject of its own paper. Note that considering the quantity p C related to the pair correlation function imposes that we have only access to |C| thus different admissible kernels may be solution to this optimization problem. It would be necessary to study the theoretical properties of convergence, the bias and the variance of the estimator, the influence of the initialization or the existence of local minima, for instance.

Determinantal patch processes: DPPs on patches

As datasets to analyze and to process keep being larger and more complex, strategies to subsample these sets or to reduce the dimension of data have recently flourished. DPP subsampling is part of these approaches, as they enable to capture the structure of data and produce a representative subset of the whole initial set, taking into account its inner diversity. DPPs have been used in different applications such as text summarization [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF], feature selection in high dimensional data [START_REF] Belhadji | A determinantal point process for column subset selection[END_REF] or approximation of a dataset adapted to a given learning task [START_REF] Tremblay | Optimized algorithms to sample determinantal point processes[END_REF]. In image processing and computer vision, DPPs have raised interest through video summarization ( [START_REF] Gong | Diverse sequential subset selection for supervised video summarization[END_REF], [START_REF] Zhang | Video summarization with long short-term memory[END_REF]). The authors of [START_REF] Gong | Diverse sequential subset selection for supervised video summarization[END_REF] introduced sequential DPP to take into account both the diversity of the frames and the chronology of the video. To represent the diversity of the frames they use a decomposition similar to the quality-diversity decomposition that is introduced in [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF] and that we recall below. Furthermore, the paper [START_REF] Zhang | Video summarization with long short-term memory[END_REF] proposes a strategy enhanced by DPPs which makes it one of the state of the art methods for video summarization. This method also uses a decomposition similar to a quality-diversity decomposition to describe the diversity in the video.

In this section, we focus on subsampling the set of patches P of an image. This procedure can be useful for compression purpose for instance. It can also be necessary in order to fit a model on the patch set using only a proportion of the set, to increase the efficiency of the algorithm. For example, several patch-based denoising methods represent the patch distribution as a Gaussian mixture model ( [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF], [START_REF] Houdard | High-dimensional mixture models for unsupervised image denoising ( HDMI)[END_REF]). These methods rely on the estimation of the parameters of such models thanks to the Expectation-Maximization (EM) algorithm. To do so, in general, they randomly and uniformly select a subset of patches. This random selection is fast but the subset needs to be large enough so that it captures the patches diversity. The size of this selection impacts the running time of the estimation process, so a smaller selection, representative of the patches of the image, would ensure a faster and more accurate estimation. DPPs offer the opportunity to select a small subset of patches that captures the whole image. For instance, the authors of [START_REF] Launay | Determinantal patch processes for texture synthesis[END_REF] have used DPPs to subsample the patches of an image and approximate the patch distribution to speed up their texture synthesis algorithm.

Tremblay et al. [START_REF] Tremblay | Determinantal point processes for coresets[END_REF] compare several k-means initializations with one using DPP sampling in the context of coresets: the authors sample an appropriate DPP to select the initial centroids for the clustering strategy. Similarly, Agarwal et al. [START_REF] Agarwal | Notes on using determinantal point processes for clustering with applications to text clustering[END_REF] propose to adapt the k-means algorithm by using a DPP initialization and they prove that this initialization compares favorably with k-means++, the most popular adaptation of the k-means algorithm, whose initialization uses another negatively correlated sampling strategy, called D 2 sampling. One advantage of an algorithm using DPPs over the second is its adaptability concerning the number of clusters. Thus, in the previous example with denoising methods, DPPs could also provide a satisfying initialization to the EM algorithm. We investigate here the possible choices of DPP kernels for such applications, in order to subsample the patch space of an image.

DPP kernels to sample in the space of image patches

When considering determinantal processes on patches, the framework is more general than before: We are no longer dealing with stationary periodic point processes. We consider a Hermitian kernel K adapted to select diverse subsets of patches from an image, as set in Definition 1. The definition of this diversity depends on the problem we want to solve: for instance, compression, reconstruction of the image or initialization of the centroids of a clustering or the EM algorithm.

First, we recall that there exists a second characterization of DPPs, using a positive semi-definite matrix L [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF].

Definition 19

We consider Y = {1, . . . , N } and L a Hermitian matrix of size N × N such that L 0, then the random set X ⊂ Y defined by

∀A ⊂ Y, P(X = A) = det(L A ) det(I + L)
is a DPP with L-ensemble kernel L. We will denote X ∼ DPP L (L).

Defined in this manner, these DPPs are called L-ensembles. Recall that the initial definition (1), using the kernel denoted by K, requires that 0 K I. This L-ensemble definition doesn't need the constraint of bounding the eigenvalues of the kernel by one. This property is convenient to define a kernel, and a diversity model adapted to a specific problem. So this characterization is increasingly used in the machine learning community. That is also the definition we use in this subsection. However, note that this definition excludes the case of projection DPPs, which have a fixed cardinality. In subsection 5.2, we consider such a projection DPP kernel. Consider the following spectral decomposition of a DPP kernel

K, K = N k=1 λ k v k v * k .
Note that the definitions using the kernels K and the L-ensemble kernel L characterize the same DPP if and only if for all k ∈ {1, . . . , N }, 0 ≤ λ n < 1 and if

K = L(L + I) -1 = I -(I + L) -1 and conversely L = K(I -K) -1 . Hence, in this case, L = N k=1 λ k 1 -λ k v k v * k .
Note that if K has any eigenvalue equal to 1, the DPP can't be associated to an L-ensemble.

In the following, consider an image u and the initial set P = {P i , i = 1, . . . , N }, the set of its patches of size (2ρ + 1) × (2ρ + 1) × d, where ρ ∈ N and d is the number of color channel. Let us present some kernels that can be used to subsample the patches of this image.

Consider an expected cardinality n, so that in average, the number of selected patches is equal to n. In practice, this target cardinality can be imposed by normalizing the kernel matrices with a constant so note that in the following formulas, the matrices are defined up to a normalization constant, precised in the experiments section. A first type of L-ensembles that are commonly used ( [START_REF] Tremblay | Determinantal point processes for coresets[END_REF], [START_REF] Launay | Determinantal patch processes for texture synthesis[END_REF]) is the class of Gaussian kernels. Let us consider a Gaussian kernel based on the intensity of the patches, that we call the intensity Gaussian kernel, defined by

∀ P i , P j ∈ P, L ij = exp - P i -P j 2 2 s 2 , ( 15
)
where s is called the bandwidth or scale parameter. This kernel depends on the squared Euclidean distance between the intensity values of pairs of patches. It is often used as a similarity measure on patches. Despite its natural limitations, this similarity measure provides good results. The value of the parameter s has an impact on how repulsive the DPP is. Notice that if s is small, due to the exponential function, L ij converges very quickly to zero as soon as i = j and the distinction between patches is not very subtle. Thus, if s is small, L is close to the identity matrix and the DPP selection of patches is similar to a random uniform selection. On the contrary, for the same reason, the larger s is, the more repulsive the DPP is. However, this scale parameter should not be set too large because this would cause high numerical instability. As noticed in [START_REF] Aggarwal | Outlier Analysis[END_REF] and [START_REF] Tremblay | Determinantal point processes for coresets[END_REF], the median of the distances of intensities between the patches is a satisfying choice for setting the value of s.

We propose to compare this kernel with another Gaussian kernel that we call the PCA kernel, which depends on the squared distance between patches in the space given by keeping only the k principal components after a Principal Component Analysis (PCA). Set A the matrix gathering all the patches of the image reshaped in column so that the size of A is d(2ρ + 1) 2 × N . We assume that A has been centered, by subtracting the average patch to all the patches. It has not been reduced, meaning that patches with high variance, for instance patches with edges, will highly influence the decomposition. Thanks to a singular value decomposition, consider U, V two unitary matrices and Σ a diagonal matrix storing the principal values of A such that A = U ΣV T . We choose to keep only k principal components and we obtain the matrix A k = V k A, where we kept only the k first rows of the matrix

V in V k of size k × d(2ρ + 1) 2 and the matrix A k = {P k i , i = 1, . . . , N } is k × N . Every initial patch P i ∈ A is associated with a projected vector P k i ∈ A k .
Thus, the PCA kernel is defined by

∀ P i , P j ∈ P, L ij = exp - P k i -P k j 2 2 s 2 . ( 16
)
This method discards principal vectors associated to small singular values and projects the patches on a low-dimensional space associated with the large singular values. This enables to find the components that best represent the variance of the patches and ignore mainly noise (depending on the number of dimension discarded). Thus, comparing patches in this low-dimensional space seems relevant to capture more precisely their dissimilarity.

A second type of common L-ensemble kernels uses a quality-diversity decomposition of the data. Kulesza and Taskar present in [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF] this decomposition that uses a given quality measure computed on each element of the set and a dissimilarity computed between pairs of elements. In this paper, each patch P i is associated with a quality measure, which is a non-negative number q i = q(P i , P) ∈ R + , depending on the patch itself and on the other patches. Each patch P i is also associated with a feature vector φ i = φ(P i ) ∈ R D , such that φ i 2 = 1, which depends only on the patch itself. The quality/diversity kernel L is defined by

∀ P i , P j ∈ P, L ij = q i φ T i φ j q j . ( 17
)
This class of kernels presents several advantages. The first advantage of this definition is its interpretability. Each patch is associated with a quality measure, that one can adapt depending on the characteristics one wants to favor. The comparison between patches is also accessible and adjustable to obtain the most adapted kernel. This decomposition has a second advantage: the kernel becomes a low-rank matrix, with a rank equal at most to D, the number of features. In case of low-rank kernels, Kulesza and Taskar [START_REF] Kulesza | Learning determinantal point processes[END_REF] propose a dual representation and a dual sampling algorithm. This sampling scheme is equivalent to the original algorithm but it takes advantage of the low-rank kernel and becomes much faster. Note that, whatever the DPP L-ensemble kernel, the cardinality of a sample generated from DPP L (L) will necessarily be lower than the rank of L. This low-rank definition imposes to sample subsets of smaller size than D, the number of features computed from the patches. Thus, this kernel is adapted when small and very small subsets of patches are needed. In these cases, it is very important to precisely control the selection process so such kernels are particularly relevant.

For the kernel that we call Qual-div kernel, we associate each patch with a feature vector given by a discrete cosine transform of the patch. Thus, each feature vector is of size d(2ρ + 1) 2 . Note that in the experiments, we use color images (with 3 color channels) and patches of size 7 × 7 (meaning that ρ = 3) so the feature vectors of length 147. We define the quality measure such that it attributes a high value to patches whose intensity is far from that of its neighbors in the pixel grid. For each pixel, we count the number of neighboring pixels whose intensity is further in Euclidean distance than a certain constant. This constant is chosen equal to half the maximum Euclidean distance between the intensities in the image. This choice gives further priority to singular patches, that can be seen as the outliers of the set of patches. As experiences will show, it highly favors textures and edges.

Minimizing the selection error

The question is to choose the best kernel, such that the sampled DPP on the patches minimizes an error computed as a distance between the selected patches and the initial set of patches P. This last kernel will be called the optimized kernel. This problem is similar to discrete optimal quantization problems [START_REF] Pagès | Introduction to vector quantization and its applications for numerics[END_REF] where the aim is to find the best subset of patches Q such that E Q∼µ (d(Q, P)) is minimal, for a given distance d. Yet, this computation is often costly and hardly tractable. In the following, we suppose that the patches are of size (2ρ + 1) × (2ρ + 1) for some positive integer ρ and we denote by ω ⊂ Z 2 the patch domain {-ρ, . . . , ρ} 2 .

First, the error, or the distance between the sample and the initial set of points, we want to minimize depends on the application. The mean square error (MSE in short) is commonly used to compare an image and its reconstruction. Here, we use a similar distance, the squared L 2 norm between the patches of the image and their nearest neighbor in the selection given by the DPP sampling on the patches. Consider Q a subset of patches. This error is defined by

E 1 = 1 N N i=1 d L 2 (P i , Q) 2 = 1 N N i=1 min Q∈Q x∈ω (P i (x) -Q(x)) 2 , ( 18
)
where ω is the patch domain. One hopes that using a DPP to generate Q will prevent from concentrating only on the most common patches and select singular patches. Given an expected cardinality n ≤ N and a kernel K n , we would consider Q ∼ DPP(K n ). The following error can be useful to verify this property:

E 2 = max i∈{1,...,N } d L 2 (P i , Q) 2 = max i∈{1,...,N } min Q∈Q x∈ω (P i (x) -Q(x)) 2 . ( 19
)
A low error value asserts that the outlier patches (non redundant) are selected. Given an expected cardinality n ≤ N and a kernel K n , we will consider Q ∼ DPP(K n ). We would like to find the DPP kernel minimizing the expectation of the errors: E Q∼DPP(Kn) (E 1 ) and E Q∼DPP(Kn) (E 2 ). Yet, this optimization problem depending on a DPP matrix K n is intractable. As in the papers by Kulesza and Taskar [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF] and Affandi [START_REF] Affandi | Learning the parameters of determinantal point process kernels[END_REF], we would like to have a closed-form minimization problem to obtain optimal parameters. These strategies are based on the quality-diversity decomposition of an L-ensemble kernel described in the previous section. Given predetermined features vectors, they determine an appropriate quality measures from the data. Here, we use a similar parametrization, using the first definition of DPPs, with a kernel matrix K. We suppose that its eigenvectors are fixed (given by features computed from the patches of the image) and we want to determine the optimal spectrum so that the associated matrix K minimizes a tractable error. Furthermore, thanks to the Campbell Formula (6), we know that the expectation of some functionals defined on point processes are tractable. That is what we use in the following.

Suppose we select a subset of patches using a DPP of kernel K: Q ∼ DPP(K). We would like to study the following measure:

R(Q) = P ∈P Q∈Q f P (Q). ( 20
)
It can be seen as a reconstruction evaluation, if the function f P involves a distance between the input patch and the patch P . With the appropriate function f P , R can represent how well a patch P ∈ P is represented by the selection Q. For instance, by considering the functions f α,P (Q) = 1 P -Q 2 ≤α or f P (Q) = e -P -Q 2 , R will return a high value if the selection Q encompasses the set of patches. Notice that if we use a function f p which depends on the L 2 distance between patches, maximizing R will favor selections similar to the ones minimizing the MSE. Thus, contrary to the previous error quantities, E 1 and E 2 , we want to generate a subset Q such that R is large. From the Campbell Formula (6) adapted to general discrete DPPs, we have

E(R(Q)) = E   P ∈P Q∈Q f P (Q)   = N j=1 E   Q∈Q f P j (Q)   = N j=1 N i=1 f P j (P i )K(P i , P i ).
Assume that K admits the eigendecomposition

K(P i , P j ) = D k=1 λ k φ k (P i )φ * k (P j ), (21)
with D ≤ N , and unknown eigenvalues (λ k ) k∈{1,••• ,D} . We suppose that the eigenvectors are fixed. In the experiments, (φ k ) k∈{1,••• ,D} are the same feature vectors as defined in the previous subsection, given by the discrete cosine transform of each patch. Then the previous expectation becomes

E(R(Q)) = D k=1 λ k N i=1 |φ k (P i )| 2 N j=1 f P j (P i ) (22)
The maximization of this quantity with respect to (λ 1 , . . . , λ D ) is a linear problem, under the linear constraints: The advantage of solving such a problem is that the solution (λ * k ) k∈{1,••• ,D} is explicit. It is on the boundary of the constraints, meaning that is a kernel K with only n non-zero eigenvalues, each one equal to 1: K is a projection DPP. Given any function f p , any integer n ≤ D, let us consider I n the set of the indices associated to the n largest coefficients of the vector ψ of size D defined by

ψ k = N i=1 |φ k (P i )| 2 N j=1 f P j (P i ).
The solution of the problem argmax

(λ k ) E (R(Q)) such that D k=1 λ k = n and ∀k, 0 ≤ λ k ≤ 1, (23)
is the set of eigenvalues (λ * k ) k=1,...,D defined by

λ * k = 1 if k ∈ I n 0 otherwise . ( 24
)
For instance, if we choose f α,P i (P j ) = 1 P i -P j 2 ≤α , then we need to maximize the function

E(R(Q)) = D k=1 λ k N i=1 |φ k (P i )| 2 N j=1 1 P i -P j ≤α = D k=1 λ k N i=1 |φ k (P i )| 2 |B(P i , α)|, ( 25
)
where B(P, α) is the ball with center P and radius α for the Euclidean distance between patch intensities, and |A| is the cardinality of the subset A. Thus, |B(P i , α)| denotes the number of patches in the image that are within a distance of P i smaller than α. In the experiments, we use this function and we choose α to be half the median of interdistances between patches. Note that this maximization problem generates a kernel called the optimized kernel that will favor patches similar to many others. This creates an interesting compromise: the DPP will tend to select diverse subsets of redundant patches. As anticipated, we will see in the experiments that this method tends to miss singular patches.

Experiments

The following figures present some results of subsampling in the space of image patches, for different cardinality. First notice that the cardinality is fixed for the uniform sampling.

It is also fixed for the last optimized kernel, as we obtain a projection kernel from the maximization problem. Concerning the three other kernels, they are defined using the L-ensemble definition in Equations ( 15), ( 16) and ( 17). We used a common normalization strategy, formalized in [START_REF] Barthelmé | Asymptotic equivalence of fixed-size and varying-size determinantal point processes[END_REF], using a L-ensemble kernel L whose eigenvalues are denoted (λ k ) k∈{1,...,N } . Given a desired expected cardinality n, we normalize L to obtain a kernel

L c = cL, where c is chosen such that N k=1 cλ k 1 + cλ k = n.
Note also that the Qual-div kernel and the optimized kernel are low-rank, with a rank equal at most to the number of features that we use to defined the kernels. In these experiments, the feature vector associated to each patch (φ in Equations ( 17) and (5.2)) is obtained from the discrete cosine transform of the patch. Note that a DPP kernel can't generate samples with more items than its rank and in the following experiments, we use patches of size 7 × 7 × 3. Thus, the rank of the two previous kernels is 147 and we can observe the results, with a step of 50, up to a cardinality equal to 100 in Figure 12.

Figures 10 and 11 show images reconstructed using the associated selected patches presented below the reconstruction. Each patch in the initial image is replaced by its nearest neighbor in the DPP selection. The final image is obtained by average: given a pixel, all the overlapping patches containing this pixel are averaged. This is a common strategy to aggregate the patches. Several other methods are proposed in the literature, such as using a weighted average [START_REF] Dabov | Image denoising by sparse 3-D transform-domain collaborative filtering[END_REF][START_REF] Salmon | From patches to pixels in non-local methods: Weighted-average reprojection[END_REF] or implicitly including the reconstruction in a global variational problem [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]. An average considering uniform weights on all the patches is often used as it does not require any other computation or information to store. Thus, after subsampling the set of patches, the initial image can be represented by its size N 1 × N 2 = N , the small set of patches of size (2ρ + 1) × (2ρ + 1) × d and a vector of indices of length N , associating each initial patch to its nearest neighbor in the selection. Figure 12 compares the errors E 1 (18), E 2 [START_REF] Philip | Circulant matrices[END_REF] and the peak signal-to-noise ratio (PSNR) of the reconstruction images generated from samples given by the different kernels. The PSNR is a metric commonly used to evaluate the quality of the reconstruction of an image. Consider an initial image I 0 and a reconstruction I 1 , both having d color channels and N pixels with a value between 0 Card Unif. sample Intens. kernel PCA kernel Qual-div. kern. Optim. kern. Note that the curves associated to the "Qual/div" and the optimized kernels stop at an expected cardinality equal to 100 selected patches since, for these two kernels, the number of features chosen to describe the patches (equal to 147 in the experiments) determines the maximal size of the generated subset of patches. .

First, as expected, a uniform sampling can produce samples which contains many similar patches. The first image (Pool) has several large and regular regions that could be represented by a few patches and these regions are often over-represented in the results. Note that when we compare the kernels using the error E 1 , in particular for the second image (Parrot), the uniform selection provides satisfying results. On the contrary, small and rare details are often missed by the uniform sampling, and the graph of Figure 12 shows that this sampling strategy compares badly with the others when considering this criteria. Furthermore, the graph presenting the PSNR results illustrates how this uniform strategy provides overall poorer reconstructed images.

Note that the optimized kernel, making a compromise between the diversity induced by DPPs and the redundancy imposed by maximizing the chosen reconstruction error [START_REF] Descamps | Automatic flamingo detection using a multiple birth and death process[END_REF], produces quantitative results similar to a uniform sampling. When observing the patches selected by this kernel in Figures 10 and11, one can see that this kernel tends to select slightly more diverse patches than a uniform sampling.

Second, the PCA kernel and the Qual-div kernel behave rather similarly. They tend to favor singular patches and patches containing edges, even sometimes over-representing them. Thus, they provide good results when looking at the second error measuring the distance between the selection and the furthest patch, especially the PCA kernel. Yet, they can provide even worst results than the uniform selection when we look at the average distance between the selection and the initial set of patches (Error E 1 (18)).

Finally, the Intensity kernel, using only the squared Euclidean distance between intensities, seems to be the most stable kernel. It provides small average error and tends to include singular patches in the selection. For both images, whatever the expected cardinality, the samples generated by this kernel produce visually satisfying reconstructions.

Thus, the choice of subsampling strategy in the patch space of an image highly depends on the purpose of the generated selection. The most stable strategy seems to be using the Intensity kernel [START_REF] Chen | Dimensionality reduction based on determinantal point process and singular spectrum analysis for hyperspectral images[END_REF], which provides a selection close in average to the initial patches and which selects also singular patches. If the priority of the application is efficiency, the best strategy may remain to use a uniform selection with a high number of selected patches. Yet, if the size of the selection needs to be low or if the selection needs to contain mainly structure and texture information, the good choice may be to use a PCA kernel or a kernel using the quality-diversity decomposition. Moreover, note that when the kernel is low-rank, it is also possible to use a dual representation of the kernel and to speed up the simulation of the process [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF][START_REF] Gartrell | Low-rank factorization of determinantal point processes[END_REF]. In this paper, we present several low-rank kernels: the quality-diversity kernel and the optimized kernel. In these cases, the computational time can be reduced.

Conclusion

In this paper, we presented several applications of determinantal point processes to image processing. We defined these models on the pixels of an image and we call such point processes determinantal pixel processes (DPixP). In this setting, we showed that the only possible hard core repulsion for DPixP is directional. Given a direction, it is possible to impose to select at most one pixel on any discrete line in the image, but any further hard core constraint leads to a degenerate kernel. We studied shot noise models based on DPixP as a method to sample homogeneous micro-textures and adapted the choice of DPixP kernel in function of a given spot function of the shot noise and of the homogeneity we look for. We also developed an algorithm to infer the Fourier coefficients of a DPixP kernel from a sample or from a set of samples. As a future work, we plan to investigate the estimation, from a texture image, of the spot function and of the DPixP kernel associated to a shot noise that could have generated the texture. We also considered DPPs defined on the set of patches of an image, in order to subsample this set while retaining the information necessary to appreciate the image. It would be interesting to use more features to describe the patches of an image, in order to improve the results of the Quality-diversity kernel and nuance the selection it generates, which for now tends to over-represent patches with textures or edges. We also would like to explore different functions used in the minimization of the selection error, depending on the purpose of the subsampling of the set of patches. At last, we intend to apply these selection strategies to patch-based methods and evaluate the gain in quality and computation time generated by the use of DPPs.

Example 1 Consider Ω = {0, 1} × {0, 1, 2}, with N 1 = 2 and N 2 = 3. The image domain is Ω = (0,0) (0,1) (0,2) (1,0) (1,1) [START_REF] Affandi | Learning the parameters of determinantal point process kernels[END_REF][START_REF] Agarwal | Notes on using determinantal point processes for clustering with applications to text clustering[END_REF] .

This image domain is associated to the following ordering of pixels Ω = 1 3 5 2 4 6 .

A kernel matrix K defined on such image domain will be defined using the previous ordering. Consider the function C defined on Ω and extended by periodicity to Z 2 as ∀x, y ∈ Ω, K(x, y) = C(x -y).

Moreover, due to the periodicity assumption, we know for instance that C((-1, -2)) = C((1, 1)) or C((-1, -1)) = C((1, 2)). We obtain the following kernel, which is indeed a 6 × 6 circulant matrix with 2 × 2 circulant blocks: It can be entirely characterized by its first column.

K =         

A.2 Diagonalization of a DPixP kernel

The following results on the diagonalization of circulant matrices and its generalization on block-circulant matrices with circulant blocks can be found in the book of Davis, Circulant matrices [START_REF] Philip | Circulant matrices[END_REF]. See [START_REF] Philip | Circulant matrices[END_REF][START_REF] Olson | Circulant matrices and their application to vibration analysis[END_REF] for more details and illustrations. First, it is well known that any circulant matrix K = circ(c 0 , . . . , c N -1 ) is diagonalized in the Fourier basis and its eigenvalues λ 0 , . . . , λ N -1 are given by ∀ ξ = 0, . . . , N -1,

λ ξ = N -1 x=0 c x e -2iπ xξ N .
Moreover, if the matrix K is Hermitian, its eigenvalues are real and most of them are constrained [START_REF] Philip | Circulant matrices[END_REF]: λ 0 is distinct but the remaining eigenvalues are such that λ ξ = λ N -ξ , if N is odd. If N is even, λ N/2 is also distinct.

One can obtain a similar result concerning block circulant matrices with circulant blocks. Consider ω N = e 2iπ N , the primitive N -th roots of unity, and the following N × N unitary matrix . . . ω

U N = 1 √ N         
(N -1) 2 N          . ( 27
)
This matrix is associated to the inverse discrete Fourier transform.

Proposition 20 (Diagonalization of a DPixP kernel) Consider a DPixP defined on Ω with kernel K, the block-circulant matrix with circulant blocks whose first column is C : Ω → C. Then, K is diagonalized in the Fourier basis, in the sense that

(U * N 2 ⊗ U * N 1 )K(U N 2 ⊗ U N 1 ) = Λ,
where ⊗ is the Kronecker product, and where Λ is diagonal such that

∀ξ ∈ Ω, Λ ξ,ξ = C(ξ) = x∈Ω C(x)e -2iπ x 1 ξ 1 N 1 + x 2 ξ 2 N 2
.

Proof Let us prove that we have

K(U N 2 ⊗ U N 1 ) = (U N 2 ⊗ U N 1 )Λ. Note that ∀x, ξ ∈ Ω, (U N 2 ⊗ U N 1 ) x,ξ = 1 √ N 1 N 2 e 2iπ x 1 ξ 1 N 1 + x 2 ξ 2 N 2 = 1 √ N e 2iπ
x,ξ , with the notation • that we use throughout the paper. Given x, ξ ∈ Ω, we have

K(U N 2 ⊗ U N 1 )
x,ξ First, let us study the 1D case, meaning that K is a kernel defined on the points of Y = {0, . . . , N -1} (to be consistant with our 2D representation) and it is circulant. Therefore, for all i, j ∈ Y, K(i, j) = c j-i = c i-j . As K is not irreducible, there exist i, j ∈ Y, such that K(i, j) = c j-i = 0. Let us denote k = inf{l > 0 such that c l = 0}, hence

c 1 = • • • = c k-1 = 0 = c -1 = • • • = c -k+1
. Notice that k is necessarily larger or equal to 2, otherwise by the circulant property, K would not have any zero coefficient and it would be irreducible. Similarly, k necessarily divides N and the only non-zero coefficients c m are multiples of k, as otherwise the non-zero elements of K would be located such that it would be possible to access to any index from any other by traveling only through non-zero coefficients: K would be irreducible. Then, if we define l such that N = k × l, there are k cycles of size l in the graph associated to K, each block with the same l coefficients {c k , c 2k , . . . , c lk }, or equivalently, ∀i 0 = 0, . . . , N -1, K(i 0 , j) = c kp , if j = kp + i 0 mod N, with p = 0, . . . , l -1, 0, otherwise.

Thus it is possible to define the permutation matrix P which gathers the cycles, and which associates K with a block diagonal matrix: ∀p = 0, . . . , l -1, ∀r = 0, . . . , k -1, P (p + lr, r + pk) = 1 meaning that P moves the coefficient K r+pk to the coordinate p + lr. Moreover, these blocks (B r ) r∈{0,...,k-1} are circulant: for all r = 0, . . . , k -1, for all i, i = 0, . . . , l -1, B r (i, i ) = K(r + ik, r + i k) = c (i-i )k , for all τ ∈ Y such that (i + τ mod N ) and (i + τ mod N ) are in the r-th cycle, B r (i + τ, i + τ ) = K(r + (i + τ )k, r + (i + τ )k) = c (i-i )k = B r (i, i ).

To conclude, K is permutation similar to a block-diagonal matrix, which is the repetition of one irreducible, circulant and Hermitian block. Now let us consider the 2D case, when K is a kernel matrix defined on Ω = {0, . . . , N 1 -1} × {0, . . . , N 2 -1} and assume that K is Hermitian, block-circulant with circulant blocks and completely reducible. Define C the function such that for all (i, j), (i , j ) ∈ Ω, K ((i, j), (i , j )) = C(i -i, j -j). As in the 1D case, define (e 1 , e 2 ) ∈ Z 2 ∩ Ω the two generating vectors such that C(r, s) = 0, ∀(r, s) inside the elementary cell generated by (e 1 , e 2 ). These two vectors generate a subgroup of Z 2 and it contains Z(0, N 2 ) + Z(N 1 , 0), as K is not irreducible and similarly as in the 1D case. Then e 1 divides N 1 , e 2 divides N 2 . As before, the only non-zero coefficients of C belong to {Ze 1 + Ze 2 } ∩ Ω. The size of the elementary cell determines the number of cycles (and future blocks) and l = {Ze 1 + Ze 2 } ∩ Ω defines the size of each cycle. It is possible to define the permutation matrix that transforms K into a blockdiagonal matrix with irreducible blocks. For all (i, j) ∈ Ω, let us define (r, s) its representative element in the elementary cell such that there exists p, q such that (i, j) = (pe 1 + qe 2 ) + (r, s) mod (N 1 , N 2 ). Thus, P associates the coefficient (i, j) = (pe 1 + qe 2 ) + (r, s) mod(N 1 , N 2 ) with the coordinate (p, q) + (r, s) (block (r, s), coefficient (p, q)). As before, the blocks

Proposition 1

 1 The cardinality |X| of the DPP X is distributed as the sum of N independent Bernoulli random variables: |X| ∼ x∈Y Ber(λ x ), where the Bernoulli variables take the value 1 with probability λ x . Then E(|X|) = x∈Y λ x = Tr(K) and Var(|X|) = x∈Y λ x (1 -λ x ).

  clearly, for all ξ ∈ Ω, C(ξ) is real and 0 ≤ C(ξ) ≤ 1. C is an admissible kernel. Let us consider two fundamental examples of DPixPs. The first one is the Bernoulli process. It corresponds to the discrete analogous of the Poisson point process: points are

Figure 1 :

 1 Figure 1: Comparison between samples (both have 148 points) of two DPixPs: a Bernoulli process (first line) and a projection DPixP defined by the kernel C such that C is the indicator function of a discrete circle (second line). For both DPixPs, from left to right, the Fourier coefficients C of the kernel, the real part of the kernel Re(C) and one associated sample.

Figure 2 :

 2 Figure 2: Example of a kernel associated with hard core repulsion in the horizontal direction.From left to right, the Fourier coefficients of C, the real part of the kernel C, a capture of the conditional density during the simulation, the associated final sample.

Figure 2

 2 Figure 2 illustrates this proposition: all non-zero Fourier coefficients are vertically aligned.The third figure presents a capture of the conditional density while the simulation is in progress, after 15 pixels already sampled. In each pixel, the probability that it is the next point selected is represented by the gray scale: the lighter a pixel is, the greater its probability of being the next point sampled. One can see that as soon as a pixel x is sampled, all the pixels belonging to the horizontal line passing through x have a zero probability of being sampled next. Note that this proposition is not limited to the horizontal or vertical cases.

Figure 3 :

 3 Figure 3: Realizations of the shot noise model driven by several spot functions and the most repulsive DPixP adapted to this spot. From left to right: the spot function, the Fourier coefficients obtained by our glutton algorithm, the real part of the associated kernel C, a sample of this most repulsive DPixP and a sample of the associated shot noise model and finally a Bernoulli shot noise model, both having the same expected number of points (n = 80).

Figure 4 :

 4 Figure 4: Determinantal shot noise realizations S M as defined in Theorem 12 with various M = 1, 2, 3, 6 and a comparison with their associated limit Gaussian random field N (0, Σ(C)) (f). The shot noise is based on the spot (a) and the DPixP with kernel C whose Fourier coefficients form an isotropic 2D Gaussian function (Figure 1, bottom).

Figure 5 :

 5 Figure 5: Three DPixP kernel functions, defined by their Fourier coefficients, generating the same DPixP.

Proposition 14 (

 14 Equivalence classes of real symmetric kernels[START_REF] Kulesza | Learning with Determinantal Point Processes[END_REF]) Let K 1 and K 2 be two real positive symmetric N × N matrices with eigenvalues bounded by 1. Then DPP(K 1 ) = DPP(K 2 ) if and only if there exists a N × N diagonal matrix D such that K 2 = D -1 K 1 D, where the coefficients of D are either 1 or -1.

Figure 6 :

 6 Figure 6: Two examples of initialization of our estimation algorithm. From left to right: the Fourier coefficients of the original kernel to retrieve (a), the initialization from one realization (b), the resulting estimated kernel (c), the initialization from 10 realizations (d), the resulting estimation from 10 realizations (e).

Figure 7 :Algorithm 2

 72 Figure 7: Experiments on several projection kernels. From left to right: the target Fourier coefficients of the kernel we want to recover, one realization of this DPixP, the estimation of the Fourier coefficients from one, from 10 and from 100 realizations, with k max = 2000.

Figure 8 :

 8 Figure 8: Experiments on general DPixP kernels. From left to right: the target Fourier coefficient of the kernel we want to recover, one realization of this DPixP, the estimation of the Fourier coefficients from one, from 100 and from 800 realizations, with k max = 2000.

λ

  k = n, and for all k ∈ {1, . . . , D}, 0 ≤ λ k ≤ 1.

Figure 9 :

 9 Figure 9: Original images considered in Figures 10 to 12.

  Figure 9: Original images considered in Figures 10 to 12.

5 25 100Figure 10 :Figure 12 :

 5251012 Figure 10: Image reconstruction comparing different expected cardinality and the DPP kernels presented in the previous subsections. For each cardinality, the first row presents the reconstruction of the image using only the patches selected by the corresponding kernel, given in the second row.

  (i, c) -I 1 (i, c)) 2

  C((0, 0)) C((1, 0)) C((0, 2)) C((1, 2)) C((0, 1)) C((1, 1)) C((1, 0)) C((0, 0)) C((1, 2)) C((0, 2)) C((1, 1)) C((0, 1)) C((0, 1)) C((1, 1)) C((0, 0)) C((1, 0)) C((0, 2)) C((1, 2)) C((1, 1)) C((0, 1)) C((1, 0)) C((0, 0)) C((1, 2)) C((0, 2)) C((0, 2)) C((1, 2)) C((0, 1)) C((1, 1)) C((0, 0)) C((1, 0)) C((1, 2)) C((0, 2)) C((1, 1)) C((0, 1)) C((1, 0)) C((0, 0))

Example 2 .

 2 -y)e 2iπ y,ξ = 1 √ N y∈Ω C(y)e 2iπ x-y,ξ = (U N 2 ⊗ U N 1 ) x,ξ C(ξ) = (U N 2 ⊗ U N 1 )Λ x,ξ, where Λ = diag ( C(ξ), ξ ∈ Ω) .A necessary and sufficient condition for the DPixP to be defined is that the eigenvalues of the kernel are in [0, 1]. In particular, the Fourier coefficients of C are real. Thus we have∀x ∈ Z 2 , C(-x) = C(x). Let us use the same illustration as in Example 1. We can deduce from the previous remark and from the periodic boundaries that the kernel is of the shape0)) C((1, 0)) C((0, 1)) C((1, 1)) C((0, 1)) C((1, 1)) C((1, 0)) C((0, 0)) C((1, 1)) C((0, 1)) C((1, 1)) C((0, 1)) C((0, 1)) C((1, 1)) C((0, 0)) C((1, 0)) C((0, 1)) C((1, 1)) C((1, 1)) C((0, 1)) C((1, 0)) C((0, 0)) C((1, 1)) C((0, 1)) C((0, 1)) C((1, 1)) C((0, 1)) C((1, 1)) C((0, 0)) C((1, 0)) C((1, 1)) C((0, 1)) C((1, 1)) C((0, 1)) C((1, 0)) C((0, 0))Note that any DPixP kernel defined on a grid Ω of size 2 × 3 is entirely characterized by four coefficients.

  Proof Define K 1 and K 2 two admissible DPP kernels, such that K 1 verifies the hypothesis of Theorem 16. By definition, DPP(K 1 ) = DPP(K 2 ) is equivalent to K 1 and K 2 having equal principal minors. In the papers[START_REF] Hartfiel | On matrices having equal corresponding principal minors[END_REF] (Theorem 7) and[START_REF] Loewy | Principal and diagonal similarity of matrices[END_REF] (Theorem 1), Hartfiel and Loewy prove that if K 1 is irreducible and for every partition of Y into two subsets, α and β such that |α| ≥ 2 and |β| ≥ 2, rank (K 1 ) α×β ≥ 2, then K 1 and K 2 have equal principal minors if and only if there exists a diagonal matrix D such that K 2

  Notice that when we consider two equivalent DPixP kernels C 1 and C 2 , generating the block-circulant matrices K 1 and K 2 , there are three possible configurations. The first one is when K 1 verifies the assumptions of Theorem 16, it leads to Proposition 17. In the second case, K 1 is irreducible, but N = N 1 N 2 ≥ 4 and there exists a partition α, β of Y such that |α| ≥ 2, |β| ≥ 2 and rank (K 1 ) α×β < 2. In the third case, K 1 is not irreducible. Assume that K 1 is irreducible, N = N 1 N 2 ≥ 4 and that there exists a partition α, β of Y such that |α| ≥ 2, |β| ≥ 2 and rank(K 1 ) α×β

	Let us characterize the second and third cases. It appears that these configurations are "rare"
	in practice.
	Case 2:

the proof is identical.

Remark 18

Acknowledgments

This work was supported by grants from Région Ile-de-France. We thank the reviewers for their valuable comments and suggestions that helped us to improve the paper.

Appendix A. DPixP kernels as block-circulant matrices with circulant blocks

Recall that a N × N circulant matrix K is entirely characterized by its generating elements c 0 , . . . , c N -1 such that K = circ(c 0 , . . . , c N -1 ), in the sense that K(x, y) = c (x-y mod N ) :

Similarly, a

. If each generating matrix K i is circulant, K is called a block-circulant matrix with circulant blocks. In that case, the matrix K can be characterized by the generating elements of all the generating matrices.

In our framework, the pixels are considered column by column so that the ordered index

Moreover, we suppose that the process is stationary and periodic. Hence we can define a function C : Ω → C such that ∀x, y ∈ Ω, K(x, y) = C(x -y), extended by periodicity to Z 2 . In this configuration, the matrix K is block-circulant with circulant blocs. To simplify the notations in the paper, we chose to identify the pixel's position in the image and its ordering.

A.1 Illustration on a small image domain

Let us present a simple illustration.

Appendix B. Identifiability of a DPixP

B.1 Case 2: K 1 is irreducible and doesn't verify the rank hypothesis (Theorem 16)

Let us study the equivalence class of a DPixP of kernel C 1 such that its associated matrix K 1 is irreducible and there exists a partition α, β of Y such that rank(K 1 ) α×β = 1, in the case Ω of size 1 × 5: Once again, we know that C 1 (0) = C 2 (0) = C 0 and thanks to the equality of principal minors of size 2, the modulus are equal so r 21 = r 11 = r. One of the principal minors of size 3 for C 1 is equal to Finally, let us assume we are in the first case, K 1 can be written

kπ , e i 4 5 kπ , e -i 4 5 kπ , e -i 2 5 kπ , which corresponds to a translation of the Fourier coefficients of C of k pixels. The second case yields to K 1 = DK 2 D -1 which corresponds to the symmetry and the translation of k pixels of the Fourier coefficients of C.

B.2 Case 3: K 1 is not irreducible

In this appendix, we consider a Hermitian block-circulant matrix K of size N × N that is completely reducible, meaning that it is permutation similar to a block diagonal matrix with irreducible blocks. We want to prove that in that case, the blocks are identical, that is they are of equal size and they are composed of the same coefficients. Moreover, we prove that these blocks are not only irreducible but also Hermitian and circulant.

(B (r,s) ) have the same size and have an identical structure. Let us consider the block (r, s), consider (i, j), (i , j ) ∈ Ω, B (r,s) ((i, j), (i , j )) = K((pe 1 + qe 2 ) + (r, s) mod(N 1 , N 2 ), (p e 1 + q e 2 ) + (r, s) mod(N 1 , N 2 )) = C((p -p)e 1 + (q -q)e 2 ) Let (τ 1 , τ 2 ) ∈ Ω be such that (i + τ 1 , j + τ 2 ), (i + τ 1 , j + τ 2 ) both belong to the cycle (r, s).

= K((pe 1 + qe 2 ) + (r, s) + (t 1 e 1 + t 2 e 2 ) mod(N 1 , N 2 ), (p e 1 + q e 2 ) + (r, s) + (t 1 e 1 + t 2 e 2 ) mod(N 1 , N 2 )) = C((p -p)e 1 + (q -q)e 2 ) = B (r,s) ((i, j), (i , j )).

Thus, for all (r, s), the associated bloc B (r,s) is block circulant with circulant blocks. Similarly, it is Hermitian. To conclude, K is permutation similar to a block diagonal matrix with one repeated irreducible, circulant, Hermitian block.