An Efficient Linearly Convergent Regularized Proximal Point Algorithm for Fused Multiple Graphical Lasso Problems

TL;DR

This paper introduces a new regularized proximal point algorithm with a superlinearly convergent semismooth Newton method for efficiently solving fused multiple graphical Lasso problems, improving convergence speed and robustness.

Contribution

It develops an efficient second-order method leveraging explicit Jacobian expressions, outperforming existing first-order approaches in fused graphical model estimation.

Findings

The proposed algorithm converges faster than state-of-the-art methods.

It demonstrates robustness and efficiency on synthetic and real datasets.

Explicit Jacobian derivation enhances the semismooth Newton method's performance.

Abstract

Nowadays, analysing data from different classes or over a temporal grid has attracted a great deal of interest. As a result, various multiple graphical models for learning a collection of graphical models simultaneously have been derived by introducing sparsity in graphs and similarity across multiple graphs. This paper focuses on the fused multiple graphical Lasso model which encourages not only shared pattern of sparsity, but also shared values of edges across different graphs. For solving this model, we develop an efficient regularized proximal point algorithm, where the subproblem in each iteration of the algorithm is solved by a superlinearly convergent semismooth Newton method. To implement the semismooth Newton method, we derive an explicit expression for the generalized Jacobian of the proximal mapping of the fused multiple graphical Lasso regularizer. Unlike those widely used first order methods, our approach has heavily exploited the underlying second order information through the semismooth Newton method. This can not only accelerate the convergence of the algorithm, but also improve its robustness. The efficiency and robustness of our proposed algorithm are demonstrated by comparing with some state-of-the-art methods on both synthetic and real data sets. Supplementary materials for this article are available online.

Tables

Table 1. Table 1 : Performances of rPPA, ADMM, and MGL on stock price data. Tolerance ε = 𝜀 absent \varepsilon= 1e-6.

Problem	$({\lambda }_1,{\lambda }_2)$	Density	Iteration			Time			Error
$(p,L)$			P	A	M	P	A	M	P	A	M
	(1e-04,1e-05)	0.039	25	3701	6	06	33	33	6.0e-07	9.5e-07	1.1e-06
SPX3a	(5e-05,5e-06)	0.144	24	3701	9	08	33	43	9.9e-07	9.5e-07	4.0e-06
(100,3)	(2e-05,2e-06)	0.241	26	5359	19	11	49	09	7.0e-07	1.0e-06	1.2e-05
	(1e-04,1e-05)	0.025	24	3301	9	15	01:14	01:26	8.1e-07	8.2e-07	2.3e-06
SPX3a	(5e-05,5e-06)	0.086	24	3301	17	22	01:15	03:22	6.8e-07	8.3e-07	4.7e-06
(200,3)	(2e-05,2e-06)	0.150	26	5920	44	30	02:17	03:59	6.4e-07	1.0e-06	8.0e-06
	(5e-04,5e-05)	0.028	24	3701	8	18	01:23	02:40	9.6e-07	1.0e-06	3.2e-06
SPX11b	(1e-04,1e-05)	0.126	24	3701	111	27	01:55	04:59	9.3e-07	1.0e-06	4.9e-06
(100,11)	(5e-05,5e-06)	0.206	24	3710	388	30	01:59	13:16	9.4e-07	1.0e-06	5.6e-06
	(5e-04,5e-05)	0.017	22	3501	28	54	03:58	46:47	6.2e-07	9.9e-07	1.7e-06
SPX11b	(1e-04,1e-05)	0.081	24	3601	477	01:19	05:13	01:34:40	9.5e-07	8.8e-07	3.5e-06
(200,11)	(5e-05,5e-06)	0.134	24	3573	1076	01:38	05:21	03:00:00	9.8e-07	1.0e-06	4.0e-05

Table 2. Table 2 : Distribution of documents of classes Student, Faculty, Course, and Project.

Class	Student	Faculty	Course	Project	Total
#train docs	1097	750	620	336	2803
#test docs	544	374	310	168	1396

Table 3. Table 3 : Performances of rPPA, ADMM, and MGL on webpages data. Tolerance ε = 𝜀 absent \varepsilon= 1e-6.

Problem	$({\lambda }_1,{\lambda }_2)$	Density	Iteration			Time			Error
$(n,L)$			P	A	M	P	A	M	P	A	M
	(1e-02,1e-03)	0.015	16	2401	4	06	25	05	5.8e-07	9.9e-07	1.2e-07
Webtest	(5e-03,5e-04)	0.047	16	2401	6	07	27	10	5.6e-07	9.9e-07	3.1e-07
(100,4)	(1e-03,1e-04)	0.219	15	701	38	05	09	49	5.7e-07	6.1e-07	9.0e-07
	(1e-02,1e-03)	0.008	18	2101	7	11	01:03	59	7.3e-07	8.5e-07	8.6e-07
Webtest	(5e-03,5e-04)	0.025	18	2101	8	12	01:04	01:05	6.8e-07	8.5e-07	4.7e-07
(200,4)	(1e-03,1e-04)	0.156	18	2101	72	12	01:12	07:31	5.5e-07	7.1e-07	9.3e-07
	(5e-03,5e-04)	0.016	18	2101	9	25	02:12	03:44	5.9e-07	8.3e-07	3.7e-07
Webtest	(1e-03,1e-04)	0.119	17	2101	258	49	02:16	39:58	7.9e-07	8.6e-07	1.0e-06
(300,4)	(5e-04,5e-05)	0.244	19	2901	1393	01:18	03:07	02:22:23	5.6e-07	8.0e-07	1.0e-06
	(1e-02,1e-03)	0.011	24	20000	3	12	03:11	04	8.4e-07	7.0e-06	3.4e-06
Webtrain	(5e-03,5e-04)	0.030	24	20000	5	13	03:35	08	8.1e-07	7.0e-06	8.9e-07
(100,4)	(1e-03,1e-04)	0.162	24	20000	22	14	03:52	01:06	7.3e-07	7.0e-06	3.5e-06
	(5e-03,5e-04)	0.015	24	20000	5	46	10:53	17	8.5e-07	6.7e-06	1.2e-06
Webtrain	(1e-03,1e-04)	0.105	24	15227	33	44	08:23	03:02	7.2e-07	1.0e-06	2.6e-06
(200,4)	(5e-04,5e-05)	0.210	24	20000	95	52	10:41	06:34	6.9e-07	6.6e-06	2.7e-06
	(5e-03,5e-04)	0.010	24	20000	7	01:31	21:31	02:07	8.1e-07	5.9e-06	-1.2e-08
Webtrain	(1e-03,1e-04)	0.077	24	20000	52	01:33	22:06	20:47	6.3e-07	5.9e-06	1.8e-06
(300,4)	(5e-04,5e-05)	0.168	24	20000	155	01:51	21:56	18:09	6.1e-07	5.9e-06	1.8e-06

Table 4. Table 1 : Partition of 20 newsgroups by topics

comp.graphics	rec.autos	sci.crypt
comp.os.ms-windows.misc	rec.motorcycles	sci.electronics
comp.sys.ibm.pc.hardware	rec.sport.baseball	sci.med
comp.sys.mac.hardware	rec.sport.hockey	sci.space
comp.windows.x
	talk.politics.misc	talk.religion.misc
misc.forsale	talk.politics.guns	alt.atheism
	talk.politics.mideast	soc.religion.christian

Table 5. Table 2 : Performances of rPPA, ADMM, and MGL on newsgroups data. Tolerance ε = 𝜀 absent \varepsilon= 1e-6.

Problem	$({\lambda }_1,{\lambda }_2)$	Density	Iteration			Time			Error
$(n,L)$			P	A	M	P	A	M	P	A	M
NGcomp	(5e-03,5e-04)	0.020	19	4201	35	53	05:27	41:26	5.9e-07	9.2e-07	6.2e-07
test	(1e-03,1e-04)	0.094	16	1543	720	01:20	02:08	01:54:03	8.8e-07	1.0e-06	1.0e-06
(300,5)	(5e-04,5e-05)	0.194	17	1391	1240	59	01:56	03:00:00	7.9e-07	1.0e-06	9.6e-06
NGrec	(5e-03,5e-04)	0.004	25	20000	4	01:16	21:18	04:43	7.1e-07	8.4e-06	8.3e-07
test	(1e-03,1e-04)	0.061	25	20000	13	01:20	22:10	04:38	6.3e-07	8.4e-06	4.6e-07
(300,4)	(5e-04,5e-05)	0.134	24	20000	37	01:23	22:05	07:55	9.8e-07	8.4e-06	1.9e-06
NGsci	(5e-03,5e-04)	0.006	22	16244	6	57	15:42	06:00	6.0e-07	1.0e-06	2.7e-07
test	(1e-03,1e-04)	0.074	21	16230	25	01:06	17:54	10:22	7.9e-07	1.0e-06	1.3e-06
(300,4)	(5e-04,5e-05)	0.156	21	16230	100	01:06	17:52	20:37	7.3e-07	1.0e-06	1.5e-06
NGtalk	(5e-03,5e-04)	0.026	17	4179	14	52	04:23	14:33	7.4e-07	1.0e-06	2.7e-07
test	(1e-03,1e-04)	0.111	17	1018	79	35	01:08	16:52	4.9e-07	1.0e-06	9.8e-07
(300,3)	(5e-04,5e-05)	0.228	17	922	434	37	01:00	52:12	6.8e-07	1.0e-06	9.9e-07
NGcomp	(5e-03,5e-04)	0.016	20	6023	13	36	07:36	26:32	8.6e-07	1.0e-06	5.2e-07
train	(1e-03,1e-04)	0.077	20	6393	153	40	08:59	36:18	6.0e-07	1.0e-06	9.7e-07
(300,5)	(5e-04,5e-05)	0.142	19	5861	662	52	07:58	01:52:34	8.3e-07	1.0e-06	1.1e-06
NGrec	(5e-03,5e-04)	0.004	26	8842	5	01:30	07:17	03:34	6.4e-07	1.0e-06	4.7e-07
train	(1e-03,1e-04)	0.067	24	8737	17	01:58	09:31	06:32	7.6e-07	1.0e-06	1.7e-06
(300,4)	(5e-04,5e-05)	0.119	24	8625	66	02:05	09:30	14:16	7.1e-07	1.0e-06	2.1e-06
NGsci	(5e-03,5e-04)	0.011	21	11166	10	41	11:15	11:32	7.9e-07	1.0e-06	2.1e-08
train	(1e-03,1e-04)	0.085	20	11137	41	01:03	12:17	13:45	9.4e-07	1.0e-06	1.5e-06
(300,4)	(5e-04,5e-05)	0.146	20	11405	226	01:02	12:41	28:27	9.1e-07	1.0e-06	1.6e-06
NGtalk	(5e-03,5e-04)	0.026	22	20000	12	01:29	21:04	13:18	9.4e-07	2.0e-06	1.4e-06
train	(1e-03,1e-04)	0.101	22	20000	74	41	22:10	12:00	6.5e-07	1.7e-06	1.6e-06
(300,3)	(5e-04,5e-05)	0.193	22	20000	402	40	23:15	45:12	6.2e-07	1.6e-06	1.6e-06

Equations

$$\begin{array}[]{cl}\min\limits_{\Theta}&\displaystyle\sum^{L}_{l=1}\left(-\log\det\,\Theta^{(l)}+\langle S^{(l)},\Theta^{(l)}\rangle\right)+\mathcal{P}(\Theta),\end{array}$$

$$\begin{array}[]{l}\mathcal{P}(\Theta)\displaystyle=\lambda_{1}\sum^{L}_{l=1}\sum_{i\neq j}|\Theta^{(l)}_{ij}|+\lambda_{2}\sum^{L}_{l=2}\sum_{i\neq j}|{\Theta^{(l)}_{ij}}-{\Theta^{(l-1)}_{ij}}|.\end{array}$$

$$\begin{array}[]{l}\Psi_{\Xi}(u):=\min_{u^{\prime}}\left\{\Xi(u^{\prime})+\frac{1}{2}\|u^{\prime}-u\|^{2}\right\},\,\,\forall u\in\mathcal{E}.\end{array}$$

$$\begin{array}[]{l}{\rm Prox}_{\Xi}(u):=\arg\min_{u^{\prime}}\left\{\Xi(u^{\prime})+\frac{1}{2}\|u^{\prime}-u\|^{2}\right\},\,\,\forall\,u\in\mathcal{E}.\end{array}$$

$$\begin{array}[]{l}\nabla\Psi_{\Xi}(u)=u-{\rm Prox}_{\Xi}(u),\,\,\forall\,u\in\mathcal{E}.\end{array}$$

$$\vartheta(A)=\left\{\begin{array}[]{ll}-{\log\det}(A),&\hbox{if $A\in\mathbb{S}^{p}_{++}$};\\ +\infty,&\hbox{otherwise}.\end{array}\right.$$

$$\begin{array}[]{l}\phi^{+}_{\beta}(x):=(\sqrt{x^{2}+4\beta}+x)/{2},\,\,\phi^{-}_{\beta}(x):=(\sqrt{x^{2}+4\beta}-x)/{2},\,\,\forall\,x\in\mathbb{R}.\end{array}$$

$${\varphi }_{\beta }^{+}(A):=Q\mathrm{Diag}({\varphi }_{\beta }^{+}(d_1),\dots ,{\varphi }_{\beta }^{+}(d_p))Q^T,{\varphi }_{\beta }^{-}(A):=Q\mathrm{Diag}({\varphi }_{\beta }^{-}(d_1),\dots ,{\varphi }_{\beta }^{-}(d_p))Q^T$$

$$({\varphi }_{\beta }^{+}{)}^{\prime }(A)[B]=Q[\Gamma \odot (Q^{T}BQ)]Q^T,$$

$${\Gamma }_{ij}=({\varphi }_{\beta }^{+}(d_i)+{\varphi }_{\beta }^{+}(d_j))/(\sqrt{d_i^2+4\beta }+\sqrt{d_j^2+4\beta }),i,j=1,2,\dots ,p.$$

$$\displaystyle\left\{\begin{array}[]{l}\mathcal{W}\in\widehat{\partial}{\rm Prox}_{\mathcal{P}}(X)\hbox{ if and only if there exist $M^{(ij)}\in\widehat{\partial}{\rm Prox}_{\varphi}(X_{[ij]})$, $i<j$}\\ \hbox{such that}\ (\mathcal{W}[Y])_{[ij]}=\left\{\begin{array}[]{ll}M^{(ij)}Y_{[ij]},&\hbox{if $i<j$},\\ Y_{[ii]},&\hbox{if $i=j$},\\ M^{(ji)}Y_{[ij]},&\hbox{if $j<i$},\end{array}\right.\,\,i,j=1,\ldots,p,\,\,\,\forall\,Y\in\mathcal{Y},\end{array}\right.$$

$${\mathrm{Prox}}_{\mathcal{P}}(Y)-{\mathrm{Prox}}_{\mathcal{P}}(X)-\mathcal{W}[Y-X]=0,\forall \mathcal{W}\in \hat{\partial }{\mathrm{Prox}}_{\mathcal{P}}(Y).$$

$$\min\limits_{\Theta,\,\Omega}\Big{\{}f(\Theta,\Omega):=\sum^{L}_{l=1}\left(\vartheta(\Omega^{(l)})+\langle S^{(l)},\Theta^{(l)}\rangle\right)+\mathcal{P}(\Theta)\,\big{|}\,\Theta-\Omega=0\Big{\}}.$$

$$\mathcal{L}(\Theta ,\Omega ,X)=f(\Theta ,\Omega )-⟨\Theta -\Omega ,X⟩,\forall (\Theta ,\Omega ,X)\in \mathcal{X}\times \mathcal{X}\times \mathcal{Y}.$$

$$\begin{array}[]{cl}\max\limits_{X}&\displaystyle\sum^{L}_{l=1}\left(-\vartheta(X^{(l)})+p\right)-\mathcal{P}^{*}(X-S).\end{array}$$

$$\Theta -{\mathrm{Prox}}_{\mathcal{P}}(\Theta +X-S)=0,{\Omega }^{(l)}-{\mathrm{Prox}}_{\vartheta }({\Omega }^{(l)}-X^{(l)})=0,l=1,\dots ,L,\Theta -\Omega =0.$$

$$\begin{array}[]{l}(\Theta^{k+1},\Omega^{k+1})\approx\arg\underset{\Theta,\,\Omega}{\min}\Big{\{}f(\Theta,\Omega)+\frac{1}{2\sigma_{k}}(\|\Theta-\Theta^{k}\|^{2}+\|\Omega-\Omega^{k}\|^{2})\,|\,\Theta-\Omega=0\Big{\}},\end{array}$$

$$\begin{array}[]{l}\sup\limits_{X}\Big{\{}\Phi_{k}(X):=\inf\limits_{\Theta,\Omega}\big{\{}\mathcal{L}(\Theta,\Omega,X)+\frac{1}{2\sigma_{k}}(\|\Theta-\Theta^{k}\|^{2}+\|\Omega-\Omega^{k}\|^{2})\big{\}}\Big{\}}.\end{array}$$

$$\begin{array}[]{l}\quad\Phi_{k}(X)\\ =\inf\limits_{\Theta}\big{\{}\mathcal{P}(\Theta)+\langle\Theta,S-X\rangle+\frac{1}{2\sigma_{k}}\|\Theta-{\Theta}^{k}\|^{2}\big{\}}\\[5.69054pt] ~{}~{}+{\sum^{L}_{l=1}}\inf_{\Omega^{(l)}}\big{\{}\vartheta(\Omega^{(l)})+\langle\Omega^{(l)},X^{(l)}\rangle+\frac{1}{2\sigma_{k}}\|\Omega^{(l)}-({\Omega}^{(l)})^{k}\|^{2}\big{\}}\\[5.69054pt] =\frac{1}{\sigma_{k}}\Psi_{\sigma_{k}\mathcal{P}}({\Theta}^{k}+\sigma_{k}(X-S))+\sum^{L}_{l=1}\frac{1}{\sigma_{k}}\Psi_{\sigma_{k}\vartheta}(({\Omega}^{k})^{(l)}-\sigma_{k}X^{(l)})\\[5.69054pt] ~{}~{}-\frac{1}{2\sigma_{k}}\|{\Theta}^{k}+\sigma_{k}(X-S)\|^{2}+\frac{1}{2\sigma_{k}}\|{\Theta}^{k}\|^{2}-{\sum^{L}_{l=1}}\big{(}\frac{1}{2\sigma_{k}}\|({\Omega}^{k})^{(l)}-\sigma_{k}X^{(l)}\|^{2}-\frac{1}{2\sigma_{k}}\|({\Omega}^{(l)})^{k}\|^{2}\big{)}.\end{array}$$

$$\left\{\begin{array}[]{l}\Theta^{k+1}={\rm Prox}_{\sigma_{k}\mathcal{P}}(\Theta^{k}+\sigma_{k}(X^{k+1}-S)),\\ (\Omega^{(l)})^{k+1}={\rm Prox}_{\sigma_{k}\vartheta}\big{(}(\Omega^{(l)})^{k}-\sigma_{k}(X^{(l)})^{k+1}\big{)},\,\,l=1,2,\dots,L,\end{array}\right.$$

$$\begin{array}[]{c}X^{k+1}\approx\arg\max\limits_{X}~{}\Big{\{}\widehat{\Phi}_{k}(X):=\Phi_{k}(X)-\frac{1}{2\sigma_{k}}\|X-X^{k}\|^{2}\Big{\}}.\end{array}$$

$$(\Omega^{(l)})^{k+1}={\rm Prox}_{\sigma_{k}\vartheta}\big{(}(\Omega^{(l)})^{k}-\sigma_{k}(X^{(l)})^{k+1}\big{)}=\phi_{\sigma_{k}}^{+}\big{(}(\Omega^{(l)})^{k}-\sigma_{k}(X^{(l)})^{k+1}\big{)}.$$

$$\begin{array}[]{rl}{\rm(A)}&\|\nabla\widehat{\Phi}_{k}(X^{k+1})\|\;\leq\;\varepsilon_{k}/\sigma_{k},\,\varepsilon_{k}\geq 0,\,\sum_{k=0}^{\infty}\varepsilon_{k}<\infty;\\[7.0pt] {\rm(B)}&\|\nabla\widehat{\Phi}_{k}(X^{k+1})\|\;\leq\;(\delta_{k}/\sigma_{k})\|(\Theta^{k+1},\Omega^{k+1})-(\Theta^{k},\Omega^{k})\|,\,\,\delta_{k}\geq 0,\sum_{k=0}^{\infty}\delta_{k}<\infty.\end{array}$$

$$\Vert ({\Theta }^{k+1},{\Omega }^{k+1},X^{k+1})-(\overline{\Theta },\overline{\Omega },\overline{X})\Vert \le {\mu }_k\Vert ({\Theta }^k,{\Omega }^k,X^k)-(\overline{\Theta },\overline{\Omega },\overline{X})\Vert ,$$

$$1>{\mu }_k:=[\kappa ({\kappa }^2+{\sigma }_k^2{)}^{-1/2}+{\delta }_k]/(1-{\delta }_k)\to {\mu }_{\infty }=\kappa ({\kappa }^2+{\sigma }_{\infty }^2{)}^{-1/2}({\mu }_{\infty }=0if{\sigma }_{\infty }=\infty )$$

$$\nabla{\Phi_{k}}(X)=-{\rm Prox}_{\sigma_{k}\mathcal{P}}\big{(}U_{k}(X)\big{)}+\big{(}\phi_{\sigma_{k}}^{+}(W^{(1)}_{k}(X)),\ldots,\phi_{\sigma_{k}}^{+}(W^{(L)}_{k}(X))\big{)},$$

$$\nabla {\hat{\Phi }}_k(X)=\nabla {\Phi }_k(X)-(X-X^k)/{\sigma }_k=0.$$

$$\left\{\begin{array}[]{l}\mathcal{V}\in\widehat{\partial}(\nabla\Phi_{k})(X)\hbox{ if and only if there exists $\mathcal{G}\in\widehat{\partial}{\rm Prox}_{\mathcal{P}}(U_{k}(X)/\sigma_{k})$ such that}\\ \mathcal{V}[D]=-\sigma_{k}\mathcal{G}[D]\\ ~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\sigma_{k}\big{(}(\phi^{+}_{\sigma_{k}})^{\prime}(W_{{}_{k}}^{(1)}(X))[D^{(1)}],\dots,(\phi^{+}_{\sigma_{k}})^{\prime}(W_{k}^{(L)}(X))[D^{(L)}]\big{)},\,\,\forall D\in\mathcal{Y}.\end{array}\right.$$

$$({\mathcal{V}}_j-{\sigma }_k^{-1}I)[D]=-\nabla {\hat{\Phi }}_k(X^j)$$

$${\hat{\Phi }}_k(X^j+{\rho }^m{D}^j)\ge {\hat{\Phi }}_k(X^j)+\mu {\rho }^m⟨\nabla {\hat{\Phi }}_k(X^j),D^j⟩.$$

Full text

Keywords: Fast linear convergence $\cdot$ Network estimation $\cdot$ Semismooth Newton method $\cdot$ Sparse Jacobian

1 Introduction

Undirected graphical models have been especially popular for learning conditional independence structures among a large number of variables where the observations are drawn independently and identically from the same distribution. The Gaussian graphical model is one of the most widely used undirected graphical models. In the high-dimensional and low-sample-size settings, it is always assumed that the conditional independence structure or the precision matrix is sparse in a certain sense. In other words, its corresponding undirected graph is expected to be sparse. To promote sparsity, there has been a great deal of interest in using the $\ell_{1}$ norm penalty in statistical applications (Banerjee et al. 2008; Friedman et al. 2008; Rothman et al. 2008). In many conventional applications, a single Gaussian graphical model is typically enough to capture the conditional independence structure of the random variables. However, due to the heterogeneity or similarity of the data involved, it is increasingly appealing to fit a collection of such models jointly, such as inferring the time-varying networks and finding the change-points (Ahmed and Xing 2009; Monti et al. 2014; Gibberd and Nelson 2017; Hallac et al. 2017; Yang and Peng 2018) and estimating multiple precision matrices simultaneously for variables from distinct but related classes (Guo et al. 2011; Danaher et al. 2014; Yang et al. 2015).

Multiple graphical models refer to the models that can estimate a collection of precision matrices jointly. Specifically, let $\Delta^{(l)}$ be $L$ random vectors (from different classes or over a temporal grid) drawn independently from different distributions $\mathcal{N}_{p}(\mu^{(l)},\Sigma^{(l)}),\,l=1,2,\dots,L,\,L\geq 2$. Assume that the multivariate random variable $\Delta^{(l)}$ has $N_{l}$ observations $\delta^{(l)}_{1},\delta^{(l)}_{2},\dots,\delta^{(l)}_{N_{l}}$, for each $l\in\{1,2,\dots,L\}$. Then the sample means are $\bar{\mu}^{(l)}=\frac{1}{N_{l}}\sum_{i=1}^{N_{l}}\delta^{(l)}_{i}$ and the sample covariance matrices are $S^{(l)}=\frac{1}{N_{l}-1}\sum_{i=1}^{N_{l}}(\delta^{(l)}_{i}-\bar{\mu}^{(l)})(\delta^{(l)}_{i}-\bar{\mu}^{(l)})^{T}$, $l=1,2,\dots,L$. The multiple graphical model for estimating the precision matrices $(\Sigma^{(l)})^{-1},\,l=1,2,\dots,L$ jointly is the model with the variable $\Theta=(\Theta^{(1)},\dots,\Theta^{(L)})\in\mathbb{S}^{p}\times\cdots\times\mathbb{S}^{p}$:

[TABLE]

where $\mathcal{P}$ is a penalty function, which usually promotes sparsity in each $\Theta^{(l)}$ and similarities among different $\Theta^{(l)}$’s. Various penalties have been considered in the literature (Ahmed and Xing 2009; Guo et al. 2011; Danaher et al. 2014; Monti et al. 2014; Yang et al. 2015; Gibberd and Nelson 2017).

In this paper, we focus on the following fused graphical Lasso (FGL) regularizer which was used by Ahmed and Xing (2009) and Yang et al. (2015):

[TABLE]

We refer to problem (1) with the FGL regularizer $\mathcal{P}$ in (2) as the FGL problem. The FGL regularizer is in some sense a generalized fused Lasso regularizer (Tibshirani et al. 2005). It applies the $\ell_{1}$ penalty to all the off-diagonal elements of the $L$ precision matrices and the consecutive differences of the elements of successive precision matrices. Many elements with the same indices in the estimated matrices $\Theta^{(1)},\dots,\Theta^{(L)}$ will be close or even identical when the parameter $\lambda_{2}$ is large enough. Therefore, the FGL regularizer encourages not only shared pattern of sparsity, but also shared values across different graphs.

Existing algorithms for solving the FGL problem are quite limited in the literature. One of the most extensively used algorithms for solving this class of problems is the alternating direction method of multipliers (ADMM) (Danaher et al. 2014; Hallac et al. 2017; Gibberd and Nelson 2017). Besides, a proximal Newton-type method (Hsieh et al. 2011; Lee et al. 2014) was implemented by Yang et al. (2015) for solving the FGL problem. As we know, ADMM could be a practical first order method for finding approximate solutions of low or moderate accuracy. However, ADMM hardly utilizes any second order information, which generally must be used in order to obtain highly accurate solutions. Although the proximal Newton-type method does incorporate some forms of second order information, a complicated quadratic approximation problem has to be solved in each iteration, and this computation is usually time-consuming. It is worth mentioning that the regularizers are often introduced to promote certain structures in the estimated precision matrices, and the trade-off between biases and variances in the resulting estimators is controlled by the regularization parameters (Fan and Lv 2010). But in practice, it is extremely hard to find the optimal regularization parameters. Therefore, a sequence of regularization parameters is applied in practice, and consequently, a sequence of corresponding optimization problems must be solved (Fan and Tang 2013). Under such a circumstance, a highly efficient and robust algorithm for solving the FGL model becomes particularly important.

In this paper, we will design a semismooth Newton (SSN) based regularized proximal point algorithm (rPPA) for solving the FGL problem, which is inspired by Li et al. (2018b), where they have convincingly demonstrated the superior numerical performance of the SSN based augmented Lagrangian method (ALM), known as Ssnal, for solving the fused Lasso problem (Tibshirani et al. 2005). Thanks to the fact that the FGL problem has close connections to the fused Lasso problem, many of the virtues and theoretical insights of the Ssnal  for solving the fused Lasso problem can be observed in our approach. However, we should emphasize that solving the FGL problem is much more challenging than solving the fused Lasso problem. Specifically, the difficulties are mainly due to the log-determinant function $\log\det\,(\cdot)$ and the matrix variables, as described below.

• (a)

Unlike the simple quadratic functions in the fused Lasso problem, the function $\log\det\,(\cdot)$ is defined on the space of positive definite matrices. Therefore, the FGL model requires the positive definiteness of their solutions. This greatly increases the difficulty and complexity of theoretical analysis and numerical implementation.

• (b)

Li et al. (2018b) constructed an efficiently computable element in the generalized Jacobian of the proximal mapping of the fused Lasso regularizer, which is an essential step for solving the fused Lasso problem. Based on the constructions, we could obtain an efficiently computable generalized Jacobian of the proximal mapping of the FGL regularizer. However, this process needs more complicated manipulations of coordinates for a collection of matrix variables, unlike the vector case of the fused Lasso problem.

The key issue in the implementation of rPPA for solving the FGL model is the computation of the solution of the subproblem in each rPPA iteration. For this purpose, we will design an SSN method to solve those subproblems. We note that the numerical performance of the SSN method relies critically on the efficient calculation of the generalized Jacobian of the proximal mapping of the FGL regularizer and that of the log-determinant function. Fortunately, the generalized Jacobian of the proximal mapping of the FGL regularizer can be constructed efficiently based on that of the proximal mapping of the fused Lasso regularizer given by Li et al. (2018b). As a result, the generalized Jacobian of the proximal mapping of the FGL regularizer would inherit the structured sparsity (referred to as second order sparsity) from that of the fused Lasso regularizer. Due to the structured sparsity, the computation of a matrix-vector product in the SSN method is reasonably cheap and thus the SSN method is quite efficient for solving each subproblem. To summarize, it can be proven that our rPPA for solving the FGL problem has a linear convergent guarantee, and the convergence rate can be arbitrarily fast by choosing a sufficiently large proximal penalty parameter. Moreover, the SSN method for solving each of rPPA subproblems can be shown to be superlinearly convergent. Thus, based on these excellent convergent properties and the novel exploitation of the second order sparsity, we can expect the SSN based rPPA for solving the FGL problem to be highly efficient. Indeed, our numerical experiments have confirmed the high efficiency and robustness of the proposed algorithm for solving the FGL problems accurately.

The remaining parts of this paper are as follows. Section 2 presents some definitions and preliminary results. In section 3, we present a semismooth Newton based regularized proximal point algorithm for solving the FGL problem and its convergence properties. The numerical performance of our proposed algorithm on time-varying stock prices data sets and categorical text data sets are evaluated in section 4. Section 5 gives the conclusion.

Notations. $\mathbb{S}^{p}_{+}$ ($\mathbb{S}^{p}_{++}$) denotes the cone of positive semidefinite (definite) matrices in the space of $p\times p$ real symmetric matrices $\mathbb{S}^{p}$. For any $A,\,B\in\mathbb{S}^{p}$, we denote $A\succeq B$ $(A\succ B)$ if $A-B\in\mathbb{S}^{p}_{+}$ ($A-B\in\mathbb{S}^{p}_{++}$). In particular, $A\succeq 0$ $(A\succ 0)$ indicates $A\in\mathbb{S}^{p}_{+}$ ($A\in\mathbb{S}^{p}_{++})$. We let $\mathcal{X}:=\mathbb{S}^{p}_{+}\times\cdots\times\mathbb{S}^{p}_{+}$ and $\mathcal{Y}=:\mathbb{S}^{p}\times\cdots\times\mathbb{S}^{p}$ to be the Cartesian product of $L$ positive semidefinite cones $\mathbb{S}^{p}_{+}$ and that of $L$ spaces of symmetric matrices $\mathbb{S}^{p}$, respectively. $\mathbb{R}^{n}$ denotes the $n$-dimensional Euclidean space, and $\mathbb{R}^{m\times n}$ denotes the set of all $m\times n$ real matrices. For any $x\in\mathbb{R}^{n}$, $\|x\|_{1}=\sum_{i=1}^{n}|x_{i}|$, and $\|x\|=\sqrt{\sum_{i=1}^{n}|x_{i}|^{2}}$. We use the Matlab notation $[A;B]$ to denote the matrix obtained by appending $B$ below the last row of $A$, when the number of columns of $A$ and $B$ is identical. For any matrix $A\in\mathbb{R}^{m\times n}$, $A_{ij}$ denotes the $(i,j)$-th element of $A$. For any $X:=(X^{(1)},\dots,X^{(L)})\in\mathcal{Y}$, $X_{[ij]}:=[X^{(1)}_{ij};\ldots;X^{(L)}_{ij}]\in\mathbb{R}^{L}$ denotes the column vector obtained by taking out the $(i,j)$-th elements across all $L$ matrices $X^{(l)},\,l=1,\dots,L$. ${\rm Diag}(D_{1},\ldots,D_{n})$ denotes the block diagonal matrix whose $i$-th diagonal block are the matrix $D_{i}$, $i=1,\dots,n$. $I_{n}$ denotes the $n\times n$ identity matrix, and $I$ denotes an identity matrix or map when the dimension is clear from the context. The function composition is denoted by $\circ$, that is, for any functions $f$ and $g$, $(f\circ g)(\cdot):=f(g(\cdot))$. The Hadamard product is denoted by $\odot$.

2 Preliminaries

Let $\mathcal{E}$ be a finite-dimensional real Hilbert space, and $\Xi:\,\mathcal{E}\rightarrow{\mathbb{R}}\cup\{+\infty\}$ be a proper and closed convex function. The Moreau-Yosida regularization (Moreau 1965; Yosida 1964) of $\Xi$ is defined by

[TABLE]

The proximal mapping associated with $\Xi$ is the unique minimizer of (3) defined by

[TABLE]

Moreover, $\Psi_{\Xi}(\cdot)$ is a continuously differentiable convex function (Lemaréchal and Sagastizábal 1997; Rockafellar and Wets 2009), and its gradient is given by

[TABLE]

For notational convenience, define $\vartheta:\mathbb{S}^{p}\rightarrow\mathbb{R}\cup\{+\infty\}$ by

[TABLE]

Let $\beta>0$ be given. Define two scalar functions as follows:

[TABLE]

In addition, the matrix counterparts of these two scalar functions can be defined by

[TABLE]

for any $A\in\mathbb{S}^{p}$ with its eigenvalue decomposition $A=Q{\rm Diag}(d_{1},d_{2},\ldots,d_{p})Q^{T}$, where $d_{1}\geq d_{2}\geq\cdots\geq d_{p}$. It is easy to show that $\phi^{+}_{\beta}$ and $\phi^{-}_{\beta}$ are well-defined. Moreover, $\phi^{+}_{\beta}(A)$ and $\phi^{-}_{\beta}(A)$ are positive definite for any $A\in\mathbb{S}^{p}$.

Proposition 2.1**.**

(Wang et al. 2010, Lemma 2.1 (b))* The function $\phi^{+}_{\beta}:\,\mathbb{S}^{p}\to\mathbb{S}^{p}$ is continuously differentiable, and its directional derivative $(\phi^{+}_{\beta})^{\prime}(A)[B]$ at $A$ for any $B\in\mathbb{S}^{p}$ is given by*

[TABLE]

where $A$ admits the eigenvalue decomposition $A=Q{\rm Diag}(d_{1},d_{2},\dots,d_{p})Q^{T},\,d_{1}\geq d_{2}\geq\cdots\geq d_{p}$, and $\Gamma\in\mathbb{S}^{p}$ is defined by

[TABLE]

Proposition 2.2**.**

(Yang et al. 2013, Proposition 2.3)* For any $A\in\mathbb{S}^{p}$, it holds that \begin{array}[]{l}{\rm Prox}_{\beta\vartheta}(A)=\phi^{+}_{\beta}(A)\,\,\hbox{and}\,\,\Psi_{\beta\vartheta}(A)=-\beta\log\det(\phi^{+}_{\beta}(A))+\frac{1}{2}\|\phi^{-}_{\beta}(A)\|^{2}.\end{array}*

2.1 Surrogate Generalized Jacobian of ${\rm Prox}_{\mathcal{P}}$

In this section, we analyse the proximal mapping of the regularizer $\mathcal{P}$ defined by (2). For any $\Theta\in\mathcal{Y}$, one might observe that the penalty term $\mathcal{P}(\Theta)$ merely penalizes the off-diagonal elements, and it is the same fused Lasso regularizer that acts on each vector $\Theta_{[ij]}\in\mathbb{R}^{L},\,i\neq j$. It holds that $\mathcal{P}(\Theta)=\sum_{i\neq j}\varphi(\Theta_{[ij]})\,\,\hbox{with}\,\,\varphi(x)=\lambda_{1}\|x\|_{1}+\lambda_{2}\|Bx\|_{1},\,\forall\,x\in\mathbb{R}^{L}.$ The function $\varphi$ is the fused Lasso regularizer, and the matrix $B\in\mathbb{R}^{(L-1)\times L}$ is defined by $Bx=[x_{1}-x_{2};\dots;x_{L-1}-x_{L}]$, $\forall\,x\in\mathbb{R}^{L}$. The formula for the generalized Jacobian of ${\rm Prox}_{\varphi}$ has been derived by Li et al. (2018b) and will be used in our subsequent algorithmic design. Define the surrogate generalized Jacobian $\widehat{\partial}{\rm Prox}_{\mathcal{P}}(X):\mathcal{Y}\rightrightarrows\mathcal{Y}$ of ${\rm Prox}_{\mathcal{P}}$ at $X$ as follows:

[TABLE]

where $\widehat{\partial}{\rm Prox}_{\varphi}(\cdot)$ is the surrogate generalized Jacobian of ${\rm Prox}_{\varphi}$ (Li et al. 2018b, Equation 22) and will be given in the supplementary materials. From Theorem 1 by Li et al. (2018b), one can obtain the following theorem, which justifies why $\widehat{\partial}{\rm Prox}_{\mathcal{P}}(X)$ in (11) can be used as the surrogate generalized Jacobian of ${\rm Prox}_{\mathcal{P}}$ at $X$.

Theorem 2.1**.**

The surrogate generalized Jacobian $\widehat{\partial}{\rm Prox}_{\mathcal{P}}(\cdot)$ defined in (11) is a nonempty compact valued, upper semicontinuous multifunction. Given any $X\in\mathcal{Y}$, any element in the set $\widehat{\partial}{\rm Prox}_{\mathcal{P}}(X)$ is self-adjoint and positive semidefinite. Moreover, there exists a neighborhood $\mathcal{U}_{X}$ of $X$ such that for all $Y\in\mathcal{U}_{X}$,

[TABLE]

2.2 Lipschitz Continuity of the Solution Mapping

By introducing an auxiliary variable $\Omega=(\Omega^{(1)},\ldots,\Omega^{(L)})\in\mathcal{Y}$, we can rewrite problem (1) equivalently as

[TABLE]

The Lagrangian function of the above problem is given by

[TABLE]

The dual problem of (12) takes the following form (Borwein and Lewis 2010, Theorem 3.3.5):

[TABLE]

The Karush-Kuhn-Tucker (KKT) optimality conditions (Han et al. 2018) for (12) are given as follows:

[TABLE]

We make the following assumption on the existence of solutions to the KKT system.

Assumption 2.1**.**

The solution set to the KKT system (14) is nonempty.

Define an operator $\mathcal{T}_{\mathcal{L}}$ by $\mathcal{T}_{\mathcal{L}}(\Theta,\Omega,X):=\{(\Theta^{\prime},\Omega^{\prime},X^{\prime})\,|\,(\Theta^{\prime},\Omega^{\prime},-X^{\prime})\in\partial\mathcal{L}(\Theta,\Omega,X)\}.$ Since the function $\vartheta(\cdot)$ is strictly convex, under Assumption 2.1, we can see that the KKT system (14) has a unique KKT point, denoted by $(\overline{\Theta},\overline{\Omega},\overline{X})$, and $\mathcal{T}^{-1}_{\mathcal{L}}(0)=\{(\overline{\Theta},\overline{\Omega},\overline{X})\}$.

Proposition 2.3**.**

There exists a nonnegative scalar $\kappa$ such that for some $\varrho>0$ it holds that $\|(\Theta,\Omega,X)-(\overline{\Theta},\overline{\Omega},\overline{X})\|\leq\kappa\|\Delta\|,\,\,\forall\,\Delta\in\mathcal{T}_{\mathcal{L}}((\Theta,\Omega,X))\,\,\hbox{and}\,\,\|\Delta\|\leq\varrho.$

Proof.

Note that the regularizer $\mathcal{P}$ defined by (2) is a positive homogeneous function. Therefore, it follows from Example 11.4(a) by Rockafellar and Wets (2009) that the conjugate function $\mathcal{P}^{*}$ is an indicator function of a nonempty convex polyhedral set. This, together with Theorem 2.7 by Li et al. (2018a) and Proposition 6 by Cui et al. (2018), proves the required result. ∎

3 Regularized Proximal Point Algorithm

In this section, we present a regularized proximal point algorithm (rPPA) for solving the problem (12) with the FGL regularizer defined by (2). Given a sequence of positive scalars $\sigma_{k}\uparrow\sigma_{\infty}\leq\infty$, the $k$-th iteration of PPA for solving (12) is given by

[TABLE]

where $k\geq 0$ and $f$ is the objective function of the problem (12).

There are many ways to solve (15). Inspired by recent progresses in solving large scale convex optimization problems (Yang et al. 2013; Li et al. 2018a, b; Zhang et al. 2019), we shall adopt the approach of solving (15) via employing a sparse SSN method to its dual. The dual of (15) takes the following form:

[TABLE]

By the definition of the Moreau-Yosida regularization (3), we can write $\Phi_{k}(\cdot)$ explicitly as follows:

[TABLE]

Therefore, by Proposition 2.2 and the definition of the proximal mapping (4), the $k$-th iteration of PPA (15) can be written as

[TABLE]

where $X^{k+1}$ approximately solves the following problem: \begin{array}[]{l}X^{k+1}\approx\arg\max\limits_{X}~{}\Phi_{k}(X).\end{array} Since $\Phi_{k}(\cdot)$ is not strongly concave in general, we consider the following rPPA.

Since in practice the inner subproblem (16) can only be solved inexactly, we will use the following standard stopping criteria studied by Rockafellar (1976):

[TABLE]

The reason for using the above stopping criteria is due to the fact that Algorithm 1 is equivalent to the primal-dual PPA in the sense of Rockafellar (1976). Moreover, we have the following convergence results.

Theorem 3.1**.**

Let $\{(\Theta^{k},\Omega^{k},X^{k})\}$ be an infinite sequence generated by Algorithm 1 under stopping criterion (A). Then the sequence $\{(\Theta^{k},\Omega^{k})\}$ converges to the unique solution $(\overline{\Theta},\overline{\Omega})$ of (12), and the sequence $\{X^{k}\}$ converges to the unique solution $\overline{X}$ of (13). Furthermore, if the criterion ${\rm(B)}$ is also executed in Algorithm 1, there exists $\bar{k}\geq 0$ such that for all $k\geq\bar{k}$,

[TABLE]

where the convergence rate

[TABLE]

and the parameter $\kappa$ is from Proposition 2.3.

Proof.

The global convergence of Algorithm 1 can be obtained from Theorem 1 by Rockafellar (1976) and the uniqueness of the KKT point. The linear rate of convergence can be derived from Proposition 2.3 and Theorem 2 by Rockafellar (1976). ∎

3.1 Semismooth Newton Method for Solving Subproblem (16)

From (5), Proposition 2.2, and Theorem 31.5 by Rockafellar (2015), we know that $\widehat{\Phi}_{k}$ is a continuously differentiable, strongly concave function and

[TABLE]

where $U_{k}(X):={\Theta}^{k}+\sigma_{k}(X-S)$ and $W^{(l)}_{k}(X):=(\Omega^{k})^{(l)}-\sigma_{k}X^{(l)}$,  $l=1,2,\ldots,L$. Therefore, one can obtain the unique solution to problem (16) by solving the nonsmooth system

[TABLE]

Recall that $\phi_{\sigma_{k}}^{+}(\cdot)$ is differentiable and its derivative is given by Proposition 2.1. Thus, the surrogate generalized Jacobian $\widehat{\partial}(\nabla\Phi_{k})(X)$ of $\nabla\Phi_{k}$ at $X$ is defined as follows:

[TABLE]

With the generalized Jacobian of $\nabla\Phi_{k}$, we are ready to solve equation (17) by the SSN method, where the Newton systems are solved inexactly by the conjugate gradient (CG) method.

Next, we derive the convergence result of the SSN method (Algorithm 2).

Theorem 3.2**.**

Let $\{X^{j}\}$ be the infinite sequence generated by Algorithm 2. Then $\{X^{j}\}$ converges to the unique optimal solution $\widehat{X}$ of (17) and $\|X^{j+1}-\widehat{X}\|=\mathcal{O}(\|X^{j}-\widehat{X}\|^{1+\tau})$.

Proof.

Since the proximal mapping ${\rm Prox}_{\mathcal{P}}$ is piecewise linear and Lipschitz continuous, we know from Theorem 7.5.17 by Facchinei and Pang (2007) that ${\rm Prox}_{\mathcal{P}}$ is directionally differentiable. This, together with Theorem 2.1, implies that ${\rm Prox}_{\mathcal{P}}$ is strongly semismooth with respect to the multifunction $\widehat{\partial}{\rm Prox}_{\mathcal{P}}$ (for its definition, see e.g., Definition 1 by Li et al. (2018b)). Therefore, the conclusion follows from the strong convexity of $\widehat{\Phi}_{k}(\cdot)$, Proposition 2.1, and Proposition 7 & Theorem 3 by Li et al. (2018b). ∎

4 Numerical Experiments

In this section, we compare the performance of our algorithm rPPA with the alternating direction method of multipliers (ADMM) and the proximal Newton-type method implemented by Yang et al. (2015) (referred to as MGL here) for which the solver is available at http://senyang.info/.

The following paragraph describes the measurement of the accuracy of an approximate optimal solution and the stopping criteria of the three methods. Since both rPPA and ADMM can generate primal and dual approximate solutions, we can assess the accuracy of their solutions by the relative KKT residuals. Unlike the primal-dual method, MGL merely gives the primal solution and the KKT residual of a solution generated by MGL is not available. Instead, we measure the relative error of the objective value obtained by MGL with respect to that computed by rPPA. Based on the KKT optimality condition (14), the accuracy of an approximate optimal solution $(\Theta,\Omega,X)$ generated by rPPA (Algorithm 1) is measured by defining the following relative residuals:

[TABLE]

Likewise, the accuracy of an approximate optimal solution $(\Theta,X,Z)$ generated by ADMM is measured by the relative KKT residual $\eta_{A}$ (defined in the supplementary material) that is analogous to $\eta_{P}$.

In our numerical experiments, we terminate rPPA if it satisfies the condition $\eta_{P}<\varepsilon$ for a given accuracy tolerance $\varepsilon$; similarly for ADMM with the stopping condition $\eta_{A}<\varepsilon.$ Note that the terminating condition for MGL is different. Let “${\rm pobj}_{P}$” and “${\rm pobj}_{M}$” be the primal objective function values computed by rPPA and MGL, respectively. MGL will be terminated when the relative difference of its objective value with respect to that obtained by rPPA is smaller than the given tolerance $\varepsilon$, i.e.,

[TABLE]

We adopt a warm-start strategy to initialize rPPA. That is, we first run ADMM (with identity matrices as the starting point) for a fixed number of iterations to generate a good initial point to warm-start rPPA. We also stop ADMM as soon as the relative KKT residual of the computed iterate is less than $100\varepsilon$. Note that such a warm-starting strategy is sound since in the initial phase of rPPA where the iterates are not close to the optimal solution (as measured by the associated relative KKT residual), it is computationally wasteful to use the more expensive rPPA iteration when the fast local linear convergence behavior of the algorithm has yet to kick in. Under such a scenario, naturally one would use cheaper iterations such as those of ADMM to generate the approximate solution points until the relative KKT residual has been sufficiently reduced.

For the tuning parameters $\lambda_{1}$ and $\lambda_{2}$, for each test instance we select three pairs that lead to reasonable sparsity. In the following tables, “P” stands for rPPA; “A” stands for ADMM; “M” stands for MGL; “nnz” denotes the number of nonzero entries in the solution $\Theta$ obtained by rPPA using the estimation: ${\rm nnz}:=\min\{k\,|\,\sum_{i=1}^{k}|\hat{x}_{i}|\geq 0.999\|\hat{x}\|_{1}\},$ where $\hat{x}\in\mathbb{R}^{p^{2}L}$ is the vector obtained via sorting all elements in $\Theta$ by magnitude in a descending order; “density” denotes the quantity nnz$/(p^{2}L)$. The time is displayed in the format of “hours:minutes:seconds”, and the fastest method in terms of running time is highlighted in red. The errors presented in the tables are the relative KKT residuals $\eta_{P}$ for rPPA and $\eta_{A}$ for ADMM; while the error for MGL is $\Delta_{M}$ in (19).

4.1 Nearest-neighbour Networks

In this section, we assess the effectiveness of the FGL model on a simulated network: nearest-neighbour network. The nearest-neighbour network is generated by modifying the data generation mechanism described by Li and Gui (2006). We set $p=500$ and $L=3$. For each $l=1,2,\dots,L$, we generate 10,000 independently and identically distributed observations from a multivariate Gaussian distribution $\mathcal{N}_{p}(0,(\Omega^{(l)})^{-1})$, where $\Omega^{(l)}$ is the precision matrix of the $l$-th class. The details of the generation of $\Omega^{(l)}$ are as follows. First of all, $p$ points are randomly generated on a unit square, their pairwise distances are calculated, and $m$-nearest neighbours of each point in terms of distance are found. The nearest-neighbour network is obtained by linking any two points that are $m$-nearest neighbours of each other. The integer $m$ controls the degree of sparsity of the network, and we set $m=5$ in our simulation. Subsequently, we add heterogeneity to the common structure by further creating individual links as follows: for each $\Omega^{(l)}$, a pair of symmetric zero elements is randomly selected and replaced with a value uniformly drawn from the interval $[-1,-0.5]\cup[0.5,1]$. This procedure is repeated ceil$(M/4)$ (the nearest integer greater than or equal to $M/4$) times, where $M$ is the number of edges in the nearest-neighbour graph. In our simulation, the true number of edges in the three networks is $3690$.

There is a pair of tuning parameters $\lambda_{1}$ and $\lambda_{2}$ which must be specified. In the FGL model, $\lambda_{1}$ drives sparsity and $\lambda_{2}$ drives similarity, and we say that $\lambda_{1}$ and $\lambda_{2}$ are the sparsity and similarity control parameters respectively. In order to show the diversity of sparsity in our experiments, we choose a series of $\lambda_{1}$ for the FGL model with $\lambda_{2}$ fixed. Figure 1 shows the relative ability of the FGL model to recover the network structures and to detect the change-points.

Figure 1(a) displays the number of true positive edges selected (i.e., TP edges) against the number of false edges selected (i.e., FP edges) for the FGL model. We say that an edge $(i,j)$ in the $l$-th network is selected in the estimate $\widehat{\Theta}^{(l)}$ if $\widehat{\Theta}^{(l)}_{ij}\neq 0$, and we say that the edge is true in the precision matrix $(\Sigma^{(l)})^{-1}$ if $((\Sigma^{(l)})^{-1})_{ij}\neq 0$ and false if $((\Sigma^{(l)})^{-1})_{ij}=0$. We can see from the figure that the FGL model with $\lambda_{2}=0.005$ can recover almost all of the true positive edges without false positive edges. Figure 1(a) also shows that for the FGL model the similarity control parameter $\lambda_{2}=0.005$ is much better than $\lambda_{2}=0.05$ in terms of the ability of true edges detection. When $\lambda_{2}=0.05$, the FGL model can merely detect about 3000 true positive edges while the the number of false positive edges is increased to over 600. One possible reason is that $\lambda_{2}=0.05$ is too large compared with the underlying optimal one in this case.

Figure 1(b) illustrates the sum of squared errors between estimated edge values and true edge values, i.e., \sum_{l=1}^{L}\sum_{i<j}\big{(}\widehat{\Theta}^{(l)}_{ij}-((\Sigma^{(l)})^{-1})_{ij}\big{)}^{2}. When the number of the total edges selected is increasing (i.e., the sparsity control parameter is decreasing), the error is decreasing and finally reaches a fairly low value.

Figure 1(c) plots the number of true positive differential edges against false positive differential edges. A differential edge is an edge that differs between classes and thus corresponds to a change-point. We say that the $(i,j)$ edge is estimated to be differential between the $l$-th and the $(l+1)$-th networks if $|\widehat{\Theta}^{(l)}_{ij}-\widehat{\Theta}^{(l+1)}_{ij}|>10^{-6}$, and we say that it is truly differential if $|((\Sigma^{(l)})^{-1})_{ij}-((\Sigma^{(l+1)})^{-1})_{ij}|>10^{-6}$. The number of differential edges is computed for all successive pairs of networks. The best point in Figure 1(c) is the red one which has approximately 2700 true positive differential edges and almost no false one. We can also see from Figure 1(c) that all the blue points have no false positive differential edge and small numbers of true positive differential edges. This might be caused by the larger similarity control parameter $\lambda_{2}=0.05$ which forces an excessive number of edges across $L$ networks to be similar.

4.2 Standard & Poor’s 500 Stocks

In this section, we compare rPPA, ADMM, and MGL on the Standard & Poor’s 500 stock price data sets. The stock price data sets contain daily returns of 500 stocks over a long period, and can be downloaded from the link www.yahoo.com. The dependency structures of different stocks vary over time. But it appears that the dependency networks change smoothly over time. Therefore, the FGL model might be able to find the interactions among these stocks and how they evolve over time.

We first consider a relatively short three-year time period from January 2004 to December 2006. During this period, there are totally 755 daily returns of 370 stocks. We call this data set SPX3a. For each year, it contains approximately 250 daily returns of each stock. Considering the limited number of observations in each year and the interpretation of the results, we choose to analyse random smaller subsets of all involved stocks, whose sizes are chosen to be $p=100$ and $p=200$, over $L=3$ periods.

In addition to the above data set over three years, a relatively long period from January 2004 to December 2014 is also considered in the experiments, which is referred to as SPX11b. Since the time period is longer than the previous one, the number of stocks becomes smaller as some stocks might disappear. During the 11-year time period, there are 2769 daily closing prices of 272 stocks. We can set a relatively large parameter $L=11$ according to years from January 2004 to December 2014. Again, we choose to analyse two random subsets of all existing stocks, of which the sizes are selected to be $p=100$ and $p=200$.

Table 1 shows the comparison of rPPA, ADMM, and MGL on the stock price data sets SPX3a and SPX11b with $100$ and $200$ selected stocks. One outstanding observation from the table is that rPPA outperforms ADMM and MGL by an obvious margin and rPPA is faster than the other two methods except for one instance. For the exceptional instance, rPPA is still faster than ADMM and merely two seconds slower than MGL. In addition, we find that both rPPA and ADMM succeeded in solving all instances; while MGL failed to solve two of them within one hour. This might imply that MGL is not robust for solving the FGL model when applied to the stock price data sets. The numerical results show convincingly that our algorithm rPPA can solve the FGL problems efficiently and robustly. The superior performance of rPPA can mainly be attributed to our ability to extract and exploit the sparsity structure (in the surrogate generalized Jacobian of ${\rm Prox}_{\mathcal{P}}$) within the SSN method to solve each rPPA subproblem very efficiently.

Figure 2 displays the sparse patterns of $11$ estimated precision matrices from year 2004 to year 2014 on the data set SPX11b with $(\lambda_{1},\lambda_{2})=$ (1e-4,1e-4) in application of the FGL model. It should be noted that this period covers the 2008 financial crisis. We manually split the time points into three stages (one stage corresponds to one row in Figure 2) to aid the interpretation of the results. Each red pattern in the left panel presents the common structure across the estimated precision matrices in its stage. And each blue pattern visualizes the individual edges specific to its own precision matrix. Generally, one can hardly expect a meaningful common structure across all the $11$ time points, and thus we provide here the common structure across parts of nearby precision matrices. One can clearly see that more edges are detected in the middle stage, and the number of the common edges across year 2008, 2009, 2010, 2011 is correspondingly larger than that in the earlier and later stages. The increased amount of interactions among the stocks over this period is likely due to the 2008 global financial crisis and its sustained effects. Another observation is that the number of edges had a drastic increase in 2007, kept at a high level during the 2008 global financial crisis and a certain period after that, and then went down to a level still higher than that of the pre-crisis period (year 2004, 2005, 2006). The sudden increase in 2007 might be seen as a prediction of the oncoming financial crisis in 2008. The increased amount of interactions among stocks after the financial crisis compared to that in the pre-crisis period may indicate some essential changes of the financial landscape. To a certain degree, the observations agree well with those observed by Yang and Peng (2018).

4.3 University Webpages

Here we evaluate the numerical performances of rPPA, ADMM, and MGL on the data set university webpages, which is provided by Cardoso-Cachopo (2007) and available at http://ana.cachopo.org/datasets-for-single-label-text-categorization. The original pages were collected from computer science departments of various universities in 1997. We selected four largest and meaningful classes in our experiment: Student, Faculty, Staff, and Department. For each class, the collection contains pages from four universities: Cornell, Texas, Washington, Wisconsin, and other miscellaneous pages from other universities. Furthermore, the original text data have been preprocessed by stemming techniques, that is, reducing words to their morphological roots. The preprocessed data sets downloaded from the link above contain two files: two thirds of the pages were randomly chosen as training set (Webtrain) and the remaining third as testing set (Webtest). Table 2 presents the distribution of documents per class.

Next, we apply the FGL model to the ${\it Webtest}$ data set for the purpose of interpreting the data. We choose tuning parameters that enforce high sparsity and similarity. In our experiment, we set $\lambda_{1}=0.005$ and $\lambda_{2}=0.003$. The resulting common structure is displayed in Figure 3. The thickness of an edge is proportional to the magnitude of the associated average partial correlation. Figure 3 shows that some standard phrases in computer science, such as program-languag, oper-system, distribut-system, softwar-engin, possess high partial correlations among their constituent words in all four classes. It successfully demonstrates the effectiveness of the FGL model for exploring the similarity across related classes. Table 3 shows the comparison of the three methods rPPA, ADMM, and MGL on the webpages data sets with data dimension $p=100$, $p=200$, and $p=300$. As can be seen, rPPA outperforms ADMM and MGL by a large margin for most of the tested webpages data sets.

5 Conclusion

We have designed an efficient and globally convergent regularized proximal point algorithm for solving the primal formulation of the fused graphical Lasso problem. From a theoretical perspective, we established the Lipschitiz continuity of the solution mapping and consequently obtained that the primal and dual sequences are locally linearly convergent. This lays the foundation for the efficiency of the proposed algorithm. Moreover, the second order information was also fully exploited, which further leads to the high efficiency of the proposed algorithm. Numerically, we demonstrated the superior efficiency and robust performance of the proposed method by comparing it with the extensively used alternating direction method of multipliers and the proximal Newton-type method (Yang et al. 2015) on both synthetic and real data sets. In summary, the proposed semismooth Newton based regularized proximal point algorithm is a highly efficient method for solving the fused graphical Lasso problems.

6 Supplementary Materials

Supplementary material:

It contains technical details (generalized Jacobian of the proximal mapping of the fused Lasso regularizer and implementation of ADMM) and numerical results (on data University Webpages and 20 Newsgroups). (pdf file)

Appendix 1 Generalized Jacobian of the Proximal Mapping of the Fused Lasso Regularizer

In this section, we recall the characterization of the generalized Jacobian of the fused Lasso regularizer (Tibshirani et al. 2005), which will be used to derive the explicit expression of the generalized Jacobian of the fused graphical Lasso (FGL) regularizer.

The fused Lasso regularizer is defined by $\varphi(x)=\lambda_{1}\|x\|_{1}+\lambda_{2}\|Bx\|_{1},\,\forall\,x\in\mathbb{R}^{L},$ where the matrix $B\in\mathbb{R}^{(L-1)\times L}$ is defined by $Bx=[x_{1}-x_{2};\dots;x_{L-1}-x_{L}]$. Denote the proximal mapping of $\lambda_{2}\|B\cdot\|_{1}$ by $x_{\lambda_{2}}(v):=\arg\min_{x}\left\{\lambda_{2}\|Bx\|_{1}+\frac{1}{2}\|x-v\|^{2}\right\},\,\,\forall\,v\in\mathbb{R}^{L}.$

Lemma 1.1**.**

(Friedman et al. 2007, Proposition 1)* Given $\lambda_{1},\lambda_{2}\geq 0$, it holds that ${\rm Prox}_{\varphi}(v)={\rm Prox}_{\lambda_{1}\|\cdot\|_{1}}(x_{\lambda_{2}}(v))={\rm sign}(x_{\lambda_{2}}(v))\odot\max(|x_{\lambda_{2}}(v)|-\lambda_{1},0),\,\,\forall\,v\in\mathbb{R}^{L}.$*

Lemma 1.2**.**

(Li et al. 2018, Lemma 1)* Given $\lambda_{2}\geq 0$, it holds that $x_{\lambda_{2}}(v)=v-B^{T}z_{\lambda_{2}}(Bv),\,\,\forall\,v\in\mathbb{R}^{L},$ where $z_{\lambda_{2}}(u):=\arg\min\limits_{z}\left\{\frac{1}{2}\|B^{T}z\|^{2}-\langle z,\,u\rangle~{}|\,\,\|z\|_{\infty}\leq\lambda_{2}\right\},\,\,\forall\,u\in\mathbb{R}^{L-1}.$*

Given $v\in\mathbb{R}^{L}$, consider the following sets:

[TABLE]

Define the multifunction $\mathcal{Q}_{z}:\mathbb{R}^{L}\rightrightarrows\mathbb{R}^{(L-1)\times(L-1)}$ by

[TABLE]

where $\Sigma_{K}={\rm Diag}(\sigma_{K})\in\mathbb{R}^{(L-1)\times(L-1)}$ with $(\sigma_{K})_{i}=0$, if $i\in K$, $(\sigma_{K})_{i}=1$, if $i\notin K$. Define also the multifunction $\mathcal{Q}_{x}:\,\mathbb{R}^{L}\rightrightarrows\mathbb{R}^{L\times L}$ by

[TABLE]

It has been shown by Li et al. (2018) that the surrogate generalized Jacobian of ${\rm Prox}_{\varphi}$ at $v$ is the set

[TABLE]

where $\partial_{B}{\rm Prox}_{\lambda_{1}\|\cdot\|_{1}}$ denotes the B-subdifferential of ${\rm Prox}_{\lambda_{1}\|\cdot\|_{1}}$ (Qi 1993, Equation (2.12)).

Appendix 2 Implementation of ADMM

In this part, we briefly describe the alternating direction method of multipliers (ADMM) for solving the dual problem of the FGL problem:

[TABLE]

This can be rewritten equivalently as follows:

[TABLE]

The augmented Lagrangian function associated with (1), given $\sigma>0$, is defined by

[TABLE]

The KKT optimality conditions are as follows:

[TABLE]

where $\vartheta(X)=-\log\det\,X$ if $X\in\mathbb{S}^{p}_{++}$ and $\vartheta(X)=+\infty$ otherwise. Due to its separable structure in terms of the variables $X$ and $Z$, ADMM is often considered as a natural choice for solving (1). The classic ADMM was first proposed by Glowinski and Marroco (1975); Gabay and Mercier (1976), and later extended by Fazel et al. (2013); Chen et al. (2017). The iteration scheme of ADMM for (1) can be described as follows: given $\tau\in(0,(1+\sqrt{5})/2)$, and an initial point $(X^{0},Z^{0},\Theta^{0})$, the $(k+1)$-th iteration is given by

[TABLE]

Here, $X^{k+1}=((X^{(1)})^{k+1},\dots,(X^{(L)})^{k+1})$ can be updated by

[TABLE]

In our implementation, we tune the parameter $\sigma$ wisely according to the progress of primal and dual feasibilities (Lam et al. 2018, Section 4.4). We also use a larger step-length $\tau$ of $1.618$, which has been demonstrated in various works to perform better than the simple case with $\tau=1$. Based on the KKT optimality conditions (2), the accuracy of an approximate optimal solution $(\Theta,X,Z)$ generated by ADMM is measured by the following relative residual:

[TABLE]

Thus far, we have provided an easily implementable framework of ADMM for which each iteration requires the computation of the proximal mapping of the log-determinant function and that of the FGL regularizer.

Appendix 3 Numerical Experiment: University Webpages

This section presents the procedure of processing the university webpages data (available at http://ana.cachopo.org/datasets-for-single-label-text-categorization) and generating sample covariance matrices which is similar to the process used by Guo et al. (2011). Actually, the previous work (Guo et al. 2011) on the data set Webtest applied a different penalty term to estimate multiple graphical models jointly. For given integer $p$, the sample covariance matrices $S^{(l)},\,l=1,2,3,4$ were constructed from the data set Webtest in the following way:

• (i)

Choose $p$ words with highest frequency which appear in each class at least once. Namely, the words we analyse are a subset of all involved words.

• (ii)

Obtain $X^{(1)}\in\mathbb{R}^{544\times p}$ from class Student, where the $(i,j)$-th element $X^{(1)}_{ij}$ denotes the number of times the $j$-th term appears in the $i$-th page of class Student. In the same way, $X^{(2)}\in\mathbb{R}^{374\times p}$, $X^{(3)}\in\mathbb{R}^{310\times p}$, and $X^{(4)}\in\mathbb{R}^{168\times p}$ can be obtained from class Faculty, Course, and Project, respectively. Denote their vertical concatenation by a new matrix $X=[X^{(1)};X^{(2)};X^{(3)};X^{(4)}]\in\mathbb{R}^{1396\times p}$.

• (iii)

The matrix $P$ is obtained by normalizing $X$ along each column: $P_{ij}=X_{ij}/\sum_{i}X_{ij}$. Then, the log-entropy weight of the $j$-th word is defined as $e_{j}=1+\sum_{i}P_{ij}(\ln P_{ij})/\ln 1396$.

• (iv)

Compute $\overline{X}$ as follows: $\overline{X}_{ij}=e_{j}\ln(1+X_{ij})$, and split $\overline{X}$ by columns accordingly: $\overline{X}=[\overline{X}^{(1)};\overline{X}^{(2)};\overline{X}^{(3)};\overline{X}^{(4)}]$.

• (v)

Generate sample covariance matrices $S^{(l)}$ from $\overline{X}^{(l)}$: $S^{(1)}={\rm cov}(\overline{X}^{(1)}),\,l=1,2,3,4$.

Following the procedure described above, we can also generate sample covariance matrices from the data set Webtrain.

In section 4.3 of the paper, we have successfully demonstrated the effectiveness of the FGL model for exploring the similarity across related classes (Figure 3). On the other hand, we believe that the model can also detect the heterogeneity among different classes. As an example, supplementary Figures 1 and 2 illustrate the differences between the Course (Figure 1) and Project (Figure 2) classes. One can see that some course related terms, such as class and assign, are of high degree in supplementary Figure 1; whereas they are not even connected in supplementary Figure 2. Besides, some teaching related terms are linked only in the Course category, such as class-assign, assign-problem, class-project. Overall, it is likely that the FGL model is capable of identifying the common and individual structures of the webpages among related classes.

Supplementary Figure 3 presents the performance profiles of rPPA, ADMM, and MGL for all 18 tested problems. The meaning of the performance profiles is given as follows: a point $(x,y)$ is on the performance curve of a particular method if and only if this method can solve up to desired accuracy $(100y)\%$ of all the tested instances within at most $x$ times of the fastest method for each instance. As can be seen, rPPA outperforms ADMM and MGL by a large margin for the most of the tested webpages data sets. In particular, focusing on $y=40\%$, we can see that rPPA is around $3\sim 5$ times faster in comparison with ADMM and MGL for over $60\%$ of the tested instances.

Appendix 4 Numerical Experiment: 20 Newsgroups

This section compares rPPA, ADMM, and MGL on newsgroups data set, which is a popular text data set. The 20 newsgroups data set is a collection of newsgroup documents, partitioned nearly evenly across $20$ different newsgroups. Different newsgroups correspond to different topics, and some of the newsgroups are closely related to each other (e.g., comp.sys.ibm.pc.hardware/comp.sys.mac.hardware), while others are highly unrelated (e.g., misc.forsale/soc.religion.christian). Supplementary Table 1 lists the 20 newsgroups (from the website: http://qwone.com/~jason/20Newsgroups/), partitioned according to subject matter. According to the topics, our numerical experiments were conducted on the four subgroups, which are likely to possess common semantic structures. The four subgroups are highlighted in supplementary Table 1 and named as NGcomp, NGrec, NGsci, and NGtalk accordingly. There are several classes in each subgroup, and we apply the FGL model to estimating jointly the precision matrices of different classes in each subgroup.

A processed version of the 20 newsgroups data set which is easy to read into Matlab can be downloaded from Jason’s page http://qwone.com/~jason/20Newsgroups/, and the downloaded data contains a training data set and a testing data set. We also adopted the procedure of generating sample covariance matrices described in the previous section 3 with a series of problem dimensionality $p=100$, $p=200$, and $p=300$.

Supplementary Table 2 shows the comparison of rPPA, ADMM, and MGL on the testing and training data sets of four subgroups with parameters $p=300$. The results for $p=200$ and $p=100$ are not shown by tables here for lack of space. Instead, we summarize all conducted instances (with different dimensionality $p=100$, $p=200$, and $p=300$, with various tuning parameters $(\lambda_{1},\lambda_{2})$) in Supplementary Figure 4. One can clearly see that rPPA outperforms ADMM and MGL by an obvious margin. It truly suggests that our proposed algorithm is efficient for solving the FGL problems.