# Divisive Hierarchical Clustering of Variables Identified by Singular Vectors

Jan O. Bauer  Department of Econometrics and Data Science, Vrije Universiteit
Amsterdam, Amsterdam, The Netherlands

###### Abstract

In this work, we present a novel method for divisive hierarchical variable
clustering. A cluster is a group of elements that exhibit higher similarity
among themselves than to elements outside this cluster. The correlation
coefficient serves as a natural measure to assess the similarity of variables.
This means that in a correlation matrix, a cluster is represented by a block
of variables with greater internal than external correlation. Our approach
provides a nonparametric solution to identify such block structures in the
correlation matrix using singular vectors of the underlying data matrix. When
divisively clustering $p$ variables, there are $2^{p-1}$ possible splits.
Using the singular vectors for cluster identification, we can effectively
reduce these number to at most $p(p-1)$, thereby making it computationally
efficient. We elaborate on the methodology and outline the incorporation of
dissimilarity measures and linkage functions to assess distances between
clusters. Additionally, we demonstrate that these distances are ultrametric,
ensuring that the resulting hierarchical cluster structure can be uniquely
represented by a dendrogram, with the heights of the dendrogram being
interpretable. To validate the efficiency of our method, we perform simulation
studies and analyze real world data on personality traits and cognitive
abilities. Supplementary materials for this article can be accessed online.

###### keywords:

Clustering , Dependence , Dissimilarity measure , Hierarchical variable
clustering , Singular value decomposition

##  1 Introduction

The study of multidimensional phenomena often requires the exploration of
hierarchical structures to represent the underlying concepts. In the field of
psychometric research, several examples of such phenomena have been described.
For instance, Carroll’s three-stratum theory (Carroll, 1993) and the Cattell-
Horn-Carroll theory, which combines Carroll’s three-stratum theory with Horn-
Cattell’s theory (McGrew, 2009), serve to structure cognitive abilities.
Likewise, personality traits, such as the Big Five, are organized into groups
(Digman, 1990).

As a motivating example, we look at well-known data on cognitive abilities
(Thurstone and Thurstone, 1941). In this data set, six cognitive abilities
(associative memory (M), verbal comprehension (V), reasoning (R), number (N),
word fluency (W), and spatial ability (S)) were measured by 17 different
ability tests on children. The correlation matrix between these tests given in
Figure 1 visually represents the cognitive abilities by its approximate block
diagonal shape. This shape is particularly noticeable for the variables Flags,
Figures and Cards, which together represent spatial ability (S): The three
variables show strong internal correlations and weaker correlations with all
other variables. This distinct pattern is visible through a block of these
variables on the lower right diagonal of the correlation matrix. While the
other cognitive abilities are also represented by blocks along the diagonal,
their shapes are comparatively more masked than those of spatial ability (S).

Figure 1: Correlation matrix of the cognitive abilities data, which contains
17 ability tests for six cognitive abilities (associative memory (M), verbal
comprehension (V), reasoning (R), number (N), word fluency (W), and spatial
ability (S)) on 212 children.

Hierarchical variable clustering can be applied to identify these six
cognitive abilities. In addition, the purpose of the method extends to answer
questions such as “How closely are the cognitive abilities verbal
comprehension (V) and word fluency (W) related?” or “What is the relationship
between the tests Flags, Figures, and Cards within the spatial ability (S)?”.
Essentially, hierarchical variable clustering uncovers the complete
interdependency structure between all variables.

Several methodologies have been proposed to model and analyze hierarchical
structures in data in order to answer such questions. Holzinger (1944) and
Jöreskog (1969) derived factor models for cases where the number of factors
and relations between groups are known. Cavicchia and Vichi (2022) introduced
a latent factor model which aims to model an unknown hierarchical structure of
order two. Further, Cavicchia et al. (2022) proposed a hierarchical principal
component analysis model consisting of disjoint groups of observed variables,
while Vichi and Saporta (2009) presented a constrained principal component
analysis for simultaneous clustering of objects and partitioning of variables.
However, all these methods are limited to a hierarchical structure of second
order and do not reveal the complete hierarchical structure.

In general, there are two main approaches for hierarchical clustering:
agglomerative (”bottom-up”) and divisive (”top-down”) clustering. For
agglomerative clustering, each element is considered as one cluster (bottom)
and the elements are clustered iteratively until all elements are contained in
a single cluster (top). Contrarily, divisive clustering starts with a single
cluster containing all elements (top) and iteratively splits one cluster into
two, until each element represents one cluster (down).

Vigneau and Qannari (2003); Vigneau et al. (2015) proposed an agglomerative
approach for grouping variables based on correlation coefficients, and
Kojadinovic (2010) derived a corpula-based clustering method. For applications
in portfolio management and neuroscience, Mantegna (1999), Liu et al. (2012),
and de Prado (2016) transformed the correlation matrices into dissimilarity
matrices for agglomerative hierarchcical variable clustering. This way,
observation-clustering methods can be applied to identify patterns and
groupings using this distance matrix. Revelle (1979) suggested a clustering
approach based on internal consistency reliability, where variables enter a
cluster if this increases reliability. However, agglomorative methods are
sensitive to misspecification by construction (Kaufman and Rousseeuw, 1990).
Therefore, other approaches are required.

Building on the motivational example provided above, as illustrated by the
correlation matrix of the Bechtoldt data in Figure 1, clustered variables are
characterized by an approximate block diagonal shape of the correlation
matrix. Cavicchia et al. (2020) exploited this characteristic and identified
hierarchical clusters by finding a so-called ultrametric correlation matrix
that fits the underlying data. Their approach is neither agglomerative nor
divisive, as it aims to determine the hierarchical structure directly through
the ultrametric correlation matrix.

However, apart from their work, no efforts have been made for non-
agglomerative hierarchical variable clustering. To address this gap, we
propose a nonparametric divisive clustering method. Our approach is based on
the singular vectors of the data matrix. It holds that these singular vectors
mirror the structure of the sample correlation matrix and can therefore be
exploited to detect the variable clusters.

The rest of this article is organized as follows: In Section 2, after
introducing basic notations and definitions, we present the concept for
divisive hierarchical variable clustering based on singular vectors. We
provide simulation studies in Section 3, and real data examples in Section 4.
Section 5 contains a discussion.

##  2 Singular vectors for identifying variable clusters

###  2.1 Setup

Let $\bm{X}$ be an $n\times p$ matrix of data that contains $n$ realizations
of each variable $X_{1},\ldots,X_{p}$. The singular value decomposition (SVD)
of the data matrix can be written as follows:

$\bm{X}=\bm{U}\bm{\Sigma}\bm{V}^{T},\;\bm{U}^{T}\bm{U}=\bm{I}_{n},\;\bm{V}^{T}\bm{V}=\bm{I}_{p}\,,$

where $\bm{\Sigma}$ contains the singular values of $\bm{X}$ on the diagonal,
and $\bm{U}=(\bm{u}_{1},\ldots,\bm{u}_{n})$ and
$\bm{V}=(\bm{v}_{1},\ldots,\bm{v}_{p})$ are the left and right singular
vectors respectively. Without loss of generality, assume that $\bm{X}$ is
standardized so that each column has an overall mean of zero and a standard
deviation of one. This implies that the right singular vectors of the data
matrix are the eigenvectors of the sample correlation matrix.

In divisive hierarchical clustering, a cluster is recursively split into two
disjoint subclusters, starting from the set of all objects. Divisive
hierarchical clustering of the variables $X_{1},\ldots,X_{p}$ is therefore a
top-down approach that can be represented by a sequence of $p$ partitions
$P^{(1)},\ldots,P^{(p)}$. Here, $P^{(1)}=\\{\\{X_{1},\ldots,X_{p}\\}\\}$
contains all variables as a single cluster (top), and
$P^{(p)}=\\{X_{1},\ldots,X_{p}\\}$ is the partition in which each variable
forms its separate cluster (down). Each partition $P^{(h)}$ is a set of
disjoint clusters $G^{(1)},\ldots,G^{(h)}$ with $G^{(i)}\neq\emptyset$,
$\bigcup_{i}G^{(i)}=\\{X_{1},\ldots,X_{p}\\}$, and $G^{(i)}\cap
G^{(j)}=\emptyset$ for all $i\in\\{1,\ldots,h\\}$ and $i\neq j$. Further, the
number of variables contained in each cluster $p_{i}=|G^{(i)}|$ is given by
its cardinality. For divisive clustering, $P^{(h+1)}$ is constructed by
splitting one cluster in $P^{(h)}$. More precisely, at each stage,
hierarchical variable clustering using singular vectors (HC-SVD) splits
according to a criterion any cluster $G^{(i)}$ contained in the partition
$P^{(h)}$ into two disjoint clusters $G^{(j)}$ and $G^{(j^{\prime})}$ to
construct the new partition $P^{(h+1)}$ that contains $G^{(j)}$ and
$G^{(j^{\prime})}$ instead of $G^{(i)}$. As the clusters are disjoint, it is
irrelevant which cluster is selected for the split.

In order to decide how a cluster should be split, a measure of dissimilarity
is required. For variable clustering, a natural measure that describes
dissimilarity between variables is given by the correlation coefficient, and
we therefore use correlation-based measures in this work. We split a cluster
$G^{(i)}$ in such a way that the distance between the variables contained in
the two splits $G^{(j)}$ and $G^{(j^{\prime})}$ is maximized.

###  2.2 Basic idea

For now, we only consider splitting one cluster $G^{(i)}$ into $G^{(j)}$ and
$G^{(j^{\prime})}$. Without loss of generality, we consider $G^{(i)}=P^{(1)}$,
i.e., we consider the first split of the divisive approach where all variables
$X_{1},\ldots,X_{p}$ are split into two disjoint clusters, meaning the
partition of the data matrix $\bm{X}$ into two distinct submatrices
$\bm{X}_{j}$ and $\bm{X}_{j^{\prime}}$ where rows are the observations and
columns are the variables contained in the clusters $G^{(j)}$ and
$G^{(j^{\prime})}$ respectively. We assume the convenient ordering
$\bm{X}=(\bm{X}_{j},\bm{X}_{j^{\prime}})$ since this order can be obtained
using column permutation.

The sample correlation matrix of $\bm{X}$ can be expressed as

$\bm{R}=\begin{pmatrix}\bm{R}_{j}&\bm{R}_{jj^{\prime}}\\\
\bm{R}_{jj^{\prime}}^{T}&\bm{R}_{j^{\prime}}\end{pmatrix}\;,$ (1)

where $\bm{R}_{j}$ is the correlation matrix of $\bm{X}_{j}$,
$\bm{R}_{j^{\prime}}$ is the correlation matrix of $\bm{X}_{j^{\prime}}$, and
$\bm{R}_{jj^{\prime}}$ contains the correlations between these two
submatrices. The objective is to find the clusters $G^{(j)}$ and
$G^{(j^{\prime})}$ with submatrices $\bm{X}_{j}$ and $\bm{X}_{j^{\prime}}$ in
such a way that the distance between these clusters is maximized, i.e., in
such a way that the correlation between the variables contained in these
clusters is minimal.

Let us assume for the moment that $\bm{R}_{jj^{\prime}}=\bm{0}$. Then
$\mathbf{R}$ has a block-diagonal structure and the clusters $G^{(j)}$ and
$G^{(j^{\prime})}$ are well separated because they are uncorrelated. Bauer
(2024) notes that the right singular vectors of $\bm{X}$ mirror the block
diagonal structure of the correlation matrix. For the purpose of completeness,
we provide the corresponding result.

######  Corollary 1.

The correlation matrix $\bm{R}$ is a block diagonal matrix if and only if the
right singular vectors of $\bm{X}$ exhibit the structure

$\bm{V}\bm{P}_{\pi}=\begin{pmatrix}\bm{V}_{j}&\bm{0}\\\
\bm{0}&\bm{V}_{j^{\prime}}\end{pmatrix}\,,$

where $\bm{V}_{j}$ are the right singular vectors of $\bm{X}_{j}$,
$\bm{V}_{j^{\prime}}$ are the right singular vectors of $\bm{X}_{j^{\prime}}$,
and $\bm{P}_{\pi}$ is a column permutation matrix leading to this block
diagonal structure.

Therefore, the singular vectors identify the two uncorrelated clusters.
Consequently, by analyzing the structure of the singular vectors, it is
possible to identify the clusters.

However, when $\bm{R}_{jj^{\prime}}\neq\bm{0}$, we conclude from the Davis-
Kahan theorem that the singular vectors are perturbed by
$\bm{R}_{jj^{\prime}}$ (Yu et al., 2015). Still, we can use the singular
vectors to find the underlying block diagonal structure of the correlation
matrix compromised by $\bm{R}_{jj^{\prime}}$ to detect the clusters. To do
this, we propose the use of sparse singular vectors. These so called sparse
loadings are sparse approximations, i.e., having many zero entries, of the
singular vectors. The intuition of using them is that by calculating sparse
loadings, the zero value components of the singular vectors perturbed by
$\bm{R}_{jj^{\prime}}$ are shrunk back to zero. Calculation of sparse loadings
is a well-established area in the literature for which numerous approaches
have been proposed (see, e.g., Zou et al. (2006), Shen and Huang (2008),
Witten et al. (2009), Yang et al. (2014), and Gataric et al. (2020) among
others). These methods typically employ a form of regularization, such as the
$\mathcal{L}_{1}$-type lasso constraint (Tibshirani, 1996), the hard
thresholding penalty (Donoho and Johnstone, 1994), or the smoothly clipped
absolute deviation penalty (Fan and Li, 2001). We consider the former to
obtain sparse loadings according to the following optimization problem:

$\min\limits_{\bm{u},\bm{v}}\|\bm{X}-\bm{u}\bm{v}^{T}\|_{F}^{2}\;\;{\rm
subject\;to}\;\;\|\bm{v}\|_{2}^{2}=1,\;\|\bm{v}\|_{1}\leq\alpha\,,$ (2)

which calculates the first sparse loading $\check{\bm{v}}_{1}$ with
$\|\check{\bm{v}}_{1}\|_{2}^{2}=1$ and $\mathcal{L}_{1}$ imposing
regularization parameter $\alpha$, and a vector $\check{\bm{u}}$ which can be
interpreted as containing the corresponding quasi singular value and quasi
left singular vector. Note that throughout this paper, $\|\bm{v}\|_{p}$
denotes the $\mathcal{L}_{p}$ norm for any vector $\bm{v}$. The remaining
sparse loadings for $i>1$ can then be calculated iteratively. For this the
data matrix $\bm{X}$ must be replaced by the residual matrices of the
sequential matrix approximations. We note that all blocks can be revealed
using only a few singular vectors since each block can be mirrored by a single
singular vector. Therefore, it is not necessary to calculate all $p$ sparse
singular vectors. This will be elaborated later in this work.

The steps for splitting any cluster $G^{(i)}$, i.e., of detecting the
variables belonging to the splits $G^{(j)}$ and $G^{(j^{\prime})}$ are
outlined in Algorithm 1, which illustrates the method for variable clustering
using singular vectors. Let therefore $s$ be the degree of sparsity for all
calculated sparse loadings, i.e. the number of non-zero components in each
sparse loading. Each sparse loading represents a split into two clusters: one
is represented by the non-zero components and the other by the components
being equal to zero. Therefore, $s$ lies between $1$ and $p_{i}-1$ where
$|G^{(i)}|=p_{i}$ are the number of variables contained in the cluster we want
to split.

Algorithm 1 Procedure to split any cluster $C^{(i)}$ into the clusters
$C^{(j)}$ and $C^{(j^{\prime})}$ using sparse loadings.

1:Input Cluster $G^{(i)}$ and corresponding standardized data matrix
$\bm{X}_{i}$.

2:for $s\in\\{1,\ldots,p_{i}-1\\}$ do

3: Calculate the first $k_{i}$ sparse loadings of the $n\times p_{i}$ data
matrix $\bm{X}_{i}$ which contains only the variables cointained in the
cluster $G^{(i)}$. Calculate each with degree of sparsity $s$.

4: Compute the $k_{i}$ distances between the potential splits indicated by the
zero and non-zero values of each sparse loading.

5:end for

6:Choose the clusters $G^{(j)}$ and $G^{(j^{\prime})}$ with the largest
distance between them among the $k_{i}(p-1)$ calculated distances.

7:Output Splits $G^{(j)}$ and $G^{(j^{\prime})}$ of $G^{(i)}$.

In order to split a cluster $G^{(i)}$ optimally, one has to select $G^{(j)}$
and $G^{(j^{\prime})}$ amongst $2^{p_{i}-1}$ possible partitions of $G^{(i)}$
which are too many to be feasible in application. Using sparse loadings for
preselection, however, we can see in Algorithm 1 that this number is reduced.

######  Remark 1.

The number of splits to check when splitting a cluster $G^{(i)}$ with
$|G^{(i)}|=p_{i}$ variables using sparse loadings as described in Algorithm 1
is reduced from $2^{p_{i}-1}$ to $k_{i}(p_{i}-1)$ with $k_{i}\leq p_{i}$.

Strictly speaking, HC-SVD decreases the number of potential splits to check
only if $p_{i}>6$. However, for small numbers such as $p_{i}\leq 5$, checking
all these potential splits is computationally feasible regardless.

###  2.3 Number of sparse loadings

Using sparse loadings, we reduce the number of splits to check. We now want to
find the smallest number $k_{i}$ of sparse loadings required for hierarchical
variable clustering. To motivate our approach, we again consider the split of
$G^{(i)}=P^{(1)}$, where the data matrix $\bm{X}$ is partitioned into two
distinct submatrices. Further, we assume that the sample correlation matrix of
$\bm{X}$ equals a block diagonal matrix with $b$ blocks. Consequently, the
right singular vectors $\bm{V}\bm{P}_{\pi}$ follow the same block diagonal
structure when the columns are permuted using a permutation matrix
$\bm{P}_{\pi}$ (see Corollary 1 in Bauer (2024)):

$\bm{R}=\begin{pmatrix}\bm{R}_{1}&\bm{0}&\bm{0}\\\ \bm{0}&\ddots&\bm{0}\\\
\bm{0}&\bm{0}&\bm{R}_{b}\end{pmatrix}\,,\qquad\bm{V}\bm{P}_{\pi}=\begin{pmatrix}\bm{V}_{1}&\bm{0}&\bm{0}\\\
\bm{0}&\ddots&\bm{0}\\\ \bm{0}&\bm{0}&\bm{V}_{b}\end{pmatrix}\,.$

Additionally, we assume that the optimal split of $G^{(i)}$ is that $G^{(j)}$
contains the variables associated with the submatrix $\bm{R}_{1}$, and
$G^{(j^{\prime})}$ the remaining variables. The first singular vector
contained in $\bm{V}_{1}$, given as the first column vector of
$\bm{V}\bm{P}_{\pi}$ above, identifies the optimal split. However, since we do
not know the permutation matrix $\bm{P}_{\pi}$, and since the singular vectors
are usually ordered according to the magnitude of the respective singular
values, we do not know the index of the first singular vector contained in
$\bm{V}_{1}$ that reflects the cluster $G^{(j)}$ among all ordered singular
vectors. Still, we can assess the minimal number of singular vectors required
to identify each block and therefore also the optimal split.

######  Remark 2.

Let the standardized data matrix $\bm{X}=\bm{U}\bm{\Sigma}\bm{V}^{T}$ have a
block diagonal correlation matrix
$\bm{R}=\operatorname{diag}(\bm{R}_{1},\ldots,\bm{R}_{b})$ with $b$ blocks,
and denote the ordered eigenvalues of $\bm{R}$ by
$\lambda_{1}\geq\cdots\geq\lambda_{p}$. If

$k=\operatorname*{arg\,min}_{j\in\\{1,\ldots,p\\}}\\{\lambda_{j}:\lambda_{j}\geq
1\\}\,,$

each block $\bm{R}_{1},\ldots,\bm{R}_{b}$ is reflected by at least one of the
singular vectors $\bm{v}_{1},\ldots,\bm{v}_{k}$.

For any cluster $G^{(i)}$ represented by a submatrix $\bm{X}_{i}$, we denote
by $k_{i}$ the number of considered singular vectors according to Remark 2. In
practice, it is unlikely that the sample correlation matrix $\bm{R}$ equals a
block diagonal matrix exactly. Still, the intuition that each block is linked
to at least one singular vector with eigenvalue larger or equal to one remains
reasonable. Therefore, we suggest to follow the intuition provided in Remark 2
and use the number of eigenvalues larger or equal to one as a reference for
the number of sparse loadings for HC-SVD. It is worth noting that even if all
$p_{i}$ sparse loadings are selected for splitting a cluster $G^{(}{i})$, the
number of possible splits is reduced from $2^{p_{i}-1}$ to the quadratic
number $p_{i}(p_{i}-1)$, which is feasible for data with a moderate number of
variables.

We note that the recommended number of sparse loadings in Remark 2 is
equivalent to the Kaiser criterion (Kaiser, 1960), which can be used to
determine the number of factors in factor analysis. Despite this equivalence,
the underlying considerations for these two concepts differ: The Kaiser
criterion aims to select all singular vectors with a certain explanatory power
in terms of variance, while our criterion selects as many singular vectors as
necessary to ensure a representation of each cluster by at least one singular
vector.

###  2.4 Distance for hierarchical variable clustering

We recap that the optimal splits $G^{(j)}$ and $G^{(j^{\prime})}$ of a cluster
$G^{(i)}$ are those which have the smallest distance between them. Therefore,
we require a measures to assess distances.

There are several measures to assess dissimilarities between sets of random
variables. The RV-coefficient (Escoufier, 1973; Robert and Escoufier, 1976) is
an extension of the squared correlation coefficient to random vectors and
therefore measures closeness between them. The RV-coefficient has been
adjusted in recent years (see, e.g., Smilde et al. (2008) and Mordant and
Segers (2022) among others). Székely et al. (2007) proposed the distance
correlation that measures the dependency between two random vectors and equals
zero if and only if the random vectors are independent. Further, there are
several dependency measures for the ”$p>n$” case (see, e.g., Székely and Rizzo
(2013), Pan et al. (2020), and Zhu et al. (2020) among others).

The RV-coefficient proved to be the most successful during our research and
benefits from a low computational cost. Therefore, we only discuss this
dissimilarity measure in what follows. In addition, we can incorporate linkage
functions traditionally used in clustering observations into the scope of
clustering variables.

######  Theorem 1.

Let $\bm{X}_{i}$ be a standardized data matrix containing the variables of the
cluster $G^{(i)}$ as columns, and let the correlation matrix be denoted
according to the splits $G^{(j)}$ and $G^{(j^{\prime})}$ as in (1). The
function $d_{jj^{\prime}}$ induced by the RV-coefficient

$d_{jj^{\prime}}=1-\frac{\|\bm{R}_{jj^{\prime}}\|_{F}^{2}}{\|\bm{R}_{j}\|_{F}\|\bm{R}_{j^{\prime}}\|_{F}}\quad\text{(RV-
coefficient)}$

defines a semidistance. If there is no perfect collinearity between two
variables, $d_{jj^{\prime}}$ induced by average or complete linkage

$\displaystyle
d_{jj^{\prime}}=\begin{cases}1-\dfrac{\|\operatorname*{vec}(\bm{R}_{jj^{\prime}})\|_{1}}{p_{j}p_{j^{\prime}}}&j\neq
j^{\prime}\\\ 0&j=j^{\prime}\end{cases}\quad\text{(average linkage)}$
$\displaystyle
d_{jj^{\prime}}=1-\|\operatorname*{vec}(\bm{R}_{jj^{\prime}})\|_{\infty}\quad\text{(single
linkage)}$

defines a semidistance.

Intuitively, clustering perfectly collinear variables, i.e., variables that
have a correlation of one in absolute terms, does not yield additional
insights and can be addressed during data preprocessing. From a technical
perspective, the positivity property for correlation based distances is
violated in the presence of perfect collinearity.

###  2.5 Hierarchical variable clustering using singular vectors

For the complete procedure of HC-SVD, it is left to obtain a $p\times p$
distance matrix $\bm{M}=(m_{ij})$ that contains the distances between
variables. This is obtained by applying Algorithm 1 iteratively on each
partition $P^{(h)}$, as described in Algorithm 2.

Algorithm 2 Procedure for HC-SVD on a data matrix $\bm{X}$

1:Input Standardized data matrix $\bm{X}$ of observed random variables
$X_{1},\ldots,X_{p}$.

2:Let $P^{(1)}=\\{\\{X_{1},\ldots,X_{p}\\}\\}$ be the partition containing all
random variables as a single cluster, and let $\bm{M}$ be a $p\times p$ matrix
with zeros on the diagonal.

3:for $h\in\\{1,\ldots,p-1\\}$ do

4: Obtain $P^{(h+1)}$ by splitting any cluster $G^{(i)}$ in $P^{(h)}$
according to Algorithm 1

5: into $G^{(j)}$ and $G^{(j^{\prime})}$.

6: Then assign the distances between the variables of $G^{(j)}$ and
$G^{(j^{\prime})}$ by $d_{jj^{\prime}}$:

7: $m_{kl}=m_{lk}=d_{jj^{\prime}}$ for all $X_{k}\in G^{(j)}$ and $X_{l}\in
G^{(j^{\prime})}$.

8:end for

9:Output Distance matrix $\bm{M}$ of HC-SVD.

For the very first split when obtaining two clusters $G^{(j)}$ and
$G^{(j^{\prime})}$ from $P^{(1)}$ to obtain the partition $P^{(2)}$, the
distances between the variables contained in $G^{(j)}$ and $G^{(j^{\prime})}$
are assigned by $d_{jj^{\prime}}$. In the next step, when the cluster, say,
$G^{(j)}$ is split into $G^{(l)}$ and $G^{(l^{\prime})}$, the distances
between the corresponding variables is updated to $d_{ll^{\prime}}$ while the
distance to the variables contained in $G^{(j^{\prime})}$ remains unchanged.
The dissimilarity between variables is therefore refined iteratively as the
algorithm proceeds to reveal the hierarchical structure.

It is desirable that the resulting distance matrix is an ultrametric distance
matrix because this ensures that the distance between clusters decreases for
each split, making the visualization of the resulting dendrogram interpretable
(Johnson, 1967; Hartigan, 1967). A distance matrix $\bm{M}=(m_{ij})$ is
ultrametric if

$m_{ij}\leq\max\\{m_{il},m_{lj}\\}\,,$

for all $i$, $j$, and $l$.

######  Remark 3.

Calculating the distances induced by the RV-coefficient, average linkage, or
single linkage as in Theorem 1, the distance matrix $\bm{M}$ resulting from
HC-SVD is ultrametric if the procedure found the correct hierarchical cluster
structure.

Given that HC-SVD only considers $k_{i}(p_{i}-1)$ splits instead of all
$2^{p_{i}-1}$ possibilities for a cluster $G^{(i)}$, there is a chance of
overlooking the true underlying cluster structure. This trade-off is made in
favor of computational feasibility, while acknowledging that satisfaction of
the ultrametric property cannot be guaranteed in such cases.

##  3 Simulation study

###  3.1 Simulation design

In this section, we evaluate the performance of HC-SVD as outlined in
Algorithm 2 for the induced distance functions derived from Theorem 1. The
sparse loadings are computed using the method of Shen and Huang (2008) which
is implemented in the irlba package (Baglama et al., 2022) within the
statistical software R (R Core Team, 2022).

We simulate data matrices containing $n$ observations from a $p$-multivariate
normal distribution $N(\bm{0},\bm{\mathrm{P}})$, where $\bm{\mathrm{P}}$ is
the known population correlation matrix that exhibits a hierarchical
structure. The data matrices are simulated $100$ times and clustered for each
simulation. We consider two different designs to evaluate the performance of
the method: one characterized by a small number of clusters with deep nesting
(a), and the other with many clusters and a flatter structure (b).

  1. (a)

We construct a design where all cluster compositions are known for $5$, $10$,
$15$, and $20$ clusters. The population correlation matrix exhibits a block
diagonal structure
$\bm{\mathrm{P}}=\operatorname{diag}(\bm{\mathrm{P}}_{1},\ldots,\bm{\mathrm{P}}_{b})$
with $b=5$ blocks, each of size $20\times 20$. In our terminology, these
blocks are referred to as clusters. Each of these clusters contains 20
variables and exhibits a stepwise correlation structure with correlation steps
$0.2+\eta_{1}$, $0.4+\eta_{2}$, $0.6+\eta_{3}$, and $0.8+\eta_{4}$ for
subclusters of four variables, while the main block diagonal maintains a
correlation of $0.95$. Here, $\eta_{1},\ldots,\eta_{4}\sim U(-0.05,0.05)$.
Figure 2 provides a visual representation of the dendrogram to illustrate the
hierarchical structure. In this simulation design, we consider
$p\in\\{100,200\\}$ and $n=300$.

  2. (b)

We construct a design where all cluster compositions are known for $20$ and
$40$ clusters. The population correlation matrix exhibits a block diagonal
structure
$\bm{\mathrm{P}}=\operatorname{diag}(\bm{\mathrm{P}}_{1},\ldots,\bm{\mathrm{P}}_{b})$
with $b=20$ blocks, each of size $3\times 3$. All these blocks exhibit a
symmetric structure with ones on the diagonal, upper triangle components given
as $\rho_{ij}=(-1)^{j}\eta_{l}^{(j-1)^{2}}$ when $i<j$, with $\eta_{l}\sim
U(0.8,0.9)$ for $l\in\\{1,\ldots,b\\}$, and lower triangle matrix following
the same pattern. Figure 2 provides a visual representation of the dendrogram
to illustrate the hierarchical structure. In this simulation design, we
consider $p\in\\{60,120\\}$ and $n=180$.

Figure 2: The dendrograms corresponding to the population correlation matrix
$\bm{\mathrm{P}}$ for design (a) (left) with $p=100$ variables, and for design
(b) (right) with $p=60$ variables. In design (a), the variable grouping is
specified by construction for 5, 10, 15, and 20 clusters. Similarly, in design
(b), the variable grouping is known for 20 and 40 clusters.

As outlined in the introduction of this paper, Cavicchia et al. (2020) have,
to our knowledge, introduced the only non-agglomerative variable clustering
approach to date. Therefore, we compare HC-SVD with their ultrametric
correlation (ULT-COR) method. Additionally, we implement divisive analysis
clustering (DIANA), introduced by Kaufman and Rousseeuw (1990) in the context
of observation clustering, for variable clustering. This can be accomplished
by treating the correlation matrix with absolute values as a distance matrix.
DIANA is available in the cluster package (Maechler et al., 2022) within the
statistical software R. The performance of the methods is evaluated using the
adjusted Rand index (Rand, 1971; Hubert and Arabie, 1985). All computational
results were obtained using R 4.1.3 on a PC running macOS version 14.0.

An intuitive question is the usage of principal component analysis (PCA) for
identification of variable clusters, given its natural connection to HC-SVD.
However, Bauer (2024) demonstrated that PCA already performs poorly for
detecting uncorrelated variable clusters, making it not suitable for
hierarchical variable clustering. A detailed discussion on this matter is
provided in the referenced work.

###  3.2 Simulation results

The results for simulation designs (a) and (b) are illustrated in Figure 3 and
Figure 4 respectively. Exact numerical results are provided in the
supplementary material. Notably, the comparison with ULT-COR was limited to
cases where $p\leq 100$ due to the long running time of the method in higher
dimensions. In fact, the method did not converge in our simulations when
$p>100$.

HC-SVD with any of the three induced distance functions average linkage,
single linkage, or RV coefficient, outperforms the other two methods across
the simulation designs. It is worth noting that single linkage shows the best
performance, although only slightly, in comparison to the other induced
distance functions.

Figure 3: Simulation design (a): Performance of HC-SVD induced by average
linkage, single linkage and the RV coefficient, DIANA, and ULT-COR for
$p\in\\{100,200\\}$ and $n=300$ measured by ARI. ULT-COR did not converge for
$p=200$ and is therefore excluded from the comparison. Figure 4: Simulation
design (b): Performance of HC-SVD induced by average linkage, single linkage
and the RV coefficient, DIANA, and ULT-COR for $p\in\\{60,120\\}$ and $n=180$
measured by ARI. ULT-COR did not converge for $p=120$ and is therefore
excluded from the comparison.

###  3.3 Computation time

The average computation times across the $100$ simulations for each simulation
study are illustrated in Figure 5. HC-SVD and DIANA exhibit comparable running
times. The computational efficiency of HC-SVD can be further enhanced when the
calculations of the distance between clusters detected by the sparse loadings
in Algorithm 1 are performed in parallel across multiple cores. It is
important to note that, in this simulation study, the calculations were
executed on a single core only.

Figure 5: Average computation time of the models for all $100$ simulations of
the four simulation studies.

##  4 Real data examples

###  4.1 Big Five personality traits

In this section, we demonstrate the application of HC-SVD on personality
traits. Specifically, we analyze a data set comprising 25 personality items,
categorized into five groups: neuroticism (N), extraversion (E), openness to
experience (O), agreeableness (A), and conscientiousness (C), according to the
NEO personality inventory for personality psychology (Costa and McCrae, 1978;
McCrae et al., 2005). These personality traits align with the Big Five
personality traits (Roccas et al., 2002). The data set contains self-
assessments from 190 university students who indicate the extent to which they
embody the personalities described by the 25 items. Data can be accessed
online from Osaka University at http://bstat.jp/en material/.

Figure 6: Sample correlation matrix of the Big Five personality traits data
with corresponding hierarchical cluster structure, which was identified by HC-
SVD using single linkage and is illustrated as a dendrogram.

Figure 6 presents both the sample correlation matrix of the data and the
hierarchical cluster structure identified by the HC-SVD using single linkage,
as illustrated by the accompanying dendrogram. HC-SVD successfully identifies
the five personality traits, which can be seen in the arrangement of the
groups agreeableness (A), conscientiousness (C), extraversion (E), neuroticism
(N) and openness to experience (O) from left to right in the dendrogram. In
addition, nuances of the five personality traits can be obtained. For example,
neuroticism (N) decomposes into the personality items {Unrest}, {Pessimistic,
Careful}, and {Worry, Sensitive}.

In hierarchical clustering, the typical assignment for the value $h_{i}$ of a
cluster $G^{(i)}$, representing the height of this cluster in the associated
dendrogram, is the distance between its splits $G^{(j)}$ and
$G^{(j^{\prime})}$, given by $h_{i}=d_{jj^{\prime}}$. In psychometric
applications, however, the focus may shift to internal consistency within
clusters rather than correlation between them. The evaluation of internal
consistency involves reliability measures such as Cronbach’s alpha (Cronbach,
1951), Revelle’s beta, or McDonald’s omega (see, e.g., Zinbarg et al. (2005)
for an overview). The evaluation of the internal cluster correlation follows a
similar rationale to the evaluation of the correlation between clusters, and
to examine the internal consistency, the value $h_{i}$ of a cluster can be
assessed using a reliability measure.

Figure 7 (right) shows the identical hierarchical variable structure
identified by HC-SVD as shown in Figure 6. In this illustration, however, the
cluster values $h_{i}$ are additionally assigned by
$h_{i}=1-\|\bm{R}_{i}\|_{2}/p_{i}$ (left). We choose the ”average” spectral
norm, as it naturally serves as an indicator for the internal consistency of a
matrix due to the characteristic of the singular value decomposition.
Specifically, the larger the common variance among variables, the more the
correlation matrix stretches in the direction of the first singular vector
(Bauer and Drabant, 2021). This selection is further motivated by the
satisfaction of the ultrametric property under mild assumptions. For
completeness, detailed discussions on the ultrametric property are provided in
the supplementary material. To our knowledge, other reliability measures do
not exhibit the ultrametric property. In addition, these measures have
limitations in terms of their underlying assumptions and calculation aspects,
as discussed by Sijtsma (2009), Dunn et al. (2014), Hayes and Coutts (2020),
and others.

Figure 7: Dendrograms of the Big Five personality traits with different
assignments for the cluster height: $h_{i}=1-\|\bm{R}_{i}\|_{2}/p_{i}$ (left)
and single linkage (right).

Using internal consistency as the cluster height shows that the data exhibits
five clusters with comparable internal consistency, each representing a Big
Five personality trait.

###  4.2 Cognitive abilities

In this section, we apply HC-SVD on the the cognitive abilities presented in
the introduction of this work. In this data set, six cognitive abilities
(associative memory (M), verbal comprehension (V), reasoning (R), number (N),
word fluency (W), and spatial ability (S)) were measured by 17 different
ability tests on 212 children (Thurstone and Thurstone, 1941; Bechtoldt,
1961). The corresponding sample correlation matrix is available in the CRAN
contributed package psych (Revelle, 2023). We visualize the outcome using a
dendrogram along with the sample correlation matrix in Figure 8. Single
linkage was chosen as the distance measure in HC-SVD.

Figure 8: Sample correlation matrix of the cognitive abilities data with
corresponding hierarchical cluster structure, which was identified by HC-SVD
using single linkage and is illustrated as a dendrogram.

The method successfully identifies the six cognitive abilities. Building on
the motivation from the previous example, we assigned the value of each
cluster using the spectral norm, i.e., $h_{i}=1-\|\bm{R}_{i}\|_{2}/p_{i}$. The
corresponding dendrogram is given in Figure 9. Cutting the dendrogram into six
clusters (left), we identify the six cognitive abilities, all of which have
comparable internal consistency, with the exception of verbal comprehension
(V), which has greater internal consistency, i.e., a lower level in the
dendrogram. However, we note that the internal consistency of both verbal
comprehension (V) and reasoning (R) together is large which is also indicated
by the strong correlation between these variables in the sample correlation
matrix (Figure 8). Consequently, an alternative interpretation of the data
emerges: the ability tests for verbal comprehension and reasoning in this
sample indicate a common cognitive ability rather than two distinct ones.
Cutting the dendrogram in five clusters (right) visualizes these five
abilities, with four abilities matching the original interpretation and an
additional fifth common ability for verbal comprehension and reasoning.

Figure 9: Dendrogram for the cognitive ability data identified by HC-SVD using
single linkage with two different horizontal cut-points. These cut-points
either identify the six cognitive abilities (left), or identify reasoning (R)
and verbal comprehension (V) as one common cluster (right).

##  5 Discussion

In this paper, we introduce a nonparametric method for divisive hierarchical
variable clustering. Our approach relies on sparse loadings to effectively
reduce the number of possible splits for a divisive approach, enabling fast
computation even for large data sets, as demonstrated in our simulation
studies.

We propose three distinct measures for dissimilarity to induce semidistances
between clusters. Furthermore, these semidistances satisfy the ultrametric
property if our proposed method identifies the true underlying hierarchical
variable structure. For completeness, we want to mention that in the course of
our research we also experimented with other dissimilarity measures, but they
either proved computationally costly or exhibited poor performance.

In two examples, we showcase the applicability of hierarchical variable
clustering using singular vectors on real data. In psychometric research, the
interest often lies in the internal reliability of variable groups. We
therefore illustrate that internal reliability, calculated for each cluster by
an averaged spectral norm of the respective sample correlation matrix,
provides interpretable values for the clusters.

Future research in the area of factor analysis can utilize the cluster values
assigned by internal reliability. In factor analysis, the aim is to identify
latent variables that account for most of the variance in the data and are a
combination of the observed variables. These factors are easier to interpret
if they are disjoint. Since the clusters in hierarchical clustering are
disjoint and their values are determined by a function of the largest
eigenvalues, which imply the explained variance, there is a natural link
between HC-SVD and factor analysis.

## Supplementary Material

The online supplementary material contains the appendix with derivations and
proofs, discussions about the ultrametric property of
$h_{i}=1-\|\bm{R}_{i}\|_{2}/p_{i}$, the exact numerical results together with
the number of sparse loadings
$k_{i}=\operatorname*{arg\,min}_{j\in\\{1,\ldots,p\\}}\\{\lambda_{j}:\lambda_{j}\geq
1\\}$ from Remark 2 used for the simulation studies presented in in this work,
and R code to perform all simulation studies.

## Acknowledgements

I am grateful to Carlo Cavicchia for providing R code for the ultrametric
correlation model.

## References

  * Baglama et al. (2022) Baglama, J., Reichel, L., Lewis, B.W., 2022.  irlba: Fast Truncated Singular Value Decomposition and Principal Components Analysis for Large Dense and Sparse Matrices.  R package version 2.3.5.1.
  * Bauer (2024) Bauer, J.O., 2024.  High-dimensional block diagonal covariance structure detection using singular vectors.  Preprint .
  * Bauer and Drabant (2021) Bauer, J.O., Drabant, B., 2021.  Principal loading analysis.  J. Multivar. Anal. 184.
  * Bechtoldt (1961) Bechtoldt, H.P., 1961.  An empirical study of the factor analysis stability hypothesis.  Psychometrika 26, 405–432.
  * Carroll (1993) Carroll, J.B., 1993.  Human Cognitive Abilities: A Survey of Factor-Analytic Studies.  Cambridge University Press, New York, NY, USA.
  * Cavicchia and Vichi (2022) Cavicchia, C., Vichi, M., 2022.  Second-order disjoint factor analysis.  Psychometrika 87, 289–309.
  * Cavicchia et al. (2020) Cavicchia, C., Vichi, M., Zaccaria, G., 2020.  The ultrametric correlation matrix for modelling hierarchical latent concepts.  AStA Adv. Stat. Anal. 14, 837–853.
  * Cavicchia et al. (2022) Cavicchia, C., Vichi, M., Zaccaria, G., 2022.  Hierarchical disjoint principal component analysis.  AStA Adv. Stat. Anal. .
  * Costa and McCrae (1978) Costa, P.T., McCrae, R.R., 1978.  Objective personality assessment, in: Storandt, M., Siegler, I.C., Elias, M.F. (Eds.), The Clinical Psychology of Aging. Springer US, Boston, MA, pp. 119–143.
  * Cronbach (1951) Cronbach, L., 1951.  Coefficient alpha and the internal structure of tests.  Psychometrika 16, 297–334.
  * de Prado (2016) de Prado, M.L., 2016.  Building diversified portfolios that outperform out of sample.  J. Portf. Manag. 42, 59–69.
  * Digman (1990) Digman, J.M., 1990.  Personality structure: Emergence of the five-factor model.  Annu. Rev. Psychol. 41, 417–440.
  * Donoho and Johnstone (1994) Donoho, D.L., Johnstone, I.M., 1994.  Ideal spatial adaptation by wavelet shrinkage.  Biometrika 81, 425–455.
  * Dunn et al. (2014) Dunn, T.J., Baguley, T.S., Brunsden, V., 2014.  From alpha to omega: A practical solution to the pervasive problem of internal consistency estimation.  Br. J. Psychol. 105 3, 399–412.
  * Escoufier (1973) Escoufier, Y., 1973.  Le traitement des variables vectorielles.  Biometrics 29, 751–760.
  * Fan and Li (2001) Fan, J., Li, R., 2001.  Variable selection via nonconcave penalized likelihood and its oracle properties.  J. Am. Stat. Assoc. 96, 1348–1360.
  * Gataric et al. (2020) Gataric, M., Wang, T., Samworth, R.J., 2020.  Sparse principal component analysis via axis-aligned random projections.  J. R. Stat. Soc. B 82, 329–359.
  * Hartigan (1967) Hartigan, J.A., 1967.  Representation of similarity matrices by trees.  J. Am. Stat. Assoc. 62, 1140–1158.
  * Hayes and Coutts (2020) Hayes, A.F., Coutts, J.J., 2020.  Use omega rather than Cronbach’s alpha for estimating reliability. But….  Commun. Methods Meas. 14, 1–24.
  * Holzinger (1944) Holzinger, K.J., 1944.  A simple method of factor analysis.  Psychometrika 9, 257–262.
  * Hubert and Arabie (1985) Hubert, L., Arabie, P., 1985.  Comparing partitions.  J. Classif. 2, 193–218.
  * Johnson (1967) Johnson, S.C., 1967.  Hierarchical clustering schemes.  Psychometrika 32, 241–254.
  * Jöreskog (1969) Jöreskog, K.G., 1969.  A general approach to confirmatory maximum-likelihood factor analysis.  Psychometrika 34, 183–202.
  * Kaiser (1960) Kaiser, H.F., 1960.  The application of electronic computers to factor analysis.  Educ. Psychol. Meas. 20, 141–151.
  * Kaufman and Rousseeuw (1990) Kaufman, L., Rousseeuw, P.J., 1990.  Finding Groups in Data: An Introduction to Cluster Analysis.  John Wiley, New York.
  * Kojadinovic (2010) Kojadinovic, I., 2010.  Hierarchical clustering of continuous variables based on the empirical copula process and permutation linkages.  Comput. Stat. Data Anal. 54, 90–108.
  * Liu et al. (2012) Liu, X., Zhu, X.H., Qiu, P., Chen, W., 2012.  A correlation-matrix-based hierarchical clustering method for functional connectivity analysis.  J. Neurosci. Methods 211, 94–102.
  * Maechler et al. (2022) Maechler, M., Rousseeuw, P., Struyf, A., Hubert, M., Hornik, K., 2022.  cluster: Cluster Analysis Basics and Extensions.  R package version 2.1.4.
  * Mantegna (1999) Mantegna, R.N., 1999.  Hierarchical structure in financial markets.  Eur. Phys. J. B 11, 193–197.
  * McCrae et al. (2005) McCrae, R.R., Costa, P.T., Jr., Martin, T.A., 2005.  The NEO–PI–3: A more readable revised NEO personality inventory.  J. Pers. Assess. 84, 261–270.
  * McGrew (2009) McGrew, K.S., 2009.  CHC theory and the human cognitive abilities project: Standing on the shoulders of the giants of psychometric intelligence research.  Intelligence 37, 1–10.
  * Mordant and Segers (2022) Mordant, G., Segers, J., 2022.  Measuring dependence between random vectors via optimal transport.  J. Multivar. Anal. 189.
  * Pan et al. (2020) Pan, W., Wang, X., Zhang, H., Zhu, H., Zhu, J., 2020.  Ball covariance: A generic measure of dependence in banach space.  J. Am. Stat. Assoc. 115, 307–317.
  * R Core Team (2022) R Core Team, 2022.  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing. Vienna, Austria.
  * Rand (1971) Rand, W.M., 1971.  Objective criteria for the evaluation of clustering methods.  J. Am. Stat. Assoc. 66, 846–850.
  * Rauf Ahmad (2019) Rauf Ahmad, M., 2019.  A significance test of the RV coefficient in high dimensions.  Comput. Stat. Data Anal. 131, 116–130.
  * Revelle (1979) Revelle, W., 1979.  Hierarchical cluster analysis and the internal structure of tests.  Multivar. Behav. Res. 14, 57–74.
  * Revelle (2023) Revelle, W., 2023.  psych: Procedures for Psychological, Psychometric, and Personality Research.  R package version 2.3.3.
  * Robert and Escoufier (1976) Robert, P., Escoufier, Y., 1976.  A unifying tool for linear multivariate statistical methods: The RV-coefficient.  J. R. Stat. Soc. C 25, 257–265.
  * Roccas et al. (2002) Roccas, S., Sagiv, L., Schwartz, S.H., Knafo, A., 2002.  The Big Five personality factors and personal values.  Pers. Soc. Psychol. Bull 28, 789–801.
  * Shen and Huang (2008) Shen, H., Huang, J.Z., 2008.  Sparse principal component analysis via regularized low rank matrix approximation.  J. Multivar. Anal. 99, 1015–1034.
  * Sijtsma (2009) Sijtsma, K., 2009.  On the use, the misuse, and the very limited usefulness of Cronbach’s alpha.  Psychometrika 74, 107–120.
  * Smilde et al. (2008) Smilde, A.K., Kiers, H.A.L., Bijlsma, S., Rubingh, C.M., van Erk, M.J., 2008.  Matrix correlations for high-dimensional data: The modified RV-coefficient.  Bioinformatics 25, 401–405.
  * Stewart and Sun (1990) Stewart, G.W., Sun, J., 1990.  Matrix Perturbation Theory.  Academic Press, Boston.
  * Székely and Rizzo (2013) Székely, G.J., Rizzo, M.L., 2013.  The distance correlation t-test of independence in high dimension.  J. Multivar. Anal. 117, 193–213.
  * Székely et al. (2007) Székely, G.J., Rizzo, M.L., Bakirov, N.K., 2007.  Measuring and testing dependence by correlation of distances.  Ann. Stat. 35, 2769–2794.
  * Thurstone and Thurstone (1941) Thurstone, L.L., Thurstone, T.G., 1941.  Factorial studies of intelligence.  University of Chicago Press, Chicago.
  * Tibshirani (1996) Tibshirani, R., 1996.  Regression shrinkage and selection via the lasso.  J. R. Stat. Soc. B 58, 267–288.
  * Vichi and Saporta (2009) Vichi, M., Saporta, G., 2009.  Clustering and disjoint principal component analysis.  Comput. Stat. Data Anal. 53, 3194–3208.
  * Vigneau et al. (2015) Vigneau, E., Chen, M., Qannari, M.E., 2015.  ClustVarLV: An R package for the clustering of variables around latent variables.  R J. 7, 134.
  * Vigneau and Qannari (2003) Vigneau, E., Qannari, E.M., 2003.  Clustering of variables around latent components.  Comm. Statist. Simulation Comput. 32, 1131–1150.
  * Witten et al. (2009) Witten, D.M., Tibshirani, R., Hastie, T.A., 2009.  A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis.  Biostatistics 10, 515–534.
  * Yang et al. (2014) Yang, D., Ma, Z., Buja, A., 2014.  A sparse singular value decomposition method for high-dimensional data.  J. Comput. Graph. Stat. 23, 923–942.
  * Yu et al. (2015) Yu, Y., Wang, T., Samworth, R.J., 2015.  A useful variant of the Davis—Kahan theorem for statisticians.  Biometrika 102, 315–323.
  * Zhu et al. (2020) Zhu, C., Zhang, X., Yao, S., Shao, X., 2020.  Distance-based and RKHS-based dependence metrics in high dimension.  Ann. Stat. 48, 3366–3394.
  * Zinbarg et al. (2005) Zinbarg, R.E., Revelle, W., Yovel, I., Li, W., 2005.  Cronbach’s $\alpha$, Revelle’s $\beta$, and McDonald’s $\omega_{h}$: Their relations with each other and two alternative conceptualizations of reliability.  Psychometrika 70, 123–133.
  * Zou et al. (2006) Zou, H., Hastie, T., Tibshirani, R., 2006.  Sparse principal component analysis.  J. Comput. Graph. Stat. 15, 265–286.