Asymptotic properties of maximum likelihood estimators for determinantal point processes

Yaozhong Hu Supported by an NSERC Discovery grant and a centennial fund from University of Alberta. Email: [email protected] Department of Mathematical and Statistical Sciences, University of Alberta, Edmonton, AB, Canada. Haiyi Shi Corresponding author. Email: [email protected] Department of Statistics and Actuarial Science, Simon Fraser University, Burnaby, BC, Canada.

Abstract

We obtain the almost sure strong consistency and the Berry-Esseen type bound for the maximum likelihood estimator $\tilde{L}_{n}$ of the ensemble $L^{\star}$ for determinantal point processes (DPPs), strengthening and completing previous work initiated in Brunel, Moitra, Rigollet, and Urschel [BMRU17]. Numerical algorithms of estimating DPPs are developed and simulation studies are performed. Lastly, we give explicit formula and a detailed discussion for the maximum likelihood estimator for blocked determinantal matrix of two by two submatrices and compare it with the frequency method.

1 Introduction

Determinantal point processes (DPPs) arise from random matrix theory [Gin65] and are first introduced to give the probability distribution of Fermionic system in thermal equilibrium in quantum physics [Mac75]. Since then, DPPs have been found in various aspects of mathematics, including for example, loop-free Markov chains [BDF10], edges of uniformly spanning trees [BP93] and in particular, recently to image processing [LDG21].

In the seminal work [KT12], Kulesza and Taskar show that DPPs demonstrate the unique characteristics comparing to various other probabilistic models in the sense that they capture the global repulsive behavior between items, give polynomial-time algorithms for statistical inference, and have geometrical intuition. Due to these advantages DPPs have played very important roles in machine learning, especially in subset selection problems, such as documentary summarization, image search, and pose determination [KT12], and so on. These real world applications necessitate the estimation of parameters of determinantal point process models. In this context, maximum likelihood estimator $\tilde{L}_{n}$ of the true parameter $L^{\star}$ is a natural choice, which in general leads to a non-convex optimization problem in our situation. Along this direction, Kulesza and Taskar split DPPs model into diversity part and quality part and only learn the quality part while the first part is fixed. They conjecture that the problem of learning the likelihood of DPPs is NP-hard, which has been proven by [GJWX22] a decade later. Brunel, Moitra, Rigollet, and Urschel [BMRU17] first study the local geometry of the expected maximum likelihood estimation of DPPs, that is, the curvature of likelihood function around its maximum. Then they prove that the maximum likelihood estimator converges to true values in probability and establish the corresponding central limit theorem. Motivated by these works, the present paper contains three further contributions.

(i)

For the consistency we improve the convergence in probability of $\tilde{L}_{n}$ to $L^{\star}$ obtained in [BMRU17] to almost sure convergence. This almost sure convergence property is important in applications since in reality we only observe a single sample path. Theoretically, it is well-known that convergence almost surely is always a hard problem. It is also the case here. We get around this difficulty by making use of the Wald’s consistency theorem in our situation.
(ii)

The asymptotic normality of the maximum estimator is also obtained in [BMRU17]. This means that $\sqrt{n}(\tilde{L}_{n}-L^{\star})$ is approximately a normal random variable. We study the rate of convergence associated this normality approximation. Namely, we obtain a Berry-Esseen type bound for the maximum likelihood estimator $\tilde{L}_{n}$ .
(iii)

Due to non-identifiability of DPPs observed in [Kul12], there is an involvement of $\hat{D}_{n}$ in the definition of $\tilde{L}_{n}$ . In general case this involvement of $\hat{D}_{n}$ seems unavoidable. To get rid of this $\hat{D}_{n}$ we consider special cases of blocked two by two matrix ensembles, where we obtain explicit form of the maximum likelihood estimator, and consequently, no involvement of $\hat{D}_{n}$ in these cases. This explicit form makes our computations much efficient. When we apply various numerical methods such as Newton-Raphson, Stochastic Gradient Descent method (SGD) to general maximum likelihood estimator we can hardly obtain satisfactory results.

The paper is organized as follows. In Section 2 we introduce some basic definitions and recall some properties of DPPs. We present a proof of a well-known result about the characterization of the determinantal process (Theorem 2.2 in next section), whose proof we cannot find explicitly in literature. In Section 3 we present our main results for the almost sure consistency and the Berry-Esseen type theorem for the maximum likelihood estimator. Newton-Raphson and stochastic gradient descent (SGD) algorithm are suggested to find the maximum likelihood estimator in general case in Section 4, and the results are reported although they are not satisfactory. In Section 5, we discuss the explicit maximum likelihood estimator for the two by two blocked ensembles and since the explicit form is available we get very satisfactory numerical results. Some concluding remarks are given in Section 6.

2 Preliminary

In this paper, we focus on the finite state point process. Let the ground set $\mathcal{Y}=[N]:=\{1,2,\cdots,N\}$ be a set of $N$ points. The set of all subsets of $\mathcal{Y}$ is denoted by $\mathbb{Y}$ . Let $(\Omega,\mathcal{F},\mathbb{P})$ be a complete probability space on which the expectation is denoted by $\mathbb{E}$ . A (finite) point process on a ground set $\mathcal{Y}$ is a random variable $\mathbf{Y}:\Omega\to\mathbb{Y}$ . We also call the probability measure $\mathcal{P}$ for $\mathbf{Y}$ the point process over the subsets of $\mathcal{Y}$ . Random subsets drawn from the point process $\mathcal{P}$ can be any subset between null set and full set $\mathcal{Y}$ .

Definition 2.1.

A $\mathbb{Y}$ -valued point process $\mathbf{Y}$ is called a determinantal point process if there is a symmetric and positive definite $N\times N$ symmetric matrix matrix $K$ such that for every fixed set $A\subseteq\mathcal{Y}$ ,

\mathbb{P}(A\subseteq\mathbf{Y})=\det(K_{A})\,,

(2.1)

where $K_{A}$ denotes the restriction of $K$ to the subset A, that is, $K_{A}\mathrel{\mathop{:}}=[K_{i,j}]_{i,j\in A}$ .

To our best knowledge the following sufficient and necessary condition for a symmetric kernel $K$ to define a determinant point process is often mentioned without directly proving it so we provide a proof as follows.

Theorem 2.2.

The sufficient and necessary condition for a symmetric kernel $K$ to define a determinant point process is $0\preceq K\preceq I$ .

Proof.

Since the marginal probability of empty set is the total probability space, $\mathbb{P}(\Omega)=\mathbb{P}(\emptyset\subseteq\mathbf{Y})=1$ . We set $\det(K_{\emptyset})=1$ . Because $\mathbb{P}$ is a probability measure, all principal minors of $K$ , i.e. $\det(K_{A})$ must be nonnegative, and thus K itself must be positive semidefinite, that is, $K\succeq 0$ .

Furthermore, from $\mathbb{P}(\emptyset=\mathbf{Y})+\mathbb{P}(\bigcup_{i=1}^{N}\{i\in\mathbf{Y}\})=1$ and using inclusion–exclusion principle we get

$\displaystyle\mathbb{P}(\bigcup_{i=1}^{N}\{i\in\mathbf{Y}\})$	$\displaystyle=$	$\displaystyle\sum_{i\in[N]}\mathbb{P}(i\in\mathbf{Y})-\sum_{\{i,j\}\subset[N]}\mathbb{P}(\{i,j\}\subseteq\mathbf{Y})+\cdots$	(2.2)
		$\displaystyle\qquad\cdots+(-1)^{N-1}\mathbb{P}([N]\subseteq\mathbf{Y})$
	$\displaystyle=$	$\displaystyle\sum_{\|A\|=1}\det(K_{A})-\sum_{\|A\|=2}\det(K_{A})+\cdots$
		$\displaystyle\qquad\cdots+(-1)^{N-1}\det(K)$
	$\displaystyle=$	$\displaystyle 1-\det(I-K)\,.$

The above last equality follows from the characteristic polynomial. Equation (2.2) also means

\mathbb{P}(\emptyset=\mathbf{Y})=\det(I-K)\geq 0.

(2.3)

Similarly, we are able to show that $\mathbb{P}(\emptyset=\mathbf{Y}\cap A)=\det(I_{A}-K_{A})\geq 0$ for any subset $A\subseteq[N]$ and hence $K\preceq I$ . So the necessary condition for a symmetric matrix to give a determinantal process is $0\preceq K\preceq I$ .

This condition turns out to be sufficient: any $0\preceq K\preceq I$ defines a DPP. To prove this, it’s sufficient to show that for every $A\subseteq[N]$ , the atomic probability is well-defined, that is, $0\leq\mathbb{P}(A=\mathbf{Y})\leq 1$ . The probability being less or equal to 1 holds since $K\preceq I$ . For the other inequality, we assume $K_{A}$ is invertible.¹¹1if $K_{A}$ is not invertible, we immediately get $\mathbb{P}(A=\mathbf{Y})=0$ . Then using Schur complement and characteristic polynomial, we have

$\displaystyle\mathbb{P}(A=\mathbf{Y})$	$\displaystyle=$	$\displaystyle\mathbb{P}(A\subseteq\mathbf{Y})-\mathbb{P}(\bigcup_{i\in\bar{A}}\{A\cup\{i\}\subseteq\mathbf{Y}\})$	(2.4)
	$\displaystyle=$	$\displaystyle\det(K_{A})-\sum_{i\in\bar{A}}\det(K_{A\cup\{i\}})+\sum_{\{i,j\}\subseteq\bar{A}}\det(K_{A\cup\{i,j\}})+$
		$\displaystyle\qquad\cdots+(-1)^{\|\bar{A}\|}\det(K)$
	$\displaystyle=$	$\displaystyle\det(K_{A})-\sum_{i\in\bar{A}}\det(K_{A})\det(K_{ii}-K_{\{i\},A}K^{-1}_{A}K_{A,\{i\}})$
		$\displaystyle\qquad+\sum_{\{i,j\}\subseteq\bar{A}}\det(K_{A})\det(K_{\{i,j\}}-K_{\{i,j\},A}K^{-1}_{A}K_{A,\{i,j\}})+$
		$\displaystyle\qquad\cdots+(-1)^{\|\bar{A\|}}\det(K_{A})\det(K_{\bar{A}}-K_{\bar{A},{A}}K^{-1}_{A}K_{A,\bar{A}})$
	$\displaystyle=$	$\displaystyle(-1)^{\|\bar{A}\|}\det(K_{A})\det((K_{\bar{A}}-K_{\bar{A},{A}}K^{-1}_{A}K_{A,\bar{A}})-I_{\bar{A}})$
	$\displaystyle=$	$\displaystyle(-1)^{\|\bar{A}\|}\det(K-I_{\bar{A}}),$

where $K_{A,B}$ denotes the matrix obtained from $K$ by keeping only those entries whose rows belong to $A$ and whose columns belong to $B$ (if $A=B$ we simply have $K_{A}$ .), $|A|$ denotes the cardinality of subset $A$ , and $\bar{A}$ the complement of set $A$ . Here we use a slight abuse of notation of $I_{\bar{A}}$ . We refer it to an N × N matrix whose restriction to $\bar{A}$ is $I_{\bar{A}}$ and has zeros everywhere else. Since $0\preceq K\preceq I$ , $\mathbb{P}(A=\mathbf{Y})=|\det(K-I_{\bar{A}})|\geq 0$ . This completes the proof of the theorem. ∎

It is quite inconvenient to work with marginal kernels since their eigenvalues should be bounded by 0 and 1, and the marginal probability is not very appropriate to describe real world data. Here we introduce a slightly smaller class of DPPs called L-ensembles, where L only needs to be positive semi-definite.

Definition 2.3.

A point process is called an L-ensemble if it is defined through a real, symmetric, positive semi-definite $N\times N$ matrix $L$ :

\mathbb{P}_{L}(A=\mathbf{Y})\propto\det(L_{A}),

(2.5)

where $A\subseteq\mathcal{Y}$ is a fixed subset.

The proportional coefficient is given by [KT12]:

\frac{1}{\sum_{A\subseteq\mathcal{Y}}\det(L_{A})}=\frac{1}{\det{(L+I)}}.

Moreover, [Mac75] shows that L-ensembles are indeed DPPs, where marginal kernels K is given by $L(L+I)^{-1}=I-(L+I)^{-1}.$

3 Maximum likelihood Estimator of DPPs

Let $\mathcal{S}_{[N]}^{++}$ denote the set of all $N$ by $N$ positive definite matrices. Suppose that $Z_{1},...,Z_{n}$ is a sample of size $n$ from a ${\text{DPP}(L^{\star})}$ for some given but unknown $L^{\star}\in\mathcal{S}_{[N]}^{++}$ . This means that $Z_{1},\cdots,Z_{n}$ are independent $\mathbb{Y}$ valued random subsets whose probability law is given by the ensemble $L^{\star}$ . The goal of this paper is to estimate this $L^{\star}$ from this given sample $Z_{1},\cdots,Z_{n}$ . The most popular one is the maximum likelihood estimator $\hat{N}_{n}$ . We shall recall its construction this our specific setting and study its strong consistency, asymptotic normality, and associated Berry-Esseen theorem.

Let $P_{L}(Z_{i})$ denote the probability of the L-ensemble being $Z_{i}$ . The scaled log-likelihood function associated with the sample $Z_{1},\cdots,Z_{n}$ is

\hat{\Phi}_{n}(L)=\frac{1}{n}\sum_{i=1}^{n}\log P_{L}(Z_{i})=\sum_{J\subseteq[N]}\hat{p}{(J)}\log\det(L_{J})-\log\det(L+I),

(3.1)

where $L\in\mathcal{S}_{[N]}^{++}$ and

\hat{p}{(J)}=\frac{1}{n}\sum_{i=1}^{n}\mathbb{I}(Z_{i}=J),

with $\mathbb{I}(\cdot)$ standing for the indicator function. It is also useful to define the expected log-maximum likelihood function given the real kernel $L^{\star}$

\Phi_{L^{\star}}(L)=\sum_{J\subseteq[N]}{p_{L^{\star}}(J)}\log\det(L_{J})-\log\det(L+I)\,,

(3.2)

where

p_{L^{\star}}(J)=\mathbb{E}{(\hat{p}{(J)})}=\frac{\det{(L^{\star}_{J})}}{\det{(L^{\star}+I)}}.

In the sequel $L^{\star}$ is always fixed. To simplify writing we let $\hat{p}_{J}$ denote $\hat{p}(J)$ , $p^{*}_{J}$ denote $p_{L^{*}}(J)$ and $\Phi$ denote $\Phi_{L^{\star}}$ .

Definition 3.1.

A maximum likelihood estimator $\hat{L}_{N}$ is a global maximizer of (3.1):

\hat{L}_{n}=\arg\!\!\max_{L\in\mathcal{S}_{[N]}^{++}}\hat{\Phi}_{n}(L)\,.

(3.3)

We notice that $\Phi{(L^{\star})}-\Phi{(L)}$ is the Kullback-Leibler divergence $\mathrm{KL}\big{(}\mathrm{DPP}(L^{\star}),\mathrm{DPP}(L)\big{)}$ , which measures the difference between distributions of $\mathrm{DPP}(L^{\star})$ and of $\mathrm{DPP}(L)$ . KL divergence is always non-negative:

\mathrm{KL}\big{(}\mathrm{DPP}(L^{\star}),\mathrm{DPP}(L)\big{)}=\Phi{(L^{\star})}-\Phi(L)\geq 0,\,\forall L\in\mathcal{S}_{[N]}^{++}.

As a consequence $L^{\star}$ is the global maxima of the expected maximum function $\Phi(L)$ .

The next Lemma gives the differential of the likelihood function $\hat{\Phi}_{n}(L)$ .

Lemma 3.2.

The gradient of log-likelihood function $\hat{\Phi}_{n}(L)$ defined in (3.1) exists and is given by

\mathop{}\!\mathrm{d}\hat{\Phi}_{n}(L)=\sum_{J\subseteq[N]}\hat{p}_{J}L_{J}^{-1}-(L+I)^{-1}.

(3.4)

Proof.

We regard determinant as a multivariate function of $N\times N$ variables and then the directional derivative of $\det(L+I)$ along direction $H$ is given by

$\displaystyle\mathop{}\!\mathrm{d}\det(L+I)(H)=$	$\displaystyle\lim_{t\to 0}\frac{\det(L+I+tH)-\det(L+I)}{t}$
$\displaystyle=$	$\displaystyle\lim_{t\to 0}\det(L+I)\Big{[}\frac{\det(I+t(L+I)^{-1}H)-1}{t}\Big{]}$
$\displaystyle=$	$\displaystyle\lim_{t\to 0}\det(L+I)\Big{[}\frac{1+t\operatorname{Tr}((L+I)^{-1}H)+o(t^{2})-1}{t}\Big{]}$
$\displaystyle=$	$\displaystyle\det(L+I)\operatorname{Tr}((L+I)^{-1}H),$	(3.5)

where the third equality follows from the power series representation of $\det(I+A)$ . Then the directional derivative of $\hat{\Phi}_{n}(L)$ along direction $H$ is

\mathop{}\!\mathrm{d}\hat{\Phi}_{n}(L)(H)=\sum_{J\subseteq[N]}\hat{p}_{J}\operatorname{Tr}(L_{J}^{-1}H_{J})-\operatorname{Tr}((L+I)^{-1}H).

(3.6)

In matrix form, the above equation becomes

\mathop{}\!\mathrm{d}\hat{\Phi}_{n}(L)=\sum_{J\subseteq[N]}\hat{p}_{J}L_{J}^{-1}-(L+I)^{-1}.

(3.7)

This proves the lemma ∎

3.1 Strong consistency

One critical issue for the maximum likelihood estimation is its consistency. Due to non-identifiability of DPPs illustrated in Theorem 4.1 in [Kul12] , any element in $\{DL^{\star}D:D\in\mathcal{D}\}$ defines the same DPP, where $\mathcal{D}$ is the collection of all diagonal matrices whose entry is either 1 or -1. Moreover, if $L_{1},L_{2}\in\mathcal{S}_{[N]}^{++}$ defines the same DPP, then there is a $D\in\mathcal{D}$ such that $L_{2}=DL_{1}D$ . Thus, we cannot expect that the maximum likelihood estimator $\hat{L}_{n}$ will converge to $L^{\star}$ . Instead, it should converge to $DL^{\star}D$ for some $D\in\mathcal{D}$ . Hence, for the consistency problem we study the limit

\ell(\hat{L}_{n},L^{\star}):=\min_{D\in\mathcal{D}}\|\hat{L}_{n}-DL^{\star}D\|_{F}.

(3.8)

[BMRU17] proves that this distance converges to zero in probability. We shall prove a stronger version: the convergence also holds almost surely.

Theorem 3.3.

Let $Z_{1},...,Z_{n}$ be a sample of size $n$ from $Z\sim{\text{DPP}(L^{\star})}$ . Let $\hat{L}_{n}$ be the maximum likelihood estimator of $L^{\star}$ defined by (3.3). Then $\ell(\hat{L}_{n},L^{\star})$ converges to zero almost surely.

Proof.

The proof is to combine the convergence in probability result of [BMRU17, Theorem 14] with Wald’s consistency theorem [Wal49] and is divided into the following five steps.

Step 1: For $0<\alpha<\beta<1$ , define a set $E_{\alpha,\beta}$

E_{\alpha,\beta}=\big{\{}L\in\mathcal{S}^{++}_{[N]}:K=L(I+L)^{-1}\in\mathcal{S}^{[\alpha,\beta]}_{[N]}\big{\}},

where $\mathcal{S}^{[\alpha,\beta]}_{[N]}$ is the set of all N by N symmetric matrices whose eigenvalues are between $\alpha$ and $\beta$ . Observe that $E_{\alpha,\beta}$ is compact since it’s bounded and closed in $\mathbb{R}^{N\times N}$ . We first show that $\ell(\hat{L}_{n},L^{\star})$ converges to zero almost surely when parameters of matrices are restricted to $L^{\star}\in E_{\alpha,\beta}$ for appropriate $0<\alpha<\beta<1$ . In order to complete this proof, we let

\Delta\hat{\Phi}(L)=\hat{\Phi}(L)-\hat{\Phi}(L^{\star})=\frac{1}{n}\sum_{i=1}^{n}\log\frac{P_{L}(Z_{i})}{P_{L^{\star}}(Z_{i})}

and

\Delta{\Phi}(L)={\Phi}(L)-{\Phi}(L^{\star})=\mathbb{E}_{L^{\star}}\big{(}\log\frac{P_{L}(Z)}{P_{L^{\star}}(Z)}\big{)}.

$\Delta{\Phi}(L)$ is the Kullback-Leibler divergence between DPP( $L^{\star}$ ) and DPP( $L$ ). By Jensen’s inequality, $\Delta{\Phi}(L)\leq 0$ for all $L$ and ${\Phi}(L)={\Phi}(L^{\star})$ if and only if $P_{L}(Z)=P_{L^{\star}}(Z)$ for all $Z\in[N]$ , which means $L=DL^{\star}D$ for some $D\in\mathcal{D}$ . In the sequel let $E$ denote $E_{L^{*}}$

For each $L\in E_{\alpha,\beta}$ , the strong law of large numbers implies

\Delta{\hat{\Phi}}(L)\xrightarrow[]{a.s.}\Delta{\Phi}(L).

(3.9)

Step 2: However, the above convergence (3.9) doesn’t imply the convergence of maximum likelihood estimator to the true values. Thus the Wald’s integrability condition is needed: for every $L\in E_{\alpha,\beta}$ , there exists $\epsilon>0$ such that,

\mathbb{E}\sup_{\begin{subarray}{c}N\in E_{\alpha,\beta}\\ \ell(L,N)<\epsilon\end{subarray}}\log\frac{P_{N}(Z)}{P_{L^{\star}}(Z)}<\infty.

(3.10)

Since $L\mapsto\log\frac{P_{L}(Z)}{P_{L^{\star}}(Z)}$ is continuous (the determinant function is continuous), for any $\delta>0$ there exists $\epsilon>0$ , when $\ell(L,N)<\epsilon$

(1-\delta)\frac{P_{L}(Z)}{P_{L^{\star}}(Z)}<\frac{P_{N}(Z)}{P_{L^{\star}}(Z)}<(1+\delta)\frac{P_{L}(Z)}{P_{L^{\star}}(Z)}.

Then the Wald’s integrability condition is satisfied. Now for every sequence $\{L_{n}\}$ converging to L, we show that $\Delta\Phi(L_{n})$ is upper semicontinuous:

$\displaystyle\limsup_{n\to\infty}\Delta\Phi(L_{n})$	$\displaystyle=$	$\displaystyle\limsup_{n\to\infty}\mathbb{E}\log\frac{P_{L_{n}}(Z)}{P_{L^{\star}}(Z)}$
	$\displaystyle\leq$	$\displaystyle\mathbb{E}\limsup_{n\to\infty}\log\frac{P_{L_{n}}(Z)}{P_{L^{\star}}(Z)}$
	$\displaystyle=$	$\displaystyle\mathbb{E}\log\frac{P_{L}(Z)}{P_{L^{\star}}(Z)}$
	$\displaystyle=$	$\displaystyle\Delta\Phi(L).$

The second inequality follows from the Fatou’s lemma and the third identity is the consequence of continuity of the function $\log\frac{P_{L_{n}}(Z)}{P_{L^{\star}}(Z)}$ . For every $\eta>0$ we define the set $K_{\eta}$

	$\displaystyle K_{\eta}$	$\displaystyle=\big{\{}L\in E_{\alpha,\beta}:\ell(L,L^{\star})\geq\eta\big{\}}$
		$\displaystyle=\bigcap_{D\in\mathcal{D}}\big{\{}L\in E_{\alpha,\beta}:\\|L-DL^{\star}D\\|_{F}\geq\eta\big{\}},$		(3.11)

which is closed and hence compact.

Since $\Delta\Phi(L)$ is an upper semicontinuous function, it achieves maximum over the compact set $K_{\eta}$ (since $E_{\alpha,\beta}$ is compact). We denote the maximum by $m(\eta)$ . And we cannot have $m(\eta)=0$ because that would imply there is an $L\in K_{\eta}$ such that $L=DL^{\star}D$ for some $D\in\mathcal{D}$ . The strong law of large numbers implies

	$\displaystyle\sup_{\begin{subarray}{c}N\in E_{\alpha,\beta}\\ \ell(L,N)<\epsilon\end{subarray}}\Delta\hat{\Phi}(N)$	$\displaystyle\leq$	$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\sup_{\begin{subarray}{c}N\in E_{\alpha,\beta}\\ \ell(L,N)<\epsilon\end{subarray}}\log\frac{P_{N}(Z_{i})}{P_{L^{\star}}(Z_{i})}$		(3.12)
		$\displaystyle\xrightarrow[]{a.s.}$	$\displaystyle\mathbb{E}\sup_{\begin{subarray}{c}N\in E_{\alpha,\beta}\\ \ell(L,N)<\epsilon\end{subarray}}\log\frac{P_{N}(Z)}{P_{L^{\star}}(Z)}.$		(3.12)

By continuity,

\lim_{\epsilon\to 0}\sup_{\begin{subarray}{c}N\in E_{\alpha,\beta}\\ \ell(L,N)<\epsilon\end{subarray}}\log\frac{P_{N}(Z))}{P_{L^{\star}}(Z)}=\log\frac{P_{L}(Z)}{P_{L^{\star}}(Z)}

and $\sup_{\epsilon}\log\frac{P_{N}}{P_{L^{\star}}}$ is a decreasing function with respect to $\epsilon$ because supremum over a smaller subset is smaller than that over a bigger subset. And by (3.10) it is integrable for all small enough $\epsilon$ . Hence by the dominated convergence theorem,

\lim_{\epsilon\to 0}\mathbb{E}\sup_{\begin{subarray}{c}N\in E_{\alpha,\beta}\\ \ell(L,N)<\epsilon\end{subarray}}\log\frac{P_{N}(Z))}{P_{L^{\star}}(Z)}=\mathbb{E}\log\frac{P_{L}(Z)}{P_{L^{\star}}(Z)}=\Delta\Phi(L).

Thus for any $L\in K_{\eta}$ and any $\gamma>0$ there exists an $\epsilon_{L}$ such that

\mathbb{E}\sup_{\begin{subarray}{c}N\in E_{\alpha,\beta}\\ \ell(L,N)<\epsilon_{L}\end{subarray}}\log\frac{P_{N}(Z)}{P_{L^{\star}}(Z)}<m(\eta)+\gamma.

(3.13)

Step 3. For each $L\in K_{\eta}$ , we define the open set:

V_{L}=\{N\in E_{\alpha,\beta}:\ell(N,L)<\epsilon_{L}\}

and then the family $\{V_{L}:L\in K_{\eta}\}$ is an open cover of $K_{\eta}$ and hence has a finite subcover: $V_{L_{1}},V_{L_{2}},....,V_{L_{d}}$ . On every $V_{L_{i}}$ we use strong law of large numbers again to obtain

	$\displaystyle\limsup_{n\to\infty}\sup_{N\in V_{L_{i}}}\Delta\hat{\Phi}(N)$	$\displaystyle\leq$	$\displaystyle\limsup_{n\to\infty}\frac{1}{n}\sum_{i=1}^{n}\sup_{N\in V_{L_{i}}}\log\frac{P_{N}(Z_{i})}{P_{L^{\star}}(Z_{i})}$		(3.14)
		$\displaystyle=$	$\displaystyle\mathbb{E}\sup_{N\in V_{L_{i}}}\log\frac{P_{N}(Z)}{P_{L^{\star}}(Z)}.$		(3.14)

From (3.13) we get

\limsup_{n\to\infty}\sup_{N\in V_{L_{i}}}\Delta\hat{\Phi}(N)<m(\eta)+\gamma\qquad i=1,2,...,d.

Since $\{V_{L_{i}}:i=1,2...,d\}$ cover $K_{\eta}$ we have

\limsup_{n\to\infty}\sup_{N\in K_{\eta}}\Delta\hat{\Phi}(N)<m(\eta)+\gamma

which, since $\gamma$ is arbitrary, implies

\limsup_{n\to\infty}\sup_{L\in K_{\eta}}\Delta\hat{\Phi}(L)<\sup_{L\in K_{\eta}}\Delta{\Phi}(L)=m(\eta).

(3.15)

Notice that $m(\eta)<0$ . From (3.15) there exists a constant $N_{1}$ such that

\sup_{L\in K_{\eta}}\Delta\hat{\Phi}(L)<\frac{m(\eta)}{2},\qquad n>N_{1}.

But

\Delta\hat{\Phi}(\hat{L}_{n})=\sup_{L\in E_{\alpha,\beta}}\Delta\hat{\Phi}({L})\geq\Delta\hat{\Phi}(L^{\star})\xrightarrow[]{a.s.}\Delta{\Phi}(L^{\star})=0,

so there exists a constant $N_{2}$ such that

\Delta\hat{\Phi}(\hat{L}_{n})\geq\frac{m(\eta)}{2},\qquad n>N_{2}

which implies that $\hat{L}_{n}\notin K_{\eta}$ , that is, $\ell(\hat{L}_{n},L)<\epsilon$ .

Step 4. Now we can remove the compactness condition. First we show that the event $\{\hat{L}_{n}\in E_{\alpha,\beta}\}$ holds almost surely. We adopt the proof from [BMRU17]. Let $\delta=\min_{J\subset[N]}P_{L^{\star}}{(J)}$ . For simplicity, we denote $P_{L^{\star}}{(J)}$ by $p^{\star}_{J}$ . Since $L^{\star}$ is positive definite, $\delta>0$ . Define the event $\mathcal{A}_{n}$ by

\mathcal{A}_{n}=\bigcap_{J\subset[N]}\big{\{}p^{\star}_{J}\leq 2\hat{p}_{J}\leq 3p^{\star}_{J}\big{\}}.

Observe that $\Phi(L^{\star})<0$ and we can find $\alpha<\exp(3\Phi(L^{\star})/\delta)$ and $\beta>1-\exp(3\Phi(L^{\star})/\delta)$ such that $0<\alpha<\beta<1$ . Then using [BMRU17, Theorem 14] we know that on the event $\mathcal{A}_{n}$ , $\hat{L}_{n}\in E_{\alpha,\beta}$ , that is,

P(\hat{L}_{n}\in E_{\alpha,\beta})\geq P(\mathcal{A}_{n}).

Because

\hat{p}_{J}=\frac{1}{n}\sum_{i=1}^{n}\mathbb{I}(Z_{i}=J)\xrightarrow[]{a.s.}P_{L^{\star}}(Z=J)=p^{\star}_{J}\,,

the event $\mathcal{A}_{n}$ holds almost surely when $n$ goes to infinity and hence $\{\hat{L}_{n}\in E_{\alpha,\beta}\}$ holds almost surely.

Step 5. Let $\mathbb{I}_{E_{n}}$ denote the indicator function of the event $\{\hat{L}_{n}\in E_{\alpha,\beta}\}$ . Then

$\displaystyle\mathbb{P}\big{(}\lim_{n\to\infty}\ell(\hat{L}_{n},L^{\star})=0\big{)}$	$\displaystyle=$	$\displaystyle\mathbb{P}\big{(}\lim_{n\to\infty}\ell(\hat{L}_{n},L^{\star})=0,\lim_{n\to\infty}\mathbb{I}_{E_{n}}=1\big{)}$
		$\displaystyle+\mathbb{P}\big{(}\lim_{n\to\infty}\ell(\hat{L}_{n},L^{\star})=0,\lim_{n\to\infty}\mathbb{I}_{E_{n}}\neq 1\big{)}$
	$\displaystyle=$	$\displaystyle\mathbb{P}\big{(}\lim_{n\to\infty}\ell(\hat{L}_{n},L^{\star})=0,\lim_{n\to\infty}\mathbb{I}_{E_{n}}=1\big{)}$
	$\displaystyle=$	$\displaystyle\mathbb{P}\big{(}\lim_{n\to\infty}\ell(\hat{L}_{n},L^{\star})=0\big{\|}\lim_{n\to\infty}\mathbb{I}_{E_{n}}=1\big{)}\mathbb{P}\big{(}\lim_{n\to\infty}\mathbb{I}_{E_{n}}=1\big{)}$
	$\displaystyle=$	$\displaystyle\mathbb{P}\big{(}\lim_{n\to\infty}\ell(\hat{L}_{n},L^{\star})=0\big{\|}\lim_{n\to\infty}\mathbb{I}_{E_{n}}=1\big{)}$
	$\displaystyle=$	$\displaystyle 1.$

The last equality follows from the fact that $\hat{L}_{n}\in E_{\alpha,\beta}$ almost surely proved in Step 4. ∎

3.2 Berry-Esseen theorem

An $N$ by $N$ matrix $[A_{ij}]_{N\times N}$ can also be viewed as an $N\times N$ dimensional column vector: $(A_{11},A_{12},...,A_{1N},A_{21},...,A_{N1},...A_{NN})^{T}$ . With this identification the Frobenius norm of the matrix $||\cdot||_{F}$ is just the $\mathcal{L}^{2}$ norm for its corresponding column vector. In the following we shall regard the matrix as the corresponding column vector.

Because of non-identifiability of DPPs, the maximum likelihood estimators are not unique. We choose the estimator $\tilde{L}_{n}$ which is closest to the given true value $L^{*}$ . In fact, let ${\hat{L}_{n}}$ be one maximal likelihood estimator. Let $\hat{D}_{n}\in\mathcal{D}$ be such that

\lVert\hat{D}_{n}\hat{L}_{n}\hat{D}_{n}-L^{\star}\rVert_{F}=\min_{{{D}\in\mathcal{D}}}\lVert{D}\hat{L}_{n}{D}-L^{\star}\rVert_{F}

(3.16)

and set

\tilde{L}_{n}=\hat{D}_{n}\hat{L}_{n}\hat{D}_{n}\,.

(3.17)

Then the strong consistency of $\tilde{L}_{n}$ immediately follows from Theorem 3.3.

We impose irreducibility assumption introduced in [BMRU17] on $L^{\star}$ , of which the assumption is equivalent to $L^{\star}$ being a one block matrix, that is, it cannot be partitioned into two block matrices. And then according to [BMRU17, Theorem 8], $\mathrm{d}^{2}{\Phi}(L^{\star})$ is negative definite and hence invertible. Let $V(L^{\star})$ denote its inverse:

V(L^{\star}):=\left[\mathrm{d}^{2}{\Phi}(L^{\star})\right]^{-1}\,.

(3.18)

Here if we vectorize ${L}$ then $\mathrm{d}^{2}{\Phi}(L^{\star})$ is an $(N\times N)\times(N\times N)$ Hessian matrix. By [VdV00, Theorem 5.41],

	$\displaystyle\sqrt{n}(\tilde{L}_{n}-L^{\star})$	$\displaystyle=$	$\displaystyle-({\mathbb{E}{(\mathrm{d}^{2}\log{P_{L^{\star}}(Z)})}})^{-1}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\mathrm{d}(\log P_{L^{\star}}(Z_{i}))+o_{P}(1)$
		$\displaystyle=$	$\displaystyle-V(L^{\star})\frac{1}{\sqrt{n}}\sum_{i=1}^{n}((L^{\star}_{Z_{i}})^{-1}-(I+L^{\star})^{-1})+o_{P}(1).$		(3.19)

In particular, [VdV00, Theorem 5.41] states that the sequence $\sqrt{n}(\tilde{L}_{n}-L^{\star})$ is asymptotically normal with mean $\bm{0}$ and covariance matrix $-V(L^{\star})$ . Hence we get the following theorem for the asymptotic normality of the maximum estimator $\tilde{L}_{n}$ .

Theorem 3.4.

Let $L^{\star}$ be irreducible. Then, $\tilde{L}_{n}$ is asymptotically normal:

\sqrt{n}(\tilde{L}_{n}-L^{\star})\stackrel{{\scriptstyle d}}{{\rightarrow}}\mathcal{N}(\bm{0},-V(L^{\star})),\quad\hbox{ as \ $n\longrightarrow\infty$}

(3.20)

where $V(L^{\star})$ is given by (3.18) and the above convergence $\stackrel{{\scriptstyle d}}{{\rightarrow}}$ holds in distribution.

Having established the asymptotic normality of the maximum estimator $\tilde{L}_{n}$ , we want to take one step further. Namely, we want to find an upper error bound on the rate of convergence of the distribution of $(-V(L^{\star}))^{-\frac{1}{2}}\sqrt{n}(\tilde{L}_{n}-L^{\star})$ to standard multidimensional normal distribution $Z\sim\mathcal{N}(\bm{0},{I})$ . In other words, we are going to obtain a Berry-Esseen type theorem. We argue that when $\tilde{L}_{n}\in E_{\alpha,\beta}$ , the bound of the maximal error is of order $n^{-\frac{1}{4}}$ . This restriction is not a big deal. Indeed, since $\alpha$ and $\beta$ can be chosen arbitrarily close to $0$ and $1$ respectively so that $E_{\alpha,\beta}$ converges to $\mathcal{S}_{[N]}^{++}$ . What’s more, since from Theorem 3.3, $\hat{L}\in E_{\alpha,\beta}$ almost surely, $\hat{D}\hat{L}\hat{D}=\tilde{L}_{n}\in E_{\alpha,\beta}$ almost surely.

Theorem 3.5.

Let the maximum likelihood estimator $\tilde{L}_{n}$ be defined by (3.17) and also belong to $E_{\alpha,\beta}$ for some $0<\alpha<\beta<1$ and let Z be an $N\times N$ standard Gaussian matrix. Then for every $x\in\mathbb{R}^{N\times N}$ ,

|\mathbb{P}((-V(L^{\star}))^{-\frac{1}{2}}\sqrt{n}(\tilde{L}_{n}-L^{\star})<x)-\mathbb{P}(Z<x)|\leq C\frac{1}{\sqrt[4]{n}},

where C is a sufficiently large constant, which is irrelevant to $x$ , subject to $\alpha,\beta$ and proportional to $N^{2}$ . Here for two $N$ by $N$ matrices $A$ and $B$ , $A<B$ means the components of $A$ is less than the corresponding components of $B$ .

Proof.

We divide the proof into four steps.

Step 1. According to (3.19), $(-V(L^{\star}))^{-\frac{1}{2}}\sqrt{n}(\tilde{L}_{n}-L^{\star})$ can be decomposed into a sum

X_{n}=\sum_{i=1}^{n}\xi_{i}:=(-V(L^{\star}))^{\frac{1}{2}}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}((L^{\star}_{Z_{i}})^{-1}-(I+L^{\star})^{-1})

(3.21)

and a term $\rho_{n}=(-V(L^{\star}))^{-\frac{1}{2}}o_{P}(1)$ whose Frobenius norm converges to zero in probability.

	$\displaystyle\|\mathbb{P}(X_{n}+\rho_{n}<{x})-\mathbb{P}(Z<{x})\|$
	$\displaystyle=\|\mathbb{P}(X_{n}+\rho_{n}<{x},\\|\rho_{n}\\|_{F}\geq k_{n})+\mathbb{P}(X_{n}+\rho_{n}<{x},\\|\rho_{n}\\|_{F}<k_{n})-\mathbb{P}(Z<{x})\|$
	$\displaystyle\leq\mathbb{P}(\\|\rho_{n}\\|_{F}\geq k_{n})+\|\mathbb{P}(X_{n}+\rho_{n}<{x},\\|\rho_{n}\\|_{F}<k_{n})-\mathbb{P}(Z<{x})\|$
	$\displaystyle\leq\mathbb{P}(\\|\rho_{n}\\|_{F}\geq k_{n})$
	$\displaystyle\qquad+\|\mathbb{P}(X_{n}+{k_{n}}{\mathbbm{1}}<x,\\|\rho_{n}\\|_{F}<k_{n})-\mathbb{P}(Z<{x})\|$
	$\displaystyle\qquad+\|\mathbb{P}(X_{n}-{k_{n}}{\mathbbm{1}}<x,\\|\rho_{n}\\|_{F}<k_{n})-\mathbb{P}(Z<{x})\|$
	$\displaystyle=:I_{1}+I_{2}+I_{3}\,,$		(3.22)

where $\{k_{n}\}$ is an arbitrary sequence of positive real number and $\mathbbm{1}$ is the $N\times N$ matrix whose entries are all 1.

Step 2. The Estimation of ( $I_{1}$ ). We claim $\mathbb{P}(\|\rho_{n}\|\geq k_{n})\leq\frac{C_{4}}{\sqrt[4]{n}},$ where $k_{n}=n^{-\frac{1}{4}}$ and $C_{4}$ is a constant.

In fact, from the proof of [VdV00, Theorem 5.41], $\rho_{n}$ has the following expression

	$\displaystyle\rho_{n}=\sqrt{n}(-V(L^{\star}))^{\frac{1}{2}}\bigg{(}$	$\displaystyle\operatorname{\mathbf{d}}^{2}\hat{\Phi}_{n}(L^{\star})-\mathbb{E}(\operatorname{\mathbf{d}}^{2}\hat{\Phi}_{n}(L^{\star}))$
		$\displaystyle+\frac{1}{2}(\tilde{L}_{n}-L^{\star})^{T}\operatorname{\mathbf{d}}^{3}\hat{\Phi}_{n}(L_{n})\bigg{)}(\tilde{L}_{n}-L^{\star}),$		(3.23)

where $L_{n}$ is a point on the line segment between $\tilde{L}_{n}$ and $L^{\star}$ . To simplify notation, let $\theta$ denote

\bigg{(}\operatorname{\mathbf{d}}^{2}\hat{\Phi}_{n}(L^{\star})-\mathbb{E}(\operatorname{\mathbf{d}}^{2}\hat{\Phi}_{n}(L^{\star}))+\frac{1}{2}(\tilde{L}_{n}-L^{\star})^{T}\operatorname{\mathbf{d}}^{3}\hat{\Phi}_{n}(L_{n})\bigg{)}(\tilde{L}_{n}-L^{\star}).

Then

$\displaystyle\mathbb{E}{\\|\rho_{n}\\|_{F}}$	$\displaystyle=$	$\displaystyle\mathbb{E}\\|\sqrt{n}(-V(L^{\star}))^{\frac{1}{2}}\theta\\|_{F}$	(3.24)
	$\displaystyle=$	$\displaystyle\sqrt{n}\mathbb{E}\\|(-V(L^{\star}))^{\frac{1}{2}}\theta\\|_{F}$
	$\displaystyle\leq$	$\displaystyle\sqrt{n}\mathbb{E}\\|(-V(L^{\star}))^{\frac{1}{2}}\\|_{op}\\|\theta\\|_{2}$
	$\displaystyle=$	$\displaystyle\sqrt{n\cdot\Lambda_{max}(-V)}\cdot\mathbb{E}\\|\theta\\|_{2}.$

$\|\cdot\|_{op}$ denotes the operator norm induced by $\mathcal{L}^{2}$ norm and $\Lambda_{max}$ denotes the largest eigenvalue. For the first inequality, we regard $\theta$ as an $N\times N$ column vector and $(-V(L^{\star}))^{\frac{1}{2}}$ is an $(N\times N)\times(N\times N)$ matrix.

$\displaystyle\mathbb{E}\\|\phi\\|_{2}=$	$\displaystyle\mathbb{E}\big{\\|}\big{(}\operatorname{\mathbf{d}}^{2}\hat{\Phi}_{n}(L^{\star})-\mathbb{E}(\operatorname{\mathbf{d}}^{2}\hat{\Phi}_{n}(L^{\star}))+\frac{1}{2}(\tilde{L}_{n}-L^{\star})^{T}\operatorname{\mathbf{d}}^{3}\hat{\Phi}_{n}(L_{n})\big{)}(\tilde{L}_{n}-L^{\star})\big{\\|}_{2}$
$\displaystyle\leq$	$\displaystyle\mathbb{E}\big{\\|}\big{(}\operatorname{\mathbf{d}}^{2}\hat{\Phi}_{n}(L^{\star})-\mathbb{E}(\operatorname{\mathbf{d}}^{2}\hat{\Phi}_{n}(L^{\star}))\big{)}(\tilde{L}_{n}-L^{\star})\big{\\|}_{2}$	(3.25)
$\displaystyle+$	$\displaystyle\mathbb{E}\\|\frac{1}{2}(\tilde{L}_{n}-L^{\star})^{T}\operatorname{\mathbf{d}}^{3}\hat{\Phi}_{n}(L_{n})(\tilde{L}_{n}-L^{\star})\big{\\|}_{2}\,.$	(3.26)

Using Cauchy-Schwartz inequality to estimate (3.25) we see

	$\displaystyle I_{1}-1$	$\displaystyle\leq$	$\displaystyle\mathbb{E}^{\frac{1}{2}}\big{\\|}\operatorname{\mathbf{d}}^{2}\hat{\Phi}_{n}(L^{\star})-\mathbb{E}(\operatorname{\mathbf{d}}^{2}\hat{\Phi}_{n}(L^{\star}))\big{\\|}^{2}_{op}\mathbb{E}^{\frac{1}{2}}\\|\tilde{L}_{n}-L^{\star}\\|^{2}_{2}$		(3.27)
		$\displaystyle\leq$	$\displaystyle\frac{N^{2}}{\sqrt{n}}\max_{i,j}({L^{\star}}^{-1})_{ij}^{2}\mathbb{E}^{\frac{1}{2}}\\|\tilde{L}_{n}-L^{\star}\\|^{2}_{2}.$		(3.27)

Let $h(x)$ be a multivariate function:

	$\displaystyle h:\quad\mathbb{R}^{N\times N}\longrightarrow{}$	$\displaystyle\,\mathbb{R}$
	$\displaystyle(x_{1},x_{2},...,x_{NN})\longmapsto$	$\displaystyle\,x_{1}^{2}+x_{2}^{2}+\cdots+x_{NN}^{2}$

Then $h$ is a continuous function. What’s more almost surely $\tilde{L}_{n}\in E_{\alpha,\beta}$ , which is a compact and convex set. Using Theorem 3.4 and portmanteau lemma we have

\mathbb{E}\big{(}h(\sqrt{n}(\tilde{L}_{n}-L^{\star}))\big{)}=n\mathbb{E}\|\tilde{L}_{n}-L^{\star}\|^{2}_{F}\longrightarrow\mathbb{E}\|\tilde{Z}\|^{2}_{F},

(3.28)

where $\tilde{Z}\sim\mathcal{N}(\bm{0},-V(L^{\star}))$ . $\mathbb{E}\|\tilde{Z}\|^{2}_{F}$ is equal to $\mathbb{E}(\tilde{Z}_{11}^{2}+\cdots+\tilde{Z}_{1n}^{2}+\tilde{Z}_{21}^{2}+\cdots+\tilde{Z}_{nn}^{2})=\mathrm{Tr}(-V(L^{\star}))$ . Then there exists a constant $C_{1}$ subject to $\alpha,\beta$ such that

\mathbb{E}^{\frac{1}{2}}\|\tilde{L}_{n}-L^{\star}\|^{2}_{2}\leq C_{1}\frac{1}{\sqrt{n}}.

(3.29)

As a result,

\text{\eqref{rho1}}\leq C_{2}N^{2}\frac{1}{n}\,,

(3.30)

where $C_{2}$ is a suitable constant.

Next, we estimate the second part, that is (3.26):

\mathbb{E}\|\frac{1}{2}(\tilde{L}_{n}-L^{\star})^{T}\operatorname{\mathbf{d}}^{3}\hat{\Phi}_{n}(L_{n})\big{)}(\tilde{L}_{n}-L^{\star})\big{\|}_{2}.

Here $\operatorname{\mathbf{d}}^{3}\hat{\Phi}_{n}(L_{n})$ is an $N\times N$ dimensional column vector whose entries are $N\times N$ matrices. Since $\hat{\Phi}_{n}(L)$ is infinitely many time differentiable, $L_{n}$ is on the line segment between $\tilde{L}_{n}$ and $L^{\star}$ , and $E_{\alpha,\beta}$ is a convex and compact set, we conclude that every entry of $\operatorname{\mathbf{d}}^{3}\hat{\Phi}_{n}(L_{n})$ is bounded. Hence there exists a constant $C_{3}\geq 0$ subject to $\alpha$ and $\beta$ such that

	$\displaystyle\mathbb{E}\\|\frac{1}{2}(\tilde{L}_{n}-L^{\star})^{T}\operatorname{\mathbf{d}}^{3}\hat{\Phi}_{n}(L_{n})\big{)}(\tilde{L}_{n}-L^{\star})\big{\\|}_{2}\leq$	$\displaystyle C_{3}\mathbb{E}\\|\tilde{L}_{n}-L^{\star}\\|^{2}_{2}$
	$\displaystyle\leq$	$\displaystyle\frac{C_{1}^{2}C_{3}}{n}.$		(3.31)

Now let $k_{n}=n^{-\frac{1}{4}}$ . Using Chebyshev’s inequality we get:

\mathbb{P}(\|\rho_{n}\|_{F}\geq k_{n})\leq\frac{\mathbb{E}\|\rho_{n}\|_{F}}{k_{n}}=\frac{C_{4}}{\sqrt[4]{n}}

(3.32)

for a suitable constant $C_{4}$ .

Step 3. Our next goal is to estimate ( $I_{2}$ ) as follows. Let $k_{n}$ be $\frac{1}{\sqrt[4]{n}}$ . Then

I_{2}\leq\frac{C_{7}}{\sqrt[4]{n}}

for some constant $C_{7}$ .

Because

	$\displaystyle\mathbb{P}(X_{n}+k_{n}\mathbbm{1}<x)-\mathbb{P}(Z<x)$
	$\displaystyle\geq\mathbb{P}(X_{n}+{k_{n}}\mathbbm{1}<x,\\|\rho_{n}\\|_{F}<k_{n})-\mathbb{P}(Z<{x})$
	$\displaystyle=\big{(}\mathbb{P}(X_{n}+k_{n}\mathbbm{1}<x)-\mathbb{P}(X_{n}+{k_{n}}\mathbbm{1}<x,\\|\rho_{n}\\|_{F}\geq k_{n})\big{)}-\mathbb{P}(Z<{x})$
	$\displaystyle\geq\mathbb{P}(X_{n}+k_{n}\mathbbm{1}<x)-\mathbb{P}(\\|\rho_{n}\\|_{F}>k_{n})-\mathbb{P}(Z<{x}),$

we have

$\displaystyle I_{2}$	$\displaystyle\leq\|\mathbb{P}(X_{n}+k_{n}\mathbbm{1}<x)-\mathbb{P}(Z<x)\|$
	$\displaystyle\quad+\|\mathbb{P}(X_{n}+k_{n}\mathbbm{1}<x)-\mathbb{P}(\|\rho_{n}\\|_{F}\geq k_{n})-\mathbb{P}(Z<{x})\|$
	$\displaystyle\leq 2\|\mathbb{P}(X_{n}+k_{n}\mathbbm{1}<x)-\mathbb{P}(Z<x)\|+\mathbb{P}(\\|\rho_{n}\\|_{F}\geq k_{n})$
	$\displaystyle=2\|\mathbb{P}(X_{n}+k_{n}\mathbbm{1}<x)-\mathbb{P}(Z+k_{n}\mathbbm{1}<x)$
	$\displaystyle\quad+\mathbb{P}(Z+k_{n}\mathbbm{1}<x)-\mathbb{P}(Z<x)\|+\mathbb{P}(\\|\rho_{n}\\|_{F}\geq k_{n})$
	$\displaystyle\leq 2\|\mathbb{P}(X_{n}+k_{n}\mathbbm{1}<x)-P(Z+k_{n}\mathbbm{1}<x)\|$
	$\displaystyle\quad+2\|\mathbb{P}(Z+k_{n}\mathbbm{1}<x)-P(Z<x)\|+\mathbb{P}(\\|\rho_{n}\\|_{F}\geq k_{n})$
	$\displaystyle=I_{21}+I_{22}+I_{23}.$	(3.33)

By multidimensional Berry-Esseen theorem in [Ben05],

I_{21}\leq C_{5}\cdot\sqrt{N}\cdot n\cdot\mathbb{E}{\|\xi_{1}\|_{2}^{3}}

(3.34)

where $C_{5}$ is a constant and $\xi_{1}$ is defined in (3.21):

	$\displaystyle\mathbb{E}{\\|\xi_{1}\\|}^{3}$	$\displaystyle=$	$\displaystyle\mathbb{E}\\|\frac{1}{\sqrt{n}}(-V(L^{\star}))^{-\frac{1}{2}}\big{(}(L^{\star}_{Z_{i}})^{-1}-(I+L^{\star})^{-1}\big{)}\\|_{2}^{3}$		(3.35)
		$\displaystyle\leq$	$\displaystyle(\frac{1}{\sqrt{n}})^{3}\mathbb{E}\\|(-V(L^{\star}))^{-\frac{1}{2}}\big{(}(L^{\star}_{Z_{i}})^{-1}-(I+L^{\star})^{-1}\big{)}\\|^{3}_{2}.$		(3.35)

Since $\mathbb{E}\|(-V(L^{\star}))^{-\frac{1}{2}}\big{(}(L^{\star}_{Z_{i}})^{-1}-(I+L^{\star})^{-1}\big{)}\|_{2}^{3}$ is a constant we get

I_{21}\leq C_{6}\sqrt{\frac{N}{n}}

(3.36)

For $I_{22}$ , since Z can be viewed as a standard Gaussian random vector, we have

$\displaystyle I_{22}$	$\displaystyle=$	$\displaystyle 2\|\mathbb{P}(x-k_{n}I<Z_{n}<x)\|$	(3.37)
	$\displaystyle\leq$	$\displaystyle 2\sum_{i,j=1}^{N}\mathbb{P}(x_{ij}-k_{n}\leq(Z_{n})_{ij}\leq x_{ij})$
	$\displaystyle=$	$\displaystyle\frac{2N^{2}}{\sqrt{2\pi}}k_{n}$

Combining (3.36), (3.37) with previous bound (3.32) to treat $I_{23}$ , where we take $k_{n}=n^{-\frac{1}{4}}$ we conclude that

I_{2}\leq\frac{C_{7}}{\sqrt[4]{n}},

where $C_{5}$ is a constant.

Step 4. As for $I_{3}$ we can use the same argument as above and conclude that $I_{3}$ is bounded by $C_{8}\cdot n^{-\frac{1}{4}}$ for some constant $C_{8}$ . This completes the proof of the theorem. ∎

4 Simulation study

In this section we test numerically the mle for the DPPs. Since in general there is no analytic expression for the mle (3.16) we need to use numerical algorithms to find the mle. We try the well-known Newton-Raphson algorithm and compare it with stochastic gradient descent algorithm.

Algorithms. First, we explain the formula of the Newton-Raphson algorithm applied to our mle of DPPs. Recall that the gradient of the maximum likelihood function is given by (3.4):

\mathop{}\!\mathrm{d}\hat{\Phi}_{n}(L)=\sum_{J\subseteq[N]}\hat{p}_{J}L_{J}^{-1}-(L+I)^{-1}.

To use Newton-Raphson algorithm we think of $\mathop{}\!\mathrm{d}\hat{\Phi}_{n}(L)$ as an $N\times N$ vector-valued function $(\mathop{}\!\mathrm{d}\hat{\Phi}_{n}(L)_{11},\mathop{}\!\mathrm{d}\hat{\Phi}_{n}(L)_{12},...,\mathop{}\!\mathrm{d}\hat{\Phi}_{n}(L)_{1N},\mathop{}\!\mathrm{d}\hat{\Phi}_{n}(L)_{21},...,\mathop{}\!\mathrm{d}\hat{\Phi}_{n}(L)_{N1},...\mathop{}\!\mathrm{d}\hat{\Phi}_{n}(L)_{NN})$ of $N\times N$ variables $(L_{11},L_{12},...,L_{1N},L_{21},...,L_{N1},...L_{NN})$ . Since

	$\displaystyle\frac{\mathop{}\!\mathrm{d}L^{-1}}{\mathop{}\!\mathrm{d}L_{ij}}$	$\displaystyle=-L^{-1}\cdot\frac{\mathop{}\!\mathrm{d}L}{\mathop{}\!\mathrm{d}L_{ij}}\cdot L^{-1}$
		$\displaystyle=-[L^{-1}_{ki}L^{-1}_{jl}]_{kl},$		(4.1)

the Jacobian matrix of $\hat{\Phi}_{n}(L)$ is given by an $N\times N$ by $N\times N$ matrix

	$\displaystyle\mathop{}\!\mathrm{d}^{2}\hat{\Phi}_{n}(L)=$	$\displaystyle-\sum_{J\subseteq[N]}\hat{p}_{J}\big{[}(L^{-1}_{ki}L^{-1}_{jl})\mathbb{I}\{i\in J,j\in J,k\in J,l\in J\}\big{]}_{ij,kl}$
		$\displaystyle+\big{[}(L+I)^{-1}_{ki}(L+I)^{-1}_{jl}\big{]}_{ij,kl},$		(4.2)

where the index $ij$ and $kl$ have order $(1,1),...(1,N),(2,1),...(2,N),...(N,1),...(N,N).$

Then Newton-Raphson algorithm updated at $(m+1)$ -th iteration is given by

L_{m+1}=L_{m}-(\mathop{}\!\mathrm{d}^{2}\hat{\Phi}_{n}(L_{m}))^{-1}\cdot\mathop{}\!\mathrm{d}\hat{\Phi}_{n}(L_{m}).

(4.3)

The matrix $\mathop{}\!\mathrm{d}^{2}\hat{\Phi}_{n}(L_{m})$ is $N^{2}\times N^{2}$ and the computation of its inverse may need some time when $N$ is large.

Newton-Raphson algorithm is time consuming in particular when $N$ is large since we need to compute the inverse of $\mathop{}\!\mathrm{d}^{2}\hat{\Phi}_{n}(L_{m})$ (which is an $N^{2}\times N^{2}$ matrix) in each iteration step. For this reason we try another method called stochastic gradient descent method. This method depends on a step size which we choose to be $\eta=0.1$ as well as initial kernel $L_{0}$ . At $(m+1)$ -th iteration we do simple random sampling in our data set to pick one sample $Z_{\rho}$ and update our kernel by

L_{m+1}=L_{m}+\eta\big{(}({L_{{m}_{Z_{\rho}}}})^{-1}-(L_{m}+I)^{-1}\big{)}.

(4.4)

In this way, we reduce high computational burden of Newton-Raphson algorithm (4.3) at each iteration but the algorithm has lower convergence rate.

Simulation results. In simulation study we consider two $3$ by $3$ matrix kernels $H_{1},H_{2}$ and one $2$ by $2$ matrix kernel $H_{3}$ . To see the trend of convergence, for each kernel we generate 300, 1000, 3000, 10000, and 30000 synthetic samples using the algorithm proposed in [HKPV06]. The simulation results show that Newton-Raphson method outperforms SGD. For the first DPP, Newton-Raphson algorithm shows a very slow trend of convergence to the true kernel $H_{1}$ whereas SGD only estimates the diagonal elements well. For the second independent DPP, both Newton-Raphson and SGD methods converge to the true kernel $H_{2}$ pretty well. For the last one, Newton-Raphson converges to a wrong critical point while SGD is numerically unstable.

Our simulation results of n = 30000 are summarized in Table 1, which suggests learning the dependence and the repulsive behavior between items seems beyond the scope of general numerical methods. This unsatisfactory result may be because on one hand both algorithms involve computing the inverse of matrices at each iteration, which could yield large values and hence may not be stable. On the other hand, the likelihood function of a DPP is usually non-concave and has at least exponentially many saddle points [BMRU17]. This unsatisfactory performance of the two well-known algorithms motivates us to find an effective method suitable to the maximum of a likelihood function of a DPP. In the next section, this issue is resolved for the two by two matrix kernel since we can obtain the explicit maximizer in this case.

True kernel	Initial kernel	${\textbf{Newton-Raphson}\atop k=100\ \ \hbox{iterations}}$	${\textbf{SGD}\atop k=60000\ \ \hbox{iterations}}$
$\begin{pmatrix}1&0.2&0\\ 0.2&2&0.3\\ 0&0.3&3\end{pmatrix}$	$\begin{pmatrix}1&0.1&0\\ 0.1&1&0.1\\ 0&0.1&1\end{pmatrix}$	$\begin{pmatrix}0.994&0.167&-0.048\\ 0.167&2.013&0.241\\ -0.048&0.241&3.029\end{pmatrix}$	$\begin{pmatrix}0.926&0.015&0.001\\ 0.015&2.323&0.025\\ 0.001&0.025&3.068\end{pmatrix}$
$\begin{pmatrix}7&0&0\\ 0&5&0\\ 0&0&9\end{pmatrix}$	$\begin{pmatrix}1&0&0\\ 0&1&0\\ 0&0&1\end{pmatrix}$	$\begin{pmatrix}7.009&0&0\\ 0&5.070&0\\ 0&0&9.0348\end{pmatrix}$	$\begin{pmatrix}7.002&0&0\\ 0&5.156&0\\ 0&0&9.035\end{pmatrix}$
$\begin{pmatrix}1&1\\ 1&2\\ \end{pmatrix}$	$\begin{pmatrix}0.5&0.1\\ 0.1&0.5\\ \end{pmatrix}$	$\begin{pmatrix}0.657&0\\ 0&1.499\\ \end{pmatrix}$	Numerically unstable

Table 1: We test the performance of algorithms for three different true kernels when the sample size is 30000. We obtain the result after a specific number

k

of iterations. In the experiment, we observe that both algorithms become unstable as iteration increases to a certain certain point so we choose

k=100

for Newton-Raphson and

k=60000

for SGD since for the later one each iteration needs less time to complete.

5 Two-by-two block kernel

Due to non-identifiability of DPPs illustrated in [Kul12, Theorem 4.1] that we mentioned earlier in Section 3, we have to use $\tilde{L}_{n}$ defined by (3.17) as our maximum likelihood estimator. This formula involves a choice of $\hat{D}_{n}$ which may also be random. However, the way to choose $\hat{D}_{n}$ is not available to us in general. To get rid of this choice of $\hat{D}_{n}$ we show in this section that if the kernels of determinantal point processes are two-by-two symmetric positive semi-definite matrices, we can obtain the desired results without involving the choice of $\hat{D}_{n}$ . Our result can also be immediately extended to any diagonally blocked two by two matrices. Since the commonly used iteration methods presented in the previous section are not so effective in finding the maximum likelihood estimator numerically this explicit analytic form we obtained is also very useful in compute the estimator numerically in practice. On the other hand, Our method effective to two by two matrices is difficult to apply to more general higher dimensional kernels.

Let $Z\sim\text{DPP}(L^{\star})$ , where $L^{\star}=\bigg{(}\begin{matrix}a^{*}&b^{*}\\ b^{*}&c^{*}\end{matrix}\bigg{)}$ , and the ground set be $\mathcal{Y}=[2]$ . To ensure $L^{\star}\in\mathcal{S}_{[2]}^{++}$ , we assume

a^{*},c^{*}>0

and

a^{*}c^{*}-{b^{*}}^{2}\geq 0.

We can always assume $b$ is non-negative since by identifiability of DPPs, $\bigg{(}\begin{matrix}a&b\\ b&c\end{matrix}\bigg{)}$ and $\bigg{(}\begin{matrix}a&-b\\ -b&c\end{matrix}\bigg{)}$ give the same DPP. For ease of notation, let $\hat{p}_{0},\hat{p}_{1},\hat{p}_{2},\hat{p}_{3}$ denote the empirical probability of the subset $\{\emptyset\}$ , $\{1\}$ , $\{2\}$ , $\{1,2\}$ respectively and let $p_{0},p_{1},p_{2},p_{3}$ denote the theoretical probability respectively. The relationship between $(a,b,c)$ and $(p_{0},p_{1},p_{2},p_{3})$ are given by

({a},{b},{c})=\Big{(}\frac{p_{1}}{p_{0}},\frac{\sqrt{p_{1}p_{2}-p_{0}p_{3}}}{p_{0}},\frac{p_{2}}{p_{0}}\Big{)},

and inversely

\begin{split}p_{0}=&\frac{1}{(a+1)(c+1)-b^{2}}\,,\qquad p_{1}=\frac{a}{(a+1)(c+1)-b^{2}}\,,\\ p_{2}=&\frac{c}{(a+1)(c+1)-b^{2}}\,,\qquad p_{3}=\frac{ac-b^{2}}{(a+1)(c+1)-b^{2}}\,.\end{split}

The likelihood function defined in (3.1) becomes now

	$\displaystyle\hat{\Phi}_{n}(L)$	$\displaystyle=\sum_{J\in[2]}\hat{p}_{J}\log(L_{J})-\log\det(L+I)$
		$\displaystyle=\hat{p}_{1}\log a+\hat{p}_{2}\log c+\hat{p}_{3}\log(ac-b^{2})-\log[(a+1)(c+1)-b^{2}]$		(5.1)

To find the critical point of (5) we first let the partial derivative of $\hat{\Phi}_{n}(L)$ with respect to $b$ equal zero and get

\displaystyle\frac{\partial\hat{\Phi}_{n}(L)}{\partial b}=

\displaystyle\frac{2\hat{p}_{3}b}{ac-b^{2}}+\frac{2b}{(a+1)(c+1)-b^{2}}=0.

(5.2)

Then we have $b$ is either equal to 0 or

b^{2}=\frac{ac-(a+1)(c+1)\hat{p}_{3}}{1-\hat{p}_{3}}.

(5.3)

If $b=0$ , then by setting the partial derivative with respect to $a$ and $c$ to zero and notice that $\hat{p}_{0}+\hat{p}_{1}+\hat{p}_{2}+\hat{p}_{3}=1$ we get the first critical point

(\hat{a},\hat{b},\hat{c})=\Bigg{(}\frac{\hat{p}_{1}+\hat{p}_{3}}{\hat{p}_{0}+\hat{p}_{2}},0,\frac{\hat{p}_{2}+\hat{p}_{3}}{\hat{p}_{0}+\hat{p}_{1}}\Bigg{)}.

(5.4)

We remark that this is the point which Newton-Raphson method converges to in the last simulation study. This critical point exists only if $\hat{p}_{0}+\hat{p}_{2}$ and $\hat{p}_{0}+\hat{p}_{1}$ is nonzero. Since empirical probability converges to its corresponding theoretical probability almost surely and $p_{0}>0$ , the strong law of large numbers implies the critical point exists almost surely when $n$ is sufficiently large.

If $b\not=0$ , then we can use (5.3) to estimate $\hat{b}$ once $\hat{a},\hat{c}$ are obtained:

\hat{b}=\sqrt{\frac{\hat{a}\hat{c}-(\hat{a}+1)(\hat{c}+1)\hat{p}_{3}}{1-\hat{p}_{3}}}.

(5.5)

To find the maximum likelihood estimators $\hat{a}$ and $\hat{c}$ of $a$ and $c$ we plug (5.3) into $\hat{\Phi}_{n}(L)$ to obtain

\hat{\Phi}_{n}(L)=\hat{p}_{1}\log a+\hat{p}_{2}\log c+(\hat{p}_{3}-1)\log(a+c+1)-(\hat{p}_{3}-1)\log\frac{\hat{p}_{3}}{1-\hat{p}_{3}}+\log\hat{p}_{3}.

(5.6)

Letting $\frac{\partial\hat{\Phi}_{n}(L)}{\partial a}$ and $\frac{\partial\hat{\Phi}_{n}(L)}{\partial c}$ equal zero yields

\displaystyle\begin{cases}\frac{\partial\hat{\Phi}_{n}(L)}{\partial a}=\frac{\hat{p}_{1}}{a}+\frac{\hat{p}_{3}-1}{a+c+1}=0\\ \frac{\partial\hat{\Phi}_{n}(L)}{\partial c}=\frac{\hat{p}_{2}}{c}+\frac{\hat{p}_{3}-1}{a+c+1}=0.\\ \end{cases}

(5.7)

The above system of function equations can be explicitly solved and combining it together with (5.5) yields

(\hat{a},\hat{b},\hat{c})=\Big{(}\frac{\hat{p}_{1}}{\hat{p}_{0}},\frac{\sqrt{\hat{p}_{1}\hat{p}_{2}-\hat{p}_{0}\hat{p}_{3}}}{\hat{p}_{0}},\frac{\hat{p}_{2}}{\hat{p}_{0}}\Big{)},

(5.8)

from which we have this critical point exists only if $\hat{p}_{0}>0$ and $\hat{p}_{1}\hat{p}_{2}-\hat{p}_{0}\hat{p}_{3}\geq 0$ . Again by strong laws of large numbers, the second critical point also exists and converges to the true value almost surely. In fact, we have almost surely,

\frac{\hat{p}_{1}}{\hat{p}_{0}}\to\frac{p_{1}}{p_{0}}=a^{*},\quad\frac{\sqrt{\hat{p}_{1}\hat{p}_{2}-\hat{p}_{0}\hat{p}_{3}}}{\hat{p}_{0}}\to\frac{\sqrt{p_{1}p_{2}-p_{0}p_{3}}}{p_{0}}=b^{*},\quad\frac{\hat{p}_{2}}{\hat{p}_{0}}\to c^{*}.

We test the estimator of the 2 by 2 kernel in earlier simulation study by taking $n=300,3000,10000,30000,$ and find the estimates respectively are

\begin{pmatrix}0.8598&0.8661\\ 0.8661&1.9159\\ \end{pmatrix},\begin{pmatrix}0.9454&0.9184\\ 0.9184&1.9868\\ \end{pmatrix},\begin{pmatrix}1.0247&1.0253\\ 1.0253&2.0241\\ \end{pmatrix},\begin{pmatrix}0.9942&1.0056\\ 1.0056&2.0015\\ \end{pmatrix},

which shows the estimator (5.8) is consistent numerically to the true kernel.

Furthermore, we can establish the central limit theorem for the estimator (5.8), which corresponds to the result in Theorem 3.4.

Theorem 5.1.

Assume $b>0$ , then the estimator $(\hat{a},\hat{b},\hat{c})$ in (5.8) is asymptotically normal,

\sqrt{n}((\hat{a},\hat{b},\hat{c})-(a^{*},b^{*},c^{*}))\xrightarrow[n\longrightarrow\infty]{}\mathcal{N}(\bm{0},-V({a^{*},b^{*},c^{*}})),

(5.9)

where the convergence holds in distribution and $V({a^{*},b^{*},c^{*}})$ is the inverse of the Hessian matrix of the expected maximum likelihood function $\Phi(a,b,c)=p_{1}\log a+p_{2}\log c+p_{3}\log(ac-b^{2})-\log[(a+1)(c+1)-b^{2}]$ .

Proof.

Let $Z_{1},...,Z_{n}$ be n independent subsets of $Z\sim\text{DPP}(L^{*})$ , where $L^{*}=\bigg{(}\begin{matrix}a^{*}&b^{*}\\ b^{*}&c^{*}\end{matrix}\bigg{)}$ . Let $X_{i}$ be the random vector $(\mathbb{I}_{\{Z_{i}=\emptyset\}},\mathbb{I}_{\{Z_{i}=\{1\}\}},\mathbb{I}_{\{Z_{i}=\{2\}\}},\mathbb{I}_{\{Z_{i}=\{1,2\}\}})^{T}$ , where $\mathbb{I}_{\{\cdot\}}$ stands for the indicator random variable. Then $X_{i}$ has mean $\bm{\mu}=(p_{0},p_{1},p_{2},p_{3})^{T}$ and covariance matrix

\mathbf{\Sigma}=\begin{pmatrix}p_{0}-p_{0}^{2}&-p_{0}p_{1}&-p_{0}p_{2}&-p_{0}p_{3}\\ -p_{0}p_{1}&p_{1}-p_{1}^{2}&-p_{1}p_{2}&-p_{1}p_{3}\\ -p_{0}p_{2}&-p_{1}p_{2}&p_{2}-p_{2}^{2}&-p_{2}p_{3}\\ -p_{0}p_{3}&-p_{1}p_{3}&-p_{2}p_{3}&p_{3}-p_{3}^{2}\end{pmatrix}\,.

By central limit theorem, $\sqrt{n}(\overline{X}_{n}-\bm{\mu})$ converges to a multivariate normal distribution with mean $\bm{0}$ and covariance $\bm{\Sigma}$ . Let a function $g:\mathbb{R}^{4}\to\mathbb{R}^{3}$ be defined by

g(x_{1},x_{2},x_{3},x_{4})=(\frac{x_{2}}{x_{1}},\frac{\sqrt{x_{2}x_{3}-x_{1}x_{4}}}{x_{1}},\frac{x_{3}}{x_{1}}).

Its Jacobi matrix $\dot{g}(\bm{x})=\big{[}\frac{\partial g_{i}}{\partial x_{j}}\big{]}_{3\times 4}$ is given by

\begin{pmatrix}-\frac{x_{2}}{x_{1}^{2}}&\frac{1}{x_{1}}&0&0\\ -\frac{x_{4}}{2x_{1}\sqrt{x_{2}x_{3}-x_{1}x_{4}}}-\frac{\sqrt{x_{2}x_{3}-x_{1}x_{4}}}{x_{1}^{2}}&\frac{x_{3}}{2x_{1}\sqrt{x_{2}x_{3}-x_{1}x_{4}}}&\frac{x_{2}}{2x_{1}\sqrt{x_{2}x_{3}-x_{1}x_{4}}}&-\frac{1}{2\sqrt{x_{2}x_{3}-x_{1}x_{4}}}\\ -\frac{x_{3}}{x_{1}^{2}}&0&\frac{1}{x_{1}}&0\end{pmatrix}.

Now we are in the position to apply Delta method [VdV00] to get

\sqrt{n}\big{(}(\hat{a},\hat{b},\hat{c})-(a^{*},b^{*},c^{*})\big{)}=\sqrt{n}\big{(}g(\overline{X}_{n})-g(\bm{\mu})\big{)}\xrightarrow[]{d}\mathcal{N}(\bm{0},\dot{g}(\bm{\mu})\bm{\Sigma}\dot{g}(\bm{\mu})^{\prime}).

After tedious matrix computations, $\dot{g}(\bm{\mu})\bm{\Sigma}\dot{g}(\bm{\mu})^{\prime}$ is found to be

D\begin{pmatrix}(a^{*}+{a^{*}}^{2})&\sigma_{12}&\sigma_{13}\\ \sigma_{12}&\frac{\frac{a^{*}c^{*}}{{b^{*}}^{2}}-1}{4}D+\frac{a^{*}+c^{*}+4a^{*}c^{*}}{4}&\sigma_{23}\\ \sigma_{13}&\sigma_{23}&c^{*}+{c^{*}}^{2}\end{pmatrix},

where

\left\{\begin{split}D=&(a^{*}+1)(c^{*}+1)-{b^{*}}^{2}\,;\\ \sigma_{12}=&(\frac{a^{*}c^{*}}{2b^{*}}+a^{*}b^{*}+\frac{a^{*}}{2b^{*}}(a^{*}c^{*}-{b^{*}}^{2}))\,;\\ \sigma_{13}=&a^{*}c^{*}\,;\\ \sigma_{23}=&\frac{a^{*}c^{*}}{2b^{*}}+b^{*}c^{*}+\frac{c^{*}}{2b^{*}}(a^{*}c^{*}-{b^{*}}^{2})\,.\\ \end{split}\right.

It is straightforward to verify the above matrix is the inverse of the Hessian matrix of the expected maximum likelihood function $\Phi(L)$ , that is, $-V(a^{*},b^{*},c^{*})$ , which in turn verifies Theorem 3.4 in this special case. However, in this two-by-two case, our maximum likelihood estimator is unique without the maneuver of the identifiability (3.17). ∎

This idea can be extended to blocked ensemble with the two by two block submatrices. If $L^{\star}$ is a matrix with k two-by-two blocks $J_{1},...,J_{k}$

\begin{pmatrix}J_{1}&&&\\ &J_{2}&&\\ &&\ddots&\\ &&&J_{k}\end{pmatrix}=\begin{pmatrix}a_{1}&b_{1}&&&&\\ b_{1}&c_{1}&&&&\\ &&a_{2}&b_{2}&&\\ &&b_{2}&c_{2}&&\\ &&&&\ddots\\ &&&&&a_{k}&b_{k}\\ &&&&&b_{k}&c_{k}\end{pmatrix},

(5.10)

where for each $1\leq i\leq k$ , $a_{i}$ , $b_{i}$ , $c_{i}>0$ and $a_{i}c_{i}-b_{i}^{2}\geq 0$ . Let ground set $\mathcal{Y}$ of this DPP be $\{J_{1}^{1},J_{1}^{2},J_{2}^{1},J_{2}^{2},...,J_{k}^{1},J_{k}^{2}\}$ and for each $1\leq i\leq k$ ,

$\displaystyle\hat{p}_{J_{i}}^{0}=$	$\displaystyle\frac{1}{n}\sum_{m=1}^{n}\mathbb{I}\{J_{i}^{1}\notin Z_{m},J_{i}^{2}\notin Z_{m}\}$	(5.11)
$\displaystyle\hat{p}_{J_{i}}^{1}=$	$\displaystyle\frac{1}{n}\sum_{m=1}^{n}\mathbb{I}\{J_{i}^{1}\in Z_{m},J_{i}^{2}\notin Z_{m}\}$	(5.12)
$\displaystyle\hat{p}_{J_{i}}^{2}=$	$\displaystyle\frac{1}{n}\sum_{m=1}^{n}\mathbb{I}\{J_{i}^{1}\notin Z_{m},J_{i}^{2}\in Z_{m}\}$	(5.13)
$\displaystyle\hat{p}_{J_{i}}^{3}=$	$\displaystyle\frac{1}{n}\sum_{m=1}^{n}\mathbb{I}\{J_{i}^{1}\in Z_{m},J_{i}^{2}\in Z_{m}\},$	(5.14)

where $Z_{1},...,Z_{n}$ are n independent subsets drawn from $\text{DPP}(L^{\star})$ . By the construction (5.10) we see that $Z\cap J_{1},...,Z\cap J_{k}$ are mutually independent. Then the result of critical point for two by two matrix can be applied:

(\hat{a}_{i},\hat{b}_{i},\hat{c}_{i})=\Bigg{(}\frac{\hat{p}_{J_{i}}^{1}}{\hat{p}_{J_{i}}^{0}},\frac{\sqrt{\hat{p}_{J_{i}}^{1}\hat{p}_{J_{i}}^{2}-\hat{p}_{J_{i}}^{0}\hat{p}_{J_{i}}^{3}}}{\hat{p}_{J_{i}}^{0}},\frac{\hat{p}_{J_{i}}^{2}}{\hat{p}_{J_{i}}^{0}}\Bigg{)},

(5.15)

for every $1\leq i\leq k$ .

However the above method encounters some difficulties when the kernel has dimension higher than 2 with general form. For example, if the kernel is a $3\times 3$ matrix

\begin{pmatrix}a&d&e\\ d&b&f\\ e&f&c\end{pmatrix},

the letting the gradient of likelihood function $\hat{\Phi}_{n}(L)$ equal zero will yield

\displaystyle\mathop{}\!\mathrm{d}\hat{\Phi}_{n}(L)=

\displaystyle\sum_{J\subseteq[3]}\hat{p}_{J}L_{J}^{-1}-(L+I)^{-1}=0\,.

Computing $L^{-1}$ and $(L+I)^{-1}$ could be troublesome. For example, $L^{-1}$ is:

\frac{1}{a(bc-f^{2})-d(cd-ef)+e(df-be)}\begin{pmatrix}bc-f^{2}&-cd+ef&-be+df\\ -cd+ef&ac-e^{2}&de-af\\ -be+df&de-af&ab-d^{2}\end{pmatrix}

which is difficult to use to obtain explicit maximum likelihood estimator.

6 Conclusion

In this paper, we study the maximum likelihood estimator for the ensemble matrix associated with determinantal process. Brunel et al show that the expected likelihood function $\Phi(L)$ is locally strongly concave around true value $L^{\star}$ if and only if $L^{\star}$ is irreducible, since the Hessian matrix of $\Phi{(L)}$ at $L^{\star}$ is negative definite. Then they prove the maximum likelihood estimator (MLE) is consistent with respect to the convergence in probability and when $L^{\star}$ is irreducible they also obtained the central limiting theorem for the MLE. We improve their results to show that the MLE is also strongly consistent with respect to the almost sure convergence. Moreover, we obtain the Berry-Esseen type result for the central limiting theorem and find the $n^{-\frac{1}{4}}$ rate of convergence of the MLE to normality. Last, we obtain the explicit form of the MLE where $L^{\star}$ is a two by two block matrix, or a block matrix whose blocks are two by two matrices. This explicit form enables us to avoid the use of $\hat{D}_{n}$ in (3.17) which is inconvenient in applications. The strong consistency and central limit theorem follows from these explicit forms, which demonstrates the general strong consistency and central limit theorem proved earlier. It would be interesting to find the explicit form of some particular higher dimensional DPPs. However, as the learning of maximum likelihood of DPPs is proven to be NP-hard, the explicit form for general ensembles, even if was found, would be very difficult to compute.

In addition to the maximum likelihood estimator there are also other approaches in lieu of MLE. Let us only mention one alternative approach. For all J such that $|J|\leq 1$ , we let

\frac{\det(L_{J})}{\det(L+I)}=\hat{p}_{J},

(6.1)

where the left hand side is the theoretical probability of $J$ and the right hand side is the empirical probability of $J$ . Taking $J=\{i\}$ suggests us the following estimator for $L_{ii}$ .

\hat{L}_{ii}=\frac{\hat{p}_{i}}{\hat{p}_{0}}.

(6.2)

Using equations (6.1) for $|J|=2$ again we are able to determine the off-diagonal elements up to the sign

L^{2}_{ij}=\frac{\hat{p}_{i}\hat{p}_{j}-\hat{p}_{\emptyset}\hat{p}_{\{i,j\}}}{\hat{p}^{2}_{\emptyset}},

(6.3)

where $i\neq j$ . Notice that this is the maximum likelihood estimator when $L$ is two dimensional. There is a question on how to choose the sign for $L_{ij}$ in (6.3), which has been resolved by [UBMR17] with graph theory. We shall not pursue this direction in current paper.

References

[BDF10] Alexei Borodin, Persi Diaconis, and Jason Fulman. On adding a list of numbers (and other one-dependent determinantal processes). Bulletin of the American Mathematical Society, 47(4):639–670, 2010.
[Ben05] Vidmantas Bentkus. A lyapunov-type bound in $R^{d}$ . Theory of Probability and Its Applications, 49(2):311–323, 2005.
[BMRU17] Victor-Emmanuel Brunel, Ankur Moitra, Philippe Rigollet, and John Urschel. Maximum likelihood estimation of determinantal point processes. arXiv preprint arXiv:1701.06501, 2017.
[BP93] Robert Burton and Robin Pemantle. Local characteristics, entropy and limit theorems for spanning trees and domino tilings via transfer-impedances. The Annals of Probability, pages 1329–1371, 1993.
[Gin65] Jean Ginibre. Statistical ensembles of complex, quaternion, and real matrices. Journal of Mathematical Physics, 6(3):440–449, 1965.
[GJWX22] Elena Grigorescu, Brendan Juba, Karl Wimmer, and Ning Xie. Hardness of maximum likelihood learning of DPPs. In Proceedings of Thirty Fifth Conference on Learning Theory, volume 178, pages 3800–3819. PMLR, 2022.
[HKPV06] J Ben Hough, Manjunath Krishnapur, Yuval Peres, and Bálint Virág. Determinantal processes and independence. Probability Survey, 3:206–229, 2006.
[KT12] Alex Kulesza and Ben Taskar. Determinantal point processes for machine learning. Foundations and Trends in Machine Learning, 5(2–3):123–286, 2012.
[Kul12] John A Kulesza. Learning with determinantal point processes. University of Pennsylvania, 2012.
[LDG21] Claire Launay, Agnès Desolneux, and Bruno Galerne. Determinantal point processes for image processing. SIAM J. Imaging Sci., 14(1):304–348, 2021.
[Mac75] O. Macchi. The coincidence approach to stochastic point processes. Advances in Applied Probability, 7(1):83–122, 1975.
[UBMR17] John Urschel, Victor-Emmanuel Brunel, Ankur Moitra, and Philippe Rigollet. Learning determinantal point processes with moments and cycles. In International Conference on Machine Learning, pages 3511–3520. PMLR, 2017.
[VdV00] Aad W Van der Vaart. Asymptotic statistics, volume 3. Cambridge university press, 2000.
[Wal49] Abraham Wald. Note on the consistency of the maximum likelihood estimate. The Annals of Mathematical Statistics, 20(4):595–601, 1949.

	$\displaystyle\|\mathbb{P}(X_{n}+\rho_{n}<{x})-\mathbb{P}(Z<{x})\|$
	$\displaystyle=\|\mathbb{P}(X_{n}+\rho_{n}<{x},\\|\rho_{n}\\|_{F}\geq k_{n})+\mathbb{P}(X_{n}+\rho_{n}<{x},\\|\rho_{n}\\|_{F}<k_{n})-\mathbb{P}(Z<{x})\|$
	$\displaystyle\leq\mathbb{P}(\\|\rho_{n}\\|_{F}\geq k_{n})+\|\mathbb{P}(X_{n}+\rho_{n}<{x},\\|\rho_{n}\\|_{F}<k_{n})-\mathbb{P}(Z<{x})\|$
	$\displaystyle\leq\mathbb{P}(\\|\rho_{n}\\|_{F}\geq k_{n})$
	$\displaystyle\qquad+\|\mathbb{P}(X_{n}+{k_{n}}{\mathbbm{1}}<x,\\|\rho_{n}\\|_{F}<k_{n})-\mathbb{P}(Z<{x})\|$
	$\displaystyle\qquad+\|\mathbb{P}(X_{n}-{k_{n}}{\mathbbm{1}}<x,\\|\rho_{n}\\|_{F}<k_{n})-\mathbb{P}(Z<{x})\|$
	$\displaystyle=:I_{1}+I_{2}+I_{3}\,,$		(3.22)

	$\displaystyle\mathbb{P}(X_{n}+k_{n}\mathbbm{1}<x)-\mathbb{P}(Z<x)$
	$\displaystyle\geq\mathbb{P}(X_{n}+{k_{n}}\mathbbm{1}<x,\\|\rho_{n}\\|_{F}<k_{n})-\mathbb{P}(Z<{x})$
	$\displaystyle=\big{(}\mathbb{P}(X_{n}+k_{n}\mathbbm{1}<x)-\mathbb{P}(X_{n}+{k_{n}}\mathbbm{1}<x,\\|\rho_{n}\\|_{F}\geq k_{n})\big{)}-\mathbb{P}(Z<{x})$
	$\displaystyle\geq\mathbb{P}(X_{n}+k_{n}\mathbbm{1}<x)-\mathbb{P}(\\|\rho_{n}\\|_{F}>k_{n})-\mathbb{P}(Z<{x}),$

$\displaystyle I_{2}$	$\displaystyle\leq\|\mathbb{P}(X_{n}+k_{n}\mathbbm{1}<x)-\mathbb{P}(Z<x)\|$
	$\displaystyle\quad+\|\mathbb{P}(X_{n}+k_{n}\mathbbm{1}<x)-\mathbb{P}(\|\rho_{n}\\|_{F}\geq k_{n})-\mathbb{P}(Z<{x})\|$
	$\displaystyle\leq 2\|\mathbb{P}(X_{n}+k_{n}\mathbbm{1}<x)-\mathbb{P}(Z<x)\|+\mathbb{P}(\\|\rho_{n}\\|_{F}\geq k_{n})$
	$\displaystyle=2\|\mathbb{P}(X_{n}+k_{n}\mathbbm{1}<x)-\mathbb{P}(Z+k_{n}\mathbbm{1}<x)$
	$\displaystyle\quad+\mathbb{P}(Z+k_{n}\mathbbm{1}<x)-\mathbb{P}(Z<x)\|+\mathbb{P}(\\|\rho_{n}\\|_{F}\geq k_{n})$
	$\displaystyle\leq 2\|\mathbb{P}(X_{n}+k_{n}\mathbbm{1}<x)-P(Z+k_{n}\mathbbm{1}<x)\|$
	$\displaystyle\quad+2\|\mathbb{P}(Z+k_{n}\mathbbm{1}<x)-P(Z<x)\|+\mathbb{P}(\\|\rho_{n}\\|_{F}\geq k_{n})$
	$\displaystyle=I_{21}+I_{22}+I_{23}.$	(3.33)