Efficient Proximal Gradient Algorithms for Joint Graphical Lasso

Abstract

We consider learning as an undirected graphical model from sparse data. While several efficient algorithms have been proposed for graphical lasso (GL), the alternating direction method of multipliers (ADMM) is the main approach taken concerning joint graphical lasso (JGL). We propose proximal gradient procedures with and without a backtracking option for the JGL. These procedures are first-order methods and relatively simple, and the subproblems are solved efficiently in closed form. We further show the boundedness for the solution of the JGL problem and the iterates in the algorithms. The numerical results indicate that the proposed algorithms can achieve high accuracy and precision, and their efficiency is competitive with state-of-the-art algorithms.

1. Introduction

Graphical models are widely used to describe the relationships among interacting objects [1]. Such models have been extensively used in various domains, such as bioinformatics, text mining, and social networks. The graph provides a visual way to understand the joint distribution of an entire set of variables.

In this paper, we consider learning Gaussian graphical models that are expressed by undirected graphs, which represent the relationship among continuous variables that follow a joint Gaussian distribution. In an undirected graph, $\mathcal{G}=(V,E)$, and edge set E represents the conditional dependencies among the variables in vertex set V.

Let $X_1,\dots ,X_p$ ($p\ge 1$) be Gaussian variables with covariance matrix $\Sigma \in {\mathbb{R}}^{p\times p}$, and $\Theta :={\Sigma }^{-1}$ be the precision matrix, if it exists. We remove the edges so that the variables $X_i$, $X_j$ are conditionally independent given the other variables if and only if the $(i,j)$-th element ${\theta }_{i,j}$ in $\Theta$ is 0:

$$\{i,j\}\notin E⟺{\theta }_{i,j}=0⟺X_i\perp \perp X_j|X_{V\backslash \{i,j\}},$$

where each edge is expressed as a set of two elements in $\{1,\dots ,p\}$. In this sense, constructing a Gaussian graphical model is equivalent to estimating a precision matrix.

Suppose that we estimate the undirected graph from data consisting of n tuples of p variables and that dimension p is much higher than sample size n. For example, if we have expression data of $p=\mathrm{20,000}$ genes for $n=100$ case/control patients, how can we construct a gene regulatory network structure from the data? It is almost impossible to estimate the locations of the nonzero elements in $\Theta$ by obtaining the inverse of sample covariance matrix $\mathit{S}\in {\mathbb{R}}^{p\times p}$, which is the unbiased estimator of $\Sigma$. In fact, if $p>n$, then no inverse ${\mathit{S}}^{-1}$ exists because the rank of $\mathit{S}\in {\mathbb{R}}^{p\times p}$ is, at most, n.

In order to address this situation, two directions are suggested:

• Sequentially find the variables on which each variable depends via regression so that the quasilikelihood is maximized [2].
• Find the locations in $\Theta$, the values of which are zeros, so that the ${\ell }_1$ regularized log-likelihood is maximized [3,4,5,6].

We follow the second approach because we assume Gaussian variables, also known as graphical lasso (GL). The  ${\ell }_1$ regularized log-likelihood is defined by:

$$\begin{array}{c}\hfill \underset{\Theta }{\mathrm{maximize}}\left\{logdet\Theta -\mathrm{trace}(\mathit{S}\Theta )-{\lambda \left|\right|\Theta \left|\right|}_1\right\},\end{array}$$

(1)

where tuning parameter $\lambda$ controls the amount of sparsity, and  ${\left|\right|\Theta \left|\right|}_1$ denotes the sum of the absolute value of the off-diagonal elements in $\Theta$. Several optimization techniques [4,7,8,9,10,11,12] have been studied for the optimization problem of (1).

In particular, we consider a generalized version of the abovementioned GL. For example, suppose that the gene regulatory networks of thirty case and seventy control patients are different. One might construct a gene regulatory network separately for each of the two categories. However, estimating each on its own does not provide an advantage if a common structure is shared. Instead, we use 100 samples to construct two networks simultaneously. Intuitively speaking, using both types of data improves the reliability of the estimation by increasing the sample size for the genes that show similar values between case and control patients, while using only one type of data leads to a more accurate estimate for genes that show significantly different values. Ref. [13] proposed a joint graphical lasso (JGL) model by including an additional convex penalty (fused or group lasso penalty) to the graphical lasso objective function for K classes. For example, K is equal to two for the case/control patients in the example. JGL includes fused graphical lasso with fused lasso penalty, which encourages sparsity and the similarity of the value of edges across K classes, and group graphical lasso with group lasso penalty, which promotes similar sparsity structure across K graphs. Although there are several approaches to handling the multiple graphical models, such as those of [14,15,16,17], the JGL is considered the most promising.

The main topic of this paper is efficiency improvement in terms of solving the JGL problem. For the GL, relatively efficient solving procedures exist. If we differentiate the ${\ell }_1$ regularized log-likelihood (1) by $\Theta$, then we have an equation to solve [4]. Moreover, several improvements have been considered for the GL, such as proximal Newton [12] and proximal gradient [10] procedures. However, for the JGL, even if we derive such an equation, we have no efficient way of handling it.

Instead, the alternating direction method of multipliers (ADMM) [18], which is a procedure for solving convex optimization problems for general purposes, has been the main approach taken [13,19,20,21]. However, ADMM does not scale well concerning the feature dimension p and number of classes K. It usually takes time to converge to a high-accuracy solution [22].

For efficient procedures to solve the JGL problem, ref. [23] proposed a method based on the proximal Newton method only when the penalty term is expressed by fused lasso. The existing method requires expensive computations for the Hessian matrix and Newton directions, which means that it would be expensive to use for high-dimensional problems.

In this paper, we propose efficient proximal-gradient-based algorithms to solve the JGL problem by extending the procedure in [10] and employing the step-size selection strategy proposed in [24]. Moreover, we provide the theoretical analysis of both methods for the JGL problem.

In our proximal gradient methods for the JGL problem, the proximal operator in each iteration is quite simple, which eases the implementation process and requires very little computation and memory at each step. Simulation experiments are used to justify our proposed methods over the existing ones.

Our main contributions are as follows:

• We propose efficient algorithms based on the proximal gradient method to solve the JGL problem. The algorithms are first-order methods and quite simple, and the subproblems can be solved efficiently with a closed-form solution. The numerical results indicate that the methods can achieve high accuracy and precision, and the computational time is competitive with state-of-art algorithms.
• We provide the boundedness for the solution to the JGL problem and the iterates in algorithms, which is related to the convergence rate of the algorithms. With the boundedness, we can guarantee that our proposed method converges linearly.

Table 1. Efficient JGL procedures.

Model	ADMM	Proximal Newton	Proximal Gradient
GL [4]	[8]	[12]	[10]
JGL [13]	[13]	[23]	Current Paper
		(for fused penalty)	(for fused and group penalties)

summarizes the relationship between the proposed and existing methods.

The remaining parts of this paper are as follows. In Section 2, we first provide the background of our proposed methods and introduce the joint graphical lasso problem. In Section 3, we illustrate the detailed content of the proposed algorithms and provide a theoretical analysis. In Section 4, we report some numerical results of the proposed approaches, including comparisons with efficient methods and performance evaluations. Finally, we draw some conclusions in Section 5.

Notation: In this paper, ${\left|\right|\mathit{x}\left|\right|}_p$ denotes the ${\ell }_p$ norm of a vector $\mathit{x}\in {\mathbb{R}}^d$, ${\left|\right|\mathit{x}\left|\right|}_p:=({\sum }_{i=1}^d|x_i{{|}^p)}^{\frac{1}{p}}$ for $p\in [1,\infty )$, and ${\left|\right|\mathit{x}\left|\right|}_{\infty }:={max}_i\left|x_i\right|$. For a matrix $\mathit{X}\in {\mathbb{R}}^{p\times q}$, ${\left|\right|\mathit{X}\left|\right|}_F$ denotes the Frobenius norm, ${\left|\right|\mathit{X}\left|\right|}_2$ denotes the spectral norm, ${\left|\right|\mathit{X}\left|\right|}_{\infty }:={max}_{i,j}\left|x_{i,j}\right|$, and ${\left|\right|\mathit{X}\left|\right|}_1:={\sum }_{i=1}^p{\sum }_{j=1}^q\left|x_{i,j}\right|$ if not specified. The inner product is defined by $\langle \mathit{X},\mathit{X}\rangle :=\mathrm{trace}\left({\mathit{X}}^T\mathit{X}\right)$.

2. Preliminaries

This section first reviews the graphical lasso (GL) problem and introduces the graphical iterative shrinkage-thresholding algorithm (G-ISTA) [10] to solve it. Then, we introduce the step-size selection strategy that we apply to the joint graphical lasso (JGL) in Section 3.2.

2.1. Graphical Lasso

Let ${\mathit{x}}_1,\dots ,{\mathit{x}}_n\in {\mathbb{R}}^p$ be $n\ge 1$ observations of dimension $p\ge 1$ that follow the Gaussian distribution with mean $\mathit{\mu }\in {\mathbb{R}}^p$ and covariance matrix $\Sigma \in {\mathbb{R}}^{p\times p}$, where without loss of generality, we assume $\mathit{\mu }=\mathbf{0}$. Let $\Theta :={\Sigma }^{-1}$, and the empirical covariance matrix $\mathit{S}:=\frac{1}{n}{\sum }_{i=1}^n{\mathit{x}}_i^T{\mathit{x}}_i$. Given penalty parameter $\lambda >0$, the graphical lasso (GL) is the procedure to find the positive definite $\Theta \in {\mathbb{R}}^{p\times p}$ such that:

$$\underset{\Theta }{\mathrm{minimize}}\left\{-logdet\Theta +\mathrm{trace}\left(\mathit{S}\Theta \right)+{\lambda \parallel \Theta \parallel }_1\right\},$$

(2)

where ${\left|\right|\Theta \left|\right|}_1={\sum }_{j\ne k}\left|{\theta }_{j,k}\right|.$ If we regard $V:=\{1,\cdots ,p\}$ as a vertex set, then we can construct an undirected graph with edge set $\left\{\{j,k\}\right|{\theta }_{j,k}\ne 0\}$, where set $\{j,k\}$ denotes an undirected edge that connects the nodes $j,k\in V$.

If we take the subgradient of (2), then we find that the optimal solution ${\Theta }_{*}$ satisfies the condition:

$$-{\Theta }_{*}^{-1}+\mathit{S}+\lambda \Phi \ni 0,$$

(3)

where $\Phi =\left({\Phi }_{j,k}\right)$ is

$${\Phi }_{j,k}=\left\{\begin{array}{cc}1,\hfill & {\theta }_{j,k}^{*}>0\hfill \\ [-1,1],\hfill & {\theta }_{j,k}^{*}=0\hfill \\ -1,\hfill & {\theta }_{j,k}^{*}<0\hfill \end{array}\right..$$

2.2. ISTA for Graphical Lasso

In this subsection, we introduce the method for solving the GL problem (2) by the iterative shrinkage-thresholding algorithm (ISTA) proposed by [10], which is a proximal gradient method usually employed in dealing with nondifferentiable composite optimization problems.

Specifically, the general ISTA solves the following composite optimization problem:

$$\underset{\mathit{x}}{\mathrm{minimize}}F\left(\mathit{x}\right):=f\left(\mathit{x}\right)+g\left(\mathit{x}\right),$$

(4)

where f and g are convex, with f differentiable and g possibly being nondifferentiable.

For the GL problem (2), we denote $f,g:{\mathbb{R}}^{p\times p}\to \mathbb{R}$ as

$$f(\Theta ):=-logdet\Theta +\mathrm{trace}\left(\mathit{S}\Theta \right),$$

and

$$g(\Theta ):={\lambda \parallel \Theta \parallel }_1.$$

If we define the quadratic approximation $Q_{\eta }:{\mathbb{R}}^{p\times p}\times {\mathbb{R}}^{p\times p}\to \mathbb{R}$ w.r.t. $f(\Theta )$ and $\eta >0$ by

$$\begin{array}{c}\hfill Q_{\eta }({\Theta }^{\prime },\Theta ):=f(\Theta )+\langle {\Theta }^{\prime }-\Theta ,\nabla f(\Theta )\rangle +\frac{1}{2\eta }\left|\right|{\Theta }^{\prime }-\Theta {\left|\right|}_F^2,\end{array}$$

(5)

then we can describe the ISTA as a procedure that iterates

$$\begin{array}{cc}\hfill {\Theta }_{t+1}& =arg\underset{\Theta }{min}\{Q_{{\eta }_t}(\Theta ,{\Theta }_t)+g(\Theta )\}\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & ={\mathrm{prox}}_{{\eta }_{t}g}({\Theta }_t-{\eta }_t\nabla f\left({\Theta }_t\right)),\hfill \end{array}$$

(7)

given initial value ${\Theta }_0$, where the value of step size ${\eta }_t>0$ may change at each iteration $t=1,2,\cdots ,$ for efficient convergence purpose, and we use the proximal operator:

$${\mathrm{prox}}_g\left(\mathit{z}\right):=arg\underset{\mathit{\theta }}{min}\{\frac{1}{2}{\parallel \mathit{z}-\mathit{\theta }\parallel }_2^2+g\left(\mathit{\theta }\right)\}.$$

(8)

Note that the proximal operator of function $g={\lambda \left|\right|\Theta \left|\right|}_1$ is the soft-thresholding operator: the absolute value $|{\theta }_{i,j}|$ of each off-diagonal element ${\theta }_{i,j}$ with $i\ne j$ becoming either ${\theta }_{i,j}-\mathrm{sgn}\left({\theta }_{i,j}\right)\lambda$ or zero (if $|{\theta }_{i,j}|<\lambda$). We use the following function for the operator in Section 3:

$${\left[{\mathcal{S}}_{\lambda }(\Theta )\right]}_{i,j}=\mathrm{sgn}\left({\theta }_{i,j}\right)\left(\right|{\theta }_{i,j}{|-\lambda )}_{+}$$

(9)

where ${\left(x\right)}_{+}:=max(x,0)$.

Definition 1.

A differentiable function $f:{\mathbb{R}}^{n\times p}\to \mathbb{R}$ is said to have a Lipschitz-continuous gradient if there exists $L>0$ (Lipschitz constant) such that

$$\left|\right|\nabla f\left(\mathit{X}\right)-{\nabla f\left(\mathit{Y}\right)\left|\right|}_F\le L\left|\right|\mathit{X}-{\mathit{Y}\left|\right|}_F,\forall \mathit{X},\mathit{Y}\in {\mathbb{R}}^{n\times p}.$$

(10)

It is known that if we choose ${\eta }_t=\frac{1}{L}$ for each step in the ISTA that minimizes $F(·)$, then the convergence rate is, at most, as follows [25]:

$$F\left({\Theta }_t\right)-F\left({\Theta }_{*}\right)=O\left(\frac{1}{t}\right)$$

(11)

However, for the GL problem (2), we know neither the exact value of the Lipschitz constant L nor any nontrivial upper bound. [10] implement a backtracking line search option in Step 1 of Algorithm 1 below to handle this issue.

The backtracking line search enables us to compute the ${\eta }_t$ value for each time $t=1,2,\dots$ by repeatedly multiplying ${\eta }_t$ by a constant $c\in (0,1)$ until ${\Theta }_{t+1}\succ 0$ ($\Theta$ is positive definite) and

$$f\left({\Theta }_{t+1}\right)\le Q_{{\eta }_t}({\Theta }_{t+1},{\Theta }_t),$$

(12)

for the ${\Theta }_{t+1}$ in (7). Additionally, (12) is a sufficient condition for (11), which was derived in [25] (see the relationship between Lemma 2.3 and Theorem 3.1 in [25]).

The whole procedure is given in Algorithm 1.

Algorithm 1 G-ISTA for problem (2).

Input: $\mathit{S}$, tolerance $ϵ>0$, backtracking constant $0<c<1$, initial value ${\eta }_0$, ${\Theta }_0$, $t=0$. While$t<t_{max}$ (until convergence) do

1: Backtracking line search: Continue to multiply ${\eta }_t$ by c until $${\Theta }_{t+1}\succ 0\mathrm{and}f\left({\Theta }_{t+1}\right)\le Q_{{\eta }_t}({\Theta }_{t+1},{\Theta }_t)$$ for ${\Theta }_{t+1}:={prox}_{{\eta }_{t}g}({\Theta }_t-{\eta }_t\nabla f\left({\Theta }_t\right))$.     2: Update iterate: ${\Theta }_{t+1}\leftarrow {prox}_{{\eta }_{t}g}({\Theta }_t-{\eta }_t\nabla f\left({\Theta }_t\right))$.     3: Set next initial step size ${\eta }_{t+1}$ by the Barzilai—Borwein method.     4: $t\leftarrow t+1$

end Output: $ϵ$-optimal solution to problem (2), ${\Theta }_{*}={\Theta }_{t+1}$.

2.3. Composite Self-Concordant Minimization

The notion of the self-concordant function was proposed in [26,27,28]. In the following, we say a convex function f is self-concordant with parameter $M\ge 0$ if

$$|f^{‴}\left(\mathit{x}\right)|\le M{f}^{″}{\left(\mathit{x}\right)}^{3/2},\mathrm{for}\mathrm{all}\mathit{x}\in \mathrm{dom}f.$$

where $\mathrm{dom}f$ is the domain of f.

Reference [24] considered a composite version of self-concordant function minimization and provided a way to efficiently calculate the step size for the proximal gradient method for the GL problem without relying on the Lipschitz gradient assumption in (10). They proved that

$$f(\Theta ):=-logdet\Theta +\mathrm{trace}(S\Theta )$$

in (2) is self-concordant and considers the following minimization:

$$F^{*}:=\underset{\mathit{x}}{\mathrm{minimize}}\{F\left(\mathit{x}\right):=f\left(\mathit{x}\right)+g\left(\mathit{x}\right)\},$$

where f is convex, differentiable, and self-concordant, and g is convex and possibly nondifferentiable. As for Algorithm 1, without using the backtracking line search, we can compute direction ${\mathit{d}}_t$ with initial step size ${\eta }_t$ as follows:

$${\mathit{d}}_t:={prox}_{{\eta }_{t}g}({\Theta }_t-{\eta }_t\nabla f\left({\Theta }_t\right))-{\Theta }_t,$$

(13)

where the operator $prox$ is defined by (8). Then, we use the modified step size ${\alpha }_t$ to update ${\Theta }_{t+1}:={\Theta }_t+{\alpha }_t{\mathit{d}}_t$, which can be determined by the direction ${\mathit{d}}_t$. After defining two parameters related to the direction: ${\beta }_t:={\eta }_t^{-1}\left|\right|{\mathit{d}}_t{\left|\right|}_F^2$ and ${\lambda }_t:={\left(\langle {\nabla }^{2}f\left({\Theta }_t\right){\mathit{d}}_t,{\mathit{d}}_t\rangle \right)}^{1/2}$, the modified step size can be obtained by

$${\alpha }_t:=\frac{{\beta }_t}{{\lambda }_t({\lambda }_t+{\beta }_t)}.$$

(14)

By Lemma 12 in [24], if the modified step size ${\alpha }_t\in (0,1]$, then it can ensure a decrease in the objective function and guarantee convergence in the proximal gradient scheme. From (14), if ${\lambda }_t\ge 1$, then the condition ${\alpha }_t\in (0,1]$ is satisfied. Therefore, we only need to check the case when ${\lambda }_t<1$. If the condition ${\alpha }_t\in (0,1]$ does not hold, we can change the value of the initial ${\eta }_t$ (such as the bisection method) to influence the value of ${\mathit{d}}_t$ in (13) until the condition is satisfied.

2.4. Joint Graphical Lasso

Let $N\ge 1,p\ge 1,K\ge 2$, and $({\mathit{x}}_1,y_1),\dots ,({\mathit{x}}_N,y_N)\in {\mathbb{R}}^p\times \{1,\dots ,K\}$, where each ${\mathit{x}}_i$ is a row vector. Let $n_k$ be the number of occurrences in $y_1,\dots ,y_N$ such that $y_i=k$, so that ${\sum }_{k=1}^K{n}_k=N$.

For each $k=1,\dots ,K$, we define the empirical covariance matrix ${\mathit{S}}^{\left(k\right)}\in {\mathbb{R}}^{p\times p}$ of the data ${\mathit{x}}_i$ as follows:

$${\mathit{S}}^{\left(k\right)}:=\frac{1}{n_k}\sum _{i:y_i=k}{\mathit{x}}_i^T{\mathit{x}}_i.$$

Given the penalty parameters ${\lambda }_1>0$ and ${\lambda }_2>0$, the joint graphical lasso (JGL) is the procedure to find the positive definite matrix ${\Theta }^{\left(k\right)}\in {\mathbb{R}}^{p\times p}$ for $k=1,\dots ,K$, such that:

$$\underset{\Theta }{\mathrm{minimize}}\left\{-\sum _{k=1}^K{n}_k\{logdet{\Theta }^{\left(k\right)}-\mathrm{trace}\left({\mathit{S}}^{\left(k\right)}{\Theta }^{\left(k\right)}\right)\}+{\lambda }_1\sum _{k=1}^K\sum _{i\ne j}\left|{\theta }_{k,i,j}\right|+P(\Theta )\right\},$$

(15)

where $P(\Theta )$ penalizes $\Theta :={[{\Theta }^{\left(1\right)},\dots ,{\Theta }^{\left(K\right)}]}^T$. For example, ref. [13] suggested the following fused and group lasso penalties:

$$P_F(\Theta ):={\lambda }_2\sum _{k<l}\sum _{i,j}|{\theta }_{k,i,j}-{\theta }_{l,i,j}|$$

and

$$P_G(\Theta ):={\lambda }_2\sum _{i\ne j}{\left\{\sum _{k=1}^K{\theta }_{k,i,j}^2\right\}}^{1/2},$$

where ${\theta }_{k,i,j}$ is the $(i,j)$-th element of ${\Theta }^{\left(k\right)}\in {\mathbb{R}}^{p\times p}$ for $k=1,\cdots ,K$.

Unfortunately, there is no equation like (3) for the JGL to find the optimum ${\Theta }_{*}$. [13] considered the ADMM to solve the JGL problem. However, ADMM is quite time consuming for large-scale problems.

3. The Proposed Methods

In this section, we propose two efficient algorithms for solving the JGL problem. One is an extended ISTA based on the G-ISTA in Section 2.2, and the other is based on the step-size selection strategy introduced in Section 2.3.

3.1. ISTA for the JGL Problem

To describe the JGL problem, we define $f,g:{\mathbb{R}}^{K\times p\times p}\to \mathbb{R}$ by

$$\begin{array}{cc}\hfill f(\Theta )& :=-\sum _{k=1}^K{n}_k\left\{logdet{\Theta }^{\left(k\right)}-\mathrm{trace}\left({\mathit{S}}^{\left(k\right)}{\Theta }^{\left(k\right)}\right)\right\},\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill g(\Theta )& :={\lambda }_1\sum _{k=1}^K\sum _{i\ne j}\left|{\theta }_{k,i,j}\right|+P(\Theta ).\hfill \end{array}$$

(17)

Then, the problem (15) reduces to:

$$\underset{\Theta }{\mathrm{minimize}}F(\Theta ):=f(\Theta )+g(\Theta ),$$

where the function f is convex and differentiable, and g is convex and nondifferentiable. Therefore, the ISTA is available for solving the JGL problem (15).

The main difference between the G-ISTA and the proposed method is that the latter needs to simultaneously consider K categories of graphical models in the JGL problem (15). What is more, there are two combined penalties in $g(\Theta )$, which complicate the proximal operator in the ISTA procedure. Consequently, the operator for the proposed method is not a simple soft thresholding operator, as it is for the G-ISTA method.

If we define the quadratic approximation $Q_{{\eta }_t}:{\mathbb{R}}^{K\times p\times p}\to \mathbb{R}$ of $f\left({\Theta }_t\right)$ by:

$$\begin{array}{c}\hfill Q_{{\eta }_t}(\Theta ,{\Theta }_t):=f\left({\Theta }_t\right)+\sum _{k=1}^K〈{\Theta }^{\left(k\right)}-{\Theta }_t^{\left(k\right)},\nabla f\left({\Theta }_t^{\left(k\right)}\right)〉+\frac{1}{2{\eta }_t}\sum _{k=1}^K\left|\right|{\Theta }^{\left(k\right)}-{\Theta }_t^{\left(k\right)}{\left|\right|}_F^2,\end{array}$$

then the update iteration is simplified as:

$$\begin{array}{cc}\hfill {\Theta }_{t+1}& =\underset{\Theta }{\mathrm{argmin}}\left\{Q_{{\eta }_t}(\Theta ,{\Theta }_t)+g(\Theta )\right\}\hfill \\ \hfill & ={\mathrm{prox}}_{{\eta }_{t}g}({\Theta }_t-{\eta }_t\nabla f\left({\Theta }_t\right)).\hfill \end{array}$$

Nevertheless, the Lipschitz gradient constant of $f(\Theta )$ is unknown over the whole domain in the JGL problem. Therefore, our approach needs a backtracking line search to calculate step size ${\eta }_t$. We show the details in Algorithm 2.

Algorithm 2 ISTA for problem (15).

Input: $\mathit{S}$, tolerance $ϵ>0$, backtracking constant $0<c<1$, initial step size ${\eta }_0$, initial iterate ${\Theta }_0$. For$t=0,1,\cdots ,$ (until convergence) do

1: Backtracking line search: Continue to multiply ${\eta }_t$ by c until $$\begin{array}{c}\hfill f\left({\Theta }_{t+1}\right)\le Q_{{\eta }_t}({\Theta }_{t+1},{\Theta }_t)\mathrm{and}{\Theta }_{t+1}^{\left(k\right)}\succ 0\mathrm{for}k=1,\cdots ,K.\end{array}$$ for ${\Theta }_{t+1}:={\mathrm{prox}}_{{\eta }_{t}g}({\Theta }_t-{\eta }_t\nabla f\left({\Theta }_t\right))$.     2: Update iterate: ${\Theta }_{t+1}\leftarrow {\mathrm{prox}}_{{\eta }_{t}g}({\Theta }_t-{\eta }_t\nabla f\left({\Theta }_t\right))$.     3: Set next initial step size ${\eta }_{t+1}$. See details in Section 3.3.

end Output: optimal solution to problem (15), ${\Theta }_{*}={\Theta }_{t+1}$.

In the update of ${\Theta }_{t+1}$, we need to compute the proximal operators for the fused and group lasso penalties. In the following, for each of them, the problem can be divided into the fused lasso problems [29] and group lasso problems [30,31] for ${\mathit{\theta }}_{i,j}\in {\mathbb{R}}^K$, $i,j=1,\dots ,p$. We apply the solutions given by (20) and (21) below.

3.1.1. Fused Lasso Penalty $P_F$

By the definition of the proximal operator in the update step, we have:

$$\begin{array}{cc}\hfill {\Theta }_{t+1}=arg\underset{\Theta }{min}\left\{\frac{1}{2}\sum _{k=1}^K\left|\right|{\Theta }^{\left(k\right)}-{\Theta }_t^{\left(k\right)}+{\eta }_t\nabla f\left({\Theta }_t^{\left(k\right)}\right){\left|\right|}_F^2\right.& +{\eta }_t{\lambda }_1\sum _{k=1}^K\sum _{i\ne j}\left|{\theta }_{k,i,j}\right|\hfill \\ \hfill & \left.+{\eta }_t{\lambda }_2\sum _{k<l}\sum _{i,j}|{\theta }_{k,i,j}-{\theta }_{l,i,j}|\right\}.\hfill \end{array}$$

(18)

Problem (18) is separable with respect to the elements ${\theta }_{k,i,j}$ in ${\Theta }^{\left(k\right)}\in {\mathbb{R}}^{p\times p}$; hence, the proximal operator can be computed in a componentwise manner: Let $\mathit{A}={\Theta }_t-{\eta }_t\nabla f\left({\Theta }_t\right)$; then, problem (18) reduces to the following for $i=1,\cdots ,p$, $j=1,\cdots ,p$:

$$\begin{array}{c}\hfill \underset{{\theta }_{1,i,j},\cdots ,{\theta }_{K,i,j}}{\mathrm{argmin}}\left\{\frac{1}{2}\sum _{k=1}^K{({\theta }_{k,i,j}-a_{k,i,j})}^2+{\eta }_t{\lambda }_1{1}_{i\ne j}\sum _{k=1}^K|{\theta }_{k,i,j}|+{\eta }_t{\lambda }_2\sum _{k<l}|{\theta }_{k,i,j}-{\theta }_{l,i,j}\left|\right\}\right\},\end{array}$$

(19)

where $1_{i\ne j}$ is an indicator function, the value of which is 1 only when $i\ne j$.

The problem (19) is known as the fused lasso problem [29,32] given $a_{k,i,j}$ for $k=1,\cdots ,K$. In particular, let $\alpha :={\eta }_t{\lambda }_1{1}_{i\ne j}$ and $\beta :={\eta }_t{\lambda }_2$. When $i\ne j$, $\alpha \ne 0$ and $\beta >0$, the solution to (19) can be obtained through the soft thresholding operator based on the solution when $\alpha =0$ by the following Lemma.

Lemma 1.

([33]) Denote the solution to parameters α and β as ${\theta }_i(\alpha ,\beta )$, and then the solution ${\theta }_i(\alpha ,\beta )$ of the fused lasso problem:

$$\frac{1}{2}\sum _{i=1}^n{(y_i-{\theta }_i)}^2+\alpha \sum _{i=1}^n|{\theta }_i|+\beta \sum _{i=1}^{n-1}|{\theta }_i-{\theta }_{i+1}|$$

(20)

is given by ${\left[S_{\alpha }\left(\theta (0,\beta )\right)\right]}_i$ when $y_1,\cdots ,y_n\in \mathbb{R}$ are given for $n\ge 1$.

Additionally, rather efficient algorithms are available for solving the fused lasso problem (20) when $\alpha =0$ (i.e., $\theta (0,\beta )$ ) such as [32,34,35].

3.1.2. Group Lasso Penalty $P_G$

By definition, the update of ${\Theta }_{t+1}$ for the group lasso penalty $P_G(\Theta )$ is as follows:

$$\begin{array}{cc}\hfill {\Theta }_{t+1}=arg\underset{\Theta }{min}\left\{\frac{1}{2}\sum _{k=1}^K\left|\right|{\Theta }^{\left(k\right)}-{\Theta }_t^{\left(k\right)}+{\eta }_t\nabla f\left({\Theta }_t^{\left(k\right)}\right){\left|\right|}_F^2\right.& +{\eta }_t{\lambda }_1\sum _{k=1}^K\sum _{i\ne j}\left|{\theta }_{k,i,j}\right|\hfill \\ \hfill & \left.+{\eta }_t{\lambda }_2\sum _{i\ne j}{(\sum _{k=1}^K{\theta }_{k,i,j}^2)}^{1/2}\right\}.\hfill \end{array}$$

Similarly, let $\mathit{A}={\Theta }_t-{\eta }_t\nabla f\left({\Theta }_t\right)$; then, the problem becomes the following for $i=1,\cdots ,p$, $j=1,\cdots ,p$:

$$\begin{array}{c}\hfill \underset{{\theta }_{1,i,j},\cdots ,{\theta }_{K,i,j}}{\mathrm{argmin}}\left\{\frac{1}{2}\sum _{k=1}^K{({\theta }_{k,i,j}-a_{k,i,j})}^2+{\eta }_t{\lambda }_1{1}_{i\ne j}\sum _{k=1}^K\left|{\theta }_{k,i,j}\right|+{\eta }_t{\lambda }_2{1}_{i\ne j}{(\sum _{k=1}^K{\theta }_{k,i,j}^2)}^{1/2}\right\}.\end{array}$$

We have ${\theta }_{k,i,j}=a_{k,i,j}$ for $i=j$. In addition, for $i\ne j$, the solution [31,36,37] is given by

$${\theta }_{k,i,j}={\mathcal{S}}_{{\eta }_t{\lambda }_1}\left(a_{k,i,j}\right){\left(1-\frac{{\eta }_t{\lambda }_2}{\sqrt{{\sum }_{k=1}^K{\mathcal{S}}_{{\eta }_t{\lambda }_1}{\left(a_{k,i,j}\right)}^2}}\right)}_{+}.$$

(21)

3.2. Modified ISTA for JGL

Thus far, we have seen that $f(\Theta )$ in the JGL problem (15) is not globally Lipschitz gradient continuous. The ISTA may not be efficient enough for the JGL case because it includes the backtracking line search procedure for this case, which needs to evaluate the objective function and the positive definiteness of ${\Theta }_{t+1}$ in Step 1 of Algorithm 2 and is inefficient when the evaluation is expensive.

In this section, we modify Algorithm 2 to Algorithm 3 based on the step-size selection strategy in Section 2.3, which takes advantage of the properties of the self-concordant function. The self-concordant function does not rely on the Lipschitz gradient assumption on the function $f(\Theta )$ [24], and we can eliminate the need for the backtracking line search.

Lemma 2.

([38]) Self-concordance is preserved by scaling and addition: if f is a self-concordant function and a constant $a\le 1$, then $af$ is self-concordant. If $f_1,f_2$ are self-concordant, then $f_1+f_2$ is self-concordant.

By Lemma 2, the function $f(\Theta )$ in (16) is self-concordant. In Algorithm 3, for the initial step size of ${\eta }_t$ in each iteration, we use the Barzilai–Borwein method [39]. We apply the step-size mechanism in Section 2.3, which is employed in Steps 3–5 of Algorithm 3.

Algorithm 3 Modified ISTA (M-ISTA).

Input: $\mathit{S}$, tolerance $ϵ>0$, initial step size ${\eta }_0$, initial iterate ${\Theta }_0$. For$t=0,1,\cdots ,$ (until convergence) do

1: Initialize ${\eta }_t$.     2: Compute $$\begin{array}{c}\hfill {\mathit{d}}_t:={\mathrm{prox}}_{{\eta }_{t}g}({\Theta }_t-{\eta }_t\nabla f\left({\Theta }_t\right))-{\Theta }_t.\end{array}$$     3: Compute $$\begin{array}{c}\hfill {\beta }_t:={\eta }_t^{-1}\left|\right|{\mathit{d}}_t{\left|\right|}_F^2\end{array}$$ and $$\begin{array}{c}\hfill {\lambda }_t:=\sum _{k=1}^K\sqrt{n_k}\left|\right|{\left({\Theta }_t^{\left(k\right)}\right)}^{-1}{\mathit{d}}_t^{\left(k\right)}{\left|\right|}_F.\end{array}$$     4: Determine the step size ${\alpha }_t:=\frac{{\beta }_t}{{\lambda }_t({\lambda }_t+{\beta }_t)}$.     5: If ${\alpha }_t>1$, then set ${\eta }_t:={\eta }_t/2$ and go back to Step 2.     6: Update ${\Theta }_{t+1}:={\Theta }_t+{\alpha }_t{\mathit{d}}_t$.

end Output: optimal solution to problem (15), ${\Theta }_{*}={\Theta }_{t+1}$.

There is no backtracking procedure in this algorithm that guarantees the positive definiteness of ${\Theta }_{t+1}$, as in Step 1 of Algorithm 2. We next show how to ensure the positive definiteness of ${\Theta }_{t+1}$ in the iterations of Algorithm 3.

Lemma 3.

([40], Theorem 2.1.1) Let f be a self-concordant function, and let $\mathit{x}\in \mathrm{dom}f$. Additionally, if

$$W\left(\mathit{x}\right)=\{\mathit{y}|{\left(\langle {\nabla }^{2}f\left(\mathit{x}\right)(\mathit{y}-\mathit{x}),\mathit{y}-\mathit{x}\rangle \right)}^{1/2}\le 1\},$$

then $W\left(\mathit{x}\right)\subset domf$.

In Algorithm 3, because we know ${\alpha }_t:=\frac{{\beta }_t}{{\lambda }_t({\lambda }_t+{\beta }_t)}<1$ with ${\beta }_t>0$ and ${\lambda }_t>0$ by Steps 3–5. Thus, we have ${\alpha }_t{\lambda }_t<1$:

$${\alpha }_t{\lambda }_t:={\alpha }_t{\left(\langle {\nabla }^{2}f\left({\Theta }_t\right){\mathit{d}}_t,{\mathit{d}}_t\rangle \right)}^{1/2}<1,$$

which implies,

$${\left(\langle {\nabla }^{2}f\left({\Theta }_t\right)({\Theta }_{t+1}-{\Theta }_t),{\Theta }_{t+1}-{\Theta }_t\rangle \right)}^{1/2}<1.$$

Hence, from Lemma 3, we see that ${\Theta }_{t+1}$ stays in the domain and maintains positive definiteness.

3.3. Theoretical Analysis

For multiple Gaussian graphical models, Honorio and Samaras [14] and Hara and Washio [17] provided lower and upper bounds for the optimal solution ${\Theta }_{*}$. However, the models they considered are different than the JGL. To the best of our knowledge, no related research has provided the bounds of the optimal solution ${\Theta }_{*}$ for the JGL problem (15).

In the following section, we show the bounds of the optimal solution ${\Theta }_{*}$ for the JGL and the iterates ${\Theta }_t$ generated by Algorithms 2 and 3, which are applied to both fused and group lasso-type penalties.

Proposition 1.

The optimal solution ${\Theta }_{*}$ of the problem (15) satisfies

$$\underset{1\le k\le K}{max}\frac{n_k}{p{\lambda }_c+n_k\left|\right|{\mathit{S}}^{\left(k\right)}{\left|\right|}_2}\le \left|\right|{\Theta }_{*}^{\left(k\right)}{\left|\right|}_2\le \frac{Np}{{\lambda }_1}+\sum _{k=1}^K\sum _{i=1}^p{\left(s_{k,i,i}\right)}^{-1},$$

where ${\lambda }_c:=\sqrt{K{\lambda }_1^2+2K{\lambda }_1{\lambda }_2+{\lambda }_2^2}$, and $s_{k,i,i}$ is the i-th diagonal element of ${\mathit{S}}^{\left(k\right)}$.

For the proof, see Appendix A.1.

Note that the objective function value $F(\Theta )$ always decreases with the increase in iteration in both algorithms due to [25] (Remark 3.1) and Lemma 12 in [24]. Therefore, the following inequality holds for Algorithms 2 and 3:

$$F\left({\Theta }_{t+1}\right)\le F\left({\Theta }_t\right)\mathrm{for}t=0,1,\cdots .$$

(22)

Then, based on the condition (22), we provide the explicit bounds of iterates ${\left\{{\Theta }_t\right\}}_{t=0,1\cdots }$ in Algorithms 2 and 3 for the JGL problem (15).

Proposition 2.

Sequence ${\left\{{\Theta }_t\right\}}_{t=0,1,\cdots ,}$ generated by Algorithms 2 and 3 can be bounded:

$$m\le \left|\right|{\Theta }_t{\left|\right|}_2\le M,$$

where $M:=\left|\right|{\Theta }_0{\left|\right|}_F+\frac{2Np}{{\lambda }_1}+2{\displaystyle \sum _{k=1}^K}{\displaystyle \sum _{i=1}^p}{s_{k,i,i}}^{-1}$, $m:=e^{-\frac{C_1}{n_m}}{M}^{(1-Kp)}$, $n_m=\underset{k}{max}{n}_k$, and constant $C_1:=F\left({\Theta }_0\right)$.

For the proof, see Appendix A.2.

With the help of Propositions 1 and 2, and the following Lemma, we can obtain the range of the step size that ensures the linear convergence rate of Algorithm 2.

Lemma 4.

Let ${\Theta }_t$ be t-th iterate in Algorithm 2. Denote ${\lambda }_{min}$ and ${\lambda }_{max}$ as the minimum and maximum eigenvalues of the corresponding matrix, respectively. Define

$$a_k:=min\{{\lambda }_{min}\left({\Theta }_t^{\left(k\right)}\right),{\lambda }_{min}\left({\Theta }_{*}^{\left(k\right)}\right)\},b_k:=max\{{\lambda }_{max}\left({\Theta }_t^{\left(k\right)}\right),{\lambda }_{max}\left({\Theta }_{*}^{\left(k\right)}\right)\}$$

and $n_l={min}_{k=1\cdots ,K}{n}_k$, $n_m={max}_{k=1\cdots ,K}{n}_k$, $a_l={min}_{k=1\cdots ,K}{a}^{\left(k\right)}$, and $b_m={max}_{k=1\cdots ,K}{b}^{\left(k\right)}$. The sequence ${\left\{{\Theta }_t\right\}}_{t=0,1,\cdots }$ generated by Algorithm 2 satisfies

$$\left|\right|{\Theta }_{t+1}-{\Theta }_{*}{\left|\right|}_F\le {\gamma }_t\left|\right|{\Theta }_t-{\Theta }_{*}{\left|\right|}_F$$

with the convergence rate ${\gamma }_t:=max\{\frac{{\eta }_t{n}_m}{a_l^2}-1,1-\frac{{\eta }_t{n}_l}{b_m^2}\}$.

Proof. 

It can be easily extended by Lemma 3 in [10]. □

Lemma 4 implies that to obtain the convergence rate ${\gamma }_t<1$, we require:

$$0<{\eta }_t<\frac{2a_l^2}{n_m}.$$

(23)

After using Propositions 1 and 2, we can obtain the bounds of $a_l$. Further, we can obtain the step size ${\eta }_t$ that satisfies (23) and guarantee s the linear convergence rate $({\gamma }_t<1)$. However, the step size is quite conservative in practice. Hence, we consider the Barzilai–Borwein method for implementation and regard the step size ${\eta }_t$ that satisfies (23) as a safe choice. When the number of backtracking iterations in Step 1 of Algorithm 2 exceeds the given maximum number to fulfill the backtracking line search condition, we can use the safe step size ${\eta }_t$ for the subsequent calculations. In Section 4.2.3, we confirm the linear convergence rate of the proposed ISTA by experiment.

4. Experiments

In this section, we evaluate the performance of the proposed methods on both synthetic and real datasets, and we compare the following algorithms:

• ADMM: the general ADMM method proposed by [13].
• FMGL: the proximal Newton-type method proposed by [23].
• ISTA: the proposed method in Algorithm 2.
• M-ISTA: the proposed method in Algorithm 3.

We perform all the tests in R Studio on a Macbook Air with 1.6 GHz Intel Core i5 and 8 GB memory. The wall times are recorded as the run times for the four algorithms.

4.1. Stopping Criteria and Model Selection

In the experiments, we consider two stopping criteria for the algorithms.

1. Relative error stopping criterion:

$$\frac{{\sum }_{k=1}^K\left|\right|{\Theta }_{t+1}^{\left(k\right)}-{\Theta }_t^{\left(k\right)}{\left|\right|}_F}{max\left\{{\sum }_{k=1}^K\right|\left|{\Theta }_t^{\left(k\right)}\right|{|}_F,1\}}\le ϵ.$$

2. Objective error stopping criterion:

$$F\left({\Theta }_t\right)-F\left({\Theta }_{*}\right)\le ϵ.$$

$ϵ$ is a given accuracy tolerance; we terminate the algorithm if the above error is smaller than $ϵ$ or the maximum number of iterations exceeds 1000. We use the objective error for convergence rate analysis and the relative error for the time comparison.

The JGL model is affected by regularized parameters ${\lambda }_1$ and ${\lambda }_2$. For selecting the parameters, we use the V-fold crossvalidation method. First, the dataset is randomly split into V segments of equal size, a single subset (test data), estimated by the other $V-1$ subsets (training data), is evaluated, and the subset is changed for the test to repeat V times so that each subset is used.

Let ${\mathit{S}}_v^{\left(k\right)}$ be the sample covariance matrix of the v-th ( $v=1,\dots ,V$) segment for class $k=1,\cdots ,K$. We estimate the inverse covariance matrix by the remaining $V-1$ subsets ${\widehat{\Theta }}_{\lambda ,-v}^{\left(k\right)}$ and choose ${\lambda }_1$ and ${\lambda }_2$, which minimize the average predictive negative log-likelihood as follows:

$$CV({\lambda }_1,{\lambda }_2)=\sum _{v=1}^V\sum _{k=1}^K\left\{n_k\mathrm{trace}\left({\mathit{S}}_v^{\left(k\right)}{\widehat{\Theta }}_{\lambda ,-v}^{\left(k\right)}\right)-\mathrm{logdet}{\widehat{\Theta }}_{\lambda ,-v}^{\left(k\right)}\right\}$$

4.2. Synthetic Data

The performance of the proposed methods was assessed on synthetic data in terms of the number of iterations, the execution time, the squared error, and the receiver operating characteristic (ROC) curve. We follow the data generation mechanism described in [41] with some modifications for the JGL model. We put the details in Appendix B.

4.2.1. Time Comparison Experiments

We vary $p,N,K\mathrm{and}{\lambda }_1$ to compare the execution time of our proposed methods with that of the existing methods. We consider only the fused penalty in our proposed method for a fair comparison in the experiments because the FMGL algorithm applies only to the fused penalty. First, we compare the performance among different algorithms under various dimensions p, which are shown in Figure 1.

Figure 1 shows that the execution time of the FMGL and ADMM increases rapidly as p increases. In particular, we observe that the M-ISTA significantly outperforms when p exceeds 200. The ISTA shows better performance than the three methods when p is less than 200, but it requires more time as p grows, compared to the M-ISTA. It is reasonable to consider that evaluating the objective function in the backtracking line search at every iteration increases the computational burden, especially when p increases, which means that the M-ISTA is a good choice for these cases. Furthermore, the ISTA can be a good candidate when the evaluation is inexpensive.

Table 2. Computational time under different settings.

Parameters Setting						Computational Time
p	K	N	${\mathit{\lambda }}_1$	${\mathit{\lambda }}_2$	precision$\mathit{ϵ}$	ADMM	FMGL	ISTA	M-ISTA
20	2	60	0.1	0.05	0.00001	10.506 s	1.158 s	2.174 s	1.742 s
	3					1.879 min	4.267 s	3.357 s	3.668 s
	5		1	0.5		1.123 min	10.556 s	4.216 s	2.874 s
30	2	120	0.1	0.05	0.0001	10.095 s	5.259 s	2.690 s	4.857 s
	3					2.014 min	38.562 s	14.722 s	31.870 s
	5		1	0.5		2.447 min	15.819 s	22.431 s	12.113 s
50	2	600	0.02	0.005	0.0001	6.427 s	10.228 s	7.213 s	4.625 s
			0.03			6.240 s	8.925 s	6.645 s	4.023 s
			0.04			7.025 s	9.381 s	6.144 s	3.993 s
200	2	400	0.09	0.05	0.0001	4.050 min	1.874 min	2.289 min	35.038 s
			0.1			4.569 min	1.137 min	1.340 min	24.852 s
			0.12			3.848 min	1.881 min	1.443 min	18.367 s

summarizes the performance of the four algorithms under different parameter settings to achieve a given precision, $ϵ$, of the relative error. The results presented in Table 2 reveal that when we increase the number of classes K, all the algorithms require more time than usual. Moreover, the execution time of ADMM becomes huge among them. When we vary ${\lambda }_1$, the algorithms become more efficient as the value increases. For most instances, the M-ISTA and ISTA outperform the existing methods, such as ADMM and FMGL. For the exceptional cases ($p=20,k=2,N=60,{\lambda }_1=0.1$ and ${\lambda }_2=0.05$), the M-ISTA and ISTA are still comparable with the FMGL and faster than ADMM.

4.2.2. Algorithm Assessment

We generate the simulation data as described in Appendix B and regard the synthetic inverse covariance matrices ${\Theta }^{\left(k\right)}$ as the true values for our assessment experiments.

First, we assessed our proposed method by drawing an ROC curve, which displays the number of true positive edges (i.e., TP edges) selected compared to the number of false positive edges (i.e., FP edges) selected. We say that an edge $(i,j)$ in the k-th class is selected in estimate ${\widehat{\Theta }}^{\left(k\right)}$ if element ${\widehat{\theta }}_{k,i,j}\ne 0$, and the edges are true positive edges selected if the precision matrix element ${\theta }_{k,i,j}\ne 0$ and false positive edges selected if the precision matrix element ${\theta }_{k,i,j}=0$, where the two quantities are defined by

$$TP=\sum _{k=1}^K\sum _{i,j}1({\theta }_{k,i,j}\ne 0)·1({\widehat{\theta }}_{k,i,j}\ne 0)$$

and

$$FP=\sum _{k=1}^K\sum _{i,j}1({\theta }_{k,i,j}=0)·1({\widehat{\theta }}_{k,i,j}\ne 0),$$

where $1(·)$ is the indicator function.

To confirm the validity of the proposed methods, we compare the ROC figures of the fused penalty and group penalty. We fix the parameters ${\lambda }_2$ for each curve and change the ${\lambda }_1$ value to obtain various numbers of selected edges because the sparsity penalty parameter ${\lambda }_1$ can control the number of selected total edges.

We show the ROC curves for fused and group lasso penalties in Figure 2a,b respectively. From the figures, we observe that both penalties show highly accurate predictions for the edge selections. The result of ${\lambda }_2=0.0166$ in the fused penalty case is better than that in ${\lambda }_2=0.05$. Additionally, the result of ${\lambda }_2=0.0966$ in the group penalty case is better than that in ${\lambda }_2=0.09$, which means that if we select the tuning parameters properly, then we can obtain precise results while simultaneously meeting our different model demands.

Then, Figure 3a,b display the mean squared error (MSE) between the estimated values and true values.

$$MSE=\frac{2}{Kp(p-1)}\sum _{k=1}^K\sum _{i<j}{({\widehat{\theta }}_{k,i,j}-{\theta }_{k,i,j})}^2,$$

where ${\widehat{\theta }}_{k,i,j}$ is the value estimated by the proposed method, and ${\theta }_{k,i,j}$ is the true precision matrix value we used in the data generation.

The figures illustrate that when the total number of edges selected increases, the errors decrease and finally achieve relatively low values.

Overall, the proposed method shows competitive efficiency not only in computational time but also in accuracy.

4.2.3. Convergence Rate

This section shows the convergence rate of the ISTA for solving the JGL problem (15) in practice, with ${\lambda }_1=0.1$, $0.09$ and $0.08$. We recorded the number of iterations to achieve the different tolerance of $F\left({\Theta }_t\right)-F\left({\Theta }_{*}\right)$ in Figure 4 and ran it on a synthetic dataset, with $p=200$, $K=2$, ${\lambda }_2=0.05$, and $N=400$. The figure reveals that as ${\lambda }_1$ decreases, more iterations are needed to converge to the specified tolerance. Moreover, the figure shows the linear convergence rate of the proposed ISTA method, which corroborate the theoretical analysis in Section 3.3.

4.3. Real Data

In this section, we use two different real datasets to demonstrate the performance of our proposed method and visualize the result.

Firstly, we used the presidential speeches dataset in [42] for the experiment to jointly estimate common links across graphs and show the common structure. The dataset contains 75 most-used words (features) from several big speeches of the 44 US presidents (samples). In addition, we used the clustering result in [42], where the authors split the 44 samples into two groups with similar features, and then we obtained two classes of samples $(K=2)$.

We used Cytoscape [43] to visualize the results when ${\lambda }_1=1.9$ and ${\lambda }_2=0.16$. We chose these relatively large tuning parameters for better interpretation of the network figure. Figure 5 shows the relationship network graph of the high-frequency words identified by the JGL model with the proposed method. As shown in the figure, each node represents a word, and the edges demonstrate the relationships between words.

We use different colors to show various structures. The black edges are a common structure between the two classes, the red edges are the specific structures for the first class $(k=1)$, and the green edges are for the second class $(k=2)$. Figure 5 shows a subnetwork on the top with red edges, meaning there are relationships among those words, and the connections only exist in the first group.

We compared the time cost among four algorithms and show the results in

Table 3. Time comparison result of two real datasets.

Dataset	Parameters Setting			Computational Time
	${\mathit{\lambda }}_1$	${\mathit{\lambda }}_2$	Precision $\mathit{ϵ}$	ADMM	FMGL	ISTA	M-ISTA
Speeches	0.1	0.05	0.0001	19.969 s	4.977 min	11.829 s	12.867 s
	0.2	0.1		4.661 min	3.209 min	11.560 s	12.682 s
	0.5	0.25		5.669 min	1.490 min	11.043 s	12.788 s
Breast cancer	0.1	0.0166	0.0001	3.809 min	7.937 min	1.305 min	1.158 min
	0.2	0.02		6.031 min	5.198 min	1.503 min	1.230 min
	0.3	0.03		5.499 min	2.265 min	1.188 min	1.061 min

. We used the crossvalidation method ($V=6$) described in Section 4.1 to select the optimal tuning parameters (${\lambda }_1=0.1$, ${\lambda }_2=0.05$). In addition, we manually chose the other two pairs of parameters for more comparisons.

Table 3 shows that ISTA outperforms the other three algorithms, and our proposed methods offer stable performance when varying the parameters, while ADMM is the slowest in most cases.

Secondly, we use a breast cancer dataset [44] for time comparison. There are 250 samples and 1000 genes in the dataset, with 192 control samples and 58 case samples $(K=2)$. Furthermore, we extract 200 genes with the highest variances among the original genes. The tuning parameter pair $({\lambda }_1=0.01$, ${\lambda }_2=0.0166)$ was chosen by the crossvalidation method. Table 3 exhibits that our proposed methods (ISTA and M-ISTA) outperform ADMM and FMGL, and M-ISTA shows the best performance in the breast cancer dataset.

5. Discussion

We propose two efficient proximal gradient descent procedures with and without the backtracking line search option for the joint graphical lasso. The first (Algorithm 2) does not require extra variables, unlike ADMM, which needs manual tuning the Lagrangian penalty parameters $\rho$ in [13] and storing and calculating dual variables. Moreover, we reduce the update iterate step to subproblems that can be solved efficiently and precisely by lasso-type problems. Based on Algorithm 2, we modified the step-size selection by extending the strategy in [24] to the second one (Algorithm 3), which does not rely on the Lipschitz assumption. Additionally, the second does not require a backtracking line search, significantly reducing the computation time needed to evaluate objective functions.

From the theoretical perspective, we reach the linear convergence rate for the ISTA. Furthermore, we derive the lower and upper bounds of the solution to the JGL problem and the iterates in the algorithms, guaranteeing that the ISTA converges linearly. Numerically, the methods are demonstrated on simulated and real datasets to illustrate their robust and efficient performance over state-of-the-art algorithms.

For further computational improvement, the most expensive step in the algorithms is to calculate the inversion of matrices required by the gradient of $f(\Theta )$. Both algorithms have a complexity of $O\left(K{p}^3\right)$ per iteration. Moreover, we can solve the matrix inversion problem with more efficient algorithms with lower complexity. In addition, we can also use the faster computation procedure in [13] to decompose the optimization problem for the proposed methods and regard it as preprocessing. Overall, the proposed methods are highly efficient for the joint graphical lasso problem.