arXiv CSopen access

A Unified Detection Framework for AI-Related Content and Artifacts

A Unified Detection Framework for AI-Related Content and Artifacts Xifeng Zhanga , Tao Hua , Yijie Pengb , Wan Tianc,d,∗ a

School of Mathematical Sciences, Capital Normal University, Beijing 100048, China School of Management and Engineering, Nanjing University, Nanjing 210008, China c Advanced Institute of Information Technology, Peking University d Wangxuan Institute of Computer Technology, Peking University, Beijing 100871, China b

arXiv:2607.07527v1 [stat.ML] 8 Jul 2026

Abstract Artificial intelligence (AI) is a double-edged sword: while it has achieved remarkable success across a wide range of domains, its deployment also calls for effective oversight and regulation, for which the detection of AI-related content and artifacts is perhaps the most direct and cost-effective approach. To this end, we propose a unified detection framework based on Mahalanobis distance scores (MDS), applicable to several important settings, including the detection of large language model (LLM) generated text, hallucination, watermark, and adversarial examples. A key component of the proposed method is to accurately characterize the positive class—such as human-generated text, factual statements, unwatermarked text, or non-adversarial samples—which requires an efficient and robust estimator of the covariance matrix of deep representations of positive samples before computing the MDS. Since the positive samples typically consist of multiple classes, and these classes may exhibit both homogeneity and heterogeneity, we develop joint estimation methods for both the casewise and cellwise minimum covariance determinant (MCD) estimators. We provide efficient optimization algorithms for both estimators and prove their convergence. We provide a reasonable definition of the breakdown point for the joint estimators and prove their corresponding high breakdown point properties. Empirical evaluations confirm the effectiveness of the proposed detection framework. Keywords: AI oversight; LLM; MDS; Robust covariance estimation; MCD; Joint estimation 1. Introduction AI, particularly with the advent of foundation models, has rapidly become a generalpurpose technology with broad scientific, industrial, and societal impact [1]. 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 [2]. 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 [3, 4, 5]. 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 ∗

Corresponding author. Email addresses: [email protected] (Xifeng Zhang), [email protected] (Tao Hu), [email protected] (Yijie Peng), [email protected] (Wan Tian)

to internal model information, such as parameters, gradients, logits, or decoding probabilities, and can exploit model-specific signals or embedded mechanisms [5]. 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 [6]. 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 [7, 8]; hallucination detection seeks to find outputs inconsistent with factual, contextual, or self-consistency-based evidence [9, 6]; adversarial example detection aims to recognize abnormal inputs or representations [3, 4, 10, 11]; and watermark detection uncovers deliberately embedded statistical signals in generated content [5]. 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 [7, 8]. Classifier-based methods train supervised or contrastive models on learned representations, with applications to LLM-generated text detection, hallucination detection, and adversarial detection [12]. Watermark-specific methods embed and subsequently detect statistical signatures in generated outputs [5]. However, these methods are often taskspecific, 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 [13] and cellwise MCD [14] 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 Lee et al. [15], who applied classconditional 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. 2

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. ——————

Single-Class ——————

1) Estimating

——————

Multi-Class

——————

1) Estimating Latent Representations

Latent Representations Pretrained Model

Casewise ...

Cellwise

MCD CellMCD

Pretrained Model

Casewise

Cellwise

...

...

...

... ... Assets

Assets

Assets Assets

，

Assets

2) Scoring

...

Joint casewise/cellwise MCD estimator

...

Assets Assets Assets

2) Scoring Latent Representations

Latent Representations Pretrained Model

Casewise ...

Cellwise

Mahalanobis Distance

Pretrained Model

Casewise

Cellwise

...

...

...

... ... Assets

Assets

Sample Under Test

Assets

MD-based Score

Assets Assets

Feature Extraction

Detector

...

Mahalanobis Distance

...

Assets Assets Assets

MD-based Score

Positive / Negative

Figure 1: Overview of the proposed framework for detecting AI-related content and artifacts under single-class and multi-class settings.

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 D1 = {(xi , Yi = 1)}ni=1 , and let the set of negative samples be S 2 (k) (k) (k) (k) k D2 = K = −1)}m i=1 , k = 1, 2, . . . , K2 . The complete dataset is k=1 D2 , where D2 = {(x̃i , Yi denoted by D = D1 ∪ D2 . Let M be a pretrained model trained on positive samples and used to extract latent representations. For each input x ∈ Rp , we obtain a representation z ∈ Rq from a chosen layer of M. The corresponding representation sets are denoted by Z1 = {(z i , Yi = 1)}ni=1 ,

Z2 =

K2 [ k=1

3

(k)

Z2 ,

(k)

(k)

(k)

k where Z2 = {(z̃ i , Yi = −1)}m i=1 , k = 1, 2, . . . , K2 . To characterize the positive class, we estimate its mean vector µ and covariance matrix Σ based on the deep representations of positive samples in Z1 . 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 b and estimation of µ and Σ can be unstable. Therefore, we use robust estimators, denoted by µ b whose construction will be detailed in Section 3. Using these estimates, the Mahalanobis Σ, distance scores for positive and negative samples are defined as

b −1 (z i − µ b ), b )⊤ Σ di = (z i − µ

(k) (k) b −1 (z̃ (k) − µ b ), b )⊤ Σ d˜i = (z̃ i − µ i

k = 1, 2, . . . , K2 .

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 S1 = {(di , Yi = 1)}ni=1 ,

(k) (k) (k) k S2 = {(d˜i , Yi = −1)}m i=1 ,

k = 1, 2, . . . , K2 ,

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 perturbed and negative sample sets be D1 = SK1 adversarially SK2 (k)counterparts. Let the positive (k) (k) (k) (k) k D and D = D , respectively, where D = {(x = 1)}ni=1 , k = 1, 2, . . . , K1 2 1 i , Yi k=1 1 k=1 2 (k) (k) (k) mk and D2 = {(x̃i , Yi = −1)}i=1 , k = 1, 2, . . . , K2 . As in the single-class setting, a pretrained model M is used to extract latent representations. The corresponding representation sets are denoted by K1 K2 [ [ (k) (k) Z1 = Z1 , Z2 = Z2 , k=1 (k)

(k)

(k)

k=1 (k)

(k)

(k)

k k and Z2 = {(z̃ i , Yi = −1)}m where Z1 = {(z i , Yi = 1)}ni=1 i=1 . To characterize the posi(k) tive classes, we estimate class-specific mean vectors µ and covariance matrices Σ(k) based on (k) the deep representations of positive samples in Z1 , k = 1, 2, . . . , K1 . 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 µ(k) and Σ(k) can be unstable. Therefore, we use joint robust estimators, 1 b (k) K1 denoted by {b µ(k) }K k=1 and {Σ }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: n o (k) (k) (k) (j) ⊤ b (j) −1 (j) b ) (Σ ) (z i − µ b ) , k = 1, 2, . . . , K1 , di = min (z i − µ 1≤j≤K1 n o (k) (k) (j) b (j) )−1 (z̃ (k) − µ b (j) )⊤ (Σ b d˜i = min (z̃ i − µ ) , k = 1, 2, . . . , K2 . i

1≤j≤K1

To ensure that the Mahalanobis distance is well-defined in high-dimensional settings, a positivedefiniteness 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 4

the dimensionality of the representations before covariance estimation. Details are provided in Appendix A.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(k) (k) (k) (k) (k) (k) k dimensional score sets S1 = {(di , Yi = 1)}ni=1 , k = 1, 2, . . . , K1 , and S2 = {(d˜i , Yi = mk −1)}i=1 , k = 1, 2, . . . , K2 , 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 Z = (z 1 , z 2 , . . . , z n )⊤ ∈ Rn×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  1 representations, 1 z i ∼ Nq (µ, Σ), i = 1, 2, . . . , n. The density is f (z i ; µ, Σ) = (2π)q/2 |Σ|1/2 exp − 2 MD2 (z i , µ, Σ) , where MD2 (z i , µ, Σ) = (z i − µ)⊤ Σ−1 (z i − µ). The negative log-likelihood, up to a constant factor, is ℓ(µ, Σ; z i ) = ln |Σ| + MD2 (z i , µ, Σ) + q ln(2π). P If all observations are reliable, minimizing ni=1 ℓ(µ, Σ; 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 H = (w1 , w2 , . . . , wnP )⊤ , where wi ∈ {0, 1}, i ∈ {1, 2, . . . , n}, indicates whether the i-th observation is retained, with ni=1 wi = h. The casewise MCD estimator solves min

µ,Σ,H

n X

wi ℓ(µ, Σ; z i ),

s.t. wi ∈ {0, 1},

i=1

n X

wi = h.

(1)

i=1

For fixed weights H, the estimators are n

b= µ

1X wi z i , h i=1

n 1X b b )(z i − µ b )⊤ . Σ= wi (z i − µ h i=1

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 Rousseeuw and Van Driessen [16], 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 = ⌊(n + q + 1)/2⌋, which attains the highest breakdown point, while larger values such as b b and Σ h = ⌊0.75n⌋ are often used when higher efficiency is desired [13, 16]. The resulting µ are then used in the single-class MDS detector. 5

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 W = (wij )1≤i≤n,1≤j≤q ∈ {0, 1}n×q be a cellwise indicator matrix, where wij = 1 means that the j-th entry of z i is retained, and wij = 0 means that the entry is flagged as con(w ) taminated. Denote by wi the i-th row of W . For each observation, let z i i be the subvector of z i containing only retained entries. Similarly, µ(wi ) is the corresponding subvector of µ, and P Σ(wi ) is the corresponding principal submatrix of Σ. Let q (wi ) = qj=1 wij 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 z i is      1 1 (wi ) (wi ) (wi ) (wi ) 2 (wi ) (wi ) , exp − MD z i , µ , Σ f zi ; µ , Σ = (w ) 2 (2π)q i /2 |Σ(wi ) |1/2 where the partial Mahalanobis distance is 2

MD



(w ) z i i , µ(wi ) , Σ(wi )



 −1  ⊤   (wi ) (wi ) (wi ) (wi ) (wi ) . zi − µ Σ = zi − µ

The corresponding negative log-likelihood is     (w ) (w ) ℓ µ(wi ) , Σ(wi ) ; z i i = ln |Σ(wi ) | + q (wi ) ln(2π) + MD2 z i i , µ(wi ) , Σ(wi ) . Unlike ordinary missing-data problems, the cellwise indicator matrix 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 ∥W ·j ∥0 ≥ h, j = 1, 2, . . . , q, where ∥W ·j ∥0 is the number of retained entries in the j-th column. We also impose λmin (Σ) ≥ a, where λmin (Σ) denotes the smallest eigenvalue of Σ, for some a > 0, thereby ensuring that the covariance matrix is nonsingular. The cellwise MCD estimator is then defined as the solution to min

µ,Σ,W

q n  X  X (wi ) (wi ) (wi ) + ℓ µ , Σ ; zi bj ∥1n − W ·j ∥0 i=1

s.t. W ∈ {0, 1}n×q ,

(2)

j=1

∥W ·j ∥0 ≥ h,

j = 1, 2, . . . , q,

λmin (Σ) ≥ a.

Here 1n denotes the n-dimensional all-ones vector, ∥1n − W ·j ∥0 counts the number of flagged cells in the j-th feature dimension, and bj > 0 controls the penalty for flagging cells in that dimension. If an observation has no retained entries, i.e., q (wi ) = 0, we set ln |Σ(wi ) | = 0 and (w ) MD2 (z i i , µ(wi ) , Σ(wi ) ) = 0. It is worth noting that directly solving (2) is computationally challenging, because the objective involves both continuous parameters (µ, Σ) and a high-dimensional binary matrix W . Following Raymaekers and Rousseeuw [14], 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 bj 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 bj jointly control the robustness–efficiency 6

trade-off: more aggressive flagging improves robustness to cellwise contamination, whereas reb are b and Σ taining more cells improves efficiency when contamination is mild. The resulting µ 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 (k) (k) within each class. Let Z = {Z (1) , . . . , Z (K) }, where Z (k) = (z 1 , . . . , z nk )⊤ ∈ Rnk ×q denotes the representation matrix of the k-th positive class, k = 1, 2, . . . , K. As a working model, we assume (k) z i ∼ Nq (µ(k) , Σ(k) ), i = 1, 2, . . . , nk , k = 1, 2, . . . , K, where µ(k) ∈ Rq and Σ(k) ∈ Rq×q denote the class-specific mean vector and covariance matrix. Following the reparameterization idea in the joint estimation of multiple graphical models [17], we introduce the following decomposition to explicitly model both shared and class-specific dependence structures: (3)

Σ(k) = Θ ⊙ Γ(k) , (k)

where ⊙ denotes the Hadamard product. Equivalently, for Σ(k) = (σij )1≤i,j≤q , ( (k)

σij =

(k)

θij γij , i ̸= j, (k) i = j, γii ,

with θii = 1 for i = 1, 2, . . . , q. Here Θ = (θij )1≤i,j≤q represents the shared covariance pattern (k) across classes, whereas Γ(k) = (γij )1≤i,j≤q captures class-specific variation. This reparameterization separates shared dependence (Θ) from class-specific deviations (Γ(k) ). We use the Root Mean Square (RMS) correlation normalization rule to ensure the identifiability of the decomposition, and details are given in Appendix Appendix A.2. (k) For k = 1, 2, . . . , K, let wi ∈ {0, 1} indicate whether the i-th observation in class k is (k) (k) retained. We write H (k) = (w1 , . . . , wnk )⊤ for the vector of casewise weights in class k with Pnk (k) = hk . The joint casewise MCD estimator is defined as the solution to i=1 wi nk K X X

min {µ(k) ,H (k) ,Σ(k) }K k=1

(k)

(k)

wi ℓ̃(µ(k) , Σ(k) ; z i ) + P(Θ, {Γ(k) }K k=1 ) (4)

k=1 i=1 (k)

s.t. wi

∈ {0, 1},

nk X

(k)

wi

= hk ,

Σ(k) = Θ ⊙ Γ(k) ,

λmin (Σ(k) ) ≥ ak ,

i=1

where

(k)

(k)

ℓ̃(µ , Σ

(k) ; z i ) = ln |Σ(k) | + q ln(2π) + MD2



(k) z i , µ(k) , Σ(k)



  (k) (k) (k) with MD2 z i , µ(k) , Σ(k) = (z i − µ(k) )⊤ (Σ(k) )−1 (z i − µ(k) ), k = 1, 2, . . . , K. Following

7

Guo et al. [17], we impose the penalty P(Θ, {Γ(k) }K k=1 ) = λ1 ∥Θ∥1,off + λ2

K X

(k)

∥Γ ∥1,off := λ1

k=1

X

|θij | + λ2

i̸=j

K X X

(k)

(5)

|γij |,

k=1 i̸=j

where λ1 , λ2 > 0 control the global and class-specific sparsity levels, respectively. Here, ∥ · ∥1,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 Θ: when θ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 Γ(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 λmin (Σ(k) ) ≥ ak > 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 (k) (k) W (k) = (wij )1≤i≤nk ,1≤j≤q ∈ {0, 1}nk ×q be the cellwise indicator matrix, where wij = 1 means (k) (k) that the j-th entry of z i is retained, and wij = 0 means that the entry is flagged as conP (k) (k) (k) taminated. Let wi denote the i-th row of W (k) , and let q (wi ) = qj=1 wij be the number (k)

of retained entries in observation i of class k. For any vector or matrix, the superscript (wi ) denotes the subvector or principal submatrix corresponding to the retained entries. Thus, (w

(k)

,k)

(k)

(k)

(k)

zi i , µ(wi ,k) , Σ(wi ,k) denote the retained subvector of z i , the corresponding subvector of (k) µ , and the corresponding principal submatrix of Σ(k) , respectively. Using the same covariance decomposition in (3), the joint cellwise MCD estimator is defined by minimizing   nk K X X (k) (k) (k) (wi ,k) (wi ,k) (wi ,k) ℓ̃ µ ,Σ ; zi

min {µ(k) ,W (k) ,Σ(k) }K k=1

k=1 i=1 q K X X

+

(k)

bj

(k)

1nk − W ·j

0

k=1 j=1 (k)

s.t. W (k) ∈ {0, 1}nk ×q , Σ(k) = Θ ⊙ Γ(k) , (w

(k)

Let MD2 (z i i

,k)

(k)

(k)

(w

W ·j

0

+ P(Θ, {Γ(k) }K k=1 )

≥ hk ,

j = 1, 2, . . . , q,

λmin (Σ(k) ) ≥ ak > 0, (k)

, µ(wi ,k) , Σ(wi ,k) ) = (z i i

,k)

k = 1, 2, . . . , K.

(k)

(k)

(w

(k)

(w

(k)

,k)

(k)

(k)

(k)

− µ(wi ,k) )⊤ (Σ(wi ,k) )−1 (z i i (k)

(w

(k)

be the partial Mahalanobis distance, we have ℓ̃(µ(wi ,k) , Σ(wi ,k) ; z i i (k)

(6)

,k)

,k)

(k)

− µ(wi ,k) ) (k)

) = ln |Σ(wi ,k) | +

q (wi ) ln(2π)+MD2 (z i i , µ(wi ,k) , Σ(wi ,k) ). Here, 1nk is the nk -dimensional all-ones vector, (k) and ∥1nk − W ·j ∥0 counts the number of flagged cells in the j-th feature dimension of class k. (k) (k) The constants bj penalize cellwise flagging, while the constraints ∥W ·j ∥0 ≥ hk prevent too many entries from being removed in any feature dimension. If all entries of an observation are (k) (k) flagged, we adopt the empty-pattern convention q (wi ) = 0, ln |Σ(wi ,k) | = 0, and MD2 = 0, so that the objective function remains well-defined. λmin (Σ(k) ) ≥ ak > 0 guarantees that the covariance matrix is nonsingular. The penalty P(Θ, {Γ(k) }K k=1 ) is the same as in (5), 8

encouraging homogeneity through the shared component Θ while preserving heterogeneity by the class-specific components Γ(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 Σ(k) = Θ ⊙ Γ(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 ∈ N as a superscript on the estimator to represent the result obtained in the t-th iteration, where 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 H (k,0) of size hk . 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: Pnk (k,t) (k) Pnk (k,t) (k) w zi zi (k,t+1) i=1 wi b = i=1 i , k = 1, 2, . . . , K. µ = Pn (k,t) k hk i=1 wi b (t+1) and the group-specific heterogeneity compoWe next update the shared structure Θ b (k,t+1) via alternating optimization: nents Γ (t)

b • We first fix Θ

b and update each Γ

b (k,t+1) = arg min Γ Γ

nk X

(k,t+1)

(k,t)

wi

:

b ℓ̃(b µ(k,t+1) , Θ

(t)

(k)

⊙ Γ; z i ) + λ2

i=1

X

|γij |.

i̸=j

b (k,t+1) and update the shared structure Θ b (t+1) : • We then fix all Γ b (t+1) = arg min Θ Θ

nk K X X

(k,t)

wi

b (k,t+1) ; z (k) ) + λ1 ℓ̃(b µ(k,t+1) , Θ ⊙ Γ i

k=1 i=1

X

|θij |.

i̸=j

b (k,t+1) = Θ b (t+1) ⊙ Γ b (k,t+1) , k = 1, 2, . . . , K. After The covariance matrices are updated by Σ 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 hk samples with the smallest Mahalanobis distances are selected as unflagged samples, forming the updated subset 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 [16]. 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 1.

9

Algorithm 1 Fast-MCD Algorithm for Joint Casewise MCD Estimator (JCASEMCD) Require: Data Z, subset sizes {hk }K k=1 , penalties λ1 , λ2 , error threshold ε, eigenvalue threshK olds {ak }k=1 . b (0) , Γ b (k,0) , Σ b (k,0) , and initial subsets b (k,0) , Θ Initialize t = 0 and generate initial estimates µ H (k,0) for k = 1, 2, . . . , K. while not converged do Update group means: (k,t+1)

b µ

(k,t) (k) zi i=1 wi

Pnk ←

hk

,

k = 1, 2, . . . , K.

Update structure parameters: b (t) and update Γ b (k,t+1) : 1. For k = 1, 2, . . . , K, fix Θ b (k,t+1) ← arg min Γ Γ

b 2. Fix all Γ

(k,t+1)

nk X

Θ



(k,t+1)

b µ

 X (t) (k) b |γij |. + λ2 , Θ ⊙ Γ; z i

i=1

b and update Θ

b (t+1) ← arg min Θ

(k,t) wi ℓ̃

i̸=j

(t+1)

nk K X X

:

(k,t) wi ℓ̃



(k,t+1)

b µ

 X (k,t+1) (k) b + λ1 |θij |. ,Θ ⊙ Γ ; zi

k=1 i=1

i̸=j

Update covariance matrices: b (k,t+1) ← Θ b (t+1) ⊙ Γ b (k,t+1) , Σ

k = 1, 2, . . . , K.

Update subsets H (k,t+1) with the hk samples having smallest Mahalanobis distances. b (k,t+1) − Σ b (k,t) ∥F < ε and maxk ∥b b (k,t) ∥2 < ε then if maxk ∥Σ µ(k,t+1) − µ Converged, exit loop. else t ← t + 1. end if end while (k,final) Reweighting: Update weights wi based on Mahalanobis distances and recompute (k,final) (k,final) b (k,final) b (final) b µ ,Σ ,Θ , Γ̂ . (k,final) (k,final) b (k,final) b (final) b return µ ,Σ ,Θ , Γ̂

10

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 W and then updating the parameters given the missing entries. For ease of discussion, we introduce notation and partitions. In the t-th (k) (k) (k)⊤ (k)⊤ iteration, for the k-th group, the i-th row is divided into three parts z i = (zij , z i,o , z i,−o )⊤ (k) after a suitable permutation of coordinates, where zij is the element in row i and column j, and (k) (k) (k,t) (k,t) z i,o and z i,−o correspond to entries marked as normal (wiτ = 1) and abnormal (wiτ = 0), (k,t) b (k) τ = 1, 2, . . . , q, respectively, with o = {τ ̸= j : wiτ = 1}. Accordingly, the mean vector µ b (k) are partitioned as and covariance matrix Σ 



(k) µj  (k)  b (k) = µ µ bo  , b (k) µ −o



(k)

σ b  jj(k) (k) b =Σ Σ b o,j

 b (k) Σ j,−o (k)  . b Σ 

b (k) Σ j,o (k) b Σ o,o

o,−o

b (k) Σ b (k) Σ b (k) Σ −o,j −o,o −o,−o

b (k,0) , We first initialize each group separately by computing initial location estimates µ (k,0) b covariance estimates Σ , and weight matrices W (k,0) , ensuring that each column contains at least hk non-missing entries not flagged as outliers. An iterative process then alternates b (k,t) between updating cellwise weights and updating parameters. In the t-th iteration, with µ b (k,t) fixed, the objective function used for updating the weight of entry (i, j) in group k is and Σ (k) (k) (zij − zbij )2 (k) (k) (k) − bj , ∆ij = ln(Cij ) + ln(2π) + (k) Cij

(7)

(k) b (k,0) )−1 ]jj ), 0}, with χ2 (1) representing the chiwhere bj = max{χ21−α (1) + ln(2π) + ln(1/([(Σ 1−α square 1 − α quantile with 1 degree of freedom. Moreover, the conditional mean and variance are (k)

(k,t)

bj zbij = µ

(k,t)

(k,t)

(k)

−1 b b b (k,t) ), +Σ o j,o (Σo,o ) (z i,o − µ

(k)

(k,t)

bjj Cij = σ

(k,t)

(k,t)

(k,t)

−1 b ⊤ b b −Σ j,o (Σo,o ) (Σj,o ) .

(8)

(k)

The index set o is determined by the current weight matrix W (k,t) . All cells with ∆ij < 0 are retained. If fewer than hk cells are retained in the j-th column, the smallest remaining (k) values of ∆ij are additionally retained until hk cells are retained [14]. Then, the Expectationb (k,t+1) , and Θ b (t+1) . b (k,t+1) , Γ Maximization (EM) algorithm [18] updates µ E-step. Given W (k,t+1) , the conditional expectation and covariance for the i-th observation (k,t+1) (k,t+1) (k) in group k are computed. Let Oi = {τ : wiτ = 1}, Ui = {τ : wiτ = 0}, and let z i,Oi be (k,t)

a subvector indexed by Oi . The conditional expectation is mi (k,t)

(k)

mi,Oi = z i,Oi ,

(k,t)

(k,t)

(k,t)

(k,t) ⊤

(k,t) b (k,t) b (k,t) −1 (k) b (k,t) b (k,t) mi,Ui = µ Ui + ΣUi ,Oi (ΣOi ,Oi ) (z i,Oi − µ Oi ),

and the conditional covariance Vi Vi,Oi ,Oi = O,

(k,t) ⊤

= (mi,Oi , mi,Ui )⊤ , where

(k,t)

(9)

satisfies

Vi,Oi ,Ui = Vi,Ui ,Oi = 0,

(k,t) b (k,t) − Σ b (k,t) (Σ b (k,t) )−1 Σ b (k,t) , Vi,Ui ,Ui = Σ Ui ,Ui Ui ,Oi Oi ,Oi Oi ,Ui

(10)

where O and 0 represent a matrix and a vector containing only zeros, respectively. When (k,t) b (k,t) . The conditional second-order b (k,t) and Vi(k,t) = Σ Oi = ∅, we use the convention mi = µ

11

central moment is then h i (k) (k) (k) (k,t) (k,t) (k,t) E (z i − µ(k) )(z i − µ(k) )⊤ z i,Oi = Vi + (mi − µ(k) )(mi − µ(k) )⊤ . P k (k,t) b (k,t+1) : µ b (k,t+1) = n1k ni=1 mi , and M-step. We first update the location parameters µ compute nk h i 1 X (k,t) (k,t) (k,t) (k,t+1) ⊤ b (k,t+1) )(m(k,t) b Scell = Vi + (mi − µ − µ ) . i nk i=1 b (t) , the group-specific structure Γ b (k,t+1) is updated by solving Next, fixing Θ h i X (k,t+1) (t) (t) −1 (k,t) b b b Γ = arg min nk ln |Θ ⊙ Γ| + tr((Θ ⊙ Γ) Scell ) + λ2 |γij |. Γ

i̸=j

b (t+1) is updated by Finally, the shared structure Θ b (t+1) = arg min Θ Θ

K X

h i X b (k,t+1) | + tr((Θ ⊙ Γ b (k,t+1) )−1 S (k,t) ) + λ1 nk ln |Θ ⊙ Γ |θij |. cell

k=1

i̸=j

b (k,t+1) = Θ b (t+1) ⊙ The covariance matrices for each group are then reconstructed via Σ b (k,t+1) , and the algorithm iterates until convergence. The full procedure is summarized in Γ Algorithm 2. Algorithm 2 Computation procedure for joint cellwise MCD estimator (JCELLMCD) K Require: Z = (Z (1) , . . . , Z (K) ), {hk }K k=1 , {ak }k=1 , λ1 , λ2 . b (k,0) , W (k,0) , Θ b (0) , Γ b (k,0) and evaluate b(k) for k = 1, 2, . . . , K. b (k,0) , Σ Initialize t = 0, µ j while not converged do For k = 1, 2, . . . , K, update cellwise weights (Algorithm 3): (k,t)

b b (k,t) , Σ W (k,t+1) = CellwiseWeightsUpdater(Z (k) , µ

(k)

, W (k,t) , hk , bj ).

Update parameters via EM (Algorithm 4): b (b µ(k,t+1) , Θ

(t+1)

b ,Γ

(k,t+1) K b (t) , Γ b (k,t) )K , {ak }K , λ1 , λ2 ). )k=1 = EM(Z, {W (k,t+1) }K µ(k,t) , Θ k=1 , (b k=1 k=1

b (k,t+1) = Θ b (t+1) ⊙ Γ b (k,t+1) , k = 1, 2, . . . , K. Reconstruct covariance matrices: Σ Check convergence and update t = t + 1. end while b (k,final) , Θ b (final) , Γ b (k,final) , W (k,final) b (k,final) , Σ return µ

12

Algorithm 3 CellwiseWeightsUpdater step for Algorithm 2 (k,t)

(k)

b b (k,t) ,Σ Require: Z (k) , µ , W (k,t) , hk , bj W (k,t+1) = W (k,t) for j = 1, . . . , q do for i = 1, . . . , nk do (k,t+1) o = {τ ̸= j : wiτ = 1}. Compute conditional mean and variance with Equation (8). (k) Compute objective function ∆ij with Equation (7). end for (k) Let N− = #{∆ij < 0, i = 1, 2, . . . , nk }, and update weights: (k,t+1)

wij

(k)

(k)

= I{N− ≥ hk }I{∆ij < 0} + I{N− < hk }I{∆ij among hk smallest}

end for return W (k,t+1)

Algorithm 4 EM step for Algorithm 2 b (t) , Γ b (k,t) )K , {ak }K , λ1 , λ2 . µ(k,t) , Θ Require: Z, {W (k,t+1) }K k=1 k=1 k=1 , (b for k = 1, . . . , K do for i = 1, . . . , nk do (k) (k) Partition indices: Oi = {τ : wiτ = 1}, Ui = {τ : wiτ = 0}. (k,t) with Equation (9). Compute conditional mean mi (k,t) with Equation (10). Compute conditional covariance Vi end for P k (k,t) b (k,t+1) = n1k ni=1 Update location: µ mi   P (k,t) (k,t) (k,t) k b (k,t+1) )⊤ . b (k,t+1) )(m(k,t) −µ Compute Scell = n1k ni=1 + (mi − µ Vi i Update group-specific structures: h i X (k,t+1) (t) (t) −1 (k,t) b b b Γ = arg min nk ln |Θ ⊙ Γ| + tr((Θ ⊙ Γ) Scell ) + λ2 |γij |. Γ

i̸=j

end for Update shared structure: b (t+1) = arg min Θ Θ

b return (b µ(k,t+1) , Θ

K X

h i X b (k,t+1) | + tr((Θ ⊙ Γ b (k,t+1) )−1 S (k,t) ) + λ1 nk ln |Θ ⊙ Γ |θij |. cell

k=1

(t+1)

i̸=j (k,t+1)

b ,Γ

)K k=1

b and {Γ b (k) }K , which correspond to the Fast-MCD algorithm for the The updates for Θ k=1 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 Appendix A.4. Details such as algorithm initialization, the identifiability of the covariance decomposition, and dimensionality reduction in the implementation are also given in Appendix A. 13

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 Appendix B. 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 6.1. Let Ω = (Ω(1) , . . . , Ω(K) ), with Ω(k) = (µ(k) , Σ(k) ) and Σ(k) = Θ ⊙ Γ(k) for k = 1, 2, . . . , K. Define the regularized joint negative log-likelihood objective function as Lcase (Ω, H) =

nk K X X

(k)

(k)

wi ℓ̃(µ(k) , Σ(k) ; z i ) + P(Θ, {Γ(k) }K k=1 )

k=1 i=1 K K X X (k) := Qcase (Ω, H) + λ1 ∥Θ∥1,off + λ2 ∥Γ(k) ∥1,off , k=1 k=1

P k (k) (k) wi = hk . Suppose λmin (Σ(k) ) ≥ ak > 0, λ1 > 0, and λ2 > 0. For where wi ∈ {0, 1} and ni=1 (t) (t) each iteration t, let LΓ,k and LΘ be the Lipschitz constants of the gradients of the smooth parts of the Γ(k) - and Θ-subproblems, respectively. The proximal gradient step sizes are chosen with (s) a strict descent margin, namely, for some constants δΓ , δΘ ∈ (0, 1), LΓ,k ≥ sups≥0 LΓ,k > 0 and (s)

LΘ ≥ sups≥0 LΘ > 0, k = 1, 2, . . . , K, (s)

(t)

0 < inf ηΓ,k ≤ ηΓ,k ≤ s≥0

1 − δΓ , LΓ,k

(s)

(t)

0 < inf ηΘ ≤ ηΘ ≤ s≥0

1 − δΘ . LΘ

Then the pre-reweighting iterations of Algorithm 1 satisfy: 1. Lcase (Ω(t+1) , H (t+1) ) ≤ Lcase (Ω(t+1) , H (t) ) ≤ Lcase (Ω(t) , H (t) ); 2. The sequence Lcase (Ω(t) , H (t) ) converges to a finite value L∗case ; 3. Any interior accumulation point (Ω∗ , H ∗ ), namely one satisfying λmin (Θ∗ ⊙Γ(k)∗ ) > ak > 0, k = 1, 2, . . . , K is blockwise stationary with respect to the continuous variables, and H ∗ is optimal for the discrete H-block given Ω∗ . Theorem 6.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.

14

Theorem 6.2. Let Ω = (Ω(1) , . . . , Ω(K) ), with Ω(k) = (µ(k) , Σ(k) ) and Σ(k) = Θ ⊙ Γ(k) for k = 1, 2, . . . , K. Denote the cellwise objective function as   nk K X X (k) (k) (k) (wi ,k) (wi ,k) (wi ,k) Lcell (Ω, W ) = ℓ̃ µ ,Σ ; zi k=1 i=1 q K X X

+

(k) (k) bj ∥1nk − W ·j ∥0 + λ1 ∥Θ∥1,off + λ2

k=1 j=1

K X

∥Γ(k) ∥1,off ,

k=1 (k)

(k)

subject to λmin (Σ(k) ) ≥ ak > 0, ∥W ·j ∥0 ≥ hk , bj ≥ 0 for all j and k, λ1 > 0, and λ2 > 0. The proximal gradient step sizes are chosen with a strict descent margin, namely, for some constants δΓ , δΘ ∈ (0, 1), ∀ t ≥ 0 and k = 1, 2, . . . , K, (s)

(t)

0 < inf ηΓ,k ≤ ηΓ,k ≤ s≥0

1 − δΓ , LΓ,k

(s)

(t)

0 < inf ηΘ ≤ ηΘ ≤ s≥0

1 − δΘ , LΘ

where LΓ,k > 0, LΘ > 0 are Lipschitz constants. Algorithm 2 satisfies: 1. Lcell (Ω(t+1) , W (t+1) ) ≤ Lcell (Ω(t) , W (t+1) ) ≤ Lcell (Ω(t) , W (t) ); 2. The sequence Lcell (Ω(t) , W (t) ) converges to a finite value; 3. Any interior limit point (Ω∗ , W ∗ ), namely one satisfying λmin (Θ∗ ⊙ Γ(k)∗ ) > ak > 0, k = 1, 2, . . . , K is stationary with respect to the continuous variables for fixed W ∗ , and W ∗ is coordinatewise optimal for the discrete cellwise weight update given Ω∗ . Theorem 6.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 Z ⋆ = ⋆(k) ⋆(k) {Z ⋆(1) , . . . , Z ⋆(K) }, where Z ⋆(k) = (z 1 , . . . , z nk )⊤ ∈ Rnk ×q denotes the clean representation b = (b b = ({b b (k) K matrix of the k-th positive class, k = 1, 2, . . . , K. Let Ω µ, Σ) µ(k) }K k=1 , {Σ }k=1 ) denote the joint estimator. For any vector r = (r1 , r2 , . . . , rK ), define rk . 1≤k≤K nk

β(r) := max

Definition 6.1 (Joint Casewise Contamination). Let Ccase (r) be the set of all contaminated (1) (K) (k) samples Z̃ = {Z̃ , . . . , Z̃ }, where for each group k ∈ {1, 2, . . . , K}, Z̃ is obtained by replacing at most rk rows of Z ⋆(k) .

15

Definition 6.2 (Joint Cellwise Contamination). Let Ccell (r) be the set of all contaminated (k) samples Z̃, where for each group k ∈ {1, 2, . . . , K} and each column j ∈ {1, 2, . . . , q}, Z̃ is obtained by replacing at most rk cells in the j-th column of Z ⋆(k) . Definition 6.3 (Finite-Sample Joint Breakdown Points). The joint location breakdown points are ( ) b Z ⋆ ) := inf β(r) : sup max ∥b b (k) (Z ⋆ )∥2 = ∞ , (Ω, δ∗ µ(k) (Z̃) − µ case,µ

Z̃∈Ccase (r)

1≤k≤K

( b Z ⋆ ) := inf ε∗cell,µ (Ω,

)

β(r) :

sup Z̃∈Ccell (r)

b (k) (Z ⋆ )∥2 = ∞ . max ∥b µ(k) (Z̃) − µ

1≤k≤K

The joint scatter explosion and implosion breakdown points are ( ) (k) b Z ⋆ ) := inf β(r) : sup max λmax (Σ b (Z̃)) = ∞ , δcase,+ (Ω, Z̃∈Ccase (r)

1≤k≤K

( b Z ⋆ ) := inf εcell,+ (Ω,

β(r) :

sup Z̃∈Ccell (r)

(k)

)

b (Z̃)) = ∞ , max λmax (Σ

1≤k≤K



 (k) ⋆ b b δcase,− (Ω, Z ) := inf β(r) : inf min λmin (Σ (Z̃)) = 0 , Z̃∈Ccase (r) 1≤k≤K   (k) ⋆ b Z ) := inf β(r) : inf b (Z̃)) = 0 . εcell,− (Ω, min λmin (Σ Z̃∈Ccell (r) 1≤k≤K

If no admissible contamination fraction in [0, 1] causes breakdown, the breakdown fraction is reported as 1. 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 6.3 (Joint Cellwise Breakdown Point). Consider the joint cellwise estimator mini(k) mizing the joint objective subject to λmin (Σ(k) ) ≥ ak > 0, ∥W ·j ∥0 ≥ hk , and a non-negative penalty K X X X (k) (k) K P(Θ, {Γ }k=1 ) = λ1 |θij | + λ2 |γij |. i̸=j

k=1 i̸=j

Assume each clean group Z ⋆(k) is in general position and jn k k + 1, k = 1, 2, . . . , K. hk ≥ 2 Let

nk − h k + 1 . 1≤k≤K nk

β0 := min

16

Then, b Z ⋆ ) = 1, εcell,− (Ω,

b Z ⋆ ) ≥ β0 , εcell,+ (Ω,

b Z ⋆ ) ≥ β0 . ε∗cell,µ (Ω,

If nk = n and hk = h for all groups, these simplify to εcell,− = 1,

εcell,+ ≥

n−h+1 , n

ε∗cell,µ ≥

n−h+1 . n

Theorem 6.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 β0 is sharp as a common guarantee for simultaneous boundedness of the location and scatter estimates under the stated contamination construction. Theorem 6.4 (Joint Casewise Breakdown Point). Suppose each clean group is in general position and hk ≥ ⌊nk /2⌋ + 1, k = 1, 2, . . . , K. For the penalized joint casewise estimator with λmin (Σ(k) ) ≥ ak > 0 and P(Θ, {Γ(k) }K k=1 ) ≥ 0, let nk − h k + 1 . 1≤k≤K nk

β0 := min Then, b Z ⋆ ) = 1, δcase,− (Ω,

b Z ⋆ ) ≥ β0 , δcase,+ (Ω,

∗ b Z ⋆ ) ≥ β0 . δcase,µ (Ω,

n−h+1 , n

n−h+1 . n

If nk = n and hk = h for all groups, δcase,− = 1,

δcase,+ ≥

∗ δcase,µ ≥

This establishes a provable high breakdown lower bound for the penalized joint casewise estimator. Theorem 6.4 similarly establishes that the joint casewise MCD estimator maintains bounded estimates even under substantial casewise contamination. The derived lower bounds on δcase,− , ∗ δcase,+ , and δcase,µ 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: LLMgenerated 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 LLMgenerated ones. We use the Human ChatGPT Comparison Corpus (HC3), a benchmark dataset widely adopted for AI-generated text detection, which contains paired human and ChatGPT3.5 responses across multiple domains [19, 20, 21]. 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 [22, 23], ChatGPT-5.4 mini [24], Gemini-3 Flash Preview [25, 26], Gemini-3.1 Flash Lite Preview [27], and Gemini-3.1 Pro Preview [28]. This yields five HC3based datasets, each with 4,842 instances. We apply the proposed unified detection framework 17

using RoBERTa-base as the pretrained encoder [29], and compare it with likelihood-based [12], rank-based [7], entropy-based [30], DetectGPT [8], and RoBERTa baseline methods. For DetectGPT, we report both the raw perturbation discrepancy score, DetectGPT-d, and its normalized version, DetectGPT-z. The results in Table 1 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 hypothesistesting procedures specifically designed for the corresponding watermarking mechanism [31]. 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 [32]. 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 inversetransform 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 Li et al. [31], which assumes access to the underlying watermarking mechanism. The results are reported in Table 2. The results in Table 2 show that the hypothesis-testing baselines of Li et al. [31] 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 Li et al. [31] depends on knowledge of the watermarking rule, its effectiveness may degrade substantially outside the matched watermark setting, as shown in Table 3. 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. 18

Table 1: Performance comparison for LLM-generated text detection. Best results are shown in bold, and secondbest results are marked with †. All metrics are computed with respect to the LLM-generated samples.

LLM

Detector

Entropy Likelihood LogRank Rank RoBERTa-base GPT-5.4 RoBERTa-large DetectGPT-d DetectGPT-z JCASEMCD JCELLMCD Entropy Likelihood LogRank Rank GPT-5.4 RoBERTa-base mini RoBERTa-large DetectGPT-d DetectGPT-z JCASEMCD JCELLMCD Entropy Likelihood LogRank Rank Gemini-3.1 RoBERTa-base flash-lite RoBERTa-large preview DetectGPT-d DetectGPT-z JCASEMCD JCELLMCD Entropy Likelihood LogRank Rank Gemini-3 RoBERTa-base flash RoBERTa-large preview DetectGPT-d DetectGPT-z JCASEMCD JCELLMCD Entropy Likelihood LogRank Rank Gemini-3.1 RoBERTa-base pro RoBERTa-large preview DetectGPT-d DetectGPT-z JCASEMCD JCELLMCD

ROC AUC 0.022 0.956 0.954 0.881 0.231 0.235 0.974 0.958 0.969 0.963 0.027 0.958 0.954 0.854 0.251 0.264 0.967 0.952 0.951 0.945 0.177 0.838 0.848 0.860 0.358 0.295 0.714 0.719 0.798 0.737 0.180 0.823 0.834 0.853 0.344 0.274 0.690 0.691 0.823 0.763 0.287 0.698 0.740 0.710 0.275 0.179 0.518 0.506 0.813 0.748

PR AUC 0.309 0.937 0.934 0.843 0.351 0.353 0.976 0.949 0.970 0.957 0.309 0.937 0.934 0.825 0.357 0.363 0.970 0.940 0.955 0.939 0.338 0.845 0.855 0.832 0.410 0.372 0.807 0.804 0.770 0.688 0.338 0.832 0.843 0.828 0.401 0.364 0.785 0.776 0.796 0.715 0.377 0.722 0.707 0.732 0.365 0.338 0.666 0.653 0.781 0.695

F1

Specificity Accuracy Precision Recall

0.667 0.928 0.932 0.840 0.667 0.667 0.928 0.921 0.915 0.913 0.667 0.929 0.929 0.799 0.667 0.667 0.919 0.916 0.892 0.887 0.667 0.790 0.798 0.804 0.669 0.667 0.713 0.720 0.737 0.700 0.667 0.776 0.787 0.787 0.668 0.667 0.691 0.686 0.768 0.717 0.667 0.669 0.699 0.682 0.667 0.667 0.667 0.667 0.754 0.717

19

0.000 0.897 0.894 0.792 0.000 0.000 0.907 0.893 0.973 0.905 0.000 0.886 0.893 0.746 0.000 0.000 0.904 0.885 0.925 0.873 0.000 0.838 0.867 0.751 0.012 0.000 0.878 0.860 0.571 0.584 0.000 0.814 0.814 0.751 0.007 0.000 0.803 0.827 0.628 0.616 0.000 0.015 0.476 0.734 0.007 0.000 0.001 0.000 0.606 0.600

0.500 0.926 0.930 0.833 0.500 0.500 0.926 0.919 0.920 0.912 0.500 0.926 0.926 0.790 0.500 0.500 0.918 0.914 0.895 0.885 0.500 0.798 0.809 0.795 0.505 0.500 0.750 0.751 0.702 0.673 0.500 0.783 0.792 0.780 0.504 0.500 0.717 0.720 0.742 0.695 0.500 0.507 0.647 0.694 0.502 0.500 0.501 0.500 0.725 0.691

0.500 0.902 0.901 0.808 0.500 0.500 0.911 0.898 0.969 0.906 0.500 0.895 0.900 0.767 0.500 0.500 0.906 0.891 0.920 0.876 0.500 0.824 0.850 0.771 0.503 0.500 0.836 0.821 0.660 0.647 0.500 0.802 0.805 0.765 0.502 0.500 0.762 0.779 0.697 0.668 0.500 0.503 0.610 0.711 0.501 0.500 0.500 0.500 0.682 0.662

1.000 0.955 0.966 0.874 1.000 1.000 0.945 0.945 0.867 0.919 1.000 0.967 0.960 0.835 1.000 1.000 0.932 0.943 0.866 0.898 1.000 0.758 0.752 0.840 0.999 1.000 0.622 0.641 0.834 0.762 1.000 0.753 0.770 0.810 1.000 1.000 0.632 0.613 0.855 0.773 1.000 0.998 0.818 0.654 0.998 1.000 1.000 1.000 0.844 0.783

Mean 0.428 0.929 0.930 0.839 0.464 0.465 0.938† 0.926 0.941 0.925 0.429 0.928† 0.928† 0.802 0.468 0.471 0.931 0.920 0.915 0.900 0.454 0.813† 0.825 0.808 0.494 0.476 0.760 0.759 0.725 0.685 0.455 0.798† 0.806 0.796 0.489 0.472 0.726 0.727 0.758 0.707 0.476 0.587 0.671 0.702† 0.474 0.455 0.550 0.547 0.744 0.699

Table 2: Performance comparison across different watermark detection tasks. Best results are shown in bold, and second-best results are marked with †. All metrics are computed with respect to the watermarked samples.

Task Human vs Gumbel Human vs Inverse Null vs Gumbel Null vs Inverse

Method Gumbel-Test JCASEMCD JCELLMCD Inverse-Test JCASEMCD JCELLMCD Gumbel-Test JCASEMCD JCELLMCD Inverse-Test JCASEMCD JCELLMCD

F1 0.630 0.575 0.591† 0.625 0.575 0.589† 0.629 0.590† 0.585 0.624 0.590† 0.581

Specificity 0.947 0.607† 0.576 0.948 0.606† 0.571 0.943 0.597 0.601† 0.945 0.590 0.596†

Accuracy 0.716 0.585 0.586† 0.713 0.584† 0.584† 0.714 0.592† 0.589 0.712 0.590† 0.585

Precision 0.901 0.589† 0.585 0.902 0.588† 0.582 0.895 0.593† 0.592 0.898 0.590† 0.587

Recall 0.485 0.562† 0.597 0.479 0.562† 0.597 0.485 0.586 0.578† 0.479 0.590 0.574†

Table 3: Performance comparison for detecting human-generated text and OPT-generated text without watermark. Best results are shown in bold, and second-best results are marked with †. All metrics are computed with respect to the OPT-generated samples.

Method Gumbel-Test Inverse-Test JCASEMCD JCELLMCD

F1 0.196 0.099 0.788 0.735†

Specificity 0.930† 0.948 0.728 0.662

Accuracy 0.523 0.501 0.777 0.720†

Precision 0.624 0.511 0.752 0.697†

Recall 0.116 0.055 0.827 0.778†

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 [33], a large-scale benchmark used for reference-free hallucination detection, and the biography-domain “wiki_bio_gpt3_hallucination” dataset [6], which contains 1,908 sentence-level samples with Wikipedia evidence and hallucination labels, enabling reference-based evaluation. Table 4: Performance comparison for hallucination detection. Best results are shown in bold, and second-best results are marked with †. All metrics are computed with respect to the hallucinated text.

Method DeBERTa-v3-large-1 DeBERTa-v3-large-2 JCASEMCD(w/o ref) JCELLMCD(w/o ref) JCASEMCD(w ref) JCELLMCD(w ref)

ROC AUC 0.623 0.644† 0.608 0.751 0.554 0.557

PR AUC 0.707 0.724† 0.489 0.734 0.640 0.653

F1

Specificity

Accuracy

Precision

Recall

0.734 0.734 0.568 0.642 0.663† 0.652

0.325 0.338 0.518† 0.577 0.441 0.447

0.631 0.634† 0.576 0.648 0.581 0.571

0.665† 0.668 0.501 0.565 0.661 0.657

0.818 0.814† 0.656 0.744 0.666 0.648

As baselines, we use two DeBERTa-based NLI models, DeBERTa-v3-large-mnli-fever-anliling-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 4 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. 20

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 [34], and the Mahalanobis distance-based detector proposed by Lee et al. [15]. We study adversarial example detection in image classification tasks. On CIFAR-10 and CIFAR-100 [35], adversarial examples are generated using four representative attack methods: FGSM [3], BIM [36], CWL2 [4], and DeepFool [37]. We also consider ImageNet [38] adversarial examples generated by Diff-PGD [39]. The target classifiers used to generate these adversarial examples are mainly ResNet [40] and DenseNet [41] 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 Lee et al. [15]. The results are reported in Figure 2. The proposed framework matches or exceeds baselines in simpler settings (e.g., CIFAR10 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. Model

Dataset

1: JCASEMCD 2: JCELLMCD 3: LID Specificity Accuracy

F1

Attack Type

4: Mahalanobis Precision

Recall

BIM CIFAR-10

CWL2 DeepFool FGSM

DenseNet

BIM CIFAR-100

0.8

CWL2 DeepFool

0.6

FGSM BIM CIFAR-10

CWL2

0.4

DeepFool FGSM

ResNet

0.2

BIM CIFAR-100

CWL2 DeepFool FGSM

ResNet-50

ImageNet

Diff-PGD

1

2

3

4

1

2

3

4

1

2

3

4

1

2

3

4

1

2

3

4

Figure 2: Heatmap of different metrics. All metrics are computed with respect to the adversarial examples.

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 21

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.

22

Appendix A. Details of implementation and optimization Appendix A.1. Initialization for JCASEMCD and JCELLMCD To initialize the estimation procedure, we employ a robust univariate filtering strategy to mitigate the adverse effects of severe cellwise outliers. Recall that Z (k) is the nk × q data matrix of the k-th group, k = 1, 2, . . . , K. For each variable j ∈ {1, 2, . . . , q}, we compute (k) (k) the coordinate-wise median medj = median(z j ) and the scaled median absolute deviation (k) (k) (k) (MAD): MADj = 1.4826 · median1≤i≤nk |zij − medj |, where the factor 1.4826 = 1/Φ−1 (0.75) makes the scaled MAD equal to the standard deviation under a normal distribution. We then (k) (k) (k) calculate the absolute MAD-scaled outlier score |zij − medj |/MADj . An observation cell (k) zij is flagged as a potential outlier and assigned a binary weight wij = 0 if its score exceeds a threshold; otherwise, wij = 1. Following the suggestion of Iglewicz and Hoaglin [42], we set the threshold to 3.5. Alternative threshold values, such as 2.5 and 3, can also be considered. To guarantee numerical stability, we enforce a minimum subset size requirement, denoted as hk . If the number of unflagged cells in any given column falls below hk , we retain the hk observations with the smallest robust distances and flag the rest. Subsequently, the flagged cellwise outliers are temporarily imputed using their respective column medians, yielding a robustly imputed data matrix Z ♯(k) . b (k,0) are computed b (k,0) and scatter matrix Σ For JCELLMCD, the initial location vector µ as the empirical mean and empirical covariance matrix of Z ♯(k) , respectively, followed by a b (k,0) is strictly symmetric positive definite. For JCASEMCD, regularization step to ensure Σ the above cellwise-cleaned estimates are used as preliminary robust estimates. Specifically, we compute the Mahalanobis distances of the original observations in Z (k) using these preliminary b (k,0) and estimates, retain the hk observations with the smallest distances, and then recompute µ b (k,0) from the retained original observations, again followed by a regularization step to ensure Σ b (0) and Γ b (k,0) are generated through the Hadamard symmetric positive definiteness. Finally, Θ product decomposition, and we present a scheme in Appendix Appendix A.2 to guarantee the uniqueness of the decomposition. This procedure provides a highly resistant starting point for the subsequent iterative algorithms. Appendix A.2. Identifiability of covariance decomposition Although the covariance decomposition is defined through the Hadamard product, such a decomposition is not unique. For example, multiplying the shared structure matrix by a positive constant while dividing the heterogeneous structure by the same constant leaves their Hadamard product unchanged. To remove this scaling ambiguity and ensure identifiability, we impose an RMS-correlation normalization rule on the decomposition. (k) Specifically, let Σ(k) = (σij )1≤i,j≤q , k = 1, 2, . . . , K. For each off-diagonal entry i ̸= j, we first define the group-specific correlation coefficient (k) σij (k) , ρij = q (k) (k) σii σjj

k = 1, 2, . . . , K.

We then define the shared structure matrix Θ by ( θij =

K

1 X  (k) 2 ρ K k=1 ij

)1/2 ,

i ̸= j,

and set θii = 1, i = 1, 2, . . . , q. Under this rule, Θ captures the common correlation strength 23

across groups rather than the signed average correlation, and therefore avoids the cancellation of positive and negative correlations across groups. Given Θ, the heterogeneous structure matrices {Γ(k) }K k=1 are uniquely determined as follows. For the diagonal entries, we set (k)

(k)

γii = σii ,

i = 1, 2, . . . , q,

k = 1, 2, . . . , K.

For the off-diagonal entries, if θij > 0, we define (k) σij (k) , γij = θij

i ̸= j,

k = 1, 2, . . . , K. (k)

If θij = 0, then by the definition of θij , we have ρij = 0 for all k = 1, 2, . . . , K, and hence (k) (k) σij = 0 for all groups. In this case, we set γij = 0 for all k, which avoids division by zero while preserving the exact reconstruction. Therefore, the covariance matrices satisfy Σ(k) = Θ ⊙ Γ(k) ,

k = 1, 2, . . . , K.

This normalization fixes the scaling ambiguity of the Hadamard product decomposition and yields a unique decomposition. It also admits a natural interpretation: the shared structure Θ captures the common dependence strength across groups, whereas the heterogeneous structures {Γ(k) }K k=1 retain group-specific scale, sign, and deviation information. Appendix A.3. Gaussian random projection dimensionality reduction Since the deep representations extracted by pretrained models are usually of high dimension, directly estimating the covariance matrix in the original representation space leads to two practical issues. First, when the feature dimension is large while the sample size is relatively limited, the sample covariance matrix and its robust alternatives may become numerically unstable or even nearly singular, which in turn affects the accuracy of Mahalanobis distance calculation as well as the stability of the decomposition into shared and heterogeneous structures in the joint estimation procedure. Second, the proposed joint casewise/cellwise MCD methods require repeated updates of the mean vectors, covariance matrices, weight matrices, shared parameters, and group-specific parameters across multiple groups. These steps involve matrix inversion, determinant evaluation, and iterative optimization. As the representation dimension increases, both the computational complexity and memory cost grow substantially, resulting in reducing efficiency during both training and detection. Therefore, in practical implementation, it is necessary to perform an appropriate dimension reduction on the deep representations, while preserving their structural information as much as possible, so as to improve the stability and computational efficiency of the subsequent robust joint estimation procedure. Motivated by these considerations, we employ Gaussian Random Projection as a preprocessing step for the deep representations in our algorithm implementation. Specifically, let the original deep representation be z i ∈ Rq , where q may be large. We generate a random matrix R ∈ Rd×q whose entries are independently drawn from the Gaussian distribution N (0, 1/d), and project the original representation into a lower-dimensional space by Rz i ∈ Rd , d ≪ q, where d is the dimension of the projected variable. All subsequent robust mean estimation, covariance estimation, and Mahalanobis distance computation are then carried out in the projected low-dimensional space. The use of Gaussian random projection has several advantages. First, it does not require additional training or complex modeling of the original representations, and can therefore be easily incorporated into the existing detection pipeline. Second, according to the Johnson–Lindenstrauss type results, random projection preserves the Euclidean geometric structure among samples with 24

high probability, which helps maintain the relative distance relationships between normal and negative samples after dimension reduction. Third, by substantially reducing the dimension of the covariance matrix, it alleviates the numerical instability in matrix inversion and determinant computation, while also significantly reducing the computational burden of the joint optimization procedure. Appendix A.4. Optimization details for Θ and Γ(k) In this appendix, we provide the detailed numerical procedures for updating the parameters Θ and Γ(k) under ℓ1 regularization. These updates correspond to the parameter update steps described in Sections 5.1 and 5.2, namely the Fast-MCD algorithm for the joint casewise MCD estimator and the EM-based procedure for the joint cellwise MCD estimator. The appendix includes the gradient derivations, proximal gradient updates, and the definition of the softthresholding operator used in the iterative optimization. When updating Θ and Γ(k) , the algorithm requires solving optimization subproblems under an ℓ1 penalty. We employ the proximal gradient method [43] for these updates. Recall that the reparameterized form of the covariance matrix is ( (k) θij γij , i ̸= j, (k) (k) (k) Σ = (σij )1≤i,j≤q , σij = (k) i = j. γii , For the JCASEMCD, we solve the alternating optimization problem: b (t) and update each Γ b (k,t+1) to decrease the objective for k = 1, 2, . . . , K: • Fix Θ b (k,t+1) = arg min Γ Γ

nk X

(k,t) wi ℓ̃



(k,t+1)

b µ

 X (t) (k) b + λ2 |γij |. , Θ ⊙ Γ; z i

i=1

i̸=j

b (k,t+1) and update the shared structure Θ b (t+1) : • Fix all Γ b (t+1) = arg min Θ Θ

Define

nk K X X

(k,t)

wi

  X b (k,t+1) ; z (k) + λ1 b (k,t+1) , Θ ⊙ Γ |θij |. ℓ̃ µ i

k=1 i=1

i̸=j

nk X 1 (k,t) (k) (k,t) b (k,t+1) )(z (k) b (k,t+1) )⊤ wi (z i − µ Scase = Pn i −µ (k,t) k i=1 i=1 wi nk X 1 (k,t) (k) b (k,t+1) )(z (k) b (k,t+1) )⊤ . = w (z i − µ i −µ hk i=1 i

For notational clarity, define b (k,t) = Θ b (t) ⊙ Γ b (k,t) , Σ Γ

b (k,t) = Θ b (t) ⊙ Γ b (k,t+1) . Σ Θ

P PK (t) (k,t) (t) (k,t) The smooth parts of the subproblems are Qcase,Γ = K k=1 Qcase,Γ and Qcase,Θ = k=1 Qcase,Θ with h i (k,t) (k,t) −1 (k,t) (k,t) b b Qcase,Γ = hk ln |ΣΓ | + tr((ΣΓ ) Scase ) , and

h i (k,t) b (k,t) | + tr((Σ b (k,t) )−1 S (k,t) ) . Qcase,Θ = hk ln |Σ Θ Θ case

25

The non-smooth part, corresponding to ℓ1 regularization, is handled by proximal operators. From matrix calculus: h (k,t) i (k,t) −1 b (k,t) )−1 , b b (k,t) )−1 S (k,t) (Σ ∇Σb (k,t) Qcase,Γ = hk (Σ ) − ( Σ case Γ Γ Γ Γ

and

h (k,t) i (k,t) −1 b b (k,t) )−1 S (k,t) (Σ b (k,t) )−1 . ∇Σb (k,t) Qcase,Θ = hk (Σ ) − ( Σ Θ Θ case Θ Θ

b (k,t) : Applying the chain rule, for Γ h i (k,t) ∇Γb (k,t) Qcase,Γ =

( (t) (k,t) θ̂ij [∇Σb (k,t) Qcase,Γ ]ij , i ̸= j, Γ

(k,t)

[∇Σb (k,t) Qcase,Γ ]ii ,

ij

i = j.

Γ

At iteration t, we perform gradient descent: (k,t+1)

γ̃ij

(k,t)

= γ̂ij

(t)

(k,t)

(k,t+1)

, ηΓ,k λ2 ), i ̸= j, i = j,

− ηΓ,k [∇Γb (k,t) Qcase,Γ ]ij ,

followed by the proximal update: ( (k,t+1) = γ̂ij

Soft(γ̃ij (k,t+1) , γ̃ii

(t)

b (t+1) : where Soft(x, ρ) = sign(x) · max{|x| − ρ, 0}. Similarly, for Θ ( (t+1) (t) Soft(θ̃ij , ηΘ λ1 ), i ̸= j, (t+1) = θ̂ij 1, i = j, with

(t+1)

(t)

(t)

(t)

=θ̂ij − ηΘ [∇Θ b (t) Qcase,Θ ]ij , ( PK (k,t+1) (k,t) [∇Σb (k,t) Qcase,Θ ]ij , i ̸= j, (t) k=1 γ̂ij Θ [∇Θ b (t) Qcase,Θ ]ij = 0, i = j. θ̃ij

For the JCELLMCD, the two subproblems are: b (t) and µ b (k,t+1) : b (k,t+1) , update Γ • Fix Θ h i X b (k,t+1) = arg min nk ln |Θ b (t) ⊙ Γ| + tr((Θ b (t) ⊙ Γ)−1 S (k,t) ) + λ2 Γ |γij |. cell Γ

b • Fix {Γ

(k,t+1)

b b (k,t+1) , update Θ }K k=1 and µ

b (t+1) = arg min Θ Θ

K X

i̸=j (t+1)

:

h i X b (k,t+1) | + tr((Θ ⊙ Γ b (k,t+1) )−1 S (k,t) ) + λ1 nk ln |Θ ⊙ Γ |θij |. cell

k=1

i̸=j

For notational clarity, define b (t) ⊙ Γ b (k,t) , b (k,t) = Θ Σ Γ

b (k,t) = Θ b (t) ⊙ Γ b (k,t+1) . Σ Θ

26

The ℓ1 regularization is handled via the proximal operator. The smooth gradient satisfies: h (k,t) i (k,t) −1 b b (k,t) )−1 S (k,t) (Σ b (k,t) )−1 , ) − ( Σ ∇Σb (k,t) Qsmooth,Γ = nk (Σ Γ Γ Γ cell Γ

and

h (k,t) i (k,t) −1 b b (k,t) )−1 S (k,t) (Σ b (k,t) )−1 , ∇Σb (k,t) Qsmooth,Θ = nk (Σ ) − ( Σ Θ Θ Θ cell Θ

with

h i (k,t) b (k,t) | + tr((Σ b (k,t) )−1 S (k,t) ) Qsmooth,Γ = nk ln |Σ Γ Γ cell

and

(k,t) Qsmooth,Θ = nk

h

i (k,t) (k,t) −1 (k,t) b b ln |ΣΘ | + tr((ΣΘ ) Scell ) .

The chain rule gives (k,t)

[∇Γb (k,t) Qsmooth,Γ ]ij =

( (t) (k,t) θ̂ij [∇Σb (k,t) Qsmooth,Γ ]ij , i ̸= j, Γ

(k,t)

[∇Σb (k,t) Qsmooth,Γ ]ii ,

i = j.

Γ

Gradient descent and proximal update: ( (k,t+1) (k,t) (t) (k,t) γ̃ij = γ̂ij − ηΓ,k [∇Γb (k,t) Qsmooth,Γ ]ij ,

b Update for Θ

(t+1)

(t+1)

θ̂ij

=

: (

(t+1)

Soft(θ̃ij 1,

(k,t+1) γ̂ij =

(k,t+1)

Soft(γ̃ij (k,t+1) , γ̃ii

(t)

, ηΓ,k λ2 ), i ̸= j, i = j.

(t)

, ηΘ λ1 ), i ̸= j, i = j,

(t) [∇Θ b (t) Qsmooth,Θ ]ij =

(t+1)

θ̃ij

(t)

(t)

(t)

= θ̂ij − ηΘ [∇Θ b (t) Qsmooth,Θ ]ij ,

(k,t+1) (k,t) [∇Σb (k,t) Qsmooth,Θ ]ij , k=1 γ̂ij

( PK

i ̸= j,

Θ

0,

i = j.

b (k) ) ≥ ak > 0, where Σ b (k) = Θ b ⊙Γ b (k) , a backtracking To enforce the constraint λmin (Σ (t) (t) line search is employed to determine the step size ηΘ and ηΓ,k at each proximal gradient update. A candidate update is accepted only if the corresponding unprojected covariance matrices satisfy the minimum eigenvalue constraint and the penalized objective function is non-increasing; otherwise, the step size is reduced and the update is retried. In addition, to improve numerical stability, we use minimum eigenvalue clipping as a numerical safeguard when reconstructing covariance matrices and in linear algebra operations such as Cholesky factorization and log-determinant computation. Appendix B. Proofs for results Appendix B.1. Proof of Theorem 6.1 Proof. First, according to the algorithm, the update of µ corresponds to solving the maximum likelihood estimate, which minimizes the above objective function and therefore does not increase it. Specifically, fixing H (k,t) and Σ(k,t) , where Σ(k,t) = Θ(t) ⊙ Γ(k,t) , the terms related to µ(k) in the k-th group are Lµ(k) =

nk X

(k,t)

wi

(k)

(k)

(z i − µ(k) )⊤ (Σ(k,t) )−1 (z i − µ(k) ).

i=1

27

Taking the gradient with respect to µ(k) yields: ∇µ(k) Lµ(k) = −2(Σ(k,t) )−1

nk X

(k,t)

wi

(k)

(z i − µ(k) ).

i=1

Setting the gradient to zero gives nk X

(k,t)

wi

(k)

(z i − µ(k) ) = 0,

i=1

which implies (k,t+1)

µ

nk (k,t) (k) zi 1 X (k,t) (k) i=1 wi wi z i . Pnk (k,t) = h k i=1 i=1 wi

Pnk =

Moreover, the Hessian matrix is ∇2µ(k) Lµ(k) = 2hk (Σ(k,t) )−1 ≻ 0. Since this is a strictly convex quadratic function, the above provides the unique global minimum. Therefore, Lcase (µ(t+1) , Σ(t) , H (t) ) ≤ Lcase (µ(t) , Σ(t) , H (t) ), (t) (k,t) where µ(t+1) = {µ(k,t+1) }K = {Σ(k,t) }K = Θ(t) ⊙ Γ(k,t) . k=1 , Σ k=1 and Σ For the update of Γ(k) , at iteration t, fixing µ(k,t+1) , H (k,t) , and Θ(t) , define the smooth part of the current subproblem as   (k,t) (t) (k) (k,t) (k,t+1) . Qcase,Γ (Γ(k) ) = Q(k) µ , Θ ⊙ Γ , H case

The Hessian of this expression consists solely of matrix products involving (Σ(k) )−1 and the (k,t) fixed sample covariance Scase , given µ(k,t+1) , H (k,t) and Θ(t) . Moreover, the relationship be(k) (k) tween Γ and Σ is linear, since Σ(k) = Θ(t) ⊙ Γ(k) . Within the feasible region, the constraint λmin (Σ(k) ) ≥ a > 0 guarantees that the spectral norm of the inverse matrix is strictly bounded, (k,t) i.e., ∥(Σ(k) )−1 ∥2 ≤ 1/a. Since ∥Scase ∥2 and the entries of Θ(t) are fixed for the current subproblem, all Hessian components are bounded locally. Therefore, there exists a finite local (t) (k,t) Lipschitz constant LΓ,k such that ∇Γ(k) Qcase,Γ is Lipschitz continuous in the current subproblem (Lemma 1.2.2 in Nesterov [44]). Hence, for any X, Y in the current feasible region, (t)

(k,t)

(k,t)

(k,t)

Qcase,Γ (Y ) ≤ Qcase,Γ (X) + ⟨∇Γ(k) Qcase,Γ (X), Y − X⟩ +

LΓ,k ∥Y − X∥2F . 2

Similarly, for the update of Θ, fixing µ(k,t+1) , H (k,t) , and Γ(k,t+1) , define (t) Qcase,Θ (Θ) =

K X

Q(k) case

  (k,t+1) (k,t) (k,t+1) µ ,Θ ⊙ Γ ,H .

k=1

Since the mapping Θ 7→ Θ ⊙ Γ(k,t+1) is linear for each fixed Γ(k,t+1) , and since Γ(k,t+1) is fixed in the current Θ-subproblem, the same argument implies that there exists a finite local Lipschitz (t) (t) constant LΘ such that ∇Θ Qcase,Θ is Lipschitz continuous in the current subproblem. Let (t)

(k,t)

Fk (Γ(k) ) = Qcase,Γ (Γ(k) ) + λ2 ∥Γ(k) ∥1,off

28

and

(t)

F (t) (Θ) = Qcase,Θ (Θ) + λ1 ∥Θ∥1,off . (t)

(t)

By the proximal gradient lemma [45], when 0 < ηΓ,k ≤ (1 − δΓ )/LΓ,k , applying this inequality with X = Γ(k,t) , Y = Γ(k,t+1) and using the optimality of the proximal gradient update yields ! (t) L 1 Γ,k (t) (t) Fk (Γ(k,t+1) ) ≤ Fk (Γ(k,t) ) − ∥Γ(k,t+1) − Γ(k,t) ∥2F . − (t) 2 2ηΓ,k (t)

(t)

Similarly, when 0 < ηΘ ≤ (1 − δΘ )/LΘ , we have F (t) (Θ(t+1) ) ≤ F (t) (Θ(t) ) −

(t)

1

L − Θ (t) 2 2ηΘ

! ∥Θ(t+1) − Θ(t) ∥2F .

Since the coefficients of the squared norm terms are nonnegative under the above local step-size conditions, the updates of Γ(k) and Θ do not increase the objective function. In summary, Lcase (Ω(t+1) , H (t) ) ≤ Lcase (Ω(t) , H (t) ). Finally, the update rule for H ensures that the objective function does not increase. Given Ω, the terms related to H (k) in the k-th group are LH (k) =

nk X

(k) wi

h i (k) (k) (k) 2 ln |Σ | + q ln(2π) + MD (z i , Ω ) ,

i=1

P k (k) (k) (k) (k) wi = hk is fixed, where MD2 (z i , Ω(k) ) = (z i − µ(k) )⊤ (Σ(k) )−1 (z i − µ(k) ). Since ni=1 ln |Σ(k) | + q ln(2π) acts as a constant. Therefore, minimizing L(k) is equivalent to minimizing P k (k) Pnk (k) (k) (k) (k) 2 wi = hk . The solution to this 0-1 linear ∈ {0, 1} and ni=1 i=1 wi MD (z i , Ω ) with wi programming problem is to select the hk samples with the smallest Mahalanobis distances. Algorithm 1 performs exactly this. Hence, Lcase (Ω(t+1) , H (t+1) ) ≤ Lcase (Ω(t+1) , H (t) ). Thus,

Lcase (Ω(t+1) , H (t+1) ) ≤ Lcase (Ω(t+1) , H (t) ) ≤ Lcase (Ω(t) , H (t) ),

establishing the first property. Moreover, let a = min1≤k≤K {ak }. For any k, given λmin (Σ(k) ) ≥ Q a, we have |Σ(k) | = qj=1 λj (Σ(k) ) ≥ aq , which implies ln |Σ(k) | ≥ q ln a. The Mahalanobis distance is non-negative: (k)

(k)

(z i − µ(k) )⊤ (Σ(k) )−1 (z i − µ(k) ) ≥ 0. The ℓ1 regularization term is also non-negative. Hence, Lcase (Ω, H) ≥

K X k=1

Q(k) case ≥

K X

hk [q ln a + q ln(2π)] := Lmin > −∞.

k=1

By the monotone convergence theorem, since the update sequence of Lcase is monotonically non-increasing and bounded below, it converges to a finite limit L∗ , establishing the second property.

29

For the third property, we first show the iterative sequence has a limit point. Since Lcase (Ω(t) , H (t) ) ≤ Lcase (Ω(0) , H (0) ) := C0 , all iterates are contained in the initial sublevel set. Moreover, the Mahalanobis-distance terms and the penalty term are nonnegative. Therefore, for each k, X hk ln |Σ(k) | ≤ C0 − hk q ln(2π) − hl · q{ln a + ln(2π)} =: Ck . l̸=k

Thus, |Σ(k) | ≤ exp(Ck /hk ). Since all eigenvalues of Σ(k) are bounded below by a, we have (k)

|Σ | =

q Y

λj (Σ(k) ) ≥ aq−1 λmax (Σ(k) ).

j=1

Consequently, λmax (Σ(k) ) ≤ a−(q−1) exp(Ck /hk ), which shows that {Σ(k,t) }t is uniformly bounded for every group k. Together with the lower eigenvalue constraint, this also implies that {(Σ(k,t) )−1 }t is uniformly bounded. We next claim that the location estimates are bounded. Let (k)

Rk = max ∥z i ∥2 < ∞. 1≤i≤nk

For any selected subset H (k,t) of size hk , nk X

(k,t)

wi

(k)

(k)

(z i − µ(k,t) )⊤ (Σ(k,t) )−1 (z i − µ(k,t) ) ≥

i=1

nk X 1 (k) (k,t) wi ∥z i − µ(k,t) ∥22 . (k,t) λmax (Σ ) i=1

(k)

Using ∥z i − µ(k,t) ∥2 ≥ |∥µ(k,t) ∥2 − Rk |, we obtain nk X

(k,t)

wi

(k)

(k)

(z i − µ(k,t) )⊤ (Σ(k,t) )−1 (z i − µ(k,t) ) ≥

i=1

1 (∥µ(k,t) ∥2 − Rk )2 . (k,t) λmax (Σ )

Since the left-hand side is bounded above on the same sublevel set and λmax (Σ(k,t) ) is uniformly bounded, it follows that {µ(k,t) }t is bounded. Finally, we prove the boundedness of the reparameterized variables. The lower bound on the likelihood part implies that the penalty term is uniformly bounded on the sublevel set: X λ1 ∥Θ(t) ∥1,off + λ2 ∥Γ(k,t) ∥1,off ≤ C0 − Lmin < ∞, k

Hence the off-diagonal entries of Θ(t) and Γ(k,t) are bounded. The diagonal entries of Θ(t) are fixed at one, while the diagonal entries of Γ(k,t) coincide with those of Σ(k,t) and are therefore bounded by the uniform bound on λmax (Σ(k,t) ). Thus Θ(t) and Γ(k,t) are bounded. Therefore, the sequence {Ω(t) } is bounded. Furthermore, since H (t) only takes a finite number of 0-1 combinations, by the Bolzano–Weierstrass theorem [46], there exists a convergent subsequence (Ω(tr ) , H (tr ) ) → (Ω∗ , H ∗ ), where tr is the subsequence index, r ∈ N+ . Moreover, since {Ω(t) } is bounded and the feasible region satisfies λmin (Σ(k,t) ) ≥ ak > 0, all iterates lie in a compact subset of the feasible region. On this compact set, the gradients of the smooth parts of the Γ(k) - and Θ-subproblems are globally Lipschitz continuous. Hence,

30

there exist finite constants LΓ,k and LΘ such that, for all iterations t, (t)

(t)

LΓ,k ≤ LΓ,k ,

LΘ ≤ LΘ .

When the step sizes are chosen with a uniform margin, i.e., (t)

0 < ηΓ,k ≤

1 − δΓ , LΓ,k

(t)

0 < ηΘ ≤

1 − δΘ , LΘ

for some δΓ , δΘ ∈ (0, 1), the proximal gradient descent inequalities imply the uniform sufficient decrease property. Specifically, there exist constants cΓ,k > 0 and cΘ > 0 such that (t)

(t)

Fk (Γ(k,t+1) ) ≤ Fk (Γ(k,t) ) − cΓ,k ∥Γ(k,t+1) − Γ(k,t) ∥2F , and

F (t) (Θ(t+1) ) ≤ F (t) (Θ(t) ) − cΘ ∥Θ(t+1) − Θ(t) ∥2F .

Summing these inequalities over t and using the fact that Lcase is bounded below, we obtain X X ∥Θ(t+1) − Θ(t) ∥2F < ∞. ∥Γ(k,t+1) − Γ(k,t) ∥2F < ∞, t

t

Hence,

Γ(k,t+1) − Γ(k,t) → 0,

Θ(t+1) − Θ(t) → 0.

Moreover, since the step sizes are bounded away from zero, i.e., (t)

(s)

ηΓ,k ≥ inf ηΓ,k > 0, s

it follows that

(t)

(s)

ηΘ ≥ inf ηΘ > 0, s

 1  (k,t+1) (k,t) Γ −Γ → 0, (t)

ηΓ,k

1  (t)

Θ

(t+1)

t ≥ 0.

−Θ

(t)



→ 0.

ηΘ

For the proximal optimality condition of the Γ-update: (k,t)

0 ∈ ∇Γ(k) Qcase,Γ (Γ(k,t) ) +

1 (t) ηΓ,k

(Γ(k,t+1) − Γ(k,t) ) + λ2 ∂∥Γ(k,t+1) ∥1,off .

Let (Ω∗ , H ∗ ) be an accumulation point satisfying the interior condition λmin (Θ∗ ⊙ Γ(k)∗ ) > ak ,

k = 1, 2, . . . , K.

By the continuity of the minimum eigenvalue, the eigenvalue constraints are inactive in a sufficiently small neighborhood of this accumulation point. Hence, along the convergent subsequence, the accepted proximal gradient updates coincide locally with the unconstrained proximal gradient updates for the corresponding penalized block subproblems. Taking the limit along the convergent subsequence and using the continuity of the gradient, the vanishing of the scaled difference term, and the closedness of the subdifferential of the ℓ1 norm yields ∗ ∗ (k)∗ 0 ∈ ∇Γ(k) Q(k) ∥1,off . case (Ω , H ) + λ2 ∂∥Γ

31

By the same reasoning for the Θ-update: 0 ∈ ∇Θ

K X

∗ ∗ ∗ Q(k) case (Ω , H ) + λ1 ∂∥Θ ∥1,off .

k=1

Meanwhile, since the µ-update is an exact minimization: 0 = ∇µ(k) Lcase (Ω∗ , H ∗ ). Finally, since the H-update selects the hk samples with the smallest Mahalanobis distances given the current parameter estimates, H ∗ belongs to the set of minimizers of the discrete H-subproblem given Ω∗ ; in the presence of ties, the conclusion is understood under the same tie-breaking rule used in the algorithm. Hence, any interior accumulation point (Ω∗ , H ∗ ) is blockwise stationary with respect to the continuous variables, and H ∗ is optimal for the discrete block. This completes the proof of Theorem 6.1. Appendix B.2. Proof of Theorem 6.2 (k) Proof. First, given Ω, for the k-th group, i-th row, and j-th column, let o = {τ ̸= j : wiτ = 1} denote the indices of the currently retained cells. By changing the status of cell j from flagged to retained, the Schur complement yields: (k)

(k)

|Σo∪{j} | = |Σ(k) oo | · Cij , (k)

(k)

(k)

(k)

−1 ⊤ where Cij = Σjj − Σj,o (Σ(k) oo ) (Σj,o ) . Meanwhile, the Mahalanobis distance can be decomposed as: (k) (k) (zij − zbij )2 2 2 MDo∪{j} = MDo + , (k) Cij

where

(k)

(k)

(k)

(k)

(k) −1 zbij = µj + Σj,o (Σ(k) oo ) (z i,o − µo ). (k)

(k)

(k)

(k)

2 When o = ∅, we use the convention that |Σ(k) bij = µj , and Cij = Σjj . oo | = 1, MDo = 0, z Therefore, compared to keeping it “flagged”, the increment in the objective function by retaining cell (i, j) is exactly: (k) (k) (zij − zbij )2 (k) (k) (k) ∆ij = ln Cij + ln(2π) + − bj . (k) Cij

This matches the quantity defined in Algorithm 3. Consequently, fixing the other parameters, the optimization over the j-th column in the k-th group is equivalent to: min

nk X

(k) (k) w1j ,...,wn j ∈{0,1} i=1 k

(k) (k) wij ∆ij

s.t.

nk X

(k)

wij ≥ hk .

i=1 (k)

This 0–1 linear programming problem admits the following solution: all cells with ∆ij < 0 should be retained, and if their number is less than hk , the remaining cells are selected in (k) ascending order of ∆ij until the columnwise lower bound hk is satisfied. The thresholding (k) rule in Algorithm 3 is precisely this optimal solution, where ∆(hk )j is the hk -th smallest value (k)

k in {∆ij }ni=1 . Thus, each columnwise weight update does not increase the objective function. Applying this argument sequentially over all columns j = 1, . . . , q and all groups k = 1, . . . , K

32

gives

Lcell (Ω(t) , W (t+1) ) ≤ Lcell (Ω(t) , W (t) ).

Given W (t+1) , the update of Ω can be viewed as a generalized EM step with penalty. We treat the flagged cells as missing variables and the retained cells as observed variables. Let Lsmooth (Ω, W (t+1) ) denote the observed negative log-likelihood part of Lcell (Ω, W (t+1) ), excludcell ing the cellwise flagging penalty and the ℓ1 regularization terms. For the current parameter Ω(t) , the EM majorization relationship gives (Ω(t) , W (t+1) ) Lsmooth (Ω, W (t+1) ) − Lsmooth cell cell ≤Qsmooth (Ω | Ω(t) , W (t+1) ) − Qsmooth (Ω(t) | Ω(t) , W (t+1) ), where terms independent of Ω are omitted from the surrogate. Define the penalized EM surrogate (t)

Qpen (Ω | Ω , W

(t+1)

(t)

) = Qsmooth (Ω | Ω , W

(t+1)

) + λ1 ∥Θ∥1,off + λ2

K X

∥Γ(k) ∥1,off .

k=1

Since the cellwise flagging penalty is fixed when W (t+1) is fixed, adding the regularization terms to both sides yields Lcell (Ω, W (t+1) ) − Lcell (Ω(t) , W (t+1) ) ≤ Qpen (Ω | Ω(t) , W (t+1) ) − Qpen (Ω(t) | Ω(t) , W (t+1) ). Therefore, it is sufficient to find a feasible update Ω(t+1) satisfying Qpen (Ω(t+1) | Ω(t) , W (t+1) ) ≤ Qpen (Ω(t) | Ω(t) , W (t+1) ), which guarantees

Lcell (Ω(t+1) , W (t+1) ) ≤ Lcell (Ω(t) , W (t+1) ).

The update of µ is the exact minimizer of the smooth EM surrogate with respect to the location parameters. Specifically, it is given by the sample mean of the conditional completedata means. Hence, the µ-step does not increase Qpen . Next, fixing µ(t+1) , W (t+1) , and Θ(t) , the Γ(k) -subproblem has the form (k,t)

min Qsmooth,Γ (Γ(k) ) + λ2 ∥Γ(k) ∥1,off , Γ

(k)

(k,t)

(k,t)

where Scase in the proof of Theorem 6.1 is replaced by Scell . According to the proximal gradient update with backtracking described in Appendix A.4, the accepted update is feasible, satisfies λmin (Θ(t) ⊙ Γ(k,t+1) ) ≥ ak and yields (k,t)

(k,t)

Qsmooth,Γ (Γ(k,t+1) ) + λ2 ∥Γ(k,t+1) ∥1,off ≤ Qsmooth,Γ (Γ(k,t) ) + λ2 ∥Γ(k,t) ∥1,off . Similarly, fixing µ(t+1) , W (t+1) , and {Γ(k,t+1) }K k=1 , the Θ-subproblem is min Θ

K X

(k,t)

Qsmooth,Θ (Θ) + λ1 ∥Θ∥1,off .

k=1

The accepted proximal gradient update with backtracking is feasible, satisfies λmin (Θ(t+1) ⊙

33

Γ(k,t+1) ) ≥ ak , k = 1, 2, . . . , K, and gives K X

(k,t) Qsmooth,Θ (Θ(t+1) ) + λ1 ∥Θ(t+1) ∥1,off ≤

k=1

K X

(k,t)

Qsmooth,Θ (Θ(t) ) + λ1 ∥Θ(t) ∥1,off .

k=1

Combining the exact µ-update with the accepted proximal gradient updates for Γ(k) and Θ, we obtain Qpen (Ω(t+1) | Ω(t) , W (t+1) ) ≤ Qpen (Ω(t) | Ω(t) , W (t+1) ). Therefore,

Lcell (Ω(t+1) , W (t+1) ) ≤ Lcell (Ω(t) , W (t+1) ).

Combining this with the W -update from the first step establishes the monotonicity of the overall algorithm: Lcell (Ω(t+1) , W (t+1) ) ≤ Lcell (Ω(t) , W (t) ). (k)

For any nonempty observed pattern A ⊂ {1, . . . , q}, ΣA is a principal submatrix of Σ(k) . Let a = min1≤k≤K {ak }. By Cauchy’s interlace theorem (Theorem 4.3.28, [47]), (k)

λmin (ΣA ) ≥ λmin (Σ(k) ) ≥ a. Therefore,

(k)

ln |ΣA | ≥ |A| ln a > −∞. (k)

For the empty pattern A = ∅, we use the convention ln |ΣA | = 0 and MD2A = 0. The partial Mahalanobis distance is non-negative, the discrete penalty satisfies q K X X (k) (k) bj ∥1 − W ·j ∥0 ≥ 0, k=1 j=1

and the ℓ1 regularization terms are also non-negative. Consequently, Lcell has a global lower bound. Since it decreases monotonically, the sequence Lcell (Ω(t) , W (t) ) converges to a finite value, establishing the second property. We next show that the generated sequence has a convergent subsequence. Since Lcell (Ω(t) , W (t) ) ≤ Lcell (Ω(0) , W (0) ) =: C0 , all iterates lie in the initial sublevel set. Fix a group k and a variable j. The constraint (k,t) ∥W ·j ∥0 ≥ hk implies that at least one retained pattern contains the variable j. For any such (k,t) nonempty pattern A with j ∈ A, the corresponding log-determinant term ln |ΣA | is bounded above on the sublevel set, because all other terms in the objective are bounded below. Together (k,t) with λmin (ΣA ) ≥ a, we obtain (k,t) |ΣA | =

|A| Y

(k,t)

(k,t)

λl (ΣA ) ≥ a|A|−1 λmax (ΣA ),

l=1 (k,t)

(k,t)

which implies that λmax (ΣA ) is uniformly bounded. Hence the diagonal element σjj is uniformly bounded. Since this holds for every j, all diagonal elements of Σ(k,t) are uniformly bounded. As Σ(k,t) is positive semidefinite, every 2 × 2 principal submatrix is also q positive (k,t) (k,t)

(k,t)

(k,t)

(k,t) (k,t)

semidefinite. Hence, for any j ̸= l, σjj σll −(σjl )2 ≥ 0, which implies |σjl | ≤ σjj σll . Therefore, the boundedness of all diagonal entries implies the boundedness of all off-diagonal 34

entries. Thus {Σ(k,t) }t is uniformly bounded for every k. The lower eigenvalue constraint also implies that {(Σ(k,t) )−1 }t is uniformly bounded. We now prove that the location estimates are bounded. Let (k)

Rk = max ∥z i ∥2 < ∞. 1≤i≤nk

For each coordinate j, choose a retained pattern A containing j. The corresponding partial Mahalanobis term satisfies (k)

(k)

(k)

(k)

(k)

(z i,A − µA )⊤ (ΣA )−1 (z i,A − µA ) ≥

1 (k) λmax (ΣA )

(k)

(k)

∥z i,A − µA ∥22 .

(k)

Since λmax (ΣA ) is uniformly bounded and the objective is bounded above on the sublevel set, (k,t) µj cannot diverge. Hence {µ(k,t) }t is bounded for every k. Finally, the lower bound on the likelihood part implies that the penalty term is uniformly bounded on the sublevel set. Thus the off-diagonal entries of Θ(t) and Γ(k,t) are bounded. The diagonal entries of Θ(t) are fixed at one, while the diagonal entries of Γ(k,t) coincide with those of Σ(k,t) and are therefore bounded. Hence Θ(t) , Γ(k,t) , Σ(k,t) , and µ(k,t) are bounded. Since W (t) only takes values from a finite set of 0-1 matrices, by the Bolzano–Weierstrass theorem there exists a convergent subsequence (Ω(tr ) , W (tr ) ) → (Ω∗ , W ∗ ), where tr ∈ N+ is the index of the subsequence. By the same sufficient-decrease argument as in the proof of (k,t) (k,t) Theorem 6.1, with Scase replaced by Scell , the accepted proximal gradient updates imply X X ∥Γ(k,t+1) − Γ(k,t) ∥2F < ∞, ∥Θ(t+1) − Θ(t) ∥2F < ∞. t

Therefore,

t

Γ(k,t+1) − Γ(k,t) → 0,

Θ(t+1) − Θ(t) → 0.

The proximal optimality condition for the Γ(k) -update gives  1  (k,t+1) (k,t) (k,t) (k,t) −Γ + λ2 ∂∥Γ(k,t+1) ∥1,off . 0 ∈ ∇Γ(k) Qsmooth,Γ (Γ ) + (t) Γ ηΓ,k Let (Ω∗ , W ∗ ) be an accumulation point satisfying the interior condition λmin (Θ∗ ⊙ Γ(k)∗ ) > ak ,

k = 1, 2, . . . , K.

By the continuity of the minimum eigenvalue, the eigenvalue constraints are inactive in a sufficiently small neighborhood of this accumulation point. Hence, along the convergent subsequence, the accepted proximal gradient updates coincide locally with the unconstrained proximal gradient updates for the corresponding penalized block subproblems. Taking the limit along the convergent subsequence and using the continuity of the gradient, the vanishing of the scaled difference term, and the closedness of the subdifferential of the ℓ1 norm, we obtain (k) 0 ∈ ∇Γ(k) Qsmooth (Ω∗ | Ω∗ , W ∗ ) + λ2 ∂∥Γ(k)∗ ∥1,off . Similarly, 0 ∈ ∇Θ

K X

(k)

Qsmooth (Ω∗ | Ω∗ , W ∗ ) + λ1 ∂∥Θ∗ ∥1,off .

k=1

35

Since the µ-update is an exact minimization of the smooth EM surrogate, 0 = ∇µ(k) Qsmooth (Ω∗ | Ω∗ , W ∗ ). At Ω = Ω∗ , the EM surrogate and the observed smooth objective have the same first-order derivatives with respect to the continuous variables. Hence the above conditions are the firstorder stationarity conditions of the penalized observed objective with W ∗ fixed. Moreover, the W -update is coordinatewise optimal for the discrete variables. Therefore, W ∗ is coordinatewise optimal given Ω∗ , with ties interpreted according to the same tie-breaking convention used in the algorithm. Hence every interior accumulation point (Ω∗ , W ∗ ) is coordinatewise optimal with respect to W and stationary with respect to the continuous variables. This completes the proof of Theorem 6.2. Appendix B.3. Proof of Theorem 6.3 Proof. First, we prove the existence of a clean benchmark with a finite objective value. Suppose that the number of contaminated cells in each column of each group does not exceed nk − hk . Then, in every column of group k, at least hk cells remain uncontaminated. Hence, for each (k) group k, we can construct a feasible benchmark mask W 0 by setting all contaminated cells (k) to 0 and all uncontaminated cells to 1. This mask satisfies ∥W 0,·j ∥0 ≥ hk , j = 1, 2, . . . , q, k = 1, 2, . . . , K. Choose a fixed feasible joint covariance parameter, for instance Θ0 = I q , (k) (k) (k) (k) Γ0 = ck I q , ck ≥ ak , Σ0 = Θ0 ⊙ Γ0 = ck I q . Also take any finite location vector µ0 , for (k) example µ0 = 0. Since the benchmark mask retains only clean cells, the observed likelihood (k) (k) K (k) K terms evaluated at ({µ0 }K k=1 , Θ0 , {Γ0 }k=1 , {W 0 }k=1 ) depend only on the original clean data (k) and on fixed constants. Moreover, P(Θ0 , {Γ0 }K k=1 ) < ∞ and the cellwise penalties are finite for (k) the fixed nonnegative constants bj . Therefore, there exists a constant M < ∞, independent of the magnitudes of the contaminated cells, such that the optimal joint cellwise objective value is at most M . This is because the objective function value corresponding to the optimal solution will not exceed the objective function value corresponding to the above feasible solutions. Thus, the required finite benchmark is provided. Second, we prove the scatter implosion breakdown point. Since every feasible solution of   (k) b ≥ ak > 0, k = 1, 2, . . . , K, the scatter matrix the joint cellwise estimator satisfies λmin Σ of any group cannot implode for any admissible contaminated sample. Hence, b Z ⋆ ) = 1. εcell,− (Ω, Third, we prove the scatter explosion breakdown point by contradiction. Suppose that there exists a sequence of contaminated samples, with at most nk − hk contaminated cells in each column of each group, such that  (k0 )  b λmax Σ →∞ for some group k0 . For the optimal cellwise mask in this group, define Ai =

n o (k ) j:w bij 0 = 1 ,

b (k0 ) for the corresponding principal submatrix. When Ai = ∅, we use qi = |Ai |, and write Σ Ai ,Ai b (k0 ) | = 0. Let j∗ ∈ arg max1≤j≤q σ̂ (k0 ) . Since Σ b (k0 ) is positive the standard convention ln |Σ Ai ,Ai jj semidefinite,  (k0 )  (k ) b λmax Σ ≤ q max σ̂jj 0 . 1≤j≤q

36

Hence, (k )

σ̂j∗ j0∗ ≥

 (k0 )  b λmax Σ .

q

(k0 )

(k )

c ∥0 ≥ hk0 , there exists at least one row i∗ such that w By the feasibility constraint ∥W bi∗ j0∗ = 1, ·j∗ or equivalently j∗ ∈ Ai∗ . Therefore,  (k0 )  b  (k0 )  λmax Σ (k0 ) b λmax ΣAi∗ ,Ai∗ ≥ σ̂j∗ j∗ ≥ . q Moreover, by the Cauchy interlacing theorem,  (k0 )   (k0 )  b b λmin Σ ≥ λ ≥ ak0 min Σ Ai ,Ai whenever Ai ̸= ∅. Thus, ln

b (k0 ) Σ Ai∗ ,Ai∗

=

qi∗ X

 (k0 )  b ln λl Σ Ai∗ ,Ai∗

l=1

 (k0 )  b ≥ ln λmax Σ Ai∗ ,Ai∗ + (qi∗ − 1) ln ak0  (k0 )  b λmax Σ + (qi∗ − 1) ln ak0 . ≥ ln q For every other row i, we also have (k0 )

b ln Σ Ai ,Ai ≥ qi ln ak0 with the above convention when qi = 0. Consequently, for a finite constant Ck0 = (nk0 q − 1) min{ln ak0 , 0}, depending only on nk0 , q, and ak0 , we obtain nk0 X

(k0 )

b ln Σ Ai ,Ai ≥ ln

 (k0 )  b λmax Σ q

i=1

+ Ck0 .

Since the partial Mahalanobis distances, the cellwise penalty, and the ℓ1 penalty are all nonb (k0 ) ) → ∞. This contradicts negative, the whole joint cellwise objective tends to +∞ if λmax (Σ the finite benchmark M established in the first step. Therefore, b Z ⋆ ) ≥ β0 . εcell,+ (Ω, Finally, we prove the location breakdown point. Suppose again that the number of contaminated cells per column in each group does not exceed nk − hk . Consider a sequence of such contaminated samples and assume, for contradiction, that b (k0 ) µ for some group k0 .

→∞ 2

From the scatter explosion result proved above, along this sequence 37

b λmax (Σ

(k0 )

) cannot diverge. Hence, there exists a finite constant CΣ such that  (k0 )  b λmax Σ ≤ CΣ

along the sequence. Let M0 be a finite constant such that the absolute values of all clean cells in group k0 are bounded by M0 . For each coordinate j, at least hk0 cells in column j are clean, while every feasible mask retains at least hk0 cells in the same column. Since hk0 > nk0 /2, these two sets must intersect. (k0 ) Thus, for each j = 1, 2, . . . , q, there exists an index i(j) such that the cell zi(j)j is both clean and retained. Therefore, ⋆(k ) (k0 ) zi(j)j0 ≤ M0 , w bi(j)j = 1. Using the lower bound on the partial Mahalanobis terms, we have nk0 X

nk

2

MD



(k0 ) ⋆(k ) 0) b b (k z i,Ai0 , µ Ai , ΣAi ,Ai



i=1

0 2 1 X ⋆(k ) 0) b (k z i,Ai0 − µ ≥ Ai CΣ i=1 2 q 2 1 X  ⋆(k0 ) (k0 ) ≥ z −µ bj CΣ j=1 i(j)j   2 1 1 (k0 ) 2 b ≥ µ − qM0 . CΣ 2 2

b (k0 ) ) ≥ The log-determinant terms are bounded from below by a finite constant because λmin (Σ ak0 . Hence, if ∥b µ(k0 ) ∥2 → ∞, the joint cellwise objective must tend to +∞, contradicting the finite benchmark M . Therefore, b Z ⋆ ) ≥ β0 . ε∗cell,µ (Ω, To show that the common lower bound is sharp, consider k∗ ∈ arg min k

nk − hk + 1 . nk

In this group, fix a column, say j∗ = 1, replace exactly nk∗ − hk∗ + 1 cells in this column by a value c, and let c → ∞, leaving all other cells unchanged. Let I be the set of rows whose first-column cell is contaminated. Since |I| = nk∗ − hk∗ + 1, any feasible mask satisfying (k ) ∥W ·1 ∗ ∥0 ≥ hk∗ must retain at least one contaminated c-cell in the first column, that is, n o  (k∗ ) # i:w bi1 = 1 ∩ I ≥ 1.

38

At convergence of the EM µ-block, the first coordinate satisfies nk

∗ 1 X (k ) m b i1 ∗ nk∗ i=1 X (k ) 1 1 = z̃i1 ∗ + nk∗ (k∗ ) nk∗

(k ) µ b1 ∗ =

i:w bi1

=

=1

 (k∗ ) −1   X  (k ) (k∗ ) (k∗ ) (k∗ ) ∗ b b b Oi µ b1 + Σ1,Oi ΣOi ,Oi z̃ i,Oi − µ (k )

i:w bi1 ∗ =0

n o  X 1 c (k ) ⋆(k ) # i:w bi1 ∗ = 1 ∩ I + zi1 ∗ nk∗ nk∗ (k ) {i:w bi1 ∗ =1}∩I c   (k∗ ) −1   X (k∗ ) 1 (k∗ ) ⋆(k ) (k ) ∗ ∗ b b b Oi + µ b1 + Σ z i,Oi − µ , 1,Oi ΣOi ,Oi nk∗ (k∗ ) i:w bi1

=0

(k )

where Oi = {j : w bij ∗ = 1, j ̸= 1}, and the usual empty-set convention is used if Oi = ∅. If b (k∗ ) remained bounded as c → ∞, then the left-hand side and all terms on b (k∗ ) and Σ both µ the right-hand side except the first one would remain bounded. The inverse submatrices are b (k∗ ) ) ≥ ak∗ . However, the first term on the right-hand side diverges to bounded because λmin (Σ +∞, a contradiction. Therefore, at contamination fraction (nk∗ − hk∗ + 1)/nk∗ = β0 , either the location or the scatter matrix, or both, must break down. This proves that the common lower bound β0 for scatter explosion and location breakdown is sharp. Appendix B.4. Proof of Theorem 6.4 Proof. The proof runs in parallel with the cellwise case, replacing the cellwise masks W (k) with (k) rowwise weights wi . First, we prove the existence of a clean benchmark with a finite objective value. Suppose that the number of casewise contaminated rows in group k does not exceed (k) nk − hk . Then there exists a subset H 0 of size hk consisting entirely of clean rows. Choose, for example, Θ0 = I q ,

(k)

Γ0 = ck I q ,

ck ≥ ak ,

and (k)

µ0 =

(k)

(k)

Σ0 = Θ0 ⊙ Γ0 = ck I q ,

1 X ⋆(k) zi . hk (k) i∈H 0 (k)

(k)

(k)

K K Evaluating the joint casewise objective at ({µ0 }K k=1 , Θ0 , {Γ0 }k=1 , {H 0 }k=1 ) gives a finite value depending only on the clean data and fixed constants, not on the magnitudes of the (k) contaminated rows. Also, P(Θ0 , {Γ0 }K k=1 ) < ∞. Hence, there exists a finite constant M < ∞ such that, under at most nk − hk contaminated rows in every group, the optimal joint casewise objective value is at most M . Second, the implosion breakdown value equals 1. Indeed, every feasible solution satisfies  (k)  b λmin Σ ≥ ak > 0, k = 1, 2, . . . , K.

Therefore, no scatter matrix can implode. Under the standard convention that a breakdown fraction is reported within [0, 1], this gives b Z ⋆ ) = 1. δcase,− (Ω, Third, we prove the scatter explosion breakdown point. Suppose, for the sake of contradic39

tion, that there exists a sequence of contaminated samples, with at most nk − hk contaminated rows in each group, such that  (k0 )  b λmax Σ →∞ b for some group k0 . Since λmin (Σ b (k0 ) = ln Σ

q X

(k0 )

) ≥ ak0 , we have

 (k0 )   (k0 )  b b ln λj Σ ≥ ln λmax Σ + (q − 1) ln ak0 .

j=1

Therefore, the log-determinant contribution of group k0 in the joint casewise objective is bounded from below by i h  (k0 )  b + (q − 1) ln ak0 . hk0 ln λmax Σ Since the Mahalanobis terms and the ℓ1 penalty are nonnegative, the joint casewise objective b (k0 ) ) → ∞. This contradicts the finite benchmark M established in the tends to +∞ as λmax (Σ first step. Hence, b Z ⋆ ) ≥ β0 . δcase,+ (Ω, Finally, we prove the location breakdown point. Suppose, for contradiction, that under at most nk −hk contaminated rows in each group, there exists a sequence of contaminated samples such that b (k0 ) → ∞ µ 2

b (k0 ) ) cannot diverge for some group k0 . From the scatter explosion result proved above, λmax (Σ along this sequence. Hence, there exists a finite constant CΣ such that  (k0 )  b ≤ CΣ . λmax Σ Any feasible subset of size hk0 must contain at least one clean row, because at most nk0 − hk0 ⋆(k ) rows are contaminated and hk0 > nk0 /2. Let z i0 0 be such a retained clean row, and let ⋆(k0 )

R0 = max z i i

2

over the clean rows of group k0 . Then ⋆(k )

z i0 0

2

≤ R0 .

The Mahalanobis term corresponding to this retained clean row satisfies  2 (k0 )  ⊤  (k0 ) −1   b µ − R0 ⋆(k ) ⋆(k ) 2 b b (k0 ) b (k0 ) ≥ Σ z i0 0 − µ z i0 0 − µ . CΣ Thus, if ∥b µ(k0 ) ∥2 → ∞, the joint casewise objective must tend to +∞, contradicting the finite benchmark M . Therefore, ∗ b Z ⋆ ) ≥ β0 . (Ω, δcase,µ To show that the common lower bound is sharp, consider the most vulnerable group k∗ ∈ arg min k

nk − hk + 1 . nk

40

In this group, replace the first coordinate of exactly nk∗ − hk∗ + 1 rows by a value c, and let c → ∞, leaving all other entries unchanged. Let I be the set of contaminated rows. Since |I| = nk∗ − hk∗ + 1, any feasible subset of size hk∗ must contain at least one row from I. At any exact minimizer, the µ-block satisfies the weighted mean equation nk

(k∗ )

b µ

∗ 1 X (k ) ⋆(k ) w bi ∗ z i ∗ . = hk∗ i=1

Hence, for the first coordinate, n o  c 1 (k ) (k ) µ b1 ∗ = # i:w bi ∗ = 1 ∩ I + hk∗ hk∗

X

⋆(k )

zi1 ∗ .

(k )

{i:w bi ∗ =1}∩I c

The first term on the right-hand side diverges to +∞, while the second term is bounded because it only involves clean data. Therefore, at contamination fraction (nk∗ −hk∗ +1)/nk∗ = β0 , either the location estimator diverges, or the scatter estimator has already broken down. This proves that the common lower bound β0 is sharp. Appendix C. Simulation studies We conducted a comprehensive set of simulation studies to evaluate the proposed joint estimation algorithm, and compared its performance with the non-joint MCD [13, 16] and cellMCD [14] methods. For each configuration, the simulation was repeated 100 times, and the reported results are the means and standard deviations of the evaluation metrics over these 100 replications. For the joint method, in each simulation replication we generated K groups of q-dimensional observations, where the k-th group contained nk samples. Clean observations were generated from a multivariate normal distribution Nq (µ(k) , Σ(k) ). The group mean vectors were set as µ(k) = (k − 1)∆ · 1q , where ∆ controls the mean shift across groups. The covariance matrices followed Σ(k) = Θ ⊙ Γ(k) , where Θ denotes the correlation structure shared by all groups, Γ(k) is a group-specific positive definite matrix for the k-th group, and ⊙ denotes the Hadamard product. The baseline correlation structures for Θ and Γ(k) were generated by two mechanisms: one was an AR(1)-type correlation structure ρ|i−j| with ρ = 0.9, 1 ≤ i, j ≤ q, referred to as setting A09; the other was a random correlation matrix with a prespecified condition number of 100, referred to as setting ALYZ [48]. Additional coordinate-wise scale perturbations were then introduced into Γ(k) to increase heterogeneity across groups. We considered two contamination schemes, namely casewise contamination and cellwise contamination. In the casewise setting, for each group we randomly selected a proportion εout of entire observations and replaced them with outliers located along the eigenvector corresponding to the smallest eigenvalue of Σ(k) , with a fixed Mahalanobis magnitude relative to µ(k) . In the cellwise setting, for each variable we independently contaminated a proportion εout of cells, and imposed anomalous shifts along the smallest-eigenvalue direction within the corresponding subspace. During estimation, the trimming subset size for the k-th group was set to hk = max{⌊(1 − εout )nk ⌋, q + 1}. Specifically, we set K = 5, q ∈ {70, 20}, nk ∈ {83, 100, 500}, ∆ = 2, and εout ∈ {0.25, 0.5}. The detailed combinations are shown in Table C.5. For each configuration, we summarized across independent replications the mean squared errors (MSE) of the mean vectors and the relative errors (RE) of covariance matrices, and also reported the Kullback-Leibler (KL) discrepancy for covariance estimation. For non-joint competitors, separate estimations were performed using covMCD from robustbase [49, 50] for casewise contamination, and cellMCD from 41

cellWise [14] for cellwise contamination. Table C.5 shows the proposed joint framework’s advantage heavily depends on the contamination level and mode. Table C.5: Simulation results under different settings. “-” indicates that the metric exceeds 106 . The values in parentheses represent the standard deviation of the 100 replications. ↓ means the smaller the better. Mode

casewise

cellwise

k

nk

εout

q

1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5

500 0.500 500 0.500 500 0.500 500 0.500 500 0.500 500 0.250 500 0.250 500 0.250 500 0.250 100 0.250 100 0.250 100 0.250 100 0.250 100 0.250 100 0.250 500 0.500 500 0.500 500 0.500 500 0.500 100 0.500 500 0.250 500 0.250 500 0.250 500 0.250 100 0.250 100 0.250 100 0.250 100 0.250 100 0.250 83 0.250

70 70 70 70 70 70 70 70 70 70 20 20 20 20 20 70 70 70 70 70 70 70 70 70 70 20 20 20 20 20

MSE(µ) ↓ Joint Non-joint 0.067 (0.087) 0.004 (0.001) 0.067 (0.114) 0.004 (0.001) 0.076 (0.115) 0.004 (0.001) 0.058 (0.089) 0.019 (0.003) 0.073 (0.116) 0.052 (0.080) 0.004 (0.001) 0.026 (0.009) 0.004 (0.001) 0.031 (0.010) 0.004 (0.001) 0.067 (0.087) 0.004 (0.001) 0.066 (0.090) 0.018 (0.003) 0.063 (0.084) 0.019 (0.007) 0.110 (0.126) 0.018 (0.006) 0.102 (0.134) 0.018 (0.006) 0.120 (0.134) 0.018 (0.005) 0.106 (0.109) 0.019 (0.006) 0.110 (0.129) 0.004 (0.001) 0.052 (0.067) 0.004 (0.001) 0.051 (0.086) 0.004 (0.001) 0.058 (0.086) 0.004 (0.001) 0.045 (0.068) 0.027 (0.012) 0.057 (0.089) 0.003 (0.001) 0.059 (0.069) 0.003 (0.001) 0.063 (0.106) 0.003 (0.001) 0.083 (0.079) 0.003 (0.001) 0.007 (0.006) 0.015 (0.003) 0.004 (0.001) 0.027 (0.009) 0.114 (0.124) 0.027 (0.007) 0.028 (0.009) 0.026 (0.009) 0.027 (0.009) 0.025 (0.009) 0.027 (0.009) 0.031 (0.010) 0.028 (0.008)

KL(Σ) ↓ Joint Non-joint 2053.184 (7.963) 32.368 (20.637) 2053.486 (10.050) 37.130 (22.755) 2052.545 (9.889) 34.097 (18.656) 2054.223 (7.897) 64787.089 (6307.088) 2052.555 (10.530) 6.933 (0.199) 78.172 (133.667) 6.936 (0.210) 575.151 (879.586) 6.947 (0.214) 6.886 (0.195) 576.000 (136.053) 58.688 (1.642) 600.097 (164.398) 3.160 (0.291) 123.034 (88.941) 3.150 (0.297) 135.536 (119.664) 3.114 (0.310) 118.094 (86.670) 3.074 (0.298) 139.030 (89.511) 3.123 (0.331) 121.324 (75.300) 35.450 (9.590) 35.181 (8.702) 35.135 (8.327) 35.586 (9.128) 102.690 (26.680) 14.332 (0.808) 576.208 (157.316) 14.441 (0.764) 583.249 (182.543) 14.434 (0.767) 14.367 (0.797) 12984.516 (26078.985) 55.587 (2.579) 31.899 (18.133) 7.399 (1.149) 101.185 (58.588) 7.380 (1.197) 178.290 (457.589) 7.428 (1.153) 66.127 (114.598) 7.369 (0.993) 93.969 (180.187) 8.967 (1.378) 78.029 (137.850)

RE(Σ) ↓ Joint Non-joint 112.888 (6.122) 0.567 (0.020) 113.360 (5.760) 0.569 (0.018) 113.080 (5.571) 0.569 (0.021) 112.871 (6.092) 1.006 (0.040) 112.747 (6.004) 129.392 (6.887) 0.412 (0.008) 0.916 (0.085) 0.413 (0.007) 0.991 (0.109) 0.413 (0.007) 75.739 (17.417) 0.411 (0.007) 67.064 (3.762) 0.936 (0.023) 67.766 (3.720) 0.488 (0.030) 17.695 (2.431) 0.488 (0.034) 17.592 (2.468) 0.486 (0.033) 17.911 (2.318) 0.483 (0.032) 17.442 (2.560) 0.485 (0.033) 17.872 (2.467) 0.884 (0.203) 129.312 (7.089) 0.863 (0.181) 129.760 (6.841) 0.881 (0.187) 129.531 (6.701) 0.879 (0.197) 129.240 (6.917) 2.371 (1.322) 129.114 (7.002) 0.474 (0.010) 67.310 (3.780) 0.474 (0.010) 67.110 (4.034) 0.475 (0.009) 109.445 (5.974) 0.473 (0.010) 0.656 (0.177) 0.733 (0.016) 0.570 (0.020) 0.693 (0.070) 17.660 (2.353) 0.690 (0.071) 0.935 (0.096) 0.681 (0.065) 0.923 (0.084) 0.689 (0.064) 0.922 (0.084) 0.720 (0.072) 0.920 (0.096)

Under a high outlier ratio (εout = 0.500) with casewise contamination, the joint method trades some estimation accuracy for enhanced stability, preventing the severe numerical breakdown often observed in the non-joint approach. Conversely, under a moderate outlier ratio (εout = 0.250), the benefits of the joint method are evident: it consistently achieves lower estimation errors across different dimensions in casewise scenarios and improves covariance estimation stability under cellwise contamination.

42

References [1] R. Bommasani, D. A. Hudson, E. Adeli, R. Altman, S. Arora, S. von Arx, M. S. Bernstein, J. Bohg, A. Bosselut, E. Brunskill, et al., On the opportunities and risks of foundation models, 2021. ArXiv:2108.07258. [2] S. Bubeck, V. Chandrasekaran, R. Eldan, J. Gehrke, E. Horvitz, E. Kamar, P. Lee, Y. T. Lee, Y. Li, S. Lundberg, et al., Sparks of artificial general intelligence: Early experiments with GPT-4, 2023. ArXiv:2303.12712. [3] I. J. Goodfellow, J. Shlens, C. Szegedy, Explaining and harnessing adversarial examples, in: International Conference on Learning Representations, 2015. [4] N. Carlini, D. A. Wagner, Towards evaluating the robustness of neural networks, in: 2017 IEEE Symposium on Security and Privacy, IEEE, 2017, pp. 39–57. [5] J. Kirchenbauer, J. Geiping, Y. Wen, J. Katz, I. Miers, T. Goldstein, A watermark for large language models, in: Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, PMLR, 2023, pp. 17061–17084. [6] P. Manakul, A. Liusie, M. J. F. Gales, SelfCheckGPT: Zero-resource black-box hallucination detection for generative large language models, in: Proceedings of the 2023 Conference on Empirical Methods in Natural Language Processing, 2023, pp. 9004–9017. [7] S. Gehrmann, H. Strobelt, A. M. Rush, GLTR: Statistical detection and visualization of generated text, in: Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics: System Demonstrations, 2019, pp. 111–116. [8] E. Mitchell, Y. Lee, A. Khazatsky, C. D. Manning, C. Finn, DetectGPT: Zero-shot machine-generated text detection using probability curvature, in: Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, PMLR, 2023, pp. 24950–24962. [9] Z. Ji, N. Lee, R. Frieske, T. Yu, D. Su, Y. Xu, E. Ishii, Y. Bang, A. Madotto, P. Fung, Survey of hallucination in natural language generation, ACM Computing Surveys 55 (2023) 1–38. [10] J. H. Metzen, T. Genewein, V. Fischer, B. Bischoff, On detecting adversarial perturbations, in: International Conference on Learning Representations, 2017. [11] W. Tian, L. Zhang, H. Cui, Abnormal sample detection based on robust Mahalanobis distance estimation in adversarial machine learning, Statistics and Its Interface 17 (2024) 91–106. [12] I. Solaiman, M. Brundage, J. Clark, A. Askell, A. Herbert-Voss, J. Wu, A. Radford, G. Krueger, J. W. Kim, S. Kreps, et al., Release strategies and the social impacts of language models, 2019. ArXiv:1908.09203. [13] M. Hubert, M. Debruyne, Minimum covariance determinant, Wiley Interdisciplinary Reviews: Computational Statistics 2 (2010) 36–43. [14] J. Raymaekers, P. J. Rousseeuw, The cellwise minimum covariance determinant estimator, Journal of the American Statistical Association 119 (2024) 2610–2621. 43

[15] K. Lee, K. Lee, H. Lee, J. Shin, A simple unified framework for detecting out-of-distribution samples and adversarial attacks, in: Advances in Neural Information Processing Systems, volume 31, 2018. [16] P. J. Rousseeuw, K. Van Driessen, A fast algorithm for the minimum covariance determinant estimator, Technometrics 41 (1999) 212–223. [17] J. Guo, E. Levina, G. Michailidis, J. Zhu, Joint estimation of multiple graphical models, Biometrika 98 (2011) 1–15. [18] A. P. Dempster, N. M. Laird, D. B. Rubin, Maximum likelihood from incomplete data via the EM algorithm, Journal of the Royal Statistical Society: Series B (Methodological) 39 (1977) 1–22. [19] B. Guo, X. Zhang, Z. Wang, M. Jiang, J. Nie, Y. Ding, J. Yue, Y. Wu, How close is ChatGPT to human experts? comparison corpus, evaluation, and detection, 2023. ArXiv:2301.07597. [20] L. Ouyang, J. Wu, X. Jiang, D. Almeida, C. L. Wainwright, P. Mishkin, C. Zhang, S. Agarwal, K. Slama, A. Ray, J. Schulman, J. Hilton, F. Kelton, L. Miller, M. Simens, A. Askell, P. Welinder, P. Christiano, J. Leike, R. Lowe, Training language models to follow instructions with human feedback, in: Advances in Neural Information Processing Systems, volume 35, 2022, pp. 27730–27744. [21] OpenAI, Introducing ChatGPT, OpenAI, 2022. Accessed April 27, 2026. [22] OpenAI, GPT-5.4 Thinking system card, OpenAI, 2026. Accessed April 27, 2026. [23] OpenAI, Introducing GPT-5.4, OpenAI, 2026. Accessed April 27, 2026. [24] OpenAI, Introducing GPT-5.4 mini and nano, OpenAI, 2026. Accessed April 27, 2026. [25] Gemini Team, Gemini 1.5: Unlocking multimodal understanding across millions of tokens of context, 2024. ArXiv:2403.05530. [26] Google, Gemini 3 Flash preview, Google AI for Developers, 2026. Accessed April 27, 2026. [27] Google, Gemini 3.1 Flash-Lite preview, Google AI for Developers, 2026. Accessed April 27, 2026. [28] Google, Gemini 3.1 Pro preview, Google Cloud Vertex AI documentation, 2026. Accessed April 27, 2026. [29] Y. Liu, M. Ott, N. Goyal, J. Du, M. Joshi, D. Chen, O. Levy, M. Lewis, L. Zettlemoyer, V. Stoyanov, RoBERTa: A robustly optimized BERT pretraining approach, 2019. ArXiv:1907.11692. [30] D. Ippolito, D. Duckworth, C. Callison-Burch, D. Eck, Automatic detection of generated text is easiest when humans are fooled, in: Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics, 2020, pp. 1808–1822. [31] X. Li, F. Ruan, H. Wang, Q. Long, W. J. Su, A statistical framework of watermarks for large language models: Pivot, detection efficiency and optimal rules, Annals of Statistics 53 (2025) 322–351.

44

[32] S. Zhang, S. Roller, N. Goyal, M. Artetxe, M. Chen, S. Chen, C. Dewan, M. Diab, X. Li, X. V. Lin, T. Mihaylov, M. Ott, S. Shleifer, K. Shuster, D. Simig, P. S. Koura, A. Sridhar, T. Wang, L. Zettlemoyer, OPT: Open pre-trained transformer language models, 2022. ArXiv:2205.01068. [33] J. Li, X. Cheng, X. Zhao, J.-Y. Nie, J.-R. Wen, HaluEval: A large-scale hallucination evaluation benchmark for large language models, in: Proceedings of the 2023 Conference on Empirical Methods in Natural Language Processing, Association for Computational Linguistics, Singapore, 2023, pp. 6449–6464. [34] X. Ma, B. Li, Y. Wang, S. M. Erfani, S. Wijewickrema, G. Schoenebeck, M. E. Houle, D. Song, J. Bailey, Characterizing adversarial subspaces using local intrinsic dimensionality, in: International Conference on Learning Representations, 2018. [35] A. Krizhevsky, G. Hinton, Learning multiple layers of features from tiny images, Technical Report, University of Toronto, Toronto, 2009. [36] A. Kurakin, I. J. Goodfellow, S. Bengio, Adversarial examples in the physical world, in: Artificial Intelligence Safety and Security, Chapman and Hall/CRC, 2018, pp. 99–112. [37] S.-M. Moosavi-Dezfooli, A. Fawzi, P. Frossard, DeepFool: A simple and accurate method to fool deep neural networks, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2016, pp. 2574–2582. [38] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, L. Fei-Fei, ImageNet: A large-scale hierarchical image database, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2009, pp. 248–255. [39] H. Xue, A. Araujo, B. Hu, Y. Chen, Diffusion-based adversarial sample generation for improved stealthiness and controllability, in: Advances in Neural Information Processing Systems, volume 36, 2023, pp. 2894–2921. [40] K. He, X. Zhang, S. Ren, J. Sun, Deep residual learning for image recognition, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2016, pp. 770–778. [41] G. Huang, Z. Liu, L. van der Maaten, K. Q. Weinberger, Densely connected convolutional networks, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2017, pp. 4700–4708. [42] B. Iglewicz, D. C. Hoaglin, How to Detect and Handle Outliers, ASQC Quality Press, Milwaukee, WI, 1993. [43] A. Beck, M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM Journal on Imaging Sciences 2 (2009) 183–202. [44] Y. Nesterov, Introductory Lectures on Convex Optimization: A Basic Course, 1st ed., Springer, New York, 2004. [45] A. Beck, First-Order Methods in Optimization, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2017. [46] R. G. Bartle, D. R. Sherbert, Introduction to Real Analysis, 4th ed., Wiley, New York, 2018. 45

[47] R. A. Horn, C. R. Johnson, Matrix Analysis, 2nd ed., Cambridge University Press, Cambridge, 2012. [48] C. Agostinelli, A. Leung, V. J. Yohai, R. H. Zamar, Robust estimation of multivariate location and scatter in the presence of cellwise and casewise contamination, TEST 24 (2015) 441–461. [49] M. Maechler, P. Rousseeuw, C. Croux, V. Todorov, A. Ruckstuhl, M. Salibian-Barrera, T. Verbeke, M. Koller, E. L. T. Conceicao, M. A. di Palma, robustbase: Basic Robust Statistics, 2026. R package version 0.99-7. [50] V. Todorov, P. Filzmoser, An object-oriented framework for robust multivariate analysis, Journal of Statistical Software 32 (2009) 1–47.

46