# The Influence Function of Graphical Lasso Estimators

Gaëtan Louvet Corresponding author. Email<EMAIL_ADDRESS>postal
address: Rue de Bruxelles, 61; 5000 Namur, Belgium Department of Mathematics
and Namur Institute for Complex Systems (Naxys), Université de Namur, Belgium
Jakob Raymaekers  Department of Quantitative Economics, Maastricht University,
The Netherlands  Germain Van Bever  Department of Mathematics and Namur
Institute for Complex Systems (Naxys), Université de Namur, Belgium  Ines
Wilms  Department of Quantitative Economics, Maastricht University, The
Netherlands

###### Abstract

The precision matrix that encodes conditional linear dependency relations
among a set of variables forms an important object of interest in multivariate
analysis. Sparse estimation procedures for precision matrices such as the
graphical lasso (Glasso) gained popularity as they facilitate
interpretability, thereby separating pairs of variables that are conditionally
dependent from those that are independent (given all other variables). Glasso
lacks, however, robustness to outliers. To overcome this problem, one
typically applies a robust plug-in procedure where the Glasso is computed from
a robust covariance estimate instead of the sample covariance, thereby
providing protection against outliers. In this paper, we study such estimators
theoretically, by deriving and comparing their influence function, sensitivity
curves and asymptotic variances.

Keywords: Asymptotic variance, Glasso, Gross-error sensitivity, Influence
function, Outliers, Robustness

##  1 Introduction

The estimation of precision (inverse covariance) matrices is indispensable in
multivariate analysis with applications ranging from economics and finance to
genomics, neuroscience and social networks (Fan et al., 2014). Let ${\bf
X}=(X_{1},\ldots,X_{p})^{\top}$ be a $p$-variate random vector with finite
second order moments, with mean $\boldsymbol{\mu}$ and covariance matrix
$\boldsymbol{\Sigma}$. We focus on estimating the precision matrix
$\boldsymbol{\Omega}:=\boldsymbol{\Sigma}^{-1}$ whose entries capture the
conditional linear dependencies between the components of ${\bf X}$.
Importantly, when ${\bf X}\sim N(\boldsymbol{\mu},\boldsymbol{\Sigma})$, the
variables $X_{i}$ and $X_{j}$ are conditionally independent given all other
variables if and only if $(\boldsymbol{\Omega})_{i,j}=0$. In that case,
recovering the sparsity pattern of the precision matrix is equivalent to
recovering the graph structure of the Gaussian graphical model $G=(V,E)$ where
the vertex set $V=\\{1,\ldots,p\\}$ and the edge set $E$ consists of the pairs
$(i,j)$ that are connected by an edge (thus having
$(\boldsymbol{\Omega})_{i,j}\neq 0$).

A common method for estimating sparse precision matrices is the graphical
lasso (Glasso) (Yuan and Lin, 2007; Banerjee et al., 2008; Rothman et al.,
2008; Friedman et al., 2008), defined as

$\widehat{\boldsymbol{\Omega}}=\operatorname*{arg\,min}_{\boldsymbol{\Omega}}\\{-\text{logdet}(\boldsymbol{\Omega})+\text{tr}({\bf
S}\boldsymbol{\Omega})+\lambda\|\boldsymbol{\Omega}^{-\text{diag}}\|_{1}\ \
\text{s.t.}\
\boldsymbol{\Omega}=\boldsymbol{\Omega}^{\top},\boldsymbol{\Omega}\succ 0\\},$
(1)

where $\text{logdet}(\boldsymbol{\Omega})=\log(|\boldsymbol{\Omega}|)$,
$\text{tr}(\cdot)$ denotes the trace, $\bf S$ is the sample covariance matrix,
$\cdot\succ 0$ denotes a positive definite matrix, and $\lambda>0$ is a
penalty parameter controlling the degree of sparsity. Throughout, we refer to
the estimator in equation (1), with the sample covariance matrix as input, as
the “standard Glasso”. By adding an $\ell_{1}$-penalty (in this case, the
$\ell_{1}$-norm of the off-diagonal elements of ${\bf\Omega}$) to the negative
log-likelihood of a sample of multivariate normal random variables, a sparse
precision matrix estimate is obtained thereby also permitting estimation in
high-dimensional settings with more parameters than observations.

While sparse precision matrix estimation forms an active area of research, far
fewer research has looked into the performance of such methods in the presence
of outliers. Nonetheless, in many multivariate settings, the occurrence of
(complex) outliers is to be expected. The standard Glasso, however, performs
poorly in presence of outliers as it relies on the sample covariance matrix
${\bf S}$ in equation (1), which is sensitive to even a single outlying data
point. Various robustification approaches of the Glasso have been proposed.
Finegold and Drton (2011) build upon the standard Glasso but use the
likelihood of the multivariate $t$-distribution instead of a normal one
thereby providing better protection against data contamination. Öllerer and
Croux (2015); Croux and Öllerer (2016); Tarr et al. (2016); Lafit et al.
(2022) consider high-dimensional precision matrix estimation under cellwise
contamination and propose to plug-in a pairwise robust covariance or
correlation estimator instead of the (non-robust) sample covariance into
equation (1). Applications of such robust sparse precision matrices to, for
instance, discriminant analysis are studied in Aerts and Wilms (2017).
However, most studies above do not provide a theoretical analysis on robust
Glasso estimators. Exceptions are Öllerer and Croux (2015); Croux and Öllerer
(2016) who derive a finite sample breakdown point of their robust Glasso under
cellwise contamination, and Loh and Tan (2018) who analyze the statistical
consistency of robust high-dimensional precision matrix estimators under
cellwise contamination.

In this paper, we study the robustness of sparse Glasso-based precision matrix
estimators theoretically by deriving their influence function. The influence
function is an important tool to measure the robustness of a statistical
functional (Hampel et al., 1986); it quantifies the influence of a small
amount of contamination placed at a given value $\boldsymbol{z}$ on a
statistical functional $T$, for data coming from a distribution $F$. More
precisely, the influence function of the statistical functional $T$ at the
distribution $F$ is, when it exists, given by

$\operatorname*{IF}(\boldsymbol{z};T,F)=\underset{\varepsilon\rightarrow
0}{\text{lim}}\frac{T(F_{\varepsilon,\boldsymbol{z}})-T(F)}{\varepsilon}=\frac{\partial}{\partial\varepsilon}T(F_{\varepsilon,\boldsymbol{z}})\Bigr{|}_{\varepsilon=0},$
(2)

where the contaminated distribution
$F_{\varepsilon,z}=(1-\varepsilon)F+\varepsilon\Delta(\boldsymbol{z})$ and
$\Delta(\boldsymbol{z})$ is the probability distribution that puts all its
mass at $\boldsymbol{z}$.

The remainder of the article is structured as follows. In Section 2, we
present the influence function of the Glasso for any plug-in scatter
functional, then discuss the special case of the standard Glasso and show its
influence function is unbounded. Section 3 discusses how the Glasso influence
function can be bounded and expressions of the gross-error sensitivity are
derived. We consider different robust Glasso estimators and compare them in
terms of their influence functions and sensitivity curves. In Section 4, we
compute the asymptotic variances and compare the robust Glasso estimators in
terms of their statistical efficiencies. Section 5 concludes. All proofs,
together with precise descriptions of the robust correlation measures used in
Section 3, are collected in the Appendix.

##  2 Influence Function of the Glasso

We first derive the general expression of the influence function of the Glasso
for any plug-in scatter functional (Section 2.1). Next, we discuss the case of
the standard Glasso and establish its lack of robustness (Section 2.2).

###  2.1 The Influence Function

We start with deriving the influence function of the associated functional
representation of the Glasso. The following notation will be used throughout.
For any $p\times p$ matrix $\boldsymbol{\Omega}$, let
$\boldsymbol{\omega}:=\operatorname*{vec}(\boldsymbol{\Omega})$ denote the
$p^{2}$-dimensional vector resulting from stacking the columns of
$\boldsymbol{\Omega}$. Let $\operatorname*{mat}$ denote the inverse operator
of $\operatorname*{vec}$, that is,
$\operatorname*{mat}(\operatorname*{vec}(\boldsymbol{\Omega}))=\boldsymbol{\Omega}$.
For a vector $\boldsymbol{\omega}\in\mathbb{R}^{p^{2}}$ with $s\leq p^{2}$
nonzero elements, let ${\bf D}_{{\omega}}$ denote a $p^{2}\times p^{2}$
symmetric permutation matrix putting the nonzero elements first. In
particular, $({\bf D}_{{\omega}}\boldsymbol{\omega})_{i}=0$ for all $i>s$ and
${\bf D}^{2}={\bf I}_{p}$, the $p\times p$ identity matrix. Let
$\|\boldsymbol{\omega}\|_{1}=\|\boldsymbol{\Omega}\|_{1}$ denote the usual
$\ell_{1}$-norm and let
$\|\boldsymbol{\omega}\|^{-\text{diag}}_{1}=\|\operatorname*{mat}(\boldsymbol{\omega})^{-\text{diag}}\|_{1}$.
Let $\otimes$ denote the Kronecker product, $C^{\infty}(\mathbb{R}^{p^{2}})$
the set of continuous and infinitely differentiable functions on
$\mathbb{R}^{p^{2}}$ and let $W^{2,2}(\mathbb{R}^{p^{2}})$ be the Sobolev
space of functions with square integrable second weak derivatives.

For any distribution $F$ on $\mathbb{R}^{p}$ and any _scatter functional_
${\bf S}(F)$, we derive the influence function of the matrix-valued functional

$\boldsymbol{T}_{\Omega}(F)=\operatorname*{arg\,min}_{\boldsymbol{\Omega}\in\mathbb{R}^{p\times
p}}\\{-\text{logdet}(\boldsymbol{\Omega})+\text{tr}({\bf
S}(F)\boldsymbol{\Omega})+\lambda\|\boldsymbol{\Omega}^{-\text{diag}}\|_{1}\ \
\text{s.t.}\
\boldsymbol{\Omega}=\boldsymbol{\Omega}^{\top},\boldsymbol{\Omega}\succ 0\\}.$
(3)

To do so, we rather provide the influence function of the equivalent vector-
valued version

$\boldsymbol{T}_{\omega}(F)=\operatorname*{arg\,min}_{\boldsymbol{\omega}\in\mathbb{R}^{p^{2}}}\\{-\text{logdet}(\operatorname*{mat}(\boldsymbol{\omega}))+(\operatorname*{vec}({\bf
S}(F)^{\top}))^{\top}\boldsymbol{\omega}+\lambda\|\boldsymbol{\omega}\|^{-\text{diag}}_{1}\\}.$
(4)

######  Remark 1.

The equivalence between both definitions stems from the fact that, for any
matrices ${\bf S}$ and $\boldsymbol{\Omega}$, $\operatorname*{tr}({\bf
S}\boldsymbol{\Omega})=(\operatorname*{vec}({\bf
S}^{\top}))^{\top}\operatorname*{vec}(\boldsymbol{\Omega})$. Moreover, in (4),
we do not impose $\boldsymbol{\omega}$ to be the vectorisation of a symmetric,
positive definite matrix. The solution to the minimisation problem, without
loss of generality, always respect those constraints. This comes from the fact
that, when looking at the first order conditions, the gradient stays symmetric
and that the logdet function acts as a barrier to the positive definite cone
(see Srinivasan and Panda, 2022, for example). Therefore, it ensures that the
solution (in matrix form) of the first order conditions is positive definite.
Hence, it holds that
$\boldsymbol{T}_{{\Omega}}(F)=\operatorname*{mat}(\boldsymbol{T}_{{\omega}}(F))$.

######  Remark 2.

The choice of the plug-in scatter functional ${\bf S}(F)$ will allow to
distinguish between robust and non-robust Glasso alternatives. For example,
the estimator provided in (1) is the empirical version of (3) obtained by
setting

${\bf S}(F)=\mathbb{E}_{F}\big{[}({\bf X}-\mathbb{E}_{F}[{\bf X}])({\bf
X}-\mathbb{E}_{F}[{\bf X}])^{\top}\big{]}\quad\textrm{and taking}\quad
F=F_{n}.$

Other choices of robust scatter functionals will be discussed in Section 3.

The theorem below gives the influence function of
$\boldsymbol{T}_{{\omega}}(F)$. The proof is provided in Appendix A.

######  Theorem 1.

Let $F$ be a distribution over $\mathbb{R}^{p}$ and let ${\bf S}(F)$ be a
scatter functional. Let $\boldsymbol{z}\in\mathbb{R}^{p}$. Assume that
$\boldsymbol{T}_{{\omega}}:=\boldsymbol{T}_{{\omega}}(F)$ has $s$ nonzero
components. Let $\boldsymbol{D}:=\boldsymbol{D}_{\boldsymbol{T}_{{\omega}}}$.
Then, the influence function of the graphical lasso satisfies
$\boldsymbol{D}\operatorname*{IF}(\boldsymbol{z};\boldsymbol{T}_{{\omega}},F)=$

$-\left(\begin{array}[]{c}\big{(}\boldsymbol{D}(\operatorname*{mat}(\boldsymbol{T}_{{\omega}})^{-1}\otimes\operatorname*{mat}(\boldsymbol{T}_{{\omega}})^{-1})\boldsymbol{D}\big{)}^{-1}_{1:s,1:s}(\boldsymbol{D}\operatorname*{IF}(\boldsymbol{z};\operatorname*{vec}({\bf
S}(F)),F))_{1:s}\\\ \bf{0}\end{array}\right).$

Theorem 1 above shows that the influence of $\boldsymbol{z}$ on null
components of $\boldsymbol{T}_{{\omega}}$ is $0$. This translates to the fact
that, under infinitesimal contamination, null values in the precision matrix
remain zero, so that the contamination does not affect the sparsity of the
solution obtained. This finding is similar to penalised regression models
where the influence function is zero for the null components of the parameter
vector (see, for example, Öllerer et al., 2015; Avella-Medina, 2017).

The influence of $\boldsymbol{z}$ on the non-zero components of
$\boldsymbol{T}_{{\omega}}$ depends on the influence of $\boldsymbol{z}$ on
the plug-in scatter functional. A somewhat surprising conclusion at first
sight is that the elements of the covariance matrix estimator corresponding to
the null elements of the precision matrix do not appear in the influence
function of the latter. This phenomenon originates from the fact that the
contamination only appears through the covariance matrix estimator in the
trace, more precisely, from the term

$\operatorname*{tr}({\bf
S}(F_{\varepsilon,\boldsymbol{z}})\boldsymbol{\Omega})=(\operatorname*{vec}({\bf
S}(F_{\varepsilon,\boldsymbol{z}})^{\top}))^{\top}\boldsymbol{\omega}.$

In turn, if
$(\boldsymbol{T}_{\Omega}(F))_{i,j}=(\boldsymbol{T}_{\Omega}(F_{\varepsilon,\boldsymbol{z}})))_{i,j}=0$,
putting contamination in the corresponding elements of ${\bf
S}(F_{\varepsilon,\boldsymbol{z}})$ will not impact the scalar product, nor
the estimation of the precision matrix.

The proof of Theorem 1 hinges on the implicit function theorem. However, due
to the non-differentiability of the penalty term, some care is required.
Similarly to Avella-Medina (2017) for regression problems, we circumvent the
problem by computing the influence function of the functional

$\boldsymbol{T}_{\omega,p_{m}}(F)=\operatorname*{arg\,min}_{\boldsymbol{\omega}\in\mathbb{R}^{p^{2}}}\\{-\text{logdet}(\operatorname*{mat}(\boldsymbol{\omega}))+(\operatorname*{vec}({\bf
S}(F)^{\top}))^{\top}\boldsymbol{\omega}+\lambda
p_{m}(\boldsymbol{\omega})\\},$ (5)

where $p_{m}\in C^{\infty}(\mathbb{R}^{p^{2}})$ is a convex and differentiable
norm which converges to $\|\cdot\|_{1}^{-\text{diag}}$ in
$W^{2,2}(\mathbb{R}^{p^{2}})$ as $m\rightarrow\infty$. Two examples of such
sequences are given by
$p_{m}(\boldsymbol{\omega})=\|\operatorname*{mat}(\boldsymbol{\omega})^{-\text{diag}}\|_{k_{m}}$
or $p^{\prime}_{m}(\boldsymbol{\omega})=\sum_{i\neq
j}|(\operatorname*{mat}(\boldsymbol{\omega}))_{i,j}|^{k_{m}}$, with
$1<k_{m}\in\mathbb{R}$ for all $m$ and $k_{m}\searrow 1$. We then define the
influence function of $\boldsymbol{T}_{\omega}$ as the limit for $m\to\infty$
of the influence function of $\boldsymbol{T}_{\omega,p_{m}}$. The result is
obtained provided the limit is independent of the chosen $p_{m}$ sequence,
which we state in the following lemma.

######  Lemma 1.

Let $F$ be a distribution over $\mathbb{R}^{p}$ and let ${\bf S}(F)$ be a
scatter functional. For any two sequences $p_{m}$ and $p^{\prime}_{m}$
converging to $\|\cdot\|_{1}^{-\text{diag}}$ in $W^{2,2}(\mathbb{R}^{p^{2}})$
and any $\boldsymbol{z}\in\mathbb{R}^{p}$, it holds that

$\lim_{m\rightarrow\infty}\big{(}IF(\boldsymbol{z};\boldsymbol{T}_{\omega,p_{m}},F)-IF(\boldsymbol{z};\boldsymbol{T}_{\omega,p^{\prime}_{m}},F)\big{)}=0.$

Moreover, if
$\sup_{\boldsymbol{z}\in\mathbb{R}^{p}}\|\operatorname*{vec}(IF(\boldsymbol{z};S(F),F)\|<\infty$,
that is, if the gross error sensitivity of $S(F)$ is bounded, then the
convergence holds true uniformly:

$\lim_{m\rightarrow\infty}\sup_{\boldsymbol{z}\in\mathbb{R}^{p}}\big{(}IF(\boldsymbol{z};\boldsymbol{T}_{\omega,p_{m}},F)-IF(\boldsymbol{z};\boldsymbol{T}_{\omega,p^{\prime}_{m}},F)\big{)}=0.$

###  2.2 Non-robustness of the Standard Glasso

Setting ${\bf S}(F)$ in equation (4) to be the classical covariance functional
yields the standard Glasso, which we now discuss. First, the standard Glasso
has an unbounded influence function, hence is not robust to outliers, due to
its reliance on the sample covariance. We illustrate its lack of robustness by
computing the influence function and sensitivity curve for a setting with
$p=3$ variables whose covariance matrix and corresponding sparse precision
matrix are respectively given by

$\boldsymbol{\Sigma}=\left(\begin{array}[]{ccc}1&1/2&1/4\\\ 1/2&1&1/2\\\
1/4&1/2&1\end{array}\right)\quad\text{and
}\quad\boldsymbol{\Omega}=\left(\begin{array}[]{ccc}4/3&-2/3&0\\\
-2/3&5/3&-2/3\\\ 0&-2/3&4/3\end{array}\right).$ (6)

We introduce contamination at point $\boldsymbol{z}=(z_{1},z_{2},0)$ with
values of $z_{1}$ and $z_{2}$ ranging from -6 to 6.

The left panel of Figure 1 displays
$\|{\operatorname*{IF}(\cdot;\boldsymbol{T}_{{\Omega}},F)}\|_{F}$, the
Frobenius norm of the influence function of the standard Glasso. The penalty
parameter is set to $\lambda=8\cdot 10^{-4}$ to ensure the correct sparsity
pattern is recovered. The influence function of the standard Glasso is
unbounded, hence it is not robust. This is to be expected as the influence
function of the standard Glasso depends on that of the classical covariance,
which is unbounded. The value of the norm increases fastest along the
direction $z_{1}=-z_{2}$.

Figure 1: Frobenius norm of the standard Glasso influence function (left) and
sensitivity curve (right) for contamination at $(z_{1},z_{2},0)$.

Next, we compute the sensitivity curve introduced by Tukey (1977). The
sensitivity curve can be thought of as a finite-sample equivalent of the
influence function, as in many cases the sensitivity curve converges to the
influence function when the sample size goes to infinity. Suppose we have a
sample
$\boldsymbol{X}_{n-1}=\\{\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{n-1}\\}$ of
size $n-1$ of $p$-variate observations. To this sample we add a variable point
$\boldsymbol{z}$ to obtain
$\boldsymbol{X}_{z}=\\{\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{n-1},\boldsymbol{z}\\}$.
We now define the sensitivity curve by

$\operatorname*{SC}{{}_{n}}({\boldsymbol{z}};\widehat{\boldsymbol{\Omega}};{\boldsymbol{X}}_{n-1})=n\cdot(\widehat{\boldsymbol{\Omega}}(\boldsymbol{X}_{z})-\widehat{\boldsymbol{\Omega}}(\boldsymbol{X}_{n-1})).$
(7)

It measures the change of the estimated precision matrix caused by adding the
contaminated point $\boldsymbol{z}$ to the clean sample, standardised by
$1/n$, the amount of contamination. The right panel of Figure 1 displays the
Frobenius norm of the sensitivity curve for the standard Glasso (again with
$\lambda=8\cdot 10^{-4}$) for $n=1000$. The results are averaged over 50
random samples of $n-1$ clean observations drawn from a multivariate normal
distribution with the covariance matrix from equation (6) and an additional
contaminated data point at location $\boldsymbol{z}=(z_{1},z_{2},0)$. The
sensitivity curve of the standard Glasso is unbounded, in line with the
obtained results of the influence function.

Next, we consider the gross-error sensitivity

$\gamma^{*}(\boldsymbol{T}_{{\Omega}},F)=\underset{\boldsymbol{z}\in\mathbb{R}^{p}}{\textnormal{sup}}\|{\operatorname*{IF}(\boldsymbol{z};\boldsymbol{T}_{{\Omega}},F)}\|_{F},$

which summarises the influence function in a single index by evaluating the
maximal influence an observation may have. In the particular case of the
standard Glasso, the gross-error sensitivity diverges since the influence
function is unbounded. It still remains of interest, however, to determine in
which direction the influence will increase the most, that is, to determine

$\gamma^{*}_{1}(\boldsymbol{T}_{{\Omega}},F):=\operatorname*{arg\,max}_{\boldsymbol{z}\in\mathcal{S}^{p-1}}\|{\operatorname*{IF}(\boldsymbol{z};\boldsymbol{T}_{{\Omega}},F)}\|_{F},$

where $\mathcal{S}^{p-1}$ denotes the hypersphere in $\mathbb{R}^{p}$. The
following lemma provides the result for the unpenalised case. The solution
turns out to be the eigenvector corresponding with either the smallest or the
largest eigenvalue, depending on a condition on the eigenvalues themselves.

######  Lemma 2.

Let $F$ be a distribution over $\mathbb{R}^{p}$ and let ${\bf S}(F)$ be the
classical covariance functional. Fix $\lambda=0$,
$g:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}:x\mapsto g(x)=x^{4}-2x^{3}$ and let
$\lambda_{1},\dots,\lambda_{p}$ (resp. ${\bf v}_{1},\dots,{\bf v}_{p}$) be the
eigenvalues (resp. eigenvectors) of $\boldsymbol{\Omega}$. It holds that

$\gamma^{*}_{1}(\boldsymbol{T}_{{\Omega}},F)=\left\\{\begin{array}[]{cl}{\bf
v}_{1}&\textrm{ if }g(\lambda_{1})\geq g(\lambda_{p})\\\ {\bf
v}_{p}&\textrm{otherwise}\end{array}\right.,$

and

$\max_{\boldsymbol{z}\in\mathcal{S}^{p-1}}\|{\operatorname*{IF}(\boldsymbol{z};\boldsymbol{T}_{{\Omega}},F)}\|_{F}=\sum_{i=1}^{p}\lambda_{i}+\max(g(\lambda_{1}),g(\lambda_{p})).$

Note that, in the penalised case, it is not possible to obtain an equivalent
result to Lemma 2 above. Indeed, in that situation, there is no guarantee that
$\boldsymbol{T}_{{\Omega}}=\boldsymbol{\Omega}$, so that the influence
function might depend on the precision matrix in more complex forms. It is,
however, easy to show that Lemma 2 also holds when $\lambda$ is taken large
enough so that $\boldsymbol{T}_{\Omega}$ is diagonal.

##  3 Influence Function of Robust Glasso Estimators

In Section 3.1 we discuss how the Glasso influence function can be bounded.
Section 3.2 then presents several plug-in covariance candidates to obtain a
robust Glasso and compares them in terms of their Glasso influence function
(Section 3.3) and sensitivity curve (Section 3.4).

###  3.1 Bounding the Influence Function via Robust Plug-in Estimators

In Section 2, we show that the influence function of the Glasso depends on the
influence function of the plug-in scatter functional. To bound the influence
function of a Glasso-type of estimator, one could plug in a robust covariance
estimator instead of the sample covariance in equation (1). In Section 3.2, we
discuss several candidates for robust covariance plug-ins. Such robust Glasso
estimators would indeed result in a bounded influence function since the
gross-error sensitivity can be bounded as

$\gamma^{*}(\boldsymbol{T}_{{\omega}},F)\leq\left\|(\boldsymbol{D}(\operatorname*{mat}(\boldsymbol{T}_{{\omega}})^{-1}\otimes\operatorname*{mat}(\boldsymbol{T}_{{\omega}})^{-1})\boldsymbol{D}\big{)}^{-1}_{1:s,1:s}\right\|\gamma^{*}((\boldsymbol{D}\operatorname*{vec}({\bf
S}))_{1:s},F),$

where the norm used is the operator norm. Note that, should
$(\boldsymbol{T}_{{\omega}}(F))_{i,j}\neq 0$ everywhere (implying
$\boldsymbol{D}=\mathbf{I}_{p}$), this bound simplifies to

$\gamma^{*}(\boldsymbol{T}_{{\omega}},F)\leq\|\boldsymbol{T}_{\Omega}\|^{2}\gamma^{*}((\operatorname*{vec}({\bf
S}))_{1:s},F).$

Before turning to the discussion of the robust covariance plug-ins, note that
the Glasso is not orthogonal equivariant for $\lambda>0$, and this regardless
of the initial covariance estimator used. Hence, the numerical results on the
Glasso influence function and sensitivity curve presented in Sections 3.3 and
3.4 are dependent on the choice of covariance $\boldsymbol{\Sigma}$. This lack
of equivariance is, however, a price one needs to pay for obtaining a
regularised estimate of the precision matrix. Nonetheless, we believe that the
results presented below provide valuable insights into the trade-off between
robustness and statistical efficiency of the various Glasso estimators.

###  3.2 Robust Covariance Plug-ins

As candidates for the robust covariance estimator to be plugged into the
Glasso objective function, we consider the popular, traditional row-wise
robust Minimum Covariance Determinant (MCD) estimator (Rousseeuw, 1984, 1985;
Rousseeuw and Van Driessen, 1999) as well as a robust covariance estimator
based on pairwise robust correlation estimators (see e.g. Öllerer and Croux,
2015; Croux and Öllerer, 2016). Unlike the MCD, the latter are frequently used
in the literature on Glasso robustifications as they are available in high-
dimensional settings and protect against cellwise outliers (Alqallaf et al.,
2009).

The influence function of the Glasso based on the MCD can simply be obtained
by plugging in the expression of the MCD influence function (see Croux and
Haesbroeck, 1999) into the expression of the Glasso influence function
provided in Theorem 1. To arrive at the Glasso influence function based on a
pairwise robust correlation estimator, we need to introduce additional
notation related to the latter’s functional representation.

Denote with $S$ a scale functional and $R$ a correlation functional. The
definition of the covariance matrix functional ${\bf T}_{\Sigma}$ based on
pairwise robust correlation estimators is given by

$({\bf T}_{\Sigma}(F))_{j,k}=S(F_{j})S(F_{k})R(F_{j,k}),$ (8)

where $F_{j}$ and $F_{k}$ denote the marginal distributions of the $j$th and
$k$th components respectively, and $F_{j,k}$ denotes the bivariate (marginal)
distribution of the $j$th and $k$th components.

As a scale estimator, we use the $Q_{n}$ estimator (Rousseeuw and Croux,
1993), whose functional version is given by
$S(F)=cH^{-1}\left(\frac{1}{4}\right)$, where $c$ is a constant ensuring
Fisher-consistency and $H$ is the distribution function of the absolute
difference of two independent random variables, each with the same
distribution $F$.

As robust correlation estimator, we consider Spearman’s rank correlation
(Spearman, 1904), Kendall’s correlation coefficient (Kendall, 1938), the
Gaussian rank correlation (Boudt et al., 2012) and the Quadrant correlation
(Blomqvist, 1950). As our main focus is the influence function, we will
consider the functional version of these estimators below. For the sake of
completeness, we give the more common finite-sample definitions of these
estimators in Appendix B.

The functional form of Spearman’s rank correlation coefficient $R_{S}$ is
given by

$R_{S}(F_{j,k})=12\operatorname{\mathbb{E}}_{F_{j,k}}\left[F_{j}\left(X_{j}\right)F_{k}\left(X_{k}\right)\right]-3,$

where $(X_{j},X_{k})\sim F_{j,k}$ and $F_{j}$ and $F_{k}$ are the marginal
distributions of $X_{j}$ and $X_{k}$ as before. Kendall’s correlation
coefficient $R_{K}$ is given in functional form by

$R_{K}(F_{j,k})=\operatorname{\mathbb{E}}_{F_{j,k}}\left[\operatorname{sign}((X_{j_{1}}-X_{j_{2}})(X_{k_{1}}-X_{k_{2}}))\right],$

where $(X_{j_{1}},X_{k_{1}})$ and $(X_{j_{2}},X_{k_{2}})$ are independent
bivariate random variables, each with distribution $F_{j,k}$. The Gaussian
rank correlation coefficient is given in functional form by

$R_{G}(F_{j,k})=\operatorname{\mathbb{E}}_{F_{j,k}}\left[\Phi^{-1}\left(F_{j}\left(X_{j}\right)\right)\Phi^{-1}\left(F_{k}\left(X_{k}\right)\right)\right],$

where $\Phi^{-1}$ denotes the quantile function of the standard normal
distribution. Note that if $F_{j,k}$ is bivariate normal, we obtain the
functional of the Pearson correlation coefficient. The Quadrant correlation
coefficient is given in functional form by

$R_{Q}(F_{j,k})=\operatorname{\mathbb{E}}_{F_{j,k}}\left[\operatorname{sign}((X_{j}-\operatorname{median}(X_{j}))(X_{k}-\operatorname{median}(X_{k})))\right].$

For comparability of the influence functions of the Glasso based on the
different plug-in covariances, it is important that all estimate the same
population quantity, i.e. are Fisher consistent. Hence, we use the Fisher
consistent versions of the Spearman, Kendall and Quadrant correlation at the
normal model by applying the non-linear transformations as given, for
instance, in Croux and Dehon (2010). In particular, we use
$\tilde{R}_{S}=2\sin\left(\frac{\pi}{6}R_{S}\right)$,
$\tilde{R}_{K}=\sin\left(\frac{\pi}{2}R_{K}\right)$ and
$\tilde{R}_{Q}=\sin\left(\frac{\pi}{2}R_{Q}\right)$. No transformation for the
Gaussian rank correlation is needed as it is Fisher consistent at the normal
model. Note that in finite samples, these transformations can destroy the
positive semidefiniteness (PSD) of the resulting covariance matrix. This is
typically dealt with by projecting the estimate onto the space of PSD
matrices, for example with respect to the Frobenius norm (as in Öllerer and
Croux, 2015) or the elementwise $\ell_{\infty}$ norm (as in Loh and Tan,
2018). This does not affect the discussion on the influence function.

The influence function of the resulting covariance estimator is then given by

$\displaystyle\operatorname*{IF}(\boldsymbol{z};(\boldsymbol{T}_{\Sigma}(F))_{j,k},F)=$
$\displaystyle\operatorname*{IF}(z_{j};S,F_{j})S(F_{k})R(F_{j,k})$
$\displaystyle+$ $\displaystyle
S(F_{j})\operatorname*{IF}(z_{k};S,F_{k})R(F_{j,k})$ $\displaystyle+$
$\displaystyle S(F_{j})S(F_{k})\operatorname*{IF}(\boldsymbol{z};R,F_{j,k}).$

To compute this influence function, the expression of the influence function
of the scale estimator, namely the $Q_{n}$, and the considered correlation
estimators are thus needed. The former is available in Rousseeuw and Croux
(1993), the latter in Croux and Dehon (2010) for Kendall, Spearman and
Quadrant, and Boudt et al. (2012) for the Gaussian rank.

In the next subsections, we present the influence functions and sensitivity
curves of the Glasso estimators based on these different robust covariances.

###  3.3 Influence Functions

We consider the same example as in Section 2.2 and obtain the influence
function of the Glasso (with the same fixed penalty parameter as before)
computed from five different robust covariance plug-ins, namely the MCD, and
the pairwise correlation estimators based on Spearman, Kendall, Gaussian rank
and Quadrant.

Classical

MCD

Gaussian rank

Kendall

Spearman

Quadrant

Figure 2: Norm of the Glasso influence function based on several covariance
matrix estimators for contamination at $(z_{1},z_{2},0)$.

Figure 2 displays the norm of the Glasso influence functions for the various
covariance estimators, and includes the one for the classical covariance (from
Figure 1) to facilitate comparisons. When using the MCD as initial estimator,
we see that the norm of the influence function resembles the one for the
standard Glasso around the center, while the influence becomes bounded for
points further away from it. The Glasso influence function based on the MCD as
initial estimator is thus B-robust (see Hampel et al., 1986). However, while
being highly robust, it should be noted that it is computationally demanding
to obtain and, in practice, only available when $p$ is small relative to $n$.

Next, we turn to the Glasso influence function based on the Gaussian rank. At
first sight, it largely resembles the results for the classical covariance, as
its influence function is also unbounded. This is not surprising, as the
influence function of the Gaussian rank correlation is unbounded (Boudt et
al., 2012). Nonetheless, the values of the norm are smaller, thereby
indicating a somewhat larger resistance to single outliers than the classical
covariance. This is most likely due to the use of $Q_{n}$ as an estimator of
the marginal scales. This unboundedness of the influence function is somewhat
misleading however, and our results on the sensitivity curve (to be discussed
in Section 3.4) indicate that the Gaussian rank yields a substantially more
robust estimator than the standard Glasso in finite samples.

The Glasso influence functions based on Kendall and Spearman are very similar:
they are bounded and smooth, hence the Glasso computed from these initial
covariance matrices is robust to outliers. Finally, the Glasso influence
function based on the Quadrant is also bounded but not smooth: it has jumps at
the coordinate axes due to the occurrence of the sign function and the median
in its definition. Hence, small changes in data points close to the median of
one of the marginals may lead to relatively large changes in its Glasso.
Nonetheless, usage of the Quadrant as initial estimator results in the overall
smallest values of the norm of the Glasso influence function. This is in line
with findings in the setting of correlation estimation where the Quadrant is
known to be the “most B-robust”, i.e. it has the lowest gross-error
sensitivity among a class of correlation estimators based on product moments
(hence not including the MCD) (Raymaekers and Rousseeuw, 2021).

We summarise the performance of the Glasso estimators based on the different
initial covariances in Figure 3. We introduce contamination in the direction
of the eigenvector corresponding to the smallest eigenvalue of
$\boldsymbol{\Sigma}$, which, according to Lemma 2, maximises the norm of the
influence function of the classical covariance. Figure 3 displays the norm of
the Glasso influence functions for increasing values of the outlier’s norm. In
line with our previous discussion, we see that the influence function for the
standard Glasso and the Glasso based on the Gaussian rank are unbounded. The
Glasso influence function based on the MCD, Quadrant, Spearman and Kendall are
all bounded, with the former two being most robust for more severe outliers.

Figure 3: Norm of the Glasso influence function based on several covariance
matrix estimators for contamination in the direction of the eigenvector
related to the smallest eigenvalue of $\boldsymbol{\Sigma}$.

###  3.4 Sensitivity Curves

We compute the Glasso sensitivity curves based on pairwise correlation
estimators computed from the Gaussian rank, Spearman and Kendall correlation,
analogously to the procedure described in Section 2.2 for the sample
covariance111We omit results for the MCD and Quadrant as the former is
computationally demanding and the latter’s discontinuity leads to unstable
behavior in finite samples.. Results are available in Figures 4 and 5 of
Appendix C.

The Glasso sensitivity curves for the Spearman and Kendall correlation are
bounded and alike, in line with the results of Section 3.3. The Glasso
computed from the Gaussian rank is considerably more robust than the standard
Glasso in finite samples. Nonetheless, its influence function does not reflect
this as it is unbounded (see Figure 2). The same was established for the
Gaussian rank correlation estimator by Boudt et al. (2012). Still, the Glasso
computed from the Spearman and Kendall correlations are considerably more
robust.

##  4 Asymptotic Variances

Apart from robustness, it is desirable that an estimator has a high
statistical efficiency. To measure the statistical efficiency of the different
Glasso estimators, we first provide the expressions of their asymptotic
variances, then compare the asymptotic variances to the standard Glasso.

Should a Bahadur representation of $\boldsymbol{T}_{\Omega}(F)$ exist, one
would then conclude to asymptotic normality and Fréchet differentiability of
the precision matrix estimator. The asymptotic variances would then be given
by

$\text{ASV}((\boldsymbol{T}_{\Omega})_{i,j},F)=\int_{\mathbb{R}^{p}}\text{IF}^{2}(\boldsymbol{z};(\boldsymbol{T}_{\Omega})_{i,j},F)dF(\boldsymbol{z}),$

for any $(i,j)$ such that $(\boldsymbol{T}_{\Omega}(F))_{i,j}\neq 0$. In
matrix form, plugging in the expression of the influence function derived in
Theorem 1 yields

$\text{ASV}((\boldsymbol{D}\boldsymbol{T}_{\omega})_{1:s},F)=\boldsymbol{A}(\boldsymbol{D}\text{ASV}(\boldsymbol{z};\operatorname*{vec}({\bf
S}(F)),F)\boldsymbol{D})_{1:s}\boldsymbol{A}^{\top},$ (9)

with

$\boldsymbol{A}=\big{(}\boldsymbol{D}(\operatorname*{mat}(\boldsymbol{T}_{{\omega}})^{-1}\otimes\operatorname*{mat}(\boldsymbol{T}_{{\omega}})^{-1})\boldsymbol{D}\big{)}^{-1}_{1:s,1:s}.$

A formal verification of the Bahadur representation seems to be challenging
and is still an open question (see, for example, Ren et al., 2015).

We now turn to the numerical example introduced in Section 2.2, and compute
the asymptotic variances for the Glasso based on the MCD, and pairwise
correlation estimators. Table 1 summarises the asymptotic efficiencies
relative to the standard Glasso for three different components— diagonal
elements (1,1), (2,2) and off-diagonal element (2,1) —of the precision matrix.
Results on the other components are similar and therefore omitted. The values
from Table 1 were obtained from (9) using numerical integration. The
asymptotic variances obtained above are consistent with simulation results
conducted by the authors.

The Glasso based on Gaussian rank returns, overall, the highest efficiencies
but is closely followed by the Glasso based on the Kendall and Spearman
correlations. Their asymptotic efficiencies remain above 80%. We suspect that
the efficiency of around 80% for the diagonal elements stems mainly from the
fact that we are using the $Q_{n}$ estimator to estimate the diagonal elements
of our initial covariance matrix. The $Q_{n}$ has an efficiency of 82% at
Gaussian data (Rousseeuw and Croux, 1993). The Glasso based on Quadrant
correlation and MCD, on the other hand, suffer from a severe loss of
statistical efficiencies with drops to around 30%-40% for the former and to
around 10% for the MCD.

Table 1: Asymptotic efficiencies of the Glasso estimator based on several robust covariance matrix estimators relative to the standard Glasso and this for three components of the precision matrix. |  | Glasso with covariance estimator |
---|---|---|---
Component |  | Gaussian rank | Kendall | Spearman | Quadrant | MCD
(1,1) |  | 0.8210 | 0.8150 | 0.8097 | 0.4866 | 0.1422
(2,2) |  | 0.8085 | 0.8091 | 0.8027 | 0.4187 | 0.1017
(2,1) |  | 0.9563 | 0.8725 | 0.8491 | 0.3004 | 0.1029

##  5 Conclusion

The standard Glasso is one of the most often used sparse estimators for
precision matrices yet it may seriously be affected by the presence of even a
single outlier. In this paper, we study the robustness of the Glasso by
deriving expressions of its influence function for any plug-in scatter
functional. We show that the influence function of the standard Glasso is
unbounded, thereby proving its lack of robustness. Nonetheless, the Glasso
influence function can be easily bounded by plugging in a robust covariance
matrix estimate instead of the sample covariance.

We consider several robust Glasso estimators relying on different covariance
matrix plug-ins (namely the MCD or pairwise correlation estimators based on
Spearman, Kendall, Gaussian rank or Quadrant), and study their trade-off in
robustness versus statistical efficiency. When plugging in a pairwise robust
correlation estimate based on Kendall or Spearman, the resulting Glasso
combines the attractive properties of (i) a bounded and smooth influence
function, and (ii) high statistical efficiency at the normal model. Opting for
the Glasso based on the Gaussian rank leads to a higher statistical
efficiency, better protection against single outliers than the standard Glasso
in finite samples but a price is paid in terms of robustness compared to the
Glasso based on Spearman or Kendall. Opting for the Glasso based on the MCD or
Quadrant correlation, on the other hand, results in a very strong resistance
against outliers but low statistical efficiency. Moreover, the MCD is
computationally demanding and only available for $p<n$, thereby limiting its
applicability to low-dimensional settings. An appropriate choice of initial
covariance matrix estimate to be plugged-into the Glasso thus depends on the
trade-off between robustness and statistical efficiency the practitioner is
willing to make.

## References

  * Aerts and Wilms (2017) Aerts, S. and Wilms, I. (2017), “Cellwise robust regularized discriminant analysis,” Statistical Analysis and Data Mining: The ASA Data Science Journal, 10, 436–447.
  * Alqallaf et al. (2009) Alqallaf, F.; Van Aelst, S.; Yohai, V. J. and Zamar, R. H. (2009), “Propagation of outliers in multivariate data,” The Annals of Statistics, 37, 311–331.
  * Avella-Medina (2016) Avella-Medina, M. (2016), “Robust penalized M-estimators for generalized linear and additive models,” Ph.D. thesis, Geneva University.
  * Avella-Medina (2017) — (2017), “Influence functions for penalized ${M}$-estimators,” Bernoulli, 23, 3178–3196.
  * Banerjee et al. (2008) Banerjee, O.; El Ghaoui, L. and d’Aspremont, A. (2008), “Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data,” Journal of Machine Learning Research, 9, 485–516.
  * Blomqvist (1950) Blomqvist, N. (1950), “On a measure of dependence between two random variables,” The Annals of Mathematical Statistics, 21, 593–600.
  * Boudt et al. (2012) Boudt, K.; Cornelissen, J. and Croux, C. (2012), “The Gaussian rank correlation estimator: Robustness properties,” Statistics and Computing, 22, 471–483.
  * Croux and Dehon (2010) Croux, C. and Dehon, C. (2010), “Influence functions of the Spearman and Kendall correlation measures,” Statistical methods & applications, 19, 497–515.
  * 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, 161–190.
  * Croux and Öllerer (2016) Croux, C. and Öllerer, V. (2016), “Robust and sparse estimation of the inverse covariance matrix using rank correlation measures,” in Recent Advances in Robust Statistics: Theory and Applications, Springer, pp. 35–55.
  * Fan et al. (2014) Fan, J.; Han, F. and Liu, H. (2014), “Challenges of big data analysis,” National Science Review, 1, 293–314.
  * Finegold and Drton (2011) Finegold, M. and Drton, M. (2011), “Robust graphical modeling of gene networks using classical and alternative $t$-distributions,” The Annals of Applied Statistics, 5, 1057–1080.
  * Friedman et al. (2008) Friedman, J.; Hastie, T. and Tibshirani, R. (2008), “Sparse inverse covariance estimation with the graphical lasso,” Biostatistics, 9, 432–441.
  * Hampel et al. (1986) Hampel, F. R.; Ronchetti, E. M.; Rousseeuw, P. J. and Stahel, W. A. (1986), Robust statistics: the approach based on influence functions, John Wiley & Sons.
  * Kendall (1938) Kendall, M. G. (1938), “A new measure of rank correlation,” Biometrika, 30, 81–93.
  * Kollo and von Rosen (2005) Kollo, T. and von Rosen, D. (2005), Advanced Multivariate Statistics with Matrices, Springer.
  * Lafit et al. (2022) Lafit, G.; Nogales, F.; Ruiz, M. and Zamar, R. (2022), “Robust graphical lasso based on multivariate Winsorization,” Working paper, arXiv:2201.03659, 1–33.
  * Loh and Tan (2018) Loh, P.-L. and Tan, X. L. (2018), “High-dimensional robust precision matrix estimation: Cellwise corruption under $\epsilon$-contamination,” Electronic Journal of Statistics, 12, 1429–1467.
  * Öllerer and Croux (2015) Öllerer, V. and Croux, C. (2015), “Robust high-dimensional precision matrix estimation,” in Modern nonparametric, robust and multivariate methods, Springer, pp. 325–350.
  * Öllerer et al. (2015) Öllerer, V.; Croux, C. and Alfons, A. (2015), “The influence function of penalized regression estimators,” Statistics, 49, 741–765.
  * Raymaekers and Rousseeuw (2021) Raymaekers, J. and Rousseeuw, P. J. (2021), “Fast robust correlation for high-dimensional data,” Technometrics, 63, 184–198.
  * Ren et al. (2015) Ren, Z.; Sun, T.; Zhang, C. and Zhou, H. (2015), “Asymptotic normality and optimalities in estimation of large gaussian graphical models,” The Annals of Statistics, 43, 991–1026.
  * Rothman et al. (2008) Rothman, A. J.; Bickel, P. J.; Levina, E. and Zhu, J. (2008), “Sparse permutation invariant covariance estimation,” Electronic Journal of Statistics, 2, 494–515.
  * Rousseeuw (1984) Rousseeuw, P. J. (1984), “Least median of squares regression,” Journal of the American statistical association, 79, 871–880.
  * Rousseeuw (1985) — (1985), “Multivariate estimation with high breakdown point,” Mathematical statistics and applications, 8, 37.
  * Rousseeuw and Croux (1993) Rousseeuw, P. J. and Croux, C. (1993), “Alternatives to the median absolute deviation,” Journal of the American Statistical association, 88, 1273–1283.
  * Rousseeuw and Van Driessen (1999) Rousseeuw, P. J. and Van Driessen, K. (1999), “A fast algorithm for the minimum covariance determinant estimator,” Technometrics, 41, 212–223.
  * Spearman (1904) Spearman, C. (1904), “General intelligence, objectively determined and measured,” The American Journal of Psychology, 15, 201–292.
  * Srinivasan and Panda (2022) Srinivasan, S. and Panda, N. (2022), “What is the gradient of a scalar function of a symmetric matrix?” Indian Journal of Pure and Applied Mathematics.
  * Tarr et al. (2016) Tarr, G.; Müller, S. and Weber, N. C. (2016), “Robust estimation of precision matrices under cellwise contamination,” Computational Statistics & Data Analysis, 93, 404–420.
  * Tukey (1977) Tukey, J. (1977), Exploratory Data Analysis, Massachusetts: Reading (Addison-Wesley).
  * Yuan and Lin (2007) Yuan, M. and Lin, Y. (2007), “Model selection and estimation in the Gaussian graphical model,” Biometrika, 94, 19–35.

##  Appendix A Proofs

We collect all proofs from the paper in this first appendix.

Proof of Theorem 1: Let $p_{m}$ be a sequence of convex norms in
$C^{\infty}(\mathbb{R}^{p^{2}})$ which converges to
$\|\cdot\|_{1}^{-\text{diag}}$ in $W^{2,2}(\mathbb{R}^{p^{2}})$ as
$m\rightarrow\infty$ and let $\boldsymbol{T}_{\omega,p_{m}}(F)$ be defined as
in (5). Without loss of generality, assume that the Hessian matrix
$\nabla^{2}p_{m}(\cdot)$ is everywhere diagonal222This is the case for a
relaxation proposed in Avella-Medina (2016), page 57-58. The proof of Lemma 1
below shows that the limit is independent of the chosen sequence.. Throughout,
we write $\boldsymbol{\Omega}=\operatorname*{mat}(\boldsymbol{\omega})$ and
$\boldsymbol{T}_{{\omega},p_{m}}(F)=\boldsymbol{T}_{{\omega},p_{m}}$ when
there is no ambiguity to the underlying distribution. Since the function to be
minimised in $\boldsymbol{T}_{\omega,p_{m}}$ is differentiable, the first
order condition writes

$-\operatorname*{vec}(\operatorname*{mat}(\boldsymbol{T}_{\omega,p_{m}}(F))^{-\top})+\operatorname*{vec}({\bf
S}(F))+\lambda\nabla p_{m}(\boldsymbol{T}_{\omega,p_{m}}(F))=0.$

Note that the first term comes from the fact that (see, Kollo and von Rosen
(2005), page 134)

$\frac{d}{d\omega_{i,j}}(-\text{logdet}(\operatorname*{mat}(\boldsymbol{\omega})))=-\frac{1}{\text{det}(\operatorname*{mat}(\boldsymbol{\omega}))}\frac{d}{d\omega_{i,j}}\text{det}(\operatorname*{mat}(\boldsymbol{\omega}))=-\big{(}\operatorname*{mat}(\boldsymbol{\omega})^{-1}\big{)}_{j,i},$

so that
$\nabla\big{(}-\text{logdet}(\operatorname*{mat}(\boldsymbol{\omega}))\big{)}=-\operatorname*{vec}(\operatorname*{mat}(\boldsymbol{\omega})^{-\top}).$
Define

$f(\varepsilon,\boldsymbol{\omega})=-\operatorname*{vec}(\operatorname*{mat}(\boldsymbol{\omega})^{-1})+\operatorname*{vec}({\bf
S}(F_{\varepsilon,\boldsymbol{z}}))+\lambda\nabla p_{m}(\boldsymbol{\omega}).$

In particular, it holds $f(0,\boldsymbol{T}_{\omega,p_{m}}(F))=0$. The partial
derivatives at $(0,\boldsymbol{T}_{\omega,p_{m}}(F))$ are given by

$\displaystyle\frac{df(\varepsilon,\boldsymbol{\omega})}{d\varepsilon}\Bigr{|}_{\varepsilon=0}$
$\displaystyle=$ $\displaystyle\operatorname*{vec}(\text{IF}({\bf z};{\bf
S}(F),F)\textrm{, \qquad and}$
$\displaystyle\frac{df(\varepsilon,\boldsymbol{\omega})}{d\boldsymbol{\omega}}\Bigr{|}_{\boldsymbol{\omega}={\boldsymbol{T}}_{\omega,p_{m}}}$
$\displaystyle=$
$\displaystyle\lambda\nabla^{2}p_{m}(T_{\boldsymbol{\omega},p_{m}})+\boldsymbol{U}(T_{\boldsymbol{\omega},p_{m}}),$

where (see Kollo and von Rosen (2005), pages 127-131)

$\boldsymbol{U}(\boldsymbol{\omega}):=\frac{d(-\operatorname*{vec}(\boldsymbol{\Omega}^{-1}))}{d\operatorname*{vec}(\boldsymbol{\Omega})}=-\frac{d(\boldsymbol{\Omega}^{-1})}{d\boldsymbol{\Omega}}=\boldsymbol{\Omega}^{-1}\otimes\boldsymbol{\Omega}^{-1}=\operatorname*{mat}(\boldsymbol{\omega})^{-1}\otimes\operatorname*{mat}(\boldsymbol{\omega})^{-1}.$

The implicit function theorem yields

$\text{IF}(\boldsymbol{z};\boldsymbol{T}_{{\omega},p_{m}},F)=-\big{(}\lambda\nabla^{2}p_{m}(\boldsymbol{T}_{\omega,p_{m}})+\boldsymbol{U}(\boldsymbol{T}_{\omega,p_{m}})\big{)}^{-1}\operatorname*{vec}(\text{IF}(\boldsymbol{z};{\bf
S}(F),F)).$

Note that this result applies here since the derivative with respect to
$\boldsymbol{\omega}$ is nowhere vanishing, a consequence of the positive
semi-definiteness of $\nabla^{2}p_{m}(\boldsymbol{\omega})$ and of the
positive definiteness of $\boldsymbol{U}(\boldsymbol{\omega})$ (see,
Srinivasan and Panda (2022)). Taking now the limit in $m$,
$\lim_{m\rightarrow\infty}\text{IF}(\boldsymbol{z};\boldsymbol{T}_{{\omega},p_{m}},F)$
rewrites

$\displaystyle\lim_{m\rightarrow\infty}-(\lambda\nabla^{2}p_{m}(\boldsymbol{T}_{{\omega},p_{m}})+\boldsymbol{U}(\boldsymbol{T}_{{\omega},p_{m}}))^{-1}\operatorname*{vec}(\operatorname*{IF}(\boldsymbol{z};{\bf
S}(F),F))$
$\displaystyle=\lim_{m\rightarrow\infty}-\boldsymbol{D}\boldsymbol{D}^{-1}(\lambda\nabla^{2}p_{m}(\boldsymbol{T}_{{\omega},p_{m}})+\boldsymbol{U}(\boldsymbol{T}_{{\omega},p_{m}}))^{-1}\boldsymbol{D}^{-\top}\boldsymbol{D}^{\top}\operatorname*{vec}(\operatorname*{IF}(\boldsymbol{z};{\bf
S}(F),F))$
$\displaystyle=\lim_{m\rightarrow\infty}-\boldsymbol{D}(\boldsymbol{D}^{\top}\lambda\nabla^{2}p_{m}(\boldsymbol{T}_{{\omega},p_{m}})\boldsymbol{D}+\boldsymbol{D}^{\top}\boldsymbol{U}(\boldsymbol{T}_{{\omega},p_{m}})\boldsymbol{D})^{-1}\boldsymbol{D}^{\top}\operatorname*{vec}(\operatorname*{IF}(\boldsymbol{z};{\bf
S}(F),F))$
$\displaystyle=\lim_{m\rightarrow\infty}-\boldsymbol{D}(\lambda\nabla^{2}p_{m}(\boldsymbol{D}\boldsymbol{T}_{{\omega},p_{m}})+\boldsymbol{D}^{\top}\boldsymbol{U}(\boldsymbol{T}_{{\omega},p_{m}})\boldsymbol{D})^{-1}\boldsymbol{D}^{\top}\operatorname*{vec}(\operatorname*{IF}(\boldsymbol{z};{\bf
S}(F),F))$
$\displaystyle=-\boldsymbol{D}\Big{(}\lim_{m\rightarrow\infty}(\lambda\nabla^{2}p_{m}(\boldsymbol{D}\boldsymbol{T}_{{\omega},p_{m}})+\boldsymbol{D}^{\top}\boldsymbol{U}(\boldsymbol{T}_{{\omega},p_{m}})\boldsymbol{D})^{-1}\Big{)}\boldsymbol{D}\operatorname*{vec}(\operatorname*{IF}(\boldsymbol{z};{\bf
S}(F),F)).$

We now show that the remaining limit converges to

$\left(\begin{array}[]{cc}(\boldsymbol{D}\operatorname*{mat}(\boldsymbol{T}_{{\omega}}(F))^{-1}\otimes\operatorname*{mat}(\boldsymbol{T}_{{\omega}}(F))^{-1})\boldsymbol{D})^{-1}_{1:s,1:s}&0\\\
0&0\end{array}\right),$

which allows to conclude (using ${\bf D}^{2}={\bf I}_{p}$). To show the
convergence, partition

$\lambda\nabla^{2}p_{m}(\boldsymbol{D}\boldsymbol{T}_{{\omega},p_{m}})+\boldsymbol{D}^{\top}\boldsymbol{U}(\boldsymbol{T}_{{\omega},p_{m}})\boldsymbol{D}=\binom{A_{11}\
A_{12}}{A_{21}\ A_{22}}$

into blocks of $s$ (resp. $(p^{2}-s)$) lines/columns. Provided all inverses
exist, the inverse is given by

$\left(\begin{array}[]{cc}A^{11}&A^{12}\\\
A^{21}&A^{22}\end{array}\right)=\left(\begin{array}[]{cc}(A_{11}-A_{12}A_{22}^{-1}A_{21})^{-1}&-A^{11}A_{12}A_{22}^{-1}\\\
-A_{22}^{-1}A_{21}A^{11}&(A_{22}-A_{21}A_{11}^{-1}A_{12})^{-1}\end{array}\right).$

Notice that, since the inverse function is continuous and the Hessian matrix
$\nabla^{2}p_{m}(\boldsymbol{T}_{{\omega},p_{m}})$ is diagonal, it holds
$(\nabla^{2}p_{m}(\boldsymbol{D}\boldsymbol{T}_{\omega}(F)))_{jj}\rightarrow
0$ for $j=1,...,s$ and
$(\nabla^{2}p_{m}(\boldsymbol{D}\boldsymbol{T}_{\omega}(F)))_{jj}$ diverges to
$+\infty$ for $j=s+1,...,p^{2}$. Therefore, $A_{22}^{-1}\rightarrow 0$,
$(A_{22}-A_{21}A_{11}^{-1}A_{12})^{-1}\rightarrow 0$ and
$(A_{11}-A_{12}A_{22}^{-1}A_{21})^{-1}\rightarrow(\boldsymbol{D}^{\top}\boldsymbol{U}(\boldsymbol{T}_{{\omega},p_{m}})\boldsymbol{D})^{-1}_{1:s,1:s},$
allowing to conclude to the claimed limit.

Note that we formally showed that
$\lim_{m\rightarrow\infty}\text{IF}(\boldsymbol{z};\boldsymbol{T}_{{\omega},p_{m}},F)$
converges to a limit that might be distinct from the influence function should
the derivative not exist. In our situation, however, the components of the
limiting result are either $0$ (showing that the component of the estimated
precision matrix stays at zero) or the derivative taken at a positive value
(which then is differentiable), which therefore concludes the proof. ∎

Proof of Lemma 1: We prove first that
$\lim_{m\rightarrow\infty}\boldsymbol{T}_{{\omega},p_{m}}(F)$ is unique. Using
Lemma 2 in Avella-Medina (2016) with
$\Lambda_{\lambda}(\boldsymbol{\omega};F,p)=-\text{logdet}(\boldsymbol{\omega})+\operatorname*{vec}({\bf
S})^{\top}\boldsymbol{\omega}+\lambda p_{m}(\boldsymbol{\omega})$ readily
gives

$\lim_{m\rightarrow\infty}\boldsymbol{T}_{{\omega},p_{m}}(F)=\boldsymbol{T}_{{\omega}}(F).$
(10)

We can now prove the independence of
$\lim_{m\rightarrow\infty}\operatorname*{IF}(\boldsymbol{z};\boldsymbol{T}_{{\omega},p_{m}},F)$
of $p_{m}$. Let $\\{p_{m}\\}_{m\geq 1}$ and $\\{p^{\prime}_{m}\\}_{m\geq 1}$
be two sequences in $C^{\infty}(\mathbb{R}^{p^{2}})$ converging to $p$ in
$W^{2,2}(\mathbb{R}^{p^{2}})$. We have

$\displaystyle\operatorname*{IF}(\boldsymbol{z};\boldsymbol{T}_{{\omega},p_{m}},F)-$
$\displaystyle\operatorname*{IF}(\boldsymbol{z};\boldsymbol{T}_{{\omega},p^{\prime}_{m}},F)=\big{(}-(\lambda\nabla^{2}p_{m}(\boldsymbol{T}_{{\omega},p_{m}}(F))+\boldsymbol{U}(\boldsymbol{T}_{{\omega},p_{m}}(F))^{-1}$
$\displaystyle+(\lambda\nabla^{2}p^{\prime}_{m}(\boldsymbol{T}_{{\omega},p^{\prime}_{m}}(F))+\boldsymbol{U}(\boldsymbol{T}_{{\omega},p^{\prime}_{m}}(F))^{-1}\big{)}\operatorname*{vec}(\operatorname*{IF}(\boldsymbol{z};S(F),F)).$

Since the inverse of a matrix is continuous and $\nabla^{2}p_{m}$ and
$\nabla^{2}p^{\prime}_{m}$ are sequences in $C^{\infty}(\mathbb{R}^{p^{2}})$,
we find, using (10),

$\lim_{m\rightarrow\infty}(\operatorname*{IF}(\boldsymbol{z};\boldsymbol{T}_{{\omega},p_{m}},F)-\operatorname*{IF}(\boldsymbol{z};\boldsymbol{T}_{{\omega},p^{\prime}_{m}},F))=0.$

Trivially, the convergence becomes uniform if
$\sup_{\boldsymbol{z}\in\mathbb{R}^{p}}\|\operatorname*{vec}(\operatorname*{IF}(\boldsymbol{z};S(F),F))\|<\infty$,
since it is the only term depending on $\boldsymbol{z}$. ∎

Proof of Lemma 2: It is easy to see that
$\operatorname*{IF}(\boldsymbol{z};{\bf
S},F)=\boldsymbol{z}\boldsymbol{z}^{\top}-\boldsymbol{\Sigma}$. Hence, since
$\textrm{tr}({\bf A}{\bf B}{\bf C})=\textrm{tr}({\bf B}{\bf C}{\bf A})$ and
using the spectral decomposition $\boldsymbol{\Omega}={\bf
O}\boldsymbol{\Lambda}{\bf O}^{\top}$, it follows

$\displaystyle\|\operatorname*{IF}(\boldsymbol{z};\boldsymbol{T}_{{\omega}},F)\|^{2}_{F}$
$\displaystyle=$
$\displaystyle\|\boldsymbol{\Omega}\operatorname*{IF}(\boldsymbol{z};{\bf
S},F)\boldsymbol{\Omega}\|^{2}_{F}=\|\boldsymbol{\Omega}(\boldsymbol{z}\boldsymbol{z}^{\top}-\boldsymbol{\Sigma})\boldsymbol{\Omega}\|^{2}_{F}=\|\boldsymbol{\Omega}^{2}(\boldsymbol{z}\boldsymbol{z}^{\top}-\boldsymbol{\Sigma})\|^{2}_{F}$
$\displaystyle=$
$\displaystyle\|\boldsymbol{\Omega}^{2}\boldsymbol{z}\boldsymbol{z}^{\top}-\boldsymbol{\Omega}\|^{2}_{F}=\|\boldsymbol{\Lambda}^{2}\boldsymbol{x}\boldsymbol{x}^{\top}-\boldsymbol{\Lambda}\|^{2}_{F}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{p}\lambda_{i}^{2}+\sum_{i=1}^{p}x_{i}^{2}(\lambda_{i}^{4}-2\lambda_{i}^{3}),$

with $\boldsymbol{x}=(x_{1},\dots,x_{p})^{\top}={\bf O}^{\top}\boldsymbol{z}$.
Since $x_{i}={\bf v}_{i}^{\top}\boldsymbol{z}$ and using
$\boldsymbol{z}=\sum_{j=1}^{p}\alpha_{j}{\bf v}_{j}$, we conclude

$\|\operatorname*{IF}(\boldsymbol{z};\boldsymbol{T}_{{\omega}},F)\|^{2}_{F}=\sum_{i=1}^{p}\lambda_{i}^{2}+\sum_{i=1}^{p}\alpha_{i}^{2}(\lambda_{i}^{4}-2\lambda_{i}^{3})\leq\sum_{i=1}^{p}\lambda_{i}^{2}+\max_{i\in\\{1,...,p\\}}(\lambda_{i}^{4}-2\lambda_{i}^{3}),$

which is reached for the eigenvector corresponding to the eigenvalue which
maximises $(\lambda^{4}-2\lambda^{3})$. Since $g$ is strictly decreasing on
the interval $[0,1.5]$ then strictly increasing, the influence function is
maximised at either $\lambda_{1}$ or $\lambda_{p}$, which completes the proof.
∎

##  Appendix B Robust Correlation Estimators

To make the paper self-contained, we hereby include an overview of the finite-
sample definitions of the different correlation estimators we consider.

First, Spearman’s rank correlation $r_{S}$ is the sample correlation of the
ranks of the observations. Second, Kendall’s correlation between variable $j$
and $k$ is given by

$r_{K}=\frac{2}{n(n-1)}\sum_{i<\ell}\textrm{sign}((x_{ij}-x_{\ell
j})(x_{ik}-x_{\ell k})).$

Third, the Gaussian rank correlation between variable $j$ and $k$ is the
sample correlation estimated from the Van Der Waerden scores of the data, as
given by

$r_{G}=\frac{\sum_{i=1}^{n}\Phi^{-1}\left(\dfrac{\text{rank}(x_{ij})}{n+1}\right)\Phi^{-1}\left(\dfrac{\text{rank}(x_{ik})}{n+1}\right)}{\sum_{i=1}^{n}\left(\Phi^{-1}\left(\dfrac{i}{n+1}\right)\right)^{2}},$

where $\text{rank}(x_{ij})$ denotes the rank of $x_{ij}$ among all components
of the $j$th variable. Finally, the Quadrant correlation between variable $j$
and $k$ is given by

$r_{Q}=\frac{1}{n}\sum_{i=1}^{n}\textrm{sign}\left\\{(x_{ij}-\textrm{med}_{\ell}(x_{\ell
j}))(x_{ik}-\textrm{med}_{\ell}(x_{\ell k}))\right\\},$

defined as the difference between, on the one hand, the frequency of the
centered observations in the first and third quadrant and, on the other hand,
the frequency of the centered observations in the second and fourth quadrant.

##  Appendix C Additional Figures

Classical

Gaussian rank

Kendall

Spearman

Figure 4: Norm of the Glasso sensitivity curve based on several covariance
matrix estimators for contamination at $(z_{1},z_{2},0)$. Figure 5: Norm of
the Glasso sensitivity curve based on several covariance matrix estimators for
contamination in the direction of the eigenvector related to the smallest
eigenvalue of $\boldsymbol{\Sigma}$.