# Clusterpath Gaussian Graphical Modeling

Daniel J.W. Touw1, Andreas Alfons1, Patrick J.F. Groenen1, and Ines Wilms2
1Department of Econometrics, Erasmus University Rotterdam
2Department of Quantitative Economics, Maastricht University

###### Abstract

Graphical models serve as effective tools for visualizing conditional
dependencies between variables. However, as the number of variables grows,
interpretation becomes increasingly difficult, and estimation uncertainty
increases due to the large number of parameters relative to the number of
observations. To address these challenges, we introduce the Clusterpath
estimator of the Gaussian Graphical Model (CGGM) that encourages variable
clustering in the graphical model in a data-driven way. Through the use of a
clusterpath penalty, we group variables together, which in turn results in a
block-structured precision matrix whose block structure remains preserved in
the covariance matrix. We present a computationally efficient implementation
of the CGGM estimator by using a cyclic block coordinate descent algorithm. In
simulations, we show that CGGM not only matches, but oftentimes outperforms
other state-of-the-art methods for variable clustering in graphical models. We
also demonstrate CGGM’s practical advantages and versatility on a diverse
collection of empirical applications.

Keywords: clusterpath, graphical modeling, hierarchical clustering, precision
matrix, unsupervised learning

##  1 Introduction

Gaussian graphical models (GGM) have become an increasingly popular tool for
summarizing conditional dependencies among a set of $p$ variables. In
particular, a GGM is a conditional dependency network where one refers to the
$p$ variables in the graphical model as the nodes and the edges represent the
conditional dependency relations among each pair of variables. Estimating GGMs
is statistically challenging when the number of parameters ($p(p+1)/2$) is
large relative to the number of observations ($n$), leading to large
variability in estimation. Yet, such a setting arises across many and diverse
scientific areas, ranging from, amongst others, gene regulatory networks in
biology to contagion networks in finance and brain activation networks in
neuroscience. This has led to a flourishing area on regularized estimation
methods for GGMs (e.g., Meinshausen Bühlmann, 2006; Peng ., 2009; Yuan, 2010;
Cai ., 2011) that impose simplifying structures on the estimated parameters.
While much of the existing literature focuses on reducing the estimation
variability through edge sparsity, our approach branches into a different
direction by using node aggregation. Specifically, we integrate convex
clustering into the GGM framework to reduce estimation variability through the
pooling of estimates for similar variables.

Edge sparsity in GGMs through, for instance, $\ell_{1}$-regularization has
long formed the predominant choice to reduce dimensionality in GGMs (e.g.,
Yuan Lin, 2007; Banerjee ., 2008; Friedman ., 2008; Rothman ., 2008). The
focus here is to sparsely estimate the precision matrix (i.e., the inverse of
the covariance matrix) since conditional independencies between variable-pairs
can be directly obtained from its sparsity patterns, or equivalently the
absence of edges in the estimated GGM, offering an appealing interpretability
advantage. Several recent studies, however, point to a number of important
drawbacks when solely relying on edge sparsity as simplifying structure.
First, when the number of nodes grows large compared to a small sample size,
conditional dependencies among the variables may become too weak to detect
(Eisenach ., 2020). Such a situation is nowadays likely to occur more
frequently since (often noisy) data across many scientific areas are measured
at higher resolutions (Wilms Bien, 2022). Second, sparsely estimated GGMs with
many variables may still contain a large number of edges, thereby hindering
interpretability (Grechkin ., 2015). Third, many real world networks exhibit
more complex structures than mere edge sparsity implies. For instance, it is
not the most suitable simplicity structure for networks containing hub nodes,
that is, a group of nodes with similar and often strong dependency structures
(Hosseini Lee, 2016; Tarzanagh Michailidis, 2018; Tarzanagh ., 2021), or
bipartite structures where some variables depend linearly on others
(Heinävaara ., 2016).

When GGMs face such challenges, one is often not interested in estimating
conditional dependencies among the many observed variables but instead among a
smaller number of (unobserved) clusters (also called communities) of original
variables that share the same behavior. For instance, biologists estimate
large-scale gene regulatory networks and cluster genes into pathways to
unravel dependencies among them (e.g., Marlin Murphy, 2009; Grechkin ., 2015;
Kumar ., 2020), neuroscientists analyzing fMRI data routinely cluster voxels
into regions of interest to learn interaction networks of brain activation
(e.g., Pircalabelu Claeskens, 2020), financial analysts may cluster company
stocks into industry sectors to study how shocks can spread over the market by
contagion (e.g., Cardoso ., 2020). Cluster analysis is one of the most popular
unsupervised learning methods to discover the underlying group structures in
data. Variable clustering in the context of GGMs not only offers simple,
interpretable dependency networks but may also boost the dependency signals
(Eisenach ., 2020).

To estimate GGMs with clustered variables, a recent yet growing interest arose
in node-based dimensionality reduction. Initial proposals assume the clusters
to be known a priori and incorporate this information in the regularization
framework to encourage within-cluster dependencies over cross-cluster
dependencies (e.g., Duchi ., 2008; Schmidt ., 2009). Subsequent works learn
the node clustering either by decoupling the clustering task from the
estimation of the GGM (e.g., Ambroise ., 2009; Marlin Murphy, 2009; Tan .,
2015; Eisenach ., 2020; Brownlees ., 2022; Shi ., 2024) or by simultaneously
conducting both tasks (e.g., Hosseini Lee, 2016; Nie ., 2016; Tarzanagh
Michailidis, 2018; Kumar ., 2020; Pircalabelu Claeskens, 2020; Wilms Bien,
2022). Yet, to the best of our knowledge, the potential to leverage penalty
structures popular in the literature on convex clustering (Pelckmans ., 2005;
Hocking ., 2011; Lindsten ., 2011) to combine node clustering jointly with the
estimation of the GGM is left largely unexplored, a notable exception being
Yao Allen (2019) which we further contrast to our proposal in Section 2.4.

In this paper, we develop a novel regularizer, called the Clusterpath
estimator of the Gaussian Graphical Model, CGGM in short, to estimate GGMs
that are node-clustered where neither the number of clusters nor their
composition needs to be known in advance. To this end, we propose a novel
penalty on the distances between the variables in the precision matrix and
embed this in a convex optimization framework for which we offer a
computationally efficient cyclic block coordinate descent algorithm. The
resulting estimated precision matrix has a block structure where all variables
belonging to the same cluster share the same within- as well as cross-cluster
dependencies. A unique property of our approach for clustering the precision
matrix is that its inverse retains the same block structure, a property not
shared by other approaches.

We evaluate the performance of CGGM against state-of-the-art methods for
estimating node-clustered GGMs through a comprehensive simulation study
covering a wide range of graph structures. Our results indicate that CGGM
frequently surpasses the benchmark methods in both estimation accuracy and
clustering performance. CGGM performs competitively in scenarios with clear
sparsity patterns, even without a dedicated sparsity-inducing penalty. While
the main focus of our study is on estimating clustered precision matrices to
create graphical models, we also demonstrate that CGGM can be easily extended
to estimate clustered covariance matrices. In fact, when a block structure in
the covariance matrix (instead of the precision matrix) is the object of
interest, we demonstrate that directly estimating a clustered covariance
matrix has practical advantages over inverting a clustered precision matrix
estimate, even though the latter also yields a block structure in the
covariance matrix. Finally, we illustrate the effectiveness and versatility of
CGGM on three practical applications involving (i) stock market data from the
S&P 100, (ii) OECD well-being indicators, and (iii) survey data on
participants’ humor styles.

The remainder of the paper is organized as follows. Section 2 introduces the
clusterpath estimator for GGMs, details the cyclic block coordinate descent
algorithm, and relates our proposal to existing work. Next, we describe the
simulation study and present its results in Section 3. Section 4 explores the
use of CGGM for directly estimating clustered covariance matrices instead of
precision matrices. Section 5 demonstrates the application of CGGM to
empirical data, and Section 6 concludes the paper.

##  2 The Clusterpath Estimator for the Gaussian Graphical Model

We begin in Section 2.1 by discussing GGMs with clustered variables, followed
by the introduction to the clusterpath estimator in Section 2.2. Section 2.3
details the cyclic block coordinate descent algorithm used to compute our
estimator. Finally, Section 2.4 relates our proposal to existing work.

###  2.1 Clustered Gaussian Graphical Models

Let $\mathbf{X}$ be an $n\times p$ matrix of $n$ multivariate normal
observations each of dimension $p$, with sample mean $\mathbf{m}$ and sample
covariance matrix $\mathbf{S}$. Denoting the population covariance matrix by
$\mathbf{\Sigma}$, our target of estimation is the precision matrix
$\mathbf{\Theta}=\mathbf{\Sigma}^{-1}$. Under the assumption of a multivariate
normal distribution, the precision matrix can be equivalently expressed in a
graph where the nodes represent the $p$ variables and the edge weights are
given by the entries in $\mathbf{\Theta}$ which represent the conditional
dependencies among the variables.

Our goal is to estimate the precision matrix $\mathbf{\Theta}$ (and hence the
graph structure of the GGM) by encouraging clustering of the $p$ nodes in the
graph. The $K$ new cluster variables $\xi_{1},\ldots,\xi_{K}$ then represent
the average of the variables that belong to that cluster, with $K$ typically
much smaller than $p$ to achieve dimensionality reduction in the parameters
defining $\mathbf{\Theta}$. The edge weights among the clustered variables
represent the conditional dependencies among the $K$ new cluster variables.
Clustering of the variables in the graph corresponds to a block structure in
the rows and columns of $\mathbf{\Theta}$, see the so-called _G-block_ format
introduced in Bunea . (2020) and discussed by Wilms Bien (2022) in the context
of GGMs. The block-structure of the precision matrix can then be written as

$\displaystyle\mathbf{\Theta}=$
$\displaystyle\begin{bmatrix}(a_{11}-r_{11})\mathbf{I}&\mathbf{0}&\ldots&\mathbf{0}\\\
\mathbf{0}&(a_{22}-r_{22})\mathbf{I}&\ldots&\mathbf{0}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{0}&\mathbf{0}&\ldots&(a_{KK}-r_{KK})\mathbf{I}\\\
\end{bmatrix}+\begin{bmatrix}r_{11}\mathbf{11}^{\top}&r_{12}\mathbf{11}^{\top}&\ldots&r_{1K}\mathbf{11}^{\top}\\\
r_{21}\mathbf{11}^{\top}&r_{22}\mathbf{11}^{\top}&\ldots&r_{2K}\mathbf{11}^{\top}\\\
\vdots&\vdots&\ddots&\vdots\\\
r_{K1}\mathbf{11}^{\top}&r_{K2}\mathbf{11}^{\top}&\ldots&r_{KK}\mathbf{11}^{\top}\\\
\end{bmatrix},$ (1)

where $\mathbf{A}=\text{diag}(a_{11},\ldots,a_{KK})$ is a $K\times K$ diagonal
matrix, $\mathbf{R}=(r_{k\ell})_{1\leq k,\ell\leq K}$ is a $K\times K$
symmetric matrix, and $\mathbf{I}$ is the identity matrix of appropriate
dimension, similarly for $\mathbf{1}$ denoting a column-vector of ones and for
$\mathbf{0}$ denoting a matrix of zeros. The within-cluster conditional
variances $a_{kk}$ and covariances $r_{kk}$ are the same for all $p_{k}$
variables within cluster $k$. These $p_{k}$ variables in cluster $k$ also have
the same conditional covariance $r_{k\ell}$ with all $p_{\ell}$ variables
belonging to another cluster $\ell$. This clustered GGM is visualized for a
toy example with $p=8$ variables and $K=3$ clusters in Figure 1.

An important implication of our assumed block structure that is not shared by
other approaches is that the clustering structure of $\mathbf{\Theta}$ is
retained when taking the inverse (see, e.g., Gower Groenen, 1991; Archakov
Hansen, 2022). This property is lost in the approaches proposed by Yao Allen
(2019), Pircalabelu Claeskens (2020), and Wilms Bien (2022). Section 4
explores the estimation of block-structured covariance matrices using CGGM
through a numerical experiment.

Figure 1: Toy example of a graph representing the clustered precision matrix
with $K=3$ clusters constructed from $p=8$ original variables. Cluster
variable $\xi_{1}$ is the average of the $p_{1}=4$ variables
$X_{1},X_{2},X_{3}$, and $X_{4}$ having within cluster conditional covariance
$r_{11}$. Cluster variable $\xi_{2}$ is a singleton ($p_{2}=1$) and equal to
the original variable $X_{5}$. Cluster variable $\xi_{3}$ is the average of
the $p_{3}=3$ variables $X_{6},X_{7}$, and $X_{8}$ having within cluster
conditional covariance $r_{33}$. The three cluster variables have conditional
covariances $r_{12},r_{13},$ and $r_{23}$.

###  2.2 Clusterpath Estimator

To estimate precision matrices with a block structure corresponding to
clustered variables, we use a convex penalization method of the form

$\displaystyle\widehat{\mathbf{\Theta}}$
$\displaystyle=\operatorname*{argmin}_{\mathbf{\Theta}}\
L(\mathbf{\Theta})\qquad\text{s.t.}\
\mathbf{\Theta}=\mathbf{\Theta}^{\top},\mathbf{\Theta}\succ 0,$ (2)
$\displaystyle L(\mathbf{\Theta})$
$\displaystyle=-\log|\mathbf{\Theta}|+\operatorname*{tr}\mathbf{S\Theta}+\lambda\mathcal{P}(\mathbf{\Theta}),$
(3)

where $\log|\cdot|$ denotes the logarithm of the determinant,
$\operatorname*{tr}(\cdot)$ is the trace, $\cdot\succ 0$ denotes a positive
definite matrix, $\lambda\geq 0$ is a tuning parameter controlling the degree
of penalization, and $\mathcal{P}(\mathbf{\Theta})$ is the convex clusterpath
penalty penalizing differences between columns $\bm{\theta}_{j}$ and
$\bm{\theta}_{j^{\prime}}$ of the precision matrix $\mathbf{\Theta}$. We take

$\mathcal{P}(\mathbf{\Theta})=\sum_{j<j^{\prime}}w_{jj^{\prime}}d_{jj^{\prime}}(\mathbf{\Theta})\qquad\text{with}\qquad
d_{jj^{\prime}}(\mathbf{\Theta})=\sqrt{(\theta_{jj}-\theta_{j^{\prime}j^{\prime}})^{2}+\sum_{\begin{subarray}{c}m=1\\\
m\notin\\{j,j^{\prime}\\}\end{subarray}}^{P}(\theta_{jm}-\theta_{j^{\prime}m})^{2}},$

and weights $w_{jj^{\prime}}$ for every pair $1\leq j,j^{\prime}\leq p$. These
weights determine the fundamental attraction between two variables and are
specified beforehand based on, for example, domain knowledge or the pairwise
distances based on the sample precision matrix $\mathbf{S}^{-1}$. We further
discuss the choice of weights in Section 2.3.

The objective function is convex since each of the terms in (3) is convex in
$\mathbf{\Theta}$ and $\mathbf{\Theta}$ lies in the convex cone of positive
definite matrices. For $\lambda=0$, one obtains the inverse of the sample
covariance matrix as estimate, hence there is no block structure in the
precision matrix, nor a clustering of nodes in the GGM. If
$d_{jj^{\prime}}(\bm{\Theta})=0$ and $\lambda$ is positive, then two vectors
${\bm{\theta}}_{j}$ and ${\bm{\theta}}_{j^{\prime}}$ have exactly the same
elements while ignoring ${\theta}_{jj^{\prime}}$ and ${\theta}_{j^{\prime}j}$.
Note that the distance $d_{jj^{\prime}}(\mathbf{\Theta})$ disregards the
difference between $\theta_{jj^{\prime}}$ and $\theta_{j^{\prime}j}$ as these
are identical due to the symmetry of $\mathbf{\Theta}$. The penalty term
$d_{jj^{\prime}}(\bm{\Theta})$ can thus be interpreted as a group lasso
penalty (Yuan Lin, 2006), where each pair of variables $j<j^{\prime}$ forms a
group whose elements are the differences between the corresponding entries in
the columns $\bm{\theta}_{j}$ and $\bm{\theta}_{j^{\prime}}$ of the precision
matrix. If all respective differences are put to zero, the estimated entries
are identical, which in turn effectively blocks columns $j$ and $j^{\prime}$
in the precision matrix, or equivalently clusters nodes $j$ and $j^{\prime}$
in the GGM. As $\lambda$ increases, estimated GGMs are obtained in which more
and more variables are clustered, thereby resulting in a clusterpath from the
$p$ original nodes (no clustering) until one clustered node (full clustering)
in the GGM.

The notion of applying a distance-based penalty to pairs of objects stems from
convex clustering (Hocking ., 2011; Lindsten ., 2011; Pelckmans ., 2005). In
convex clustering, a copy of the data matrix is estimated while penalizing the
distances between the rows, thereby facilitating the clustering of
observations. Building upon this framework, convex biclustering (Chi ., 2017)
extends the conventional model by introducing a penalty on the distances
between columns, which allows for the clustering of variables as well. The
alternative provided by clustering variables using the precision matrix offers
insights into the underlying probabilistic dependencies within the data.

Next, we re-express the objective function $L(\mathbf{\Theta})$ in terms of
$\mathbf{A}$ and $\mathbf{R}$ (thereby assuming $1\leq K\leq p$ clusters)
instead of $\mathbf{\Theta}$ to permit more efficient expressions and updates
for solving optimization problem (2). To this end, we start by re-writing
equation (1) as

$\bm{\Theta}=\mathbf{UR}^{*}\mathbf{U}^{\top}+\begin{bmatrix}(a_{11}-r_{11})\mathbf{J}_{1}&\mathbf{0}&\cdots&\mathbf{0}\\\
\mathbf{0}&(a_{22}-r_{22})\mathbf{J}_{2}&\cdots&\mathbf{0}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{0}&\mathbf{0}&\cdots&(a_{KK}-r_{KK})\mathbf{J}_{K}\\\ \end{bmatrix},$
(4)

where $\mathbf{U}$ is the $p\times K$ cluster membership matrix with
$u_{jk}=1$ if variable $j$ belongs to cluster $k$ and zero otherwise,
$\mathbf{R}^{*}=\mathbf{R}+\text{diag}((a_{11}-r_{11})/p_{1},\ldots,(a_{KK}-r_{KK})/p_{K})$,
and $\mathbf{J}_{k}$ is the $p_{k}\times p_{k}$ centering matrix. Now,
equation (4) presents an orthogonal decomposition where the precision matrix
is split in two parts, the first part is a linear space that depends on
$\mathbf{A}$ and all elements of $\mathbf{R}$, the second part is an
orthogonal linear space that only depends on the diagonal blocks of the
precision matrix, namely $\mathbf{A}$ and the diagonal elements of
$\mathbf{R}$.

Using decomposition (4), we can equivalently express the objective function in
(3) as

$\displaystyle L(\mathbf{\Theta})$ $\displaystyle=L(\mathbf{A,R})$ (5)
$\displaystyle=-\log|\mathbf{P}^{1/2}\mathbf{R}^{*}\mathbf{P}^{1/2}|-\sum_{k=1}^{K}(p_{k}-1)\log(a_{kk}-r_{kk})+\operatorname*{tr}\mathbf{SURU}^{\top}+\sum_{\ell=1}^{K}(a_{\ell\ell}-r_{\ell\ell})\operatorname*{tr}\mathbf{S}_{\ell}+\lambda\mathcal{P}(\mathbf{\mathbf{A,R}}),$

where the diagonal matrix
$\mathbf{P}=\mathbf{U}^{\top}\mathbf{U}=\text{diag}(p_{1},\ldots,p_{K})$
contains the number of variables $p_{k}$ per cluster on the diagonal,
$\mathbf{P}^{1/2}$ simply contains their square roots, $\mathbf{S}_{\ell}$ is
the sample covariance matrix computed from the $p_{\ell}$ variables in cluster
$\ell$, and the clusterpath penalty is reformulated as

$\mathcal{P}(\mathbf{\mathbf{A,R}})=\sum_{k<\ell}\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{\ell}d_{k\ell}(\mathbf{A},\mathbf{R}),$

where $\mathbf{W}$ is the symmetric $p\times p$ weight matrix containing the
individual weights $w_{jj^{\prime}}$, $\mathbf{u}_{k}$ represents the $k$th
column of the membership matrix $\mathbf{U}$ and

$d_{k\ell}^{2}(\mathbf{A},\mathbf{R})=(a_{kk}-a_{\ell\ell})^{2}+(p_{k}-1)(r_{kk}-r_{k\ell})^{2}+(p_{\ell}-1)(r_{\ell\ell}-r_{k\ell})^{2}+\sum_{\begin{subarray}{c}m=1\\\
m\notin\\{k,\ell\\}\end{subarray}}^{K}p_{m}(r_{km}-r_{\ell m})^{2}.$

Full details of this derivation are provided in Appendix A.

We are now ready to optimize the objective function $L(\mathbf{A,R})$, which
can be done efficiently—especially as the number of clusters $K$ decreases.

###  2.3 Cyclic Block Coordinate Descent Algorithm

We propose a cyclic block coordinate descent algorithm to minimize objective
function (5) under the constraints from (2), where the blocks are formed by
the clusters of variables. To permit cyclic block updates, the objective
$L(\mathbf{A},\mathbf{R})$ needs to be separated into those parts that depend
on cluster $k$, namely $a_{kk}$, $\mathbf{r}_{k}$, and $r_{kk}$, and those
that do not. The $K-1$ vector $\mathbf{r}_{k}$ contains the elements of the
$k$th column of $\mathbf{R}$ excluding $r_{kk}$. We distinguish
$\mathbf{r}_{k}$ and $r_{kk}$ because the diagonal elements of $\mathbf{R}$
interact differently with the objective function compared to other elements
related to the $k$th cluster. Additionally, we partition the objective
function into four distinct components, as given by

$L(a_{kk},\mathbf{r}_{k},r_{kk})=L_{\text{det}}(a_{kk},\mathbf{r}_{k},r_{kk})+L_{\text{cov}}(a_{kk},\mathbf{r}_{k},r_{kk})+L_{\text{cpath}}(a_{kk},\mathbf{r}_{k},r_{kk})+C,$
(6)

where $L_{\text{det}}$ ($L_{\text{cov}}$) represents the log-determinant
(trace) part of the objective function that depends on cluster $k$,
$L_{\text{cpath}}$ is the clusterpath penalty that depends on cluster $k$, and
finally $C$ is a constant collecting parts that are independent of cluster
$k$. Details of these expressions are given in Appendix B.

The proposed algorithm now involves two key steps. For each cluster $k$, it
first assesses whether there is an eligible fusion candidate. In the instances
where no suitable candidate is identified, the algorithm proceeds to the
second step, wherein the parameters associated with cluster $k$ are updated
using Newton’s method. The reason we choose a cyclic block coordinate descent
algorithm lies with the use of Newton’s method since the computation of a
Newton descent direction for all parameters in $\mathbf{A}$ and $\mathbf{R}$
simultaneously would result in a computationally intractable algorithm even
for moderately sized data sets.

In Section 2.3.1, we provide details on the two key steps of the optimization
algorithm. We further discuss two important components related to the
clusterpath penalty, namely the weight matrix (Section 2.3.2) and the tuning
parameter (Section 2.3.3). Finally, in Section 2.3.4, we discuss a refitting
step to reduce bias in the CGGM estimate.

####  2.3.1 Outline of the Cyclic Block Coordinate Descent Algorithm

We outline the two-step algorithm minimizing the objective function in (6)
iteratively for cluster $1\leq k\leq K$ for a fixed value for $\lambda$.

Cluster fusions. The first step is to verify whether cluster $k$ is close
enough to another cluster to warrant their fusion. Cluster $k$ is fused with
cluster $\ell$ if $d_{k\ell}(\mathbf{A},\mathbf{R})$ is smaller than some
user-defined threshold $\varepsilon_{f}$ representing the minimum required
similarity for fusion. To account for the scale of the data, we take
$\varepsilon_{f}=\tau\operatorname*{median}_{j,j^{\prime}}(d_{jj^{\prime}}(\mathbf{\Theta})),$
where we typically set $\tau=10^{-3}$ in our experiments. In case of multiple
candidates for which this applies, cluster $\ell$ is chosen as
$\ell=\operatorname*{argmin}_{\ell^{\prime}}d_{k\ell^{\prime}}(\mathbf{A},\mathbf{R})$.
If a fusion is performed, the algorithm does not proceed to the next step as
cluster $k$ has ceased to exist, instead it advances to the next cluster. If
no eligible fusion candidate is found in the first step, we proceed to the
second step.

Descent direction. The second step is to compute a descent direction. We opt
for a Newton descent direction as it offers efficient convergence properties.
Let $\nabla L(a_{kk},\mathbf{r}_{k},r_{kk})$ denote the gradient and
$\nabla^{2}L(a_{kk},\mathbf{r}_{k},r_{kk})$ the Hessian of
$L(a_{kk},\mathbf{r}_{k},r_{kk})$, see Appendix B for their derivations. The
descent direction is then computed as
$\bm{\delta}_{k}=-\nabla^{2}L(a_{kk},\mathbf{r}_{k},r_{kk})^{-1}\nabla
L(a_{kk},\mathbf{r}_{k},r_{kk})$.

Finally, we determine the optimal step size $s^{*}$ given the descent
direction using a line search. The reason for augmenting the Newton descent
direction with a line search is twofold. First, the minimization of
$L(a_{kk},\mathbf{r}_{k},r_{kk})$ is constrained by the restriction that
$\mathbf{\Theta}$ should be positive definite. Hence, a step size should be
chosen that ensures the update to adhere this restriction. Assuming that
$\mathbf{\Theta}$ is positive definite for the current values of
$[a_{kk},\mathbf{r}_{k},r_{kk}]$ and the number of clusters is larger than
one, positive definiteness of $\bm{\Theta}$ is preserved by a step size $s$
that satisfies

$\displaystyle(a_{kk}+s\delta_{a_{kk}})+(p_{k}-1)(r_{kk}+s\delta_{r_{kk}})-p_{k}(\mathbf{r}_{k}+s\bm{\delta}_{\mathbf{r}_{k}})^{\top}\mathbf{V}(\mathbf{r}_{k}+s\bm{\delta}_{\mathbf{r}_{k}})$
$\displaystyle>0$ $\displaystyle
a_{kk}+s\delta_{a_{kk}}-r_{kk}-s\delta_{r_{kk}}$ $\displaystyle>0,$ (7)

where $\delta_{a_{kk}}$, $\bm{\delta}_{\mathbf{r}_{k}}$, and $\delta_{r_{kk}}$
represent the descent directions for $a_{kk}$, $\mathbf{r}_{k}$, and $r_{kk}$.
In case there is a single cluster left, the first inequality reduces to

$(a_{kk}+s\delta_{a_{kk}})+(p_{k}-1)(r_{kk}+s\delta_{r_{kk}})>0.$ (8)

A second reason to opt for a line search is the fact that the Hessian may not
be well-behaved due to the presence of the $\ell_{2}$-norm in the clusterpath
penalty. If, for the current values of $\mathbf{A}$ and $\mathbf{R}$, at least
one of the distances $d_{k\ell}(\mathbf{A},\mathbf{R})$ is close to zero,
$L(a_{kk},\mathbf{r}_{k},r_{kk})$ is not locally smooth. As a result, the
Hessian is not well-behaved and the favorable convergence properties of
Newton’s method do not hold. A line search for the optimal step size
guarantees that the update $s^{*}\bm{\delta}_{k}$ does not increase the
objective function.

The pseudo-code for the algorithm is presented in Algorithm 1. To initialize
$\mathbf{A}$ and $\mathbf{R}$, one can use the inverse of the sample
covariance matrix if it is available; then
$\mathbf{A}^{(1)}=\text{diag}(\mathbf{S}^{-1})$ and
$\mathbf{R}^{(1)}=\mathbf{S}^{-1}$. In other cases, one can start from a
regularized estimator such as the graphical lasso (Yuan Lin, 2007; Banerjee .,
2008; Rothman ., 2008; Friedman ., 2008).

Algorithm 1 Pseudo-code for the algorithm that minimizes the CGGM objective
function for fixed $\lambda$.

Input Initial estimates for $\mathbf{A}^{(1)}$ and $\mathbf{R}^{(1)}$, sample
covariance matrix $\mathbf{S}$, weight matrix $\mathbf{W}$, tuning parameter
$\lambda$, fusion threshold $\varepsilon_{f}$, convergence threshold
$\varepsilon_{c}$, maximum number of iterations $t_{\text{max}}$
Output $\widehat{\mathbf{A}}$ and $\widehat{\mathbf{R}}$ that minimize
$L(\mathbf{A},\mathbf{R})$ with $\widehat{\mathbf{U}}$ that describes cluster
membership

1:$\mathbf{U}^{(1)}\leftarrow\mathbf{I}$

2:$L^{(1)}\leftarrow L(\mathbf{A}^{(1)},\mathbf{R}^{(1)})$

3:$L^{(0)}\leftarrow(1+2\varepsilon_{c})L^{(1)}$

4:$t\leftarrow 1$

5:while $L^{(t-1)}/L^{(t)}-1>\varepsilon_{c}$ and $t\leq t_{\text{max}}$ do

6: $t\leftarrow t+1$

7: $\mathbf{A}^{(t)}\leftarrow\mathbf{A}^{(t-1)}$

8: $\mathbf{R}^{(t)}\leftarrow\mathbf{R}^{(t-1)}$

9: $\mathbf{U}^{(t)}\leftarrow\mathbf{U}^{(t-1)}$

10: for $k=1,\ldots,K$ do

11:
$\ell\leftarrow\operatorname*{argmin}_{\ell^{\prime}}d_{k\ell^{\prime}}(\mathbf{A}^{(t)},\mathbf{R}^{(t)})$

12: if $d_{k\ell}(\mathbf{A}^{(t)},\mathbf{R}^{(t)})\leq\varepsilon_{f}$ then

13: Fuse clusters $k$ and $\ell$ by modifying $\mathbf{A}^{(t)}$,
$\mathbf{R}^{(t)}$, $\mathbf{U}^{(t)}$

14: else

15:
$\bm{\delta}_{k}\leftarrow-\nabla^{2}L\bigl{(}a^{(t)}_{kk},\mathbf{r}^{(t)}_{k},r^{(t)}_{kk}\bigr{)}^{-1}\nabla
L\bigl{(}a^{(t)}_{kk},\mathbf{r}^{(t)}_{k},r^{(t)}_{kk}\bigr{)}$

16: Compute maximum step size $s_{\text{max}}$ using equations (2.3.1) and (8)

17: Select optimal step size $s^{*}\in[0,s_{\text{max}})$

18:
$\bigl{[}a^{(t)}_{kk},\mathbf{r}^{(t)}_{k},r^{(t)}_{kk}\bigr{]}\leftarrow\bigl{[}a^{(t)}_{kk},\mathbf{r}^{(t)}_{k},r^{(t)}_{kk}\bigr{]}+s^{*}\bm{\delta}_{k}$

19:$\widehat{\mathbf{A}}\leftarrow\mathbf{A}^{(t)}$

20:$\widehat{\mathbf{R}}\leftarrow\mathbf{R}^{(t)}$

21:$\widehat{\mathbf{U}}\leftarrow\mathbf{U}^{(t)}$

####  2.3.2 Weight Matrix

An important element of the clusterpath estimator is the weight matrix
$\mathbf{W}$ that contains information about the preferences of clustering
variables $j$ and $j^{\prime}$ through the weight $w_{jj^{\prime}}$. A large
value for $w_{jj^{\prime}}$ incentivizes clustering of the corresponding
variables, but it does not guarantee this. In the convex clustering
literature, it is standard practice to set the weights based on the distances
between the original data points (see, e.g., Chi ., 2017; Hocking ., 2011; Sun
., 2021). Furthermore, clustering performance increases when using $k$-nearest
neighbors to assign zero weights to pairs of variables that are beyond each
other’s $k$ nearest neighbors (Chi Lange, 2015). Following these standard
practices, we take

$w_{jj^{\prime}}=w_{j^{\prime}j}=\begin{cases}\exp(-\phi
d^{2}_{jj^{\prime}}(\mathbf{S}^{-1}))&\text{if
$(j,j^{\prime})\in\mathcal{E}$}\\\ 0&\text{otherwise},\end{cases}$

where $\phi$ is a tuning parameter and $\mathcal{E}$ is the set of variable-
pairs $(j,j^{\prime})$ that should be assigned a nonzero weight.

A potential issue with relying solely on $k$-nearest neighbors for the
construction of $\mathcal{E}$ is the occurrence of disconnected subgroups. To
alleviate this issue, we add variable-pairs to $\mathcal{E}$ using the
procedure described in Kruskal (1956) until all disconnected subgroups are
eliminated. This approach effectively integrates the $k$-nearest neighbors
structure with a minimum spanning tree to obtain a weight matrix corresponding
to a connected graph (Touw ., 2022).

####  2.3.3 Tuning Parameters

The tuning parameter $\lambda$ forms another key ingredient of the clusterpath
estimator as it controls the degree of clustering. First, note that in the
current form of objective function (3), the tuning parameter $\lambda$ is
sensitive to the number of variables $p$. To reduce this sensitivity, the
different terms of the objective function may be scaled. We suggest to scale
the terms pertaining to the log likelihood by $p^{-1}$ and the penalty by
$\kappa=((p-1)^{1/2}\sum_{j<j^{\prime}}w_{jj^{\prime}})^{-1}$ to reduce this
dependence on $p$. These scaling factors can then easily be absorbed into a
rescaled tuning parameter $\gamma=p\kappa\lambda$, which then replaces
$\lambda$ in the objective function.

To set a sequence of tuning parameters, we implement an automated procedure.
First, $\lambda$ is initialized at 0.5 and undergoes iterative increments of
50% until the minimum number of clusters, as determined by the weight matrix,
is attained. This yields an increasing series of values $\\{\lambda_{q}:1\leq
q\leq Q\\}$ alongside their corresponding solutions
$\\{\widehat{\mathbf{\Theta}}_{q}:1\leq q\leq Q\\}$. To ensure smoothness in
the trajectory, additional values of $\lambda$ are inserted whenever the
difference between consecutive solutions, as measured by
$\|\widehat{\mathbf{\Theta}}_{q-1}-\widehat{\mathbf{\Theta}}_{q}\|/\|\widehat{\mathbf{\Theta}}_{q-1}\|$,
exceeds 0.01. Consequently, a continuum of solutions for $\mathbf{\Theta}$ is
obtained, transitioning smoothly from minimal to maximal regularization, with
a hierarchical clustering structure. Throughout this iterative process, the
algorithm leverages existing solutions for the precision matrix as warm starts
to speed up finding solutions for new values of $\lambda$.

Apart from the tuning parameter $\lambda$ that controls the degree of
clustering, one also needs to select the tuning parameters $k$ and $\phi$ that
are used to compute the weight matrix. Typical values of interest for $k$
range from 1 to $p^{1/2}$. The choice for $\phi$ depends heavily on the scale
of the data. The smallest value of $\phi$ that is of interest ensures that
most nonzero weights lie between 0.9 and 1.0, whereas the largest value
results in a more uniform distribution of the nonzero weights.

To select the tuning parameters $\lambda$, $k$, and $\phi$ we adopt a cross-
validation procedure. We form a grid of tuning parameters ($\lambda$, $k$,
$\phi$) and find the combination of the tuning parameters that minimizes the
cross-validated likelihood-based score given by

$\frac{1}{G}\sum_{g=1}^{G}-\log|\widehat{\mathbf{\Theta}}_{-\mathcal{F}_{g}}|+\operatorname*{tr}\mathbf{S}_{\mathcal{F}_{g}}\widehat{\mathbf{\Theta}}_{-\mathcal{F}_{g}},$
(9)

where $\widehat{\mathbf{\Theta}}_{-\mathcal{F}_{g}}$ represents the precision
matrix estimated on the sample excluding the observations in the $g$th fold
($1\leq g\leq G$, one can take for instance $G=5$), and
$\mathbf{S}_{\mathcal{F}_{g}}$ denotes the sample covariance matrix computed
solely on the samples within the $g$th fold. The density of the grid for $k$
and $\phi$ depends on $p$; low-dimensional data allows for shorter computation
times and denser grids. In this paper, we use grids of five to eight evenly
spaced values for both tuning parameters.

####  2.3.4 Refitted CGGM

While minimization of the objective function in (5) yields a clustered
$\widehat{\mathbf{\Theta}}$, the CGGM estimate may be biased due to the
shrinkage of the distances between different clusters. Following Wilms Bien
(2022), we therefore also consider a refitted version of CGGM. First, CGGM is
used to obtain a clustering in $\widehat{\mathbf{\Theta}}$, as defined by the
cluster membership estimate $\widehat{\mathbf{U}}$. We then re-estimate the
precision matrix $\mathbf{\Theta}$ constrained by the obtained clustering by
maximizing the likelihood subject to the clustering constraint. This is
equivalent to minimizing the objective in (5) with initial clustered estimates
for $\mathbf{A}$ and $\mathbf{R}$ and $\lambda=0$. By re-estimating the
precision matrix, we aim to capitalize on the clustering capabilities of CGGM
while refining parameter estimates for improved accuracy.

###  2.4 Connections to Related Work

Clustering in the context of GGMs arises either at the level of the
observations111A separate literature stream thus considers GGMs where the
observations come from different subpopulations (i.e. clusters) that are
either a priori known (e.g., Danaher ., 2014; Mohan ., 2014) or learned from
the data by unifying clustering with GGM estimation into one optimization
framework (e.g., Gao ., 2016; Hao ., 2018). or at the level of the variables
(i.e., the nodes). Here, we focus on GGMs with clustered nodes. Contributions
to this strand of the literature can be divided into three main streams.

A first stream assumes one can a priori cluster the variables into known
groups which share similar dependencies (e.g., Duchi ., 2008; Schmidt ., 2009;
Grechkin ., 2015; Millington Niranjan, 2019). Domain-knowledge may, however,
not always be available to impose a grouping, thereby still calling for
unsupervised clustering procedures. A second stream of papers therefore uses
two-step sequential procedures that first learn how to cluster the variables
and then turn to the estimation of the GGM (e.g., Ambroise ., 2009; Marlin
Murphy, 2009; Tan ., 2015; Eisenach ., 2020; Brownlees ., 2022). A two-step
procedure that decouples the clustering step from the GGM estimation may,
however, lead to suboptimal results, thereby giving rise to a third stream
that considers both tasks jointly.

In this third stream, Nie . (2016), Kumar . (2020), Cardoso . (2020), and
Gheche Frossard (2020), amongst others, study $k$-connected component graphs
permitting variable clustering in $k$ groups due to the spectral properties of
the Laplacian matrix. Besides, penalized likelihood procedures form a popular
alternative to encourage various forms of node-clustering (e.g., Hosseini Lee,
2016; Tarzanagh Michailidis, 2018). To this end, our proposal is most closely
related to Yao Allen (2019), Pircalabelu Claeskens (2020), and Wilms Bien
(2022).

The first to explore the use of a clusterpath penalty in precision matrix
estimation are Yao Allen (2019). However, a key distinction with our proposal
lies in the distance metric used in the penalty: their approach does not
account for the diagonal elements of $\mathbf{\Theta}$ whereas ours does. This
is an important distinction in the context of GGMs since variables clustered
by CGGM share identical values on the diagonal of the precision matrix. As a
consequence, unlike for the proposal of Yao Allen (2019), the clustering in
$\bm{\widehat{\Theta}}$ is retained in its inverse (see, e.g., Gower Groenen,
1991; Archakov Hansen, 2022). Pircalabelu Claeskens (2020) use a penalized
likelihood procedure with a group penalty to identify communities of nodes,
their procedure—unlike ours—does not result in a precision matrix with block
structure. Wilms Bien (2022), finally, use a penalized likelihood procedure
with a tree-based penalty to cluster nodes, hence their choice of penalty
differs from ours as it requires side-information on the similarity of nodes
to guide their clustering. We include the proposals of Pircalabelu Claeskens
(2020) and Wilms Bien (2022) as benchmarks in our simulation study since their
software implementations are publicly available.

##  3 Simulations

To investigate the behavior of the proposed method CGGM, we perform an
extensive simulation study. We discuss the data generating processes in
Section 3.1, the benchmark methods and evaluation criteria in Section 3.2, and
the results in Section 3.3.

###  3.1 Simulation Designs

For most of the simulations, we follow the designs of Wilms Bien (2022)
together with some additional variations. In all settings, we sample data from
a multivariate normal distribution with mean zero and covariance matrix
$\mathbf{\Sigma}=\mathbf{\Theta}^{-1}$, and we repeat the process 100 times.

Baseline simulation designs. We simulate $n=120$ observations on $p=15$
variables using four different structures for the precision matrix
$\mathbf{\Theta}$ (see Figure 2, top row):

  * •

_Random_ : The clusters are of equal size with $p_{1}=p_{2}=p_{3}=5$; one pair
of clusters is selected at random to be connected via an edge.

  * •

_Chain_ : The clusters are of equal size with $p_{1}=p_{2}=p_{3}=5$; adjacent
clusters are connected via an edge.

  * •

_Unbalanced_ : Same as the chain design, but the clusters are of unequal size
with $p_{1}=3$, $p_{2}=5$, and $p_{3}=7$.

  * •

_Unstructured_ : Each variable forms its own cluster; edges between variables
are drawn with probability $\pi=0.1$.

The first three designs thus employ a block structure with $K=3$ clusters,
while the fourth does not exhibit any variable aggregation structure. In all
four designs, the diagonal elements of $\mathbf{\Theta}$ are set to 1, the
elements within a cluster of variables to 0.5, and the non-zero elements
between clusters to 0.25.

Figure 2: Structure in the precision matrix $\mathbf{\Theta}$ in the four
baseline simulation designs taken from Wilms Bien (2022) (top), as well as in
the designs with clustering structure on the diagonal and with blockdiagonal
structure, respectively, using balanced and unbalanced clusters sizes
(bottom). The color shade indicates the magnitude of the elements. Diagonal
elements are on a separate color scale than the off-diagonal ones to highlight
their differing roles in CGGM in comparison to the benchmark methods.

Increasing the number of variables. We focus on the chain design as in Wilms
Bien (2022) and vary the number of variables $p\in\\{15,30,60\\}$. We set the
corresponding number of observations $n\in\\{120,240,480\\}$ such that the
ratio $n/p$ is constant. The number of clusters is kept fixed at $K=3$.

Increasing the number of clusters. We again focus on the chain design and vary
the number of clusters $K\in\\{3,5,6,10\\}$ while keeping the number of
observations and variables fixed at $n=240$ and $p=30$, respectively, as in
Wilms Bien (2022).

Approximate block structure. We take the four baseline designs and modify the
structure in $\mathbf{\Theta}$ by uniformly drawing elements within a cluster
of variables from the interval $[0.4,0.6]$ and non-zero elements between
clusters from the interval $[0.2,0.3]$. The assumptions of CGGM are not
(fully) met in the resulting designs, as the clusters do not correspond to
blocks of equal values in $\mathbf{\Theta}$ but blocks of approximately
similar values.

Clustering structure on diagonal. Using the same dimensions and number of
clusters as in the baseline designs, all off-diagonal elements in
$\mathbf{\Theta}$ are set to 0.5. The diagonal elements are set to 1 for the
variables in the first cluster, to 2 for the variables in the second cluster,
and to 3 for the variables in the third cluster. We investigate both a
_balanced_ design with $p_{1}=p_{2}=p_{3}=5$, and an _unbalanced_ design with
$p_{1}=3$, $p_{2}=5$, and $p_{3}=7$ (see left two panels in the bottom row of
Figure 2).

Blockdiagonal structure. In these designs, the precision matrix
$\mathbf{\Theta}$ does not exhibit a variable clustering structure, but a
blockdiagonal structure inspired by the simulation design of Pircalabelu
Claeskens (2020). We use the same dimensions as in the baseline design and
$K=3$ blocks on the diagonal. Outside of the diagonal blocks, edges between
variables are drawn with probability $\pi=0.1$. The diagonal values of
$\mathbf{\Theta}$ are set to 1, the off-diagonal elements in the diagonal
blocks are set to 0.5, and non-zero elements outside the diagonal blocks are
set to 0.25. If such a randomly drawn precision matrix is not positive semi-
definite, we repeat this process until a positive semi-definite precision
matrix is obtained. We consider both a _balanced_ design with
$p_{1}=p_{2}=p_{3}=5$, and an _unbalanced_ design with $p_{1}=3$, $p_{2}=5$,
and $p_{3}=7$ (see the right two panels in the bottom row of Figure 2).

###  3.2 Methods and Evaluation Criteria

We apply CGGM with and without the parameter re-estimation step described in
Section 2.3.4, which we refer to as _CGGM-raw_ and _CGGM-refit_. We thereby
set the convergence tolerance $\varepsilon_{c}=10^{-7}$ and the maximum number
of iterations to $t_{\text{max}}=5000$. For comparison, we include the tree-
aggregated graphical lasso (_TAGL_) of Wilms Bien (2022), the community-based
group graphical lasso (_ComGGL_) of Pircalabelu Claeskens (2020), and the
inverse of the sample covariance matrix ($\mathbf{S}^{-1}$). TAGL performs
node aggregation based on side-information in the form of a tree-based
variable hierarchy, which needs to be specified a priori. We generate an
_ideal_ and a _realistic_ tree hierarchy as described in Wilms Bien (2022)
(see Figure 5 in their paper for an illustration). As both trees contain the
true variable clustering, we also generate a _misspecified_ tree in the same
manner as the realistic tree, except that with probability 0.1, each variable
is incorrectly assigned to the subsequent cluster (or the first cluster if the
variable in fact belongs to the last cluster). ComGGL does not perform node
aggregation but community detection in the form of a blockdiagonal structure.
As opposed to CGGM, both TAGL and ComGGL also apply lasso-regularization for
edge sparsity.

Concerning tuning parameter selection, we apply 3-fold cross-validation. For
the weight matrix in CGGM, we use candidate values
$\phi\in\\{1,1.5,2,2.5,3\\}$ and numbers of neighbors $k\in\\{1,2,3,4,5\\}$.
For each combination of $\phi$ and $k$, we select candidate values for the
regularization parameter $\lambda$ as described in Section 2.3.3. For TAGL, we
first perform a binary search for the smallest value of the aggregation
parameter that aggregates all variables into one cluster (for each of the tree
hierarchies) while keeping the sparsity parameter at 0, and we determine the
smallest value of the sparsity parameter that removes all edges while keeping
the aggregation parameter at 0 (corresponding to the graphical lasso of
Friedman ., 2008). The candidate values for the aggregation and sparsity
parameters are then obtained as fractions of those aforementioned values,
where the fractions vary from 0 to 1 in 10 steps, with each step being twice
as large as the previous step. For ComGGL, we similarly set the candidate
values of the grouping and sparsity parameters as fractions of certain maximum
values. While we fix the maximum value of the grouping parameter at 1, we
again set the maximum value of the sparsity parameter as the smallest value
that removes all edges in the graphical lasso. As in Pircalabelu Claeskens
(2020), we set the balancing parameter to 1.

We evaluate the methods in terms of their estimation accuracy and aggregation
performance. Regarding estimation accuracy, we compute the Frobenius norm
$\|\mathbf{\Theta}-\mathbf{\hat{\Theta}}\|_{F}$. For aggregation performance,
we report the estimated number of clusters $\hat{K}$ as well as the adjusted
Rand index (ARI) (Hubert Arabie, 1985) of the obtained clustering of
variables. Note that $\hat{K}$ directly follows from the obtained clustering
for the optimal values of the tuning parameters from cross-validation.

###  3.3 Results

In Figures 3–5, we report the evaluation criteria for different simulation
designs. As CGGM and TAGL outperform ComGGL and the inverse of the sample
covariance matrix ($\mathbf{S}^{-1}$) in almost all cases, we focus on the
former two methods in the following discussion of the results.

Baseline simulation designs. In terms of estimation accuracy and aggregation,
CGGM-raw and CGGM-refit perform similar to TAGL-ideal, outperform TAGL-
realistic in some designs, and outperform and TAGL-misspecified in all designs
(see Figure 3). The refitting step of CGGM is advantageous in the three block
designs but disadvantageous in the unstructured design. For the latter, the
refitting step often aggregates all variables into one cluster, likely due to
the high degree of sparsity in the unstructured design. This pooling of
estimates for the large number of zero elements reduces noise, potentially
improving the out-of-sample cross-validation score in equation (9), despite
introducing some bias in the relatively small number of nonzero elements.
Conversely, CGGM-raw allows for shrinkage between elements without
aggregation, thus improving estimation performance for designs in which many
off-diagonal elements share the same value without belonging to the same
cluster. While TAGL-misspecified often selects the correct number of clusters,
the ARI reveals that it struggles to accurately classify the variables within
those clusters due to the misspecified tree. This pattern persists across most
of the simulation designs.

Figure 3: Results for the baseline simulation designs (in separate columns).
The top row contains boxplots of the Frobenius norm together with black
diamonds representing the average. In the middle and bottom rows, diamonds
display the average of the estimated number of clusters and the average ARI,
respectively. In addition, the size of the grey dots represent the frequency
of different values across the replications. Reference lines are drawn for the
true number of clusters and the maximum possible ARI value. Aggregation
performance is omitted for $\mathbf{S}^{-1}$ as it does not perform any
aggregation. In the unstructured design, a misspecified tree for TAGL does not
exist since any tree hierarchy contains the true clustering (each variable
being its own cluster). Figure 4: Results for the chain simulation design with
an increasing number of clusters (in separate columns). See Figure 3 for
explanatory notes.

Increasing the number of variables or clusters. Given that CGGM-refit clearly
outperforms CGGM-raw in the baseline chain design, we omit the results for the
latter in the variations with increasing number of variables or clusters.
First, the results remain stable when increasing the number of variables and
are similar to those in the baseline setting (see Figure 10 in Appendix C),
likely due to maintaining a fixed ratio of $n$ to $p$. CGGM-refit performs
comparably to TAGL-ideal and TAGL-realistic, and clearly surpasses TAGL-
misspecified. Second, as the number of clusters increases, the aggregation
performance of CGGM-refit remains stable, although the estimation error in
terms of the Frobenius norm increases (see Figure 4). In the chain design, the
sparsity of the precision matrix increases in the number of clusters. Hence, a
method that enforces sparsity—such as TAGL—is expected to perform well,
particularly when up to 56% of the elements are zero (as in the design with
$K=10$). While the median estimation error of CGGM-refit increases in
comparison to TAGL-ideal and TAGL-realistic, it remains considerably lower
than that of TAGL-misspecified. To investigate estimation performance further,
we also considered the negative log-likelihood of the estimated precision
matrix, which remains similar for CGGM-refit and TAGL-ideal as the number of
clusters increases (results not shown). Finally, CGGM-refit in some settings
exhibits lower variability in the estimation error compared to TAGL-realistic.

Approximate block structure. The results for estimation accuracy are highly
similar to those of the baseline simulation designs (see Figure 11 in Appendix
C). In terms of aggregation performance, CGGM-raw struggles to find the block
structure when it is no longer exact, but the refitting step overcomes this
issue. TAGL-realistic also performs worse in terms of aggregation, though the
effect is less pronounced than with CGGM-raw. Overall, aggregation performance
of CGGM-refit and TAGL-ideal remains very similar to the baseline designs. It
is important to note that, technically, the identified block structure is
incorrect at the population level. However, the noise reduction in finite
samples is still beneficial.

Figure 5: Results for the simulation designs with the clustering structure on
the diagonal (left two columns) and with a blockdiagonal structure (right two
columns). See Figure 3 for explanatory notes.

Clustering structure on diagonal. CGGM-refit clearly outperforms the
alternative methods in terms of aggregation performance (see left two columns
in Figure 5). All variants of TAGL tend to aggregate all variables into one
cluster, as they exclude the elements on the diagonal from the aggregation
penalty. As a result, the variability in the estimation error of TAGL is
considerably larger compared to CGGM. Despite the better aggregation
performance, the estimation error is slightly higher for CGGM-refit than for
CGGM-raw. A likely explanation of this phenomenon is that CGGM-refit no longer
applies shrinkage to the difference between elements from different clusters.
In this design, this behavior does not align with the structure of the
underlying precision matrix, in which all off-diagonal elements share the same
value. Consequently, by shrinking the off-diagonal elements from different
clusters towards each other, CGGM-raw achieves superior estimation accuracy
even without finding the true cluster structure.

Blockdiagonal structure. In evaluating aggregation performance, CGGM-refit
demonstrates a sizeable advantage over alternative methods (see right two
columns in Figure 5). However, it exhibits a slightly higher estimation error
compared to CGGM-raw, TAGL-ideal, and TAGL-realistic. As in the designs with
an approximate block structure, it is worth noting that the identified block
structure is technically incorrect at the population level. Nonetheless,
identifying groups remains beneficial, considering the limited number of
nonzero elements outside the diagonal blocks. Furthermore, the performance of
CGGM-raw is noteworthy, as it surpasses all alternatives in terms of
estimation accuracy with minimal aggregation.

Conclusions from the simulations. CGGM-refit exhibits excellent performance
similar to TAGL-ideal while not requiring auxiliary information. Additionally,
the simulations reveal that a misspecification of the tree that encodes the
aggregation is detrimental to the performance of TAGL, positioning CGGM as a
particularly strong alternative when accurate information on the aggregation
structure is unavailable. Notably, CGGM consistently keeps up with—and even
surpasses—alternatives in simulation designs characterized by a clear sparsity
structure, despite the absence of a penalty term specifically aimed at
identifying such patterns. In general, the refitting step is beneficial, but
may lead to overaggregation in unstructured settings. Correspondingly, in some
scenarios, we find that there is a trade-off between CGGM-raw and CGGM-refit
in terms of estimation accuracy and aggregation performance.

##  4 Estimation of a Clustered Covariance Matrix

While this paper so far has been focused on estimating a clustered precision
matrix in order to simplify the partial dependency structure among the
variables, it can also be of interest to obtain a clustered covariance matrix.
For instance, a block structure in the covariance matrix implies that the
variables of a given block have equal loadings in latent factor models. In
this context, a clustered covariance matrix may reduce uncertainty in
estimating the latent factors.

Contrary to other approaches for finding block structures in the precision
matrix, such as those of Yao Allen (2019), Pircalabelu Claeskens (2020), and
Wilms Bien (2022), CGGM retains the found block structure when taking the
inverse of the obtained estimate (see Section 2.1). If model (1) holds for the
precision matrix $\mathbf{\Theta}$, it may therefore be reasonable to obtain a
clustered estimate of the covariance matrix $\mathbf{\Sigma}$ by taking
$\mathbf{\hat{\Sigma}}=\mathbf{\hat{\Theta}}^{-1}$, where
$\mathbf{\hat{\Theta}}$ denotes the estimate of the precision matrix from
CGGM. However, adherence to an exact block structure as in model (1) may be a
rather strong assumption in practice. If this model only holds approximately
with somewhat similar values within the blocks of $\mathbf{\Theta}$, we have
seen in the simulations that estimating a clustered precision matrix with CGGM
is nevertheless beneficial. But if an approximate block structure is present
in the covariance matrix $\mathbf{\Sigma}$, the structure in $\mathbf{\Theta}$
is much more noisy with larger variability between the elements in the
corresponding blocks, as illustrated in Figure 6. In such settings, first
estimating a clustered precision matrix with CGGM and then taking the inverse
may not succeed in finding the correct structure.

Figure 6: Illustration of how an exact block structure is retained between the
covariance matrix $\mathbf{\Sigma}$ and the precision matrix $\mathbf{\Theta}$
(left), whereas an approximate block structure in $\mathbf{\Sigma}$
corresponds to a much more noisy structure in $\mathbf{\Theta}$ (right).

To estimate a clustered covariance matrix, one could attempt to adapt the
clusterpath algorithm to be based on the objective function

$-\log|\mathbf{\Sigma}^{-1}|+\operatorname*{tr}\mathbf{S}\mathbf{\Sigma}^{-1}+\lambda\mathcal{P}(\mathbf{\Sigma}).$

An algorithm optimizing this objective function must deal with an additional
layer of complexity, as evaluation of $\mathbf{\Sigma}$ is required in the
penalty part and evaluation of $\mathbf{\Sigma}^{-1}$ in the likelihood part
of the objective function.

As an alternative, consider a random vector $\bm{X}$ that follows a normal
distribution with covariance matrix $\mathbf{\Sigma}$. Then there exists an
affine transformation matrix $\mathbf{A}=\mathbf{A}(\mathbf{\Sigma})$ so that
$\bm{Y}=\mathbf{A}^{\top}\bm{X}$ follows a normal distribution with covariance
matrix $\mathbf{\Theta}=\mathbf{\Sigma}^{-1}$. That is, $\mathbf{\Sigma}$ is
the precision matrix for the random vector $\bm{Y}$, which implies that we can
apply CGGM to minimize the objective function $L(\mathbf{\Sigma})$ from (5)
with the following minor modification. As we do not observe realizations of
$\bm{Y}$, we cannot compute the sample covariance matrix of those
realizations. Instead, we take the inverse of the sample covariance matrix
$\mathbf{S}^{-1}$ of the observed realizations of $\bm{X}$ as input for the
CGGM algorithm.

Figure 7: Results for the simulation designs regarding the estimation of a
clustered covariance matrix (in separate panels per column). Compare Figure 3
for explanatory notes. Aggregation performance is omitted for the sample
covariance matrix $\mathbf{S}$ as it does not perform any aggregation.

To study this further, we conduct simulations using the same settings as in
the baseline chain design—which employs an exact block structure—and the
approximate block structure chain design (see Section 3.1), but with the block
structure in the covariance matrix $\mathbf{\Sigma}$ rather than the precison
matrix $\mathbf{\Theta}$. We apply CGGM for obtaining a clustered precision
matrix followed by taking the inverse (denoted by
CGGM-$\mathbf{\hat{\Theta}}^{-1}$) and the modification for obtaining a
clustered covariance matrix (denoted by CGGM-$\mathbf{\hat{\Sigma}}$), as well
as the sample covariance matrix $\mathbf{S}$. Both variants of CGGM include
the parameter re-estimation step from Section 2.3.4. The results are shown in
Figure 7. Clearly, CGGM-$\mathbf{\hat{\Sigma}}$ succeeds in finding the
relevant structure in most replications whereas
CGGM-$\mathbf{\hat{\Theta}}^{-1}$ struggles to do so, even for the setting
with an exact block structure. While both variants of CGGM improve upon the
sample covariance matrix $\mathbf{S}$ in terms of the estimation error,
CGGM-$\mathbf{\hat{\Sigma}}$ exhibits the lowest error. Although a more
thorough evaluation is beyond the scope of this paper, our findings are an
indication that applying CGGM-$\mathbf{\hat{\Sigma}}$ may be preferable over
CGGM-$\mathbf{\hat{\Theta}}^{-1}$ when a clustered covariance matrix is of
primary interest. We therefore recommend users to first decide on the object
of interest, a block-structured covariance matrix or precision matrix, and
then use CGGM accordingly.

##  5 Applications

We demonstrate the practical usefulness and versatility of CGGM on three
applications: a finance application using stock data in Section 5.1, an
application with country-level well-being indicators in Section 5.2, and a
survey-based behavioral science application in Section 5.3.

###  5.1 S&P 100 Stocks

Consider a financial data set containing price information of $p=101$ stocks
of the companies that constituted the S&P 100 on September 18th, 2023. Daily
stock price information was collected for the period ranging from January 3rd,
2023, until December 29th, 2023 ($n=250$).222The S&P 100 contains 100
companies that may have one or more classes of stocks listed. During the
timeframe of our analysis, Alphabet was the sole company with more than one
class, distinguished by the symbols GOOG and GOOGL, therefore resulting in
$p=101$ instead of $p=100$ variables in our analysis. We compute daily
realized ranges as given by

$r_{tj}=\frac{(\log H_{tj}-\log L_{tj})^{2}}{4\log 2},$

where $H_{tj}$ and $L_{tj}$ are the high and low prices for stock $j$
($j=1,\ldots,p$) during trading day $t$ (Parkinson, 1980) and study the
conditional dependency structure of the stocks’ realized ranges. As a
preprocessing step, we first fit the popular heterogeneous autoregressive
(HAR) model of Corsi (2009) to the individual daily realized range series to
capture time dependencies, and then apply CGGM to the standardized residual
series to learn the conditional dependency structure among the stocks; our
procedure is in line with Wilms Bien (2022).

We compare the performance of CGGM to that of TAGL. While the former learns
how to cluster the variables in an unsupervised data-driven way, the latter
requires side-information in the form of a tree that encodes the similarity
between the variables stocks to do so. To this end, we use the Global Industry
Classification Standard (GICS). The tree then consists of the $p=101$ stocks
(leaves), 11 sectors as middle layer where each stock/company belongs to one
industry sector, and one root node that aggregates all sectors. For both
methods, we use 5-fold cross-validation to select the tuning parameters ($k$,
$\phi$, and $\lambda$ for CGGM; $\lambda_{1}$ and $\lambda_{2}$ for TAGL) and
refit the precision matrix subject to the obtained variable clustering. For
the weight matrix in CGGM, we use a grid of $\phi\in\\{2,6,10,14,18,22\\}$ and
$k\in\\{1,3,5,7,9,11,13,15\\}$ and for TAGL the tuning parameter grid is
determined via the same procedure as described in Section 3.2.

The cluster solution returned by both procedures differs considerably. CGGM
groups the stocks into $K=3$ clusters: the first cluster contains BAC (issued
by Bank of America Corp.), the second cluster contains GOOG and GOOGL (both
issued by Alphabet), and the last one contains the remaining stocks. That is,
CGGM puts the Bank of America and the technology company Alphabet in the
spotlight as opposed to the rest of the market. TAGL, in contrast, results in
$K=9$ clusters where the industry sectors are assigned to the different
clusters on an almost one-to-one basis. In the left panel of Figure 8, we
present the distribution of the sectors over the clusters for both CGGM and
TAGL. The unbalanced cluster solution of the former versus the more balanced
one of the latter directly stands out. The role of the stocks BAC, GOOG, and
GOOGL in our analysis is in line with the central, market-wide role of banks
(e.g., Martins, 2023) and technology companies (e.g., Cardoso ., 2020) since
their influence on prices affects stocks across the whole market.

Given the considerable difference in clusters returned by both methods, a
natural question is how both perform in capturing the conditional dependency
structure. To this end, we conduct an out-of-sample exercise to compare their
performance. An additional outer loop for 10-fold cross-validation is used to
compute out-of-sample errors on each of the $G=10$ test samples according to
the likelihood-based score in equation (9). These errors on the test data are
visualized in the right panel of Figure 8 with CGGM on the horizontal axis and
TAGL on the vertical axis. For each of the test samples (dots), CGGM has a
lower error thereby indicating that the more unbalanced cluster structure
forms a better description of the conditional dependency structure in this
context.

Figure 8: Finance application. Distribution of the sectors over the $K=3$
cluster solution of CGGM and the $K=9$ cluster solution of TAGL (left). Out-
of-sample errors across the ten replications (dots) for CGGM versus TAGL
(right).

###  5.2 OECD Well-Being Indicators

The second data set we analyze is one on OECD well-being indicators, collected
during 2018. These data were previously analyzed in Cavicchia . (2022), and we
acknowledge the authors for providing us with data access. The data contain
$p=11$ variables related to well-being: education, jobs, income, safety,
health, environment, civic engagement, accessibility to services, housing,
community, and life satisfaction, on which $n=36$ countries are given a score
ranging from 0 to 10. Note, however, that we split the sample of countries
into two groups since Cavicchia . (2022) have found evidence for two groups of
countries ($n_{1}=21$ and $n_{2}=14$) with distinct clustering
structures.333The first group consists of the countries Australia, Austria,
Belgium, Canada, Denmark, Finland, France, Germany, Iceland, Ireland, Italy,
Japan, Luxembourg, Netherlands, New Zealand, Norway, Spain, Sweden,
Switzerland, United Kingdom, United States, the second group consists of the
countries Chile, Czech Republic, Estonia, Greece, Hungary, Israel, Korea,
Latvia, Lithuania, Mexico, Poland, Portugal, Slovak Republic, Slovenia,
Turkey. In our analysis, we omit Lithuania from the second group due to a
missing value, hence $n_{1}=21$ and $n_{2}=14$.

This application highlights the capabilities of CGGM in estimating a cluster
hierarchy on the standardized well-being data. We thus report the whole
clusterpath solution obtained by applying CGGM to estimate the precision
matrix with different values of the tuning parameter $\lambda$ for the two
country groups. Moreover, due to the small sample sizes in both groups, we do
not use cross-validation to select the tuning parameters $\phi$ and $k$ for
the weight matrix but use fixed values based on existing conventions instead,
namely $\phi=0.5$ and $k=3$ (see e.g., Chi Lange, 2015; Sun ., 2021; Wang .,
2018). Since our primary interest is in the clusterpaths, the refitting step
described in Section 2.3.4 is not required as it would not affect the obtained
paths. Figure 9 visualizes the dendrograms for the two groups of countries as
obtained from CGGM’s complete clusterpath solution ranging from $p=11$
clusters to one.

Figure 9: Well-being application. Clusterpath dendrograms obtained by CGGM on
the two groups of countries.

The analysis of well-being indicators across the two groups of countries
reveals distinct clustering patterns that reflect their differing levels of
development. In the first group (e.g., Australia, France, Switzerland), the
aggregation of civic engagement, life satisfaction, community, environment,
housing, accessibility to services, and jobs appears to reflect the social and
environmental dimensions of well-being. Conversely, in the second group (e.g.,
Chile, Hungary, Mexico), safety and health are combined with the variables
representing social and environmental dimensions, possibly suggesting a
greater reliance on the community for these services. This contrasts with more
developed countries where such services may be perceived as more self-evident.
Additionally, in the second group, income, jobs, and education are grouped
together, reflecting their importance to socioeconomic progress. Finally, note
that we use CGGM with the precision matrix as object of interest; its block
structure directly transfers to the covariance matrix thereby making our
results comparable to those obtained in Cavicchia . (2022), who estimate the
clustering structure on the covariance matrix.

###  5.3 Humor Styles Questionnaire

In surveys in the behavioral sciences, multiple rating-scale items are
typically used to measure each latent construct of interest. For further
analysis, it is common practice to compute mean scores across the respective
items as measurements of each latent construct, which implies the assumption
that the covariance matrix of the questionnaire items follows a block
structure. The humor styles questionnaire (HSQ) developed by Martin . (2003)
measures four latent constructs corresponding to different styles of humor:
_affiliative_ , _self-enhancing_ , _aggressive_ , and _self-defeating_ , with
each construct being measured by eight items. An overview of the items is
presented in Table LABEL:tab:HSQ in Appendix D. If we can retrieve these four
clusters of items in the covariance matrix, the assumption behind computing
mean scores can be considered reasonable. We use responses to the HSQ on a
five-point rating scale (anchored by 1 = “never or very rarely true” to 5 =
“very often or always true”), which we obtained from
https://openpsychometrics.org/_rawdata/. In addition to the aforementioned
$p=32$ items, participants were asked at the end of the questionnaire to
indicate the accuracy of their responses. We restrict our analysis to
participants who reported that their responses are fully accurate, and after
removing 9 observations with missing responses, we retain $n=182$ respondents.

Using 5-fold cross-validation, we apply the CGGM algorithm for clustering the
covariance matrix (including the parameter re-estimation step) with a grid of
$\phi\in\\{1,1.5,2,2.5,3\\}$ and $k\in\\{1,2,3,4,5\\}$ for the weight matrix.
We obtain four clusters of items that largely overlap with the grouping of
Martin . (2003), as shown in Table 1. Looking at the wording of the
questionnaire items, the clustering found by CGGM is intuitive. All items of
the _affiliative_ group, which CGGM places in Cluster 1, refer to humor in
social settings. In addition, the two items of the _self-enhancing_ group that
complete Cluster 1—items 6 and 30—address being alone as opposed to being with
people. On the other hand, the remaining six items of the _self-enhancing_
group, which CGGM puts in Cluster 2, all refer to using humor for dealing with
feeling depressed, upset, unhappy, sad, or having problems. Item 28 from the
_self-defeating_ group, which completes Cluster 2, also addresses having
problems or feeling unhappy. In contrast, the other seven items of this group,
which form Cluster 4 in the CGGM solution, refer to laughing or letting others
laugh at oneself. Finally, all items from the _aggressive_ group address
teasing and offensive or inappropriate humor, and are placed in Cluster 3 by
CGGM.

|  Cluster 1 |  Cluster 2 |  Cluster 3 |  Cluster 4
---|---|---|---|---
Affiliative |  {1, 5, 9, 13, 17, 21, 25, 29} |  |  |
Self-enhancing |  {6, 30} |  {2, 10, 14, 18, 22, 26} |  |
Aggressive |  |  |  {3, 7, 11, 15, 19, 23, 27, 31} |
Self-defeating |  |  {28} |  |  {4, 8, 12, 16, 20, 24, 32}
Table 1: Humor styles questionnaire application. Allocation of the items
developed by Martin . (2003) for measuring the four humor styles (in rows) to
the four clusters found by CGGM for clustering the covariance matrix (in
columns).

Although the above interpretation of the obtained clusters is subjective, our
findings suggest that the survey design of the HSQ may require further
finetuning. Moreover, if subsequent analysis is based on mean scores, a
slightly different allocation of the items to the latent constructs—the humor
styles—may be more in line with the implicit assumption of a block structure
in the covariance matrix.

##  6 Conclusion

We introduce a novel method to estimate Gaussian graphical models (GGMs)
subject to node-clustering. Our method, which we call CGGM, uses a clusterpath
penalty to produce a hierarchical clustering of variables (nodes in the graph)
without relying on pre-existing notions of the cluster composition. An R
package called CGGMR implements the proposed method and is available on the
GitHub page of the first author.

In a comprehensive simulation study covering a wide range of graph structures,
we compare CGGM to similar benchmarks such as TAGL and ComGGL. CGGM oftentimes
surpasses the benchmarks both in terms of estimation accuracy and clustering
performance, though the latter requires a refitting step particularly in noisy
simulation designs. Through a diverse set of applications, we also demonstrate
the versatility of CGGM in (i) learning cluster hierarchies in a fully data-
driven way as compared to other benchmarks such as TAGL that require
additional side-information for this task, (ii) delivering accurate cluster
hierarchies that are transferable between the precision matrix and covariance
matrix, and (iii) directly estimating clustered covariance matrices if these
are the primary object of interest.

Concerning future work, CGGM may be extended to partial correlation matrices
instead of precision matrices as the object of study. In contrast to the
graphical lasso (glasso; Friedman ., 2008), Carter . (2024) apply the sparsity
penalty to the elements of the partial correlation matrix and call their
method the partial correlation glasso. As a result, the sparsity structure
estimated by this method is invariant to scalar multiplication of the data.
CGGM can also be modified to accommodate clustering based on the partial
correlation matrix, but this comes at the cost of the convexity of the
objective function as with the partial correlation glasso.

## Computational Details

All computations were performed with R (R Core Team, 2024). _Replication files
will be made publicly available upon acceptance and a link will be included
here._

## Acknowledgments

Andreas Alfons is supported by a grant from the Dutch Research Council (NWO),
research program Vidi, grant number VI.Vidi.195.141. Ines Wilms is supported
by the same funding agency under grant number VI.Vidi.211.032.

## References

  * Ambroise . (2009) Ambroise2009InferringSGAmbroise, C., Chiquet, J. Matias, C.  2009\.  Inferring sparse Gaussian graphical models with latent structure Inferring sparse Gaussian graphical models with latent structure.  Electronic Journal of Statistics3205–238.  https://doi.org/10.1214/08-EJS314
  * Archakov Hansen (2022) archakov2022canonicalArchakov, I. Hansen, PR.  2022\.  A canonical representation of block matrices with applications to covariance and correlation matrices A canonical representation of block matrices with applications to covariance and correlation matrices.  The Review of Economics and Statistics1–39.  https://doi.org/10.1162/rest_a_01258
  * Banerjee . (2008) banerjee2008model_selectionBanerjee, O., El Ghaoui, L. d’Aspremont, A.  2008\.  Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data.  Journal of Machine Learning Research9485–516.
  * Brownlees . (2022) brownlees2022communityBrownlees, C., Guðmundsson, GS. Lugosi, G.  2022\.  Community detection in partial correlation network models Community detection in partial correlation network models.  Journal of Business & Economic Statistics401216–226.  https://doi.org/10.1080/07350015.2020.1798241
  * Bunea . (2020) bunea2020modelBunea, F., Giraud, C., Luo, X., Royer, M. Verzelen, N.  2020\.  Model assisted variable clustering: minimax-optimal recovery and algorithms Model assisted variable clustering: minimax-optimal recovery and algorithms.  Annals of Statistics481111–137.  https://doi.org/10.1214/18-AOS1794
  * Cai . (2011) cai2011constrainedCai, T., Liu, W. Luo, X.  2011\.  A constrained $\ell_{1}$ minimization approach to sparse precision matrix estimation A constrained $\ell_{1}$ minimization approach to sparse precision matrix estimation.  Journal of the American Statistical Association106494594–607.  https://doi.org/10.1198/jasa.2011.tm10155
  * Cardoso . (2020) cardoso2020algorithmsCardoso, JVdM., Ying, J. Palomar, DP.  2020\.  Algorithms for learning graphs in financial markets Algorithms for learning graphs in financial markets.  arXiv preprint arXiv:2012.15410.  https://doi.org/10.48550/arXiv.2012.15410
  * Carter . (2024) carter2024partialCarter, JS., Rossell, D. Smith, JQ.  2024\.  Partial correlation graphical lasso Partial correlation graphical lasso.  Scandinavian Journal of Statistics51132-63.  https://doi.org/10.1111/sjos.12675
  * Cavicchia . (2022) cavicchia2022gaussianCavicchia, C., Vichi, M. Zaccaria, G.  2022\.  Gaussian mixture model with an extended ultrametric covariance structure Gaussian mixture model with an extended ultrametric covariance structure.  Advances in Data Analysis and Classification162399–427.  https://doi.org/10.1007/s11634-021-00488-x
  * Chi . (2017) chi2017cvxbiclustChi, EC., Allen, GI. Baraniuk, RG.  2017\.  Convex Biclustering Convex biclustering.  Biometrics73110–19.  https://doi.org/10.1111/biom.12540
  * Chi Lange (2015) chi2015amaChi, EC. Lange, KL.  2015\.  Splitting methods for convex clustering Splitting methods for convex clustering.  Journal of Computational and Graphical Statistics244994-1013.  https://doi.org/10.1080/10618600.2014.948181
  * Corsi (2009) corsi2009simpleCorsi, F.  2009\.  A simple approximate long-memory model of realized volatility A simple approximate long-memory model of realized volatility.  Journal of Financial Econometrics72174–196.  https://doi.org/10.1093/jjfinec/nbp001
  * Danaher . (2014) danaher2014jointDanaher, P., Wang, P. Witten, DM.  2014\.  The joint graphical lasso for inverse covariance estimation across multiple classes The joint graphical lasso for inverse covariance estimation across multiple classes.  Journal of the Royal Statistical Society: Series B762373.  https://doi.org/10.1111/rssb.12033
  * Duchi . (2008) duchi2008projectedDuchi, J., Gould, S. Koller, D.  2008\.  Projected subgradient methods for learning sparse Gaussians Projected subgradient methods for learning sparse Gaussians.  Proceedings of the Twenty-Fourth Conference on Uncertainty in Artificial Intelligence Proceedings of the twenty-fourth conference on uncertainty in artificial intelligence ( 153–160).
  * Eisenach . (2020) eisenach2020highEisenach, C., Bunea, F., Ning, Y. Dinicu, C.  2020\.  High-dimensional inference for cluster-based graphical models High-dimensional inference for cluster-based graphical models.  Journal of Machine Learning Research21531-55.
  * Friedman . (2008) friedman2008glassoFriedman, J., Hastie, T. Tibshirani, R.  2008\.  Sparse Inverse Covariance Estimation with the Graphical Lasso Sparse inverse covariance estimation with the graphical lasso.  Biostatistics93432–441.  https://doi.org/10.1093/biostatistics/kxm045
  * Gao . (2016) gao2016estimationGao, C., Zhu, Y., Shen, X. Pan, W.  2016\.  Estimation of multiple networks in Gaussian mixture models Estimation of multiple networks in Gaussian mixture models.  Electronic Journal of Statistics1011133–1154.  https://doi.org/10.1214/16-EJS1135
  * Gheche Frossard (2020) gheche2020multilayerGheche, ME. Frossard, P.  2020\.  Multilayer clustered graph learning Multilayer clustered graph learning.  arXiv preprint arXiv:2010.15456.  https://doi.org/10.48550/arXiv.2010.15456
  * Gower Groenen (1991) gower1991structuredGower, JC. Groenen, PJF.  1991\.  Applications of the modified Leverrier-Faddeev algorithm for the construction of explicit matrix decompositions and inverses Applications of the modified Leverrier-Faddeev algorithm for the construction of explicit matrix decompositions and inverses.  Utilitas Mathematica4051–64.
  * Grechkin . (2015) grechkin2015pathwayGrechkin, M., Fazel, M., Witten, D. Lee, SI.  2015\.  Pathway graphical lasso Pathway graphical lasso.  Proceedings of the AAAI Conference on Artificial Intelligence Proceedings of the AAAI Conference on Artificial Intelligence ( 29).  https://doi.org/10.1609/aaai.v29i1.9636
  * Hao . (2018) hao2018simultaneousHao, B., Sun, WW., Liu, Y. Cheng, G.  2018\.  Simultaneous clustering and estimation of heterogeneous graphical models Simultaneous clustering and estimation of heterogeneous graphical models.  Journal of Machine Learning Research182171–58.
  * Heinävaara . (2016) heinavaara2016inconsistencyHeinävaara, O., Leppä-Aho, J., Corander, J. Honkela, A.  2016\.  On the inconsistency of $\ell_{1}$-penalised sparse precision matrix estimation On the inconsistency of $\ell_{1}$-penalised sparse precision matrix estimation.  BMC Bioinformatics171699–107.  https://doi.org/10.1186/s12859-016-1309-x
  * Hocking . (2011) hocking2011convex_clusteringHocking, TD., Joulin, A., Bach, F. Vert, JP.  2011\.  Clusterpath: An algorithm for clustering using convex fusion penalties Clusterpath: An algorithm for clustering using convex fusion penalties.  The 28th International Conference on Machine Learning. The 28th International Conference on Machine Learning.  Bellevue, Washington.
  * Hosseini Lee (2016) hosseini2016learningHosseini, MJ. Lee, SI.  2016\.  Learning sparse Gaussian graphical models with overlapping blocks Learning sparse Gaussian graphical models with overlapping blocks.  Advances in Neural Information Processing Systems Advances in Neural Information Processing Systems ( 29).
  * Hubert Arabie (1985) hubert1985Hubert, L. Arabie, P.  1985\.  Comparing partitions Comparing partitions.  Journal of Classification21193–218.  https://doi.org/10.1007/BF01908075
  * Kruskal (1956) kruskal1956mstKruskal, JB.  1956\.  On the shortest spanning subtree of a graph and the traveling salesman problem On the shortest spanning subtree of a graph and the traveling salesman problem.  Proceedings of the American Mathematical society7148–50.
  * Kumar . (2020) kumar2020unifiedKumar, S., Ying, J., De M. Cardoso, JV. Palomar, DP.  2020\.  A unified framework for structured graph learning via spectral constraints A unified framework for structured graph learning via spectral constraints.  Journal of Machine Learning Research2111–60.
  * Lindsten . (2011) lindsten2011convex_clusteringLindsten, F., Ohlsson, H. Ljung, L.  2011\.  Just relax and come clustering!: A convexification of $k$-means clustering Just relax and come clustering!: A convexification of $k$-means clustering .  Linköping, SwedenDepartment of Electrical Engineering, Linköping University.
  * Marlin Murphy (2009) marlin2009sparseMarlin, BM. Murphy, KP.  2009\.  Sparse Gaussian graphical models with unknown block structure Sparse Gaussian graphical models with unknown block structure.  Proceedings of the 26th Annual International Conference on Machine Learning Proceedings of the 26th Annual International Conference on Machine Learning ( 705–712).  https://doi.org/10.1145/1553374.1553465
  * Martin . (2003) martin2003HSQMartin, RA., Puhlik-Doris, P., Larsen, G., Gray, J. Weir, K.  2003\.  Individual differences in uses of humor and their relation to psychological well-being: Development of the Humor Styles Questionnaire Individual differences in uses of humor and their relation to psychological well-being: Development of the humor styles questionnaire.  Journal of Research in Personality37148–75.  https://doi.org/10.1016/S0092-6566(02)00534-2
  * Martins (2023) martins2023stockMartins, AM.  2023\.  Stock market effects of silicon valley bank and credit suisse failure: Evidence for a sample of European listed banks Stock market effects of silicon valley bank and credit suisse failure: Evidence for a sample of European listed banks.  Finance Research Letters58104296.  https://doi.org/10.1016/j.frl.2023.104296
  * Meinshausen Bühlmann (2006) meinshausen2006highMeinshausen, N. Bühlmann, P.  2006\.  High-dimensional graphs and variable selection with the lasso High-dimensional graphs and variable selection with the lasso.  The Annals of Statistics3431436–1462.  https://doi.org/10.1214/009053606000000281
  * Millington Niranjan (2019) millington2019quantifyingMillington, T. Niranjan, M.  2019\.  Quantifying influence in financial markets via partial correlation network inference Quantifying influence in financial markets via partial correlation network inference.  2019 11th International Symposium on Image and Signal Processing and Analysis (ISPA) 2019 11th International Symposium on Image and Signal Processing and Analysis (ISPA) ( 306–311).  https://doi.org/10.1109/ISPA.2019.8868437
  * Mohan . (2014) mohan2014nodeMohan, K., London, P., Fazel, M., Witten, D. Lee, SI.  2014\.  Node-based learning of multiple Gaussian graphical models Node-based learning of multiple Gaussian graphical models.  Journal of Machine Learning Research151445–488.
  * Nie . (2016) nie2016constrainedNie, F., Wang, X., Jordan, M. Huang, H.  2016\.  The constrained Laplacian rank algorithm for graph-based clustering The constrained Laplacian rank algorithm for graph-based clustering.  Proceedings of the AAAI Conference on Artificial Intelligence Proceedings of the AAAI Conference on Artificial Intelligence ( 30).  https://doi.org/10.1609/aaai.v30i1.10302
  * Parkinson (1980) parkinson1980extremeParkinson, M.  1980\.  The extreme value method for estimating the variance of the rate of return The extreme value method for estimating the variance of the rate of return.  Journal of Business53161–65.
  * Pelckmans . (2005) pelckmans2005convex_clusteringPelckmans, K., De Brabanter, J., Suykens, JAK. de Moor, B.  2005\.  Convex clustering shrinkage Convex clustering shrinkage.  PASCAL Workshop on Statistics and Optimization of Clustering Workshop. PASCAL Workshop on Statistics and Optimization of Clustering Workshop.
  * Peng . (2009) peng2009partialPeng, J., Wang, P., Zhou, N. Zhu, J.  2009\.  Partial correlation estimation by joint sparse regression models Partial correlation estimation by joint sparse regression models.  Journal of the American Statistical Association104486735–746.  https://doi.org/10.1198/jasa.2009.0126
  * Pircalabelu Claeskens (2020) pircalabelu2020communityPircalabelu, E. Claeskens, G.  2020\.  Community-based group graphical lasso Community-based group graphical lasso.  Journal of Machine Learning Research2111–32.
  * R Core Team (2024) RR Core Team.  2024\.  R: A Language and Environment for Statistical Computing R: A language and environment for statistical computing [].  Vienna, Austria.  https://www.R-project.org/
  * Rothman . (2008) rothman2008covar_estimationRothman, AJ., Bickel, PJ., Levina, E. Zhu, J.  2008\.  Sparse permutation invariant covariance estimation Sparse permutation invariant covariance estimation.  Electronic Journal of Statistics2494–515.  https://doi.org/10.1214/08-EJS176
  * Schmidt . (2009) schmidt2009optimizingSchmidt, M., Berg, E., Friedlander, M. Murphy, K.  2009\.  Optimizing costly functions with simple constraints: A limited-memory projected quasi-newton algorithm Optimizing costly functions with simple constraints: A limited-memory projected quasi-newton algorithm.  Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics ( 5, 456–463).
  * Shi . (2024) shi2024simultaneousShi, D., Wang, T. Ying, Z.  2024\.  Simultaneous Identification of Sparse Structures and Communities in Heterogeneous Graphical Models Simultaneous identification of sparse structures and communities in heterogeneous graphical models.  arXiv preprint arXiv:2405.09841.  https://doi.org/10.48550/arXiv.2405.09841
  * Sun . (2021) sun2021ssnalSun, D., Toh, KC. Yuan, Y.  2021\.  Convex clustering: Model, theoretical guarantee and efficient algorithm Convex clustering: Model, theoretical guarantee and efficient algorithm.  Journal of Machine Learning Research2291–32.
  * Tan . (2015) tan2015clusterTan, KM., Witten, D. Shojaie, A.  2015\.  The cluster graphical lasso for improved estimation of Gaussian graphical models The cluster graphical lasso for improved estimation of Gaussian graphical models.  Computational Statistics & Data Analysis8523–36.  https://doi.org/10.1016/j.csda.2014.11.015
  * Tarzanagh . (2021) tarzanagh2021fairTarzanagh, DA., Balzano, L. Hero, AO.  2021\.  Fair structure learning in heterogeneous graphical models Fair structure learning in heterogeneous graphical models.  arXiv preprint arXiv:2112.05128.  https://doi.org/10.48550/arXiv.2112.05128
  * Tarzanagh Michailidis (2018) tarzanagh2018estimationTarzanagh, DA. Michailidis, G.  2018\.  Estimation of graphical models through structured norm minimization Estimation of graphical models through structured norm minimization.  Journal of Machine Learning Research1811–48.
  * Touw . (2022) touw2023convexTouw, DJW., Groenen, PJF. Terada, Y.  2022\.  Convex clustering through MM: An efficient algorithm to perform hierarchical clustering Convex clustering through MM: An efficient algorithm to perform hierarchical clustering.  arXiv preprint arXiv:2211.01877.  https://doi.org/10.48550/arXiv.2211.01877
  * Wang . (2018) wang2018sparseWang, B., Zhang, Y., Sun, WW. Fang, Y.  2018\.  Sparse convex clustering Sparse convex clustering.  Journal of Computational and Graphical Statistics272393–403.  https://doi.org/10.1080/10618600.2017.1377081
  * Wilms Bien (2022) wilms2021tag_lassoWilms, I. Bien, J.  2022\.  Tree-based Node Aggregation in Sparse Graphical Models Tree-based node aggregation in sparse graphical models.  Journal of Machine Learning Research232431–36.
  * Yao Allen (2019) yao2019clusteredYao, T. Allen, GI.  2019\.  Clustered Gaussian graphical model via symmetric convex clustering Clustered Gaussian graphical model via symmetric convex clustering.  2019 IEEE Data Science Workshop (DSW) 2019 IEEE Data Science Workshop (DSW) ( 76–82).  https://doi.org/10.1109/DSW.2019.8755774
  * Yuan (2010) yuan2010highYuan, M.  2010\.  High dimensional inverse covariance matrix estimation via linear programming High dimensional inverse covariance matrix estimation via linear programming.  Journal of Machine Learning Research11792261–2286.
  * Yuan Lin (2006) yuan2006groupYuan, M. Lin, Y.  2006\.  Model selection and estimation in regression with grouped variables Model selection and estimation in regression with grouped variables.  Journal of the Royal Statistical Society: Series B68149–67.  https://doi.org/10.1111/j.1467-9868.2005.00532.x
  * Yuan Lin (2007) yuan2007graphical_modelYuan, M. Lin, Y.  2007\.  Model selection and estimation in the Gaussian graphical model Model selection and estimation in the Gaussian graphical model.  Biometrika94119–35.  https://doi.org/10.1093/biomet/asm018

##  Appendix A Derivations Clusterpath Estimator

The main purpose of re-expressing the objective function $L(\mathbf{\Theta})$
in terms of $\mathbf{A}$ and $\mathbf{R}$ is to obtain a more compact
expression that permits straightforward updates for solving optimization
problem (1). To this end, we begin by defining the key variables and then work
our way through the different terms of the objective function (3).

For any number of clusters $1\leq K\leq p$, the precision matrix contains the
block-structure

$\mathbf{\Theta}=\begin{bmatrix}(a_{11}-r_{11})\mathbf{I}+r_{11}\mathbf{11}^{\top}&r_{12}\mathbf{11}^{\top}&\cdots&r_{1K}\mathbf{11}^{\top}\\\
r_{21}\mathbf{11}^{\top}&(a_{22}-r_{22})\mathbf{I}+r_{22}\mathbf{11}^{\top}&\cdots&r_{2K}\mathbf{11}^{\top}\\\
\vdots&\vdots&\ddots&\vdots\\\
r_{K1}\mathbf{11}^{\top}&r_{K2}\mathbf{11}^{\top}&\cdots&(a_{KK}-r_{KK})\mathbf{I}+r_{KK}\mathbf{11}^{\top}\\\
\end{bmatrix}.$

Without loss of generality, we assume that the variables have been reordered
such that those in cluster 1 form the first block, those in cluster 2 the
second block, and so on. Let $\mathbf{U}$ be the $p\times K$ membership matrix
with $u_{jk}=1$ if variable $j$ belongs to cluster $k$ and zero otherwise.
Then, $\mathbf{P}=\mathbf{U}^{\top}\mathbf{U}$ is the diagonal matrix with the
number of variables $p_{k}$ per cluster as diagonal elements and
$\mathbf{p}=\mathbf{1}^{\top}\mathbf{U}$ is the $K\times 1$ vector with
elements $p_{k}$. It can be verified that
$\mathbf{X}\mathbf{U}\mathbf{P}^{-1}$ computes the cluster averages of the
original variables in $\mathbf{X}$.

First, we focus on the first term in objective function (3), the log
determinant of $\mathbf{\Theta}$, to obtain a convenient expression in
$\mathbf{A}$ and $\mathbf{R}$. We show that $\mathbf{\Theta}$ can be split in
a linear space that depends on $\mathbf{A}$ and $\mathbf{R}$ and an orthogonal
linear space that is only dependent on the diagonal blocks of
$\mathbf{\Theta}$ and thus on $\mathbf{A}$ and the diagonal elements of
$\mathbf{R}$. To do so, we first consider what happens if we take cluster
averages, that is,

$\displaystyle\mathbf{P}^{-1}\mathbf{U}^{\top}\bm{\Theta}\mathbf{U}\mathbf{P}^{-1}$
$\displaystyle=\mathbf{P}^{-1}\mathbf{U}^{\top}\left(\mathbf{URU}^{\top}+\begin{bmatrix}(a_{11}-r_{11})\mathbf{I}&\mathbf{0}&\cdots&\mathbf{0}\\\
\mathbf{0}&(a_{22}-r_{22})\mathbf{I}&\cdots&\mathbf{0}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{0}&\mathbf{0}&\cdots&(a_{KK}-r_{KK})\mathbf{I}\\\
\end{bmatrix}\right)\mathbf{U}\mathbf{P}^{-1}$
$\displaystyle=\mathbf{R}+\mathbf{P}^{-1}\begin{bmatrix}(a_{11}-r_{11})&0&\cdots&0\\\
0&(a_{22}-r_{22})&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&(a_{KK}-r_{KK})\\\ \end{bmatrix}$
$\displaystyle=\mathbf{R}+\begin{bmatrix}(a_{11}-r_{11})/p_{1}&0&\cdots&0\\\
0&(a_{22}-r_{22})/p_{2}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&(a_{KK}-r_{KK})/p_{KK}\\\ \end{bmatrix}$
$\displaystyle=\mathbf{R}^{*}.$

Then,

$\mathbf{UR}^{*}\mathbf{U}^{\top}=\begin{bmatrix}r_{11}^{*}\mathbf{11}^{\top}&r_{12}\mathbf{11}^{\top}&\cdots&r_{1K}\mathbf{11}^{\top}\\\
r_{21}\mathbf{11}^{\top}&r_{22}^{*}\mathbf{11}^{\top}&\cdots&r_{2K}\mathbf{11}^{\top}\\\
\vdots&\vdots&\ddots&\vdots\\\
r_{K1}\mathbf{11}^{\top}&r_{K2}\mathbf{11}^{\top}&\cdots&r_{KK}^{*}\mathbf{11}^{\top}\\\
\end{bmatrix},$

where $r_{kk}^{*}=(a_{kk}-r_{kk})/p_{k}+r_{kk}$ so that
$\bm{\Theta}-\mathbf{UR}^{*}\mathbf{U}^{\top}$ has off-diagonal blocks equal
to zero and diagonal blocks equal to
$(a_{kk}-r_{kk})(\mathbf{I}-p_{k}^{-1}\mathbf{11}^{\top})=(a_{kk}-r_{kk})\mathbf{J}_{k}$
with $\mathbf{J}_{k}$ the $p_{k}\times p_{k}$ centering matrix. This yields

$\mathbf{\Theta}=\mathbf{UR}^{*}\mathbf{U}^{\top}+\begin{bmatrix}(a_{11}-r_{11})\mathbf{J}_{1}&\mathbf{0}&\cdots&\mathbf{0}\\\
\mathbf{0}&(a_{22}-r_{22})\mathbf{J}_{2}&\cdots&\mathbf{0}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{0}&\mathbf{0}&\cdots&(a_{KK}-r_{KK})\mathbf{J}_{K}\\\ \end{bmatrix},$
(10)

as also displayed in (4). It can be verified that this decomposition is indeed
orthogonal because post-multiplying the final term with $\mathbf{U}$ yields
zero due to the centering matrices.

Next, it is well known that the determinant is equal to the product of the
eigenvalues of the matrix. Therefore, we require the eigenvalues of both parts
of decomposition (10). Let the eigendecomposition of $\mathbf{\Theta}$ be
given by

$\displaystyle\mathbf{\Theta}$
$\displaystyle=\mathbf{Q}\mathbf{\Gamma}\mathbf{Q}^{\top}$
$\displaystyle=\mathbf{Q}_{1}\mathbf{\Gamma}_{1}\mathbf{Q}_{1}^{\top}+\mathbf{Q}_{2}\mathbf{\Gamma}_{2}\mathbf{Q}_{2}^{\top}$
$\displaystyle=\mathbf{UR}^{*}\mathbf{U}^{\top}+\begin{bmatrix}(a_{11}-r_{11})\mathbf{J}_{1}&\mathbf{0}&\cdots&\mathbf{0}\\\
\mathbf{0}&(a_{22}-r_{22})\mathbf{J}_{2}&\cdots&\mathbf{0}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{0}&\mathbf{0}&\cdots&(a_{KK}-r_{KK})\mathbf{J}_{K}\\\ \end{bmatrix}.$

We now need to find $\mathbf{\Gamma}_{1}$ in the eigendecomposition
$\mathbf{UR}^{*}\mathbf{U}^{\top}=\mathbf{Q}_{1}\mathbf{\Gamma}_{1}\mathbf{Q}_{1}^{\top}$.
Pre- and post-multiplying with an orthonormal matrix does not change the
eigenvalues. A convenient orthonormal matrix for pre-multiplication is
therefore $\mathbf{P}^{-1/2}\mathbf{U}$ giving

$\displaystyle\mathbf{UR}^{*}\mathbf{U}^{\top}$
$\displaystyle=\mathbf{Q}_{1}\mathbf{\Gamma}_{1}\mathbf{Q}_{1}^{\top}$
$\displaystyle\mathbf{P}^{-1/2}\mathbf{U}^{\top}\mathbf{UR}^{*}\mathbf{U}^{\top}\mathbf{U}\mathbf{P}^{-1/2}$
$\displaystyle=\mathbf{P}^{-1/2}\mathbf{U}^{\top}\mathbf{Q}_{1}\mathbf{\Gamma}_{1}\mathbf{Q}_{1}^{\top}\mathbf{U}\mathbf{P}^{-1/2}$
$\displaystyle\mathbf{P}^{1/2}\mathbf{R}^{*}\mathbf{P}^{1/2}$
$\displaystyle=\mathbf{Q}_{1}^{*}\mathbf{\Gamma}_{1}{\mathbf{Q}_{1}^{*}}^{\top}$
$\displaystyle\begin{bmatrix}p_{1}r_{11}^{*}&(p_{1}p_{2})^{1/2}r_{12}&\cdots&(p_{1}p_{K})^{1/2}r_{1K}\\\
(p_{1}p_{2})^{1/2}r_{21}&p_{2}r_{22}^{*}&\cdots&(p_{2}p_{K})^{1/2}r_{2K}\\\
\vdots&\vdots&\ddots&\vdots\\\
(p_{1}p_{K})^{1/2}r_{K1}&(p_{2}p_{K})^{1/2}r_{K2}&\cdots&p_{K}r_{KK}^{*}\\\
\end{bmatrix}$
$\displaystyle=\mathbf{Q}_{1}^{*}\mathbf{\Gamma}_{1}{\mathbf{Q}_{1}^{*}}^{\top},$

which shows that in the case of $K$ clusters, the eigenvalues corresponding to
the space spanned by $\mathbf{R}^{*}$ in decomposition (10) can be obtained by
considering the eigenvalues of a $K\times K$ matrix.

The second part of the space in decomposition (10) is spanned by the diagonal
centering matrices $(a_{kk}-r_{kk})\mathbf{J}_{k}$. The centering matrix is a
projector matrix that has $p_{k}-1$ eigenvalues of 1, so that
$(a_{kk}-r_{kk})\mathbf{J}_{k}$ has $p_{k}-1$ eigenvalues of $a_{kk}-r_{kk}$.

Combining both parts, we can express

$-\log|\bm{\Theta}|=-\log|\mathbf{P}^{1/2}\mathbf{R}^{*}\mathbf{P}^{1/2}|-\sum_{k=1}^{K}(p_{k}-1)\log(a_{kk}-r_{kk}).$

Second, we address the second term in objective function (3), the trace of the
matrix product $\mathbf{S\Theta}$. It is straightforward to express this term
in terms of $\mathbf{A}$ and $\mathbf{R}$, namely

$\operatorname*{tr}\mathbf{S\Theta}=\operatorname*{tr}\mathbf{SURU}^{\top}+\sum_{\ell=1}^{K}(a_{\ell\ell}-r_{\ell\ell})\operatorname*{tr}\mathbf{S}_{\ell},$

where $\mathbf{S}_{\ell}$ is the sample covariance matrix computed from the
$p_{\ell}$ variables in cluster $\ell$.

Third, the clusterpath penalty term in objective function (3) can be written
as

$\mathcal{P}(\mathbf{\Theta})=\sum_{j<j^{\prime}}w_{jj^{\prime}}d_{jj^{\prime}}(\bm{\Theta})=\sum_{k<\ell}\sum_{j\in\mathcal{C}_{k}}\sum_{j^{\prime}\in\mathcal{C}_{\ell}}w_{jj^{\prime}}d_{jj^{\prime}}(\mathbf{\Theta}),$

where $\mathcal{C}_{k}$ denotes the set of variables belonging to cluster $k$.
Because $d_{jj^{\prime}}(\bm{\Theta})$ evaluates to the same value for all
$j\in\mathcal{C}_{k}$ and $j^{\prime}\in\mathcal{C}_{\ell}$, the penalty can
be further reduced into

$\mathcal{P}(\mathbf{\Theta})=\sum_{k<\ell}\left(\sum_{j\in\mathcal{C}_{k}}\sum_{j^{\prime}\in\mathcal{C}_{\ell}}w_{jj^{\prime}}\right)d_{k\ell}(\mathbf{A},\mathbf{R})=\sum_{k<\ell}\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{\ell}d_{k\ell}(\mathbf{A},\mathbf{R}),$

where $\mathbf{W}$ is the symmetric $p\times p$ weight matrix containing the
individual weights $w_{jj^{\prime}}$, $\mathbf{u}_{k}$ denotes the $k$th
column of the membership matrix $\mathbf{U}$, and

$d_{k\ell}^{2}(\mathbf{A},\mathbf{R})=(a_{kk}-a_{\ell\ell})^{2}+(p_{k}-1)(r_{kk}-r_{k\ell})^{2}+(p_{\ell}-1)(r_{\ell\ell}-r_{k\ell})^{2}+\sum_{\begin{subarray}{c}m=1\\\
m\notin\\{k,\ell\\}\end{subarray}}^{K}p_{m}(r_{km}-r_{\ell m})^{2}.$

Note that the subscripts $k\ell$ in $d_{k\ell}(\mathbf{A},\mathbf{R})$ always
refer to the cluster formulation of the Euclidean distance.

Finally, combining all decomposed terms we obtain expression (5) for the
objective function $L(\mathbf{\Theta})$ in terms of $\mathbf{A}$ and
$\mathbf{R}$, namely

$\displaystyle L(\mathbf{\Theta})={}$ $\displaystyle L(\mathbf{A,R})$
$\displaystyle={}$
$\displaystyle-\log|\mathbf{P}^{1/2}\mathbf{R}^{*}\mathbf{P}^{1/2}|-\sum_{k=1}^{K}(p_{k}-1)\log(a_{kk}-r_{kk})+\operatorname*{tr}\mathbf{SURU}^{\top}+\sum_{\ell=1}^{K}(a_{\ell\ell}-r_{\ell\ell})\operatorname*{tr}\mathbf{S}_{\ell}$
$\displaystyle+\lambda\sum_{k<\ell}\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{\ell}d_{k\ell}(\mathbf{A},\mathbf{R}).$

##  Appendix B Derivations Cyclic Block Coordinate Descent Algorithm

We provide detailed derivations of the parameterization of the objective
function in Appendix B.1 and its gradient and Hessian in Appendix B.2.

###  B.1 Parameterization

In Section 2.2, efficient expressions have been obtained that make use of the
block structure of $\mathbf{\Theta}$. What remains to be done is to separate
$L(\mathbf{A},\mathbf{R})$ into those parts that depend on cluster $k$ and
those that do not. In the following, we make the assumption, without loss of
generality, that the variables are arranged to consistently position cluster
$k$ as the last cluster. Let the elements of $\mathbf{\Theta}$ that pertain to
cluster $k$ be separated from the others (denoted by the subscript 0), so that

$\displaystyle\mathbf{\Theta}$
$\displaystyle=\left[\begin{array}[]{ccc|c}(a_{11}-r_{11})\mathbf{I}+r_{11}\mathbf{11}^{\top}&r_{12}\mathbf{11}^{\top}&\cdots&r_{13}\mathbf{11}^{\top}\\\
r_{21}\mathbf{11}^{\top}&(a_{22}-r_{22})\mathbf{I}+r_{22}\mathbf{11}^{\top}&\cdots&r_{2k}\mathbf{11}^{\top}\\\
\vdots&\vdots&\ddots&\vdots\\\ \hline\cr
r_{k1}\mathbf{11}^{\top}&r_{k2}\mathbf{11}^{\top}&\cdots&(a_{kk}-r_{kk})\mathbf{I}+r_{kk}\mathbf{11}^{\top}\\\
\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{c|c}\mathbf{\Theta}_{0}&\mathbf{\Theta}_{0k}\\\
\hline\cr\mathbf{\Theta}_{0k}^{\top}&\mathbf{\Theta}_{kk}\end{array}\right].$

Consequently, updating all parameters in $\mathbf{A}$ and $\mathbf{R}$ that
are associated with cluster $k$ updates an entire block of variables in
$\mathbf{\Theta}$. To separate the term
$-\log|\mathbf{P}^{1/2}\mathbf{R}^{*}\mathbf{P}^{1/2}|$ in objective function
(5), it helps to write $\mathbf{P}^{1/2}\mathbf{R}^{*}\mathbf{P}^{1/2}$ as the
$2\times 2$ block matrix

$\mathbf{P}^{1/2}\mathbf{R}^{*}\mathbf{P}^{1/2}=\left[\begin{array}[]{c|c}\mathbf{P}_{0}^{1/2}\mathbf{R}_{0}^{*}\mathbf{P}_{0}^{1/2}&p_{k}^{1/2}\mathbf{P}_{0}^{1/2}\mathbf{r}_{k}\\\
\hline\cr
p_{k}^{1/2}(\mathbf{P}_{0}^{1/2}\mathbf{r}_{k})^{\top}&p_{k}r_{kk}^{*}\\\
\end{array}\right]$

with $\mathbf{r}_{k}$ the $K-1$ vector containing all elements in the
$k^{\text{th}}$ column of $\mathbf{R}$ apart from $r_{kk}$, and

$\mathbf{P}=\left[\begin{array}[]{cc}\mathbf{P}_{0}&\mathbf{0}\\\
\mathbf{0}^{\top}&p_{k}\\\ \end{array}\right].$

From linear algebra, we have the following expression for the log determinant
of the $2\times 2$ block matrix

$-\log|\mathbf{P}^{1/2}\mathbf{R}^{*}\mathbf{P}^{1/2}|=-\log|\mathbf{P}_{0}^{1/2}\mathbf{R}_{0}^{*}\mathbf{P}_{0}^{1/2}|-\log(a_{kk}+(p_{k}-1)r_{kk}-p_{k}\mathbf{r}_{k}^{\top}(\mathbf{R}_{0}^{*})^{-1}\mathbf{r}_{k}).$

This allows us to rewrite the negative log determinant of $\mathbf{\Theta}$ as

$\displaystyle-\log|\mathbf{\Theta}|={}$
$\displaystyle-\log|\mathbf{P}_{0}^{1/2}\mathbf{R}_{0}^{*}\mathbf{P}_{0}^{1/2}|-\log(a_{kk}+(p_{k}-1)r_{kk}-p_{k}\mathbf{r}_{k}^{\top}(\mathbf{R}_{0}^{*})^{-1}\mathbf{r}_{k})$
$\displaystyle-(p_{k}-1)\log(a_{kk}-r_{kk})-\sum_{\begin{subarray}{c}\ell=1\\\
\ell\neq k\end{subarray}}^{K}(p_{\ell}-1)\log(a_{\ell\ell}-r_{\ell\ell}).$
(11)

Using the re-expression from objective function (5), the trace of
$\mathbf{S}\mathbf{\Theta}$ separates readily into

$\operatorname*{tr}\mathbf{S\Theta}=2\mathbf{u}_{k}^{\top}\mathbf{S}\mathbf{U}_{0}\mathbf{r}_{k}+\mathbf{u}_{k}^{\top}\mathbf{S}\mathbf{u}_{k}r_{kk}+(a_{kk}-r_{kk})\operatorname*{tr}\mathbf{S}_{k}+\operatorname*{tr}\mathbf{S}\mathbf{U}_{0}\mathbf{R}_{0}\mathbf{U}_{0}^{\top}+\sum_{\begin{subarray}{c}\ell=1\\\
\ell\neq
k\end{subarray}}^{K}(a_{\ell\ell}-r_{\ell\ell})\operatorname*{tr}\mathbf{S}_{\ell}.$

It can be verified that the last two terms are constant in $a_{kk}$,
$\mathbf{r}_{k}$, and $r_{kk}$, and the remaining terms are linear in these
unknowns. Finally, due to the symmetry of $\mathbf{R}$, the clusterpath
penalty is not separable into parts dependent and independent of cluster $k$.

Combining the foregoing results, $L(\mathbf{A},\mathbf{R})$ can be expressed
in terms of parts $a_{kk}$, $\mathbf{r}_{k}$, and $r_{kk}$, and we obtain
$L(a_{kk},\mathbf{r}_{k},r_{kk})$ in equation (6) as

$\displaystyle L(a_{kk},\mathbf{r}_{k},r_{kk})={}$ $\displaystyle
L_{\text{det}}(a_{kk},\mathbf{r}_{k},r_{kk})+L_{\text{cov}}(a_{kk},\mathbf{r}_{k},r_{kk})+L_{\text{cpath}}(a_{kk},\mathbf{r}_{k},r_{kk})+C$
$\displaystyle={}$
$\displaystyle-\log(a_{kk}+(p_{k}-1)r_{kk}-p_{k}\mathbf{r}_{k}^{\top}(\mathbf{R}_{0}^{*})^{-1}\mathbf{r}_{k})-(p_{k}-1)\log(a_{kk}-r_{kk})$
$\displaystyle+2\mathbf{u}_{k}^{\top}\mathbf{S}\mathbf{U}_{0}\mathbf{r}_{k}+\mathbf{u}_{k}^{\top}\mathbf{S}\mathbf{u}_{k}r_{kk}+(a_{kk}-r_{kk})\operatorname*{tr}\mathbf{S}_{k}$
$\displaystyle+\lambda\left({\sum_{k<\ell}\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{\ell}d_{k\ell}(\mathbf{A},\mathbf{R})}\right)+~{}C$

with

$\displaystyle C={}$
$\displaystyle-\log|\mathbf{P}_{0}^{1/2}\mathbf{R}_{0}^{*}\mathbf{P}_{0}^{1/2}|-\sum_{\begin{subarray}{c}\ell=1\\\
\ell\neq k\end{subarray}}^{K}(p_{\ell}-1)\log(a_{\ell\ell}-r_{\ell\ell})$
$\displaystyle+\operatorname*{tr}\mathbf{S}\mathbf{U}_{0}\mathbf{R}_{0}\mathbf{U}_{0}^{\top}+\sum_{\begin{subarray}{c}\ell=1\\\
\ell\neq
k\end{subarray}}^{K}(a_{\ell\ell}-r_{\ell\ell})\operatorname*{tr}\mathbf{S}_{\ell}.$

###  B.2 Derivations of the Gradient and Hessian

In the derivations that follow, we obtain expressions for the gradient and
Hessian of (6). To shorten notation, let
$\mathbf{V}=(\mathbf{R}_{0}^{*})^{-1}$ and

$h(a_{kk},\mathbf{r}_{k},r_{kk})=a_{kk}+(p_{k}-1)r_{kk}-p_{k}\mathbf{r}_{k}^{\top}\mathbf{V}\mathbf{r}_{k}.$

For convenience, the elements $r_{km}$ of the vector $\mathbf{r}_{k}$ are
treated individually. Consequently, we should note that $1\leq m\leq K$,
$m\neq k$. For the gradient, we obtain

$\displaystyle\frac{\partial
L_{\text{det}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial a_{kk}}={}$
$\displaystyle-\frac{1}{h(a_{kk},\mathbf{r}_{k},r_{kk})}-\frac{p_{k}-1}{a_{kk}-r_{kk}}$
$\displaystyle\frac{\partial
L_{\text{det}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial r_{km}}={}$
$\displaystyle\frac{2p_{k}}{h(a_{kk},\mathbf{r}_{k},r_{kk})}\mathbf{v}_{m}^{\top}\mathbf{r}_{k}$
$\displaystyle\frac{\partial
L_{\text{det}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial r_{kk}}={}$
$\displaystyle-\frac{p_{k}-1}{h(a_{kk},\mathbf{r}_{k},r_{kk})}+\frac{p_{k}-1}{a_{kk}-r_{kk}}$
$\displaystyle\frac{\partial
L_{\text{cov}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial a_{kk}}={}$
$\displaystyle\operatorname*{tr}\mathbf{S}_{k}$ $\displaystyle\frac{\partial
L_{\text{cov}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial r_{km}}={}$
$\displaystyle 2\mathbf{u}_{m}^{\top}\mathbf{S}\mathbf{u}_{k}$
$\displaystyle\frac{\partial
L_{\text{cov}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial r_{kk}}={}$
$\displaystyle\mathbf{u}_{k}^{\top}\mathbf{S}\mathbf{u}_{k}-\operatorname*{tr}\mathbf{S}_{k}$
$\displaystyle\frac{\partial
L_{\text{cpath}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial a_{kk}}={}$
$\displaystyle\lambda\sum_{\begin{subarray}{c}\ell=1\\\ \ell\neq
k\end{subarray}}^{K}\frac{\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{\ell}}{d_{k\ell}(\mathbf{A},\mathbf{R})}(a_{kk}-a_{\ell\ell})$
$\displaystyle\frac{\partial
L_{\text{cpath}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial r_{km}}={}$
$\displaystyle\lambda\sum_{\begin{subarray}{c}\ell=1\\\
\ell\notin\\{k,m\\}\end{subarray}}^{K}\left(\frac{\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{\ell}}{d_{k\ell}(\mathbf{A},\mathbf{R})}p_{m}(r_{km}-r_{m\ell})+\frac{\mathbf{u}_{m}^{\top}\mathbf{W}\mathbf{u}_{\ell}}{d_{m\ell}(\mathbf{A},\mathbf{R})}p_{k}(r_{km}-r_{k\ell})\right)$
$\displaystyle+\lambda\frac{\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{m}}{d_{km}(\mathbf{A},\mathbf{R})}((p_{k}-1)(r_{km}-r_{kk})+(p_{m}-1)(r_{km}-r_{mm}))$
$\displaystyle\frac{\partial
L_{\text{cpath}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial r_{kk}}={}$
$\displaystyle\lambda(p_{k}-1)\sum_{\begin{subarray}{c}\ell=1\\\ \ell\neq
k\end{subarray}}^{K}\frac{\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{\ell}}{d_{k\ell}(\mathbf{A},\mathbf{R})}(r_{kk}-r_{k\ell}).$

In the derivations for the Hessian, we require an additional index
$m^{\prime}$ to define the off-diagonal elements in the block
$\partial^{2}L(a_{kk},\mathbf{r}_{k},r_{kk})/\partial\mathbf{r}_{k}^{2}$ that
satisfies $1\leq m^{\prime}\leq K$, $m^{\prime}\notin\\{k,m\\}$. For the
Hessian, we obtain

$\displaystyle\frac{\partial^{2}L_{\text{det}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial
a_{kk}^{2}}={}$
$\displaystyle\frac{1}{h^{2}(a_{kk},\mathbf{r}_{k},r_{kk})}+\frac{p_{k}-1}{(a_{kk}-r_{kk})^{2}}$
$\displaystyle\frac{\partial^{2}L_{\text{det}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial
a_{kk}\partial r_{km}}={}$
$\displaystyle-\frac{2p_{k}}{h^{2}(a_{kk},\mathbf{r}_{k},r_{kk})}\mathbf{v}_{m}^{\top}\mathbf{r}_{k}$
$\displaystyle\frac{\partial^{2}L_{\text{det}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial
a_{kk}\partial r_{kk}}={}$
$\displaystyle\frac{p_{k}-1}{h^{2}(a_{kk},\mathbf{r}_{k},r_{kk})}-\frac{p_{k}-1}{(a_{kk}-r_{kk})^{2}}$
$\displaystyle\frac{\partial^{2}L_{\text{det}}(a_{kk},\mathbf{r}_{k},r_{km})}{\partial
r_{km}^{2}}={}$
$\displaystyle\frac{2p_{k}}{h(a_{kk},\mathbf{r}_{k},r_{kk})}v_{mm}+\frac{4p_{k}^{2}}{h^{2}(a_{kk},\mathbf{r}_{k},r_{kk})}\mathbf{v}_{m}^{\top}\mathbf{r}_{k}\mathbf{r}_{k}^{\top}\mathbf{v}_{m}$
$\displaystyle\frac{\partial^{2}L_{\text{det}}(a_{kk},\mathbf{r}_{k},r_{km})}{\partial
r_{km}\partial r_{km^{\prime}}}={}$
$\displaystyle\frac{2p_{k}}{h(a_{kk},\mathbf{r}_{k},r_{kk})}v_{mm^{\prime}}+\frac{4p_{k}^{2}}{h^{2}(a_{kk},\mathbf{r}_{k},r_{kk})}\mathbf{v}_{m}^{\top}\mathbf{r}_{k}\mathbf{r}_{k}^{\top}\mathbf{v}_{m^{\prime}}$
$\displaystyle\frac{\partial^{2}L_{\text{det}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial
r_{km}\partial r_{kk}}={}$
$\displaystyle-\frac{2p_{k}(p_{k}-1)}{h^{2}(a_{kk},\mathbf{r}_{k},r_{kk})}\mathbf{v}_{m}^{\top}\mathbf{r}_{k}$
$\displaystyle\frac{\partial^{2}L_{\text{det}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial
r_{kk}^{2}}={}$
$\displaystyle\frac{(p_{k}-1)^{2}}{h^{2}(a_{kk},\mathbf{r}_{k},r_{kk})}+\frac{p_{k}-1}{(a_{kk}-r_{kk})^{2}}$
$\displaystyle\frac{\partial^{2}L_{\text{cpath}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial
a_{kk}^{2}}={}$ $\displaystyle\lambda\sum_{\begin{subarray}{c}\ell=1\\\
\ell\neq
k\end{subarray}}^{K}\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{\ell}\left(\frac{1}{d_{k\ell}(\mathbf{A},\mathbf{R})}-\frac{(a_{kk}-a_{\ell\ell})^{2}}{d^{3}_{k\ell}(\mathbf{A},\mathbf{R})}\right)$
$\displaystyle\frac{\partial^{2}L_{\text{cpath}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial
a_{kk}\partial r_{km}}={}$
$\displaystyle-\lambda\sum_{\begin{subarray}{c}\ell=1\\\
\ell\notin\\{k,m\\}\end{subarray}}^{K}\frac{\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{\ell}}{d^{3}_{k\ell}(\mathbf{A},\mathbf{R})}p_{m}(a_{kk}-a_{\ell\ell})(r_{km}-r_{m\ell})$
$\displaystyle-\lambda\frac{\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{m}}{d^{3}_{km}(\mathbf{A},\mathbf{R})}(a_{kk}-a_{mm})((p_{k}-1)(r_{km}-r_{kk})+(p_{m}-1)(r_{km}-r_{mm}))$
$\displaystyle\frac{\partial^{2}L_{\text{cpath}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial
a_{kk}\partial r_{kk}}={}$
$\displaystyle-\lambda(p_{k}-1)\sum_{\begin{subarray}{c}\ell=1\\\ \ell\neq
k\end{subarray}}^{K}\frac{\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{\ell}}{d^{3}_{k\ell}(\mathbf{A},\mathbf{R})}(a_{kk}-a_{\ell\ell})(r_{kk}-r_{k\ell})$
$\displaystyle\frac{\partial^{2}L_{\text{cpath}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial
r_{km}^{2}}={}$ $\displaystyle\lambda\sum_{\begin{subarray}{c}\ell=1\\\
\ell\notin\\{k,m\\}\end{subarray}}^{K}p_{m}\frac{\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{\ell}}{d_{k\ell}(\mathbf{A},\mathbf{R})}\left(1-\frac{p_{m}(r_{km}-r_{m\ell})^{2}}{d^{2}_{k\ell}(\mathbf{A},\mathbf{R})}\right)$
$\displaystyle+\lambda\sum_{\begin{subarray}{c}\ell=1\\\
\ell\notin\\{k,m\\}\end{subarray}}^{K}p_{k}\frac{\mathbf{u}_{m}^{\top}\mathbf{W}\mathbf{u}_{\ell}}{d_{m\ell}(\mathbf{A},\mathbf{R})}\left(1-\frac{p_{k}(r_{km}-r_{k\ell})^{2}}{d^{2}_{m\ell}(\mathbf{A},\mathbf{R})}\right)$
$\displaystyle+\lambda\frac{\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{m}}{d_{km}(\mathbf{A},\mathbf{R})}\left(p_{k}+p_{m}-2-\frac{((p_{k}-1)(r_{km}-r_{kk})+(p_{m}-1)(r_{km}-r_{mm}))^{2}}{d^{2}_{km}(\mathbf{A},\mathbf{R})}\right)$
$\displaystyle\frac{\partial^{2}L_{\text{cpath}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial
r_{km}\partial r_{km^{\prime}}}={}$
$\displaystyle-\lambda\sum_{\begin{subarray}{c}\ell=1\\\
\ell\notin\\{k,m,m^{\prime}\\}\end{subarray}}^{K}\frac{\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{\ell}}{d^{3}_{k\ell}(\mathbf{A},\mathbf{R})}p_{m}p_{m^{\prime}}(r_{km}-r_{m\ell})(r_{km^{\prime}}-r_{m^{\prime}\ell})$
$\displaystyle-\lambda\frac{\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{m^{\prime}}}{d^{3}_{km^{\prime}}(\mathbf{A},\mathbf{R})}((p_{k}-1)(r_{km^{\prime}}-r_{kk})+(p_{m}-1)(r_{km^{\prime}}-r_{m^{\prime}m^{\prime}}))p_{m}(r_{km}-r_{mm^{\prime}})$
$\displaystyle+\lambda
p_{k}\frac{\mathbf{u}_{m}^{\top}\mathbf{W}\mathbf{u}_{m^{\prime}}}{d_{mm^{\prime}}(\mathbf{A},\mathbf{R})}\left(1-\frac{p_{k}(r_{km^{\prime}}-r_{km})^{2}}{d^{2}_{mm^{\prime}}(\mathbf{A},\mathbf{R})}\right)$
$\displaystyle-\lambda\frac{\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{m}}{d_{km}(\mathbf{A},\mathbf{R})}p_{m^{\prime}}(r_{km^{\prime}}-r_{mm^{\prime}})((p_{k}-1)(r_{km}-r_{kk})+(p_{m}-1)(r_{km}-r_{mm}))$
$\displaystyle\frac{\partial^{2}L_{\text{cpath}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial
r_{km}\partial r_{kk}}={}$
$\displaystyle-\lambda\sum_{\begin{subarray}{c}\ell=1\\\
\ell\notin\\{k,m\\}\end{subarray}}^{K}\frac{\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{\ell}}{d^{3}_{k\ell}(\mathbf{A},\mathbf{R})}p_{m}(p_{k}-1)(r_{kk}-r_{k\ell})(r_{km}-r_{m\ell})$
$\displaystyle-\lambda(p_{k}-1)\frac{\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{m}}{d_{km}(\mathbf{A},\mathbf{R})}\left(1-\frac{(p_{k}-1)(r_{km}-r_{kk})}{d^{2}_{km}(\mathbf{A},\mathbf{R})}(r_{km}-r_{kk})\right.$
$\displaystyle+\left.\frac{(p_{m}-1)(r_{km}-r_{mm})}{d^{2}_{km}(\mathbf{A},\mathbf{R})}(r_{km}-r_{kk})\right)$
$\displaystyle\frac{\partial^{2}L_{\text{cpath}}(a_{kk},\mathbf{r}_{k},r_{kk})}{\partial
r_{kk}^{2}}={}$
$\displaystyle\lambda(p_{k}-1)\sum_{\begin{subarray}{c}\ell=1\\\ \ell\neq
k\end{subarray}}^{K}\mathbf{u}_{k}^{\top}\mathbf{W}\mathbf{u}_{\ell}\left(\frac{1}{d_{k\ell}(\mathbf{A},\mathbf{R})}-\frac{(p_{k}-1)(r_{kk}-r_{k\ell})^{2}}{d^{3}_{k\ell}(\mathbf{A},\mathbf{R})}\right),$

where the terms relating to $L_{\text{cov}}(a_{kk},\mathbf{r}_{k},r_{kk})$ are
left out as they evaluate to zero. It should be noted that if the size of the
$k$th cluster is one, the variable $r_{kk}$ is effectively not a part of the
objective function. Consequently, the $k$th element of the gradient is zero
and the $k$th row and column of the Hessian are also filled with zeros.

##  Appendix C Additional Simulation Results

This appendix provides supplementary simulation results from the experiments
discussed in Section 3. Figure 10 illustrates the results for the chain
simulation design with varying numbers of variables. Figure 11 shows the
results for designs featuring an approximate block structure. The performance
of the methods is evaluated based on estimation accuracy (Frobenius norm) and
clustering quality (number of clusters and ARI).

Figure 10: Results for the chain simulation design with an increasing number
of variables (in separate columns). The top row contains boxplots of the
Frobenius norm together with black diamonds representing the average. In the
middle and bottom rows, diamonds display the average of the estimated number
of clusters and the average ARI, respectively. In addition, the size of the
grey dots represent the frequency of different values across the replications.
Reference lines are drawn for the true number of clusters and the maximum
possible ARI value. Aggregation performance is omitted for $\mathbf{S}^{-1}$
as it does not perform any aggregation.  Figure 11: Results for the simulation
designs with an approximate block structure (in separate columns). The top row
contains boxplots of the Frobenius norm together with black diamonds
representing the average. In the middle and bottom rows, diamonds display the
average of the estimated number of clusters and the average ARI, respectively.
In addition, the size of the grey dots represent the frequency of different
values across the replications. Reference lines are drawn for the true number
of clusters and the maximum possible ARI value. Aggregation performance is
omitted for $\mathbf{S}^{-1}$ as it does not perform any aggregation. In the
unstructured design, a misspecified tree for TAGL does not exist since any
tree hierarchy contains the true clustering (each variable being its own
cluster).

##  Appendix D Humor Styles Questionnaire

This appendix contains information on the HSQ developed by Martin . (2003),
which is analyzed in Section 5.3. The $p=32$ items of the HSQ, grouped by the
humor style being measured, are listed in Table LABEL:tab:HSQ together with
their position in the survey.

Humor style | Items
---|---
Affiliative | 1. |  I usually don’t laugh or joke around much with other people.∗
| 5. |  I don’t have to work very hard at making other people laugh—I seem to be a naturally humorous person.
| 9. |  I rarely make other people laugh by telling funny stories about myself.∗
| 13. |  I laugh and joke a lot with my closest friends.
| 17. |  I usually don’t like to tell jokes or amuse people.∗
| 21. |  I enjoy making people laugh.
| 25. |  I don’t often joke around with my friends.∗
| 29. |  I usually can’t think of witty things to say when I’m with other people.∗
Self-enhancing | 2. |  If I am feeling depressed, I can usually cheer myself up with humor.
| 6. |  Even when I’m by myself, I’m often amused by the absurdities of life.
| 10. |  If I am feeling upset or unhappy I usually try to think of something funny about the situation to make myself feel better.
| 14. |  My humorous outlook on life keeps me from getting overly upset or depressed about things.
| 18. |  If I’m by myself and I’m feeling unhappy, I make an effort to think of something funny to cheer myself up.
| 22. |  If I am feeling sad or upset, I usually lose my sense of humor.∗
| 26. |  It is my experience that thinking about some amusing aspect of a situation is often a very effective way of coping with problems.
| 30. |  I don’t need to be with other people to feel amused—I can usually find things to laugh about even when I’m by myself.
Aggressive | 3. |  If someone makes a mistake, I will often tease them about it.
| 7. |  People are never offended or hurt by my sense of humor.∗
| 11. |  When telling jokes or saying funny things, I am usually not very concerned about how other people are taking it.
| 15. |  I do not like it when people use humor as a way of criticizing or putting someone down.∗
| 19. |  Sometimes I think of something that is so funny that I can’t stop myself from saying it, even if it is not appropriate for the situation.
| 23. |  I never participate in laughing at others even if all my friends are doing it.∗
| 27. |  If I don’t like someone, I often use humor or teasing to put them down.
| 31. |  Even if something is really funny to me, I will not laugh or joke about it if someone will be offended.∗
Self-defeating | 4. |  I let people laugh at me or make fun at my expense more than I should.
| 8. |  I will often get carried away in putting myself down if it makes my family or friends laugh.
| 12. |  I often try to make people like or accept me more by saying something funny about my own weaknesses, blunders, or faults.
| 16. |  I don’t often say funny things to put myself down.∗
| 20. |  I often go overboard in putting myself down when I am making jokes or trying to be funny.
| 24. |  When I am with friends or family, I often seem to be the one that other people make fun of or joke about.
| 28. |  If I am having problems or feeling unhappy, I often cover it up by joking around, so that even my closest friends don’t know how I really feel.
| 32. |  Letting others laugh at me is my way of keeping my friends and family in good spirits.
Table 2: List of the $p=32$ items of the humor styles questionnaire of Martin
. (2003), grouped by the humor style being measured and reported together with
their position in the survey. Items marked with an asterisk are reverse-keyed,
i.e., the responses on the five-point rating scale are reversed prior to
analysis.