Differentially private analysis of networks with covariates via a generalized $\beta$-model

Let $G_{n}$ be an undirected graph on $n\geq 2$ nodes labeled by “$1,\ldots,n$”. Let $A=(a_{ij})_{n\times n}$ be the adjacency matrix of $G_{n}$, where $a_{ij}$ is an indicator denoting whether node $i$ is connected to node $j$. That is, $a_{ij}=1$ if there is a link between $i$ and $j$; $a_{ij}=0$ otherwise. We do not consider self-loops here, i.e., $a_{ii}=0$. Let $d_{i}=\sum_{j\neq i}a_{ij}$ be the degree of node $i$ and $d=(d_{1},\ldots,d_{n})^{\top}$ be the degree sequence of the graph $G_{n}$. We also observe a $p$-dimensional vector $z_{ij}$, the covariate information attached to the edge between nodes $i$ and $j$. The covariate $z_{ij}$ can be formed according to the similarity or dissimilarity between nodal attributes $z_{i}$ and $z_{j}$ for nodes $i$ and $j$. Specifically, $z_{ij}$ can be represented through a symmetric function $g(\cdot,\cdot)$ with $z_{i}$ and $z_{j}$ as its arguments. As an example if $z_{i}$ is an indicator of genders (e.g., $1$ for male and $-1$ for female), then we could use $z_{ij}=z_{i}\cdot z_{j}$ to denote the similarity or dissimilarity measurement between $i$ and $j$.

The $\beta$-model with covariates [Graham, (2017); Yan et al., (2019)] assumes that the edge $a_{ij}$ between $i$ and $j$ conditional on the unobserved degree effects and observed covariates has the following probability:

$$\mathbb{P}(a_{ij}=a|\beta,z_{ij})=\frac{e^{(\beta_{i}+\beta_{j}+z_{ij}^{\top}\gamma)a}}{1+e^{\beta_{i}+\beta_{j}+z_{ij}^{\top}\gamma}},~{}~{}a\in\{0,1\},$$

(1)

independent of other edges. The parameter $\beta_{i}$ is the intrinsic individual effect that reflects the node heterogeneity to participate in network connection. The common parameter $\gamma$ is exogenous, measuring the homophily or heterophilic effect. A larger homophily component $z_{ij}^{\top}\gamma$ means a larger homophily effect. We will refer to $\gamma$ as the homophily parameter hereafter although it could represent heterophilic measurement. Hereafter, we call model (1) the covariate-adjusted $\beta$-model.

The log-likelihood function is

$$\ell(\beta,\gamma)=\sum_{i=1}^{n}\beta_{i}d_{i}+\sum_{i<j}a_{ij}z_{ij}^{\top}\gamma-\sum_{i<j}\log\left(1+\frac{e^{\beta_{i}+\beta_{j}+z_{ij}^{\top}\gamma}}{1+e^{\beta_{i}+\beta_{j}+z_{ij}^{\top}\gamma}}\right).$$

(2)

The maximum likelihood equations are

$$\begin{array}[]{rcl}d_{i}&=&\sum_{j\neq i}\frac{e^{\beta_{i}+\beta_{j}+z_{ij}^{\top}\gamma}}{1+e^{\beta_{i}+\beta_{j}+z_{ij}^{\top}\gamma}},\\ \sum_{i<j}a_{ij}z_{ij}&=&\sum_{i<j}\frac{z_{ij}e^{\beta_{i}+\beta_{j}+z_{ij}^{\top}\gamma}}{1+e^{\beta_{i}+\beta_{j}+z_{ij}^{\top}\gamma}}.\end{array}$$

(3)

The R language provides a standard package “glm” to solve (3), which implements an iteratively reweighted least squares method for generalized linear models [McCullagh and Nelder, (1989)].

Given an original database $D$ with records of $n$ persons, we consider a randomized data releasing mechanism $Q$ that takes $D$ as input and outputs a sanitized database $S=(S_{1},\ldots,S_{\ell})$ for public use. As an illustrated example, the additive noise mechanism returns the answer $f(D)+z$ to the query $f(D)$, where $z$ is a random noise. Let $\epsilon$ be a positive real number and $\mathcal{S}$ denote the sample space of $Q$. The data releasing mechanism $Q$ is $\epsilon$-differentially private if for any two neighboring databases $D_{1}$ and $D_{2}$ that differ on a single element (i.e., the data of one person), and all measurable subsets $B$ of $\mathcal{S}$ [Dwork et al., (2006)],

$$Q(S\in B|D_{1})\leq e^{\epsilon}\times Q(S\in B|D_{2}).$$

This says the probability of an output $S$ given the input $D_{1}$ is less than that given the input $D_{2}$ multiplied by a privacy factor $e^{\epsilon}$. The privacy parameter $\epsilon$ is chosen according to the privacy policy, which controls the trade-off between privacy and utility. It is generally public. Smaller value of $\epsilon$ means more privacy protection.

Differential privacy requires that the distribution of the output is almost the same whether or not an individual’s record appears in the original database. We illustrate why it protects privacy with an example. Suppose a hospital wants to release some statistics on the medical records of their patients to the public. In response, a patient may wish to make his record omitted from the study due to a privacy concern that the published results will reveal something about him personally. Differential privacy alleviates this concern because whether or not the patient participates in the study, the probability of a possible output is almost the same. From a theoretical point, any test statistic has nearly no power for testing whether an individual’s data is in the original database or not [Wasserman and Zhou, (2010)].

What is being protected in the differential privacy is precisely the difference between two neighboring databases. Within network data, depending on the definition of the graph neighbor, differential privacy is divided into $k$-node differential privacy [Hay et al., (2009)] and $k$-edge differential privacy [Nissim et al., (2007)]. Two graphs are called neighbors if they differ in exactly $k$ edges, then differential privacy is $k$-edge differential privacy. The special case with $k=1$ is generally referred to as edge differential privacy. Analogously, we can define $k$-node differential privacy by letting graphs be neighbors if one can be obtained from the other by removing $k$ nodes and its adjacent edges. Edge differential privacy protects edges not to be detected, whereas node differential privacy protects nodes together with their adjacent edges, which is a stronger privacy policy. However, it may be infeasible to design algorithms that are both node differential privacy and have good utility since it generally needs a large noise [e.g., Hay et al., (2009)]. Following Hay et al., (2009) and Karwa and Slavković, (2016), we use edge differential privacy here.

Let $\delta(G,G^{\prime})$ be the harming distance between two graphs $G$ and $G^{\prime}$, i.e., the number of edges on which $G$ and $G^{\prime}$ differ. The formal definition of $(\epsilon,k)$-edge differential privacy is as follows.

Let $f:\mathcal{G}\rightarrow\mathbb{R}^{\ell}$ be a function. The global sensitivity [Dwork et al., (2006)] of the function $f$, denoted $\Delta f$, is defined below.

The global sensitivity measures the largest change for the query function $f$ in terms of the $L_{1}$-norm between any two neighboring graphs. The magnitude of noises added in the differentially private algorithm $Q$ crucially depends on the global sensitivity. If the outputs are the network statistics, then a simple algorithm to guarantee edge differential privacy is the Laplace Mechanism [e.g., Dwork et al., (2006)] that adds the Laplacian noise proportional to the global sensitivity of $f$.

When $f(G)$ is integer, one can use a discrete Laplace random variable as the noise, where it has the probability mass function:

$$\mathbb{P}(X=x)=\frac{1-\lambda}{1+\lambda}\lambda^{|x|},~{}~{}x\in\{0,\pm 1,\ldots\},\lambda\in(0,1).$$

(4)

Lemma 1 still holds if the continuous Laplace distribution is replaced by the discrete version and the privacy parameter is chosen by $\epsilon=-\Delta(f)\log\lambda$; see Karwa and Slavković, (2016).

We introduce a nice property on differential privacy: any function of a differentially private mechanism is also differentially private, as stated in the lemma below.

From the log-likelihood function (2), we know that $(d,\sum_{i<j}a_{ij}z_{ij})$ is the sufficient statistic. Thus, the private information in the covariates-adjusted $\beta$–model is $(d,\sum_{i<j}a_{ij}z_{ij})$. We use the continuous Laplace mechanism in Lemma 1 and its discrete version to release the covariate statistic $\sum_{i<j}z_{ij}a_{ij}$ and the degree sequence $d$ under $(\epsilon_{n}/2,k_{n})$-edge differential privacy, respectively. The joint mechanism satisfies $(\epsilon_{n},k_{n})$-edge differential privacy. The subscript $n$ means that $k_{n}$ and $\epsilon_{n}$ are allowed to depend on $n$. If we add $k_{1}$ or remove $k_{2}=k_{n}-k_{1}$ edges in $G_{n}$ and denote the induced graph as $G_{n}^{\prime}$, then

$$\|d-d^{\prime}\|_{1}=2k_{n},~{}~{}\|\sum_{i<j}(a_{ij}-a^{\prime}_{ij})z_{ij}\|_{1}\leq pk_{n}z_{*}$$

where $d^{\prime}$ is the degree sequence of $G_{n}^{\prime}$, $a^{\prime}_{ij}$ is the value of edge $(i,j)$ in $G_{n}^{\prime}$ and $z_{*}=\max_{ijk}|z_{ijk}|$. So the global sensitivity is $2k_{n}$ for $d$ and $pk_{n}z_{*}$ for $\sum_{i<j}a_{ij}z_{ij})$. We release the sufficient statistics $d$ and $y:=\sum_{i<j}a_{ij}z_{ij}$ as follows:

$\displaystyle\begin{array}[]{rcl}\tilde{d}_{i}&=&d_{i}+\xi_{i},~{}~{}i=1,\ldots,n,\\ \tilde{y}_{t}&=&\sum_{i<j}z_{ijt}a_{ij}+\eta_{t},~{}~{}t=1,\ldots,p,\end{array}$

(7)

where $\xi_{i}$, $i=1,\ldots,n$, are independently generated from the discrete Laplace distribution with $\lambda_{n1}=e^{-\epsilon_{n}/(4k_{n})}$, and $\eta_{t}$, $t=1,\ldots,p$, are independently generated from the Laplace distribution with $\lambda_{n2}=2pk_{n}z_{*}/\epsilon_{n}$.

Write $\mu(x)=e^{x}/(1+e^{x})$. Define

$$\pi_{ij}:=z_{ij}^{\top}\gamma+\beta_{i}+\beta_{j}.$$

(8)

It is clear that $\mu(\pi_{ij})$ is the expectation of $a_{ij}$. When we emphasize the arguments $\beta$ and $\gamma$ in $\mu(\cdot)$, we write $\mu_{ij}(\beta,\gamma)$ instead of $\mu(\pi_{ij})$. To estimate model parameters, we directly replace $d$ and $y$ in maximum likelihood equations (3) with their noisy observed values $\tilde{d}$ and $\tilde{y}$:

$$\begin{array}[]{rcl}\tilde{d}_{i}&=&\sum_{j\neq i}\mu_{ij}(\beta,\gamma),~{}~{}i=1,\ldots,n,\\ \tilde{y}&=&\sum_{i=1}^{n}\sum_{j=1,j<i}^{n}z_{ij}\mu_{ij}(\beta,\gamma).\end{array}$$

(9)

Because the expectations of the noises are zero, the above equation are the same as the moment equations.

Let $(\widehat{\beta},\widehat{\gamma})$ be the solution to the equations (9). Since $(\tilde{d},\tilde{y})$ satisfies $(\epsilon_{n},k_{n})$-edge differential privacy, $(\widehat{\beta},\widehat{\gamma})$ is also $(\epsilon_{n},k_{n})$-edge differentially private according to Lemma 2. A two-step iterative algorithm by alternating between solving the first equation in (9) via the fixed point method in Chatterjee et al., (2011) for a given $\gamma$ and solving the second equation in (9) via the Newton method or the gradient descent method, can be employed to obtain the solution.

To derive consistency of the differentially private estimator, let us first define a system of functions based on the estimating equations (9). Define

$$F_{i}(\beta,\gamma)=\sum\limits_{j=1,j\neq i}^{n}\mu_{ij}(\beta,\gamma)-\tilde{d}_{i},~{}~{}i=1,\ldots,n,$$

(12)

and $F(\beta,\gamma)=(F_{1}(\beta,\gamma),\ldots,F_{n}(\beta,\gamma))^{\top}$. Further, we define $F_{\gamma,i}(\beta)$ as the value of $F_{i}(\beta,\gamma)$ for an arbitrarily fixed $\gamma$ and $F_{\gamma}(\beta)=(F_{\gamma,1}(\beta),\ldots,F_{\gamma,n}(\beta))^{\top}$. Let $\widehat{\beta}_{\gamma}$ be a solution to $F_{\gamma}(\beta)=0$. Correspondingly, we define two functions for exploring the asymptotic behaviors of the estimator of the homophily parameter:

	$\displaystyle Q(\beta,\gamma)=\sum_{j<i}z_{ij}\mu_{ij}(\beta,\gamma)-\tilde{y},$		(13)
	$\displaystyle Q_{c}(\gamma)=\sum_{j<i}z_{ij}\mu_{ij}(\widehat{\beta}_{\gamma},\gamma)-\tilde{y}.$		(14)

$Q_{c}(\gamma)$ could be viewed as the profile function of $Q(\beta,\gamma)$ in which the degree parameter $\beta$ is profiled out. It is clear that

$$F(\widehat{\beta},\widehat{\gamma})=0,~{}~{}F_{\gamma}(\widehat{\beta}_{\gamma})=0,~{}~{}Q(\widehat{\beta},\widehat{\gamma})=0,~{}~{}Q_{c}(\widehat{\gamma})=0.$$

By the compound function derivation law, we have

	$\displaystyle 0=\frac{\partial F_{\gamma}(\widehat{\beta}_{\gamma})}{\partial\gamma^{\top}}=\frac{\partial F(\widehat{\beta}_{\gamma},\gamma)}{\partial\beta^{\top}}\frac{\partial\widehat{\beta}_{\gamma}}{\gamma^{\top}}+\frac{\partial F(\widehat{\beta}_{\gamma},\gamma)}{\partial\gamma^{\top}},$		(15)
	$\displaystyle\frac{\partial Q_{c}(\gamma)}{\partial\gamma^{\top}}=\frac{\partial Q(\widehat{\beta}_{\gamma},\gamma)}{\partial\beta^{\top}}\frac{\partial\widehat{\beta}_{\gamma}}{\gamma^{\top}}+\frac{\partial Q(\widehat{\beta}_{\gamma},\gamma)}{\partial\gamma^{\top}}.$		(16)

By solving $\partial\widehat{\beta}_{\gamma}/\partial\gamma^{\top}$ in (15) and substituting it into (16), we get the Jacobian matrix $Q_{c}^{\prime}(\gamma)$ $(=\partial Q_{c}(\gamma)/\partial\gamma)$:

$\displaystyle\frac{\partial Q_{c}(\gamma)}{\partial\gamma^{\top}}=\frac{\partial Q(\widehat{\beta}_{\gamma},\gamma)}{\partial\gamma^{\top}}-\frac{\partial Q(\widehat{\beta}_{\gamma},\gamma)}{\partial\beta^{\top}}\left[\frac{\partial F(\widehat{\beta}_{\gamma},\gamma)}{\partial\beta^{\top}}\right]^{-1}\frac{\partial F(\widehat{\beta}_{\gamma},\gamma)}{\partial\gamma^{\top}},$

(17)

where

$$\frac{\partial Q(\widehat{\beta}_{\gamma},\gamma)}{\partial\gamma^{\top}}:=\frac{\partial Q(\beta,\gamma)}{\partial\gamma^{\top}}{\Big{|}}_{\beta=\widehat{\beta}_{\gamma},\gamma=\gamma},~{}~{}\frac{\partial F(\widehat{\beta}_{\gamma},\gamma)}{\partial\beta^{\top}}:=\frac{\partial F(\beta,\gamma)}{\partial\beta^{\top}}{\Big{|}}_{\beta=\widehat{\beta}_{\gamma},\gamma=\gamma}.$$

The asymptotic behavior of $\widehat{\gamma}$ crucially depends on the Jacobian matrix $Q_{c}^{\prime}(\gamma)$. Since $\widehat{\beta}_{\gamma}$ does not have a closed form, conditions that are directly imposed on $Q_{c}^{\prime}(\gamma)$ are not easily checked. To derive feasible conditions, we define

$$H(\beta,\gamma)=\frac{\partial Q(\beta,\gamma)}{\partial\gamma}-\frac{\partial Q(\beta,\gamma)}{\partial\beta}\left[\frac{\partial F(\beta,\gamma)}{\partial\beta}\right]^{-1}\frac{\partial F(\beta,\gamma)}{\partial\gamma},$$

(18)

which is a general form of $\partial Q_{c}(\gamma)/\partial\gamma$. Note that $H(\widehat{\beta}_{\gamma},\gamma)$ is the Fisher information matrix of the concentrated likelihood function $\ell_{c}(\gamma)$, where the degree parameter $\beta$ is profiled out. When $\beta\in B(\beta^{*},\epsilon_{n1})$ and $b_{n}^{2}\kappa_{n}^{2}\epsilon_{n1}=o(1)$, we have the approximation:

$$\frac{1}{n^{2}}H(\beta,\gamma^{*})=\frac{1}{n^{2}}H(\beta^{*},\gamma^{*})+o(1),$$

whose proof is given in the supplementary material. We assume that there exists a number $\rho_{n}$ such that

$$\sup_{\beta\in B(\beta^{*},\epsilon_{n1})}\|H^{-1}(\beta,\gamma^{*})\|_{\infty}\leq\frac{\rho_{n}}{n^{2}}.$$

Note that the dimension of $H(\beta,\gamma)$ is fixed and every its entry is a sum $n(n-1)/2$ of terms. If $n^{-2}H(\beta,\gamma^{*})$ converges to a constant matrix, then $\rho_{n}$ is bounded. Moreover, if $n^{-2}H(\beta,\gamma^{*})$ is positively definite, then

$$\rho_{n}=\sqrt{p}\times\sup_{\beta\in B(\beta^{*},\epsilon_{n1})}1/\lambda_{\min}(\beta),$$

where $\lambda_{\min}(\beta)$ is the smallest eigenvalue of $n^{-2}H(\beta,\gamma^{*})$.

We use a two-stage Newton iterative sequence to show consistency. In the first stage, we obtain an upper bound of the error between $\widehat{\beta}_{\gamma}$ and $\beta$ in terms of the $\ell_{\infty}$ norm for a given $\gamma$. This is done by verifying the well-known Newton-Kantororich conditions, under which the optimal error bounds are established. Then we derive the upper bound of the error between $\widehat{\gamma}$ and $\gamma$ by using a profiled function $Q_{c}(\gamma)$ constructed from estimating equations. Now we formally state the consistency result.

The scalar factor $\tau_{n}$ appears due to that the magnitude of $\|Q_{c}(\gamma^{*})\|_{\infty}$ is $O(\tau_{n}n\log n)$. Note that the error bound of the parameter estimator in Theorem 3 of Karwa and Slavković, (2016) does not depend on the privacy parameter. Our result here characterizes how the error bound varies on $\epsilon_{n}$. In view of the above theorem, we present the consistency conditions and the error bounds under two special cases. The first case is that the parameters and covariates are bounded. The second is that $k_{n}/\epsilon_{n}$ goes to zero.

The asymptotic distribution of $\widehat{\beta}$ depends crucially on the inverse of the Fisher information matrix $V$ of $\beta$. Given $m,M>0$, we say an $n\times n$ matrix $V=(v_{ij})$ belongs to the matrix class $\mathcal{L}_{n}(m,M)$ if $V$ is a diagonally balanced matrix with positive elements bounded by $m$ and $M$, i.e.,

$$\begin{array}[]{l}v_{ii}=\sum_{j=1,j\neq i}^{n}v_{ij},~{}~{}i=1,\ldots,n,\\ 0<m\leq v_{ij}\leq M,~{}~{}i,j=1,\ldots,n;i\neq j.\end{array}$$

Clearly, $F^{\prime}_{\gamma}(\beta)$ belongs to the matrix class $\mathcal{L}(b_{n}^{-1},1/4)$ when $\beta\in B(\beta^{*},\epsilon_{n1})$ and $\gamma\in B(\gamma^{*},\epsilon_{n2})$. We will obtain the asymptotic distribution of the estimator $\widehat{\beta}$ through obtaining its asymptotic expression, which depends on the inverse of $F^{\prime}_{\gamma}(\beta)$. However, its inverse does not have a closed form. Yan et al., (2015) proposed to approximate the inverse $V^{-1}$ of $V$ by a diagonal matrix

$$S=\mathrm{diag}(1/v_{11},\ldots,1/v_{nn}),$$

(19)

and obtained the upper bound of the approximate error.

Note that $v_{ii}=\sum_{j\neq i}\mathrm{Var}(a_{ij})$. By the central limit theorem for the bounded case, as in Loéve, (1977, p.289), if $b_{n}=o(n^{1/2})$, then $v_{ii}^{-1/2}\{d_{i}-\mathbb{E}(d_{i})\}$ converges in distribution to the standard normal distribution. When considering the asymptotic behaviors of the vector $(d_{1},\ldots,d_{r})$ with a fixed $r$, one could replace the degrees $d_{1},\ldots,d_{r}$ by the independent random variables $\tilde{d}_{i}=d_{i,r+1}+\ldots+d_{in}$, $i=1,\ldots,r$. Therefore, we have the following proposition.

Part (2) follows from part (1) and the fact that

$$\lim_{r\to\infty}\limsup_{t\to\infty}\mathrm{Var}\left(\sum_{k=r+1}^{n}c_{i}\frac{d_{i}-\mathbb{E}(d_{i})}{\sqrt{v_{ii}}}\right)=0$$

by Theorem 4.2 of Billingsley, (1995). To prove the above equation, it suffices to show that the eigenvalues of the covariance matrix of $(d_{i}-\mathbb{E}(d_{i}))/v_{ii}^{1/2}$, $i=r+1,\ldots,n$ are bounded by 2 (for all $r<n$). This is true by the well-known Perron-Frobenius theory: if $A$ is a symmetric matrix with nonnegative elements, then its largest eigenvalue is less than the largest value of row sums.

We apply a second order Taylor expansion to $\sum_{j\neq i}\mu_{ij}(\widehat{\beta},\widehat{\gamma})$ to derive the asymptotic expression of $\widehat{\beta}$. In the expansion, the first order term is the sum of $V(\widehat{\beta}-\beta)$ and $V_{\gamma\beta}(\widehat{\gamma}-\gamma)$, where $V_{\gamma\beta}=\partial F(\beta,\gamma)/\partial\gamma^{\top}$. Since $V^{-1}$ does not have a closed form, we work with $S$ defined at (19) to approximate it. By Theorem 1, $\widehat{\gamma}$ has a $n^{-1}$ convergence rate up to a logarithm factor. This makes that the term $V_{\gamma\beta}(\widehat{\gamma}-\gamma)$ is an remainder term. The second order term in the expansion is also asymptotically neglect. Then we represent $\widehat{\theta}-\theta$ as the sum of $S(d-\mathbb{E}d)$ and a remainder term. The central limit theorem is proved by establishing the asymptotic normality of $S(d-\mathbb{E}d)$ and showing the remainder is negligible. We formally state the central limit theorem as follows.

Let $T_{ij}$ be an $n$-dimensional column vector with $i$th and $j$th elements ones and other elements zeros. Define

$$\begin{array}[]{c}V=\frac{\partial F(\beta^{*},\gamma^{*})}{\partial\beta^{\top}},~{}~{}V_{\gamma\beta}=\frac{\partial Q(\beta^{*},\gamma^{*})}{\partial\beta^{\top}},\\ s_{ij}(\beta,\gamma)=(a_{ij}-\mathbb{E}a_{ij})(z_{ij}-V_{\gamma\beta}V^{-1}T_{ij}).\end{array}$$

Note that $s_{ij}(\beta,\gamma)$, $i<j$, are independent vectors. By direct calculations,

$$(V_{\gamma\beta})_{kj}=\frac{\partial Q_{k}(\beta^{*},\gamma^{*})}{\partial\beta_{j}}=\sum_{i\neq j}z_{ijk}\mu^{\prime}(\pi^{*}_{ij}),$$

such that

$$\|V_{\gamma\beta}\|_{\infty}\leq\frac{z_{*}}{4}(n-1).$$

By Lemma 4, we have

$$\|V_{\gamma\beta}V^{-1}\|_{\infty}\leq\|V_{\gamma\beta}\|_{\infty}\|V^{-1}\|_{\infty}\leq\frac{z_{*}(n-1)}{4}\cdot\frac{b_{n}(3n-4)}{2(n-1)(n-2)}<b_{n}z_{*}.$$

Therefore, all $s_{ij}(\beta^{*},\gamma^{*})$ are bounded. Let $\bar{Q}=Q(\beta^{*},\gamma^{*})-\eta$ and $\bar{F}=F(\beta^{*},\gamma^{*})-\xi$. Note that

$$\mathrm{Cov}(\sum\nolimits_{i<j}s_{ij}(\beta^{*},\gamma^{*}))=\mathrm{Cov}(\bar{Q}-V_{\beta\gamma}^{\top}V^{-1}\bar{F})=H(\beta^{*},\gamma^{*}),$$

where $H(\beta,\gamma)$ is defined at (18). By the central limit theorem for the bounded case, as in Loéve, (1977, p.289), we have the following proposition.