# Post-clustering Inference under Dependence

Javier González-Delgado, Mathis Deronzier, Juan Cortés, Pierre Neuvial

► To cite this version:

Javier González-Delgado, Mathis Deronzier, Juan Cortés, Pierre Neuvial. Post-clustering Inference under Dependence. 2025. hal-04250364v2

HAL Id: hal-04250364

https://laas.hal.science/hal-04250364v2

Preprint submitted on 15 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Post-clustering Inference under Dependence

Javier González-Delgado¹,²,³, Mathis Deronzier², Juan Cortés³ and Pierre Neuvial²

¹Université de Rennes, ENSAI, CNRS, CREST-UMR 9194, F-35000 Rennes, France.

²Institut de Mathématiques de Toulouse, UMR5219, Université de Toulouse, CNRS, UPS, F-31062 Toulouse Cedex 9, France.

³LAAS-CNRS, Université de Toulouse, CNRS, F-31400 Toulouse, France.

## Abstract

Recent work by Gao et al. [21] has laid the foundations for post-clustering inference, establishing a theoretical framework allowing to test for differences between means of estimated clusters. Additionally, they studied the estimation of unknown parameters while controlling the selective type I error. However, their theory was developed for independent observations identically distributed as p-dimensional Gaussian variables, where the parameter estimation could only be performed for spherical covariance matrices. Here, we aim at extending this framework to a more convenient scenario for practical applications, where arbitrary dependence structures between observations and features are allowed. We establish sufficient conditions for extending the setting presented in [21] to the general dependence framework. Moreover, we assess theoretical conditions allowing the compatible estimation of a covariance matrix. The theory is developed for hierarchical agglomerative clustering algorithms with several types of linkages, and for the k-means algorithm. We illustrate our method with synthetic data and real data of protein structures.

## 1 Introduction

Post-selection inference has gained substantial attention in recent years due to its potential to address practical problems in diverse fields. The issue of using data to answer a question that has been chosen based on the same data was formalized in [20], where the basis of selective hypothesis testing was rigorously set with the definition of the selective type I error. This paved the way to perform selective testing when null hypotheses are chosen through clustering algorithms, bypassing the naive data splitting that reveals unsuitable in this context. However, their proposed approach, referred to as data carving, as well as more recent approaches like data fission [31] are difficult to implement in practice because they require knowledge of the covariance structure between variables. Moreover, they often involve the non-trivial calibration of a tuning parameter that controls the proportion of information allocated for model selection and for inference. The seminal work by Gao et al. [21] established a theoretical framework allowing selective testing after clustering using all the information in the data set. Their method is defined for independent observations identically distributed as p-dimensional Gaussian random variables with a spherical covariance matrix. This corresponds to the following matrix normal model [23]:

$$X \sim \mathcal{MN}_{n \times p}(\mu, I_n, \sigma^2 I_p), \qquad (\text{ind-MN})$$

where $\mu \in \mathcal{MN}_{n \times p}(\mathbb{R})$ and $\sigma > 0$. Under (ind-MN), the authors in [21] defined a p-value that controls the selective type I error when testing for a difference in means between a pair of estimated clusters. This p-value can be efficiently computed for hierarchical clustering algorithms with common linkage functions. Moreover, the authors in [21] made another remarkable contribution by addressing the estimation of $\sigma$ while controlling the selective type I error, which had not been addressed in previous works [31, 41] despite its major importance in applications. They showed that if $\sigma$ is asymptotically over-estimated, the p-value is asymptotically super-uniform under the null, and provided an estimator $\hat{\sigma}$ that can be used in practice. They also proposed an extension of their testing procedure to known arbitrary covariance
<figure><img src="image_2.png" /><figcaption>Figure 1: Empirical cumulative distribution functions (ECDF) of *p*-values defined in [21] testing for the difference in cluster means after performing a hierarchical clustering algorithm (HAC) with average linkage. The ECDF were computed from *M* = 2000 realizations of a matrix normal model with **μ** = **0**<sub>*n*×*p*</sub> and non-diagonal **U** and **Σ**. For each realization, the test compared the means of two randomly selected clusters after setting the HAC to choose three clusters. We set *n* = 100 and *p* ∈ {5, 20, 50}. In (a), dependence between observations is ignored. In (b), the covariance between features is assumed to be spherical to allow its estimation using the approach in [21].</figcaption></figure>

structures between features, still assuming i.i.d. observations. However, the estimation of the covariance
between features remained unaddressed.

Despite the notable contribution of [21], model (ind-MN) is somewhat limited in view of more complex
applications. In real problems, features describing observations are unlikely to be independent with
identical variance, but rather present more general covariance structures Σ. In the same way, observations
may present non-negligible dependence structures when, for instance, they can be drawn from time series
models or simulated with physical models involving time evolution. Note that ignoring dependence
between features and observations implies the loss of selective type I error control. This can be illustrated
by a simple simulation scenario based on matrix normal samples with non-diagonal covariance matrices,
accounting for the dependence structures between observations and/or features. If model assumptions
are not satisfied, the approach presented in [21] does not control the selective type I error, as illustrated
in Figure 1 by the fact that the distribution of the corresponding p-values is above the diagonal (which
corresponds to uniform p-values). This deviation from uniformity increases with the dimension of the
feature space. Details about the corresponding simulation are given in Appendix D.1.

The practical motivation of the present work is to perform inference after clustering protein conforma-
tions. Protein structures are non-static and their conformational variability is essential to understand the
relationship between sequence, structural properties and function [28]. Due to the high complexity of the
conformational space, clustering techniques have emerged as powerful tools to characterize the structural
variability of proteins, by extracting families of representative states [3, 11, 39, 43]. Usually, Euclidean
distances between pairs of amino acids are considered as *p*-dimensional descriptors of protein conforma-
tions [6, 11, 30]. These distances are highly correlated and hardly match the model (ind-MN). Moreover,
protein data is often simulated with Molecular Dynamics approaches that simulate the time-evolution of
the protein according to physical models [2]. In that case, independence between observations cannot be
assumed.
<figure><img src="image_3.png" /><figcaption>Figure 2: (a): Sample drawn from (gen-MN) with $n = 100$, $p = 2$ and non-diagonal covariance matrices $\mathbf{U}$ and $\mathbf{\Sigma}$. The mean matrix $\boldsymbol{\mu}$ is divided into three clusters. Observations are classified into three groups by the $k$-means algorithm and the difference between cluster means is tested using the approach proposed in this work. Classification using hierarchical agglomerative clustering (HAC) with average linkage yielded the same partition. (b): Data in (a) whitened and classified into three groups by the $k$-means algorithm. The differences between cluster means are tested assuming (ind-MN) and using the approach presented in [9]. (c): Data in (a) whitened and classified into three groups using HAC with average linkage. The differences between cluster means are tested assuming (ind-MN) and using the approach proposed in [21]. In all panels, $p_{ij}$ denotes the $p$-value for the difference between the means of clusters $i$ and $j$, for $i, j = 1, 2, 3$.</figcaption></figure>

Accordingly, our aim is to go one step further and extend the framework introduced in [21] to a more general setting where arbitrary dependence structures between both observations and features are admitted, allowing for the estimation of one of them. We present a generalization of [21] where the model (ind-MN) is extended to

$$ \mathbf{X} \sim \mathcal{MN}_{n \times p}(\boldsymbol{\mu}, \mathbf{U}, \mathbf{\Sigma}), \qquad (\text{gen-MN}) $$

where $\boldsymbol{\mu} \in \mathcal{MN}_{n \times p}(\mathbb{R})$, $\mathbf{U} \in \mathcal{MN}_{n \times n}(\mathbb{R})$ and $\mathbf{\Sigma} \in \mathcal{MN}_{p \times p}(\mathbb{R})$. Our techniques follow the same reasoning steps as the ones in [21], establishing sufficient conditions that allow an extension to (gen-MN).

The reader might wonder whether it is necessary to develop a new framework for (gen-MN) given that the matrix normal data can be whitened to fit into the Gao et al. model. Indeed, as we have

$$ \mathbf{X} \sim \mathcal{MN}_{n \times p}(\boldsymbol{\mu}, \mathbf{U}, \mathbf{\Sigma}) \Leftrightarrow \mathrm{vec}(\mathbf{X}) \sim \mathcal{N}_{np}(\mathrm{vec}(\boldsymbol{\mu}), \mathbf{\Sigma} \otimes \mathbf{U}), $$

the transformed random vector $(\mathbf{\Sigma} \otimes \mathbf{U})^{-\frac{1}{2}} \mathrm{vec}(\mathbf{X})$ has covariance matrix $\mathbb{I}_{np}$ and can be de-vectorized to fit (ind-MN). However, clustering the original and whitened data often leads to different partitions, and thus to different null hypotheses. In some cases, de-correlating the observations and features of $\mathbf{X}$ might yield a misleading impression of the underlying class structure. This is illustrated in Figure 2, where we show that whitening a sample drawn from (gen-MN) and performing a selective test defined for (ind-MN) might substantially alter the significance of the differences between cluster means, as well as the overall clustering partition. Details on this numerical analysis are provided in Appendix D.2. Note also that whitening a $n \times p$ matrix normal sample involves the inversion of a $np \times np$ matrix. These considerations, together with the unsuitability of the whitening approach when any of the covariance matrices is unknown, justifies the need of developing a new framework for the general model (gen-MN).

The paper is organized as follows:
* Section 2 presents our extension of [21] to the general model (gen-MN) when both covariance matrices are known.

* In Section 3, we explore the scenarios that allow the asymptotic over-estimation of either $\mathbf{U}$ or $\mathbf{\Sigma}$ while respecting the asymptotic control of the selective type I error. We provide an estimator that can be used in several common practical scenarios.

* Section 4 illustrates all the results through numerical experiments on synthetic data, and evaluates the robustness of the presented approach to model misspecification. Finally, Section 5 shows how this theory can be applied to perform inference after clustering protein structures.

# 2 Selective inference for clustering under general dependence

In [21], the authors consider the problem of selective inference after hierarchical clustering in the case of independent observations and features (with an extension to arbitrary known dependence between features). Here, we aim to extend the method to allow for general dependence structures between both observations and features. We consider $n$ observations of $p$ features drawn from the matrix normal distribution (gen-MN), where $\mathbf{U}$ and $\mathbf{\Sigma}$ are required to be positive definite. Each row of $\mathbf{X}$ is a vector of features in $\mathbb{R}^p$. The dependence between such features is given by $\mathbf{\Sigma}$, and $\mathbf{U}$ encodes the dependence between observations. If observations are independent with unit variance, we have $\mathbf{U} = \mathbf{I}_n$, and if features are independent with equal variance we can write $\mathbf{\Sigma} = \sigma^2\mathbf{I}_p$ for a given $\sigma > 0$. These two assumptions define the model in [21], which we aim to extend to the most general $\mathbf{U}$ and $\mathbf{\Sigma}$.

## 2.1 Problem setting and Gao et al.'s approach

Let us first recall the setting originally introduced in [21]. We will denote by $X_i$ (resp. $\mu_i$) the $i$-th row of $\mathbf{X}$ (resp. $\mu$) and, for a group of observations $\mathcal{G} \subseteq [n] = \{1, \dots, n\}$, $X_{\mathcal{G}}$ will denote the submatrix of $\mathbf{X}$ with rows $X_i$ for $i \in \mathcal{G}$. We also consider the mean of $\mathcal{G}$ in $\mathbf{X}$ and its empirical counterpart, denoted respectively by

$$ \bar{\mu}_{\mathcal{G}} = \frac{1}{|\mathcal{G}|} \sum_{i \in \mathcal{G}} \mu_i \quad \text{and} \quad \bar{X}_{\mathcal{G}} = \frac{1}{|\mathcal{G}|} \sum_{i \in \mathcal{G}} X_i. \qquad (1) $$

Letting

$$ C_{[n]} = \{(\mathcal{G}_1, \mathcal{G}_2), \mathcal{G}_1, \mathcal{G}_2 \subset [n] : \mathcal{G}_1 \cap \mathcal{G}_2 = \emptyset\}, \qquad (2) $$

be the set of all pairs of non-overlapping groups of observations, for any $(\mathcal{G}_1, \mathcal{G}_2) \in C_{[n]}$ we can define the column vector $\nu_{\mathcal{G}_1, \mathcal{G}_2}$ having as components

$$ [\nu_{\mathcal{G}_1, \mathcal{G}_2}]_i = 1\{i \in \mathcal{G}_1\} / |\mathcal{G}_1| - 1\{i \in \mathcal{G}_2\} / |\mathcal{G}_2|, \qquad (3) $$

for $i \in [n]$. This allows the difference between the (empirical) group means to be written compactly as

$$ \bar{\mu}_{\mathcal{G}_1} - \bar{\mu}_{\mathcal{G}_2} = \mu^T \nu_{\mathcal{G}_1, \mathcal{G}_2}, \quad \text{and} \quad \bar{X}_{\mathcal{G}_1} - \bar{X}_{\mathcal{G}_2} = \mathbf{X}^T \nu_{\mathcal{G}_1, \mathcal{G}_2}. \qquad (4) $$

For the sake of a clearer notation, we will simply write $\nu = \nu_{\mathcal{G}_1, \mathcal{G}_2}$ when the context is clear. Let $\mathcal{C}$ be a clustering algorithm, $\mathbf{x}$ a realization of $\mathbf{X}$ and $\mathcal{G}_1, \mathcal{G}_2$ an arbitrary pair of clusters in $\mathcal{C}(\mathbf{x})$. The goal of post-clustering inference is to assess the null hypothesis

$$ H_0^{\{\mathcal{G}_1, \mathcal{G}_2\}} : \mu^T \nu_{\mathcal{G}_1, \mathcal{G}_2} = 0, \qquad (H0) $$
by controlling the selective type I error for clustering at level $\alpha$, i.e. by ensuring that

$$ \mathbb{P}_{H_0^{\{G_1, G_2\}}} \left( \text{reject } H_0^{\{G_1, G_2\}} \text{ based on } \mathbf{X} \text{ at level } \alpha \mid \mathcal{G}_1, \mathcal{G}_2 \in \mathcal{C}(\mathbf{X}) \right) \le \alpha \quad \forall \alpha \in (0, 1). \quad (5) $$

If the inequality in the previous equation can be replaced by an equality, we will say that the selective type I error is controlled exactly at level $\alpha$. The ideal scenario to define a p-value for (H0) satisfying (5) would be to only condition on the event $\{G_1, G_2 \in \mathcal{C}(\mathbf{X})\}$, which is the broader conditioning set that allows selective type I error control. However, making the p-value analytically tractable often needs the refinement of the conditioning set by adding extra technical events (see also Appendix B).

The idea in [21] is to decompose $\mathbf{X}$ using the projection onto the orthogonal complement of $\nu$, that is, $\pi_{\nu}^{\perp} = \mathbf{I}_n - \nu\nu^{T}/||\nu||_2^2$. This naturally brings out the difference between empirical cluster means $\mathbf{X}^T\nu$, which can be used as a test statistic to evaluate (H0):

$$ \mathbf{X} = \pi_{\nu}^{\perp} \mathbf{X} + (\mathbf{I}_n - \pi_{\nu}^{\perp} \mathbf{X}) = \pi_{\nu}^{\perp} \mathbf{X} + \left( \frac{||\mathbf{X}^T \nu||_2}{||\nu||_2^2} \right) \nu \operatorname{dir}(\mathbf{X}^T \nu)^T, \quad (6) $$

where $\operatorname{dir}(\nu) = v/||v||_2 \mathbf{1}\{v \ne 0\}$ for all $v \in \mathbb{R}^p$. The previous decomposition depends on the quantities $\pi_{\nu}^{\perp} \mathbf{X}$ and $\operatorname{dir}(\mathbf{X}^T \nu)$, whose null distributions remain unknown. As a consequence, the authors in [21] condition on their values for a realization $\mathbf{x}$ of $\mathbf{X}$, defining the following quantity:

$$ \begin{gathered} p(\mathbf{x}; \{G_1, G_2\}) = \mathbb{P}_{H_0^{\{G_1, G_2\}}} \left( \begin{array}{l} ||\mathbf{X}^T \nu||_2 \ge ||\mathbf{x}^T \nu||_2 \\ \pi_{\nu}^{\perp} \mathbf{X} = \pi_{\nu}^{\perp} \mathbf{x}, \operatorname{dir}(\mathbf{X}^T \nu) = \operatorname{dir}(\mathbf{x}^T \nu) \end{array} \mid G_1, G_2 \in \mathcal{C}(\mathbf{X}) \right), \quad (\text{p-GBW}) \\ \pi_{\nu}^{\perp} \mathbf{X} = \pi_{\nu}^{\perp} \mathbf{x}, \operatorname{dir}(\mathbf{X}^T \nu) = \operatorname{dir}(\mathbf{x}^T \nu) \end{gathered} \end{equation} $$

as a p-value for (H0) [21, Theorem 1].

The key challenge in proposing (p-GBW) is finding an efficient characterization suitable for practical application. The idea in [21] involves two steps. The first is the definition of a test statistic based on the norm induced by the null covariance matrix of $\mathbf{X}^T \nu$ (up to a positive multiplicative factor). More precisely, if $\mathbf{A}$ is the covariance matrix of a non-degenerated, centered p-dimensional Gaussian vector $y$, then $||y||_2^2 = y^T \mathbf{A}^{-1} y$ follows a $\chi_p^2$ distribution. This implies that $||\mathbf{X}^T \nu||_2$ follows a $\sigma||\nu||_2 \cdot \chi_p$ distribution under (H0), thereby justifying the choice of the $\ell^2$-norm. The second step is to show that $||\mathbf{X}^T \nu||_2$ is independent of both the direction and the projection in (p-GBW). Consequently, the p-value (p-GBW) can be expressed in terms of a $\chi_p$ distribution truncated to a set $\hat{\mathcal{S}}$ that accounts for the event $\{C_1, C_2 \in \mathcal{C}(\mathbf{X})\}$. If $\mathcal{C}$ is a hierarchical clustering algorithm, the set $\hat{\mathcal{S}}$-and thus (p-GBW)- can be efficiently computed for several types of linkages. Otherwise, it can be approximated with a Monte Carlo procedure.

## 2.2 Extension to the general matrix normal model

### 2.2.1 Feasibility of a straightforward extension of Gao et al.

Here, we aim at extending (p-GBW) for the general model (gen-MN), following the same strategy to ensure the tractability of the p-value. Noticing that, under (H0), $\mathbf{X}^T \nu \sim N_p(0, \mathbf{V}_{G_1, G_2})$, where

$$ \mathbf{V}_{G_1, G_2} = \nu^T \mathbf{U} \nu \mathbf{\Sigma}, \quad (7) $$
a natural extension corresponds to replace $|| \cdot ||_2$ by the more general norm

$$ ||v||_{V_{G_1,G_2}} = \sqrt{v^T V_{G_1,G_2}^{-1} v}, \quad \forall v \in \mathbb{R}^p, \qquad (8) $$

which satisfies $||\mathbf{X}^T \nu||_{V_{G_1,G_2}} \sim \chi_p$ under the null. This choice leads us to consider the quantity

$$ p_{V_{G_1,G_2}}(\mathbf{x}; \{G_1, G_2\}) = \mathbb{P}_{H_0^{\{G_1, G_2\}}} \left( \begin{array}{l} ||\mathbf{X}^T \nu||_{V_{G_1, G_2}} \ge ||\mathbf{x}^T \nu||_{V_{G_1, G_2}} \\ \mathbf{G}_1, \mathbf{G}_2 \in \mathcal{C}(\mathbf{X}), \\ \pi_{\nu}^{\perp} \mathbf{X} = \pi_{\nu}^{\perp} \mathbf{x}, \mathrm{dir}_{V_{G_1, G_2}}(\mathbf{X}^T \nu) = \mathrm{dir}_{V_{G_1, G_2}}(\mathbf{x}^T \nu) \end{array} \right), \quad (\text{p-gen}) $$

as a candidate p-value to extend Theorem 1 in [21], where $\mathrm{dir}_{V_{G_1, G_2}}(v) = v/||v||_{V_{G_1, G_2}} 1\{v \ne 0\}$ for all $v \in \mathbb{R}^p$.

A straightforward generalization of [21] to (gen-MN) can be obtained when the test statistic is independent of the projection and direction in (p-gen), using the same argument as in the second step of the approach of [21]. However, the following result shows that the independence on the extra conditioning events holds if and only if the norm (8) is chosen to define the test statistic and $\mathbf{U}$ belongs to the class of positive definite compound symmetry matrices:

$$ \mathcal{CS}(n) = \left\{ (a-b)\mathbf{I}_n + b\mathbf{1}_{n \times n} : a \ge 0, -\frac{a}{n-1} < b < a \right\}, \qquad (9) $$

where $\mathbf{1}_{n \times n}$ is a $n \times n$ matrix of ones. Note that $\mathcal{CS}(n)$ is the set of covariance matrices of the vectors $(y_1+\epsilon, \dots, y_n+\epsilon)$, where the $y_i$ are centered i.i.d. Gaussian variables and $\epsilon$ is a centered noise independent of the $y_i$.

**Proposition 2.1.** Let $\mathcal{C}$ be a clustering algorithm and $\mathbf{x}$ a realization of $\mathbf{X} \sim \mathcal{MN}_{n \times p}(\mu, \mathbf{U}, \mathbf{\Sigma})$. For any $p \times p$ symmetric positive definite matrix $\mathbf{A}$, let $\mathrm{dir}_{\mathbf{A}}(v) = v/||v||_{\mathbf{A}} 1\{v \ne 0\}$ for all $v \in \mathbb{R}^p$. Then,

$$ (i) \quad \mathbf{U} \in \mathcal{CS}(n) \Leftrightarrow \mathbf{X}^T \nu_{G_1, G_2} \perp \pi_{\nu_{G_1, G_2}}^{\perp} \mathbf{X} \quad \text{for all } (G_1, G_2) \in \mathcal{C}_{[n]}, $$

$$ (ii) \quad \text{For any } (G_1, G_2) \in \mathcal{C}_{[n]}, \quad \mathbf{A} = c\mathbf{V}_{G_1, G_2} \text{ for some } c > 0 \quad \begin{array}{l} H_0^{\{G_1, G_2\}} \\ \Leftrightarrow \\ ||\mathbf{X}^T \nu_{G_1, G_2}||_{\mathbf{A}} \perp \mathrm{dir}_{\mathbf{A}}(\mathbf{X}^T \nu_{G_1, G_2}). \end{array} $$

The previous result is proved in Appendix A.1.1. The first equivalence is established by showing that both conditions are simultaneously equivalent to $\nu_{G_1, G_2}$ being an eigenvector of $\mathbf{U}$. The second follows from the fact that a Gaussian vector is independent of its direction if and only if it is centered with spherical covariance, a condition that we rewrite in terms of the matrix $\mathbf{A}$.

Proposition 2.1 shows that the choice of the norm (8) to define the test statistic not only ensures the tractability of its null distribution but also its independence with respect to the direction in (p-gen). Furthermore, the independence $\mathbf{X}^T \nu \perp \pi_{\nu}^{\perp} \mathbf{X}$ is equivalent to $\mathbf{U} \in \mathcal{CS}(n)$. In other words, the direct extension of the strategy in [21] to the general model (gen-MN) imposes a compound symmetry constraint on the dependence between observations. We develop the framework $\mathbf{U} \in \mathcal{CS}(n)$ in Section 2.2.2. In Section 2.2.3, we explore the extension of the same strategy to arbitrary $\mathbf{U}$, focusing on the characterization of quantities of the form (p-GBW) when the extra conditioning events are not independent of the test statistic.
### 2.2.2 Compound symmetry dependence between observations

If the dependence between observations **U** has a compound symmetry structure, the quantity (p-gen) can be efficiently written in terms of a truncated $\chi_p$ distribution, and used as a p-value for (H0). This is stated in the next result, that extends Theorem 1 in [21].

**Theorem 2.2.** Let $\mathcal{C}$ be a clustering algorithm and $\mathbf{x}$ a realization of $\mathbf{X} \sim \mathcal{MN}_{n \times p}(\mu, \mathbf{U}, \mathbf{\Sigma})$ with $\mathbf{U} \in \mathcal{CS}(n)$. Then,

$$p_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\}) = 1 - \mathbb{F}_p (\|\mathbf{x}^T \nu\|_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}, S_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\})), \quad (\text{p-tract})$$

where $\mathbb{F}_p(t, \mathcal{S})$ is the cumulative distribution function of a $\chi_p$ random variable truncated to the set $\mathcal{S}$ and

$$S_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\}) = \left\{ \phi \ge 0 : \mathcal{G}_1, \mathcal{G}_2 \in \mathcal{C} \left( \pi_{\nu}^{\perp} \mathbf{x} + \phi \frac{\nu}{\|\nu\|_2^2} \mathrm{dir}_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}^T \nu) \right) \right\}, \quad (10)$$

for any $(\mathcal{G}_1, \mathcal{G}_2) \in \mathcal{C}_{[n]}$. Furthermore, $p_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\})$ is a p-value for (H0) that controls the selective type I error for clustering (5) exactly at level $\alpha$.

The proof of the previous result is presented in Appendix A.1.2. One can easily verify that setting $\mathbf{U} = \mathbf{I}_n$ and $\mathbf{\Sigma} = \sigma^2 \mathbf{I}_p$ in Theorem 2.2 yields exactly Theorem 1 in [21]. In this general version, the information about the variance has been extracted from the statistic null distribution, which now remains the same independently of $\mathbf{U}, \mathbf{\Sigma}$, and moved it into the test statistic itself by making it dependent on the scale matrices. More precisely, $\|\mathbf{X}^T \nu\|_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}$ is the Mahalanobis distance [34] between the empirical group means with respect to the null distribution of their difference. This distance generalizes to multiple dimensions the idea of quantifying how many standard deviations away a point is from the mean of its distribution, and therefore integrates the dependence structure between columns and rows in $\mathbf{X}$.

Following (p-tract), the computation of (p-gen) only depends on the characterization of the one-dimensional set

$$\hat{S}_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}} = S_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\}) = \left\{ \phi \ge 0 : \mathcal{G}_1, \mathcal{G}_2 \in \mathcal{C}(\mathbf{x}'_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\phi)) \right\}, \quad (11)$$

where

$$\mathbf{x}'_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\phi) = \pi_{\nu}^{\perp} \mathbf{x} + \phi \frac{\nu}{\|\nu\|_2^2} \mathrm{dir}_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}^T \nu). \quad (12)$$

The data set (12) is analogous to $\mathbf{x}'(\phi)$ in [21, Equation (13)] for the norm (8), and its interpretation is equivalent. Indeed, we can rewrite both $\mathbf{x}'(\phi)$ and (12) as

$$\mathbf{x}'(\phi) = \mathbf{x} + \frac{\nu}{\|\nu\|_2^2} (\phi - \|\mathbf{x}^T \nu\|_2) \mathrm{dir}(\mathbf{x}^T \nu), \quad (13)$$

$$\mathbf{x}'_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\phi) = \mathbf{x} + \frac{\nu}{\|\nu\|_2^2} (\phi - \|\mathbf{x}^T \nu\|_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}) \mathrm{dir}_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}^T \nu). \quad (14)$$

Consequently, we can interpret (12) as a perturbed version of $\mathbf{x}$, but where the perturbation is based on the norm $\| \cdot \|_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}$ instead of $\| \cdot \|_2$. Thus, (11) is the set of non-negative $\phi$ for which applying the clustering algorithm $\mathcal{C}$ to the perturbed data set $\mathbf{x}'_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\phi)$ yields $\mathcal{G}_1$ and $\mathcal{G}_2$. As shown in [21], the set

$$\hat{S} = \{\phi \ge 0 : \mathcal{G}_1, \mathcal{G}_2 \in \mathcal{C}(\mathbf{x}'(\phi))\}, \quad (15)$$
can be explicitly characterized for hierarchical agglomerative clustering with several types of linkages. The next Lemma shows that we do not need to re-adapt the work in [21] to the set (11), as its points are given by a scale transformation of the points in $\hat{S}$.

**Lemma 2.3.** Let $\mathbf{X}$ be a realization of $\mathbf{X}$ and $\mathcal{G}_1, \mathcal{G}_2$ an arbitrary pair of clusters in $\mathcal{C}(\mathbf{x})$. Let $\hat{S}$ denote the set (15) defined in [21, Equation (12)]. Then,

$$ \hat{S}_{\mathcal{V}_{\mathcal{G}_1, \mathcal{G}_2}} = \frac{\|\mathbf{x}^T \nu\|_{\mathcal{V}_{\mathcal{G}_1, \mathcal{G}_2}}}{\|\mathbf{x}^T \nu\|_2} \hat{S}, \qquad (16) $$

where $\hat{S}_{\mathcal{V}_{\mathcal{G}_1, \mathcal{G}_2}}$ is defined in (11).

Consequently, the work in [21, Section 3] can be applied here to characterize the set (11) and, therefore, to compute the p-value defined in (p-gen). An explicit characterization of (11) is possible when $\mathcal{C}$ is a hierarchical clustering algorithm with squared Euclidean distance, along with either single linkage or a linkage satisfying a linear Lance-Williams update [21, Equation 20], e.g. average, weighted, Ward, centroid or median linkage. The efficient computation of (11) can also be extended to k-means clustering using the work in [9], as shown in Appendix B. Otherwise, the p-value (p-gen) can be approximated with a Monte Carlo procedure, adapting the importance sampling approach presented in [21, Section 4.1]. Following the same notation, we sample $\omega_1, \dots, \omega_N \sim \mathcal{N}(\|\mathbf{x}^T \nu\|_{\mathcal{V}_{\mathcal{G}_1, \mathcal{G}_2}}, 1)$ i.i.d. and approximate (p-gen) as

$$ p_{\mathcal{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\}) \approx \frac{\sum_{i=1}^{N} \pi_i \mathbf{1}\{\omega_i \ge \|\mathbf{x}^T \nu\|_{\mathcal{V}_{\mathcal{G}_1, \mathcal{G}_2}}, \mathcal{G}_1, \mathcal{G}_2 \in \mathcal{C}(\mathbf{x}'(\omega_i))\}}{\sum_{i=1}^{N} \pi_i \mathbf{1}\{\mathcal{G}_1, \mathcal{G}_2 \in \mathcal{C}(\mathbf{x}'(\omega_i))\}}, \qquad (17) $$

for $\pi_i = f_1(\omega_i)/f_2(\omega_i)$, where $f_1$ is the density of a $\chi_p$ random variable, and $f_2$ is the density of a $\mathcal{N}(\|\mathbf{x}^T \nu\|_{\mathcal{V}_{\mathcal{G}_1, \mathcal{G}_2}}, 1)$ random variable.

### 2.2.3 General dependence between observations

The goal of this section is to study whether a quantity of the form (p-GBW) can be (i) efficiently characterized in terms of a known distribution and/or (ii) used as a p-value for (H0), in the case where the restriction $\mathbf{U} \in \mathcal{C}\mathcal{S}(n)$ is not necessarily satisfied. Following from Proposition 2.1(i), this means that the projection $\pi_\nu^\perp \mathbf{X}$ is not independent of $\mathbf{X}^T \nu$ in general and, therefore, that the distribution of interest for the definition of the test statistic is not that of $\mathbf{X}^T \nu$, but rather that of the conditioned vector:

$$ \bar{\mathbf{X}}_\nu(\mathbf{x}) := \mathbf{X}^T \nu \mid \{\pi_\nu^\perp \mathbf{X} = \pi_\nu^\perp \mathbf{x}, \mathrm{dir}(\mathbf{X}^T \nu) = \pm \mathrm{dir}(\mathbf{x}^T \nu)\}, \quad \text{for } \mathbf{x} \in \mathbb{R}^{n \times p}. \qquad (18) $$

Adding the $\pm$ symbol allows to express the conditioning set as a linear constraint, which is more suitable for Gaussian processes. Then, the distribution of $\mathbf{X}^T \nu$ conditioned on the original conditioning set can be recovered by truncating the density function of $\bar{\mathbf{X}}_\nu(\mathbf{x})$ to the half space $\{\mathbf{y} \in \mathbb{R}^p : \langle \mathbf{y}, \mathbf{x}^T \nu \rangle \ge 0\}$. The null distribution of (18) is derived in the next result. In what follows, we will denote by $\mathbf{A}^\dagger$ the Moore-Penrose pseudo-inverse of a matrix $\mathbf{A}$.

**Theorem 2.4.** Let $\mathcal{C}$ be a clustering algorithm and $\mathbf{x}$ a realization of $\mathbf{X} \sim \mathcal{MN}_{n \times p}(\mu, \mathbf{U}, \mathbf{\Sigma})$. Then, under (H0),

$$ \bar{\mathbf{X}}_\nu(\mathbf{x}) \sim \mathcal{N}_p(0, \mathbf{\Gamma}_\nu), \qquad (19) $$

with

$$ \mathbf{\Gamma}_\nu = (\mathbf{I}_p \otimes \nu^T)(\mathbf{\Sigma} \otimes \mathbf{U} - (\mathbf{\Sigma} \otimes \mathbf{U})\mathbf{A}_\nu^T(\mathbf{A}_\nu(\mathbf{\Sigma} \otimes \mathbf{U})\mathbf{A}_\nu^T)^\dagger\mathbf{A}_\nu(\mathbf{\Sigma} \otimes \mathbf{U}))(\mathbf{I}_p \otimes \nu), \qquad (20) $$
where $\mathbf{A}_x$ is a $2np \times p$ matrix given by:

$$ \mathbf{A}_x = \begin{bmatrix} \pi_{x\nu}^{\perp} (\mathbf{I}_p \otimes \pi_{\nu}) \\ \mathbf{I}_p \otimes \pi_{\nu}^{\perp} \end{bmatrix}, \qquad (21) $$

with $\pi_{\nu} = \mathbf{I}_n - \pi_{\nu}^{\perp}$, $\mathbf{x}_{\nu} = \text{vec}(\pi_{\nu}\mathbf{x})$ and $\pi_{x\nu}^{\perp} = \mathbf{I}_{np} - \mathbf{x}_{\nu}^{T}\mathbf{x}_{\nu}/\|\mathbf{x}_{\nu}\|_2^2$.

The proof of Theorem 2.4, presented in Appendix A.1.3, proceeds in two main steps. First, we express the conditioning set in (18) as a linear constraint, which allows us to apply Proposition 3.13 in [19], characterizing the distribution of Gaussian vectors $z$ conditioned to events of the form $\{\mathbf{A}z = y\}$. The structure of the conditioned covariance and mean matrices motivates a detailed analysis of some specific matrix families (Lemma A.4 and Corollary A.5). The corresponding results allow us to prove that $\bar{\mathbf{X}}_{\nu}(\mathbf{x})$ is centered under (H0) for all $\mathbf{x} \in \mathbb{R}^{n \times p}$.

Following from (19), in order to define a quantity of the form (p-GBW), the same reasoning as in the previous sections leads us to consider the following norm based on the covariance matrix (20):

$$ \|v\|_{\Gamma_x} = \sqrt{v^T \Gamma_x^\dagger v}, \quad \forall v \in \mathbb{R}^p, \qquad (22) $$

for any $\mathbf{x} \in \mathbb{R}^{n \times p}$, where we have considered the generalized inverse of (20) as this matrix is not full-rank. This leads us to define the quantity:

$$ \begin{aligned} p_{\Gamma}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\}) = \mathbb{P}_{H_0^{\{\mathcal{G}_1, \mathcal{G}_2\}}} & \left( \|\mathbf{X}^T \nu\|_{\Gamma_x} \ge \|\mathbf{x}^T \nu\|_{\Gamma_x} \mid \mathcal{G}_1, \mathcal{G}_2 \in \mathcal{C}(\mathbf{X}), \right. \\ & \left. \pi_{\nu}^{\perp} \mathbf{X} = \pi_{\nu}^{\perp} \mathbf{x}, \text{ dir}(\mathbf{X}^T \nu) = \pm \text{dir}(\mathbf{x}^T \nu) \right), \end{aligned} \qquad (\text{p-Gamma}) $$

as a candidate $p$-value for (H0) under (gen-MN). The previous quantity has an unusual form for a $p$-value, since the test statistic $\|\mathbf{X}^T \nu\|_{\Gamma_x}$ depends on the realization $\mathbf{x}$. However, the following result shows that its distribution under (H0) is almost surely independent of $\mathbf{x}$, yielding an efficient characterization of (p-Gamma). Its proof is presented in Appendix A.1.3.

**Proposition 2.5.** In the conditions of Theorem 2.4, the quantity $\|\bar{\mathbf{X}}_{\nu}(\mathbf{x})\|_{\Gamma_x}$ follows $\mathbf{x}$-a.s. a $\chi_1$ distribution under (H0). Moreover, the quantity (p-Gamma) can be written as:

$$ p_{\Gamma}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\}) = 1 - \mathbb{F}_1 \left( \|\mathbf{x}^T \nu\|_{\Gamma_x}, S_{\Gamma_x}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\}) \right), \qquad (23) $$

for any $(\mathcal{G}_1, \mathcal{G}_2) \in \mathcal{C}_{[n]}$, where $\mathbb{F}_1(t, S)$ is the cumulative distribution function of a $\chi_1$ random variable truncated to the set $S$ and

$$ S_{\Gamma_x}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\}) = \frac{\|\mathbf{x}^T \nu\|_{\Gamma_x}}{\|\mathbf{x}^T \nu\|_2} \{\phi \ge 0 : \mathcal{G}_1, \mathcal{G}_2 \in \mathcal{C}(\mathbf{x}'(\pm\phi))\}, \qquad (24) $$

where $\mathbf{x}'(\phi)$ is the perturbed data set defined in (14).

The previous result allows the efficient computation of (p-Gamma) in terms of a $\chi_1$ distribution. Equation (23) is the counterpart of [21, Equation (9)] and (p-gen) for the most general case, where (gen-MN) holds with arbitrary $\mathbf{U}$. However, although the quantity (p-Gamma) is the natural extension of (p-GBW) in this context, and can be efficiently characterized via (23), assessing whether it controls the selective type I error is a challenging problem. More precisely, the null distribution of the conditioned vector (18) depends on the realization $\mathbf{x}$ and, consequently, the norm required to ensure that the test statistic is distribution-free is also dependent on $\mathbf{x}$. As a consequence, (p-Gamma) compares two quantities that
behave differently under (H0). Indeed, in order to assess whether (5) is satisfied, it is necessary to understand the behavior of the null distribution of $||\bar{\mathbf{X}}_\nu(\mathbf{x})||^2_{\Gamma_\mathbf{x}} = \bar{\mathbf{X}}_\nu(\mathbf{x})^T \Gamma_\mathbf{x}^\dagger \bar{\mathbf{X}}_\nu(\mathbf{x})$, which is a nontrivial problem. Nevertheless, since Proposition 2.5 allows the computation of (p-Gamma) in practice, we are able to illustrate numerically that the quantity (23) does not control the selective type I error for several $\mathbf{U} \notin \mathcal{CS}(n)$ structures. We present these simulations in Appendix D.3.

As we further discuss in Section 6, the analyses presented above suggest that defining a tractable p-value of the form (p-GBW) that ensures the selective type I error control requires the conditioning on events that are *independent* of the test statistic. Following from Proposition 2.1, this is ensured if and only if the covariance structure between observations has a compound symmetry structure and the norm (8) is used to define the p-value.

## 3 Unknown dependence structures

The selective inference framework introduced for (gen-MN) in Section 2 assumes that both scale matrices $\mathbf{U}$ and $\mathbf{\Sigma}$ are known, which is a quite unrealistic scenario. Under the independence assumption made in [21], where $\mathbf{\Sigma} = \sigma^2 \mathbf{I}_p$ and $\mathbf{U} = \mathbf{I}_n$, the authors showed in Theorem 4 that over-estimating $\sigma$ yields asymptotic control of the selective type I error, and provided such an estimator $\hat{\sigma}$ that can be used in practice.

The simultaneous estimation of $\mathbf{U}$ and $\mathbf{\Sigma}$ from a single copy of $\mathbf{X}$ is a challenging task due to the intrinsic limitations of the matrix normal model. The non-identifiability of both matrices under (gen-MN) makes their existing estimators interdependent. Besides, multiple realizations of $\mathbf{X}$ are needed to ensure their existence and uniqueness [17]. The same goes for the estimation of $\mathbf{U} \otimes \mathbf{\Sigma}$, that fully determines the covariance structure of $\mathbf{X}$ [14–16, 44], even when $\mathbf{U}$ is restricted to the class $\mathcal{CS}(n)$ [1]. All of this hinders the estimation of both covariance matrices from a single copy of $\mathbf{X}$. Furthermore, we not only require any estimator of $\mathbf{U}$ and $\mathbf{\Sigma}$ but one which is compatible with the selective type I error control. Consequently, we opt to investigate the situation where only one of the scale matrices is known, and assess theoretical conditions that allow asymptotic control of the selective type I error when estimating the other one. We also provide an estimator that satisfies these conditions for some common dependence models.

Let us recall that, for the model (gen-MN), we have

$$ \mathbf{X} \sim \mathcal{MN}_{n \times p}(\mu, \mathbf{U}, \mathbf{\Sigma}) \Leftrightarrow \mathbf{X}^T \sim \mathcal{MN}_{p \times n}(\mu^T, \mathbf{\Sigma}, \mathbf{U}). \quad (25) $$

Therefore, the methods presented in this section can be equally applied to estimate $\mathbf{U}$ or $\mathbf{\Sigma}$ when the other is known, by transposing $\mathbf{X}$ if needed. From now on, we assume that the dependence structure between observations $\mathbf{U}$ is known, and study under which conditions we can suitably estimate $\mathbf{\Sigma}$. In Section 3.1, we focus on the case where a computationally tractable p-value can be defined according to Theorem 2.2, assessing the applicability of (p-tract) when $\mathbf{\Sigma}$ is estimated with $\mathbf{U} \in \mathcal{CS}(n)$. Since the robustness of (p-tract) to $\mathbf{U} \notin \mathcal{CS}(n)$ will be numerically studied, in Section 3.2 we explore the theoretical guarantees that can be provided in that case regarding the estimation of $\mathbf{\Sigma}$.

### 3.1 Compound symmetry covariance between observations

Let $\hat{\mathbf{\Sigma}}(\mathbf{x})$ be an estimate of $\mathbf{\Sigma}$ for a given realization $\mathbf{x}$ of $\mathbf{X}$. Following from Theorem 2.2, the p-value (p-gen) has the closed form (p-tract) if $\mathbf{U} \in \mathcal{CS}(n)$. In that case, the estimation of $\mathbf{\Sigma}$ comes down
to studying under which conditions the p-value

$$p_{\hat{V}_{G_1, G_2}}(\mathbf{x}; \{G_1, G_2\}) = 1 - \mathbb{F}_p (\|\mathbf{x}^T \nu\|_{\hat{V}_{G_1, G_2}}; S_{\hat{V}_{G_1, G_2}}(\mathbf{x}; \{G_1, G_2\})), \quad (\text{hat-p-tract})$$

where $\hat{V}_{G_1, G_2} = \nu^T U \nu \hat{\Sigma}(\mathbf{x})$, controls the selective type I error. Theorem 3.1 below generalizes Theorem 4 in [21] for the estimation of $\Sigma$ under the model (gen-MN) by relying on the Loewner partial order, defined below. The proof is included in Appendix A.2.

**Definition 3.1** (Definition 7.7.1 in [26]). For two square matrices of equal size A, B, we write $A \succeq B$ if and only if A, B are Hermitian and $A - B$ is positive semidefinite. This binary relation between square matrices is called the Loewner partial order.

**Theorem 3.1.** For $n \in \mathbb{N}$, let $\mathbf{X}^{(n)} \sim \mathcal{MN}_{n \times p}(\mu^{(n)}, U^{(n)}, \Sigma)$ with $U^{(n)} = (a-b)\mathbf{I}_n + b\mathbf{1}_{n \times n}$ for some $a > b > 0$. Let $\mathbf{x}^{(n)}$ be a realization of $\mathbf{X}^{(n)}$ and $G_1^{(n)}, G_2^{(n)}$ a pair of clusters estimated from $\mathbf{x}^{(n)}$. If $\hat{\Sigma}(\mathbf{X}^{(n)})$ is a positive definite estimator of $\Sigma$ such that

$$\lim_{n \to \infty} \mathbb{P}_{H_0}^{\{\sigma_1^{(n)}, \sigma_2^{(n)}\}} \left( \hat{\Sigma}(\mathbf{X}^{(n)}) \succeq \Sigma \mid G_1^{(n)}, G_2^{(n)} \in \mathcal{C}(\mathbf{X}^{(n)}) \right) = 1, \quad (\text{over-est})$$

then,

$$\limsup_{n \to \infty} \mathbb{P}_{H_0}^{\{\sigma_1^{(n)}, \sigma_2^{(n)}\}} \left( p_{\hat{V}_{G_1^{(n)}, G_2^{(n)}}}(\mathbf{X}^{(n)}; \{G_1^{(n)}, G_2^{(n)}\}) \le \alpha \mid G_1^{(n)}, G_2^{(n)} \in \mathcal{C}(\mathbf{X}^{(n)}) \right) \le \alpha, \quad (26)$$

for any $\alpha \in [0, 1]$.

Note that the Loewner partial order is a natural extension to Hermitian matrices of the usual order in $\mathbb{R}$. If we replace $\Sigma$ by $\sigma^2 \mathbf{I}_p$ in Theorem 3.1, the condition $\hat{\Sigma} \succeq \Sigma$ becomes $\hat{\sigma} \ge \sigma$, as in [21, Theorem 4]. We aim now at providing an estimator of $\Sigma$ satisfying condition (over-est). The asymptotic properties of such an estimator strongly depend on the asymptotic dependence structure between observations, given by the sequence of matrices $\{U^{(n)}\}_{n \in \mathbb{N}}$ of Theorem 3.1. First, let us consider

$$\hat{\Sigma} = \hat{\Sigma}(\mathbf{X}) = \frac{1}{n-1} (\mathbf{X} - \bar{\mathbf{X}})^T U^{-1} (\mathbf{X} - \bar{\mathbf{X}}), \quad (\text{hat-Sigma})$$

where $\bar{\mathbf{X}}$ is a $n \times p$ matrix having as rows the mean across rows of $\mathbf{X}$, i.e.

$$\bar{\mathbf{X}} = \mathbf{1}_n \otimes \frac{1}{n} \sum_{k=1}^{n} X_k, \quad (27)$$

where $\mathbf{1}_n$ is a column $n$-vector of ones. Note that (hat-Sigma) is constructed by first de-correlating the observations using $U$, then subtracting off the column means and finally taking the sample covariance matrix. Following [23, Corollary 2.3.10.2], subtracting off the true mean matrix $\mu$ instead of $\bar{\mathbf{X}}$ would lead to a consistent estimator without making any assumption on $U$, as the rows of $U^{-\frac{1}{2}}(\mathbf{X} - \mu)$ are $n$ i.i.d. copies of a $p$-dimensional centered Gaussian vector of covariance matrix $\Sigma$. However, $\mu$ needs to be considered unknown in the context of clustering analysis. Note also that the estimator $\hat{\Sigma}$ is a positive definite matrix if the matrix $\mathbf{X} - \bar{\mathbf{X}}$ has full rank. In order to ensure that (hat-Sigma) satisfies condition (over-est), some additional assumptions regarding the asymptotic behavior of the matrices $\mu^{(n)}$ are required.
**Assumption 3.1** (Assumptions 1 and 2 in [21]). For all $n \in \mathbb{N}$, there are exactly $K^*$ distinct mean vectors among the first $n$ observations, i.e.

$$ \{\mu_i^{(n)}\}_{i=1,...,n} = \{\theta_1,...,\theta_{K^*}\}. \qquad (28) $$

Moreover, the proportion of the first $n$ observations that have mean vector $\theta_k$ converges to $\pi_k > 0$, i.e.

$$ \lim_{n \to \infty} \frac{1}{n} \sum_{i=1}^{n} 1\{\mu_i^{(n)} = \theta_k\} = \pi_k, \qquad (29) $$

for all $k \in \{1,..., K^*\}$, where $\sum_{k=1}^{K^*} \pi_k = 1$.

If observations are independent, Assumption 3.1 is the only requirement for (hat-Sigma) to asymptotically over-estimate $\Sigma$ in the sense of Theorem 3.1. For non-diagonal $U^{(n)}$, the following condition on $\{\mu^{(n)}\}_{n \in \mathbb{N}}$ needs to be assumed.

**Assumption 3.2.** If $U^{(n)}$ is non-diagonal for all $n \in \mathbb{N}$, for any $k, k' \in \{1,..., K^*\}$, the proportion of the first $n$ observations at distance $r \ge 1$ in $X^{(n)}$ having means $\theta_k$ and $\theta_{k'}$ converges, and its limit converges to $\pi_k \pi_{k'}$ when the lag $r$ tends to infinity. More precisely,

$$ \lim_{n \to \infty} \frac{1}{n} \sum_{i=1}^{n-r} 1\{\mu_i^{(n)} = \theta_k\} 1\{\mu_{i+r}^{(n)} = \theta_{k'}\} = \pi_{kk'}^r \underset{r \to \infty}{\longrightarrow} \pi_k \pi_{k'}. \qquad (30) $$

Note that we are requiring the proportion of pairs of observations having a given a pair of means to approach the product of individual proportions (29) when both observations are far away in $X^{(n)}$. Assumption 3.2 can be alternatively formulated in terms of strong mixing of measure-preserving dynamical systems [29, Chapter 20]. This is proved in Appendix A.2.

If $U^{(n)}$ is compound symmetry for fixed $a > b > 0$ and Assumptions 3.1 and 3.2 hold for a given sequence $\{\mu^{(n)}\}_{n \in \mathbb{N}}$, the following result ensures that $\hat{\Sigma}$ asymptotically over-estimates (in the sense of the Loewner partial order) the dependence structure $\Sigma$ between features.

**Proposition 3.2.** Let $X^{(n)} \sim \mathcal{MN}_{n \times p}(\mu^{(n)}, U^{(n)}, \Sigma)$, where $U^{(n)} = (a-b)I_n + b1_{n \times n}$ for some $a > b > 0$ and $\mu^{(n)}$ satisfies Assumptions 3.1 and 3.2 for some $K^* > 1$. Let $\hat{\Sigma}$ be the estimator defined in (hat-Sigma). Then,

$$ \lim_{n \to \infty} \mathbb{P} \left( \hat{\Sigma} \left( \mathbf{X}^{(n)} \right) \succeq \Sigma \right) = 1. \qquad (31) $$

Finally, it suffices to estimate $\Sigma$ using an independent and identically distributed copy of $X^{(n)}$ to have (over-est) provided (31) holds. Such a copy is sometimes available in practical applications, as the one we present in Section 5. Combining this observation with Proposition 3.2, we obtain our final result:

**Proposition 3.3.** Let $X^{(n)} \sim \mathcal{MN}_{n \times p}(\mu^{(n)}, U^{(n)}, \Sigma)$, where $U^{(n)} = (a-b)I_n + b1_{n \times n}$ for some $a > b > 0$ and $\mu^{(n)}$ satisfies Assumptions 3.1 and 3.2 for some $K^* > 1$. Let $\mathbf{x}^{(n)}$ be a realization of $X^{(n)}$ and $G_1^{(n)}$, $G_2^{(n)}$ a pair of clusters estimated from $\mathbf{x}^{(n)}$. Let $\mathbf{Y}^{(n)}$ be an independent and identically distributed copy of $X^{(n)}$. Then, the estimator $\hat{\Sigma}(\mathbf{Y}^{(n)})$ defined in (hat-Sigma) satisfies the conditions of Theorem 3.1, i.e.

$$ \lim_{n \to \infty} \mathbb{P}_{H_0}^{\{\sigma_1^{(n)}, \sigma_2^{(n)}\}} \left( \hat{\Sigma}(\mathbf{Y}^{(n)}) \succeq \Sigma \mid G_1^{(n)}, G_2^{(n)} \in \mathcal{C}(\mathbf{X}^{(n)}) \right) = 1. \qquad (32) $$

Following from the previous result and from Theorem 3.1, if Assumptions 3.1 and 3.2 hold and $U \in \mathcal{CS}(n)$, selective type I error is asymptotically controlled when using (hat-Sigma) to estimate $\Sigma$.
This constitutes an extension of the over-estimation framework presented in [21], which holds under model (ind-MN), to the more general (gen-MN) with compound symmetry dependence between observations.

## 3.2 Arbitrary covariance between observations

In Section 3.1, we proved that p-values (hat-p-tract) are asymptotically super-uniform under (H0) if, besides Assumptions 3.1 and 3.2, the following conditions hold:

(a) The p-value (p-tract) (for known $\Sigma$) is uniformly distributed under (H0) or, equivalently, $U \in CS(n)$,

(b) The estimator (hat-Sigma) satisfies (over-est).

However, as it will be numerically illustrated in Section 4.4, the null uniformity of (p-tract) is robust to $U$ structures that do not fit in $CS(n)$. Consequently, the null super-uniformity of (hat-p-tract) will be robust to $U \notin CS(n)$ as long as (b) is satisfied. In this section, we investigate the theoretical conditions that need to be imposed to an arbitrary sequence $\{U^{(n)}\}_{n \in \mathbb{N}}$ so that the estimator (hat-Sigma) satisfies (over-est). To that end, besides Assumptions 3.1 and 3.2, the quantities

$$ \frac{1}{n} \sum_{l,s=1}^{n} (U^{(n)})_{ls}^{-1} \mathbf{1}\{\mu_l^{(n)} = \theta_k\} \mathbf{1}\{\mu_s^{(n)} = \theta_{k'}\} \quad (33) $$

are also required to converge. Furthermore, we need to know their limit explicitly to assess (over-est). Below, we state sufficient conditions on the sequence $\{U^{(n)}\}_{n \in \mathbb{N}}$ that -together with Assumptions 3.1 and 3.2- ensure the convergence of (33) to a tractable limit. Note that these technical assumptions can be difficult to verify for a given model of dependence, and other unknown sufficient conditions might guarantee that (hat-Sigma) asymptotically over-estimates $\Sigma$. This point is investigated numerically in Section 4.4.2.

**Assumption 3.3.** Let $\{U^{(n)}\}_{n \in \mathbb{N}}$ be a sequence of real positive definite matrices, and let $(U^{(n)})_{ij}^{-1}$ denote the $i, j$ entry of $(U^{(n)})^{-1}$ for any $n \in \mathbb{N}$. Then, every superdiagonal of $(U^{(n)})^{-1}$ defines asymptotically a convergent sequence, whose limits sum up to a real value. More precisely, for any $i \in \mathbb{N}$ and any $r \ge 0$,

$$ \lim_{n \to \infty} (U^{(n)})_{i i+r}^{-1} = \Lambda_{i i+r}, \quad \text{where} \quad \lim_{i \to \infty} \Lambda_{i i+r} = \lambda_r \quad \text{and} \quad \sum_{r=0}^{\infty} \lambda_r = \lambda \in \mathbb{R}. \quad (34) $$

Moreover, for each $r \ge 0$ any of the following conditions are satisfied:

(i) It exists a sequence $\{\alpha_i\}_{i=1}^\infty \in \ell_1$ such that $|(U^{(n)})_{i i+r}^{-1} - \Lambda_{i i+r}| \le \alpha_i$ for all $n \in \mathbb{N}$,

(ii) For each $i \in \mathbb{N}$, the sequence $\{(U^{(n)})_{i i+r}^{-1}\}_{n \in \mathbb{N}}$ is non-decreasing or non-increasing.

Note that Assumptions 3.2 and 3.3 implicitly require an ordering of the observations in $X$. More precisely, they require the existence of a permutation of the rows in $X$ such that their conditions are satisfied. The following result generalizes Proposition 3.2 to arbitrary sequences $\{U^{(n)}\}_{n \in \mathbb{N}}$.

**Proposition 3.4.** Let $X^{(n)} \sim MN_{n \times p}(\mu^{(n)}, U^{(n)}, \Sigma)$, where $\mu^{(n)}$ and $U^{(n)}$ satisfy Assumptions 3.1, 3.2 and 3.3 for some $K^* > 1$. Let $\hat{\Sigma}$ be the estimator defined in (hat-Sigma). Then,

$$ \lim_{n \to \infty} \mathbb{P} \left( \hat{\Sigma} \left( \mathbf{X}^{(n)} \right) \succeq \Sigma \right) = 1. \quad (35) $$
As a consequence, Proposition 3.3 directly holds for arbitrary {$U^{(n)}$}$_{n \in \mathbb{N}}$ if Assumption 3.3 is added to its hypotheses. Our proof of Proposition 3.4 relies on the following Lemma, which makes use of Assumptions 3.1, 3.2 and 3.3 explicitly. Both results are proved in Appendix A.2.

**Lemma 3.5.** Let $\mathbf{X}^{(n)} \sim \mathcal{MN}_{n \times p}(\mu^{(n)}, \mathbf{U}^{(n)}, \mathbf{\Sigma})$, where $\mu^{(n)}$ and $\mathbf{U}^{(n)}$ satisfy Assumptions 3.1, 3.2 and 3.3 for some $K^* > 1$. Then,

$$ \lim_{n \to \infty} \frac{1}{n} \sum_{l,s=1}^{n} (U^{(n)})_{ls}^{-1} \mathbf{1}\{\mu_l^{(n)} = \theta_k\} \mathbf{1}\{\mu_s^{(n)} = \theta_{k'}\} = 2(\lambda - \lambda_0)\pi_k \pi_{k'} + \lambda_0 \pi_k \delta_{kk'}, \quad (36) $$

for any $k, k' \in \{1, \dots, K'\}$, and where $\pi_k, \pi_{k'}$ and $\lambda_0, \lambda$ are defined in Assumptions 3.1 and 3.3 respectively.

Assessing whether a model of dependence satisfies Assumption 3.3 is not trivial as it requires full knowledge of how the inverse matrices $(U^{(n)})^{-1}$ grow up when the sample size increases. However, we are able to show that Assumption 3.3 is satisfied for some specific dependence models and, consequently, that selective type I error can be controlled when $\mathbf{\Sigma}$ is over-estimated in such cases. The following remarks are proved in Appendix A.2.

**Remark 3.1** (Compound symmetry). Let $\mathbf{U}^{(n)} = (a-b)\mathbf{I}_n + b\mathbf{1}_{n \times n}$ for some $a > b > 0$. Then, {$\mathbf{U}^{(n)}$}$_{n \in \mathbb{N}}$ satisfies Assumption 3.3.

The compatibility of compound symmetry structures with the over-estimation of $\mathbf{\Sigma}$ can be explained within this more general framework: Remark 3.1 and Proposition 3.4 imply Proposition 3.2. Therefore, we do not provide a direct proof of the latter result. We can also consider the case of independent observations with different variances along features. Note that, if the matrix $\mathbf{X}$ is transposed, any general dependence structure between observations $\mathbf{U}$ can be estimated if independent features with known variances are provided, which is already an important generalization of [21].

**Remark 3.2** (Diagonal). Let $\mathbf{U}^{(n)} = \text{diag}(\lambda_1, \dots, \lambda_n)$. If the sequence {$\lambda_n$}$_{n \in \mathbb{N}}$ is convergent, the sequence {$\mathbf{U}^{(n)}$}$_{n \in \mathbb{N}}$ satisfies Assumption 3.3.

We can extend the complexity of $\mathbf{U}^{(n)}$ to auto-regressive covariance structures of any lag. This is mainly thanks to the fact that the inverses of such matrices are tractable and banded, i.e. their non-zero entries are confined to a centered diagonal band. Under model (gen-MN), assuming that $\mathbf{U}^{(n)}$ is the covariance matrix of an auto-regressive process of order $P$ means that

$$ \frac{1}{\sqrt{\Sigma_{jj}}} X_{ij}^{(n)} = \frac{1}{\sqrt{\Sigma_{jj}}} \sum_{s=1}^{P} \beta_s X_{i-s j}^{(n)} + \varepsilon_i, \quad \forall j \in \{1, \dots, p\}, \quad (37) $$

where {$\varepsilon_i$}$_{i=1, \dots, n}$ are i.i.d univariate centered normal variables and {$\beta_s$}$_{s=1, \dots, P} \subset \mathbb{R}$ are the model coefficients. Then, for any $j \in \{1, \dots, p\}$, the entries of $\mathbf{U}^{(n)}$ would be given by

$$ U_{ii'} = \text{Cov} \left( \frac{X_{ij}}{\sqrt{\Sigma_{jj}}}, \frac{X_{i'j}}{\sqrt{\Sigma_{jj}}} \right), \quad \forall i, i' \in [n], \forall j \in \{1, \dots, p\}. \quad (38) $$

If model (37) is assumed, the covariance matrix $\mathbf{U}^{(n)}$ and its inverse have a tractable structure. For example, for the simplest auto-regressive process where $P = 1$, and the $i$-th observation depends linearly only on the $(i-1)$-th one, the entries of $\mathbf{U}^{(n)}$ have the form $U_{ij}^{(n)} = \sigma^2 \rho^{|i-j|}$, for $\sigma > 0$. To ensure the
the positive definiteness of $\mathbf{U}^{(n)}$, we need $|\rho| < 1$ (see the form of eigenvalues in [45]). This is equivalent to ask the process to be stationary. Then, the inverse of $\mathbf{U}^{(n)}$ is a tridiagonal matrix of the form

$$ (\mathbf{U}^{(n)})^{-1} = \frac{1}{\sigma^2(1 - \rho^2)} \begin{pmatrix} 1 & -\rho & & & \\ -\rho & 1 + \rho^2 & -\rho & & \\ & -\rho & \ddots & \ddots & \\ & & \ddots & 1 + \rho^2 & -\rho \\ & & & -\rho & 1 \end{pmatrix}. \quad (39) $$

The super and sub-diagonals trivially satisfy condition (i) in Assumption 3.3 with $\lambda_{\pm 1} = -\rho/(1 - \rho^2)$. Then, the entries of the main diagonal define the sequences

$$ \sigma^2(1 - \rho^2) \left\{ \left( U^{(n)} \right)_{ii}^{-1} \right\}_{n \in \mathbb{N}} = \begin{cases} \{1, 1, \dots\} & \text{if } i = 1, \\ \{\xi_1, \dots, \xi_{i-1}, 1, 1 + \rho^2, 1 + \rho^2, \dots\} & \text{if } i > 1, \end{cases} $$

for every $i \in \mathbb{N}$, where the entries $\sigma^2(1 - \rho^2) (U^{(n)})_{ii}^{-1} = \xi_n$ for $i > n$ can be chosen as needed. Note that these sequences do not satisfy condition (i) in Assumption 3.3, but they are non-decreasing (choosing appropriately the $\xi_k$). Consequently, Assumption 3.3 holds and we have $\Lambda_{11} = 1/(\sigma^2((1 - \rho^2)))$, $\Lambda_{ii} = \lambda_0 = (1 + \rho^2)/(\sigma^2((1 - \rho^2)))$ for all $i > 1$ and, finally, $\lambda = (1 - \rho^2)/(\sigma^2((1 - \rho^2)))$. For any $P \ge 1$, the inverse matrices are banded with $2P + 1$ non-zero diagonals and we can follow the same reasoning. However, for $P > 2$, we need to require the coefficients $\beta_1, \dots, \beta_P$ to have the same sign.

**Remark 3.3 (Auto-regressive).** Let $\mathbf{U}^{(n)}$ be the covariance matrix of an auto-regressive process of order $P \ge 1$ such that, if $P > 2$, $\beta_k\beta_{k'} \ge 0$ for all $k, k' \in \{1, \dots, P\}$. Then, the sequence $\{\mathbf{U}^{(n)}\}_{n \in \mathbb{N}}$ satisfies Assumption 3.3.

The above remarks imply that (hat-Sigma) satisfies (over-est) in the above-studied compound symmetry, diagonal and auto-regressive models. Consequently, the asymptotic null super-uniformity of (hat-p-tract) will be robust to $\mathbf{U}$ being diagonal or auto-regressive as long as the null uniformity of (p-tract) is robust to $\mathbf{U}$ belonging to such models (and Assumptions 3.1 and 3.2 hold).

# 4 Numerical experiments

In this section, we assess the numerical performance of the proposed approach in several scenarios simulated with synthetic data. We start by simulating settings that satisfy condition (ii) in Theorem 2.2, that is, choosing $\mathbf{U} \in \mathcal{CS}(n)$ and using the p-value (p-tract). The following three cases are considered for the scale matrices $\mathbf{U}$ and $\mathbf{\Sigma}$:

(D1) $\mathbf{U} = \mathbf{I}_n$ and $\mathbf{\Sigma}$ is the covariance matrix of an AR(1) model, i.e. $\Sigma_{ij} = \sigma^2 \rho^{|i-j|}$, with $\sigma = 1$ and $\rho = 0.5$.

(D2) $\mathbf{U} = b\mathbf{1}_{n \times n} + (a - b)\mathbf{I}_n$, with $a = 0.5$ and $b = 1$. $\mathbf{\Sigma}$ is a Toeplitz matrix, i.e. $\Sigma_{ij} = t(|i - j|)$, with $t(s) = 1 + 1/(1 + s)$ for $s \in \mathbb{N}$.

(D3) $\mathbf{U} = b\mathbf{1}_{n \times n} + (a - b)\mathbf{I}_n$, with $a = 0.2$ and $b = 2$. $\mathbf{\Sigma}$ is a diagonal matrix with diagonal entries given by $\Sigma_{ii} = 1 + 1/i$.
<figure><img src="image_4.png" /><figcaption>Figure 3: Empirical cumulative distribution functions (ECDF) of *p*-values (p-gen) with *C* being a hierarchical clustering algorithm with average linkage. The ECDF were computed from *M* = 2000 realizations of (gen-MN) under the three dependence settings (*D1*), (*D2*) and (*D3*) with *μ* = *0*<sub>*n*×*p*</sub>, *n* = 100 and *p* ∈ {5, 20, 50}.</figcaption></figure>

We simulated matrix normal data in settings (D1), (D2) and (D3) and performed *k*-means and hierarchical agglomerative clustering (HAC) with average, centroid, single and complete linkages. In Section 4.1 we illustrate the uniformity of the *p*-values (p-gen) under a global null hypothesis, assuming that both scale matrices are known. In Section 4.2, we consider the case where **U** is known and the covariance between features **Σ** is estimated. We show, as proved in Section 3.1, that *p*-values are super-uniform for large enough sample sizes. In Section 4.3, we assess the relative efficiency of the considered algorithms in terms of power, for the three dependence scenarios. Finally, in Section 4.4, we study the robustness of the proposed approach to model misspecification.

## 4.1 Uniform *p*-values under a global null hypothesis

To illustrate the null distribution of *p*-values, we followed the same steps as in [21, Section 5.1]. For *n* = 100 and *p* ∈ {5, 20, 50}, we simulated *M* = 2000 samples drawn from model (gen-MN) in settings (D1), (D2) and (D3) with $\mu = 0_{n \times p}$ a zero matrix, so that the null hypothesis (H0) holds for any pair of clusters in $C(\mathbf{X})$. For each simulated sample, we used *k*-means and HAC to estimate three clusters and tested (H0) for two randomly selected clusters. Results for HAC with average linkage are displayed in Figure 3, where the empirical cumulative distribution functions (ECDF) of the simulated *p*-values are shown. The results for *k*-means and HAC with centroid, single and complete linkage are analogous to those for average linkage and we present them in Appendix D.4. The *p*-values for HAC with complete linkage were computed as their Monte Carlo approximation (17) with *N* = 2000 iterations. In all cases, the *p*-values follow a uniform distribution when the null hypothesis (H0) holds.

## 4.2 Super-uniform *p*-values for unknown **Σ**

In this section, we illustrate that *p*-values (hat-p-tract) are asymptotically super-uniform under (H0) when **Σ** is asymptotically over-estimated in the sense of Loewner partial order, as proved in Theorem 3.1. We use the estimator (hat-Sigma) that asymptotically over-estimates **Σ** for $\mathbf{U} \in \mathcal{CS}(n)$ if Assumptions 3.1 and 3.2 hold. The estimate is computed using an independent and identically distributed copy of the sample where the clustering was performed, following Proposition 3.3.
<figure><img src="image_5.png" /><figcaption>Figure 4: Empirical cumulative distribution functions (ECDF) of p-values (hat-p-tract) with C being a hierarchical clustering algorithm with average linkage. The ECDF are computed from M = 5000 realizations of (gen-MN) under the three dependence settings (D1), (D2) and (D3) with n = 100, p = 5 and μ given by (40) with δ ∈ {6, 8}. Only samples for which the null hypothesis held were kept, as described in Section 4.2.</figcaption></figure>

We follow the same steps as in [21, Section D.1]. For n = 100 and p = 5, we simulate M = 5000 samples drawn from (gen-MN) in settings (D1), (D2) and (D3) with μ being divided into two clusters:

$$ \mu_{ij} = \begin{cases} \frac{\delta}{j} & \text{if } i \le \frac{n}{2}, \\ -\frac{\delta}{j} & \text{otherwise,} \end{cases} \quad \forall i \in [n], \forall j \in \{1, \dots, p\}, \qquad (40) $$

with δ ∈ {6, 8}. For k-means and HAC with average, centroid, single and complete linkage we set C to choose three clusters. The samples for which (H0) held when comparing two randomly selected clusters are kept. Results for HAC with average linkage are presented in Figure 4. The results for k-means and HAC with centroid, single and complete linkage are analogous and we present them in Appendix D.4. All simulations illustrate the asymptotic super-uniformity of p-values (p-gen) under the null hypothesis, when Σ is asymptotically over-estimated using (hat-Sigma). Moreover, as the distance between clusters δ decreases, the over-estimation is less severe and the null distribution of p-values approaches the one of a uniform random variable.

It is important to remark that Figure 4 serves only to illustrate the validity of Theorem 3.1, but in no way to interpret the conservativeness of p-values when Σ is over-estimated. The deviation from uniformity of the null distribution of (hat-p-tract) or, equivalently, the power of the corresponding test, depends on the measure of the conditioning set, which in Figure 4 is determined by the frequency of iterations satisfying (H0).

## 4.3 Power analysis

We now assess the relative efficiency of the five clustering algorithms considered in terms of power, as well as their power loss when one of the scale matrices is estimated using (hat-Sigma). As in [21, Section 5.2], we consider the conditional power of the p-value (p-gen), which is the probability of rejecting the null (H0) for a randomly selected pair of clusters given that they are different. To estimate the conditional power, we simulate M = 5000 samples drawn from (gen-MN) under the three settings (D1), (D2) and
(D3) with $\mu$ dividing the $n = 200$ observations into three true clusters:

$$ \mu_i = \begin{cases} (-\frac{\delta}{2}, 0, \dots, 0) & \text{if } i \le \lfloor \frac{n}{3} \rfloor, \\ (0, \dots, 0, \frac{\sqrt{3}\delta}{2}) & \text{if } \lfloor \frac{n}{3} \rfloor < i \le \lfloor \frac{2n}{3} \rfloor, \quad \forall i \in [n], \\ (\frac{\delta}{2}, 0, \dots, 0) & \text{otherwise,} \end{cases} \quad (41) $$

for $p = 5$ and for 13 evenly-spaced values of $\delta \in [4, 10]$. Then, we estimate the conditional power as the proportion of rejections at level $\alpha = 0.05$ among the samples for which the null hypothesis (H0) did not hold (which were above the 90% of $n$ in all settings). The conditional power as a function of $\delta$ is shown in Figure 5(a-c) for the three scenarios (D1), (D2) and (D3) and the five considered clustering algorithms. The $p$-values for HAC with complete linkage are estimated using the approximation (17) with $N = 2000$ iterations.

Figure 5(a-c) shows that, in all cases, conditional power increases with the distance between true clusters. Regarding HAC, we observe that average linkage presents the best relative efficiency among the four considered linkages in all the dependence settings, followed closely by complete linkage, which seems to weaken in (D2). This might suggest that conditional power depends on the scale matrices and some scenarios might strongly differ from the overall observed behavior. Indeed, the qualitative difference between average or complete linkage and centroid or single linkage that is observed in (D1) and (D3) considerably lessens in (D2). In (D1) and (D3), the performance of single linkage is undoubtedly the lowest, and large differences between clusters are required to attain satisfactory levels of conditional power. However, single linkage achieves one of the best performances in (D2).

The relative efficiency of the $k$-means algorithm in terms of conditional power is the best in (D2), but one of the worst among all the considered algorithms in (D1) and (D3) settings. These unsatisfactory performances might be explained by the behavior already pointed out by the authors in [9], who referred to the fact that conditioning on too much information entails a loss of power [7, 20, 27, 33]. Recall that the truncation set for $k$-means post-clustering inference defined in [7] is non-maximal to allow its efficient computation (see Appendix B and [9, Equation (9)]). This approach, although respecting the selective type I error as shown in Theorem B.1, might sacrifice the efficiency in terms of power of the corresponding test, as illustrated in Figure 5(a,c).

Next, we evaluate the loss of power entailed by estimating one of the scale matrices using (hat-Sigma). Recall that, following Theorem 3.1, the $p$-values (hat-p-tract) are asymptotically super-uniform under the null, so conditional power is expected to decrease due to both the estimation of unknown parameters and the conservativeness of the testing approach. We repeat the previously described analysis but replacing $\Sigma$ by its estimate (hat-Sigma), and calculate the counterparts of the curves in Figure 5(a-c) for $p$-values (hat-p-tract). They are shown in Figure 5(d-f). In Figure 5(g-i), we depict the loss of power in estimation, defined as the absolute difference of the conditional power computed with known and over-estimated $\Sigma$, for every fixed clustering algorithm and value of $\delta$.

Figure 5(g-i) illustrates how power loss varies substantially across settings (D1), (D2) and (D3). Overall, average and centroid linkages exhibit the slightest loss, falling below 10% for $\delta > 6$ in (D1) and (D2). A greater separation between clusters is required to achieve a reasonable power loss under (D3). The power loss curve of complete linkage closely resembles that of average and centroid linkages in (D1) and (D3), but takes substantially higher values in (D2). Conversely, single linkage shows a similar behavior to centroid and average linkages in (D2) but differs notably in (D1) and (D3). Once again, we find that the $k$-means algorithm exhibits the worst relative efficiency in terms of power loss, especially in (D1) and (D3). A similar behavior was observed in [9] for $k$-means clustering when over-estimating
<figure><img src="image_6.png" /><figcaption>Figure 5: (a-f): conditional power for the test proposed in Section 2 under model (gen-MN) with the three dependence settings (D1), (D2) and (D3) and the mean matrix defined in (41). The conditional power is estimated as the proportion of rejection at level $\alpha = 0.05$ among the subset of the $M = 5000$ realizations of (gen-MN) for which the null hypothesis (H0) holds. In (a-c), $\Sigma$ is known and in (d-f) it is over-estimated using (hat-Sigma). (g-i): power loss in estimation defined as the absolute difference of the curves in (a-c) and (d-f).</figcaption></figure>

$\sigma$ under (ind-MN) using the estimator proposed in [21]. This suggests that the unsatisfactory efficiency of post-k-means inference is intrinsic to the $p$-value defined in [9], and that the extension proposed here inherits that drawback. An alternative approach would be to explore the use of consistent estimators of $\Sigma$ under (gen-MN), which would reduce power loss as demonstrated in [9] for the simpler model (ind-MN). Following all panels in Figure 5, we can conclude that HAC with average linkage exhibits the highest relative efficiency and lower power loss when $\Sigma$ is estimated, making it the most suitable algorithm in practice. Note that substantial power loss in the estimation of unknown parameters was similarly observed in the methods proposed in [21, 50], as demonstrated in [50] for HAC algorithms under (ind-MN).
## 4.4 Robustness to model misspecification

We conclude the numerical simulations on synthetic data by studying the robustness of the proposed approach to model misspecification. We particularly evaluate settings where the theoretical constraints on the dependence between observations given by $\mathbf{U}$ are not satisfied or known. First, in Section 4.4.1, we analyze how $p$-values (p-gen) behave when the covariance matrix $\mathbf{U}$ is not compound symmetry, but is compatible with the over-estimation of $\mathbf{\Sigma}$. Then, in Section 4.4.2, we explore the setting where $\mathbf{U}$ does not fit into $\mathcal{CS}(n)$ nor belongs to any of the models stated in Remarks 3.2 or 3.3. Finally, in Section 4.4.3, we evaluate the validity of the method when $\mathbf{U} \neq \mathbf{I}_n$ is unknown and observations are assumed to be independent.

### 4.4.1 Non-compound-symmetry $\mathbf{U}$ structures

In this section we evaluate the robustness of $p$-values (p-gen) and (hat-p-tract) to $\mathbf{U} \notin \mathcal{CS}(n)$. We choose three dependence settings that satisfy Assumption 3.3, so that (hat-Sigma) satisfies (over-est). In all cases, $\mathbf{\Sigma}$ is a diagonal matrix with entries $\Sigma_{ii} = 1 + 1/i$. The dependence structure between observations is given by the three following settings:

(D4) $\mathbf{U}$ is a diagonal matrix with entries $U_{ii} = 1 + 1/i$.

(D5) $\mathbf{U}$ is the covariance matrix of an AR(1) model with $\sigma = 1$ and $\rho = 0.1$.

(D6) $\mathbf{U}$ is the covariance matrix of an AR(2) model with $\sigma = 1$, $\beta_1 = 0.4$ and $\beta_2 = 0.1$.

We start by simulating the distribution of $p$-values (p-gen) under the global null hypothesis, repeating the numerical experience described in Section 4.1. The counterpart of Figure 3 for (D4), (D5) and (D6) is presented in Figure 6. The empirical distribution of $p$-values does not markedly deviate from uniformity in settings (D4) and (D5), especially for $p \in \{5, 10\}$. This was expected since the $\mathbf{U}$ matrices in both cases do not deviate substantially from the compound symmetry structure. In (D6), the entries of $\mathbf{U}$ decay more slowly to zero along the columns, which makes this structure to deviate more from $\mathcal{CS}(n)$. This results in a greater departure from uniformity of the $p$-value distribution, as seen in Figure 6(c). However, this deviation occurs within the super-uniformity regime, meaning that the $p$-values still maintain statistical guarantees, despite the power loss. The corresponding results for $k$-means and HAC with centroid, single and complete linkages are analogous. We include them in Appendix D.4.

The previous analysis suggests that $p$-values (p-gen) are robust to small deviations from $\mathbf{U} \in \mathcal{CS}(n)$. As discussed in Section 3.2, if the over-estimate condition (over-est) of Theorem 3.1 is satisfied, this would mean that $p$-values (hat-p-tract) are equally robust in that setting. Following from Remarks 3.2 and 3.3, settings (D4), (D5) and (D6) are compatible with the asymptotic over-estimation of $\mathbf{\Sigma}$ using (hat-Sigma). Consequently, we reproduce the analyses of Section 4.2 for such dependence structures to assess whether the previously illustrated robustness is maintained with estimation. Results are presented in Figure 7 for HAC with average linkage and in Appendix D.4 for the remaining clustering algorithms. In all cases, the empirical null distribution of $p$-values is super-uniform, confirming the robustness of (hat-p-tract) to small deviations from $\mathbf{U} \in \mathcal{CS}(n)$.

### 4.4.2 Non-admissible $\mathbf{U}$ for the over-estimation of $\mathbf{\Sigma}$

Let us recall that Assumption 3.3 is a sufficient condition for the sequence $\{\mathbf{U}^{(n)}\}_{n \in \mathbb{N}}$ to ensure that (hat-Sigma) satisfies (over-est). As discussed in Section 3, proving that a given dependence model satisfies
<figure><img src="image_7.png" /><figcaption>Figure 6: Empirical cumulative distribution functions (ECDF) of p-values (p-gen) with C being a hierarchical clustering algorithm with average linkage. The ECDF were computed from $M = 2000$ realizations of (gen-MN) under the three dependence settings (D4), (D5) and (D6) with $\mu = 0_{n \times p}$, $n = 100$ and $p \in \{5, 20, 50\}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Empirical cumulative distribution functions (ECDF) of p-values (hat-p-tract) with C being a hierarchical clustering algorithm with average linkage. The ECDF are computed from $M = 5000$ realizations of (gen-MN) under the three dependence settings (D4), (D5) and (D6) with $n = 50$, $p = 5$ and $\mu$ given by (40) with $\delta \in \{4, 6, 8\}$. Only samples for which the null hypothesis held were kept, as described in Section 4.2.</figcaption></figure>

this Assumption is non-trivial in most cases. In Remarks 3.2, 3.1 and 3.3, we showed that Assumption 3.3 is satisfied by three common dependence structures, but other sequences ${U^{(n)}}_{n \in \mathbb{N}}$ might also satisfy the same sufficient condition or other unknown hypotheses that ensure that (hat-Sigma) asymptotically over-estimates $\Sigma$. In this section, we repeat the simulations of Section 4.2 under three settings that do not fit Remarks 3.1, 3.2 or 3.3:

$$ (D7) \quad U \text{ is a Toeplitz matrix with } U_{ij} = 1 + 1/(1 + |i-j|). $$

$$ (D8) \quad U \text{ is the covariance matrix of an AR(3) model with } \sigma = 1, \beta_1 = 0.4, \beta_2 = -0.2 \text{ and } \beta_3 = 0.1. $$

$$ (D9) \quad U \text{ is a banded matrix with } U_{ii} = 1, U_{ii+1} = 0.6, U_{ii+2} = 0.5, U_{ii+3} = 0.2 \text{ and } U_{ii+r} \text{ for all } r > 3. $$

In all cases, we chose $\Sigma$ as a diagonal matrix with entries $\Sigma_{ii} = 1 + 1/i$. We also set $n = 50$, $p = 5$ and $\delta \in \{4, 6, 8\}$. Results are presented in Figure 8 for HAC with average linkage, and in Appendix D.4 for the rest of clustering algorithms. The simulated p-values are super-uniform in all settings, suggesting that
<figure><img src="image_9.png" /><figcaption>Figure 8: Empirical cumulative distribution functions (ECDF) of p-values (hat-p-tract) with C being a hierarchical clustering algorithm with average linkage. The ECDF are computed from $M = 5000$ realizations of (gen-MN) under the three dependence settings (D7), (D8) and (D9) with $n = 50$, $p = 5$ and $\mu$ given by (40) with $\delta \in \{4, 6, 8\}$. Only samples for which the null hypothesis held were kept, as described in Section 4.2.</figcaption></figure>

(hat-Sigma) might asymptotically over-estimate $\Sigma$ for further models of dependence between observations. Note that, in particular, results corresponding to (D8) suggest that the requirement $\beta_k\beta_{k'} \ge 0$ for $P > 2$ in Remark 3.3 is not very restrictive.

These results might also motivate further theoretical inspection on Toeplitz and banded structures to verify whether they satisfy Assumption 3.3. Extensive work has been done on the asymptotic behavior of continuous functions of Toeplitz matrices [22]. However, it mainly concerns their average behavior rather than their element-wise one. Notably, in [22], it is proved that the mean of the eigenvalues of $(U^{(n)})^{-1}$ converges when $n$ tends to infinity, if the sequence $\{U_{1n}\}_{n \in \mathbb{N}}$ is absolutely summable. This implies that the mean of the sequence $\{(U^{(n)})_{ii}^{-1}\}_{i=1,...,n}$ also converges with $n$. However, this is insufficient to state convergence of (33) and the asymptotic behavior of the individual entries need to be studied. If we impose $U^{(n)}$ to be banded, the entry-wise convergence of the elements $(U^{(n)})_{i+r}^{-1}$ has been demonstrated in [12] for the tridiagonal case. This, together with the exponential decay of the entries of banded matrices [13], is enough to prove the first part of Assumption 3.3 for tridiagonal Toeplitz matrices. Unfortunately, the existing results do not ensure that any of the conditions (i) or (ii) in Assumption 3.3 hold. Assessing that remaining step is mathematically very challenging and it is left for future work.

### 4.4.3 Ignoring weak dependence between observations

In real applications, it might be common that the practitioner lacks knowledge of both dependence structures between observations and variables. As discussed in Section 3, simultaneous estimation of both matrices $U$ and $\Sigma$ is unfeasible under the matrix normal model (gen-MN) when only one or few copies of $X$ are available. Consequently, even ignoring the control of statistical guarantees, we are unable to simultaneously consider a pair of estimators $\hat{U}$, $\hat{\Sigma}$ (or one of the Kronecker product $U \otimes \Sigma$) in the context of this work. In practice, a common alternative strategy is to assume weak dependence between observations, and ignore this dependence by considering $U = I_n$ in the method. In this section, we study the robustness of the proposed approach when observations are supposed independent but it is known that $U \ne I_n$.
<figure><img src="image_10.png" /><figcaption>Figure 9: Empirical cumulative distribution functions (ECDF) of p-values (hat-p-tract) with C being a hierarchical clustering algorithm with average linkage. The ECDF are computed from M = 5000 realizations of (gen-MN) as described in Section 4.4.3 with n = 50, p = 5 and μ given by (40) with $\delta \in \{4, 6, 8\}$. Only samples for which the null hypothesis held were kept, as described in Section 4.2.</figcaption></figure>

We consider $\mathbf{X}$ drawn from (gen-MN) with $\Sigma$ a diagonal matrix having as entries $\Sigma_{ii} = 1 + 1/i$ as in the previous section. The dependence between observations is encoded by the covariance matrix of an AR(1) model, that is, $U_{ij} = \sigma^2 \rho^{|i-j|}$, with $\sigma = 1$ and $\rho \in \{0.1, 0.2, 0.3, 0.4, 0.5\}$. Once again, we repeated the simulations described in Section 4.2 and computed the p-values (hat-p-tract) using (hat-Sigma) to estimate $\Sigma$ and assuming $\mathbf{U} = \mathbf{I}_n$. Results for HAC with average linkage are presented in Figure 9, and in Appendix D.4 for the rest of clustering algorithms. In all cases, the simulated p-values do not substantially deviate from the super-uniform regime. Besides, if we take a closer look at $[0, 0.1]$, we see that the simulated ECDF strictly lie below the diagonal for small values of $\rho$. In other words, when the dependence between observations is weak, the proposed test is robust to departures from the assumption $\mathbf{U} = \mathbf{I}_n$, and the estimation of $\Sigma$ using (hat-Sigma) yields p-values that asymptotically control the selective type I error.

# 5 Application to clustering of protein structures

Proteins are essential molecules in all living organisms. Many of their numerous functions are closely related to their non-static structure, which exhibits high variability within numerous protein families [18, 32, 37]. The characterization of such intrinsic structural complexity represents a highly active area of research in the field of Structural Biology. In this pursuit, clustering methods applied to protein conformations have provided valuable insights into this challenging problem [3, 11]. One of the most commonly-chosen descriptors to characterize a protein conformation is the set of pairwise Euclidean distances between every pair of amino acids along the sequence [30, 35, 40], usually referred to as distance maps. As these distances are strongly correlated, assuming a constant diagonal covariance matrix as in [21] seems very unrealistic. Instead, we opt for the more convenient model

$$ \mathbf{X} \sim \mathcal{MN}_{n \times p}(\mu, \mathbf{I}_n, \Sigma), \qquad (42) $$

where $\Sigma$ can be estimated using (hat-Sigma). Each row of $\mathbf{X}$ corresponds to a protein conformation, featured by a vector of Euclidean distances between every pair of amino acids, which constitute the
<figure><img src="image_11.png" /><figcaption>Figure 10: Average pairwise distances between every pair of amino acids across the conformations of each cluster. The clusters were found after performing hierarchical clustering with average linkage on the protein data presented in Section 5.</figcaption></figure>

columns of X. We perform hierarchical agglomerative clustering with average linkage (as it showed the best relative efficiency in Section 4.3) to estimate $k = 6$ clusters among $n = 2000$ conformations of a disordered protein called Histatin-5 (Hst5). The number of clusters was chosen arbitrarily. The corresponding sequence is 24 amino acids long, so $p = 23 \cdot 24/2 = 276$. The conformations were generated using Flexible-Meccano [5,38] and refined using previously reported small-angle X-ray scattering (SAXS) data [42]. Note that Flexible-Meccano is a sampling algorithm that generates an independent conformation at each iteration, contrary to Molecular Dynamics simulation techniques that present temporal dependence between samples. This justifies our choice of $U = I_n$. Moreover, we had access to an independent replica of the simulated ensemble that we used to estimate $\Sigma$, as it is usual for generated protein ensembles. The obtained estimate $\hat{\Sigma}$ substantially deviated from the spherical structure. Figure 10 shows the average distance map across all conformations in a given cluster or, in other words, the empirical cluster means as defined in (1). Table 1 presents the p-values corresponding to every pair of clusters, corrected for multiple testing using the Bonferroni-Holm adjustment [25].

<table><caption>Table 1: p-values (p-gen) computed under model (42) retrieved after testing (H0) on the protein data presented in Section 5. The hierarchical clustering algorithm was set to find six clusters using average linkage. In blue, adjusted p-values for which the null is not rejected at level $\alpha = 0.05$.</caption><thead><tr><th>Cluster</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>2</td><td>2.187589·10<sup>-4</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>3</td><td>3.039844·10<sup>-11</sup></td><td>1.41·10<sup>-3</sup></td><td></td><td></td><td></td></tr><tr><td>4</td><td>1.070993·10<sup>-10</sup></td><td>0.300540</td><td>2.98464·10<sup>-4</sup></td><td></td><td></td></tr><tr><td>5</td><td>3.038979·10<sup>-16</sup></td><td>0.093018</td><td>6.015797·10<sup>-5</sup></td><td>0.105446</td><td></td></tr><tr><td>6</td><td>1.729616·10<sup>-6</sup></td><td>0.010612</td><td>9.290826·10<sup>-9</sup></td><td>2.105·10<sup>-3</sup></td><td>5.624624·10<sup>-5</sup></td></tr></tbody></table>
The p-values presented in Table 1 show significant differences between the most part of the average distance maps depicted in Figure 10. The non-rejecting pairs of clusters at level $\alpha = 0.05$, marked in blue in Table 1, suggest that clusters 2, 4 and 5 could be merged into a single group. Indeed, when looking at the corresponding empirical means in Figure 10, we appreciate that these three clusters are characterized by large distances between pairs of amino acids that are far apart in the sequence, which indicates a lack of interactions between the sequence termini and a more extended structure of the corresponding conformations. This feature appears as an exclusive and prominent characteristic of clusters 2, 4 and 5, which might explain the non-rejection of the corresponding nulls. For the rest of rejecting pairs of clusters, clear differences in distance patterns are retrieved in Figure 10, accounting for significant changes on Hst5 structure between the corresponding groups. The results presented in Table 1 are coherent with the HAC dendrogram, presented in Figure C.1, showing that clusters 2, 4, and 5 form a subgroup that is promptly separated from the rest.

# 6 Discussion

The seminal work by Gao et al. [21] has laid the foundation for selective inference after clustering by introducing a theoretical framework allowing to test differences between cluster means, conditioning on having estimated those clusters. Furthermore, the authors have tackled the problem of estimating unknown parameters while controlling the selective type I error, which had been overlooked in previous works [31,41], but which is crucial for the practical application of this theory. Their contribution motivates extensions of post-clustering inference to more general frameworks that arise in complex real applications, where observations or features present non-negligible dependence structures. To generalize the model considered in [21] to the more general (gen-MN), we consider a p-value of the form (p-GBW), choosing a test statistic based on $X^T\nu$ and conditioning on both its direction and the projection $\pi_\nu^\perp X$, as done in [21]. In that setting, we prove that the strategy of [21] can be extended to (gen-MN) if and only if the dependence structure between observations $U$ is compound symmetry. Otherwise, we show that the natural generalization of (p-GBW) to arbitrary $U$ yields a quantity that can be efficiently characterized, but whose statistical guarantees are difficult to assess. Numerically, we illustrate that the control of the selective type I error is not ensured in that setting. We also generalize the estimation of one covariance matrix compatible with the selective type I error control when $U \in \mathcal{CS}(n)$. These extensions, presented in Sections 2 and 3 respectively, and numerically illustrated in Sections 4 and 5, represent the main contributions of this work.

The theoretical framework presented in Section 2 limits the use of p-values of the form (p-GBW) to structures $U \in \mathcal{CS}(n)$. Following from the analyses that we present in Section 2.2.3, generalizing the family of admissible $U$ is a complex problem in this context and would require exploring p-values with alternative conditioning sets. As we have suggested, such a strategy would require the definition of extra conditioning events that are *independent* of the test statistic. According to Proposition 2.1, this would mean to replace the projection $\pi_\nu^\perp$ by one that is independent of $X^T\nu$ for any $U$. If the projection is taken with respect to the scalar product defined by $U$, that is,

$$ \pi_{U;\nu}^{\perp} X = X - \frac{X^T U \nu}{\nu^T U \nu} \nu, \qquad (43) $$

the independence $\pi_{U;\nu}^{\perp} X \perp X^T \nu$ follows from the Cochran Theorem [10]. However, replacing (43) in (6) and proceeding with the same reasoning would mean to consider a test statistic based on $X^T U \nu$, that would account for a less interpretable null hypothesis of the form $\mu^T U \nu = 0$. Besides, maintaining both the projection (43) and the null hypothesis (H0) would substantially complicate the derivation of
a tractable *p*-value. For this reason, we believe that extending the conditional post-clustering inference approaches to arbitrary structures **U** would require a substantial shift in framework and should follow alternative paths to the strategy initiated in [21]. Recall, however, that the method proposed here has been shown to be robust to **U** *notin* CS(*n*) in several scenarios.

The estimation of unknown parameters, which is essential for practical applications, inherently leads to a loss of power in any hypothesis test. This has been illustrated in Section 4.3 for the method proposed here. A relevant avenue for future work would be the exploration of alternative scenarios where the power loss in estimation could be mitigated. One possibility would be to develop a framework inspired by the work of Yun and Foygel Barber [50], in which they consider, under model (ind-MN), a test statistic that does not depend on the unknown parameter $\sigma$. This results in a method that is relatively more efficient than the one proposed in [21] in some settings. Adapting this idea to the general model (gen-MN) would require an appropriate test statistic that does not depend on the unknown parameters **U** and $\Sigma$. However, the direct adaptation of [50] to (gen-MN) presents a non-trivial theoretical challenge while offering limited practical advantages compared to the extension presented here. Indeed, an efficient computation is proposed only for binary partitions of the data. An alternative approach would be the definition of consistent estimators of **U** or $\Sigma$ that are compatible with the selective type I error control. This was studied in [9] in the context of *k*-means clustering under (ind-MN), where the authors showed that considering a median-based consistent estimator of $\sigma$ yields better performances than the over-estimation strategy proposed in [21].

Clustering is a multidimensional method that incorporates information from *p* descriptors to classify *n* observations. However, the estimated groups are often distinguished by a subset of variables, whose determination is essential in various fields of application [36, 46]. The framework presented in [21] has also been adapted to feature-level post-clustering inference [8, 24], testing for the difference of the *g*-th coordinate of cluster means, for a fixed $g \in \{1, \dots, p\}$. In that case, clustering is performed on the complete data set **X** but inference is carried out on the *g*-th column, modeled by a *n*-dimensional Gaussian. In a recent contribution [8], the covariance matrix is let arbitrary and *p*-values can be efficiently computed following a similar reasoning as in [21]. Nevertheless, none of these works deal with the estimation of unknown parameters. The extension of the over-estimation strategy presented in Section 3 to this framework is non-trivial, and would represent a very relevant line for future research.

As discussed in Appendix B, performing analytically tractable post-clustering inference requires the addition of technical events to the conditioning set, which implies a reduction in power. Investigating whether these conditions might be relaxed is an interesting path for future research. The problem of power loss due to extra conditioning is not exclusive to this method. Techniques like data fission [31] need to calibrate the conditioning information and consequences in terms of power are analogous. However, it is still unknown whether power loss is more drastic in one method or the other. An interesting contribution would be to establish a framework allowing for a proper comparison of this effect when performing post-clustering inference using data fission and the approach proposed in [21]. Nevertheless, extending this comparison to practical applications would be unfeasible as long as the estimation of the covariance structure with statistical guarantees cannot be carried out in both methods.

## Code availability

The methods introduced in the present work were implemented in the R package *PCIdep*, available at https://github.com/gonzalez-delgado/PCIdep. All the numerical experiments on synthetic and real data
can be reproduced with the code available at https://github.com/gonzalez-delgado/PCIdep-experiments. The dataset of protein structures used in Section 5 can be downloaded at https://doi.org/10.5281/zenodo.10021202.

# Acknowledgments

We thank Amin Sagar and Pau Bernadó for providing protein structure data.

This work was supported by the French National Research Agency (ANR) under grants: ANR-11-LABX-0040 (LabEx CIMI) within the French State Program “Investissements d’Avenir”, ANR-22-CE45-0003 (CORNFLEX project) and ANR-21-CE40-0007 (GAP project).

# A Proofs

## A.1 Proofs of Section 2

### A.1.1 Proof of Proposition 2.1

We begin by recalling a useful established result. We then state and prove Lemma A.2, which is essential for the proof of Proposition 2.1, presented at the end of the section.

**Lemma A.1 (Proposition 3.4 in [19]).** Let $y \sim N(0, S)$ be a $p$-dimensional non-degenerated Gaussian vector and $F \subset \mathbb{R}^p$ a vector subspace. We denote by $\mathbf{P}_F$ the orthogonal projection on $F$ and by $\mathbf{P}_F^\perp$ the orthogonal projection on $F^\perp$. Then, $\mathbf{SP}_F = \mathbf{P}_F S$ if and only if the Gaussian vectors $\mathbf{P}_F y$ and $\mathbf{P}_F^\perp y$ are independent.

**Lemma A.2.** Let $\mathbf{T}$ be a $n \times n$ positive definite symmetric matrix. Then, $\mathbf{T} \in \mathcal{CS}(n)$ if and only if $\nu_{\mathcal{G}_1, \mathcal{G}_2}$ is an eigenvector of $\mathbf{T}$ for all $(\mathcal{G}_1, \mathcal{G}_2) \in \mathcal{C}_{[n]}$.

*Proof of Lemma A.2.* Let $\mathbf{T} = (a-b)\mathbf{I}_n + b\mathbf{1}_{n \times n} \in \mathcal{CS}(n)$. Then, $\mathbf{T}\nu_{\mathcal{G}_1, \mathcal{G}_2} = (a-b)\nu_{\mathcal{G}_1, \mathcal{G}_2}$ as $\mathbf{1}_{n \times n}\nu_{\mathcal{G}_1, \mathcal{G}_2} = \mathbf{0}_n$ for any $(\mathcal{G}_1, \mathcal{G}_2) \in \mathcal{C}_{[n]}$. To prove the reciprocal implication, we first define the set

$$ \mathcal{C}_P = \{(\mathcal{G}_1, \mathcal{G}_2) | \mathcal{G}_1, \mathcal{G}_2 \subset \mathcal{P}, \mathcal{G}_1 \cap \mathcal{G}_2 = \emptyset\}, $$

for any $\mathcal{P} \subset [n]$. Then, we prove the following proposition by induction over $k \ge 2$:

For any $\mathcal{P} \subset [n]$ with $2 \le |\mathcal{P}| \le k$, if $\nu_{\mathcal{G}_1, \mathcal{G}_2}$ is an eigenvector of $\mathbf{T}$ for all $(\mathcal{G}_1, \mathcal{G}_2) \in \mathcal{C}_P$, then the restriction of $\mathbf{T}$ on $F_P := \operatorname{span}\{\nu_{\mathcal{G}_1, \mathcal{G}_2} : (\mathcal{G}_1, \mathcal{G}_2) \in \mathcal{C}_P\}$ is a uniform scaling. (H$_k$)

* Initialization ($k = 2, 3$). If $\mathcal{P} = \{p_1, p_2\}$, ($H_2$) holds as $F_P = \operatorname{span}\{\nu_{\{p_1\}, \{p_2\}}\}$ and $\nu_{\{p_1\}, \{p_2\}}$ is an eigenvector of $\mathbf{T}$. The same strategy yields ($H_3$).

* Induction. Let $(H_k)$ be true for $3 < k < n$. Let $\mathcal{P} \subset [n]$ with $|\mathcal{P}| = k + 1$ and assume that $\nu_{\mathcal{G}_1, \mathcal{G}_2}$ is an eigenvector of $\mathbf{T}$ for any $(\mathcal{G}_1, \mathcal{G}_2) \in \mathcal{C}_P$. Consider also $\mathcal{P}_1, \mathcal{P}_2 \subset \mathcal{P}$ with $|\mathcal{P}_1| = |\mathcal{P}_2| = k$ and $\mathcal{P}_1 \ne \mathcal{P}_2$. Note that from previous assumptions we have $\mathcal{P}_1 \cup \mathcal{P}_2 = \mathcal{P}$. Now, since $\mathcal{C}_{\mathcal{P}_1}$ and $\mathcal{C}_{\mathcal{P}_2}$ are subsets of $\mathcal{C}_P$, property $(H_k)$ ensures that the restrictions of $\mathbf{T}$ on $F_{\mathcal{P}_1}$ and $F_{\mathcal{P}_2}$ are uniform scalings, that is,

$$ \mathbf{T}_{|F_{\mathcal{P}_1}} = \lambda_{\mathcal{P}_1} \mathbf{I}_n \quad \text{and} \quad \mathbf{T}_{|F_{\mathcal{P}_2}} = \lambda_{\mathcal{P}_2} \mathbf{I}_n \quad \text{for some} \quad \lambda_{\mathcal{P}_1}, \lambda_{\mathcal{P}_2} \in \mathbb{R}. $$
Moreover, as $|\mathcal{P}_1 \cap \mathcal{P}_2| = k - 1 \ge 2$, there exist two distinct elements $i_1$ and $i_2$ in the intersection $\mathcal{P}_1 \cap \mathcal{P}_2$. Then, $\nu_{\{i_1\},\{i_2\}} \in F_{\mathcal{P}_1} \cap F_{\mathcal{P}_2}$. Since $F_{\mathcal{P}_1}$ and $F_{\mathcal{P}_2}$ share a non-zero element, we have $\lambda_{\mathcal{P}_1} = \lambda_{\mathcal{P}_2}$. We conclude the induction step noticing that $F_{\mathcal{P}} = F_{\mathcal{P}_1} + F_{\mathcal{P}_2}$ (by inclusion and dimensional argument).

* Conclusion. The property $(H_k)$ is initialized and inductive, then true for any $2 \le k \le n$.

Following from the previous reasoning, $(H_n)$ is true. Then, if $\nu_{\mathcal{G}_1, \mathcal{G}_2}$ is an eigenvector of $\mathbf{T}$ for all $(\mathcal{G}_1, \mathcal{G}_2) \in \mathcal{C}_{[n]}$, the restriction $\mathbf{T}|_{F_{[n]}}$ is a uniform scaling of parameter $\lambda$. Moreover, as $\mathbf{T}$ is symmetric, both $F_{[n]}$ and its orthogonal $F_{[n]}^\perp$ are stable under $\mathbf{T}$. It can be easily shown that $F_{[n]}^\perp = \text{span}\{\mathbf{1}_n\}$. Then, $\mathbf{1}_n$ is an eigenvector of $\mathbf{T}$, whose associated eigenvalue will be denoted by $\beta$. Noting that $n^{-1}\mathbf{1}_{n \times n} = n^{-1}\mathbf{1}_n \cdot \mathbf{1}_n^T$ is the orthogonal projection over $\text{span}\{\mathbf{1}_n\}$ and $\mathbf{I}_n - n^{-1}\mathbf{1}_{n \times n}$ is the orthogonal projection over $F_{[n]}$, we can write:

$$ \begin{aligned} \mathbf{T} &= \mathbf{T}(n^{-1}\mathbf{1}_{n \times n} + \mathbf{I}_n - n^{-1}\mathbf{1}_{n \times n}) = \mathbf{T}|_{F_{[n]}^\perp} n^{-1}\mathbf{1}_{n \times n} + \mathbf{T}|_{F_{[n]}}(\mathbf{I}_n - n^{-1}\mathbf{1}_{n \times n}) \\ &= \beta n^{-1}\mathbf{1}_{n \times n} + \lambda(\mathbf{I}_n - n^{-1}\mathbf{1}_{n \times n}) = (\lambda\mathbf{I}_n + n^{-1}(\beta - \lambda)\mathbf{1}_{n \times n}) \in \mathcal{CS}(n), \end{aligned} $$

concluding the proof.

<figure><img src="image_12.png" /></figure>

*Proof of Proposition 2.1.* We start showing the first equivalence in Proposition 2.1. Let us denote by $\Lambda \subset \mathbb{R}^{n \times p}$ the kernel of the linear mapping $\nu_{\mathcal{G}_1, \mathcal{G}_2}^T : \mathbf{M} \in \mathbb{R}^{n \times p} \mapsto \nu_{\mathcal{G}_1, \mathcal{G}_2}^T \mathbf{M}$ and by $\Lambda^\perp$ its orthogonal complement. We omit their dependence on $\mathcal{G}_1, \mathcal{G}_2$ for the sake of a simpler notation. Next, we denote by $\mathbf{\Pi}_\Lambda := \pi_{\nu_{\mathcal{G}_1, \mathcal{G}_2}}^\perp$ and $\mathbf{\Pi}_{\Lambda^\perp} := \mathbf{I} - \mathbf{\Pi}_\Lambda = \nu_{\mathcal{G}_1, \mathcal{G}_2}^T \nu_{\mathcal{G}_1, \mathcal{G}_2} / \|\nu_{\mathcal{G}_1, \mathcal{G}_2}\|^2$ the orthogonal projections on $\Lambda$ and $\Lambda^\perp$, respectively. Then, for all $(\mathcal{G}_1, \mathcal{G}_2) \in \mathcal{C}_{[n]}$, we have:

$$ \begin{aligned} \mathbf{X}^T \nu_{\mathcal{G}_1, \mathcal{G}_2} \perp \pi_{\nu_{\mathcal{G}_1, \mathcal{G}_2}}^\perp \mathbf{X} &\iff \mathbf{\Pi}_{\Lambda^\perp} \mathbf{X} \perp \mathbf{\Pi}_\Lambda \mathbf{X} \\ (\text{By Lemma A.1}) &\iff (\mathbf{\Sigma} \otimes \mathbf{U})(\mathbf{I}_p \otimes \mathbf{\Pi}_{\Lambda^\perp}) = (\mathbf{I}_p \otimes \mathbf{\Pi}_{\Lambda^\perp})(\mathbf{\Sigma} \otimes \mathbf{U}) \\ &\iff \mathbf{\Sigma} \otimes \mathbf{U} \mathbf{\Pi}_{\Lambda^\perp} = \mathbf{\Sigma} \otimes \mathbf{\Pi}_{\Lambda^\perp} \mathbf{U} \\ (\text{By injectivity of the mapping } \mathbf{M} \mapsto \mathbf{\Sigma} \otimes \mathbf{M}) &\iff \mathbf{U} \mathbf{\Pi}_{\Lambda^\perp} = \mathbf{\Pi}_{\Lambda^\perp} \mathbf{U} \\ &\iff \text{The eigenspaces of } \mathbf{\Pi}_{\Lambda^\perp} \text{ are stable under } \mathbf{U} \\ &\iff \text{span } \nu_{\mathcal{G}_1, \mathcal{G}_2} \text{ and } \nu_{\mathcal{G}_1, \mathcal{G}_2}^\perp \text{ are stable under } \mathbf{U} \\ &\iff \nu_{\mathcal{G}_1, \mathcal{G}_2} \text{ is an eigenvector of } \mathbf{U}. \end{aligned} $$

In the last equivalence, we consider the matrix $\mathbf{\Pi}_{\Lambda^\perp}$ as a linear operator on $\mathbb{R}^{n \times 1}$. As this holds for every $(\mathcal{G}_1, \mathcal{G}_2)$ in $\mathcal{C}_{[n]}$, equivalence (i) in Proposition 2.1 follows directly from Lemma A.2.

The second equivalence in Proposition 2.1 is a consequence of the following well-known result. For any $p$-dimensional Gaussian vector $z \sim \mathcal{N}(\mu, \mathbf{A})$,

$$ \|z\|_2 \perp \text{dir}(z) \iff \mu = 0 \text{ and } \mathbf{A} = \lambda \mathbf{I}_p \text{ for some } \lambda > 0. \quad (44) $$
Let $y \sim \mathcal{N}(0, S)$ be a $p$-dimensional Gaussian vector and consider $z = \sqrt{\mathbf{A}^{-1}}y$, for any $p \times p$ positive definite matrix $\mathbf{A}$. Then, $||y||_{\mathbf{A}} = ||z||_2$ and $z \sim \mathcal{N}(0, \sqrt{\mathbf{A}^{-1}\mathbf{S}\sqrt{\mathbf{A}^{-1}}})$. Consequently, we have:

$$ ||y||_{\mathbf{A}} \perp \mathrm{dir}_{\mathbf{A}}(y) \iff ||z||_2 \perp \frac{y}{||z||_2} $$

$$ (\mathbf{M} \mapsto \sqrt{\mathbf{A}^{-1}}\mathbf{M} \text{ is a one-to-one mapping}) \iff ||z||_2 \perp \mathrm{dir}(z) $$

$$ (\text{Equivalence (44)}) \iff \sqrt{\mathbf{A}^{-1}\mathbf{S}}\sqrt{\mathbf{A}^{-1}} = \lambda\mathbf{I}, \text{ for some } \lambda > 0. $$

$$ (\mathbf{M} \mapsto \sqrt{\mathbf{M}^{-1}\mathbf{S}}\sqrt{\mathbf{M}^{-1}} \text{ is one-to-one on positive matrices}) \iff \mathbf{A} = \lambda\mathbf{S}, \text{ for some } \lambda > 0. $$

Setting $\mathbf{S} = \mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}$ yields the result. $\square$

### A.1.2 Proofs of Section 2.2.2

*Proof of Theorem 2.2.* We follow the steps of the proof of Theorem 1 in [21]. We begin by deriving the null distribution of the test statistic $||\mathbf{X}^T\nu||_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}$ under (H0). First, from [23, Theorem 2.3.10], we have:

$$ \mathbf{X}^T \nu \stackrel{H_0^{\{\mathcal{G}_1, \mathcal{G}_2\}}}{\sim} \mathcal{N}_p(0, \mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}), \qquad (45) $$

which yields

$$ ||\mathbf{X}^T \nu||_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}} \stackrel{H_0^{\{\mathcal{G}_1, \mathcal{G}_2\}}}{\sim} \chi_p, \qquad (46) $$

where the norm $||\cdot||_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}$ is defined in (8). Let us now build the $p$-value for $H_0^{\{\mathcal{G}_1, \mathcal{G}_2\}}$, by slightly adapting the reasoning in [21]. On one hand, for any $\nu \in \mathbb{R}^n$, we have

$$ \mathbf{X} = \pi_{\nu}^{\perp} \mathbf{X} + (\mathbf{I}_n - \pi_{\nu}^{\perp} \mathbf{X}) = \pi_{\nu}^{\perp} \mathbf{X} + \left( \frac{||\mathbf{X}^T \nu||_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}}{||\nu||_2^2} \right) \nu \mathrm{dir}_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{X}^T \nu)^T. \qquad (47) $$

On the other hand, from Proposition 2.1 we have $\pi_{\nu}^{\perp} \mathbf{X} \perp \mathbf{X}^T \nu$, which implies $||\mathbf{X}^T \nu||_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}} \perp \pi_{\nu}^{\perp} \mathbf{X}$, and $||\mathbf{X}^T \nu||_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}} \perp \mathrm{dir}_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{X}^T \nu)$. We can now plug (47) in the definition of (p-gen) and, taking into account the previous independence relationships, we can write:

$$ p_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\}) = \mathbb{P}_{H_0^{\{\mathcal{G}_1, \mathcal{G}_2\}}} \left( \begin{aligned} & (||\mathbf{X}^T \nu||_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}} \ge ||\mathbf{x}^T \nu||_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}) \\ & ||\mathbf{X}^T \nu||_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}} \in S_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\}) \end{aligned} \right), \qquad (48) $$

where the set $S_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\})$ is defined in (10). Consequently, if we denote by $\mathbb{F}_p(t, S)$ the cumulative distribution function of a $\chi_p$ random variable truncated to the set $S$, from (48) and (46) we have

$$ p_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\}) = 1 - \mathbb{F}_p(||\mathbf{x}^T \nu||_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}, S_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\})), \qquad (49) $$

which proves the first statement (p-tract). The control of selective type I error is proved identically to the reasoning in the proof of [21, Theorem 1]. $\square$

*Proof of Lemma 2.3.* Let us first show that the perturbed data sets $\mathbf{x}'(\phi)$, defined in [21, Equation (13)] and $\mathbf{x}'_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\phi)$, defined in (12) are the same up to a scale transformation, i.e. that

$$ \mathbf{x}'_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\phi) = \mathbf{x}' \left( \frac{||\mathbf{x}^T \nu||_2}{||\mathbf{x}^T \nu||_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}} \right) \quad \forall \phi \ge 0. \qquad (50) $$
Note first that we can write

$$ \left( \frac{\|\mathbf{x}^T \nu\|_2}{\|\mathbf{x}^T \nu\|_{\mathbf{V}_{G_1, G_2}}} \phi - \|\mathbf{x}^T \nu\|_2 \right) \mathrm{dir}(\mathbf{x}^T \nu) = \left( \phi - \|\mathbf{x}^T \nu\|_{\mathbf{V}_{G_1, G_2}} \right) \mathrm{dir}_{\mathbf{V}_{G_1, G_2}}(\mathbf{x}^T \nu), \quad (51) $$

where $\mathrm{dir}(u) = u/\|u\|_2 \mathbf{1}\{u \neq 0\}$. Replacing (51) in (14), we have (50). Finally, it suffices to remark that

$$ \begin{aligned} \hat{\mathbf{S}}_{\mathbf{V}_{G_1, G_2}} &= \{\phi \ge 0 : G_1, G_2 \in \mathcal{C}(\mathbf{x}'_{\mathbf{V}_{G_1, G_2}}(\phi))\} = \{\phi \ge 0 : G_1, G_2 \in \mathcal{C}\left(\mathbf{x}'\left(\frac{\|\mathbf{x}^T \nu\|_2}{\|\mathbf{x}^T \nu\|_{\mathbf{V}_{G_1, G_2}}} \phi\right)\right)\} \\ &= \left\{\frac{\|\mathbf{x}^T \nu\|_{\mathbf{V}_{G_1, G_2}}}{\|\mathbf{x}^T \nu\|_2} \phi : G_1, G_2 \in \mathcal{C}(\mathbf{x}'(\phi))\right\} = \frac{\|\mathbf{x}^T \nu\|_{\mathbf{V}_{G_1, G_2}}}{\|\mathbf{x}^T \nu\|_2} \hat{\mathbf{S}}, \end{aligned} $$

which concludes the proof.

<figure><img src="image_13.png" /></figure>

### A.1.3 Proofs of Section 2.2.3

We start by stating some technical results that are needed for the proof of Theorem 2.4. In what follows, we will use the notation $\mathbf{S}$ to denote both a $d \times d$ real matrix and its associated linear mapping, that is, the map $\mathbf{S} : \mathbb{R}^d \to \mathbb{R}^d$ such that $\mathbf{S}(y) = \mathbf{S}y$ for all $y \in \mathbb{R}^d$. For any vector subspace $F \subset \mathbb{R}^d$, we will denote by $\mathbf{\Pi}_F$ the orthogonal projection onto $F$.

**Theorem A.3 (Proposition 3.13 in [19]).** Let $\mathbf{X} \sim \mathcal{N}(\mathbf{\mu}, \mathbf{\Sigma})$ be a Gaussian vector in $\mathbb{R}^n$, let $\mathbf{A}$ be a matrix in $M_{p,n}(\mathbb{R})$. Then, the conditional vector $(\mathbf{X}|\mathbf{A}\mathbf{X} = \mathbf{y})$ is a Gaussian vector satisfying

$$ (\mathbf{X}|\mathbf{A}\mathbf{X} = \mathbf{y}) \sim \mathcal{N}(\mathbf{\mu} + \mathbf{\Sigma}\mathbf{A}^T(\mathbf{A}\mathbf{\Sigma}\mathbf{A}^T)^{\dagger}(\mathbf{y} - \mathbf{A}\mathbf{\mu}), \mathbf{\Sigma} - \mathbf{\Sigma}\mathbf{A}^T(\mathbf{A}\mathbf{\Sigma}\mathbf{A}^T)^{\dagger}\mathbf{A}\mathbf{\Sigma}), \quad (52) $$

where $(\mathbf{A}\mathbf{\Sigma}\mathbf{A}^T)^{\dagger}$ is the Moore-Penrose pseudoinverse of the matrix $\mathbf{A}\mathbf{\Sigma}\mathbf{A}^T$.

**Lemma A.4.** Let $F, G$ be two orthogonal subspaces of $\mathbb{R}^d$. For any full-rank symmetric matrix $\mathbf{S} \in \mathbb{R}^{d \times d}$, let $\mathbf{S}_{F,G}$ be the $2d \times 2d$ matrix:

$$ \mathbf{S}_{F,G} := (\mathbf{\Pi}_F, \mathbf{\Pi}_G)^T \mathbf{S} (\mathbf{\Pi}_F, \mathbf{\Pi}_G) = \begin{bmatrix} \mathbf{\Pi}_F \mathbf{S} \mathbf{\Pi}_F & \mathbf{\Pi}_F \mathbf{S} \mathbf{\Pi}_G \\ \mathbf{\Pi}_G \mathbf{S} \mathbf{\Pi}_F & \mathbf{\Pi}_G \mathbf{S} \mathbf{\Pi}_G \end{bmatrix}. \quad (53) $$

Then, the range of the linear mapping associated to $\mathbf{S}_{F,G}$ is given by:

$$ \mathrm{Range}(\mathbf{S}_{F,G}) = \tilde{F} \times \tilde{G}, \quad \text{with} \quad \tilde{F} = F \cap \mathbf{S}(F \oplus G) \quad \text{and} \quad \tilde{G} = G \cap \mathbf{S}(F \oplus G). \quad (54) $$

Moreover, the restriction of $\mathbf{S}_{F,G}$ to its range is a one-to-one mapping whose inverse is given by:

$$ \mathbf{S}_{F,G}^{-1}(u, v)^T = (\mathbf{\Pi}_{\tilde{F}} \circ \mathbf{S}^{-1}(u + v), \mathbf{\Pi}_{\tilde{G}} \circ \mathbf{S}^{-1}(u + v))^T, \quad \forall (u, v) \in \mathrm{Range}(\mathbf{S}_{F,G}). \quad (55) $$

**Proof of Lemma A.4.** We show (54) using double inclusion. The following reasoning shows that $\mathrm{Range}(\mathbf{S}_{F,G}) \subset \tilde{F} \times \tilde{G}$.

$$ \begin{aligned} \mathbf{S}_{F,G}(\mathbb{R}^d \times \mathbb{R}^d) &\subset (\mathbf{\Pi}_F \mathbf{S} \mathbf{\Pi}_F)(\mathbb{R}^d) + (\mathbf{\Pi}_F \mathbf{S} \mathbf{\Pi}_G)(\mathbb{R}^d) \times (\mathbf{\Pi}_G \mathbf{S} \mathbf{\Pi}_F)(\mathbb{R}^d) + (\mathbf{\Pi}_G \mathbf{S} \mathbf{\Pi}_G)(\mathbb{R}^d) \\ &\subset \mathbf{\Pi}_F \mathbf{S}(F \oplus G) \times \mathbf{\Pi}_G \mathbf{S}(F \oplus G) \\ &\subset F \cap \mathbf{S}(F \oplus G) \times G \cap \mathbf{S}(F \oplus G) = \tilde{F} \times \tilde{G}. \end{aligned} $$
We show now the reciprocal inclusion. Letting $(u, v) \in \tilde{F} \times \tilde{G}$, we have:

$$
\begin{align*}
\mathbf{S}_{F,G} (\mathbf{S}^{-1}(u+v), \mathbf{S}^{-1}(u+v))^T &= (\mathbf{\Pi}_F \mathbf{S} \mathbf{\Pi}_F \mathbf{S}^{-1}(u+v) + \mathbf{\Pi}_F \mathbf{S} \mathbf{\Pi}_G \mathbf{S}^{-1}(u+v), \\
&\qquad \mathbf{\Pi}_F \mathbf{S} \mathbf{\Pi}_F \mathbf{S}^{-1}(u+v) + \mathbf{\Pi}_F \mathbf{S} \mathbf{\Pi}_G \mathbf{S}^{-1}(u+v)) \\
(\text{As } F \perp G : \mathbf{\Pi}_F + \mathbf{\Pi}_G = \mathbf{\Pi}_{F \oplus G}) &= (\mathbf{\Pi}_F \mathbf{S} \mathbf{\Pi}_{F \oplus G} \mathbf{S}^{-1}(u+v), \mathbf{\Pi}_G \mathbf{S} \mathbf{\Pi}_{F \oplus G} \mathbf{S}^{-1}(u+v)) \\
(\text{As } u+v \in \mathbf{S}(F \oplus G)) &= (\mathbf{\Pi}_F \mathbf{S}(\mathbf{S}^{-1}(u+v)), \mathbf{\Pi}_G \mathbf{S}(\mathbf{S}^{-1}(u+v))) \\
&= (\mathbf{\Pi}_F(u+v), \mathbf{\Pi}_G(u+v)) \\
(\text{As } u \in F \text{ and } v \in G) &= (u, v).
\end{align*}
$$

Consequently, we have $\tilde{F} \times \tilde{G} \subset \text{Range}(\mathbf{S}_{F,G})$. We conclude by showing (55). Following from the fact that the range of the linear mapping associated to any symmetric matrix is orthogonal to its kernel, we have that $\mathbf{S}_{F,G} = \mathbf{S}_{F,G} \circ \mathbf{\Pi}_{\text{Range}(\mathbf{S}_{F,G})} = \mathbf{S}_{F,G} \circ \mathbf{\Pi}_{\tilde{F} \times \tilde{G}}$. This, together with the fact that $\mathbf{\Pi}_{\tilde{F} \times \tilde{G}} = (\mathbf{\Pi}_{\tilde{F}}, \mathbf{\Pi}_{\tilde{G}})$, yields:

$$
\mathbf{S}_{F,G} (\mathbf{S}^{-1}(u+v), \mathbf{S}^{-1}(u+v))^T = \mathbf{S}_{F,G} (\mathbf{\Pi}_{\tilde{F}} \circ \mathbf{S}^{-1}(u+v), \mathbf{\Pi}_{\tilde{G}} \circ \mathbf{S}^{-1}(u+v))^T,
$$

for all $(u, v) \in \text{Range}(\mathbf{S}_{F,G})$, which concludes the proof. $\square$

**Lemma A.5.** Let $F, F', G, G'$ be subspaces of $\mathbb{R}^d$ such that $F' \subset F$, $G' \subset G$ and $F \perp G$. For any symmetric matrix $\mathbf{S} \in \mathbb{R}^{d \times d}$, let $\mathbf{S}_{F,G}$ be the one defined in (53) and let $\mathbf{S}_{F',G'}$ be defined analogously. Then, the following inclusions hold:

$$
(i) \operatorname{Range}(\mathbf{S}_{F', G'}) \subset \operatorname{Range}(\mathbf{S}_{F, G}),
$$

$$
(ii) \quad \mathbf{S}_{F', G'}^{\dagger}(0_d \times \mathbb{R}^d) \subset \mathbf{S}_{F, G}^{\dagger}(0_d \times \mathbb{R}^d),
$$

$$
(iii) \quad \mathbf{S}_{F', G'}^{\dagger}(\mathbb{R}^d \times 0_d) \subset \mathbf{S}_{F, G}^{\dagger}(\mathbb{R}^d \times 0_d),
$$

where $\mathbf{A}^{\dagger}(\cdot)$ is the linear mapping associated to the Moore-Penrose pseudo-inverse of a matrix $\mathbf{A}$.

*Proof of Lemma A.5.* We start by showing (i). As, by hypothesis, we have:

$$
F' \cap \mathbf{S}(F' \oplus G') \subset F \cap \mathbf{S}(F \oplus G) \quad \text{and} \quad G' \cap \mathbf{S}(F' \oplus G') \subset G \cap \mathbf{S}(F \oplus G),
$$

Equation (54) yields $\text{Range}(\mathbf{S}_{F',G'}) \subset \text{Range}(\mathbf{S}_{F,G})$. Let us show (ii). In the following, we will write $\text{Range}(\mathbf{S}_{F',G'}) = \tilde{F}' \times \tilde{G}'$ and $\text{Range}(\mathbf{S}_{F,G}) = \tilde{F} \times \tilde{G}$, as in (54). Inclusion (i) implies:

$$
\tilde{F}' \subset \tilde{F} \quad \text{and} \quad \tilde{G}' \subset \tilde{G}. \tag{56}
$$

As the pseudo-inverse of a symmetric matrix can be written as the composition of the orthogonal pro-
jection onto its range with its inverse on its range, Equation (55) yields:

$$
\begin{align*}
\mathbf{S}_{F', G'}^{\dagger} (0_d \times \mathbb{R}^d) &= \mathbf{S}_{F', G'}^{-1} \circ \mathbf{\Pi}_{\tilde{F}' \times \tilde{G}'} (0_d \times \mathbb{R}^d) \\
&= \mathbf{S}_{F, G'}^{-1} (0_d \times \tilde{G}') \\
&= (\mathbf{\Pi}_{\tilde{F}'}, \mathbf{\Pi}_{\tilde{G}'}) (\mathbf{S}^{-1} \tilde{G}').
\end{align*}
$$

Following the same reasoning we can show that $\mathbf{S}_{F,G}^{\dagger}(0_d \times \mathbb{R}^d) = (\mathbf{\Pi}_{\tilde{F}}, \mathbf{\Pi}_{\tilde{G}})(\mathbf{S}^{-1}\tilde{G})$, so if we prove that

$$
(\Pi_{\tilde{F}'}, \Pi_{\tilde{G}'}) (\mathbf{S}^{-1} \tilde{G}') \subset (\Pi_{\tilde{F}}, \Pi_{\tilde{G}}) (\mathbf{S}^{-1} \tilde{G}), \quad (57)
$$
inclusion (ii) will follow. Let $(h_{\tilde{F}'}, h_{\tilde{G}'}) \in (\Pi_{\tilde{F}'}, \Pi_{\tilde{G}'}) (\mathbf{S}^{-1}\tilde{G}')$ and let $h = h_{\tilde{F}'} + h_{\tilde{G}'}$. Thus, $h \in (\tilde{F}' \oplus \tilde{G}') \cap (\mathbf{S}^{-1}\tilde{G}') \subset (\tilde{F} \oplus \tilde{G}) \cap (\mathbf{S}^{-1}\tilde{G})$. From the unicity of the decomposition in $\tilde{F} \oplus \tilde{G}$ and (56), we have $(\Pi_{\tilde{F}}, \Pi_{\tilde{G}})(h) = (h_{\tilde{F}'}, h_{\tilde{G}'})$, which yields (57). The reasoning to prove (iii) is identical. $\square$

We are now ready to prove Theorem 2.4 and Proposition 2.5. Throughout the following proofs we will manipulate two intrinsically similar vector spaces, $\mathbb{R}^{n \times p}$ and $\mathbb{R}^{np \times 1}$, that are identified through the isometry $\text{vec} : \mathbb{R}^{n \times p} \to \mathbb{R}^{np \times 1}$. For the sake of a simpler notation, we will write $\tilde{F} = \text{vec}(F)$ for any vector space $F \subset \mathbb{R}^{n \times p}$. Then, the orthogonal projections onto $\tilde{F}$ and $F$ are identified via the equality $\Pi_{\tilde{F}} = \mathbf{I}_p \otimes \Pi_F$.

**Proof of Theorem 2.4.** We prove Theorem 2.4 in two steps. First, we show that the conditioned vector (18) has a p-dimensional normal distribution under (H0), explicitly deriving its mean and covariance matrix. Then, we will show that such distribution is centered. To shed light on the objects introduced in this proof, we keep the notation of the proof of Proposition 2.1. In particular, we denote by $\Lambda \subset \mathbb{R}^{n \times p}$ the kernel of the linear mapping $\nu^T : \mathbf{M} \in \mathbb{R}^{n \times p} \mapsto \nu^T \mathbf{M}$ and by $\Lambda^\perp$ its orthogonal complement. This means that $\Pi_\Lambda = \pi_\nu^\perp$ and $\Pi_{\Lambda^\perp} = \pi_\nu$, respectively. The idea is to find a matrix $\mathbf{A}_x$ and a vector $y_x$ such that the conditioned vector

$$ \text{vec}(\mathbf{X}) | \{ \Pi_\Lambda \mathbf{X} = \Pi_\Lambda \mathbf{x}, \text{dir}(\nu^T \mathbf{X}) = \pm \text{dir}(\nu^T \mathbf{x}) \} \quad (58) $$

can be rewritten as $\text{vec}(\mathbf{X}) | \{ \mathbf{A}_x \text{vec}(\mathbf{X}) = y_x \}$. Then, applying Theorem A.3 would yield an explicit Gaussian distribution for

$$ (\mathbf{I}_p \otimes \nu) \text{vec}(\mathbf{X}) | \{ \mathbf{A}_x \text{vec}(\mathbf{X}) = y_x \} = \nu^T \mathbf{X} | \{ \mathbf{A}_x \text{vec}(\mathbf{X}) = y \} = \bar{\mathbf{X}}_\nu(\mathbf{x}). \quad (59) $$

We start by rewriting the condition $\text{dir}(\nu^T \mathbf{X}) = \pm \text{dir}(\nu^T \mathbf{x})$ as follows. First, we have:

$$ \begin{align*} \text{dir}(\nu^T \mathbf{X}) = \pm \text{dir}(\nu^T \mathbf{x}) &\iff \nu^T \mathbf{X} \in \text{span}(\nu^T \mathbf{x}) \\ &\iff \mathbf{X} \in V_x := (\Lambda \oplus \text{span}(\mathbf{x})) \\ &\iff \text{vec}(\mathbf{X}) \in \tilde{V}_x = \tilde{\Lambda} \oplus \text{span}(\text{vec}(\mathbf{x})) \\ &\iff \Pi_{\tilde{V}_x^\perp} \text{vec}(\mathbf{X}) = 0. \end{align*} $$

Writing $\mathbf{x}_\nu := \text{vec}(\Pi_{\Lambda^\perp} \mathbf{x})$, we have that

$$ \tilde{V}_x^\perp = (\tilde{\Lambda} \oplus \text{span}(\text{vec}(\mathbf{x})))^\perp = (\tilde{\Lambda} \oplus \text{span}(\mathbf{x}_\nu))^\perp = \tilde{\Lambda}^\perp \cap \mathbf{x}_\nu^\perp, $$

where $\mathbf{x}_\nu^\perp$ denotes the orthogonal complement of $\mathbf{x}_\nu$. Since $\mathbf{x}_\nu \in \tilde{\Lambda}^\perp$, we can write $\Pi_{\tilde{V}_x^\perp} = \Pi_{\mathbf{x}_\nu^\perp} \circ \Pi_{\tilde{\Lambda}^\perp}$. This yields:

$$ \text{dir}(\nu^T \mathbf{X}) = \pm \text{dir}(\nu^T \mathbf{x}) \iff (\Pi_{\mathbf{x}_\nu^\perp} \circ \Pi_{\tilde{\Lambda}^\perp}) \text{vec}(\mathbf{X}) = \Pi_{\mathbf{x}_\nu^\perp} (\mathbf{I}_p \otimes \Pi_{\Lambda^\perp}) \text{vec}(\mathbf{X}) = 0. \quad (60) $$

Finally, using that

$$ \Pi_\Lambda \mathbf{X} = \Pi_\Lambda \mathbf{x} \iff (\mathbf{I}_p \otimes \Pi_\Lambda) \text{vec}(\mathbf{X}) = (\mathbf{I}_p \otimes \Pi_\Lambda) \text{vec}(\mathbf{x}), \quad (61) $$

we can characterize the conditioning set in (58) as follows:

$$ \text{dir}(\nu \mathbf{X}) = \pm \text{dir}(\nu \mathbf{x}) \text{ and } \Pi_\Lambda \mathbf{X} = \Pi_\Lambda \mathbf{x} \iff \mathbf{A}_x \text{vec}(\mathbf{x}) = y_x, \quad (62) $$
where $\mathbf{A}_x$ and $y_x$ are defined as

$$ \mathbf{A}_x = \begin{bmatrix} \Pi_{x_\nu}^\perp (\mathbf{I}_p \otimes \Pi_{\Lambda^\perp}) \\ \mathbf{I}_p \otimes \Pi_\Lambda \end{bmatrix}, \quad y_x = \begin{bmatrix} 0_{np} \\ (\mathbf{I}_p \otimes \Pi_\Lambda) \mathrm{vec}(\mathbf{x}) \end{bmatrix}, \qquad (63) $$

and $\Pi_{x_\nu}^\perp$ corresponds to the object $\pi_{x_\nu}^\perp$ defined in Theorem 2.4. Finally, using Theorem A.3 and the properties of the multivariate Gaussian distribution, we have that

$$ \bar{\mathbf{X}}_\nu(\mathbf{x}) \sim \mathcal{N}_p(\bar{\mu}_\nu(\mathbf{x}), \Gamma_\mathbf{x}), $$

where

$$ \bar{\mu}_\nu(\mathbf{x}) = (\mathbf{I}_p \otimes \nu^T)(\mathrm{vec}(\boldsymbol{\mu}) + (\boldsymbol{\Sigma} \otimes \mathbf{U})\mathbf{A}_{\mathbf{x}}^T(\mathbf{A}_{\mathbf{x}}(\boldsymbol{\Sigma} \otimes \mathbf{U})\mathbf{A}_{\mathbf{x}}^T)^\dagger(y_{\mathbf{x}} - \mathbf{A}_{\mathbf{x}} \mathrm{vec}(\boldsymbol{\mu}))), \quad (64) $$

and $\Gamma_\mathbf{x}$ is defined in (20).

We conclude by showing that $\bar{\mu}_\nu(\mathbf{x}) = 0_p$ under (H0) for all $\mathbf{x} \in \mathbb{R}^{n \times p}$. In what follows, we assume that (H0) holds. First, note that (H0) implies

$$ (\mathbf{I}_p \otimes \nu^T) \mathrm{vec}(\boldsymbol{\mu}) = 0_{np} \quad \text{and} \quad \mathbf{A}_{\mathbf{x}} \mathrm{vec}(\boldsymbol{\mu}) = (0_{np}, \mathrm{vec}(\boldsymbol{\mu}))^T, \qquad (65) $$

yielding

$$ y_{\mathbf{x}} - \mathbf{A}_{\mathbf{x}} \mathrm{vec}(\boldsymbol{\mu}) = (0_{np}, (\mathbf{I}_p \otimes \Pi_{\Lambda}) \mathrm{vec}(\mathbf{x}) - \mathrm{vec}(\boldsymbol{\mu})). \qquad (66) $$

Consequently, proving $\bar{\mu}_\nu(\mathbf{x}) = 0$ comes down to show that $0_{np} \times \mathbb{R}^{np}$ is included in the kernel of the linear operator defined by the matrix

$$ (\mathbf{I}_p \otimes \nu^T)(\boldsymbol{\Sigma} \otimes \mathbf{U})\mathbf{A}_{\mathbf{x}}^T(\mathbf{A}_{\mathbf{x}}(\boldsymbol{\Sigma} \otimes \mathbf{U})\mathbf{A}_{\mathbf{x}}^T)^\dagger, \qquad (67) $$

or, equivalently, in the kernel of the linear operator associated to

$$ \Pi_{\tilde{\Lambda} \perp} (\boldsymbol{\Sigma} \otimes \mathbf{U}) \mathbf{A}_{\mathbf{x}}^T (\mathbf{A}_{\mathbf{x}} (\boldsymbol{\Sigma} \otimes \mathbf{U}) \mathbf{A}_{\mathbf{x}}^T)^\dagger. \qquad (68) $$

Let us consider the matrix $\mathbf{A} = (\Pi_{\tilde{\Lambda} \perp}, \Pi_{\tilde{\Lambda}})^T$. Then, if the following statements hold:

$$ (S1) \quad \Pi_{\tilde{\Lambda} \perp} (\boldsymbol{\Sigma} \otimes \mathbf{U}) \mathbf{A}^T (\mathbf{A} (\boldsymbol{\Sigma} \otimes \mathbf{U}) \mathbf{A}^T)^\dagger (0_{np} \times \mathbb{R}^{np}) = 0_{np}, $$

$$ (S2) \quad \Pi_{\tilde{\Lambda} \perp} (\boldsymbol{\Sigma} \otimes \mathbf{U}) \mathbf{A}_{\mathbf{x}}^T (\mathbf{A} (\boldsymbol{\Sigma} \otimes \mathbf{U}) \mathbf{A}_{\mathbf{x}}^T)^\dagger (0_{np} \times \mathbb{R}^{np}) \subset \Pi_{\tilde{\Lambda} \perp} (\boldsymbol{\Sigma} \otimes \mathbf{U}) \mathbf{A}^T (\mathbf{A}_{\mathbf{x}} (\boldsymbol{\Sigma} \otimes \mathbf{U}) \mathbf{A}^T)^\dagger (0_{np} \times \mathbb{R}^{np}), $$

the subspace $0_{np} \times \mathbb{R}^{np}$ is included in the kernel of (68) and the result follows.

Since $\Pi_{\tilde{\Lambda} \perp}$ is a sub-block of $\mathbf{A}$, (S1) is equivalent to the equality:

$$ \mathbf{A} (\boldsymbol{\Sigma} \otimes \mathbf{U}) \mathbf{A}^T (\mathbf{A} (\boldsymbol{\Sigma} \otimes \mathbf{U}) \mathbf{A}^T)^\dagger (0_{np} \times \mathbb{R}^{np}) = (0_{np}, V)^T, \qquad (69) $$

for a subspace $V \subset \mathbb{R}^{np}$. From the properties of the Moore-Penrose pseudo-inverse, we have that

$$ \mathbf{A} (\boldsymbol{\Sigma} \otimes \mathbf{U}) \mathbf{A}^T (\mathbf{A} (\boldsymbol{\Sigma} \otimes \mathbf{U}) \mathbf{A}^T)^\dagger = \Pi_{\mathrm{Range}(\mathbf{A} (\boldsymbol{\Sigma} \otimes \mathbf{U}) \mathbf{A}^T)}. $$

This, together with Lemma A.4, yields (69).

To prove (S2), it suffices to show that:

$$ \mathbf{A}_{\mathbf{x}}^T (\mathbf{A}_{\mathbf{x}} (\boldsymbol{\Sigma} \otimes \mathbf{U}) \mathbf{A}_{\mathbf{x}}^T)^\dagger (0_{np} \times \mathbb{R}^{np}) \subset \mathbf{A}^T (\mathbf{A} (\boldsymbol{\Sigma} \otimes \mathbf{U}) \mathbf{A}^T)^\dagger (0_{np} \times \mathbb{R}^{np}). $$
Inclusion (ii) in Lemma A.5 yields:

$$ (\mathbf{A}_{\mathbf{x}}(\boldsymbol{\Sigma} \otimes \mathbf{U})\mathbf{A}_{\mathbf{x}}^{\mathrm{T}})^{\dagger}(0_{np} \times \mathbb{R}^{np}) \subset (\mathbf{A}(\boldsymbol{\Sigma} \otimes \mathbf{U})\mathbf{A}^{\mathrm{T}})^{\dagger}(0_{np} \times \mathbb{R}^{np}). $$

Finally, following the same strategy as in the proof of Lemma A.5, we can show that the previous inclusion is stable when composed by $\mathbf{A}_{\mathbf{x}}$ on the left side and $\mathbf{A}$ on the right side, which yields (S2). □

**Proof of Proposition 2.5.** We keep the notation of the proof of Theorem 2.4. To show (23), the key idea is to prove that the rank of $\mathbf{\Gamma}_{\mathbf{x}}$ is $\mathbf{x}$-a.s. constant equal to one. From the condition $\mathrm{dir}(\nu^{\mathrm{T}}\mathbf{X}) = \pm\mathrm{dir}(\nu^{\mathrm{T}}\mathbf{x})$, we have clearly that the rank of $\mathbf{\Gamma}_{\mathbf{x}}$ is upper bounded by one. Moreover, since the mapping $\nu: \mathbb{R}^{1 \times p} \to \mathbb{R}^{n \times p}$ defined by $z \mapsto \nu z$ is injective, the rank of the covariance matrix of $\bar{\mathbf{X}}_{\nu}(\mathbf{x})$ is the same as the rank of $\nu \bar{\mathbf{X}}_{\nu}(\mathbf{x})$ and, from the proof of Theorem 2.4, the same as the rank of the matrix

$$ \mathbf{\Pi}_{\Lambda \perp} (\mathbf{X} | \{\mathbf{\Pi}_{\Lambda} \mathbf{X} = \mathbf{\Pi}_{\Lambda} \mathbf{x}, \mathbf{\Pi}_{\mathbf{x}_{\nu}^{\perp}} \mathbf{\Pi}_{\Lambda^{\perp}} \mathbf{X} = \mathbf{0}\}). $$

Following the steps of the proof of Theorem A.3 (Proposition 3.13 in [19]), we can decompose $\mathbf{\Pi}_{\Lambda \perp} \mathbf{X}$ as the sum of two independent Gaussian vectors $\mathbf{Y}$ and $\mathbf{Z}$, with with $\mathbf{Y} = \mathbf{\Pi}_{\Lambda} \mathbf{X} + \mathbf{\Pi}_{\mathbf{x}_{\nu}^{\perp}} \mathbf{\Pi}_{\Lambda^{\perp}} \mathbf{X}$. Thus, $\mathbf{Z}$ must be non-zero since otherwise $\mathbf{\Pi}_{\mathbf{x}_{\nu}} \mathbf{X} = 0$, and $\mathbf{X}$ is non degenerated. As $\mathbf{\Gamma}_{\mathbf{x}}$ is the covariance matrix of $\mathbf{Z}$, its rank is $\mathbf{x}$-a.s. equal to one. This implies that $||\bar{\mathbf{X}}_{\nu}(\mathbf{x})||_{\mathbf{\Gamma}_{\mathbf{x}}} \sim \chi_1$ $\mathbf{x}$-a.s. under (H0), where $|| \cdot ||_{\mathbf{\Gamma}_{\mathbf{x}}}$ is defined in (22). Following the same steps as in the proofs of Theorem 2.2 and Lemma 2.3, we have (23) and (24). □

## A.2 Proofs of Section 3

**Proof of Theorem 3.1.** We follows the steps of the proof of Theorem 4 in [21]. For simplicity, we use $\hat{p}_n$ to denote $p_{\hat{\mathbf{V}}_{\mathcal{G}_1^{(n)}, \mathcal{G}_2^{(n)}}}(\mathbf{X}^{(n)}; \{\mathcal{G}_1^{(n)}, \mathcal{G}_2^{(n)}\})$, $p_n$ to denote $p_{\mathbf{V}_{\mathcal{G}_1^{(n)}, \mathcal{G}_2^{(n)}}}(\mathbf{X}^{(n)}; \{\mathcal{G}_1^{(n)}, \mathcal{G}_2^{(n)}\})$, $\hat{\mathbf{V}}_n$ to denote $\hat{\mathbf{V}}_{\mathcal{G}_1^{(n)}, \mathcal{G}_2^{(n)}}$, $\mathbf{V}_n$ to denote $\mathbf{V}_{\mathcal{G}_1^{(n)}, \mathcal{G}_2^{(n)}}$ and $\nu_n$ to denote $\nu_{\mathcal{G}_1^{(n)}, \mathcal{G}_2^{(n)}}$. We will also write the difference of cluster means as the row vector $\nu_n^{\mathrm{T}}\mathbf{X}^{(n)}$ for the sake of a clearer notation. If we show that

$$ \hat{\mathbf{\Sigma}}(\mathbf{X}^{(n)}) \succeq \mathbf{\Sigma} \Rightarrow \hat{p}_n \geq p_n, \qquad (70) $$

then the result follows using the same reasoning as in the proof of [21, Theorem 4], replacing the usual order $\geq$ in $\mathbb{R}$ by the Loewner partial order $\succeq$ between matrices. Consequently, we only need to prove (70). First note that, as the Kronecker product is distributive, we have

$$ \hat{\mathbf{\Sigma}}(\mathbf{X}^{(n)}) \succeq \mathbf{\Sigma} \Rightarrow \hat{\mathbf{V}}_n \succeq \mathbf{V}_n. \qquad (71) $$

Next, by Corollary 7.7.4(a) and Theorem 7.7.2(a) in [26], we can write

$$ \begin{aligned} & \hat{\mathbf{V}}_n \succeq \mathbf{V}_n \Leftrightarrow \mathbf{V}_n^{-1} \succeq \hat{\mathbf{V}}_n^{-1} \Rightarrow (\nu_n^{\mathrm{T}}\mathbf{X}^{(n)}) \mathbf{V}_n^{-1} (\nu_n^{\mathrm{T}}\mathbf{X}^{(n)})^T \\ & \geq (\nu_n^{\mathrm{T}}\mathbf{X}^{(n)}) \hat{\mathbf{V}}_n^{-1} (\nu_n^{\mathrm{T}}\mathbf{X}^{(n)})^T \Leftrightarrow ||\nu_n^{\mathrm{T}}\mathbf{X}^{(n)}||_{\mathbf{V}_n} \geq ||\nu_n^{\mathrm{T}}\mathbf{X}^{(n)}||_{\hat{\mathbf{V}}_n}. \end{aligned} \qquad (72) $$

Let us then state that, if $\mathbb{F}_p(t, c, S)$ denotes the cumulative distribution function of a $c \cdot \chi_p$ distribution truncated to the set $S$, for $c > 0$, it follows that

$$ \mathbb{F}_p(t, c, aS) = \mathbb{F}_p\left(\frac{t}{a}, \frac{c}{a}, S\right), \qquad (73) $$
for any $a > 0$. We prove (73) as a technical lemma after the proof. Consequently, taking

$$a = \frac{\|\nu_n^T \mathbf{X}^{(n)}\|_{\hat{\mathbf{V}}_n}}{\|\nu_n^T \mathbf{X}^{(n)}\|_{\mathbf{V}_n}} \le 1, \qquad (74)$$

we have

$$\begin{aligned} 1 - \hat{p}_n &= \mathbb{F}_p (\|\nu_n^T \mathbf{X}^{(n)}\|_{\hat{\mathbf{V}}_n}, \mathcal{S}_{\hat{\mathbf{V}}_n}) = \mathbb{F}_p (\|\nu_n^T \mathbf{X}^{(n)}\|_{\hat{\mathbf{V}}_n}, a \mathcal{S}_{\hat{\mathbf{V}}_n}) \\ &= \mathbb{F}_p \left( \frac{1}{a} \|\nu_n^T \mathbf{X}^{(n)}\|_{\hat{\mathbf{V}}_n}, \frac{1}{a}, \mathcal{S}_{\hat{\mathbf{V}}_n} \right) = \mathbb{F}_p \left( \|\nu_n^T \mathbf{X}^{(n)}\|_{\mathbf{V}_n}, \frac{1}{a}, \mathcal{S}_{\mathbf{V}_n} \right) \\ &\le \mathbb{F}_p (\|\nu_n^T \mathbf{X}^{(n)}\|_{\mathbf{V}_n}, 1, \mathcal{S}_{\mathbf{V}_n}) = 1 - p_n, \end{aligned} \qquad (75)$$

where the last inequality follows from Lemma A.3 in [21]. This shows (70). □

**Lemma A.6.** For $c > 0$ and $\emptyset \neq S \subset \mathbb{R}$, let $\mathbb{F}_p(t, c, S)$ denote the cumulative distribution function of a $c \cdot \chi_p$ distribution truncated to $S$. Then, for any $a > 0$, it holds

$$\mathbb{F}_p(t, c, a S) = \mathbb{F}_p \left( \frac{t}{a}, \frac{c}{a}, S \right).$$

*Proof of Lemma A.6.* First, if we denote by $f(t, c, S)$ the probability density function of a $c \cdot \chi_p$ distribution truncated to the set $S$, we have

$$f(t, c, a S) = \frac{1}{a} f \left( \frac{t}{a}, \frac{c}{a}, S \right). \qquad (76)$$

Indeed, following the first lines of the proof of [21, Lemma A.3], we can rewrite $f(t, c, a S)$ as

$$f(t, c, a S) = \frac{t^{p-1} \mathbb{1}\{t \in a S\}}{\int u^{p-1} \exp(-\frac{u^2}{2c^2}) \mathbb{1}\{t \in a S\} du} \exp\left(-\frac{t^2}{2c^2}\right), \qquad (77)$$

that we can easily express in terms of $t/a$ as

$$f(t, c, a S) = \frac{\left(\frac{t}{a}\right)^{p-1} \mathbb{1}\{\frac{t}{a} \in S\}}{\int \left(\frac{u}{a}\right)^{p-1} \exp\left(-\frac{\left(\frac{u}{a}\right)^2}{2(c/a)^2}\right) \mathbb{1}\{\frac{t}{a} \in S\} du} \exp\left(-\frac{(t/a)^2}{2(c/a)^2}\right) = \frac{1}{a} f\left(\frac{t}{a}, \frac{c}{a}, S\right), \qquad (78)$$

where the last equality follows from taking the variable change $y = u/a$ in the integral. Finally, we have

$$\mathbb{F}_p(t, c, a S) = \int_0^t f(x, c, a S) dx = \frac{1}{a} \int_0^t f\left(\frac{x}{a}, \frac{c}{a}, S\right) dx = \int_0^t f\left(u, \frac{c}{a}, S\right) du = \mathbb{F}_p\left(\frac{t}{a}, \frac{c}{a}, S\right),$$

which concludes the proof. □

*Alternative formulation of Assumption 3.2.* Assumption 3.2 can be formulated in terms of strong mixing of measure-preserving dynamical systems [29, Chapter 20]. To show this, let us consider the sets $A_k = \{i \in \mathbb{N} : \mu_i^{(i)} = \theta_k\}$ for any $k = 1, \dots, K^*$. This makes the family $\mathcal{F} = \mathcal{P}(A)$ with $\mathcal{A} = \{A_k\}_{k=1}^{K^*}$ a $\sigma$-algebra on $\mathbb{N}$. Next, let $P_n$ denote the measure defined by $P_n(A) = \frac{1}{n}|A \cap [n]|$ for any $A \in \mathcal{F}$ and $P$ denote the measure defined by $P(\cup_{s \in S} A_s) = \sum_{s \in S} \pi_s$ for any $S \in \mathcal{P}(\{1, \dots, K^*\})$. Note that the pair $(\mathbb{N}, \mathcal{F})$ can be provided with either $P_n$ or $P$ to form a measure space. Besides, Assumption 3.1 states the setwise convergence of $P_n$ to $P$ when $n \to \infty$. Finally, for any $k = 1, \dots, K^*$, we can define the
transformation $T(A_k) = \{i \in \mathbb{N} : \mu_{i-1}^{(i-1)} = \theta_k\}$, which is measure-preserving on $(\mathbb{N}, \mathcal{F}, P)$ [29, Definition 20.6]. Then, Equation (30) can be rewritten as:

$$P_n(T^{-r}(A_k) \cap A_{k'}) \underset{n \to \infty}{\longrightarrow} P(T^{-r}(A_k) \cap A_{k'}) \underset{r \to \infty}{\longrightarrow} P(A_k) P(A_{k'}), \quad (79)$$

for any $k, k' \in \{1, \dots, K^*\}$. The first limit in (79) follows from Assumption 3.1, whereas the second one is equivalent to state that the measure-preserving dynamical system $(\mathbb{N}, \mathcal{F}, P, T)$ is (strong) mixing (see [29, Definition 20.04]). □

**Proof of Remark 3.1.** Let $U^{(n)} = b\mathbf{1}_{n \times n} + (a-b)\mathbf{I}_n$. As $b \in (-\frac{a}{n-1}, a)$ if and only if $U^{(n)}$ is positive definite, condition $0 < b < a$ is needed to ensure positive definiteness for all $n \in \mathbb{N}$. Following the Sherman-Morrison formula [4], we can derive an explicit expression for the sequence of inverse matrices:

$$\left(U^{(n)}\right)^{-1} = \frac{1}{a-b} \mathbf{I}_n + \frac{-b}{(a-b)(nb+a-b)}, \quad \forall n \in \mathbb{N}. \quad (80)$$

Consequently, for every $r \ge 0$ and every $i \in \mathbb{N}$, we have

$$\left(U^{(n)}\right)^{-1}_{i i+r} = \begin{cases} \frac{1}{a-b} + \frac{-b}{(a-b)(nb+a-b)} & \text{if } r=0, \\ \frac{-b}{(a-b)(nb+a-b)} & \text{if } r>0, \end{cases}$$

which are monotone, so condition (ii) in Assumption 3.3 is satisfied. Then, we have

$$\Lambda_{ii+r} = \begin{cases} \frac{1}{a-b} & \text{if } r=0, \\ 0 & \text{if } r>0, \end{cases}$$

for all $i \in \mathbb{N}$, $\lambda_0 = 1/(a-b)$ and $\lambda_r = 0$ for $r > 0$. Consequently, Assumption 3.3 holds. □

**Proof of Remark 3.2.** The case of diagonal matrices is straightforward as both $U^{(n)}$ and $(U^{(n)})^{-1}$ are defined by a sequence $\{\lambda_i\}_{i \in \mathbb{N}}$. Every diagonal entry of the inverse satisfies $(U^{(n)})^{-1}_{ii} = \frac{1}{\lambda_i}$ for all $n \in \mathbb{N}$ and, as we asked the $\lambda_i$ to converge to $\lambda$, which is strictly positive due to the positive definiteness of $U^{(n)}$, Assumption 3.3 is satisfied. □

**Proof of Remark 3.3.** The inverse of an auto-regressive covariance matrix of lag $P \ge 1$ is banded with $2P-1$ non-zero diagonals. Its explicit form is derived in [47] for a stationary process of any lag, and the cases $P \le 3$ are discussed in detail in [48]. From these results we can derive the behavior of the sequences $\{(U^{(n)})^{-1}_{ii+r}\}$ as $n$ increases. The diagonal elements define the sequences

$$\sigma^2 \left\{ \left( U^{(n)} \right)^{-1}_{ii} \right\}_{n \in \mathbb{N}} = \begin{cases} \{1 + \sum_{k=1}^{i-1} \beta_k^2, 1 + \sum_{k=1}^{i-1} \beta_k^2, \dots \} & \text{if } i \le P+1, \\ \{0, \dots, 0, 1, 1 + \beta_1^2, 1, 1 + \beta_1^2 \beta_2^2, \dots, 1 + \sum_{k=1}^P \beta_k^2, 1 + \sum_{k=1}^P \beta_k^2, \dots \} & \text{if } i > P+1, \end{cases}$$

where the sums are taken as zero if the upper limit of summation is zero. Note that these sequences do not satisfy condition (i) in Assumption 3.3 as, even if each sequence reaches its limit after a finite number of terms, the index of the term where the limit is reached diverges with $i$. In other words, we can dominate the sequence, but not by a summable one. However, for all $i \in \mathbb{N}$ the series are non-decreasing
so condition (ii) is satisfied and we have

$$ \sigma^2 \Lambda_{ii} = \begin{cases} 1 + \sum_{k=1}^{i-1} \beta_k^2 & \text{if } i \le P+1 \\ 1 + \sum_{k=1}^{P} \beta_k^2 & \text{if } i > P+1. \end{cases} $$

Then, $\sigma^2 \lambda_0 = 1 + \sum_{k=1}^{P} \beta_k^2$. The sequences outside the main diagonal show a similar behavior, but they are not positive in general. As, following the same reasoning, they do not satisfy condition (i) in Assumption 3.3, we force them to satisfy condition (ii). For any $0 < r \le P$, we have

$$ \sigma^2 \left\{ \left( U^{(n)} \right)_{ii+r}^{-1} \right\}_{n \in \mathbb{N}} = \begin{cases} \{-\beta_r + \sum_{k=1}^{i-(r+1)} \beta_k \beta_{k+r}, -\beta_r + \sum_{k=1}^{i-(r+1)} \beta_k \beta_{k+r}, \dots \} & \text{if } i \le P+1, \\ \{0, i \ne 1, 0, -\beta_r + \beta_1 \beta_{1+r}, -\beta_r + \beta_1 \beta_{1+r} + \beta_2 \beta_{2+r}, \dots, \\ -\beta_r + \sum_{k=1}^{P-r} \beta_k \beta_{k+r}, -\beta_r + \sum_{k=1}^{P-r} \beta_k \beta_{k+r}, \dots \} & \text{if } i > P+1. \end{cases} \quad (81) $$

For these sequences to satisfy condition (ii) we need them to be non-decreasing or non-increasing. For $P \le 2$ this is always satisfied but, for $P > 2$, we need to require all the $\beta_k$ to have the same sign. In that case, condition (ii) holds and we have

$$ \sigma^2 \Lambda_{ii+r} = \begin{cases} -\beta_r + \sum_{k=1}^{i-(r+1)} \beta_k \beta_{k+r} & \text{if } i \le P+1, \\ -\beta_r + \sum_{k=1}^{P-r} \beta_k \beta_{k+r} & \text{if } i > P+1, \end{cases} $$

and, consequently, $\sigma^2 \lambda_r = -\beta_r + \sum_{k=1}^{P-r} \beta_k \beta_{k+r}$. As the sequence $\{\lambda_r\}_{r=1}^\infty$ is non-zero for a finite number of terms (due to the bandedness of the inverse matrix), its sum converges and Assumption 3.3 is satisfied. $\square$

*Proof of Lemma 3.5.* We start by rewriting the sum in (36) as a sum along each diagonal. Using the symmetry of $(U^{(n)})^{-1}$ we have,

$$ \begin{aligned} & \lim_{n \to \infty} \frac{1}{n} \sum_{l,s=1}^{n} (U^{(n)})_{ls}^{-1} \mathbf{1}\{\mu_l^{(n)} = \theta_k\} \mathbf{1}\{\mu_s^{(n)} = \theta_{k'}\} \\ & = \lim_{n \to \infty} \frac{1}{n} \sum_{r=1}^{n-1} \sum_{i=1}^{n-r} (U^{(n)})_{ii+r}^{-1} \mathbf{1}\{\mu_i^{(n)} = \theta_k\} \mathbf{1}\{\mu_{i+r}^{(n)} = \theta_{k'}\} \end{aligned} \quad (82) $$

$$ + \lim_{n \to \infty} \frac{1}{n} \sum_{r=1}^{n-1} \sum_{i=1}^{n-r} (U^{(n)})_{ii+r}^{-1} \mathbf{1}\{\mu_{i+r}^{(n)} = \theta_k\} \mathbf{1}\{\mu_i^{(n)} = \theta_{k'}\} \quad (83) $$

$$ + \lim_{n \to \infty} \frac{1}{n} \sum_{i=1}^{n} (U^{(n)})_{ii}^{-1} \mathbf{1}\{\mu_i^{(n)} = \theta_k\} \mathbf{1}\{\mu_i^{(n)} = \theta_{k'}\}, \quad (84) $$

where (82),(83) and (84) are respectively the sums along all the superdiagonals, subdiagonals and along the main diagonal. Let us detail the general reasoning that we use to show that the three quantities converge. Let $\{a_i^{(n)}\}_{i \in \mathbb{N}}$ be a double sequence such that $\lim_{n \to \infty} a_i^{(n)} = a_i \in \mathbb{R}$, and let $\{b_i^{(n)}\}_{i \in \mathbb{N}}$ be a binary Cesàro summable double sequence, i.e. such that $\lim_{n \to \infty} \frac{1}{n} \sum_{i=1}^{n} b_i^{(n)} = b$ and $b_i^{(n)} \in \{0, 1\}$ for all $i, n \in \mathbb{N}$. Let us first show that, if $\{a_i^{(n)}\}_{n \in \mathbb{N}}$ satisfies any of the conditions (i) or (ii), and the sequence $\{a_i^{(1)} - a_i\}_{i=1}^\infty \in \ell_1$, we can write

$$ \lim_{n \to \infty} \frac{1}{n} \sum_{i=1}^{n} a_i^{(n)} b_i^{(n)} = \lim_{n \to \infty} \frac{1}{n} \sum_{i=1}^{n} a_i b_i^{(n)}. \quad (85) $$
First, note that

$$ \lim_{n \to \infty} \frac{1}{n} \sum_{i=1}^{n} a_i^{(n)} b_i^{(n)} = \lim_{n \to \infty} \frac{1}{n} \sum_{i=1}^{n} (a_i^{(n)} - a_i) b_i^{(n)} + \lim_{n \to \infty} \frac{1}{n} \sum_{i=1}^{n} a_i b_i^{(n)}. \quad (86) $$

Therefore, it suffices to show that the first term in (86) is zero to have (85). Using Hölder's inequality, we have

$$ \begin{aligned} \lim_{n \to \infty} \frac{1}{n} \left| \sum_{i=1}^{n} (a_i^{(n)} - a_i) b_i^{(n)} \right| & \le \lim_{n \to \infty} \frac{1}{n} \sum_{i=1}^{n} \left| (a_i^{(n)} - a_i) b_i^{(n)} \right| \\ & \le \lim_{n \to \infty} \left( \sum_{i=1}^{n} (a_i^{(n)} - a_i)^2 \right)^{\frac{1}{2}} \lim_{n \to \infty} \frac{1}{n} \left( \sum_{i=1}^{n} b_i^{(n)} \right)^{\frac{1}{2}}. \end{aligned} $$

On one hand,

$$ \lim_{n \to \infty} \frac{1}{n} \left( \sum_{i=1}^{n} b_i^{(n)} \right)^{\frac{1}{2}} = \lim_{n \to \infty} \frac{1}{\sqrt{n}} \lim_{n \to \infty} \left( \frac{1}{n} \sum_{i=1}^{n} b_i^{(n)} \right)^{\frac{1}{2}} = 0. $$

On the other hand, let us show that

$$ \lim_{n \to \infty} \sum_{i=1}^{n} (a_i^{(n)} - a_i)^2 = 0 \quad (87) $$

if {$a_i^{(n)}$}$_{n \in \mathbb{N}}$ satisfies any of the conditions (i) or (ii). If {$a_i^{(n)}$}$_{n \in \mathbb{N}}$ satisfies (i), the sequence {($a_i^{(n)}$ - $a_i$)$^2$}$_{n \in \mathbb{N}}$ is dominated by the sequence {$\alpha_i^2$}$_{i \in \mathbb{N}}$, which is summable as $\ell_1 \subset \ell_2$. Then, (85) holds following the Dominated Convergence Theorem [49, Theorem 9.20]. If {$a_i^{(n)}$}$_{n \in \mathbb{N}}$ is non-increasing, then $a_i^{(n+1)} - a_i \le a_i^{(n)} - a_i$ implies $(a_i^{(n+1)} - a_i)^2 \le (a_i^{(n)} - a_i)^2$ and $\tilde{a}_i^{(n)} := (a_i^{(n)} - a_i)^2$ is a non-increasing and non-negative sequence. Similarly, if {$a_i^{(n)}$}$_{n \in \mathbb{N}}$ is non-decreasing, then $a_i^{(n+1)} - a_i \ge a_i^{(n)} - a_i$ implies $(a_i^{(n+1)} - a_i)^2 \le (a_i^{(n)} - a_i)^2$ and $\tilde{a}_i^{(n)}$ is again a non-increasing and non-negative sequence. Then, the sequence $z_i^{(n)} := \tilde{a}_i^{(1)} - \tilde{a}_i^{(n)}$ is non-negative and non-decreasing. Thus, following the Monotone Convergence Theorem [49, Theorem 8.5], we have

$$ \lim_{n \to \infty} \sum_{i=1}^{n} z_i^{(n)} = \lim_{n \to \infty} \sum_{i=1}^{n} (a_i^{(1)} - a_i)^2, \quad (88) $$

which implies (87) if the limit in the right side of (88) exists and is finite. This is guaranteed if we ask the sequence {$a_i^{(1)} - a_i$}$_{i=1}^\infty$ to be summable. This always holds in our case as we can arbitrarily define the entries $(U^{(n)})_{ii+r}^{-1}$ for $i > n$. Consequently, if we write {${(U^{(1)})_{i+r}^{-1}}$}$_{i=1}^\infty = {((U^{(1)})_{1+r}^{-1}, \Lambda_{2+1+r}, \Lambda_{3+1+r}, \dots}$, the sequence {${(U^{(1)})_{i+r}^{-1} - \Lambda_{i+r}}$}$_{i=1}^\infty$ is trivially summable. This proves (85).

Now, if we have that $\lim_{i \to \infty} a_i = a$, let us show that

$$ \lim_{n \to \infty} \frac{1}{n} \sum_{i=1}^{n} a_i b_i^{(n)} = ab. \quad (89) $$

First, let separate the sum in (89) as

$$ \frac{1}{n} \sum_{i=1}^{n} a_i b_i^{(n)} = \frac{1}{n} \sum_{i=1}^{n} (a_i - a) b_i^{(n)} + \frac{a}{n} \sum_{i=1}^{n} b_i^{(n)}. \quad (90) $$
The right term tends to ab when $n \to \infty$. Let's show that the first term tends to zero. For any $i_0 \in \mathbb{N}$, we can write

$$ \left| \frac{1}{n} \sum_{i=1}^{n} (a_i - a) b_i^{(n)} \right| \le \left| \frac{1}{n} \sum_{i=1}^{i_0-1} (a_i - a) b_i^{(n)} \right| + \left| \frac{1}{n} \sum_{i=i_0}^{n} (a_i - a) b_i^{(n)} \right| \quad (91) $$

$$ \le \sup_{i<i_0} |a_i - a| \frac{1}{n} \sum_{i=1}^{i_0-1} b_i^{(n)} + \sup_{i \ge i_0} |a_i - a| \frac{1}{n} \sum_{i=i_0}^{n} b_i^{(n)} \le \frac{C}{n} + \sup_{i \ge i_0} |a_i - a| \frac{1}{n} \sum_{i=i_0}^{n} b_i^{(n)}, \quad (92) $$

where C is a real constant. Then, following the definition of limit, when can choose $i_0$ as the one such that for all $i \ge i_0$ we have $|a_i - a| \le \frac{1}{n}$. Therefore,

$$ \left| \frac{1}{n} \sum_{i=1}^{n} (a_i - a) b_i^{(n)} \right| \le \frac{C}{n} + \frac{1}{n^2} \sum_{i=i_0}^{n} b_i^{(n)}, \quad (93) $$

which tends to zero when $n \to \infty$ using that $\{b_i^{(n)}\}_{i \in \mathbb{N}}$ has Cesàro sum b. Thus, we have (89). As the sequences $(U^{(n)})_{i+r}^{-1}$ have limits $\Lambda_{i+r}$ when $i \to \infty$, following Assumption 3.2, and the products of indicator functions are Cesàro summable thanks to Assumptions 3.1 and 3.2, we can use (85) and (89) to rewrite the three limits in (82), (83), (84) as

$$ \begin{aligned} & \lim_{n \to \infty} \frac{1}{n} \sum_{l,s=1}^{n} (U^{(n)})_{ls}^{-1} \mathbf{1}\{\mu_l^{(n)} = \theta_k\} \mathbf{1}\{\mu_s^{(n)} = \theta_{k'}\} \\ & = \lim_{n \to \infty} \sum_{r=1}^{n-1} \lambda_r (\pi_{kk'}^r + \pi_{k'k}^r) + \lambda_0 \pi_k \delta_{kk'} = 2(\lambda - \lambda_0)\pi_k \pi_{k'} + \lambda_0 \pi_k \delta_{kk'}, \end{aligned} \quad (94) $$

where the last limit is derived following the same reasoning as to prove (89). This concludes the proof. □

*Proof of Proposition 3.4.* We start by proving the element-wise convergence in probability of (hat-Sigma). More precisely, we show that

$$ \hat{\Sigma}_{ij}^{(n)} \xrightarrow{p} \Sigma_{ij} + \lambda_0 \sum_{k=1}^{K^*} \pi_k (\theta_{ki} - \tilde{\theta}_i) (\theta_{kj} - \tilde{\theta}_j), \quad (95) $$

for all $i,j \in \{1, \dots, p\}$, where $\hat{\Sigma}_{ij}^{(n)}$ is the $ij$ entry of $\hat{\Sigma}(X^{(n)})$, that is,

$$ \hat{\Sigma}_{ij} = \frac{1}{n-1} \sum_{l,s=1}^{n} (X_{li} - \bar{X}_i) (U^{-1})_{ls} (X_{sj} - \bar{X}_j), \quad \forall i,j \in \{1, \dots, p\}, \quad (96) $$

where $\bar{X}_i = \frac{1}{n} \sum_{k=1}^{n} X_{ki}$, and we have defined $\tilde{\theta}_i = \sum_{k=1}^{K^*} \pi_k \theta_{ki}$. Recall that all the quantities in (95) have been defined in Assumptions 3.1 and 3.3. To prove (95), it suffices to show, following the same reasoning as in the proof of [21, Lemma C.1], that

$$ \lim_{n \to \infty} \mathbb{E} (\hat{\Sigma}_{ij}^{(n)}) = \Sigma_{ij} + \lambda_0 \sum_{k=1}^{K^*} \pi_k (\theta_{ki} - \tilde{\theta}_i) (\theta_{kj} - \tilde{\theta}_j) \quad \text{and} \quad \mathrm{Var} (\hat{\Sigma}_{ij}^{(n)}) = 0. \quad (97) $$
Indeed, (97) implies convergence in mean of $\hat{\Sigma}_{ij}^{(n)}$ towards the limit of its expectation and, following Markov's inequality, convergence in probability. Let start by rewriting $\hat{\Sigma}_{ij}^{(n)}$. Following (96), we can write

$$ \hat{\Sigma}_{ij}^{(n)} = \frac{1}{n-1} \sum_{l,s=1}^{n} X_{li}^{(n)} X_{js}^{(n)} (U^{(n)})_{ls}^{-1} - \frac{1}{n-1} \bar{X}_{j}^{(n)} \sum_{l,s=1}^{n} X_{li}^{(n)} (U^{(n)})_{ls}^{-1} \\ - \frac{1}{n-1} \bar{X}_{i}^{(n)} \sum_{l,s=1}^{n} X_{sj}^{(n)} (U^{(n)})_{ls}^{-1} + \frac{1}{n-1} \bar{X}_{i}^{(n)} \bar{X}_{j}^{(n)} \sum_{l,s=1}^{n} (U^{(n)})_{ls}^{-1}. \quad (98) $$

For simplicity, we denote as $A_{ij}^{(n)}$, $B_{ij}^{(n)}$, $C_{ij}^{(n)}$ and $D_{ij}^{(n)}$ the four terms in (98) respectively. First, let us derive their asymptotic expectations.

$$ \begin{aligned} \mathbb{E}(A_{ij}^{(n)}) &= \frac{1}{n-1} \sum_{l,s=1}^{n} (U^{(n)})_{ls}^{-1} \mathbb{E}(X_{li}^{(n)} X_{sj}^{(n)}) \\ &= \frac{1}{n-1} \sum_{l,s=1}^{n} (U^{(n)})_{ls}^{-1} \mu_{li}^{(n)} \mu_{sj}^{(n)} + \frac{\Sigma_{ij}}{n-1} \sum_{l,s=1}^{n} (U^{(n)})_{ls}^{-1} U_{sl}^{(n)} \\ &= \sum_{k,k'=1}^{K^*} \frac{1}{n-1} \sum_{l,s=1}^{n} (U^{(n)})_{ls}^{-1} \mathbb{1}\{\mu_l^{(n)} = \theta_k\} \mathbb{1}\{\mu_s^{(n)} = \theta_{k'}\} \theta_{ki} \theta_{k'j} + \frac{n}{n-1} \Sigma_{ij}. \end{aligned} $$

Using Lemma 3.5, we have

$$ \lim_{n \to \infty} \mathbb{E}(A_{ij}^{(n)}) = 2(\lambda - \lambda_0) \sum_{k=1}^{K^*} \pi_k \theta_{ki} \sum_{k=1}^{K^*} \pi_k \theta_{kj} + \lambda_0 \sum_{k=1}^{K^*} \pi_k \theta_{ki} \theta_{kj} + \Sigma_{ij}. \quad (99) $$

Then,

$$ \begin{aligned} \mathbb{E}(B_{ij}^{(n)}) &= \frac{1}{n(n-1)} \sum_{l,s,r=1}^{n} (U^{(n)})_{ls}^{-1} \mathbb{E}(X_{li}^{(n)} X_{rj}^{(n)}) \\ &= \frac{1}{n(n-1)} \sum_{l,s,r=1}^{n} (U^{(n)})_{ls}^{-1} \mu_{li}^{(n)} \mu_{rj}^{(n)} + \frac{\Sigma_{ij}}{n-1} \\ &= \frac{1}{n} \sum_{r=1}^{n} \mu_{rj}^{(n)} \frac{1}{n-1} \sum_{l,s}^{n} (U^{(n)})_{ls}^{-1} \mu_{li}^{(n)} + \frac{\Sigma_{ij}}{n-1} \\ &= \sum_{k=1}^{K^*} \frac{1}{n} \sum_{r=1}^{n} \mathbb{1}\{\mu_r^{(n)} = \theta_k\} \theta_{kj} \sum_{k=1}^{K^*} \frac{1}{n-1} \sum_{l,s=1}^{n} (U^{(n)})_{ls}^{-1} \mathbb{1}\{\mu_l^{(n)} = \theta_k\} \theta_{ki} + \frac{\Sigma_{ij}}{n-1}. \end{aligned} $$

Using the same reasoning as to prove Lemma 3.5, we have

$$ \lim_{n \to \infty} \frac{1}{n-1} \sum_{l,s=1}^{n} (U^{(n)})_{ls}^{-1} \mathbb{1}\{\mu_l^{(n)} = \theta_k\} = (2(\lambda - \lambda_0) + \lambda_0)\pi_k. $$

This, together with Assumption 3.1, yields

$$ \lim_{n \to \infty} \mathbb{E}(B_{ij}^{(n)}) = \lim_{n \to \infty} \mathbb{E}(C_{ij}^{(n)}) = (2(\lambda - \lambda_0) + \lambda_0) \sum_{k=1}^{K^*} \pi_k \theta_{kj} \sum_{k=1}^{K^*} \pi_k \theta_{ki}, \quad (100) $$
where $B_{ij}^{(n)}$ and $C_{ij}^{(n)}$ have the same expectation by symmetry. Finally,

$$ 
\begin{aligned}
\mathbb{E}(D_{ij}^{(n)}) &= \frac{1}{n^2(n-1)} \sum_{l,s=1}^{n} (U^{(n)})_{ls}^{-1} \sum_{r,r'=1}^{n} \mathbb{E}(X_{ri}^{(n)} X_{r'j}^{(n)}) \\
&= \frac{1}{n-1} \sum_{l,s=1}^{n} (U^{(n)})_{ls}^{-1} \left[ \frac{1}{n^2} \sum_{r,r'=1}^{n} \mu_{ri}^{(n)} \mu_{r'j}^{(n)} + \frac{\Sigma_{ij}}{n^2} \sum_{r,r'=1}^{n} U_{rr'}^{(n)} \right].
\end{aligned}
$$

Using the same reasoning as to prove Lemma 3.5, we have

$$ \lim_{n \to \infty} \frac{1}{n-1} \sum_{l,s=1}^{n} (U^{(n)})_{ls}^{-1} = 2(\lambda - \lambda_0) + \lambda_0. \qquad (101) $$

Moreover, we state that

$$ \lim_{n \to \infty} \frac{1}{n^2} \sum_{l,s=1}^{n} U_{ls}^{(n)} = 0. \qquad (102) $$

We prove (102) at the end of the proof. This claim, together with (101) and Assumption 3.1, yields

$$ \lim_{n \to \infty} \mathbb{E}(D_{ij}^{(n)}) = (2(\lambda - \lambda_0) + \lambda_0) \sum_{k=1}^{K^*} \pi_k \theta_{ki} \sum_{k=1}^{K^*} \pi_k \theta_{kj}. \qquad (103) $$

Consequently, following (99), (100) and (103), we have

$$ 
\begin{aligned}
\lim_{n \to \infty} \mathbb{E}(\hat{\Sigma}_{ij}^{(n)}) &= \Sigma_{ij} + \lambda_0 \left[ \sum_{k=1}^{K^*} \pi_k \theta_{ki} \theta_{kj} - \sum_{k=1}^{K^*} \pi_k \theta_{ki} \sum_{k=1}^{K^*} \pi_k \theta_{kj} \right] \\
&= \Sigma_{ij} + \lambda_0 \sum_{k=1}^{K^*} \pi_k (\theta_{ki} - \tilde{\theta}_i) (\theta_{kj} - \tilde{\theta}_j).
\end{aligned}
\qquad (104) $$

This is the first statement in (97). To prove the second one, we show that the variance of each term in (98) tends to zero. To do so, we need the explicit form of the non-centered 4-th moments of a Gaussian distribution. More precisely, if $X_1, \dots, X_4$ are four Gaussian random variables with $\mathbb{E}(X_i) = \mu_i$ and $\text{Cov}(X_i, X_j) = \sigma_{ij}$, for $i, j \in \{1, \dots, 4\}$, we need the explicit form of the quantity

$$ \mathbb{E}(X_1 X_2 X_3 X_4) - \mathbb{E}(X_1 X_2) \mathbb{E}(X_3 X_4). \qquad (105) $$

The first term can be derived using the moment generating function of a 4-dimensional normal distribution

$$ M_{(X_1, \dots, X_4)}(t_1, \dots, t_4) = \exp \left( \sum_{i=1}^{4} \mu_i t_i + \frac{1}{2} \sum_{i,j=1}^{n} \sigma_{ij} t_i t_j \right), $$

and computing

$$ \mathbb{E}(X_1 X_2 X_3 X_4) = \frac{\partial M_{(X_1, \dots, X_4)}(t_1, \dots, t_4)}{\partial t_1 \dots \partial t_4} \Big|_0. $$

Doing so, and using $\mathbb{E}(X_i X_j) = \mu_i \mu_j + \sigma_{ij}$, we can derive

$$ \mathbb{E}(X_1 X_2 X_3 X_4) - \mathbb{E}(X_1 X_2) \mathbb{E}(X_3 X_4) = \sigma_{13}\sigma_{24} + \sigma_{14}\sigma_{23} + \mu_1\mu_4\sigma_{23} + \mu_1\mu_3\sigma_{24} + \mu_2\mu_3\sigma_{14} + \mu_2\mu_4\sigma_{13}. \qquad (106) $$
We are ready to prove that $\mathrm{Var}(\hat{\Sigma}_{ij}^{(n)})$ tends to zero. First, using $\mathrm{Var}(X) = \mathbb{E}(X^2) - \mathbb{E}(X)^2$, we have

$$ \mathrm{Var} (A_{ij}^{(n)}) = \frac{1}{(n-1)^2} \sum_{l,s,k,r=1}^{n} (U^{(n)})_{sl}^{-1} (U^{(n)})_{kr}^{-1} [\mathbb{E}(X_{li} X_{sj} X_{ri} X_{kj}) - \mathbb{E}(X_{li} X_{sj}) \mathbb{E}(X_{ki} X_{rj})]. \quad (107) $$

Using (106), we can separate (107) into the following six terms:

$$ \mathrm{Var} (A_{ij}^{(n)}) = \frac{\Sigma_{ii} \Sigma_{jj}}{(n-1)^2} \sum_{l,s,k,r=1}^{n} (U^{(n)})_{ls}^{-1} (U^{(n)})_{kr}^{-1} U_{lk}^{(n)} U_{sr}^{(n)} \quad (108) $$

$$ + \frac{\Sigma_{ij}^2}{(n-1)^2} \sum_{l,s,k,r=1}^{n} (U^{(n)})_{ls}^{-1} (U^{(n)})_{kr}^{-1} U_{lr}^{(n)} U_{sk}^{(n)} \quad (109) $$

$$ + \frac{\Sigma_{jj}}{(n-1)^2} \sum_{l,s,k,r=1}^{n} (U^{(n)})_{ls}^{-1} (U^{(n)})_{kr}^{-1} U_{sr}^{(n)} \mu_{li}^{(n)} \mu_{ki}^{(n)} \quad (110) $$

$$ + \frac{\Sigma_{ij}}{(n-1)^2} \sum_{l,s,k,r=1}^{n} (U^{(n)})_{ls}^{-1} (U^{(n)})_{kr}^{-1} U_{sk}^{(n)} \mu_{li}^{(n)} \mu_{rj}^{(n)} \quad (111) $$

$$ + \frac{\Sigma_{ij}}{(n-1)^2} \sum_{l,s,k,r=1}^{n} (U^{(n)})_{ls}^{-1} (U^{(n)})_{kr}^{-1} U_{lr}^{(n)} \mu_{ki}^{(n)} \mu_{sj}^{(n)} \quad (112) $$

$$ + \frac{\Sigma_{ii}}{(n-1)^2} \sum_{l,s,k,r=1}^{n} (U^{(n)})_{ls}^{-1} (U^{(n)})_{kr}^{-1} U_{lk}^{(n)} \mu_{sj}^{(n)} \mu_{rj}^{(n)}. \quad (113) $$

Each of these terms tend to zero when $n \to \infty$. For (108), we have

$$ \frac{\Sigma_{ii} \Sigma_{jj}}{(n-1)^2} \sum_{l,s,k,r=1}^{n} (U^{(n)})_{ls}^{-1} (U^{(n)})_{kr}^{-1} U_{lk}^{(n)} U_{sr}^{(n)} = \frac{\Sigma_{ii} \Sigma_{jj}}{(n-1)^2} \sum_{l,s,r=1}^{n} (U^{(n)})_{ls}^{-1} U_{sr}^{(n)} \delta_{lr} \\ = \frac{\Sigma_{ii} \Sigma_{jj}}{(n-1)^2} \sum_{l,s=1}^{n} (U^{(n)})_{ls}^{-1} U_{sl}^{(n)} = \frac{\Sigma_{ii} \Sigma_{jj}}{(n-1)^2} \sum_{l=1}^{n} \delta_{ll} = \frac{n}{(n-1)^2} \Sigma_{ii} \Sigma_{jj} \underset{n \to \infty}{\longrightarrow} 0. $$

Identically we can show that (109) tends to zero. For (110), we have

$$ \frac{\Sigma_{jj}}{(n-1)^2} \sum_{l,s,k,r=1}^{n} (U^{(n)})_{ls}^{-1} (U^{(n)})_{kr}^{-1} U_{sr}^{(n)} \mu_{li}^{(n)} \mu_{ki}^{(n)} \\ = \frac{\Sigma_{jj}}{(n-1)^2} \sum_{l,k,r=1}^{n} (U^{(n)})_{kr}^{-1} \delta_{lr} \mu_{li}^{(n)} \mu_{ki}^{(n)} = \frac{\Sigma_{jj}}{(n-1)^2} \sum_{l,k=1}^{n} (U^{(n)})_{kl}^{-1} \mu_{li}^{(n)} \mu_{ki}^{(n)} \\ = \sum_{r,r'=1}^{K^*} \frac{\Sigma_{jj}}{(n-1)^2} \sum_{l,k=1}^{n} (U^{(n)})_{kl}^{-1} \mathbb{1}\{\mu_l^{(n)} = \theta_r\} \mathbb{1}\{\mu_k^{(n)} = \theta_{r'}\} \mu_{li}^{(n)} \mu_{ki}^{(n)} \theta_{ri} \theta_{r'i} \underset{n \to \infty}{\longrightarrow} 0, $$

where the limit is derived using Lemma 3.5. The same reasoning is used to show that (111), (112) and (113) tend to zero when $n \to \infty$. Therefore, we have $\lim_{n \to \infty} \mathrm{Var}(A_{ij}^{(n)}) = 0$. The same strategy, together with (101) and (102), is used to show that $\lim_{n \to \infty} \mathrm{Var}(B_{ij}^{(n)}) = \lim_{n \to \infty} \mathrm{Var}(C_{ij}^{(n)}) = \lim_{n \to \infty} \mathrm{Var}(D_{ij}^{(n)}) = 0$. Thus, we have (95). Note that the sum in (95) can be written as the *ij* term of a matrix. Indeed, we have

$$ \hat{\Sigma}_{ij}^{(n)} - \Sigma_{ij} \xrightarrow{p} \lambda_0 (\Theta^T \mathrm{diag}(\pi_1, \dots, \pi_{K^*}) \Theta)_{ij}, \quad (114) $$
where $\Theta$ is a $p \times K^*$ matrix having as entries $\Theta_{ij} = \theta_{ij} - \tilde{\theta}_j$. As $\lambda_0, \pi_1, \dots, \pi_{K^*} \ge 0$, the matrix $\lambda_0(\Theta^T \text{diag}(\pi_1, \dots, \pi_{K^*}) \Theta)$ is positive semi-definite, so the entries of $\hat{\Sigma} (\mathbf{X}^{(n)}) - \Sigma$ converge in probability to the entries of a positive semi-definite matrix. Note that, as both $\hat{\Sigma} (\mathbf{X}^{(n)})$ and $\Sigma$ are positive definite, the eigenvalues of their difference are real. Finally, since the eigenvalues depend continuously on the entries of the matrix, the eigenvalues of $\hat{\Sigma} (\mathbf{X}^{(n)}) - \Sigma$ converge in probability to the eigenvalues of a positive semi-definite matrix, which are non-negative. Therefore, we have (35).

Let us conclude by showing (102). To do so, note that we can write,

$$1 = \frac{1}{n} \sum_{k,l,s=1}^{n} (U^{(n)})_{lk}^{-1} U_{ks}^{(n)} = \frac{2}{n} \sum_{s=1}^{n} \sum_{r=1}^{n-1} \sum_{i=1}^{n-r} (U^{(n)})_{ii+r}^{-1} U_{i+r s}^{(n)} + \frac{1}{n} \sum_{s,i=1}^{n} (U^{(n)})_{ii}^{-1} U_{is}^{(n)}.$$

Using the same reasoning as in the proof of Lemma 3.5, we have

$$1 = 2 \lim_{n \to \infty} \sum_{r=1}^{n-1} \lambda_r \left( \lim_{n \to \infty} \frac{1}{n} \sum_{i,s=1}^{n} U_{i+r s}^{(n)} \right) + \lambda_0 \lim_{n \to \infty} \frac{1}{n} \sum_{i,s=1}^{n} U_{is}^{(n)},$$

which diverges unless the third limit is finite, which implies (102). □

## B Non-maximal conditioning sets

The methodology presented in Section 2.2.2 sets up the framework to perform selective inference after hierarchical clustering. Exploring its adaptation to further clustering algorithms involves, as shown in [9], the redefinition of p-values by constraining the conditional event that define (p-GBW) and (p-gen). In this section, we revisit the procedure of post-clustering inference introduced in Section 2.2.2 and rewrite it in a more general form that allows its straightforward adaptation to the scenario where more conditioning is imposed.

When defining a p-value for (H0) that controls the selective type I error (5), one may think of conditioning only on having selected the pair of clusters that define the null hypothesis, i.e. on the event

$$M_{G_1, G_2}(\mathbf{X}) = \{G_1, G_2 \in \mathcal{C}(\mathbf{X})\}. \qquad (115)$$

However, this is generally not enough to ensure the analytical tractability of the p-value. When considering a matrix normal distribution for the p-dimensional observations, two further conditions are imposed as shown in [21]. Following Section 2.2.2, this corresponds to conditioning on the event

$$M_{G_1, G_2}(\mathbf{X}) \cap \{\pi_{\nu}^{\perp} \mathbf{X} = \pi_{\nu}^{\perp} \mathbf{x}, \text{dir}_{V_{G_1, G_2}}(\mathbf{X}^T \nu) = \text{dir}_{V_{G_1, G_2}}(\mathbf{x}^T \nu)\}, \qquad (116)$$

which is the maximal event for which any analytically tractable p-value has been shown to control (5) under the general model (gen-MN). If we denote by $T_{G_1, G_2}(\mathbf{X}, \mathbf{x})$ the second set in (116), we can rewrite (p-gen) as

$$p_{V_{G_1, G_2}}(\mathbf{x}; \{G_1, G_2\}) = \mathbb{P}_{H_0^{(G_1, G_2)}} \left( \|\mathbf{X}^T \nu\|_{V_{G_1, G_2}} \ge \|\mathbf{x}^T \nu\|_{V_{G_1, G_2}} \mid M_{G_1, G_2}(\mathbf{X}) \cap T_{G_1, G_2}(\mathbf{X}, \mathbf{x}) \right). \qquad (117)$$

Then, from Theorem 2.2 and its proof we can rewrite the truncation set in (p-tract) as

$$S_{V_{G_1, G_2}}(\mathbf{x}; \{G_1, G_2\}) = \{\phi \in \mathbb{R} : M_{G_1, G_2}(\mathbf{x}'_{V_{G_1, G_2}}(\phi))\}, \qquad (118)$$
where $\mathbf{x}'_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\phi)$ is defined in (12). Consequently, in the conditions of Theorem 2.2, (p-gen) is analytically tractable as

$$ p_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\}) = 1 - \mathbb{F}_p \left( \left\| \mathbf{x}^T \nu \right\|_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}, \left\{ \phi \ge 0 : M_{\mathcal{G}_1, \mathcal{G}_2} \left( \mathbf{x}'_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\phi) \right) \right\} \right), \quad (119) $$

where $\mathbb{F}_p$ is defined in Theorem 2.2. Uncoupling $M_{\mathcal{G}_1, \mathcal{G}_2}(\mathbf{X})$ and $T_{\mathcal{G}_1, \mathcal{G}_2}(\mathbf{X}, \mathbf{x})$ in (117) allows us to characterize the null distribution of the p-value in terms of the conditioning event (115). This is useful to study the scenarios where, for technical reasons, subsets of (115) are chosen to define the p-value for (H0). This is the case in [9], where the framework of [21] under model (ind-MN) has been adapted to perform selective inference after k-means clustering. To allow the efficient computation of their truncation set, the authors condition on $T_{\mathcal{G}_1, \mathcal{G}_2}(\mathbf{X}, \mathbf{x})$ but also on all the intermediate clustering assignments for the n observations [9, Equation (9)], which is a subset of (115). In accordance with (118) and (119), this more restrictive conditioning yielded the same p-value (p-GBW) as in [21] except from a different truncation set, based on the finer conditioning event. The following result characterizes this framework under our general model (gen-MN) and for an arbitrary non-maximal conditioning event. As such, it is a generalization of Theorem 2.2.

**Theorem B.1.** In the conditions of Theorem 2.2, let $\emptyset \ne E_{\mathcal{G}_1, \mathcal{G}_2}(\mathbf{X}) \subset M_{\mathcal{G}_1, \mathcal{G}_2}(\mathbf{X})$ for any $(\mathcal{G}_1, \mathcal{G}_2) \in C_{[n]}$. Then, the quantity

$$ p_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\}; E_{\mathcal{G}_1, \mathcal{G}_2}) = \mathbb{P}_{H_0^{\{\mathcal{G}_1, \mathcal{G}_2\}}} \left( \left\| \mathbf{X}^T \nu \right\|_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}} \ge \left\| \mathbf{x}^T \nu \right\|_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}} \middle| E_{\mathcal{G}_1, \mathcal{G}_2}(\mathbf{X}) \cap T_{\mathcal{G}_1, \mathcal{G}_2}(\mathbf{X}, \mathbf{x}) \right) \quad (120) $$

is a p-value for (H0) that controls the selective type I error for clustering (5) at level $\alpha$. Furthermore, it satisfies

$$ p_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\}; E_{\mathcal{G}_1, \mathcal{G}_2}) = 1 - \mathbb{F}_p \left( \left\| \mathbf{x}^T \nu \right\|_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}, \left\{ \phi \ge 0 : E_{\mathcal{G}_1, \mathcal{G}_2} \left( \mathbf{x}'_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\phi) \right) \right\} \right), \quad (121) $$

where $\mathbb{F}_p(t, S)$ is the cumulative distribution function of a $\chi_p$ random variable truncated to the set $S$ and $\mathbf{x}'_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\phi)$ is defined in (12).

**Proof of Theorem B.1.** We omit the proof of (121) as it is identical to the one of (p-tract). Here, we show that the p-values defined using a non-maximal conditioning set $E(\mathbf{X}) \subset M(\mathbf{X})$ as (120) control the selective type I error for clustering (5). First, note that we have

$$ \mathbb{P}_{H_0^{\{\mathcal{G}_1, \mathcal{G}_2\}}} \left( p_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\}; E) \le \alpha \middle| E(\mathbf{X}) \cap T(\mathbf{X}) \right) = \alpha \quad (122) $$

following (120), for any $\alpha \in (0, 1)$. For simplicity, we will denote

$$ A = 1 \{ p_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\}; E) \le \alpha \}. \quad (123) $$
Then, following a similar reasoning as in the proof of [21, Theorem 1] and the tower property of conditional expectation, we can write

$$ \mathbb{P}_{H_0^{\{\mathcal{G}_1, \mathcal{G}_2\}}} \left( p_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\mathbf{x}; \{\mathcal{G}_1, \mathcal{G}_2\}; E) \le \alpha \mid M(\mathbf{X}) \right) = \mathbb{E}_{H_0^{\{\mathcal{G}_1, \mathcal{G}_2\}}} \left( A \mid M(\mathbf{X}) \right) \quad (124) $$

$$ = \mathbb{E}_{H_0^{\{\mathcal{G}_1, \mathcal{G}_2\}}} \left[ \mathbb{E}_{H_0^{\{\mathcal{G}_1, \mathcal{G}_2\}}} \left( A \mid M(\mathbf{X}) \cap E(\mathbf{X}) \cap T(\mathbf{X}) \right) \mid M(\mathbf{X}) \right] \quad (125) $$

$$ = \mathbb{E}_{H_0^{\{\mathcal{G}_1, \mathcal{G}_2\}}} \left[ \mathbb{E}_{H_0^{\{\mathcal{G}_1, \mathcal{G}_2\}}} \left( A \mid E(\mathbf{X}) \cap T(\mathbf{X}) \right) \mid M(\mathbf{X}) \right] = \mathbb{E}_{H_0^{\{\mathcal{G}_1, \mathcal{G}_2\}}} \left[ \alpha \mid M(\mathbf{X}) \right] = \alpha, \quad (126) $$

where the third equality follows from the fact $E(\mathbf{X}) \subset M(\mathbf{X})$ and the last equality follows from (122). □

Note that, following (119), replacing $E_{\mathcal{G}_1, \mathcal{G}_2}(\mathbf{X})$ by $M_{\mathcal{G}_1, \mathcal{G}_2}(\mathbf{X})$ yields exactly Theorem 2.2. Once again, the efficient computation of (121) depends on the efficient computation of the truncation set $E_{\mathcal{G}_1, \mathcal{G}_2}(\mathbf{x}'_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\phi))$. As shown for the maximal conditioning event in Lemma 2.3, it suffices to characterize the truncation set when the perturbed data set $\mathbf{x}'$ is defined with respect to any norm.

**Lemma B.2.** Let $\mathbf{x}$ be a realization of $\mathbf{X}$ and $\mathcal{G}_1, \mathcal{G}_2$ an arbitrary pair of clusters in $\mathcal{C}(\mathbf{x})$. Let $\mathbf{x}'$ denote the set (14) defined in [21, Equation (12)]. Then,

$$ E_{\mathcal{G}_1, \mathcal{G}_2} \left( \mathbf{x}'_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}(\phi) \right) = \frac{\|\mathbf{x}^T \nu\|_{\mathbf{V}_{\mathcal{G}_1, \mathcal{G}_2}}}{\|\mathbf{x}^T \nu\|_2} E_{\mathcal{G}_1, \mathcal{G}_2} \left( \mathbf{x}'(\phi) \right). \quad (127) $$

The proof of Lemma B.2 is omitted as it is identical to that of Lemma 2.3. In [9], the authors characterized $E_{\mathcal{G}_1, \mathcal{G}_2}(\mathbf{x}'(\phi))$ when $E_{\mathcal{G}_1, \mathcal{G}_2}$ corresponds to all intermediate clustering assignments of a k-means algorithm. Therefore, we can benefit from their efficient computation procedure and compute the truncation set under model (gen-MN) using Lemma B.2. As such, we are able to perform selective inference after k-means clustering for $\mathbf{U} \in \mathcal{CS}(n)$ and arbitrary $\mathbf{\Sigma}$. The estimation procedure presented in Section 3 remains identical for this case.

## C Supplementary Figures

<figure><img src="image_14.png" /><figcaption>Figure C.1: HAC dendrogram for the Hst5 protein ensemble data, with the six estimated clusters marked with colored rectangles.</figcaption></figure>
# D Additional numerical simulations

In this section we describe the numerical experiments illustrated in Figures 1 and 2 and present the results of the simulations described in Section 4 when C is a k-means or a hierarchical agglomerative clustering (HAC) algorithm with centroid, single and complete linkages.

## D.1 Numerical simulation of Figure 1

Figure 1 simulates the null distribution of p-values defined in [21] when data present dependence structures between observations and features, and p-values are computed assuming (ind-MN). We consider the general matrix normal model $X \sim MN_{n \times p}(\mu, U, \Sigma)$, where we set $\mu = 0_{n \times p}$, that is, the global null hypothesis. The matrices $U \in MN_{n \times n}(\mathbb{R})$ and $\Sigma \in MN_{p \times p}(\mathbb{R})$ encode the dependence structure between observations and features respectively. We choose U the covariance matrix of a stationary auto-regressive process of first order, AR(1), whose entries are given by $U_{ij} = \phi \rho^{|i-j|}$, for $\phi > 0$ and $|\rho| < 1$. The dependence between features is given by a Toeplitz matrix with entries $\Sigma_{ij} = 1 + 1/|i-j|$. We choose $\phi = 1$, $\rho = 0.2$ and generate $M = 2000$ realizations of $X$. For each one, we set the HAC algorithm with average linkage to choose three clusters and test for the difference in means of a pair of randomly selected clusters. The p-values are computed using the approach defined in [21] assuming that $X$ follows (ind-MN) with $\sigma^2 = 2$, that is, neglecting the off-diagonal entries of the covariance matrices $U$ and $\Sigma$.

## D.2 Numerical simulation of Figure 2

Figure 2 illustrates the effect of whitening matrix normal data with dependent observations and features and performing post-clustering inference assuming (ind-MN) afterwards. Data were first simulated from the general model (gen-MN) with $n = 100$, $p = 2$. We set $U$ as the covariance matrix of an AR(1) process, that is, $U_{ij} = \phi \rho^{|i-j|}$ for $\phi > 0$ and $|\rho| < 1$. We chose $\phi = 1$ and $\rho = 0.2$ The dependence between features was encoded by a Toeplitz matrix $\Sigma$ with entries $\Sigma_{ij} = 1 + 1/|i-j|$. The mean matrix $\mu$ divided the observations into three clusters and its entries were given by:

$$ \mu_i = \begin{cases} (-5, 0, \dots, 0) & \text{if } i \le \lfloor \frac{n}{3} \rfloor, \\ (0, \dots, 0, 5\sqrt{3}) & \text{if } \lfloor \frac{n}{3} \rfloor < i \le \lfloor \frac{2n}{3} \rfloor, \quad \forall i \in [n]. \\ (5, 0, \dots, 0) & \text{otherwise,} \end{cases} $$

The sample drawn from this model is presented in Figure 2(a). Its observations are classified into three groups using the k-means algorithms and compared using the p-values (p-gen) presented in this work, that account for the dependence structures $U$ and $\Sigma$. In panels (b,c), data is whitened by taking the transformation $(\Sigma \otimes U)^{-\frac{1}{2}} \text{vec}(X)$ and de-vectorizing the resulting random vector into a $n \times p$ matrix. Then, observations are classified into three groups using k-means (b) and HAC with average linkage (c) algorithms and the differences between cluster means are tested using the approaches proposed in [9] (b) and [21] (c), that assume model (ind-MN).

## D.3 Numerical analysis of (p-Gamma)

In this Section, we simulate the distribution of (p-Gamma) under a global null hypothesis, that is, setting $\mu = 0_{n \times p}$. Following Proposition 2.5, the quantity (p-Gamma) has the closed form (23), allowing its implementation in practice. We follow the same strategy as in Section D.1, generating $M = 2000$
<figure><img src="image_15.png" /><figcaption>Figure D.1: Empirical cumulative distribution functions (ECDF) of quantities (23) with $\mathcal{C}$ being a hierarchical agglomerative clustering algorithm (HAC) with complete linkage. The ECDF were computed from $M = 2000$ realizations of (gen-MN) under the three dependence settings (D4), (D5) and (D6) with $\boldsymbol{\mu} = \mathbf{0}_{n \times p}$, $n = 20$ and $p = 5$.</figcaption></figure>

realizations of $\mathbf{X} \sim \mathcal{MN}_{n \times p}(\mathbf{0}_{n \times p}, \mathbf{U}, \mathbf{\Sigma})$, setting the HAC algorithm to choose three clusters and computing (23) for a pair of randomly selected groups. We choose $\mathbf{\Sigma}$ to be a diagonal matrix with entries $\mathbf{\Sigma}_{ii} = 1 + 1/i$, and repeat the simulation under the following three settings:

(D4) $\mathbf{U}$ is a diagonal matrix with entries $U_{ii} = 1 + 1/i$.

(D5) $\mathbf{U}$ is the covariance matrix of an AR(1) model with $\sigma = 1$ and $\rho = 0.1$.

(D6) $\mathbf{U}$ is the covariance matrix of an AR(2) model with $\sigma = 1$, $\beta_1 = 0.4$ and $\beta_2 = 0.1$.

Note that the truncation set in (23) has slightly changed with respect to (11), due to the relaxation of the direction equality in (p-Gamma), that now includes the event $\{\text{dir}(\mathbf{X}^T\nu) = -\text{dir}(\mathbf{x}^T\nu)\}$. As shown in Proposition 2.5, this yields a broader truncation set (24) including also perturbations in the sense of $-\mathbf{x}^T\nu$. Adapting the efficient characterization of (11) to this setting is not straightforward. However, this is immediate under a Monte Carlo computation of (23), as we only need to replace $\mathcal{C}(\mathbf{x}'(\omega_i))$ by $\mathcal{C}(\mathbf{x}'(\pm\omega_i))$ in (17). As this is sufficient for the purpose of this analysis, we limit this experience to HAC clustering with complete linkage. Results, showing that selective type I error is not controlled in any of the previous settings, are presented in Figure D.1.

## D.4 Additional numerical simulations of Section 4

In this section, we present the counterparts of Figures 3, 4, 6, 7, 8 and 9 for $k$-means and HAC with centroid, single and complete linkage. In Figures D.6 and D.7, the simulation for $k$-means was performed for $\delta \in \{6, 8, 10\}$, as the proportion of samples for which the null hypothesis held was very low for $\delta = 4$.
<figure><img src="image_16.png" /><figcaption>Figure D.2: Empirical cumulative distribution functions (ECDF) of *p*-values (p-gen) with *C* being a hierarchical agglomerative clustering algorithm (HAC) with centroid (a-c), single (d-f) and complete (g-i) linkage and a *k*-means algorithm (j-l). The ECDF were computed from *M* = 2000 realizations of (gen-MN) under the three dependence settings (*D1*), (*D2*) and (*D3*) with *μ* = *0*<sub>*n*×*p*</sub>, *n* = 100 and *p* ∈ {5, 20, 50}.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure D.3: Empirical cumulative distribution functions (ECDF) of *p*-values (hat-p-tract) with *C* being a HAC algorithm with centroid (a-c), single (d-f) and complete (g-i) linkage and a *k*-means algorithm (j-l). The ECDF were computed from *M* = 5000 realizations of (gen-MN) under the three dependence settings (*D1*), (*D2*) and (*D3*) with *n* = 100, *p* = 5 and *μ* given by (40). Only samples for which the null hypothesis held were kept, as described in Section 4.2.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure D.4: Empirical cumulative distribution functions (ECDF) of *p*-values (p-gen) with *C* being a hierarchical agglomerative clustering algorithm (HAC) with centroid (a-c), single (d-f) and complete (g-i) linkage and a *k*-means algorithm (j-l). The ECDF were computed from *M* = 2000 realizations of (gen-MN) under the three dependence settings (*D4*), (*D5*) and (*D6*) with *μ* = *0*<sub>*n*×*p*</sub>, *n* = 100 and *p* ∈ {5, 20, 50}.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure D.5: Empirical cumulative distribution functions (ECDF) of *p*-values (hat-p-tract) with *C* being a HAC algorithm with centroid (a-c), single (d-f) and complete (g-i) linkage and a *k*-means algorithm (j-l). The ECDF were computed from *M* = 5000 realizations of (gen-MN) under the three dependence settings (*D4*), (*D5*) and (*D6*) with *n* = 100, *p* = 5 and *μ* given by (40). Only samples for which the null hypothesis held were kept, as described in Section 4.2.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure D.6: Empirical cumulative distribution functions (ECDF) of *p*-values (hat-*p*-tract) with *C* being a HAC algorithm with centroid (a-c), single (d-f) and complete (g-i) linkage and a *k*-means algorithm (j-l). The ECDF were computed from *M* = 5000 realizations of (gen-MN) under the three dependence settings (*D7*), (*D8*) and (*D9*) with *n* = 50, *p* = 5 and *μ* given by (40). Only samples for which the null hypothesis held were kept, as described in Section 4.4.2.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure D.7: Empirical cumulative distribution functions (ECDF) of *p*-values (hat-p-tract) with *C* being a HAC algorithm with centroid, single and complete linkage and a *k*-means algorithm. The ECDF were computed from *M* = 5000 realizations of (gen-MN) as described in Section 4.4.3 with *n* = 50, *p* = 5 and *μ* given by (40) with *δ* ∈ {4, 6, 8} for HAC and *δ* ∈ {6, 8, 10} for *k*-means. Only samples for which the null hypothesis held were kept, as described in Section 4.4.3.</figcaption></figure>
References

[1] S. E. Ahmed, S. Fallahpour, D. von Rosen, and T. von Rosen. Estimation of Several Intraclass Correlation Coefficients. *Comm. Statist. Simulation Comput.*, Oct. 2015.

[2] M. P. Allen and D. J. Tildesley. *Computer Simulation of Liquids*. Oxford University Press, 06 2017.

[3] R. Appadurai, J. K. Koneru, M. Bonomi, P. Robustelli, and A. Srivastava. Clustering heterogeneous conformational ensembles of intrinsically disordered proteins with t-distributed stochastic neighbor embedding. *Journal of Chemical Theory and Computation*, June 2023.

[4] M. S. Bartlett. An Inverse Matrix Adjustment Arising in Discriminant Analysis. *The Annals of Mathematical Statistics*, 22(1):107 – 111, 1951.

[5] P. Bernadó, L. Blanchard, P. Timmins, D. Marion, R. W. H. Ruigrok, and M. Blackledge. A structural model for unfolded proteins from residual dipolar couplings and small-angle x-ray scattering. *Proceedings of the National Academy of Sciences of the United States of America*, 102(47):17002–17007, 2005.

[6] A. R. Camacho-Zarco, S. Kalayil, D. Maurin, N. Salvi, E. Delaforge, S. Milles, M. R. Jensen, D. J. Hart, S. Cusack, and M. Blackledge. Molecular basis of host-adaptation interactions between influenza virus polymerase PB2 subunit and ANP32a. *Nature Communications*, 11(1), July 2020.

[7] Y. Chen, S. Jewell, and D. Witten. More powerful selective inference for the graph fused lasso. *Journal of Computational and Graphical Statistics*, 32(2):577–587, 2023.

[8] Y. T. Chen and L. L. Gao. Testing for a difference in means of a single feature after clustering. *arXiv*, 2023.

[9] Y. T. Chen and D. M. Witten. Selective inference for k-means clustering. *Journal of Machine Learning Research*, 24(152):1–41, 2023.

[10] W. G. Cochran. The distribution of quadratic forms in a normal system, with applications to the analysis of covariance. *Math. Proc. Cambridge Philos. Soc.*, 30(2):178–191, Apr. 1934.

[11] A. Conev, M. M. Rigo, D. Devaurs, A. F. Fonseca, H. Kalavadwala, M. V. de Freitas, C. Clementi, G. Zanatta, D. A. Antunes, and L. E. Kavraki. EnGens: a computational framework for generation and analysis of representative protein conformational ensembles. *Briefings in Bioinformatics*, 24(4):bbad242, 07 2023.

[12] P. M. Crespo and J. Gutierrez-Gutierrez. On the elementwise convergence of continuous functions of hermitian banded toeplitz matrices. *IEEE Transactions on Information Theory*, 53(3):1168–1176, 2007.

[13] S. Demko, W. F. Moss, and P. Smith. Decay rates for inverses of band matrices. *Mathematics of Computation*, 43:491–499, 1984.

[14] H. Derksen and V. Makam. Maximum likelihood estimation for matrix normal models via quiver representations. *SIAM Journal on Applied Algebra and Geometry*, 5(2):338–365, 2021.

[15] M. Drton, A. Grosdos, and A. McCormack. Rational Maximum Likelihood Estimators of Kronecker Covariance Matrices. *arXiv*, Jan. 2024.
[16] M. Drton, S. Kuriki, and P. Hoff. Existence and uniqueness of the Kronecker covariance MLE. *Ann. Stat.*, 49(5):2721–2754, Oct. 2021.

[17] P. Dutilleul. The MLE algorithm for the matrix normal distribution. *Journal of Statistical Computation and Simulation*, 64(2):105–123, 1999.

[18] H. J. Dyson and P. E. Wright. Intrinsically unstructured proteins and their functions. *Nature Reviews Molecular Cell Biology*, 6:197–208, 2005.

[19] M. L. Eaton. *Multivariate Statistics*. SPIE, Jan. 2007.

[20] W. Fithian, D. Sun, and J. Taylor. Optimal inference after model selection, 2017. arXiv:1410.2597.

[21] L. L. Gao, J. Bien, and D. Witten. Selective inference for hierarchical clustering. *Journal of the American Statistical Association*, 10(0):1–11, 2022.

[22] R. M. Gray. Toeplitz and circulant matrices: A review. *Foundations and Trends® in Communications and Information Theory*, 2(3):155–239, 2006.

[23] A. Gupta and D. Nagar. *Matrix Variate Distributions*. Chapman and Hall/CRC, 2018.

[24] B. Hivert, D. Agniel, R. Thiébaut, and B. P. Hejblum. Post-clustering difference testing: Valid inference and practical considerations with applications to ecological and biological data. *Comput. Statist. Data Anal.*, 193:107916, May 2024.

[25] S. Holm. A simple sequentially rejective multiple test procedure. *Scandinavian Journal of Statistics*, 6(2):65–70, 1979.

[26] R. Horn and C. Johnson. *Matrix Analysis*. Matrix Analysis. Cambridge University Press, 2013.

[27] S. Jewell, P. Fearnhead, and D. Witten. Testing for a change in mean after changepoint detection. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 84(4):1082–1104, Apr. 2022.

[28] A. Kessel and N. Ben-Tal. *Introduction to Proteins*. Chapman and Hall/CRC, Mar. 2018.

[29] A. Klenke. *Probability Theory*. Springer International Publishing, Cham, Switzerland, 2020.

[30] T. Lazar, M. Guharoy, W. Vranken, S. Rauscher, S. J. Wodak, and P. Tompa. Distance-based metrics for comparing conformational ensembles of intrinsically disordered proteins. *Biophysical Journal*, 118(12):2952–2965, 2020.

[31] J. Leiner, B. Duan, L. Wasserman, and A. Ramdas. Data Fission: Splitting a Single Data Point. *J. Am. Stat. Assoc.*, 2023.

[32] A. Liljas, L. Liljas, J. Piskur, G. Lindblom, P. Nissen, and M. Kjeldgaard. *Textbook Of Structural Biology*. World Scientific Publishing, Singapore, 2009.

[33] K. Liu, J. Markovic, and R. Tibshirani. More powerful post-selection inference, with application to the lasso, 2018. arXiv:1801.09037.

[34] P. Mahalanobis. On the generalized distance in statistics. *Proceedings of the National Academy of Sciences, India (Calcutta)*, 2(1):44–55, 1936.

[35] K. Nishikawa, T. Ooi, Y. Isogai, and N. Saitô. Tertiary structure of proteins. i. representation and computation of the conformations. *Journal of the Physical Society of Japan*, 32(5):1331–1337, 1972.
[36] V. Ntranos, L. Yi, P. Melsted, and L. Pachter. A discriminative learning approach to differential expression analysis for single-cell RNA-seq. *Nature Methods*, 16(2):163–166, Jan. 2019.

[37] C. J. Oldfield and A. K. Dunker. Intrinsically disordered proteins and intrinsically disordered protein regions. *Annual Review of Biochemistry*, 83(1):553–584, 2014.

[38] V. Ozenne, F. Bauer, L. Salmon, J.-r. Huang, M. R. Jensen, S. Segard, P. Bernadó, C. Charavay, and M. Blackledge. Flexible-meccano: a tool for the generation of explicit ensemble descriptions of intrinsically disordered proteins and their associated experimental observables. *Bioinformatics*, 28(11):1463–1470, 2012.

[39] R. Pearce and Y. Zhang. Deep learning techniques have significantly impacted protein structure prediction and protein design. *Current Opinion in Structural Biology*, 68:194–207, June 2021.

[40] D. Phillips. British biochemistry, past and present. In *London Biochemical Society Symposia*, page 11. Academic Press, 1970.

[41] D. G. Rasines and G. A. Young. Splitting strategies for post-selection inference. *Biometrika*, 12 2022. asac070.

[42] A. Sagar, C. M. Jeffries, M. V. Petoukhov, D. I. Svergun, and P. Bernadó. Comment on the optimal parameters to derive intrinsically disordered protein conformational ensembles from small-angle X-ray scattering data using the ensemble optimization method. *Journal of Chemical Theory and Computation*, 17(4):2014–2021, 2021.

[43] J. Shao, S. W. Tanner, N. Thompson, and T. E. Cheatham. Clustering molecular dynamics trajectories: 1. characterizing the performance of different clustering algorithms. *Journal of Chemical Theory and Computation*, 3(6):2312–2334, Oct. 2007.

[44] I. Soloveychik and D. Trushin. Gaussian and robust Kronecker product covariance estimation: Existence and uniqueness. *J. Multivariate Anal.*, 149:92–113, July 2016.

[45] W. F. Trench. Asymptotic distribution of the spectra of a class of generalized kac–murdock–szegö matrices. *Linear Algebra and its Applications*, 294(1):181–192, 1999.

[46] A. Vandenbon and D. Diez. A clustering-independent method for finding differentially expressed genes in single-cell transcriptome data. *Nature Communications*, 11(1), Aug. 2020.

[47] A. P. Verbyla. A note on the inverse covariance matrix of the autoregressive process1. *Australian Journal of Statistics*, 27(2):221–224, 1985.

[48] J. Wise. The autocorrelation function and the spectral density function. *Biometrika*, 42(1/2):151–159, 1955.

[49] J. Yeh. *Real Analysis*. World Scientific, 3rd edition, 2014.

[50] Y.-J. Yun and R. Foygel Barber. Selective inference for clustering with unknown variance. *Electronic Journal of Statistics*, 17(2):1923–1946, 2023.