Tests of Missing Completely At Random based on sample covariance matrices
Alberto Bordino and Thomas B. Berrett
Department of Statistics, University of Warwick

We study the problem of testing whether the missing values of a potentially high-dimensional dataset are Missing Completely at Random (MCAR). We relax the problem of testing MCAR to the problem of testing the compatibility of a sequence of covariance matrices, motivated by the fact that this procedure is feasible when the dimension grows with the sample size. Tests of compatibility can be used to test the feasibility of positive semi-definite matrix completion problems with noisy observations, and thus our results may be of independent interest.

Our first contributions are to define a natural measure of the incompatibility of a sequence of correlation matrices, which can be characterised as the optimal value of a Semi-definite Programming (SDP) problem, and to establish a key duality result allowing its practical computation and interpretation.

By studying the concentration properties of the natural plug-in estimator of this measure, we introduce novel hypothesis tests that we prove have power against all distributions with incompatible covariance matrices.

The choice of critical values for our tests rely on a new concentration inequality for the Pearson sample correlation matrix, which may be of interest more widely. By considering key examples of missingness structures, we demonstrate that our procedures are minimax rate optimal in certain cases.

We further validate our methodology with numerical simulations that provide evidence of validity and power, even when data are heavy tailed.

§ INTRODUCTION

Incomplete data are a common occurrence in almost all areas of statistical application, and the mechanisms leading to such data are diverse. For example, subjects in a survey may choose not to respond to certain questions, leading to missing values, or a practitioner may wish to combine data collected in different studies, where different variables were recorded in each. With incomplete data, traditional approaches become unreliable or even inapplicable, leading to a significant effect on the conclusions that can be drawn from the data. The most common approaches to dealing with missing values are to remove any incomplete observations, and thus to perform a complete-case analysis, or to replace any missing entry with a representative value, using an imputation method <cit.>. However, the validity of such procedures, and the choice of an appropriate one, depends crucially on the mechanism that determines the missingness. Mechanisms have traditionally been classified as Missing Completely At Random (MCAR), Missing At Random (MAR) and Missing Not At Random (MNAR) <cit.> according to the dependence structure between the variables themselves and their missingness, with such assumptions being required to link observations to targets of inference.

The typical formal setting is to suppose that we observe independent and identically distributed copies of a random object $\mathbf{X} \circ \Omega$, where $\mathbf{X}$ takes values in some product space $\mathcal{X} = \prod_{j=1}^d \mathcal{X}_j$, where $\Omega$ takes values in $\{0,1\}^d$ and where we define the operator $\circ$ by

\[(x \circ \omega)_j =
\begin{cases}
x_j \quad \text{ if } \omega_j = 1, \\
\text{NA} \quad \text{ if } \omega_j = 0. \\
\end{cases}
\]

The assumptions named above then control the dependence between the data $\mathbf{X}$ and the missingness indicator $\Omega$. The simplest case of MCAR is when these are independent, denoted $\mathbf{X} \perp \!\!\! \perp \Omega$, which essentially means that the data we observe is representative of the population, even if it is incomplete. For example, consider the simple problem of estimating $\mathbb{E}[XY]$ from complete data $(X_1, Y_1), \ldots, (X_N, Y_N)$ and incomplete data $(X_{N+1}, \text{NA}), \ldots, $ $(X_n, \text{NA})$. It is easy to see that under MCAR the complete case estimator $\hat{\mu}^{CC} := N^{-1} \sum_{i=1}^N X_i Y_i$ is an unbiased estimator of $\mathbb{E}[XY]$. When MCAR or similar assumptions hold, we can often employ statistical methodologies that are easy to interpret and make good use of all incomplete data, and solid theoretical guarantees have been developed in various modern statistical problems such as high-dimensional regression [Loh and Wainwright, 2012],
high-dimensional or sparse principal component analysis [Zhu et al., 2022, Elsener and van de Geer, 2018], classification [Cai and Zhang, 2019, Sell et al., 2023], and precision matrix and changepoint estimation [Follain et al., 2022]. However, if MCAR does not hold, which is common in practice, alternative methods may be required.

Hypothesis tests can be used to guide practitioners in deciding whether or not missingness assumptions are reasonable. The goal of this work is to study the problem of testing the hypothesis of MCAR, which has been the subject of much research in the missing data literature. Most prior work has been developed within the context of parametric models. For example, [Little, 1988], works under the hypothesis that the data is Gaussian in the setting that all pairs of variable are observed together (see Section <ref> for further details). [Fuchs, 1982] considers discrete data in the setting that a large number of complete cases are available. In both cases the methods are likelihood ratio tests, with the MLEs calculated using the EM algorithm [Dempster et al., 1977] and validity and power guarantees based on classical asymptotics. More recently, Berrett and Samworth, 2023 provided a nonparametric formulation of the problem and methodology that was proved to be widely powerful under minimal assumptions.
The key insight of Berrett and Samworth, 2023 is to relate the problem of testing MCAR to the problem of testing compatibility. For $S \subseteq [d]:=\{1,\ldots,d\}$ denote by $\{\Omega=\mathbbm{1}_S\}$ the event that $X_j$ is observed if and only if $j \in S$, and write $\mathbb{S}=\{S : \mathbb{P}(\Omega=\mathbbm{1}_S) >0\}$ for the set of all possible observation patterns. Then, under MCAR, the distribution $P_S$ of the observation $X_S|\{\Omega = \mathbbm{1}_S \}$ is equal to the marginal distribution of the population distribution $ \mathcal{L}(\mathbf{X})$ on $\mathcal{X}_S := \prod_{j \in S} \mathcal{X}_j$. Hence, if $P_{\mathbb{S}}:=(P_S : S \in \mathbb{S})$ is incompatible, in the sense that there is no distribution $P$ on $\mathcal{X}$ with marginal distribution $P_S$ on $\mathcal{X}_S$ for all $S \in \mathbb{S}$, then the data cannot be MCAR. In fact, it is shown that this reasoning is tight in that it is not possible to rule out MCAR based on observations of $\mathbf{X} \circ \Omega$ if $P_\mathbb{S}$ is compatible. In general, fully testing the compatibility of a sequence of distributions requires us to look at complex interactions between the distributions, and methods for doing so will have sample complexity that is exponential in the dimension $d$. Our work aims to provide methods that are valid and powerful without strong assumptions while being effective as the dimension $d$ grows.

Our methodology will be based on testing the compatibility of sequences of covariance matrices, which can be estimated consistently even for large $d$. Earlier studies have employed the covariance matrix to assess MCAR. As briefly discussed above, [Little, 1988] studied a likelihood ratio test of MCAR, effectively examining the homogeneity of means and covariances under the assumption of normality. However, Little expressed scepticism about its effectiveness unless the sample size is exceptionally large and the assumption of normality holds. This scepticism was further validated in simulations by [Kim and Bentler, 2002], who also developed a test for consistency of means and covariances based on generalised least squares. Both of these approaches work by comparing the sample covariance matrix associated to a given missingness pattern to the corresponding submatrix of an estimated complete covariance matrix. More recently, [Jamshidian and Jalal, 2010] developed $k$-sample tests of the equality of covariance matrices, given complete data, based on Hawkins' test [Hawkins, 1981]. Using empirical evidence, they then argued that these tests could be combined with imputation techniques to test the homogeneity of covariance matrices calculated using incomplete data. These methodologies can be effective when the corresponding assumptions are met and when a complete covariance matrix can be consistently estimated.

Our method works by directly checking the compatibility of the observed sample covariance matrices, making no assumptions on the form of the underlying distributions and not requiring the estimation of a complete covariance matrix. In particular, this second point means that our test can be applied with any collection $\mathbb{S}$ of missingness patterns. More precisely, at the population level, we will consider $\Sigma_{\mathbb{S}} = (\Sigma_S:S \in \mathbb{S})$, a sequence of suitably-normalised covariance matrices $\Sigma_S$ associated to the law of $\mathbf{X}_S|\{\Omega = \mathbbm{1}_S\}$, and design a statistical test to check if $\Sigma_{\mathbb{S}}$ is compatible, meaning that each $\Sigma_S$ can be obtained by marginalising a general $d \times d$ covariance matrix $\Sigma$, i.e. $(\Sigma)_S = \Sigma_S$. If MCAR holds then for each $S \in \mathbb{S}$ we must have $\mathrm{Cov}(\mathbf{X}_S |\Omega = \mathbbm{1}_S) = (\mathrm{Cov}(\mathbf{X}))_S$, the block of the covariance matrix of $\mathbf{X}$ corresponding to the variables in $S$, so that the sequence $(\Sigma_S : S \in \mathbb{S})$ must be compatible. Hence, if we can reject the hypothesis $H_0: \Sigma_{\mathbb{S}} \textit{ compatible}$, then we can reject the hypothesis of MCAR. See Figure <ref> for a pictorial summary of the key concepts so far.

More generally, one can consider the problem of testing the compatibility of moments of order $p \geq 1$ and, if it is found that these moments are incompatible, one can reject MCAR. For $p = 1$, this problem reduces to testing the compatibility of mean vectors, which essentially boils down to testing the equality of means. This has been studied in the statistical literature for nearly a century, and we refer to existing methods for solving this problem <cit.>. In order to have power against a wider range of alternatives, while limiting the complexity of the testing procedure, we restrict attention in this work to the natural $p=2$ problem. Here there are still various ways in which compatibility can fail. For example, we can rule out $H_0$ if $\Sigma_{\mathbb{S}}$ is inconsistent, in the sense that that there are two observation patterns $S_1, S_2 \in \mathbb{S}$ for which $(\Sigma_{S_1})_{S_1 \cap S_2} \neq (\Sigma_{S_2})_{S_1 \cap S_2}$, meaning that there exists a pair of variables whose covariance takes different values in different observation patterns. Testing the consistency of covariance matrices reduces to testing the equality of smaller covariance matrices, which has again been previously studied <cit.>. The corresponding nonparametric problem of testing the consistency of distributions was studied by [Li and Yu, 2015, Spohn et al., 2021] in the context of testing MCAR. However, there exist $\Sigma_{\mathbb{S}}$ that are consistent but not compatible. As a concrete example, consider the case where $d=3$, where $\mathbb{S}=\{\{1,2\},\{1,3\},\{2,3\} \}$, and where

\[
\Sigma_{ij} = \begin{pmatrix} 1 & \rho_{ij} \\ \rho_{ij} & 1 \end{pmatrix}
\]

with $\rho_{23}=\rho_{13}=-\rho_{12}=\rho$. Then $\Sigma_{\mathbb{S}}$ is compatible if and only if $\rho \leq 1/2$, even though it is always consistent.

Our framework. We relax the methodology in [Berrett and Samworth, 2023] and consider covariance matrices instead of full distributions. The price we pay is to create an extra ring (red area) in which we cannot detect departure from $H_0$ just by looking at $\Sigma_\mathbb{S}$ or $\sigma^2_\mathbb{S}$, the sequence of variances. For example, if the sequence of third-moment tensors were inconsistent, but $\Sigma_\mathbb{S}$ were compatible and $\sigma^2_\mathbb{S}$ consistent, we would not be able to reject MCAR, although $P_\mathbb{S}$ would be incompatible.

The above example is relatively simple because any pair of variables is observed together so that the full covariance matrix can be estimated. However, we can characterise compatibility for any $\mathbb{S}$ (see Proposition <ref>). While the compatibility of distributions can be characterised using linear programming <cit.>, characterising the compatibility of covariance matrices requires ideas from semi-definite programming (SDP), which studies linear optimisation problems over spectrahedra <cit.>. If $\Sigma_{\mathbb{S}}$ is consistent, compatibility is equivalent to the feasibility of a positive semi-definite matrix completion problem, where we observe a partial symmetric matrix $A = (a_{ij})$ for positions $(i, j)$ in a certain set of edges, and aim to construct a positive semi-definite completion of $A$. This problem is extensively studied owing to its widespread applications in diverse fields such as probability, statistics, systems engineering and geophysics; see, for example, [Laurent, 2009] and the references therein for an introduction to the topic. Statistical questions associated with such problems are relatively under-explored, though we mention the recent work [Waghmare and Panaretos, 2022] that provides estimated completions of covariance operators in settings where completions always exist. A distinct but related problem that has received more attention in the statistics literature is low-rank matrix completion. In particular, significant contributions [Candès and Recht, 2009, Candès and Tao, 2010, Recht, 2011] have been made in the realm of convex optimization, where a low-rank matrix is recovered from partial observations after introducing a nuclear norm penalty. In our work we make no low-rank assumptions and our main interest is in answering the question of whether or not a positive semi-definite completion exists.

We now briefly outline our main contributions. In Section <ref> we define a numerical measure $R(\Sigma_\mathbb{S})$ of the incompatibility of a sequence of correlation matrices $\Sigma_\mathbb{S}$ and establish some key properties, including an interpretable and useful dual representation (Proposition <ref>). We then combine this index of incompatibility with a measure $V(\sigma^2_\mathbb{S})$ of the inconsistency of variances (e.g. Proposition <ref>), in order to introduce a numerical measure of the incompatibility of sequences of covariance matrices. In Section <ref> we turn to the empirical estimation of such indices and the introduction of data-driven testing procedures. Under the non-singularity assumption that $\Sigma_S \succeq c I_S$ for all $S \in \mathbb{S}$, we first introduce an oracle test that relies on knowledge of $c>0$ and give a result on its validity and power (Theorem <ref>). Next, we avoid this restriction by introducing a sample-splitting test that we prove to be valid (Proposition <ref>). The analysis of both tests is based on a novel concentration inequality for the spectral norm of the difference between the Pearson sample correlation matrix and its population version (Proposition <ref>). In Section <ref> we study the performance of our oracle test in various examples and show that its separation rate is near-minimax optimal in some cases, while studying properties of the associated semi-definite programs. In Section <ref> we validate our methodology in numerical experiments. Section <ref> contains the proofs of our main results. The Appendix contains background and auxiliary results.

We conclude the introduction with some notation that is used throughout the paper. For $d \in \mathbb{N}$, we write $[d] := \{1, \cdots, d\}$. Given $a, b \geq 0$, we write $a \lesssim b$
to mean that there exists a universal constant $C > 0$ such that $a \leq Cb$.
We use $a \wedge b$ for $\min\{a,b\}$, and $a \vee b$ for $\max\{a,b \}$. We will denote with $\boldsymbol{0}_d$ the null vector of dimension $d$, with $\boldsymbol{1}_d$ the all-one vector, with $\boldsymbol{e}_j$ the $j$-th element of the canonical basis of $\mathbb{R}^d$, with $\boldsymbol{O}_d$ the zero matrix of dimension $d$, and with $I_d$ the identity matrix of dimension $d$. We will omit the subscript with the dimension $d$ when it is clear from the context. For symmetric matrices $A, B$ of dimension $d$, we write $A \succeq 0$ to mean that $A$ is positive semi-definite, write $A \succeq B$ to mean that $A - B \succeq 0$, write $\operatorname{diag}(A)$ to indicate the diagonal matrix having the same diagonal as $A$, and $\operatorname{diag}(\boldsymbol{v})$ for a vector $\boldsymbol{v}=(v_1,\ldots,v_d)$ to indicate a diagonal matrix with diagonal elements equal to $v_i$. We will indicate the trace of $A$ with $\operatorname{tr}(A)$, the determinant with either $|A|$ or $\operatorname{det}(A)$, and the minimum and maximum eigenvalues of $A$ with $\lambda_{\mathrm{max}}(A)$ and $\lambda_{\mathrm{min}}(A)$, respectively. We will use  $\|\cdot\|_*$ for nuclear norm, or Schatten-1 norm of a matrix, $\|\cdot\|_2$ for the spectral norm, and $\|\cdot\|_F$ for the Frobenius norm. For random elements $X, Y$, we write $X \perp \!\!\! \perp Y$ to mean that $X$ and $Y$ are independent. For $\sigma>0$, a random variable $X$ with mean $\mu=\mathbb{E}[X]$ is said to be $\sigma$-subgaussian if

\[
\mathbb{E}\left[e^{\lambda(X-\mu)}\right] \leq e^{\sigma^2 \lambda^2 / 2} \quad \text { for all } \lambda \in \mathbb{R},
\]

while, for $(\nu,\alpha) \in (0,\infty)^2$, it is said to be $(\nu,\alpha)$-subexponential if

\mathbb{E}\left[e^{\lambda(X-\mu)}\right] \leq e^{\frac{\nu^2 \lambda^2}{2}} \quad \text { for all }|\lambda|<\frac{1}{\alpha}.

A random vector $\mathbf{X}$ in $\mathbb{R}^n$ is said to be $\sigma$-subgaussian if every one-dimensional projection, i.e. $v^T\mathbf{X}$ with $v \in \mathbb{R}^n$ and $\|v\|=1$, is $\sigma$-subgaussian in the sense defined above.

§ MEASURE OF INCOMPATIBILITY FOR COVARIANCE MATRICES

In this section we develop our index of the incompatibility of population covariance matrices. We check the compatibility of covariance matrices by checking the consistency of the variances and the compatibility of the correlation matrices separately. Standardising the covariance matrices is necessary to have a well-posed problem, and we choose to work with correlation matrices because the resulting compatibility measure is more tractable. Other standardisation are possible, though, each leading to a different measure of incompatibility, with different properties. In Appendix <ref> we introduce another measure of incompatibility based on a different standardisation, analyse its properties, and derive a test of MCAR based on its estimation from data.
Here, we standardise the sequence of covariance matrices, resulting in a sequence of correlation matrices $\Sigma_\mathbb{S}$, introduce a measure of incompatibility for $\Sigma_\mathbb{S}$, and develop some of its basic properties. We then define a measure of the inconsistency of variances, and combine the two for an overall measure of the incompatibility of sequences of covariance matrices. In order to do this, we must first introduce some basic algebraic objects for sequences of symmetric and positive semi-definite matrices. Our key notation is collected in Table <ref> below.

Notation     Definition     Meaning

$\mathbb{S}$     A subset of the power set of $[d]$     Set of all missingness patterns

$\mathbb{S}_j$     $\{S \in \mathbb{S} : j \in S\}$     Set of all patterns that contain $j$

$\mathbb{S}_{jj'}$     $\{S \in \mathbb{S} : j, j' \in S\}$     Set of all patterns that contain $(j,j')$

$\mathcal{M} \equiv \mathcal{M}_d$      $\{ X \in \mathbb{R}^{d \times d} : X = X^T\}$      Space of symmetric matrices

$\mathcal{M}_\mathbb{S}$     $\{ (X_S : S \in \mathbb{S}) : X_S \in \mathcal{M}_{|S|} \text{ for all } S \in \mathbb{S}\}$     Space of sequences of symmetric matrices

$\langle X,Y \rangle$      $ \mathrm{tr}(X Y)$ for $X,Y \in \mathcal{M}$     Frobenius inner product

$\langle X_\mathbb{S}, Y_\mathbb{S} \rangle_\mathbb{S}$     $\sum_{S \in \mathbb{S}} \mathrm{tr}(X_S Y_S)$ for $X_\mathbb{S}, Y_\mathbb{S} \in \mathcal{M}_\mathbb{S}$     Sum of Frobenius inner products

$\mathcal{P}^*$     $\{ \Sigma \in \mathcal{M} : \Sigma \succeq 0 \}$     Cone of PSD matrices

$\mathcal{P}$     $ \{ \Sigma \in \mathcal{P}^* : \mathrm{diag}(\Sigma)=I_d \}$     Set of correlation matrices

$\Sigma_\mathbb{S} \succeq_\mathbb{S} 0$     $\Sigma_\mathbb{S} \succeq_\mathbb{S} 0$ if and only if $\Sigma_S \succeq 0$ for all $S \in \mathbb{S}$     Loewner order for sequences of matrices

$\mathcal{P}_\mathbb{S}^*$     $\{ \Sigma_\mathbb{S} \in \mathcal{M}_\mathbb{S} : \Sigma_\mathbb{S} \succeq_\mathbb{S} 0 \}$     Sequences of PSD matrices

$\mathcal{P}_\mathbb{S}$     $\{ \Sigma_\mathbb{S} \in \mathcal{P}_\mathbb{S}^* : \mathrm{diag}(\Sigma_S)=I_{|S|} \, \text{for all } S \in \mathbb{S} \}$     Sequence of correlation matrices

$A$     $A : \mathcal{M} \rightarrow \mathcal{M}_\mathbb{S}$ with $ (A X)_S = (X_{jj'})_{j,j' \in S}$     Marginalisation operator on matrices

$\mathcal{P}_\mathbb{S}^{0,*}$     $ \{ A \Sigma : \Sigma \in \mathcal{P}^*\}$      Set of compatible sequences of PSD matrices

$\mathcal{P}_\mathbb{S}^{0}$     $\{ A \Sigma : \Sigma \in \mathcal{P}\}$      Set of compatible sequences of correlation matrices

$\mathcal{Y}$     $\{\operatorname{diag}(v) : v \in \mathbb{R}^d \text{ and } \sum_{j = 1}^d v_j = 0\}$     Space of diagonal matrices with null trace

$\boldsymbol{O}_\mathbb{S}$     $(O_{|S|} : S \in \mathbb{S})$     Sequence of zero matrices

$I_{\mathbb{S}}$      $(I_{|S|} : S \in \mathbb{S})$     Sequence of identity matrices

Table of definitions commonly used in the main text.

Crucially, we say that an element of $\mathcal{P}_\mathbb{S}^*$ is compatible if and only if it is an element of $\mathcal{P}_{\mathbb{S}}^{0,*}$.

In order to characterise compatibility, we first state a basic property of the linear operator $A$ defined in Table <ref> above.

The adjoint operator $A^* : \mathcal{M}_\mathbb{S} \rightarrow \mathcal{M}$ of $A$ is given by

\[
(A^* X_\mathbb{S})_{jj'} = \sum_{S \in \mathbb{S}} \mathbbm{1}_{\{j,j' \in S\}} (X_S)_{jj'} = \sum_{S \in \mathbb{S}_{jj'}} (X_S)_{jj'},
\]

where we recall that $\mathbb{S}_{jj'} = \{S \in \mathbb{S} : j,j' \in S\} = \mathbb{S}_j \cap \mathbb{S}_{j'}$.

Now, the following proposition fully characterises compatibility in terms of the non-negativity of a collection of linear functionals.

For $\Sigma_\mathbb{S} \in \mathcal{P}_\mathbb{S}^*$ we have $\Sigma_\mathbb{S} \in \mathcal{P}_\mathbb{S}^{0,*}$ if and only if

\[
\langle X_\mathbb{S}, \Sigma_\mathbb{S} \rangle_\mathbb{S} \geq 0 \quad \text{for all } X_\mathbb{S} \in \mathcal{M}_\mathbb{S} \text{ satisfying } A^* X_\mathbb{S} \succeq 0.
\]

The proof can be found in Section <ref>. This is an extension of well-known characterisation of the feasibility of positive semi-definite matrix completion <cit.>. Indeed, when $\Sigma_\mathbb{S}$ is consistent, we can show that $\langle X_\mathbb{S}, \Sigma_\mathbb{S} \rangle_\mathbb{S} = \langle A^* X_\mathbb{S}, \Sigma \rangle$, where $\Sigma$ is the $d \times d$ symmetric matrix with $\Sigma_{jj'}=(\Sigma_S)_{jj'}$ for all $S \in \mathbb{S}_{jj'}$, if $\mathbb{S}_{jj'} \neq \emptyset$, and $\Sigma_{jj'}=0$ if $\mathbb{S}_{jj'} = \emptyset$. Here $\Sigma$ can be thought of as a partial matrix that is padded with zeros in unobserved positions. Since $A^* X_\mathbb{S}$ is also zero in these positions, the value of $\Sigma$ there is arbitrary. Now our characterisation reduces to checking that $\langle A^* X_\mathbb{S}, \Sigma \rangle \geq 0$ for all $X_\mathbb{S} \in \mathcal{M}_\mathbb{S}$ satisfying $A^* X_\mathbb{S} \succeq 0$, which is equivalent to checking $\langle X, \Sigma \rangle \geq 0$ for all $X \in \mathcal{M}$ with $X_{jj'}=0$ if $S_{jj'} = \emptyset$, which coincides with (4) in [Laurent, 2009].

Proposition <ref> provides a characterisation of compatibility, but in order to assess the significance of departures from the null hypothesis and thus to define hypothesis tests, we will need a numerical measure of incompatibility. A natural way to do this is to minimise $\langle X_\mathbb{S}, \Sigma_\mathbb{S} \rangle_\mathbb{S}$ over a subset of $\{ X_\mathbb{S} \in \mathcal{M}_\mathbb{S} : A^* X_\mathbb{S} \succeq 0\}$ that still characterises compatibility, but gives finite minimal values. First, observe that checking the compatibility of covariance matrices is equivalent to checking the consistency of the variances $\sigma^2_S = \operatorname{diag}(\mathrm{Cov}(\mathbf{X}_S|\Omega=\mathbbm{1}_S)$ for $S \in \mathbb{S}$, and the compatibility of the correlation matrices $\operatorname{Corr}(\mathbf{X}_S|\Omega = \mathbbm{1}_S)$ for $S \in \mathbb{S}$.
Now, whenever $\Sigma_\mathbb{S}$ is a sequence of correlation matrices we define

\begin{align}
\label{Eq:Primal}
R(\Sigma_\mathbb{S}) &:= \sup\biggr\{-\frac{1}{d} \langle \Sigma_\mathbb{S}, X_\mathbb{S} \rangle_\mathbb{S} : X_\mathbb{S}+X_\mathbb{S}^0 \succeq_\mathbb{S} 0,  A^* X_\mathbb{S} + Y \succeq 0 \text{ for some } Y \in \mathcal{Y} \biggr\} \nonumber \\
&= 1- \frac{1}{d}\inf\{ \langle \Sigma_\mathbb{S}, Y_\mathbb{S} \rangle_\mathbb{S} : Y_\mathbb{S} \succeq_\mathbb{S} 0,  A^* Y_\mathbb{S} + Y \succeq I_d \text{ for some } Y \in \mathcal{Y} \}
\end{align}

where $X_\mathbb{S}^0 = (X_S^0 : S \in \mathbb{S})$, where $X_S^0 = \operatorname{diag}\left(1 /\left|\mathbb{S}_j\right|: j \in S\right)$, and where $\mathcal{Y}$ is the set of diagonal $d \times d$ matrices with trace equal to zero. The objective function of this optimisation problem is a one-to-one mapping of the linear functional appearing in our characterisation of compatibility. Moreover, by choosing $Y=\boldsymbol{O}$ and noting that $X_\mathbb{S}^0 \succ_\mathbb{S} 0$, we can see that for any $X_\mathbb{S}$ that satisfies $A^* X_\mathbb{S} \succeq 0$, the sequence $\epsilon X_\mathbb{S}$ is feasible for $\epsilon>0$ sufficiently small. Thus, by Proposition <ref> we have that $R(\Sigma_\mathbb{S})>0$ whenever $\Sigma_\mathbb{S}$ is incompatible. On the other hand, when $\Sigma_\mathbb{S} = A\Sigma$ is compatible and $(X_\mathbb{S},Y)$ is feasible we have $\langle \Sigma_\mathbb{S}, X_\mathbb{S} \rangle_\mathbb{S} = \langle \Sigma, A^* X_\mathbb{S} \rangle =  \langle \Sigma, A^* X_\mathbb{S} + Y \rangle \geq 0$, where the second equality holds because $\Sigma$ has a constant diagonal. Combining this with the observation that $X_\mathbb{S} = \boldsymbol{O}_\mathbb{S}$ is feasible, we see that $R(\Sigma_\mathbb{S})=0$ when $\Sigma_\mathbb{S}$ is compatible.

In the above argument we did not use the specific form of the lower bound $X_\mathbb{S} \succeq_\mathbb{S} -X_\mathbb{S}^0$ anywhere, and it would also have been possible to optimise over the restricted set of $X_\mathbb{S}$ that are feasible with $Y=\boldsymbol{O}$. The specific choice of the feasible set in the definition of $R(\Sigma_\mathbb{S})$ was made because it leads to an interpretable dual formulation. While there exist semi-definite programs for which strong duality does not hold, Slater's condition (see Appendix <ref> for an introduction to the theory of semi-definite programming) is satisfied in our problem, so we do not encounter such issues. This is formalised in the result below.

For $\Sigma_\mathbb{S} \in \mathcal{P}_\mathbb{S}$ we have

\begin{align}
\label{Eq:Dual}
R(\Sigma_{\mathbb{S}}) &=\inf \{ \epsilon \in [0,1] : \Sigma_\mathbb{S} \in (1-\epsilon) \mathcal{P}_\mathbb{S}^0 + \epsilon \mathcal{P}_\mathbb{S} \} \nonumber \\
&= 1-\frac{1}{d} \sup\{ \mathrm{tr}(\Sigma) : A\Sigma \preceq_\mathbb{S} \Sigma_\mathbb{S}, \, \Sigma_{11}=\ldots=\Sigma_{dd}, \, \Sigma \succeq 0 \}.
\end{align}

This result shows that our measure of incompatibility $R(\Sigma_\mathbb{S})$

can be interpreted as the smallest amount of perturbation $\epsilon$ that a compatible sequence of correlation matrices must be corrupted by to result in the input sequence $\Sigma_{\mathbb{S}}$. It is immediate from this representation that $R(\Sigma_\mathbb{S})$ takes values in $[0,1]$. Moreover, it follows from Slater's condition that the optimal value of the dual problem is attained. Thus, writing $\lambda^* = 1- R(\Sigma_{\mathbb{S}})$, we can write $\Sigma_{\mathbb{S}}$ as

\begin{align}\label{eq:dual_decomposition}
\Sigma_{\mathbb{S}} = \lambda^*A\Sigma + (1-\lambda^*)\Sigma_{\mathbb{S}}',
\end{align}

where $\Sigma \in \mathcal{P}$ and $\Sigma_{\mathbb{S}}' \in \mathcal{P}_\mathbb{S}$. By the maximality of $\lambda^* = 1-R(\Sigma_\mathbb{S})$, it must be the case that $R(\Sigma_{\mathbb{S}}')=1$.
Indeed, if this were not the case, it would be possible to write $\Sigma_{\mathbb{S}}' = \lambda'A\Sigma' + (1-\lambda')\Sigma_{\mathbb{S}}''$ for some $\lambda' \in (0,1]$, which would contradict the fact that $\lambda^*$ is optimal.

This argument shows that, whenever $\mathbb{S}$ is such that there exists an incompatible sequence $\Sigma_\mathbb{S}$, the maximal value $R(\Sigma_\mathbb{S})=1$ is attainable, so that the quantity $R(\Sigma_\mathbb{S})$ is on an interpretable scale between compatibility at one extreme and maximal incompatibility at the other.

We remark that this dual interpretation of $R(\cdot)$ aligns with a similar representation of the incompatiblity of sequences of distributions defined by [Berrett and Samworth, 2023]. In this earlier work it is shown that the incompatibility of sequences of distributions can be understood through linear programming techniques. Our work here, however, shows that we must consider the more complex problem of semi-definite programming to understand the incompatibility of sequences of covariance matrices. Despite this additional complexity, since Slater's condition is satisfied for our problem, the primal-dual interior point method has a computational complexity which is polynomial in the number of constraints and the dimension of the unknown square matrix (Section 6.4.1. of [Nesterov and Nemirovskii, 1994], Section 5.7. of [Vandenberghe and Boyd, 1996]). This ensures that $R(\Sigma_\mathbb{S})$ can be always computed efficiently without additional assumptions.

We conclude with some basic properties of $R(\cdot)$.

The following hold:

(i) $R(\cdot)$ is convex.

(ii) $R(\cdot)$ is continuous.

(iii) Suppose $\mathbb{S} \subseteq \mathbb{S}'$ and $\Sigma_{\mathbb{S}} \subseteq_{\mathbb{S}} \Sigma_{\mathbb{S}'}$, where the inclusion $\subseteq_{\mathbb{S}}$ means that every correlation matrix in $\Sigma_{\mathbb{S}}$ is also in $\Sigma_{\mathbb{S}'}$. Then $R(\Sigma_{\mathbb{S}}) \leq R(\Sigma_{\mathbb{S}'})$.

It is interesting to observe that property (iii) says that $R$ is monotone with respect to the inclusion operator, so that additional information can only make a sequence less compatible.

Having established the main properties of our measures of the incompatiblity of correlation matrices, we now turn to the simpler problem of measuring the inconsistency of individual variances.

There could be cases for which variances are not consistent across different patterns, but a test based on $R(\cdot)$ alone would fail to reject the null hypothesis. In order to take into account such deviations from the null, we define an analogous test statistic for the consistency of the variances. To this aim, writing $\sigma_\mathbb{S}^2=(\sigma_S^2 : S \in \mathbb{S})$ for the collection of individual variances, we fix our units of measurement such that

\[
\bar{\operatorname{av}}_j(\sigma^2_{\mathbb{S}}) :=  |\mathbb{S}_j|^{-1} \sum_{S \in \mathbb{S}_j} \sigma_{S,j}^2 = 1,
\]

for all $j \in [d]$. This is a natural constraint, analogous to the standardisation of variables in complete-data problems, that does not remove information that may be present in the individual variances. For such $\sigma_\mathbb{S}^2$, define

\[
V(\sigma_{\mathbb{S}}^2) := 1 - \min_{j \in [d]}\min_{S \in \mathbb{S}_j} \sigma_{S,j}^2 = \max_{j \in [d]} \max_{S \in \mathbb{S}_j} \left(1 - \sigma_{S,j}^2\right).
\]

Under the hypothesis $\bar{\operatorname{av}}_j(\sigma^2_{\mathbb{S}}) = 1$ for all $j \in [d]$, it is clear that $V(\sigma_{\mathbb{S}}^2) = 0$ if and only if $\sigma^2_{S,j} = 1$ for all $j \in [d]$, $S \in \mathbb{S}_j$. On the other hand, we have $V(\sigma_{\mathbb{S}}^2) > 0$, if and only if there exists at least one variance strictly less than $1$. It is clear from the definition that $V$ is bounded by one, and that this extreme value is attainable when $\mathbb{S}$ is non-trivial and there exists $j$ such $\sigma_{S,j}^2=0$ for some $S \in \mathbb{S}_j$. The following result gives a dual representation for $V(\sigma^2_\mathbb{S})$, providing justification for our specific measure of inconsistency.

If $\bar{\operatorname{av}}_j(\sigma^2_{\mathbb{S}}) = 1$ for all $j \in [d]$, then

\[
V(\sigma_{\mathbb{S}}^2) = \inf\left\{ \epsilon \in [0,1] : \sigma_{\mathbb{S}}^2 = (1-\epsilon)A_{V} \boldsymbol{1}_d + \epsilon {\sigma'}_{\mathbb{S}}^2 \text{ with } \bar{\operatorname{av}}_j(\sigma'^2_{\mathbb{S}}) = 1 \text{ for all } j \in [d] \right\},
\]

where $(A_V \sigma^2)_S = (\sigma^2_k)_{k \in S}$.

This result gives a dual representation for $V(\sigma^2_\mathbb{S})$, which is analogous to Proposition <ref> and leads to a similar interpretation, based on the idea of finding the smallest perturbation to make the sequence consistent.

Combining our measures of the incompatibility of correlation matrices and the inconsistency of variances, we define an overall measure of the incompatibility of sequences of covariance matrices by
\[
T =  R(\Sigma_\mathbb{S}) + V(\sigma_{\mathbb{S}}^2).
\]

It is easy to see that $T=0$ if and only if the sequence of correlation matrices $\Sigma_\mathbb{S}$ is compatible and the sequence of variances $\sigma_{\mathbb{S}}^2$ is consistent, and that $T \leq 2$.

§ STATISTICAL TESTS

Having introduced our population-level measure of incompatibility, in this section we design a testing procedure based on the estimation of $T$ to test compatibility with finite amounts of data. For each $S \in \mathbb{S}$ we assume that we have access to an independent sample

\[
X_{S,1},\ldots,X_{S,n_S} \overset{\text{i.i.d.}}{\sim} P_S
\]

for some sample size $n_S$ and some distribution $P_S$ with correlation matrix $\Sigma_S$ and vector of variances $\sigma_S^2$. We write $\hat{\Sigma}_S$ and $\hat{\sigma}_S^2$ for the sample correlation matrix and sample variances, respectively, of $X_{S,1},\ldots,X_{S,n_S}$ for each $S \in \mathbb{S}$, and write $\hat{\Sigma}_\mathbb{S} = (\hat{\Sigma}_S : S \in \mathbb{S})$ and $\hat{\sigma}_\mathbb{S}^2 = (\hat{\sigma}_S^2 : S \in \mathbb{S})$ for the sequences collecting these estimators. Our approach is based on analysing the concentration properties of the plug-in estimates  $ R(\hat{\Sigma}_\mathbb{S})$ and $V(\hat{\sigma}_{\mathbb{S}}^2)$ and thus introducing suitable critical values.

The analysis of $R(\hat{\Sigma}_\mathbb{S})$, in particular, is challenging, as it is defined as the optimal value of a semi-definite program with an unbounded feasible set. In fact, without further assumptions, it is not possible to restrict attention to a compact feasible set. On the other hand, most statistical techniques for the analysis of suprema of empirical processes require feasible sets to be totally bounded so that, for example, covering arguments can be applied.

Fortunately, under the assumption that $\Sigma_\mathbb{S} \succ_\mathbb{S} 0$, our dual problem (<ref>) is strictly feasible and hence Slater's condition implies that the optimal value is attained in the primal problem (<ref>). This assumption is reasonable in many areas of application, and similar assumptions of invertibility have been used frequently in the literature [Meinshausen and Bühlmann, 2006, Cai et al., 2011]. In fact, if we assume the stronger condition that $\Sigma_\mathbb{S} \succeq_\mathbb{S} c I_\mathbb{S}$ for some $c>0$, we will see that the optimal value is always attained in a compact set whose size depends on $c$.

Indeed, the strict feasibility of the dual problem (<ref>) implies that there exists $X_{\mathbb{S}} \in \mathcal{M}_\mathbb{S}$ such that $X_\mathbb{S} + X_\mathbb{S}^0 \succeq_\mathbb{S} 0$ and $R(\Sigma_{\mathbb{S}}) = -d^{-1} \langle X_{\mathbb{S}}, \Sigma_{\mathbb{S}} \rangle_{\mathbb{S}}$. This in turn implies that

\[
\sum_{S \in \mathbb{S}}\|X_S + X_S^0\|_*
= \langle X_{\mathbb{S}} + X_{\mathbb{S}}^0, I_{\mathbb{S}} \rangle_{\mathbb{S}} \leq \frac{1}{c} \langle X_{\mathbb{S}} + X_{\mathbb{S}}^0, \Sigma_{\mathbb{S}} \rangle_{\mathbb{S}} = \frac{d}{c}\{1-R(\Sigma_\mathbb{S})\} \leq \frac{d}{c},
\]

so that we have a bound on the sum of the nuclear norms of the matrices in the sequence $X_\mathbb{S} + X_\mathbb{S}^0$. In finding the optimal value of the primal problem (<ref>), then, we may restrict attention to

\begin{equation}
\label{Eq:CompactFeasibleSet}
\mathcal{F}_c := \biggl\{ X_\mathbb{S} \in \mathcal{M}_\mathbb{S} :  X_{\mathbb{S}} + X_{\mathbb{S}}^0 \succeq_{\mathbb{S}} 0 , \sum_{S \in \mathbb{S}}\|X_S + X_S^0\|_* \leq d/c,  A^* X_\mathbb{S} + Y \succeq 0 \text{ for some } Y \in \mathcal{Y} \biggr\},
\end{equation}

which is compact.

Before moving on to describe how to construct a statistical test under this new assumption, we give a brief discussion of the norm on $\mathcal{M}_\mathbb{S}$ defined by

\[
\|X_{\mathbb{S}}\|_{*, \mathbb{S}} := \sum_{S \in \mathbb{S}}\|X_S\|_*,
\]

which reduces to $\sum_{S \in \mathbb{S}} \operatorname{tr}(X_S)$ in case that $X_{\mathbb{S}} \succeq_{\mathbb{S}} 0$. For each $S \in \mathbb{S}$, the nuclear norm $\|X_S\|_*$ can be thought of as the $\ell_1$ norm applied to the eigenvalues of $X_S$. As these are then summed to give $\|X_{\mathbb{S}}\|_{*, \mathbb{S}}$, we see that $\| \cdot \|_{*, \mathbb{S}}$ can be thought of as an $\ell_1$ norm on $\mathcal{M}_\mathbb{S}$. It is easy to see that the dual norm of $\|\cdot\|_{*, \mathbb{S}}$ with respect to the inner product $\langle \cdot, \cdot \rangle_{\mathbb{S}}$ is

\[
\|X_{\mathbb{S}}\|_{2, \mathbb{S}} := \max_{S \in \mathbb{S}}\|X_S\|_2,
\]

where $\|\cdot\|_2$ is the usual spectral norm of a matrix. This follows after writing the sequence of matrices in block-diagonal form, and allows us to derive the following generalisation of Holder's inequality in the space of matrix sequences,

\begin{equation}
\label{Eq:Holder}
|\langle X_{\mathbb{S}}, Y_{\mathbb{S}} \rangle_{\mathbb{S}}| \leq \|X_{\mathbb{S}}\|_{*, \mathbb{S}} \|Y_{\mathbb{S}}\|_{2, \mathbb{S}} = \sum_{S \in \mathbb{S}}\|Y_S\|_* \cdot \max_{S \in \mathbb{S}}\|X_S\|_2.
\end{equation}

This inequality will be used in the proof of the following result, which provides valid critical values for the test statistic $\hat{T} = R(\hat{\Sigma}_{\mathbb{S}}) + V(\hat{\sigma}_{\mathbb{S}}^2)$ and gives conditions under which the resulting test has large power.

Suppose we observe $\mathbf{X}_{S,1}, \ldots, \mathbf{X}_{S, n_S} \overset{\text{i.i.d.}}{\sim} P_S$ for each $S \in \mathbb{S}$ independently, where each $P_S$ is $\nu$-subgaussian, with the sequence of variances $\sigma^2_\mathbb{S}$ satisfying  $\bar{\operatorname{av}}_j(\sigma^2_{\mathbb{S}}) = 1$ for all $j \in [d]$, and the sequence of population correlation matrices $\Sigma_{\mathbb{S}}$ satisfying $\Sigma_{\mathbb{S}} \succeq_\mathbb{S} cI_\mathbb{S}$ for $c >0$.

Define $C_{\alpha} = \max\{C_{\alpha}^{(R)} , C_{\alpha}^{(V)}\}$, with

\[
C_{\alpha}^{(V)} = C_0 \nu^2 \sqrt{\frac{\log\left( \sum_{j \in [d]} |\mathbb{S}_j|/\alpha\right)}{\min_{S \in \mathbb{S}} n_S}},
\]

\begin{align*}
C_{\alpha}^{(R)} = \frac{1}{c}&\max_{S \in \mathbb{S}}\Biggr\{ C_1 \frac{\nu^2}{\sigma_\mathrm{min}^2} \sqrt{\frac{|S| + \log(|\mathbb{S}|/\alpha)}{n_S}} \vee \frac{|S| + \log(|\mathbb{S}|/\alpha)}{n_S} + C_2 \frac{\nu^4}{\sigma^4_{\mathrm{min}}} \sqrt{\frac{|S| \log(|S\|\mathbb{S}|/\alpha)}{n_S}} \\
& \quad \quad \quad \quad  + C_3 \frac{\nu^4}{\sigma^4_{\mathrm{min}}} \left( \sqrt{\frac{|S| + \log(|\mathbb{S}|/\alpha)}{n_S}} \vee \frac{|S| + \log(|\mathbb{S}|/\alpha)}{n_S} \right) \sqrt{\frac{|S| \log(|S|/\alpha)}{n_S}}\Biggr\},
\end{align*}

for universal constants $C_0, C_1, C_2, C_3 >0$. Then, for $C_0, C_1,C_2,C_3$ and universal constant $K>0$ chosen sufficiently large, for all $\alpha \in (0,1)$ such that $\sigma^4_{\mathrm{min}} n \geq K \nu^4 \log(\max_{S \in \mathbb{S}}|S|/\alpha)$, the test that rejects $H_0$ if and only if $\hat{T} = R(\hat{\Sigma}_{\mathbb{S}}) + V(\hat{\sigma}_{\mathbb{S}}^2) \geq C_{\alpha}$ has Type I error bounded by $\alpha$. Moreover, for $\beta \in (0,1)$ such that $\sigma^4_{\mathrm{min}} n \geq K \nu^4 \log(\max_{S \in \mathbb{S}}|S|/\beta)$, if $T > C_{\alpha} + C_{\beta}$, then $\mathbb{P}(\hat{T} \leq C_{\alpha}) \leq \beta.$

In proving this result we give concentration inequalities for the random quantities $R(\hat{\Sigma}_\mathbb{S})$ and $V(\hat{\sigma}_\mathbb{S}^2)$. The analysis of $R(\hat{\Sigma}_\mathbb{S})$ is crucially based on the fact that, under $H_0$ and in light of the inequality (<ref>), we can control the oscillation $|R(\hat{\Sigma}_\mathbb{S}) - R(\Sigma_\mathbb{S})|$ using $\max_{S \in \mathbb{S}} \|\hat{\Sigma}_{S} - \Sigma_{S}\|_2$, where the $\Sigma_{S}$ are the Pearson population correlation matrices and $\hat{\Sigma}_{S}$ are the corresponding Pearson sample correlation matrices. To this end, we derive the following tail bound for the spectral norm $\|\hat{P} - P\|_2$, where $P$ is the population correlation matrix and $\hat{P}$ is the sample correlation matrix of complete data, which may be of independent interest.

Suppose we observe an i.i.d sample $\mathbf{X}_1, \cdots, \mathbf{X}_n \sim \mathbf{X}$, where $\mathbf{X}$ is a $\nu$-subgaussian random vector in $\mathbb{R}^d$, and let $\Sigma$ and $\hat{\Sigma}$ be the population and sample covariance matrices, respectively. Let $P = D^{-1/2}\Sigma D^{-1/2}$ be the population correlation matrix, where $D = \operatorname{diag}(\Sigma)$, and $\hat{P} = \hat{D}^{-1/2} \hat{\Sigma} \hat{D}^{-1/2}$ be the sample correlation matrix, where $\hat{D} = \operatorname{diag}(\hat{\Sigma})$. Then, there exist universal constants $ C_1, C_2, C_3, K > 0$ such that, for every $t \in [0,1]$ such that $\sigma^4_{\mathrm{min}} n \geq K \nu^4 \log(d/t)$ we have

\begin{align*}
\|\hat{P} - P \|_2 & \leq C_1 \frac{\nu^2}{\sigma_\mathrm{min}^2} \left( \sqrt{\frac{d + \log(1/t)}{n}} \vee \frac{d + \log(1/t)}{n}\right) + C_2 \frac{\nu^4}{\sigma^4_{\mathrm{min}}} \sqrt{\frac{d \log(d/t)}{n}}  \\
& \quad \quad + C_3 \frac{\nu^4}{\sigma^4_{\mathrm{min}}} \left( \sqrt{\frac{d + \log(1/t)}{n}} \vee \frac{d + \log(1/t)}{n} \right) \sqrt{\frac{d \log(d/t)}{n}}
\end{align*}

with probability $\geq 1 - t$, where $\sigma_\mathrm{min}^2 := \min\limits_{j \in [d]} \Sigma_{jj}$.

First, observe that the dependence on $1/\sigma_\mathrm{min}^2$ is reasonable, as the smaller the minimum variance the more problematic the normalisation matrix $D^{-1/2}$. Fortunately, since we work under the assumption that $\Sigma_{\mathbb{S}} \succeq_{\mathbb{S}} cI_\mathbb{S}$, we have that $\sigma_\mathrm{min}^2 \geq c$ when we apply this result. Second, observe that if $t > 2\exp\{ -n + 2d \log3\}$, i.e. $n \gtrsim d$, the subgaussian regime prevails, and we obtain that

\[
\|\hat{P} - P\|_2 \lesssim \sqrt{\frac{d \log d}{n}}
\]

in probability. Similar rates, with logarithmic factors, were found in high-dimensional covariance matrix estimation with missing observations [Lounici, 2014], sample covariance matrix estimator of reduced effective rank population matrices [Bunea and Xiao, 2015], concentration of the adjacency matrix and of the Laplacian in random graphs [Oliveira, 2010], and in the statistical analysis of latent generalized correlation matrix estimation in transelliptical distribution [Han and Liu, 2017]. In particular, using the additional assumption that the data is generated according to a transelliptical distribution, [Han and Liu, 2017] gave an estimator $\hat{K}$ based on Kendall's tau and proved that

\[
\|\hat{K} - P\|_2 \lesssim \sqrt{\frac{r(P) \log(d)}{n}},
\]

where $r(P) := \operatorname{tr}(P)/\|P\|_2$ is the effective dimension of $P$. This is analogous to the bound given in Proposition <ref>, where we have an extra factor of $\nu^2/\sigma_\mathrm{min}^2$, which can be interpreted as the condition number and might lead to a suboptimal bound when it is large, and the ambient dimension $d$ in place of the intrinsic dimension $r(P)$. This would improve the bound sensibly in the case of an approximately low-rank correlation matrix, but in the worst case the bounds have the same rates.

As well as providing a critical value for our test, Theorem <ref> also gives upper bounds on the minimax separation rate for this testing problem. When $c>0$ and $\nu>0$ are fixed, our analysis gives an upper bound on the minimax rate of the order

\[
2C_{1/4}\lesssim \max_{S \in \mathbb{S}}\sqrt{\frac{|S|\log(|S\|\mathbb{S}|)}{n_S}}.
\]

whenever $n_S \gtrsim |S|$ for all $S \in \mathbb{S}$. This is our main regime of interest, and we see in our examples in Section <ref> below that reliable testing is only possible when sample sizes are large compared with dimensions, up to logarithmic factors.

We now illustrate the behaviour of this bound in certain examples where the expression for $C_\alpha$ can be simplified. The corresponding upper bounds on the minimax separation rate will be complemented by lower bounds in Section <ref> to follow.

In the $d$-cycle example, with $d \geq 3$ and $\mathbb{S}= \{ \{1,2\} , \{2,3\} ,\ldots, \{d-1,d\}, \{d,1\}\}$, we have $|S|=2$ for all $S \in \mathbb{S}$ and $|\mathbb{S}| = d$ so that

\[
C_{\alpha} \lesssim  \max\limits_{S \in \mathbb{S}} \sqrt{\frac{\log(d/\alpha)}{n_S}} =   \sqrt{\frac{\log(d/\alpha)}{\min\limits_{S \in \mathbb{S}} n_S}}.
\]

If $n_S = n$ for all $S \in \mathbb{S}$, this reduces to

\[
C_{\alpha} \lesssim  \sqrt{\frac{\log(d/\alpha)}{n}},
\]

and by considering the sub-problem of testing the consistency of variances, we will see in Theorem <ref> that this upper bound is optimal in this case. Our results reveal that, in this specific example, testing compatibility is no harder than testing consistency, up to constant factors.

Consider the block-$3$-cycle $\mathbb{S} = \{[2d], [d]\cup([3d]\setminus[2d]),[3d]\setminus[d]\}$, with $d \geq 1$. Then

\[
C_{\alpha} \lesssim  \max\limits_{S \in \mathbb{S}} \sqrt{\frac{d\log(d/\alpha)}{n_S}} =  \sqrt{\frac{d\log(d/\alpha)}{\min\limits_{S \in \mathbb{S}}  n_S}}.
\]

As before, if $n_S = n$ for all $S \in \mathbb{S}$, this reduces to

\[
C_{\alpha} \lesssim  \sqrt{\frac{d\log(d/\alpha)}{n}}.
\]

We prove the minimax optimality, up to logarithmic factors, of this rate in Theorem <ref>. In particular, this shows that the optimal separation rates for this testing problem are not signficantly faster than the optimal rates for the estimation of $\Sigma_\mathbb{S}$ with operator norm loss.

Consider the example $\mathbb{S} = \{S_{(-1)}, \cdots, S_{(-d)}\}$, with $d \geq 2$, where $S_{(-i)} = \{1, \ldots, i-1, i+1, \ldots, d\}$. This corresponds to the setting where all observations have a single missing value. In this case, $|S_{(-i)}|=d-1$, and $|\mathbb{S}| = d$, so that

\[
C_{\alpha} \lesssim  \max\limits_{S \in \mathbb{S}} \sqrt{\frac{d\log(d/\alpha)}{n_S}} =  \sqrt{\frac{d\log(d/\alpha)}{\min\limits_{S \in \mathbb{S}} n_S}}.
\]

The test we developed in Theorem <ref> depends on the unknown quantity $c > 0$, which quantifies the level of positive-definiteness of $\Sigma_{\mathbb{S}}$ by requiring $\Sigma_{\mathbb{S}} \succeq_{\mathbb{S}} cI_{\mathbb{S}}$. As a result, the practical implementation of this method would be difficult. In the remainder of this section we relax this assumption using the well-known technique of sample splitting <cit.> to calculate a critical value that depends on the data. We randomly divide the data into two non-overlapping sets, $\mathcal{X}_1$ and $\mathcal{X}_2$, and then perform the so-called two steps of hunting and testing: sample $\mathcal{X}_1$ is used to select a test statistic from a large class of potential test statistics, which is then applied to $\mathcal{X}_2$ to produce
the final test statistic. This has the advantage of being easier to calibrate, due to the fact that $\mathcal{X}_1 \perp \! \! \! \perp \mathcal{X}_2$. This approach has been used for a variety of problems, such as testing the location of multiple samples [Cox, 1975], constructing conformal prediction intervals [Lei et al., 2018, Solari and Djordjilović, 2022], goodness-of-fit testing [Janková et al., 2020], conditional mean independence testing
[Scheidegger et al., 2022, Lundborg et al., 2022], and conducting inference that is agnostic to the asymptotic regime [Kim and Ramdas, 2024].

In our problem we use sample splitting as follows. Using the data in $\mathcal{X}_1$ we form the sample correlation matrices $\hat{\Sigma}_\mathbb{S}^{(1)}$ then calculate the value of our incompatibility index $R(\hat{\Sigma}_\mathbb{S}^{(1)})$ and retain the optimal primal solution $\hat{X}^{(1)}_{\mathbb{S}}$, which satisfies $R(\hat{\Sigma}_\mathbb{S}^{(1)}) = -d^{-1}\langle \hat{X}^{(1)}_{\mathbb{S}}, \hat{\Sigma}_{\mathbb{S}}^{(1)} \rangle_{\mathbb{S}}$. Now, using the data in $\mathcal{X}_2$ we form the sample correlation matrices $\hat{\Sigma}_\mathbb{S}^{(2)}$ and calculate the statistic

\[
\hat{R}_2 := - \frac{1}{d} \langle \hat{X}^{(1)}_{\mathbb{S}}, \hat{\Sigma}_{\mathbb{S}}^{(2)} \rangle_{\mathbb{S}}.
\]

We then take our final test statistic to be $V(\hat{\sigma}_{\mathbb{S}}^2) + \hat{R}_2$, where $\hat{\sigma}_\mathbb{S}^2$ is calculated using the entire sample. Analytically, the advantage of using this approach is that the statistic $\hat{R}_2$ is much easier to understand than the statistic $R(\hat{\Sigma}_\mathbb{S})$ used for the oracle test above. Indeed, the independence of $\hat{X}^{(1)}_{\mathbb{S}}$ and $ \hat{\Sigma}_{\mathbb{S}}^{(2)}$ means that it is enough to give concentration inequalities for linear functionals of $ \hat{\Sigma}_{\mathbb{S}}^{(2)}$, rather than for the supremum $R(\hat{\Sigma}_\mathbb{S})$ over a class of linear functionals. The following result introduces a critical value for this test statistic, which is proved to result in a valid test.

Suppose we observe $\mathbf{X}_{S,1}, \ldots, \mathbf{X}_{S, n_S} \overset{\text{i.i.d.}}{\sim} P_S$ for each $S \in \mathbb{S}$ independently, where each $P_S$ is $\nu$-subgaussian, with the sequence of variances $\sigma^2_\mathbb{S}$ satisfying  $\bar{\operatorname{av}}_j(\sigma^2_{\mathbb{S}}) = 1$ for all $j \in [d]$, and the sequence of population correlation matrices $\Sigma_{\mathbb{S}}$ satisfying $\Sigma_{\mathbb{S}} \succ_\mathbb{S} 0$.

Define $C_{\alpha}(\mathcal{X}_1) = \max\{C_{\alpha}^{(R)}(\mathcal{X}_1) , C_{\alpha}^{(V)}\}$, with

\[
C_{\alpha}^{(V)} = C_0 \nu^2 \sqrt{\frac{\log\left(\sum_{j \in [d]} |\mathbb{S}_j|/\alpha\right)}{\min_{S \in \mathbb{S}} n_S}},
\]

\begin{align*}
C_{\alpha}^{(R)}(\mathcal{X}_1) = \frac{\| \hat{X}^{(1)}_{\mathbb{S}}\|_{*, \mathbb{S}}}{d}&\max_{S \in \mathbb{S}}\Biggr\{ C_1 \frac{\nu^2}{\sigma_\mathrm{min}^2} \sqrt{\frac{|S| + \log(|\mathbb{S}|/\alpha)}{n_S}} \vee \frac{|S| + \log(|\mathbb{S}|/\alpha)}{n_S} + C_2 \frac{\nu^4}{\sigma^4_{\mathrm{min}}} \sqrt{\frac{|S| \log(|S| |\mathbb{S}|/\alpha)}{n_S}} \\
& \quad \quad \quad \quad  + C_3 \frac{\nu^4}{\sigma^4_{\mathrm{min}}} \left( \sqrt{\frac{|S| + \log(|\mathbb{S}|/\alpha)}{n_S}} \vee \frac{|S| + \log(|\mathbb{S}|/\alpha)}{n_S} \right) \sqrt{\frac{|S| \log(|S|/\alpha)}{n_S}}\Biggr\},
\end{align*}

for some universal constants $C_0, C_1, C_2, C_3 >0$. Then, for $C_0, C_1,C_2,C_3$ and universal constant $K>0$ chosen sufficiently large, for all $\alpha \in (0,1)$ such that $\sigma^4_{\mathrm{min}} n \geq K \nu^4 \log(\max_{S \in \mathbb{S}}|S|/\alpha)$, the test that rejects $H_0$ if and only if $V(\hat{\sigma}_{\mathbb{S}}^2) +\hat{R}_2 \geq C_{\alpha}(\mathcal{X}_1)$ has Type I error bounded by $\alpha$.

This test gives a random threshold $C_{\alpha}(\mathcal{X}_1)$ which does not depend on $c$, and it is of the order of

\[
C_{\alpha}(\mathcal{X}_1) \lesssim \frac{\| \hat{X}^{(1)}_{\mathbb{S}}\|_{*, \mathbb{S}}}{d} \max_{S \in \mathbb{S}}\sqrt{\frac{|S|\log(|S\|\mathbb{S}|/\alpha)}{n_S}}
\]

whenever $C_\alpha^{(R)} \geq C_\alpha^{(V)}$. Comparing this to the critical value used for the oracle test of Theorem <ref>, observe that $\| \hat{X}^{(1)}_{\mathbb{S}}\|_{*, \mathbb{S}}/d$ takes the place of $1/c$. In this regard, notice that for $X_{\mathbb{S}} \in \mathcal{F}_c$ we have
\[
\frac{\| X_{\mathbb{S}} + X_{\mathbb{S}}^0\|_{*, \mathbb{S}}}{d} \leq \frac{1}{c},
\]

by definition of $\mathcal{F}_c$. This provides some justification for thinking of $\| \hat{X}^{(1)}_{\mathbb{S}}\|_{*, \mathbb{S}}/d$ as an estimated lower bound on $1/c$ without the need of estimating the spectrum of the $\Sigma_S$.

To conclude, we should point out that these procedures suffer from the problem of replicability, due to the fact that $C_{\alpha}(\mathcal{X}_1)$ is random, implying that, for example,  splitting the data in two different ways could potentially lead to contradictory conclusions. Moreover, it is likely that splitting the sample reduces the power of our tests, since we are not making full use of the information in the data. These issues were tackled in [Guo and Shah, 2023], where rank-transformed subsampling was introduced to increase the power of $k$-fold sample splitting.

§ OPTIMALITY AND EXAMPLES

In this section, we assess the optimality of the oracle test given in Theorem <ref>, under the settings of Examples <ref> and <ref>, i.e. when $\mathbb{S}= \{ \{1,2\} , \{2,3\} ,\ldots, \{d-1,d\}, \{d,1\}\}$, the $d$-cycle, and when $\mathbb{S} = \{[2d], [d]\cup([3d]\setminus[2d]),[3d]\setminus[d]\}$, the block $3$-cycle. These collections $\mathbb{S}$ provide examples where our methodology is provably near rate-optimal. For a given dimension $d$, these two examples further demonstrate the range of optimal rates that can arise for different collections $\mathbb{S}$. Assuming for simplicity that $n_S=n$ for all $S \in \mathbb{S}$, we will see that the optimal rate in the $d$-cycle case is $\{\log(d)/n\}^{1/2}$, while for the block 3-cycle it is $(d/n)^{1/2}$ up to logarithmic factors. Together, these results show that the structure of $\mathbb{S}$ can have a significant effect on the difficulty of the problem.

We will characterise the optimality of a testing procedure using the minimax framework, where we aim at finding the smallest separation between the null and the alternative hypotheses such that there exists a test that can distinguish between $H_0$ and $H_1$ up to a given level of error. More precisely, given $\rho>0$, we are interested in testing

\[
H_0: R(\Sigma_\mathbb{S}) + V(\sigma_{\mathbb{S}}^2) = 0 \quad \text{ vs. } \quad H_1: R(\Sigma_\mathbb{S}) + V(\sigma_{\mathbb{S}}^2) > \rho,
\]

and our goal is to find the smallest value of $\rho$ such that there exists a test with uniform error control. Write $\Psi \equiv \Psi_\mathbb{S}(n_\mathbb{S})$ for the set of all tests, that is measurable functions of the data $(X_{S,i} : S \in \mathbb{S}, i \in [n_S])$ taking values in $\{0,1\}$. Write $\mathcal{P}_\mathbb{S}(0)$ for the set of all sequences of distributions on $(\mathbb{R}^S : S \in \mathbb{S})$ such that the associated correlation matrices and variances satisfy $R(\Sigma_\mathbb{S}) + V(\sigma_{\mathbb{S}}^2) = 0$, and write $\mathcal{P}_\mathbb{S}(\rho)$ for the set of all sequences of distributions on $(\mathbb{R}^S : S \in \mathbb{S})$ such that the associated correlation matrices and variances satisfy $R(\Sigma_\mathbb{S}) + V(\sigma_{\mathbb{S}}^2) > \rho$. Given a  sequence of distributions $P_\mathbb{S} = (P_S : S \in \mathbb{S})$ on $(\mathbb{R}^S : S \in \mathbb{S})$ and a sequence of sample sizes $n_\mathbb{S}=(n_S : S \in \mathbb{S})$, we write $P_\mathbb{S}^{\otimes n_\mathbb{S}}$ for the distribution of the entire dataset $(X_{S,i} : S \in \mathbb{S}, i \in [n_S])$ when each observation is independent and $X_{S,i} \sim P_S$ for each $i \in [n_S]$ and $S \in \mathbb{S}$. For a fixed $\eta \in (0,1)$ we may then define the minimax separation to be

\[
\rho^* \equiv \rho_\mathbb{S}^*(n_\mathbb{S},\eta) := \inf\left\{\rho > 0: \inf_{\varphi \in \Psi} \biggl( \sup_{P_{\mathbb{S},0} \in \mathcal{P}_{\mathbb{S}}(0)} P_{\mathbb{S},0}^{\otimes n_{\mathbb{S}}}(\varphi = 1) + \sup_{P_{\mathbb{S},1} \in \mathcal{P}_{\mathbb{S}}(\rho)} P_{\mathbb{S},1}^{\otimes n_{\mathbb{S}}} (\varphi = 0) \biggr) \leq \eta \right\}
\]

In our analysis we take $\eta=3/4$, but this is an arbitrary choice and any constant value in $(0,1)$ would result in the same qualitative behaviour. In common with previous work on minimax testing, we prove lower bounds on $\rho^*$ by constructing suitable (prior) distributions $\mu_0, \mu_1$ whose support is contained in $\mathcal{P}_{\mathbb{S}}(0), \mathcal{P}_{\mathbb{S}}(\rho)$, respectively. In our proofs it will be sufficient to consider mean-zero Gaussian distributions with suitable priors over their covariance matrices. Having chosen these priors we can bound the minimal error probability by writing

\begin{align*}
\sup_{P_{\mathbb{S},0} \in \mathcal{P}_{\mathbb{S}}(0)} P_{\mathbb{S},0}^{\otimes n_{\mathbb{S}}}(\varphi = 1)  + \sup_{P_{\mathbb{S},1} \in \mathcal{P}_{\mathbb{S}}(\rho)} P_{\mathbb{S},1}^{\otimes n_{\mathbb{S}}} (\varphi = 0) &\geq \mathbb{E}_{\mu_0} P_{\mathbb{S},0}^{\otimes n_{\mathbb{S}}}(\varphi = 1) + \mathbb{E}_{\mu_1} P_{\mathbb{S},1}^{\otimes n_{\mathbb{S}}}(\varphi = 0) \\
& \geq 1 - \operatorname{TV}\bigl( \mathbb{E}_{\mu_0}P_{\mathbb{S},0}^{\otimes n_{\mathbb{S}}} , \mathbb{E}_{\mu_1}P_{\mathbb{S},1}^{\otimes n_{\mathbb{S}}} \bigr),
\end{align*}

$\mathbb{E}_{\mu_i}P_{\mathbb{S},i}^{\otimes n_{\mathbb{S}}}$ denotes the mixture distribution of the dataset resulting from generating $P_{\mathbb{S},i} \sim \mu_i$ then, conditionally on  $P_{\mathbb{S},i}$, generating the data. Then, the idea behind this method of finding a lower bound on $\rho^*$ is to find priors $\mu_0, \mu_1$  with the largest separation possible such that no test can successfully distinguish between $\mathbb{E}_{\mu_0}P_{\mathbb{S},0}^{\otimes n_{\mathbb{S}}}$ and $\mathbb{E}_{\mu_1}P_{\mathbb{S},1}^{\otimes n_{\mathbb{S}}}$.

§.§ Cycles

Graph associated to the $d$-cycle $\mathbb{S}_d=\{\{1,2\},\ldots, \{d,1\}\}$.

Recall that we refer to $\mathbb{S}$ as a $d$-cycle when $\mathbb{S} = \mathbb{S}_d:=\{\{1,2\},\ldots, \{d,1\}\}$. This can be illustrated by Figure <ref>, where an edge in the graph represents an element of $\mathbb{S}$, so that two nodes are connected if and only if the corresponding variables are simultaneously observed. In this subsection, additions in subscripts of the form $(j, j+1)$ for $j \in [d]$ are intended modulo $d$, where $d$ in the size of the cycle, so that $(0,1)$ and $(d,d+1)$ are equivalent to $(d,1)$. Our main statistical result in this subsection is  following minimax lower bound.

Let $\mathbb{S} = \mathbb{S}_d$ for $d \geq 3$, with sample sizes $n_{\mathbb{S}} = (n_1, \ldots, n_d)$. There exists a universal constant $c_1 > 0$ such that

\[
\rho^* \geq c_1 \sqrt{\frac{\log d}{\min_{j \in [d]}n_j}}.
\]

This result shows that the power guarantees for the oracle test given in Theorem <ref> are optimal up to constants, in the case of a $d$-cycle. Combined with this upper bound, our construction of the prior distributions $\mu_0,\mu_1$ in the proof of Theorem <ref> shows that testing the consistency of the variances $\sigma_\mathbb{S}^2$, i.e. testing $H_{0,j} : \sigma_{\{j-1,j\},j}^2 = \sigma_{\{j,j+1\},j}^2$ for each $j=1,\ldots,d$, captures the essential statistical difficulty of the problem.

In the remainder of this subsection we explore the properties of $R(\cdot)$ in the cycle example. In particular, we provide explicit expression in simple cases, we discuss the meaning of maximal incompatibility and we prove a result showing that $R(\cdot)$ is bounded below by the maxima of suitable linear functions if $\Sigma_{\mathbb{S}_d}$ is not too singular. Here we write $\Sigma_{\mathbb{S}_d}:= (\Sigma_{\{1,2\}}, \cdots, \Sigma_{ \{d,1\}})$ for a sequence of $2 \times 2$ correlation matrices with

\[
\Sigma_{j,j+1} = \begin{pmatrix} 1 & \rho_{j,j+1} \\ \rho_{j,j+1} & 1 \end{pmatrix}.
\]

Our next result shows that singular matrices can be removed from $\Sigma_{\mathbb{S}_d}$ when $d\geq4$, without affecting the value of $R(\cdot)$, reducing the length of the cycle.

Fix $d \geq 3$ and $k \geq 1$. Let $\Sigma_{\mathbb{S}_{d+k}}$ be a $(d+k)$-cycle with correlations $(\rho_{\{1,2\}},\ldots,\rho_{\{d+k,1\}})$ such that  $|\rho_{j,j+1}| = 1$ for all $j \in \{d+1, \cdots, d+k \}$ and let $\Sigma_{\mathbb{S}_{d}}$ represent a $d$-cycle with correlations $(\bar{\rho}_{\{1,2\}},\ldots,\bar{\rho}_{\{d,1\}})$ such that $\bar{\rho}_{j,j+1} = \rho_{j, j+1}, \text{ for all } j \in [d-1]$ and

\[
\begin{cases}
\bar{\rho}_{d,1} = \rho_{d, d+1} \text{ if } \prod\limits_{j = d+1}^{d+k} \rho_{j,j+1} = 1\\
\bar{\rho}_{d,1} = -\rho_{d, d+1} \text{ if } \prod\limits_{j = d+1}^{d+k} \rho_{j,j+1} = -1.
\end{cases}
\]

Then we have $R(\Sigma_{\mathbb{S}_{d+k}}) = R(\Sigma_{\mathbb{S}_{d}})$.

This reduction applies when the correlations associated to an edge belonging to the path from node $d+1$ to node $1$ are either $+1$ or $-1$. In this setting, we are allowed to identify node $1$ with node $d+1$ in such a way that the incompatibility measure of the red $d$-cycle $\Sigma_{\mathbb{S}_{d}}$ in Figure <ref> is the same as the one of the original $(d+k)$-cycle $\Sigma_{\mathbb{S}_{d+k}}$. This is to be expected, as $\rho_{j,j+1}=\pm 1$ means that variables $j$ and $j+1$ can be identified, up to change in scale, and the dimensionality of the problem can be reduced. Clearly, this result is invariant under cyclic permutations of the nodes' labels.

Illustration of Proposition <ref>. All the correlations associated to an edge belonging to the path from node $d+1$ to node $1$ are either $+1$ or $-1$, and are such that $\prod_{i = d+1}^{d+k} \rho_{i,i+1} = \pm 1$. The new $d$-cycle $\Sigma_{\mathbb{S}_{d}}$ in red has the same correlations as the $(d+k)$-cycle ${\Sigma}_{\mathbb{S}_{d+k}}$, except for $\bar{\rho}_{d,1}$ corresponding to the new edge $\{ 1, d\}$, which is equal to $\pm \rho_{d, d+1}$.

We now give some explicit expressions for $R(\cdot)$ in special cases and discuss a case for which $R(\Sigma_{\mathbb{S}}) = 1$, meaning that $\Sigma_{\mathbb{S}}$ is maximally incompatible. It will be convenient for the rest of the subsection to reparametrise the correlations as $\rho_j = \cos\theta_j$, with $\theta_j \in [0, \pi]$.

If there are $\theta_1,\theta_2 \in [0,\pi]$ such that $\theta_1 \geq \theta_2$ and

\[
\Sigma_{\mathbb{S}_3} = \Bigg\{\begin{pmatrix}
1 & \cos\theta_1 \\
\cos\theta_1 & 1
\end{pmatrix},\begin{pmatrix}
1 & \cos\theta_2 \\
\cos\theta_2 & 1
\end{pmatrix}, \begin{pmatrix}
1 & 1 \\
1 & 1
\end{pmatrix} \Bigg\},
\]

then $R(\Sigma_{\mathbb{S}_3}) = (\cos\theta_2 - \cos\theta_1)/2$. In particular, setting $\theta_2=0$ we see that if

\[
\Sigma_{\mathbb{S}_3} = \Bigg\{\begin{pmatrix}
1 & \cos\theta_1 \\
\cos\theta_1 & 1
\end{pmatrix}, \begin{pmatrix}
1 & 1 \\
1 & 1
\end{pmatrix}, \begin{pmatrix}
1 & 1 \\
1 & 1
\end{pmatrix}\Bigg\},
\]

then $R(\Sigma_{\mathbb{S}_3}) = \sin^2(\theta_1/2)$. Moreover, assuming without loss of generality that at most one correlation is negative, as justified in Proposition <ref> below, for a general $3$-cycle $\Sigma_{\mathbb{S}_3}$ we have $R(\Sigma_{\mathbb{S}_3}) = 1$ if and only if
$$\Sigma_{{\mathbb{S}}_3} := \left\{\begin{pmatrix}
1 & -1 \\ -1 & 1
\end{pmatrix}, \begin{pmatrix}
1 & 1 \\ 1 & 1
\end{pmatrix}, \begin{pmatrix}
1 & 1 \\ 1 & 1
\end{pmatrix} \right\}.$$

These results can be extended to a general $d$ using Proposition <ref>.

Related to the last claim of Example <ref>, another important property for a $d$-cycle is that we can always assume without loss of generality that at most one $\theta_i$ is larger than $\pi/2$, which is equivalent to having at most one negative $\rho_i$, without changing the value of $R(\cdot)$. This is shown in the following result.

Consider the $d$-cycle with $\mathbb{S}_d=\{\{1,2\},\ldots, \{d,1\}\}$ and $\Sigma_{\mathbb{S}_d}:= (\Sigma_{\{1,2\}}, \cdots, \Sigma_{ \{d,1\}})$, where the correlations $\rho_j = \cos\theta_j$ are uniquely determined by $0 \leq \theta_1, \ldots, \theta_d \leq \pi$. Then, there exists another sequence of angles $0 \leq \tilde{\theta}_1, \ldots, \tilde{\theta_d} \leq \pi$ with at most one $\tilde{\theta}_i$ larger than $\pi/2$ such that the corresponding $d$-cycle $\tilde{\Sigma}_{\mathbb{S}_d}:= (\tilde{\Sigma}_{\{1,2\}}, \cdots, \tilde{\Sigma}_{ \{d,1\}})$ satisfies $R(\tilde{\Sigma}_{\mathbb{S}_d}) = R(\Sigma_{\mathbb{S}_d})$.

As the proof in Section <ref> shows, this new cycle $\tilde{\Sigma}_{\mathbb{S}_d}$ is obtained after changing some $X_i$ into $-X_i$ in such a way that  at most $\cos{\theta_1}$ is negative.

The last result we present on the $d$-cycle gives an explicit lower bound for $R$ in the case that $\Sigma_{\mathbb{S}_d}$ is incompatible. This is related to the results of [Barrett et al., 1993] characterising exactly when the partial correlation matrix

\[\Sigma_{partial} =
\begin{pmatrix} 1 & \cos \theta_1 & \ast & \cdots & \cos \theta_d \\ \cos \theta_1 & 1 & \cos \theta_2 & \cdots & \ast \\ \ast & \cos \theta_2 & 1 & \ldots & \ast \\ \vdots & \vdots & \vdots & \ddots &\vdots \\ \cos \theta_d & \ast & \ast & \cdots & 1 \end{pmatrix}
\]

has a positive semi-definite completion. [Barrett et al., 1993] shows that this is the case if and only if

\begin{align}
\label{eq:barrett_charcterisation}
\sum\limits_{j \in K} \theta_j \leq (|K| -1 ) \pi + \sum\limits_{j \not\in K} \theta_j
\end{align}

for all $K \subseteq [d]$ with $|K|$ odd. Remarkably, this shows that the parametrisation $\rho_{j,j+1}=\cos \theta_j$ allows us to characterise the feasibility of positive semi-definite matrix completion in terms of a finite number of linear inequalities. If we know that $0 \leq \theta_d \leq \theta_{d-1} \leq \ldots \leq \theta_1 \leq \pi$ then this reduces to checking

\[
\sum_{j=1}^k \theta_j \leq (k-1) \pi + \sum_{j=k+1}^d \theta_j
\]

for all odd $k \in [d]$. Furthermore, if $0 \leq \theta_1, \ldots, \theta_d \leq \pi$ with at most one $\theta_j$ larger than $\pi/2$, then $\Sigma_{partial}$ has a positive semi-definite completion if and only if
\[
2\max_{j \in [d]}\theta_j \leq \sum_{j=1}^d \theta_j.
\]

Proposition <ref> shows that we can always work under this setting, so that the problem of whether $\Sigma_{partial}$ has a PSD completion or not is determined by one condition only, namely $2\max_{j \in [d]}\theta_j \leq \sum_{j=1}^d \theta_j$. This is a novel contribution per se, since it is not present in [Barrett et al., 1993]. It is interesting to see that Barrett's characterisation enables us to derive a more explicit expression for $R(\Sigma_{\mathbb{S}_d})$ in the case of a $d$-cycle.

Let $\Sigma$ be the optimum solution to the dual problem (<ref>), and let $\boldsymbol{\varphi}^* = (\varphi_1^*, \ldots, \varphi_d^*)$ be such that $\Sigma_{j,j+1} = \cos\varphi_j^*$ for each $j \in [d]$. Then:

(i) $R(\Sigma_{\mathbb{S}_d}) = |\rho_j - (1-R(\Sigma_{\mathbb{S}_d}))\cos\varphi_j^*|, \text{ for all } j \in [d]$;

(ii)$\boldsymbol{\varphi}^* = (\varphi_1^*, \ldots, \varphi_d^*)$ is unique, and $\boldsymbol{\varphi}^*(\theta_1, \ldots, \theta_d)$ is continuous for varying $(\theta_1, \ldots, \theta_d) \in [0, \pi]^d$;

(iii) if $\theta_1 = \max_{j \in [d]} \theta_j$, with $\theta_2,\ldots,\theta_d \leq \pi/2$, then
\[
1-R(\Sigma_{\mathbb{S}_d}) = \frac{1-\epsilon_j \cos\theta_j}{1-\epsilon_j \cos\varphi_j^*}, \quad \text{ for all } j \in [d],
\]

where $\boldsymbol{\epsilon}_d = (\epsilon_1, \ldots, \epsilon_d) = (-1, +\boldsymbol{1}_{d-1})$. Also, $\varphi_1^* = \sum_{j = 2}^d \varphi_j^*$.

Observe that part (ii) only says that the entries of $\Sigma$ corresponding to the cycle pattern are unique, not the whole $\Sigma$ itself. Indeed, given the unique optimal $\boldsymbol{\varphi}^*$, there may exist infinitely many positive semi-definite completions. In fact, if $A$ is a partial symmetric matrix admitting a positive semi-definite completion, then there exists a unique positive semi-definite completion with maximum determinant (see Theorem 2 in [Grone et al., 1984]). Also, as a sanity check, observe that the optimal choice of the signs in part (iii) makes $\Sigma_{\mathbb{S}}'$ in (<ref>) as incompatible as possible, in accordance with the dual representation given in Proposition <ref>. For a general sequence of angles $(\theta_1, \ldots, \theta_d)$, it is sufficient to use the transformation given in  Proposition <ref> to reduce to the case where at most one angle is larger than $\pi/2$, choose $\boldsymbol{\epsilon}_d$ as outlined above, and perform the inverse transformation to obtain the signs for the original $(\theta_1, \ldots, \theta_d)$. As an immediate corollary of this, it is easy to see that, under the same set of hypotheses, we have
\[
1-R(\Sigma_{\mathbb{S}_d}) = \frac{1+ \cos\theta_1}{1+ \cos\varphi_1^*},
\]

where $\varphi_1^*$ is the solution of
\[
\begin{cases}
\varphi_1^* = \sum_{j = 2}^d \varphi_j^* \\
\\
\cos \varphi_j^* = 1 - \frac{1-\cos\theta_j}{1+\cos\theta_1}(1+\cos\varphi_1^*), \quad \text{ for all }  j \in \{2, \ldots, d\}.
\end{cases}
\]

This is a relatively explicit expression for $R(\cdot)$ for a general $d$-cycle.

Our final result in this subsection shows that, provided not too many of our input matrices are close to being singular, $R(\Sigma_{\mathbb{S}_d})$ can be bounded below by a finite maximum of linear functionals that is zero if and only if $\Sigma_{\mathbb{S}_d}$ is compatible. This lower bound constitutes another sanity check for our measure $R(\cdot)$, since the quantities appearing in the lower bound are a natural quantitive version of the qualitative conditions given in [Barrett et al., 1993] to check whether the partial matrix $\Sigma_{partial}$ defined above admits a PSD completion.

Consider the $d$-cycle with $\mathbb{S}=\mathbb{S}_d=\{\{1,2\},\ldots, \{d,1\}\}$ and suppose that

\[
\Sigma_{\{j,j+1\}} = \begin{pmatrix} 1 & \cos\theta_j \\ \cos\theta_j & 1 \end{pmatrix}.
\]

Assume further that there exist $c>0$ and two indices $k,j \in [d]$ such that $1-\rho_j^2 \geq c, 1-\rho_k^2 \geq c$, so that $\Sigma_{\{j,j+1\}}$ and $\Sigma_{\{k,k+1\}}$ are bounded away from singularity. Then, whenever $\Sigma_{\mathbb{S}_d}$ is incompatible, we have

\[
R(\Sigma_{\mathbb{S}_d}) \geq c' \mathop{\max\limits_{K \subseteq [d]}}\limits_{|K| \text{ odd}}\left(\sum\limits_{i \in K}\theta_i - (|K| -1) \pi - \sum\limits_{i \in K^c}\theta_i\right),
\]

where $c'>0$ depends only on $c$.

The proof of Proposition <ref> can be found in Section <ref>. First, observe that this lower bound reduces to

\[
c' \left(\theta_1 - \sum_{i = 2}^d \theta_i \right)_+,
\]

in the case that $\theta_1 = \max_{j \in [d]} \theta_j$ and $\theta_2,\ldots,\theta_d \leq \pi/2$, which we have already argued that we may assume without loss of generality. Furthermore, as a sanity check, the simple explicit expressions found in Example <ref>, in which we have seen that

\[
1 & \cos\theta_1 \\
\cos\theta_1 & 1
\end{pmatrix},\begin{pmatrix}
1 & \cos\theta_2 \\
\cos\theta_2 & 1
\end{pmatrix}, \begin{pmatrix}
1 & 1 \\
1 & 1
\end{pmatrix}, \ldots, \begin{pmatrix}
1 & 1 \\
1 & 1
\end{pmatrix}\Bigg\}\right) = (\cos\theta_2 - \cos\theta_1)/2,
\]

is in accordance with Proposition <ref>, since $\cos\theta_2 - \cos\theta_1 \gtrsim_c \theta_1 - \theta_2$ when $\theta_1, \theta_2$ are bounded away from $\{0, \pi \}$.

§.§ Block cycles

So far, we studied with particular care the case of a $d$-cycle, which is a relatively simple high-dimensional setting, since it is a collection of $d$ two-dimensional distributions. We now describe an evolution of this setting, where we consider a block-matrix version of the $3$-cycle. In this case the number of variables per missingness pattern is large and we will see that the minimax separation rates are correspondingly much larger than in the $d$-cycle, though the number of variables is of the same order.

Let $\mathbb{S} = \{[2d], [d]\cup([3d]\setminus[2d]),[3d]\setminus[d]\}$ for some $d \geq 1$. Writing $n_{\mathbb{S}} = (n_1, n_2, n_3)$ for the sample sizes within each pattern, there exists a universal constant $c_1 > 0$ such that

\[
\rho^* \geq c_1 \sqrt{\frac{d}{(n_1 \wedge n_2) \log^4 (ed)}}
\]

whenever $n_1 \wedge n_2 \geq d/2$.

This result shows that, up to logarithmic factors in $d$, the minimax separation rates for this testing problem are the same as the minimax estimation rates for estimating $\Sigma_\mathbb{S}$ in the operator norm distance. This is related to the fact that $R(\Sigma_\mathbb{S})$ is a non-smooth functional of $\Sigma_\mathbb{S}$. Indeed, the following result shows that we can construct examples of $\Sigma_\mathbb{S}$ such that $R(\Sigma_\mathbb{S})$ can be bounded below using the function $x \mapsto \max(0,x)$; see below for more discussion of the relevant literature.

Consider the set of patterns $\mathbb{S} = \{[2d], [d]\cup([3d]\setminus[2d]),[3d]\setminus[d]\}$ for some $d \geq 1$, and suppose that
\[
\Sigma_{\mathbb{S}} = \Bigg\{\begin{pmatrix}
I_d & P \\
P^T & I_d
\end{pmatrix},\begin{pmatrix}
I_d & -P \\
-P^T & I_d
\end{pmatrix}, \begin{pmatrix}
I_d & \beta I_d \\
\beta I_d & I_d
\end{pmatrix}\Bigg\},
\]

for some $P \in \mathbb{R}^{d \times d}$ such that $\|P\|_2 \leq 1$ and some $\beta \in [0,1]$. Then:

(i) $R(\Sigma_{\mathbb{S}}) = 0$ if and only if $\|P\|_2^2 \leq \frac{1-\beta}{2}$,

(ii) $R(\Sigma_{\mathbb{S}}) \geq \frac{3}{4d}\sum_{j = 1}^d (\sigma_j^2(P) -  \frac{1-\beta}{2})_+$, where $\sigma_j(P)$ is the $j$-th singular value of $P$.

This shows that, for $\Sigma_\mathbb{S}$ of the form above, we can relate our testing problem to the problem of testing whether the vector of squared singular values of $P$ belongs to the orthant $(-\infty,0]^d$, or is separated from it in the $\ell_1$ distance. In a Gaussian location model a similar problem, measuring separation with the $\ell_2$ distance, was considered by Blanchard et al., 2018, and part of our lower bound construction is inspired by this work. However, the consideration of singular values of matrices rather than Gaussian means means  that new technical tools are required. In this regard, the techniques of [Thépaut and Verzelen, 2021], who consider the estimation of quantities of the form $\sum_{j=1}^d \sigma_j(P)^q$ for $q >0$, are useful. We also mention that such problem are related to the estimation of $\ell_1$ distances, for which good references include Cai and Low, 2011 and Jiao et al., 2016.

§ NUMERICAL STUDIES

The tests introduced in Theorem <ref> and Proposition <ref> give finite-sample Type-I error control over $\mathcal{P}_\mathbb{S}^0$, the parameter space associated to the null hypothesis. However, these procedures can be too conservative in some examples of interest, and have the further downside of depending on the unknown subgaussian variance proxy $\nu$. For these reasons, we propose Monte-Carlo versions of our tests, that can be applied without any knowledge of unknown parameters, and we compare it with Little's test [Little, 1988]. In Section <ref> we introduce a bootstrap method that checks the incompability of sample correlation matrices, while in Section <ref> we extend this method to also check the consistency of samples means and variances.

Little's test can be applied when all pairs of variables are observed together, so that the EM algorithm [Dempster et al., 1977] can be applied to find estimators $\hat{\mu}$ and $\hat{\Lambda}$ of the mean and covariance matrix of the data under the null hypothesis of MCAR. Little's test is a generalised likelihood ratio test whose validity is is based on the assumption that the data $(X_{S,i} : S \in \mathbb{S}, i \in [n_S])$ are Gaussian. Writing $\tilde{\Lambda} = n_{\mathrm{tot}}\hat{\Lambda}/(n_{\mathrm{tot}}-1)$, where $n_{\mathrm{tot}} = \sum_{S \in \mathbb{S}} n_S$, define
\[
d^2 = \sum_{S \in \mathbb{S}} n_S (\bar{X_S} - \hat{\mu}_{|S})\tilde{\Lambda}_{|S}^{-1}(\bar{X_S} - \hat{\mu}_{|S}),
\]

\[
d^2_{\mathrm{cov}} =  \sum_{S \in \mathbb{S}} n_S[\operatorname{tr}(\hat{\Sigma}_S \hat{\Lambda}_{|S}^{-1})-|S|-\log|\hat{\Sigma}_S| + \log|\hat{\Lambda}_{|S}|],
\]

\[
d^2_{\mathrm{aug}} = d^2 + d^2_{\mathrm{cov}}.
\]

Then, under MCAR, $d^2_{\mathrm{aug}}$ converges in law to a $\chi^2$-distribution with
\[
f = \sum_{S \in \mathbb{S}} \frac{1}{2}|S|(|S|+3) - \frac{1}{2}d(d+3),
\]

degrees of freedom by Wilks' theorem. Based on these asymptotic results, Little's test rejects MCAR if and only if $d^2_{\mathrm{aug}} > \chi^2_f(1-\alpha)$,
where $\chi^2_f(1-\alpha)$ is such that $\mathbb{P}(W \geq \chi^2_f(1-\alpha)) = \alpha$, and where $W$ is $\chi^2$-distributed with $f$ degrees of freedom. Similarly, one can define a test based on $d^2_{\mathrm{cov}}$ that ignores the means and only considers the partial covariance matrices, which converges to a $\chi^2$-distribution with
\[
f' = \sum_{S \in \mathbb{S}} \frac{1}{2}|S|(|S|+1) - \frac{1}{2}d(d+1)
\]

degrees of freedom by Wilks' theorem. In the next section, we will compare our bootstrap method with these two versions of Little's test, one based on $d^2_{\mathrm{aug}}$, the other on $d^2_{\mathrm{cov}}$. This is because in Section <ref> departures from the null are due to an incompatible sequence of correlation matrices, while the means are assumed to be constant, hence we do not want to give ourselves a clear advantage by comparing only with $d^2_{\mathrm{aug}}$.

§.§ Correlation matrices and simulations for $d$-cycles

In this section, we design a bootstrap version of our test, and compare it with Little's test in detecting departures from MCAR due to an incompatible sequence of correlation matrices $\Sigma_{\mathbb{S}}$. Recall that we write $\hat{\Sigma}_\mathbb{S} = (\hat{\Sigma}_S)_{S \in \mathbb{S}}$, where $\hat{\Sigma}_S$ is the sample correlation matrix of the data $X_S:=(X_{S,i},i\in [n_S])$ for $S \in \mathbb{S}$. From Proposition <ref> we can write $\hat{\Sigma}_\mathbb{S} = (1-R(\hat{\Sigma}_\mathbb{S}))\hat{Q}_\mathbb{S} + R(\hat{\Sigma}_\mathbb{S})\hat{\Sigma}'_\mathbb{S}$, where $\hat{Q}_\mathbb{S}$ can be thought as the closest compatible sequence of correlation matrices to $\hat{\Sigma}_\mathbb{S}$, and can be computed at the same time as the test statistic $R^{(0)} := R(\hat{\Sigma}_\mathbb{S})$. We then transform the original data by calculating $\tilde{X}_{S} := \hat{Q}_S^{1/2}\hat{\Sigma}_S^{-1/2}\operatorname{scale}(X_{S})$ for all $S \in \mathbb{S}$, where $\mathrm{scale}(\cdot)$ represents the step of standardising the original data. This transformation means that the sample correlation matrices of $\tilde{X}_\mathbb{S}:=(\tilde{X}_S : S \in \mathbb{S})$ are compatible. Fixing $B \in \mathbb{N}$, for each $b \in [B]$ and $S \in \mathbb{S}$ we generate $\tilde{X}_{S}^{(b)}$ as a nonparametric bootstrap sample from $\tilde{X}_{S}$ and calculate the sample correlation matrix $\hat{\Sigma}_{S,b} = \operatorname{Cor}(\tilde{X}_{S}^{(b)})$. Then, for each $b \in [B]$ we compute the corresponding test statistic $R^{(b)} := R(\hat{\Sigma}_{S,b} : S \in \mathbb{S})$. Finally, we reject $H_0$ at a significance level $\alpha \in (0,1)$ if and only if $1 + \sum_{i=1}^B \mathbbm{1}\{R^{(b)} \geq R^{(0)}\} \geq \alpha(1+B)$.

We compare this test with Little's procedure in the settings given in Theorem <ref>, namely in the case of a $d$-cycle. For our first settings, we set $n_\mathbb{S} = (n_S)_{S \in \mathbb{S}} = (200, \ldots, 200)$, and simulate $X_{\{j,j+1\},i} \overset{\text{i.i.d.}}{\sim} N(\boldsymbol{0}_2, \Sigma_{\{j,j+1 \}})$ for $i \in [200]$ and $j \in [d]$, where
\[
\Sigma_{\mathbb{S}_d} = \left\{\begin{pmatrix}
1 & \cos\theta_1 \\
\cos\theta_1 & 1
\end{pmatrix}, \ldots,  \begin{pmatrix}
1 & \cos\theta_d \\
\cos\theta_d & 1
\end{pmatrix}\right\},
\]

for certain values of $\theta_1, \ldots, \theta_d \in [0, \pi]$, and compare our bootstrap test with Little's procedure. Here, we repeat the the experiment $M = 200$ times, and report the average decision as an estimate of the power function. This makes sense only for $d=3$, while for $d \geq 4$ there exists at least one pair of variables that are never observed together, making the EM algorithm to estimate $\hat{\Lambda}$ inapplicable. As for the case $d = 3$, Figure <ref> shows two different simulations, with different values of $(\theta_1, \theta_2, \theta_3)$, where our test performs very similarly to Little's tests.

Simulation of the power functions of our method with $B = 99$ (blue), Little's method based on $d^2_{\mathrm{aug}}$ (green), and Little's method based on $d^2_{\mathrm{cov}}$ (orange), with Gaussian data. In each example, we vary $\theta_1 \in [\theta_2 + \theta_3, (\theta_2 + \theta_3 + \pi)/2]$, with $(\theta_2, \theta_3)$ equal to $(\pi/3, \pi/6)$ (left), $(\pi/4, \pi/4)$ (right). For each of this setting, we repeat the experiment $M = 200$ times, and report the average decision. The nominal level $\alpha = 0.05$ in red.

As stated above, for $d \geq 4$ Little's test cannot be applied, while our test remains valid since there are no constraints on $\mathbb{S}$. In Figure <ref> below, we show the power function of our bootstrap test in the case of a $d$-cycle, with $d \in \{100, 200\}$, with $\theta_2 = \ldots = \theta_d = \frac{\pi}{2(d-1)}$, and varying $\theta_1$ in $[\pi/2, 5\pi/8]$. We repeat the procedure $M = 100$ times, and report the average decision as an estimate of the power function.

Simulation of the power functions of our method with $B=99$(blue) for $d = 100$ (left), and $d = 200$ (right), with Gaussian data. In each example, we fix $\theta_2 = \ldots = \theta_d = \frac{\pi}{2(d-1)}$, and vary $\theta_1$ in $[\pi/2, 5\pi/8]$. For each of this setting, we repeat the experiment $M = 100$ times, and report the average decision. The nominal level $\alpha = 0.05$ in red.

Our simulations so far have used Gaussian data, so that Little's test is valid. We now repeat our simulations with a heavy-tailed data distribution in order to assess the robustness of the methods. To this aim, we consider again a $3$-cycle, and generate $X_{\{j,j+1\},i} \overset{\text{i.i.d.}}{\sim} \log N(\boldsymbol{0}_2, \Sigma_{\{j,j+1 \}})$ for all $i \in [200], j \in [3]$, where $\log N(\boldsymbol{0}_2, \Sigma_{\{j,j+1 \}})$ stands for the log-normal distribution, meaning that if $Y \sim \log N(\boldsymbol{0}_2, \Sigma_{\{j,j+1 \}})$ then $Y_i = e^{Z_i}$, with $Z \sim N(\boldsymbol{0}_2, \Sigma_{\{j,j+1 \}})$. Figure <ref> below shows the analogue of Figure <ref>, in the sense that the parameters $(\theta_1, \theta_2, \theta_3)$ are the same, but we generated artificial data from a multivariate log-normal distribution rather than a Gaussian distribution. It is interesting to see that Little's test does not have Type-I error control.

On the other hand, our test succeeds in controlling the Type-I error and, although being conservative, its power increases as $\Sigma_{\mathbb{S}_3}$ becomes more incompatible.

Simulation of the power functions of our method with $B =99$(blue), Little's method based on $d^2_{\mathrm{aug}}$ (green), and Little's method based on $d^2_{\mathrm{cov}}$ (orange), with log-normal data. In each example, we vary $\theta_1 \in [\theta_2 + \theta_3, (\theta_2 + \theta_3 + \pi)/2]$, with $(\theta_2, \theta_3)$ equal to $(\pi/3, \pi/6)$ (left), $(\pi/4, \pi/4)$ (right). For each of this setting, we repeat the experiment $M = 200$ times, and report the average decision. The nominal level $\alpha = 0.05$ in red.

§.§ Omnibus approach

We now aim at designing a test that
is able to detect departures from MCAR due to inconsistent means and variances as well as incompatible correlation matrices. To this aim, we define a new bootstrap test that checks both compatibility of $\Sigma_{\mathbb{S}}$, and consistency of $\sigma^2_\mathbb{S}$ and $\mu_\mathbb{S}$, the sequence of means along the different patterns in $\mathbb{S}$. We will use as test statistics estimators of $T = R(\Sigma_\mathbb{S}) + V(\sigma_{\mathbb{S}}^2) + M(\mu_{\mathbb{S}})$, where
\[
M(\mu_{\mathbb{S}}) = \max_{j \in [d]} \max_{S_1, S_2 \in \mathbb{S}_j} |\mu_{S_1,j} - \mu_{S_2,j}|.
\]

Observe that $M(\mu_{\mathbb{S}}) = 0$ if and only if $\mu_{\mathbb{S}}$ is consistent. This results in Algorithm <ref>, which is implemented in the R-package  [Berrett et al., 2022].
MCAR bootstrap test checking compatibility of correlation matrices, and consistency of means and variances
[1]
Given data $X_\mathbb{S}$, compute $\hat{\sigma}^2_\mathbb{S} = \operatorname{Var}X_\mathbb{S}$, and rescale it such that $\operatorname{av}_j(\hat{\sigma}^2_{\mathbb{S}}) = 1$ for all $j \in [d]$; i.e. replace $X_j$ with $X_j/\sqrt{\operatorname{av}_j(\hat{\sigma}^2_{\mathbb{S}})}$ for all $j \in [d]$.
Compute $\hat{\mu}_\mathbb{S} = \mathbb{E}X_\mathbb{S}$, $\hat{M}_\mathbb{S} = \operatorname{Cov}X_\mathbb{S} = \operatorname{diag}( \sigma_{\mathbb{S}}^2)^{1/2} \cdot \Sigma_{\mathbb{S}} \cdot \operatorname{diag}(\sigma_{\mathbb{S}}^2)^{1/2}$.
Compute $T^{(0)} = R(\hat{\Sigma}_\mathbb{S}) + V(\hat{\sigma}_{\mathbb{S}}^2) + M(\hat{\mu}_{\mathbb{S}})$, and compute at the same time the dual decomposition $\hat{\Sigma}_\mathbb{S} = (1-R(\hat{\Sigma}_\mathbb{S}))\hat{Q}_\mathbb{S} + R(\hat{\Sigma}_\mathbb{S})\hat{\Sigma}'_\mathbb{S}$.
Rotate the original data $X_\mathbb{S}$, i.e. for all $S \in \mathbb{S}$, for all $i \in [n_S]$ do $\tilde{X}_{S,i} = \hat{Q}_S^{1/2}\hat{M}_S^{-1/2}(X_{S,i}-\hat{\mu}_S + \hat{\mu}_{|S})$, where $\hat{\mu}_j = |\mathbb{S}_j|^{-1}\sum_{S\in \mathbb{S}_j} \mu_{S,j}$.
$b \in [B]$
For all $S \in \mathbb{S}$, let $\tilde{X}_{S,i}^{(b)}$ be a nonparametric bootstrap sample from $\tilde{X}_{S,i}$, for $i \in [n_S]$.
Compute $\hat{\mu}_{\mathbb{S},b} = \mathbb{E}X_\mathbb{S}^{(b)}$, $\hat{M}_{\mathbb{S},b} = \operatorname{Cov}X_\mathbb{S}^{(b)} = \operatorname{diag}( \hat{\sigma}_{\mathbb{S},b}^2)^{1/2} \cdot \hat{\Sigma}_{\mathbb{S},b} \cdot \operatorname{diag}( \hat{\sigma}_{\mathbb{S},b}^2)^{1/2}$.
Compute $T^{(b)} = R(\hat{\Sigma}_{\mathbb{S},b}) + V(\hat{\sigma}_{\mathbb{S},b}^{2}) + M(\hat{\mu}_{\mathbb{S},b})$.
Reject $H_0$ if and only if $1 + \sum_{i=1}^B \mathbbm{1}\{T^{(b)} \geq T^{(0)}\} \geq \alpha(1+B)$.

As before, we compare this algorithm with Little's test based on $d^2_{\mathrm{aug}}$. In this section we generate complete artificial data according to various distributions, and then delete entries using the R package  [Rockel, 2020]. MCAR data are generated with the function , where each entry of the data matrix is deleted independently of the others with probability $p \in (0,1)$. Deviations from the null are generated by partitioning the columns in two groups, group A where the missing values are generated, and group B which determines the missingness mechanism, with two different mechanisms being considered. First,  sets threshold values, splits the rows into two further groups depending on whether columns in group B have values greater or smaller than the threshold, and deletes some entries in columns in group A in a such a way that the probability for a value to be missing in group A divided by the
probability for a value to be missing in group B equals 1 divided by , with  to be specified as an input parameter. Second,  deletes each entry in a column of group A with probability proportional to the rank of the same row in the corresponding column of group B. For further details on these functions, and other methods to generate MCAR, MAR, MNAR data, refer to [Santos et al., 2019]. These three functions were also chosen in the numerical analysis of a test of MCAR based on U-statistics in Aleksić, 2023.

Type-I error under MCAR data generate with , for varying probability $p$ of having a missing value. Data from Clayton copula with parameter $1$ and log-normal margins. $B = 99, M =300$.

Power function under MAR data generate with , for varying probability $p$ of having a missing value. Data from the same Clayton coupla. $B = 99, M =300$.

Power function under MAR data generate with , for varying probability $p$ of having a missing value. Data from the same Clayton coupla. $B = 99, M =300$.

For Figures <ref>, <ref>, <ref>, we generated $5$-dimensional datasets of sample size $n = 1000$ distributed according to a Clayton copula, with parameter $1$ and log-normal margins, using the function  from the R-package  [Hofert et al., 2020]. For Figure <ref> we deleted the first two variables with  for different values of $p \in \{0.03, 0.06, \ldots, 0.3\}$, in order to get an artificial setting coming from the null. For each $p$, we repeat the simulation $300$ times, and report the average Type-I error. Alternatives to the null were generated using , with $x = 9$, for Figure <ref>, and  for Figure <ref>. Again, for each $p$, we repeat the simulations $300$ times, and report the average power. The simulations show that Little's test is not able to recognise MCAR in this setting, and rejects the null hypothesis with high probability. On the other hand, Algorithm <ref> has good control of the Type-I error, although being a little conservative, and its power increases as the missingness probability $p$ gets bigger and the effective sample sizes for the incomplete patterns increase. Our test performs slightly better in the case of a $3$-dimensional dataset of sample size $n = 1000$ distributed according to a Clayton copula, with parameter $1$ and chi-squared margins. The results of these simulations are shown in Figures <ref>, <ref>, <ref>. In this case, we delete the first two columns, while the third one is always complete. Here our method retains Type I error control and is more powerful than in the previous settings, while Little's test does not have good control of the Type I error.

Type-I error under MCAR data generate with , for varying probability $p$ of having a missing value. Data from Clayton copula with parameter $1$ and chi-squared margins. $B = 99, M =300$.

Power function under MAR data generate with , for varying probability $p$ of having a missing value. Data from the same Clayton coupla. $B = 99, M =300$.

Power function under MAR data generate with , for varying probability $p$ of having a missing value. Data from the same Clayton coupla. $B = 99, M =300$.

§ PROOFS

§.§ Proofs for Section <ref>

For any $X \in \mathcal{M}$ and $X_\mathbb{S} \in \mathcal{M}_\mathbb{S}$ we have

\begin{align*}
\langle AX, X_\mathbb{S} \rangle_\mathbb{S} &= \sum_{S \in \mathbb{S}} \sum_{j,j' \in S} ((AX)_S)_{jj'} (X_S)_{jj'} = \sum_{S \in \mathbb{S}} \sum_{j,j' \in S} X_{jj'} (X_S)_{jj'} \\
&= \sum_{j,j'=1}^d X_{jj'} \sum_{S \in \mathbb{S}} \mathbbm{1}_{j,j' \in S} (X_S)_{jj'} = \langle X , A^* X_\mathbb{S} \rangle,
\end{align*}

as claimed.

The strategy is to use a semi-definite programming version of Farkas' lemma. This is well known in the relevant literature, but we provide a statement and short proof for completeness; see Proposition <ref> in Appendix <ref>. First, rewrite the matrix completion problem

\[
\text{find } \Sigma \in \mathcal{M} \text{ such that } \begin{cases}
\Sigma_{jj'} = (\Sigma_S)_{jj'}, \forall S \in \mathbb{S}_{jj'} \\
\Sigma \succcurlyeq 0
\end{cases}
\]

\begin{align}\label{eq:SDP_completion}
\text{find } \Sigma \in \mathcal{M} \text{ such that } \begin{cases}
\langle \Sigma, E_{jj'} \rangle = (\Sigma_S)_{jj'}, \forall S \in \mathbb{S}_{jj'} \\
\Sigma \succcurlyeq 0
\end{cases}
\end{align}

where $E_{jj'} = (\boldsymbol{e}_{j} \boldsymbol{e}_{j'}^T + \boldsymbol{e}_{j'} \boldsymbol{e}_{j}^T)/2$ and $\boldsymbol{e}_j$ is the $j$-th column vector of the standard orthonormal basis of $\mathbb{R}^d$. In order to apply the semi-definite version of Farkas' lemma we transform our problem so that the equality constraints have zero on the right-hand side. To this end, define

\[
H_j := \begin{pmatrix}
0 & \boldsymbol{e}_j^T/2 \\
\boldsymbol{e}_j/2 & \boldsymbol{O}
\end{pmatrix}
\text{ and }
G_{S, jj'} := \begin{pmatrix}
-(\Sigma_S)_{jj'} & \boldsymbol{0}^T \\
\boldsymbol{0} & E_{jj'}
\end{pmatrix},
\]

and consider the completion problem
\begin{align}\label{eq:SDP_completion2}
\text{find } \tilde{\Sigma} \in \mathcal{M} \text{ such that } \begin{cases}
\langle \tilde{\Sigma}, H_j \rangle = 0, \forall j \in [d] \\
\langle \tilde{\Sigma}, G_{S, jj'} \rangle = 0, \forall S \in \mathbb{S}_{jj'} \\
\tilde{\Sigma} \succcurlyeq 0.
\end{cases}
\end{align}

The condition $\langle \tilde{\Sigma}, H_j \rangle = 0, \forall j \in [d]$ forces $\tilde{\Sigma}$ to be in block diagonal form

\[
\tilde{\Sigma} := \begin{pmatrix}
\gamma_{0,0} & \boldsymbol{0}^T \\
\boldsymbol{0} & \Sigma
\end{pmatrix}.
\]

Now, observe that (<ref>) has a solution if and only if (<ref>) has a non-zero solution. Indeed, for every solution $\Sigma_0$ of (<ref>), then $\operatorname{diag}(1, \Sigma_0)$ is a solution of (<ref>). On the other hand, suppose that $\tilde{\Sigma}_0 = \operatorname{diag}(\gamma_{0,0}, \Sigma_0) \neq \boldsymbol{O}$ is a solution of (<ref>). This implies that $\gamma_{0,0} \neq 0$, otherwise $0 = \langle \tilde{\Sigma}_0, G_{S, jj'} \rangle = -\gamma_{0,0}(\Sigma_S)_{jj'} + \Sigma_{jj'} = \Sigma_{jj'}$, which would imply $\tilde{\Sigma}_0 = \boldsymbol{O}$. Being $\gamma_{0,0} \neq 0$, we can rescale the bigger block in $\tilde{\Sigma}_0$ by $\gamma_{0,0}$, i.e. $\Sigma_0 =: \gamma_{0,0} G$, and get $0 = \langle \tilde{\Sigma}_0, G_{S, jj'} \rangle = -\gamma_{0,0}(\Sigma_S)_{jj'} + \gamma_{0,0} G_{jj'} = -(\Sigma_S)_{jj'} + G_{jj'}$, which shows that $G$ is a solution of (<ref>). This further implies that we can assume without loss of generality that $\gamma_{0,0} = 1$ when (<ref>) admits a non-zero solution. Now, by Proposition <ref>, we know that (<ref>) has a non-zero solution $\tilde{\Sigma} = \operatorname{diag}(1, \Sigma)$ if and only if

\[
\sum\limits_{S \in \mathbb{S}}\sum\limits_{j,j' \in S} (X_S)_{jj'}G_{S, jj'} = \sum\limits_{S \in \mathbb{S}}\begin{pmatrix}
- \langle \Sigma_S, X_S \rangle & \boldsymbol{0}^T \\
\boldsymbol{0} & \frac{1}{2}X_S
\end{pmatrix} = \begin{pmatrix}
- \langle \Sigma_{\mathbb{S}}, X_{\mathbb{S}} \rangle & \boldsymbol{0}^T \\
\boldsymbol{0} & \frac{1}{2} A^* X_{\mathbb{S}}
\end{pmatrix} \nsucc 0,
\]

for all sequences of matrices $X_{\mathbb{S}}$, not necessarily PSD. Now, this block matrix is positive definite if and only if both $A^* X_{\mathbb{S}} \succ 0$ and $\langle \Sigma_{\mathbb{S}}, X_{\mathbb{S}} \rangle < 0$. Hence, (<ref>) has a non-zero solution if and only if $\langle \Sigma_{\mathbb{S}}, X_{\mathbb{S}} \rangle \geq 0$ for all $X_\mathbb{S}$ such that $A^* X_{\mathbb{S}} \succ 0$, and the claim follows.

Weak duality, i.e. LHS $\leq$ RHS, always holds for SDPs (see [Blekherman et al., 2012]), but we include a short proof for the sake of completeness. In fact, for any $\Sigma_\mathbb{S} \in \mathcal{P}_\mathbb{S}$, we can rewrite
\begin{align}\label{eq:RHS_dual}
\inf \{ \epsilon \in [0,1] : \Sigma_\mathbb{S} \in (1-\epsilon) \mathcal{P}_\mathbb{S}^0 + \epsilon \mathcal{P}_\mathbb{S} \}
\end{align}

\begin{align*}
\inf\{  \epsilon \in [0,1] : \Sigma_\mathbb{S} \in (1-\epsilon)  \mathcal{P}_\mathbb{S}^0 + \epsilon \mathcal{P}_\mathbb{S}\} &= 1 - \sup\{ \epsilon \in [0,1] : \Sigma_\mathbb{S} \in \epsilon \mathcal{P}_\mathbb{S}^0 + (1-\epsilon) \mathcal{P}_\mathbb{S} \} \\
& = 1 - \frac{1}{d} \sup\{ \mathrm{tr}(\Sigma) : \Sigma \in \mathcal{P}^*, \Sigma_\mathbb{S} - A \Sigma \succeq_\mathbb{S} 0, \Sigma_{11} = \ldots = \Sigma_{dd} \}.
\end{align*}

Now, for any $Y_\mathbb{S} \in \mathcal{P}_\mathbb{S}^*$ such that $A^* Y_\mathbb{S} + Y \succeq I_d \text{ for some } Y \in \mathcal{Y}$, and any $\Sigma \in \mathcal{P}^*$ such that $\Sigma_\mathbb{S} - A \Sigma \succeq_\mathbb{S} 0$, we have

\begin{align*}
\mathrm{tr}(\Sigma) &= \langle I_d, \Sigma \rangle = -\langle A^*Y_\mathbb{S} + Y - I_d, \Sigma \rangle + \langle A^*Y_\mathbb{S} + Y, \Sigma \rangle \leq \langle A^*Y_\mathbb{S}, \Sigma \rangle + \langle Y, \Sigma \rangle\\
&= \langle A^*Y_\mathbb{S}, \Sigma \rangle = \langle Y_\mathbb{S}, A\Sigma \rangle_\mathbb{S} = \langle Y_\mathbb{S}, \Sigma_\mathbb{S} \rangle_\mathbb{S} - \langle Y_\mathbb{S}, \Sigma_\mathbb{S} - A\Sigma \rangle_\mathbb{S} \leq \langle Y_\mathbb{S}, \Sigma_\mathbb{S} \rangle_\mathbb{S}.
\end{align*}

This shows that (<ref>)
is lower bounded by

\begin{align}\label{eq:LHS_dual}
1 - \frac{1}{d} \inf\{ \langle Y_\mathbb{S}, \Sigma_\mathbb{S} \rangle_\mathbb{S} : Y_\mathbb{S} \in \mathcal{P}_\mathbb{S}^*, A^*Y_\mathbb{S} + Y \succeq I_d \}.
\end{align}

Weak duality follows upon noting that $A^* X_\mathbb{S}^0 = I_d$ and $\langle X_\mathbb{S}^0, \Sigma_\mathbb{S} \rangle_\mathbb{S} =d$ and setting $X_\mathbb{S}=Y_\mathbb{S}-X_\mathbb{S}^0$. This is not surprising, as we already mentioned that weak duality always holds for SDP problems.

We will now prove strong duality for this problem. Our strategy is to write our primal and dual problems in standard form and check Slater's condition for the primal problem (<ref>). We already mentioned that (<ref>) can be written as

\[
1 - \frac{1}{d} \sup\{ \mathrm{tr}(\Sigma) : \Sigma \in \mathcal{P}^*, \Sigma_{11} = \ldots = \Sigma_{dd}, \Sigma_\mathbb{S} - A\Sigma \succeq_\mathbb{S} 0\}.
\]

We now write this maximisation problem in standard form by introducing variables $(Z_S : S \in \mathbb{S}) = \Sigma_\mathbb{S} - A \Sigma \in \mathcal{P}_\mathbb{S}^*$.  Enumerating $\mathbb{S}$ as $\{S_1,\ldots,S_m\}$, we instead optimise over block-diagonal matrices of the form

\[
X  = \begin{pmatrix} \Sigma & 0 & \cdots & 0 \\ 0 & Z_{S_1} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & Z_{S_m} \end{pmatrix}
\]

For such $X$ our constraints are equivalent to $X \succeq 0$,

\[
\langle E_{jj} - E_{11}, X \rangle = 0 \quad \text{ for } j =2,\ldots,d
\]

\[
\langle E_{jj'} + E_{S,jj'} , X \rangle = (\Sigma_S)_{jj'} \quad \text{ for } S \in \mathbb{S} \text{ and } j,j' \in S,
\]

where $E_{jj'} = (\boldsymbol{e}_j\boldsymbol{e}_{j'}^T +\boldsymbol{e}_{j'}\boldsymbol{e}_j^T)/2$ is the binary symmetric matrix of the same dimension as $X$ with its only non-zero entries being in the $(j,j')$-th and $(j',j)$-th positions of the top left block, and where $E_{S,jj'} = (\boldsymbol{e}_{S,j}\boldsymbol{e}_{S,j'}^T +\boldsymbol{e}_{S,j'}\boldsymbol{e}_{S,j}^T)/2$ is the binary symmetric matrix of the same dimension as $X$ with its only non-zero entries being in the $(j,j')$-th and $(j',j)$-th positions of the block occupied by $Z_S$ in $X$. Write $C$ for the diagonal matrix of the same dimension as $X$ with $I_d$ in the top left block, and all other entries equal to zero. It is now possible to write

\begin{align}
\label{Eq:StandardDual}
\sup&\{ \mathrm{tr}(\Sigma) : \Sigma \in \mathcal{P}^*,\quad \Sigma_{11} = \ldots = \Sigma_{dd}, \quad\Sigma_\mathbb{S} - A\Sigma \succeq_\mathbb{S} 0\} \nonumber \\
&= \sup\{ \langle C, X \rangle : X \succeq 0 \text{ is block diagonal}, \langle E_{jj} - E_{11}, X \rangle = 0 \text{ for } j =2,\ldots,d, \nonumber \\
& \hspace{150pt} \langle E_{jj'} + E_{S,jj'} , X \rangle = (\Sigma_S)_{jj'} \quad \text{ for } S \in \mathbb{S} \text{ and } j,j' \in S \},
\end{align}

so that our dual problem (<ref>) is now in standard form. Our primal problem (<ref>) is put into standard form by writing

\begin{align}
\label{Eq:StandardPrimal}
\inf \{ &\langle \Sigma_\mathbb{S}, Y_\mathbb{S} \rangle : A^* Y_\mathbb{S} + Y \succeq I_d, \quad Y_\mathbb{S} \succeq_\mathbb{S} 0, \quad Y_\mathbb{S} \in \mathcal{M}_\mathbb{S}, \quad Y \in \mathcal{Y} \} \nonumber \\
& =\inf\biggl\{ \sum_{S \in \mathbb{S}} \sum_{j,j' \in S} (\Sigma_S)_{jj'} y_{S,jj'} : \sum_{S \in \mathbb{S}} \sum_{j,j' \in S} y_{S,jj'} E_{jj'} + \sum_{j=2}^d y_{jj}(E_{jj}-E_{11}) \succeq I_d, \nonumber \\
& \hspace{100pt} \sum_{S \in \mathbb{S}} \sum_{j,j' \in S} y_{S,jj'} E_{S,jj'} \succeq 0 \text{ for all } S \in \mathbb{S}, \quad y_{jj},y_{S,jj'} \in \mathbb{R} \text{ for all } S,j,j' \biggr\}.
\end{align}

With the problems written in standard form, it is now clear that (<ref>) is the dual problem associated to (<ref>); see  Theorem 3.1 in [Vandenberghe and Boyd, 1996]. Observe further that the primal problem is strictly feasible since $Y_{\mathbb{S}} = X_{\mathbb{S}}^0$ satisfies the linear constraints with  $Y$ equal to the zero matrix. Hence, by standard duality results (Theorem 2.15 in [Blekherman et al., 2012], Theorem 3.1 in [Vandenberghe and Boyd, 1996]), we have that

\begin{align*}
%\sup&\{ \mathrm{tr}(\Sigma) : \Sigma \in \mathcal{P}^*,\quad \Sigma_{11} = \ldots = \Sigma_{dd}, \quad\Sigma_\mathbb{S} - A\Sigma \succeq_\mathbb{S} 0\} \\
& \sup\{ \langle C, X \rangle : X \succeq 0 \text{ is block diagonal}, \langle E_{jj} - E_{11}, X \rangle = 0 \text{ for } j =2,\ldots,d, \\
& \hspace{200pt} \langle E_{jj'} + E_{S,jj'} , X \rangle = (\Sigma_S)_{jj'} \quad \text{ for } S \in \mathbb{S} \text{ and } j,j' \in S \} \\
& = \inf\biggl\{ \sum_{S \in \mathbb{S}} \sum_{j,j' \in S} (\Sigma_S)_{jj'} y_{S,jj'} : \sum_{S \in \mathbb{S}} \sum_{j,j' \in S} y_{S,jj'} E_{jj'} + \sum_{j=2}^d y_{jj}(E_{jj}-E_{11}) \succeq I_d, \\
& \hspace{100pt} \sum_{S \in \mathbb{S}} \sum_{j,j' \in S} y_{S,jj'} E_{S,jj'} \succeq 0 \text{ for all } S \in \mathbb{S}, \quad y_{jj},y_{S,jj'} \in \mathbb{R} \text{ for all } S,j,j' \biggr\}, %\\
%& = \inf \{ \langle \Sigma_\mathbb{S}, Y_\mathbb{S} \rangle : A^* Y_\mathbb{S} + Y \succeq I_d, \quad Y_\mathbb{S} \succeq_\mathbb{S} 0, \quad Y_\mathbb{S} \in \mathcal{M}_\mathbb{S}, \quad Y \in \mathcal{Y} \},
\end{align*}

and the result follows.

(i) Convexity follows easily from basic properties of the supremum. Indeed, consider $\tilde{\Sigma}_\mathbb{S}: = \lambda \Sigma_\mathbb{S}^{(1)} + (1 - \lambda) \Sigma_\mathbb{S}^{(2)}$ with $\lambda \in [0,1]$.   Observe that $R$ is well defined at $\tilde{\Sigma}_\mathbb{S}$, as the convex combination of correlation matrices is still a correlation matrix. Then, for all $\lambda \in [0,1]$,

\begin{align*}
R(\tilde{\Sigma}_\mathbb{S}) &= \sup\biggr\{ -\frac{1}{d} \langle X_{\mathbb{S}}, \tilde{\Sigma}_{\mathbb{S}} \rangle : X_\mathbb{S} + X_\mathbb{S}^0 \succeq_\mathbb{S} 0, A^* X_\mathbb{S} + Y \succeq 0 \text{ for some } Y \in \mathcal{Y} \biggr\} \\
&= \sup\biggr\{ -\frac{1}{d} \langle X_{\mathbb{S}}, \lambda \Sigma_\mathbb{S}^{(1)} + (1 - \lambda) \Sigma_\mathbb{S}^{(2)} \rangle : X_\mathbb{S} + X_\mathbb{S}^0 \succeq_\mathbb{S} 0, A^* X_\mathbb{S} + Y \succeq 0 \text{ for some } Y \in \mathcal{Y} \biggr\} \\
&\leq \lambda \sup\biggr\{ -\frac{1}{d} \langle X_{\mathbb{S}},  \Sigma_\mathbb{S}^{(1)} \rangle: X_\mathbb{S} + X_\mathbb{S}^0 \succeq_\mathbb{S} 0, A^* X_\mathbb{S} + Y \succeq 0 \text{ for some } Y \in \mathcal{Y} \biggr\} \\
&  \hspace{50pt} + (1-\lambda) \sup\biggr\{ -\frac{1}{d} \langle X_{\mathbb{S}}, \Sigma_\mathbb{S}^{(2)} \rangle : X_\mathbb{S} + X_\mathbb{S}^0 \succeq_\mathbb{S} 0, A^* X_\mathbb{S} + Y \succeq 0 \text{ for some } Y \in \mathcal{Y} \biggr\} \\ &= \lambda R(\Sigma_\mathbb{S}^{(1)}) + (1-\lambda) R(\Sigma_\mathbb{S}^{(2)}),
\end{align*}

and the convexity of $R(\cdot)$ follows.

(ii) $R$ acts on $\mathcal{P}_{\mathbb{S}}$, which is the space of correlation matrices over the patterns $\mathbb{S}$. Now, the spectrahedron of all correlation matrices of dimension $p$,
\[\mathcal{E}_p = \left\{\left(x_1, \ldots, x_{\binom{p}{2}}\right) \in \mathbb{R}^{\binom{p}{2}} : \Sigma_{\boldsymbol{x}} = \begin{pmatrix}
1 & x_1 & \cdots & x_{p-1}\\
x_1 & 1 & \cdots & x_{2p-3}\\
\vdots & \vdots & \ddots & \vdots\\
x_{p-1} & x_{2p-3} & \cdots & 1\\
\end{pmatrix} \succeq 0 \right\},\]

is called the elliptope, and identifies a closed subset of $\mathbb{R}^{\binom{p}{2}}$. This follows from the fact that the symmetry condition $\Sigma_{\boldsymbol{x}} = \Sigma_{\boldsymbol{x}}^T$ defines a linear subspace of $\mathbb{R}^p$ of dimension $\binom{p}{2}$, while the PSD condition $v^T \Sigma_{\boldsymbol{x}} v \geq 0$ for all $v \in \mathbb{R}^p$ defines a closed subset of $\mathbb{R}^{\binom{p}{2}}$, which is a convex cone. For further insights, refer to [Laurent and Poljak, 1996]. This implies that, for every pattern $\mathbb{S}$,  $\mathcal{P}_{\mathbb{S}}$ can be identified with a closed subspace of $\mathbb{R}^s$, where $s = \sum_{S \in \mathbb{S}} \binom{|S|}{2}$. The continuity of $R$ follows from the fact that every convex function that is finite on  $\mathbb{R}^s$ is necessarily continuous (see Corollary 10.1.1. in [Rockafellar, 1970]).

To prove (iii), we will make use of the fact that the dual characterisation allows us to express $R(\Sigma_{\mathbb{S}'})$ as
\[
1 - \frac{1}{d'} \sup\{ \mathrm{tr}(\Sigma) : \Sigma \in \mathcal{P}^*, \Sigma_{11} = \ldots = \Sigma_{d'd'}, \Sigma_{\mathbb{S}'} - A_{\mathbb{S}'} \Sigma \succeq_{\mathbb{S}'} 0\},
\]

where $d'=|\cup_{S \in \mathbb{S}'}S|$. Now, let $\tilde{\Sigma}$ be an optimal feasible matrix for $\Sigma_{\mathbb{S}'}$, where all the diagonal elements of $\tilde{\Sigma}$ are the equal to each other by definition of $R$. Then, if we consider the restriction of $\tilde{\Sigma}$ on $\cup_{S \in \mathbb{S}}S$, call it  $\tilde{\Sigma}_{|\mathbb{S}}$, it is clear that $ \Sigma_{\mathbb{S}} - A_{\mathbb{S}} \tilde{\Sigma}_{|\mathbb{S}} \succeq_{\mathbb{S}} 0$, since $ \Sigma_{\mathbb{S}'} - A_{\mathbb{S}'} \Sigma \succeq_{\mathbb{S}'} 0$ and $\Sigma_{\mathbb{S}} \subseteq \Sigma_{\mathbb{S}'}$ by hypothesis, while $\tilde{\Sigma}_{|\mathbb{S}} \succeq 0$ follows again by Cauchy's interlacing theorem. Hence, calling $d=|\cup_{S \in \mathbb{S}}S|$, for every $\tilde{\Sigma}$ that is optimal for $\Sigma_{\mathbb{S}'}$, we can construct a feasible $\tilde{\Sigma}_{|\mathbb{S}}$ for $\Sigma_{\mathbb{S}}$ such that $1-\operatorname{tr}(\tilde{\Sigma}_{|\mathbb{S}})/d = R(\Sigma_{\mathbb{S}'})$, which completes the proof.

Let $\sigma_\mathbb{S}^2$ be a nonnegative sequence such that $\bar{\operatorname{av}}_j(\sigma^2_{\mathbb{S}}) = 1 \text{ for all } j \in [d]$. Using, for the third equality, the facts that $A_{V} \boldsymbol{1}_d$ also satisfies these properties and that $\bar{\text{av}}$ is linear, we have that

\begin{align*}
1 - V(\sigma_{\mathbb{S}}^2) &= 1 - \inf\left\{ \epsilon \in [0,1] : \sigma_{\mathbb{S}}^2 = (1-\epsilon)A_{V} \boldsymbol{1}_d + \epsilon {\sigma'}_{\mathbb{S}}^2 \text{ with } \bar{\operatorname{av}}_j(\sigma'^2_{\mathbb{S}}) = 1 \text{ for all } j \in [d] \right\} \\
& = \sup\left\{ \epsilon \in [0,1] : \sigma_{\mathbb{S}}^2 = \epsilon A_{V} \boldsymbol{1}_d + (1-\epsilon) {\sigma'}_{\mathbb{S}}^2 \text{ with } \bar{\operatorname{av}}_j(\sigma'^2_{\mathbb{S}}) = 1 \text{ for all } j \in [d] \right\} \\
& =  \sup\left\{ \epsilon \in [0,1] : \epsilon \leq \min_{j \in [d]} \min_{S \in \mathbb{S}_j} \sigma_{S,j}^2  \right\} = \min_{j \in [d]} \min_{S \in \mathbb{S}_j} \sigma_{S,j}^2,
\end{align*}

as claimed.

§.§ Proofs for Section <ref>

We are interested in finding $C_{\alpha} \in (0,1)$ such that $\forall \alpha \in (0,1)$

\[ \mathbb{P}_{H_0} \left(\hat{T} \geq C_{\alpha}  \right) = \mathbb{P}_{H_0} \left(R(\hat{\Sigma}_{\mathbb{S}}) + V(\hat{\sigma}_{\mathbb{S}}^2) \geq C_{\alpha}  \right) \leq \alpha. \]

First, observe that
\[
\mathbb{P}_{H_0} \left(R(\hat{\Sigma}_{\mathbb{S}}) + V(\hat{\sigma}_{\mathbb{S}}^2) \geq C_{\alpha}  \right) \leq \mathbb{P}_{H_0}\left(R(\hat{\Sigma}_{\mathbb{S}}) \geq C_{\alpha}/2  \right) + \mathbb{P}_{H_0}\left(V(\hat{\sigma}_{\mathbb{S}}^2) \geq C_{\alpha}/2  \right),
\]

and let us concentrate on the first term. For simplicity we replace our assumption that $\Sigma_{\mathbb{S}} \succeq_\mathbb{S} cI_\mathbb{S}$ by the assumption that $\Sigma_{\mathbb{S}} \succeq_\mathbb{S} 2cI_\mathbb{S}$. We have

\begin{align}
\label{Eq:NonSingularityDecomposition}
\mathbb{P}_{H_0}\left(R(\hat{\Sigma}_{\mathbb{S}}) \geq C_{\alpha}/2  \right) & \leq  \mathbb{P}_{H_0}\left(R(\hat{\Sigma}_{\mathbb{S}}) \geq C_{\alpha}/2,\, \hat{\Sigma}_{\mathbb{S}} \succeq_\mathbb{S} cI_\mathbb{S} \right) + 1 - \mathbb{P}_{H_0}\left(\hat{\Sigma}_{\mathbb{S}} \succeq_\mathbb{S} cI_\mathbb{S} \right).
\end{align}

Since $\Sigma_{\mathbb{S}} \succeq_\mathbb{S} 2cI_\mathbb{S}$ by assumption, we may bound the second part of (<ref>) by writing

\begin{align*}
1 & = \mathbb{P}_{H_0}\left(\Sigma_{\mathbb{S}} \succeq_\mathbb{S} 2cI_\mathbb{S} \right) = \mathbb{P}_{H_0}\left(\Sigma_{\mathbb{S}} - \hat{\Sigma}_{\mathbb{S}} + \hat{\Sigma}_{\mathbb{S}}  \succeq_\mathbb{S} 2cI_\mathbb{S} \right) \\
& \leq \mathbb{P}_{H_0}\left(\Sigma_{\mathbb{S}} - \hat{\Sigma}_{\mathbb{S}}  \succeq_\mathbb{S} cI_\mathbb{S} \right) + \mathbb{P}_{H_0}\left( \hat{\Sigma}_{\mathbb{S}}  \succeq_\mathbb{S} cI_\mathbb{S} \right) \\
& \leq \mathbb{P}_{H_0}\left(\|\hat{\Sigma}_{\mathbb{S}} - \Sigma_{\mathbb{S}}\|_{2,\mathbb{S}} \geq c \right) + \mathbb{P}_{H_0}\left( \hat{\Sigma}_{\mathbb{S}}  \succeq_\mathbb{S} cI_\mathbb{S} \right).
\end{align*}

This implies that

\[
1 - \mathbb{P}_{H_0}\left(\hat{\Sigma}_{\mathbb{S}} \succeq_\mathbb{S} cI_\mathbb{S} \right) \leq \mathbb{P}_{H_0}\left(\|\hat{\Sigma}_{\mathbb{S}} - \Sigma_{\mathbb{S}}\|_{2,\mathbb{S}} \geq c \right) .
\]

Now, define
\[
\bar{\mathrm{tr}}(X_\mathbb{S}) = \sum_{j=1}^d |\mathbb{S}_j|^{-1} \sum_{S \in \mathbb{S}_j} (X_S)_{jj}
\]

and observe that, for all $X_\mathbb{S} \in \mathcal{M}_\mathbb{S}$, we have $\langle X_\mathbb{S}^0, X_\mathbb{S} \rangle_\mathbb{S} = \bar{\mathrm{tr}}(X_\mathbb{S})$. See Proposition <ref> in Appendix <ref> for a proof of this fact. Using the arguments leading up to (<ref>) above, the first term on the right-hand side of (<ref>) can be written as

\[
\mathbb{P}_{H_0}\left(R(\hat{\Sigma}_{\mathbb{S}}) \geq C_{\alpha}/2, \,\hat{\Sigma}_{\mathbb{S}} \succeq_\mathbb{S} cI_\mathbb{S} \right) = \mathbb{P}_{H_0}\left(\sup\limits_{X_{\mathbb{S}} \in \mathcal{F}_c} -\frac{1}{d} \langle X_{\mathbb{S}}, \hat{\Sigma}_{\mathbb{S}} \rangle_{\mathbb{S}}  \geq C_{\alpha}/2, \,\hat{\Sigma}_{\mathbb{S}} \succeq_\mathbb{S} cI_\mathbb{S} \right),
\]

where $\mathcal{F}_c =  \{X_{\mathbb{S}} + X_{\mathbb{S}}^0 \succeq_{\mathbb{S}} 0 , A^* X_\mathbb{S} + Y \succeq 0 \text{ for some } Y \in \mathcal{Y}, \langle X_{\mathbb{S}} + X_{\mathbb{S}}^0, 2cI_{\mathbb{S}} \rangle_{\mathbb{S}} \leq d\}$. Discarding the condition $A^* X_\mathbb{S} + Y \succeq 0 \text{ for some } Y \in \mathcal{Y}$ and enlarging our feasible to   $\mathcal{\tilde{F}}_c := \{X_{\mathbb{S}} + X_{\mathbb{S}}^0 \succeq_{\mathbb{S}} 0 ,  \langle X_{\mathbb{S}} + X_{\mathbb{S}}^0, 2cI_{\mathbb{S}} \rangle_{\mathbb{S}} \leq d\}$, we have

\begin{align*}
\mathbb{P}_{H_0} \left(R(\hat{\Sigma}_{\mathbb{S}}) \geq C_{\alpha}/2, \,\hat{\Sigma}_{\mathbb{S}} \succeq_\mathbb{S} cI_\mathbb{S} \right) &= \mathbb{P}_{H_0} \left(\hat{R} - R \geq C_{\alpha}/2, \,\hat{\Sigma}_{\mathbb{S}} \succeq_\mathbb{S} cI_\mathbb{S} \right) \leq  \mathbb{P}_{H_0} \left(|\hat{R} - R| \geq C_{\alpha}/2, \,\hat{\Sigma}_{\mathbb{S}} \succeq_\mathbb{S} cI_\mathbb{S} \right) \\
& \leq \mathbb{P}_{H_0} \left(\sup\limits_{X_{\mathbb{S}} \in \mathcal{\tilde{F}}_c} \left| - \frac{1}{d} \langle X_{\mathbb{S}}, \hat{\Sigma}_{\mathbb{S}} - \Sigma_{\mathbb{S}}\rangle_{\mathbb{S}} \right| \geq C_{\alpha}/2\right) \\
&= \mathbb{P}_{H_0} \left(\sup\limits_{X_{\mathbb{S}} \in \mathcal{\tilde{F}}_c} \left| \langle X_{\mathbb{S}} + X_{\mathbb{S}}^0, \hat{\Sigma}_{\mathbb{S}} - \Sigma_{\mathbb{S}}\rangle_{\mathbb{S}} - (\operatorname{\bar{tr}}(\hat{\Sigma}_{\mathbb{S}}) - d)\right| \geq d \cdot C_{\alpha}/2 \right) \\
&\leq \mathbb{P}_{H_0} \left(\|\hat{\Sigma}_{\mathbb{S}} - \Sigma_{\mathbb{S}}\|_{2,\mathbb{S}} \cdot \sup\limits_{X_{\mathbb{S}} \in \mathcal{\tilde{F}}_c} \|X_{\mathbb{S}} + X_{\mathbb{S}}^0\|_{*, \mathbb{S}} \geq d \cdot C_{\alpha}/2 \right) \\
& \leq \mathbb{P}_{H_0} \left(\|\hat{\Sigma}_{\mathbb{S}} - \Sigma_{\mathbb{S}}\|_{2,\mathbb{S}} \cdot d/c \geq d \cdot C_{\alpha}/2 \right) \\
& = \mathbb{P}_{H_0} \left(\|\hat{\Sigma}_{\mathbb{S}} - \Sigma_{\mathbb{S}}\|_{2,\mathbb{S}} \geq c \cdot C_{\alpha}/2 \right),
\end{align*}

where we used Holder's inequality for sequences of matrices, and the fact that $\operatorname{\bar{tr}}(\hat{\Sigma}_{\mathbb{S}}) = d$, since $\hat{\Sigma}_{\mathbb{S}}$ is a sequence of sample correlation matrices. Putting all the pieces together, we have
\begin{align*}
\mathbb{P}_{H_0}\left(\hat{T}\geq C_{\alpha}  \right) & \leq \mathbb{P}_{H_0} \left(\|\hat{\Sigma}_{\mathbb{S}} - \Sigma_{\mathbb{S}}\|_{2,\mathbb{S}} \geq c \cdot C_{\alpha}/2 \right) + \mathbb{P}_{H_0} \left(\|\hat{\Sigma}_{\mathbb{S}} - \Sigma_{\mathbb{S}}\|_{2,\mathbb{S}} \geq c \right) + \mathbb{P}_{H_0}\left(V(\hat{\sigma}_{\mathbb{S}}^2) \geq C_{\alpha}/2  \right) \\
& \leq 2 \mathbb{P}_{H_0} \left(\|\hat{\Sigma}_{\mathbb{S}} - \Sigma_{\mathbb{S}}\|_{2,\mathbb{S}} \geq c \cdot C_{\alpha}/2 \right) + \mathbb{P}_{H_0}\left(V(\hat{\sigma}_{\mathbb{S}}^2) \geq C_{\alpha}/2  \right),
\end{align*}

since $\mathbb{P}(X \geq x_1) + \mathbb{P}(X \geq x_2) \leq 2\mathbb{P}(X \geq \min\{x_1, x_2\})$. Hence, in order to bound this probability above by $\alpha$, it is sufficient to find $C_{\alpha}$ such that

\begin{equation}
\label{Eq:MaxProbs}
\max\left\{\mathbb{P}_{H_0} \left(\|\hat{\Sigma}_{\mathbb{S}} - \Sigma_{\mathbb{S}}\|_{2,\mathbb{S}} \geq c \cdot C_{\alpha}/2 \right), \mathbb{P}_{H_0}\left(V(\hat{\sigma}_{\mathbb{S}}^2) \geq C_{\alpha}/2  \right) \right\} \leq \alpha/3.
\end{equation}

As for the first term inside the maximum, we have

\begin{align*}
\mathbb{P}_{H_0} \Bigl( \|\hat{\Sigma}_{\mathbb{S}} - \Sigma_{\mathbb{S}}\|_{2,\mathbb{S}} & \geq c \cdot C_{\alpha}/2 \Bigr)  \leq \mathbb{P}_{H_0} \left(\max_{S \in \mathbb{S}} \|\hat{\Sigma}_{S} - \Sigma_{S}\|_2 \geq c \cdot C_{\alpha}/2 \right) \\
& \leq \sum_{S \in \mathbb{S}} \mathbb{P}_{H_0} \left( \|\hat{\Sigma}_{S} - \Sigma_{S}\|_2 \geq c \cdot C_{\alpha}/2  \right) \leq |\mathbb{S}| \cdot \max_{S \in \mathbb{S}} \mathbb{P}_{H_0} \left( \|\hat{\Sigma}_{S} - \Sigma_{S}\|_2 \geq c \cdot C_{\alpha}/2  \right).
\end{align*}

Hence, calling

\begin{align*}
C_{t}(S) &:= C_1 \frac{\nu^2}{\sigma_\mathrm{min}^2} \left( \sqrt{\frac{|S| + \log(1/t)}{n}} \vee \frac{|S| + \log(1/t)}{n}\right) + C_2 \frac{\nu^4}{\sigma^4_{\mathrm{min}}} \sqrt{\frac{|S| \log(|S|/t)}{n}}  \\
& \quad \quad + C_3 \frac{\nu^4}{\sigma^4_{\mathrm{min}}} \left( \sqrt{\frac{|S| + \log(1/t)}{n}} \vee \frac{|S| + \log(1/t)}{n} \right) \sqrt{\frac{|S| \log(|S|/t)}{n}}
\end{align*}

for all $S \in \mathbb{S}$, with $C_1, C_2, C_3 > 0$ sufficiently big universal constants, it is immediate to see using Proposition <ref> that it is sufficient to take

\begin{equation}\label{eq:C_alpha_R}
C_{\alpha} \geq \frac{2}{c}\max_{S \in \mathbb{S}}C_{\alpha/3|\mathbb{S}|}(S) ,
\end{equation}

in order to have $\mathbb{P}_{H_0} \left(\|\hat{\Sigma}_{\mathbb{S}} - \Sigma_{\mathbb{S}}\|_{2,\mathbb{S}} \geq c \cdot C_{\alpha}/2 \right) \leq \alpha/3$. As for the second term in the maximum in (<ref>), since $V(\sigma_{\mathbb{S}}^2) = 0$ under the null,

\begin{align*}
\mathbb{P}_{H_0}\left(V(\hat{\sigma}_{\mathbb{S}}^2) \geq C_{\alpha}/2  \right) & = \mathbb{P}_{H_0}\left(V(\hat{\sigma}_{\mathbb{S}}^2) - V(\sigma_{\mathbb{S}}^2) \geq C_{\alpha}/2  \right) \leq  \mathbb{P}_{H_0}\left(|V(\hat{\sigma}_{\mathbb{S}}^2) - V(\sigma_{\mathbb{S}}^2)| \geq C_{\alpha}/2  \right) \\
& = \mathbb{P}_{H_0}\left(|\min_{j \in [d]} \min_{S \in \mathbb{S}_j} \hat{\sigma}^2_{S,j} - \min_{j \in [d]} \min_{S \in \mathbb{S}_j} \sigma^2_{S,j}| \geq C_{\alpha}/2  \right) \\
& \leq \mathbb{P}_{H_0}\left(\max_{j \in [d]} \max_{S \in \mathbb{S}_j} |\hat{\sigma}^2_{S,j} - \sigma^2_{S,j}| \geq C_{\alpha}/2  \right) \\
&  \leq \sum_{j \in [d]} \sum_{S \in \mathbb{S}_j} \mathbb{P}_{H_0}\left( |\hat{\sigma}^2_{S,j} - \sigma^2_{S,j}| \geq C_{\alpha}/2  \right) \\
& \leq \left(\sum_{j \in [d]} |\mathbb{S}_j| \right) \max_{j \in [d]} \max_{S \in \mathbb{S}_j} \mathbb{P}_{H_0}\left( |\hat{\sigma}^2_{S,j} - \sigma^2_{S,j}| \geq C_{\alpha}/2  \right).
\end{align*}

Now,  the standard Chernoff method gives, for all $j \in [d]$, for all $S \in \mathbb{S}_j$,

\begin{align*}
\mathbb{P}_{H_0}\left( |\hat{\sigma}^2_{S,j} - \sigma^2_{S,j}| \geq C_{\alpha}/2  \right)
& \leq \exp\left\{-\frac{n_{S} C^2_{\alpha}}{256 \nu^4} \right\} \leq \exp\left\{-\frac{\min_{S \in \mathbb{S}} n_S C^2_{\alpha}}{256 \nu^4} \right\},
\end{align*}

so that $\mathbb{P}_{H_0}\left(V(\hat{\sigma}_{\mathbb{S}}^2) \geq C_{\alpha}/2  \right) \leq \alpha/3$ is satisfied if
\[
\left(\sum_{j \in [d]} |\mathbb{S}_j| \right) \exp\left\{-\frac{\min_{S \in \mathbb{S}} n_S C^2_{\alpha}}{256 \nu^4} \right\} \leq \frac{\alpha}{3}.
\]

Hence, it is sufficient to take

\begin{equation}\label{eq:C_alpha_V}
C_{\alpha} \geq 16 \nu^2 \sqrt{\frac{\log\left( 3 \sum_{j \in [d]} |\mathbb{S}_j|/\alpha\right)}{\min_{S \in \mathbb{S}} n_S}}.
\end{equation}

In order to satisfy both (<ref>) and (<ref>) at the same time, it is sufficient to take the maximum between the two right-hand sides, and the statement follows.

First, observe that a generic element of $\hat{P} - P$, can be written as

\[
\hat{P}_{ij} - P_{ij} = \frac{\hat{\sigma}_{ij}}{\hat{\sigma}_{i}\hat{\sigma}_{j}} - \frac{\sigma_{ij}}{\sigma_i \sigma_j} = \left(\frac{1}{\hat{\sigma}_{i}\hat{\sigma}_{j}} - \frac{1}{\sigma_i \sigma_j} \right)\hat{\sigma}_{ij} + \frac{\hat{\sigma}_{ij} - \sigma_{ij}}{\sigma_i \sigma_j}.
\]

This implies that $\hat{P} - P$ can be written as
\[
\hat{P} - P = W \circ \hat{\Sigma} + \tilde{W} \circ (\hat{\Sigma} - \Sigma),
\]

where $\circ$ stands for the matrix pointwise product, also known as Hadamard product, and $W, \tilde{W}$ satisfy $W_{ij} = 1/\hat{\sigma}_{i}\hat{\sigma}_{j} - 1/\sigma_i \sigma_j$ and $\tilde{W}_{ij} = 1/\sigma_i \sigma_j$. We will now bound the operator norm of this difference using the following facts about Hadamard products. First, as shown in (3.7.12) of Johnson, 1989, if $A,B \in \mathcal{M}$ and $A \succeq 0$ then we have

\[
\|A \circ B\|_2 \leq \max_{j \in [d]} |A_{jj}| \|B\|_2.
\]

Second, for arbitrary $A,B \in \mathcal{M}$ and $v=\sum_{i=1}^d v_i \mathbf{e}_i \in \mathbb{R}^d$, we have

\begin{align*}
\|(A \circ B) v \| & = \|(A \circ B)  \sum_{i = 1}^d v_i\boldsymbol{e}_i \| \leq \sum_{i = 1}^d |v_i| \|(A \circ B)e_i\| \\
& = \sum_{i = 1}^d |v_i| \|(A_{1i} B_{1i}, \ldots, A_{di} B_{di})^T\| \leq \sum_{i = 1}^d |v_i| \left(\max_{j \in [d]} |A_{ji}|\right) \|(B_{1i}, \ldots, B_{di})^T\| \\
& \leq \left(\max_{i,j \in [d]} |A_{ij}|\right) \sum_{i = 1}^d |v_i| \|B\boldsymbol{e}_i\| \leq \left(\max_{i,j \in [d]} |A_{ij}|\right) \sum_{i = 1}^d |v_i| \|B\|_2 \leq \sqrt{d} \left(\max_{i,j \in [d]} |A_{ij}|\right) \|B\|_2 \|v\|,
\end{align*}

where $\{\mathbf{e}_i\}_{i \in [d]}$ is the canonical basis for $\mathbb{R}^d$. Since $\tilde{W} \succeq 0$ these facts imply that

\begin{align*}
&\|\hat{P} - P\|_2  \leq \|W \circ \hat{\Sigma}\|_2 + \|\tilde{W} \circ (\hat{\Sigma} - \Sigma)\|_2 \leq \|W \circ \hat{\Sigma}\|_2 + \left(\max\limits_{j \in [d]} |\tilde{W}_{jj}| \right)\|\hat{\Sigma} - \Sigma\|_2 \\
& \leq \sqrt{d} \left(\max\limits_{i,j \in [d]} |W_{ij}| \right)\|\hat{\Sigma}\|_2 + \left(\max\limits_{j \in [d]} |\tilde{W}_{jj}| \right)\|\hat{\Sigma} - \Sigma\|_2 \\
& \leq \sqrt{d} \left(\max\limits_{i,j \in [d]} |W_{ij}| \right)\|\Sigma\|_2 +  \sqrt{d} \left(\max\limits_{i,j \in [d]} |W_{ij}| \right)\|\hat{\Sigma} - \Sigma\|_2 + \left(\max\limits_{j \in [d]} |\tilde{W}_{jj}| \right)\|\hat{\Sigma} - \Sigma\|_2 \\
& = \sqrt{d} \left(\max\limits_{i,j \in [d]} \left|\frac{1}{\hat{\sigma}_{i}\hat{\sigma}_{j}} - \frac{1}{\sigma_i \sigma_j} \right| \right)\|\Sigma\|_2 + \sqrt{d} \left(\max\limits_{i,j \in [d]} \left|\frac{1}{\hat{\sigma}_{i}\hat{\sigma}_{j}} - \frac{1}{\sigma_i \sigma_j} \right| \right)\|\hat{\Sigma} - \Sigma\|_2 + \left(\max\limits_{j \in [d]} \left|\frac{1}{\sigma_j^2}\right| \right)\|\hat{\Sigma} - \Sigma\|_2.
\end{align*}

For the second and final terms,  Proposition <ref> in Appendix <ref> ensures that
\[
\mathbb{P}\left(\|\hat{\Sigma}-\Sigma\|_2 \geq x\right) \leq 2 \cdot 9^d \exp \left\{-n \frac{x}{16 \nu^2} \wedge\left(\frac{x}{16 \nu^2}\right)^2  \right\},
\]

and inverting this bound, we get that
\[
\|\hat{\Sigma} - \Sigma\|_2 \leq C_1 \nu^2 \sqrt{\frac{d + \log(4/t)}{n}} \vee \frac{d + \log(4/t)}{n},
\]

with probability $\geq 1 - t/2$, for a universal constant $C_1 > 0$ sufficiently big. In this regard, we remark that universal constants might change from line to line, but we will use the same notation to ease the presentation. As for the first and second terms, writing $\sigma_{\min}^2$ for the smallest component of $(\sigma_1^2,\ldots,\sigma_d^2)$, observe that

\begin{align*}
\max\limits_{i,j \in [d]} \left|\frac{1}{\hat{\sigma}_{i}\hat{\sigma}_{j}} - \frac{1}{\sigma_i \sigma_j} \right| & = \max\limits_{i,j \in [d]} \frac{1}{\sigma_i \sigma_j} \left|\frac{\sigma_i \sigma_j}{\hat{\sigma}_{i}\hat{\sigma}_{j}} - 1 \right| \leq  \frac{1}{\sigma_\mathrm{min}^2} \max\limits_{i,j \in [d]} \left|\frac{\sigma_i \sigma_j}{\hat{\sigma}_{i}\hat{\sigma}_{j}} - 1 \right| \\
& = \frac{1}{\sigma_\mathrm{min}^2} \max\limits_{i,j \in [d]} \left|\frac{\sigma_i \sigma_j}{\hat{\sigma}_{i}\hat{\sigma}_{j}} - \frac{\sigma_i}{\hat{\sigma}_{i}} + \frac{\sigma_i}{\hat{\sigma}_{i}} - 1 \right| \\
& \leq \frac{1}{\sigma_\mathrm{min}^2} \left( \max\limits_{i,j \in [d]} \left|\frac{\sigma_j}{\hat{\sigma}_{j}}\left(\frac{\sigma_i}{\hat{\sigma}_{i}} - 1 \right)\right| + \max\limits_{i \in [d]} \left|\frac{\sigma_i}{\hat{\sigma}_{i}} - 1 \right|\right) \\
& = \frac{1}{\sigma_\mathrm{min}^2} \left( \max\limits_{i \in [d]} \left|\frac{\sigma_i}{\hat{\sigma}_{i}} - 1 \right| \max\limits_{j \in [d]} \left|\frac{\sigma_j}{\hat{\sigma}_{j}}\right| + \max\limits_{i \in [d]} \left|\frac{\sigma_i}{\hat{\sigma}_{i}} - 1 \right|\right) \\
& = \frac{1}{\sigma_\mathrm{min}^2} \left(\|D^{1/2}\hat{D}^{-1/2} - I\|_2 \|D^{1/2}\hat{D}^{-1/2}\|_2 + \|D^{1/2}\hat{D}^{-1/2} - I\|_2 \right),
\end{align*}

which implies that, in order to control $\max_{i,j \in [d]} \left|1/\hat{\sigma}_{i}\hat{\sigma}_{j} - 1/\sigma_i \sigma_j \right|$, it is enough to control $\|D^{1/2}\hat{D}^{-1/2} - I\|_2$. To this aim, first observe that, for all $t \in [0,1]$,

\[
\mathbb{P}(|\sigma_j^2/\hat{\sigma}_{j}^2 - 1| > t) \leq  \mathbb{P}(|\hat{\sigma}_j^2/\sigma_{j}^2 - 1| > t/2).
\]

\begin{align*}
\mathbb{P}\left(|\sigma_j^2/\hat{\sigma}_{j}^2 - 1| > t\right) & = \mathbb{P}\left(\left\{\sigma_j^2/\hat{\sigma}_{j}^2 > 1+t \right\} \cup \left\{\sigma_j^2/\hat{\sigma}_{j}^2 < 1-t\right\}\right) \\
& \leq \mathbb{P}\left(\sigma_j^2/\hat{\sigma}_{j}^2 > 1+t \right) + \mathbb{P}\left(\sigma_j^2/\hat{\sigma}_{j}^2 < 1-t \right) \\
& = \mathbb{P}\left(\hat{\sigma}_{j}^2/\sigma_j^2 < (1+t)^{-1} \right) + \mathbb{P}\left(\hat{\sigma}_{j}^2/\sigma_j^2 > (1-t)^{-1} \right) \\
& = \mathbb{P}\left(\hat{\sigma}_{j}^2 - \sigma_j^2 < -\sigma_j^2t(1+t)^{-1} \right) + \mathbb{P}\left(\hat{\sigma}_{j}^2 - \sigma_j^2 > \sigma_j^2t(1-t)^{-1}\right) \\
& = \mathbb{P}\left(\hat{\sigma}_{j}^2/\sigma_j^2 - 1 < -t(1+t)^{-1} \right) + \mathbb{P}\left(\hat{\sigma}_{j}^2/\sigma_j^2 - 1  > t(1-t)^{-1}\right),
\end{align*}

and since $t(1+t)^{-1}, t(1-t)^{-1} \geq t/2$ for $t \in [0,1]$, we can conclude
\[
\mathbb{P}\left(|\sigma_j^2/\hat{\sigma}_{j}^2 - 1| > t) \leq  \mathbb{P}(|\hat{\sigma}_j^2/\sigma_{j}^2 - 1| > t/2\right).
\]

This is helpful, since we know how to control $|\hat{\sigma}_j^2/\sigma_{j}^2 - 1|$. Indeed, since the $\mathbf{X}_i$ are $\nu$-subgaussian random vectors, the $X_{ij}$ for varying $i \in [n]$ are i.i.d. $\nu$-subgaussian random variables, hence the $X_{ij}^2$ for varying $i \in [n]$ are i.i.d. subexponential with parameters $(4\sqrt{2}\nu^2, 4\nu^2)$ (see Proposition <ref> in the Appendix). For every $j \in [d]$ and every $t >0$, we can therefore use the standard Chernhoff method for subexponential random variables (see Proposition <ref> in the Appendix) to see that
\[
\mathbb{P}\left(|\hat{\sigma}_j^2/\sigma_{j}^2 - 1| > t\right) = \mathbb{P}\left(\left|\frac{1}{n}\sum_{i \in [n]} (X_{ij}^2 - \mathbb{E}[X_{ij}^2])\right| > \sigma_j^2 t \right) \leq \exp{\{-n \sigma_j^4 t^2/(64 \nu^4) \}}
\]

if $0 \leq \sigma_j^2 t \leq 8 \nu^2.$ In particular, this remains valid for $0 \leq t \leq 8$, being $\sigma_j^2 \leq \nu^2$, and since we are interested in small values of $t$, we are allowed to focus just on this subgaussian regime in our analysis. Inverting the previous bound, we have that

\begin{align*}
1- \mathbb{P} & \left(|\sigma_j^2/\hat{\sigma}_{j}^2 - 1| \leq \frac{16 \nu^2}{\sigma_\mathrm{min}^2}\sqrt{\frac{\log(4d/t)}{n}}, \forall j \in [d] \right)  = \mathbb{P} \left(\exists j \in [d]: |\sigma_j^2/\hat{\sigma}_{j}^2 - 1| \leq \frac{16 \nu^2}{\sigma_\mathrm{min}^2}\sqrt{\frac{\log(4d/t)}{n}} \right) \\
& \leq \sum_{j=1}^d \mathbb{P}\left(|\sigma_j^2/\hat{\sigma}_{j}^2 - 1| \leq \frac{16 \nu^2}{\sigma_j^2}\sqrt{\frac{\log(4d/t)}{n}}\right) \leq \sum_{j=1}^d \mathbb{P}\left(|\hat{\sigma}_j^2/\sigma_{j}^2 - 1| \leq \frac{8 \nu^2}{\sigma_j^2}\sqrt{\frac{\log(4d/t)}{n}}\right) \leq t.
\end{align*}

Now, following similar lines as in [Oliveira, 2010], since for any $x \in [-3/4, 3/4]$ we have $|\sqrt{1+x} - 1| \leq x$ by the mean value theorem, if we assume $9 \sigma^4_{\mathrm{min}} n \geq 1024 \nu^4 \log(4d/t)$ and take $x = \sigma_j^2/\hat{\sigma}_{j}^2 - 1$, we also have

\[
\mathbb{P}\left(|\sigma_j/\hat{\sigma}_{j} - 1| \leq \frac{16 \nu^2}{\sigma_\mathrm{min}^2}\sqrt{\frac{\log(4d/t)}{n}}, \forall j \in [d] \right) \geq 1 - t,
\]

which leads to

\[
\mathbb{P}\left(\|D^{1/2}\hat{D}^{-1/2}-I\|_2 \leq \frac{16 \nu^2}{\sigma_\mathrm{min}^2}\sqrt{\frac{\log(4d/t)}{n}} \right) \geq 1 - t.
\]

It follows that

\begin{align*}
\|D^{-1/2}\hat{D}^{1/2} &- I\|_2 \|D^{-1/2}\hat{D}^{1/2}\|_2 + \|D^{-1/2}\hat{D}^{1/2} - I\|_2 \\
&\leq \frac{16 \nu^2}{\sigma_\mathrm{min}^2}\sqrt{\frac{\log(4d/t)}{n}}\left(1+\frac{16 \nu^2}{\sigma_\mathrm{min}^2}\sqrt{\frac{\log(4d/t)}{n}}\right) + \frac{16 \nu^2}{\sigma_\mathrm{min}^2}\sqrt{\frac{\log(4d/t)}{n}}
\end{align*}

with probability $\geq 1 - t$, and since for any $x \in [0, 1/2]$ we have $x(1+x)+x \leq 5x/2$, assuming  $\sigma^4_{\mathrm{min}} n \geq 1024 \nu^4 \log(4d/t)$, we obtain
\[
\|D^{-1/2}\hat{D}^{1/2} - I\|_2 \|D^{-1/2}\hat{D}^{1/2}\|_2 + \|D^{-1/2}\hat{D}^{1/2} - I\|_2 \leq \frac{40 \nu^2}{\sigma_\mathrm{min}^2}\sqrt{\frac{\log(4d/t)}{n}}
\]

with probability $\geq 1 - t$. This allows to conclude that

\[
\|D^{-1/2}\hat{D}^{1/2} - I\|_2 \|D^{-1/2}\hat{D}^{1/2}\|_2 + \|D^{-1/2}\hat{D}^{1/2} - I\|_2 \leq \frac{40 \nu^2}{\sigma_\mathrm{min}^2}\sqrt{\frac{\log(8d/t)}{n}} \leq C_2 \frac{\nu^2}{\sigma_\mathrm{min}^2}\sqrt{\frac{\log(d/t)}{n}}
\]

with probability $\geq 1 - t/2$, for a universal constant $C_2 > 0$ sufficiently big. This implies that

\begin{align*}
&\|\hat{P} - P\|_2 \\