Robust covariance estimation with missing values and cell-wise contamination
Authors: gregoire pacreau
Karim Lounici
Venue: NeurIPS 2023 poster (2023)
Keywords: robust statistics
missing values
cell-wise contamination

Abstract: Large datasets are often affected by cell-wise outliers in the form of missing or erroneous data. However, discarding any samples containing outliers may result in a dataset that is too small to accurately estimate the covariance matrix. Moreover, the robust procedures designed to address this problem require the invertibility of the covariance operator and thus are not effective on high-dimensional data. In this paper, we propose an unbiased estimator for the covariance in the presence of missing values that does not require any imputation step and still achieves near minimax statistical accuracy with the operator norm. We also advocate for its use in combination with cell-wise outlier detection methods to tackle cell-wise contamination in a high-dimensional and low-rank setting, where state-of-the-art methods may suffer from numerical instability and long computation times. To complement our theoretical findings, we conducted an experimental study which demonstrates the superiority of our approach over the state of the art both in low and high dimension settings.

Robust covariance estimation with missing values and cell-wise contamination
Karim Lounici
Gregoire Pacreau CMAP Ecole Polytechnique
CMAP
Ecole Polytechnique
Palaiseau, France
Palaiseau, France
karim.oumiicc@polyechniiqqeeee
gregoire.paceauupppyycchhhgeeeee
Abstract
Large datasets are often affected by cell-wise outliers in the form of missing erroneous data. However, discarding any samples containing outliers may result dataset that is too small accurately estimate the covariance matrix. Moreover, the robust procedures designed to address this problem require the invertibility nigh-dimensional data. In
the covariance operator and thus are not effective this paper, we propose an unbiased estimator for the covariance in the presence of missing values that does not require any imputation step and still achieves near minimax statistical accuracy with the operator norm. We also advocate for its in combination with cell-wise outlier detection methods tackle cell-wise contamination in high-dimensional and low-rank setting, where state-of-the-ar methods may suffer from numerical instability and long computation times. To complement our theoretical findings, we conducted an experimental study which demonstrates the superiority of our approach over the state of the art both in low and high dimension settings.
Introduction
Outliers are a common occurrence in datasets, and they can significantly affect the accuracy of data analysis. While research on outlier detection and treatment has been ongoing since the 1960s, much of it has focused on cases where entire samples are outliers (Huber's contamination model) [ 41410 13 While sample- wise contamination is common issue in many datasets, modern data analysis often involves combining data from multiple sources. For example, data may be collected from array of sensors, each with independent probability of failure, financial data may come from multiple companies, where reporting errors from one source do not necessarily impact the validity of the information from the other sources. Discarding an entire sample as outlier when only a few features are contaminated can result in the loss of valuable information, especially in high-dimensional datasets where samples are already scarce. It is important identify and address the specific contaminated features, rather than simply treating the entire sample as an outlier. In fact, if each dimension of sample has a contamination probability of e, then the probability that sample containing at least one outlier is given by 1 (1 - E)², where is the dimensionality of the sample. high dimension, this probability can quickly exceed 50%, surpassing the breakdown point of many robust estimators designed for the Huber sample-we wise contamination setting. Hence, it is crucial to develop robust methods that can handle cell-wise contaminations and still provide accurate results.
The issue of cell-wise contamination, where individual cells in a dataset may be contaminated, was first introduced in [3] However, the issue of missing data due to outliers was studied much earlier, dating back the work of 351 Although missing values in dataset are much easier detect than outliers, they can lead to errors in estimating the location and scale of the underlying distribution 25
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

172
1558
Opcrstrr crue occerttaa cilaniee
Clastedd
Computation time
Matrix
Estimator
P=50
p=100
inversion
72
tailMV DDCMV
10-3 + 10-4
10-2+ 10
n0
0.6 10-a
0.7 0.007
no
DI
6 0.5
74 1
yes
22-
TSGS
2.2
200 + 10
yes
0.30
1.5
0.11
0.15
1.30
6.25
0.30
Datwett confanttaccan race (1
Figure Left: Estimation error of the covariance matrix for n 100. p p= 500, r(Σ) = under a Dirac contamination (tailMV and DDCMV are our methods). Here = = 1 and 5 varies in (0, 1). Right: For each method, mean computation time (in seconds) over 20 repetitions and whether it uses matrix inversion. For 100, we had to raise r (E) to 10 otherwise both DT and TSGS would fail due to numerical instability.
and can negatively affect the performance of supervised learning algorithms [7]. This motivated the development of the field of data imputation. Several robust estimation methods have been proposed to handle missing data, including Expectation Maximization (EM)-based algorithms 7 maximum likelihood estimation [14] and Multiple Imputation (25). among which we can find k-nearest neighbor imputation
42 Recently, sophisticated solutions based deep (2) or Diffusion schemes [48] have been proposed to perform complex tasks like artificial data generation or image inpainting. The aforementioned references focus solely on minimising the entrywise error for imputed entries. Noticeably, our practical findings reveal that applying state-of-the-art imputation methods to complete the dataset, followed by covariance estimation on the completed dataset, does not yield satisfactory results when evaluating the covariance estimation error using the operator norm.
39 and iterative imputation
learning, GANs 47
28
9 VAE
In comparison to data missingness or its sample-wing wise counterpart, the cell-wise contamination problem is less studied. The Detection Imputation (DI) algorithm of 32 is an EM type procedure combining a robust covariance estimation method with an outlier detection method to iteratively update the covariance estimation. Other methods include adapting methodology created for Huber contamination for the cell-wise problem, such as in 6 or [2]. In high dimensional statisties, however, most of these methods fail due to high computation time and numerical instability. Or they are simply not designed to work in this regime since they are based on the Mahalanobis distance, which requires an inversion of the estimated covariance matrix. This a major issue since classical covariance matrix estimators have many eigenvalues close to zero or even exactly equal to zero in high-dimension To the best of our knowledge, no theoretical result exists concerning the statistical accuracy of these methods in the cell-wise contamination setting contrarily to the extensive literature on Huber's contamination
Contributions. In this paper we address the problem of high- dimensional covariance estimation in the presence of missing observations and cell-wise contamination. To formalize this problem, we adopt and generalize the setting introduced in [T. We propose and investigate two different strategies, the first based on filtering outliers and debiasing and the second based filtering outliers followed by imputation and standard covariance estimation. We propose novel computationally efficient and numerically stable procedures that avoid matrix inversion, making them well-suited for high-dimensional data. We derive non-asymptotic estimation bounds of the covariance with the operator norm and minimax lower bounds, which clarify the impact of the missing value rate and outlier contamination rate. Our theoretical results also improve over [26] in the MCAR and no contamination. Next, we conduct an experimental study on synthetic data, comparing our proposed methods to the state-of-the--e (SOTA) methods. Our results demonstrate that SOTA methods fail the high-dimensional regime due to matrix inversions, while our proposed methods perform well in this regime, highlighting their effectiveness. Then we demonstrate the practical utility of our approach by applying it to real-life datasets, which highlights that the use of existing estimation methods significantly alters the spectral properties of the estimated covariance matrices. This implies

---

that cell-wise contamination can significantly impact the results of dimension reduction techniques like PCA by completely altering the computed principal directions. Our experiments demonstrate that our methods are more robust to cell-wise contamination than SOTA methods and produce reliable estimates of the covariance.
2 Missing values and cell-wise contamination setting
Let ,I...., x be n i.dd. copies of a zero mean random vector X admitting unknown covariance operator Σ = E [ @ X], where @ is the outer product. Denote by x; the jth component of vector X, for any j E [p]. All our results are non-asymptotic and cover a wide range of configurations for n and p including the high- nigh-dimensional setting p >> n. In this paper, consider the following two realistic scenarios where the measurements are potentially corrupted.
Missing values. We assume that each component X(") observed independently from the others with probability ò E (0, 1). Formally, we observe the random vector Y E IR? defined as follows:
(1) = x sisn,nn s
where dyy are independent realisations of a bernoulli random variable of parameter 8. This corresponds the Missing Completely Random (MCAR) setting of [35 Our theory also covers the more general Missing at Random (MAR) setting in Theorem
Cell-wise contamination. Here we assume that some missing components x. can be replaced with probability& by some independent noise variables, representing either poisoning of the data random mistakes in measurements. The observation vector Y then satisfies:
y - = + +(1 (1- -1-di,j
1 iisiin Sspp
(2)
where < are independent erroneous measurements and ei,j are i.i.d. bernoulli random variables
with parameter e. We also assume that the variables X,, c( di,j- ei are mutually independent.
In this scenario, component x(6) is either perfectly observed with probability 8, replaced by random noise with probability e'= 8) or missing with probability (1 - 5) - =). Cell-wise contamination as introduced in [ corresponds to the case where 1, and thus g' = 8.
In both these settings, the task of estimating the mean of the random vectors X; well-understood, as it reduces to the classical Huber setting for component- simponent-wise mean estimation. One could for instance apply the Tuker median on each component separately [ However, the problem becomes more complex when we consider non-linear functions of the data, such the covariance operator. Robust covariance estimators originally designed for the Huber setting may not be suitable when applied in the presence of missing values cell-wise contaminations.
We study simple estimator based on a correction of the classical covariance estimator on Yı...., Y, as introduced in [26] for the missing values scenario. The procedure is based on the following observation, linking Σ the covariance of the data with missing values and Σ the true covariance:
Σ= -8 8-2) diag(E*) 8-22
Note that this formula assumes the knowledge of ö. In the missing values scenario, can be efficiently estimated by a simple count of the values exactly set to or equal to NaN (not a number). In the contamination setting (2) the operator n
E (Y & Y) satisfies, for A E [E & E]:
r 82 (6 82 )d)ag(E) + e(1 - 6)A.
In this setting, as one does not know the exact location and number of outliers we propose to estimate 8 by the proportion of data remaining after the application of filtering procedure.
We denote by @ the Hadamard (or term by term) product of two matrices and by @ the outer product vectors, i.e. Var y E Rdd x @ y = xy" T. We denote by ||-|| and ||-I| the operator and Frobenius norms of matrix respectively. We denote by I1-II2 the vector normm.
Notations.
3

---

3 Estimation of covariance matrices with missing values
We consider the scenario outlined in ① where the matrix E is of approximately low rank. To quantify this, we use the concept of effective rank, which provides a useful measure of the inherent complexity of matrix. Specifically, the effective rank of 2 is defined as follows
rE)) = E I |ل
tr (E) |231
(4)
We note that 0 S r(z) S ra)k(E)) Furthermore, for approximately low rank matrices with rapidly decaying eigenvalues, we have r(E) < rank(E)) This section presents a novel analysis of the estimator defined in equation (3. which yields non-asymplotic minimax optimal estimation bound in the operator norm. Our findings represent substantial enhancement over the suboptimal guarantees reported in 26 [8] Similar results could be established for the Frobenius norm using more straightforward arguments, as those in 4 or [1I We give priority to the operator norm since aligns naturally with learning tasks such PCA. See To 211 22]
and the references cited therein.
We need the notion of Orlicz norms. For any a 2 1, the da -norms of real-valued random variable V are defined as: |Vil,, inf (u > 0, E Eexp (IV) ///l///") /u²) I 2}. A random vector X E TRP is
sub-Gaussian if and only if Va E RP, xx)x))l S II(x,
Minimax lower-bound. We now provide a minimax lower bound for the covariance estimation with missing values problem. Let Sp the set of p x p symmetric semi-positive matrices. Then, define 7} the set of matrices of Sp with effective rank at most F.
== E Sp r(S) S
Theorem 1. Let p. n, be strictly positive integers such that p 2 max max 27, Let X1...., X be i.idd. random vectors in 12 with covariance matrix D E Cy. Let (dr,)
sequence of Bernoulli random variables with probability of success (lld)iiiinnss1 be an i.i.d. (0, 1]. independent from
the XI,...,,, Xn. We observe n i.i.d. vectors Yi...., Yn E Iº such that y- = d, [n]. j E [p]. Then there exists two absolute constants C and (0, 1) such that:
(15- 2 C 2 115-11200
(E)
max Py: Σ XEC,
23
(5)
where infs, represents the infimumm over all estimators E. of matrix E. based on Yi,..., Y.-
Sketch of proof. We first build a sufficiently large test set of hard-to-leamy covariance operators exploiting entropy properties of the Grassmann manifold such that the distance between any two
Next, in order to control the Kullback- Leibler divergence of the observations with missing values, we exploit in particular interlacing properties of the eigenvalues of the perturbed covariance operators 37
distinct covariance operator is at least of the order | (())
This lower bound result improves upon 26 Theorem 2] as it relaxes the hypotheses on n and F. More specifically, the lower bound in [2 requires > 2n 882 while only need the mild assumption p 2 max(n, 2r, 2r Our proof leverages the properties of the Grassmann manifold, which has been previously utilized in different settings such sparse PCA without missing values or contamination 45 and low-rank covariance estimation without missing values or contamination 23 However, tackling missing values in the Grassmann approach adds technical challenge these proofs as they modify the distribution of observations. Our proof requires several additional nontrivial arguments to control the distribution divergences, which is a crucial step in deriving the minimax lower bound.
Non-asympiotic upper-bound the operator norm. We provide an upper bound of the es- timation error in operator norm. We write Y; = de @ X;. Let Er n -" Y, @ Y, be the classical covariance estimator of the covariance of Y. When the dataset contains missing values and corruptions, ς is biased estimator of Exploiting Equation unbiased estimator of the covariance matrix Σ:
26 proposed the following
- 8 S (8-1-8 -82))daag {4
The following result is from I8 Theorem 4.2).

---

Lemma 1. Let Xi...., .n be i.i.d. sub-Gaussian random variables in R', with covariance matrix E, and let dis, i E nI. E [l. be i.i.d bernoulli random variables with probability of success 5 > O. Then there exists an absolute constant C such that, for 1 >> O, with probability at least 1 - -
Σ =|| | SCIEl (V (12))Logr(2) 8't2
r (t logr(2)) ((ggr+))))))))
V
log(n)
82n
82nn
This result uses recent unbounded version of the non-commutative Bernstein incquality, thus yielding some improvement upon the previous best known bound of [26]. Theorem [ and Lemmal ! provide some important insights on the minimax rate of estimation in the missing values setting. In
the high-dimensional regime p 2 max(n, 2 and n Z 8 (Z)) (log r(E)) log" n, we observe that
the two bounds coincide up to a logarithmic factor in r (2), hence clarifying the impact of missing data on the estimation rate via the parameter 8.
Heterogeneous missingness. We can extend the correction to the more general case where each feature has different missing value rate known the Missing Random (MAR) setting in B5 We denote by 5, E (0, 1] the probability to observe feature X 1 S j S p and we set 5 := (5, (o))) As in the MCAR setting, the probabilities (6))) can be readily estimated by tallying the number of missing entries for each feature. Hence they will be assumed to be known for the sake of brevity. Let 6mv "07"))s()) be the vector containing the inverse of the observing probabilities and Ain 8m Sw In this case, the corrected estimator becomes:
E = Ainy ④ ΣY + (diag (dim) - Ainv) C diag
Σ
Let 5 = max, and ē = min {o,} be the largest and smallest probabilities to observe a feature. Theorem 2. (i) Let X1,..., X, be i.i.d. sub-Gaussian random variables in RP, with covariance
matrix E. We consider the MAR setting described above. Then the estimator 8 satisfies, for any
t O, with probability at least -e
8 | - - S SCIEII
(V (r(E) log
VLV ((E)((+ogrr()) + (o)r(E)) log
(9)
8²
ôn
(ii) Let P, F be strictly positive integers such that 2 max(n, (n, 27 Let XI,..., Xn be i.dd random vectors Rp with covariance matrix E E Cr. Then.
(|S- 2 Σ
Σ.)
inf max Px I DEC,
2B. B.
(10)
8
I x 5 then the rates for the MCAR and MAR settings match. The proof is a straightforward adaptation of the proof in the MCAR setting.
4 Optimal estimation of covariance matrices with cell-wise contamination
In this section, we consider the cell-wise contamination setting 2 We derive both an upper bound on the operator norm error of the estimator Q and a minimax lower bound for this specific setting. Let us assume that the ..,...., are sub-Gaussian r.v. Note also that == E[g @ Ei] is diagonal in the cell-wise contamination setting (2)
Minimax lower-bound. The lower bound for missing values still applies to the contaminated case as missing values are particular case of cell-wise contamination. But we want more general lower bound that also covers the case of adversarial contaminations.
Theorem 3. Let p. n, F be strictly positive integers such that p 2 maxfnn 2 Let XI,..., ... Xn, be i.idd random vectors in Rd with covariance matrix E E C,. Let of bermoulli random variables of probability of success 8 E (0, 11, independent to the Xiy..., X
(isis, be i.i.d. sequence
We observe i...d. vectors Yı,..., Yn. RR satisfying (2) where
i.i.d. of arbitrary distribution Q. Then there exists two absolute constants C >> 0 and 3 E (0, 1) such that:
||--||
|||||
r(E)
c(1 8)
inf max maxP S EECr Q
2
N 3
(11)
8
8
5

---

where inf represents the infimum over all estimators of matrix Σ and maxq is the maximum over all contamination Q.
The proof of this theorem adapts an argument developed to derive minimax lower bounds in the Huber contamination setting. See App. G.3 for the full proof.
Non-asymptotic upper-bound in the operator norm. Note that the term c(1 6)J A the cell-wise contamination setting is negligible when ò R 1 or O. Using the DDC detection procedure of [32), we can detect the contaminations and make smaller without decreasing å too much. For simplicity, we assume from now on that the {! are i.i.d. with common variance of. Hence A = T We
further assume that the c!u) are sub-Gaussian since we observed in our experiments that filtering removed all the large-valued contaminated cells and only a few inconspicuous contaminated cells remained. Our procedure O satisfies the following result.
Theorem 4. Let the assumptions of Theorem 7 be satisfied, We assume in addition that the observa-
tions I......y, Yn satisfy 2 with E [0,1) and å E (0, 1] and i.i.d. sub-Gaussian e! 's. Then, for any > 0, with probability at least 1 - et=
S - > $ |ل3||
(VMLE) log r(E)
r(E)(1 + logr(E)) log(n) + 8
e(1- 2
V
82m
82n
82n
(1 8)E
+ 82 /1 (1 -8)E)
V
+ log(p) D()
8(( o) of PV tr (2) log(n) 1 ++ log(p)
n
n
(1-6) (1-0))(p c(p- 2 2 + o where D(5.p) V
+
(1-6)
Eot (|tr(E) - 5(p 2) + || 53
See App for the proof. As emphasized in 20 the effective rank r(E) provides a measure of the statistical complexity of the covariance learning problem in the absence of any contamination. However, when cell-wise contamination is present, the statistical complexity of the problem may increase from r(Σ) r(A) p. Fortunately, if the filtering process reduces the proportion of cell-wise contamination & such that (1 d) tr()) < 8tr(E) and æ ||А|| I 5 Σ||- Then we can effectively mitigate the impact of cell-wise contamination. Indeed, deduce from Theorem g that
{> - > S |Ell (V
r(2)
((g))ogr log r())
t
r(E)(t + logr(2)))
c(1 6) +
log(n)
82n
8m
82n
+ 8 (0(1-))t((E))
t + log(p)
(vo + log())))g(n))
t + log(p)
(12)
where we considered for convenience the reasonable scenario where 8 (pp (p-2) 2) 2 tr(E) and of 2 |E|| The combination of the upper bound (12) with the lower bound Theorem provides the first insights into the impact of cell-wise contamination covariance estimation.
5 Experiments
In our experiments, refers either to the debiased MCAR covariance estimator or to MAR extension B). The synthetic data generation is described in App- A We also performed experiments on real life datasets described in App. B All experiments were conducted on a 2020 MacBook Air
with MI processor (8 cores, 3.4 GHz).

---

Table 1: Execution time of the covariance estimation procedures (in milliseconds) with n 300 averaged over all values of the contamination rate å and 20 repetitions.
method
p 50
p 100
p = 500
MV (ours)
0.29 0.03
0.49 + 0.08
9.7 + 4.5
KNNImputer (KNN)
26 + 9.8
45 + 17
470 = 190
Iterativelmputer (II)
940 + 350
2, 2,800 + 900
3.7 x 109 1.1 x 10a
Gain
6, 900 + 480
1.1 x 10 250
8.8 x 10* 1.1 103
MIWAE
5.1 x 10a + 2.8 10s
6.7 x 10a 1 550
1.77 x 10 5.8 103
Sepestics
LIVAEE LS
MIDAA
61T
GASS
Spantal 1014
02
07
Milling oilees
Mistza valao
Figure Estimation error on a synthetic dataset with 50, 300, r (E)= = 5.
Figure 3: Estimation error on synthetic
dataset with p 500, 300, r (E) = 5.
5.1 Missing Values
We compared our method to popular imputations methods: KNNImputer (KNNI), which imputes the missing values based on the k-nearest neighbours 39 and Iterat iiee IeI (II), which inspired by the R package MICE
22 as coded sklearn [30]; and two recent GANs-based imputation methods MIWAE 28 and GA IN 47 found in the package hyperimpute [15. The deep methods were tested using the same architectures, hyperparameters and early stopping rules as their respective papers.
In Figures 2 and Table we compare our estimator MV defined in (6) to these imputation methods combined with the usual covariance estimator on synthetic data (see App. A for details of data generation) terms of statistical accuracy and execution time. First, MV beats all other methods in low-dimensional scenarios and maintains competitive edge with II in high-dimensionall situations when the missing data rate remains below 30%. Furthermore, it stands as the second-bess best choice when dealing with missing data rates exceeding 35%. Next, MV has by far the smallest execution time down several orders of magnitude while the execution time of II increases very quickly with the dimension and can become impractical (see Figure y for a dataset too large for T)). Overall, the procedures MV and II perform better than MIWAE and GAIN in this experiment. Our understanding training metrics designed to minimize the entrywise error of imputation.
is that MIWA WAE and GAIN suspect this may be why their performances for the estimation of covariance with operator norm not par with other minimax methods. An interesting direction would be to investigate whether training MI
WAE and CAIN with different metrics may improve the operator norm performance.
We refer to App. for more experiments in the MAR setting of 28% Annex 3] which led to similar conclusions. These results confirm that imputation of missing values not mandatory for accurate estimation of the covariance operator. Another viable option is to apply a debiasing correction to the empirical covariance computed on the original data containing missing values. The advantage of this approach is its low computational cost even in high high-dimension.
7

---

Table 2: We consider contaminated data following model contaminated with Dirac contamination of high intensity with and for several values of á a grid. For each ò, we average the proportion of real data 8 and contaminated data after filtering over 20 repetitions. Values are displayed in percentages (ō must be high, é low)). STD stands for standard deviation.
CONTAMINATION
TAIL CUT
DDC 99%
DDC 90%
RATE
8
STD
STD
6
STD
0.1 %
99.6
0.023
0.000
0.000
99.1
0.029
0.000
0.000
94.8
0.054
0.00
0.00
1%
98.8
0.027
0.000
0.000
98.2
0.037
0.000
0.00
94.3
0.102
0.00
0.00
5%
94.9
0.013
0.000
94.6
0.000
0.000
91.8
10%
90.0
0.004
0.000
0.000
89.9
0.016
0.00
0.000
88.2
0.109
0.000
0,000
20%
80.0
0.000
20.0
0.000
80.0
0.003
0.017
0.035
79.4
0.035
0.009
0.022
30%
70.0
0.000
30.0
0.000
70.0
0.001
3.48
2.19
69.9
0.015
2.930
2.31
5.2 Cell-wise contamination
Methods tested. Our baselines are the empirical covariance estimator applied without care for contamination and an oracle which knows the position of every outlier, deletes them and then computes the MV bias correction procedure G In view of Theorems and this oracle procedure is the best possible in the setting of cell-wise contamination. Hence, we have a practical framework to assess the performance of any procedure designed to handle cell-wise contamination.
The SOTA methods the cell-wise contamination setting are the DI DeteccionnInpuutation) method 32 and the TSGS method (Two Step Generalised S-estimator) 2 Both these methods were designed to work in the standard setting n > p but cannot handle the high-dimensions included comparisons of our methods to them in the standard setting n > p. The code for DI and TSGS are from the R packages cellwise and GSE respectively.
setting as we
already mentioned. Nevertheless,
We combine the DDC detection procedure 34 to first detect and remove outliers with several estimators developed to handle missing values. Our main estimators are DDCMV (short for Detecting Deviating Cells Missing Values), which uses first DDC and then computes the debiaised covariance estimator
6 on the filtered data, and t sailll 1MV, which detects outliers through thresholding and then uses again ( But we also proposed to combine the DDC procedure with imputation methods KNNI, II, GAIN and MIWAE and finally compute the standard covariance estimator on the completed data. Hence we define four additional novel robust procedures which call DDCKNN, DDCII, DDCGAIN and DDCMIWAE. To the best of our knowledge, neither the first approach combining filtering with debiasing nor the second alternative approach combining filtering with missing values imputation have never been tested to deal with cell- wise contamination. detailed description of each method is provided in App.
Outlier detection and estimation error under cell-wise contamination synthetic data. We showed that the error of a covariance estimator under cell-wise contamination depends on the proportion of remaining outliers after filtration. In Table 2 we investigate the filtering power of the Tail Cut and DDC methods in presence of Dirac contamination. We consider the cell-wise contamination setting in the most difficult case = 1 which means that an entry is either correctly observed or replaced by an outlier (in other words, the dataset does not contain any missing value). For each values of 5 in a grid, the quantities 8 and é are the proportions of true entries and remaining contaminations after filtering averaged over 20 repetitions. The DDC based methods are particularly efficient since the proportion of Dirac contamination drops from 1 - 5 to virtually 0 for any d 2 0.74. In Figures and 4 we see that the performance of our method is virtually the same as the oracle OracleMV as long as the filtering procedure correctly eliminates the Dirac contaminations. As soon as the filtering procedure fails, the statistical accuracy brutally collapses and our DDC based estimators no longer do better than the usual empirical covariance. In Table S and Figure 5 we repeated the same experiment but with a centered Gaussian contamination. Contrarily to the Dirac contamination scenario, we see in Figure S that the statistical accuracy of our DDC based methods slowly degrades as the contamination rate increases but their performance remains significantly better than that of the usual empirical covariance.
8

---

Oyesee 400 OFICCAA
Figure 4: Estimation function of the contamination rate for n = 500, p = 400, r(E)= 5 and Dirac contamination.
Figure 5: Estimation error as a function of the contamination rate for 500, p = 400, r(E) = 5 and Gaussian contamination
1X
1400
4.1
12
4.
1:3
DIK
2000
7.6
2.6
DIC SE
47
EOO MIRAA
11
47
TOO SIWAR
5:3
DOC WW..
BOC NOW
DOC SIOUE
300 GAR
0000
Figure 6: DI fails on ATTEND since the co- variance matrix is approximately low rank. The dataset has only features and the effec- tive rank of its covariance matrix is below 2.
Figure 7: Woolridge's CEOSAL dataset fails both TSGS and DI with its dimension of 13 and effective rank of around 2.5.
5.3 The effect cell-wise contamination on real-life datasets
We tested the methods on datasets from sklearn and Woolridge's book on econometrics A6 These low dimensional datasets (less than 20 features) representing various medical, social and economic phenomena. We also included 2 high-dimensional datasets. See App. for the list of the datasets.
One interesting observation is that the instability of Mahalanobis distance-based algorithms is not limited to nigh-dimensional datasets. Even datasets with a relatively small number features can exhibit instability. This can be seen in the performance of DI on the Attend dataset, as depicted in Figure where it fails to provide accurate results. Similarly, both TSGS and D fail to perform well on the CEOSAL2 dataset, as shown in Figure despite both datasets having fewer than 15 features.
On the Abalone dataset, once have removed obvious outliers (which are detected by both DDC and the tail procedure), all estimators reached a consensus with the non-robust classical estimator, meaning that this dataset provides a ground truth against which we can evaluate and compare the performance of robust procedures in our study. To this end, we artificially contaminate 5% of the cells random in the dataset with Dirac contamination and compare the spectral error of the different robust estimators. As expected, TSGS and all our new procedures succeed at correcting the eиror, however becomes unstable (see Table 3). DDC MIWAE close SOTA TSGS for cellwise contamination and DDC II performs better. We also performed experiments on two high- dimensional datasets, where our methods return stable estimates
the covariance (DDCMV99 and DDCMV95 are within 3% of each other) and farther away from the classical estimator (See Figures g anddy y Note also that DDCTI's I'S computation time explodes and even returns out-of-memory errors duc to the high computation cost of II that we already highlighted in Table I

---

20
1.5
5-1
97
9.5
5-6
6.1
DOC MV
DOC
DOC 119
5
DOC GAN
10N MOOO
060 002
7000
DOC MINK
DOC GAN
D0C
DOC
Figure 8: Relative spectral difference (in %) between covariance estimators on SP500 stock returns 2021 and 2022. On high- dimensional data, DDCII becomes inconsis- tent with the other procedures, maybe because IT does not scale well with dimension.
Figure 9: Relative spectral difference (in %) between between covariance estimators on NASDAQ stock returns over 2021 and 2022. Here, DDCII fails due to out-of-memory er-
Table Relative spectral difference (in %) between estimated covariance matrices on Abalone with 5% synthetic contamination (6 = 0.95, E = 1). On the cleaned dataset, all the robust estimators are very close to empirical covariance (relative differences 5%), so we consider empirical covariance matrix the truth. Here the D procedure fails probably due to numerical errors.
relative error estimator
Classical
DDC MIWAF
DDC GAIN
DDCMV99
DDCMV95
DDC
DDC KNN
TSGS
DI
Truth std
12.8
4.12
6.81
1.70
2.06
3.46
5.06
3.06
8.85
0.45
0.29
0.26
0.10
0.092
0.18
0.35
0.21
1.48
Classical
13.1
14.3
13.0
13.0
12.9
13.1
13.4
14.9
DDCMV99
2.99
2.52
2.22
1.87
2.66
5.44
8.79
DDCMV95
5.27
5.03
4.04
3.71
8.28
9,99
DDC II DDC KNN
0.465
1.88
3.49
3.27
8.28
3.19
3.46
8.15
DDC MIWAE DDC GAIN TSGS
1.70
7.22
5.97
6.71
6.94
6 Conclusion and future work
In this paper, we have extended theoretical guarantees on the spectral error of our covariance estimators robust missing data to the missing random setting. We have also derived the first theoretical guarantees in the cell-wise contamination setting. We highlighted in our numerical experimental study that in the missing value setting, our debiased estimator designed tackle missing values without imputation offers statistical accuracy similar to the SOTA Iterat ve Imputed e for dramatic computational gain. We also found that SOTA algorithms in the cell-wise contamination setting often fail in the standard setting p < n for dataset with fast decreasing eigenvalues (resulting in approximately low rank covariance), a setting which is commonly encountered in many real life applications. This is due to the fact that these methods use matrix inversion which is unstable to small eigenvalues in the covariance structure and can even fail return any estimate. In contrast, we showed that our strategy combining filtering with estimation procedures designed to tackle missing values produce far more stable and reliable results. future work, we plan to improve our theoretical upper and lower bounds in the cell-wise contamination setting to fully clarify the impact of this type of contamination covariance estimation.
This paper based upon work partially supported by the Chaire Business Analytic for Future Banking and EU Project ELIAS under grant agreement No. 101 101237.
Acknowledgements.
10

---

References
[1] Pedro Abdalla and Nikita Zhivotovskiy. Covariance Estimation: Optimal Dimension-free Guarantees for Adversarial Corruption and Heavy Tails, July 2023.
[2] Claudio Agostinelli, Andy Leung. Victor J. Yohai, and Ruben H. Zamar. Robust estimation of multivariate location and scatter in the presence of cellwise and casewise contamination, June 2014.
[3] Fatemah Alqallaf, Stefan Van Aelst, Victor J. Yohai, and Ruben H. Zamar. Propagation of outliers in multivariate data. The Annals of Statistics, 37(1):311-331, February 2009.
[4] Florentina Bunea and Luo Xiao. On the sample covariance matrix estimator of reduced effective rank population matrices, with applications to fPCA. Bernoulli, 21(2), May 2015.
[5] Mengjie Chen, Chao Gao, and Zhao Ren. Robust Covariance and Scatter Matrix Estimation under Huber's Contamination Model.
[math, stat), June 2017.
[6] Mike Danilov, Victor Yohai, and Ruben Zamar. Robust Estimation of Multivariate Location and Scatter in the Presence Missing Data. JASA. Journal of the American Statistical Association, 107, September 2012.
[7] A. Dempster, N. M. Laird, and D. B. Rubin. Maximum Likelihood from Incomplete Data via the EM Algorithm. Journal of the Royal Statistical Society, Series B (Methodological), 39(1);1-38, 1977.
[8] Sheela Devadas, Peter Haine, and Keaton Stubis. The Schur-Horn Theorem. 2015.
[9] Weinan Dong, Daniel Yee Tak Fong. Jin-sun Yoon, Eric Yuk Fai Wan, Laura Elizabeth Bedford, Eric Ho Man Tang, and Cindy Lo Kuen Lam. Generative adversarial networks for imputing missing data for big data clinical research. BMC Medical Research Methodology, 21(1):78, April 2021.
[10] Alessio Farcomeni. Robust Constrained Clustering in Presence of Entry-Wise Outliers. Techno- metrics, 56, February 2014.
[11] Peter J. Huber. Robust Estimation of Location Parameter. The Annals of Mathematiial Statistics, 35(1):73-101, 1964.
[12] Peter J. Huber and Elvezio M. Ronchetti. Robust Statistics, 2nd Edition I Wiley. Wiley Series in Probability and Statistics. John Wiley and Sons, Inc., 2009.
[13] Mia Hubert, Michiel Debruyne, and Peter J. Rousseeuw. Minimum Covariance Determinant and Extensions. WIREs Computational Statistics, 10(3), May 2018
[14] Mortaza Jamshidian and Peter M. Bentler. ML Estimation of Mean and Covariance Structures with Missing Data Using Complete Data Routines. Journal of Educational and Behavioral Statistics, 24(1):21 1999.
[15] Daniel Jarrett, Bogdan Cebere, Tennison Liu, Alicia Curth, and Mihaela van der Schaar. Hyperimpute: Generalized iterative imputation with automatic model selection. 2022.
[16] Charles R. Johnson. Matrix theory and applications (charles r. johnson, ed.). Proceedings of symposia in applied mathematics, 40, 1989.
[17] Julie Josse, Nicolas Prost, Erwan Scomnet, and Gaël Varoquaux. On the consistency of supervised learning with missing values.
[cs, math, star],Jlyy 2020.
[18] Yegor Klochkov and Nikita Zhivotovskiy. Uniform Hanson-Wrighh type concentration inequali- ties for unbounded entries via the entropy method, August 2019.
[19] Vladimir Koltchinskii and Karim Lounici. Asymptotics and concentration bounds for bilin- ear forms of spectral projectors of sample covariance. Annales de l l'Institut Henri Poincaré, Probabilités et Statistiques, 52(4):1976 2013, 2016.

---

[20] Vladimir Koltchinskii and Karim Lounici. Concentration inequalities and moment bounds for sample covariance operators. Bernoulli, 23(1): 10-133, February 2017.
[21] Vladimir Koltchinskii and Karim Lounici. New asymptotic results in principal component analysis. Sankhya A, 79(2):254 297, Aug 2017.
[22] Vladimir Koltchinskii and Karim Lounici. Normal approximation and concentration of spectral projectors sample covariance. The Annals of Statistics, 45(1):121 157, 2017.
[23] Vladimir Koltchinskii, Karim Lounici, and Alexander B. Tsybakov. Estimation of Low-Rank Covariance Function, April 2015.
[24] Andy Leung, Hongyang Zhang, and Ruben H. Zamar. Robust regression estimation and inference the presence of cellwise and casewise contamination. Computational Statistics & Data Analysis, 99:1-11, July 2016.
[25] Roderick Little and Donald Rubin. Statistical Analysis with Missing Data, Second Edition. Wiley Series in Probability and Mathematical Statistics. Probability and Mathematical Statistics. Wiley edition, 2002.
[26] Karim Lounici. High-dimensional covariance matrix estimation with missing observations. Bernouli,, 20(3):122--1058, August 2014.
[27] Chao Ma. Sebastian Tschiatschek, José Miguel Hemmández-Lobato Richard Turner, and Cheng Zhang. VAEM: A Deep Generative Model for Heterogeneous Mixed Type Data, June 2020.
[28] Pierre-Alexandee Mattei and Jes Frellsen. MIWAE: Deep Generative Modelling and Imputation of Incomplete Data Sets. In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, Proceed- ings of the 36th International Conference on Machine Learning, ICML 2019, 9-15 June 2019, Long Beach, California, USA, volume 97 of Proceedings of Machine Learning Research, pages 4413 4423. PMLR, 2019.
[29] Alain Pajor. Metric Entropy of the Grassmann Manifold. Complex Geometry Analysis, 34:181- 188, 1998.
[30] Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and Duchesnay. Scikit-learn: Machine learning Python, Journal of Machine Learning Research, 12:2825-288 2830, 2011.
[31] Nikita Puchkin, Fedor Noskov, and Vladimir Spokoiny. Sharper dimension-free bounds on the Frobenius distance between sample covariance and its expectation, August 2023.
[32] Jakob Raymaekers and Peter Rousseeuw. Handling cellwise outliers by sparse regression and robust covariance. 1912: [stat], December 2020.
[33] Jakob Raymackers and Peter Rousseeuw. Fast robust correlation for high-dimensional data. Technometrics, 63(2):184-198, April 2021.
[34] Peter J. Rousseeuw and Wannes Van den Bossche. Detecting deviating data cells. Technometrics, 6045), April 2018.
[35] Donald B. Rubin. Inference and Missing Data. Biometrika, 63(3):581-592, 1976.
[36] Eckhard Schlemm. The kearns saul inequality for bernoulli and poisson-binonie distributions. Journal of theoretical probabiliry, 29:48 62, 2016.
[37] R. C. Thompson. Principal submatrices of normal and Hermitian matrices. Illinois Journal of Mathematics, 10(2):296-308, June 1966.
[38] Joel A. Tropp. An introduction to matrix concentration inequalities. Foundations and Trends in Machine Learning, 1(1-,):1-230, 2015.
12

---

[39] Olga Troyanskaya, Michael Cantor, Gavin Sherlock, Pat Brown, Trevor Hastie, Robert Tib- shirani, David Botstein, and Russ B. Altman. Missing value estimation methods for DNA microarrays. Bioinformatics, 17(6):520-5255 520 June 2001.
[40] Alexandre B. Tsybakov. Nonparametric estimators. In Alexandre B. Tsybakov, editor, Intro- duction to Nonparametric Estimation, Springer Series Statistics, pages 1-76. Springer, New York, NY, 2009.
[41] John W. Tukey. The Ninther, Technique for Low-Effort Robust (Resistant) Location in Large Samples. In H. A. David, editor, Contributions to Survey Sampling and Applied Statistics, pages 251-257. Academic Press, January 1978.
[42] Stef van Buuren and Karin Groothuis- Oudoshoo--Odssoom... Mice: Multivariate Imputation by Chained Equations in R. Journal of Statistical Software, 45:1-67, December 2011.
[43] Roman Vershynin. Introduction to the non-asymptotic analysis of random matrices, November 2011.
[44] Roman Vershynin. High- Sige--iinnsooall Probability: An Introduction with Applications in Data Science. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, Cambridge, 2018
[45] Vincent Q. Vu and Jing Lei. Minimax sparse principal subspace estimation in high dimensions. The Annals of Statistics, (6)999905-2 December 2013.
[46] Jeffrey M. Wooldridge. Introductory econometrics modern approach Jeffrey M. Wooldridge... Cengage learning, 2016.
[47] Jinsung Yoon, James Jordon, and Mihaela van der Schaar. GAIN: Missing Data Imputation using Generative Adversarial Nets, June 2018
[48] Shuhan Zheng and Nontawat Charoenphakdee. Diffusion models for missing value imputation in tabular data, March 2023.
13


Reviews:
[{"review_id": "jfRELdKxHg", "reviewer": "Reviewer_vhsW", "summary": "This paper studies the problem of estimating a covariance matrix in the presence of missing observations and/or outliers.\nThe authors first consider the missing values framework and theoretically analyse the unbiased covariance estimator introduced in [Karim Lounici. High-dimensional covariance matrix estimation with missing observations]. Then, in the framework (missing values + outliers), they also recommend its use in combination with cell-wise outlier detection methods such as DDC.\nFinally, to support their theoretical results, they carried out some experiments to demonstrate the advantage of their approach over existing methods.", "questions": "My first question is more related to Theorem 1. However, as the proofs of the other theorems use sometimes the same techniques, they are also related to them.\n\n1\\ In the proof of Theorem 1 (Appendix E.2.3), why there is no absolute value on (1/d - 1/d^2) when this scalar moves outside of the norm? This could change the proof because then |(1/d - 1/d^2)| is no longer always less than 1/d\n\n2\\ In the section \"Experiments\", it is stated that the best oracle is the one that knows \"the position of every outlier, deletes them and then computes the MV bias correction procedure\". Isn't an oracle that knows \"the position of every outlier\" but instead of deleting them, performs a good imputation and then uses a standard estimator better?\n\n3\\ It seems that only the estimator of equation (5) is used in the experiments. Why not also use the estimator from equation (10)?\n\n================\n\nOne final remark: My score reflects my concerns about the first question. Of course, I'm willing to raise the score if there's something I'm missing.", "limitations": "-", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1\\  It is a real pleasure to read this paper. Indeed, It is very well written, which makes it easy to understand\n\n2\\ The contributions are clear and well-supported. In addition, the theoretical results (if true - see next section) are theoretically sound because there are probability bounds on the difference between the estimated covariance and the true covariance, but more importantly there are minimax bounds.\n\n3\\ The experimental part is again very clear and seems to compare well to existing methods with the error bars and calculation time indicated.", "weaknesses": "1\\ The contamination is assumed to be sub-Gaussian. This contrasts with the recent robust literature, which attempts to take into account any type of contamination.\n\n2\\ The probability of missing values d and the outlier epsilon are assumed to be known in the theoretical results (although I understand that a full analysis of the estimator when d and epsilon are estimated is certainly difficult). \n\n3\\  3. There are possible problems in the theorems which should be clarified (see Questions).\n\n================\nAdditional remarks:\n\nFigure 1 is never cited and explained.\n\nI think the formal definition of the operator norm should be given in the notation part.\n\nline 101: \"The The\"\n\nline 171: \"procedure procedure\"\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Q1. We thank the referee for pointing out the issue in the proof. We actually found it after submitting the paper and found an alternative simpler argument leading to an improved bound.\n\nIndeed Ineq. (29) in our paper should be: \n$$\n\\Vert\\widehat{\\Sigma} - \\Sigma\\Vert  \\leq \\delta^{-2} \\Vert diag \\left(\\widehat{\\Sigma}^Y - \\Sigma^Y\\right)\\Vert+ \\delta^{-2} \\Vert\\widehat{\\Sigma}^Y - \\Sigma^Y\\Vert.\n $$\nNext we use the Schur-Horn Theorem to get\n$$\n\\Vert diag \\left( \\hat{\\Sigma}^Y - \\Sigma^Y\\right) \\Vert \\leq \\Vert\\hat{\\Sigma}^Y - \\Sigma^Y\\Vert.\n$$\nThen, using Ineq (23), we get\n$$\n            \\Vert \\widehat{\\Sigma} - \\Sigma\\Vert \\lesssim  \\frac{\\Vert\\Sigma\\Vert}{\\delta}\\left(\\sqrt{\\frac{{r}(\\Sigma)}{n}} \\lor \\frac{{r}(\\Sigma)}{n} \\lor \\sqrt{\\frac{t}{n}} \\lor \\frac{t}{n}\\right).\n$$\nWe got rid of the $\\log p$ in the bound making our Theorem 1 completely dimension-free and sharp minimax optimal.\n\nWe corrected the proof. We also note that the expectation sign in the LHS of (23) was a typo that we corrected. Lemma 2 is no longer needed. We hope this clarifies this issue. Thanks again! \n\nQ2. Thanks for the interesting question. According to minimax theory, given the available information, it is impossible to estimate faster than a certain universal rate, regardless of the procedure used. Hence, if our dataset has missing values, imputation will not provide more information on the problem as we rely only on the observed data to impute the missing values. Therefore we should not expect a faster minimax rate of estimation after imputation than the universal rate provided by the minimax lower bound. The minimax lower bound in Thm 2 guarantees that no estimation method can achieve a faster rate than the universal rate in Eq (7). Hence an imputation method may attain the minimax rate but cannot do better. Additionally, we proved that our debiasing method achieves the minimax rate. This means that we can build an efficient minimax optimal procedure for this covariance estimation problem without imputation.\nNow, in the cellwise contamination setting, we observe that corrupted entries lead to a deterioration in the minimax estimation rate, as shown in Eq (14), when compared to the missing values setting. Hence the cellwise contamination problem is statistically more challenging than the missing values setting. This is why we propose a user-friendly procedure that identifies and removes the outliers before applying the MV procedure.  We hope this clarifies our choice of estimator in experiments on cellwise contaminations.\n\nQ3. We considered estimator (10) only for the theoretical exploration of cellwise contamination problem. Indeed (10) requires estimating probability $\\epsilon$ and contamination covariance $\\Lambda$. This is difficult as we may not easily tell apart contamination from signal. It seemed more convenient to us in real life applications to apply a good filtering strategy which eliminates most of the contaminations (thus making the term $\\epsilon \\|\\Lambda\\|$ negligible) followed by the MV estimator (5). This strategy is user-friendly and beats the SOTA method for covariance estimation with cell-wise contaminations ([1,4,21,22,23]) in our experimental study. \n\nWe mostly agree with all 3 weaknesses pointed out by the referee but we would like to explain our choices in light of phenomena we observed in our experimental study.\n\n1. We observe in practice that heavy-tail cellwise contaminations (high-magnitude entries) are very easily filtered out of the dataset and only inconspicuous contaminations may remain (those of the same magnitude as real data or smaller). This is why we only investigated in Thm 3 the impact of sub-gaussian contaminations on the performance of a simple procedure. In the lower bound of Thm 4, we make no assumption on the distribution of the cellwise contamination. These 2 thms reveal that even small contaminations impact the minimax rate.\n\nThe literature on robust estimation is very interesting. Some contributions consider weaker assumptions like a $L_4-L_2$ condition on the distribution of $X$ and no contamination. But this requires more complex covariance estimation procedures like one based on uniform trimming in [Abdalla and Zhivotovskiy, 2022] which are often computationally intensive. It would be interesting to explore this direction both theoretically and in practice within the cell-wise contamination setting with a focus on the computational/statistical trade-off. But this is beyond the scope of this paper.\n\nAbdalla and Zhivotovskiy (2022). Covariance Estimation: Optimal Dimension-free Guarantees for Adversarial Corruption and Heavy Tails (arXiv:2205.08494)\n\n2. We mostly agree. In a missing values setting, we can easily estimate the probability of missing values by counting the number of NA entries in a dataset and derive theoretical guarantees for it. But, in the contamination setting of Eq(2), a full theoretical analysis of our estimators including estimation of $\\delta$ and $\\epsilon$ is difficult. From a practical aspect, we have built estimation procedures for $\\delta$ and $\\epsilon$ which work reasonably well in experiments. See Tab 2 with $\\epsilon=1$ which shows that we can accurately estimate $\\delta$ in some cellwise contamination setting.\n \n3.\tPlease see our response to Q1.\n\nFinally, thanks for pointing out some typos. We have corrected them. We also added a mention to Figure 1 in our discussion on line 257 along with Figure 4.\n"}, {"review_id": "UAhZCZ3s58", "reviewer": "Reviewer_LnKP", "summary": "This work proposes an unbiased estimator for the covariance matrix under missing entries. Unlike most existing approaches, the proposed estimator does not require an imputation step to discard missing values before computing the covariance matrix. Interestingly, the proposed method employs an approach that avoids matrix inversion in high dimensions, thereby overcoming large complexity issues. \nMoreover, along with relevant outlier detection methods, the proposed estimator can be used used under the presence of cell-wise outliers to tackle such outliers in a high dimensional setting with inherent low rank structure. The effectiveness of the proposed method is demonstrated via experimental studies. ", "questions": "What is STD in Table 2?\n\nPlease provide an expression for the relative spectral difference used in Table 3.\n\n", "limitations": "NA", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Derivation of non-asymptotic estimation bounds of the covariance matrix with the operator norm and matching minimax lower bounds, enabling the quantification of impact of rates of missing values rate outliers.\n\nTraditional methods rely on matrix inversions resulting in failure in the high-dimensional setting, while the proposed avoids matrix multiplication and therefore, performs well. \n\nExperiments demonstrate that the proposed method is more robust to cell-wise contamination than comparing state-of-the-art methods, while producing reliable estimates of the covariance in low computation duration.\n\n", "weaknesses": "Missing definitions of metrics -- Please see below.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the referee for the careful reading of our paper. As recommended by the referee, we will precise the definitions of the metrics in the revision of this paper.\n\n- Q1. STD stands for standard deviation. We will precise it in the table caption.\n\n- Q2. The formula for the relative spectral difference of matrices $A$ and $B$ is :\n$$d(A,B) = \\frac{\\Vert A - B \\Vert}{max(\\Vert A \\Vert, \\Vert B \\Vert)}$$\ni.e. the ratio between the spectral difference and the largest spectral norm. We will add it to the revision of the paper.\n"}, {"review_id": "EoPNwDTBxZ", "reviewer": "Reviewer_zw68", "summary": "The objective of this paper is to mitigate the influence of outliers and missing values in covariance estimation without resorting to imputations.  A particular type of missing data that satisfies assumption (1) is considered and an unbiased estimator of the true covariance  that achieves the minimax bound with the operator norm is proposed.  For the outliers, they are further assumed to be cell-wise contaminated satisfying assumption (2) and two further assumptions (see lines 169 and  173).  Under these assumptions, an unbiased estimator of the true covariance is proposed, and its upper and lower bounds in the operator norm are established.  Numerical experiments demonstrate the advantages of the proposed approaches over existing approaches in a high-dimension and low rank setting. ", "questions": "Does equation (3) only involve missing values with no cell-wise contaminations as it does not involve any of the parameters in (2)?   \n\nThe \\xi_k in line 86 are random vectors with diagonal covariance matrix (see line 169).    Is this diagonal assumption a key to the results in Section 4? What would happen if they were not diagonal? \n\n", "limitations": "The proposed approach involves several assumptions that will limit its applicability in practice.\n\n", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This is a paper with both new methodology and theory for an interesting problem. The primary contribution is its strong theoretical results, which significantly  improve over existing results. ", "weaknesses": "The assumptions on missing values seem strong.  Usually, measurements for a subject are collected at the same time.  So, forcing the missing mechanism to be independent across the components seems to limit its applications. Also, why would all the components have the same missing probability? Some components might be harder to measure, hence more likely to be missing.  This assumption is much stronger than MCAR. \n\nThe definition of cell-wise contamination in (2) is restrictive as an outlier for a vector data does not have to be an outlier in any of the component. Also, it is further assumed in line 169 that components of cell-wise contaminations are uncorrelated. \n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the referee for the careful reading of our paper. \n\nWe discuss first the weaknesses identified by the referee.\n\n- W1. Regarding the first identified weakness, we decided to present first the classical MCAR setting for the sake of clarity. But we also covered the heterogeneous missingness case at the end of Section 3 (line 161). In that setting, each component has a different missingness probability and we provided a statistical guarantee in Eq. 8 (See also the proof in Appendix E.4). We did not give more details due to lack of space but in that case, the debiased covariance scheme is\n$$\n\\widetilde{\\Sigma} = \\mathrm{diag}(\\delta_1^{-1},\\cdots,\\delta_p^{-1}) \\odot \\mathrm{diag}(\\Sigma^Y) + A\\odot  \\left( \\Sigma^Y - \\mathrm{diag}(\\Sigma^Y) \\right),\n$$\nwhere the $p\\times p$ matrix $A$ admits entries  $A_{j,j}=0$ for any $j\\in [p]$ and off-diagonal entries\n$$\nA_{j,k} = \\frac{1}{\\delta_j\\delta_k},\\quad \\text{if $j\\neq k$}.\n$$\nHere $\\odot$ stands for the (Hadamard) entrywise product. \n\nWe will add more details on the heterogeneous missingness case in the appendix in the revision of the paper.\n\n\n- W2. Regarding the second identified weakness, we understand that the setting described by the referee may correspond to the classical Huber’s contamination setting that has been the focus of many studies in the literature. Our understanding is the cell-wise contamination setting and the Huber contamination settings are very different (one cannot be considered as a subcase of the other). As a consequence, methods developed for Huber contamination may not be adapted to cell-wise contamination and vice versa. We discussed this aspect on lines 25-31, but we would be happy to further precise this aspect in the revision if requested by the referee.\n\nWe answer now the questions raised by the referee. \n\n- Q1. Yes we consider the estimator in Eq (5) which corresponds to the case of missing values and no contamination. Later in the paper, we consider two strategies to handle the cell-wise contamination setting. The first is the estimator in Eq. (10) that we studied from a theoretical perspective in order to determine the minimax rate of estimation in the cell-wise contamination setting. The second strategy is a combination of filtering to detect and replace all contaminated entries by a $0$ and then apply the missing values estimator Eq. (5) on the filtered dataset. This second strategy is user-friendly and produces good results in our experimental study. It beats the SOTA method for covariance estimation with cell-wise contaminations ([1,4,21,22,23]). \n\n\n- Q2. Thanks for this great question. We considered in Eq (2) the setting where the cell-wise contamination happened at random independently of each other as it is common in this literature considered a realistic assumption to model sensor failures ([1,4,21,22,23]). The consequence is that the matrix $ \\Lambda $ is diagonal.\nComing back to your question, we believe that it may be possible to cover the case of non diagonal $\\Lambda$ by introducing an appropriate modification of Eq (10) and a similar (albeit more technical) analysis as the one we use to prove Theorem 3. We may not pursue this direction further as the estimator in Eq (10) was introduced solely for theory purposes to understand the minimax rate in the cell-wise contamination setting. As we mentioned above, our practical solution is rather to filter the dataset first and then apply the Missing Values estimator in Eq. (5).\n\nWe will add discussion on these aspects in the revision of our paper. Thanks!\n\n"}, {"review_id": "85KTs5mn3z", "reviewer": "Reviewer_PVnK", "summary": "The paper “Robust covariance estimation with missing values and cell-wise contamination”  presents an unbiased estimator for covariance in the presence of missing values that do not require any imputation step, and the authors claim they achieve competitive results compared to other state-of-the-art methods. The authors check their findings in synthetical and real-world datasets and compare them to other methods such as imputation using KNN or the MICE algorithm", "questions": "I still have some questions that could maybe help to make the paper more easy to understand:\n\n1. In Equation 2, why $X_i, \\psi_i, d_{i,j}, e_{i,j}$ are mutually independent? Isn’t this a rather strong assumption? \n2. The oracle using the MV method seems to always perform perfectly. More comments on this would be helpful to understand why. My first believe is that there might not be a fair comparison w.r.t the other methods. \n3. In Table 2: What are DDCMV90, DDCM95? I cannot find where these are defined. \n4. What is the so called “classical” covariance estimator? There is no definition of such method.", "limitations": "The authors don’t make any comment about the limitations of their method. However, I found a rather strong limitation in Appendix B. The method is tested on real-world datasets from sklearn, but the authors remove the categorical variables and variables that might be multimodal. This is a strong limitation since the method is tested on very smooth data distributions. I believe the authors should perform experiments in these datasets, assuming the natural artefacts that might appear.\n\nThe lack of a clear thread that makes the narrative of the paper easy to follow, the difficulty of understanding the final results from the experimental section because some figures do not appear and the need for further comparison of the proposed method with more SOTA models and other datasets drive me to believe that this paper lies below the acceptance threshold for NeurIPS conference at this point. I also hope some of my concerns can be addressed during the rebuttal period.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The notion of rigorously establishing and deriving bounds for covariance matrix estimation in the presence of missing data is intriguing, as it suggests that the bounds established for covariance estimation in complete data scenarios should extend to account for missing data. Hence, this contribution from the paper opens up new avenues of thought for researchers in the field of \"missing data.”", "weaknesses": "However, there are some substantial flaws in the paper that I would like to point out. First and probably most important, I consider the paper is not well structured and it is challenging to follow, even for experienced readers in the field of missing data handling. There is no evident thread from which the reader can follow the ideas presented throughout the paper. The connection between sections, paragraphs and ideas is often misleading: the reader might be introduced to new notations, and new concepts, or expected to remember previous concepts that were not actually introduced (Line 108 you comment $\\Sigma$ is low rank with no extra justification). Specially, after reading the paper the ultimate goal of the paper still remains unclear: is it providing a novel method for covariance estimation under missing data, imputing missing data directly? \n\nThe paper also lacks a clear section presenting the related work. These sections are often included in papers because: 1) they permit newcomers to the field to get some background, 2) they set the presented paper in the context of other historical or concurrent works, and outline the clear differences w.r.t others and 3) they allow experienced readers to corroborate if the proposed paper has properly undergone a prior deep literature review. I feel this paper lacks this important aspect: authors claim they compare to SOTA methods, but the most recent paper dates from 2021, then 2020, and after that 2018. Although MICE and KNN imputation are widely used in the missing data literature, the missing data is a rather hot topic in the machine learning field, and many works have proposed novel methods in the context of VAEs([1], [2], [3], [4]), GANs([5]) or recently diffusion models ([6]). I consider claiming this paper can outperform SOTA methods might be a strong assertion.\n\nThis, and the lack of comparison with standard datasets used in the missing data literature (tabular datasets from UCI or high-dimensional datasets such as MNIST, and CelebA, for example) makes the experimental section difficult to follow. The synthetic and real-world datasets are not described at all in the main manuscript and are referred to the Appendix, but the authors make comments about these datasets in the main text. From Figures 4 and 5 it is still not evident to me which method they propose, because it seems like they mix already proposed methods with their new covariance estimation. What’s more, the authors refer to the Figures and Tables in the Appendix and make comments about the results in the main text. I believe if comments about Figures are to be made, then these figures should be added in the main manuscript; or conversely, comments should be made in the Appendix. Something in between might make the reading a bit cumbersome. \n\n[1] Nazabal, A., Olmos, P. M., Ghahramani, Z., & Valera, I. (2020). Handling incomplete heterogeneous data using vaes. *Pattern Recognition*, *107*, 107501.\n\n[2] Mattei, P. & Frellsen, J.. (2019). MIWAE: Deep Generative Modelling and Imputation of Incomplete Data Sets. Proceedings of the 36th International Conference on Machine Learning, in Proceedings of Machine Learning Research 97:4413-4423 Available from https://proceedings.mlr.press/v97/mattei19a.html.\n\n[3] Ma, C., Tschiatschek, S., Palla, K., Hernández-Lobato, J. M., Nowozin, S., & Zhang, C. (2018). Eddi: Efficient dynamic discovery of high-value information with partial vae. *arXiv preprint arXiv:1809.11142*.Chicago\n\n[4] Ma, C., Tschiatschek, S., Turner, R., Hernández-Lobato, J. M., & Zhang, C. (2020). VAEM: a deep generative model for heterogeneous mixed type data. *Advances in Neural Information Processing Systems*, *33*, 11237-11247.\n\n[5] Yoon, J., Jordon, J., & Schaar, M. (2018, July). Gain: Missing data imputation using generative adversarial nets. In *International conference on machine learning* (pp. 5689-5698). PMLR.\n\n[6] Zheng, S., & Charoenphakdee, N. (2022). Diffusion models for missing value imputation in tabular data. Representation Learning Workshop at NeurIPS.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the referee for the detailed feedback. \n\nWe would like to address first the weaknesses pointed out by the referee.\n\n- The approximately low-rank assumption is standard in the literature on the theoretical analysis of covariance estimation. For many datasets including numerous UCI datasets, the covariances are approximately low-rank. \n\n-  _Purpose of our paper_. We certainly agree that robust covariance estimation with missing data is a different problem from data imputation. We will complete our literature discussion with the references the referee kindly provided and a message along these lines:\n\n“Data imputation is a fundamental field of research with the development of sophisticated solutions based on deep learning, GANs, VAE or Diffusion schemes [1,2,3,4,5] to perform complex tasks like artificial data generation or image inpainting. While the need for efficient imputation methods is undeniable for these applications, for the specific task of covariance estimation, we propose a simple, user-friendly debiasing scheme requiring no imputation and attaining the minimax optimal rate.”\n\nWe will be happy to complete this discussion as the referee sees fit regarding the data imputation literature.\n\nWe confirm that our goal is to study the robust covariance estimation problem when the dataset contains missing values and/or cell-wise contaminations. Hence we discussed SOTA methods for this specific problem. We do not treat in our work the fundamental problem of data imputation and we honestly believe that our paper entitled “Robust covariance estimation with missing values and cell-wise contamination” was not misleading on our purpose and contributions which we enumerated in the global response to the referees.\n\n\n\n\n- _Lack of comparison with standard datasets used in the missing data literature_. In Appendix B Table 5, we describe the 10 real datasets of various sizes and taken from diverse domains that we used in our benchmark. Note that 3 datasets taken from sklearn are actually from the UCI database. We could not place this description in the main body. Indeed, given the 9 pages limitation and our main objective to present our theoretical findings, we made the choice to place this description in Appendix B as it is quite common for theoretical papers submitted at ML conferences.\n\n- _Compared methods in Figures 4 and 5_. Our practical method for cell-wise contamination comprises 2 steps: 1st. filtering the dataset 2nd applying Estimator in Eq (5) called MV as we described on line 212. We applied different filtering procedures in combination with MV in Fig 4 and 5 as we described in lines 240-243. \n\n- _Discussion of additional material in Appendix_. We agree. In the revision, we will move the comments next to their figures in Appendix whenever required.\n\n\nWe address now the questions raised by the referee.\n\n- Q1. These assumptions are standard in the literature on cell-wise contamination like in [Agostinelli et al. (2014), Eq. (2)], especially when we investigate the theoretical properties of this problem. Note that these assumptions are relaxed in Theorem 4 where the contaminations $\\xi$ admit arbitrary distributions.\n\nAgostinelli, Leung, Yohai, and Zamar (2014). Robust estimation of multivariate location and scatter in the presence of cellwise and casewise contamination\n\n- Q2. This may not be easy to see with the naked eye, but the performance of oracleMV does degrade as the contamination parameter increases but slower than for the other methods. Hence it may appear constant for contamination probability in the range $(0, 0.4)$. In the rebuttal pdf file, we included a figure where we extended the x-axis with contamination probability ranging from $0$ to $1$. We can see more clearly that the performance of oracleMV deteriorates as the contamination probability increases as we should expect from theory.\nNote however that in the literature on robust estimation with contamination, it is proved that this problem is hopeless when the contamination probability exceeds 50% since in that case, the true observations are in minority. This explains why we did not go beyond the range (0, 0.5) in our Exps. In practice, methods capable of handling 20% or more cellwise contaminations are considered quite good and of practical interest. Remarkably, our DDCMV method is doing significantly better than the classical sample covariance estimator for contamination probabilities up to 35%. Conversely, other methods experience breakdowns at contamination levels between 15% and 18%. \n \n- Q3. As described on lines 240-243, we implemented the DDC filtering combined with MV (Eq 5). 90% or 95% correspond to the trimming parameter of the DDC filtering method.\n\n\n- Q4. This is the classical empirical covariance estimator that we introduced on line 117:\n$\\hat{\\Sigma}^Y = \\sum_{i=1}^n Y_i \\otimes Y_i/n$\n\n\nRegarding the limitation on quantitative features raised by the referee, we would like to clarify that our paper concerns robust covariance estimation, primarily with applications like PCA in mind. Therefore, we focused solely on quantitative features. For studying relationships between quantitative and categorical features (e.g., color of the eyes), alternative methods like ANOVA tests would be more appropriate, as computing the covariance between a quantitative and a categorical covariate does not make sense.\n"}]

Author rebuttal:
We wish to thank the reviewers for their insightful evaluation of our paper. We appreciate all the comments and remarks, which we will implement in our revision. We discuss here some points that were raised by several referees. 

The goal of our paper is to study the robust covariance estimation problem when the dataset contains missing values and/or cell-wise contaminations. We clarify some important aspects of our contributions:

- __For the specific problem of robust covariance estimation with missing values__, we prove (in Thms 1 and 2) that a debiasing scheme without any imputation can provide an optimal estimation procedure with theoretical minimax guarantees. This minimax result firmly establishes the fact that there cannot exist any method (even one relying on data imputation) that achieves a rate faster than the minimax rate. In other words, imputation based methods may also achieve results on par with but cannot surpass the performance of our simple method in terms of minimax optimality. Furthermore, the imputation methods we implemented incurred a higher computational cost.

- For the cell-wise contamination setting, we prove theoretical results (Thm 3 and 4) on the impact of contaminations on estimation rates. We advocate that a simple solution that performs filtering of the dataset first and then applying a simple debiasing scheme for the missing values case performs better than the SOTA methods for the robust covariance estimation problem with cell-wise contaminations ([1,4,21,22,23]). As we observed in our experiments, these SOTA methods perform poorly and are numerically unstable even on simple datasets.

- On the practical aspect, our solutions are user-friendly and computationally fast.


- As requested by one referee, we added in the attached pdf file a completed Figure 4 with missing probability in the range (0,1). We can see that the accuracy of oracle MV does degrades as missing probability increases to 1 as predicted by theory.

As requested by one referee, we added in the attached pdf file a completed Figure 4 with missing probability in the range (0,1). We can see that the accuracy of oracle MV does degrades as missing probability increases to 1 as predicted by theory.


Meta-review: This paper studies the task of covariance estimation in the presence of missing data and cell-wise contamination. It provides an unbiased estimator that can achieve minimax optimal statistical accuracy. In particular, the authors show that their provided estimator achieves a non-asymptotically optimal statistical error, provided that both the missing values and the contamination model follow a Bernoulli model. The authors then complement their theoretical results with comprehensive experiments in both low and high-dimensional regimes.

Four reviewers have reviewed the paper, and their overall assessment of the paper was positive. I agree with this assessment and believe the paper is a solid contribution with interesting results. I should also mention that the authors did a great job addressing the comments raised by the reviewers. I encourage the authors to address these comments and add the complementary simulations to the revised paper.
Decision: Accept (poster)