# Learning the hub graphical Lasso model with the structured sparsity via an
efficient algorithm

Chengjing Wang School of Mathematics, Southwest Jiaotong University, No.999,
Xian Road, West Park, High-tech Zone, Chengdu 611756, China (Email:
renascencewang@hotmail.com).  Peipei Tang Corresponding author, School of
Computer and Computing Science, Hangzhou City University, Hangzhou 310015,
China (Email: tangpp@hzcu.edu.cn).  Wenling He School of Mathematics,
Southwest Jiaotong University, No.999, Xian Road, West Park, High-tech Zone,
Chengdu 611756, China (Email: 1165649932@qq.com).  Meixia Lin Engineering
Systems and Design, Singapore University of Technology and Design, Singapore.
This author’s research is supported by The Singapore University of Technology
and Design under MOE Tier 1 Grant SKI 2021_02_08.
(Email:meixia_lin@sutd.edu.sg).

###### Abstract

Graphical models have exhibited their performance in numerous tasks ranging
from biological analysis to recommender systems. However, graphical models
with hub nodes are computationally difficult to fit, particularly when the
dimension of the data is large. To efficiently estimate the hub graphical
models, we introduce a two-phase algorithm. The proposed algorithm first
generates a good initial point via a dual alternating direction method of
multipliers (ADMM), and then warm starts a semismooth Newton (SSN) based
augmented Lagrangian method (ALM) to compute a solution that is accurate
enough for practical tasks. The sparsity structure of the generalized Jacobian
ensures that the algorithm can obtain a nice solution very efficiently.
Comprehensive experiments on both synthetic data and real data show that it
obviously outperforms the existing state-of-the-art algorithms. In particular,
in some high dimensional tasks, it can save more than 70% of the execution
time, meanwhile still achieves a high-quality estimation.

Keywords:  hub graphical Lasso, structured sparsity, augmented Lagrangian
method, semismooth Newton method

AMS subject classification: 65K05, 90C06, 90C25, 90C90

##  1 Introduction

Graphical models have been widely adopted to characterize relationships of
data variables in various applications, such as gene finding, medical
diagnosis, modeling of protein structures, and so on. To encode conditional
dependence relationships among a set of variables, the Gaussian graphical
model is employed to estimate the precision matrix $\Theta=\Sigma^{-1}$ of the
data distribution, where $\Sigma$ is the corresponding covariance matrix, and
then the sparsity pattern of $\Theta$ is used to determine the conditional
independence graph [8]. In particular, each variable in the Gaussian graphical
model represents a node in a graph, and an edge will link two nodes if the
corresponding entry of the precision matrix $\Theta$ is non-zero. A series of
variants of Gaussian graphical model have been designed to satisfy the
different requirements of practical tasks, one may refer to [43], [31], [20],
[10], [27], [21], [38], [28], [36], and [6], to name but a few.

Among these studies, one interesting problem is to construct a graph with
_hubs_ , i.e., learning a graphical model that contains a few special nodes,
which are _densely-connected_ to several other nodes. Such kind of graph is
universal in the real world. For instance, in the social network of Twitter,
only a few users (e.g., famous actors) have a large number of followers and
their tweets have strong influence, while most users do not [15]. A more
typical example is the gene-regulatory networks, which are generally thought
to approximate a hierarchical scale-free network topology [3], i.e., the
number of edges follows a power-law distribution. In addition, one can refer
to [1], [18], and [14] for more examples.

To model a graph with hubs, Tan et al. [37] proposed to estimate the precision
matrix via adding a sparse-group Lasso-type regularization to the Gaussian
negative log-likelihood. Specifically, let $S$ be an empirical covariance
matrix with its $j$-th column as $S_{j}$, denote the set of $p\times p$ square
matrices by $\mathbb{M}^{p}$ and the set of $p\times p$ symmetric matrices by
$\mathbb{S}^{p}$. The problem of learning a graph with hub nodes can be
formulated as the following _hub graphical Lasso_ (HGL) optimization problem:

${\min_{\Theta\in\mathbb{S}^{p}}}\left\\{-\log\det(\Theta)+\langle
S,\Theta\rangle+P(\Theta)\colon\Theta\succ 0\right\\},$ (1)

where $\langle S,\Theta\rangle$ is the standard matrix inner product between
$S$ and $\Theta$, $\Theta\succ 0$ means that $\Theta\in\mathbb{S}^{p}$ is
positive definite. The function $P(\Theta)$ is referred to as the _hub
penalty_ , which takes the form

$\displaystyle P(\Theta)=\min_{Z\in\mathbb{S}^{p},V\in\mathbb{M}^{p}}$
$\displaystyle\Big{\\{}\lambda_{1}\|Z-{\rm diag}(Z)\|_{1}+\lambda_{2}\|V-{\rm
diag}(V)\|_{1}$ $\displaystyle+\lambda_{3}\sum_{j=1}^{p}\|{\left(V-{\rm
diag}(V)\right)}_{j}\|\,\Big{|}\,\Theta=Z+V+V^{T}\Big{\\}}.$

For any $x\in\mathbb{R}^{p}$, we denote its q-norm as
$\|x\|_{q}=(\sum_{i=1}^{p}|x_{i}|^{q})^{1/q}$. For simplicity, we denote
$\|\cdot\|=\|\cdot\|_{2}$. For any $X\in\mathbb{M}^{p}$, we denote
$\|X\|_{1}=\sum_{i=1}^{p}\sum_{j=1}^{p}|X_{ij}|$. This penalty function
attempts to decompose the estimated precision matrix $\Theta$ into three
matrices $Z$, $V$ and $V^{T}$. The non-zero entries of $Z$ indicate the edges
between the non-hub nodes, and the non-zero columns of $V$ corresponds to the
hub nodes. The $\ell_{1}$ penalty on the off-diagonal entries of $Z$ promotes
the sparsity of the solution, which reduces the links between the non-hub
nodes. The combination of the $\ell_{1}$ and $\ell_{2}$ penalties for the
columns of $V$ induces group sparsity, which enforces the entries of each
column to be either all zeros or almost entirely nonzeros.

However, in the HGL model (1), it is assumed that there are no priors provided
for the nodes. That is, there is no _priori_ of the hubs. To further make use
of the prior knowledge of the hubs that may be available in advance, Bai and
Zhou [25] improved the HGL via dividing the hub penalty in a discriminatory
way:

$\displaystyle\widetilde{P}(\Theta)=\min_{Z\in\mathbb{S}^{p},V\in\mathbb{M}^{p}}\Big{\\{}\lambda_{1}\|Z-{\rm
diag}(Z)\|_{1}+\lambda_{2}\sum_{j\notin\mathcal{D}}\|{(V-{\rm
diag}(V))}_{j}\|_{1}$
$\displaystyle\qquad\qquad+\lambda_{3}\sum_{j\notin\mathcal{D}}\|{(V-{\rm
diag}(V))}_{j}\|+\lambda_{4}\sum_{j\in\mathcal{D}}\|{(V-{\rm
diag}(V))}_{j}\|_{1}$
$\displaystyle\qquad\qquad+\lambda_{5}\sum_{j\in\mathcal{D}}\|{(V-{\rm
diag}(V))}_{j}\|\,\Big{|}\,\Theta=Z+V+V^{T}\Big{\\}},$

where $\mathcal{D}\subseteq\\{1,\dots,p\\}$ is the index set of hubs with
prior information. Here we set $\lambda_{2}\geq\lambda_{4}$ and
$\lambda_{3}\geq\lambda_{5}$ such that the known hubs are imposed a smaller
penalty. It should be noted that the _discriminated hub graphical Lasso_
(DHGL) reduces to the HGL if $\mathcal{D}=\emptyset$, and reduces to the
classical _graphical Lasso_ (GL) as
$\lambda_{2},\lambda_{3},\lambda_{4},\lambda_{5}\to\infty$ [43].

Now we directly focus on the DHGL optimization problem. For convenience, we
use $Q(Z)$ and $R(V)$ to denote the penalty function of $Z$ and $V$, i.e.,

$Q(Z)=\lambda_{1}\|Z-{\rm diag}(Z)\|_{1},\quad\forall\,Z\in\mathbb{S}^{p},$

and

$\displaystyle R(V)$ $\displaystyle=$
$\displaystyle\lambda_{2}\sum_{j\not\in\mathcal{D}}\|{(V-{\rm
diag}(V))}_{j}\|_{1}+\lambda_{3}\sum_{j\not\in\mathcal{D}}\|{(V-{\rm
diag}(V))}_{j}\|$ $\displaystyle+\lambda_{4}\sum_{j\in\mathcal{D}}\|{(V-{\rm
diag}(V))}_{j}\|_{1}+\lambda_{5}\sum_{j\in\mathcal{D}}\|{(V-{\rm
diag}(V))}_{j}\|,\quad\forall\,V\in\mathbb{M}^{p}.$

In addition, we define that $\log 0:=-\infty$. Then, the optimization problem
can be written as

${\min_{\Theta,Z\in\mathbb{S}^{p},V\in\mathbb{M}^{p}}}\left\\{\langle
S,\Theta\rangle-\log\det(\Theta)+Q(Z)+R(V)\,\Big{|}\,\Theta=Z+\mathcal{T}V,\Theta\succeq
0\right\\},$ (P)

where $\Theta\succeq 0$ indicates that $\Theta\in\mathbb{S}^{p}$ is positive
semidefinite, and $\mathcal{T}\colon\mathbb{M}^{p}\to\mathbb{S}^{p}$ is a
linear operator which is defined as $\mathcal{T}V:=V+V^{T}$, for any
$V\in\mathbb{M}^{p}$.

To solve the optimization problem HGL or DHGL, Tan et al. [37] applied the
_alternating direction method of multipliers_ (ADMM) to the problem (1), and
Bai and Zhou [25] applied the ADMM to solve the problem (P). (For short, we
call this kind of framework as pADMM.) Unfortunately, as the scale of the
problems becomes large, the efficiency of the pADMM decreases. On one hand,
the ADMM often cannot achieve a solution with high accuracy, which may ruin
the performance of the HGL or DHGL model. On the other hand, the computational
cost of the ADMM is extremely high and nearly cannot be acceptable when the
dimension $p$ of the data is very high. Hence, an efficient algorithm to solve
the problem (P) accurately is needed. For the classical graphical Lasso
problem (where $P(\Theta)=\lambda\|\Theta-{\rm diag}{\Theta}\|_{1}$), a number
of optimization algorithms were adopted to solve it, e.g., the ADMM [35],
Nesterov’s smooth gradient method [9, 26], the block coordinate ascent method
[2, 13], the interior point method [43, 24], the Newton method [30, 22], and
the dual spectral projected gradient method [29]. Besides, many efficient
algorithms have been presented towards other forms of penalty functions. For
example, Hsieh et al. [23] employed an iterative quadratic approximation to
estimate the precision matrix of the reweighted graphical Lasso model. And the
similar idea was also taken by [41] when considering the fused multiple
graphical Lasso problem. In addition, Bybee and Atchadé [5] introduced an
approximate majorize-minimize algorithm to compute the Gaussian graphical
model with change points. Tarzanagh and Michailidis [38] applied a linearized
multi-block ADMM to learn the Gaussian graphical model with structured sparse
components. More recently, Zhang et al. [45] proposed a proximal point
algorithm to solve the group graphical Lasso problem with the inner subproblem
solved by the dual-based semismooth Newton method. Nevertheless, none of these
well-designed algorithms can be directly used to solve the DHGL problem.
Moreover, in order to solve this kind of large-scale problems, the sparsity
structure should be further explored delicately.

In this paper, we develop a two-phase algorithm to compute a high-quality
solution of the DHGL problem efficiently. Specifically, in the first phase, we
apply the ADMM to the dual of the problem (P) (for short, we call it dADMM),
to generate a good initial point to warm start the second phase. In the second
phase, we apply the augmented Lagrangian method (ALM) with the subproblems
solved by the semismooth Newton (SSN) method to obtain a highly accurate
solution. When using the SSN method, we need to thoroughly exploit the
sparsity structure of the generalized Jacobian to save the computational cost.
The systematic empirical evaluations show that the algorithm is promising and
it outperforms the existing pADMM by a large margin in all data settings.

In summary, the contributions of this paper are twofold:

  * (1)

We design a highly efficient two-phase algorithm to solve the discriminated
hub graphical Lasso problem. In particular, we design a surrogate generalized
Jacobian of the proximal mapping of the hub penalty, which paves the way for
the SSN method to efficiently solve the inner subproblems of the ALM.

  * (2)

We analyze the convergence of the dADMM, the ALM, and the SSN method. We
further conduct systematic experiments to evaluate the proposed algorithm
empirically, and the experimental results on both synthetic and real data
confirm the general efficacy and the efficiency of our approach.

The rest of this paper is organized as follows. In Section 2, we present the
dual of the problem (P) and some necessary preliminaries. In Section 3, we
devote to present the two-phase algorithm. In Section 4, we describe the
numerical issues for solving the involved subproblems. In Section 5, we
implement the numerical experiments to demonstrate the efficacy and efficiency
of the proposed algorithm. In Section 6, we conclude this paper.

###  1.1 Additional notations

Let $\mathcal{X}$, $\mathcal{Y}$ and $\mathcal{Z}$ be real finite dimensional
Euclidean spaces. We denote $x\circ y$ as the Hadamard product of
$x,y\in\mathcal{X}$. Let $f:\mathcal{X}\rightarrow\mathcal{Y}$ and
$g:\mathcal{Y}\rightarrow\mathcal{Z}$, then $g\circ f$ denotes the composition
of the two functions. Given $r>0$, we denote
$\mathcal{B}_{2}^{r}:=\\{x\in\mathbb{R}^{p}\,|\,\|x\|\leq r\\}$ and
$\mathcal{B}_{\infty}^{r}:=\\{x\in\mathbb{R}^{p}\,|\,\|x\|_{\infty}\leq r\\}$.
For a given closed convex set $C$ and a vector $x$, we denote the Euclidean
projection of $x$ onto $C$ by $\Pi_{C}(x):=\mathop{\mbox{argmin}}_{y\in
C}{\|x-y\|}$.

##  2 Preliminaries

In this section, firstly, we formulate the dual problem of the DHGL problem.
Then we briefly introduce the proximal mapping and the Moreau envelope, which
play important roles in the ALM. Finally, we calculate the proximal mapping
and Moreau envelope for the negative log-determinant function
$-\log\det(\cdot)$, the regularization terms $Q(\cdot)$ and $R(\cdot)$ in the
DHGL problem.

###  2.1 Duality and the optimality condition

The Lagrangian function associated with the problem (P) is

$\displaystyle L(\Theta,Z,V;Y):=$ $\displaystyle\left\langle
S,\Theta\right\rangle-\log\det(\Theta)+Q(Z)+R(V)+\delta_{\mathbb{S}_{+}^{p}}(\Theta)$
$\displaystyle-\left\langle
Y,Z+\mathcal{T}V-\Theta\right\rangle,\quad\forall\left(\Theta,Z,V,Y\right)\in\mathbb{S}^{p}\times\mathbb{S}^{p}\times\mathbb{M}^{p}\times\mathbb{S}^{p},$

where $\delta_{\mathbb{S}^{p}_{+}}(\cdot)$ is an indicator function whose
effective domain consists of p-dimensional positive semidefinite matrices.
Furthermore, we can obtain the corresponding dual problem:

$\begin{array}[]{cl}{\underset{Y,\Omega,\Lambda\in\mathbb{S}^{p}}{\max}}&{\log\det(\Omega)-Q^{*}(Y)-R^{*}(\Lambda)+p}\\\
{\text{s.t.}}&{S-\Omega+Y=0,}\\\ &{\Lambda-\mathcal{T}^{*}Y=0,}\\\
&{\Omega\succeq 0,}\end{array}$ (D)

where $\mathcal{T}^{*}\colon\mathbb{S}^{p}\to\mathbb{M}^{p}$ is the adjoint
operator of $\mathcal{T}$, $Q^{*}(\cdot)$ and $R^{*}(\cdot)$ are the Fenchel
conjugate functions of $Q(\cdot)$ and $R(\cdot)$, defined by

$\begin{array}[]{l}Q^{*}(Y)=\sup\left\\{\left\langle Z,Y\right\rangle-Q(Z)\mid
Z\in\mathbb{S}^{p}\right\\},\quad\forall\,{Y}\in\mathbb{S}^{p},\\\
R^{*}(\Lambda)=\sup\left\\{\langle V,\Lambda\rangle-R(V)\mid
V\in\mathbb{M}^{p}\right\\},\quad\forall\,\Lambda\in\mathbb{M}^{p}.\end{array}$

Note that for any $Y\in\mathbb{S}^{p}$, we have $\mathcal{T}^{*}Y=2Y$. The KKT
condition for the problems (P) and (D) is

$\left\\{\begin{array}[]{l}{Z+\mathcal{T}V=\Theta,}\\\ {S-\Omega+Y=0,}\\\
{\Theta\Omega=I,}\\\ {Y\in\partial Q(Z),}\\\ {\mathcal{T}^{*}Y\in\partial
R(V),}\\\ {\Theta\succeq 0,\quad\Omega\succeq 0.}\\\ \end{array}\right.$ (2)

###  2.2 The proximal mapping and the Moreau envelope

In a real finite dimensional Euclidean space $\mathcal{H}$, which is equipped
with an inner product $\langle\cdot,\cdot\rangle$ and its induced norm
$\|\cdot\|$, given a closed convex function
$f:\mathcal{H}\rightarrow\mathbb{R}$, the Moreau envelope of $f$ is defined by

$e_{f}(x):=\min_{y\in\mathcal{H}}\left\\{f(y)+\frac{1}{2}\|y-x\|^{2}\right\\}.$

The unique minimizer of the above minimization problem is called the proximal
mapping of $f$ at $x$, and we denote it by $\textrm{prox}_{f}(x)$. That is,

$\textrm{prox}_{f}(x):=\mathop{\arg\min}_{y\in\mathcal{H}}\left\\{f(y)+\frac{1}{2}\|y-x\|^{2}\right\\}.$

It is proven in [34, Theorem 2.26] that the envelope function $e_{f}(\cdot)$
is convex and continuously differentiable with

$\nabla e_{f}(x)=x-\textrm{prox}_{f}(x)=\textrm{prox}_{f^{*}}(x),$

where the last equality is derived by the Moreau identity, which states that
the following equation holds for any $t>0$ and $x\in\mathcal{H}$,

$\textrm{prox}_{tf}(x)+t\textrm{prox}_{f^{*}/t}(x/t)=x.$

###  2.3 Properties of the negative log-determinant function

For convenience, we denote $h(X)=-\log\det(X)$ for any $X\succeq 0$. The
proximal mapping of $h(\cdot)$ and its corresponding Hessian is well studied
in [39, Lemma 2.1, Lemma 3.1]. We state them in the following two
propositions.

######  Proposition 2.1.

For any $X\in\mathbb{S}^{p}$, let its eigendecomposition be
$X=P\textrm{Diag}(d)P^{T}$, where $d$ is the vector of eigenvalues in the
descending order and $P$ is an orthonormal matrix whose columns are the
corresponding eigenvectors. Given $\gamma>0$, we define two scalar functions

$\displaystyle\phi_{\gamma}^{+}(x):=\frac{\sqrt{x^{2}+4\gamma}+x}{2}$ and
$\displaystyle\phi_{\gamma}^{-}(x):=\frac{\sqrt{x^{2}+4\gamma}-x}{2},\,\forall\,x\in\mathbb{R}.$

We also define their vector counterparts
$\phi_{\gamma}^{+}(\cdot):\mathbb{R}^{p}\to\mathbb{R}^{p}$ and
$\phi_{\gamma}^{-}(\cdot):\mathbb{R}^{p}\to\mathbb{R}^{p}$ as

$\displaystyle\phi_{\gamma}^{+}(d):=\begin{pmatrix}\phi_{\gamma}^{+}(d_{1})\\\
\vdots\\\ \phi_{\gamma}^{+}(d_{p})\\\ \end{pmatrix}$ and
$\displaystyle\phi_{\gamma}^{-}(x):=\begin{pmatrix}\phi_{\gamma}^{-}(d_{1})\\\
\vdots\\\ \phi_{\gamma}^{-}(d_{p})\\\ \end{pmatrix},\ \forall
d\in\mathbb{R}^{p},$

and define their matrix counterparts
$\phi_{\gamma}^{+}(\cdot)\colon\mathbb{S}^{p}\to\mathbb{S}^{p}$ and
$\phi_{\gamma}^{-}(\cdot):\mathbb{S}^{p}\to\mathbb{S}^{p}$ as

$\displaystyle\phi_{\gamma}^{+}(X):=P\operatorname{Diag}(\phi_{\gamma}^{+}(d))P^{T}$
and
$\displaystyle\phi_{\gamma}^{-}(X):=P\operatorname{Diag}(\phi_{\gamma}^{-}(d))P^{T},\
\forall X\in\mathbb{S}^{p}.$

Now, we have the following two properties:

  * (i)

Both of $\phi_{\gamma}^{+}(X)$ and $\phi_{\gamma}^{-}(X)$ are positive
definite for any $X\in\mathbb{S}^{p}$, and

$X=\phi_{\gamma}^{+}(X)-\phi_{\gamma}^{-}(X),\quad\langle\phi_{\gamma}^{+}(X),\phi_{\gamma}^{-}(X)\rangle=\gamma
I.$

  * (ii)

The mapping $\phi_{\gamma}^{+}(\cdot)$ is continuously differentiable and its
Fréchet derivative
$(\phi_{\gamma}^{+})^{\prime}(X)$ is

$(\phi_{\gamma}^{+})^{\prime}(X)[H]=P(M\circ(P^{T}HP))P^{T},\quad\forall
H\in\mathbb{S}^{p},$

where $M\in\mathbb{S}^{p}$ is defined by

$M_{ij}=\frac{\phi_{\gamma}^{+}(d_{i})+\phi_{\gamma}^{+}(d_{j})}{\sqrt{d_{i}^{2}+4\gamma}+\sqrt{d_{j}^{2}+4\gamma}},\quad
i,j=1,\dots,p.$

######  Proposition 2.2.

Given $G\in\mathbb{S}^{p}$ and $\sigma>0$, we have

$\displaystyle\min_{\Omega\succeq
0}\Big{\\{}\frac{1}{2\sigma}\|G-\Omega\|^{2}-\log\det(\Omega)\Big{\\}}=\frac{1}{2\sigma}\|\phi_{\sigma}^{-}(G)\|^{2}-\log\det(\phi_{\sigma}^{+}(G)),$

where the minimizer is attained at $\Omega=\phi_{\sigma}^{+}(G)$. That is,
${\rm prox}_{\sigma h}(G)=\phi_{\sigma}^{+}(G)$.

###  2.4 Properties of the regularization term $Q(\cdot)$

The regularization term $Q(Z)$ is called the graphical Lasso regularizer [13]:

$Q(Z)=\lambda_{1}\|Z-{\rm
diag}(Z)\|_{1}=\sum_{j=1}^{p}\lambda_{1}\|{\left(Z-{\rm
diag}(Z)\right)}_{j}\|_{1}.$

For $j=1,\dots,p$, let
$Z_{[j]}:=[Z_{1j},\dots,Z_{j-1,j},0,Z_{j+1,j},\dots,Z_{pj}]^{T}\in\mathbb{R}^{p}$,
and define $\mathcal{P}_{j}:\mathbb{M}^{p}\to\mathbb{R}^{p}$ as a linear map
satisfying $\mathcal{P}_{j}Z=Z_{[j]}$. For $x\in\mathbb{R}^{p}$, the adjoint
$\mathcal{P}_{j}^{*}:\mathbb{R}^{p}\to\mathbb{M}^{p}$ at $x$ is defined as
$\mathcal{P}_{j}^{*}x$ with placing $x$ in the j-th column of a $p\times p$
all-zero matrix $X$, and then setting $X_{jj}$ to be zero. Additionally,
define an operator $\mathcal{P}:\mathbb{M}^{p}\to\mathbb{R}^{p^{2}}$ as
$\mathcal{P}:=[\mathcal{P}_{1};\mathcal{P}_{2};\dots;\mathcal{P}_{p}]$, such
that $\mathcal{P}X$ conducts the operation $\mathcal{P}_{j}$ on the j-th
column of $X$ for all $j=1,\dots,p$, which actually sets the diagonal elements
of $X$ to be zeros and then stacks each column of the corresponding matrix to
be a column vercor. Thus
$Q(Z)=\sum_{j=1}^{p}\lambda_{1}\|\mathcal{P}_{j}Z\|_{1}$. Then we can compute
the closed-form expression of the proximal mapping $\textrm{prox}_{Q}(\cdot)$,
which is stated in the following proposition. We ignore its proof since it is
very similar to that of Proposition 2.4.

######  Proposition 2.3.

For any $Y\in\mathbb{S}^{p}$, the proximal operator of $Q(\cdot)$ at $Y$ can
be computed as

$\operatorname{prox}_{Q}(Y)=\mathop{\arg\min}_{Z\in\mathbb{S}^{p}}\left\\{\frac{1}{2}\|Z-Y\|^{2}+Q(Z)\right\\}=Y-\mathcal{P}^{*}\Pi_{\mathcal{B}^{p,\lambda_{1}}_{\infty}}(\mathcal{P}Y),$

where
$\mathcal{B}^{p,\lambda_{1}}_{\infty}:=\underbrace{\mathcal{B}_{\infty}^{\lambda_{1}}\times\dots\times\mathcal{B}_{\infty}^{\lambda_{1}}}_{p}$,
and
$\Pi_{\mathcal{B}^{p,\lambda_{1}}_{\infty}}(\cdot):\mathbb{R}^{p^{2}}\to\mathbb{R}^{p^{2}}$
is a block-wise projection onto the $\ell_{\infty}$ norm ball with the radius
$\lambda_{1}$. Specifically, the j-th block of
$\Pi_{\mathcal{B}_{\infty}^{p,\lambda_{1}}}(\mathcal{P}Y)$ is

$[\Pi_{\mathcal{B}_{\infty}^{p,\lambda_{1}}}(\mathcal{P}Y)]_{(j-1)p+1:jp}=\Pi_{\mathcal{B}_{\infty}^{\lambda_{1}}}(\mathcal{P}_{j}Y)=\operatorname{sign}(\mathcal{P}_{j}Y)\circ\min(|\mathcal{P}_{j}Y|,\lambda_{1}),\,j=1,\ldots,p.$

Next, for any $Y\in\mathbb{S}^{p}$, we define an alternative for the Clarke
generalized Jacobian of $\textrm{prox}_{Q}(\cdot)$ at $Y$, i.e.,
$\partial\textrm{prox}_{Q}(Y)$, as below

$\widehat{\partial}\textrm{prox}_{Q}(Y):=\\{\mathcal{I}-\mathcal{P}^{*}\Sigma\mathcal{P}\mid\Sigma\in\partial\Pi_{\mathcal{B}_{\infty}^{p,\lambda_{1}}}(\mathcal{P}Y)\\},$
(3)

where $\partial\Pi_{\mathcal{B}_{\infty}^{p,\lambda_{1}}}(\mathcal{P}Y)$ is
the Clarke generalized Jacobian of
$\Pi_{\mathcal{B}_{\infty}^{p,\lambda_{1}}}(\cdot)$ at $\mathcal{P}Y$, which
takes the form

$\partial\Pi_{\mathcal{B}_{\infty}^{p,\lambda_{1}}}(\mathcal{P}Y)=\\{\textrm{Diag}(\Sigma_{1},\dots,\Sigma_{p})\mid\Sigma_{j}\in\partial\Pi_{\mathcal{B}_{\infty}^{\lambda_{1}}}(\mathcal{P}_{j}Y),\
j=1,\cdots,p\\}.$

Particularly, for any $y\in\mathbb{R}^{p}$, we have

$\partial\Pi_{\mathcal{B}_{\infty}^{\lambda_{1}}}(y)=\left\\{{\rm
Diag}(u)\middle|u\in\mathbb{R}^{p},\
u_{i}\in\left\\{\begin{array}[]{cl}\\{1\\}&\text{if}~{}|y_{i}|<\lambda_{1}\\\
\\{t\mid 0\leq t\leq 1\\}&\text{if}~{}|y_{i}|=\lambda_{1}\\\
\\{0\\}&\text{if}~{}|y_{i}|>\lambda_{1}\end{array}\right.,\
i=1,\cdots,p\right\\}.$

Note that it has been proven in [19, Example 2.5] that for any
$H\in\mathbb{S}^{p}$,

$\widehat{\partial}\textrm{prox}_{Q}(Y)[H]=\partial\textrm{prox}_{Q}(Y)[H].$

###  2.5 Properties of the regularization term $R(\cdot)$

The regularization term $R(\cdot)$ is essentially the weighted group graphical
Lasso regularizer:

$\displaystyle R(V)=$
$\displaystyle\lambda_{2}\sum_{j\notin\mathcal{D}}\|{\left(V-{\rm
diag}{(V)}\right)}_{j}\|_{1}+\lambda_{3}\sum_{j\notin\mathcal{D}}\|{\left(V-{\rm
diag}{(V)}\right)}_{j}\|$
$\displaystyle+\lambda_{4}\sum_{j\in\mathcal{D}}\|{\left(V-{\rm
diag}{(V)}\right)}_{j}\|_{1}+\lambda_{5}\sum_{j\in\mathcal{D}}\|{\left(V-{\rm
diag}{(V)}\right)}_{j}\|,\quad\forall\,V\in\mathbb{M}^{p}.$

For simplicity, for any $V\in\mathbb{M}^{p}$, we denote

$\psi(V)=\sum_{j=1}^{p}w_{1,j}\|{\left(V-{\rm
diag}{(V)}\right)}_{j}\|_{1},\quad\varphi(V)=\sum_{j=1}^{p}w_{2,j}\|{\left(V-{\rm
diag}{(V)}\right)}_{j}\|,$

where

$w_{1,j}=\left\\{\begin{array}[]{cl}\lambda_{2},&\text{if}~{}j\notin\mathcal{D},\\\
\lambda_{4},&\text{if}~{}j\in\mathcal{D}.\end{array}\right.\quad
w_{2,j}=\left\\{\begin{array}[]{cl}\lambda_{3},&\text{if}~{}j\notin\mathcal{D},\\\
\lambda_{5},&\text{if}~{}j\in\mathcal{D}.\end{array}\right.$

Then, borrowing the linear operators $\mathcal{P}_{j},j=1,\dots,p$, we can
write

$\displaystyle
R(V)=\psi(V)+\varphi(V)=\sum_{j=1}^{p}(w_{1,j}\|\mathcal{P}_{j}V\|_{1}+w_{2,j}\|\mathcal{P}_{j}V\|).$
(4)

The next proposition provides a crucial property. It reveals that the proximal
mapping $\textrm{prox}_{R}(\cdot)$ of $R=\psi+\varphi$ can be decomposed into
the composition of the proximal mappings $\textrm{prox}_{\psi}(\cdot)$ and
$\textrm{prox}_{\varphi}(\cdot)$.

######  Proposition 2.4.

For any $Y\in\mathbb{M}^{p}$ and $\sigma>0$, the proximal operator of
$R(\cdot)$ at $Y$ is

$\displaystyle\operatorname{prox}_{R}(Y)$ $\displaystyle=$
$\displaystyle\operatorname{prox}_{\varphi}\circ\operatorname{prox}_{\psi}(Y)=X-\mathcal{P}^{*}\Pi_{\mathcal{B}^{w_{2}}_{2}}(\mathcal{P}X),$

where

$\displaystyle
X=\operatorname{prox}_{\psi}(Y)=Y-\mathcal{P}^{*}\Pi_{\mathcal{B}^{w_{1}}_{\infty}}(\mathcal{P}Y),$
$\displaystyle\mathcal{B}^{w_{2}}_{2}:=\mathcal{B}_{2}^{w_{2,1}}\times\dots\times\mathcal{B}_{2}^{w_{2,p}},\quad\mathcal{B}^{w_{1}}_{\infty}:=\mathcal{B}_{\infty}^{w_{1,1}}\times\dots\times\mathcal{B}_{\infty}^{w_{1,p}},$

and $\Pi_{\mathcal{B}^{w_{1}}_{\infty}}(\cdot)$ and
$\Pi_{\mathcal{B}^{w_{2}}_{2}}(\cdot)$ are block-wise projections onto
$\ell_{\infty}$ and $\ell_{2}$ norm ball, respectively. Specifically, for
$j=1,\ldots,p$, the j-th block of
$\Pi_{\mathcal{B}^{w_{1}}_{\infty}}(\mathcal{P}Y)$ and
$\Pi_{\mathcal{B}^{w_{2}}_{2}}(\mathcal{P}X)$ are

$\displaystyle\begin{aligned}
[\Pi_{\mathcal{B}^{w_{1}}_{\infty}}(\mathcal{P}Y)]_{(j-1)p+1:jp}&=\Pi_{\mathcal{B}_{\infty}^{w_{1,j}}}(\mathcal{P}_{j}Y)=\operatorname{sign}(\mathcal{P}_{j}Y)\circ\min(|\mathcal{P}_{j}Y|,w_{1,j}),\end{aligned}$

and

$\displaystyle\begin{aligned}
[\Pi_{\mathcal{B}^{w_{2}}_{2}}(\mathcal{P}X)]_{(j-1)p+1:jp}&=\Pi_{\mathcal{B}_{2}^{w_{2,j}}}(\mathcal{P}_{j}X)=\left\\{\begin{array}[]{cl}w_{2,j}\frac{\mathcal{P}_{j}X}{\|\mathcal{P}_{j}X\|}&\text{if}~{}\|\mathcal{P}_{j}X\|>w_{2,j},\\\
\mathcal{P}_{j}X&\text{otherwise,}\end{array}\right.\end{aligned}$

respectively.

###### Proof.

See Appendix A. ∎

Note that even though the proximal mapping $\textrm{prox}_{R}(\cdot)$ can be
obtained by the composition of $\textrm{prox}_{\psi}(\cdot)$ and
$\textrm{prox}_{\varphi}(\cdot)$, the classical chain rule usually cannot work
in calculating the generalized Jacobian of $\textrm{prox}_{R}(\cdot)$ since
$\textrm{prox}_{\psi}(\cdot)$ and $\textrm{prox}_{\varphi}(\cdot)$ are non-
differentiable [7]. Therefore, we define a surrogate generalized Jacobian of
$\textrm{prox}_{R}(\cdot)$ as follows.

######  Definition 2.1.

For any $Y\in\mathbb{M}^{p}$,
$\widehat{\partial}\operatorname{prox}_{R}(Y):\mathbb{M}^{p}\rightrightarrows\mathbb{M}^{p}$
is defined as

$\widehat{\partial}\operatorname{prox}_{R}(Y):=\\{(\mathcal{I}-\mathcal{P}^{*}\Delta\mathcal{P})(\mathcal{I}-\mathcal{P}^{*}\Xi\mathcal{P})\,|\,\Delta\in\partial\Pi_{\mathcal{B}^{w_{2}}_{2}}(\mathcal{P}X),\
\Xi\in\partial\Pi_{\mathcal{B}^{w_{1}}_{\infty}}(\mathcal{P}Y)\\},$

where $X=Y-\mathcal{P}^{*}\Pi_{\mathcal{B}^{w_{1}}_{\infty}}(\mathcal{P}Y)$,
and

$\displaystyle\partial\Pi_{\mathcal{B}^{w_{1}}_{\infty}}(\mathcal{P}Y)=\\{\operatorname{Diag}(\Xi_{1},\dots,\Xi_{p})\mid\Xi_{j}\in\partial\Pi_{\mathcal{B}_{\infty}^{w_{1,j}}}(\mathcal{P}_{j}Y),\
j=1,\cdots,p\\},$
$\displaystyle\partial\Pi_{\mathcal{B}^{w_{2}}_{2}}(\mathcal{P}X)=\\{\operatorname{Diag}(\Delta_{1},\dots,\Delta_{p})\mid\Delta_{j}\in\partial\Pi_{\mathcal{B}_{2}^{w_{2,j}}}(\mathcal{P}_{j}X),\
j=1,\cdots,p\\}.$

For any $y\in\mathbb{R}^{p}$ and $j=1,\ldots,p$, we have

$\displaystyle\begin{aligned}
&\partial\Pi_{\mathcal{B}_{\infty}^{w_{1,j}}}(y)\\!=\\!\left\\{{\rm
Diag}(u)\middle|u\\!\in\\!\mathbb{R}^{p},u_{i}\in\left\\{\begin{array}[]{cl}\\{1\\}&\text{if}~{}|y_{i}|<w_{1,j}\\\
\\{t\mid 0\leq t\leq 1\\}&\text{if}~{}|y_{i}|=w_{1,j}\\\
\\{0\\}&\text{if}~{}|y_{i}|>w_{1,j}\end{array}\right.,j=1,\cdots,p\right\\},\\\
&\partial\Pi_{\mathcal{B}^{w_{2,j}}_{2}}(y)=\left\\{\begin{array}[]{cl}\left\\{\frac{w_{2,j}}{\|y\|}(I-\frac{yy^{T}}{\|y\|^{2}})\right\\}&\text{if}~{}\|y\|>w_{2,j}\\\
\Big{\\{}I-t\frac{yy^{T}}{(w_{2,j})^{2}}\,\Big{|}\,0\leq t\leq
1\Big{\\}}&\text{if}~{}\|y\|=w_{2,j}\\\
\\{I\\}&\text{if}~{}\|y\|<w_{2,j}\end{array}\right..\end{aligned}$

Similarly, by [19, Example 2.5] we have that for any $H\in\mathbb{M}^{p}$,
$\widehat{\partial}\textrm{prox}_{R}(Y)[H]=\partial\textrm{prox}_{R}(Y)[H]$.

######  Proposition 2.5.

$\operatorname{prox}_{Q}(\cdot)$ and $\operatorname{prox}_{R}(\cdot)$ are
strongly semismooth. That is, let $X,Y\in\mathbb{S}^{p}$, we have

$\operatorname{prox}_{Q}(Y)-\operatorname{prox}_{Q}(X)-\mathcal{W}^{Q}(Y-X)=O(\|Y-X\|^{2}),\quad\forall\,\mathcal{W}^{Q}\in\widehat{\partial}\operatorname{prox}_{Q}(Y),$
(5)

and for any $U,V\in\mathbb{M}^{p}$, we have

$\operatorname{prox}_{R}(V)-\operatorname{prox}_{R}(U)-\mathcal{W}^{R}(V-U)=O(\|V-U\|^{2}),\quad\forall\,\mathcal{W}^{R}\in\widehat{\partial}\operatorname{prox}_{R}(V).$
(6)

###### Proof.

Since $\operatorname{prox}_{Q}(\cdot)$ is piecewise affine, according to [12,
Proposition 7.4.7] we have (5). Similar to [44, Theorem 3.1], we can easily
obtain that (6) holds. ∎

##  3 A two-phase algorithm

In this section, we propose a two-phase algorithm to solve the DHGL problem.
In Phase I, we adopt the dADMM to generate a nice initial point to warm start
the Phase II algorithm. It is worth noting that though the dADMM can be used
alone to solve the problem, it is not efficient enough to achieve an accurate
solution needed in practical applications. In Phase II, we apply the
superlinearly convergent ALM to compute a more accurate solution. Furthermore,
in order to fully leverage the sparsity structure of the Jacobians, we use the
SSN method to solve the subproblems.

###  3.1 Phase I: dADMM

The minimization form of the dual problem (D) is

$\begin{array}[]{cl}\underset{Y,\Omega,\Gamma,\Lambda\in\mathbb{S}^{p}}{\min}&-\log\det(\Omega)+Q^{*}(\Gamma)+R^{*}(\Lambda)+\delta_{\mathbb{S}^{p}_{+}}(\Omega)-p\\\
\text{s.t.}&S-\Omega+Y=0\color[rgb]{0,0,0.9},\\\
&{\Gamma-Y=0}\color[rgb]{0,0,0.9},\\\
&{\Lambda-\mathcal{T}^{*}Y=0}\color[rgb]{0,0,0.9}.\end{array}$

Its corresponding Lagrangian function is given by

$\displaystyle\mathcal{L}(Y,\Omega,\Gamma,\Lambda;\Theta,Z,V):=$
$\displaystyle\,{-\log\det(\Omega)+Q^{*}(\Gamma)+R^{*}(\Lambda)+\delta_{\mathbb{S}_{+}^{p}}(\Omega)-p}$
(7) $\displaystyle\qquad{-\left\langle\Theta,S-\Omega+Y\right\rangle-\langle
Z,\Gamma-Y\rangle-\langle V,\Lambda-\mathcal{T}^{*}Y\rangle,}$
$\displaystyle{\qquad\,\forall\,\Omega,Y,\Theta,Z,\Gamma,\Lambda\in\mathbb{S}^{p},V\in\mathbb{M}^{p}.}$

For $\sigma>0$, the associated augmented Lagrangian function is

$\displaystyle\mathcal{L}_{\sigma}(Y,\Omega,\Gamma,\Lambda;\Theta,Z,V)=$
$\displaystyle\,\mathcal{L}(Y,\Omega,\Gamma,\Lambda;\Theta,Z,V)+\frac{\sigma}{2}\|S-\Omega+Y\|^{2}$
(8)
$\displaystyle\qquad+\frac{\sigma}{2}\|\Gamma-Y\|^{2}+\frac{\sigma}{2}\|\Lambda-\mathcal{T}^{*}Y\|^{2}.$

Based on the augmented Lagrangian function (8), we can adopt the dADMM. We
summarize the dADMM for solving the problem (D) in Algorithm 1. One may refer
to Section 2 for the detailed computations of the related proximal mappings.

Algorithm 1 dADMM for DHGL

0:
$\Omega^{0},\Theta^{0}\in\mathbb{S}_{++}^{p},Z^{0},\Gamma^{0},\Lambda^{0}\in\mathbb{S}^{p},V^{0}\in\mathbb{M}^{p},\sigma>0$,
$\tau\in(0,(1+\sqrt{2})/2)$, (a typical choice is $\tau=1.618$). Let $k=0$,
iterate as follows:

1: Compute
$Y^{k+1}=\frac{1}{6\sigma}(\Theta^{k}-Z^{k}-V^{k}-{V^{k}}^{T})+\frac{1}{6}(\Gamma^{k}+2\Lambda^{k}-(S-\Omega^{k})).$

2: Compute
$\left\\{\begin{array}[]{l}\Omega^{k+1}=\phi_{1/\sigma}^{+}(S+Y^{k+1}-\Theta^{k}/\sigma),\\\
\Gamma^{k+1}=\textrm{prox}_{Q^{*}/\sigma}(Y^{k+1}+Z^{k}/\sigma),\\\
\Lambda^{k+1}=\textrm{prox}_{R^{*}/\sigma}(\mathcal{T}^{*}Y^{k+1}+V^{k}/\sigma).\end{array}\right.$

3: Compute
$\left\\{\begin{array}[]{l}\Theta^{k+1}=\Theta^{k}-\tau\sigma(S-\Omega^{k+1}+Y^{k+1}),\\\
Z^{k+1}=Z^{k}-\tau\sigma(\Gamma^{k+1}-Y^{k+1}),\\\
V^{k+1}=V^{k}-\tau\sigma(\Lambda^{k+1}-\mathcal{T}^{*}Y^{k+1}).\end{array}\right.$

4: $k\leftarrow k+1$, go to Step 1.

Now we give the convergence result of the dADMM for solving the DHGL problem.
For clarity and brevity, we simply use $X_{P}$ and $X_{D}$ to denote the
primal and dual variables, respectively, i.e.,
$X_{P}:=(\Theta,Z,V)\in\mathbb{S}^{p}_{+}\times\mathbb{S}^{p}\times\mathbb{M}^{p}$
and
$X_{D}:=(Y,\Omega,\Gamma,\Lambda)\in\mathbb{S}^{p}\times\mathbb{S}^{p}_{+}\times\mathbb{S}^{p}\times\mathbb{S}^{p}$.
Based on [17, Theorem 2], we have the following convergence result.

######  Theorem 3.1.

Let $\tau\in(0,(1+\sqrt{5})/2)$. Let $\\{(X_{D}^{k},X_{P}^{k})\\}$ be the
sequence generated by the dADMM. Then
$\\{X_{D}^{k}:=(Y^{k},\Omega^{k},\Gamma^{k},\Lambda^{k})\\}$ converges to an
optimal solution of the problem (D), and
$\\{X_{P}^{k}:=(\Theta^{k},Z^{k},V^{k})\\}$ converges to an optimal solution
of the problem (P).

###  3.2 Phase II: ALM

For the ALM, we still work on the augmented Lagrangian function. At the k-th
iteration of the ALM for solving the problem (D), the most important task is
to solve the following subproblem

$\min_{Y,\Omega,\Gamma,\Lambda\in\mathbb{S}^{p}}\left\\{\mathcal{L}_{\sigma_{k}}(Y,\Omega,\Gamma,\Lambda;\Theta^{k},Z^{k},V^{k})\right\\},$
(9)

where $\sigma_{k}$ is a given parameter, $(\Theta^{k},Z^{k},V^{k})$ are given
points. Since $\Omega,\Gamma,\Lambda$ are independent of each other, we can
solve them separately as follows

$\begin{array}[]{l}\Omega^{\prime}=\phi_{1/\sigma_{k}}^{+}(S+Y-\frac{1}{\sigma_{k}}\Theta^{k}),\\\
\Gamma^{\prime}=\textrm{prox}_{Q^{*}/\sigma_{k}}(Y+\frac{1}{\sigma_{k}}Z^{k}),\\\
\Lambda^{\prime}=\textrm{prox}_{R^{*}/\sigma_{k}}(\mathcal{T}^{*}Y+\frac{1}{\sigma_{k}}V^{k}).\end{array}$

Then we can substitute them to the problem (9), we actually solve the
subproblem

$\displaystyle\min_{Y\in\mathbb{S}^{p}}$ $\displaystyle\Phi_{k}(Y),$ (10)

where

$\displaystyle\Phi_{k}(Y)=\min_{\Omega,\Gamma,\Lambda\in\mathbb{S}^{p}}\left\\{\mathcal{L}_{\sigma_{k}}(Y,\Omega,\Gamma,\Lambda;\Theta^{k},Z^{k},V^{k})\right\\}$
$\displaystyle=~{}\sigma_{k}\left[(e_{h/\sigma_{k}})(S+Y-\frac{1}{\sigma_{k}}\Theta^{k})+(e_{Q^{*}/\sigma_{k}})(Y+\frac{1}{\sigma_{k}}Z^{k})+(e_{R^{*}/\sigma_{k}})(\mathcal{T}^{*}Y+\frac{1}{\sigma_{k}}V^{k})\right]$
$\displaystyle\quad-\frac{1}{2\sigma_{k}}\left(\|\Theta^{k}\|^{2}+\|Z^{k}\|^{2}+\|V^{k}\|^{2}\right)-p.$

For clarity, we present the framework of the ALM in Algorithm 2.

Algorithm 2 ALM for DHGL

0:
$\Theta^{0}\in\mathbb{S}^{p}_{++},Z^{0}\in\mathbb{S}^{p},V^{0}\in\mathbb{M}^{p};\sigma_{0}>0$.
Let $k=0$, iterate as follows:

1: Compute  $\displaystyle
Y^{k+1}\approx\mathop{\arg\min}_{Y\in\mathbb{S}^{p}}\Phi_{k}(Y).$

2: Compute
$\displaystyle\Omega^{k+1}=\phi_{1/\sigma_{k}}^{+}\big{(}S+Y^{k+1}-\frac{1}{\sigma_{k}}\Theta^{k}\big{)},$
$\displaystyle\Gamma^{k+1}=\textrm{prox}_{Q^{*}/\sigma_{k}}\big{(}Y^{k+1}+\frac{1}{\sigma_{k}}Z^{k}\big{)},$
$\displaystyle\Lambda^{k+1}=\textrm{prox}_{R^{*}/\sigma_{k}}\big{(}\mathcal{T}^{*}Y^{k+1}+\frac{1}{\sigma_{k}}V^{k}\big{)}.$

3: Compute
$\displaystyle\Theta^{k+1}=\Theta^{k}-\sigma_{k}\big{(}S-\Omega^{k+1}+Y^{k+1}\big{)},$
$\displaystyle Z^{k+1}=Z^{k}-\sigma_{k}\big{(}\Gamma^{k+1}-Y^{k+1}\big{)},$
$\displaystyle
V^{k+1}=V^{k}-\sigma_{k}\big{(}\Lambda^{k+1}-\mathcal{T}^{*}Y^{k+1}\big{)}.$

4: Update $\sigma_{k+1}\uparrow\sigma_{\infty}$ and $k\leftarrow k+1$, go to
Step 1.

####  3.2.1 An implementable stopping criteria for the ALM subproblems

Since the subproblem (9) often has no analytical solution, the stopping
criteria for the subproblem in the algorithm needs to be discussed. For
$X_{D}:=(Y,\Omega,\Gamma,\Lambda)\in\mathbb{S}^{p}\times\mathbb{S}^{p}_{+}\times\mathbb{S}^{p}\times\mathbb{S}^{p}$
and
$X_{P}:=(\Theta,Z,V)\in\mathbb{S}^{p}_{+}\times\mathbb{S}^{p}\times\mathbb{M}^{p}$,
we define functions $f_{k}$ and $g_{k}$ as follows

$\displaystyle f_{k}(X_{D})$
$\displaystyle:=\mathcal{L}_{\sigma_{k}}(X_{D};X_{P}^{k})=\sup_{X_{P}}\Big{\\{}\mathcal{L}(X_{D};X_{P})-\frac{1}{2\sigma_{k}}\|X_{P}-X_{P}^{k}\|^{2}\Big{\\}},$
$\displaystyle g_{k}(X_{P})$
$\displaystyle:=\inf_{X_{D}}\Big{\\{}\mathcal{L}(X_{D};X_{P})-\frac{1}{2\sigma_{k}}\|X_{P}-X_{P}^{k}\|^{2}\Big{\\}}$
$\displaystyle=\log\det(\Theta)-\langle S,\Theta\rangle-Q(Z)-R(V)$
$\displaystyle\quad-\frac{1}{2\sigma_{k}}(\|\Theta-\Theta^{k}\|^{2}+\|Z-Z^{k}\|^{2}+\|V-V^{k}\|^{2})\color[rgb]{0,0,0.9}-\delta_{F_{P}}(\Theta,Z,V),$

where $F_{P}:=\left\\{(\Theta,Z,V)\mid\Theta-Z-\mathcal{T}V=0\right\\}$.

To obtain $X_{D}^{k+1}$, the following criteria on the approximate computation
of the subproblem (9) were proposed in [33]:

$\displaystyle f_{k}(X_{D}^{k+1})-\inf
f_{k}(X_{D})\leq\varepsilon_{k}^{2}/(2\sigma_{k}),\quad\varepsilon_{k}\geq
0,\quad\sum_{k=0}^{\infty}\varepsilon_{k}<\infty,$ (A) $\displaystyle
f_{k}(X_{D}^{k+1})-\inf
f_{k}(X_{D})\leq\delta_{k}^{2}/(2\sigma_{k})\|X_{P}^{k+1}-X_{P}^{k}\|^{2},\quad
0\leq\delta_{k}\leq 1,\quad\sum_{k=0}^{\infty}\delta_{k}<\infty,$ (B)

The convergence results when executing the stopping criteria (A) and (B) for
(9) can be obtained directly by adopting [33, Theorem 4, 5]. However, the
stopping criteria (A) and (B) are usually not in a direct usage due to the
unknown value of $\inf f_{k}(X_{D})$. To obtain implementable stopping
criteria, we need to find an upper bound for $f_{k}(X_{D}^{k+1})-\inf
f_{k}(X_{D})$. Since

$\inf f_{k}(X_{D})=\sup g_{k}(X_{P})\geq g_{k}(X_{P}^{k+1}),$

we have

$f_{k}(X_{D}^{k+1})-\inf f_{k}(X_{D})\leq
f_{k}(X_{D}^{k+1})-g_{k}(X_{P}^{k+1}).$

Hence, we terminate the subproblem (9) if $(X_{D}^{k+1},X_{P}^{k+1})$
satisfies the following conditions: given nonnegative summable sequences
$\\{\varepsilon_{k}\\}$ and $\\{\delta_{k}\\}$ such that $\delta_{k}<1$ for
all $k\geq 0$,

$\displaystyle
f_{k}(X_{D}^{k+1})-g_{k}(X_{P}^{k+1})\leq\varepsilon_{k}^{2}/(2\sigma_{k}),\quad\varepsilon_{k}\geq
0,\quad\sum_{k=0}^{\infty}\varepsilon_{k}<\infty,$ (A′) $\displaystyle
f_{k}(X_{D}^{k+1})-g_{k}(X_{P}^{k+1})\leq\delta_{k}^{2}/(2\sigma_{k})\|X_{P}^{k+1}-X_{P}\|^{2},\quad
0\leq\delta_{k}\leq 1,\quad\sum_{k=0}^{\infty}\delta_{k}<\infty.$ (B′)

####  3.2.2 Convergence results of the ALM

The global convergence and asymptotically superlinear convergence of the ALM
have been carefully explored by [33], thus we can directly obtain the result
without much effort. Firstly, by denoting $f(X_{P})$ as the objective function
of (P)

$f(X_{P})=\langle
S,\Theta\rangle-\log\det(\Theta)+Q(Z)+R(V)+\delta_{F_{P}}(\Theta,Z,V),$

we can define the maximal monotone operator $\mathcal{T}_{f}:=\partial f$ and
its inverse operator

$\mathcal{T}^{-1}_{f}(U_{P})=\underset{X_{P}}{\arg\min}\Big{\\{}f_{k}(X_{P})-\langle
U_{P},X_{P}\rangle\Big{\\}}.$

Now we are ready to present the convergence result of the ALM for solving the
dual problem.

######  Theorem 3.2.

(i) Let $\\{(X_{D}^{k},X_{P}^{k})\\}$ be the sequence generated by the ALM
with the stopping criterion (A′). Then $\\{(X_{D}^{k},X_{P}^{k})\\}$ is
bounded, $\\{X_{D}^{k}\\}$ converges to an optimal solution of (D), and
$\\{X_{P}^{k}\\}$ converges to an optimal solution of (P).

(ii) Let $\rho$ be a positive number such that
$\rho>\sum_{k=0}^{\infty}\varepsilon_{k}$. Assume that there exists $\kappa>0$
such that

$\mathrm{dist}\big{(}X_{P},\mathcal{T}_{f}^{-1}(0)\big{)}\leq\kappa\mathrm{dist}\big{(}0,\mathcal{T}_{f}(X_{P})\big{)}$

for all $X_{P}$ satisfying
$\mathrm{dist}\big{(}X_{P},\mathcal{T}_{f}^{-1}(0)\big{)}\leq\rho$. Suppose
that the initial point $X_{P}^{0}$ satisfies

$\mathrm{dist}\big{(}X_{P}^{0},\mathcal{T}_{f}^{-1}(0)\big{)}\leq\rho-\sum_{k=0}^{\infty}\varepsilon_{k}.$

Let $\\{(X_{D}^{k},X_{P}^{k})\\}$ be the sequence generated by the ALM with
the stopping criteria (A′) and (B′). Then for $k\geq 0$, it holds that

$\mathrm{dist}\big{(}X_{P}^{k+1},\mathcal{T}_{f}^{-1}(0)\big{)}\leq\varrho_{k}\mathrm{dist}\big{(}X_{P}^{k},\mathcal{T}_{f}^{-1}(0)\big{)},$

where

$\displaystyle\varrho_{k}$ $\displaystyle:=$
$\displaystyle\left[(1+\delta_{k})\kappa\left(\kappa^{2}+\sigma_{k}^{2}\right)^{-1/2}+\delta_{k}\right]/\left(1-\delta_{k}\right)\rightarrow\varrho_{\infty}:=\kappa\left(\kappa^{2}+\sigma_{\infty}^{2}\right)^{-1/2},$
$\displaystyle\left(\varrho_{\infty}=0\text{ if
}\sigma_{\infty}=\infty\right).$

####  3.2.3 The SSN for solving the ALM subproblem

Solving the subproblem (10) is the most critical part for the efficiency of
the algorithm. By ignoring the superscript, the subproblem at each iteration
can be formulated as

$\displaystyle\min_{Y\in\mathbb{S}^{p}}\Big{\\{}\Phi(Y):=$
$\displaystyle~{}\sigma\left[(e_{h/\sigma})(S+Y-\frac{1}{\sigma}\overline{\Theta})+(e_{Q^{*}/\sigma})(Y+\frac{1}{\sigma}\overline{Z})+(e_{R^{*}/\sigma})(2Y+\frac{1}{\sigma}\overline{V})\right]$
(13)
$\displaystyle\qquad-\frac{1}{2\sigma}\left(\|\overline{\Theta}\|^{2}+\|\overline{Z}\|^{2}+\|\overline{V}\|^{2}\right)-p\Big{\\}},$

where $\overline{\Theta},\overline{Z},\overline{V}$ are fixed.

Note that $\Phi(\cdot)$ is strictly convex and continuously differentiable,
and for any $Y\in\mathbb{S}^{p}$,

$\displaystyle\nabla\Phi(Y)=$
$\displaystyle\,-\phi^{+}_{\sigma}(-\sigma(S+Y)+\overline{\Theta})+\textrm{prox}_{\sigma
Q}(\sigma Y+\overline{Z})+2\textrm{prox}_{\sigma R}(2\sigma Y+\overline{V}).$

Thus, the unique solution $\overline{Y}$ of (13) can be obtained by solving
the following nonsmooth system of equations

$\nabla\Phi(Y)=0.$

We adopt the SSN method to solve this nonsmooth system. However, it is
difficult to characterize the structure of $\partial(\nabla\Phi)(Y)$ exactly.
Therefore, we also define a multifunction $\mathcal{V}$ as a surrogate of
$\partial(\nabla\Phi)(Y)$. The multifunction
$\mathcal{V}\colon\mathbb{S}^{p}\rightrightarrows\mathbb{S}^{p}$ is defined as
follows

$\displaystyle\mathcal{V}(Y):=$
$\displaystyle~{}\sigma\Big{[}(\phi^{+}_{\sigma})^{\prime}(-\sigma(S+Y)+\overline{\Theta})+\widehat{\partial}\textrm{prox}_{\sigma
Q}(\sigma Y+\overline{Z})+4\widehat{\partial}\textrm{prox}_{\sigma R}(2\sigma
Y+\overline{V})\Big{]}.$

The SSN method is presented in Algorithm 3.

Algorithm 3 SSN for the ALM subproblem

0: $\eta\in(0,1),\beta\in(0,1]$, $\mu\in(0,1/2)$, and $\delta\in(0,1)$. Choose
$Y^{0}\in\mathbb{S}^{p}$ and $j=0$.

1: Select $\mathcal{W}^{Q}\in\widehat{\partial}\textrm{prox}_{\sigma Q}(\sigma
Y^{j}+\overline{Z})$,
$\mathcal{W}^{R}\in\widehat{\partial}\textrm{prox}_{\sigma R}(2\sigma
Y^{j}+\overline{V})$. Let
$W=\sigma[(\phi^{+}_{\sigma})^{\prime}(-\sigma(S+Y^{j})+\overline{\Theta})+\mathcal{W}^{Q}+4\mathcal{W}^{R}]$.
Solve the following linear system by the Preconditioned Conjugate Gradient
(PCG) method  $W[H]=-\nabla\Phi(Y^{j})$ (14) to find $H^{j}\in\mathbb{S}^{p}$
such that
$\|W[H^{j}]+\nabla\Phi(Y^{j})\|\leq\min(\eta,\|\nabla\Phi(Y^{j})\|^{1+\beta})$.

2: Set $\alpha_{j}=\delta^{m_{j}}$ , where $m_{j}$ is the smallest nonnegative
integer $m$ for which
$\Phi(Y^{j}+\delta^{m}H^{j})\leq\Phi(Y^{j})+\mu\delta^{m}\langle\nabla\Phi(Y^{j}),H^{j}\rangle.$

3: Set $Y^{j+1}=Y^{j}+\alpha_{j}H^{j}$.

4: $j\leftarrow j+1$, go to Step 1.

Similar to the proof of [40, Proposition 2], we can prove that $\nabla\Phi$ is
strongly semismooth in combination with Proposition 2.1 and Proposition 2.5,
based on which we can get the following convergence result of the SSN method.

######  Theorem 3.3.

Let $\\{Y^{j}\\}$ be the sequence generated by Algorithm 3, then $\\{Y^{j}\\}$
converges to the unique optimal solution $\overline{Y}$ of (13), and the
convergence rate is at least superlinear:

$\|Y^{j+1}-\overline{Y}\|=\mathcal{O}(\|Y^{j}-\overline{Y}\|^{1+\beta}),\quad\beta\in(0,1].$

##  4 Numerical issues for solving the subproblem (14)

The most important part in implementing Algorithm 3 is how to efficiently
solve the subproblem (14). The Newton system we need to solve is

$\Big{[}(\phi^{+}_{\sigma})^{\prime}(-\sigma(S+Y)+\overline{\Theta})[H]+\mathcal{W}^{Q}(H)+4\mathcal{W}^{R}(H)\Big{]}=-\sigma^{-1}\nabla\Phi(Y),$
(15)

where $\mathcal{W}^{Q}\in\widehat{\partial}\textrm{prox}_{\sigma Q}(\sigma
Y+\overline{Z})$ and
$\mathcal{W}^{R}\in\widehat{\partial}\textrm{prox}_{\sigma R}(2\sigma
Y+\overline{V})$. Here, we can choose $\mathcal{W}^{Q}$ and $\mathcal{W}^{R}$
such that

$\displaystyle\mathcal{W}^{Q}(H)=H-U\circ
H,\quad\text{where}~{}U_{ij}=\left\\{\begin{array}[]{cl}0,&\text{if}~{}|(\sigma
Y+\overline{Z})_{ij}|>\sigma\lambda_{1}~{}\text{or}~{}i=j,\\\
1,&\text{otherwise},\end{array}\right.$

and

$\displaystyle\mathcal{W}^{R}(H)=H-\Sigma^{\psi}\circ
H-\mathcal{P}^{*}\Sigma^{\varphi}\mathcal{P}(H-\Sigma^{\psi}\circ H),$

where

$\displaystyle\Sigma^{\psi}_{ij}=\left\\{\begin{array}[]{cl}0,&\text{if}~{}|(2\sigma
Y+\overline{V})_{ij}|>\sigma w_{1,j}~{}\text{or}~{}i=j,\\\
1,&\text{otherwise},\end{array}\right.$

and

$\displaystyle\Sigma^{\varphi}=\textrm{Diag}(\Sigma^{\varphi}_{1},\ldots,\Sigma^{\varphi}_{p}),$

for $j=1,\ldots,p$,

$\displaystyle\begin{aligned}
&\Sigma^{\varphi}_{j}=\left\\{\begin{array}[]{cl}\frac{\sigma
w_{2,j}}{\|z_{j}\|}(I-\frac{z_{j}z_{j}^{T}}{\|z_{j}\|^{2}}),&\text{if}~{}\|z_{j}\|>\sigma
w_{2,j},\\\ I,&\text{if}~{}\|z_{j}\|\leq\sigma
w_{2,j},\end{array}\right.\end{aligned}$

where $z_{j}=\mathcal{P}_{j}X$ with $X=2\sigma
Y+\overline{V}-\mathcal{P}^{*}\Pi_{\mathcal{B}^{w_{1}}_{\infty}}(\mathcal{P}(2\sigma
Y+\overline{V}))$.

We need to apply the PCG method to solve the system (15). To achieve faster
convergence, we hope to devise an easy-to-compute preconditioner. Similar to
[39], we define an operator
$\mathcal{T}^{\phi}:\mathbb{S}^{p}\rightarrow\mathbb{S}^{p}$ as
$\mathcal{T}^{\phi}(H):=P(M\circ(P^{T}HP))P^{T}$. As we know, the standard
basis in $\mathbb{S}^{p}$ is given by
$\\{E_{ij}:=\alpha_{ij}(e_{i}e_{j}^{T}+e_{j}e_{i}^{T}):1\leq i\leq j\leq
p\\}$, where $e_{i}$ is the i-th unit vector in $\mathbb{R}^{p}$ and
$\alpha_{ij}=1/\sqrt{2}$ if $i\neq j$ and $\alpha_{ij}=1/2$ otherwise. Let
$\mathbf{T}^{\phi},\mathbf{W}^{Q}$ and $\mathbf{W}^{R}$ denote the matrix
representations of the linear operators $\mathcal{T}^{\phi},\mathcal{W}^{Q}$
and $\mathcal{W}^{R}$, respectively. Then the diagonal element of
$\textbf{T}^{\phi}$ with respect to the basis element $E_{ij}$ is

$\displaystyle{\bf T}^{\phi}_{(ij),(ij)}$ $\displaystyle=$
$\displaystyle\langle E_{ij},\mathcal{T}^{\phi}(E_{ij})\rangle=\langle
P^{T}E_{ij}P,M\circ(P^{T}E_{ij}P)\rangle$ (20) $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}((P\circ P)M(P\circ
P)^{T})_{ij}+\langle v^{(ij)},Mv^{(ij)}\rangle,&\mbox{if
$i\not=j$},\\\\[5.0pt] ((P\circ P)M(P\circ
P)^{T})_{ij},&\mbox{otherwise},\end{array}\right.$

where $v^{(ij)}=P_{i}\circ P_{j}$ and $P_{i},P_{j}$ are the $i$th and $j$th
rows of $P$, respectively. The diagonal element of $\textbf{W}^{Q}$ with
respect to the basis element $E_{ij}$ is

$\displaystyle\textbf{W}^{Q}_{(ij),(ij)}=\langle
E_{ij},\mathcal{W}^{Q}(E_{ij})\rangle=\langle E_{ij},E_{ij}-U\circ
E_{ij}\rangle=1-U_{ij}.$ (21)

The diagonal element of $\textbf{W}^{R}$ with respect to the basis element
$E_{ij}$ is

$\displaystyle\textbf{W}^{R}_{(ij),(ij)}$ $\displaystyle=$
$\displaystyle\langle E_{ij},\mathcal{W}^{Q}(E_{ij})\rangle=\langle
E_{ij},E_{ij}-\Sigma^{\psi}\circ
E_{ij}-\mathcal{P}^{*}\Sigma^{\varphi}\mathcal{P}(E_{ij}-\Sigma^{\psi}\circ
E_{ij})\rangle$ (22) $\displaystyle=$ $\displaystyle
1-\Sigma^{\psi}_{ij}-\langle
E_{ij},\mathcal{P}^{*}\Sigma^{\varphi}\mathcal{P}(E_{ij}-\Sigma^{\psi}\circ
E_{ij})\rangle.$

In (20), we note that to compute all the diagonal elements of
$\textbf{T}^{\phi}$, the computational cost of the terms $\\{\langle
v^{(ij)},Mv^{(ij)}\rangle\\}_{i,j}$ is $O(p^{4})$ flops. Fortunately, the term
$((P\circ P)M(P\circ P)^{T})_{ij}$ is a very good approximation of ${\bf
T}^{\phi}_{(ij),(ij)}$, whose computational cost is only $O(p^{3})$ flops for
all $i,j$. Similarly, in (22), the computational cost of the terms $\\{\langle
E_{ij},\mathcal{P}^{*}\Sigma^{\varphi}\mathcal{P}(E_{ij}-\Sigma^{\psi}\circ
E_{ij})\rangle\\}_{i,j}$ is also $O(p^{4})$ flops, and $1-\Sigma^{\psi}_{ij}$
is a good approximation to $\textbf{W}^{R}_{(ij),(ij)}$. Hence, we propose a
preconditioner $\mathcal{V}_{D}:\mathbb{S}^{p}\rightarrow\mathbb{S}^{p}$, such
that

$\displaystyle\mathcal{V}_{D}(Y)[H]=D\circ H,$

where $D\in\mathbb{S}^{p}$ and $D_{ij}=((P\circ P)M(P\circ
P)^{T})_{ij}+1-U_{ij}+4(1-\Sigma^{\psi}_{ij})$, for $1\leq i,j\leq p$.

##  5 Numerical experiments

In this section, we demonstrate the performance of our two-phase algorithm.
The pADMM and dADMM serve as benchmark algorithms. We evaluate the efficacy
and efficiency of our algorithm on both synthetic and real-world data. The
algorithms are implemented in Matlab R2019a. All the numerical experiments in
this paper are carried out on a laptop computer based on 64-bit $\rm
Windows$10 system. The computer is configured as: $\rm Intel(R)$$\rm Core(TM)\
i5-6200U\ CPU\ @2.30GHz$$\rm 2.40GHz$, 4G running memory.

###  5.1 Experimental setup

####  5.1.1 Stopping criterion

In our experiments, we define the $R_{P}$, $R_{D}$, and $R_{C}$ to quantify
the infeasibilities of the primal and dual problems, and the complementarity
condition, where

$\displaystyle R_{P}$ $\displaystyle:=\frac{\|\Theta-Z-
V-V^{T}\|}{1+\|\Theta\|},\quad R_{D}:=\frac{\|S-\Omega+Y\|}{1+\|S\|},$
$\displaystyle R_{C}$
$\displaystyle:=\max\left\\{\frac{\|\Theta\Omega-I\|}{1+\|\Theta\|+\|\Omega\|},\frac{\|Z-\textrm{prox}_{Q}(Y+Z)\|}{1+\|Z\|},\frac{\|V-\textrm{prox}_{R}(V+\mathcal{T}^{*}Y)\|}{1+\|V\|}\right\\}.$

We stop the algorithm when $R_{P}$, $R_{D}$, and $R_{C}$ achieve a sufficient
precision:

$\max\left\\{R_{P},R_{D},R_{C}\right\\}<\texttt{Tol}=10^{-6}.$

We also use these three precision metrics to determine when to start Phase II
of the algorithm. Particularly, we switch to phase II if
$\max\\{R_{P},R_{D},R_{C}\\}<10^{-4}$. Moreover, we fix the maximum iteration
numbers of the algorithms. Specifically, in our two-phase algorithm, the Phase
I or Phase II algorithm is terminated if the number of the iterations reaches
200. In the benchmark algorithms, we confine the maximum number of iterations
to 10,000.

####  5.1.2 Performance measures

Given examples
$\textbf{x}_{1},\dots,\textbf{x}_{n}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}\mathcal{N}(\textbf{0},\Sigma)$
with $p$ features, the goal of the DHGL is to estimate the inverse of the
matrix $\Sigma$. Let $\widehat{\Theta}$ be the final estimation of the inverse
empirical covariance matrix, we directly regard the entry of
$\widehat{\Theta}$ whose value is less than a threshold $\epsilon$ as $0$,
such that we can analyze its sparsity pattern. In our experiments, the
threshold $\epsilon$ is fixed to be $10^{-5}$. We refer to the index set of
the estimated hub nodes by $\widehat{\mathcal{H}}_{r}$, and the rule of
identifying a node as a hub is that it contains more than $r$ edges, i.e.,
$\widehat{\mathcal{H}}_{r}=\big{\\{}i\mid\sum_{j=1,j\neq
i}^{p}(1_{\\{|\Omega_{ij}|>\epsilon\\}})>r\big{\\}}$. In our experiments, we
simply set $r$ to be $p/5$.

(1) On synthetic data. The true covariance matrix $\Sigma$ is known. Thus the
true hubs in the graph is available, and we denote the index set of the true
hub nodes by $\mathcal{H}$. Hence, we can use the following four most commonly
used measurements to evaluate the efficacy of the algorithms.

  * (a)

Correctly estimated number of edges:

$\sum_{j<j^{\prime}}(1_{\\{|\Theta_{jj^{\prime}}|>\epsilon~{}\mathrm{and}~{}|\Sigma^{-1}_{jj^{\prime}}|\neq
0\\}}).$

  * (b)

Proportion of correctly estimated hub nodes:

$\frac{|\mathcal{H}_{r}\cap\mathcal{H}|}{|\mathcal{H}|}.$

  * (c)

Proportion of correctly estimated hub edges:

$\frac{\sum_{j\in\mathcal{H},j^{\prime}\neq
j}\left(1_{\left\\{|\Theta_{jj^{\prime}}|>\epsilon~{}\mathrm{and}~{}|\Sigma^{-1}_{jj^{\prime}}|\neq
0\right\\}}\right)}{\sum_{j\in\mathcal{H},j^{\prime}\neq
j}\left(1_{\left\\{|\Sigma^{-1}_{jj^{\prime}}|\neq 0\right\\}}\right)}.$

  * (d)

The mean squared error:

$\mathcal{M}_{\rm
MSE}:=\frac{\sum_{j<j^{\prime}}{\left(\Theta_{jj^{\prime}}-\Sigma^{-1}_{jj^{\prime}}\right)}^{2}}{p(p-1)}.$

(2) On real-world data. In the case of real data, the actual covariance matrix
$\Sigma$ is unknown and cannot be evaluated using the previous four
measurements, so we only compare the computatational time for different
algorithms.

####  5.1.3 Hyperparameter tuning

Similar to other machine learning models, the hyperparameter tuning also plays
a critical role in the DHGL. Tan et al. [37] proposed to select the
hyperparameters by optimizing a Bayesian information criterion (BIC)-type
quantity, but this quantity tends to yield unnecessarily dense graphs when the
dimension $p$ is large [11], and it also introduces two additional
hyperparamters that are too sensitive to the quality of the solution. As a
result, instead of using the BIC-type quantity, we adopt the conventional
k-fold cross-validation [4], and set $k=5$ in our experiments.

For each hyperparameter, we set different search scopes, respectively. In our
experiments, we fix $\lambda_{1}=0.4$, and consider the cases where
$\lambda_{3}\in\\{1,1.5,2\\}$ and $\lambda_{5}\in\\{0.5,1\\}$. For the tuning
of $\lambda_{2}$ and $\lambda_{4}$, we use a relatively finer grid of their
values. For details, we select $\lambda_{2}$ from $[0.1,0.5]$ and
$\lambda_{4}$ from $[0.05,0.15]$.

Another problem in the hyperparameter tuning is to specify the index set
$\mathcal{D}$, which indicates the possible hubs. Li, Bai, and Zhou [25]
recommended employing the GL or HGL before applying the DHGL to identify the
hubs previously, and designed two individual strategies towards various cases
where the hubs are known or unknown. For completeness, we include their
proposed strategies in Algorithms 4 and 5. Furthermore, we note that the
proposed strategies also alleviate the difficulty in hyperparameter selections
because the setting of the hyperparameters in the HGL can be directly adapted
to the DHGL.

Algorithm 4 DHGL with known hubs

1: Using the HGL with cross-validation to obtain the estimated hubs, denoted
by $\mathcal{H}_{\rm HGL}$.

2: Set $\mathcal{D}=\mathcal{K}\setminus\mathcal{H}_{\rm HGL}$, where
$\mathcal{K}$ is the set of known hubs.

3: If $\mathcal{D}\neq\emptyset$, get the estimation $\Theta$ by solving the
DHGL, and denote the estimated hubs by $\mathcal{H}_{\rm DHGL}$, where
$\lambda_{1},\lambda_{2},\lambda_{3}$ are directly set by the same values as
those in the HGL and $\lambda_{4},\lambda_{5}$ are selected based on the
cross-validation. Then let $\mathcal{H}_{r}=\mathcal{H}_{\rm
HGL}\cup\mathcal{H}_{\rm DHGL}$ as the set of the estimated hubs. If
$\mathcal{D}=\emptyset$, use the estimation $\Theta$ of the HGL as the final
estimation result, and let $\mathcal{H}_{r}=\mathcal{H}_{\rm HGL}$.

Algorithm 5 DHGL with unknown hubs

1: Using the HGL with cross-validation to obtain the estimated hubs, denoted
by $\mathcal{H}_{\rm HGL}$.

2: (Prior Information Construction) Tune the regularization parameter
$\lambda$ of the GL from large to small until $|\mathcal{H}_{\rm
GL}\setminus\mathcal{H}_{\rm HGL}|>0$ and $|\mathcal{H}_{\rm
GL}\cup\mathcal{H}_{\rm HGL}|\leq\max\\{|\mathcal{H}_{\rm
HGL}+a,b|\mathcal{H}_{\rm HGL}|\\}$, where $a,b$ are typically set to be $2$
and $1.1$.

3: Set $\mathcal{D}=\mathcal{H}_{\rm GL}\setminus\mathcal{H}_{\rm HGL}$ which
is nonempty.

4: Use the DHGL to estimate $\Theta$, where
$\lambda_{1},\lambda_{2},\lambda_{3}$ remain the same values as those in the
HGL, $\lambda_{4}$ is equally set to be $\lambda_{2}$, and $\lambda_{5}$ is
selected using the cross-validation.

###  5.2 Experiments on synthetic data

We consider two cases (the hubs are known or unknown) and apply the
corresponding strategies. We apply three different algorithms to solve the
core optimization problems (i.e., HGL and DHGL).

Now, we introduce three simulation settings used for the evaluation. We first
randomly generate a $p\times p$ adjacency matrix on the basis of the following
three set-ups.

  * I -

Network with hub nodes. We randomly select $|\mathcal{H}|$ nodes as the hubs.
For any $i<j$, we set the probability of $A_{ij}=1$ to be $0.7$ if the
corresponding node is the hub and to be $0.02$ otherwise. That is, for any
$i<j$, we have

$\mathrm{Prob}(A_{ij}=1)=\left\\{\begin{array}[]{ll}0.7,&\text{if}~{}i\in\mathcal{H},\\\
0.02,&\text{otherwise}.\end{array}\right.$

Finally, the adjacency matrix $A$ is completed by setting $A_{ji}=A_{ij}$ for
$i=1,\dots,p$.

  * II -

Network consisting of two connected subnetworks with hub nodes. The adjacency
matrix $A=\begin{pmatrix}A_{1}&0\\\ 0&A_{2}\end{pmatrix}$, where $A_{1}$ and
$A_{2}$ are generated in accordance with the set-up I.

  * III -

Scale-free network. For any given node, the probability that it has $k$ edges
is proportional to $k^{-\alpha}$. For this simulation setting, we use the
graph library NetworkX [16] to generate the adjacency matrix.

Next, we use the adjacency matrix $A$ to create a matrix $E$. To be specific,
for any $(i,j)$ such that $A_{ij}\neq 0$, let $E_{ij}$ follow a uniform
distribution with the interval $[-0.75,-0.25]\cup[0.25,0.75]$, i.e., in other
words,

$E_{ij}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\left\\{\begin{array}[]{ll}0,&\text{if}~{}A_{ij}=0,\\\
\operatorname{Unif}([-0.75,-0.25]\cup[0.25,0.75]),&\text{otherwise.}\end{array}\right.$

Then, we can obtain a symmetric matrix $\overline{E}$ by setting
$\overline{E}=(E+E^{T})/2$. The inverse covariance matrix is constructed by
setting $\Sigma^{-1}=\overline{E}+(0.1-\rho_{\min}(\overline{E}))I$, where
$\rho_{\min}(\overline{E})$ indicates the smallest eigenvalue of
$\overline{E}$. Now, we can generate the data matrix
$X=(\textbf{x}_{1};\dots;\textbf{x}_{n})$, where each row is sampled from the
normal distribution $\mathcal{N}(0,\Sigma)$. Finally, we conduct an additional
standardization, such that the standard deviation of each feature is equal to
one. To clearly illustrate the generated inverse covariance matrix
$\Sigma^{-1}$, we also give an example in Figure 1 for each set-up.

(a) Set-up I (b) Set-up II (c) Set-up III

Figure 1: Examples of the inverse covariance matrix in three set-ups.

For the generated data matrix $X\in\mathbb{R}^{n\times p}$, we considered the
following seven cases:

$\displaystyle(n,p)$ $\displaystyle\in$
$\displaystyle\big{\\{}(100,300),(300,500),(500,800),(800,1000),(1000,1500),(1500,2000),$
$\displaystyle(2000,2500)\big{\\}}.$

We also need to set the number of the true hub nodes (i.e., the cardinality of
the hub index set $|\mathcal{H}|$) in set-up I and II. In our experiments, we
set $|\mathcal{H}|$ to be $5$, $5$, $10$, $10$, $30$, $30$ and $30$,
respectively for each case. Besides, the cardinality $|\mathcal{H}|$ in set-up
III is dependent on the generated adjacency matrix rather than our setting.
In this simulation setting, a node will be identified as a hub if it has more
than $r$ edges, i.e., $\mathcal{H}=\\{i\mid\sum_{j=1,j\neq
i}^{p}A_{ij}>r,i=1,\dots,p\\}$, and we also set $r=p/5$.

####  5.2.1 Experiments with known hub nodes

When the hub nodes are known, for the DHGL problem, the numerical performances
of the pADMM, dADMM and two-phase algorithms on the synthetic data of
simulation settings I, II and III are compared. The specific numerical
performances are listed in the tables below, where “P”, “D”, “T”, “T(I)” and
“T(II)” denote the pADMM, dADMM, two-phase algorithm, Phase I algorithm and
Phase II algorithm, respectively, and the computational time is recorded in
the format of “hours:minutes:seconds”. As can be seen from the numerical test
results in Tables 1-3, for the synthetic data of known hubs, the two-phase
algorithm is more efficient than the pADMM and dADMM in all the cases.

Table 1 In simulation setting I, the test results for the DHGL problem in
which the hub nodes are known  $(n,p)$ max{$R_{P},R_{D},R_{C}$} iter time P D
T P D T(I) T(II) P D T $(100,300)$ 9.84e-07 9.90e-07 1.29e-07 471 594 200
12(19) 0:00:11 00:00:18 00:00:10 $(300,500)$ 9.34e-07 9.89e-07 1.14e-07 518
329 126 7(8) 00:00:53 00:00:39 00:00:19 $(500,800)$ 9.33e-07 9.93e-07 1.91e-07
781 339 135 5(6) 00:04:13 00:02:25 00:01:09 $(800,1000)$ 9.81e-07 9.70e-07
9.80e-07 955 366 126 4(5) 00:08:42 00:06:02 00:01:35 $(1000,1500)$ 9.90e-07
9.94e-07 8.87e-07 1378 900 200 37(45) 00:33:21 00:27:34 00:19:28 $(1500,2000)$
1.00e-06 9.83e-07 5.01e-07 1609 646 200 12(13) 01:19:47 00:39:53 00:18:16
$(2000,2500)$ 9.99e-07 9.96e-07 5.12e-07 2857 1693 200 31(44) 04:07:22
03:09:18 01:06:12

Table 2 In simulation setting II, the test results for the DHGL problem in
which the hub nodes are known  $(n,p)$ max{$R_{P},R_{D},R_{C}$} iter time P D
T P D T(I) T(II) P D T $(100,300)$ 9.99e-07 9.89e-07 5.62e-07 490 317 200
12(16) 00:00:12 00:00:10 00:00:07 $(300,500)$ 9.96e-07 9.84e-07 8.21e-07 683
372 200 12(14) 00:01:15 00:00:44 00:00:31 $(500,800)$ 9.99e-07 9.89e-07
3.85e-07 781 484 200 13(22) 00:04:16 00:03:11 00:02:15 $(800,1000)$ 9.96e-07
9.84e-07 7.19e-07 1056 551 200 13(20) 00:09:30 00:06:11 00:03:41 $(1000,1500)$
9.95e-07 9.91e-07 4.71e-07 1293 908 200 15(23) 00:31:02 00:27:12 00:11:21
$(1500,2000)$ 9.98e-07 9.96e-07 7.95e-07 2012 1112 200 15(16) 01:50:30
01:09:27 00:17:12 $(2000,2500)$ 9.98e-07 9.99e-07 5.13e-07 1764 1367 200
18(34) 02:55:01 02:36:06 01:04:21

Table 3 In simulation setting III, the test results for the DHGL problem in
which the hub nodes are known  $(n,p)$ max{$R_{P},R_{D},R_{C}$} iter time P D
T P D T(I) T(II) P D T $(100,300)$ 9.89e-07 9.77e-07 3.59e-07 449 338 200
13(16) 00:00:10 00:00:11 00:00:09 $(300,500)$ 9.99e-07 9.82e-07 2.67e-07 714
359 200 12(13) 00:01:15 00:00:49 00:00:35 $(500,800)$ 9.93e-07 9.90e-07
5.01e-07 688 505 200 12(12) 00:03:43 00:03:54 00:01:52 $(800,1000)$ 9.96e-07
9.84e-07 3.60e-07 1058 546 200 13(13) 00:09:52 00:09:19 00:03:13 $(1000,1500)$
9.93e-07 9.97e-07 3.13e-07 778 606 200 14(20) 00:20:18 00:21:09 00:09:50
$(1500,2000)$ 9.98e-07 9.96e-07 4.09e-07 1401 968 200 14(14) 01:15:24 01:22:27
00:21:02 $(2000,2500)$ 9.95e-07 9.94e-07 5.07e-07 2609 1402 200 17(21)
04:30:37 03:48:58 01:11:18

####  5.2.2 Results with unknown hub nodes

When the hub nodes are unknown, for the DHGL problem, the numerical
performances of the pADMM, dADMM and two-phase algorithms on the synthetic
data of simulation Settings I, II and III are compared. The specific
numerical performances are shown in the tables below. As can be seen from the
numerical test results in Tables 4-6, for the synthetic data of unknown hubs,
the two-phase algorithm is the most efficient algorithm in all the cases.

Table 4 In simulation setting I, the test results for the DHGL problem of the
hub node are unknown  $(n,p)$ max{$R_{P},R_{D},R_{C}$} iter time P D T P D
T(I) T(II) P D T $(100,300)$ 9.92e-07 9.98e-07 2.34e-07 737 483 200 12(21)
00:00:19 00:00:15 00:00:11 $(300,500)$ 9.93e-07 9.92e-07 5.44e-07 1038 636 200
12(12) 00:01:50 00:01:22 00:00:31 $(500,800)$ 9.98e-07 9.92e-07 3.94e-07 1443
954 200 14(14) 00:07:55 00:06:55 00:01:50 $(800,1000)$ 9.96e-07 8.96e-07
7.50e-07 1610 976 200 14(14) 00:15:02 00:12:10 00:04:17 $(1000,1500)$ 1.00e-06
9.95e-07 2.52e-07 2238 1579 200 19(24) 00:55:42 00:51:59 00:12:12
$(1500,2000)$ 9.98e-07 9.95e-07 4.43e-07 2784 2003 200 17(17) 02:27:29
02:51:18 00:21:21 $(2000,2500)$ 9.99e-07 1.00e-06 4.96e-07 3954 2505 200
19(53) 06:31:37 05:45:03 02:23:26

Table 5 In simulation setting II, the test results for the DHGL problem of
the hub node are unknown  $(n,p)$ max{$R_{P},R_{D},R_{C}$} iter time P D T P D
T(I) T(II) P D T $(100,300)$ 9.87e-07 9.93e-07 4.54e-07 658 418 200 11(11)
00:00:16 00:00:12 00:00:08 $(300,500)$ 9.95e-07 9.84e-07 4.40e-07 861 372 200
12(19) 00:01:29 00:00:45 00:00:36 $(500,800)$ 9.92e-07 9.89e-07 2.06e-07 929
484 200 13(13) 00:04:57 00:03:21 00:01:38 $(800,1000)$ 9.97e-07 9.84e-07
3.79e-07 1043 551 200 13(13) 00:12:52 00:09:35 00:02:48 $(1000,1500)$ 9.94e-07
9.91e-07 2.51e-07 1531 908 200 15(15) 00:39:33 00:32:09 00:07:34 $(1500,2000)$
9.96e-07 9.96e-07 4.25e-07 1860 1112 200 15(15) 01:43:56 01:15:23 00:15:25
$(2000,2500)$ 9.98e-07 9.99e-07 2.77e-07 2446 1367 200 18(18) 04:01:29
02:59:30 00:34:40

Table 6 In simulation setting III, the test results for the DHGL problem of
the hub node are unknown  $(n,p)$ max{$R_{P},R_{D},R_{C}$} iter time P D T P D
T(I) T(II) P D T $(100,300)$ 9.90e-07 9.99e-07 5.18e-07 438 259 200 10(10)
00:00:10 00:00:09 00:00:08 $(300,500)$ 9.92e-07 9.95e-07 4.67e-07 591 380 200
10(16) 00:01:01 00:00:53 00:00:42 $(500,800)$ 9.92e-07 9.83e-07 2.34e-07 1051
485 200 12(12) 00:05:39 00:04:18 00:01:53 $(800,1000)$ 9.96e-07 9.94e-07
3.42e-07 959 557 200 12(12) 00:08:53 00:06:57 00:03:29 $(1000,1500)$ 9.97e-07
9.96e-07 3.18e-07 1341 937 200 14(14) 00:35:25 00:39:04 00:09:12 $(1500,2000)$
9.98e-07 9.92e-07 4.34e-07 1470 853 200 15(15) 01:22:28 01:14:32 00:22:33
$(2000,2500)$ 9.98e-07 9.89e-07 2.80e-07 2596 793 200 16(16) 04:06:21 01:30:01
00:29:11

In order to explain the graph recovery of the three algorithms more
intuitively, we select an example with dimension $p=800$ for the simulation
setting II , and draw the estimated adjacency matrix.

Figure 2: The adjacency matrix of graph estimation by using three algorithms.
The four subfigures are the true adjacency matrix (i.e., ground truth), the
estimation results of the pADMM, dADMM, and two-phase algorithm, respectively.
The number of the non-zero elements is denoted by nz.

It can be seen from Figure 2 that the estimations via the pADMM and dADMM
contain a large number of wrong edges. However, the two-phase algorithm still
converges well and accurately identifies all the hub nodes.

Meanwhile, in order to illustrate the efficacy of the two-phase algorithm, we
draw the relevant data of four criteria (the same as those in [37]) to
evaluate the efficacy of the algorithm in seven different dimensions for the
simulation setting I when the hub nodes are known. As can be seen from Figure
3, the two-phase algorithm performs well in the correct estimation of the
number of edges, the proportion of correctly estimated hub nodes, the
proportion of correctly estimated hub edges, and the mean squared error, which
also fully demonstrate the efficacy of the two-phase algorithm.

Figure 3: Results of synthetic data efficacy measures

###  5.3 Experiments on real-world data

For the DHGL problem, the numerical performances of the pADMM, dADMM and two-
phase algorithm on the following five real data are compared. We list the
numerical performance in the table below. As can be seen from the numerical
test results in Table 7, the two-phase algorithm can obtain the desired
solutions with the given accuracy for all data sets. In contrast, for most
cases, the pADMM and dADMM reach the maximal number of iterations and cannot
obtain the corresponding results with the given accuracy. Meanwhile, the
computational time of the two-phase algorithm is much less than that of the
pADMM and dADMM.

Table 7 For real data, the test results of the DHGL problem name$(n,p)$
max{$R_{P},R_{D},R_{C}$} iter time P D T P D T(I) T(II) P D T
Iconix$(10455,255)$ 3.27e-03 1.56e-02 5.78e-07 10000 10000 200 28(44) 00:04:00
00:03:31 00:00:10 Rosetta$(6316,301)$ 6.35e-03 2.88e-02 6.67e-07 10000 10000
200 29(35) 00:04:22 00:05:05 00:00:13 animals$(33,102)$ 5.12e-05 9.99e-07
9.43e-07 10000 6545 200 14(14) 00:00:42 00:00:33 00:00:02 cancer$(20531,801)$
6.11e-03 9.74e-02 9.23e-07 10000 10000 200 25(45) 00:59:37 01:07:35 00:02:41
zoo$(100,17)$ 9.97e-07 1.00e-06 5.39e-07 2935 1347 200 14(14) 00:00:03
00:00:02 00:00:01

##  6 Conclusion

In this paper, we have developed a two-phase algorithm to solve the hub
graphical Lasso model with the structured sparsity. Specifically, we design
the dADMM in Phase I to generate a good initial point to warm start Phase II
of the ALM. The SSN method is applied to solve the inner subproblems of the
ALM. We take full advantage of the sparsity structure of the generalized
Jacobian and make the performances of the SSN and ALM efficiently. Numerical
experiments on both synthetic data and real data have demonstrated the
efficacy and efficiency of the proposed algorithm.

##  Appendix A Proof of Proposition 2.4

###### Proof.

Firstly, we calculate the proximal operators of $\varphi$ and $\psi$. Since
for any $V\in\mathbb{M}^{p}$,

$\varphi(V)=\sum_{j=1}^{p}w_{2,j}\|\mathcal{P}_{j}V\|,$

by denoting
$\widehat{\mathcal{P}}_{j}:\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}$ as
$\widehat{\mathcal{P}}_{j}y=[y_{1};\cdots;y_{j-1};0;y_{j+1};\cdots;y_{p}]$ for
any $y\in\mathbb{R}^{p}$, we have that

$\displaystyle\varphi^{*}(Y)$ $\displaystyle=$
$\displaystyle\sup_{V\in\mathbb{M}^{p}}\Big{\\{}\langle
Y,V\rangle-\sum_{j=1}^{p}w_{2,j}\|\mathcal{P}_{j}V\|\Big{\\}}$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{p}\sup_{V_{j}\in\mathbb{R}^{p}}\Big{\\{}\langle
Y_{j},V_{j}\rangle-w_{2,j}\|\widehat{\mathcal{P}}_{j}V_{j}\|\Big{\\}}$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{p}\sup_{V_{j}\in\mathbb{R}^{p}}\Big{\\{}Y_{jj}V_{jj}+\langle\mathcal{P}_{j}Y,\widehat{\mathcal{P}}_{j}V_{j}\rangle-
w_{2,j}\|\widehat{\mathcal{P}}_{j}V_{j}\|\Big{\\}}.$

Note that

$\displaystyle
Y_{jj}V_{jj}+\langle\mathcal{P}_{j}Y,\widehat{\mathcal{P}}_{j}V_{j}\rangle-
w_{2,j}\|\widehat{\mathcal{P}}_{j}V_{j}\|\leq
Y_{jj}V_{jj}+\|\widehat{\mathcal{P}}_{j}V_{j}\|(\|\mathcal{P}_{j}Y\|-w_{2,j}).$

Therefore,

$\displaystyle\varphi^{*}(Y)$ $\displaystyle=$ $\displaystyle\delta_{C}(Y),$

where $C:=\Big{\\{}Y\in\mathbb{M}^{p}\,\Big{|}\,\|\mathcal{P}_{j}Y\|\leq
w_{2,j}\textrm{ and }Y_{jj}=0,j=1,2,\cdots,p\Big{\\}}$. Then by the Moreau
identity, we have

$\displaystyle\operatorname{prox}_{\varphi}(Y)=Y-\operatorname{prox}_{\varphi^{*}}(Y)=Y-\Pi_{C}(Y)=Y-\mathcal{P}^{*}\Pi_{\mathcal{B}^{w_{2}}_{2}}(\mathcal{P}Y),$

where

$[\Pi_{\mathcal{B}^{w_{2}}_{2}}(\mathcal{P}Y)]_{(j-1)p+1:jp}=\Pi_{\mathcal{B}_{2}^{w_{2,j}}}(\mathcal{P}_{j}Y)=\left\\{\begin{array}[]{cl}w_{2,j}\frac{\mathcal{P}_{j}Y}{\|\mathcal{P}_{j}Y\|}&\text{if}~{}\|\mathcal{P}_{j}Y\|>w_{2,j},\\\
\mathcal{P}_{j}Y&\text{otherwise.}\end{array}\right.$ (23)

Similarly, we can calculate that

$\displaystyle\operatorname{prox}_{\psi}(Y)=Y-\mathcal{P}^{*}\Pi_{\mathcal{B}^{w_{1}}_{\infty}}(\mathcal{P}Y),$

where

$\displaystyle[\Pi_{\mathcal{B}^{w_{1}}_{\infty}}(\mathcal{P}Y)]_{(j-1)p+1:jp}$
$\displaystyle=\Pi_{\mathcal{B}_{\infty}^{w_{1,j}}}(\mathcal{P}_{j}Y)=\operatorname{sign}(\mathcal{P}_{j}Y)\circ\min(|\mathcal{P}_{j}Y|,w_{1,j}).$

In the following, we prove the result about the proximal mapping of
$R(\cdot)$. We can see from (4) that the function $R(V)$ has a separable
structure, thus we only need to prove that

$\displaystyle\operatorname{prox}_{w_{1,j}\|\mathcal{P}_{j}\cdot\|_{1}+w_{2,j}\|\mathcal{P}_{j}\cdot\|}(Y)=\operatorname{prox}_{w_{2,j}\|\mathcal{P}_{j}\cdot\|}\circ\operatorname{prox}_{w_{1,j}\|\mathcal{P}_{j}\cdot\|_{1}}(Y),\
j=1,2,\cdots,p.$

By [42, Theorem 1], it suffices to show that for each $j=1,\cdots,p$,

$\displaystyle\partial(w_{1,j}\|\mathcal{P}_{j}Y\|_{1})\subseteq\partial(w_{1,j}\|\mathcal{P}_{j}Z\|_{1}),$

where $Z={\rm
prox}_{\varphi}(Y)=Y-\mathcal{P}^{*}\Pi_{\mathcal{B}_{2}^{w_{2}}}(\mathcal{P}Y)$.
We prove this relation in two cases:

(1) If $\|\mathcal{P}_{j}Y\|\leq w_{2,j}$, from (23), we have
$\mathcal{P}_{j}Z=0$. Then according to [32, Theorem 23.9], we obtain

$\displaystyle\partial(w_{1,j}\|\mathcal{P}_{j}Z\|_{1})=\Big{\\{}\mathcal{P}_{j}^{*}v\,\Big{|}\,v\in[-w_{1,j},w_{1,j}]^{p}\Big{\\}},$

which contains $\partial(w_{1,j}\|\mathcal{P}_{j}Y\|_{1})$.

(2) If $\|\mathcal{P}_{j}Y\|>w_{2,j}$, from (23), we have
$\mathcal{P}_{j}Z=(1-\frac{w_{2,j}}{\|\mathcal{P}_{j}Y\|})(\mathcal{P}_{j}Y)$,
which implies
$\textrm{sign}(\mathcal{P}_{j}Y)=\textrm{sign}(\mathcal{P}_{j}Z)$. Therefore,
$\partial(w_{1,j}\|\mathcal{P}_{j}Y\|_{1})\subseteq\partial(w_{1,j}\|\mathcal{P}_{j}Z\|_{1})$.

Now we have finished all the proof. ∎

## References

  * [1] R. Albert, Scale-free networks in cell biology, Journal of Cell Science, 118 (2005), pp. 4947–4957.
  * [2] O. Banerjee, L. El Ghaoui, and A. d’Aspremont, Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data, Journal of Machine Learning Research, 9 (2008), pp. 485–516.
  * [3] A.-L. Barabasi and Z. N. Oltvai, Network biology: Understanding the cell’s functional organization, Nature Reviews Genetics, 5 (2004), pp. 101–113.
  * [4] P. J. Bickel and E. Levina, Regularized estimation of large covariance matrices, The Annals of Statistics, 36 (2008), pp. 199–227.
  * [5] L. Bybee and Y. Atchadé, Change-point computation for large graphical models: A scalable algorithm for Gaussian graphical models with change-points, Journal of Machine Learning Research, 19 (2018), pp. 1–38.
  * [6] A. Chaturvedi and J. Scarlett, Learning Gaussian graphical models via multiplicative weights, in International Conference on Artificial Intelligence and Statistics, PMLR, 2020, pp. 1104–1114.
  * [7] F. Clarke, Optimization and Nonsmooth Analysis, John Wiley and Sons, New York, 1983.
  * [8] D. Cox and N. Wermuth, Multivariate Dependencies: Models, Analysis and Interpretation, Chapman & Hall/CRC Monographs on Statistics & Applied Probability, CRC Press, 2014.
  * [9] A. d’Aspremont, O. Banerjee, and L. El Ghaoui, First-order methods for sparse covariance selection, SIAM Journal on Matrix Analysis and Applications, 30 (2008), pp. 56–66.
  * [10] A. Defazio and T. Caetano, A convex formulation for learning scale-free networks via submodular relaxation, in Advances in Neural Information Processing Systems, vol. 25, Curran Associates, Inc., 2012.
  * [11] M. Drton and M. H. Maathuis, Structure learning in graphical modeling, Annual Review of Statistics and Its Application, 4 (2017), pp. 365–393.
  * [12] F. Facchinei and J.-S. Pang, Finite-Dimensional Variational Inequalities and Complementarity Problems, vol. II, Springer, 2003.
  * [13] J. Friedman, T. Hastie, and R. Tibshirani, Sparse inverse covariance estimation with the graphical lasso, Biostatistics, 9 (2008), pp. 432–441.
  * [14] X. Gao, W. Shen, C.-M. Ting, S. C. Cramer, R. Srinivasan, and H. Ombao, Modeling brain connectivity with graphical models on frequency domain, arXiv preprint arXiv:1810.03279, (2018).
  * [15] N. Hafiene, W. Karoui, and L. B. ROMDHANE, Influential nodes detection in dynamic social networks: A survey, Expert Systems with Applications, (2020), p. 113642.
  * [16] A. A. Hagberg, D. A. Schult, and P. J. Swart, Exploring network structure, dynamics, and function using networkx, in Proceedings of the 7th Python in Science Conference, Pasadena, CA USA, 2008, pp. 11–15.
  * [17] D. Han, D. F. Sun, and L. Zhang, Linear rate convergence of the alternating direction method of multipliers for convex composite programming, Mathematics of Operations Research, 43 (2018), pp. 622–637.
  * [18] I. Himelboim, M. A. Smith, L. Rainie, B. Shneiderman, and C. Espina, Classifying twitter topic-networks using social network analysis, Social Media+ Society, 3 (2017), p. 2056305117691545.
  * [19] J.-B. Hiriart-Urruty, J.-J. Strodiot, and V. H. Nguyen, Generalized Hessian matrix and second-order optimality conditions for problems with $c^{1,1}$ data, Applied Mathematics and Optimization, 11 (1984), pp. 43–56.
  * [20] J. Honorio, D. Samaras, N. Paragios, R. Goldstein, and L. E. Ortiz, Sparse and locally constant Gaussian graphical models, in Advances in Neural Information Processing Systems, 2009, pp. 745–753.
  * [21] M. J. Hosseini and S.-I. Lee, Learning sparse Gaussian graphical models with overlapping blocks, in Advances in Neural Information Processing Systems, vol. 29, Curran Associates, Inc., 2016, pp. 3808–3816.
  * [22] C.-J. Hsieh, M. A. Sustik, I. S. Dhillon, and P. Ravikumar, Sparse inverse covariance matrix estimation using quadratic approximation, in Advances in Neural Information Processing Systems, vol. 24, Curran Associates Inc., 2011, pp. 2330–2338.
  * [23] C.-J. Hsieh, M. A. Sustik, I. S. Dhillon, and P. Ravikumar, QUIC: Quadratic approximation for sparse inverse covariance estimation, Journal of Machine Learning Research, 15 (2014), pp. 2911–2947.
  * [24] L. Li and K.-C. Toh, An inexact interior point method for $\ell_{1}$-regularized sparse covariance selection, Mathematical Programming Computation, 2 (2010), pp. 291–315.
  * [25] Z. Li, J. Bai, and W. Zhou, Learning Gaussian graphical models using discriminated hub graphical lasso, in 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE, 2018, pp. 2471–2475.
  * [26] Z. Lu, Smooth optimization approach for sparse covariance selection, SIAM Journal on Optimization, 19 (2009), pp. 1807–1827.
  * [27] Z. Meng, B. Eriksson, and A. Hero, Learning latent variable Gaussian graphical models, in Proceedings of the 31st International Conference on Machine Learning, vol. 32, PMLR, 2014, pp. 1269–1277.
  * [28] A. J. Molstad and A. J. Rothman, Shrinking characteristics of precision matrix estimators, Biometrika, 105 (2018), pp. 563–574.
  * [29] T. Nakagaki, M. Fukuda, S. Kim, and M. Yamashita, A dual spectral projected gradient method for log-determinant semidefinite problems, Computational Optimization and Applications, (2020), pp. 33–68.
  * [30] F. Oztoprak, J. Nocedal, S. Rennie, and P. A. Olsen, Newton-like methods for sparse inverse covariance estimation, in Advances in Neural Information Processing Systems, vol. 25, Curran Associates, Inc., 2012, pp. 755–763.
  * [31] P. Ravikumar, G. Raskutti, M. J. Wainwright, and B. Yu, Model selection in Gaussian graphical models: High-dimensional consistency of $\ell_{1}$-regularized MLE, in Proceedings of the 21st International Conference on Neural Information Processing Systems, Curran Associates Inc., 2008, pp. 1329–1336.
  * [32] R. T. Rockafellar, Convex Analysis, Princeton University Press, 1970\.
  * [33] R. T. Rockafellar, Augmented Lagrangians and applications of the proximal point algorithm in convex programming, Mathematics of Operations Research, 1 (1976), pp. 97–116.
  * [34] R. T. Rockafellar and R. J.-B. Wets, Variational Analysis, vol. 317, Springer Science & Business Media, 2009.
  * [35] K. Scheinberg, S. Ma, and D. Goldfarb, Sparse inverse covariance selection via alternating linearization methods, in Advances in Neural Information Processing Systems, vol. 23, Curran Associates, Inc., 2010, pp. 2101–2109.
  * [36] D. W. Soh and S. Tatikonda, Learning unfaithful k-separable Gaussian graphical models., Journal of Machine Learning Research, 20 (2019), pp. 1–30.
  * [37] K. M. Tan, P. London, K. Mohan, S.-I. Lee, M. Fazel, and D. Witten, Learning graphical models with hubs, Journal of Machine Learning Research, 15 (2014), pp. 3297–3331.
  * [38] D. A. Tarzanagh and G. Michailidis, Estimation of graphical models through structured norm minimization, Journal of Machine Learning Research, 18 (2018), pp. 7692–7739.
  * [39] C. Wang, D. F. Sun, and K.-C. Toh, Solving log-determinant optimization problems by a Newton-CG primal proximal point algorithm, SIAM Journal on Optimization, 20 (2010), pp. 2994–3013.
  * [40] C. Wang and P. Tang, A dual semismooth Newton based augmented Lagrangian method for large-scale linearly constrained sparse group square-root lasso problems, Journal of Scientific Computing, (2023).
  * [41] S. Yang, Z. Lu, X. Shen, P. Wonka, and J. Ye, Fused multiple graphical lasso, SIAM Journal on Optimization, 25 (2015), pp. 916–943.
  * [42] Y.-L. Yu, On decomposing the proximal map, in Advances in Neural Information Processing Systems, Curran Associates, Inc., 2013, pp. 91–99.
  * [43] M. Yuan and Y. Lin, Model selection and estimation in the Gaussian graphical model, Biometrika, 94 (2007), pp. 19–35.
  * [44] Y. Zhang, N. Zhang, D. F. Sun, and K.-C. Toh, An efficient Hessian based algorithm for solving large-scale sparse group lasso problems, Mathematical Programming, 179 (2020), pp. 223–263.
  * [45] , A proximal point dual Newton algorithm for solving group graphical lasso problems, SIAM Journal on Optimization, 30 (2020), pp. 2197–2220.