# Variable selection and covariance structure identification using sparse
loadings

Jan O. Bauer<EMAIL_ADDRESS>Faculty of Business and Economics, University
of Basel
Peter Merian-Weg 6, 4002 Basel, Switzerland.  Baden-Wuerttemberg Cooperative
State University Mannheim, Coblitzallee 1-9, 68163 Mannheim, Germany

###### Abstract

We provide sparse principal loading analysis which is a new concept that
reduces dimensionality of cross sectional data and identifies the underlying
covariance structure. Sparse principal loading analysis selects a subset of
existing variables for dimensionality reduction while variables that have a
small distorting effect on the covariance matrix are discarded. Therefore, we
show how to detect these variables and provide methods to assess their
magnitude of distortion. Sparse principal loading analysis is twofold and can
also identify the underlying block diagonal covariance structure using sparse
loadings. This is a new approach in this context and we provide a required
criterion to evaluate if the found block-structure fits the sample. The method
uses sparse loadings rather than eigenvectors to decompose the covariance
matrix which can result in a large loss of information if the loadings of
choice are too sparse. However, we show that this is no concern in our new
concept because sparseness is controlled by the aforementioned evaluation
criterion. Further, we show the advantages of sparse principal loading
analysis both in the context of variable selection and covariance structure
detection, and illustrate the performance of the method with simulations and
on real datasets. Supplementary material for this article is available online.
_Keywords:_ Covariance Structure, Dimensionality Reduction, Principal Loading
Analysis, Sparse Principal Component Analysis, Sparsity, Variable Selection
_2020 Mathematics Subject Classification:_ 62H25, 62F07, 65F15, 65F50.

##  1 Introduction

Principal loading analysis (PLA) is a method developed by Bauer and Drabant
(2021) to reduce dimensionality using variable selection. The method chooses a
subset of variables based on the impact of linked eigenvectors on the
covariance matrix. PLA relies on hard-thresholding to link the eigenvectors to
variables. Bauer and Drabant (2021) and Bauer (2021) provide recommendations
for the threshold based on simulations in their work. However, since simulated
results have their limitations, it is natural to search for an extension that
does not rely on simulations. One approach by Bauer and Drabant (2023) is to
analyze PLA from a regression point of view to obtain threshold values in line
with coefficient significant tests of multivariate linear regression.

In this work, we provide the new concept of sparse principal loading analysis
(SPLA) which gives two major contributions. First, SPLA does not rely on
thresholding and therefore overcomes the concerns discussed above. The
proposed concept is based on sparse loadings rather than eigenvectors to make
the usage of a threshold redundant. Transforming a random vector using sparse
loadings yields variables that are potentially correlated. Therefore, we
propose the usage of a measure based on linear regression to correct
calculation of the explained variance. Second, SPLA identifies the block
covariance structure of the underlying random variables using sparse loadings.
This is a new approach in this context. We contribute a criterion to evaluate
the detected structure which is necessary to find the one that fits the
underlying block-structure.

Calculation of sparse loadings is not new and the respective literature is
broad. Jolliffe et al. (2003) proposed a method that constrains the loadings
using lasso regularization (see Tibshirani (1996)). Zou et al. (2006)
formulated the calculation of the eigenvectors as a regression problem with
quadratic penalization. The eigenvectors are then sparsed using lasso
regularization as well. Both regularizations together yield the elastic net
(see Zou and Hastie (2005)). Qi et al. (2013) consider regularization by an
elasticnet-type norm. On the other hand, Shen and Huang (2008) utilized the
connection between the eigendecomposition of the covariance matrix and the
singular value decomposition of the sample matrix to derive sparse loadings.
Their approach is based on a regularized singular value decomposition which
has also been used by Witten et al. (2009) for their method.

Our work is twofold since SPLA also identifies the underlying block-diagonal-
structure of the covariance matrix. There are recent parametric (see, e.g.,
Srivastava and Reid (2012) and Yamada et al. (2017)) and non-parametric works
(see, e.g., Pavlenko et al. (2012) and Devijver and Gallopin (2018))
especially in the high-dimensional context that have the same goal. The reason
is that various techniques exist to estimate covariance matrices that follow a
sparse structure (see, e.g., Ledoit and Wolf (2004), Karoui (2008), Bickel and
Levina (2008), and many more). Further, detection of the block-diagonal
structure prior to network inferences for Gaussian graphical models can
improve performance (Tan et al., 2015).

This article is organized as follows: Section 2 provides notation needed for
the remainder of this work. In Section 3, we motivate the new concept of SPLA.
We discuss how SPLA is used to detect an underlying block-diagonal covariance
structure in Section 4. In Section 5, we provide calculations of the explained
total variance of SPLA which is required for application. Sparseness of the
loadings to identify an underlying block-structure is essential for SPLA. We
provide a criterion to evaluate the found block covariance structure in
Section 6. Further, we contribute the algorithm to conduct SPLA in practice.
Although a sparser model is easier to interpret, it is problematic to work
with a model that is sparser than the underlying covariance matrix since this
results in too little explanatory power. Therefore, we also discuss this
variance-interpretability-trade-off and show respective solutions. A
simulation study to evaluate the performance for block-diagonal covariance
detection is given in Section 7. In Section 8, we provide real data examples
for variable selection and block covariance structure detection, and we give
an example for block covariance structure detection for synthetic data. SPLA
will also be compared to PLA. We take a resume and suggest extensions in
Section 9. Proofs and more detailed results are deferred to the supplementary
material.

##  2 Setup

We first state some notation in this section. Afterwards, we give three
assumptions used throughout this work and elaborate their respective
objective.

We use the $\mathcal{L}_{2}$ norm $\|\bm{v}\|_{2}$ and the
$\mathcal{L}_{\infty}$ norm $\|\bm{v}\|_{\infty}$ for column vectors. For any
matrix $\bm{A}$, $\|\bm{A}\|$ denotes the matrix norm and we will use the
spectral norm $\|\bm{A}\|_{2}$ and the infinity norm $\|\bm{A}\|_{\infty}$ as
induced norms, and the Frobenius norm $\|\bm{A}\|_{F}$. Further, ${\rm
tr}(\bm{A})$ denotes the trace and $\bm{A}^{+}$ denotes the Moore–Penrose
inverse of a matrix $\bm{A}$. $\mathcal{O}_{p}(\cdot)$ denotes stochastic
boundedness.

We consider a random vector
$\bm{X}=(X_{1},\ldots,X_{M})\equiv(\bm{X}_{\mathcal{K}},\bm{X}_{\mathcal{D}})\in\mathbb{R}^{M}$
consisting of two random vectors
$\bm{X}_{\mathcal{K}}^{\top}\in\mathbb{R}^{K}$ and
$\bm{X}_{\mathcal{D}}^{\top}\in\mathbb{R}^{D}$ such that $K+D=M$. The
covariance matrix of $\bm{X}$ is given by
$\tilde{\bm{\Sigma}}=(\tilde{\sigma}_{i,j})$ for $i,j\in\\{1,\ldots,M\\}$, and
we consider a perturbation matrix $\bm{\mathrm{E}}$ such that

$\tilde{\bm{\Sigma}}\equiv\begin{pmatrix}\tilde{\bm{\Sigma}}_{1}&\tilde{\bm{\Sigma}}_{12}\\\
\tilde{\bm{\Sigma}}_{12}^{\top}&\tilde{\bm{\Sigma}}_{2}\end{pmatrix}\equiv\begin{pmatrix}\bm{\Sigma}_{1}&\bm{0}\\\
\bm{0}&\bm{\Sigma}_{2}\end{pmatrix}+\begin{pmatrix}\bm{0}&\bm{\mathrm{E}}_{12}\\\
\bm{\mathrm{E}}_{12}^{\top}&\bm{0}\end{pmatrix}\equiv\bm{\Sigma}+\bm{\mathrm{E}}\;,$
(1)

can be written as a sum of matrices. The covariance matrices of
$\bm{X}_{\mathcal{K}}$ and $\bm{X}_{\mathcal{D}}$ are given by
$\tilde{\bm{\Sigma}}_{1}\equiv\bm{\Sigma}_{1}$ and
$\tilde{\bm{\Sigma}}_{2}\equiv\bm{\Sigma}_{2}$ respectively. $\bm{\Sigma}$ as
well as $\bm{\mathrm{E}}$ are technical constructions to facilitate the
distinction of the covariance matrices of the two random vectors
$\bm{X}_{\mathcal{K}}$ and $\bm{X}_{\mathcal{D}}$ from the covariances
$\bm{\mathrm{E}}$ between the two random vectors. The covariance matrix
follows the eigendecomposition
$\tilde{\bm{\Sigma}}=\bm{V}\bm{\Lambda}\bm{V}^{\top}$. $\bm{\Lambda}$ is
hereby a diagonal matrix containing the eigenvalues
$\lambda_{1}>\cdots>\lambda_{M}$ of $\tilde{\bm{\Sigma}}$ in descending order,
and $\bm{V}=(\bm{v}_{1},\ldots,\bm{v}_{M})$ is an orthonormal matrix
containing the respective eigenvectors where
$\bm{v}_{i}^{\top}=(v_{i}^{(1)},\ldots,v_{i}^{(M)})$ for all
$i\in\\{1,\ldots,M\\}$. Further, we will introduce sparse representations of
the eigenvectors throughout this work where
$\bm{U}=(\bm{u}_{1},\ldots,\bm{u}_{M})$ denotes the matrix containing these
sparse loadings.

The independent and identically distributed (IID) sample
$\bm{x}=(\bm{x}_{1},\ldots,\bm{x}_{M})\equiv(\bm{x}_{\mathcal{K}},\bm{x}_{\mathcal{D}})$
$\in\mathbb{R}^{N\times M}$ contains $N$ observations of the random vector
$\bm{X}$. In general, the sample counterparts are indicated using a hat. For
example,
$\hat{\tilde{\bm{\Sigma}}}\equiv\hat{\bm{\Sigma}}+\hat{\bm{\mathrm{E}}}$ is
the sample covariance matrix of $\bm{x}$.

Variables or blocks of variables are linked to sparse loadings in SPLA. If we
want to reduce dimensionality, we look for a set of variables $\mathcal{D}$ to
discard while keeping the remaining variables $\mathcal{K}$. If we want to
detect the underlying block covariance structure, our interest is to identify
if two disjoint blocks of variables $\mathcal{D}$ and $\mathcal{K}$ do exist.
Hence, we define $\mathcal{D}\subset\\{1,\ldots,M\\}$ with ${\rm
card}(\mathcal{D})=D$ such that $1\leq D<M$. $\mathcal{D}$ will contain the
indices of a variable-block $\\{X_{d}\\}\equiv\\{X_{d}\\}_{d\in\mathcal{D}}$,
and we introduce the quasi-complement
$\mathcal{K}\equiv\\{1,\ldots,M\\}\backslash\mathcal{D}$ with ${\rm
card}(\mathcal{K})=K$ to label the remaining block of variables
$\\{X_{k}\\}\equiv\\{X_{k}\\}_{k\in\mathcal{K}}$. In an analogous manner,
$\mathit{\Delta}$ with ${\rm card}(\mathit{\Delta})=D$ will be used to index
eigenvectors with respective eigenvalues linked to the $\\{X_{d}\\}$.

For convenience purposes, we consider only the two blocks $\\{X_{d}\\}$ and
$\\{X_{k}\\}$ throughout this work except for the examples in Section 8.
Further, since $\bm{X}\equiv(\bm{X}_{\mathcal{K}},\bm{X}_{\mathcal{D}})$, the
index sets are given by $\mathcal{D}=\\{K+1,\ldots,M\\}$ and
$\mathcal{K}=\\{1,\ldots,K\\}$ respectively. However, our results remain true
for the trivial extension to any number of blocks and for the general case if
$\mathcal{D}\subset\\{1,\ldots,M\\}$. All matrices following an analogous
block-structure as $\tilde{\bm{\Sigma}}$ given in (1) are then to be replaced
by their general form such as
$\tilde{\bm{\Sigma}}_{1}=\tilde{\bm{\Sigma}}_{[\mathcal{K},\mathcal{K}]}$,
$\tilde{\bm{\Sigma}}_{2}=\tilde{\bm{\Sigma}}_{[\mathcal{D},\mathcal{D}]}$,
$\tilde{\bm{\Sigma}}_{12}=\tilde{\bm{\Sigma}}_{[\mathcal{K},\mathcal{D}]}$,
and
$\tilde{\bm{\Sigma}}_{12}^{\top}=\tilde{\bm{\Sigma}}_{[\mathcal{D},\mathcal{K}]}$.
For our purposes, this matrix-notation is only needed for the proofs to which
we refer for a more elaborate explanation.

Given the notation above, the following assumptions are made:

######  Assumption 1.

$\mathbb{E}(\bm{X})=\bm{0}$.

We assume the random variables to have mean zero for convenience purposes to
simplify the covariance matrix to
$\mathbb{E}(\bm{X}^{\top}\bm{X})=\tilde{\bm{\Sigma}}$.

######  Assumption 2.

$\tilde{\bm{\Sigma}}$ is positive definite.

We want to extract the information contained in the $\\{X_{d}\\}$ that can be
explained by the $\\{X_{k}\\}$. This is done by regression projection where
$\tilde{\bm{\Sigma}}$ is required to be invertible. This corresponds to the
full rank assumption of $\bm{x}$ in the sample case. Since
$\tilde{\bm{\Sigma}}=\bm{\Sigma}$ if $\bm{\mathrm{E}}=\bm{0}$ by construction,
we implicitly assume that $\bm{\Sigma}$ is positive definite as well.

######  Assumption 3.

$\bm{U}$ is a block-diagonal matrix with $\bm{U}^{\top}\bm{U}=\bm{I}$.

Orthogonality of the sparse loadings is needed to obtain results regarding the
explained total variance in this work. For convenience purposes, we therefore
assume orthogonality throughout the whole work which holds without loss of
generality: if $\bm{U}$ was not orthogonal, we would orthogonalize the
loadings by using the singular value decomposition of $\bm{U}$ with singular
values replaced by ones as it is done in a Procrustes problem (see, e.g.,
Mardia et al. (1979)). This can be done due to the block-diagonal shape of
$\bm{U}$ which will be introduced in (3), and because the exact values of the
non-zero components of $\bm{U}$ are not of concern. Both will be discussed
elaborately in Section 3.

##  3 Methodology and Motivation

Firstly, we recap PLA in this section. The method relies on choosing a
threshold value which, however, can be difficult in practice. Hence, we
motivate the usage of SPLA based on sparse loadings which does not require
hard-thresholding anymore.

PLA is a concept for dimensionality reduction where a subset of existing
variables is selected while the other variables are discarded. We refer to
Bauer and Drabant (2021) for an elaborate explanation. However, the intuition
is that blocks of variables are discarded which distort the covariance matrix
only slightly. We scan for a block $\bm{X}_{\mathcal{D}}$ by checking if all
$K$ elements $k\in\mathcal{K}$ of $D$ eigenvectors
$\\{v_{\delta}^{(k)}\\}_{\delta\in\mathit{\Delta}}$ of the covariance matrix
$\tilde{\bm{\Sigma}}$ are smaller in absolute terms than a certain cut-off
value $\tau$. That is to say, we check if $|v_{\delta}^{(k)}|\leq\tau$ for all
$(k,\delta)\in\mathcal{K}\times\mathit{\Delta}\equiv\\{(k,\delta):k\in\mathcal{K},\delta\in\mathit{\Delta}\\}$.
For $\bm{X}=(\bm{X}_{\mathcal{K}},\bm{X}_{\mathcal{D}})$, we therefore scan if
the eigenvectors $\bm{V}$ are of shape

$\displaystyle\begin{split}\bm{v}_{\delta^{c}}&=\begin{pmatrix}\bm{\ast}_{K}\\\
\bm{\varepsilon}_{D}\end{pmatrix}\;\text{ for
}\delta^{c}\in\mathit{\Delta}^{c}\equiv\\{1,\ldots,M\\}\backslash\mathit{\Delta}\;\text{
with }|\mathit{\Delta}^{c}|=K\\\
\bm{v}_{\delta}&=\begin{pmatrix}\bm{\varepsilon}_{K}\\\
\bm{\ast}_{D}\end{pmatrix}\;\text{ for }\delta\in\mathit{\Delta}\;\text{ with
}|\mathit{\Delta}|=D\end{split}$ (2)

and hence represent the blocks $\bm{x}_{\mathcal{K}}$ and
$\bm{x}_{\mathcal{D}}$ respectively. $\bm{\ast}_{K}\in\mathbb{R}^{K}$,
$\bm{\ast}_{D}\in\mathbb{R}^{D}$ are hereby vectors containing components
larger than $\tau$ in absolute terms, and
$\bm{\varepsilon}_{K}\in\mathbb{R}^{K}$,
$\bm{\varepsilon}_{D}\in\mathbb{R}^{D}$ are vectors containing components that
are all smaller than $\tau$ in absolute terms respectively. The intuition
comes from the special case when $\tau=0$, i.e. when $\bm{\mathrm{E}}=\bm{0}$,
since then the $\\{\bm{v}_{\delta^{c}}\\}_{\delta^{c}\in\mathit{\Delta}^{c}}$
and $\\{\bm{v}_{\delta}\\}_{\delta\in\mathit{\Delta}}$ represent the exclusive
distortion of the covariance matrix by the $\\{X_{k}\\}$ and the $\\{X_{d}\\}$
respectively.

######  Lemma 1.

The eigenvectors $\bm{V}\equiv(\bm{v}_{1},\ldots,\bm{v}_{M})$ of the
covariance matrix $\tilde{\bm{\Sigma}}\equiv\bm{\Sigma}+\bm{\mathrm{E}}$ from
(1) are of shape $\bm{v}_{\delta^{c}}=(\bm{\ast}_{K},\bm{0}_{D})^{\top}$ and
$\bm{v}_{\delta}=(\bm{0}_{K},\bm{\ast}_{D})^{\top}$ if and only if
$\bm{\mathrm{E}}=\bm{0}$.

In practice, we are more likely to face the situation that
$\bm{\mathrm{E}}\neq\bm{0}$. However, if the perturbation $\bm{\mathrm{E}}$ is
small, then the components of $\bm{\varepsilon}_{K}$ and
$\bm{\varepsilon}_{D}$ are small as well because they are bounded by
$\bm{\mathrm{E}}$.

######  Corollary 1.

For $\bm{\varepsilon}_{K}\in\mathbb{R}^{K}$ and
$\bm{\varepsilon}_{D}\in\mathbb{R}^{D}$ from (2) it holds that
$\|\bm{\varepsilon}_{K}\|_{\infty}=\mathcal{O}_{p}(\|\bm{\mathrm{E}}\|_{l})$
and
$\|\bm{\varepsilon}_{D}\|_{\infty}=\mathcal{O}_{p}(\|\bm{\mathrm{E}}\|_{l})$
for $l\in\\{2,\infty\\}$.

Therefore, if the eigenvectors are of shape as in (2) with
$\bm{\varepsilon}_{K}$, $\bm{\varepsilon}_{D}$ having all components smaller
than $\tau$ in absolute terms, then most distortion of the variables
$\\{X_{d}\\}$ is geometrically represented by the eigenvectors
$\\{\bm{v}_{\delta}\\}_{\delta\in\mathit{\Delta}}$. The explained variance of
the $\\{X_{d}\\}$ can then be evaluated either by
$\sum_{d\in\mathcal{D}}\tilde{\sigma}_{d,d}$ or approximated by
$(\sum_{m}\lambda_{m})^{-1}(\sum_{\delta\in\mathit{\Delta}}\lambda_{\delta})$
where $\lambda_{i}$ is the eigenvalue corresponding to the eigenvector
$\bm{v}_{i}$. If the explained variance is sufficiently small for the
underlying purpose of application, the variables $\\{X_{d}\\}$ are discarded.
The same holds for the block $\\{X_{k}\\}$ in an analogous manner.

The concern in PLA is the proper choice of the threshold $\tau$ which is
elementary for the detection of the blocks $\\{X_{d}\\}$ and $\\{X_{k}\\}$
respectively. There are no theoretical derivations for the choice of the
threshold yet. Bauer and Drabant (2021) provided simulation results and Bauer
and Drabant (2023) built a connection between $\tau$ and hypothesis testing in
multivariate linear regression. However, if we use sparse loadings instead of
eigenvectors, we do not rely on a threshold anymore. This is the intuition of
the new concept of SPLA: Instead of checking if the eigenvectors follow a
shape according to (2), we rather construct sparse loadings of shape

$\displaystyle\begin{split}\bm{u}_{\delta^{c}}&=\begin{pmatrix}\bm{\ast}_{K}\\\
\bm{0}_{D}\end{pmatrix}\;\text{ for
}\delta^{c}\in\mathit{\Delta}^{c}\equiv\\{1,\ldots,M\\}\backslash\mathit{\Delta}\;\text{
with }|\mathit{\Delta}^{c}|=K\\\ \bm{u}_{\delta}&=\begin{pmatrix}\bm{0}_{K}\\\
\bm{\ast}_{D}\end{pmatrix}\;\text{ for }\delta\in\mathit{\Delta}\;\text{ with
}|\mathit{\Delta}|=D\end{split}$

which represent the underlying covariance matrix. While eigenvectors follow a
natural order due to the size of their respective eigenvalues, for our
purposes such an order is not required for the sparse loadings. Instead, we
will always consider $\bm{U}$ to be a matrix of block diagonal shape

$\bm{U}\equiv\begin{pmatrix}\bm{\ast}_{K\times K}&\bm{0}_{K\times D}\\\
\bm{0}_{D\times K}&\bm{\ast}_{D\times
D}\end{pmatrix}\equiv\begin{pmatrix}\bm{U}_{1}&\bm{0}\\\
\bm{0}&\bm{U}_{2}\end{pmatrix}\;,$ (3)

which reduces $\mathit{\Delta}^{c}=\\{1,\ldots,K\\}$ and
$\mathit{\Delta}=\\{K+1,\ldots,M\\}$. If the sparse loadings do not follow
this shape by construction, we will transform $\bm{U}$ to a matrix of block-
diagonal shape $\bm{P}_{1}\bm{U}\bm{P}_{2}$ using permutation matrices
$\bm{P}_{1}$ and $\bm{P}_{2}$ for row and column permutations respectively. We
note that the exact component-values in $\bm{\ast}_{K\times K}$ and
$\bm{\ast}_{D\times D}$ are not of concern, since we are only interested if an
underlying structure representing the blocks $\\{X_{k}\\}$ and $\\{X_{d}\\}$
does exist. Therefore, we can always orthogonalize $\bm{U}$ as described in
Assumption 3 despite of its changes of the component-values.

##  4 Block Covariance Structure Identification

In this section, we discuss how SPLA identifies if the random vector $\bm{X}$
follows a block covariance structure.

According to Lemma 1, it is reasonable to consider the eigenvectors for
identification of the shape of the covariance matrix. Further, it holds that
the covariance matrix follows a block-diagonal structure if all eigenvectors
are either close to $(\bm{\ast}_{K},\bm{0}_{D})^{\top}$ or close to
$(\bm{0}_{K},\bm{\ast}_{D})^{\top}$.

######  Lemma 2.

For $\bm{\varepsilon}_{K}\in\mathbb{R}^{K}$ and
$\bm{\varepsilon}_{D}\in\mathbb{R}^{D}$ from (2) it holds that
$\|\bm{\mathrm{E}}\|_{l}=\mathcal{O}_{p}(\|(\bm{\varepsilon}_{K},\bm{\varepsilon}_{D})^{\top}\|_{\infty})$
for $l\in\\{2,F,\infty\\}$.

Instead of checking if the eigenvectors are close to the aforementioned shape
in the $\mathcal{L}_{1}$ norm, we rather suggest to check if sparse loadings
exist that are not only of shape
$\bm{u}_{\delta^{c}}=(\bm{\ast}_{K},\bm{0}_{D})^{\top}$ and
$\bm{u}_{\delta}=(\bm{0}_{K},\bm{\ast}_{D})^{\top}$ but also represent the
underlying covariance matrix adequately. A criterion to evaluate if the sparse
loadings are representative for the covariance matrix and therefore for the
random vector are provided in Section 6.

Note that using sparse loadings instead of eigenvectors loses the bounds on
the perturbation. When using hard-thresholding on the eigenvectors in PLA, the
magnitude of $\bm{\mathrm{E}}$ relates to the identification as a block-
diagonal covariance matrix as follows:

######  Corollary 2.

Let $\tau$ be the threshold for PLA and $\lambda_{1}>\cdots>\lambda_{M}$ be
the eigenvalues of $\tilde{\bm{\Sigma}}$ with $\lambda_{0}\equiv\infty$ and
$\lambda_{M+1}\equiv-\infty$. If the covariance matrix does not follow a
block-structure in PLA-sense, then

$\|\bm{\mathrm{E}}\|_{l}>\tau\cdot
2^{-3/2}\min\limits_{j\in\\{1,\ldots,M+1\\}}(\lambda_{j-1}-\lambda_{j})\;,$

for $l\in\\{2,\infty\\}$.

This bound on $\bm{\mathrm{E}}$ does not hold for the concept of SPLA when
using sparse loadings however. The reason is that SPLA considers block
covariance structure identification based on the covariance matrix and on the
precision within a block, and not exclusively on the covariances between
blocks. Hence, SPLA takes into account the relation between $\bm{\mathrm{E}}$
and $\bm{\Sigma}_{1}^{-1}$. This relation is important to evaluate if the
found covariance structure is reasonable and will be considered more
elaborately in Section 6.

##  5 Explained variance

Transforming the random vector using sparse loadings yields variables that are
potentially correlated. This is a concern when calculating the explained
variance for each loading or for each block of variables. Solutions to that
issue will be addressed in this section.

Since $\mathbb{E}[(\bm{X}\bm{v}_{i})^{\top}\bm{X}\bm{v}_{j}]=0$ and
$\mathbb{E}[(\bm{X}\bm{v}_{i})^{\top}\bm{X}\bm{v}_{i}]=\lambda_{j}$ for
$i,j\in\\{1,\ldots,M\\}$ and $i\neq j$, it is well known that the eigenvalues
contain the magnitude of dispersion of $\bm{X}$ along the eigenvectors. Hence,
the overall distortion of the covariance matrix is given by
$\sum_{i}\lambda_{i}$. In the sparse case, however, the projected variables
$\\{\bm{X}\bm{u}_{i}\\}\equiv\\{\bm{X}\bm{u}_{i}\\}_{i\in\\{1,\ldots,M\\}}$
might be correlated such that
$\mathbb{E}[(\bm{X}\bm{u}_{i})^{\top}\bm{X}\bm{u}_{j}]\neq 0$ for
$i,j\in\\{1,\ldots,M\\}$. Therefore, quasi-eigenvalues
$\mathbb{E}[(\bm{X}\bm{u}_{i})^{\top}\bm{X}\bm{u}_{i}]=\bm{u}_{i}^{\top}\tilde{\bm{\Sigma}}\bm{u}_{i}$
overvalue the explained total variance because the $\\{\bm{X}\bm{u}_{i}\\}$
contain similar contributions.

To find a different measure for the explained variance, we recap that the
partial covariance matrix $\tilde{\bm{\Sigma}}_{2\cdot 1}$ contains the
information left in $\bm{X}_{\mathcal{D}}$ after eliminating the effects of
$\bm{X}_{\mathcal{K}}$ by regression projection.

######  Remark 1.

Let $\bm{X}\equiv(\bm{X}_{\mathcal{K}},\bm{X}_{\mathcal{D}})$ with
$\mathbb{E}[\bm{X}^{\top}\bm{X}]=\tilde{\bm{\Sigma}}$. It holds that

$\mathbb{E}[(\bm{X}_{\mathcal{D}}-\bm{X}_{\mathcal{K}}\bm{\beta})^{\top}(\bm{X}_{\mathcal{D}}-\bm{X}_{\mathcal{K}}\bm{\beta})]=\bm{\Sigma}_{2}-\bm{\mathrm{E}}_{12}^{\top}\bm{\Sigma}_{1}^{-1}\bm{\mathrm{E}}_{12}\equiv\tilde{\bm{\Sigma}}_{2\cdot
1}\;,$ (4)

where $\bm{\beta}$ contains the respective regression coefficients.

Hence, an initial proposal is to use a measure for the explained variance that
is based on the partial covariance matrix $\tilde{\bm{\Sigma}}_{2\cdot 1}$ to
overcome the concern of correlation among the $\\{\bm{X}\bm{u}_{i}\\}$. McCabe
(1984) discussed several measures and we will use one of them for our examples
in Section 8. However, in order to find the variables we want to discard as
well as the number of variables we want to discard, each partial covariance
matrix for any possible combination of variables $\sum_{m=1}^{M-1}{M\choose
m}=2^{M}-2$ has to be calculated which is not feasible in practice. Further,
calculating the partial covariance matrix is computationally demanding.
Therefore, we suggest to firstly use a different evaluation of the explained
variance to find possible blocks of variables we consider to discard.

Zou et al. (2006) proposed to correct the explained variance using regression
projection as well. For calculating the explained variance of
$\bm{X}\bm{u}_{i}$, they suggest to correct upwardly for
$\\{\bm{X}\bm{u}_{j}\\}_{j\in\\{1,\ldots,i-1\\}}$ and therefore to use

$\mathbb{E}[(\bm{X}\bm{u}_{i}-\bm{X}\bm{U}_{i-1}\bm{\beta})^{\top}(\bm{X}\bm{u}_{i}-\bm{X}\bm{U}_{i-1}\bm{\beta})]\equiv
r_{i,i}^{2}\;,$ (5)

where $\bm{U}_{i-1}\equiv(\bm{u}_{1},\ldots,\bm{u}_{i-1})$, and $\bm{\beta}$
contains the respective regression coefficients. The advantage is the
following computationally more efficient measure for application in the sample
case using a QR decomposition.

######  Remark 2.

Let $\bm{x}\hat{\bm{U}}=\hat{\bm{Q}}\hat{\bm{R}}$ with
$\hat{\bm{R}}=(\hat{r}_{i,j})$ be the QR decomposition of the sample
$\bm{x}\hat{\bm{U}}$. It holds that

$\|\bm{x}\hat{\bm{u}}_{i}-\bm{x}\hat{\bm{U}}_{i-1}\hat{\bm{\beta}}\|_{2}^{2}=\hat{r}_{i,i}^{2}\;,$
(6)

where $\hat{\bm{U}}_{i-1}\equiv(\hat{\bm{u}}_{1},\ldots,\hat{\bm{u}}_{i-1})$,
and $\hat{\bm{\beta}}$ contains the respective regression coefficients.
Further, $(N-1)^{-1}\hat{r}_{i,i}^{2}\overset{p}{\longrightarrow}r_{i,i}^{2}$
is the sample counterpart of $r_{i,i}^{2}$ for $i\in\\{1,\ldots,M\\}$.

Hence, in the sample case the explained variance of $\bm{X}\bm{u}_{i}$
corrected by $\\{\bm{X}\bm{u}_{j}\\}_{j\in\\{1,\ldots,i-1\\}}$ can be
estimated by $(N-1)^{-1}\hat{r}_{i,i}^{2}$, and the explained total variance
can be estimated by $(N-1)^{-1}\sum_{i}\hat{r}_{i,i}^{2}$ respectively.
However, while (6) gives a computationally efficient measure, there are also
two downsides: Firstly, the evaluation of blocks of variables is not taken
into account since the measure considers each sparse loading as an $1\times 1$
block. Secondly, the measure is biased because we do not correct for all
sparse loadings. Since we correct for more loadings when $i$ increases, it
holds that this bias decreases the larger $i$. However, (6) still provides a
good first intuition regarding the explained variance which benefits from the
efficient calculation. Therefore, we propose that SPLA is based on (6) for a
first analysis to reduce the number of block-combinations that have to be
evaluated as well as on the partial covariance matrix (4) for the final
evaluation.

##  6 Identifying the Underlying Block-Structure

When calculating sparse loadings, we might create sparse loadings even if the
underlying eigenvector structure is not sparse-ish, or we might create
loadings that are sparser than the underlying eigenvector structure. In
Section 6.1, we therefore discuss how to evaluate if the obtained sparse
loadings are reasonable or not. Further, we recap methods to calculate sparse
loadings. Afterwards, we have all pieces together to provide the algorithm for
SPLA in Section 6.2. Additionally, we show in Section 6.3 that the explained
total variance of the sparse model is controlled by the introduced evaluation
in Section 6.1.

###  6.1 Sparseness Evaluation Criterion

The sparser the loadings, the easier it is to interpret the underlying
structure of the sample. Of course, it is only reasonable to enforce sparse
loadings in application if the underlying population eigenvectors are sparse-
ish. However, since we have only access to the sample eigenvectors in
practice, we need a criterion to decide if the found sparseness is reasonable
or not.

Let therefore $\bm{u}_{\delta^{\ast}}$ be the sparse loading associated with
the block $\\{X_{d}\\}$ with the smallest index. As a sparseness criterion, we
propose to evaluate the ratio between the corrected explained variance
$r_{{\delta^{\ast}},{\delta^{\ast}}}$ of $\bm{u}_{\delta^{\ast}}$ and its
uncorrected variance
$\mathbb{E}[(\bm{X}\bm{u}_{\delta^{\ast}})^{2}]=\bm{u}_{{\delta^{\ast}}}^{\top}\tilde{\bm{\Sigma}}\bm{u}_{{\delta^{\ast}}}$.
We use $\bm{u}_{\delta^{\ast}}$ because if we used another $\bm{u}_{d}$ with
$d<{\delta^{\ast}}$, then this sparse loading would be corrected by
$\bm{u}_{\delta^{\ast}}$ i.e. it would be corrected by variables included in
the same block. It holds the following result for the corrected explained
variance given by the largest sparse loading $\bm{u}_{\delta^{\ast}}$:

######  Theorem 1.

Let $1\notin\mathit{\Delta}$ and let
$\\{\bm{u}_{\delta}\\}_{\delta\in\mathit{\Delta}}$ be the sparse loadings
associated with the block $\\{X_{d}\\}$. It holds for
$r_{{\delta^{\ast}},{\delta^{\ast}}}^{2}\equiv\mathbb{E}[(\bm{X}\bm{u}_{\delta^{\ast}}-\bm{X}\bm{U}_{{\delta^{\ast}}-1}\bm{\beta})^{\top}(\bm{X}\bm{u}_{\delta^{\ast}}-\bm{X}\bm{U}_{{\delta^{\ast}}-1}\bm{\beta})]$
that

$r_{{\delta^{\ast}},{\delta^{\ast}}}^{2}=\bm{u}_{{\delta^{\ast}}}^{\top}\tilde{\bm{\Sigma}}\bm{u}_{{\delta^{\ast}}}-\mathcal{O}_{p}(\|\bm{\mathrm{E}}\bm{u}_{{\delta^{\ast}}}\|_{2}^{2}\|\bm{\Sigma}_{1}^{-1}\|_{2})\;,$

where ${\delta^{\ast}}\equiv\min\mathit{\Delta}$ is the smallest index of the
sparse loadings associated with the $\\{X_{d}\\}$,
$\bm{U}_{{\delta^{\ast}}-1}\equiv(\bm{u}_{1},\ldots,\bm{u}_{{\delta^{\ast}}-1})$,
and $\bm{\beta}$ contains the respective regression coefficients.

The bound depends on the covariances between the $\\{X_{d}\\}$ and the
$\\{X_{k}\\}$, and on the precision matrix of the $\\{X_{k}\\}$. The latter
one follows the indication by Bauer and Drabant (2021) that blocks become more
distinguishable the stronger the relation is within the blocks. In case the
underlying covariance matrix follows a perfect block-diagonal structure, i.e.
if $\bm{\mathrm{E}}=\bm{0}$, we can conclude that
$r_{{\delta^{\ast}},{\delta^{\ast}}}^{2}=\bm{u}_{{\delta^{\ast}}}^{\top}\tilde{\bm{\Sigma}}\bm{u}_{{\delta^{\ast}}}$
reduces to the largest sparse quasi-eigenvalue. We note, however, that the
bound depends not on the plain covariances between the blocks but rather on a
weighted covariance $\bm{\mathrm{E}}\bm{u}_{{\delta^{\ast}}}$ with weights
$\bm{u}_{{\delta^{\ast}}}$. This is problematic because different loadings of
same shape yield different evaluation outcomes. Additionally, variables with
zero-components in the corresponding loading are not considered in the
criterion. However, we can fix this issue by giving equal weights to all
variables. Therefore, the closer the ratio of the following criterion is to
one, the more likely we detected the correct underlying sparse structure.

######  Remark 3.

Let $1\notin\mathit{\Delta}$ and let
${\delta^{\ast}}\equiv\min\mathit{\Delta}$ be the smallest index of the sparse
loadings associated with the $\\{X_{d}\\}$. Let further
$\bm{w}_{{\delta^{\ast}}}\equiv(\bm{0}_{K},\bm{1}_{D})^{\top}$. We replace the
${\delta^{\ast}}$th loading by
$\bm{u}_{{\delta^{\ast}}}^{\top}\equiv\bm{w}_{{\delta^{\ast}}}/\|\bm{w}_{{\delta^{\ast}}}\|_{2}=D^{-1/2}\bm{w}_{{\delta^{\ast}}}$
and we replace $\\{\bm{u}_{i}\\}_{i\neq{\delta^{\ast}}}$ such that Assumption
3 is satisfied. The block evaluation criterion (EC) for $\\{X_{d}\\}$ is given
by

${\rm
EC}\equiv\frac{r_{{\delta^{\ast}},{\delta^{\ast}}}^{2}}{\bm{u}_{{\delta^{\ast}}}^{\top}\tilde{\bm{\Sigma}}\bm{u}_{{\delta^{\ast}}}}\in(0,1]$
(7)

and it holds that

${\rm
EC}=\frac{\bm{u}_{{\delta^{\ast}}}^{\top}\tilde{\bm{\Sigma}}\bm{u}_{{\delta^{\ast}}}-\mathcal{O}_{p}(\|D^{-1/2}\bm{\mathrm{E}}_{12}\bm{1}_{D}\|_{2}^{2}\|\bm{\Sigma}_{1}^{-1}\|_{2})}{\bm{u}_{{\delta^{\ast}}}^{\top}\tilde{\bm{\Sigma}}\bm{u}_{{\delta^{\ast}}}}$

due to Theorem 1.

Clearly, (7) equals one for the extreme case when $\bm{\mathrm{E}}=\bm{0}$ and
therefore when the two random vectors $\bm{X}_{\mathcal{K}}$ and
$\bm{X}_{\mathcal{D}}$ are perfectly uncorrelated. Further, the criterion is
positive by construction because it is a scaled variance term with
$\tilde{\bm{\Sigma}}$ being positive definite (Assumption 2).

The vector $\bm{w}_{{\delta^{\ast}}}$ weighs all covariances between
$\\{X_{d}\\}$ and $\\{X_{k}\\}$ equally making the EC a corrected criterion.
When replacing a loading $\bm{u}_{{\delta^{\ast}}}$, we have to change all
remaining loadings $\\{\bm{u}_{i}\\}_{i\neq{\delta^{\ast}}}$ to maintain
orthogonality of $\bm{U}$. However, such a change is quickly implemented and a
respective procedure is given in Supplementary Material 1.

For completion, we shall mention that we obtain similar bounds to Theorem 1
also for the case when the block-diagonal shape of the loadings given in (3)
is not satisfied. However, in this case it is not feasible to correct for all
$K$ loadings associated with the block $\\{X_{k}\\}$ since we would correct
for variables contained in the same block otherwise.

######  Lemma 3.

Let $1\notin\mathit{\Delta}$ and let
$\\{\bm{u}_{\delta}\\}_{\delta\in\mathit{\Delta}}$ be the sparse loadings
associated with the block $\\{X_{d}\\}$. If the block-diagonal condition of
$\bm{U}$ in Assumption 3 is violated, it holds for
$r_{{\delta^{\ast}},{\delta^{\ast}}}^{2}\equiv\mathbb{E}[(\bm{X}\bm{u}_{\delta^{\ast}}-\bm{X}\bm{U}_{{\delta^{\ast}}-1}\bm{\beta})^{\top}(\bm{X}\bm{u}_{\delta^{\ast}}-\bm{X}\bm{U}_{{\delta^{\ast}}-1}\bm{\beta})]$
that

$r_{{\delta^{\ast}},{\delta^{\ast}}}^{2}=\bm{u}_{\delta^{\ast}}^{\top}\tilde{\bm{\Sigma}}\bm{u}_{\delta^{\ast}}-\mathcal{O}_{p}(\|\bm{\mathrm{E}}\bm{u}_{\delta^{\ast}}\|_{2}^{2}\|\bm{\Sigma}_{1}^{-1}\|_{2})\;,$

where ${\delta^{\ast}}\equiv\min\mathit{\Delta}$ is the smallest index of the
sparse loadings associated with the $\\{X_{d}\\}$,
$\bm{U}_{{\delta^{\ast}}-1}\equiv(\bm{u}_{\kappa_{1}},\bm{u}_{\kappa_{2}},\ldots)$
with $\kappa_{1},\kappa_{2},\ldots<{\delta^{\ast}}$ and ${\rm
card}(\\{\kappa_{1},\kappa_{2},\ldots\\})<K$, and $\bm{\beta}$ contains the
respective regression coefficients.

We evaluate the sparseness only in an upward manner in the EC due to the
nature of the measure $r_{i,i}^{2}$ from $(\ref{eq:expvarzouPopulation})$ as
discussed in Section 5 in relation to Remark 2. Further, following Theorem 1
we evaluate a block with respect to its largest associated sparse loading. If
we used a different loading $\bm{u}_{\delta}$ with
$\delta\neq{\delta^{\ast}}\equiv\min\mathit{\Delta}$, the bound on the quasi-
eigenvalue would not only depend on $\bm{\mathrm{E}}_{12}$ and
$\bm{\Sigma}_{1}^{-1}$, but also on $\bm{\Sigma}_{2}$.

Note that ${\delta^{\ast}}=1$, i.e. using the first loading $\bm{u}_{1}$, is
not feasible for the evaluation criterion EC. We have to keep this in mind for
the general case with more than two blocks. Despite ${\delta^{\ast}}=1$ is
mathematically possible in (6), no additional information is provided since
loadings are corrected only in an upward direction. Therefore, the first
loading is not corrected for the effect of other loadings and therefore the EC
always equals one in this case.

######  Corollary 3.

Let $\\{\bm{u}_{\delta}\\}_{\delta\in\mathit{\Delta}}$ be the sparse loadings
associated with the block $\\{X_{d}\\}$. If
${\delta^{\ast}}\equiv\min\mathit{\Delta}=1$, then the block-structure
evaluation criterion EC from (7) reduces to one.

If this situation arises, however, one evaluates the block $\\{X_{k}\\}$
instead of $\\{X_{d}\\}$ by using its largest loading $\bm{u}_{\kappa^{\ast}}$
with $\kappa^{\ast}=\min\\{1,\ldots,M\\}\backslash\mathit{\Delta}$. Clearly,
verifying the block $\\{X_{k}\\}$ implies verification of $\\{X_{d}\\}$ as a
block. The condition that $1\notin\mathit{\Delta}$ in Theorem 1 is fulfilled
by construction in this work since we assumed for convenience that
$\mathcal{D}=\mathit{\Delta}=\\{K+1,\ldots,M\\}$ and therefore that
${\delta^{\ast}}\equiv\min\mathit{\Delta}=K+1$. However, in a more general
setting $1\in\mathit{\Delta}$ is possible as can be seen in the examples in
Section 8.

We note that other evaluation criteria based on the partial covariance matrix
can also be used.

######  Lemma 4.

Let $\tilde{\bm{\Sigma}}_{2\cdot
1}\equiv\bm{\Sigma}_{2}-\bm{\mathrm{E}}_{12}^{\top}\bm{\Sigma}_{1}^{-1}\bm{\mathrm{E}}_{12}$
be the partial covariance matrix. It holds that

$\frac{\|\tilde{\bm{\Sigma}}_{2\cdot
1}\|}{\|\bm{\Sigma}_{2}\|}=\frac{\|\bm{\Sigma}_{2}\|-\mathcal{O}_{p}(\|\bm{\mathrm{E}}_{12}\|^{2}\|\bm{\Sigma}_{1}^{-1}\|)}{\|\bm{\Sigma}_{2}\|}\in(0,1]\;.$

As mentioned in Section 5, however, calculation of the partial covariance
matrix is computationally costly.

###  6.2 Sparse Principal Loading Analysis

In Section 5, we proposed that a measure for the explained variance in SPLA
should be based on the QR decomposition (6) for a first analysis as well as on
the partial covariance matrix (4) for the final evaluation. To decide if the
found block-structure is reasonable, we additionally check the evaluation
criterion.

The missing piece is the calculation of the sparse loadings. Guerra-Urzola et
al. (2021) compared the performance of several methods and assessed that the
one by Shen and Huang (2008) performs best. Since it is a special case of the
approach by Witten et al. (2009), we use their method for the majority of this
work. For some of our examples, however, we applied the procedure of Zou et
al. (2006). The latter one allows sparsening of each loading individually by
penalization which allows to create any composition of sparse loadings and has
been used for the examples in Section 8.1 and Section 8.2 respectively. Bot
methods are recapped in Supplementary Material 1, however, we refer to the
original works for elaborate explanations.

We can now contribute an algorithm for SPLA. Discard $\\{X_{d}\\}$ according
to SPLA proceeds as follows:

Algorithm 1 Sparse principal loading analysis

1:Calculate the sparse loadings $\\{\bm{u}_{i}\\}$.

2:Assess if the detected sparseness is meaningful by checking the EC from (7).
Repeat step 1-2 as long as the sparseness is reasonable.

3:Evaluate the detected blocks of variables using their explained variances
from (6) and select potential blocks for discarding accordingly.

4:Check the choice in step 3 using the partial covariance matrices from (4).
Discard the blocks if selection is verified.

In case we are only interested in finding the underlying block covariance
structure, the first three steps 1 to 3 of Algorithm 1 are sufficient.
Further, in step 2 we decide that the detected sparseness is not reasonable
when the evaluation criterion falls below a fixed parameter $c_{\rm EC}$.

###  6.3 Variance-Interpretability-Trade-Off

In this section, we discuss the trade-off between the explained total variance
and interpretability, and elaborate the importance of a large explained total
variance. We will show that the variance-interpretability-trade-off is
addressed appropriately when the evaluation of the underlying sparseness
discussed in Section 6.1 is considered.

For applications, we have to be aware of the trade-off between the explained
total variance and interpretability of SPLA. We discussed in Section 6.1 that
the more sparse the loadings, the easier it is to interpret the underlying
sample. However, this comes at the price of a small explained total variance
which, on the other hand, might make interpretation not reasonable anymore
because the sparse loadings simply provide too little information.

In application, the explained total variance of the sample sparse loadings
according to (6) relative to the total variance of $\bm{x}$ is given by
$(N-1)^{-1}\sum_{i}\hat{r}_{i,i}^{2}/{\rm tr}(\hat{\tilde{\bm{\Sigma}}})$. Due
to the perturbation
$\hat{\bm{\mathrm{E}}}\equiv\hat{\tilde{\bm{\Sigma}}}-\hat{\bm{\Sigma}}$, we
obtain the following result:

######  Theorem 2.

Let $\varepsilon(\bm{x})$ be a perturbation of the sample such that
$\check{\bm{x}}\equiv\bm{x}-\varepsilon(\bm{x})$ has covariance matrix
$\hat{\bm{\Sigma}}$. It holds that

$\|\varepsilon(\bm{x})\|=\mathcal{O}_{p}(\|\hat{\bm{\mathrm{E}}}\|)\;.$

If further
$(N-1)^{-1}\|\hat{\bm{\Sigma}}^{-1}\|_{2}^{1/2}\|\varepsilon(\bm{x})\|_{2}<1$,
it holds that

$(N-1)^{-1}\sum\limits_{i}\hat{r}_{i,i}^{2}={\rm
tr}(\hat{\tilde{\bm{\Sigma}}})-(N-1)^{-1}\mathcal{O}_{p}(\|\hat{\bm{\mathrm{E}}}\|_{F}^{2})\;.$

Theorem 2 contains two implications. Firstly, since
$\|\varepsilon(\bm{x})\|=\mathcal{O}_{p}(\|\hat{\bm{\mathrm{E}}}\|)$, the
condition that
$(N-1)^{-1}\|\hat{\bm{\Sigma}}^{-1}\|_{2}^{1/2}\|\varepsilon(\bm{x})\|_{2}<1$
evaluates the relation within blocks in comparison to the relation between
blocks as discussed in Section 6.1. Secondly, the explained total variance
drops in case we erroneously assign a block. In turn, the explained total
variance is large in case we correctly detect a block. However, recap that we
already control for correct specification of blocks in SPLA due to the EC from
(7). We can conclude that the variance-interpretability-trade-off is of no
concern when we follow Algorithm 1 and therefore evaluate the sparseness
according to the evaluation criterion.

##  7 Simulation Study

In this section, we deepen the understanding of SPLA based on simulations. In
Section 7.1, we analyse the behaviour of the EC from Section 6 for block-
diagonal covariance detection when the correlation among blocks varies.
Afterwards, we evaluate the percentage of identified block-diagonal covariance
matrices (identification rate) of SPLA in Section 7.2. We further discuss
considerations regarding our simulation design in Section 7.3. Code to
replicate all simulation results is provided in Supplementary Material 2.

###  7.1 Corrected Evaluation Criterion Analysis

In this section, we simulate a sample following a block-diagonal correlation
matrix of six blocks to demonstrate the evaluation using the EC.

Figure 1: Sample correlation matrices of $\bm{x}^{[0.1]}$, $\bm{x}^{[0.2]}$,
$\bm{x}^{[0.4]}$, and $\bm{x}^{[0.5]}$ respectively.

Let therefore $Y$ and $Z_{i}$ for $i\in\\{1,\ldots,5\\}$ be IID $N(0,10)$
distributed, and let $W$ be $N(0,1)$ distributed. We consider a random vector
$\bm{X}^{[\rho]}=(X_{1},\ldots,X_{14})$ with
$X_{j}=\sqrt{1-\rho}Z_{i}+\sqrt{\rho}Y+W$ for
$(j,i)\in\\{(1,1),(2,1),(3,2),(4,2),\ldots,(13,6),(14,6)\\}$ such that
$\bm{X}^{[\rho]}$ has a covariance matrix with six $2\times 2$ blocks on the
diagonal. Hereby, $\rho$ reflects the approximate correlation between the
blocks. We simulate a sample $\bm{x}^{[\rho]}$ by firstly drawing $N=100$
observations from $\bm{X}^{[\rho]}$ and secondly standardizing the sample.
Therefore, the sample covariance matrix equals the sample correlation matrix
which is illustrated in Figure 1 for $\rho\in\\{0.1,0.2,0.4,0.5\\}$.
Afterwards, we redraw $s\in\\{1,\ldots,100\\}$ samples $\bm{x}^{[\rho|s]}$
with $N=100$ observations and calculate the respective EC according to
$(\ref{eq:CEC})$ for block two, four, and six. The results are illustrated as
boxplots in Figure 2.

Figure 2: EC for block two, four, and six for each of the
$s\in\\{1,\ldots,100\\}$ samples $\bm{x}^{[\rho|s]}$ drawn from
$\bm{X}^{[\rho]}$ for $\rho\in\\{0.1,0.2,0.3,0.4,0.5\\}$ respectively.

The EC decreases with an increase of the correlation among the blocks and,
therefore, with an increase of $\rho$. Further, the EC decreases when
evaluating more blocks. This is due to the construction of the corrected
measure for the explained variance in $(\ref{eq:expvarzou})$ because we
control in an upward direction.

According to the EC for the sixth block in Figure 2, $c_{\rm EC}=0.6$ seems
appropriate because it confirms the block for $\rho<0.3$ and rejects the block
for $\rho>0.5$. This is meaningful, since strong correlation among variables
tells us that they are not separated into different blocks. We therefore use
$c_{\rm EC}=0.6$ for the remainder of this work. For completion, we shall
emphasize that choosing the appropriate sparseness depends on application.
Choosing a too small $c_{\rm EC}$, however, is no concern because erroneously
detected blocks are corrected in step 4 of Algorithm 1 using the partial
covariance matrix.

###  7.2 Performance

In this section, we compare the performance of block-diagonal covariance
structure detection. Recently, Devijver and Gallopin (2018) derived
nonasymptotic approaches that select the number of blocks either using a slope
heuristic robust regression (SHRR) or using a slope heuristic dimension jump
(SHDJ). We refer to their work for an elaborate explanation. These approaches
outperform existing methods when the number of blocks is unknown a priori and
we therefore compare SPLA to them.

Similar to Section 7.1, we simulate $s\in\\{1,\ldots,100\\}$ samples
$\bm{x}^{[\rho|s]}$ for $\rho\in\\{0,0.1,\ldots,0.9\\}$ with number of
observations $N\in\\{1000,500,100,50,20\\}$. We compare the performance of the
methods according to the identification rate and the results are illustrated
in Figure 3.

Figure 3: Performance of the block-diagonal covariance structure detection by
SHDJ, SHRR, and SPLA with $c_{\rm EC}=0.6$ measured by the identification rate
(the percentage of samples where all blocks where identified) on
$s\in\\{1,\ldots,100\\}$ samples $\bm{x}^{[\rho|s]}$ for
$\rho\in\\{0,0.1,\ldots,0.7\\}$ with number of observations
$N\in\\{1000,500,100,50,20\\}$. Further, the EC for each block is illustrated
for the case when $N=100$.

SPLA performs best for large sample sizes. The identification rate decreases
for larger $\rho$ since the covariance matrix can no longer be considered to
be a block-diagonal matrix due to the increasing correlation among variables.
This is illustrated by the decreasing EC. Still, we could increase
identification rate by assuming a small EC to be reasonable as well. This is
illustrated in Figure 4, where the identification rate continuous to be high
for an increasing $\rho$. However, this demonstration of adapting $c_{\rm EC}$
is only given for completion since it might be useful for a different context.
One should be aware that a small EC might not reflect a reasonable block-
structure as discussed in Section 6.

Figure 4: Performance of the block-diagonal covariance structure detection by
SPLA for $c_{\rm EC}\in\\{0.1,0.2,0.3,0.4,0.5,0.6\\}$ measured by the
identification rate (the percentage of samples where all blocks where
identified) on $s\in\\{1,\ldots,100\\}$ samples $\bm{x}^{[\rho|s]}$ for
$\rho\in\\{0,0.1,\ldots,0.9\\}$ with number of observations
$N\in\\{500,100,20\\}$. Identification rates for SHDJ are drawn to make the
graphics comparable to Figure 3.

SHDJ and SHRR increase performance the smaller the sample size because they
are nonasymptotic approaches. Further, they can also be applied for the case
when $M>N$ which is not discussed for SPLA in this work.

###  7.3 Considerations about the Simulation Design

Simulations always face limitations due to their artificialness. However,
there are further considerations due to the nature of SPLA. In the previous
sections, we constructed a spiked covariance model (Johnstone, 2001) similar
to Zou et al. (2006) and Shen and Huang (2008) respectively to cause strong
correlations within each block. The reason is that SPLA might rightfully split
an intended block into several blocks if the within-block-correlation is
small, making it difficult to analyse the obtained results. As an
illustration, we simulated $1000$ times the matrix
$\bm{A}\in\mathbb{R}^{3\times 3}$ by drawing each element $a_{i,j}$ from a
uniform $U(0,1)$ distribution. The covariance matrix
$\hat{\tilde{\bm{\Sigma}}}=\bm{A}^{\top}\bm{A}$ was constructed accordingly
and we consider the respective correlation matrix. By construction, the
correlation matrix reflects three correlated variables. However, in Figure 5
we see that SPLA detected one, two, or three blocks respectively for different
realisations. The EC is small when SPLA detected one block (i.e. no split of
the variables) and large for two or three blocks respectively which indicates
that the found block-structure was reasonable.

Figure 5: Number of detected blocks by SPLA in relation to the EC for
simulated samples with correlation matrices obtained from the simulated
covariance matrix $\hat{\tilde{\bm{\Sigma}}}=\bm{A}^{\top}\bm{A}$.

Therefore, our simulations follow a spiked covariance model having a large
within-block-correlation to ensure that an intended block cannot be
interpreted as several blocks from a SPLA perspective.

##  8 Examples

We provide two examples for SPLA on real data and one example on synthetic
data in this section. The first example in Section 8.1 mainly serves to
illustrate detection of the underlying block covariance structure and results
are compared to the ones obtained by PLA. In the second example in Section
8.2, the principal focus lies on variable selection by SPLA. Afterwards, we
discuss a synthetic example from the sparse principal component analysis
context in Section 8.3. Code to replicate all examples is provided in
Supplementary Material 3.

###  8.1 OECD Sample

For the first example, we use data from Mankiw et al. (1992) that consist of
$N=22$ OECD countries with population larger than one million. The sample
contains the real GDP in 1960 (Y60) and in 1985 (Y85) per person of working
age respectively. Further, it consists of the average annual ratios of real
investment to GDP (I/Y), the percentage of working-age population that is in
secondary school (SCH), the annual population growth from 1960 to 1985 (POP),
and the average of annual ratios of gross domestic expenditure on research and
development to nominal GDP (RD). Data are available in the CRAN contributed
package AER by Kleiber and Zeileis (2008). We standardized the sample i.e. we
consider the correlation matrix rather than the covariance matrix because the
variables have different scales.

The detected blocks by both PLA and SPLA are summarized in Table 1. A large
number of blocks $K$ eases interpretation while a large EC indicates that the
block-structure represents the underlying sample. Therefore, splitting the
OECD sample into the four blocks {I/Y}, {SCH}, {POP}, and {RD, Y85, Y60} with
EC = 0.84 is an appropriate choice since this provides a reasonable as well as
more detailed representation of the underlying sample. The corresponding
sparse loadings which we obtain when performing SPLA according to Algorithm 1
are given in Table 2. In general, we continuously increase the sparseness of
the loadings to find the underlying structure until the EC indicates that
sparseness is too large. This procedure is indicated in Table 1, where we stop
with four blocks since five blocks result in a small EC.

Table 1: Number of blocks $K$ of the OECD sample detected by PLA and SPLA. For PLA, we either provide the threshold $\tau$ or we denote $--$ if the structure cannot be identified by any threshold. For SPLA, we give the smallest EC of all detected blocks to evaluate if the respective structure is reasonable. $K$ | Blocks | PLA | SPLA
---|---|---|---
2 | {I/Y, POP}, {SCH, RD, Y85, Y60} | $--$ | $0.96$
2 | {I/Y, POP, SCH}, {RD, Y85, Y60} | $0.40$ | $0.84$
3 | {I/Y}, {POP}, {SCH, RD, Y85, Y60} | $--$ | $0.96$
3 | {RD}, {Y85, Y60}, {I/Y, SCH, POP} | $0.50$ | $0.53$
4 | {I/Y}, {POP}, {SCH}, {RD, Y85, Y60} | $--$ | $0.84$
5 | {I/Y}, {POP}, {SCH}, {RD}, {Y85, Y60} | $--$ | $0.45$

Note that PLA does not detect as many combinations as SPLA due to its nature
of hard-thresholding. Further, PLA does not provide information to decide
which number of blocks to choose. This is a concern because the detected
blocks with $\tau=0.5$ could be chosen which are a poor representation of the
underlying sample according to the EC.

Table 2: Sample sparse loadings $\hat{\bm{u}}_{i}$ with $i\in\\{1,\ldots,6\\}$ of the covariance matrix of the standardized OECD sample. All loading components equal to zero are left blank. ${\rm SV}\equiv\hat{r}_{i,i}^{2}/{\rm tr}(\hat{\tilde{\bm{\Sigma}}})$ denotes the share of explained total variance (in percent), CV denotes the cumulative share of explained total variance for each block (in percent). $\hat{\bm{u}}_{i}$ | $\hat{\bm{u}}_{3}$ | $\hat{\bm{u}}_{4}$ | $\hat{\bm{u}}_{2}$ | $\hat{\bm{u}}_{1}$ | $\hat{\bm{u}}_{5}$ | $\hat{\bm{u}}_{6}$
---|---|---|---|---|---|---
I/Y | 1.00 |  |  |  |  |
SCH |  | 1.00 |  |  |  |
RD |  |  | 1.00 |  |  |
POP |  |  |  | 0.50 | -0.87 | 0.00
Y85 |  |  |  | 0.64 | 0.36 | 0.68
Y60 |  |  |  | 0.59 | 0.34 | -0.73
EC |  | 0.96 | 0.93 | 0.84
SV | 16.67 | 16.04 | 15.57 | 40.26
CV | 16.67 | 32.71 | 48.28 | 88.54

Variable selection appears not to be feasible in this case, because no block
has a small cumulative share of explained total variance and, therefore, there
is no block that distorts the covariance matrix only a little. As discussed in
Section 5 however, we have to be aware that the explained variance given by
$\bm{R}$ is biased. Following Algorithm 1, we therefore check the results
using the partial covariance matrices. Instead of calculating the partial
covariance matrix for all possible combinations $\sum_{m=1}^{6-1}{6\choose
m}=62$, we are now left with calculations for the three respective blocks. The
partial covariance does not indicate that any block should be discarded either
and the respective results can be found in Table 5 for completion.

###  8.2 EXAM Sample

The second example deals with the EXAM sample from Mardia et al. (1979), where
we have access to grades of $N=88$ students in a five part exam containing the
subjects Mechanics (mec), Vectors (vec), Algebra (alg), Analysis (ana), and
Statistics (stat). Each subject has been scored separately. The sample was
taken from the CRAN contributed package bootstrap by Tibshirani (2019).
Following the procedure in Algorithm 1, we obtain the results given in Table
3.

Table 3: Sparse loadings $\bm{u}_{i}$ with $i\in\\{1,\ldots,5\\}$ of the covariance matrix of the EXAM sample. All loading components equal to zero are left blank. ${\rm SV}\equiv r_{i,i}^{2}/{\rm tr}(\tilde{\bm{\Sigma}})$ denotes the share of explained total variance (in percent), CV denotes the cumulative share of explained total variance for each block (in percent). $\bm{u}_{i}$ | $\bm{u}_{1}$ | $\bm{u}_{2}$ | $\bm{u}_{3}$ | $\bm{u}_{4}$ | $\bm{u}_{5}$
---|---|---|---|---|---
vec | 1.00 |  |  |  |
mec |  | 1.00 |  |  |
alg |  |  | -0.09 | -0.99 | 0.00
ana |  |  | -0.65 | 0.06 | -0.76
sta |  |  | -0.75 | 0.06 | 0.66
EC |  | 0.74 | 0.72
SV | 13.21 | 19.28 | 38.98
CV | 13.21 | 32.49 | 71.47

Due to the evaluation criteria of $0.74$ and $0.72$, we can conclude that the
block-structure consisting of the three blocks {vec}, {mec}, and {alg, ana,
sta} reflects the EXAM sample well. Further, the shared explained total
variance of {vec} equals 13.21%. Since a single variable out of five should on
average explain 20% and since the 13.27% are biased upwardly as discussed in
Section 5, we consider to discard this variable. We verify or falsify this
selection by checking the partial covariance matrix. Accordingly, the block
{vec} explains only 7.45% (Table 5) and discarding is verified.

###  8.3 Synthetic Example

We use an example based on Zou et al. (2006) and Shen and Huang (2008) to
illustrate the block covariance structure detection of a synthetic sample.
First, we consider a random vector $\bm{X}=(X_{1},\ldots,X_{8})$ with
$X_{j}=Z_{1}+\theta_{j}$ for $j\in\\{1,\ldots,4\\}$ and
$X_{j}=Z_{2}+\theta_{j}$ for $j\in\\{5,\ldots,8\\}$. $Z_{1}$, $Z_{2}$, and the
$\theta_{j}$ are independent and we let $Z_{1}\sim N(0,290)$, $Z_{2}\sim
N(0,300)$ and $\theta_{j}\sim N(0,1)$ for $j\in\\{1,\ldots,10\\}$. A synthetic
sample $\bm{x}$ consisting of $N=5000$ observations is generated accordingly.
By construction, the sample consists of the two blocks
$\\{X_{1},\ldots,X_{4}\\}$ and $\\{X_{5},\ldots,X_{8}\\}$. SPLA detects this
underlying block-structure which is reflected by a EC of $0.9998$ (Table 4).

We extend the synthetic example and add the two random variables
$X_{j}=-0.3\cdot Z_{1}+0.925\cdot Z_{2}+\theta_{j}$ for $j\in\\{9,10\\}$ to
the random vector $\bm{X}=(X_{1},\ldots,X_{10})$. Both Zou et al. (2006) and
Shen and Huang (2008) discuss the appropriate sparseness of the loadings.
Hereby, Zou et al. (2006) use the knowledge about the synthetic sample while
the approach of Shen and Huang (2008) is based on the explained variance of
the most important loadings in a variance-sense. Expressing their results in
SPLA-sense, they discuss if the loadings are of shape

$\begin{pmatrix}\bm{\ast}_{4\times 4}&\bm{0}&\bm{0}\\\
\bm{0}&\bm{\ast}_{4\times 4}&\bm{0}\\\ \bm{0}&\bm{0}&\bm{\ast}_{2\times
2}\end{pmatrix}\;\text{ or }\;\begin{pmatrix}\bm{\ast}_{4\times 4}&\bm{0}\\\
\bm{0}&\bm{\ast}_{6\times 6}\end{pmatrix}$

respectively and therefore if the three blocks $\\{X_{1},\ldots,X_{4}\\}$,
$\\{X_{5},\ldots,X_{8}\\}$ and $\\{X_{9},X_{10}\\}$, or if the two blocks
$\\{X_{1},\ldots,X_{4}\\}$ and $\\{X_{5},\ldots,X_{10}\\}$ represent the
random vector. From a SPLA-perspective, we clearly confirm the two-block-
structure with an EC of $0.9910$ and results can be found in Table 4. The
block $\\{X_{9},X_{10}\\}$ in the three-block-scenario has an EC of $0.0009$
and we can therefore conclude that it is not adequate to split
$\\{X_{1},\ldots,X_{10}\\}$ into $\\{X_{5},\ldots,X_{8}\\}$ and
$\\{X_{9},X_{10}\\}$.

Table 4: Block indices $j$ for the $\\{X_{j}\\}$ of the synthetic examples given by $\bm{X}=(X_{1},\ldots,X_{8})$ and $\bm{X}=(X_{1},\ldots,X_{10})$ respectively. EC is the respective criterion for the blocks according to Remark 3. $\bm{X}$ | $\bm{X}=(X_{1},\ldots,X_{8})$ | $\bm{X}=(X_{1},\ldots,X_{10})$ | $\bm{X}=(X_{1},\ldots,X_{10})$
---|---|---|---
$j$ | $\\{1,\ldots,4\\}$ | $\\{5,\ldots,8\\}$ | $\\{1,\ldots,4\\}$ | $\\{5,\ldots,10\\}$ | $\\{1,\ldots,4\\}$ | $\\{5,\ldots,8\\}$ | $\\{9,10\\}$
EC |  | 0.9998 |  | 0.9910 |  | 0.9998 | 0.0009

##  9 Concluding Remarks

We introduced SPLA as a new concept for variable selection and block
covariance structure identification based on sparse loadings. We proposed how
the sparse loadings can be used to identify the underlying block-structure,
and gave a criterion to evaluate if the detected structure fits the covariance
matrix. Regarding dimensionality reduction, calculation of the explained total
variance is a concern and we provided an eligible measure to control the loss
of explanatory power. Algorithm 1 can be used to apply SPLA for variable
selection and block covariance structure identification in practice. The
algorithm is available in the CRAN contributed package prinvars by Bauer and
Holzapfel (2023).

Some improvements have to be addressed in future research. From an applied
perspective, the choice of the method to calculate the sparse loadings needs
to be discussed. This should address the computational speed by comparing
existing methods that calculate sparse loadings. However, the methods need
also be evaluated with respect to the created block-structure. From a
theoretical perspective, considering how SPLA is related to the optimal
variable selection based on the partial covariance matrix is of interest.
While SPLA yields blocks that possibly explain only little of the overall
variance, it is not clear yet how optimal these blocks are in a variable
selection-sense.

There are methods to calculate sparse loadings that are stable in the high-
dimensional setting. This motivates block covariance structure identification
in the generalized spiked covariance model (Bai and Yao, 2012) or the spiked
population model (Johnstone, 2001) as a special case where the $M\times M$
population covariance matrix is assumed to be of form

$\begin{pmatrix}{\rm Cov}(\bm{X})&\bm{0}\\\
\bm{0}&\bm{V}_{M^{\prime}}\end{pmatrix}$

where $M^{\prime}=M-P$ and $M=M(N)$ such that $M/N\to\theta>0$ ($N\to\infty$).
A natural extension of SPLA would be the identification of this block
covariance structure in this high dimensional scenario.

Further, SPLA in time series analysis and functional data analysis is of
interest. Let therefore $\\{X_{t}\\}_{t\in\mathcal{T}}$ be a stochastic
process where $\mathcal{T}$ is discrete. In time series analysis, the
autocovariance function $\gamma_{X}(s,t)\equiv{\rm Cov}(X_{s},X_{t})$ gives
the covariance with itself at pairs of time points. Further, weak stationarity
of time series requires that $\gamma_{X}(h)\equiv{\rm Cov}(X_{t},X_{t+h})$ for
$t,h\in\mathbb{Z}$ does not depend on time and therefore assumes that the
autocovariance matrix consists of one block in a SPLA-sense. Hence, if the
autocovariance matrix consists of several blocks, SPLA detects structural
breaks where the weak stationary assumption is violated. On the other hand,
SPLA can be used for functional (time) data analysis to detect time intervals
that are less important from a variance perspective. This can be done by
extending the concept of SPLA to the eigenfunctions of the covariance function
of the underlying stochastic process.

##  Appendix A Complementary Results

This section provides complementary results about the partial covariance
matrix for both examples discussed in Section 8.1 and Section 8.2
respectively.

Table 5 contains the trace of the sample partial covariance matrix for each of
the identified blocks from Section 8. Considering the trace is in line with
one of the measures given by McCabe (1984). Further, scaling the matrices by
the overall explained variance ${\rm tr}(\hat{\tilde{\bm{\Sigma}}})$ gives us
the contribution of each block corrected by all remaining blocks in percent.

Table 5: Share of the total variance explained by the sample partial covariance matrix $\hat{\tilde{\bm{\Sigma}}}_{2\cdot 1}$ of the variables $\\{X_{d}\\}$ eliminating the effects of the remaining variables $\\{X_{k}\\}$ such that $\\{X_{d}\\}\cup\\{X_{k}\\}=\\{X_{j}\\}_{j\in\\{1,\ldots,M\\}}$ both for the OECD and for the EXAM sample respectively. $\mathcal{D}$ | ${\rm tr}(\hat{\tilde{\bm{\Sigma}}}_{2\cdot 1})/{\rm tr}(\hat{\tilde{\bm{\Sigma}}})\cdot 100$ | $\mathcal{D}$ | ${\rm tr}(\hat{\tilde{\bm{\Sigma}}}_{2\cdot 1})/{\rm tr}(\hat{\tilde{\bm{\Sigma}}})\cdot 100$
---|---|---|---
{I/Y} | 10.23 | {vec} | 7.45
{SCH} | 12.41 | {mec} | 17.97
{POP} | 12.94 | {alg, ana, sta} | 46.49
{RD, Y85, Y60} | 41.73 |  |

SUPPLEMENTARY MATERIAL

  1. 1.

Methods to calculate sparse loadings and weight matrix construction for Remark
3.

  2. 2.

Simulation results

  3. 3.

Examples

  4. 4.

Proofs

## References

  * Bai and Yao (2012) Bai, Z. and J. Yao (2012).  On sample eigenvalues in a generalized spiked population model.  Journal of Multivariate Analysis 106, 167–177.
  * Bauer (2021) Bauer, J. O. (2021).  Correlation based principal loading analysis.  4th International Conference on Mathematics and Statistics, 27–34.
  * Bauer and Drabant (2021) Bauer, J. O. and B. Drabant (2021).  Principal loading analysis.  Journal of Multivariate Analysis 184.
  * Bauer and Drabant (2023) Bauer, J. O. and B. Drabant (2023).  Regression based thresholds in principal loading analysis.  Journal of Multivariate Analysis 193.
  * Bauer and Holzapfel (2023) Bauer, J. O. and R. Holzapfel (2023).  prinvars: Principal Variables.  R package version 1.0.0.
  * Bickel and Levina (2008) Bickel, P. J. and E. Levina (2008).  Covariance regularization by thresholding.  The Annals of Statistics 36(6), 2577–2604.
  * Devijver and Gallopin (2018) Devijver, E. and M. Gallopin (2018).  Block-diagonal covariance selection for high-dimensional gaussian graphical models.  Journal of the American Statistical Association 113(521), 306–314.
  * Guerra-Urzola et al. (2021) Guerra-Urzola, R., K. Van Deun, J. C. Vera, and K. Sijtsma (2021).  A guide for sparse PCA: model comparison and applications.  Psychometrika 86, 893–919.
  * Johnstone (2001) Johnstone, I. M. (2001).  On the distribution of the largest eigenvalue in principal components analysis.  The Annals of Statistics 29(2), 295–327.
  * Jolliffe et al. (2003) Jolliffe, I. T., N. T. Trendafilov, and M. Uddin (2003).  A modified principal component technique based on the lasso.  Journal of Computational and Graphical Statistics 12(3), 531–547.
  * Karoui (2008) Karoui, N. E. (2008).  Operator norm consistent estimation of large-dimensional sparse covariance matrices.  The Annals of Statistics 36(6), 2717–2756.
  * Kleiber and Zeileis (2008) Kleiber, C. and A. Zeileis (2008).  AER: Applied Econometrics with R.  R package version 1.2-10.
  * Ledoit and Wolf (2004) Ledoit, O. and M. Wolf (2004).  A well-conditioned estimator for large-dimensional covariance matrices.  Journal of Multivariate Analysis 88(2), 365–411.
  * Mankiw et al. (1992) Mankiw, N. G., D. Romer, and D. Weil (1992).  A contribution to the empirics of economic growth.  Quarterly Journal of Economics 107(May), 407–437.
  * Mardia et al. (1979) Mardia, K., J. Kent, and J. Bibby (1979).  Multivariate Analysis (first ed.).  San Diego: Academic Press.
  * McCabe (1984) McCabe, G. P. (1984).  Principal variables.  Technometrics 26(2), 137–44.
  * Pavlenko et al. (2012) Pavlenko, T., A. Björkström, and A. Tillander (2012).  Covariance structure approximation via glasso in high-dimensional supervised classification.  Journal of Applied Statistics 39(8), 1643–1666.
  * Qi et al. (2013) Qi, X., R. Luo, and H. Zhao (2013).  Sparse principal component analysis by choice of norm.  Journal of Multivariate Analysis 114, 127–160.
  * Shen and Huang (2008) Shen, H. and J. Z. Huang (2008).  Sparse principal component analysis via regularized low rank matrix approximation.  Journal of Multivariate Analysis 99(6), 1015–1034.
  * Srivastava and Reid (2012) Srivastava, M. S. and N. Reid (2012).  Testing the structure of the covariance matrix with fewer observations than the dimension.  Journal of Multivariate Analysis 112, 156–171.
  * Tan et al. (2015) Tan, K. M., D. Witten, and A. Shojaie (2015).  The cluster graphical lasso for improved estimation of gaussian graphical models.  Computational Statistics & Data Analysis 85, 23–36.
  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society. Series B 58(1), 267–288.
  * Tibshirani (2019) Tibshirani, R. J. (2019).  bootstrap: Functions for the Book ”An Introduction to the Bootstrap”.  R package version 2019.6.
  * Witten et al. (2009) Witten, D. M., R. Tibshirani, and T. A. Hastie (2009).  A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis.  Biostatistics 10(3), 515–534.
  * Yamada et al. (2017) Yamada, Y., M. Hyodo, and T. Nishiyama (2017).  Testing block-diagonal covariance structure for high-dimensional data under non-normality.  Journal of Multivariate Analysis 155, 305–316.
  * Zou and Hastie (2005) Zou, H. and T. Hastie (2005).  Regularization and variable selection via the elastic net.  Journal of the Royal Statistical Society. Series B 67(2), 301–320.
  * Zou et al. (2006) Zou, H., T. Hastie, and R. Tibshirani (2006).  Sparse principal component analysis.  Journal of Computational and Graphical Statistics 15(2), 265–286.