# Robust covariance estimation and explainable outlier detection for matrix-
valued data

Marcus Mayrhofer
TU Wien
<EMAIL_ADDRESS>
& Una Radojičić
TU Wien
<EMAIL_ADDRESS>
& Peter Filzmoser
TU Wien
<EMAIL_ADDRESS>

###### Abstract

The minimum covariance determinant (MCD) estimator is a popular method for
robustly estimating the mean and covariance of multivariate data. We extend
the MCD to the setting where the observations are matrices rather than vectors
and introduce the matrix minimum covariance determinant (MMCD) estimators for
robust parameter estimation. These estimators hold equivariance properties,
achieve a high breakdown point, and are consistent under elliptical matrix-
variate distributions. We have also developed an efficient algorithm with
convergence guarantees to compute the MMCD estimators. Using the MMCD
estimators, we can compute robust Mahalanobis distances that can be used for
outlier detection. Those distances can be decomposed into outlyingness
contributions from each cell, row, or column of a matrix-variate observation
using Shapley values, a concept for outlier explanation recently introduced in
the multivariate setting. Simulations and examples reveal the excellent
properties and usefulness of the robust estimators.

_Keywords_ Covariance with Kronecker structure $\cdot$ Explainable artificial
intelligence $\cdot$ Image data $\cdot$ Minimum covariance determinant $\cdot$
Shapley values

##  1 Introduction

Thanks to modern data collection tools, the amount and complexity of available
information are increasing rapidly, and often, matrix-valued data are
observed. Compared to classical multivariate observations, where values for
$p$ variables are recorded for one subject, matrix-valued observations are
recorded on a grid of $p\times q$ variables. These are then naturally
represented as a matrix with $p$ rows and $q$ columns. Some examples include
image data, where $p$ and $q$ are given by the resolution of the image, or
multivariate data measured on $p$ variables, where the measurements for a
subject are available for $q$ replications (e.g., different time points,
different spatial locations, different experimental conditions, etc.).
Frequently, matrix-valued data are analyzed as classical multivariate data by
stacking the matrix columns (or rows) to a vector of length $p\cdot q$. Thus,
if $n$ observations are available, the data are arranged in a matrix of
dimension $n\times pq$. Depending on the dimensions, this can create high-
dimensional data, possibly with a sample size lower than the resulting
dimensionality, which constitutes a limitation for several multivariate
statistical methods.

As an alternative to vectorizing matrix-valued observations, we model them
under the assumption that they originate from a certain matrix-variate
distribution. As in the multivariate case, elliptical matrix-variate
distributions are studied frequently (Gupta and Varga,, 2012), with particular
focus being placed on the most prominent representative, i.e., the matrix
normal distribution (Dawid,, 1981). Formally, a random matrix
$\bm{X}\in\mathbb{R}^{p\times q}$ follows a matrix normal distribution
$\mathcal{MN}(\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}})$
with mean $\bm{\operatorname{M}}\in\mathbb{R}^{p\times q}$, row covariance
$\bm{\Sigma}^{\operatorname{row}}\in\operatorname{PSD}(p)$, and column
covariance $\bm{\Sigma}^{\operatorname{col}}\in\operatorname{PSD}(q)$, if and
only if its vectorized version $\operatorname{vec}(\bm{X})\in\mathbb{R}^{pq}$
follows a multivariate normal distribution
$\mathcal{N}(\operatorname{vec}(\bm{\operatorname{M}}),\bm{\Sigma}^{\operatorname{col}}\otimes\bm{\Sigma}^{\operatorname{row}})$
(Gupta and Nagar,, 1999). Here, $\operatorname{vec}(\cdot)$ is the
vectorization operator, stacking the columns of a matrix on top of each other,
$\otimes$ is the Kronecker product, and $\operatorname{PSD}(a)$, with
$a\in\mathbb{N}$, denotes the class of all positive definite symmetric
$a\times a$ matrices. The probability density function (pdf) of a matrix
normal random variable $\bm{X}$ is given by

$\displaystyle
f(\bm{X}|\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}})=\frac{\exp(-\frac{1}{2}\operatorname{tr}(\bm{\Omega}^{\operatorname{col}}(\bm{X}-\bm{\operatorname{M}})^{\prime}\bm{\Omega}^{\operatorname{row}}(\bm{X}-\bm{\operatorname{M}})))}{(2\pi)^{\nicefrac{{pq}}{{2}}}\det(\bm{\Sigma}^{\operatorname{col}})^{\nicefrac{{p}}{{2}}}\det(\bm{\Sigma}^{\operatorname{row}})^{\nicefrac{{q}}{{2}}}},$
(1)

with
$\bm{\Omega}^{\operatorname{row}}=(\bm{\Sigma}^{\operatorname{row}})^{-1}$ and
$\bm{\Omega}^{\operatorname{col}}=(\bm{\Sigma}^{\operatorname{col}})^{-1}$
denoting the precision matrices among the rows and columns, respectively.

Regarding the estimation of location and covariance for an i.i.d. sample
$\bm{\mathfrak{X}}=(\bm{X}_{1},\ldots,\bm{X}_{n})\in\mathbb{R}^{n\times
p\times q}$, with
$\bm{X}_{i}\sim\mathcal{MN}(\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}})$,
we can either work with the vectorized observations or directly with the
matrices. In the former setting, the existence and uniqueness of the maximum
likelihood estimator (MLE) for the covariance is guaranteed almost surely if
$n\geq pq+1$. However, this approach does not take advantage of the Kronecker
structure of the covariance matrix and instead directly estimates the entire
$pq$-dimensional matrix $\bm{\Sigma}$. In contrast, if we utilize the
knowledge of the inherent data structure, we only need to estimate the
$p$-dimensional rowwise covariance matrix $\bm{\Sigma}^{\operatorname{row}}$
and the $q$-dimensional columnwise covariance matrix
$\bm{\Sigma}^{\operatorname{col}}$. For the matrix-variate sample
$\bm{\mathfrak{X}}$, the MLEs for the mean, as well as for the rowwise and
columnwise covariance are given by (Dutilleul,, 1999):

$\displaystyle\hat{\bm{\operatorname{M}}}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\bm{X}_{i}$ (2)
$\displaystyle\hat{\bm{\Sigma}}^{\operatorname{row}}$
$\displaystyle=\frac{1}{qn}\sum_{i=1}^{n}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}})\hat{\bm{\Omega}}^{\operatorname{col}}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}})^{\prime}$
(3) $\displaystyle\hat{\bm{\Sigma}}^{\operatorname{col}}$
$\displaystyle=\frac{1}{pn}\sum_{i=1}^{n}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}})^{\prime}\hat{\bm{\Omega}}^{\operatorname{row}}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}})$
(4)

Soloveychik and Trushin, (2016) showed that for $n$ i.i.d. samples from a
continuous $p\times q$ matrix-variate distribution, there exists no unique
maximum of the matrix normal likelihood function if
$n<\max(\nicefrac{{p}}{{q}},\nicefrac{{q}}{{p}})+1$, and that a unique maximum
exists almost surely if
$n\geq\lfloor\nicefrac{{p}}{{q}}+\nicefrac{{q}}{{p}}\rfloor+2$. While there
exists no closed-form solution for the MLEs of
$\bm{\Sigma}^{\operatorname{row}}$ and $\bm{\Sigma}^{\operatorname{col}}$, an
iterative estimation procedure was proposed by Dutilleul, (1999). The idea of
the so-called _flip-flop_ algorithm is to alternate between the computation of
$\hat{\bm{\Sigma}}^{\operatorname{row}}$ and
$\hat{\bm{\Sigma}}^{\operatorname{col}}$ based on Equations (3) and (4),
respectively, until a convergence criterion is met. The algorithm is
constructed such that positive definite estimates of subsequent iterations are
nondecreasing in likelihood (Lu and Zimmerman,, 2005) and it converges almost
surely to the unique maximum from any symmetric positive definite
initialization of either $\hat{\bm{\Sigma}}^{\operatorname{row}}$ or
$\hat{\bm{\Sigma}}^{\operatorname{col}}$, if
$n\geq\lfloor\nicefrac{{p}}{{q}}+\nicefrac{{q}}{{p}}\rfloor+2$ (Soloveychik
and Trushin,, 2016).

In this paper, we are interested in a robust estimation of the parameters
$\bm{\operatorname{M}}$, $\bm{\Sigma}^{\operatorname{row}}$, and
$\bm{\Sigma}^{\operatorname{col}}$ of the matrix normal distribution. Existing
proposals include a generalization of Tyler’s M-estimator (Tyler,, 1987)
introduced by Soloveychik and Trushin, (2016), a robust estimator for
structured covariance matrices with Kronecker structure as a particular case
(Sun et al.,, 2016), distribution-free robust covariance estimation (Zhang et
al.,, 2022), and ML estimation for the matrix t-distribution (Thompson et
al.,, 2020).

We propose a new robust estimation method, the matrix-variate minimum
covariance determinant (MMCD) estimator. This is an extension of the minimum
covariance determinant (MCD) estimator (Rousseeuw,, 1985), a popular robust
estimator of location and covariance for the case of multivariate data. We
will show that the MMCD is equivariant under linear matrix transformations,
consistent under elliptical distributions, and achieves a high breakdown
point. It also achieves a higher breakdown point than the MCD when applied to
vectorized matrix-valued observations. Furthermore, an algorithm similar to
the Fast-MCD algorithm (Rousseeuw and Driessen,, 1999) is developed for the
computation of the MMCD, which makes the methodology also accessible to
higher-dimensional matrices.

The robust MMCD estimators can then be employed for outlier detection using
the Mahalanobis distances (Mahalanobis,, 1936) for matrix-valued observations.
Because it is essential to understand the reasons for the outlyingness, we
will extend the concept of Shapley values introduced in Mayrhofer and
Filzmoser, (2023) for outlier explanation in the multivariate case to the
matrix-variate setting. Shapley values (Shapley,, 1953) are well-known from
explainable AI (Lundberg and Lee,, 2017), but their computation is usually
time-consuming. Our proposal is computationally efficient, and the resulting
Shapley values preserve their attractive properties (Shapley,, 1953).

The paper is organized as follows. In Section 2, we introduce the MMCD
estimators, then proceed to derive their theoretical properties in Section 3.
Section 4 is devoted to computational details for the MMCD estimators. In
Section 5, we propose Shapley values for outlier explanation and present their
properties. In Sections 6 and 7, we illustrate the performance of the proposed
methods on numerical simulations and real-world examples. Section 8 concludes
our findings. The Supplement contains more theoretical background information
in this context, proofs, technical derivations, and additional numerical
results.

##  2 The MMCD estimators

The goal is to derive robust estimators for the parameters
$\bm{\operatorname{M}}$, $\bm{\Sigma}^{\operatorname{row}}$, and
$\bm{\Sigma}^{\operatorname{col}}$ in Equation (1) mimicking the approach of
the MCD estimator, whose objective is to identify the subset of $h$ out of $n$
samples, with $\nicefrac{{n}}{{2}}\leq h\leq n$, whose sample covariance
matrix has the smallest determinant. For the MMCD estimators, we use the
matrix MLE estimators instead of the sample covariance of the vectorized data.
Therefore, we use the connection between the MCD estimator and the MLEs for
mean and covariance of the multivariate normal distribution, as pointed out in
Raymaekers and Rousseeuw, (2023): Minimizing the determinant of the sample
covariance matrix of a data subset is equivalent to minimizing the
(log-)likelihood function for that subset. The latter corresponds to
minimizing a trimmed (log-)likelihood, where trimming is done by introducing
binary weights for the observations. This idea has been similarly used in the
context of other robust estimators (e.g. Neykov et al.,, 2007; García-Escudero
et al.,, 2010; Kurnaz et al.,, 2018). More details regarding the matrix-
variate setting and the MCD procedure are provided in Supplement A.

The MLEs given in Equations (2)-(4) are not robust against outlying (matrix-
variate) observations. Thus, we introduce weights
$\bm{w}=(w_{1},\dots,w_{n})\in\mathbb{R}^{n}$ for a given sample
$\bm{\mathfrak{X}}=(\bm{X}_{1},\ldots,\bm{X}_{n})$ that is independently drawn
from
$\mathcal{MN}(\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}})$
to formulate the weighted log-likelihood function
$l(\bm{w},\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}}|\bm{\mathfrak{X}})$
as

$\displaystyle-\frac{1}{2}\sum_{i=1}^{n}w_{i}\Bigl{(}p\ln(\det(\bm{\Sigma}^{\operatorname{col}}))+q\ln(\det(\bm{\Sigma}^{\operatorname{row}}))+\operatorname{MMD}^{2}(\bm{X})+pq\ln(2\pi)\Bigr{)},$
(5)

where $\operatorname{MMD}^{2}(\bm{X})$ denotes the squared matrix Mahlanobis
distance defined as

$\displaystyle\operatorname{MMD}^{2}(\bm{X}):=\operatorname{MMD}^{2}(\bm{X};\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}})=\operatorname{tr}(\bm{\Omega}^{\operatorname{col}}(\bm{X}-\bm{\operatorname{M}})^{\prime}\bm{\Omega}^{\operatorname{row}}(\bm{X}-\bm{\operatorname{M}})).$
(6)

Setting $w_{i}=1$ for all $i=1,\ldots,n$, yields the traditional log-
likelihood function, and its maximization yields the MLEs of Equation (2)-(4).
However, by taking binary weights, $w_{i}\in\\{0,1\\}$, with the constraint
that $\sum_{i=1}^{n}w_{i}=h$, we see that $n-h$ contributions are trimmed.
Since contributions from outliers should be trimmed, the task is to identify
the subset of regular observations $H\subset\\{1,\ldots,n\\}$ with $|H|=h$,
where $w_{i}=1$ for $i\in H$ and $0$ otherwise. The resulting constrained
optimization problem of finding the weighted MLE can be written as

$\displaystyle\begin{split}&\max_{\bm{w},\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}}}l(\bm{w},\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}}|\bm{\mathfrak{X}})\\\
&\text{subject to}\quad w_{i}\in\\{0,1\\}\text{ for all
}i=1,\dots,n\quad\text{and}\quad\sum_{i=1}^{n}w_{i}=h.\end{split}$ (7)

To improve clarity, we will use the following notation for subsamples of
$\bm{\mathfrak{X}}$ and estimators based on it: Let
$H\subseteq\\{1,\dots,n\\}$ be a subset of size $h=\left|H\right|$, then
$\bm{\mathfrak{X}}_{H}:=(\bm{X}_{i})_{i\in H}$ denotes an $h$-subset of
$\bm{\mathfrak{X}}$. An estimator for a parameter $\theta$ based on the sample
$\bm{\mathfrak{X}}$ is denoted as $\hat{\theta}_{\bm{\mathfrak{X}}}$ or simply
as $\hat{\theta}$ if it is clear on which sample the estimator is computed.
Similarly, if an estimator is based on an $h$-subset, it is denoted as
$\hat{\theta}_{H}$ or as $\hat{\theta}_{\bm{\mathfrak{X}}_{H}}$.

######  Proposition 2.0.1.

Let $\bm{\mathfrak{X}}=(\bm{X}_{1},\dots,\bm{X}_{n})$,
$\nicefrac{{n}}{{2}}\leq h\leq n$ and
$h\geq\lfloor\nicefrac{{p}}{{q}}+\nicefrac{{q}}{{p}}\rfloor+2$, be an i.i.d.
sample from
$\mathcal{MN}(\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}})$.
Maximizing the weighted log-likelihood function (7) is equivalent to
minimizing

$\displaystyle\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{col}}_{H}\otimes\hat{\bm{\Sigma}}^{\operatorname{row}}_{H}))=p\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{col}}_{H}))+q\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{row}}_{H}))$
(8)

across all subsets $H\subset\\{1,\dots,n\\}$ with $|H|=h$. In Equation (8),

$\displaystyle\hat{\bm{\operatorname{M}}}_{H}$
$\displaystyle=\frac{1}{h}\sum_{i\in H}\bm{X}_{i},$ (9)
$\displaystyle\hat{\bm{\Sigma}}^{\operatorname{row}}_{H}$
$\displaystyle=\frac{1}{qh}\sum_{i\in
H}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})\hat{\bm{\Omega}}^{\operatorname{col}}_{H}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})^{\prime}\text{,
and}$ (10) $\displaystyle\hat{\bm{\Sigma}}^{\operatorname{col}}_{H}$
$\displaystyle=\frac{1}{ph}\sum_{i\in
H}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})^{\prime}\hat{\bm{\Omega}}^{\operatorname{row}}_{H}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})$
(11)

denote the MLEs based on the observations in $H$, and
$\hat{\bm{\Omega}}^{\operatorname{row}}_{H}=(\hat{\bm{\Sigma}}^{\operatorname{row}}_{H})^{-1}$
and
$\hat{\bm{\Omega}}^{\operatorname{col}}_{H}=(\hat{\bm{\Sigma}}^{\operatorname{col}}_{H})^{-1}$
denote the corresponding precision matrices.

A proof is given in Supplement B. Based on this proposition, we obtain the
matrix-variate counterpart to the multivariate MCD estimator’s objective,
resulting in robust estimators of the parameters $\bm{\operatorname{M}}$,
$\bm{\Sigma}^{\operatorname{row}}$, and $\bm{\Sigma}^{\operatorname{col}}$.

######  Definition 2.0.1.

Let $\bm{\mathfrak{X}}=(\bm{X}_{1},\dots,\bm{X}_{n})$,
$\nicefrac{{n}}{{2}}\leq h\leq n$ and
$h\geq\lfloor\nicefrac{{p}}{{q}}+\nicefrac{{q}}{{p}}\rfloor+2$, be an i.i.d
sample of a continuous $p\times q$ matrix-variate distribution. The _raw_
matrix minimum covariance determinant (MMCD) estimators are defined as

$\displaystyle(\hat{\bm{\operatorname{M}}}_{H^{\ast}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{H^{\ast}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{H^{\ast}}):=\operatorname*{arg\,min}_{\begin{subarray}{c}\hat{\bm{\operatorname{M}}}_{H},\hat{\bm{\Sigma}}^{\operatorname{row}}_{H},\hat{\bm{\Sigma}}^{\operatorname{col}}_{H}\\\
H\subset\\{1,\dots,n\\},\left|H\right|=h\end{subarray}}p\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{col}}_{H}))+q\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{row}}_{H})),$
(12)

with $\hat{\bm{\operatorname{M}}}_{H}$,
$\hat{\bm{\Sigma}}^{\operatorname{row}}_{H}$, and
$\hat{\bm{\Sigma}}^{\operatorname{col}}_{H}$ as in Equations (9), (10), and
(11), respectively.

Because the estimators in Definition (2.0.1) are MLEs based on the
observations in the subset $H$, they almost surely exist and are positive
definite if $h\geq\lfloor\nicefrac{{p}}{{q}}+\nicefrac{{q}}{{p}}\rfloor+2$. If
$p=1$ and/or $q=1$, optimization problem (12) coincides with the optimization
problem of the MCD estimator, and one obtains the univariate or multivariate
MCD estimator, respectively (Rousseeuw,, 1985).

##  3 Properties of the MMCD estimators

As outlined in Lopuhaa and Rousseeuw, (1991), the MCD estimator has the
highest possible breakdown point an affine equivariant estimator can attain.
Now that we have generalized the MCD estimator to the matrix-variate setting,
we want to study how the breakdown point and equivariance properties translate
to this setting.

Affine equivariance in the multivariate setting can be motivated by the fact
that location and covariance estimators should transform in the same way under
affine transformations as the parameters of elliptically symmetrical unimodal
distributions (from here on called elliptical distributions), see Maronna et
al., (2019). We can define the matrix-variate analog of affine equivariance
based on the properties of elliptical distributions, which are frequently
employed to study robustness properties of normal theory under nonnormal
situations (Gupta and Nagar,, 1999).

Using the same notation and parameters as for the matrix normal distribution,
a random matrix $\bm{X}$ with an absolutely continuous distribution has an
elliptical distribution, denoted
$\mathcal{ME}(\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}})$,
if its density can be written as

$\displaystyle
f(\bm{X})=\det(\bm{\Sigma}^{\operatorname{row}})^{\nicefrac{{q}}{{2}}}\det(\bm{\Sigma}^{\operatorname{col}})^{\nicefrac{{p}}{{2}}}g(\operatorname{tr}(\bm{\Omega}^{\operatorname{col}}(\bm{X}-\bm{\operatorname{M}})^{\prime}\bm{\Omega}^{\operatorname{row}}(\bm{X}-\bm{\operatorname{M}}))),$
(13)

with any positive $g$ such that $f(\bm{X})$ integrates to one. Linear
functions of a random matrix
$\bm{X}\sim\mathcal{ME}(\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}})$
also have an elliptical distribution (Gupta and Varga,, 2012). This means that
for constant matrices $\bm{A}\in\mathbb{R}^{r\times p}$,
$\operatorname{rank}(\bm{A})=r\leq p$, $\bm{B}\in\mathbb{R}^{q\times s}$,
$\operatorname{rank}(\bm{B})=s\leq q$, and $\bm{C}\in\mathbb{R}^{r\times s}$,
the transformed random matrix $\bm{Z}=\bm{A}\bm{X}\bm{B}+\bm{C}$ has density

$\displaystyle\bm{Z}\sim
E_{r,s}(\bm{A}\bm{\operatorname{M}}\bm{B}+\bm{C},\bm{A}\bm{\Sigma}^{\operatorname{row}}\bm{A}^{\prime},\bm{B}^{\prime}\bm{\Sigma}^{\operatorname{col}}\bm{B}).$
(14)

Let $\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{X}}}$,
$\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{X}}}$, and
$\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{X}}}$ denote the
estimators based on a sample $\bm{\mathfrak{X}}=(\bm{X}_{1},\dots,\bm{X}_{n})$
generated by
$f(\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}})$.
Then the estimators of the sample
$\bm{\mathfrak{Z}}=(\bm{A}\bm{X}_{1}\bm{B}+\bm{C},\dots,\bm{A}\bm{X}_{n}\bm{B}+\bm{C})$
should transform in the same way as the parameters in Equation (14), i.e.,

$\displaystyle\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Z}}}=\bm{A}\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{X}}}\bm{B}+\bm{C},\quad\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Z}}}=\bm{A}\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{X}}}\bm{A}^{\prime},\quad\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Z}}}=\bm{B}^{\prime}\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{X}}}\bm{B}.$
(15)

Properties (15) provide a suitable generalization of affine equivariance to
the matrix-variate setting, and it is easy to verify that they hold for the
MLEs. However, they do not imply affine equivariance of the location and
covariance estimators for the vectorized observations. This would only hold
for transformations with the Kronecker structure
$\operatorname{vec}(\bm{A}\bm{X}\bm{B}+\bm{C})=(\bm{B}^{\prime}\otimes\bm{A})\operatorname{vec}(\bm{X})+\operatorname{vec}(\bm{C})$.
To avoid confounding definitions, we refer to Properties (15) as _matrix
affine equivariance_.

######  Lemma 3.0.1.

Let $\bm{\mathfrak{X}}=(\bm{X}_{1},\dots,\bm{X}_{n})$ be a sample of $p\times
q$ matrices, where
$\bm{X}_{i}\sim\mathcal{MN}(\bm{\operatorname{M}}_{\bm{\mathfrak{X}}},\bm{\Sigma}^{\operatorname{row}}_{\bm{\mathfrak{X}}},\bm{\Sigma}^{\operatorname{col}}_{\bm{\mathfrak{X}}})$,
and let $\bm{\mathfrak{Z}}=(\bm{Z}_{1},\dots,\bm{Z}_{n})$ be the affine
transformation of $\bm{\mathfrak{X}}$, i.e.,
$\bm{Z}_{i}=\bm{A}\bm{X}_{i}\bm{B}+\bm{C}$, $\bm{A}\in\operatorname{PSD}(p)$,
$\bm{B}\in\operatorname{PSD}(q)$, and $\bm{C}\in\mathbb{R}^{p\times q}$. The
following then holds:

  1. (a)

The MMCD estimators as in Definition 2.0.1 are matrix affine equivariant.

  2. (b)

$\operatorname{MMD}^{2}(\bm{Z}_{i};\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Z}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Z}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Z}}})=\operatorname{MMD}^{2}(\bm{X}_{i};\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{X}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{X}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{X}}})$,
where
$(\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Z}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Z}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Z}}})$
are matrix affine equivariant location and covariance estimators of the
transformed sample $\bm{\mathfrak{Z}}$.

Lemma 3.0.1 shows that the MMCD estimators are equivariant under linear matrix
transformations, and a proof is given in Supplement B.1.

The finite sample breakdown point of an estimator evaluates its resilience to
contamination. It refers to the largest proportion of observations that may be
arbitrarily replaced by outliers such that the estimator still contains some
information about the true parameter (Maronna et al.,, 2019). Let
$\bm{\mathfrak{X}}$ be a sample of $n$ matrix-variate observations in
$\mathbb{R}^{p\times q}$ and suppose $\bm{\mathfrak{Y}}$ is a corrupted
version, obtained by replacing $m$ samples of $\bm{\mathfrak{X}}$ by arbitrary
matrices. The finite sample breakdown point of a location estimator
$\hat{\bm{\operatorname{M}}}$ is given by

$\displaystyle\varepsilon^{\ast}(\hat{\bm{\operatorname{M}}},\bm{\mathfrak{X}})=\min_{1\leq
m\leq
n}\left\\{\frac{m}{n}:\sup_{m}\left\lVert\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{X}}}-\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}}\right\rVert=\infty\right\\}$
(16)

and the (joint) finite sample breakdown point of row and columnwise covariance
estimators $\hat{\bm{\Sigma}}^{\operatorname{row}}$ and
$\hat{\bm{\Sigma}}^{\operatorname{col}}$ is given by

$\displaystyle\varepsilon^{\ast}(\hat{\bm{\Sigma}}^{\operatorname{row}},\hat{\bm{\Sigma}}^{\operatorname{col}},\bm{\mathfrak{X}})=\min_{1\leq
m\leq
n}\left\\{\frac{m}{n}:\sup_{m}\max_{i,j}\left|\log(\lambda_{i}(\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Y}}})\lambda_{j}(\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Y}}}))-\log(\lambda_{i}(\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{X}}})\lambda_{j}(\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{X}}}))\right|=\infty\right\\}.$
(17)

While the MCD and the MMCD estimators coincide for the case that $p=1$ and/or
$q=1$, the following theorem shows that the MMCD estimators achieve a higher
breakdown point than the MCD estimators applied to the vectorized samples if
$p\geq 2$ and $q\geq 2$.

######  Theorem 3.0.1.

Let $\bm{\mathfrak{X}}$ be a collection of $n$ i.i.d. samples from a
continuous $p\times q$ matrix-variate distribution, where
$d=\lfloor\nicefrac{{p}}{{q}}+\nicefrac{{q}}{{p}}\rfloor$,
$p,q\in\mathbb{N},p\geq 2,q\geq 2$, and let $\hat{\bm{\operatorname{M}}}$,
$\hat{\bm{\Sigma}}^{\operatorname{row}}$, and
$\hat{\bm{\Sigma}}^{\operatorname{col}}$ denote the MMCD estimators, then

$\displaystyle\varepsilon^{\ast}(\hat{\bm{\operatorname{M}}},\bm{\mathfrak{X}})=\varepsilon^{\ast}(\hat{\bm{\Sigma}}^{\operatorname{row}},\hat{\bm{\Sigma}}^{\operatorname{col}},\bm{\mathfrak{X}})=\frac{1}{n}\lfloor\min(n-h+1,h-(d+1))\rfloor=:\frac{m}{n},$

with $\nicefrac{{n}}{{2}}\leq h\leq n$ and $h\geq d+2$.

A proof building on similar methodologies established in Rousseeuw, (1985) and
Lopuhaa and Rousseeuw, (1991) is given in Supplement B.1. Theorem 3.0.1
implies that the maximum breakdown point of the MMCD estimators is
$\nicefrac{{1}}{{n}}\lfloor\nicefrac{{(n-d)}}{{2}}\rfloor$ and is attained if
$h=\lfloor\nicefrac{{(n+d+2)}}{{2}}\rfloor$. This implies that the maximum
breakdown point of the MMCD covariance estimators for $p\geq 2,q\geq 2$ is
higher than the upper bound for the breakdown point of affine equivariant
covariance estimators applied to vectorized samples, which is given by
$\nicefrac{{1}}{{n}}\lfloor\nicefrac{{(n-pq+1)}}{{2}}\rfloor$ (Davies,, 1987;
Lopuhaa and Rousseeuw,, 1991). However, as mentioned earlier, affine
equivariance in the matrix-variate setting does not imply affine equivariance
in the multivariate setting. Thus, the mentioned upper bound for the
vectorized observations does not apply. In other words, since affine
equivariance (in the vectorized case) is not a requirement for matrix-variate
affine equivariance, it is possible to achieve a higher breakdown point for
the MMCD estimators than for the MCD estimators of vectorized samples.

To illustrate the advantage of respecting the inherent data structure of
matrix-variate data for the breakdown properties, we compare the maximum
breakdown points of the MCD and MMCD estimators in Figure 1 for different
combinations of $p$ and $q$, and for different sample sizes $n$. Here, the MCD
estimator is applied to the vectorized data, and the dimensionality of the
samples is $pq$, which can get large. This affects the computability of the
MCD estimator since it requires a subset size larger than the dimension.

Figure 1: Comparison of the maximum breakdown point of the MMCD estimators for
matrix-variate data with $p=1,5,10,20$ rows and $q=20$ columns, and the MCD
estimator applied to the vectorized data. When $p=1$, both estimators and
their breakdown points coincide. However, increasing the number of rows yields
better breakdown properties for the MMCD estimators, as the proportion between
the number of rows and columns
$d=\lfloor\nicefrac{{p}}{{q}}+\nicefrac{{q}}{{p}}\rfloor$ is approaching 2.

By multiplying the raw MMCD rowwise or columnwise covariance estimator by a
consistency factor $c(\alpha)$, we can obtain consistency at the normal model
as outlined in Theorem 3.0.2. The proof of this theorem shows that the
consistency factor of the MCD estimator and the MMCD estimator must coincide,
and therefore, we use the consistency factor

$\displaystyle
c(\alpha)=\frac{\alpha}{F_{\chi^{2}_{p+2}}(\chi^{2}_{\alpha;pq})}$ (18)

proposed by Croux and Haesbroeck, (1999), where $F_{\chi^{2}_{p+2}}$ denotes
the CDF of the chi-square distribution with $p+2$ degrees of freedom, and
$\chi^{2}_{pq;\alpha}$ denotes the $\alpha$ quantile of the chi-square
distribution with $pq$ degrees of freedom.

######  Theorem 3.0.2.

Let $\bm{X}_{1},\dots,\bm{X}_{n}$ be a random sample from an elliptical
matrix-variate distribution
$\mathcal{ME}(\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{row}})$
with positive definite covariances
$\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}}$, and let
$(\hat{\bm{\operatorname{M}}},\hat{\bm{\Sigma}}^{\operatorname{row}},\hat{\bm{\Sigma}}^{\operatorname{col}})$
be the corresponding MMCD estimators. Then it holds that

$\|\hat{\bm{\operatorname{M}}}-\bm{\operatorname{M}}\|\xrightarrow{a.s.}0,\quad\left|c_{1}\hat{\bm{\Sigma}}^{\operatorname{col}}\otimes\hat{\bm{\Sigma}}^{\operatorname{row}}-\bm{\Sigma}^{\operatorname{col}}\otimes\bm{\Sigma}^{\operatorname{row}}\right|\xrightarrow{a.s.}0,$

where $c_{1}$ is a distribution-specific consistency factor as in Croux and
Haesbroeck, (1999).

The proof of the consistency of the _raw_ MMCD estimators relies on the strong
consistency of the MCD estimator given in Butler et al., (1993); Cator and
Lopuhaä, (2012) and is provided in Supplement B.1.

In the multivariate case, Lopuhaa and Rousseeuw, (1991) proposed a one-step
reweighting approach to improve the efficiency of the MCD estimator and showed
that the breakdown point of the original estimator is preserved. We generalize
this approach to the matrix-variate setting to obtain the _reweighted_ MMCD
estimator.

######  Definition 3.0.1.

Let $\bm{\mathfrak{X}}$ be a collection of $n$ i.i.d. samples from a
continuous $p\times q$ matrix-variate distribution, where
$d=\lfloor\nicefrac{{p}}{{q}}+\nicefrac{{q}}{{p}}\rfloor$,
$p,q\in\mathbb{N},p\geq 2,q\geq 2$, and let $\hat{\bm{\operatorname{M}}}$,
$\hat{\bm{\Sigma}}^{\operatorname{row}}$, and
$\hat{\bm{\Sigma}}^{\operatorname{col}}$ denote the _raw_ MMCD estimators as
in Definition 2.0.1. The _reweighted_ MMCD estimators are given by

$\displaystyle\tilde{\bm{\operatorname{M}}}$
$\displaystyle=\frac{1}{\sum_{i=1}^{n}w(d(\bm{X}_{i}))}\sum_{i=1}^{n}w(d(\bm{X}_{i}))\bm{X}_{i},$
(19) $\displaystyle\tilde{\bm{\Sigma}}^{\operatorname{row}}$
$\displaystyle=\frac{1}{q\sum_{i=1}^{n}w(d(\bm{X}_{i}))}\sum_{i=1}^{n}w(d(\bm{X}_{i}))(\bm{X}_{i}-\tilde{\bm{\operatorname{M}}})\tilde{\bm{\Omega}}^{\operatorname{col}}(\bm{X}_{i}-\tilde{\bm{\operatorname{M}}})^{\prime}\text{,
and}$ (20) $\displaystyle\tilde{\bm{\Sigma}}^{\operatorname{col}}$
$\displaystyle=\frac{1}{p\sum_{i=1}^{n}w(d(\bm{X}_{i}))}\sum_{i=1}^{n}w(d(\bm{X}_{i}))(\bm{X}_{i}-\tilde{\bm{\operatorname{M}}})^{\prime}\tilde{\bm{\Omega}}^{\operatorname{row}}(\bm{X}_{i}-\tilde{\bm{\operatorname{M}}}),$
(21)

where $w:[0,\infty)\rightarrow[0,\infty)$ is a non-increasing and bounded
weight function such that $w(d(\bm{X}_{i}))>0$ for at least
$\lfloor\nicefrac{{(n+d+2)}}{{2}}\rfloor$ observations that vanishes for large
distances, i.e., $w(d(\bm{X}_{i}))=0$ if $d(\bm{X}_{i})>c_{1}>0$.

######  Theorem 3.0.3.

Let $\bm{\mathfrak{X}}$ be a collection of $n$ i.i.d. samples from a
continuous $p\times q$ matrix-variate distribution, where
$d=\lfloor\nicefrac{{p}}{{q}}+\nicefrac{{q}}{{p}}\rfloor$,
$p,q\in\mathbb{N},p\geq 2,q\geq 2$, and let
$\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{X}}}$,
$\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{X}}}$, and
$\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{X}}}$ denote the _raw_
MMCD estimators as in Definition 2.0.1 with breakdown points

$\displaystyle\varepsilon^{\ast}(\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{X}}},\bm{\mathfrak{X}})=\varepsilon^{\ast}(\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{X}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{X}}},\bm{\mathfrak{X}})=\frac{1}{n}\lfloor\min(n-h+1,h-(d+1))\rfloor=:\frac{m}{n},$

and let $\tilde{\bm{\operatorname{M}}}_{\bm{\mathfrak{X}}}$,
$\tilde{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{X}}}$, and
$\tilde{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{X}}}$ denote the
_reweighted_ estimators as in Definition 3.0.1. Then,

$\displaystyle\varepsilon^{\ast}(\tilde{\bm{\operatorname{M}}}_{\bm{\mathfrak{X}}},\bm{\mathfrak{X}})\geq\frac{m}{n}\quad\text{and}\quad\varepsilon^{\ast}(\tilde{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{X}}},\tilde{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{X}}},\bm{\mathfrak{X}})\geq\frac{m}{n}.$

A proof is given in Supplement B.1. For the algorithm used to compute the
reweighted MMCD estimators introduced in the following section, we use the
weight function $w:[0,\infty)\mapsto\\{0,1\\}$ with

$\displaystyle
w(\operatorname{MMD}^{2}(\bm{X}_{i})):=\begin{cases}1&\text{if}\quad i\in
H\vee\operatorname{MMD}^{2}(\bm{X}_{i})<\chi^{2}_{pq;0.975}\\\
0&\text{otherwise}\end{cases}.$ (22)

Note that the $h$ observations in the $h$-subset of the raw MMCD estimator
have the lowest MMDs, and the condition that all observations $i\in H$ get a
positive weight ensures that the reweighting step does not lead to an
estimator that uses fewer than $h$ samples.

##  4 Algorithm

Rousseeuw and Driessen, (1999) proposed the Fast-MCD algorithm to efficiently
compute the MCD estimator. The key idea to find the $h$-subset with the lowest
covariance determinant is based on the concentration step (C-step): after each
C-step, the objective function is smaller or equal as before, and by
repeatedly applying C-steps convergence is reached within finitely many
iterations.

###  4.1 Adapting the C-step

Adapting the structure of the Fast-MCD algorithm to the matrix-variate setting
leads to the development of the Fast-MMCD algorithm. This adaptation
necessitates a modification in the covariance estimation during the C-step to
derive suitable counterparts for computing the MMCD estimators. However, this
process encounters a challenge due to the involvement of two covariance
matrices, as depicted in Equations (10) and (11), both lacking closed-form
solutions for their estimation. To address this issue, we incorporate the
flip-flop algorithm introduced by Dutilleul, (1999) within the C-step.
Consider a matrix-variate random sample
$\bm{\mathfrak{X}}=(\bm{X}_{1},\dots,\bm{X}_{n})$, with
$\bm{X}_{i}\in\mathbb{R}^{p\times q}$, and any $h$-subset
$H_{\operatorname{old}}\subset\\{1,\ldots,n\\}$, with
$\left|H_{\operatorname{old}}\right|=h>\lfloor\nicefrac{{p}}{{q}}+\nicefrac{{q}}{{p}}\rfloor+2$.
First, the MLEs
$(\hat{\bm{\operatorname{M}}}_{H_{\operatorname{old}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{H_{\operatorname{old}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{H_{\operatorname{old}}})$
are computed based on the observations in the subset $H_{\operatorname{old}}$
using the flip-flop algorithm, which is non-decreasing in likelihood. Next,
compute the squared Mahalanobis distances
$d^{2}_{i}(H_{\operatorname{old}}):=\operatorname{MMD}^{2}(\bm{X}_{i},\hat{\bm{\operatorname{M}}}_{H_{\operatorname{old}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{H_{\operatorname{old}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{H_{\operatorname{old}}})$
for all $i=1,\ldots,n$. In Proposition 2.0.1, we showed that $\sum_{i\in
H_{\operatorname{old}}}d^{2}_{i}(H_{\operatorname{old}})=hpq$, hence only the
terms of the log-likelihood function involving the determinants change in this
step. To construct the new subset $H_{\operatorname{new}}$, sort the squared
MMDs in ascending order, resulting in a permutation $\pi$ of
$\\{1,\ldots,n\\}$ such that
$d^{2}_{\pi(1)}(H_{\operatorname{old}})\leq\ldots\leq
d^{2}_{\pi(n)}(H_{\operatorname{old}})$, and define a new $h$-subset
$H_{\operatorname{new}}=\\{\pi(1),\ldots,\pi(h)\\}$. Since the estimators do
not change in this step, the terms involving the determinant are constant and
by construction, the sum of the Mahalanobis distances either decreases or
stays constant. Hence, the reordering is non-decreasing in likelihood.
Finally, the MLEs are updated using the flip-flop algorithm based on the
observations in the subset $H_{\operatorname{new}}$, resulting in estimators
$(\hat{\bm{\operatorname{M}}}_{H_{\operatorname{new}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{H_{\operatorname{new}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{H_{\operatorname{new}}})$,
increasing the likelihood once more, and it follows that

$p\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{col}}_{H_{\operatorname{new}}}))+q\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{row}}_{H_{\operatorname{new}}}))\leq
p\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{col}}_{H_{\operatorname{old}}}))+q\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{row}}_{H_{\operatorname{old}}})).$
(23)

By repeatedly applying such C-steps, we can decrease the covariance
determinant in subsequent iterations as in Equation (23). This results in a
decreasing and non-negative sequence of determinants that must converge after
exploring finitely many $h$-subsets. Similar to the multivariate case, we
obtain equality of the determinants from one $h$-subset to the next if and
only if the estimators do not change from one to the next iteration. However,
this does not necessarily imply that we have found a global optimum. A pseudo-
code for this matrix-variate version of the C-step is given in Algorithm 1 in
Supplement C.

###  4.2 The Fast-MMCD algorithm

The Fast-MMCD algorithm is a matrix-variate extension of the Fast-MCD
procedure of Rousseeuw and Driessen, (1999), aiming to alleviate the C-steps
dependence on the initial subset by using multiple initial subsets,
iteratively conducting C-steps on each until convergence, and ultimately
selecting the solution with the lowest determinant. While this explains the
idea of the algorithm, there are more computational considerations and
adjustments in the full Fast-MMCD. A pseudo-code of the FastMMCD Algorithm 2
is given in Supplement C.

As in its multivariate counterpart, the Fast-MMCD procedure uses so-called
_elemental_ subsets to initialize the procedure. This means that we use $m$
subsets of size $d+2,d=\lfloor\nicefrac{{p}}{{q}}+\nicefrac{{q}}{{p}}\rfloor$,
instead of size $h$, to increase the probability of obtaining at least one
clean initial subset. Using $m=500$ elemental subsets by default allows for a
reasonable tradeoff between a wide variety of settings where we likely obtain
at least one clean elemental subset and the computational demands of computing
initial estimators. If either $p\ll q$ or $q\gg p$, $d$ will be large, and
using more initial subsets is recommended. Using elemental subsets increases
not only the robustness of the initial estimators but also the computational
efficiency.

Moreover, the Fast-MMCD procedure only uses 2 C-step and MLE iterations for
the initial elemental subsets to ensure even faster computation of the initial
estimators. In the MLE procedure, Werner et al., (2008) demonstrated that the
same asymptotic efficiency can be attained using only two iterations instead
of iterating until convergence. As for the C-step, Rousseeuw and Driessen,
(1999) outlined that after two iterations, subsets with the lowest covariance
determinant during the procedure can already be identified, even before
reaching convergence. Moreover, simulations show that we can identify those
initial subsets that yield robust solutions after 2 C-steps, whether we use 2
MLE iterations or iterate the flip-flop algorithm until convergence. This is
described in detail in Supplement C.1, where we also show that elemental
subsets indeed yield more robust solutions than their larger counterparts in
case of high contamination.

The initialization step of the FastMMCD procedure yields $m$ initial
estimators and we keep the $10$ estimators with the lowest covariance
determinant. Using those as initial estimators, we iterate C-steps until
convergence on the full data set $\bm{\mathfrak{X}}$. The solution with the
lowest covariance determinant then yields the raw MMCD estimators.

The raw MMCD estimators are scaled using the consistency factor $c(\alpha)$
given in Equation (18) to achieve consistency at the normal model as outlined
in Theorem 3.0.2. Based on those rescaled raw MMCD estimators, the reweighted
estimators described in Definition 3.0.1 are computed using the weights given
in Equation (22). The reweighted MMCD estimators are then scaled using
$c(\tilde{\alpha})=c(\nicefrac{{\tilde{h}}}{{n}})$, where $\tilde{h}$ denotes
the number of observations with weights one. Simulations in Section 6 show
that the finite sample efficiency of the reweighed MMCD estimators increases
to over $90\%$ as the sample size increases.

The Fast-MMCD algorithm repeatedly computes Mahalanbois distances for all $n$
samples, which is computationally expensive when $n$ gets large. To improve
the computational efficiency for settings where $n$ is large, we implemented
the subsampling approach proposed by Rousseeuw and Driessen, (1999). The idea
is to split the sample of $n$ observations into several smaller subsamples and
compute initial estimators on those subsamples before working on the large set
with $n$ observations.

##  5 Outlier detection and explainability

Given a sample $\bm{\mathfrak{X}}=(\bm{X}_{1},\dots,\bm{X}_{n})$ of matrix-
variate observations, the task for outlier detection is to identify those
observations which are “far away” from the center of the data cloud with
respect to its shape. In robust statistics, it is common to consider the
Mahalanobis distance for this purpose, assume an underlying normal
distribution of the observations, and use a quantile of the chi-square
distribution as an outlier cutoff value (Maronna et al.,, 2019). Here, we
follow the same idea: an observation $\bm{X}_{i}$ is flagged as an outlier if

$\operatorname{MMD}^{2}(\bm{X}_{i};\hat{\bm{\operatorname{M}}},\hat{\bm{\Sigma}}^{\operatorname{row}},\hat{\bm{\Sigma}}^{\operatorname{col}})>\chi^{2}_{pq;0.975}\
,$

for $i\in\\{1,\ldots,n\\}$ and the MMCD estimators
$\hat{\bm{\operatorname{M}}}$, $\hat{\bm{\Sigma}}^{\operatorname{row}}$, and
$\hat{\bm{\Sigma}}^{\operatorname{col}}$.

Even though this information is valuable in practice, it is not very useful
for understanding the reasons for outlyingness. This is the goal of outlier
explainability, where the contributions of the cells/rows/columns of the
matrix-valued observations are investigated in more detail. We will use the
concept of Shapley values for this purpose and first briefly review how this
is applied to multivariate data before extending it to the matrix-variate
case. For details, we refer to Mayrhofer and Filzmoser, (2023).

###  5.1 Shapley values for multivariate data

Let $\bm{x}=(x_{1},\ldots,x_{p})^{\prime}$ denote an observation vector from a
population with expectation vector
$\bm{\mu}=(\mu_{1},\ldots,\mu_{p})^{\prime}$ and covariance matrix
$\bm{\Sigma}$, and $P=\\{1,\ldots,p\\}$ the index set of the variables. Then
the outlyingness contributions
$\bm{\phi}(\bm{x},\bm{\mu},\bm{\Sigma})=\bm{\phi}(\bm{x})=(\phi_{1}(\bm{x}),\ldots,\phi_{p}(\bm{x}))$
based on the Shapley value assign each variable its average marginal
contribution to the squared Mahalanobis distance, i.e.,

$\displaystyle\phi_{k}(\bm{x},\bm{\mu},\bm{\Sigma})=\sum_{S\subseteq
P\setminus\\{k\\}}\frac{\left|S\right|!(p-\left|S\right|-1)!}{p!}\Delta_{k}\operatorname{MD}^{2}(\hat{\bm{x}}^{S})=(x_{k}-\mu_{k})\sum_{j=1}^{p}(x_{j}-\mu_{j})\omega_{jk},$
(24)

with marginal contributions

$\Delta_{k}\operatorname{MD}^{2}(\hat{\bm{x}}^{S}):=\operatorname{MD}^{2}(\hat{\bm{x}}^{S\cup\\{k\\}})-\operatorname{MD}^{2}(\hat{\bm{x}}^{S})\quad\text{and}\quad\hat{x}_{j}^{S}:=\begin{cases}x_{j}&\text{if
}j\in S\\\ \mu_{j}&\text{if }j\notin S\end{cases}$ (25)

as the components of $\hat{\bm{x}}^{S}$. Here, $\omega_{jk}$ is the element
$(j,k)$ of $\bm{\Omega}=\bm{\Sigma}^{-1}$. Since $\bm{\phi}(\bm{x})$ is based
on the Shapley value, it is the only decomposition of the squared Mahalanobis
distance based on Equation (25) that fulfills the following properties:

  * •

_Efficiency:_ The contributions $\phi_{j}(\bm{x})$, for $j=1,\ldots,p$, sum up
to the squared Mahalanobis distance of $\bm{x}$, hence
$\sum_{j=1}^{p}\phi_{j}(\bm{x})=\operatorname{MD}^{2}(\bm{x})$.

  * •

_Symmetry:_ If
$\operatorname{MD}^{2}(\hat{\bm{x}}^{S\cup\\{j\\}})=\operatorname{MD}^{2}(\hat{\bm{x}}^{S\cup\\{k\\}})$
holds for all subsets $S\subseteq P\setminus\\{j,k\\}$ for two coordinates $j$
and $k$, then $\phi_{j}(\bm{x})=\phi_{k}(\bm{x})$.

  * •

_Monotonicity:_ Let $\bm{\mu},\tilde{\bm{\mu}}\in\mathbb{R}^{p}$ be two
vectors and $\bm{\Sigma},\tilde{\bm{\Sigma}}\in\operatorname{PSD}(p)$ be two
matrices. If

$\displaystyle\operatorname{MD}_{\bm{\mu},\bm{\Sigma}}^{2}(\hat{\bm{x}}^{S\cup\\{j\\}})-\operatorname{MD}_{\bm{\mu},\bm{\Sigma}}^{2}(\hat{\bm{x}}^{S})\geq\operatorname{MD}_{\tilde{\bm{\mu}},\tilde{\bm{\Sigma}}}^{2}(\hat{\bm{x}}^{S\cup\\{j\\}})-\operatorname{MD}_{\tilde{\bm{\mu}},\tilde{\bm{\Sigma}}}^{2}(\hat{\bm{x}}^{S})$

holds for all subsets $S\subseteq P$, then
$\phi_{j}(\bm{x},\bm{\mu},\bm{\Sigma})\geq\phi_{j}(\bm{x},\tilde{\bm{\mu}},\tilde{\bm{\Sigma}})$.

In words, the coordinate $\phi_{k}(x)$ of the Shapley value is the average
marginal contribution of the $k$-th variable to the squared Mahlanobis
distance and is obtained by averaging over all marginal outlyingness
contributions $\Delta_{k}\operatorname{MD}^{2}(\hat{\bm{x}}^{S})$ across all
possible subsets $S\subseteq P\setminus\\{k\\}$. Although this suggests an
exponential computational complexity, which becomes costly especially if $p$
is large, the second equality in Equation (24) reveals just linear complexity;
for a proof we refer to Mayrhofer and Filzmoser, (2023). This allows for
another insight into the Shapley value by comparing it to the squared
Mahalanobis distance, which can be written as
$\sum_{j,k=1}^{p}(x_{j}-\mu_{j})(x_{k}-\mu_{k})\omega_{jk}$. While the latter
calculates an outlyingness measure by aggregating the contributions
$(x_{j}-\mu_{j})(x_{k}-\mu_{k})\omega_{jk}$ of all variables for the entire
observation, Equation (24) shows that a coordinate $\phi_{k}(x)$ of the
Shapley value only considers the contributions that are associated with the
$k$-th variable.

###  5.2 Shapley value for matrix-valued data

To define Shapley values for matrix-variate data, we can use the connection
between the matrix and multivariate Mahalanobis distance; see Equation (6).
Let $\bm{X}\in\mathbb{R}^{p\times q}$ be a matrix-variate sample with mean
$\bm{\operatorname{M}}\in\mathbb{R}^{p\times q}$ and covariance matrices
$\bm{\Sigma}^{\operatorname{row}}\in\operatorname{PSD}(p)$ and
$\bm{\Sigma}^{\operatorname{col}}\in\operatorname{PSD}(q)$. The
$pq$-dimensional vectorized observation is denoted as
$\bm{x}=\operatorname{vec}(\bm{X})$, with mean
$\bm{\mu}=\operatorname{vec}(\bm{\operatorname{M}})$ and covariance matrix
$\bm{\Sigma}=\bm{\Sigma}^{\operatorname{col}}\otimes\bm{\Sigma}^{\operatorname{row}}$.
Based on Equation (24), we can obtain outlyingness contributions for every
coordinate of $\bm{x}$ and hence for every cell of the matrix $\bm{X}$ by

$\displaystyle\phi_{a}(\bm{x})$
$\displaystyle=(x_{a}-\mu_{a})\sum_{b=1}^{pq}(x_{b}-\mu_{b})\omega_{ab}=(x_{jk}-m_{jk})\sum_{i=1}^{p}\sum_{l=1}^{q}(x_{il}-m_{il})\omega_{ij}^{\operatorname{row}}\omega_{kl}^{\operatorname{col}}=\phi_{jk}(\bm{X}),$

with $a=i+(l-1)p$ and $b=j+(k-1)p$. Analyzing the outlyingness scores for all
$pq$ cells of an observation $\bm{X}\in\mathbb{R}^{p\times q}$ might be
unpractical, especially when we want to compare multiple observations. Instead
of the cellwise scores, we can consider row- and columnwise outlyingness
scores for matrix-variate data. To obtain those scores, we could either sum up
the cellwise contributions $\phi_{jk}(\bm{X})$, with $j=1,\ldots,p$ and
$k=1,\ldots,q$, or replace an entire row $\bm{x}_{i.}$ of $\bm{X}$ with its
corresponding mean, i.e., the $i$-th row of $\bm{\operatorname{M}}$. For the
latter setting, we must adapt Equation (25) to the matrix case. The rowwise
marginal contributions to the matrix Mahalanobis distance are defined as as

$\displaystyle\Delta_{j}\operatorname{MMD}(\hat{\bm{X}}^{S})$
$\displaystyle:=\operatorname{MMD}(\hat{\bm{X}}^{S\cup\\{j\\}})-\operatorname{MMD}(\hat{\bm{X}}^{S}),$

where the $i$-th row of $\hat{\bm{X}}^{S}$ is given as $(x_{i1},\dots,x_{iq})$
if $i\in S$ and $(m_{i1},\dots,m_{iq})$ if $i\notin S$.

######  Proposition 5.2.1.

The $j$-th coordinate of the rowwise Shapley value is given by

$\displaystyle\begin{split}\phi_{j.}(\bm{X})&:=\sum_{S\subseteq
P\setminus\\{j\\}}w(\left|S\right|)\Delta_{j}\operatorname{MMD}(\hat{\bm{X}}^{S})=\sum_{i=1}^{p}\sum_{k=1}^{q}\sum_{l=1}^{q}(x_{jl}-m_{jl})(x_{ik}-m_{ik})\omega_{ij}^{\operatorname{row}}\omega_{kl}^{\operatorname{col}},\end{split}$
(26)

###### Proof.

A proof for Equation (26) can be found in Supplement D. ∎

From the above, we conclude that adding up the contributions of all cells in a
row yields identical scores as when we adapt the marginal contributions to a
rowwise replacement. Using the first approach, we can sum up the columns and
obtain columnwise outlying contributions. We can achieve the same by
transposing the matrix in the latter case. The cellwise, rowwise, and
columnwise Shapley values of all features can be computed efficiently using
matrix operations as listed in Equations (27), (28), and (29), respectively.
In Equation 27, the symbol $\circ$ refers to element-wise multiplication.

$\displaystyle\bm{\Phi}(\bm{X})$
$\displaystyle=(\bm{X}-\bm{\operatorname{M}})\circ\bm{\Omega}^{\operatorname{row}}(\bm{X}-\bm{\operatorname{M}})\bm{\Omega}^{\operatorname{col}}\in\mathbb{R}^{p\times
q}$ (27) $\displaystyle\bm{\phi}_{\operatorname{row}}(\bm{X})$
$\displaystyle=\operatorname{diag}(\bm{\Omega}^{\operatorname{row}}(\bm{X}-\bm{\operatorname{M}})\bm{\Omega}^{\operatorname{col}}(\bm{X}-\bm{\operatorname{M}})^{\prime})\in\mathbb{R}^{p}$
(28) $\displaystyle\bm{\phi}_{\operatorname{col}}(\bm{X})$
$\displaystyle=\operatorname{diag}((\bm{X}-\bm{\operatorname{M}})^{\prime}\bm{\Omega}^{\operatorname{row}}(\bm{X}-\bm{\operatorname{M}})\bm{\Omega}^{\operatorname{col}})\in\mathbb{R}^{q}$
(29)

##  6 Simulations

In the simulation studies outlined in this section, our primary focus is to
rigorously assess the performance of the MMCD estimators. We aim to validate
their demonstrated theoretical properties and compare efficiency against ML
estimators. Despite our initial intention to include various robust estimators
as mentioned in Section 1, practical constraints arose as the relevant
routines were exclusively accessible in Matlab, while our current framework
operates within R. For the sake of consistency and practical implementation,
we concentrate on comparing the efficiency of the raw and reweighted MMCD,
alongside an in-depth analysis of MLE and (reweighted) MMCD on contaminated
data. To ensure the highest possible breakdown point across all simulations
and examples discussed in this paper, we set
$h=\lfloor\nicefrac{{(n+d+2)}}{{2}}\rfloor$ for the MMCD estimators.

We generate samples from a matrix normal distribution with
$(p,q)\in\\{(5,20),(50,20)$, $(100,50)\\}$. The considered sample sizes are
$n\in\\{20,50,100,300,1000\\}$. The mean matrix has all zero entries. For the
rowwise covariance matrix we use randomly generated covariance matrices
proposed by Agostinelli et al., (2015),
$\bm{\Sigma}^{\operatorname{row}}=\bm{\Sigma}^{\operatorname{rnd}}\in\operatorname{PSD}(p)$,
generally yielding low correlations. For the columnwise covariance, we have
$\bm{\Sigma}^{\operatorname{col}}=\bm{\Sigma}^{\operatorname{mix}}(0.7)\in\operatorname{PSD}(q)$,
with entries $\sigma_{jk}^{\operatorname{mix}}(0.7)=0.7^{\left|j-k\right|}$.
Each simulation setting is replicated 100 times.

Using solely the clean samples, we assess the normal finite-sample efficiency
by comparing the ratio

$\frac{D(\hat{\bm{\Sigma}}^{\operatorname{row}}_{\operatorname{MLE}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{\operatorname{MLE}})}{D(\hat{\bm{\Sigma}}^{\operatorname{row}}_{\operatorname{MMCD}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{\operatorname{MMCD}})},$

where
$D(\hat{\bm{\Sigma}}^{\operatorname{row}},\hat{\bm{\Sigma}}^{\operatorname{col}})$
denotes the Kullback-Leiber (KL) divergence of the estimators
$\hat{\bm{\Sigma}}^{\operatorname{row}}$ and
$\hat{\bm{\Sigma}}^{\operatorname{col}}$ in the matrix normal setting
$\mathcal{MN}(\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}})$,
and is given by

$\displaystyle\begin{split}D(\hat{\bm{\Sigma}}^{\operatorname{row}},\hat{\bm{\Sigma}}^{\operatorname{col}})=&\operatorname{tr}(\bm{\Omega}^{\operatorname{row}}\hat{\bm{\Sigma}}^{\operatorname{row}})\operatorname{tr}(\bm{\Omega}^{\operatorname{col}}\hat{\bm{\Sigma}}^{\operatorname{col}})\\\
&-q\log(\det(\bm{\Omega}^{\operatorname{row}}\hat{\bm{\Sigma}}^{\operatorname{row}}))-p\log(\det(\bm{\Omega}^{\operatorname{col}}\hat{\bm{\Sigma}}^{\operatorname{col}}))-pq,\end{split}$
(30)

with
$\bm{\Omega}^{\operatorname{row}}=(\bm{\Sigma}^{\operatorname{row}})^{-1}$ and
$\bm{\Omega}^{\operatorname{col}}=(\bm{\Sigma}^{\operatorname{col}})^{-1}$.
Figure 2, shows that the efficiency of the raw MMCD estimators is below $0.5$
on average for both columnwise covariance matrices. In contrast, the
reweighted estimators’ efficiency is above $0.5$ for $n=100$ and it rises to
over $0.9$ as the sample size increases.

Figure 2: Comparison of the finite-sample efficiency of raw and reweighted
MMCD.

For the setting with contamination, we generate two scenarios where 20% of the
observations are outliers, with a mean matrix of either values of one or
three, and with the same covariance matrices as for the regular observations.
As evaluation measures, we use KL divergence (30) to analyze the quality of
the covariance estimation. Additionally, we analyze outlier detection
capabilities of the squared Mahalanobis distance based on MMCD and MLE
covariance estimators, with the $\chi^{2}_{pq,0.99}$ quantile as a detection
threshold. The performance is measured by the F-score, the harmonic mean of
precision and recall. The results are presented in Figure 3, where the KL
divergence is plotted on a log scale to enable a comparison between the robust
MMCD estimators and the ML estimators.

Figure 3: Comparison between F-Score and KL divergence of MMCD and ML
estimators for simulated data from a matrix normal distribution with $20\%$
contamination.

In the upper row of Figure 3, the boxplots depict F-scores across various
parameter configurations. Notably, the MMCD estimators exhibit improved
performance as sample sizes increase across all settings, consistently
outperforming ML estimators. However, for $(p,q)=(5,20)$, in a scenario
involving a minor mean shift, the F-scores derived from MMCD exhibit some
volatility with larger sample sizes. This situation arises due to the
proximity of outliers to regular observations, posing challenges in their
identification. Notably, a more pronounced mean shift significantly simplifies
outlier detection. The KL divergence metric distinctly favors the robust
estimators over classical ones. Yet, in the case of $(p,q)=(5,20)$ with a
small mean shift and limited sample sizes, the median results between MMCD and
ML estimators appear closer, albeit MMCD demonstrates significantly lower
variability.

In Supplement E, we expand our analysis to cover additional simulation
scenarios, incorporating diverse performance metrics such as precision,
recall, Frobenius error, and the angle between eigenvalues of the covariance
matrices.

##  7 Examples

###  7.1 Glacier weather data – Sonnblick observatory

We analyze the publicly available weather data from Austria’s highest weather
station, located in the Austrian Central Alps at an elevation of 3106 m above
sea level on top of the glaciated mountain “Hoher Sonnblick” (datasource:
GeoSphere Austria - https://data.hub.geosphere.at). The observed parameters
are monthly averages of temperature (T), precipitation (P), proportion of
solid precipitation (SP), air pressure (AP), and sunshine hours (SH). We
consider the monthly values between 1891 and 2022 and exclude five years with
missing values, yielding $n=127$ observations of $p=5$ times $q=12$
dimensional matrices. Our goal is to identify observations that show a
different weather pattern than the majority of the data and explain why the
corresponding years deviate from the majority. We did not adjust for a
possible yearly trend in this exploratory analysis as we wish to understand
long-term patterns and shifts in climate without the influence of adjustments.

Figure 4: Yearly outlyingness contributions for the glacier weather data.
Regular years are white, and years that contain missing data are gray.
Outliers are colored as follows: blue for “above average”, red for “below
average”, and color intensity proportional to rowwise Shapley value.

In total, outlier detection based on the MMCD estimators flags 23 outlying
matrices, which are indicated in Figure 4 as colored years: If the aggregated
monthly measurements are above their average, the cells are colored red;
otherwise, they are colored blue. The rowwise Shapley value is then used to
determine color brightness, i.e., the larger the outlyingness contribution,
the darker the color. Years with missing observations are grey; years with
only white cells refer to regular observations. It is visible that the outlier
frequency increases in the last period. Moreover, more recent outliers are
characterized by increased temperature, precipitation, air pressure, and a
lack of solid precipitation (e.g. snow) – a clear signal of a climate change.

Figure 5: Outlyingess contributions based on cellwise Shapley values for the
years 1895 and 2022 of the glacier weather data using the same color scheme as
in Figure 4.

In Figure 5, we use cellwise Shapley values to understand which parameters in
which months contributed most to the outlyingness of 1895 and 2022,
corresponding to the first and last outlying observation in the data set,
respectively, where the color scheme is inherited from Figure 4. The largest
outlyingness contribution is due to an unusually large amount of precipitation
in March 1895. Overall, high amounts of precipitation were observed that year,
with a high percentage of snow even in the summer months. In contrast, the
largest outlyingness contributions in 2022 are due to a very sunny March and
low percentages of snowfall in May, June, and August.

###  7.2 Darwin data

We consider the DARWIN (Diagnosis AlzheimeR WIth haNdwriting) (Cilia et al.,,
2022) data set containing handwriting samples of 174 subjects, 89 diagnosed
with Alzheimer’s disease (AD), and 85 healthy subjects (H). Each individual
completed 25 handwriting tasks on paper, and the pen movements were recorded
using a graphic tablet. The tasks are ordered in difficulty. From the raw
handwriting data, 18 features were extracted: Total Time, Air Time, Paper
Time, Mean Speed on paper, Mean Speed in air, Mean Acceleration on paper, Mean
Acceleration in air, Mean Jerk on paper, Mean Jerk in air, Pressure Mean,
Pressure Variance, Generalization of the Mean Relative Tremor (GMRT) on paper,
GMTR in air, Mean GMRT, Pendowns Number, Max X Extension, Max Y Extension, and
Dispersion Index. For a more detailed description of the data, we refer to
Cilia et al., (2018). In Cilia et al., (2022), each task was considered
separately to train a classifier, and the combination of the classifiers led
to an improvement in the classification of subjects. Our focus here lies not
in the classification task but rather in explaining the differences between AD
and H groups. We treat the observations as matrices, with the rows
representing the extracted features and the columns representing the tasks.
Because of linear dependencies, the variables Total Time and Mean GMRT were
excluded. Further, the variable Air Time had several extreme and unreliable
measurements and was thus also excluded. This yields observation matrices with
$p=15$ features and $q=25$ tasks.

We applied the MMCD procedure only on the healthy subjects and used the robust
estimators to compute MMDs for all observations. Thus, the MMDs presented in
Figure 6 left are generally smaller for the H group, whereas all observations
from the AD group exceed the outlier cutoff value. The fact that healthy
subjects also exceed the cutoff value shows the heterogeneity in this group.

Figure 6: Plot of robust MMD based on MMCD estimators for the Darwin data on
the left, and average proportional rowwise Shapley values for the H and AD
subjects on the right.

In the right panel of Figure 6, we consider the _average proportional
contributions_ of the variables to the MMDs for the H and AD groups. The
outlyingness contributions are based on the rowwise Shapley values, resulting
in 15 scores for each individual. Since those scores sum up to the squared
MMD, we can divide them by the squared MMD to get proportional contributions,
and by averaging over all individuals in the H and AD groups, respectively, we
obtain the values shown in this plot. Large differences between the AD and H
groups indicate variables that are important to distinguish between healthy
individuals and those who have Alzheimer’s disease. For example, Pressure Mean
and Paper Time are evidently higher in the AD group.

###  7.3 Video data

In this example, we examine a surveillance video of a beach sourced from Li et
al., (2004). The video comprises 633 frames, each sized at $128\times 160$
pixels; five selected frames are shown in Figure 7. The majority of the frames
depict the beach scene. Around frame 500, a man walks into the scene from the
left and partly disappears behind the tree. As he continues walking, he
reappears on the right side of the tree and remains in the video until the
end.

Figure 7: Selected frames of the video data.

For our analyis, we converted the original RGB video to a grayscale video,
applied the Fast-MMCD procedure, and obtained MMDs for all 633 frames, which
are visualized in Figure 8. The plot on the left shows the robust MMDs for all
633 frames, and the one on the right for frames 471 to 633 to better highlight
the increase in MMD when the man enters the scenery, with a short drop in MMD
when he disappears behind the tree. We indicate frames 487, 491, and 495, also
presented in Figure 9 in terms of their cellwise Shapley values. We see that
the pixels that form the contours of the man and most of the pixels of the
man’s head contribute most to the outlyingness. When the man disappears behind
the tree, there are fewer pixels with high outlyingness contributions. Since
the sum of the contributions amounts to the squared MMD of an observation,
this explains the behavior of the MMDs of the frames shown in Figure 8(b).

(a) Frames 1 to 633.

(b) Frames 471 to 633.

Figure 8: Plot of robust MMD based on MMCD estimators for the video data.

It is interesting to see a certain increase in the MMD in Figure 8(a) between
frames 400 and 450. Here, the Shapley values on the contour of the palm tree
contribute the most to the outlyingness. This could be caused by a slight
shifting of the camera or a small movement of the palm tree due to wind.

Figure 9: Outlyingness scores based on cellwise Shapley values are shown in
red, where darker colors indicate higher outlyingness contributions, and the
grayscale video frames are displayed in the background.

##  8 Summary and conclusions

Matrix-valued observations, like images or dual-factor data tables, are common
in various fields. To apply multivariate methods on matrix-valued data, the
matrices are typically converted to vectors by stacking either the rows or
columns. This disrupts the inherent data structure and increases
dimensionality, thereby complicating parameter estimation. Thus, it is often
preferable to model matrix-valued data directly with matrix-variate
distributions. In this setting, Maximum Likelihood (ML) estimation methods
exist for estimating the mean, as well as the row and column covariances,
respectively. However, these estimators are sensitive to deviations caused by
outliers among matrix-valued observations.

This work introduced the MMCD (matrix minimum covariance determinant)
estimators as a robust counterpart to the ML estimators in the matrix-variate
normal model. Several desirable properties are achieved: equivariance under
linear matrix transformations, high breakdown point, and consistency under
elliptical matrix-variate distributions. The proposed reweighted versions lead
to higher efficiency but not to any loss in terms of breakdown point. An
algorithm along the lines of the FastMCD procedure (Rousseeuw and Driessen,,
1999) allows for efficient computation of the estimators. Simulation
experiments validate the theoretical properties and advantages over robust
estimation for vectorized observations regarding breakdown and computational
efficiency.

We further extended the outlier explanation concept based on Shapley values
(Mayrhofer and Filzmoser,, 2023) to the matrix-variate setting. This allows
for an additive decomposition of the matrix-variate Mahalanobis distance of an
observation into Shapley contributions of either the rows, the columns, or the
matrix cells. The resulting Shapley values greatly aid with diagnostics,
particularly in revealing those cells (rows, columns) of the matrix with the
most substantial contributions to the outlyingness of the observation.

The efficiency of MMCD estimators in outlier detection for large sample sizes
is evident from the simulations. However, our future research aims to improve
and extend these estimators. For instance, smaller sample sizes might benefit
from integrating finite sample corrections proposed by Pison et al., (2002) to
enhance the results. Furthermore, the iterative computation of MMCD covariance
estimators, which involves inverse covariance matrices, requires data that
ensures full-rank estimates at each iteration. This requirement may be impeded
for example in image data, in case certain rows or columns maintain constant
pixel values across all observations. To solve this, regularization involving
a linear combination of the covariance matrix with a full-rank target matrix
can be used (Ledoit and Wolf,, 2004), similarly to the multivariate setting
(Boudt et al.,, 2020).

The MMCD objective can be expressed as a trimmed maximum likelihood problem,
and thus, can be extended to tensor-valued data using ML estimation for the
tensor normal distribution Manceur and Dutilleul, (2013). The framework of
Raymaekers and Rousseeuw, (2023) can be used to develop a cellwise robust
version of the MMCD. Our ongoing research focuses on extending the MMCD
estimators and outlier explanations based on Shapley values to the field of
functional data analysis. Our goal is to introduce robust estimators and
enhance interpretability for multivariate functional data. In the future, we
also plan to incorporate these robust estimators as plug-in estimators to
robustify established multivariate methodologies in the matrix-variate domain,
like principal component analysis and discriminant analysis.

Disclosure statement
The authors report there are no competing interests to declare.

Software and data availability
The R package robustmatrix available on CRAN includes a parallelized C++
implementation of the FastMMCD algorithm and a vignette to reproduce the
examples presented in this paper.

## References

  * Agostinelli et al., (2015) Agostinelli, C., Leung, A., Yohai, V. J., and Zamar, R. H. (2015).  Robust estimation of multivariate location and scatter in the presence of cellwise and casewise contamination.  Test, 24:441–461.
  * Boudt et al., (2020) Boudt, K., Rousseeuw, P. J., Vanduffel, S., and Verdonck, T. (2020).  The minimum regularized covariance determinant estimator.  Statistics and Computing, 30(1):113–128.
  * Butler et al., (1993) Butler, R., Davies, P., and Jhun, M. (1993).  Asymptotics for the minimum covariance determinant estimator.  The Annals of Statistics, pages 1385–1400.
  * Cator and Lopuhaä, (2012) Cator, E. A. and Lopuhaä, H. P. (2012).  Central limit theorem and influence function for the MCD estimators at general multivariate distributions.  Bernoulli, 18(2):520 – 551.
  * Cilia et al., (2022) Cilia, N. D., De Gregorio, G., De Stefano, C., Fontanella, F., Marcelli, A., and Parziale, A. (2022).  Diagnosing alzheimer’s disease from on-line handwriting: a novel dataset and performance benchmarking.  Engineering Applications of Artificial Intelligence, 111:104822.
  * Cilia et al., (2018) Cilia, N. D., De Stefano, C., Fontanella, F., and Di Freca, A. S. (2018).  An experimental protocol to support cognitive impairment diagnosis by using handwriting analysis.  Procedia Computer Science, 141:466–471.
  * Croux and Haesbroeck, (1999) Croux, C. and Haesbroeck, G. (1999).  Influence function and efficiency of the minimum covariance determinant scatter matrix estimator.  Journal of Multivariate Analysis, 71(2):161–190.
  * Davies, (1987) Davies, P. L. (1987).  Asymptotic behaviour of s-estimates of multivariate location parameters and dispersion matrices.  The Annals of Statistics, pages 1269–1292.
  * Dawid, (1981) Dawid, A. P. (1981).  Some matrix-variate distribution theory: notational considerations and a bayesian application.  Biometrika, 68(1):265–274.
  * Dutilleul, (1999) Dutilleul, P. (1999).  The mle algorithm for the matrix normal distribution.  Journal of Statistical Computation and Simulation, 64(2):105–123.
  * García-Escudero et al., (2010) García-Escudero, L. A., Gordaliza, A., Matrán, C., and Mayo-Iscar, A. (2010).  A review of robust clustering methods.  Advances in Data Analysis and Classification, 4:89–109.
  * Gupta and Nagar, (1999) Gupta, A. and Nagar, D. (1999).  Matrix Variate Distributions.  Monographs and Surveys in Pure and Applied Mathematics. Taylor & Francis.
  * Gupta and Varga, (2012) Gupta, A. K. and Varga, T. (2012).  Elliptically contoured models in statistics, volume 240.  Springer Science & Business Media.
  * Kurnaz et al., (2018) Kurnaz, F. S., Hoffmann, I., and Filzmoser, P. (2018).  Robust and sparse estimation methods for high-dimensional linear and logistic regression.  Chemometrics and Intelligent Laboratory Systems, 172:211–222.
  * Ledoit and Wolf, (2004) Ledoit, O. and Wolf, M. (2004).  A well-conditioned estimator for large-dimensional covariance matrices.  Journal of Multivariate Analysis, 88(2):365–411.
  * Li et al., (2004) Li, L., Huang, W., Gu, I. Y.-H., and Tian, Q. (2004).  Statistical modeling of complex backgrounds for foreground object detection.  IEEE Transactions on Image Processing, 13(11):1459–1472.
  * Lopuhaa and Rousseeuw, (1991) Lopuhaa, H. P. and Rousseeuw, P. J. (1991).  Breakdown points of affine equivariant estimators of multivariate location and covariance matrices.  The Annals of Statistics, pages 229–248.
  * Lu and Zimmerman, (2005) Lu, N. and Zimmerman, D. L. (2005).  The likelihood ratio test for a separable covariance matrix.  Statistics & probability letters, 73(4):449–457.
  * Lundberg and Lee, (2017) Lundberg, S. M. and Lee, S.-I. (2017).  A unified approach to interpreting model predictions.  In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R., editors, Advances in Neural Information Processing Systems 30, pages 4765–4774. Curran Associates, Inc.
  * Mahalanobis, (1936) Mahalanobis, P. C. (1936).  On the generalized distance in statistics.  Proceedings of the National Institute of Sciences (Calcutta), 2:49–55.
  * Manceur and Dutilleul, (2013) Manceur, A. M. and Dutilleul, P. (2013).  Maximum likelihood estimation for the tensor normal distribution: Algorithm, minimum sample size, and empirical bias and dispersion.  Journal of Computational and Applied Mathematics, 239:37–49.
  * Maronna et al., (2019) Maronna, R. A., Martin, R. D., Yohai, V. J., and Salibián-Barrera, M. (2019).  Robust statistics: Theory and Methods (with R).  John Wiley & Sons.
  * Mayrhofer and Filzmoser, (2023) Mayrhofer, M. and Filzmoser, P. (2023).  Multivariate outlier explanations using shapley values and mahalanobis distances.  Econometrics and Statistics.
  * Neykov et al., (2007) Neykov, N., Filzmoser, P., Dimova, R., and Neytchev, P. (2007).  Robust fitting of mixtures using the trimmed likelihood estimator.  Computational Statistics & Data Analysis, 52(1):299–308.
  * Pison et al., (2002) Pison, G., Van Aelst, S., and Willems, G. (2002).  Small sample corrections for lts and mcd.  Metrika, 55:111–123.
  * Raymaekers and Rousseeuw, (2023) Raymaekers, J. and Rousseeuw, P. J. (2023).  The cellwise minimum covariance determinant estimator.  Journal of the American Statistical Association, 0(0):1–12.
  * Roś et al., (2016) Roś, B., Bijma, F., de Munck, J. C., and de Gunst, M. C. (2016).  Existence and uniqueness of the maximum likelihood estimator for models with a kronecker product covariance structure.  Journal of Multivariate Analysis, 143:345–361.
  * Rousseeuw, (1985) Rousseeuw, P. (1985).  Multivariate estimation with high breakdown point.  Mathematical Statistics and Applications Vol. B, pages 283–297.
  * Rousseeuw and Driessen, (1999) Rousseeuw, P. J. and Driessen, K. V. (1999).  A fast algorithm for the minimum covariance determinant estimator.  Technometrics, 41(3):212–223.
  * Shapley, (1953) Shapley, L. S. (1953).  A value for n-person games.  Contributions to the Theory of Games, 2(28):307–317.
  * Soloveychik and Trushin, (2016) Soloveychik, I. and Trushin, D. (2016).  Gaussian and robust kronecker product covariance estimation: Existence and uniqueness.  Journal of Multivariate Analysis, 149:92–113.
  * Srivastava et al., (2008) Srivastava, M. S., von Rosen, T., and Von Rosen, D. (2008).  Models with a kronecker product covariance structure: estimation and testing.  Mathematical Methods of Statistics, 17:357–370.
  * Sun et al., (2016) Sun, Y., Babu, P., and Palomar, D. P. (2016).  Robust estimation of structured covariance matrix for heavy-tailed elliptical distributions.  IEEE Transactions on Signal Processing, 64(14):3576–3590.
  * Thompson et al., (2020) Thompson, G. Z., Maitra, R., Meeker, W. Q., and Bastawros, A. F. (2020).  Classification with the matrix-variate-t distribution.  Journal of Computational and Graphical Statistics, 29(3):668–674.
  * Tyler, (1987) Tyler, D. E. (1987).  A distribution-free m-estimator of multivariate scatter.  The Annals of Statistics, pages 234–251.
  * Werner et al., (2008) Werner, K., Jansson, M., and Stoica, P. (2008).  On estimation of covariance matrices with Kronecker product structure.  IEEE Transactions on Signal Processing, 56(2):478–491.
  * Zhang et al., (2022) Zhang, Y., Shen, W., and Kong, D. (2022).  Covariance estimation for matrix-valued data.  Journal of the American Statistical Association, pages 1–12.

Supplementary material to:

Robust covariance estimation and explainable outlier detection for matrix-
valued data

Marcus Mayrhofer, Una Radojicic, and Peter Filzmoser

##  Appendix A Preliminaries

Consider an i.i.d. sample
$\bm{\mathfrak{X}}=(\bm{X}_{1},\ldots,\bm{X}_{n})\in\mathbb{R}^{n\times
p\times q}$, with
$\bm{X}_{i}\sim\mathcal{MN}(\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}})$.
Due to the factored covariance structure of matrix normal data, the rowwise
and columnwise covariance matrices $\bm{\Sigma}^{\operatorname{row}}$ and
$\bm{\Sigma}^{\operatorname{col}}$ are only identified up to a multiplicative
constant $\kappa\neq 0$, since replacing $\bm{\Sigma}^{\operatorname{row}}$ by
$\kappa\bm{\Sigma}^{\operatorname{row}}$ and
$\bm{\Sigma}^{\operatorname{col}}$ by
$\nicefrac{{1}}{{\kappa}}\bm{\Sigma}^{\operatorname{col}}$ does not change the
pdf of $\bm{X}$. While the Kronecker product
$\bm{\Sigma}^{\operatorname{col}}\otimes\bm{\Sigma}^{\operatorname{row}}$ can
be uniquely identified, the issue of trivial non-uniqueness of
$\bm{\Sigma}^{\operatorname{row}}$ and $\bm{\Sigma}^{\operatorname{col}}$ is
commonly solved by either fixing a diagonal entry, the determinant, or the
norm of either matrix (Roś et al.,, 2016; Soloveychik and Trushin,, 2016). For
simplicity, we assume that the first diagonal entry of
$\bm{\Sigma}^{\operatorname{col}}$ is set to one. This implies that the
uniqueness of
$\bm{\Sigma}^{\operatorname{col}}\otimes\bm{\Sigma}^{\operatorname{row}}$ is
equivalent to the uniqueness of $\bm{\Sigma}^{\operatorname{col}}$ and
$\bm{\Sigma}^{\operatorname{row}}$ with the identifiability constraint
$\sigma^{\operatorname{col}}_{11}=1$. The multiplicative constant for their
estimators is also chosen such that
$\hat{\sigma}^{\operatorname{col}}_{11}=1$.

Instead of using Equations (2)-(4) for mean and covariance estimation, it is
also possible to consider the vectorized samples
$\bm{x}_{i}=\operatorname{vec}(\bm{X}_{i})\sim\mathcal{N}(\bm{\mu},\bm{\Sigma}),i=1,\dots,n$,
where $\bm{\mu}=\operatorname{vec}(\bm{M})$ and
$\bm{\Sigma}=\bm{\Sigma}^{\operatorname{col}}\otimes\bm{\Sigma}^{\operatorname{row}}$
denote the mean and covariance matrix, respectively. Then the maximum
likelihood estimators for mean and covariance are given by

$\displaystyle\hat{\bm{\mu}}=\frac{1}{n}\sum_{i=1}^{n}\bm{x}_{i}\quad\text{and}\quad\hat{\bm{\Sigma}}=\frac{1}{n}\sum_{i=1}^{n}(\bm{x}_{i}-\hat{\bm{\mu}})(\bm{x}_{i}-\hat{\bm{\mu}})^{\prime},$
(31)

respectively. The computation of the MLEs for matrix-variate samples based on
Equations (2)-(4) involves estimating $p(p+1)/2+q(q+1)/2$ parameters instead
of $pq(pq+1)/2+1$ parameters for the vectorized observations according to
Equation (31). This raises the question of whether fewer than $pq+1$
observations are sufficient for guaranteeing the existence and uniqueness of
MLEs for i.i.d. samples from a matrix normal distribution. This question was
investigated in several papers, such as Dutilleul, (1999); Lu and Zimmerman,
(2005); Srivastava et al., (2008); Roś et al., (2016); Soloveychik and
Trushin, (2016). We rely on the latter for the most recent proof of those
conditions. Note that it is not necessary to assume that the sample consists
of i.i.d. observations. In fact, the i.i.d. assumption can be relaxed to allow
for statistically dependent samples and it is not even necessary to require
identical distribution (Soloveychik and Trushin,, 2016, Remarks 2 and 6). The
critical condition for existence and uniqueness is that the sample contains at
least $n\geq\lfloor\nicefrac{{p}}{{q}}+\nicefrac{{q}}{{p}}\rfloor+2$
observations that are not collinear. The same holds for the existence and
uniqueness of the MMCD estimators, where $n$ is replaced by $h$, and for
properties like the breakdown point the assumptions could be relaxed only
requiring that the sample is in general position, i.e., no subset of $r,2\leq
r\leq\lfloor\nicefrac{{p}}{{q}}+\nicefrac{{q}}{{p}}\rfloor+2$ samples lies on
an $r-2$ dimensional subspace. However, the i.i.d. assumption is still
necessary when we consider properties like consistency.

The idea of the multivariate MCD estimator is as follows: Let
$\bm{x}_{i}=(x_{i1},\ldots,x_{ip})^{\prime}\in\mathbb{R}^{p}$ denote the
$i$-th observation of a data set in the multivariate setting, where
$i=1,\ldots,n$. The objective of the MCD estimator is to find the subset of
$h$ out of $n$ observations whose sample covariance matrix has the lowest
determinant, with $\nicefrac{{n}}{{2}}\leq h\leq n$ and $h>p$. In total, there
are $\binom{h}{n}$ possible $h$-subsets, and thus, a strategy needs to be used
to tackle the optimization problem efficiently. This has been done with the
so-called Fast-MCD algorithm (Rousseeuw and Driessen,, 1999), which internally
sorts the observations based on their Mahalanobis distances. For an
observation $\bm{x}_{i}$ from a population with mean
$\bm{\mu}\in\mathbb{R}^{p}$ and covariance
$\bm{\Sigma}\in\operatorname{PSD}(p)$ it is given by

$\displaystyle\operatorname{MD}(\bm{x}_{i},\bm{\mu},\bm{\Sigma})=\sqrt{(\bm{x}_{i}-\bm{\mu})^{\prime}\bm{\Sigma}^{-1}(\bm{x}_{i}-\bm{\mu})}.$

Since the Mahalanobis distance is vital for the computation of the MCD
estimator, it will also be crucial in a matrix-variate extension, where it can
be directly derived from the Mahalanobis distance of a vectorized matrix-
variate observation $\bm{X}$ as

$\displaystyle\operatorname{MMD}^{2}(\bm{X})=\operatorname{MMD}^{2}(\bm{X};\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}})$
$\displaystyle=\operatorname{MD}^{2}(\operatorname{vec}(\bm{X}))$
$\displaystyle=\operatorname{vec}(\bm{X}-\bm{\operatorname{M}})^{\prime}(\bm{\Omega}^{\operatorname{col}}\otimes\bm{\Omega}^{\operatorname{row}})\operatorname{vec}(\bm{X}-\bm{\operatorname{M}})$
$\displaystyle=\sum_{i=1}^{p}\sum_{j=1}^{p}\sum_{k=1}^{q}\sum_{l=1}^{q}(x_{ik}-m_{ik})(x_{jl}-m_{jl})\omega_{ij}^{\operatorname{row}}\omega_{kl}^{\operatorname{col}}$
$\displaystyle=\operatorname{tr}(\bm{\Omega}^{\operatorname{col}}(\bm{X}-\bm{\operatorname{M}})^{\prime}\bm{\Omega}^{\operatorname{row}}(\bm{X}-\bm{\operatorname{M}})),$

where $m_{ij}$, $\omega_{ij}^{\operatorname{row}}$ and
$\omega_{ij}^{\operatorname{col}}$ denote the elements $(i,j)$ of the matrices
$\bm{\operatorname{M}}$, $\bm{\Omega}^{\operatorname{row}}$ and
$\bm{\Omega}^{\operatorname{col}}$, respectively. If $\bm{X}$ has a matrix
normal distribution, then the squared matrix Mahalanbois distance has a
$\chi^{2}$ distribution with $pq$ degrees of freedom,
$\operatorname{MMD}^{2}(\bm{X})\sim\chi^{2}_{pq}$ (Gupta and Nagar,, 1999).

##  Appendix B Proofs of Section 2

###### Proof of Proposition 2.0.1.

In optimization problem (7) we want to maximize

$\displaystyle\begin{split}l(\bm{w},\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}}|\bm{\mathfrak{X}})=&-\frac{1}{2}\sum_{i=1}^{n}w_{i}\Bigl{(}p\ln(\det(\bm{\Sigma}^{\operatorname{col}}))+q\ln(\det(\bm{\Sigma}^{\operatorname{row}}))\Bigr{)}\\\
&-\frac{1}{2}\sum_{i=1}^{n}w_{i}\operatorname{MMD}^{2}(\bm{X}_{i})-hpq\ln(2\pi)\end{split}$
(32)

subject to $w_{i}\in\\{0,1\\}$ for all $i=1,\dots,n$ and
$\sum_{i=1}^{n}w_{i}=h$. In Equation (32),
$\operatorname{MMD}^{2}(\bm{X}_{i})$ is defined as in Equation (6).

For any random $h$-subset $H$ (or equivalently the corresponding set of
weights $\bm{w}$) the constrained MLEs for $\bm{\operatorname{M}}$,
$\bm{\Sigma}^{\operatorname{row}}$, and $\bm{\Sigma}^{\operatorname{col}}$ of
Equation (32) can be written as:

$\displaystyle\hat{\bm{\operatorname{M}}}_{H}$
$\displaystyle=\frac{1}{h}\sum_{i=i}^{n}w_{i}\bm{X}_{i}=\frac{1}{h}\sum_{i\in
H}\bm{X}_{i}$ $\displaystyle\hat{\bm{\Sigma}}^{\operatorname{row}}_{H}$
$\displaystyle=\frac{1}{qh}\sum_{i=i}^{n}w_{i}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})\hat{\bm{\Omega}}^{\operatorname{col}}_{H}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})^{\prime}=\frac{1}{qh}\sum_{i\in
H}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})\hat{\bm{\Omega}}^{\operatorname{col}}_{H}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})^{\prime}$
$\displaystyle\hat{\bm{\Sigma}}^{\operatorname{col}}_{H}$
$\displaystyle=\frac{1}{ph}\sum_{i=i}^{n}w_{i}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})^{\prime}\hat{\bm{\Omega}}^{\operatorname{row}}_{H}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})=\frac{1}{ph}\sum_{i\in
H}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})^{\prime}\hat{\bm{\Omega}}^{\operatorname{row}}_{H}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})$

Using those estimators to compute the sum of the Mahalanobis distances
$\operatorname{MMD}^{2}(\bm{X}_{i})$ in Equation (32) we obtain

$\displaystyle\sum_{i=1}^{n}w_{i}\operatorname{MMD}^{2}(\bm{X}_{i})=$
$\displaystyle\sum_{i\in
H}\operatorname{tr}\left(\hat{\bm{\Omega}}^{\operatorname{col}}_{H}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})^{\prime}\hat{\bm{\Omega}}^{\operatorname{row}}_{H}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})\right)$
$\displaystyle=$ $\displaystyle\sum_{i\in
H}\operatorname{tr}\left((\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})\hat{\bm{\Omega}}^{\operatorname{col}}_{H}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})^{\prime}\hat{\bm{\Omega}}^{\operatorname{row}}_{H}\right)$
$\displaystyle=$ $\displaystyle\operatorname{tr}\left(\sum_{i\in
H}((\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})\hat{\bm{\Omega}}^{\operatorname{col}}_{H}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{H})^{\prime})\hat{\bm{\Omega}}^{\operatorname{row}}_{H}\right)$
$\displaystyle=$
$\displaystyle\operatorname{tr}\left(qh\hat{\bm{\Sigma}}^{\operatorname{row}}_{H}\hat{\bm{\Omega}}^{\operatorname{row}}_{H}\right)=hpq.$

Thus, the terms in the second row of Equation (32) are all constant, and it is
sufficient to maximize only the term in the first row, which contains the
(negative) determinant of Equation (8). ∎

###  B.1 Properties of MMCD estimators

###### Proof of Lemma 3.0.1.

Ad (a): We show that the MMCD estimators are matrix affine equivariant. Let us
consider the objective of the MMCD for the transformed samples, which is to
minimize

$\displaystyle\det(\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Z}}_{H}}\otimes\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Z}}_{H}})$
$\displaystyle=\det\Bigl{(}(\bm{B}^{\prime}\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{X}}_{H}}\bm{B})\otimes(\bm{A}\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{X}}_{H}}\bm{A}^{\prime})\Bigr{)}$
$\displaystyle=\Bigl{[}\det(\bm{B}^{\prime}\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{X}}_{H}}\bm{B})\Bigr{]}^{p}\Bigl{[}\det(\bm{A}\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{X}}_{H}}\bm{A}^{\prime})\Bigr{]}^{q}$
$\displaystyle=\Bigl{[}\det(\bm{B}^{\prime})\det(\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{X}}_{H}})\det(\bm{B})\Bigr{]}^{p}\Bigl{[}\det(\bm{A})\det(\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{X}}_{H}})\det(\bm{A}^{\prime})\Bigr{]}^{q}$
$\displaystyle=4\det(\bm{B})^{p}\det(\bm{A})^{q}\det(\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{X}}_{H}})^{p}\det(\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{X}}_{H}})^{q}.$

Since $4\det(\bm{B})^{p}\det(\bm{A})^{q}$ is constant, the objective does not
change, and we obtain the same $h$-subset. Since the MMCD estimators
correspond to the trimmed MLEs and the objective is not affected by the
transformation, the matrix affine equivariance of the MMCD estimators follows
from the matrix affine equivariance of the MLEs.

Ad (b): Suppose that
$(\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Z}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Z}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Z}}})$
are matrix affine equivariant estimators of location and covariance of the
transformed sample $\bm{\mathfrak{Z}}$, then
$\operatorname{MMD}^{2}(\bm{Z}_{i};\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Z}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Z}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Z}}})\overset{!}{=}\operatorname{MMD}^{2}(\bm{X}_{i};\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{X}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{X}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{X}}})$:

$\displaystyle\operatorname{MMD}^{2}(\bm{Z}_{i};\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Z}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Z}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Z}}})$
$\displaystyle=$
$\displaystyle\operatorname{tr}(\hat{\bm{\Omega}}^{\operatorname{col}}_{\bm{\mathfrak{Z}}}(\bm{Z}_{i}-\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Z}}})^{\prime}\hat{\bm{\Omega}}^{\operatorname{row}}_{\bm{\mathfrak{Z}}}(\bm{Z}_{i}-\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Z}}}))$
$\displaystyle=$
$\displaystyle\operatorname{tr}\Big{(}\big{(}\bm{B}^{-1}\hat{\bm{\Omega}}^{\operatorname{col}}_{\bm{\mathfrak{X}}}(\bm{B}^{\prime})^{-1}(\bm{A}\bm{X}_{i}\bm{B}+\bm{C}-(\bm{A}\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{X}}}\bm{B}+\bm{C}))^{\prime}\big{)}$
$\displaystyle\qquad\big{(}(\bm{A}^{\prime})^{-1}\hat{\bm{\Omega}}^{\operatorname{row}}_{\bm{\mathfrak{X}}}\bm{A}^{-1}(\bm{A}\bm{X}_{i}\bm{B}+\bm{C}-(\bm{A}\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{X}}}\bm{B}+\bm{C}))\big{)}\Big{)}$
$\displaystyle=$
$\displaystyle\operatorname{tr}(\bm{B}^{-1}\hat{\bm{\Omega}}^{\operatorname{col}}_{\bm{\mathfrak{X}}}(\bm{B}^{\prime})^{-1}\bm{B}^{\prime}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{X}}})^{\prime}\bm{A}^{\prime}(\bm{A}^{\prime})^{-1}\hat{\bm{\Omega}}^{\operatorname{row}}_{\bm{\mathfrak{X}}}\bm{A}^{-1}\bm{A}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{X}}})\bm{B})$
$\displaystyle=$
$\displaystyle\operatorname{tr}(\hat{\bm{\Omega}}^{\operatorname{col}}_{\bm{\mathfrak{X}}}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{X}}})^{\prime}\hat{\bm{\Omega}}^{\operatorname{row}}_{\bm{\mathfrak{X}}}(\bm{X}_{i}-\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{X}}}))=\operatorname{MMD}^{2}(\bm{X}_{i};\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{X}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{X}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{X}}}).$

∎

The proofs of Theorems 3.0.1 and 3.0.3 require some definitions and properties
related to the vector space of matrices, which are introduced before the
proofs of the theorems. Since all matrices of a fixed size form a vector
space, objects such as ellipsoids or a simplex that are defined on the more
common vector spaces are also defined here. Let

$\displaystyle
E(\bm{T},\bm{U},\bm{V})=\\{\bm{X}:\operatorname{tr}(\bm{V}^{-1}(\bm{X}-\bm{T})^{\prime}\bm{U}^{-1}(\bm{X}-\bm{T}))\leq
1\\}$ (33)

be the ellipsoid containing the matrices $\bm{X}\in\mathbb{R}^{p\times q}$
with $\operatorname{MMD}^{2}(\bm{X};\bm{T},\bm{U},\bm{V})\leq 1$, where
$\bm{T}\in\mathbb{R}^{p\times q}$, $\bm{U}\in\operatorname{PSD}(p)$ and
$\bm{V}\in\operatorname{PSD}(q)$. The volume of this ellipsoid is given by

$\displaystyle\operatorname{vol}(E(\bm{T},\bm{U},\bm{V}))=\underbrace{\frac{\pi^{\nicefrac{{pq}}{{2}}}}{\Gamma(\nicefrac{{pq}}{{2}}+1)}}_{=:\beta_{pq}}\prod_{i=1}^{p}\prod_{j=1}^{q}\sqrt{\lambda_{i}(\bm{U})\lambda_{j}(\bm{V})}=\beta_{pq}\underbrace{\det(\bm{U})^{\nicefrac{{q}}{{2}}}\det(\bm{V})^{\nicefrac{{p}}{{2}}}}_{=:\det(E(\bm{T},\bm{U},\bm{V}))},$
(34)

where $\Gamma$ is the gamma function,
$0<\lambda_{p}(\bm{U})\leq\ldots\leq\lambda_{1}(\bm{U})$ and
$0<\lambda_{q}(\bm{V})\leq\ldots\leq\lambda_{1}(\bm{V})$ are the eigenvalues
of $\bm{U}$ and $\bm{V}$, respectively. Moreover, the axes have lengths
$\sqrt{\lambda_{i}(\bm{U})\lambda_{j}(\bm{V})}$.

Let $\bm{A}$ be a symmetric nonnegative definite $p\times p$ matrix, then

$\displaystyle\lambda_{1}(\bm{A})=\sup_{\bm{z}\in\mathbb{R}^{p}}\frac{\bm{z}^{\prime}\bm{A}\bm{z}}{\bm{z}^{\prime}\bm{z}}\quad\text{and}\quad\lambda_{n}(\bm{A})=\inf_{\bm{z}\in\mathbb{R}^{p}}\frac{\bm{z}^{\prime}\bm{A}\bm{z}}{\bm{z}^{\prime}\bm{z}}.$
(35)

Consider another symmetric nonnegative definite $p\times p$ matrix $\bm{B}$,
then using Equation (35) we get that

$\displaystyle\lambda_{1}(\bm{A}+\bm{B})\leq\lambda_{1}(\bm{A})+\lambda_{1}(\bm{B})\quad\text{and}\quad\lambda_{p}(\bm{A}+\bm{B})\geq\lambda_{p}(\bm{A})+\lambda_{p}(\bm{B}).$
(36)

If $\bm{A}\in\operatorname{PSD}(p)$ with eigenvalues
$0<\lambda_{p}(\bm{A})\leq\ldots\leq\lambda_{1}(\bm{A})$ then the eigenvalues
of $\bm{A}^{-1}$ are the reciprocals of the eigenvalues of $\bm{A}$, i.e.
$\lambda_{i}(\bm{A}^{-1})=\lambda^{-1}_{i}(\bm{A})$. Hence, we have that

$\displaystyle\frac{1}{\lambda_{1}(\bm{A})}=\inf_{\bm{z}\in\mathbb{R}^{p}}\frac{\bm{z}^{\prime}\bm{A}^{-1}\bm{z}}{\bm{z}^{\prime}\bm{z}},$

which implies that for any $\bm{x}\in\mathbb{R}^{p}$

$\displaystyle\frac{1}{\lambda_{1}(\bm{A})}\leq\frac{\bm{x}^{\prime}\bm{A}^{-1}\bm{x}}{\bm{x}^{\prime}\bm{x}}\Leftrightarrow\bm{x}^{\prime}\bm{x}\leq\bm{x}^{\prime}\bm{A}^{-1}\bm{x}\lambda_{1}(\bm{A}).$
(37)

Suppose $\bm{A}\in\operatorname{PSD}(p)$, $\bm{B}\in\operatorname{PSD}(q)$ and
let $\lambda(\bm{A})$ be an eigenvalue of $\bm{A}$ with corresponding
eigenvector $\bm{v}(\bm{A})$, and $\lambda(\bm{B})$ an eigenvalue of $\bm{B}$
with corresponding eigenvector $\bm{v}(\bm{B})$. Then
$\lambda(\bm{A})\lambda(\bm{B})$ is an eigenvalue of $\bm{B}\otimes\bm{A}$
with corresponding eigenvector $\bm{v}(\bm{B})\otimes\bm{v}(\bm{A})$. We
denote the sequence of eigenvalues of $\bm{A}$ and $\bm{B}$ as
$0<\lambda_{p}(\bm{A})\leq\ldots\leq\lambda_{1}(\bm{A})$ and
$0<\lambda_{q}(\bm{b})\leq\ldots\leq\lambda_{1}(\bm{b})$, respectively. It
follows that the smallest eigenvalue
$\lambda_{pq}(\bm{A},\bm{B})=\lambda_{p}(\bm{A})\lambda_{q}(\bm{B})$ and the
largest eigenvalue
$\lambda_{1}(\bm{A},\bm{B})=\lambda_{1}(\bm{A})\lambda_{1}(\bm{B})$. Moreover
note that for $\bm{Z}\in\mathbb{R}^{p\times q}$

$\displaystyle\operatorname{vec}(\bm{Z})^{\prime}(\bm{B}\otimes\bm{A})\operatorname{vec}(\bm{Z})=\operatorname{tr}(\bm{B}\bm{Z}^{\prime}\bm{A}\bm{Z}),$

as in Equation (6), which implies that

$\displaystyle\lambda_{pq}(\bm{A},\bm{B})=\inf_{\bm{Z}\in\mathbb{R}^{p\times
q}}\frac{\operatorname{tr}(\bm{B}\bm{Z}^{\prime}\bm{A}\bm{Z})}{\operatorname{tr}(\bm{Z}^{\prime}\bm{Z})}\quad\text{and}\quad\lambda_{1}(\bm{A},\bm{B})=\sup_{\bm{Z}\in\mathbb{R}^{p\times
q}}\frac{\operatorname{tr}(\bm{B}\bm{Z}^{\prime}\bm{A}\bm{Z})}{\operatorname{tr}(\bm{Z}^{\prime}\bm{Z})}.$

This leads us to the matrix-variate version of Equation (37), where for any
matrix $\bm{X}\in\mathbb{R}^{p\times q}$

$\displaystyle\left\lVert\bm{X}\right\rVert_{F}^{2}=\operatorname{tr}(\bm{X}^{\prime}\bm{X})\leq\operatorname{tr}(\bm{B}^{-1}\bm{X}^{\prime}\bm{A}^{-1}\bm{X})\lambda_{1}(\bm{A},\bm{B})=\operatorname{tr}(\bm{B}^{-1}\bm{X}^{\prime}\bm{A}^{-1}\bm{X})\lambda_{1}(\bm{A})\lambda_{1}(\bm{B})$
(38)

######  Lemma B.1.1.

Take $p,q\in\mathbb{N}$,
$d=\lfloor\nicefrac{{p}}{{q}}+\nicefrac{{q}}{{p}}\rfloor$, $d+2\leq s\leq pq$,
and matrices $\bm{X}_{1},\ldots,\bm{X}_{s}\in\mathbb{R}^{p\times q}$ that are
in general position, i.e., no subset of $r,2\leq r\leq s$ samples lies on an
$r-2$ dimensional subspace. For an ellipsoid $E(\bm{T},\bm{U},\bm{V})$ as in
Equation (33), containing the matrices $\bm{X}_{1},\ldots,\bm{X}_{s}$, it
holds that for every $C>0$ there exists a constant
$\alpha:=\alpha(\bm{X}_{1},\ldots,\bm{X}_{s})>0$ only depending on
$\bm{X}_{1},\ldots,\bm{X}_{s}$ such that
$\left\lVert\bm{T}\right\rVert_{F}=\sqrt{\operatorname{tr}(\bm{T}^{\prime}\bm{T})}>\alpha$
implies $\det(E(\bm{T},\bm{U},\bm{V}))>C$, i.e.,

$\displaystyle\forall
C>0~{}\exists\alpha>0:\left\lVert\bm{T}\right\rVert_{F}>\alpha\implies\det(E(\bm{T},\bm{U},\bm{V}))>C.$

###### Proof.

The samples $\bm{X}_{1},\ldots,\bm{X}_{s}$ are in general position, which
implies that they span a nonempty $s-1$ simplex. Since
$E(\bm{T},\bm{U},\bm{V})$ contains those samples, it also contains the simplex
spanned by those matrices. This implies that there exists a constant $a>0$
only depending on $\bm{X}_{1},\ldots,\bm{X}_{s}$, such that the length of $k$,
$s-1\leq k\leq pq$, of the $pq$ axes of the ellipsoid
$E(\bm{T},\bm{U},\bm{V})$ is at least $a$, i.e., there are $k$ out of $pq$
indices $(i,j),1\leq i\leq p,1\leq j\leq q$ such that

$\displaystyle\sqrt{\lambda_{i}(\bm{U})\lambda_{j}(\bm{V})}>a.$ (39)

In Equation (39), $\lambda_{i}(\bm{U}),i\in\\{1,\dots,p\\}$, are the
eigenvalues of $\bm{U}$, and $\lambda_{j}(\bm{V}),j\in\\{1,\dots,q\\}$, are
the eigenvalues of $\bm{V}$. For any matrix $\bm{X}$ contained in
$E(\bm{T},\bm{U},\bm{V})$, Equations (33) and (38) imply that

$\displaystyle\begin{split}\left\lVert\bm{X}-\bm{T}\right\rVert_{F}^{2}&=\operatorname{tr}((\bm{X}-\bm{T})^{\prime}(\bm{X}-\bm{T}))\\\
&\leq\operatorname{tr}(\bm{V}^{-1}(\bm{X}-\bm{T})^{\prime}\bm{U}^{-1}(\bm{X}-\bm{T}))\lambda_{1}(\bm{U})\lambda_{1}(\bm{V})\\\
&\leq\lambda_{1}(\bm{U})\lambda_{1}(\bm{V}).\end{split}$ (40)

Without loss of generality, we assume that the matrix of all zeros
$\bm{0}\in\mathbb{R}^{p\times q}$ is contained in the ellipsoid
$E(\bm{T},\bm{U},\bm{V})$, then Equation (40) implies that
$\left\lVert\bm{T}\right\rVert_{F}^{2}\leq\lambda_{1}(\bm{U})\lambda_{1}(\bm{V})$.
Take $\alpha=C/(a^{pq-1})$, then we have that

$\displaystyle\frac{C}{a^{pq-1}}<\left\lVert\bm{T}\right\rVert_{F}\leq\sqrt{\lambda_{1}(\bm{U})\lambda_{1}(\bm{V})}\Leftrightarrow
C<\sqrt{\lambda_{1}(\bm{U})\lambda_{1}(\bm{V})}a^{pq-1}$

and from Equaiton (39) it follows that

$\displaystyle\begin{split}\det(E(\bm{T},\bm{U},\bm{V})):&=\det(\bm{U})^{\nicefrac{{q}}{{2}}}\det(\bm{V})^{\nicefrac{{p}}{{2}}}\\\
&=\prod_{i=1}^{p}\prod_{j=1}^{q}\sqrt{\lambda_{i}(\bm{U})\lambda_{j}(\bm{V})}\\\
&>\sqrt{\lambda_{1}(\bm{U})\lambda_{1}(\bm{V})}a^{pq-1}>C.\end{split}$

∎

######  Proof of Theorem 3.0.1.

We show that the breakdown points of the MMCD estimators of location and
covariance defined in Equations (16) and (17), respectively, are both
$\nicefrac{{m}}{{n}}$, with $m=\lfloor\min(n-h+1,h-(d+1))\rfloor$,
$d=\lfloor\nicefrac{{p}}{{q}}+\nicefrac{{q}}{{p}}\rfloor$. First, we prove
that
$\varepsilon^{\ast}(\hat{\bm{\operatorname{M}}},\bm{\mathfrak{X}})=\varepsilon^{\ast}(\hat{\bm{\Sigma}}^{\operatorname{row}},\hat{\bm{\Sigma}}^{\operatorname{col}},\bm{\mathfrak{X}})\geq\nicefrac{{m}}{{n}}$.
Let $\bm{\mathfrak{Y}}$ be the sample obtained by replacing at most $m-1$
matrices of $\bm{\mathfrak{X}}$ by arbitrary $p\times q$ matrices. Since
$n-(m-1)\geq h$, $\bm{\mathfrak{Y}}$ contains at least $h$ matrices of the
orginial sample $\bm{\mathfrak{X}}$ and because $m-1\leq h-(d+1)-1$, every
subset of size $h$ of $\bm{\mathfrak{Y}}$ includes at least $d+2$ matrices of
the original sample $\bm{\mathfrak{X}}$. Hence, the MMCD estimators can almost
surely be computed for any $h$-subset of $\bm{\mathfrak{Y}}$. Let us consider
three ellipsoids:

  * •

Let $E_{\max}=E(\bm{0},c_{\max}\bm{I},\bm{I})$ denote the smallest sphere that
contains all samples in $\bm{\mathfrak{X}}$, where $c_{\max}$ is chosen
accordingly.

  * •

Let $E_{h}=E(\bm{0},c_{h}\bm{I},\bm{I})$ denote the smallest sphere that
contains the $h$ samples of $\bm{\mathfrak{X}}$ that are also in
$\bm{\mathfrak{Y}}$, where $c_{h}$ is chosen accordingly.

  * •

Let
$E_{\operatorname{MMCD}}=E(\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Y}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Y}}})$
denote the MMCD ellipsoid.

It follows that
$\det(E_{\operatorname{MMCD}})\leq\det(E_{h})\leq\det(E_{\max})=:\alpha$,
where for an ellipsoid $E=E(\bm{T},\bm{U},\bm{V})$, $\det(E)$ is defined in
(34). Note that $\bm{\mathfrak{X}}$ is a collection of random samples from a
continuous distribution and therefore it is in general position almost surely.
Further, $E_{\operatorname{MMCD}}$ covers at least $h$ samples, and those
include at least $d+2$ samples of $\bm{\mathfrak{X}}$, which span a nonempty
$d+1$ simplex. Lemma B.1.1 shows that there exists a constant $\alpha>0$ that
only depends on those $d+2$ samples such that, if
$\left\lVert\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}}\right\rVert_{F}>C$
it would imply $\det(E_{\operatorname{MMCD}})>\alpha$. As shown above, this is
not possible, hence
$\left\lVert\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}}\right\rVert_{F}\leq
C$.

Similarly, since $\bm{\mathfrak{Y}}$ contains at least $d+2$ matrices of the
original sample $\bm{\mathfrak{X}}$, the MMCD estimators almost surely yield
positive definite covariance estimates
$\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Y}}}$ and
$\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Y}}}$. More
specifically, let $\bm{\mathfrak{X}}_{T}$, $T\subseteq H$ be the subset of the
at least $d+2$ matrices of the original sample that are in
$\bm{\mathfrak{Y}}$. Since $\left|T\right|\geq d+2=\lfloor p/q+q/p\rfloor+2$
the MLE estimators
$(\hat{\bm{M}}_{\bm{\mathfrak{X}}_{T}},\hat{\bm{\Sigma}}_{\bm{\mathfrak{X}}_{T}}^{\operatorname{row}},\hat{\bm{\Sigma}}_{\bm{\mathfrak{X}}_{T}}^{\operatorname{row}})$
of this subsample are almost surely positive definite. Let
$E_{T}=E(\hat{\bm{M}}_{\bm{\mathfrak{X}}_{T}},\hat{\bm{\Sigma}}_{\bm{\mathfrak{X}}_{T}}^{\operatorname{row}},\hat{\bm{\Sigma}}_{\bm{\mathfrak{X}}_{T}}^{\operatorname{row}})$
denote the corresponding ellipsoid which is the smallest ellipsoid, of the
type $E=E(\textbf{T},\textbf{U},\textbf{V})$ as in Equation (33), containing
the samples $\bm{\mathfrak{X}}_{T}$ as one can think of it as the MMCD
ellipsoid for those $\left|T\right|\geq d+2$ samples with $H=T$. This further
implies that the volume of the corresponding ellipsoid $E_{T}$ is bounded from
below by a constant only depending on $\bm{\mathfrak{X}}$, i.e.
$\det(E_{T})\geq v>0$. As $E_{\operatorname{MMCD}}$ is also an ellipsoid
containing the samples $\bm{\mathfrak{X}}_{T}$,
$\det(E_{\operatorname{MMCD}})\geq\det(E_{T})\geq v>0$. Moreover, it also
means that there exists a constant $k$ depending only on $\bm{\mathfrak{X}}$,
such that $E_{T}\subseteq kE_{\operatorname{MMCD}}$, implying that there
exists a constant $\gamma>0$ depending only on $\bm{\mathfrak{X}}$, such that
$\lambda_{i}(\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Y}}})\lambda_{j}(\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Y}}})>\gamma,1\leq
i\leq p,1\leq j\leq q$. Especially,
$\lambda_{p}(\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Y}}})\lambda_{q}(\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Y}}})>\gamma$.
Since also $\det(E_{\operatorname{MMCD}})\leq\alpha$ there exists a constant
$\delta>0$, depending only on $\bm{\mathfrak{X}}$ such that
$\lambda_{i}(\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Y}}})\lambda_{j}(\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Y}}})<\delta,1\leq
i\leq p,1\leq j\leq q$.

Next we show that
$\varepsilon^{\ast}(\hat{\bm{\operatorname{M}}},\bm{\mathfrak{X}})=\varepsilon^{\ast}(\hat{\bm{\Sigma}}^{\operatorname{row}},\hat{\bm{\Sigma}}^{\operatorname{col}},\bm{\mathfrak{X}})\leq\nicefrac{{m}}{{n}}$.
If $m=n-h+1$, we replace $m=n-h+1$ matrices of $\bm{\mathfrak{X}}$ to obtain
$\bm{\mathfrak{Y}}$, then $n-m=h-1$, implying that every subset of $h$ samples
of $\bm{\mathfrak{Y}}$ contains at least one conaminated sample. Hence,
$E_{\operatorname{MMCD}}=E(\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Y}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Y}}})$
also includes at least one contaminated sample. Let
$\left\lVert\bm{X}\right\rVert_{F}\rightarrow\infty$ for all contaminated
samples $\bm{X}$, then at least one eigenvalue of $E_{\operatorname{MMCD}}$
explodes and the MMCD location and covariance estimators break down. Finally,
consider the case where $m=h-(d+1)$. To construct $\bm{\mathfrak{Y}}$, take
any $d+1$ samples of $\bm{\mathfrak{X}}$ and consider the $d$ dimensional
hyperplane $L$ they determine. Replace $h-(d+1)$ samples that are not in $L$
and replace them with matrices on $L$. Then $L$ contains $h$ points of
$\bm{\mathfrak{Y}}$ and the ellipsoid covering those points has volume zero
and hence determinant zero. Since $\bm{\mathfrak{X}}$ is in general position,
we can construct $\bm{\mathfrak{Y}}$ such that no other lower dimensional
hyperplane contains $h$ points of $\bm{\mathfrak{Y}}$. Hence,
$\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}}$ lies on $L$ and
$E_{\operatorname{MMCD}}=E(\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Y}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Y}}})$
has zero determinant. This implies that at least one eigenvalue is zero, hence
the MMCD location and covariance estimators break down. ∎

######  Proof of Theorem 3.0.2.

Let $(\bm{X}_{1},\dots,\bm{X}_{n})$ be a sample of matrix-variate observations
and $(\bm{x}_{1},\dots,\bm{x}_{n})$, $\bm{x}_{i}=\mathrm{vec}(\bm{X}_{i})$,
$i=1,\dots n$ its vectorized form. The MCD estimator can also be found as a
solution to the following maximization problem:

$\displaystyle\max_{\bm{w},\hat{\bm{\mu}},\hat{\bm{\Sigma}}}l(\bm{w},\hat{\bm{\mu}},\hat{\bm{\Sigma}}|(\bm{x}_{1},\dots,\bm{x}_{n}))=\max_{\bm{w},\hat{\bm{\mu}},\hat{\bm{\Sigma}}}-\frac{1}{2}\sum_{i=1}^{n}w_{i}(\ln(\det(\hat{\bm{\Sigma}}))+pq\ln(2\pi)+\operatorname{MD}^{2}(\bm{x}_{i},\hat{\bm{\mu}},\hat{\bm{\Sigma}}))$

subject to $w_{1},\dots,w_{n}\in\\{0,1\\}$, $\sum_{i=1}^{n}w_{i}=h$,
$\hat{\bm{\mu}}\in\mathbb{R}^{pq}$,
$\hat{\bm{\Sigma}}\in\operatorname{PSD}(pq)$; see Raymaekers and Rousseeuw,
(2023) for more insight. Similarly, the MMCD estimator is a solution to the
following maximization problem:

$\displaystyle\max_{\bm{w},\hat{\bm{\operatorname{M}}},\hat{\bm{\Sigma}}^{\operatorname{row}},\hat{\bm{\Sigma}}^{\operatorname{col}}}l(\bm{w},\hat{\bm{\operatorname{M}}},\hat{\bm{\Sigma}}^{\operatorname{row}},\hat{\bm{\Sigma}}^{\operatorname{col}}|(\bm{X}_{1},\dots,\bm{X}_{n}))$
$\displaystyle=$
$\displaystyle\max_{\bm{w},\hat{\bm{\operatorname{M}}},\hat{\bm{\Sigma}}^{\operatorname{row}},\hat{\bm{\Sigma}}^{\operatorname{col}}}-\frac{1}{2}\sum_{i=1}^{n}w_{i}\Bigl{(}p\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{col}}))+q\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{row}}))+\operatorname{MMD}^{2}(\bm{X}_{i})+pq\ln(2\pi)\Bigr{)}$

subject to $w_{1},\dots,w_{n}\in\\{0,1\\}$, $\sum_{i=1}^{n}w_{i}=h$,
$\hat{\bm{\operatorname{M}}}\in\mathbb{R}^{p\times q}$,
$\hat{\bm{\Sigma}}^{\operatorname{row}}\in\operatorname{PSD}(p)$,
$\hat{\bm{\Sigma}}^{\operatorname{col}}\in\operatorname{PSD}(q)$; see
Proposition 2.0.1.

Denote further
$(\bm{w}_{\operatorname{MCD}},\hat{\bm{\mu}}_{\operatorname{MCD}},\hat{\bm{\Sigma}}_{\operatorname{MCD}})$
and
$(\bm{w}_{\operatorname{MMCD}},\hat{\bm{\mu}}_{\operatorname{MMCD}},\hat{\bm{\Omega}}_{\operatorname{MMCD}}\otimes\hat{\bm{\Sigma}}_{\operatorname{MMCD}})$
weights, mean and covariance estimators for the vectorized sample
$(\bm{x}_{1},\dots,\bm{x}_{n})$, based on MCD and MMCD, respectively. As
$\bm{X}_{i}\sim\mathcal{MN}(\bm{\operatorname{M}},\bm{\Sigma},\bm{\Omega})$,
then $\bm{x}_{i}\sim\mathcal{N}(\bm{\mu},\bm{\Omega}\otimes\bm{\Sigma})$, with
$\bm{\mu}=\operatorname{vec}(\bm{\operatorname{M}})$ and $i=1,\dots,n$.
Moreover, the mean estimator $\hat{\bm{\mu}}_{\operatorname{MCD}}$ and
properly scaled covariance estimator $\hat{\bm{\Sigma}}_{\operatorname{MCD}}$
are strongly consistent for the population counterparts $\bm{\mu}$ and
$\bm{\Omega}\otimes\bm{\Sigma}$; see e.g. Croux and Haesbroeck, (1999) and
Cator and Lopuhaä, (2012). Especially, this implies that for every $\delta>0$
there exists $n\in\mathbb{N}$ such that

$\left\lVert\hat{\bm{\mu}}_{\operatorname{MCD}}-\bm{\mu}\right\rVert\overset{a.s.}{<}\delta,\quad\left\lVert\hat{\bm{\Sigma}}_{\operatorname{MCD}}-\bm{A}\otimes\bm{B}\right\rVert\overset{a.s.}{<}\delta,$

for some
$\bm{A}\otimes\bm{B}\in\operatorname{PSD}(p)\otimes\operatorname{PSD}(q)$. In
the following, we will drop ${a.s.}$ superscript from (in)equality signs when
it is clear from the context. For fixed weights $\bm{w}$, the log-likelihood
function
$l_{\cdot,\bm{w}|(\bm{x}_{1},\dots,\bm{x}_{n})}:(\bm{\mu},\bm{\Sigma})\mapsto
l(\bm{\mu},\bm{\Sigma}|\bm{w},(\bm{x}_{1},\dots,\bm{x}_{n}))$ is continuous in
both $\bm{\mu}$, and $\bm{\Sigma}$, and its continuity implies

$\left|l(\bm{w}_{\operatorname{MCD}},\hat{\bm{\mu}}_{\operatorname{MCD}},\hat{\bm{\Sigma}}_{\operatorname{MCD}}|(\bm{x}_{1},\dots,\bm{x}_{n}))-l(\bm{w}_{\operatorname{MCD}},\bm{\mu},\bm{A}\otimes\bm{B}|(\bm{x}_{1},\dots,\bm{x}_{n}))\right|<\varepsilon,$

for $\varepsilon=\varepsilon(\delta)>0$. Moreover, the solution
$(\bm{w}_{\operatorname{MCD}},\hat{\bm{\mu}}_{\operatorname{MCD}},\hat{\bm{\Sigma}}_{\operatorname{MCD}})$
is optimal for $l(\cdot|(\bm{x}_{1},\dots,\bm{x}_{n}))$, implying that

$0<l(\bm{w}_{\operatorname{MCD}},\hat{\bm{\mu}}_{\operatorname{MCD}},\hat{\bm{\Sigma}}_{\operatorname{MCD}}|(\bm{x}_{1},\dots,\bm{x}_{n}))-l(\bm{w}_{\operatorname{MCD}},\bm{\mu},\bm{A}\otimes\bm{B}|(\bm{x}_{1},\dots,\bm{x}_{n}))<\varepsilon.$
(41)

Similarly,
$(\bm{w}_{\operatorname{MMCD}},\hat{\bm{\mu}}_{\operatorname{MMCD}},\hat{\bm{\Omega}}_{\operatorname{MMCD}}\otimes\hat{\bm{\Sigma}}_{\operatorname{MMCD}})$
is a maximizer of $l(\cdot|(\bm{x}_{1},\dots,\bm{x}_{n}))$ in the set of all
feasible weights, means, and covariances with Kronecker product structure. As
$(\bm{w}_{\operatorname{MCD}},\bm{\mu},\bm{A}\otimes\bm{B})$ belongs to the
same set,

$l(\bm{w}_{\operatorname{MMCD}},\hat{\bm{\mu}}_{\operatorname{MMCD}},\hat{\bm{\Omega}}_{\operatorname{MMCD}}\otimes\hat{\bm{\Sigma}}_{\operatorname{MMCD}}|(\bm{x}_{1},\dots,\bm{x}_{n}))>l(\bm{w}_{\operatorname{MCD}},\bm{\mu},\bm{A}\otimes\bm{B}|(\bm{x}_{1},\dots,\bm{x}_{n})).$

Denote further
$\hat{\bm{S}}_{\operatorname{MMCD}}=\frac{1}{h}\sum_{i=1}^{n}w_{{\operatorname{MMCD}},i}(\bm{x}_{i}-\hat{\bm{\mu}}_{\mathrm{MMCD}})(\bm{x}_{i}-\hat{\bm{\mu}}_{\mathrm{MMCD}})^{\prime}$
to be the estimate of the covariance of the vectorized observations, based on
the weights (subset) produced by the MMCD algorithm. As
$\hat{\bm{S}}_{\operatorname{MMCD}}$ is optimal for $l$ given fixed weights
$\bm{w}_{\mathrm{MMCD}}$,

$\displaystyle
l(\bm{w}_{\operatorname{MCD}},\hat{\bm{\mu}}_{\operatorname{MCD}},\hat{\bm{\Sigma}}_{\operatorname{MCD}}|(\bm{x}_{1},\dots,\bm{x}_{n}))$
$\displaystyle>l(\bm{w}_{\operatorname{MMCD}},\hat{\bm{\mu}}_{\operatorname{MMCD}},\hat{\bm{S}}_{\operatorname{MMCD}}|(\bm{x}_{1},\dots,\bm{x}_{n}))$
$\displaystyle>l(\bm{w}_{\operatorname{MMCD}},\hat{\bm{\mu}}_{\operatorname{MMCD}},\hat{\bm{\Omega}}_{\operatorname{MMCD}}\otimes\hat{\bm{\Sigma}}_{\operatorname{MMCD}}|(\bm{x}_{1},\dots,\bm{x}_{n}))$
$\displaystyle>l(\bm{w}_{\operatorname{MCD}},\bm{\mu},\bm{A}\otimes\bm{B}|(\bm{x}_{1},\dots,\bm{x}_{n})).$
(42)

(41) and (B.1) now give that

$0<l(\bm{w}_{\operatorname{MCD}},\hat{\bm{\mu}}_{\operatorname{MCD}},\hat{\bm{\Sigma}}_{\operatorname{MCD}}|(\bm{x}_{1},\dots,\bm{x}_{n}))-l(\bm{w}_{\operatorname{MMCD}},\hat{\bm{\mu}}_{\operatorname{MMCD}},\hat{\bm{S}}_{\operatorname{MMCD}}|(\bm{x}_{1},\dots,\bm{x}_{n}))<\varepsilon,$

i.e., due to Raymaekers and Rousseeuw, (2023),

$0<\det(\hat{\bm{S}}_{\operatorname{MMCD}})-\det(\hat{\bm{\Sigma}}_{\operatorname{MCD}})<\varepsilon,$
(43)

for $\varepsilon=\varepsilon(n)>0$, arbitrarily small $(\varepsilon(n)\to
0,\,n\to\infty)$. As both $\hat{\bm{\Sigma}}_{\operatorname{MCD}}$ and
$\hat{\bm{S}}_{\operatorname{MMCD}}$ are weighted sample covariances for the
random sample of vectorized observations calculated using the weights
satisfying the same constraints, Corollary 4.1. in Cator and Lopuhaä, (2012)
(taking $P_{t}$ to be the empirical measure based on the sample
$(\bm{x}_{1},\dots,\bm{x}_{n})$) implies that

$\hat{\bm{\mu}}_{\operatorname{MMCD}}\xrightarrow{a.s.}\bm{\mu},\quad\hat{\bm{S}}_{\operatorname{MMCD}}\xrightarrow{a.s.}c_{1}^{-1}\bm{\Omega}\otimes\bm{\Sigma}\quad,$
(44)

where $c_{1}>0$ is a distribution-specific consistency factor of MCD for the
given in Croux and Haesbroeck, (1999).

To complete the proof consider reparametrization of
$l(\bm{w},\bm{a},\bm{A}|(\bm{x}_{1},\dots,\bm{x}_{n}))$ for fixed weights
$\bm{w}\in\mathbb{R}^{n}$, mean $\bm{a}\in\mathbb{R}^{pq}$, and covariance
$\bm{A}\in\operatorname{PSD}(pq)$ in terms of the precicion matrix
$\bm{B}=\bm{A}^{-1}$. Denote this new parametrization as
$g(\bm{B}|\bm{w},\bm{a},(\bm{x}_{1},\dots,\bm{x}_{n}))=l(\bm{w},\bm{a},\bm{B}^{-1}|\bm{x}_{1},\dots,\bm{x}_{n})$,
which is now concave in $\bm{B}$. Especially, for
$\bm{w}=\bm{w}_{\operatorname{MMCD}}$ and
$\bm{a}=\hat{\bm{\mu}}_{\operatorname{MMCD}}$, the function
$g(\bm{B}|\bm{w}_{\operatorname{MMCD}},\hat{\bm{\mu}}_{\operatorname{MMCD}},(\bm{x}_{1},\dots,\bm{x}_{n}))$
is concave in $\bm{B}$ and achieves a unique global maximum at
$\bm{B}=\hat{\bm{S}}_{\operatorname{MMCD}}$. Equations (41) and (B.1) then
give

$\displaystyle 0<$ $\displaystyle
l(\bm{w}_{\operatorname{MMCD}},\hat{\bm{\mu}}_{\operatorname{MMCD}},\hat{\bm{S}}_{\operatorname{MMCD}}|(\bm{x}_{1},\dots,\bm{x}_{n}))$
$\displaystyle-l(\bm{w}_{\operatorname{MMCD}},\hat{\bm{\mu}}_{\operatorname{MMCD}},\hat{\bm{\Omega}}_{\operatorname{MMCD}}\otimes\hat{\bm{\Sigma}}_{\operatorname{MMCD}}|(\bm{x}_{1},\dots,\bm{x}_{n}))<\varepsilon,$

further implying that

$\displaystyle 0<$ $\displaystyle
g(\hat{\bm{S}}_{\operatorname{MMCD}}^{-1}|\bm{w}_{\operatorname{MMCD}},\hat{\bm{\mu}}_{\operatorname{MMCD}},(\bm{x}_{1},\dots,\bm{x}_{n}))$
$\displaystyle-g((\hat{\bm{\Omega}}_{\operatorname{MMCD}}\otimes\hat{\bm{\Sigma}}_{\operatorname{MMCD}})^{-1}|\bm{w}_{\operatorname{MMCD}},\hat{\bm{\mu}}_{\operatorname{MMCD}},(\bm{x}_{1},\dots,\bm{x}_{n}))<\varepsilon,$

as both $\hat{\bm{S}}_{\operatorname{MMCD}}$ and
$\hat{\bm{\Omega}}_{\operatorname{MMCD}}\otimes\hat{\bm{\Sigma}}_{\operatorname{MMCD}}$
are a.s. positive definite for $n$ large enough. Concavity of $g$ and the fact
that $\hat{\bm{S}}_{\operatorname{MMCD}}^{-1}$ is its global maximum further
imply that

$\|\hat{\bm{S}}_{\operatorname{MMCD}}^{-1}-(\hat{\bm{\Omega}}_{\operatorname{MMCD}}\otimes\hat{\bm{\Sigma}}_{\operatorname{MMCD}})^{-1}\|<\delta_{1},$

for $\delta_{1}=\delta_{1}(\varepsilon)\to 0$ as $n\to\infty$. Almost sure
positive definiteness of $\hat{\bm{S}}_{\operatorname{MMCD}}$ and
$\hat{\bm{\Omega}}_{\operatorname{MMCD}}\otimes\hat{\bm{\Sigma}}_{\operatorname{MMCD}}$,
and continuity of matrix inverse imply that

$\|\hat{\bm{S}}_{\operatorname{MMCD}}-\hat{\bm{\Omega}}_{\operatorname{MMCD}}\otimes\hat{\bm{\Sigma}}_{\operatorname{MMCD}}\|<\delta,$
(45)

for $\delta=\delta(\varepsilon)\to 0$ as $n\to\infty$. Equations (44) and (45)
now complete the proof. Observe that the proof indicates that the
distribution-specific consistency factor is inherited from the MCD covariance
estimator; see Croux and Haesbroeck, (1999). ∎

######  Proof of Theorem 3.0.3.

We show that the breakdown points of the reweighted MMCD estimators are at
least as high as the breakdown points of the raw MMCD estimators. Let
$\bm{\mathfrak{Y}}$ be the sample obtained by replacing at most $m-1$ matrices
of $\bm{\mathfrak{X}}$ by arbitrary $p\times q$ matrices. Let
$\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}}$,
$\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Y}}}$, and
$\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Y}}}$ denote the _raw_
MMCD estimators and $\tilde{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}}$,
$\tilde{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Y}}}$, and
$\tilde{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Y}}}$ denote the
_reweighted_ MMCD estimators based on the corrupted sample
$\bm{\mathfrak{Y}}$. Further,
$d(\bm{Y}_{i})=\operatorname{MMD}(\bm{Y}_{i};\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Y}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Y}}})$,
$i\in N=\\{1,\dots,n\\}$, denote the matrix Mahalanbois distances of the
corrupted sample based on the _raw_ MMCD estimators. Since
$m\leq\varepsilon^{\ast}(\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{X}}},\bm{\mathfrak{X}})-1=\varepsilon^{\ast}(\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{X}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{X}}},\bm{\mathfrak{X}})-1$
it follows that there exist constants $k_{0}$, $k_{1}$, and $k_{2}$ that only
depend on $\bm{\mathfrak{X}}$, such that

$\displaystyle\begin{split}&\left\lVert\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}}\right\rVert\leq
k_{0}<\infty\quad\text{and}\\\
&0<k_{1}<\lambda_{p}(\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Y}}})\lambda_{q}(\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Y}}})\leq\lambda_{1}(\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Y}}})\lambda_{1}(\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Y}}})\leq
k_{2}<\infty.\end{split}$ (46)

Since at least $\lfloor\nicefrac{{(n+d+2)}}{{2}}\rfloor$ have a positive
weight and at most $\lfloor\nicefrac{{(n-d)}}{{2}}\rfloor-1$ observations are
replaced, there are at least $d+2$ observations of the orginal sample
$\bm{\mathfrak{X}}$ contained in $\bm{\mathfrak{Y}}$ that have a positive
weight. Let $T\subseteq N$ denote the indices of those samples, then we have
that

$\displaystyle\sum_{i=1}^{n}w(d(\bm{Y}_{i}))=\sum_{i\in N\setminus
T}w(d(\bm{Y}_{i}))+\sum_{i\in T}w(d(\bm{X}_{i}))\geq\sum_{i\in
T}w(d(\bm{X}_{i}))\geq(d+2)c_{0}>0,$ (47)

with $c_{0}:=\min_{i\in T}w(d(\bm{X}_{i}))>0$. This implies that the
denominators of $\tilde{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}}$,
$\tilde{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Y}}}$, and
$\tilde{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Y}}}$ are always
positive.

Let us now show that there exists a constant $\alpha_{0}<\infty$ only
dependent on $\bm{\mathfrak{X}}$ such that
$\left\lVert\tilde{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}}\right\rVert_{F}<\alpha_{0}$.
From Equation (38) we have that

$\displaystyle\left\lVert\bm{Y}_{i}-\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}}\right\rVert_{F}^{2}$
$\displaystyle\leq\operatorname{tr}(\hat{\bm{\Omega}}^{\operatorname{col}}_{\bm{\mathfrak{Y}}}(\bm{Y}_{i}-\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}})^{\prime}\hat{\bm{\Omega}}^{\operatorname{row}}_{\bm{\mathfrak{Y}}}(\bm{Y}_{i}-\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}}))\lambda_{1}(\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Y}}})\lambda_{1}(\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Y}}})$
$\displaystyle=d(\bm{Y}_{i})\lambda_{1}(\hat{\bm{\Sigma}}^{\operatorname{row}}_{\bm{\mathfrak{Y}}})\lambda_{1}(\hat{\bm{\Sigma}}^{\operatorname{col}}_{\bm{\mathfrak{Y}}}).$

When computing $\tilde{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}}$ we have
that $w(d(\bm{Y}_{i}))=0$ if $d(\bm{Y}_{i})>c_{1}$ and for all
$\bm{Y}_{i}\in\bm{\mathfrak{Y}}$ that are assigned positive weights, Equation
(46) yields

$\displaystyle\left\lVert\bm{Y}_{i}\right\rVert_{F}^{2}\leq\left\lVert\bm{Y}_{i}-\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}}\right\rVert_{F}^{2}+\left\lVert\hat{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}}\right\rVert_{F}^{2}\leq
c_{1}k_{2}+k_{0}^{2}.$ (48)

Since the denominator of is
$\tilde{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}}$ bounded according to
Equation (47), $w$ is non-increasing and bounded, and $k_{0}$ and $k_{2}$ are
only dependent on $\bm{\mathfrak{X}}$, there exsits a constant $\alpha_{0}$
only dependent on $\bm{\mathfrak{X}}$ such that

$\displaystyle\left\lVert\tilde{\bm{\operatorname{M}}}_{\bm{\mathfrak{Y}}}\right\rVert_{F}\leq\alpha_{0}<\infty.$
(49)

To show that the covariance does not break down, we first consider the case of
the weight function $w(d_{i})=\mathds{1}(d_{i}\leq c_{1})$, for $c_{1}>0$. Let
$S\subseteq\\{1,\dots,N\\}$ denote the subset of indices of the
$s=\left|S\right|$ samples of
$\bm{\mathfrak{Y}}=\\{\bm{Y}_{1},\dots,\bm{Y}_{n}\\}$ for which $d_{i}\leq
c_{1},i\in S$. Observe that the $h$ samples $\bm{Y}_{i}$, $i\in H$ are those
with the smallest MD, hence $T\subseteq H\subseteq S$. Let
$\hat{\bm{M}}_{\bm{\mathfrak{Y}}_{T}},\hat{\bm{\Sigma}}_{\bm{\mathfrak{Y}}_{T}},\hat{\bm{\Omega}}_{\bm{\mathfrak{Y}}_{T}}$
and
$\hat{\bm{M}}_{\bm{\mathfrak{Y}}_{S}},\hat{\bm{\Sigma}}_{\bm{\mathfrak{Y}}_{S}},\hat{\bm{\Omega}}_{\bm{\mathfrak{Y}}_{S}}$
denote the MLE estimators of $\bm{\mathfrak{Y}}_{T}=(\bm{Y}_{i})_{i\in T}$ and
$\bm{\mathfrak{Y}}_{S}=(\bm{Y}_{i})_{i\in S}$, respectively. Consider the
following three ellipsoids:

  * •

Let
$E_{T}=E(\hat{\bm{M}}_{\bm{\mathfrak{Y}}_{T}},\hat{\bm{\Sigma}}_{\bm{\mathfrak{Y}}_{T}},\hat{\bm{\Omega}}_{\bm{\mathfrak{Y}}_{T}})$
denote the ellipsoid corresponding to the MLEs of $\bm{\mathfrak{Y}}_{T}$,
i.e., the smallest ellipsoid containing those at least $d+2$ samples.

  * •

Let
$E_{S}=E(\hat{\bm{M}}_{\bm{\mathfrak{Y}}_{S}},\hat{\bm{\Sigma}}_{\bm{\mathfrak{Y}}_{S}},\hat{\bm{\Omega}}_{\bm{\mathfrak{Y}}_{S}})$
denote the ellipsoid corresponding to the MLEs of $\bm{\mathfrak{Y}}_{S}$.

  * •

Let $E_{0}=E(\bm{0},k\bm{I}_{p},\bm{I}_{q})$ denote the smallest sphere
containing the samples $\bm{\mathfrak{Y}}_{S}$, where $k=c_{1}k_{2}+k_{0}^{2}$
is as in (48).

Observe first that as $E_{T}$ is the smallest ellipsoid containing the samples
$\bm{\mathfrak{Y}}_{T}=\bm{\mathfrak{X}}_{T}$ that are also in $E_{S}$, there
exists a constant $a_{1}$ depending only on $\bm{\mathfrak{X}}_{T}$, such that
$E_{T}\subseteq
a_{1}E_{S}:=E(\hat{\bm{M}}_{\bm{\mathfrak{Y}}_{S}},a_{1}\hat{\bm{\Sigma}}_{\bm{\mathfrak{Y}}_{S}},\hat{\bm{\Omega}}_{\bm{\mathfrak{Y}}_{S}})$.
On the other hand, $E_{S}$ is the smallest ellipsoid containing
$\bm{\mathfrak{Y}}_{S}$. As these points are also in $E_{0}$, then there exist
$\alpha=\alpha(c_{1})$ such that $\det(E_{S})\leq\det(E_{0})\leq\alpha$.
Equivalent argumentation as in the proof of Theorem 3.0.1 completes the first
part of the proof.

Let now $w=w(d_{i})$ be an arbitrarily, nondecreasing, bounded weight
function, such that $w(d_{i})=0$ if $d_{i}>c_{1}$, $i=1,\dots,n$. The weighted
log-likelihood function for the sample $\bm{\mathfrak{Y}}$, with the weights
satisfying $\sum_{i=1}^{n}w_{i}=s$ is given by

$\displaystyle
l(\bm{w},\bm{\operatorname{M}},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}}|\bm{\mathfrak{Y}})=$
$\displaystyle-\frac{1}{2}\sum_{i=1}^{m}w_{i}\Bigl{(}p\ln(\det(\bm{\Sigma}^{\operatorname{col}}))+q\ln(\det(\bm{\Sigma}^{\operatorname{row}}))$
$\displaystyle+\operatorname{tr}(\bm{\Omega}^{\operatorname{col}}(\bm{Y}_{i}-\bm{\operatorname{M}})^{\prime}\bm{\Omega}^{\operatorname{row}}(\bm{Y}_{i}-\bm{\operatorname{M}}))+pq\ln(2\pi)\Bigr{)}$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\Bigl{(}s\bigl{(}p\ln(\det(\bm{\Sigma}^{\operatorname{col}}))+q\ln(\det(\bm{\Sigma}^{\operatorname{row}}))\bigr{)}$
$\displaystyle+\sum_{i=1}^{s}\operatorname{tr}(\bm{\Omega}^{\operatorname{col}}(\bm{Z}_{i}-\bm{\operatorname{M}})^{\prime}\bm{\Omega}^{\operatorname{row}}(\bm{Z}_{i}-\bm{\operatorname{M}}))+pq\ln(2\pi)\Bigr{)}$
$\displaystyle=$
$\displaystyle\,l(\tilde{\bm{w}},\bm{M},\bm{\Sigma}^{\operatorname{row}},\bm{\Sigma}^{\operatorname{col}}|\bm{\mathfrak{Z}}),$

where $\tilde{\bm{w}}=(\tilde{w}(d_{1}),\dots,\tilde{w}(d_{n}))$, the new
weight function satisfies $\tilde{w}(d_{i})=\mathds{1}(d_{1}\leq c_{1})$,
$\bm{\mathfrak{Z}}=\\{\bm{Z}_{1},\dots,\bm{Z}_{n}\\}$, and
$\bm{Z}_{i}=\sqrt{w_{i}}\bm{Y}_{i}$, $i=1,\dots,n$. To complete the proof it
is sufficient to observe the following: $\bm{Z}_{1},\dots,\bm{Z}_{h}$ contains
at least $d+2$ points of the form $\sqrt{w_{i}}\bm{X}_{i}$ and are in a
general position, as $w_{i}\geq a_{2}>0$, for some constant depending only on
$\bm{\mathfrak{X}}$. Moreover,
$\|\bm{Z}_{i}\|_{F}^{2}=w_{i}\|\bm{Y}_{i}\|_{F}^{2}\leq
w_{i}(c_{1}k_{2}+k_{0}^{2})\leq w(0)(c_{1}k_{2}+k_{0}^{2})$, $i=1,\dots,s$.
The statement now follows from the first part of the proof, observing that
assumption $\sum_{i=1}^{m}w_{i}=s$ without loss of generality, since
$0<w(0)\leq\sum_{i=1}^{n}w_{i}\leq sw(0)<\infty$. ∎

##  Appendix C MMCD algorithm

Algorithm 1 Iterative C-step procedure for the MMCD estimators

1:procedure
CSTEP($(\bm{X}_{1},\ldots,\bm{X}_{n}),H_{\operatorname{old}},\varepsilon>0$)

2:
$(\hat{\bm{\operatorname{M}}}_{H_{\operatorname{new}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{H_{\operatorname{new}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{H_{\operatorname{new}}})=\text{MLE}((\bm{X}_{i})_{i\in
H_{\operatorname{old}}})$

3: $h=\left|H_{\operatorname{old}}\right|$

4: repeat

5:
$(\hat{\bm{\operatorname{M}}}_{H_{\operatorname{old}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{H_{\operatorname{old}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{H_{\operatorname{old}}})=(\hat{\bm{\operatorname{M}}}_{H_{\operatorname{new}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{H_{\operatorname{new}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{H_{\operatorname{new}}})$

6:
$\bm{d}=(\operatorname{MMD}^{2}(\bm{X}_{1};\hat{\bm{\operatorname{M}}}_{H_{\operatorname{old}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{H_{\operatorname{old}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{H_{\operatorname{old}}}),\ldots,\operatorname{MMD}^{2}(\bm{X}_{n};\hat{\bm{\operatorname{M}}}_{H_{\operatorname{old}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{H_{\operatorname{old}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{H_{\operatorname{old}}}))$

7: $\pi(i)=\\{\\{1,\ldots,n\\}\rightarrow\\{1,\ldots,n\\}:i\mapsto
j:d_{\pi(1)}\leq\ldots\leq d_{\pi(n)}\\}$

8: $H_{\operatorname{new}}=\\{\pi(1),\pi(2),\ldots,\pi(h)\\}$

9:
$(\hat{\bm{\operatorname{M}}}_{H_{\operatorname{new}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{H_{\operatorname{new}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{H_{\operatorname{new}}})=\text{MLE}((\bm{X}_{i})_{i\in
H_{\operatorname{new}}})$

10: until
$\left|p(\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{col}}_{H_{\operatorname{old}}}))-\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{col}}_{H_{\operatorname{new}}})))+q(\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{row}}_{H_{\operatorname{old}}}))-\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{row}}_{H_{\operatorname{new}}})))\right|~{}<~{}\varepsilon$

11: return
$\hat{\bm{\operatorname{M}}}_{H_{\operatorname{new}}},\hat{\bm{\Sigma}}^{\operatorname{row}}_{H_{\operatorname{new}}},\hat{\bm{\Sigma}}^{\operatorname{col}}_{H_{\operatorname{new}}},\bm{d},H_{\operatorname{new}}$

12:end procedure

Algorithm 2 Fast _reweighted_ MMCD procedure

1:procedure MMCD($\bm{\mathfrak{X}}=(\bm{X}_{1},\ldots,\bm{X}_{n})$)

2: $h=\lfloor\nicefrac{{(n+d+2)}}{{2}}\rfloor$

3: $\alpha=\nicefrac{{h}}{{n}}$

4: $N=\\{1,\dots,n\\}$

5: for $k=1$ to $500$ do

6: $H_{k}$ = sample($N$, size = $d+2$)

7:
$(\hat{\bm{\operatorname{M}}}_{k},\hat{\bm{\Sigma}}^{\operatorname{row}}_{k},\hat{\bm{\Sigma}}^{\operatorname{col}}_{k},\bm{d}_{k},H_{k})$
= CSTEP2($\bm{\mathfrak{X}},H_{k}$) $\triangleright$ 2 MLE and C-step
iterations

8:
$\delta_{k}=p\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{col}}_{k}))+q\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{row}}_{k}))$

9: end for

10:
$\pi_{\delta}(i)=\\{\\{1,\ldots,500\\}\rightarrow\\{1,\ldots,500\\}:i\mapsto
j:\delta_{\pi_{\delta}(1)}\leq\ldots\leq\delta_{\pi_{\delta}(500)}\\}$

11: for $l\in\\{\pi_{\delta}(1),\pi_{\delta}(2),\dots,\pi_{\delta}(10)\\}$ do

12:
$(\hat{\bm{\operatorname{M}}}_{l},\hat{\bm{\Sigma}}^{\operatorname{row}}_{l},\hat{\bm{\Sigma}}^{\operatorname{col}}_{l},\bm{d}_{l},H_{l})$
= CSTEP($\bm{\mathfrak{X}},H_{l}$) $\triangleright$ Iterating C-steps until
convergence

13:
$\delta_{l}=p\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{col}}_{l}))+q\ln(\det(\hat{\bm{\Sigma}}^{\operatorname{row}}_{l}))$

14: end for

15: $j=\operatorname*{arg\,min}_{k\in N}(\delta_{k})$

16:
$(\hat{\bm{\operatorname{M}}},\hat{\bm{\Sigma}}^{\operatorname{row}},\hat{\bm{\Sigma}}^{\operatorname{col}})=(\hat{\bm{\operatorname{M}}}_{j},c(\alpha)\hat{\bm{\Sigma}}^{\operatorname{row}}_{j},\hat{\bm{\Sigma}}^{\operatorname{col}}_{j})$
$\triangleright$ Consistency scaling for raw MMCD

17:
$\bm{d}=(\operatorname{MMD}^{2}(\bm{X}_{1};\hat{\bm{\operatorname{M}}},\hat{\bm{\Sigma}}^{\operatorname{row}},\hat{\bm{\Sigma}}^{\operatorname{col}}),\ldots,\operatorname{MMD}^{2}(\bm{X}_{n};\hat{\bm{\operatorname{M}}},\hat{\bm{\Sigma}}^{\operatorname{row}},\hat{\bm{\Sigma}}^{\operatorname{col}}))$

18: $H=H_{j}\cup\\{i\in N|d_{i}<\chi^{2}_{0.975;pq}\\}$

19:
$(\hat{\bm{\operatorname{M}}},\hat{\bm{\Sigma}}^{\operatorname{row}},\hat{\bm{\Sigma}}^{\operatorname{col}})=\operatorname{MLE}(\bm{X}_{i\in