1 Introduction↩︎

AI, particularly with the advent of foundation models, has rapidly become a general-purpose technology with broad scientific, industrial, and societal impact [@bommasani2021opportunities]. Recent advances in LLMs and generative models have enabled the automated production, transformation, and interaction with diverse digital content, including natural language, code, and multimodal data [@bubeck2023sparks]. While these capabilities create substantial opportunities in education, healthcare, scientific discovery, and software development, they also introduce AI-related artifacts that are difficult to monitor at scale, such as synthetic text, hallucinated or factually inconsistent outputs, watermark traces, and adversarial perturbations [@goodfellow2015explaining; @carlini2017towards; @kirchenbauer2023watermarking]. These artifacts often manifest as observable statistical irregularities in model outputs or learned representations, for which detection provides a scalable and cost-effective mechanism for AI oversight, especially when the underlying model is inaccessible or continually changing. Existing detection methods can be broadly divided into white-box and black-box paradigms. White-box methods assume access to internal model information, such as parameters, gradients, logits, or decoding probabilities, and can exploit model-specific signals or embedded mechanisms [@kirchenbauer2023watermarking]. In contrast, black-box methods operate only on observable outputs or representations obtained from accessible feature extractors, making them more suitable for proprietary, evolving, or API-only systems [@manakul2023selfcheckgpt]. From a practical perspective, we therefore focus on black-box detection and discuss this class of methods below.

Despite being studied separately, LLM-generated text detection, hallucination detection, adversarial example detection, and watermark detection can all be viewed as determining whether a sample deviates from a reference distribution in a suitable representation space. LLM-generated text detection can be viewed as identifying distributional discrepancies between human-written and model-generated text [@gehrmann2019gltr; @mitchell2023detectgpt]; hallucination detection seeks to find outputs inconsistent with factual, contextual, or self-consistency-based evidence [@ji2023survey; @manakul2023selfcheckgpt]; adversarial example detection aims to recognize abnormal inputs or representations [@goodfellow2015explaining; @carlini2017towards; @metzen2017detecting; @tian2024abnormal]; and watermark detection uncovers deliberately embedded statistical signals in generated content [@kirchenbauer2023watermarking]. A wide range of black-box detection methods have been developed under this perspective. Likelihood-based methods exploit token- or sequence-level probability patterns assigned by language models [@gehrmann2019gltr; @mitchell2023detectgpt]. Classifier-based methods train supervised or contrastive models on learned representations, with applications to LLM-generated text detection, hallucination detection, and adversarial detection [@solaiman2019release]. Watermark-specific methods embed and subsequently detect statistical signatures in generated outputs [@kirchenbauer2023watermarking]. However, these methods are often task-specific, rely on tailored features or supervised training data, and require retraining or careful data curation when transferred to new domains or models. They also rarely provide a unified statistical treatment of heterogeneous reference samples arising from multiple classes, sources, domains, or generation mechanisms.

To address these limitations, we propose a general black-box detection framework based on MDS. The key idea is to characterize positive samples in a learned representation space and detect negative samples that deviate from this characterization. Specifically, for positive samples, we use an appropriate pretrained model to extract deep representations, estimate their location and covariance structure, and compute MDS for detection. For MDS estimation, we adopt both casewise MCD [@hubert2010minimum] and cellwise MCD [@raymaekers2024cellwise] as the basic principles. Since positive samples may consist of multiple classes, we further propose joint estimators for both the casewise and cellwise settings using a reparameterization technique, together with efficient optimization algorithms. We establish the convergence of the proposed algorithms and the high breakdown point properties of the resulting estimators. Extensive experiments on LLM-generated text detection, hallucination detection, watermark detection, and adversarial example detection demonstrate the effectiveness and generality of the proposed framework.

It is worth noting that our approach is inspired by [@lee2018simple], who applied class-conditional Gaussian modeling to adversarial and out-of-distribution detection. However, their method assumes a common covariance matrix across classes and relies on empirical covariance estimation, limiting its robustness in heterogeneous or contaminated settings. In summary, our contributions are as follows:

1. We propose a unified black-box detection framework for a broad class of AI-related content and artifacts, including LLM-generated text, hallucinations, watermarked text, and adversarial examples. The framework is flexible and broadly applicable across different detection tasks.

2. We develop joint estimation methods for multiple covariance matrices in low-dimensional settings. The proposed methods cover both casewise and cellwise contamination mechanisms and allow positive samples to contain multiple classes with shared and class-specific covariance structures.

3. We introduce a suitable definition of the breakdown point for joint estimators and establish the high-breakdown-point properties of the proposed joint estimators.

The remainder of this paper is organized as follows. Section 2 introduces the unified detection framework. Sections 3 and 4 present the proposed MCD estimators and their multi-class extensions. Section 5 develops optimization algorithms, and Section 6 establishes convergence guarantees and high-breakdown-point properties. Section 7 reports experimental results, and Section 8 concludes the paper. Implementation details and proofs are provided in the appendices.

2 Unified detection framework via Mahalanobis Distance Scores↩︎

In this section, we introduce the proposed unified framework for detecting AI-related content and artifacts. For ease of exposition, we first consider the single-class positive-sample setting in Section 2.1, and then extend the framework to the multi-class setting in Section 2.2. Figure 1 provides a schematic illustration of the overall framework.

2.1 Single-class MDS detector↩︎

We first consider the setting where positive samples are assumed to come from a single class. For example, in image adversarial example detection, the positive samples may be clean images from one class, while the negative samples are their adversarially perturbed counterparts. Let the set of positive samples be \(\mathcal{D}_{1} = \{(\boldsymbol{x}_{i}, Y_{i}=1)\}^n_{i=1},\) and let the set of negative samples be \(\mathcal{D}_{2} = \bigcup_{k=1}^{K_{2}} \mathcal{D}_{2}^{(k)},\) where \(\mathcal{D}_{2}^{(k)} = \{(\tilde{\boldsymbol{x}}_{i}^{(k)},Y_{i}^{(k)}=-1)\}_{i=1}^{m_k}, k=1,2,\ldots,K_{2}.\) The complete dataset is denoted by \(\mathcal{D} = \mathcal{D}_{1} \cup \mathcal{D}_{2}\). Let \(\mathcal{M}\) be a pretrained model trained on positive samples and used to extract latent representations. For each input \(\boldsymbol{x} \in \mathbb{R}^{p}\), we obtain a representation \(\boldsymbol{z} \in \mathbb{R}^{q}\) from a chosen layer of \(\mathcal{M}\). The corresponding representation sets are denoted by \[\mathcal{Z}_{1} = \{(\boldsymbol{z}_{i}, Y_{i}=1)\}_{i=1}^{n}, \quad \mathcal{Z}_{2} = \bigcup_{k=1}^{K_{2}}\mathcal{Z}_{2}^{(k)},\] where \(\mathcal{Z}_{2}^{(k)}=\{(\tilde{\boldsymbol{z}}_{i}^{(k)},Y_{i}^{(k)}=-1)\}_{i=1}^{m_k}\), \(k=1,2,\ldots,K_2\). To characterize the positive class, we estimate its mean vector \(\boldsymbol{\mu}\) and covariance matrix \(\boldsymbol{\Sigma}\) based on the deep representations of positive samples in \(\mathcal{Z}_{1}\). This choice is natural because positive samples are typically abundant and easier to obtain, and pretrained models are usually trained primarily on positive samples. Since the positive samples used for estimation may still contain contamination, direct empirical estimation of \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\) can be unstable. Therefore, we use robust estimators, denoted by \(\widehat{\boldsymbol{\mu}}\) and \(\widehat{\boldsymbol{\Sigma}}\), whose construction will be detailed in Section 3. Using these estimates, the Mahalanobis distance scores for positive and negative samples are defined as \[d_{i} = (\boldsymbol{z}_{i} - \widehat{\boldsymbol{\mu}})^{\top} \widehat{\boldsymbol{\Sigma}}^{-1} (\boldsymbol{z}_{i} - \widehat{\boldsymbol{\mu}}), \quad \tilde{d}_{i}^{(k)} = (\tilde{\boldsymbol{z}}_{i}^{(k)} - \widehat{\boldsymbol{\mu}})^{\top} \widehat{\boldsymbol{\Sigma}}^{-1} (\tilde{\boldsymbol{z}}_{i}^{(k)} - \widehat{\boldsymbol{\mu}}), \quad k=1,2,\ldots,K_2.\] The resulting scores measure how far each sample is from the positive class in the learned representation space. Intuitively, a larger score indicates that the sample is less consistent with the positive class and therefore more likely to be negative. Classification can then be performed on the one-dimensional score sets \[\mathcal{S}_{1} = \{(d_i,Y_i=1)\}_{i=1}^{n}, \quad \mathcal{S}_{2}^{(k)} = \{(\tilde{d}_i^{(k)},Y_i^{(k)}=-1)\}_{i=1}^{m_k}, \quad k=1,2,\ldots,K_{2},\] by training a score-based binary classifier.

2.2 Multi-class MDS detector↩︎

We next consider the multi-class setting, where positive samples may come from multiple classes or subpopulations. For example, in image adversarial example detection, the positive samples may consist of clean images from different digit classes, while the negative samples are their adversarially perturbed counterparts. Let the positive and negative sample sets be \(\mathcal{D}_{1} = \bigcup_{k=1}^{K_{1}} \mathcal{D}_{1}^{(k)}\) and \(\mathcal{D}_{2} = \bigcup_{k=1}^{K_{2}} \mathcal{D}_{2}^{(k)},\) respectively, where \(\mathcal{D}_{1}^{(k)} = \{(\boldsymbol{x}_{i}^{(k)},Y_{i}^{(k)}=1)\}^{n_k}_{i=1}, k=1,2,\ldots,K_{1}\) and \(\mathcal{D}_{2}^{(k)} = \{(\tilde{\boldsymbol{x}}_{i}^{(k)},Y_{i}^{(k)}=-1)\}_{i=1}^{m_k}, k=1,2,\ldots,K_{2}.\) As in the single-class setting, a pretrained model \(\mathcal{M}\) is used to extract latent representations. The corresponding representation sets are denoted by \[\mathcal{Z}_{1} = \bigcup_{k=1}^{K_{1}}\mathcal{Z}_{1}^{(k)}, \quad \mathcal{Z}_{2} = \bigcup_{k=1}^{K_{2}}\mathcal{Z}_{2}^{(k)},\] where \(\mathcal{Z}_{1}^{(k)}=\{(\boldsymbol{z}_{i}^{(k)},Y_{i}^{(k)}=1)\}_{i=1}^{n_k}\) and \(\mathcal{Z}_{2}^{(k)}=\{(\tilde{\boldsymbol{z}}_{i}^{(k)},Y_{i}^{(k)}=-1)\}_{i=1}^{m_k}\). To characterize the positive classes, we estimate class-specific mean vectors \(\boldsymbol{\mu}^{(k)}\) and covariance matrices \(\boldsymbol{\Sigma}^{(k)}\) based on the deep representations of positive samples in \(\mathcal{Z}_{1}^{(k)}\), \(k=1,2,\ldots,K_{1}\). This multi-class formulation allows different positive classes to have distinct geometric structures in the representation space. Since the positive samples used for estimation may still contain contamination, direct empirical estimation of \(\boldsymbol{\mu}^{(k)}\) and \(\boldsymbol{\Sigma}^{(k)}\) can be unstable. Therefore, we use joint robust estimators, denoted by \(\{\widehat{\boldsymbol{\mu}}^{(k)}\}_{k=1}^{K_{1}} \text{ and } \{\widehat{\boldsymbol{\Sigma}}^{(k)}\}_{k=1}^{K_{1}},\) whose construction will be detailed in Section 4.

Using these estimates, the Mahalanobis distance score for each sample is defined as its minimum distance to the positive classes: \[\begin{align} d_{i}^{(k)} &= \min_{1\leq j\leq K_{1}} \left\{ (\boldsymbol{z}_{i}^{(k)}-\widehat{\boldsymbol{\mu}}^{(j)})^{\top} (\widehat{\boldsymbol{\Sigma}}^{(j)})^{-1} (\boldsymbol{z}_{i}^{(k)}-\widehat{\boldsymbol{\mu}}^{(j)}) \right\}, \quad k=1,2,\ldots,K_{1}, \\ \tilde{d}_{i}^{(k)} &= \min_{1\leq j\leq K_{1}} \left\{ (\tilde{\boldsymbol{z}}_{i}^{(k)}-\widehat{\boldsymbol{\mu}}^{(j)})^{\top} (\widehat{\boldsymbol{\Sigma}}^{(j)})^{-1} (\tilde{\boldsymbol{z}}_{i}^{(k)}-\widehat{\boldsymbol{\mu}}^{(j)}) \right\}, \quad k=1,2,\ldots,K_{2}. \end{align}\] To ensure that the Mahalanobis distance is well-defined in high-dimensional settings, a positive-definiteness constraint is imposed during covariance estimation. Details are provided in Sections 3-4. In addition, to mitigate the computational burden and numerical instability associated with high-dimensional covariance estimation, Gaussian random projection is applied to reduce the dimensionality of the representations before covariance estimation. Details are provided in 9.3. This minimum-distance rule reflects the intuition that a sample should be regarded as normal if it is close to at least one positive class, whereas a sample far from all positive classes is more likely to be negative. Classification can then be performed on the one-dimensional score sets \(\mathcal{S}_{1}^{(k)} = \{(d_i^{(k)},Y_i^{(k)}=1)\}_{i=1}^{n_k}, k=1,2,\ldots,K_{1},\) and \(\mathcal{S}_{2}^{(k)} = \{(\tilde{d}_i^{(k)},Y_i^{(k)}=-1)\}_{i=1}^{m_k}, k=1,2,\ldots,K_{2},\) by training a score-based binary classifier.

3 Single-class robust estimation↩︎

In the single-class MDS detector introduced in Section 2.1, the key step is to estimate the mean vector and covariance matrix of the deep representations of positive samples. Let \(\boldsymbol{Z}=(\boldsymbol{z}_1,\boldsymbol{z}_2,\ldots,\boldsymbol{z}_n)^\top \in \mathbb{R}^{n\times q}\) denote the representation matrix of positive samples. Since these samples may still contain contaminated observations or corrupted feature entries, empirical estimates of the mean and covariance matrix can be unstable. We therefore introduce two robust estimation methods: the casewise MCD, which is designed for contaminated observations, and the cellwise MCD, which is designed for contamination at the feature-entry level.

3.1 Casewise minimum covariance determinant↩︎

We first consider casewise contamination, where entire observations may be unreliable. To motivate the estimator, consider a Gaussian working model for the positive representations, \(\boldsymbol{z}_i \sim \mathcal{N}_{q}(\boldsymbol{\mu},\boldsymbol{\Sigma}), i=1, 2, \ldots, n.\) The density is \(f(\boldsymbol{z}_i;\boldsymbol{\mu},\boldsymbol{\Sigma}) = \frac{1}{(2\pi)^{q/2}|\boldsymbol{\Sigma}|^{1/2}} \exp\left\{ -\frac{1}{2}\mathrm{MD}^2(\boldsymbol{z}_i,\boldsymbol{\mu},\boldsymbol{\Sigma}) \right\},\) where \(\mathrm{MD}^2(\boldsymbol{z}_i,\boldsymbol{\mu},\boldsymbol{\Sigma}) = (\boldsymbol{z}_i-\boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1} (\boldsymbol{z}_i-\boldsymbol{\mu}).\) The negative log-likelihood, up to a constant factor, is \[\ell (\boldsymbol{\mu},\boldsymbol{\Sigma};\boldsymbol{z}_i) = \ln|\boldsymbol{\Sigma}| + \mathrm{MD}^2(\boldsymbol{z}_i,\boldsymbol{\mu},\boldsymbol{\Sigma}) + q\ln(2\pi).\] If all observations are reliable, minimizing \(\sum_{i=1}^{n}\ell (\boldsymbol{\mu},\boldsymbol{\Sigma};\boldsymbol{z}_i)\) leads to the usual empirical mean and covariance matrix. However, this estimator is sensitive to contaminated observations.

The casewise MCD estimator addresses this issue by selecting a subset of \(h\) observations with the smallest covariance determinant. Let \(\boldsymbol{H}=(w_1,w_2,\ldots,w_n)^\top\), where \(w_i \in \{0,1\}, i\in\{1,2,\ldots,n\}\), indicates whether the \(i\)-th observation is retained, with \(\sum_{i=1}^{n} w_i = h\). The casewise MCD estimator solves \[\label{eq:casewise95mcd95obj} \min_{\boldsymbol{\mu},\boldsymbol{\Sigma},\boldsymbol{H}} \sum_{i=1}^{n} w_i \ell (\boldsymbol{\mu},\boldsymbol{\Sigma};\boldsymbol{z}_i), \quad \text{s.t. } w_i\in\{0,1\},\; \sum_{i=1}^{n}w_i=h.\tag{1}\] For fixed weights \(\boldsymbol{H}\), the estimators are \[\widehat{\boldsymbol{\mu}} = \frac{1}{h} \sum_{i=1}^{n} w_i \boldsymbol{z}_i,\quad \widehat{\boldsymbol{\Sigma}} = \frac{1}{h} \sum_{i=1}^{n} w_i (\boldsymbol{z}_i-\widehat{\boldsymbol{\mu}}) (\boldsymbol{z}_i-\widehat{\boldsymbol{\mu}})^\top .\] Substituting these estimates into 1 shows that the casewise MCD is equivalent to finding the subset of size \(h\) whose covariance matrix has the minimum determinant.

It is worth noting that directly solving 1 is computationally challenging, since it involves a combinatorial search over all subsets of size \(h\). In practice, the casewise MCD is commonly computed using the Fast-MCD algorithm proposed by [@rousseeuw1999fast], which efficiently approximates the MCD solution through concentration steps and multiple initial subsets. The tuning parameter \(h\) controls the trade-off between robustness and efficiency: a smaller \(h\) yields higher robustness against contamination, whereas a larger \(h\) uses more observations and can improve statistical efficiency when contamination is mild. A common choice is \(h=\lfloor (n+q+1)/2 \rfloor\), which attains the highest breakdown point, while larger values such as \(h=\lfloor 0.75n \rfloor\) are often used when higher efficiency is desired [@hubert2010minimum; @rousseeuw1999fast]. The resulting \(\widehat{\boldsymbol{\mu}}\) and \(\widehat{\boldsymbol{\Sigma}}\) are then used in the single-class MDS detector.

3.2 Cellwise minimum covariance determinant↩︎

The casewise MCD removes entire observations and is therefore suitable when contamination affects whole samples. In deep representations, however, contamination may occur only in a subset of feature entries. Removing an entire observation in this case may discard useful information. The cellwise MCD addresses this problem by assigning binary weights at the cell level. Let \(\boldsymbol{W}=(w_{ij})_{1\leq i \leq n, 1\leq j \leq q}\in\{0,1\}^{n\times q}\) be a cellwise indicator matrix, where \(w_{ij}=1\) means that the \(j\)-th entry of \(\boldsymbol{z}_i\) is retained, and \(w_{ij}=0\) means that the entry is flagged as contaminated. Denote by \(\boldsymbol{w}_i\) the \(i\)-th row of \(\boldsymbol{W}\). For each observation, let \(\boldsymbol{z}_i^{(\boldsymbol{w}_i)}\) be the subvector of \(\boldsymbol{z}_i\) containing only retained entries. Similarly, \(\boldsymbol{\mu}^{(\boldsymbol{w}_i)}\) is the corresponding subvector of \(\boldsymbol{\mu}\), and \(\boldsymbol{\Sigma}^{(\boldsymbol{w}_i)}\) is the corresponding principal submatrix of \(\boldsymbol{\Sigma}\). Let \(q^{(\boldsymbol{w}_i)} = \sum_{j=1}^{q} w_{ij}\) denote the number of retained entries in the \(i\)-th observation.

Under the same Gaussian working model, the observed likelihood for the retained entries of \(\boldsymbol{z}_i\) is \[f\left(\boldsymbol{z}_{i}^{(\boldsymbol{w}_{i})}; \boldsymbol{\mu}^{(\boldsymbol{w}_{i})},\boldsymbol{\Sigma}^{(\boldsymbol{w}_{i})}\right) = \frac{ 1 }{ (2\pi)^{q^{(\boldsymbol{w}_i)}/2} |\boldsymbol{\Sigma}^{(\boldsymbol{w}_i)}|^{1/2} } \exp\left\{ -\frac{1}{2} \mathrm{MD}^{2} \left( \boldsymbol{z}_{i}^{(\boldsymbol{w}_{i})}, \boldsymbol{\mu}^{(\boldsymbol{w}_{i})},\boldsymbol{\Sigma}^{(\boldsymbol{w}_{i})} \right) \right\},\] where the partial Mahalanobis distance is \[\mathrm{MD}^{2} \left( \boldsymbol{z}_{i}^{(\boldsymbol{w}_{i})}, \boldsymbol{\mu}^{(\boldsymbol{w}_{i})},\boldsymbol{\Sigma}^{(\boldsymbol{w}_{i})} \right) = \left( \boldsymbol{z}_{i}^{(\boldsymbol{w}_{i})} - \boldsymbol{\mu}^{(\boldsymbol{w}_{i})} \right)^{\top} \left( \boldsymbol{\Sigma}^{(\boldsymbol{w}_{i})} \right)^{-1} \left( \boldsymbol{z}_{i}^{(\boldsymbol{w}_{i})} - \boldsymbol{\mu}^{(\boldsymbol{w}_{i})} \right).\] The corresponding negative log-likelihood is \[\ell \left(\boldsymbol{\mu}^{(\boldsymbol{w}_{i})},\boldsymbol{\Sigma}^{(\boldsymbol{w}_{i})}; \boldsymbol{z}_{i}^{(\boldsymbol{w}_{i})}\right) = \ln |\boldsymbol{\Sigma}^{(\boldsymbol{w}_{i})}| + q^{(\boldsymbol{w}_{i})}\ln(2\pi) + \mathrm{MD}^{2} \left( \boldsymbol{z}_{i}^{(\boldsymbol{w}_{i})}, \boldsymbol{\mu}^{(\boldsymbol{w}_{i})},\boldsymbol{\Sigma}^{(\boldsymbol{w}_{i})} \right).\]

Unlike ordinary missing-data problems, the cellwise indicator matrix \(\boldsymbol{W}\) is not observed in advance but is estimated from the data. To avoid flagging too many entries, the cellwise MCD imposes a lower bound on the number of retained cells in each column and adds a penalty for flagged cells. Specifically, for a given \(h\), we require \(\|\boldsymbol{W}_{\cdot j}\|_0 \geq h, j=1,2,\ldots,q,\) where \(\|\boldsymbol{W}_{\cdot j}\|_0\) is the number of retained entries in the \(j\)-th column. We also impose \(\lambda_{\min}(\boldsymbol{\Sigma}) \geq a,\) where \(\lambda_{\min}(\boldsymbol{\Sigma})\) denotes the smallest eigenvalue of \(\boldsymbol{\Sigma}\), for some \(a>0\), thereby ensuring that the covariance matrix is nonsingular. The cellwise MCD estimator is then defined as the solution to \[\label{eq:cellwise95mcd95obj} \begin{align} \min_{\boldsymbol{\mu},\boldsymbol{\Sigma},\boldsymbol{W}} \quad & \sum_{i=1}^{n} \ell \left(\boldsymbol{\mu}^{(\boldsymbol{w}_{i})},\boldsymbol{\Sigma}^{(\boldsymbol{w}_{i})};\boldsymbol{z}_{i}^{(\boldsymbol{w}_{i})}\right) + \sum_{j=1}^{q} b_j \|\mathbf{1}_{n}-\boldsymbol{W}_{\cdot j}\|_0 \\ \mathrm{s.t.} \quad & \boldsymbol{W}\in\{0,1\}^{n\times q}, \quad \|\boldsymbol{W}_{\cdot j}\|_0 \geq h,\quad j=1,2,\ldots,q,\quad \lambda_{\min}(\boldsymbol{\Sigma})\geq a . \end{align}\tag{2}\] Here \(\mathbf{1}_{n}\) denotes the \(n\)-dimensional all-ones vector, \(\|\mathbf{1}_{n}-\boldsymbol{W}_{\cdot j}\|_0\) counts the number of flagged cells in the \(j\)-th feature dimension, and \(b_j>0\) controls the penalty for flagging cells in that dimension. If an observation has no retained entries, i.e., \(q^{(\boldsymbol{w}_i)}=0\), we set \(\ln |\boldsymbol{\Sigma}^{(\boldsymbol{w}_i)}| = 0\) and \(\mathrm{MD}^2( \boldsymbol{z}_i^{(\boldsymbol{w}_i)}, \boldsymbol{\mu}^{(\boldsymbol{w}_i)}, \boldsymbol{\Sigma}^{(\boldsymbol{w}_i)} )=0\).

It is worth noting that directly solving 2 is computationally challenging, because the objective involves both continuous parameters \((\boldsymbol{\mu},\boldsymbol{\Sigma})\) and a high-dimensional binary matrix \(\boldsymbol{W}\). Following [@raymaekers2024cellwise], the cellwise MCD can be computed through an iterative procedure that alternates between updating the cellwise indicator matrix and estimating the mean vector and covariance matrix from the retained cells. The penalty parameters \(b_j\) regulate the degree of cellwise flagging and are typically chosen so that, in the absence of contamination, only a small expected fraction of cells is flagged. The parameter \(h\) specifies the minimum number of retained entries in each feature dimension and prevents degenerate solutions that remove too many cells. Thus, \(h\) and \(b_j\) jointly control the robustness–efficiency trade-off: more aggressive flagging improves robustness to cellwise contamination, whereas retaining more cells improves efficiency when contamination is mild. The resulting \(\widehat{\boldsymbol{\mu}}\) and \(\widehat{\boldsymbol{\Sigma}}\) are then used in the single-class MDS detector.

4 Multi-class robust estimation↩︎

In the multi-class MDS detector introduced in Section 2.2, positive samples may come from multiple classes. The covariance structures of different classes may exhibit both homogeneity and heterogeneity. Estimating the covariance matrix for each class separately may fail to borrow information across classes and thus lead to insufficient statistical efficiency. Therefore, building on the casewise and cellwise MCD principles, we extend these two robust estimators to the multi-class setting and propose corresponding joint estimation methods. For notational simplicity, we write the total number of classes as \(K\) in this section.

4.1 Joint casewise MCD estimator↩︎

We first consider casewise contamination, where entire observations may be unreliable within each class. Let \(\boldsymbol{Z}=\{\boldsymbol{Z}^{(1)},\ldots,\boldsymbol{Z}^{(K)}\}\), where \(\boldsymbol{Z}^{(k)}=(\boldsymbol{z}_{1}^{(k)},\ldots,\boldsymbol{z}_{n_k}^{(k)})^\top \in \mathbb{R}^{n_k\times q}\) denotes the representation matrix of the \(k\)-th positive class, \(k=1,2,\ldots,K\). As a working model, we assume \(\boldsymbol{z}_{i}^{(k)} \sim \mathcal{N}_{q} (\boldsymbol{\mu}^{(k)}, \boldsymbol{\Sigma}^{(k)}), i=1,2,\ldots,n_k, k=1,2,\ldots,K,\) where \(\boldsymbol{\mu}^{(k)}\in\mathbb{R}^{q}\) and \(\boldsymbol{\Sigma}^{(k)}\in\mathbb{R}^{q\times q}\) denote the class-specific mean vector and covariance matrix. Following the reparameterization idea in the joint estimation of multiple graphical models [@guo2011joint], we introduce the following decomposition to explicitly model both shared and class-specific dependence structures: \[\label{eq:sigma-decomposition} \boldsymbol{\Sigma}^{(k)} = \boldsymbol{\Theta} \odot \boldsymbol{\Gamma}^{(k)},\tag{3}\] where \(\odot\) denotes the Hadamard product. Equivalently, for \(\boldsymbol{\Sigma}^{(k)}=(\sigma_{ij}^{(k)})_{1\leq i,j \leq q}\), \[\sigma_{ij}^{(k)} = \begin{cases} \theta_{ij}\gamma_{ij}^{(k)}, & i\neq j,\\ \gamma_{ii}^{(k)}, & i=j, \end{cases}\] with \(\theta_{ii}=1\) for \(i=1,2,\ldots,q\). Here \(\boldsymbol{\Theta}=(\theta_{ij})_{1\leq i,j \leq q}\) represents the shared covariance pattern across classes, whereas \(\boldsymbol{\Gamma}^{(k)}=(\gamma_{ij}^{(k)})_{1\leq i,j \leq q}\) captures class-specific variation. This reparameterization separates shared dependence (\(\boldsymbol{\Theta}\)) from class-specific deviations (\(\boldsymbol{\Gamma}^{(k)}\)). We use the Root Mean Square (RMS) correlation normalization rule to ensure the identifiability of the decomposition, and details are given in Appendix 9.2.

For \(k=1,2,\ldots,K\), let \(w_i^{(k)}\in\{0,1\}\) indicate whether the \(i\)-th observation in class \(k\) is retained. We write \(\boldsymbol{H}^{(k)} = (w_1^{(k)},\ldots,w_{n_k}^{(k)})^\top\) for the vector of casewise weights in class \(k\) with \(\sum_{i=1}^{n_k}w_i^{(k)}=h_k.\) The joint casewise MCD estimator is defined as the solution to \[\label{eq:joint95casewise95mcd} \begin{align} \min_{\{\boldsymbol{\mu}^{(k)}, \boldsymbol{H}^{(k)}, \boldsymbol{\Sigma}^{(k)}\}_{k=1}^{K}} \quad & \sum_{k=1}^{K} \sum_{i=1}^{n_k} w_i^{(k)} \tilde{\ell}(\boldsymbol{\mu}^{(k)}, \boldsymbol{\Sigma}^{(k)};\boldsymbol{z}_i^{(k)}) + \mathbf{P}(\boldsymbol{\Theta}, \{\boldsymbol{\Gamma}^{(k)}\}_{k=1}^{K}) \\ \mathrm{s.t.} \quad & w_i^{(k)}\in\{0,1\}, \quad \sum_{i=1}^{n_k}w_i^{(k)}=h_k, \quad \boldsymbol{\Sigma}^{(k)}=\boldsymbol{\Theta}\odot\boldsymbol{\Gamma}^{(k)}, \quad \lambda_{\min}(\boldsymbol{\Sigma}^{(k)})\geq a_k, \end{align}\tag{4}\] where \[\tilde{\ell}({\boldsymbol{\mu}}^{(k)}, \boldsymbol{\Sigma}^{(k)};\boldsymbol{z}_i^{(k)}) = \ln|\boldsymbol{\Sigma}^{(k)}| + q\ln(2\pi) + \mathrm{MD}^2 \left( \boldsymbol{z}_{i}^{(k)}, \boldsymbol{\mu}^{(k)}, \boldsymbol{\Sigma}^{(k)} \right)\] with \(\mathrm{MD}^2 \left( \boldsymbol{z}_{i}^{(k)}, \boldsymbol{\mu}^{(k)}, \boldsymbol{\Sigma}^{(k)} \right) = (\boldsymbol{z}_{i}^{(k)}-\boldsymbol{\mu}^{(k)})^\top (\boldsymbol{\Sigma}^{(k)})^{-1} (\boldsymbol{z}_{i}^{(k)}-\boldsymbol{\mu}^{(k)})\), \(k=1,2,\ldots,K\). Following [@guo2011joint], we impose the penalty \[\label{eq:joint95casewise95penalty} \mathbf{P}(\boldsymbol{\Theta}, \{\boldsymbol{\Gamma}^{(k)}\}_{k=1}^{K}) = \lambda_1 \|\boldsymbol{\Theta}\|_{1,\mathrm{off}} + \lambda_2 \sum_{k=1}^{K} \|\boldsymbol{\Gamma}^{(k)}\|_{1,\mathrm{off}} := \lambda_1 \sum_{i\neq j} |\theta_{ij}| + \lambda_2 \sum_{k=1}^{K} \sum_{i\neq j} |\gamma_{ij}^{(k)}|,\tag{5}\] where \(\lambda_1,\lambda_2>0\) control the global and class-specific sparsity levels, respectively. Here, \(\|\cdot\|_{1, \mathrm{off}}\) denotes the sum of the absolute values of the off-diagonal elements of a matrix. The first term promotes a sparse shared structure through \(\boldsymbol{\Theta}\): when \(\theta_{ij}\) is set to zero, the corresponding covariance entry is removed from all classes simultaneously, encouraging homogeneity across classes. The second term regularizes the class-specific components \(\boldsymbol{\Gamma}^{(k)}\), allowing different classes to retain distinct covariance magnitudes and signs when the shared component is nonzero, thereby preserving heterogeneity. Therefore, the penalty jointly encourages common structure across classes and flexible class-specific deviations. The constraint \(\lambda_{\min}(\boldsymbol{\Sigma}^{(k)})\geq a_k > 0\) ensures that each covariance matrix is nonsingular.

4.2 Joint cellwise MCD estimator↩︎

We next consider cellwise contamination, where only a subset of feature entries may be unreliable. In this setting, removing entire observations can be inefficient, since many uncontaminated entries may still contain useful information. The joint cellwise MCD estimator therefore assigns binary weights at the cell level within each class. For class \(k\), let \(\boldsymbol{W}^{(k)}=(w_{ij}^{(k)})_{1\leq i \leq n_k, 1\leq j \leq q}\in\{0,1\}^{n_k\times q}\) be the cellwise indicator matrix, where \(w_{ij}^{(k)}=1\) means that the \(j\)-th entry of \(\boldsymbol{z}_{i}^{(k)}\) is retained, and \(w_{ij}^{(k)}=0\) means that the entry is flagged as contaminated. Let \(\boldsymbol{w}_{i}^{(k)}\) denote the \(i\)-th row of \(\boldsymbol{W}^{(k)}\), and let \(q^{(\boldsymbol{w}_{i}^{(k)})} = \sum_{j=1}^{q} w_{ij}^{(k)}\) be the number of retained entries in observation \(i\) of class \(k\). For any vector or matrix, the superscript \((\boldsymbol{w}_{i}^{(k)})\) denotes the subvector or principal submatrix corresponding to the retained entries. Thus, \(\boldsymbol{z}_{i}^{(\boldsymbol{w}_{i}^{(k)}, k)}, \boldsymbol{\mu}^{(\boldsymbol{w}_{i}^{(k)}, k)}, \boldsymbol{\Sigma}^{(\boldsymbol{w}_{i}^{(k)}, k)}\) denote the retained subvector of \(\boldsymbol{z}_{i}^{(k)}\), the corresponding subvector of \(\boldsymbol{\mu}^{(k)}\), and the corresponding principal submatrix of \(\boldsymbol{\Sigma}^{(k)}\), respectively. Using the same covariance decomposition in 3 , the joint cellwise MCD estimator is defined by minimizing \[\label{eq:joint95cellwise95mcd} \begin{align} \min_{\{\boldsymbol{\mu}^{(k)}, \boldsymbol{W}^{(k)}, \boldsymbol{\Sigma}^{(k)}\}_{k=1}^{K}} \quad & \sum_{k=1}^{K} \sum_{i=1}^{n_k} \tilde{\ell} \left( \boldsymbol{\mu}^{(\boldsymbol{w}_{i}^{(k)}, k)}, \boldsymbol{\Sigma}^{(\boldsymbol{w}_{i}^{(k)}, k)}; \boldsymbol{z}_{i}^{(\boldsymbol{w}_{i}^{(k)}, k)} \right) \\ &\quad + \sum_{k=1}^{K} \sum_{j=1}^{q} b_j^{(k)} \left\| \mathbf{1}_{n_k} - \boldsymbol{W}_{\cdot j}^{(k)} \right\|_0 + \mathbf{P}(\boldsymbol{\Theta}, \{\boldsymbol{\Gamma}^{(k)}\}_{k=1}^{K}) \\ \mathrm{s.t.} \quad & \boldsymbol{W}^{(k)}\in\{0,1\}^{n_k\times q}, \quad \left\|\boldsymbol{W}_{\cdot j}^{(k)}\right\|_0\geq h_k, \quad j=1,2,\ldots,q, \\ & \boldsymbol{\Sigma}^{(k)}=\boldsymbol{\Theta}\odot\boldsymbol{\Gamma}^{(k)}, \quad \lambda_{\min}(\boldsymbol{\Sigma}^{(k)})\geq a_k > 0, \quad k=1,2,\ldots,K. \end{align}\tag{6}\] Let \(\mathrm{MD}^2 (\boldsymbol{z}_{i}^{(\boldsymbol{w}_{i}^{(k)}, k)},\boldsymbol{\mu}^{(\boldsymbol{w}_{i}^{(k)}, k)},\boldsymbol{\Sigma}^{(\boldsymbol{w}_{i}^{(k)}, k)}) = (\boldsymbol{z}_{i}^{(\boldsymbol{w}_{i}^{(k)},k)}-\boldsymbol{\mu}^{(\boldsymbol{w}_{i}^{(k)},k)})^{\top} (\boldsymbol{\Sigma}^{(\boldsymbol{w}_{i}^{(k)},k)})^{-1} (\boldsymbol{z}_{i}^{(\boldsymbol{w}_{i}^{(k)},k)}-\boldsymbol{\mu}^{(\boldsymbol{w}_{i}^{(k)},k)})\) be the partial Mahalanobis distance, we have \(\tilde{\ell}(\boldsymbol{\mu}^{(\boldsymbol{w}_{i}^{(k)}, k)},\boldsymbol{\Sigma}^{(\boldsymbol{w}_{i}^{(k)}, k)};\boldsymbol{z}_{i}^{(\boldsymbol{w}_{i}^{(k)}, k)}) = \ln |\boldsymbol{\Sigma}^{(\boldsymbol{w}_{i}^{(k)},k)}| + q^{(\boldsymbol{w}_{i}^{(k)})}\ln(2\pi) + \mathrm{MD}^2 (\boldsymbol{z}_{i}^{(\boldsymbol{w}_{i}^{(k)}, k)},\boldsymbol{\mu}^{(\boldsymbol{w}_{i}^{(k)}, k)},\boldsymbol{\Sigma}^{(\boldsymbol{w}_{i}^{(k)}, k)})\). Here, \(\mathbf{1}_{n_k}\) is the \(n_k\)-dimensional all-ones vector, and \(\| \mathbf{1}_{n_k} - \boldsymbol{W}_{\cdot j}^{(k)} \|_0\) counts the number of flagged cells in the \(j\)-th feature dimension of class \(k\). The constants \(b_j^{(k)}\) penalize cellwise flagging, while the constraints \(\|\boldsymbol{W}_{\cdot j}^{(k)}\|_0\geq h_k\) prevent too many entries from being removed in any feature dimension. If all entries of an observation are flagged, we adopt the empty-pattern convention \(q^{(\boldsymbol{w}_i^{(k)})}=0\), \(\ln |\boldsymbol{\Sigma}^{(\boldsymbol{w}_i^{(k)},k)}|=0\), and \(\mathrm{MD}^2=0\), so that the objective function remains well-defined. \(\lambda_{\min}(\boldsymbol{\Sigma}^{(k)})\geq a_k > 0\) guarantees that the covariance matrix is nonsingular. The penalty \(\mathbf{P}(\boldsymbol{\Theta}, \{\boldsymbol{\Gamma}^{(k)}\}_{k=1}^{K})\) is the same as in 5 , encouraging homogeneity through the shared component \(\boldsymbol{\Theta}\) while preserving heterogeneity by the class-specific components \(\boldsymbol{\Gamma}^{(k)}\).

5 Optimization algorithms↩︎

In this section, we develop optimization algorithms for the two joint estimators introduced in Section 4. The main difficulty is that both estimators involve discrete selection variables and continuous covariance parameters under the shared-structure reparameterization \(\boldsymbol{\Sigma}^{(k)}=\boldsymbol{\Theta}\odot\boldsymbol{\Gamma}^{(k)}\). For the joint casewise MCD estimator, the discrete variables correspond to the retained observations in each class. For the joint cellwise MCD estimator, the discrete variables correspond to retained cells in each class-specific representation matrix. We solve both problems by alternating between updating the discrete selection variables and updating the continuous parameters. We use \(t \in \mathbb{N}\) as a superscript on the estimator to represent the result obtained in the \(t\)-th iteration, where \(\mathbb{N}\) represents the set of non-negative integers, and \(t=0\) represents the initial value.

5.1 Optimization for joint casewise MCD estimator↩︎

Under casewise contamination, we first initialize, for each group \(k\), a subset \(\boldsymbol{H}^{(k,0)}\) of size \(h_k\). Starting from this initial subset, the algorithm iteratively updates both the model parameters and the subset composition until convergence. We first update the group means based on the current subset: \[\widehat{\boldsymbol{\mu}}^{(k,t+1)} = \frac{\sum_{i=1}^{n_k} w_i^{(k,t)} \boldsymbol{z}_i^{(k)}}{\sum_{i=1}^{n_k} w_i^{(k,t)}} = \frac{\sum_{i=1}^{n_k} w_i^{(k,t)} \boldsymbol{z}_i^{(k)}}{h_k}, \quad k = 1,2,\ldots, K.\]

We next update the shared structure \(\widehat{\boldsymbol{\Theta}}^{(t+1)}\) and the group-specific heterogeneity components \(\widehat{\boldsymbol{\Gamma}}^{(k,t+1)}\) via alternating optimization:

• We first fix \(\widehat{\boldsymbol{\Theta}}^{(t)}\) and update each \(\widehat{\boldsymbol{\Gamma}}^{(k,t+1)}\): \[\widehat{\boldsymbol{\Gamma}}^{(k,t+1)} = \arg\min_{\boldsymbol{\Gamma}} \sum_{i=1}^{n_k} w_i^{(k,t)} \tilde{\ell}(\widehat{\boldsymbol{\mu}}^{(k,t+1)}, \widehat{\boldsymbol{\Theta}}^{(t)} \odot \boldsymbol{\Gamma}; \boldsymbol{z}_i^{(k)}) + \lambda_2 \sum_{i \neq j} |\gamma_{ij}|.\]

• We then fix all \(\widehat{\boldsymbol{\Gamma}}^{(k,t+1)}\) and update the shared structure \(\widehat{\boldsymbol{\Theta}}^{(t+1)}\): \[\widehat{\boldsymbol{\Theta}}^{(t+1)} = \arg\min_{\boldsymbol{\Theta}} \sum_{k=1}^{K} \sum_{i=1}^{n_k} w_i^{(k,t)} \tilde{\ell}(\widehat{\boldsymbol{\mu}}^{(k,t+1)}, \boldsymbol{\Theta} \odot \widehat{\boldsymbol{\Gamma}}^{(k,t+1)}; \boldsymbol{z}_i^{(k)}) + \lambda_1 \sum_{i \neq j} |\theta_{ij}|.\]

The covariance matrices are updated by \(\widehat{\boldsymbol{\Sigma}}^{(k, t+1)} = \widehat{\boldsymbol{\Theta}}^{(t+1)} \odot \widehat{\boldsymbol{\Gamma}}^{(k,t+1)}, k = 1,2,\ldots, K\). After updating the parameters, the algorithm computes the Mahalanobis distances of all samples to their group centers using the updated mean and covariance estimates. Within each group, the \(h_k\) samples with the smallest Mahalanobis distances are selected as unflagged samples, forming the updated subset \(\boldsymbol{H}^{(k,t+1)}\). This process is repeated until convergence. To enhance robustness, a reweighting step is performed after convergence, following the Fast-MCD algorithm [@rousseeuw1999fast]. The subset is updated based on Mahalanobis distances, and the final parameter estimates are computed with respect to this reweighted subset. The complete computation procedure is presented in Algorithm 2.

5.2 Optimization for joint cellwise MCD estimator↩︎

For the cellwise setting, it is difficult to directly solve Equation 6 . The optimization proceeds by first estimating the weight matrices \(\boldsymbol{W}\) and then updating the parameters given the missing entries. For ease of discussion, we introduce notation and partitions. In the \(t\)-th iteration, for the \(k\)-th group, the \(i\)-th row is divided into three parts \(\boldsymbol{z}_i^{(k)} = (z_{ij}^{(k)}, \boldsymbol{z}_{i,o}^{(k)\top}, \boldsymbol{z}_{i,-o}^{(k)\top})^{\top}\) after a suitable permutation of coordinates, where \(z_{ij}^{(k)}\) is the element in row \(i\) and column \(j\), and \(\boldsymbol{z}_{i,o}^{(k)}\) and \(\boldsymbol{z}_{i,-o}^{(k)}\) correspond to entries marked as normal \((w_{i\tau}^{(k,t)} = 1)\) and abnormal \((w_{i\tau}^{(k,t)} = 0)\), \(\tau=1,2,\ldots,q\), respectively, with \(o = \{\tau \neq j: w_{i\tau}^{(k,t)} = 1\}\). Accordingly, the mean vector \(\widehat{\boldsymbol{\mu}}^{(k)}\) and covariance matrix \(\widehat{\boldsymbol{\Sigma}}^{(k)}\) are partitioned as \[\widehat{\boldsymbol{\mu}}^{(k)} = \begin{pmatrix} \mu_j^{(k)} \\ \widehat{\boldsymbol{\mu}}_o^{(k)} \\ \widehat{\boldsymbol{\mu}}_{-o}^{(k)} \end{pmatrix}, \quad \widehat{\boldsymbol{\Sigma}}^{(k)} = \begin{pmatrix} \widehat{\sigma}_{jj}^{(k)} & \widehat{\boldsymbol{\Sigma}}_{j,o}^{(k)} & \widehat{\boldsymbol{\Sigma}}_{j,-o}^{(k)} \\ \widehat{\boldsymbol{\Sigma}}_{o,j}^{(k)} & \widehat{\boldsymbol{\Sigma}}_{o,o}^{(k)} & \widehat{\boldsymbol{\Sigma}}_{o,-o}^{(k)} \\ \widehat{\boldsymbol{\Sigma}}_{-o,j}^{(k)} & \widehat{\boldsymbol{\Sigma}}_{-o,o}^{(k)} & \widehat{\boldsymbol{\Sigma}}_{-o,-o}^{(k)} \end{pmatrix}.\]

We first initialize each group separately by computing initial location estimates \(\widehat{\boldsymbol{\mu}}^{(k,0)}\), covariance estimates \(\widehat{\boldsymbol{\Sigma}}^{(k,0)}\), and weight matrices \(\boldsymbol{W}^{(k,0)}\), ensuring that each column contains at least \(h_k\) non-missing entries not flagged as outliers. An iterative process then alternates between updating cellwise weights and updating parameters. In the \(t\)-th iteration, with \(\widehat{\boldsymbol{\mu}}^{(k,t)}\) and \(\widehat{\boldsymbol{\Sigma}}^{(k,t)}\) fixed, the objective function used for updating the weight of entry \((i,j)\) in group \(k\) is \[\label{eq:W95delta} \Delta_{ij}^{(k)} = \ln(C_{ij}^{(k)}) + \ln(2\pi) + \frac{(z_{ij}^{(k)} - \widehat{z}_{ij}^{(k)})^2}{C_{ij}^{(k)}} - b_j^{(k)},\tag{7}\] where \(b_j^{(k)} = \max \{ \chi_{1-\alpha}^2(1) + \ln (2\pi) + \ln (1/([(\widehat{\boldsymbol{\Sigma}}^{(k,0)})^{-1}]_{jj}), 0\}\), with \(\chi_{1-\alpha}^2(1)\) representing the chi-square \(1-\alpha\) quantile with 1 degree of freedom. Moreover, the conditional mean and variance are \[\label{eq:cond95mean95var95W} \widehat{z}_{ij}^{(k)} = \widehat{\mu}_j^{(k,t)} + \widehat{\boldsymbol{\Sigma}}_{j,o}^{(k,t)} (\widehat{\boldsymbol{\Sigma}}_{o,o}^{(k,t)})^{-1} (\boldsymbol{z}_{i,o}^{(k)} - \widehat{\boldsymbol{\mu}}_o^{(k,t)}), \quad C_{ij}^{(k)} = \widehat{\sigma}_{jj}^{(k,t)} - \widehat{\boldsymbol{\Sigma}}_{j,o}^{(k,t)} (\widehat{\boldsymbol{\Sigma}}_{o,o}^{(k,t)})^{-1} (\widehat{\boldsymbol{\Sigma}}_{j,o}^{(k,t)})^\top.\tag{8}\] The index set \(o\) is determined by the current weight matrix \(\boldsymbol{W}^{(k,t)}\). All cells with \(\Delta_{ij}^{(k)}<0\) are retained. If fewer than \(h_k\) cells are retained in the \(j\)-th column, the smallest remaining values of \(\Delta_{ij}^{(k)}\) are additionally retained until \(h_k\) cells are retained [@raymaekers2024cellwise]. Then, the Expectation-Maximization (EM) algorithm [@dempster1977maximum] updates \(\widehat{\boldsymbol{\mu}}^{(k,t+1)}\), \(\widehat{\boldsymbol{\Gamma}}^{(k,t+1)}\), and \(\widehat{\boldsymbol{\Theta}}^{(t+1)}\).

E-step. Given \(\boldsymbol{W}^{(k,t+1)}\), the conditional expectation and covariance for the \(i\)-th observation in group \(k\) are computed. Let \(\mathcal{O}_i = \{\tau: w_{i\tau}^{(k,t+1)}=1\}\), \(\mathcal{U}_i = \{\tau: w_{i\tau}^{(k,t+1)}=0\}\), and let \(\boldsymbol{z}_{i, \mathcal{O}_i}^{(k)}\) be a subvector indexed by \(\mathcal{O}_i\). The conditional expectation is \(\boldsymbol{m}_i^{(k,t)} = ({\boldsymbol{m}_{i,\mathcal{O}_i}^{(k,t)}}^{\top}, {\boldsymbol{m}_{i,\mathcal{U}_i}^{(k,t)}}^{\top})^{\top}\), where \[\label{eq:cond95mean} \boldsymbol{m}_{i,\mathcal{O}_i}^{(k,t)} = \boldsymbol{z}_{i,\mathcal{O}_i}^{(k)}, \quad \boldsymbol{m}_{i,\mathcal{U}_i}^{(k,t)} = \widehat{\boldsymbol{\mu}}_{\mathcal{U}_i}^{(k,t)} + \widehat{\boldsymbol{\Sigma}}_{\mathcal{U}_i,\mathcal{O}_i}^{(k,t)} (\widehat{\boldsymbol{\Sigma}}_{\mathcal{O}_i,\mathcal{O}_i}^{(k,t)})^{-1} (\boldsymbol{z}_{i,\mathcal{O}_i}^{(k)} - \widehat{\boldsymbol{\mu}}_{\mathcal{O}_i}^{(k,t)}),\tag{9}\] and the conditional covariance \(V_{i}^{(k,t)}\) satisfies \[\label{eq:cond95var} \begin{align} V_{i,{\mathcal{O}_i,\mathcal{O}_i}}^{(k,t)} = \mathbf{O},\quad V_{i,{\mathcal{O}_i,\mathcal{U}_i}}^{(k,t)} = V_{i,{\mathcal{U}_i,\mathcal{O}_i}}^{(k,t)} = \mathbf{0},\quad V_{i,{\mathcal{U}_i,\mathcal{U}_i}}^{(k,t)} = \widehat{\boldsymbol{\Sigma}}_{\mathcal{U}_i,\mathcal{U}_i}^{(k,t)} - \widehat{\boldsymbol{\Sigma}}_{\mathcal{U}_i,\mathcal{O}_i}^{(k,t)} (\widehat{\boldsymbol{\Sigma}}_{\mathcal{O}_i,\mathcal{O}_i}^{(k,t)})^{-1} \widehat{\boldsymbol{\Sigma}}_{\mathcal{O}_i,\mathcal{U}_i}^{(k,t)}, \end{align}\tag{10}\] where \(\mathbf{O}\) and \(\mathbf{0}\) represent a matrix and a vector containing only zeros, respectively. When \(\mathcal{O}_i=\varnothing\), we use the convention \(\boldsymbol{m}_i^{(k,t)}=\widehat{\boldsymbol{\mu}}^{(k,t)}\) and \(V_i^{(k,t)}=\widehat{\boldsymbol{\Sigma}}^{(k,t)}\). The conditional second-order central moment is then \[\mathbb{E}\left[ (\boldsymbol{z}_i^{(k)}-\boldsymbol{\mu}^{(k)}) (\boldsymbol{z}_i^{(k)}-\boldsymbol{\mu}^{(k)})^\top \,\middle\vert\, \boldsymbol{z}_{i,\mathcal{O}_i}^{(k)} \right] = V_i^{(k,t)} + (\boldsymbol{m}_i^{(k,t)}-\boldsymbol{\mu}^{(k)}) (\boldsymbol{m}_i^{(k,t)}-\boldsymbol{\mu}^{(k)})^\top .\]

M-step. We first update the location parameters \(\widehat{\boldsymbol{\mu}}^{(k,t+1)}\): \(\widehat{\boldsymbol{\mu}}^{(k,t+1)} = \frac{1}{n_k} \sum_{i=1}^{n_k} \boldsymbol{m}_i^{(k,t)},\) and compute \[S_{\mathrm{cell}}^{(k,t)} = \frac{1}{n_k} \sum_{i=1}^{n_k} \left[ V_i^{(k,t)} + (\boldsymbol{m}_i^{(k,t)} - \widehat{\boldsymbol{\mu}}^{(k,t+1)})(\boldsymbol{m}_i^{(k,t)} - \widehat{\boldsymbol{\mu}}^{(k,t+1)})^\top \right].\] Next, fixing \(\widehat{\boldsymbol{\Theta}}^{(t)}\), the group-specific structure \(\widehat{\boldsymbol{\Gamma}}^{(k,t+1)}\) is updated by solving \[\widehat{\boldsymbol{\Gamma}}^{(k,t+1)} = \arg\min_{\boldsymbol{\Gamma}} n_k \left[ \ln |\widehat{\boldsymbol{\Theta}}^{(t)} \odot \boldsymbol{\Gamma}| + \mathrm{tr}((\widehat{\boldsymbol{\Theta}}^{(t)} \odot \boldsymbol{\Gamma})^{-1} S_{\mathrm{cell}}^{(k,t)}) \right] + \lambda_2 \sum_{i \neq j} |\gamma_{ij}|.\] Finally, the shared structure \(\widehat{\boldsymbol{\Theta}}^{(t+1)}\) is updated by \[\widehat{\boldsymbol{\Theta}}^{(t+1)} = \arg\min_{\boldsymbol{\Theta}} \sum_{k=1}^{K} n_k \left[ \ln |\boldsymbol{\Theta} \odot \widehat{\boldsymbol{\Gamma}}^{(k,t+1)}| + \mathrm{tr}((\boldsymbol{\Theta} \odot \widehat{\boldsymbol{\Gamma}}^{(k,t+1)})^{-1} S_{\mathrm{cell}}^{(k,t)}) \right] + \lambda_1 \sum_{i\neq j} |\theta_{ij}|.\]

The covariance matrices for each group are then reconstructed via \(\widehat{\boldsymbol{\Sigma}}^{(k,t+1)} = \widehat{\boldsymbol{\Theta}}^{(t+1)} \odot \widehat{\boldsymbol{\Gamma}}^{(k,t+1)},\) and the algorithm iterates until convergence. The full procedure is summarized in Algorithm 3.

The updates for \(\widehat{\boldsymbol{\Theta}}\) and \(\{\widehat{\boldsymbol{\Gamma}}^{(k)}\}_{k=1}^{K}\), which correspond to the Fast-MCD algorithm for the joint casewise MCD estimator in Section 5.1 and the EM-based procedure for the joint cellwise MCD estimator in Section 5.2, are essentially constrained optimization problems. We use proximal gradient descent with backtracking search to obtain approximate solutions. Detailed update formulas are provided in 9.4. Details such as algorithm initialization, the identifiability of the covariance decomposition, and dimensionality reduction in the implementation are also given in 9.

6 Theoretical properties↩︎

In this section, we discuss the main theoretical properties of the proposed joint estimators. Specifically, the joint casewise MCD estimator is optimized via the Fast-MCD algorithm, and the joint cellwise MCD estimator is optimized via the EM-based procedure. Detailed statements and proofs of these properties are provided in the subsequent theorems and 10.

6.1 Convergence guarantees↩︎

The following theorems establish convergence guarantees for the proposed joint estimators. In particular, they show that the algorithms generate non-increasing sequences of their respective objective values and that these objective values converge to finite limits. Moreover, any interior accumulation point satisfies suitable blockwise stationarity conditions for the continuous variables, together with optimality or coordinatewise optimality of the corresponding discrete selection variables.

Theorem 1 ensures that the Fast-MCD algorithm for the joint casewise MCD estimator, before the final reweighting step, produces a non-increasing sequence of the regularized negative log-likelihood values. This implies that each pre-reweighting iteration either decreases or maintains the objective value, and that the objective values converge to a finite limit. Moreover, any interior accumulation point is blockwise stationary with respect to the continuous parameters and is optimal with respect to the discrete casewise selection block. Consequently, practitioners can rely on the algorithm to produce stable estimates for casewise-contaminated data.

Theorem 2 establishes analogous convergence guarantees for the joint cellwise MCD estimator using the EM-based generalized procedure. The cellwise objective values decrease monotonically and converge to a finite limit. Moreover, any interior accumulation point satisfies the stationarity conditions for the continuous parameters with the cellwise weights fixed, while the limiting cellwise weight matrix is coordinatewise optimal for the discrete weight update. These results provide a theoretical justification for the stability of the proposed procedure under cellwise contamination and missing-entry patterns.

6.2 Finite-sample breakdown property↩︎

To provide context, the following theorems characterize the robustness of the proposed joint estimators under casewise and cellwise contamination. We establish finite-sample breakdown points, showing how much contamination the estimators can tolerate while maintaining bounded estimates. In order to discuss the breakdown point properties of the joint estimators, we introduce some definitions under the joint estimation framework. Let \(\boldsymbol{Z}^{\star}=\{\boldsymbol{Z}^{\star(1)},\ldots,\boldsymbol{Z}^{\star(K)}\}\), where \(\boldsymbol{Z}^{\star(k)}=(\boldsymbol{z}_{1}^{\star(k)},\ldots,\boldsymbol{z}_{n_k}^{\star(k)})^\top \in \mathbb{R}^{n_k\times q}\) denotes the clean representation matrix of the \(k\)-th positive class, \(k=1,2,\ldots,K\). Let \(\widehat{\boldsymbol{\Omega}} = (\widehat{\boldsymbol{\mu}}, \widehat{\boldsymbol{\Sigma}}) = (\{\widehat{\boldsymbol{\mu}}^{(k)}\}_{k=1}^K,\{\widehat{\boldsymbol{\Sigma}}^{(k)}\}_{k=1}^K)\) denote the joint estimator. For any vector \(\mathbf{r} = (r_1,r_2,\ldots, r_K)\), define \[\beta(\mathbf{r}) := \max_{1 \le k \le K} \frac{r_k}{n_k}.\]

The above definition is appropriate because, within the joint estimation framework, the estimator comprises the location and scatter parameters of multiple positive classes. Unbounded shift, covariance explosion, or covariance degradation in the estimation of any single class will lead to the failure of the overall joint estimation. Therefore, we measure the joint contamination intensity by the maximum contamination ratio across all classes and define breakdown based on the worst-case scenario. When \(K=1\), these definitions reduce to the classic casewise and cellwise breakdown definitions, maintaining consistency with the existing theoretical frameworks of MCD and cellMCD. Based on these definitions, we establish the breakdown point properties of the joint estimators. We assume the dataset is in general position: any \((q-1)\)-dimensional affine subspace contains at most \(q\) points.

Theorem 3 indicates that the joint cellwise MCD estimator is highly robust: it can tolerate a large proportion of contaminated cells without the location or covariance estimates becoming unbounded. The lower bound \(\beta_0\) is sharp as a common guarantee for simultaneous boundedness of the location and scatter estimates under the stated contamination construction.

Theorem 4 similarly establishes that the joint casewise MCD estimator maintains bounded estimates even under substantial casewise contamination. The derived lower bounds on \(\delta_{\mathrm{case},-}\), \(\delta_{\mathrm{case},+}\), and \(\delta_{\mathrm{case},\boldsymbol{\mu}}^*\) provide a clear, provable measure of the estimator’s finite-sample robustness.

7 Real data applications↩︎

We evaluate our unified MDS-based framework on four real-world AI-related tasks: LLM-generated text, watermark, hallucination, and adversarial example detection. Performance is assessed using comprehensive metrics, including area under the receiver operating characteristic curve (ROC AUC), area under the precision-recall curve (PR AUC), F1 score, specificity, accuracy, precision, and recall. Our implementation is available at https://github.com/Astringency/JointMCD.git.

7.1 LLM-generated text detection↩︎

LLM-generated text detection aims to distinguish human-written passages from LLM-generated ones. We use the Human ChatGPT Comparison Corpus (HC3), a benchmark dataset widely adopted for AI-generated text detection, which contains paired human and ChatGPT-3.5 responses across multiple domains [@guo2023close; @ouyang2022training; @openai2022chatgpt]. We sample 1,000 instances from each domain, except for the “wiki_csai” subset, where all 842 available samples are used. To assess cross-generator generalization, we further collect responses to the same questions from five recent LLMs: ChatGPT-5.4 [@openai2026gpt54T; @openai2026gpt54], ChatGPT-5.4 mini [@openai2026gpt54mininano], Gemini-3 Flash Preview [@gemini2024gemini15; @google2026gemini3flashpreview], Gemini-3.1 Flash Lite Preview [@google2026gemini31flashlitepreview], and Gemini-3.1 Pro Preview [@google2026gemini31propreview]. This yields five HC3-based datasets, each with 4,842 instances. We apply the proposed unified detection framework using RoBERTa-base as the pretrained encoder [@liu2019roberta], and compare it with likelihood-based [@solaiman2019release], rank-based [@gehrmann2019gltr], entropy-based [@ippolito2020automatic], DetectGPT [@mitchell2023detectgpt], and RoBERTa baseline methods. For DetectGPT, we report both the raw perturbation discrepancy score, DetectGPT-d, and its normalized version, DetectGPT-z.

5pt

The results in Table ¿tbl:tab:detect95LGC? demonstrate the effectiveness and robustness of the proposed methods across different LLM generators: JCASEMCD achieves ROC AUC values of 0.969 and 0.951 on GPT-5.4 and GPT-5.4 mini, respectively, comparable to the strongest DetectGPT variants while yielding the highest specificity and precision, and remains competitive on the more challenging Gemini-based datasets, achieving the best overall performance on Gemini-3.1 Pro Preview in terms of ROC AUC, PR AUC, F1 score, and accuracy.

7.2 Watermark detection↩︎

Watermark detection can be viewed as a special case of LLM-generated text detection, where the generated text may contain deliberately embedded statistical traces. When the watermarking rule of the generating model is known, such traces can be detected by hypothesis-testing procedures specifically designed for the corresponding watermarking mechanism [@li2025statistical]. In contrast, the proposed framework does not require explicit knowledge of the watermarking rule, but instead detects distributional deviations in the space of deep representations.

We conduct watermark detection experiments based on Facebook OPT-1.3B [@zhang2022opt]. Three types of generated text are considered: non-watermarked text, denoted by Null; text generated with Gumbel-Max watermarking, denoted by Gumbel; and text generated with inverse-transform watermarking, denoted by Inverse. We consider two types of detection tasks: distinguishing human-written text from watermarked text, and distinguishing non-watermarked generated text from watermarked generated text. In these tasks, we denote unwatermarked text as positive samples and watermarked text as negative samples. The proposed method is compared with the hypothesis-testing framework of [@li2025statistical], which assumes access to the underlying watermarking mechanism. The results are reported in Table ¿tbl:tab:watermark95all95metrics?.

The results in Table ¿tbl:tab:watermark95all95metrics? show that the hypothesis-testing baselines of [@li2025statistical] perform strongly when the watermarking rule is known and matched to the test procedure, achieving the best F1 score, accuracy, precision, and specificity. This is expected because these methods are tailored to the known watermarking mechanism and directly exploit its statistical signal. In contrast, the proposed method is mechanism-agnostic and does not rely on explicit watermarking rules. Although its overall performance is generally lower, it provides a more balanced classification behavior: the testing-based methods attain very high precision and specificity but relatively low recall, indicating conservative detection of the negative class, namely watermarked text, whereas the proposed method improves recall while maintaining moderate precision, leading to a more balanced precision–recall trade-off.

This distinction becomes more important when the watermarking mechanism is unknown or the generated text is watermark-free. In this setting, we treat human-written text as positive samples and watermark-free LLM-generated text as negative samples. Since the testing framework of [@li2025statistical] depends on knowledge of the watermarking rule, its effectiveness may degrade substantially outside the matched watermark setting, as shown in Table 1. These results indicate that specialized testing procedures are preferable when the watermarking rule is fully known, whereas the proposed framework offers a more general, task-agnostic detection strategy.

7.3 Hallucination detection↩︎

Hallucination detection aims to determine whether generated text contains factually incorrect or internally inconsistent content. Compared with LLM-generated text detection, hallucination detection is more directly concerned with factual reliability and logical consistency. Hallucination detection can be conducted in either a reference-free or reference-based setting, depending on whether external evidence or reference answers are available. In both settings, correct facts are treated as positive samples and hallucinatory text as negative samples. We consider two benchmark datasets: HaluEval [@li2023halueval], a large-scale benchmark used for reference-free hallucination detection, and the biography-domain “wiki_bio_gpt3_hallucination” dataset [@manakul2023selfcheckgpt], which contains 1,908 sentence-level samples with Wikipedia evidence and hallucination labels, enabling reference-based evaluation.

As baselines, we use two DeBERTa-based NLI models, DeBERTa-v3-large-mnli-fever-anli-ling-wanli-binary and DeBERTa-v3-large-zeroshot-v1.1-all-33, denoted as DeBERTa-NLI and DeBERTa-ZS, respectively; both require reference information. The results in Table 2 show that the reference-free JCELLMCD achieves the best overall performance, with ROC AUC and PR AUC values of 0.751 and 0.734, respectively, outperforming the DeBERTa baselines by alleviating their low-specificity issue. JCELLMCD also consistently outperforms JCASEMCD, suggesting that the cellwise strategy is more effective in capturing localized hallucination anomalies.

7.4 Adversarial examples detection↩︎

Adversarial examples are deliberately perturbed inputs designed to mislead a model while remaining visually or semantically similar to benign samples. In this setting, normal images are positive samples, and adversarial images are negative samples. Since such perturbations can induce abnormal patterns in deep feature representations even when the input-space changes are small, adversarial example detection provides another natural setting for evaluating the proposed MDS-based framework. Existing detection methods include Local Intrinsic Dimensionality (LID), which detects adversarial samples by characterizing the intrinsic dimensionality of their local neighborhoods [@ma2018characterizing], and the Mahalanobis distance-based detector proposed by [@lee2018simple].

We study adversarial example detection in image classification tasks. On CIFAR-10 and CIFAR-100 [@krizhevsky2009learning], adversarial examples are generated using four representative attack methods: FGSM [@goodfellow2015explaining], BIM [@kurakin2018adversarial], CWL2 [@carlini2017towards], and DeepFool [@moosavi2016deepfool]. We also consider ImageNet [@deng2009imagenet] adversarial examples generated by Diff-PGD [@xue2023diffusion]. The target classifiers used to generate these adversarial examples are mainly ResNet [@he2016deep] and DenseNet [@huang2017densely] models. We use ResNet as the feature extractor to obtain deep representations for the proposed framework, and compare its performance with LID and the Mahalanobis distance-based method of [@lee2018simple]. The results are reported in Figure 6. The proposed framework matches or exceeds baselines in simpler settings (e.g., CIFAR-10 under FGSM/BIM) by effectively capturing adversarial deviations. However, it underperforms the standard Mahalanobis baseline in complex scenarios (e.g., CIFAR-100 under CWL2/DeepFool), likely due to the instability of joint covariance estimation in high-dimensional, multi-class spaces. Thus, while effective for moderate-complexity tasks, it requires further stabilization or dimensionality reduction for more challenging settings.

8 Conclusion↩︎

This paper proposed a unified framework for detecting AI-related anomalous samples, including LLM-generated text, hallucinations, watermarked text, and adversarial examples. The framework combined pretrained deep representations, robust estimation of the mean and covariance structure of positive samples, and MDSs for detection. To handle multi-class positive samples with shared and class-specific structures, we developed joint casewise and cellwise MCD estimators and established the convergence of the corresponding optimization algorithms and the high-breakdown-point properties of the estimators. Experiments demonstrated the effectiveness and generality of the proposed framework. These results showed that combining deep representations, robust joint covariance estimation, and MDS provided a flexible and statistically principled approach for detecting diverse AI-related anomalous samples. Future work will consider relaxing the normality assumption, improving high-dimensional covariance estimation, and extending the framework to online and distribution-shift settings.