Asymptotic Normality of Log Likelihood Ratio and Fundamental Limit of the Weak Detection for Spiked Wigner Matrices

For the spiked Wigner matrices,

• •

  (Asymptotic normality - Theorem 2.3) we prove that the log-LR converges to a Gaussian, and

• •

  (Optimal error - Corollary 2.5) we compute the limit of the sum of the Type-I and Type-II errors of the LR test.

Let $\mathbb{P}_{0}$ and $\mathbb{P}_{1}$ be the probability distributions of a spiked Wigner matrix $M$ under ${\boldsymbol{H}}_{0}$ and ${\boldsymbol{H}}_{1}$, respectively. If the noise is Gaussian, from the normalization in Definition 2.1 and Assumption 2.2, the likelihood ratio of $\mathbb{P}_{1}$ with respect to $\mathbb{P}_{0}$ is given by

$$\begin{split}&{\mathcal{L}}(M;\omega)=\frac{\mathrm{d}\mathbb{P}_{1}}{\mathrm{d}\mathbb{P}_{0}}\\ &:=\frac{1}{2^{N}}\sum_{{\boldsymbol{x}}}\prod_{i<j}\frac{\exp(-\frac{N}{2}(\sqrt{N}M_{ij}-\sqrt{\omega N}x_{i}x_{j})^{2})}{\exp(-\frac{N}{2}(\sqrt{N}M_{ij})^{2})}\prod_{k}\frac{\exp(-\frac{N}{2}(\sqrt{N}M_{kk}-\sqrt{\omega N}x_{k}^{2})^{2})}{\exp(-\frac{N}{2}(\sqrt{N}M_{kk})^{2})}\\ &=\frac{1}{2^{N}}\sum_{{\boldsymbol{x}}}\prod_{i<j}\exp\left(\sqrt{\omega}N^{2}M_{ij}x_{i}x_{j}-\frac{\omega}{2}\right)\prod_{k}\exp\left(N\sqrt{\omega}M_{kk}-\frac{\omega}{2}\right),\end{split}$$

(1.1)

where we used that $x_{i}^{2}=N^{-1}$ for all $i$. The random off-diagonal term in (LABEL:eq:Gaussian_LR),

$$\frac{1}{2^{N}}\sum_{{\boldsymbol{x}}}\prod_{i<j}\exp\left(\sqrt{\omega}N^{2}M_{ij}x_{i}x_{j}\right)=\frac{1}{2^{N}}\sum_{{\boldsymbol{x}}}\exp\left(\sqrt{\omega}N^{2}\sum_{i<j}M_{ij}x_{i}x_{j}\right),$$

coincides with the free energy of the Sherrington–Kirkpatrick model for which the asymptotic normality is known. (The diagonal part is independent from the off-diagonal part and its asymptotic (log)-normality can be easily proved by the CLT.)

If the noise is non-Gaussian, however, then we instead have

$${\mathcal{L}}(M;\omega)=\frac{1}{2^{N}}\sum_{{\boldsymbol{x}}}\prod_{i<j}\frac{p(\sqrt{N}M_{ij}-\sqrt{\omega N}x_{i}x_{j})}{p(\sqrt{N}M_{ij})}\prod_{k}\frac{p_{d}(\sqrt{N}M_{kk}-\sqrt{\omega N}x_{k}^{2})}{p_{d}(\sqrt{N}M_{kk})}\,,$$

(1.2)

where $p$ and $p_{d}$ are the densities of the (normalized) off-diagonal terms and the diagonal terms, respectively. Applying the Taylor expansion,

$$\frac{p(\sqrt{N}M_{ij}-\sqrt{\omega N}x_{i}x_{j})}{p(\sqrt{N}M_{ij})}\simeq 1-\sqrt{\omega N}\,\frac{p^{\prime}(\sqrt{N}M_{ij})}{p(\sqrt{N}M_{ij})}x_{i}x_{j}+\frac{\omega}{2N}\frac{p^{\prime\prime}(\sqrt{N}M_{ij})}{p(\sqrt{N}M_{ij})}\,,$$

and after comparing it to the Taylor expansion of the exponential function, it is possible to rewrite the leading order terms of the off-diagonal part of the LR as

$$\frac{1}{2^{N}}\sum_{{\boldsymbol{x}}}\exp\left(\sum_{i<j}(A_{ij}x_{i}x_{j}+B_{ij})\right),$$

(1.3)

where $A_{ij}$ and $B_{ij}$ depend only on $p(\sqrt{N}M_{ij})$ and its derivatives. (See (4.3)-(LABEL:eq:ABC) for a precise definition.) Since $A_{ij}$ and $B_{ij}$ are not independent, the asymptotic normality of (1.3) is not obvious, and in fact, requires hard analysis.

The main technical difficulty in the analysis of the log-LR arises from that the mechanism for the Gaussian convergence of the terms containing $A_{ij}$, which we call the spin glass term, and those with $B_{ij}$, which we call the CLT term, are markedly different. While the latter can be readily proved by the CLT, the former is a result from the theory of spin glass based on the combinatorial analysis for the graph structure involving $A_{ij}$’s. To overcome the difficulty, adopting the strategy of [1], we decompose the spin glass term, $2^{-N}\sum_{{\boldsymbol{x}}}\exp(\sum_{i<j}(A_{ij}x_{i}x_{j}))$, into two parts, which are asymptotically orthogonal to each other, conditional on $(B_{ij})$. (See Proposition 4.4 for more detail.) We can then verify the part that depends on $(B_{ij})$ and express the dependence by a conditional expectation on $(B_{ij})$. Combining the $B$-dependent part of the spin glass term with the CLT term, we can prove the Gaussian convergence with a precise formula for the mean and the variance of the limiting distribution.

From the asymptotic normality, it is immediate to obtain the sum of the Type-I and Type-II errors of the LR test. The limiting error in Corollary 2.5 converges to $0$ as $\lambda\nearrow(1/F)$, which shows that the detection threshold is given by $\lambda_{c}=1/F$. While the limiting error of the LR test for the case $\lambda>1/F$ also converges to $0$, which can be indirectly deduced from the fact that the transformed PCA in [23] can reliably detect the signal, our result is not directly applicable to this case; the asymptotic normality holds only in the regime $\lambda<1/F$ and it is expected the asymptotic behavior of the log-LR changes if $\lambda>1/F$, which corresponds to the low temperature regime in the theory of spin glass.

The Rademacher prior we consider in this work is one of the most natural models for the rank-one signal-plus-noise data, and it also corresponds to the most fundamental spin configuration in the SK model. Our method relies on this special structure of the spike at least in the technical level. We believe that our main result would change with different priors, e.g., the spherical prior instead of the Rademacher prior, but the change is restricted to a subleading term (proportional to $\omega^{2}$) in the mean and the variance of the limiting Gaussian. Note that such a change would not occur with the Gaussian noise where the quadratic term (of $\omega$) vanishes, which can be made rigorously since the asymptotic behavior of the log-LR with the spherical prior coincides that with many i.i.d. priors including the Rademacher prior [13]. See the discussion at the end of Section 2.2 and also the numerical experiments in Section 3 for details.

For the spiked IID matrices, we also prove the asymptotic normality of the log-LR (Theorem 2.8) and compute the limit of the sum of the Type-I and Type-II errors of the LR test (Corollary 2.9).

Suppose that the data matrix $Y$ is of the form $Y=X+\sqrt{\lambda}{\boldsymbol{x}}{\boldsymbol{x}}^{T}$, where the signal ${\boldsymbol{x}}$ is an $N$-dimensional vector and the noise $X$ is an $N\times N$ IID matrix. (See Definition 2.7 for a precise definition.) If the noise is Gaussian, similar to the spiked Wigner case in (LABEL:eq:Gaussian_LR), we find that the likelihood ratio

$$\begin{split}{\mathcal{L}}(Y;\omega)&:=\frac{\mathrm{d}\mathbb{P}_{1}}{\mathrm{d}\mathbb{P}_{0}}:=\frac{1}{2^{N}}\sum_{{\boldsymbol{x}}}\prod_{i,j=1}^{N}\frac{\exp(-\frac{N}{2}(\sqrt{N}Y_{ij}-\sqrt{\omega N}x_{i}x_{j})^{2})}{\exp(-\frac{N}{2}(\sqrt{N}Y_{ij})^{2})}\\ &=\frac{1}{2^{N}}\sum_{{\boldsymbol{x}}}\prod_{i<j}\exp\left(\sqrt{\omega}N^{2}(Y_{ij}+Y_{ji})x_{i}x_{j}-\omega\right)\prod_{k}\exp\left(N\sqrt{\omega}Y_{kk}-\frac{\omega}{2}\right).\end{split}$$

(1.4)

Since $(Y_{ij}+Y_{ji})$ is again a Gaussian random variable whose variance is $2N^{-1}$, under ${\boldsymbol{H}}_{0}$, the exponent in the off-diagonal term $(\sqrt{\omega}N^{2}(Y_{ij}+Y_{ji})x_{i}x_{j}-\omega)$ is equal in distribution to $(\sqrt{2\omega}N^{2}Y_{ij}x_{i}x_{j}-\omega)$. Thus, when compared to the exponent in the off-diagonal term in (LABEL:eq:Gaussian_LR), we find that the SNR $\omega$ is effectively doubled, and in particular, reliable detection is impossible if $\omega<1/2$. Note that in this case it is possible to reliably detect the signal if $\omega>1/2$ by PCA for the symmetrized matrix $(Y+Y^{T})/\sqrt{2}$, whereas the largest singular value of $M$ is separated from other singular values only when $\omega>1$ (see, e.g., Section 3.1 of [6]).

If the noise is non-Gaussian, however, we face the same issue as in the spiked Wigner matrices; the log-LR can be decomposed into the spin glass term and the CLT term as in (1.3), which are not independent. The issue can be resolved by adapting the analysis for the spiked Wigner matrices, where we need to consider $A_{ij}+A_{ji}$ and $B_{ij}+B_{ji}$ in place of $A_{ij}$ and $B_{ij}$. It is then possible to prove the asymptotic normality of the log-LR with precise formulas for the mean and the variance.

The results in Theorem 2.8 and Corollary 2.9 show that reliable detection is impossible if $\lambda<1/(2F)$. On the other hand, if we transform the data matrix entrywise and then symmetrize the transformed matrix, the largest eigenvalue of the resulting matrix becomes an outlier when $\lambda>1/(2F)$, and hence reliable detection is possible by PCA. Thus, we can conclude that the threshold $\lambda_{c}=1/(2F)$ for the reliable detection. See the discussion at the end of Section 2.3 for details.

The phase transition of the largest eigenvalue of the spiked Wigner matrix has been extensively studied in random matrix theory. It is known as the BBP transition after the seminal work by Baik, Ben Arous, and Péché [2] for spiked (complex) Wishart ensembles. Corresponding results for the spiked Wigner matrix were proved under various assumptions [22, 14, 7, 5]. For the behavior of the eigenvector associated with the largest eigenvalue of the spiked Wigner matrix, we refer to [5].

The impossibility of the reliable detection for the case $\lambda<1$ with the Gaussian noise was considered by Montanari, Reichman, and Zeitouni [18] and Johnstone and Onatski [16]. The optimality of the PCA for the Gaussian noise, the sub-optimality of the PCA and the optimality of the transformed PCA for the non-Gaussian noise was proved by Perry, Wein, Bandeira, and Moitra [23] under mild assumptions on the prior. We remark that the threshold $\lambda_{c}$ may not be exactly $1$ depending on the prior [13].

The analysis on the LR test and its fundamental limit were studied by El Alaoui, Krzakala, and Jordan [13] under the assumption that the signal is an i.i.d. random vector and the noise is Gaussian. We remark that there exists a non-LR test based on the linear spectral statistics of the data matrix, which is optimal if the noise is Gaussian [9].

The detection problem for the spiked rectangular model, introduced by Johnstone [15], has also been widely studied. We refer to [19, 20, 12, 17, 10] for the weak detection under Gaussian noise, including the optimal error of the hypothesis test, and [23, 17] for the sub-optimality of the PCA under non-Gaussian noise. For the spiked IID matrices and their applications, we refer to [8] and references therein.

The SK model is one of the most fundamental models of spin glass. The model was introduced by Sherrington and Kirkpatrick [24] as a mean-field version of the Edwards–Anderson model [11]. The (deterministic) limit of the free energy of the SK model was first predicted by Parisi [21], and it was later rigorously proved by Talagrand [25]. The fluctuation of the free energy for the SK model was studied by Aizenman, Lebowitz, and Ruelle [1]; they showed the Gaussian convergence of the fluctuation of the free energy in the high temperature case, which also holds with non-Gaussian noise. To the best of our best knowledge, the fluctuation in the low temperature case, including the zero-temperature case, still remains open.

The rest of the paper is organized as follows. In Section 2, we precisely define the model and state the main results. In Section 3, we present some examples and conduct numerical experiments to compare the outcomes with the theoretical results. In Sections 4 and 5, we prove our main results on the spiked Wigner matrices. We conclude the paper in Section 6 by summarizing our results and proposing future research directions. Some details of the proof, including the proof for the results on the spiked IID matrices, can be found in Appendix.

We begin by precisely defining the model we consider in this paper. We suppose that the noise matrix is a Wigner matrix, which is defined as follows:

For the data matrix $M$, we assume that it is a spiked Wigner matrix, i.e., $M=H+\sqrt{\lambda}{\boldsymbol{x}}{\boldsymbol{x}}^{T}$, and also additionally assume the following:

Note that we assume the off-diagonal entries are identically distributed. While it is possible to consider a spiked Wigner matrix with non-identically distributed noise matrix by assuming certain conditions on their density functions, we refrain from it to avoid unnecessary complications. Similarly, we remark that some conditions in Definition 2.1 and Assumption 2.2 are due to technical reasons, especially the finite moment condition and the symmetric condition for the density of the noise, and we believe that it is possible to relax the conditions. However, for the sake of brevity, we do not attempt to relax them further in this paper.

Recall that the SNR $\omega$ for the alternative ${\boldsymbol{H}}_{1}$ is fixed. Our main result is the following asymptotic normality for the log-LR statistic when testing between ${\boldsymbol{H}}_{0}:\lambda=0$ versus ${\boldsymbol{H}}_{1}:\lambda=\omega>0$.

If the off-diagonal density $p$ is Gaussian, we find that $F=1$, $G=2$, and the result in Theorem 2.3 coincides with Theorem 2 (and Remark below it) of [13], which states that $\log{\mathcal{L}}(M;\lambda)$ converges in distribution to ${\mathcal{N}}(-\rho^{\prime},2\rho^{\prime})$ where

$$\rho^{\prime}=-\frac{1}{4}\big{(}\log(1-\omega)+\omega(1-2F_{d})\big{)}.$$

From the comparison with the Gaussian case, ignoring the term $(\omega^{2}/4)(2F^{2}-G)$ in (2.2), the result in Theorem 2.3 heuristically shows that the effective SNR is $\omega F$ when the noise is non-Gaussian.

The change of the effective SNR with non-Gaussian noise can be understood as follows: under ${\boldsymbol{H}}_{1}$, if we apply a function $q$ to a normalized entry $\sqrt{N}M_{ij}$, then by the Taylor expansions,

$$q(\sqrt{N}M_{ij})=q(\sqrt{N}H_{ij}+\sqrt{\omega N}x_{i}x_{j})\approx q(\sqrt{N}H_{ij})+\sqrt{\omega N}q^{\prime}(\sqrt{N}H_{ij})x_{i}x_{j}.$$

(2.3)

Approximating the right-hand side further by

$$\begin{split}q(\sqrt{N}H_{ij})+\sqrt{\omega N}q^{\prime}(\sqrt{N}H_{ij})x_{i}x_{j}&\approx q(\sqrt{N}H_{ij})+\sqrt{\omega N}\mathbb{E}[q^{\prime}(\sqrt{N}H_{ij})]x_{i}x_{j}\\ &=\sqrt{N}\left(\frac{q(\sqrt{N}H_{ij})}{\sqrt{N}}+\sqrt{\omega}\mathbb{E}[q^{\prime}(\sqrt{N}H_{ij})]x_{i}x_{j}\right),\end{split}$$

we find that the matrix obtained by transforming $M$ entrywise via $q$ is of the form $\sqrt{N}(Q+\sqrt{\omega^{\prime}}{\boldsymbol{x}}{\boldsymbol{x}}^{T})$, with $Q_{ij}=q(\sqrt{N}H_{ij})/\sqrt{N}$ and $\sqrt{\omega^{\prime}}=\sqrt{\omega}\mathbb{E}[q^{\prime}(\sqrt{N}H_{ij})]$, which shows that the SNR is effectively changed. The approximation can be made rigorous [23], and after optimizing $q$, it can be shown that effective SNR is $\omega F$, which is independent of the prior under mild assumptions. We remark that the optimal $q$ is given by $-p^{\prime}/p$, and it is unique up to a constant factor.

An immediate consequence of Theorem 2.3 is the following estimate for the error probability of the LR test.

The proof of Corollary 2.5 from Theorem 2.3 is standard; we refer to the proof of Theorem 2 of [9].

Note that the error $\mathrm{err}(\omega)$ converges to $0$ as $\omega F\nearrow 1$, which is consistent with the discussion below Theorem 2.3 that the effective SNR is $\omega F$.

The most non-trivial part in Theorem 2.3 and Corollary 2.5 is the third term of $\rho$ in (2.2),

$$\frac{\omega^{2}}{4}(2F^{2}-G).$$

(2.4)

We conjecture that this term is due to the Rademacher prior and it would change if we assume a different prior. Our heuristic argument for it is as follows. In [9], an algorithm was proposed for hypothesis testing based on the linear spectral statistics (LSS) of the data matrix $M$, which is optimal if the noise is Gaussian; the LSS of the matrix $M$ whose eigenvalues are $\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{N}$ is

$$\sum_{i=1}^{N}f(\lambda_{i})$$

for a function $f$ continuous on an open neighborhood of $[-2,2]$. If the noise is non-Gaussian, the algorithm in [9] can be improved by pre-transforming the data matrix entrywise as in (2.3). The error of the improved algorithm is

$$\mathrm{erfc}\left(\sqrt{\widetilde{\rho}}/2\right),$$

(2.5)

where

$$\widetilde{\rho}=-\frac{1}{4}\left(\log(1-\omega F)+\omega(F-2F_{d})+\frac{\omega^{2}}{4}(2F^{2}-\widetilde{G})\right)$$

with

$$\widetilde{G}=\left(\int_{-\infty}^{\infty}\frac{(p^{\prime}(x))^{2}p^{\prime\prime}(x)}{p(x)}\mathrm{d}x\right)^{2}\Big{/}\left(\frac{3}{2}\int_{-\infty}^{\infty}\frac{(p^{\prime}(x))^{2}p^{\prime\prime}(x)}{p(x)}\mathrm{d}x-F^{2}\right).$$

The error (2.5) differs from the error in Corollary 2.5 only in the term involving $\omega^{2}$. The hypothesis test based on the eigenvalues are spherical in nature, since it does not depend on the eigenvectors, or specific directions. It means that the prior naturally associated with the LSS-based test is the spherical prior, which is the uniform distribution on the unit sphere. We thus conjecture that the error term in (2.4) depends on the prior of the model.

The results on the spiked Wigner matrices in Theorem 2.3 can be naturally extended to the spiked IID matrices defined as follows:

For a spiked IID matrix of the form $Y=X+\sqrt{\lambda}{\boldsymbol{x}}{\boldsymbol{x}}^{T}$, we have following result.

See Appendix C for the proof of Theorem 2.8. From Theorem 2.8, we can obtain the following estimate for the error probability of the LR test.

We omit the proof of Corollary 2.9; see the proof of Theorem 2 of [9].

Comparing Theorem 2.8 with Theorem 2.3, we find that the effective SNR is $2\omega F$, which is doubled from the effective $\omega F$ for the spiked Wigner matrices. The doubling of the effective SNR can also be heuristically checked from the LR as follows. By definition,

$$\begin{split}{\mathcal{L}}(Y;\omega)=\frac{1}{2^{N}}\sum_{{\boldsymbol{x}}}\prod_{i\neq j}\frac{p(\sqrt{N}Y_{ij}-\sqrt{\omega N}x_{i}x_{j})}{p(\sqrt{N}Y_{ij})}\prod_{k}\frac{p(\sqrt{N}X_{kk}-\sqrt{\omega N}x_{k}^{2})}{p(\sqrt{N}X_{kk})}\,,\end{split}$$

(2.9)

The off-diagonal term can be further decomposed into the upper-diagonal term and the lower-diagonal term, i.e.,

$$\prod_{i\neq j}\frac{p(\sqrt{N}Y_{ij}-\sqrt{\omega N}x_{i}x_{j})}{p(\sqrt{N}Y_{ij})}=\prod_{i<j}\frac{p(\sqrt{N}Y_{ij}-\sqrt{\omega N}x_{i}x_{j})}{p(\sqrt{N}Y_{ij})}\prod_{i>j}\frac{p(\sqrt{N}Y_{ij}-\sqrt{\omega N}x_{i}x_{j})}{p(\sqrt{N}Y_{ij})}.$$

The upper-diagonal term and the lower-diagonal term are almost independent, which effectively makes the LR becomes squared and the log-LR doubled. However, these terms are not exactly independent due to the spike $x_{i}x_{j}$, and this is why the quadratic term $(\omega^{2}/2)(2F^{2}-G)$ in $\rho^{*}$ is not the same as the four times the quadratic term $(\omega^{2}/4)(2F^{2}-G)$, which is what one would get by changing $\omega$ to $2\omega$ in (2.2). Note that the linear term $\omega(F-2F_{d})$ vanishes in (2.8) due to the assumption that $F_{d}=F$, since the term $\omega F$, arising from the off-diagonal terms, gets doubled for the spiked IID matrices but $2\omega F_{d}$ does not as it originates from the diagonal terms.

Theorem 2.8 implies the impossibility of the reliable detection when $\lambda<1/(2F)$. To conclude that the threshold $\lambda_{c}=1/(2F)$ for the reliable detection, we can consider the following variant of the PCA. We first apply the function $q$ entrywise to the data matrix as in (2.3). After normalizing the transformed matrix by dividing it by $\sqrt{N}$, it is approximately a spiked IID matrix. In the next step, we symmetrize the (normalized) transformed matrix, i.e., we compute

$$\widetilde{Y}_{ij}=\frac{q(\sqrt{N}Y_{ij})+q(\sqrt{N}Y_{ji})}{\sqrt{2N}}.$$

From the Taylor expansion, we find that

$$\begin{split}\frac{q(\sqrt{N}Y_{ij})+q(\sqrt{N}Y_{ji})}{\sqrt{2N}}\approx\frac{q(\sqrt{N}X_{ij})+q(\sqrt{N}X_{ji})}{\sqrt{2N}}+\sqrt{2\omega}\mathbb{E}[q^{\prime}(\sqrt{N}X_{ij})]x_{i}x_{j},\end{split}$$

which is approximately a spiked Wigner matrix. Following the proof in [23], it is then immediate to prove that the reliable detection is possible by applying the PCA to $\widetilde{Y}$ if $\lambda_{c}>1/(2F)$.

Assume ${\boldsymbol{H}}_{0}$. Recall that

$${\mathcal{L}}(M;\omega)=\frac{1}{2^{N}}\sum_{{\boldsymbol{x}}}\prod_{i<j}\frac{p(\sqrt{N}M_{ij}-\sqrt{\omega N}x_{i}x_{j})}{p(\sqrt{N}M_{ij})}\prod_{k}\frac{p_{d}(\sqrt{N}M_{kk}-\sqrt{\omega N}x_{k}^{2})}{p_{d}(\sqrt{N}M_{kk})}\,.$$

(4.3)

We first focus on the off-diagonal part of (4.3) and assume $i<j$. Recall the definition of $P^{(s)}_{ij}$ in Lemma 4.3. Note that $P^{(s)}_{ij}$ is an odd function of $\sqrt{N}M_{ij}$ if $s$ is odd, whereas it is an even function if $s$ is even. Recall that $x_{i}^{2}=N^{-1}$. Since $P^{(s)}_{ij}={\mathcal{O}}(1)$ from Lemma 4.2, by the Taylor expansion,

$$\begin{split}\frac{p(\sqrt{N}M_{ij}-\sqrt{\omega N}x_{i}x_{j})}{p(\sqrt{N}M_{ij})}=1-\sqrt{\omega N}P^{(1)}_{ij}x_{i}x_{j}+\frac{\omega P^{(2)}_{ij}}{2N}-\frac{\omega\sqrt{\omega}}{6\sqrt{N}}P^{(3)}_{ij}x_{i}x_{j}+\frac{\omega^{2}P^{(4)}_{ij}}{24N^{2}}+{\mathcal{O}}(N^{-\frac{5}{2}}),\end{split}$$

(4.4)

uniformly on $i$ and $j$. Taking the logarithm and Taylor expanding it again, we obtain

$$\begin{split}&\log\left(\frac{p(\sqrt{N}M_{ij}-\sqrt{\omega N}x_{i}x_{j})}{p(\sqrt{N}M_{ij})}\right)\\ &=-\sqrt{\omega N}x_{i}x_{j}\left(P^{(1)}_{ij}+\frac{\omega}{6N}\Big{(}P^{(3)}_{ij}-3P^{(1)}_{ij}P^{(2)}_{ij}+2(P^{(1)}_{ij})^{3}\Big{)}\right)+\frac{\omega}{2N}\left(P^{(2)}_{ij}-(P^{(1)}_{ij})^{2}\right)\\ &\qquad+\frac{\omega^{2}}{24N^{2}}\left(P^{(4)}_{ij}-3(P^{(2)}_{ij})^{2}-4P^{(1)}_{ij}P^{(3)}_{ij}+12(P^{(1)}_{ij})^{2}P^{(2)}_{ij}-6(P^{(1)}_{ij})^{4}\right)+{\mathcal{O}}(N^{-\frac{5}{2}}).\end{split}$$