This work was supported by the Grants-in-Aid for Scientific Research (KAKENHI)
under Grant JP18H01446.

Corresponding author: Masahiro Yukawa (e-mail: yukawa@elec.keio.ac.jp).

# Learning Sparse Graph with Minimax Concave Penalty under Gaussian Markov
Random Fields

TATSUYA KOYAKUMARU1  MASAHIRO YUKAWA1  EDUARDO PAVEZ2  and ANTONIO ORTEGA.2
Department of Electronics and Electrical Engineering, Keio University,
Kanagawa 223-8522, Japan  Department of Electrical and Computer Engineering,
University of Southern California, Los Angeles, California, CA 90089, USA

###### Abstract

This paper presents a convex-analytic framework to learn sparse graphs from
data. While our problem formulation is inspired by an extension of the
graphical lasso using the so-called combinatorial graph Laplacian framework, a
key difference is the use of a nonconvex alternative to the $\ell_{1}$ norm to
attain graphs with better interpretability. Specifically, we use the weakly-
convex minimax concave penalty (the difference between the $\ell_{1}$ norm and
the Huber function) which is known to yield sparse solutions with lower
estimation bias than $\ell_{1}$ for regression problems. In our framework, the
graph Laplacian is replaced in the optimization by a linear transform of the
vector corresponding to its upper triangular part. Via a reformulation relying
on Moreau’s decomposition, we show that overall convexity is guaranteed by
introducing a quadratic function to our cost function. The problem can be
solved efficiently by the primal-dual splitting method, of which the
admissible conditions for provable convergence are presented. Numerical
examples show that the proposed method significantly outperforms the existing
graph learning methods with reasonable CPU time.

###### Index Terms:

Graph signal processing, graph learning, graphical lasso, minimax concave
penalty, primal-dual splitting method, proximity operator

=-15pt

##  I Introduction

How can we learn sparse graphs with enhanced interpretability under the
Gaussian Markov random field (GMRF) [1]? This is the central question
addressed in this article. A graph, containing a set of vertices and edges, is
a mathematical tool to represent the dependencies among components (such as
nodes of a network or pixels of an image), through the selection of pairwise
relations (edge weights) between each pair of objects (vertices). In
particular, the strength of the relation can be expressed in terms of
(nonnegative) graph weights. In the present context, graph “sparseness” is an
important property because it tends to provide better interpretability, i.e.,
relative to all possible connections between nodes, only a few edges are non-
zero and provide information about the major relationships between objects.

The problem of learning graphs from data has been studied widely in a variety
of fields including signal processing, machine learning, and statistics [2, 3,
4, 5, 6, 7, 8, 9, 10, 11, 12]. Graph learning has been considered in multiple
applications such as design of functional brain network architectures [13],
molecular biology [14], and network anomaly detection [15]. We also refer the
reader to [16, 17] for comprehensive reviews of graph learning. The graphical
model approach [2, 3, 4, 5, 6] represents dependencies with the data in graph
form and has gained significant popularity owing to two main reasons. First,
the graphical model is built upon a solid statistical foundation, so that the
edge weights have a physical meaning under certain assumptions. For instance,
if the observed data are derived from a GMRF model, the weights are based on
partial correlation coefficients [18]. Second, it provides excellent
versatility as it assumes no specific structure on the graph. A particular
example of graph learning algorithm is the graphical lasso [2, 19, 7, 8, 9],
which employs $\ell_{1}$ regularization on the edge weights to obtain the
sparse inverse covariance matrix of a GMRF model [1]. This approach has been
extended and modified in [10] to learn several types of Laplacian matrices,
including a formulation where the inverse covariance matrix has a
combinatorial graph Laplacian (CGL) structure [10]. As noted earlier, research
on sparse graph learning is motivated by the fact that sparsity enhances the
interpretability of the learned graphs [20, 10]. All those sparsity-seeking
methods exploit convex penalties (such as the $\ell_{1}$ norm) mainly due to
their mathematical tractability.

To clarify the motivation of the present study, let us turn our attention to
sparse linear regression. A plethora of nonconvex alternatives to the
$\ell_{1}$ regularization have been proposed to reduce the estimation bias
while maintaining the benefit of variance reduction [21, 22, 23, 12]. Among
them, we focus on the minimax concave (MC) penalty [21, 24] because:

  1. (i)

it saturates (i.e., it returns a constant value when the variable being
estimated exceeds a given threshold) thereby reducing the estimation biases
significantly;

  2. (ii)

it has been shown to bridge the gap between the $\ell_{0}$ and $\ell_{1}$
norms in a parametric way [25];

  3. (iii)

it is a weakly convex function [26]; more specifically, it is given by
subtracting from the $\ell_{1}$ norm its Moreau envelope.

Property (iii) is of particular importance from an optimization viewpoint
because the overall convexity of the cost function is ensured when the MC
penalty is used with strongly convex loss functions, and also because the
decomposed form in terms of a convex function and its Moreau envelope is
compatible with the efficient operator splitting methods. The MC penalty has
been used in various sparse estimation problems, e.g., feature selection with
a sparse support vector machine (SVM) [27] and gear fault diagnosis from noisy
vibration signals [28]. To the best of authors’ knowledge, the $\ell_{p}$
quasi-norm for $p\in(0,1)$ is the only function, excluding the MC penalty,
that is known to possess property (ii) above, but it lacks properties (i) and
(iii). On the other hand, the smoothly clipped absolute deviation (SCAD)
penalty [12] is similar to the MC penalty and it could be an alternative
choice, since it actually possesses property (i) as well as the weak-convexity
part of property (iii) although further investigations would be needed to
determine whether property (ii) and the other part of property (iii) also hold
for the SCAD penalty. While nonconvex alternatives to the $\ell_{1}$ penalty
have been successful in the context of sparse linear regression, their study
for graph learning has been limited, and most of the existing graph learning
methods use the $\ell_{1}$ or $\ell_{2}$ regularization. A few exceptions
include approaches using the log function [11] or the $\ell_{0}$ norm [29]. A
use of the SCAD penalty was also mentioned in [12, 30]. The recent works [30,
31] have observed that an increase of the regularization parameter of the
$\ell_{1}$ penalty in the CGL estimation framework ultimately does not lead to
a sparse solution and instead produces a dense solution associated with a
fully connected graph. Based on this observation, in [30, 31], it has been
shown that the use of the MC penalty (as well as other nonconvex penalties)
yields better performance. However, these approaches are based on a nonconvex
formulation and thus there is no guarantee that the generated sequence of
graphs converges to a global optimum. This motivates us to devise another
formulation which benefits from the weak convexity of the MC penalty to
guarantee overall convexity of the entire cost so that the generated graphs
converge to provably global optimum.

The goal of this article is to present a novel graph learning framework based
on a convex formulation involving the nonconvex (but weakly convex) MC penalty
to produce sparse graphs, and specifically sparse CGL matrices. Since CGLs are
symmetric matrices, we remove this redundancy by representing a CGL matrix
using a linear transform of the vector of graph weights corresponding to the
upper-triangular part.111Although the one-to-one linear operator for
representing the CGL was used in the structured graph learning via Laplacian
spectral constraints (SGL) method [11], the proposed method is more efficient
(as shown in Section IV) due to the proposed reformulation, which allows to
use the primal-dual splitting method [32] (as shown in Section III). Here, the
upper-triangular part represents the undirected relations among nodes and
completely characterizes the CGL matrix, so that our estimate is automatically
guaranteed to have a Laplacian structure without the need to impose any
constraints. This is in sharp contrast to the existing CGL approaches, which
typically require both a positive semi-definite constraint and a linear
constraint. Our formulation involves the nonconvex MC penalty, instead of the
$\ell_{1}$ norm, while essentially keeping the same terms (the “nonsmooth”
log-determinant term and the linear term) as the graphical lasso formulation
but with the linear operator mentioned above. Note here that the negative log-
determinant function is differentiable but with non-Lipschitz-continuous
gradient. Due to the nonconvexity and the nonsmoothness, the problem cannot be
solved directly using existing optimization methods. To circumvent the
difficulty, we invoke the classical Moreau’s decomposition and show that the
Tikhonov regularization convexifies the overall cost function, reformulating
the problem into a canonical form of the primal-dual splitting method [32]. We
present the admissible conditions under which the convergence to the global
optimal point is guaranteed by the primal-dual splitting method. Numerical
examples show that the proposed method outperforms the conventional CGL method
(its $\ell_{1}$-based counterpart) for three types of graph. Compared to the
state-of-the-art method, the structured graph learning via Laplacian spectral
constraints (SGL) [11], the proposed method achieves comparable or better
performance, depending on the type of graph, with up to 40 times shorter CPU
time. In addition, experiments with real data show that the method produces a
sparser graph than other existing methods.

New features of the present work relative to our preliminary work [33] include
detailed proofs of the mathematical results and refined experimental results
as well as additional simulation results using real data.

##  II PRELIMINARIES

We present notation, and then show some mathematical tools used in this work.
We finally present the primal-dual splitting method which is used to solve the
proposed optimization problem to be presented in Section III.

###  II-A Notation

The sets of real numbers and nonnegative real numbers are denoted by
$\mathbb{R}$ and $\mathbb{R}_{+}$, respectively. The transpose of
vector/matrix is denoted by $(\cdot)^{T}$. Given a vector
$\bm{x}:=[x_{1},x_{2},\cdots,x_{n}]^{T}$ $\in\mathbb{R}^{n}$, define the
$\ell_{1}$ and the $\ell_{2}$ norms by $\|\bm{x}\|_{1}$ := $\sum_{i=1}^{n}$
$|x_{i}|$ and $\|\bm{x}\|_{2}$ :=
$\left(\sum_{i=1}^{n}x_{i}^{2}\right)^{\frac{1}{2}}$, respectively. Similarly,
given a matrix $\bm{X}\in\mathbb{R}^{n\times n}$ with its $(i,j)$ component
denoted by $x_{i,j}$, define the $\ell_{1}$ and the Frobenius norms by
$\|\bm{X}\|_{1}:=\sum_{i,j=1}^{n}|x_{i,j}|$ and $\|\bm{X}\|_{\rm
F}:=\left(\sum_{i,j=1}^{n}x_{i,j}^{2}\right)^{\frac{1}{2}}$, respectively.
Given a pair of matrices $\bm{A}$ and $\bm{B}$, define the inner product
$\langle\bm{A},\bm{B}\rangle:=\sum_{i,j=1}^{n}a_{i,j}b_{i,j}$. Let $\bm{I}$
and $\bm{1}$ denote the identity matrix and the vector of ones, respectively,
and let $\text{diag}(\bm{x})$ represent the diagonal matrix consisting of the
components of a vector $\bm{x}$.

We consider undirected weighted graphs with nonnegative edge weights. The
graph $\mathcal{G}=(\mathcal{V},\mathcal{E},\bm{W})$ is composed of a set of
nodes $\mathcal{V}$, edges $\mathcal{E}$, and a symmetric weight matrix
$\bm{W}\in\mathbb{R}^{n\times n}$ with $w_{i,j}>0$ if $(i,j)\in\mathcal{E}$,
and $w_{i,j}=0$ if $(i,j)\not\in\mathcal{E}$, where $n=|\mathcal{V}|$ is the
number of nodes. Here, $(i,i)\not\in\mathcal{E}$ for any $i$ by convention.
CGL is defined by $\bm{\Theta}=\bm{D}-\bm{W}\in\mathbb{R}^{n\times n}$, where
$\bm{D}:=\text{diag}(\bm{W1})$ is the degree matrix. CGL has zero row-sums
with its minimum eigenvalue also zero which is simple when the graph is
connected.

###  II-B Mathematical tools

The conjugate of a function $f(\bm{w})$ is denoted by
$f^{*}(\bm{y})=\sup_{\bm{w}\in\mathbb{R}^{N}}\langle\bm{w},\bm{y}\rangle-f(\bm{w})$,
$\bm{y}\in\mathbb{R}^{N}$. The set of proper lower semicontinuous convex
functions from $\mathbb{R}^{N}$ to ($-\infty,+\infty$] is denoted by
$\Gamma_{0}(\mathbb{R}^{N})$.222 A function $f$ is proper if
$\operatorname{dom}f:=\\{\bm{w}\in\left.\mathbb{R}^{N}\mid
f(\bm{w})<+\infty\right\\}\neq\emptyset$, and lower semicontinuous at $\bm{w}$
if $f(\bm{w})\leq\lim\inf_{\bm{y}\rightarrow\bm{w}}\\!f(\bm{y})$.  The
proximity operator of $f\in\Gamma_{0}(\mathbb{R}^{N})$ of index $\gamma>0$ is
defined as follows [34]:

$\operatorname{prox}_{\gamma
f}(\bm{w}):=\underset{\bm{y}\in\mathbb{R}^{N}}{\operatorname{argmin}}\left(f(\bm{y})+\frac{1}{2\gamma}\|\bm{w}-\bm{y}\|^{2}_{2}\right).$
(1)

Uniqueness and existence of the minimizer is guaranteed by the strong
convexity and coercivity of $f+\frac{1}{2\gamma}\|\bm{w}-\cdot\|^{2}_{2}$. The
indicator function with respect to a given set $\mathcal{S}$ is denoted by

$\iota_{\mathcal{S}}(\bm{w}):=\begin{cases}0,&\mbox{if
}\bm{w}\in\mathcal{S},\\\ +\infty,&\mbox{otherwise}.\end{cases}$ (2)

It is clear by definition that
$\mathrm{prox}_{\iota_{C}}(\bm{w})=P_{C}(\bm{w}):=\underset{\bm{y}\in
C}{\operatorname{argmin}}\|\bm{w}-\bm{y}\|_{2}$. The Moreau envelope of a
function $f\in\Gamma_{0}(\mathbb{R}^{N})$ of index $\gamma>0$ is defined as
follows [34, Definition 12.20]:

${}^{\gamma}\\!f(\bm{w})$ $\displaystyle:=$
$\displaystyle\min_{\bm{y}\in\mathbb{R}^{N}}\left(f(\bm{y})+\frac{1}{2\gamma}\|\bm{w}-\bm{y}\|_{2}^{2}\right).$
(3)

Using the Moreau envelope ${}^{\gamma}\|\cdot\|_{1}$ of $\|\cdot\|_{1}$, which
is the widely known Huber function, the MC penalty [24] is defined as

$\displaystyle\phi_{\mathrm{MC}}(\bm{w})$ $\displaystyle=$
$\displaystyle\|\bm{w}\|_{1}-^{\gamma}\\!\|\cdot\|_{1}(\bm{w}).$ (4)

The nonconvex function $\phi_{\mathrm{MC}}$ here is known to induce a sparser
and less biased estimate with respect to the $\ell_{1}$ penalty.

###  II-C Primal-dual splitting method

Let $\mathcal{X}$ and $\mathcal{Y}$ be real Hilbert spaces: in the present
case, $\mathcal{X}:=\mathbb{R}^{N}$ and $\mathcal{Y}:=\mathbb{R}^{n\times n}$.
The primal-dual splitting method [32] solves convex optimization problems in
the following form:

$\min_{\bm{w}\in\mathcal{X}}[F(\bm{w})+G(\bm{w})+H(L(\bm{w}))],$ (5)

where $F\\!:\\!\mathcal{X}\rightarrow\mathbb{R}$ is a differentiable convex
function with Lipschitz continuous gradient $\nabla F$,
$G\in\Gamma_{0}(\mathcal{X})\\!$ and $H\in\Gamma_{0}(\mathcal{Y})$ are
proximable proper lower semicontinuous convex functions, and
$L:\mathcal{X}\rightarrow\mathcal{Y}$ is a bounded linear operator with its
adjoint operator denoted by $L^{*}$. Here, “proximable” means that the
proximity operator of the function can be computed easily (in a closed form in
the present case). The primal dual splitting method is given in Algorithm 1.

Algorithm 1 Primal-dual splitting method

0: Initial estimate $\bm{w}_{0}\in\mathcal{X},\bm{V}_{0}\in\mathcal{Y}$,
tolerance $\epsilon>0$, proximity parameters $\tau>0$ and $\sigma>0$,
relaxation parameters $\rho_{k}>0$.

while
$\frac{\|\bm{w}_{k+1}-\bm{w}_{k}\|^{2}_{2}}{\|\bm{w}_{k}\|^{2}_{2}}>\epsilon$
do

1\. $\tilde{\bm{w}}_{k+1}:=\operatorname{prox}_{\tau
G}\left(\bm{w}_{k}-\tau\nabla F\left(\bm{w}_{k}\right)-\tau
L^{*}\bm{V}_{k}\right)$

2\. $\tilde{\bm{V}}_{k+1}:=\operatorname{prox}_{\sigma
H^{*}}\left(\bm{V}_{k}+\sigma
L\left(2\tilde{\bm{w}}_{k+1}-\bm{w}_{k}\right)\right)$

3\.
$\left(\bm{w}_{k+1},\bm{V}_{k+1}\right)\\!:=\\!\rho_{k}\\!\left(\\!\tilde{\bm{w}}_{k+1},\tilde{\bm{V}}_{k+1}\\!\right)\\!+\\!\left(1\\!-\\!\rho_{k}\right)\\!\left(\bm{w}_{k},\bm{V}_{k}\right)$

end while

##  III PROPOSED ALGORITHMS

Due to its structure (i.e., symmetry and zero row-sums), the CGL is completely
defined by its upper (or lower) triangular part excluding the main diagonal,
or, in other words, by a length-$\frac{n(n-1)}{2}$ vector, where the CGL is of
size $n\times n$. Since all the off-diagonal components of CGL need to be
nonnegative, our variable vector is constrained to the nonnegative cone (the
nonnegative orthant) $C:=\mathbb{R}_{+}^{\frac{n(n-1)}{2}}$. Given this, we
define a specific linear operator $L:\\!C\\!\rightarrow\\!\mathbb{R}^{n\times
n}$ that maps a nonnegative vector of size $n(n-1)/2$ to its corresponding
CGL. For $n=4$, for instance, $L$ is defined as follows:

$\displaystyle L:[w_{1},w_{2},w_{3},w_{4},w_{5},w_{6}]^{T}\mapsto$
$\displaystyle\left[\begin{array}[]{cccc}\\!w_{1}\\!+\\!w_{2}\\!+\\!w_{3}\\!&-w_{1}&-w_{2}&-w_{3}\\\
-w_{1}&\\!w_{1}\\!+\\!w_{4}\\!+\\!w_{5}\\!&-w_{4}&-w_{5}\\\
-w_{2}&-w_{4}&\\!w_{2}\\!+\\!w_{4}\\!+\\!w_{6}\\!&-w_{6}\\\
-w_{3}&-w_{5}&-w_{6}&\\!w_{3}\\!+\\!w_{5}\\!+\\!w_{6}\\!\\\
\end{array}\right].$

###  III-A Problem formulation

The CGL formulation presented in [10] is a popular extension of graphical
lasso for imposing a Laplacian constraint. With a slight modification using
the linear operator $L$ introduced above, the CGL formulation is given by

$\displaystyle\mathrm{P}_{0}:\min_{\bm{w}\in
C}\;-\mathrm{logdet}(L(\bm{w})\\!+\\!\bm{J})\\!+\\!\langle\bm{S},L(\bm{w})\rangle\\!+\\!\lambda_{1}\|\bm{w}\|_{1},$
(6)

where $\bm{J}:=\frac{1}{n}\bm{11}^{T}\in\mathbb{R}^{n\times n}$,
$\bm{S}\in\mathbb{R}^{n\times n}$ stands for the sample covariance obtained
from data, and $\lambda_{1}\geq 0$ is the regularization parameter. Note here
that $L(\bm{w})+\bm{J}$ is positive definite if and only if the graph of
$L(\bm{w})$ is connected, i.e., $L(\bm{w})\circ\bm{I}-L(\bm{w})$ is an
irreducible matrix, where $\circ$ denotes the Hadamard product.

Replacing the regularization term of Problem $P_{0}$ by the MC penalty given
in (4), the problem reads as follows:333Although the formulation in P1 has
been considered in the literature [30, 31], it was also considered earlier in
the authors’ previous works [35, 33] as an intermediate step, and the present
study is independent from [30, 31].

$\displaystyle\mathrm{P}_{1}:\min_{\bm{w}\in
C}\;-\mathrm{logdet}(L(\bm{w})+\bm{J})+\langle\bm{S},L(\bm{w})\rangle$ (7)
$\displaystyle+\lambda_{1}\underbrace{\left[\|\bm{w}\|_{1}-^{\gamma}\\!\|\cdot\|_{1}(\bm{w})\right].}_{\mathrm{MC}}$

We introduce the Tikhonov regularization term
$\frac{\lambda_{2}}{2}\|\bm{w}\|_{2}^{2}$, $\lambda_{2}\geq 0$, which also
plays a role of convexification as shown below. Introducing the indicator
function $\iota_{C}(\bm{w})$ to accommodate the constraint as well, Problem
$\mathrm{P}_{1}$ is transformed into the following unconstrained optimization
problem:

$\displaystyle\mathrm{P_{2}}:$
$\displaystyle\min_{\bm{w}\in\mathbb{R}^{\frac{n(n-1)}{2}}}\;\underbrace{-\lambda_{1}\
^{\gamma}\|\cdot\|_{1}(\bm{w})+\frac{\lambda_{2}}{2}\|\bm{w}\|_{2}^{2}}_{F(\bm{w})}$
$\displaystyle+\underbrace{\iota_{C}(\bm{w})\\!+\\!\lambda_{1}\|\bm{w}\|_{1}\\!+\\!\langle\bm{S},L(\bm{w})\rangle}_{G(\bm{w})}\underbrace{-\mathrm{logdet}(L(\bm{w})\\!+\\!\bm{J})}_{H(L(\bm{w}))}.$
(8)

By using Moreau’s decomposition $\frac{1}{2\gamma}\|\cdot\|_{2}^{2}=\hskip
1.00006pt^{\gamma}f+\hskip
1.00006pt^{\gamma^{-1}}f^{*}\circ\gamma^{-1}\,\bm{I}$ [34, Theorem 14.3], the
MC penalty term can be rewritten as

$\displaystyle-^{\gamma}\|\cdot\|_{1}(\bm{w})$ $\displaystyle\\!=\\!$
${}^{\gamma^{-1}}(\|\cdot\|_{1}^{*})(\gamma^{-1}\bm{w})\\!-\\!\frac{1}{2\gamma}\|\bm{w}\|_{2}^{2}$
(9) $\displaystyle\\!=\\!$
${}^{\gamma}\iota_{B_{\infty}}(\bm{w})-\frac{1}{2\gamma}\|\bm{w}\|_{2}^{2},$
(10)

where
$B_{\infty}:=\mathrm{lev}_{\leq\gamma}\|\cdot\|_{\infty}:=\\{\bm{x}\in\mathbb{R}^{n}\mid\|\bm{x}\|_{\infty}\leq\gamma\\}$
is the $\ell_{\infty}$ ball of radius-$\gamma$ [34, Example 13.32]. The
usefulness of this decomposition of the MC penalty term has been observed also
in [36, 37, 38, 39]. Using (10), the function $F$ of P2 can be rewritten as

$F(\bm{w})=\lambda_{1}\
^{\gamma}\iota_{B_{\infty}}(\bm{w})-\frac{\lambda_{1}}{2\gamma}\|\bm{w}\|_{2}^{2}+\frac{\lambda_{2}}{2}\|\bm{w}\|_{2}^{2},$
(11)

of which the convexity is ensured clearly by choosing $\lambda_{1}$ and
$\lambda_{2}$ such that $\lambda_{2}\geq\gamma^{-1}\lambda_{1}$ (see
Proposition 4 below). The Tikhonov regularization term
$\frac{\lambda_{2}}{2}\|\bm{w}\|_{2}^{2}$ thus has a convexification property,
as mentioned above. On the other hand, the functions $G(\bm{w})$ and $H\circ
L(\bm{w})$ are convex, since the composition of a convex function with an
arbitrary affine operator is also a convex function. Hence, under the
convexity condition given above, Problem $\mathrm{P_{2}}$ takes the form of
(5), and it can be solved by the primal-dual splitting method.

###  III-B Optimization algorithm

The proposed algorithm is derived by applying the primal-dual spitting method
to Problem $\mathrm{P_{2}}$.

####  III-B1 Derivation of $\tilde{\bm{w}}_{k+1}$

We define the soft thresholding operator for a length-$\frac{n(n-1)}{2}$
positive vector
$\bm{\delta}:=[\delta_{1},\delta_{2},\cdots,\delta_{\frac{n(n-1)}{2}}]^{T}$ by

$\displaystyle[\mathrm{soft}_{\bm{\delta}}(\bm{w})]_{i}=\left\\{\begin{array}[]{lll}w_{i}-\delta_{i},&\mathrm{if}~{}w_{i}\geq\delta_{i},\\\
0,&\mathrm{if}~{}|w_{i}|<\delta_{i},\\\
w_{i}+\delta_{i},&\mathrm{if}~{}w_{i}\leq-\delta_{i},\\\ \end{array}\right.$
(15)

where $[\cdot]_{i}$ is the $i$th component of the argument. The convex
projection onto the nonnegative cone $C$ is given by

$\displaystyle[P_{C}(\bm{w})]_{i}=\left\\{\begin{array}[]{lll}w_{i},&\mathrm{if}~{}w_{i}\geq
0,\\\ 0,&\mathrm{if}~{}w_{i}<0.\\\ \end{array}\right.$ (18)

Applying Step 1 of Algorithm 1 to Problem $\mathrm{P_{2}}$ yields

$\displaystyle\tilde{\bm{w}}_{k+1}$ $\displaystyle=\operatorname{prox}_{\tau
G}\bigl{[}\bm{w}_{k}-\tau L^{*}(\bm{V}_{k})-\tau\nabla
F\left(\bm{\bm{w}}_{k}\right)\bigr{]}$ (19)
$\displaystyle=\operatorname{prox}_{\tau G}\Bigl{[}\bm{w}_{k}-\tau
L^{*}(\bm{V}_{k})-\tau\bigl{(}\gamma^{-1}\lambda_{1}\mathrm{prox}_{\|\cdot\|_{1}}(\bm{w}_{k})$
$\displaystyle\quad-\gamma^{-1}\lambda_{1}\bm{w}_{k}+\lambda_{2}\bm{w}_{k}\bigr{)}\Bigr{]}.\
$ (20)

The operators $L^{*}$ and $\mathrm{prox}_{\tau G}$ can be computed by using
the following propositions.

###### Proposition 1

Let $\bm{M}\in\mathbb{R}^{n\times n}$ be an arbitrary CGL matrix with its
$(p,q)$ component denoted by $m_{p,q}$. Then, for any
$p,q\in\\{1,2,\cdots,n\\}$ such that
$(2n-p-1)p/2+q-n\in\\{1,2,\cdots,n(n-1)/2\\}$, it holds that

$[L^{*}(\bm{M})]_{(2n-p-1)p/2+q-n}=m_{p,p}+m_{q,q}-m_{p,q}-m_{q,p}.$

Proof: See Appendix A.

###### Proposition 2

The proximity operator of
$G(\bm{w})\\!:=\\!\iota_{c}(\bm{w})\\!+\\!\lambda_{1}\|\bm{w}\|_{1}\\!+\\!\langle\bm{S},L(\bm{w})\rangle$
of index $\tau>0$ can be expressed by

$\displaystyle\mathrm{prox}_{\tau
G}(\bm{w})=P_{C}(\bm{w}-\tau(\lambda_{1}\bm{1}+L^{*}(\bm{S}))).$ (21)

Proof: See Appendix B.

By using Proposition 2 and
$\mathrm{prox}_{\|\cdot\|_{1}}=\mathrm{soft}_{\bm{1}}$, (20) can be rewritten
as

$\displaystyle\tilde{\bm{w}}_{k+1}=P_{C}[\bm{w}_{k}-\tau
L^{*}(\bm{V}_{k})-\tau(\lambda_{1}\bm{1}+L^{*}(\bm{S}))$
$\displaystyle-\tau(\gamma^{-1}\lambda_{1}\mathrm{soft}_{\bm{1}}(\bm{w}_{k})-\gamma^{-1}\lambda_{1}\bm{w}_{k}+\lambda_{2}\bm{w}_{k})].$
(22)

####  III-B2 Derivation of $\tilde{\bm{V}}_{k+1}$

Substituting Moreau’s decomposition [34, Theorem 14.3]

$\displaystyle\operatorname{prox}_{\sigma H^{*}}(\bm{u})$ $\displaystyle=$
$\displaystyle\bm{u}-\sigma\operatorname{prox}_{\sigma^{-1}H}(\sigma^{-1}\bm{u})$
(23)

with $H:=-\mathrm{logdet}(\cdot+\bm{J})$ into Step 2 of Algorithm 1 yields

$\displaystyle\tilde{\bm{V}}_{k+1}$ $\displaystyle\\!=\\!\bm{V}_{k}+\sigma
L\left(2\tilde{\bm{w}}_{k+1}-\bm{\bm{w}}_{k}\right)$ $\displaystyle-$
$\displaystyle\sigma\operatorname{prox}_{\sigma^{-1}\left(-\mathrm{logdet}(\cdot+\bm{J})\right)}\left[\sigma^{-1}\bm{\bm{V}}_{k}\\!+\\!L\left(2\tilde{\bm{w}}_{k+1}\\!-\\!\bm{\bm{w}}_{k}\right)\right].$

Here, the proximity operator
$\mathrm{prox}_{\sigma^{-1}(-\mathrm{logdet}(\hskip 1.39998pt\cdot\hskip
1.39998pt+\hskip 1.39998pt\bm{J}))}$ can be written in a closed form, as shown
in the following proposition.

###### Proposition 3

For a positive semi-definite matrix $\bm{W}\in\mathbb{R}^{n\times n}$, it
holds that

$\displaystyle\mathrm{prox}_{\sigma^{-1}(-\mathrm{logdet}(\hskip
1.39998pt\cdot\hskip 1.39998pt+\hskip 1.39998pt\bm{J}))}(\bm{W})$
$\displaystyle=\bm{Q}\,\text{diag}\left(\\!\frac{\mu_{1}\\!+\\!\sqrt{\mu_{1}^{2}\\!+\\!4\sigma^{-1}}}{2},\\!\cdots\\!,\frac{\mu_{n}\\!+\\!\sqrt{\mu_{n}^{2}\\!+\\!4\sigma^{-1}}}{2}\right)\\!\bm{Q}^{T}$
$\displaystyle\quad-\bm{J},$ (25)

where $\mu_{i}$ is the $i$th eigenvalue of $\bm{W}+\bm{J}$, and
$\bm{Q}:=[\bm{q}_{1}\bm{q}_{2}\ldots\bm{q}_{n}]$ with the eigenvectors
$\bm{q}_{i}$ of $\bm{W}+\bm{J}$.
Proof: See Appendix C.

An application of Proposition 3 to (LABEL:y_np1_op) yields

$\displaystyle\tilde{\bm{V}}_{k+1}\\!=\\!\bm{V}_{k}+\sigma
L\left(2\tilde{\bm{w}}_{k+1}-\bm{w}_{k}\right)+\sigma\bm{J}$
$\displaystyle-\\!\sigma\\!\left[\\!\bm{U}\text{diag}\\!\left(\\!\dfrac{\nu_{1}\\!+\\!\sqrt{\nu_{1}^{2}\\!+\\!4\sigma^{-\\!1}}}{2}\\!,\\!\cdots\\!,\\!\dfrac{\nu_{n}\\!+\\!\sqrt{\nu_{n}^{2}\\!+\\!4\sigma^{-\\!1}}}{2}\\!\right)\\!\bm{U}^{T}\\!\right],$
(26)

where $\nu_{i}$ is the eigenvalue of
$\sigma^{-1}\bm{\bm{V}}_{k}+L\left(2\tilde{\bm{w}}_{k+1}-\bm{\bm{w}}_{k}\right)$,
and $\bm{U}:=[\bm{u}_{1}\bm{u}_{2}\ldots\bm{u}_{n}]$ with the eigenvectors
$\bm{u}_{i}$ of
$\sigma^{-1}\bm{\bm{V}}_{k}+L\left(2\tilde{\bm{w}}_{k+1}-\bm{\bm{w}}_{k}\right)$.
The proposed algorithm is given in Algorithm 2. Our formulation based on the
MC penalty is expected to yield a sparser solution with better
interpretability than the conventional $\ell_{1}$-based methods due to the
efficient sparsity promoting property of the MC penalty. In addition, our
representation of CGL using the linear operator reduces the number of
variables approximately by half, while also transforming the positive semi-
definite constraint of the graph Laplacian to the nonnegativity constraint
$\bm{w}\in C$. However, our formulation needs $\mathcal{O}(n^{3})$ complexity
due to the need to execute matrix multiplication and eigenvalue decomposition.
This computational drawback can be mitigated by using the eigenvalue
decomposition method for symmetric matrices [40]. Note that, in the particular
case of $\gamma:=+\infty$, the proposed algorithm gives an alternative way to
solve the graphical lasso problem for CGL.

Algorithm 2 Proposed graph learning algorithm

0: Initial estimate ($\bm{w}_{0},\bm{V}_{0}$), tolerance $\epsilon>0$,
proximity parameters $\tau>0$, $\sigma>0$, covariance matrix $\bm{S}$,
regularization parameter $\lambda_{1}\geq 0,\lambda_{2}\geq 0$, minimax
concave parameter $\gamma^{-1}>0$, relaxation parameters $\rho_{k}>0.$

0: Graph Laplacian $\bm{\Theta}$

while
$\frac{\|\bm{w}_{k+1}-\bm{w}_{k}\|^{2}_{2}}{\|\bm{w}_{k}\|^{2}_{2}}>\epsilon$
do

1\. Compute the vector $\tilde{\bm{w}}_{k+1}$ by (22)

2\. Find the eigenvalues $\nu_{i}$ and the matrix
$\bm{U}=[\bm{u}_{1}\ldots\bm{u}_{n}]$ containing all the corresponding (unit-
norm) eigenvectors of
$\left(\bm{J}+\sigma^{-1}\bm{V}_{k}+L\left(2\tilde{\bm{w}}_{k+1}-\bm{w}_{k}\right)\right)$

3\. Compute the vector $\tilde{\bm{V}}_{k+1}$ by (26)

4\.
$\left(\bm{w}_{k+1},\bm{V}_{k+1}\right)\\!=\\!\rho_{k}\\!\left(\tilde{\bm{w}}_{k+1},\tilde{\bm{V}}_{k+1}\right)\\!+\\!\left(1\\!-\\!\rho_{k}\right)\\!\left(\bm{w}_{k},\bm{V}_{k}\right)$

end while

return $\bm{\Theta}=L(\bm{w}_{k})$

###  III-C Convergence conditions

Convergence is guaranteed under the following conditions.

###### Proposition 4

If $\lambda_{2}\geq\gamma^{-1}\lambda_{1}$, the function $F$ is convex. In
this case, Algorithm 1 converges to a global minimizer of $\mathrm{P}_{2}$ if
the following conditions are jointly satisfied:

  1. 1.

$\dfrac{1}{\tau}\geq 2\sigma n+\dfrac{\lambda_{2}}{2},$

  2. 2.

$0<\rho_{k}<2-\dfrac{\lambda_{2}}{2}\left(\dfrac{1}{\tau}-2\sigma
n\right)^{-1}.$

Proof:  See Appendix D.

From a theoretical side, a use of $\lambda_{2}$ satisfying the convexity
condition shown in Proposition 4 ensures convergence to a global minimizer.
From a practical side, on the other hand, a use of $\lambda_{2}$ violating the
convexity condition may yield better performance, as will be seen in Section
IV. However, we emphasize that improved performance comes with no theoretical
guarantees. A remarkable advantage of the present framework is its flexibility
due to the use of powerful convex analytic solver, which allows to extend the
presented framework in many possible directions including dynamic graph
learning [41, 42].

##  IV Numerical Examples

We show the efficacy of the proposed method through some experiments with
synthetic and real data. We first show the performances of the proposed method
for different regularization parameters. We then compare the performance of
the proposed method with CGL estimation [10] and SGL estimation [11]. 444We
used the implementations of the SGL and CGL algorithms in
spectralGraphTopology
(https://CRAN.R-project.org/package=spectralGraphTopology).

###  IV-A Experiments with synthetic data

Dataset generation: We consider three types of graph: (i) grid graph
$\mathcal{G}_{\mathrm{grid}}^{(\sqrt{n},\sqrt{n})}$ with nodes connected to
their four nearest neighbors (except the nodes at boundaries), (ii) random
modular graph (a.k.a. stochastic block model)
$\mathcal{G}_{\mathrm{M}}^{(n,0.01,0.3)}$ with four modules where the nodes
are connected across the modules and within each module with probabilities
0.01 and 0.3, respectively, and (iii) Erdös-Rényi graph
$\mathcal{G}_{\mathrm{ER}}^{(n,0,1)}$ with nodes connected to other nodes with
probability 0.1. The graph weights are randomly drawn from the uniform
distribution over the interval [0.1, 3.0], regarded as the ground-truth graph
Laplacian $\bm{\Theta}$ in this experiment. From each graph generated, data
are generated from $\mathcal{N}(0,\bm{\Theta}^{\dagger})$, where
$(\cdot)^{\dagger}$ denotes the Moore-Penrose pseudo inverse, and the
covariance matrix $\bm{S}$ is computed from data. For each type of graph, we
randomly generate $15$ graphs with $n=100$ nodes using the toolbox given in
[43].

Performance measure: The relative error (RE) and F-score (FS) are used as
performance measures:

$\displaystyle\operatorname{RE}\left(\widehat{\bm{\Theta}},\bm{\Theta}_{\star}\right):=\frac{\left\|\widehat{\bm{\Theta}}-\bm{\Theta}_{\star}\right\|^{2}_{\rm
F}}{\left\|\bm{\Theta}_{\star}\right\|^{2}_{\rm F}},$ (27)
$\displaystyle\mathrm{FS}\left(\widehat{\bm{\Theta}},\bm{\Theta}_{\star}\right):=\frac{2\mathrm{tp}}{2\mathrm{tp}+\mathrm{f}\mathrm{n}+\mathrm{fp}},$
(28)

where $\mathrm{tp,fp}$, and $\mathrm{fn}$ stand for true-positive, false-
positive, and false-negative, respectively. Here, RE indicates the discrepancy
between the ground-truth graph Laplacian
$\bm{\Theta}_{\star}\in\mathbb{R}^{n\times n}$ and its estimate
$\widehat{\bm{\Theta}}\in\mathbb{R}^{n\times n}$, while FS is a measure of
accuracy of binary classification (taking values in [0,1]), indicating whether
the sparse structures are extracted correctly.

####  IV-A1 Performance of the proposed method

We study the impacts of the parameters $\lambda_{1}$ and $\gamma^{-1}$ of the
MC penalty on the performance of the proposed method. We also tested the case
of $\lambda_{2}:=0$, which makes the entire cost function nonconvex for any
$\lambda_{1}>0$, with the other parameters tuned manually (see Table I). To
study the impact of $\lambda_{1}$, we fix $\gamma^{-1}:=2.25$ which gave a
best performance in the nonconvex case. Since $\tau$ is the algorithm
parameter and it only affects the convergence speed, we fix it to $\tau:=1.0$.
The other parameters are then set to $\lambda_{2}:=\gamma^{-1}\lambda_{1}$ and
$\sigma\approx(1/\tau-\lambda_{2}/2)/(2n)$ according to the convexity
condition (See Proposition 4). Figure 1 plots the RE and FS curves across
$m/n$ in modular graph $\mathcal{G}_{\mathrm{M}}^{(100,0.01,0.3)}$ for
different choices of $\lambda_{1}$. To study the impact of $\gamma^{-1}$, on
the other hand, we fix $\lambda_{1}:=1.0\times 10^{-4}$ and choose the other
parameters in the same way as in Fig. 1. Figure 2 plots the RE and FS curves
for different choices of $\gamma^{-1}$ under the same conditions as in Fig. 1,
with $\lambda_{1}$=$1.0\times 10^{-4}$. In Figs. 1 and 2, the proposed method
using the convexity condition attains better performance than the nonconvex
case when $m/n$ is small, while the nonconvex case is better when $m/n$ is
large. To be more specific, when $m/n$ is small, using larger $\lambda_{1}$,
or larger $\gamma^{-1}$, yields better performance. Although the
regularization parameter $\lambda_{2}$ for the Tikhonov regularization needs
to be sufficiently large to ensure the convexity of the entire objective,
using a $\lambda_{2}$ that is too large tends to yield a less sparse solution,
which means degradation of graph interpretability (cf. Section III-C).

TABLE I: Parameter settings for each graph. parameter | grid | modular | ER
---|---|---|---
| $\lambda_{1}$ | 0.005 | 0.01 | 0.01
nonconvex | $\lambda_{2}$ | 0 | 0 | 0
| $\tau$ | 1.0 | 1.0 | 1.0
| $\sigma$ | 0.05 | 0.05 | 0.01
| $\lambda_{1}$ | $1.0\\!\times\\!10^{-4}$ | $1.0\\!\times\\!10^{-4}$ | $1.0\\!\times\\!10^{-4}$
convex | $\lambda_{2}$ | $2.5\\!\times\\!10^{-4}$ | $2.5\\!\times\\!10^{-4}$ | $2.5\\!\times\\!10^{-4}$
| $\tau$ | 1.0 | 1.0 | 1.0
| $\sigma$ | $4.9\\!\times\\!10^{-3}$ | $4.9\\!\times\\!10^{-3}$ | $4.9\\!\times\\!10^{-3}$
$\gamma^{-1}$ | 2.25
$\rho_{k}$ | 1.0
Maximum iterations | 5000
Tolerance error | $1.0\times 10^{-4}$

(a) RE (b) FS

Figure 1: Experimental results of each $\lambda_{1}$ value under the convexity
condition in random modular graph $\mathcal{G}_{\mathrm{M}}^{(100,0.01,0.3)}$
estimation.

(a) RE (b) FS

Figure 2: Experimental results of each $\gamma^{-1}$ value under the convexity
condition in random modular graph $\mathcal{G}_{\mathrm{M}}^{(100,0.01,0.3)}$
estimation.

####  IV-A2 Comparisons in estimation accuracy

Parameters: The best parameters are chosen manually, see Table I. The
generated graphs depend only on $\lambda_{1}$, $\lambda_{2}$, and
$\gamma^{-1}$, while being independent of the algorithm parameters $\tau$,
$\sigma$, and $\rho_{k}$ in principle. It is very important to tune the
parameters $\lambda_{1}$, $\lambda_{2}$, and $\gamma^{-1}$ carefully for
better performance, although the generated graph changes gradually as each of
those parameters changes. As shown in Section IV-A1, for the proposed method
under the convexity condition, $\lambda_{1}=1.0\times 10^{-4}$ is used. We
mention that a smaller value of $\lambda_{1}$ tends to give a better result
for RE and FS when $m/n$ is large. Although the convergence is guaranteed
under the condition of $\lambda_{2}$ described in Section III-C,
$\lambda_{2}:=0$ (for which there is no guarantee of convergence to the global
minimizer) gave the best performance in the current experiments. Regarding the
parameters for CGL and SGL, we follow the parameter selection techniques
proposed in [10] and [11], respectively.

Results: Figure 3 shows the ground truth and the learned graph for $m/n=100$
where $m$ is the number of measurements. One can see that the proposed method
yields a more accurate graph than the other methods; in particular, the graph
obtained by the proposed method is remarkably sparse. Figures 4 – 6 show the
performances in RE and FS across $m/n$ for the grid graph
$\mathcal{G}_{\mathrm{grid}}^{(10,10)}$, the random modular graph
$\mathcal{G}_{\mathrm{M}}^{(100,0.01,0.3)}$, and the Erdös-Rényi graph
$\mathcal{G}_{\mathrm{ER}}^{(100,0.1)}$, respectively. In Fig. 4, the proposed
method (nonconvex) significantly outperforms CGL due to the use of the MC
penalty, while the proposed method under the convexity condition achieves
approximately the same RE performance as the nonconvex case with the degraded
FS performance for large $m/n$. The performance of SGL is close to that of the
proposed method in this case. The proposed method under the convexity
condition exhibits approximately the same performance as the proposed method
for $\lambda_{2}=0$ (the nonconvex case). In Fig. 5, on the other hand, the
proposed algorithm outperforms SGL considerably. The proposed method
significantly outperforms CGL in FS for large $m/n$, although those two
methods exhibit comparable performances in RE (and in FS as well for small
$m/n$). In Fig. 6, the proposed method outperforms the other methods in FS for
a wide range of $m/n$ values, while the proposed method under the convexity
condition achieves approximately the same RE performance as the nonconvex case
with its FS performance close to those of CGL and SGL for large $m/n$. We
remark that the difference between the proposed method and CGL is more notable
in Fig. 4 than in Figs. 5 and 6 because the graph in Fig. 4 is approximately
four times sparser than that of Fig. 5 and twice sparser than that of Fig. 6.
Thus, the sparsity assumption is a better match to the actual data in for the
graph from Fig. 4. We finally remark that, for CGL and SGL, a small
regularization parameter was used because use of large regularization
parameters with the $\ell_{1}$ norm leads to increased errors, which degrade
the quality (e.g., interpretability) of the learned graphs. This is the reason
why the graphs obtained by CGL and SGL in Fig. 3 are not sufficiently sparse.

(a) Ground truth (b) Proposed
---

(c) CGL (d) SGL
---

Figure 3: Visualization of estimate graphs, in which the edge colors represent
the edge weights.

(a) RE (b) FS

Figure 4: Experimental results for grid graph
$\mathcal{G}_{\mathrm{grid}.}^{(10,10)}$

(a) RE (b) FS

Figure 5: Experimental results for random modular graph
$\mathcal{G}_{\mathrm{M}.}^{(100,0.01,0.3)}$

(a) RE (b) FS

Figure 6: Experimental results for Erdös-Rényi graph
$\mathcal{G}_{\mathrm{ER}.}^{(100,0.1)}$

####  IV-A3 Comparisons in computation time

We investigate how the computation time with tolerance error $1.0\times
10^{-4}$ changes with the size of the graph for the modular graph
$\mathcal{G}_{\mathrm{modular}}^{(n,0.01,0.3)}.$ We set $m/n=5000$ and perform
graph learning 15 times as in Section IV-A. The computation time for
$n=160,240,320,400$ is summarized in Table II. It can be seen that the
proposed method is 5.38 – 43.4 times faster than SGL. Although the computation
time of the proposed method is higher due mainly to the eigenvalue
decomposition of the larger sized matrix compared to CGL, we emphasize that
the performance improvements are remarkable especially for the grid graph (See
Fig. 4). The significant advantage in CPU time is due to the fact that the
proposed method requires a few thousand iterations for approximate
convergence, while SGL requires over $10^{5}$ iterations on average with per-
iteration complexity of order $\mathcal{O}(n^{3})$.

TABLE II: A comparison of average CPU time (in seconds) for a modular graph $\mathcal{G}_{\mathrm{modular}}^{(n,0.01,0.3)}$ with $m/n=5000$. $n$ | proposed method | CGL | SGL
---|---|---|---
160 | 111.2 | 14.95 | 6958
240 | 427.8 | 51.99 | 40565
320 | 1601 | 152.6 | 69555
400 | 4484 | 381.3 | 145613

###  IV-B Experiments with real data

We test our method for the animal dataset [44], in which each node represents
each animal and the edges represent how much the animals are related to each
other. The dataset contains binary values (i.e., it is a categorical dataset)
which represent the answers to some questions such as “has lungs?” for
instance. There are 102 such questions in total, answered for 33 different
animals. The covariance matrix is created based on this data set, and the
graph is learned in the same way as in the previous experiment. The results
are shown in Fig. 7. It can be seen that the proposed method produces a
sparser graph than CGL and SGL while preserving the dominant links.

(a) Proposed (b) CGL (c) SGL

Figure 7: The relations of animals from animal dataset.

##  V Conclusion

We presented a graph learning method inserting the nonconvex MC penalty into
the extension of the graphical lasso formulation to produce sparse and
accurate graphs. With the linear-operator-based representation of CGL together
with the reformulation through the Moreau decomposition, an efficient
algorithm was derived with the primal-dual splitting method, for which an
admissible choice of parameters was presented to ensure the provable
convergence. Numerical examples showed that the proposed method significantly
outperformed CGL for the high-sparsity grid graphs, while still achieving
higher F-scores than CGL and SGL (the state-of-the-art method) for the low-
sparsity modular graphs. In addition, the comparisons in CPU time showed that
the proposed method was dramatically faster than SGL.

##  Appendix A Proof of Proposition 1

Let
$\bm{w}:=[w_{1,2},w_{1,3},\\!\cdots\\!,w_{1,n},w_{2,3},\cdots,w_{2,n},\\!\cdots\\!,w_{n-1,n}]^{T}$
$\in C$. By definition of adjoint operator, it holds that $\langle
L(\bm{w}),\bm{M}\rangle$=$\langle L^{*}(\bm{M}),\bm{w}\rangle$, of which the
left side can be expanded as
$w_{1,2}(m_{1,1}+m_{2,2}-m_{1,2}-m_{2,1})+w_{1,3}(m_{1,1}+m_{3,3}-m_{1,3}-m_{3,1})+\cdots$.
It can therefore be seen that
$[L^{*}(\bm{M})]_{1}=m_{1,1}+m_{2,2}-m_{1,2}-m_{2,1}$ and
$[L^{*}(\bm{M})]_{2}=m_{1,1}+m_{3,3}-m_{1,3}-m_{3,1}$. In general, it can be
verified that
$[L^{*}(\bm{M})]_{(2n-p-1)p/2+q-n}=m_{p,p}+m_{q,q}-m_{p,q}-m_{q,p}$.

##  Appendix B Proof of Proposition 2

The assertion can be verified by combining the basic property [34, Proposition
24.8]

$\mathrm{prox}_{\phi+\langle\cdot,\bm{u}\rangle}(\bm{x})=\mathrm{prox}_{\phi}(\bm{x}-\bm{u})$
(B.1 )

and the following fact [45, Proposition 1]:

$\operatorname{prox}_{\iota_{C}+\eta\|\cdot\|_{1}}(\bm{w})=P_{C}(\bm{w}-\eta\bm{1})$
(B.2 )

for $\eta>0$. Indeed, since
$G:=\iota_{C}+\lambda_{1}\|\cdot\|_{1}+\langle\bm{S},L(\cdot)\rangle=\iota_{C}+\lambda_{1}\|\cdot\|_{1}+\langle\cdot,L^{*}(\bm{S})\rangle$,
we have $\tau G=\iota_{C}+\tau\lambda_{1}\|\cdot\|_{1}+\langle\cdot,\tau
L^{*}(\bm{S})\rangle$. Applying (B.1) to
$\phi:=\iota_{C}+\tau\lambda_{1}\|\cdot\|_{1}$ and $u=\tau L^{*}(\bm{S})$, we
obtain

$\displaystyle\mathrm{prox}_{\tau G}(\bm{w})$ $\displaystyle=$
$\displaystyle\mathrm{prox}_{\iota_{C}+\tau\lambda_{1}\|\cdot\|_{1}+\langle\cdot,\tau
L^{*}(\bm{S})\rangle}(\bm{w})$ (B.3 ) $\displaystyle=$
$\displaystyle\mathrm{prox}_{\iota_{C}+\tau\lambda_{1}\|\cdot\|_{1}}(\bm{w}-\tau
L^{*}(\bm{S}))$ $\displaystyle=$ $\displaystyle P_{C}(\bm{w}-\tau
L^{*}(\bm{S})-\tau\lambda_{1}\bm{1}).$

##  Appendix C Proof of Proposition 3

The proximity operator of

$f(\bm{W})\\!\mapsto\\!\left\\{\begin{array}[]{ll}{\\!\\!-\log\\!\operatorname{det}(\bm{W})},&\\!{\text{
if }\bm{W}\succ 0,}\\\ \\!\\!{+\infty,}&\\!{\text{
otherwise,}}\end{array}\right.$

is given by [34, Example 24.66]

$\displaystyle\begin{aligned}
&\operatorname{prox}_{-\sigma^{-1}\mathrm{logdet}(\cdot)}(\bm{W})\\\
&=\frac{1}{2}\bm{Q}\\!\operatorname{diag}\\!\left(\mu_{1}\\!+\\!\sqrt{\mu_{1}^{2}\\!+\\!4\sigma^{-1}},\\!\cdots\\!,\mu_{n}\\!+\\!\sqrt{\mu_{n}^{2}\\!+\\!4\sigma^{-1}}\right)\\!\bm{Q}^{T}.\end{aligned}$

By (C) and the property [34, Proposition 24.8]:

$\mathrm{prox}_{f(\cdot+\bm{J})}(\bm{W})=\mathrm{prox}_{f}(\bm{W}+\bm{J})-\bm{J},$
(C.2 )

we obtain the result.

##  Appendix D Proof of Proposition 4

It is clear from (11) that $F$ is convex when
$\lambda_{2}\geq\gamma^{-1}\lambda_{1}$. The convergence condition of the
primal dual splitting method is given as follows [32]:

  1. 1.

$\frac{1}{\tau}\geq\sigma\|L\|^{2}+\frac{\beta}{2}$,

  2. 2.

$0<\rho_{k}<2-\frac{\beta}{2}\left(\frac{1}{\tau}-\sigma\|L\|^{2}\right)^{-1}$,

where $\beta$ is the Lipschitz constant of $\nabla F$ and
$\|L\|:=\sup_{\bm{w}\neq\bm{0}}\frac{\|L(\bm{w})\|_{\rm F}}{\|\bm{w}\|_{2}}$
is the operator norm. We show below the Lipschitz constant of $F$ and the
operator norm of $L$. (Although it is shown in [11] that $\|L\|=\sqrt{2n}$, we
show the proof for self-containedness.)

###  D-A Derivation of the Lipschitz constant of $\nabla$ $F$

The gradient of $F$ is given by

$\nabla
F(\bm{w})=\gamma^{-1}\lambda_{1}\mathrm{prox}_{\|\cdot\|_{1}}(\bm{w})+(\lambda_{2}-\gamma^{-1}\lambda_{1})\bm{w}.$

Hence, by the nonexpansiity of the proximity operator as well as the
triangular inequality of norm, we obtain

$\displaystyle\|\nabla F(\bm{w})-\nabla F(\bm{w}^{\prime})\|_{2}$
$\displaystyle=\left\|\left(\gamma^{-1}\lambda_{1}\mathrm{prox}_{\|\cdot\|_{1}}(\bm{w})\\!+\\!(\lambda_{2}-\gamma^{-1}\lambda_{1})\bm{w}\right)\right.$
$\displaystyle\quad\left.-\left(\gamma^{-1}\lambda_{1}\mathrm{prox}_{\|\cdot\|_{1}}(\bm{w}^{\prime})\\!+\\!(\lambda_{2}-\gamma^{-1}\lambda_{1})\bm{w}^{\prime}\right)\right\|_{2}$
(D.1 )
$\displaystyle\leq\|\gamma^{-1}\lambda_{1}\mathrm{prox}_{\|\cdot\|_{1}}(\bm{w})-\gamma^{-1}\lambda_{1}\mathrm{prox}_{\|\cdot\|_{1}}(\bm{w}^{\prime})\|_{2}$
$\displaystyle\quad+\|(\lambda_{2}-\gamma^{-1}\lambda_{1})(\bm{w}-\bm{w}^{\prime})\|_{2}$
(D.2 )
$\displaystyle\leq\gamma^{-1}\lambda_{1}\|\bm{w}-\bm{w}^{\prime}\|_{2}+(\lambda_{2}-\gamma^{-1}\lambda_{1})\|\bm{w}-\bm{w}^{\prime}\|_{2}$
(D.3 ) $\displaystyle\leq\lambda_{2}\|\bm{w}-\bm{w}^{\prime}\|_{2},$ (D.4 )

from which $\nabla F$ is $\lambda_{2}$-Lipschitz contiuous.

###  D-B Derivation of $\|L\|$

By definition of Laplacian, we have $[L(\bm{w})]_{i,i}=-\sum_{j\neq
i}[L(\bm{w})]_{i,j}$, and hence it holds that

$\displaystyle\|L\|^{2}=\sup_{\bm{w}\neq\bm{0}}\frac{\|L(\bm{w})\|_{\rm
F}^{2}}{\|\bm{w}\|_{2}^{2}}$
$\displaystyle=\sup_{\bm{w}\neq\bm{0}}\frac{\displaystyle\sum_{p=1}^{n}\left[\left(-\displaystyle\sum_{q\neq
p}[L(\bm{w})]_{p,q}\right)^{2}\right]+2\left(\sum_{p<q}[L(\bm{w})]_{p,q}^{2}\right)}{\displaystyle\sum_{p<q}[L(\bm{w})]_{p,q}^{2}}$
$\displaystyle=\sup_{\bm{w}\neq\bm{0}}\frac{\displaystyle\sum_{p=1}^{n}\left[\left(-\sum_{q\neq
p}[L(\bm{w})]_{p,q}\right)^{2}\right]}{\displaystyle\sum_{p<q}[L(\bm{w})]_{p,q}^{2}}+2.$
(D.5 )

Using the Caucy-Schwartz inequality, we have

$\displaystyle\sum_{p=1}^{n}\left[\left(-\sum_{q\neq
p}[L(\bm{w})]_{p,q}\right)^{2}\right]\\!$
$\displaystyle\\!\leq\\!\sum_{p=1}^{n}\left[(n-1)\sum_{q\neq
p}[L(\bm{w})]_{p,q}^{2}\right]$
$\displaystyle=2(n-1)\\!\sum_{p<q}[L(\bm{w})]_{p,q}^{2},$ (D.6 )

where the inequality holds with equality when $\bm{w}=\alpha\bm{1}$,
$\alpha\in\mathbb{R}$. By (D.5) and (D.6), we obtain $\|L\|^{2}=2n$, where the
upper bound is obtained when we consider the complete graph with all weights
equal to one; i.e., $\bm{w}=\bm{1}$.

## References

  * [1] H. Rue and L. Held  “Gaussian Markov random fields: theory and applications”  In _Chapman & Hall/CRC_, 2005
  * [2] Jerome Friedman, Trevor Hastie and Robert Tibshirani  “Sparse inverse covariance estimation with the graphical lasso”  In _Biostatistics_ 9.3 Oxford University Press, 2008, pp. 432–441
  * [3] Rahul Mazumder and Trevor Hastie  “The graphical lasso: New insights and alternatives”  In _Electronic Journal of Statistics_ 6 The Institute of Mathematical Statisticsthe Bernoulli Society, 2012, pp. 2125–2149  DOI: 10.1214/12-EJS740
  * [4] N. Meinshausen and P. Buhlmann  “High-dimensional graphs and variable selection with the Lasso”  In _The Annals of Statistics_ 34, 2006, pp. 1436–1462
  * [5] Onureena Banerjee, Laurent El Ghaoui and Alexandre d’Aspremont  “Model Selection Through Sparse Maximum Likelihood Estimation for Multivariate Gaussian or Binary Data”  In _Journal of Machine Learning Research_ 9.15, 2008, pp. 485–516  URL: http://jmlr.org/papers/v9/banerjee08a.html
  * [6] S. Segarra, A.. Marques, G. Mateos and A. Ribeiro  “Network Topology Inference from Spectral Templates”  In _IEEE Transactions on Signal and Information Processing over Networks_ 3.3, 2017, pp. 467–483  DOI: 10.1109/TSIPN.2017.2731051
  * [7] Hao Wang  “Bayesian Graphical Lasso Models and Efficient Posterior Computation”  In _Bayesian Analysis_ 7.4 International Society for Bayesian Analysis, 2012, pp. 867–886  DOI: 10.1214/12-BA729
  * [8] Shiliang Sun, Rongqing Huang and Ya Gao  “Network-Scale Traffic Modeling and Forecasting with Graphical Lasso and Neural Networks”  In _Journal of Transportation Engineering_ 138.11, 2012, pp. 1358–1367  DOI: 10.1061/(ASCE)TE.1943-5436.0000435
  * [9] Brenden Lake and Joshua Tenenbaum  “Discovering Structure by Learning Sparse Graphs”  In _Proceedings of 33rd Annual Conference of the Cognitive Science Society_ , 2010, pp. 778–783
  * [10] H.. Egilmez, E. Pavez and A. Ortega  “Graph Learning From Data Under Laplacian and Structural Constraints”  In _IEEE Journal of Selected Topics in Signal Processing_ 11.6, 2017, pp. 825–841  DOI: 10.1109/JSTSP.2017.2726975
  * [11] S. Kumar, J. Ying, J.. de Miranda Cardoso and D.. Palomar  “A unified framework for structured graph learning via spectral constraints”  In _Journal of Machine Learning Research_ 21, 2020, pp. 1–60
  * [12] Clifford Lam and Jianqing Fan  “Sparsistency and rates of convergence in large covariance matrix estimation”  In _Annals of Statistics_ 37.6B The Institute of Mathematical Statistics, 2009, pp. 4254–4278  DOI: 10.1214/09-AOS720
  * [13] Fabrizio Vecchio, Francesca Miraglia and Paolo Maria Rossini  “Connectome: Graph theory application in functional brain network architecture”  In _Clinical Neurophysiology Practice_ 2, 2017, pp. 206–213  DOI: https://doi.org/10.1016/j.cnp.2017.09.003
  * [14] O. Mason and M. Verwoerd  “Graph theory and networks in Biology”  In _IET Systems Biology_ 1.30 Institution of EngineeringTechnology, 2007, pp. 89–119  URL: https://digital-library.theiet.org/content/journals/10.1049/iet-syb_20060038
  * [15] M.. Bhuyan, D.. Bhattacharyya and J.. Kalita  “Network Anomaly Detection Methods, Systems and Tools”  In _IEEE Communications Surveys Tutorials_ 16.1, 2014, pp. 303–336
  * [16] G. Mateos, S. Segarra, A.. Marques and A. Ribeiro  “Connecting the Dots: Identifying Network Structure via Graph Signal Processing”  In _IEEE Signal Processing Magazine_ 36.3, 2019, pp. 16–43  DOI: 10.1109/MSP.2018.2890143
  * [17] X. Dong, D. Thanou, M. Rabbat and P. Frossard  “Learning Graphs From Data: A Signal Representation Perspective”  In _IEEE Signal Processing Magazine_ 36.3, 2019, pp. 44–63  DOI: 10.1109/MSP.2018.2887284
  * [18] Daphne Koller and Nir Friedman  “Probabilistic Graphical Models: Principles and Techniques - Adaptive Computation and Machine Learning”  The MIT Press, 2009
  * [19] P. Danaher, P. Wang and D.. Witten  “The joint graphical lasso for inverse covariance estimation across multiple classes”  In _Journal of Royal Statistical Society. Series B: Statistical Methodological_ 76.2, 2014, pp. 373–397
  * [20] P. Buhlmann and S. Geer  “Statistics for High-Dimensional Data:Methods, Theory and Applications”, 2011
  * [21] Cun Hui Zhang  “Nearly unbiased variable selection under minimax concave penalty”  In _Annals of Statistics_ 38.2 The Institute of Mathematical Statistics, 2010, pp. 894–942  DOI: 10.1214/09-AOS729
  * [22] Fei Wen, Lei Chu, Peilin Liu and Robert C. Qiu  “A Survey on Nonconvex Regularization-Based Sparse and Low-Rank Recovery in Signal Processing, Statistics, and Machine Learning”  In _IEEE Access_ 6, 2018, pp. 69883–69906  DOI: 10.1109/ACCESS.2018.2880454
  * [23] Rahul Mazumder, Jerome H. Friedman and Trevor Hastie  “SparseNet: Coordinate Descent With Nonconvex Penalties”  In _Journal of the American Statistical Association_ 106.495 Taylor & Francis, 2011, pp. 1125–1138  DOI: 10.1198/jasa.2011.tm09738
  * [24] I. Selesnick  “Sparse Regularization via Convex Analysis”  In _IEEE Transactions on Signal Processing_ 65.17, 2017, pp. 4481–4494  DOI: 10.1109/TSP.2017.2711501
  * [25] Jiro Abe, Masao Yamagishi and Isao Yamada  “Linearly involved generalized Moreau enhanced models and their proximal splitting algorithm under overall convexity condition”  In _Inverse Problems_ 36.3 IOP Publishing, 2020, pp. 1–36  DOI: 10.1088/1361-6420/ab551e
  * [26] E.A. Nurminskii  “The quasigradient method for the solving of the nonlinear programming problems”  In _Cybernetics and Systems Analysis_ 9, 1973, pp. 145–150
  * [27] Léa Laporte et al.  “Nonconvex Regularizations for Feature Selection in Ranking With Sparse SVM”  In _IEEE Transactions on Neural Networks and Learning Systems_ 25.6, 2014, pp. 1118–1130  DOI: 10.1109/TNNLS.2013.2286696
  * [28] Shibin Wang et al.  “Nonconvex Sparse Regularization and Convex Optimization for Bearing Fault Diagnosis”  In _IEEE Transactions on Industrial Electronics_ 65.9, 2018, pp. 7332–7342  DOI: 10.1109/TIE.2018.2793271
  * [29] X. Shen, W. Pan and Y. Zhu  “Likelihood-based selection and sharp parameter estimation”  In _Journal of the American Statistical Association_ 107.497, 2012, pp. 223–232  DOI: 10.1080/01621459.2011.645783
  * [30] Jiaxi Ying, José Vinı́cius Miranda Cardoso and Daniel Palomar  “Nonconvex Sparse Graph Learning under Laplacian Constrained Graphical Model”  In _Advances in Neural Information Processing Systems_ 33, 2020, pp. 7101–7113  URL: https://proceedings.neurips.cc/paper/2020/file/4ef42b32bccc9485b10b8183507e5d82-Paper.pdf
  * [31] Yangjing Zhang, Kim Chuan Toh and Defeng Sun  “Learning Graph Laplacian with MCP”  In _ArXiv e-prints_ , 2020  eprint: 2010.11559
  * [32] L Condat  “A Primal–Dual Splitting Method for Convex Optimization Involving Lipschitzian, Proximable and Linear Composite Terms”  In _Journal of Optimization Theory and Applications_ 158, 2013, pp. 460–479  DOI: 10.1007/s10957-012-0245-9
  * [33] Tatsuya Koyakumaru, Masahiro Yukawa, Eduardo Pavez and Antonio Ortega  “A Graph Learning Algorithm Based on Gaussian Markov Random Fields and Minimax Concave Penalty”  In _Proceedings of 46th International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2021, pp. 5390–5394
  * [34] Heinz H. Bauschke and Patrick L. Combettes  “Convex Analysis and Monotone Operator Theory in Hilbert Spaces”, 2017
  * [35] T.  “Graph learning based on Gaussian Markov random fields and minimax concave penalty” Bachelor Thesis (in Japanese), 2020  Keio University
  * [36] K. and M.  “Robust recovery of jointly-sparse signals using minimax concave loss function”  In _IEEE Transactions on Signal Process._ 69, 2020, pp. 669–681
  * [37] Hiroyuki Kaneko and Masahiro Yukawa  “Normalized least-mean-square algorithms with minimax concave penalty”  In _Proceedings of 45th International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2020, pp. 5440–5444
  * [38] M., H., K. and I.  “Linearly-involved Moreau-Enhanced-over-Subspace Model: Debiased Sparse Modeling and Stable Outlier-Robust Regression” submitted for publication, 2021
  * [39] Kei Komuro, Masahiro Yukawa and Renato L.. Cavalcante  “Distributed Sparse Optimization With Minimax Concave Regularization”  In _2021 IEEE Statistical Signal Processing Workshop (SSP)_ , 2021, pp. 31–35
  * [40] L. Le Magoarou, R. Gribonval and N. Tremblay  “Approximate Fast Graph fourier Transforms via Multilayer Sparse Approximations”  In _IEEE Transactions on Signal and Information Processing over Networks_ 4.2, 2018, pp. 407–420  DOI: 10.1109/TSIPN.2017.2710619
  * [41] David Hallac, Youngsuk Park, Stephen Boyd and Jure Leskovec  “Network Inference via the Time-Varying Graphical Lasso”  In _Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2017, pp. 205–213  DOI: 10.1145/3097983.3098037
  * [42] Koki Yamada, Yuichi Tanaka and Antonio Ortega  “Time-varying Graph Learning Based on Sparseness of Temporal Variation”  In _Proceedings of 44th International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2019, pp. 5411–5415  DOI: 10.1109/ICASSP.2019.8682762
  * [43] Nathanaël Perraudin et al.  “GSPBOX: A toolbox for signal processing on graphs”  In _ArXiv e-prints_ , 2014  arXiv:1408.5781 [cs.IT]
  * [44] Charles Kemp and Joshua B. Tenenbaum  “The discovery of structural form”  In _Proceedings of the National Academy of Sciences_ 105.31 National Academy of Sciences, 2008, pp. 10687–10692  DOI: 10.1073/pnas.0802631105
  * [45] Masahiro Yukawa and Hideaki Kagami  “Supervised nonnegative matrix factorization via minimization of regularized Moreau-envelope of divergence function with application to music transcription”  In _Journal of the Franklin Institute_ 355.4 Elsevier Limited, 2018, pp. 2041–2066  DOI: 10.1016/j.jfranklin.2017.12.002