# Efficient Graph Laplacian Estimation by Proximal Newton††thanks: This
research was supported by the Israel Ministry of National Infrastructure,
Energy, and Water Resources and by the ISRAEL SCIENCE FOUNDATION (grant No.
1148/22). This work was also supported in part by the Israeli Council for
Higher Education (CHE) via the Data Science Research Center, Ben-Gurion
University of the Negev, Israel.

Yakov Medvedovsky School of Electrical and Computer Engineering, Ben-Gurion
University of the Negev. Eran Treister Department of Computer Science, Ben-
Gurion University of the Negev. Tirza Routtenberg†

###### Abstract

The Laplacian-constrained Gaussian Markov Random Field (LGMRF) is a common
multivariate statistical model for learning a weighted sparse dependency graph
from given data. This graph learning problem can be formulated as a maximum
likelihood estimation (MLE) of the precision matrix, subject to Laplacian
structural constraints, with a sparsity-inducing penalty term. This paper aims
to solve this learning problem accurately and efficiently. First, since the
commonly used $\ell_{1}$-norm penalty is inappropriate in this setting and may
lead to a complete graph, we employ the nonconvex minimax concave penalty
(MCP), which promotes sparse solutions with lower estimation bias. Second, as
opposed to existing first-order methods for this problem, we develop a second-
order proximal Newton approach to obtain an efficient solver, utilizing
several algorithmic features, such as using Conjugate Gradients,
preconditioning, and splitting to active/free sets. Numerical experiments
demonstrate the advantages of the proposed method in terms of both
computational complexity and graph learning accuracy compared to existing
methods.

##  1 Introduction

Graphs are fundamental mathematical structures used in various fields to
represent data, signals, and processes. Weighted graphs naturally arise when
processing networked data applications, such as computer, social, sensor,
energy, transportation, and biological networks Newman, (2010), where the data
is inherently related to a graph associated with the underlying network. As a
result, graph learning is used in a wide range of applications, both for
virtual, data-based networks and for physical networks of infrastructures; the
latter include the brain Vecchio et al., (2017), financial analysis Giudici
and Spelta, (2016), and electrical networks Deka et al., (2018). In this
context, a fundamental graph learning problem is the estimation of the
Laplacian matrix, which plays a central role in spectral graph theory Chung,
(1997) and machine learning Von Luxburg, (2007). Moreover, in the field of
graph signal processing (GSP) Shuman et al., (2013), the Laplacian matrix is
used to extend basic signal processing operations, such as filtering Milanfar,
(2012), sampling Anis et al., (2016), and signal recovery Narang and Ortega,
(2012); Kroizer et al., (2022). Therefore, accurate and efficient Laplacian
learning is extremely valuable for data processing in networks.

Gaussian Markov Random Fields (GMRFs) Rue and Held, (2005) are probabilistic
graphical models that have been studied extensively over the last two decades
Banerjee et al., (2006); Hsieh et al., (2013, 2011). In particular, there has
been a significant effort to develop approaches for learning sparse precision
matrices under GMRF models, where this precision represents a graph topology.
The graphical LASSO (GLASSO) Friedman et al., (2008) is a widely-used approach
to estimate the sparse precision matrix using $\ell_{1}$-norm regularization.
Laplacian-constrained GMRF (LGMRF) models are a specific type of GMRF, in
which the precision is a Laplacian matrix. Laplacian constraints imply that
smooth signals, which have similar values for nodes connected by large
weights, have a higher probability Ying et al., 2020b . Smooth graph signals
are common in various network applications, such as power systems Dabush et
al., (2023). Thus, LGMRF models have been widely studied in semi-supervised
learning frameworks Zhu et al., (2003), to analyze real-world datasets Cardoso
and Palomar, (2020), and in GSP Dong et al., (2016); Egilmez et al., (2017).

The estimation of standard and Laplacian-constrained GMRF models are closely
related. However, the extension of GLASSO and similar GMRF estimation
algorithms to the case of LGMRF, which also considers the Laplacian
constraints, is not straightforward. For example, since the Laplacian matrix
is sparse in many applications, then similarly to the GLASSO approach,
$\ell_{1}$-norm regularization methods were used to promote the Laplacian
sparsity Egilmez et al., (2017); Kumar et al., (2019); Liu et al., (2019);
Zhao et al., (2019). However, it was recently shown that the $\ell_{1}$-norm
is an inappropriate penalty for promoting the sparsity of the precision matrix
under the LGMRF model, as it leads to an inaccurate recovery of the
connectivity pattern of the graph Ying et al., 2020a . Thus, to properly
address the Laplacian estimation problem, other sparsity-promoting penalty
terms should be used Ying et al., 2020a ; Ying et al., 2020b . However, the
existing methods are of first-order (i.e., gradient-based), while in general
GMRF estimation literature Oztoprak et al., (2012); Hsieh et al., (2013,
2011); Treister and Turek, (2014) it is known that second-order methods, i.e.,
proximal Newton (Lee et al.,, 2014), yield better performance in terms of run
time. In such approaches, a quadratic approximation is applied on the smooth
part, while the non-smooth regularization term remains intact. This exploits a
connection between the gradient and Hessian of the MLE objective and hence is
efficient for GMRF estimation. Therefore, there is a need for second-order
methods for Laplacian learning that utilize proper regularization (to yield
accurate sparsity patterns) and are efficient in run time.

Contribution. In this paper, we present an efficient graph estimation method
for the LGMRF model based on the proximal Newton approach. Our method involves
solving a nonconvex penalized MLE problem with MCP penalty, instead of the
commonly-used $\ell_{1}$-norm penalty. At each iteration, the smooth part of
the objective is approximated by a second-order Taylor expansion around the
current precision matrix, while the non-smooth penalty and Laplacian
constraints remain intact. The resulting Newton problem is a penalized
quadratic minimization under linear constraints to guarantee that the next
iterate is a Laplacian matrix. To the best of our knowledge, the proposed
method is the first proximal Newton method for the LGMRF model estimation.
Furthermore, our framework includes two algorithmic novelties that are
nontrivial in this context. First, the inner, constrained, Newton problem is
solved by a projected nonlinear conjugate gradient method. This yields a
significant computational speedup over gradient-based first-order methods.
Second, we introduce a diagonal preconditioner to improve the performance
further. We show some theoretical results regarding the problem and the
algorithms we propose. Furthermore, we demonstrate the effectiveness of the
proposed method in learning sparse graphs through numerical experiments.

Organization and notation. The rest of this paper is organized as follows. In
Section 2, we describe the considered Laplacian learning problem, simplify its
formulation, and discuss the sparsity-promoting penalty function. In Section
3, we develop the proposed NewGLE method. Simulations are shown in Section 5
and the paper is concluded in Section 6.

In the rest of this paper we denote vectors and matrices in boldface lowercase
and uppercase letters, respectively. The elements of the vectors ${\bf{1}}$
and ${\bf{0}}$ are ones and zeros, respectively. The notations
$\|\cdot\|_{\textbf{F}}$, $|\cdot|$, $\otimes$, $(\cdot)^{T}$, $(\cdot)^{-1}$,
$(\cdot)^{\dagger}$, $|\cdot|_{+}$, and $\text{Tr}(\cdot)$ denote the
Frobenius norm, determinant operator, Kronecker product, transpose, inverse,
Moore-Penrose pseudo-inverse, pseudo-determinant, and trace, respectively. The
sets $\mathcal{S}^{p}$ and $\mathcal{S}^{p}_{+}$ are of real symmetric and
real positive semi-definite $p\times p$ matrices, respectively. For matrix
${\bf{A}}\in\mathbb{R}^{p\times p}$, $\nabla_{{\bf{A}}}f\in\mathbb{R}^{p\times
p}$ and $\nabla_{{\bf{A}}}^{2}f\in\mathbb{R}^{p^{2}\times p^{2}}$ are the
gradient and Hessian of the scalar function $f({\bf{A}})$.

##  2 Problem Formulation

In this section, we formulate the considered Laplacian learning problem. We
start with definitions related to GSP. Then, we present the LGMRF and
formulate the graph estimation problem under the LGMRF and MLE in Subsection
2.1. The sparsity-promoting penalty function is discussed in Subsection 2.2.

We consider an undirected, connected, and weighted graph
$\mathcal{G}\left({\mathcal{V},\mathcal{E}}\right)$, where
$\mathcal{V}=\left\\{{1,\dots,p}\right\\}$ and
$\mathcal{E}=\left\\{{1,\dots,m}\right\\}$ are the set of vertices and the set
of edges, respectively. The nonnegative weighted adjacency matrix of the
graph, ${\bf{W}}$, has the $(i,j)$-th element $w_{i,j}>0$ if vertex $i$ and
$j$ are connected, and zero otherwise. From spectral graph theory Chung,
(1997), the rank of the Laplacian matrix for a connected graph with $p$ nodes
is $p-1$. Thus, the set of Laplacian matrices for connected graphs can be
defined as Ying et al., 2020b

$\displaystyle\mathcal{L}=\left\\{{\bf{L}}\in\mathcal{S}^{p}_{+}\Big{|}\begin{array}[]{ll}{\bf{L}}_{i,j}\leq
0,\forall i\neq j,~{}i,j=1,\ldots,p\\\
{\bf{L}}{\bf{1}}={\bf{0}},~{}~{}{\textrm{rank}}\left({{\bf{L}}}\right)=p-1\end{array}\right\\}.$
(3)

###  2.1 GMRF and Constrained MLE

In this paper, we use the formulation of the graph estimation problem based on
the probabilistic graphical model Koller and Friedman, (2009); Banerjee et
al., (2008). We data samples from a zero-mean Gaussian distribution
parametrized by a positive semi-definite precision matrix, ${\bf{L}}$, i.e.,
${\bf{x}}\sim\mathcal{N}\left({{\bf{0}},{\bf{L}}^{\dagger}}\right)$, where
$\mathcal{N}\left({\boldsymbol{\mu},\boldsymbol{\Sigma}}\right)$ denotes the
normal distribution with mean $\boldsymbol{\mu}$ and covariance matrix
$\boldsymbol{\Sigma}$. This defines a Laplacian-constrained improper LGMRF
model Kumar et al., (2019) with parameters $\left({{\bf{0}},{\bf{L}}}\right)$,
where the precision matrix is ${\bf{L}}\in\mathcal{L}$ and $\mathcal{L}$ is
defined in (3).

Given $n$ independent and identically distributed (i.i.d.) samples
${\bf{x}}_{1},\ldots,{\bf{x}}_{n}$ drawn from an LGMRF, our goal is to find
the constrained MLE of ${\bf{L}}$ based on these samples and under the
Laplacian constraints. The MLE can be found by solving the following
constrained minimization of the negative log-likelihood:

$\hat{{\bf{L}}}_{MLE}=\underset{{\bf{L}}\in\mathcal{L}}{\text{argmin}}\left\\{{{\rm
Tr}{\left({{\bf{L}}{\bf{S}}}\right)}-\log|{\bf{L}}|_{+}}\right\\},$ (4)

where
${\bf{S}}\stackrel{{\scriptstyle\triangle}}{{=}}\frac{1}{n}\sum_{i=1}^{n}{\bf{x}}_{i}{\bf{x}}^{T}_{i}$
is the empirical covariance matrix.

Sparsity plays an important role in high-dimensional learning, which helps
avoid over-fitting and improves the identification of the relationships among
data, especially where ${\bf{S}}$ is low rank. In particular, the Laplacian is
a sparse matrix in various applications, e.g. electrical networks Halihal and
Routtenberg, (2022). A sparse graph estimation problem under the LGMRF model
can be formulated by adding a sparse penalty function to the estimator in (4),
which results in

$\displaystyle\begin{aligned}
&\underset{{\bf{L}}\in\mathcal{L}}{\text{minimize}}&{\rm
Tr}{\left({{\bf{L}}{\bf{S}}}\right)}-\log|{\bf{L}}|_{+}+\rho\left({{\bf{L}};\lambda}\right),\end{aligned}$
(5)

where $\rho\left({{\bf{L}};\lambda}\right)$ is a general sparsity-promoting
penalty function and $\lambda$ is a tuning parameter.

The objective function in (5) (or in (4)) involves the pseudo-determinant
term, since the target matrix ${\bf{L}}$ is a singular matrix. As a result,
the optimization problem becomes hard to solve Holbrook, (2018). As a remedy
to this problem, it is demonstrated in Egilmez et al., (2017) that if
${\bf{L}}\in\mathcal{L}$, then we can use the relation

$|{\bf{L}}|_{+}=|{\bf{L}}+{\bf{J}}|,$ (6)

where ${\bf{J}}=\frac{1}{p}{\bf{1}}{\bf{1}}^{T}$. By substituting (6) in (5),
the problem is simplified. In addition, we show that the constraints in (3)
are redundant, and simplify them in the following theorem whose proof appears
in Appendix A.

######  Theorem 1.

The optimization problem

$\displaystyle\underset{{\bf{L}}\in\mathcal{F}}{\text{minimize}}$
$\displaystyle{\rm
Tr}{\left({{\bf{L}}{\bf{S}}}\right)}-\log|{\bf{L}}+{\bf{J}}|+\rho\left({{\bf{L}};\lambda}\right),$
(7)

where the feasible set is given by

$\displaystyle\mathcal{F}=\left\\{{\bf{L}}\in\mathcal{S}^{p}\Big{|}\begin{array}[]{ll}{\bf{L}}_{i,j}\leq
0,\forall i\neq j,~{}i,j=1,\ldots,p,\\\
{\bf{L}}{\bf{1}}={\bf{0}}\end{array}\right\\}.$ (10)

is equivalent to the optimization problem stated in (5).

###  2.2 Sparsity Promoting Penalty Functions

Recent works like (Egilmez et al.,, 2017) introduced the $\ell_{1}$-norm
penalized MLE under the LGMRF model to estimate a sparse graph. The choice of
$\rho$ in (7) as $\|{\bf{L}}\|_{1,\text{off}}$, the absolute sum of all off-
diagonal elements of ${\bf{L}}$, results in the same objective function as in
the well-known GLASSO problem but where the feasible set is restricted to
${\bf{L}}\in\mathcal{F}$. However, unlike in GLASSO, under Laplacian
constraints the $\ell_{1}$-norm is ineffective in promoting a sparse solution.
In particular, it is shown by Ying et al., 2020a  that substituting
$\rho\left({{\bf{L}};\lambda}\right)=\lambda\|{\bf{L}}\|_{1,\text{off}}$ in
the objective of (7) and choosing the penalty parameter $\lambda$ to be large
enough leads to a fully-connected (complete) graph as the solution. Moreover,
numerical experiments in (Ying et al., 2020a, ) show that the number of the
estimated positive edges increases as $\lambda$ increases. An intuitive
explanation for this phenomenon is that the $\ell_{1}$-norm promotes sparsity
uniformly on all off-diagonal entrances of the matrix. Hence, for $\lambda>0$,
minimizing the $\ell_{1}$-norm separately with a positive definite constraint,
as in GLASSO, leads to a diagonal matrix with a positive diagonal. However,
unlike a regular GMRF model, the Laplacian constraints,
${\bf{L}}{\bf{1}}={\bf{0}}$ and ${\bf{L}}_{i,j}\leq 0,$ for $i\neq j$, relate
the off-diagonal entries with the diagonal entries of ${\bf{L}}$. Hence,
minimizing the $\ell_{1}$-norm would lead to the zeros matrix. Therefore, as
$\lambda$ increases, the eigenvalues of ${\bf{L}}$ get closer to zero, and the
$-\log|\cdot|$ term in the objective function increases significantly. Since
$\ell_{1}$-norm penalizes all off-diagonal elements in a uniform way, to
prevent a significant increase of the objective function, the minimization
process inserts relatively small (in absolute value) negative numbers, even
for zero elements in the ground-truth Laplacian matrix. This way, both the
$\ell_{1}$-norm and the expression $-\log|\cdot|$ in the optimization problem
remain relatively small.

Alternative approaches using nonconvex penalties, such as clipped absolute
deviation (SCAD) Fan and Li, (2001) and the minimax concave penalty (MCP)
Zhang, 2010a , have been proposed to alleviate this issue. These penalty
functions penalize the elements in the precision matrix more selectively than
the $\ell_{1}$-norm, such that elements that receive significant values from
the MLE will receive insignificant penalties, while elements closer to zero
will receive a greater penalty. In this paper, we use the MCP penalty
function, which is applied on the off-diagonal entries of ${\bf{L}}$ as
follows:

$\rho_{MCP}\left({{\bf{L}};\lambda}\right)=\sum\nolimits_{i\neq
j}MCP\left({{\bf{L}}_{ij};\gamma,\lambda}\right),$ (11)

where the $MCP(\cdot)$ is a scalar function, defined as

$\displaystyle
MCP(x;\gamma,\lambda)=\left\\{\begin{array}[]{lr}\lambda|x|-\frac{x^{2}}{2\gamma}\quad\textrm{if}\quad|x|\leq\gamma\lambda\\\
\frac{1}{2}\gamma\lambda^{2}\quad\quad\quad\textrm{if}\quad|x|>\gamma\lambda\end{array}\right..$
(14)

The function is constant and has a derivative of 0 for values larger than
$\gamma\lambda$ in magnitude, hence only elements smaller than $\gamma\lambda$
are affected by $\rho_{MCP}()$.

To conclude, our goal in this paper is to develop an efficient algorithm for
finding the constrained MLE of the Laplacian matrix, which is the precision
matrix of the LGMRF model. This is done by solving the optimization problem in
(7) with
$\rho\left({{\bf{L}};\lambda}\right)=\rho_{MCP}\left({{\bf{L}};\lambda}\right)$,
defined in (11). I.e., the problem we solve is given by:

$\displaystyle\begin{aligned}
&\underset{{\bf{L}}\in\mathcal{F}}{\text{minimize}}&{\rm
Tr}{\left({{\bf{L}}{\bf{S}}}\right)}-\log|{\bf{L}}+{\bf{J}}|+\rho_{MCP}\left({{\bf{L}};\lambda}\right)\\\
\end{aligned},$ (15)

where $\mathcal{F}$ is defined in (10).

##  3 Method

In this section, we develop the proposed proximal Newton method for the graph
Laplacian estimation, named NewGLE, that aims to solve (15) efficiently using
a proximal Newton approach adopted to the problem. First, we define the
constrained Newton optimization problem for the graph learning in Subsection
3.1, keeping the Laplacian constraints when finding the Newton direction.
Then, in Subsection 3.2, we present the approach that maps a Laplacian matrix
to a vector and simplifies the set of constraints. Next, we describe the idea
of restricting the Newton direction into a “free set” in Subsection 3.3. This
constitutes the NewGLE method, which is summarized in Algorithm 1. Following
that, we present our inner solver for finding the Newton direction, which is
composed of the nonlinear projected conjugate gradient that is used together
with a diagonal preconditioner, in Subsection 3.4 and Appendix C. Finally, we
provide an algorithm (Algorithm 2 in the Appendix) for finding the Newton
direction.

###  3.1 Proximal Newton for Graph Laplacian Estimation

In the “proximal Newton” approach, a quadratic approximation is applied on the
smooth part of the objective function, while leaving the non-smooth term
intact, in order to obtain the Newton descent direction. This approach is
considered to be among the state-of-the-art methods for solving the GLASSO
problem Rolfs et al., (2012); Hsieh et al., (2013, 2011); Mazumder and Hastie,
(2012); Treister et al., (2016), which is highly related to our problem in
(15). The advantage of this method lies in the treatment of the $\log\det$
term (as we show next), which appears in both problems. However, the GLASSO
methods cannot be applied to our graph learning problems because we aim to
estimate a precision matrix that satisfies the Laplacian constraints (and is
singular). In contrast, the learned precision matrix under the GMRF model in
GLASSO is a general positive definite matrix. In addition, we consider the MCP
penalty, while the above-mentioned methods consider the $\ell_{1}$-norm
penalty, which is not suitable for our case, as discussed in Subsection 2.2.

In our proximal Newton approach we design a constrained Newton problem, where
the smooth part of the objective function in (15) is

$f\left({{\bf{L}}}\right)={\rm
Tr}{\left({{\bf{L}}{\bf{S}}}\right)}-\log{|{\bf{L}}+{\bf{J}}|},$ (16)

and the penalty function, $\rho_{MCP}\left({{\bf{L}};\lambda}\right)$, is the
non-smooth term. At the $t$-th iteration of the proximal Newton approach, the
smooth part of the objective is approximated by a second-order Taylor
expansion around the current estimation ${\bf{L}}^{(t)}$. To this end, we use
the gradient and Hessian of $f({\bf{L}})$ that are given by Hsieh et al.,
(2014); Chin et al., (2013)

$\nabla_{{\bf{L}}}f={\bf{S}}-{\bf{Q}},\quad\nabla_{{\bf{L}}}^{2}f={\bf{Q}}\otimes{\bf{Q}},$
(17)

where $\otimes$ is the Kronecker product and
${\bf{Q}}=\left({{\bf{L}}+{\bf{J}}}\right)^{-1}$.

The gradient in (17) already shows the main difficulty in solving the
optimization problem in (15): it contains ${\bf{Q}}$, the inverse of the
matrix ${\bf{L}}+{\bf{J}}$, which is expensive to compute. The advantage of
the proximal Newton approach for this problem is the low overhead: by
calculating ${\bf{Q}}$ in $\nabla_{{\bf{L}}}f({\bf{L}})$, we also get part of
the Hessian computation at the same cost.

In the same spirit of the proximal Newton approach, in addition to the penalty
term, we also keep the Laplacian constraint intact in the Newton problem. That
is, at iteration $t$, the Newton direction $\Delta^{(t)}$ solves the
constrained penalized quadratic minimization problem:

$\displaystyle\underset{\Delta\in\mathcal{S}^{p}}{\text{minimize}}$
$\displaystyle f({\bf{L}}^{(t)})+{\rm
Tr}{\left({\Delta\left({{\bf{S}}-{\bf{Q}}}\right)}\right)}+\frac{1}{2}{\rm
Tr}{\left({\Delta{\bf{Q}}\Delta{\bf{Q}}}\right)}$ (18)
$\displaystyle{+\frac{1}{2}\|\varepsilon\cdot\Delta\|_{F}^{2}}+\rho_{MCP}({\bf{L}}^{(t)}+\Delta;\lambda)$
s.t. $\displaystyle{\bf{L}}^{(t)}+\Delta\in\mathcal{F},$

where ${\bf{Q}}=\left({{\bf{L}}^{(t)}+{\bf{J}}}\right)^{-1}$. It should be
noted that for the sake of simplicity, we use the notation ${\bf{Q}}$ instead
of ${\bf{Q}}^{(t)}$, although ${\bf{Q}}$ is updated with the iteration $t$.
Note that the gradient and Hessian of $f(\cdot)$ in (17) at the iterate
${\bf{L}}^{(t)}$ are featured in the second and third terms in (18),
respectively. The first term of (18) is constant and can be ignored.
Altogether, the objective function in (18) is quadratic with the MCP sparsity-
promoting penalty. In addition, we have a constraint to guarantee that the
next update of ${\bf{L}}$ is feasible, i.e.,
$\left({{\bf{L}}^{(t)}+\Delta}\right)\in\mathcal{F}$. The fourth term, with a
small and symmetric weight matrix $\varepsilon\geq 0$ penalizing the squares
of the off-diagonals of $\Delta$. This term is used for the stabilization of
the iterative solution of (18), and to make sure that the Hessian of (18) is
positive definite, since $\rho_{MCP}$ is concave. The advantage of solving
(18) with $\rho_{MCP}$ intact, as opposed to majorizing $\rho_{MCP}$ (e.g.,
using its gradient instead) is that the underlying gradient-based iterative
solver of (18) has access to the curvature of $\rho_{MCP}$ and to its behavior
in the sub-sections in its definition (14).

Once the Newton problem is approximately solved, yielding the direction
$\Delta^{(t)}$, it is added to ${\bf{L}}^{(t)}$ employing a linesearch
procedure to sufficiently reduce the objective in (15). To this end, the
updated iterate is

${\bf{L}}^{(t+1)}={\bf{L}}^{(t)}+\alpha^{(t)}\Delta^{(t)},$ (19)

where the parameter $\alpha^{(t)}$ is obtained by a backtracking linesearch
using Armijo’s rule Armijo, (1966). We state the following theorem on such a
solution process, whose proof appears in Appendix D:

######  Theorem 2.

Let $\\{{\bf{L}}^{(t)}\\}$ be a series of points produced by a proximal
quadratic minimization, constrained by $\mathcal{F}$ with some SPD matrices
$0\prec\lambda^{{\bf{M}}}_{min}{\bf{I}}\preceq{\bf{M}}^{(t)}\preceq\lambda^{{\bf{M}}}_{max}{\bf{I}}$
as the Hessian (as in (48)), followed by the linesearch (19), starting from
${\bf{L}}^{(0)}\in\mathcal{F}$. Then any limit point $\bar{\bf{L}}$ of
$\\{{\bf{L}}^{(t)}\\}$ is a stationary point of (15).

In the next section, we simplify problem (18) to ease the treatment of the
constraints.

###  3.2 Laplacian parameterization

To simplify the Laplacian structural constraints in (18), we use a linear
operator that maps a vector ${\bf{w}}\in\mathbb{R}^{p(p-1)/2}$ to a matrix
$\mathcal{P}\left({{{\bf{w}}}}\right)\in\mathbb{R}^{p\times p}$.

######  Definition 1.

The linear operator
$\mathcal{P}:\mathbb{R}^{p(p-1)/2}\mapsto\mathbb{R}^{p\times p}$ is defined as
Ying et al., 2020a

$\displaystyle\mathcal{P}\left({{{\bf{w}}}}\right)_{i,j}=\left\\{\begin{array}[]{lr}-w_{k}&i>j\\\
\mathcal{P}\left({{{\bf{w}}}}\right)_{j,i}&i<j\\\ -\sum_{j\neq
i}\mathcal{P}\left({{{\bf{w}}}}\right)_{i,j}&i=j\end{array},\right.$ (23)

where $k=i-j+\frac{j-1}{2}\left({2p-j}\right)$ for $i>j$.

The operator $\mathcal{P}$ defines a linear mapping from a nonnegative weight
vector ${\bf{w}}\geq 0$ to a Laplacian matrix
$\mathcal{P}\left({{{\bf{w}}}}\right)\in\mathcal{F}$, defined in (10). Using
this parametrization, we rewrite the optimization problem in (18) such that
the decision variable is a vector instead of a matrix. In particular, in
Appendix B, we show that (18) is equivalent to the following simpler vector-
parameterized Newton problem that we solve

$\displaystyle\begin{aligned}
\underset{\delta\in\mathbb{R}^{p(p-1)/2}}{\text{minimize}}\quad&f_{N}(\delta)\stackrel{{\scriptstyle\triangle}}{{=}}{\rm
Tr}{\left({\mathcal{P}\left({{\delta}}\right)\left({{\bf{S}}-{\bf{Q}}}\right)}\right)}\;+{\|\tilde{\varepsilon}\odot\delta\|_{2}^{2}+}\\\
\quad&\frac{1}{2}{\rm
Tr}{\left({\mathcal{P}\left({{\delta}}\right){\bf{Q}}\mathcal{P}\left({{\delta}}\right){\bf{Q}}}\right)}+2\tilde{\rho}_{mcp}\left({{\bf{w}}+\delta;\lambda}\right)\\\
\textrm{s.t.}\quad&\delta\geq-{\bf{w}},\end{aligned}$ (24)

where $\tilde{\varepsilon}\in\mathcal{R}^{p(p-1)/2}$ is the vectorized upper
triangle of $\varepsilon$ (without the diagonal) and $\tilde{\rho}_{MCP}$ is
given by

$\tilde{\rho}_{MCP}\left({{{\bf{w}}};\lambda}\right)\stackrel{{\scriptstyle\triangle}}{{=}}\sum_{k=1}^{p\left({p-1}\right)/2}MCP\left({w_{k};\gamma,\lambda}\right)$
(25)

and $MCP(x;\gamma,\lambda)$ is defined in (14). The only remaining constraint
in (24) is an inequality constraint, which simplifies the projection onto the
feasible space.

###  3.3 Restricting the updates to free sets

To ease the minimization of (15) we limit the minimization of the Newton
problem to a “free set” of variables, while keeping the rest as zeros. This
idea was suggested in (Hsieh et al.,, 2011) for the GLASSO problem that
includes the $\ell_{1}$-norm penalty. The free set of a matrix ${\bf{L}}$ is
defined as

$\displaystyle\textbf{Free}\left({{\bf{L}}}\right)=\left\\{{k:{\bf{L}}_{i,j}\neq
0\lor\left[{{\bf{S}}-{\bf{Q}}}\right]_{i,j}>\lambda}\right\\},$ (26)

where $\lor$ denotes the logistic OR operator and the relation between $(i,j)$
and $k$ appears in Definition (1). This set comes from the subgradient of the
$\ell_{1}$-norm term, in addition to the non-positivity constraint—the zero
elements that are not in $\textbf{Free}\left({{\bf{L}}}\right)$ will remain
zero after a (projected) coordinate descent update. Because the subgradient of
the MCP penalty is identical to that of $\ell_{1}$-norm at zero, the free set
is defined in the same way as for the $\ell_{1}$-norm. As ${\bf{L}}^{(t)}$
approaches a solution $\bar{\bf{L}}$,
$\textbf{Free}\left({{\bf{L}}^{(t)}}\right)$ approaches
$\left\\{{(i,j):\bar{\bf{L}}_{i,j}=0}\right\\}$. Restricting (24) to the free
set variables improves the Hessian’s condition number and decreases the number
of iterations needed to solve (24).

So far, we discussed the outer iteration of the proposed method, which is
summarized in Algorithm 1.

Input: Empirical covariance matrix ${\bf{S}}$

Result: Estimated Laplacian matrix $\hat{{\bf{L}}}$

Initialization: Set ${\bf{L}}^{(0)}\in\mathcal{L}$

while  _Stopping criterion is not achieved_ do

  * •

Compute ${\bf{Q}}=\left({{\bf{L}}^{(t)}+{\bf{J}}}\right)^{-1}$.

* •

Define $\textbf{Free}\left({{\bf{L}}^{(t)}}\right)$ as in (26).

* •

Find the Newton direction $\delta^{(t)}$ by solving

(24) restricted to $\textbf{Free}\left({{\bf{L}}^{(t)}}\right)$. % by Alg. 2

* •

Update:
${\bf{L}}^{(t+1)}={\bf{L}}^{(t)}+\alpha\mathcal{P}\left({{\delta^{(t)}}}\right)$,
where $\alpha$

is achieved by linesearch.

end while

Algorithm 1 NewGLE

###  3.4 Solution of the Newton problem

To get the Newton direction for the $t$-th iteration, we approximately solve
the quadratic, linearly constrained and MCP-regularized problem in (24) using
an iterative method. We use the nonlinear conjugate gradient (NLCG) method and
modify it for our case. Specifically, we use the variant by Dai and Yuan,
(1999), and since (24) is linearly constrained we add projections between the
iterations. Moreover, to speed up the convergence, we use NLCG together with a
diagonal preconditioner, using the diagonal of the Hessian, as suggested in
Zibulevsky and Elad, (2010) for the LASSO problem.

The projected and preconditioned NLCG algorithm is the inner and most
computationally expensive part of the Algorithm 1. The algorithm itself and
its details are given in Appendix C. A key element of the approach is the
diagonal preconditioner ${\bf{D}}$, which is applied by replacing the gradient
$\left({\nabla_{\delta}f_{N}}\right)$ in the search direction by
${\bf{D}}^{-1}\cdot\nabla_{\delta}f_{N}$. The matrix ${\bf{D}}$ is defined as
the diagonal of the Hessian of $f_{N}(\delta)$, i.e.

${\bf{D}}_{k,k}=\left({\left[{{\bf{Q}}}\right]_{i,i}+\left[{{\bf{Q}}}\right]_{j,j}-2\left[{{\bf{Q}}}\right]_{i,j}}\right)^{2}+\tilde{\varepsilon}_{k}^{2},$
(27)

where the relation between $(i,j)$ and $k$ appears in Definition 1. For the
full derivation of the Hessian diagonal please refer to Appendix C. Note that
the operator ${\bf{D}}$ is diagonal, so calculating its inverse matrix is
trivial.

##  4 More Related Works

GMRF estimation with $\ell_{1}$ regularization. The $\ell_{1}$ regularized MLE
under the GMRF model has been extensively studied recently (see d’Aspremont et
al., (2008); Friedman et al., (2008); Rolfs et al., (2012); Hsieh et al.,
(2013, 2011); Treister et al., (2016); Shalom et al., (2022) and references
therein). Some of these methods include a proximal Newton approach. However,
these methods cannot be directly applied to the graph learning problem, as the
goal here is to learn a Laplacian precision matrix as opposed to a general
positive definite matrix.

The $\ell_{1}$ regularization bias and its alternatives. One drawback of the
mentioned approaches is that the $\ell_{1}$ penalty causes an estimation bias,
which is relevant to several variable selection problems, e.g., least squares
minimization and precision estimation under $\ell_{1}$ priors. To reduce the
$\ell_{1}$ bias, two-stage methods, nonconvex penalties, and hard-thresholding
approaches have been introduced for least squares variable selection Zhang,
2010a ; Breheny and Huang, (2011); Loh and Wainwright, (2013) and GMRF
estimation Chen et al., (2018); Lam and Fan, (2009); Shen et al., (2012);
Finder et al., (2022). Several of these methods use nonconvex penalties, such
as SCAD, MCP, and capped $\ell_{1}$-norm (Zhang, 2010b, ). Different
algorithms have been proposed to solve the resulting nonconvex problems, e.g.,
minimization-maximization (Hunter and Li,, 2005), local quadratic and linear
approximations (Fan and Li,, 2001; Zou and Li,, 2008).

LGMRF estimation. Recent works Egilmez et al., (2017); Kumar et al., (2019);
Liu et al., (2019); Zhao et al., (2019) investigate the $\ell_{1}$-norm
penalized MLE under the LGMRF model to estimate a graph Laplacian from data.
However, it has been shown recently that imposing an $\ell_{1}$-norm penalty
to the Laplacian-constrained MLE produces an unexpected behavior: the number
of nonzero graph weights grows as the regularization parameter increases Ying
et al., 2020a .

To alleviate this issue, nonconvex penalties, such as the MCP, have been
proposed. For example, the work of Koyakumaru et al., (2023) suggests adding
Tikhonov regularization to convexify the cost function and use the primal-dual
splitting method to solve the optimization problem. However, as the Tikhonov
regularization parameter increases, the bias of the estimation increases as
well. Another approach in Ying et al., 2020a ; Ying et al., 2020b ; Vieyra,
(2022) uses a nonconvex estimation method by solving a sequence of weighted
$\ell_{1}$-norm penalized subproblems, according to a majorization-
minimization framework. In each majorization step, the $\ell_{1}$-norm weights
are updated, and then in the minimization step, the $\ell_{1}$-penalized
function is minimized by the projected gradient descent (PGD) method. The
approach in Ying et al., (2021); Tugnait, (2021) requires solving two
subproblems: The first is an initial estimator of the precision matrix, e.g.,
the MLE. Then, the estimation is used to define the $\ell_{1}$ weights for the
second subproblem. Both stages are solved by PGD.

All the abovementioned methods are gradient-based methods that do not consider
Hessian information. On the other hand, the experience with the closely
related GMRF estimation problem suggests that proximal Newton approaches would
result in superior computational efficiency compared to those methods.

##  5 Experimental Results

In this section, we present experimental results for the evaluation of the
proposed method. We define the experimental setup in Subsection 5.1. Then, we
present the results in terms of graph learning performance and computational
efficiency in Subsections 5.2 and 5.3, respectively. Our MATLAB code for
replicating these experiments will be made publicly available on GitHub upon
the acceptance of this paper for publication.

###  5.1 Experimental setup

We compare the performance of the following methods:

  * •

MLE, defined in (4), where the objective function does not include any
sparsity-promoting penalty.

  * •

Nonconvex Graph Learning (NGL) with MCP Ying et al., 2020b , which solves a
sequence of weighted $\ell_{1}$-norm penalized subproblems via the
majorization-minimization framework.

  * •

Adaptive Laplacian constrained Precision matrix Estimation (ALPE) Ying et al.,
(2021), which is based on a weighted $\ell_{1}$-norm-regularized MLE. The
weights are precalculated using the MLE.

  * •

PGD for solving (15) using the Projected Gradient Descent method. We developed
this method to demonstrate the advantages of using the second-order
approximation.

  * •

The proposed method (NewGLE) for the MCP-regularized MLE, which is summarized
in Algorithm 1 and Algorithm 2 (in Appendix C). We used $\gamma=1.01$ and
$\varepsilon=0$, and to start the solution we applied a few PGD iterations,
which were taken into account in the results.

The regularization parameter, $\lambda$, is fine-tuned for each method for the
best performance. In addition, all algorithms use the convergence criterion
$\frac{\|{\bf{L}}^{(t+1)}-{\bf{L}}^{\textrm{(t)}}\|_{\textrm{F}}}{\|{\bf{L}}^{(t+1)}\|_{\textrm{F}}}\leq\epsilon$
with the tolerance $\epsilon=10^{-4}$. All the experiments were conducted on a
machine with 2 Intel Xeon E5-2660 2.0GHz processors with 28 cores and 512GB
RAM.

We use the relative error (RE) and F-score (FS) to evaluate the performance of
the algorithms, where
$\textrm{RE}=\frac{\|\bar{{\bf{L}}}-{\bf{L}}^{*}\|_{\textrm{F}}}{\|{\bf{L}}^{*}\|_{\textrm{F}}}$
in which $\bar{{\bf{L}}}$ and ${\bf{L}}^{*}$ denote the estimated and true
precision matrices, respectively, and
$\textrm{FS}=\frac{2\textrm{tp}}{2\textrm{tp}+\textrm{fp}+\textrm{fn}}$, where
the true-positive (tp), false-positive (fp), and false-negative (fn) detection
of graph edges are calculated by comparing the supports of $\bar{{\bf{L}}}$
and ${\bf{L}}^{*}$. The F-score takes values in $[0,1]$, where $1$ indicates
perfect recovery of the support.

We generated several datasets based on different graph-based models. To create
the datasets, we first created a graph, and then its associated Laplacian
matrix, ${\bf{L}}^{*}$, is used to generate independent data samples from the
Gaussian distribution
$\mathcal{N}\left({{\bf{0}},({\bf{L}}^{*})^{\dagger}}\right)$. The graph
connectivity is based on two options: 1) a random planar graph consisting of
$p=1,000$ nodes Treister et al., (2016), and 2) a random Barabasi-Albert graph
Zadorozhnyi and Yudin, (2012) of degree 2 with $p=100$ nodes. The edge weights
for both graph models are uniformly sampled from the range $[0.5,2]$. The
curves in Figs. 1 and 5 are the results of an average of 10 Monte Carlo
realizations, and in Fig. 2 are the result of an average of 50 realizations.

###  5.2 Results - graph learning performance

(a)

(b)

Figure 1: Performance comparisons under (a) RE, (b) F-score with different
sample size ratios $n/p$ for planar graphs with 1,000 nodes.

(a)

(b)

Figure 2: Performance comparisons in terms of (a) RE, and (b) F-score with
different sample size ratios $n/p$ for Barabasi-Albert graph of degree 2 with
100 nodes.

Figures 1 and 2 show the graph learning performance of different methods for
learning the topology of planar graphs with 1,000 nodes (Fig. 1) and of
Barabasi-Albert graphs of degree two with 100 nodes (Fig. 2). The performance
is presented in terms of RE and F-score versus the sample size ratio $n/p$,
where $n$ is the number of data samples used to calculate the sample
covariance matrix, ${\bf{S}}$, and $p$ is the number of nodes.

Figures 1(a) and 2(a) show that the REs of all the estimators decrease as the
sample size ratio, $n/p$, increases. Figures 1(b) and 2(b) show that the
F-score of all the methods increases as the sample size ratio $n/p$ increases.
Figures 1 and 2 show that the proposed NewGLE method outperforms all compared
methods (ALPE, NGL, PGD, and MLE) in both the RE and F-score senses. It should
be noted that in 1(b) the comparison with the MLE has been removed since it
has a significantly lower F-score than the other methods and in order to
highlight the differences between the other methods. As shown in Figs. 1(b)
and 2(b) for a small sample size ratio, the methods do not provide a perfect
F-score. This is because there may not be enough data samples to recover the
graph connectivity effectively. Still, the proposed NewGLE method consistently
outperforms all compared methods in the F-score sense, especially when the
sample size ratio is small. This outcome is significant in large-scale
problems since, in these cases, the sample size ratio is usually small.

###  5.3 Results - Computational Costs

(a)

(b)

Figure 3: Run-time performance comparisons under RE with a sample size ratio
(a) $n/p=0.5$ and (b) $n/p=15$ for random planar graphs with 1,000 nodes.

(a)

(b)

Figure 4: Run-time performance comparisons under F-score with a sample size
ratio (a) $n/p=0.5$ and (b) $n/p=15$ for random planar graphs with 1,000
nodes. Figure 5: Average convergence times over sample size ratio $n/p$ for
random planar graphs with 1,000 nodes.

Figures 3-5 compare the computational efficiency of the proposed NewGLE method
with the other methods in terms of run time. These experiments were conducted
on relatively large planar graphs of $p=1,000$ nodes so that the measuring
times less depend on minor implementation details. Figure 3 shows the
convergence of the compared methods in terms of RE for two sample size ratios:
(a) $n/p=0.5$, and (b) $n/p=15$. Notably, the proposed method converges
significantly faster than the other methods in both cases, while the
difference between the methods is more significant for the case when the
sample size ratio is smaller, $n/p=0.5$.

Figure 4 shows the F-score of the compared methods during two realizations:
(a) with sample size ratio $n/p=0.5$, and (b) with sample size ratio $n/p=15$.
Notably, the proposed method reached the highest F-score and got it
significantly faster than the run time taken for the other methods to achieve
their highest F-score values. In addition, it can be seen that in the case of
$n/p=15$, all the methods succeeded in estimating the connectivity of the
graph, and the NewGLE, PGD, and ALPE methods managed to estimate it perfectly
(i.e., F-score = 1). In addition, it can be seen that the advantage of the
proposed method is more notable for a smaller sample size ratio, i.e.,
$n/p=0.5$. Additionally, the results presented in Figs. 3 and 4 reveal that
while the RE and F-score performance of the ALPE method, as demonstrated in
Subsection (5.2), is commendable and comparable to the proposed method, the
implementation of the ALPE method necessitates a significant investment in
computational resources to effectively learn a graph, as measured by RE and
F-score metrics. The reason for this is that the ALPE method requires an
estimation of the weights, which are used for sparsity-promoting
regularization.

Finally, Fig. 5 presents the average convergence time for each method over the
sample size ratio $n/p$. It shows that the average convergence time of all
estimators decreases as the sample size ratio $n/p$ increases. In addition, it
emphasizes the significant advantage of the proposed method in terms of
computational efficiency, especially when there are few samples and the sample
size ratio is small (i.e., $n/p\leq 1$), which is typically the case in large-
scale problems.

##  6 Conclusion

In this paper, we developed a proximal Newton method for minimizing a
regularized maximum likelihood objective function for estimating a sparse
graph Laplacian precision matrix from data distributed according to an LGMRF
model. Our method is based on a second-order approximation of the smooth part
of the objective function, while keeping the regularization and constraints
intact. Furthermore, we used the nonconvex MCP function to promote the
sparsity of the resulting graph without introducing a significant bias. Beyond
the quadratic approximation itself, we used three algorithmic features to
treat the Newton problem at each iteration efficiently: 1) we used a
free/active set splitting, (2) we used a nonlinear conjugate gradient method,
and (3) we used diagonal preconditioning. This resulted in an accurate and
efficient estimation of the sparse weighted graph.

Numerical studies demonstrated the effectiveness of the proposed algorithm in
terms of accuracy and computational complexity. The proposed method
outperformed recent methods in terms of convergence time, while achieving
similar or better F-scores and relative errors, particularly for problems with
a small sample size ratio (i.e., $n/p<1$). The plausible performance of the
proposed method is particularly important in high-dimensional problems, where
computational efficiency is crucial, and sample scarcity is a common issue.

## References

  * Anis et al., (2016) Anis, A., Gadde, A., and Ortega, A. (2016).  Efficient sampling set selection for bandlimited graph signals using graph spectral proxies.  IEEE Transactions on Signal Processing, 64(14):3775–3789.
  * Armijo, (1966) Armijo, L. (1966).  Minimization of functions having Lipschitz continuous first partial derivatives.  Pacific Journal of mathematics, 16(1):1–3.
  * Banerjee et al., (2008) Banerjee, O., El Ghaoui, L., and d’Aspremont, A. (2008).  Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data.  The Journal of Machine Learning Research, 9:485–516.
  * Banerjee et al., (2006) Banerjee, O., Ghaoui, L. E., d’Aspremont, A., and Natsoulis, G. (2006).  Convex optimization techniques for fitting sparse gaussian graphical models.  In Proceedings of the 23rd international conference on Machine learning, pages 89–96.
  * Breheny and Huang, (2011) Breheny, P. and Huang, J. (2011).  Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection.  The annals of applied statistics, 5(1):232.
  * Cardoso and Palomar, (2020) Cardoso, J. V. d. M. and Palomar, D. P. (2020).  Learning undirected graphs in financial markets.  In Asilomar Conference on Signals, Systems, and Computers, pages 741–745.
  * Chen et al., (2018) Chen, J., Xu, P., Wang, L., Ma, J., and Gu, Q. (2018).  Covariate adjusted precision matrix estimation via nonconvex optimization.  In International Conference on Machine Learning, pages 922–931. PMLR.
  * Chin et al., (2013) Chin, G. M., Nocedal, J., Olsen, P. A., and Rennie, S. J. (2013).  Second order methods for optimizing convex matrix functions and sparse covariance clustering.  IEEE Transactions on Audio, Speech, and Language Processing, 21(11):2244–2254.
  * Chung, (1997) Chung, F. R. (1997).  Spectral graph theory, volume 92.  American Mathematical Soc.
  * Dabush et al., (2023) Dabush, L., Kroizer, A., and Routtenberg, T. (2023).  State estimation in unobservable power systems via graph signal processing tools.  Sensors, 23(3).
  * Dai and Yuan, (1999) Dai, Y.-H. and Yuan, Y. (1999).  A nonlinear conjugate gradient method with a strong global convergence property.  SIAM Journal on optimization, 10(1):177–182.
  * d’Aspremont et al., (2008) d’Aspremont, A., Banerjee, O., and El Ghaoui, L. (2008).  First-order methods for sparse covariance selection.  SIAM Journal on Matrix Analysis and Applications, 30(1):56–66.
  * Deka et al., (2018) Deka, D., Backhaus, S., and Chertkov, M. (2018).  Structure learning in power distribution networks.  IEEE Trans. Control of Network Systems, 5(3):1061–1074.
  * Dong et al., (2016) Dong, X., Thanou, D., Frossard, P., and Vandergheynst, P. (2016).  Learning Laplacian matrix in smooth graph signal representations.  IEEE Transactions on Signal Processing, 64(23):6160–6173.
  * Egilmez et al., (2017) Egilmez, H. E., Pavez, E., and Ortega, A. (2017).  Graph learning from data under Laplacian and structural constraints.  IEEE Journal of Selected Topics in Signal Processing, 11(6):825–841.
  * Fan and Li, (2001) Fan, J. and Li, R. (2001).  Variable selection via nonconcave penalized likelihood and its oracle properties.  Journal of the American statistical Association, 96(456):1348–1360.
  * Finder et al., (2022) Finder, S. E., Treister, E., and Freifeld, O. (2022).  Effective learning of a GMRF mixture model.  IEEE Access, 10:7289–7299.
  * Friedman et al., (2008) Friedman, J., Hastie, T., and Tibshirani, R. (2008).  Sparse inverse covariance estimation with the graphical lasso.  Biostatistics, 9(3):432–441.
  * Giudici and Spelta, (2016) Giudici, P. and Spelta, A. (2016).  Graphical network models for international financial flows.  Journal of Business & Economic Statistics, 34(1):128–138.
  * Halihal and Routtenberg, (2022) Halihal, M. and Routtenberg, T. (2022).  Estimation of the admittance matrix in power systems under Laplacian and physical constraints.  In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 5972–5976.
  * Holbrook, (2018) Holbrook, A. (2018).  Differentiating the pseudo determinant.  Linear Algebra and its Applications, 548:293–304.
  * Horn and Johnson, (2012) Horn, R. A. and Johnson, C. R. (2012).  Matrix analysis.  Cambridge university press.
  * Hsieh et al., (2011) Hsieh, C.-J., Dhillon, I., Ravikumar, P., and Sustik, M. (2011).  Sparse inverse covariance matrix estimation using quadratic approximation.  Neural information processing systems, 24.
  * Hsieh et al., (2014) Hsieh, C.-J., Sustik, M. A., Dhillon, I. S., and Ravikumar, P. (2014).  QUIC: quadratic approximation for sparse inverse covariance estimation.  Journal of Machine Learning Research, 15(1):2911–2947.
  * Hsieh et al., (2013) Hsieh, C.-J., Sustik, M. A., Dhillon, I. S., Ravikumar, P. K., and Poldrack, R. (2013).  BIG & QUIC: Sparse inverse covariance estimation for a million variables.  Advances in neural information processing systems, 26.
  * Hunter and Li, (2005) Hunter, D. R. and Li, R. (2005).  Variable selection using MM algorithms.  Annals of statistics, 33(4):1617.
  * Koller and Friedman, (2009) Koller, D. and Friedman, N. (2009).  Probabilistic graphical models: principles and techniques.  MIT press.
  * Koyakumaru et al., (2023) Koyakumaru, T., Yukawa, M., Pavez, E., and Ortega, A. (2023).  Learning sparse graph with minimax concave penalty under Gaussian Markov random fields.  IEICE Transaction on Fundamentals of Electronics, Communications and Computer Sciences, 106(1):23–34.
  * Kroizer et al., (2022) Kroizer, A., Routtenberg, T., and Eldar, Y. C. (2022).  Bayesian estimation of graph signals.  IEEE Transactions on Signal Processing, 70:2207–2223.
  * Kumar et al., (2019) Kumar, S., Ying, J., de Miranda Cardoso, J. V., and Palomar, D. (2019).  Structured graph learning via Laplacian spectral constraints.  Neural information processing systems, 32.
  * Lam and Fan, (2009) Lam, C. and Fan, J. (2009).  Sparsistency and rates of convergence in large covariance matrix estimation.  The Annals of statistics, 37(6B):4254–4278.
  * Lee et al., (2014) Lee, J. D., Sun, Y., and Saunders, M. A. (2014).  Proximal Newton-type methods for minimizing composite functions.  SIAM Journal on Optimization, 24(3):1420–1443.
  * Liu et al., (2019) Liu, T., Hoang, M. T., Yang, Y., and Pesavento, M. (2019).  A block coordinate descent algorithm for sparse gaussian graphical model inference with Laplacian constraints.  In International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), pages 236–240.
  * Loh and Wainwright, (2013) Loh, P.-L. and Wainwright, M. J. (2013).  Regularized M-estimators with nonconvexity: Statistical and algorithmic theory for local optima.  Advances in Neural Information Processing Systems, 26.
  * Mazumder and Hastie, (2012) Mazumder, R. and Hastie, T. (2012).  Exact covariance thresholding into connected components for large-scale graphical lasso.  The Journal of Machine Learning Research, 13(1):781–794.
  * Milanfar, (2012) Milanfar, P. (2012).  A tour of modern image filtering: New insights and methods, both practical and theoretical.  IEEE signal processing magazine, 30(1):106–128.
  * Narang and Ortega, (2012) Narang, S. K. and Ortega, A. (2012).  Perfect reconstruction two-channel wavelet filter banks for graph structured data.  IEEE Transactions on Signal Processing, 60(6):2786–2799.
  * Newman, (2010) Newman, M. (2010).  Networks: An Introduction.  Oxford University Press, Inc., New York, NY, USA.
  * Oztoprak et al., (2012) Oztoprak, F., Nocedal, J., Rennie, S., and Olsen, P. A. (2012).  Newton-like methods for sparse inverse covariance estimation.  Neural information processing systems, 25.
  * Rolfs et al., (2012) Rolfs, B., Rajaratnam, B., Guillot, D., Wong, I., and Maleki, A. (2012).  Iterative thresholding algorithm for sparse inverse covariance estimation.  Advances in Neural Information Processing Systems, 25.
  * Rue and Held, (2005) Rue, H. and Held, L. (2005).  Gaussian Markov random fields: theory and applications.  Chapman and Hall/CRC.
  * Shalom et al., (2022) Shalom, G., Treister, E., and Yavneh, I. (2022).  pista: preconditioned iterative soft thresholding algorithm for graphical lasso.  arXiv preprint arXiv:2205.10027.
  * Shen et al., (2012) Shen, X., Pan, W., and Zhu, Y. (2012).  Likelihood-based selection and sharp parameter estimation.  Journal of the American Statistical Association, 107(497):223–232.
  * Shuman et al., (2013) Shuman, D. I., Narang, S. K., Frossard, P., Ortega, A., and Vandergheynst, P. (2013).  The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains.  IEEE signal processing magazine, 30(3):83–98.
  * Treister and Turek, (2014) Treister, E. and Turek, J. S. (2014).  A block-coordinate descent approach for large-scale sparse inverse covariance estimation.  In Advances in neural information processing systems, pages 927–935.
  * Treister et al., (2016) Treister, E., Turek, J. S., and Yavneh, I. (2016).  A multilevel framework for sparse optimization with application to inverse covariance estimation and logistic regression.  SIAM Journal on Scientific Computing, 38(5):S566–S592.
  * Tugnait, (2021) Tugnait, J. K. (2021).  Sparse graph learning under Laplacian-related constraints.  IEEE Access, 9:151067–151079.
  * Vecchio et al., (2017) Vecchio, F., Miraglia, F., and Rossini, P. M. (2017).  Connectome: Graph theory application in functional brain network architecture.  Clinical neurophysiology practice, 2:206–213.
  * Vieyra, (2022) Vieyra, M. V. (2022).  Robust estimation of Laplacian constrained Gaussian graphical models with trimmed non-convex regularization.  In ICML 2022-Workshop on Principles of Distribution Shift.
  * Von Luxburg, (2007) Von Luxburg, U. (2007).  A tutorial on spectral clustering.  Statistics and computing, 17(4):395–416.
  * (51) Ying, J., Cardoso, J. V. d. M., and Palomar, D. P. (2020a).  Does the $\ell_{1}$-norm learn a sparse graph under Laplacian constrained graphical models?  arXiv preprint arXiv:2006.14925.
  * (52) Ying, J., de M. Cardoso, J. V., and Palomar, D. (2020b).  Nonconvex sparse graph learning under Laplacian constrained graphical model.  Advances in Neural Information Processing Systems, 33:7101–7113.
  * Ying et al., (2021) Ying, J., de Miranda Cardoso, J. V., and Palomar, D. (2021).  Minimax estimation of Laplacian constrained precision matrices.  In International Conference on Artificial Intelligence and Statistics, pages 3736–3744. PMLR.
  * Zadorozhnyi and Yudin, (2012) Zadorozhnyi, V. N. and Yudin, E. B. (2012).  Structural properties of the scale-free barabasi-albert graph.  Automation and Remote Control, 73:702–716.
  * (55) Zhang, C.-H. (2010a).  Nearly unbiased variable selection under minimax concave penalty.  The Annals of statistics, 38(2):894–942.
  * (56) Zhang, T. (2010b).  Analysis of multi-stage convex relaxation for sparse regularization.  Journal of Machine Learning Research, 11(3).
  * Zhao et al., (2019) Zhao, L., Wang, Y., Kumar, S., and Palomar, D. P. (2019).  Optimization algorithms for graph Laplacian estimation via ADMM and MM.  IEEE Transactions on Signal Processing, 67(16):4231–4244.
  * Zhu et al., (2003) Zhu, X., Lafferty, J., and Ghahramani, Z. (2003).  Combining active learning and semi-supervised learning using gaussian fields and harmonic functions.  In ICML 2003 workshop on the continuum from labeled to unlabeled data in machine learning and data mining, volume 3.
  * Zibulevsky and Elad, (2010) Zibulevsky, M. and Elad, M. (2010).  L1-l2 optimization in signal and image processing.  IEEE Signal Processing Magazine, 27(3):76–88.
  * Zou and Li, (2008) Zou, H. and Li, R. (2008).  One-step sparse estimates in nonconcave penalized likelihood models.  The Annals of statistics, 36(4):1509–1533.

APPENDIX

##  Appendix A PROOF OF THEOREM 1

In this section we prove Theorem 1 (which is also repeated below). We note
that this property was used in the works Egilmez et al., (2017); Ying et al.,
2020b  without explicit proof, and here we provide it.

Theorem 1 The optimization problem

$\displaystyle\underset{{\bf{L}}\in\mathcal{F}}{\text{minimize}}$
$\displaystyle{\rm
Tr}{\left({{\bf{L}}{\bf{S}}}\right)}-\log|{\bf{L}}+{\bf{J}}|+\rho\left({{\bf{L}};\lambda}\right),$

where the feasible set is given by

$\displaystyle\mathcal{F}=\left\\{{\bf{L}}\in\mathcal{S}^{p}\Big{|}\begin{array}[]{ll}{\bf{L}}_{i,j}\leq
0,\forall i\neq
j,~{}i,j=1,\ldots,p,\;\;{\bf{L}}{\bf{1}}={\bf{0}}\end{array}\right\\}.$ (29)

is equivalent to the optimization problem stated in (5).

###### Proof.

Since the objective functions in (5) and in (7) are identical (through (6)),
it remains to prove that the feasible sets (solution spaces) of the two
problems are equal under this objective function. It should be noted that
there are two distinctions between the feasible sets $\mathcal{L}$ and
$\mathcal{F}$ in problems (5) and (7), respectively. The first is that
$\mathcal{L}$ includes a positive semi-definite (PSD) constraint (i.e.,
${\bf{L}}\in\mathcal{S}^{p}_{+}$), and the second is that $\mathcal{L}$
includes a rank constraint (i.e., $\text{rank}\left({{\bf{L}}}\right)=p-1$).

To address the PSD constraint in problem (5), we observe that the constraints
${\bf{L}}{\bf{1}}={\bf{0}},$ and $\left[{{\bf{L}}}\right]_{i,j}\leq 0$ for
$i\neq j$ are present. Therefore, it follows that every matrix
${\bf{L}}\in\mathcal{F}$ is a diagonally dominant matrix with positive
diagonal entries. Hence, requiring symmetry is sufficient to ensure that the
minimizer of (7) is PSD as well (Horn and Johnson,, 2012).

Next, we show that for any graph Laplacian ${\bf{L}}\in\mathcal{F}$, adding
the constraint ${\textrm{rank}}\left({{\bf{L}}}\right)=p-1$ is either
equivalent to ${\textrm{rank}}\left({{\bf{L}}+{\bf{J}}}\right)=p$, or that the
matrix $\left({{\bf{L}}+{\bf{J}}}\right)$ is non-singular under the considered
objective function. To this end, we note that any ${\bf{L}}\in\mathcal{F}$
admits an eigenvalue decomposition ${\bf{U}}\Lambda_{L}{\bf{U}}^{T}$, where
${\bf{U}}$ is $p\times p$ matrix of the eigenvectors of ${\bf{L}}$, and
$\Lambda_{L}$ is a diagonal eigenvalue matrix as follows

$\Lambda_{L}=\textrm{diag}\left({\left[{0,\lambda_{2},\dots,\lambda_{p}}\right]^{T}}\right),$
(30)

where $0,\lambda_{2},\ldots,\lambda_{p}$ are the eigenvalues of ${\bf{L}}$.
The rank of ${\bf{L}}$ is the number of positive eigenvalues. In addition, the
first eigenvector of ${\bf{L}}$ is $\frac{1}{\sqrt{p}}{\mathbf{1}}$ and it
corresponds to the eigenvalue $\lambda_{1}=0$. On the other hand, an
eigendecomposition of ${\bf{J}}$ can be ${\bf{U}}\Lambda_{J}{\bf{U}}^{T}$,
where

$\Lambda_{J}=\textrm{diag}\left({\left[{1,0,\dots,0}\right]^{T}}\right)$ (31)

and the eigenvector of ${\bf{J}}$ that is associated with the single nonzero
eigenvalue is $\frac{1}{\sqrt{p}}{\mathbf{1}}$—the same one as in ${\bf{L}}$.
Hence, for the same eigenvector matrix ${\bf{U}}$, the eigenvalues of the
matrix ${\bf{L}}+{\bf{J}}$ are $1,\lambda_{2},\ldots,\lambda_{p}$, and
${\textrm{rank}}\left({{\bf{L}}+{\bf{J}}}\right)={\textrm{rank}}\left({{\bf{L}}}\right)+1$.
Hence, to conclude, if ${\bf{L}}\in\mathcal{F}$ and
$\left({{\bf{L}}+{\bf{J}}}\right)$ is a non-singular matrix, then
${\textrm{rank}}\left({{\bf{L}}}\right)=p-1$.

Now we come to the last part of the proof. Assume by contradiction that the
minimizer of (7) satisfies that the matrix $\left({{\bf{L}}+{\bf{J}}}\right)$
is singular, then the expression $-\log{\left({|{\bf{L}}+{\bf{J}}|}\right)}$
in the objective approaches infinity. Thus, any other matrix
${\bf{L}}\in\mathcal{F}$ for which the matrix
$\left({{\bf{L}}+{\bf{J}}}\right)$ is a non-singular matrix will yield a lower
(finite) objective value in (7). Hence it can be concluded that the solution
of (7) satisfies that the $\left({{\bf{L}}+{\bf{J}}}\right)$ matrix is non-
singular, i.e., ${\textrm{rank}}\left({{\bf{L}}+{\bf{J}}}\right)=p$ and hence
${\textrm{rank}}\left({{\bf{L}}}\right)=p-1$. ∎

##  Appendix B DERIVATION OF THE PARAMETERIZED NEWTON PROBLEM

First, we simplify the writing constraints of (18). We assume that the
previous iterate ${\bf{L}}^{(t)}$ satisfies the constraints, and by definition
${\bf{L}}^{(t+1)}$ will satisfy them as well. To achieve that we also choose
${\bf{L}}^{(0)}\in\mathcal{F}$. Based on the definition of $\mathcal{F}$ in
(10), under the assumption that ${\bf{L}}^{(t)}$ satisfies the constraints, in
order to obtain an update ${\bf{L}}^{(t+1)}$ that satisfy these constraints,
we require the following conditions on the variable $\Delta$:

  * •

${\bf{L}}^{(t)}+\Delta\in\mathcal{S}^{p}$. Hence, we require that
$\Delta\in\mathcal{S}^{p}$.

  * •

The off-diagonal elements should satisfy
$\left[{{\bf{L}}^{(t)}+\Delta}\right]_{i,j}\leq 0$, $\forall i\neq j$. Hence,
we require that $\Delta_{i,j}\leq-{\bf{L}}^{(t)}_{i,j}$.

  * •

$\left({{\bf{L}}^{(t)}+\Delta}\right){\bf{1}}={\bf{0}}$. Assuming that
${\bf{L}}^{(t)}{\bf{1}}={\bf{0}}$, then we only need to require
$\Delta{\bf{1}}={\bf{0}}$.

By substituting these conditions in (18) and removing the constant term in the
objective, we can rewrite the Newton problem from (18) as follows:

$\displaystyle\underset{\Delta\in\mathcal{S}^{p}}{\text{minimize}}$
$\displaystyle{\rm
Tr}{\left({\Delta\left({{\bf{S}}-{\bf{Q}}}\right)}\right)}+\frac{1}{2}{\rm
Tr}{\left({\Delta{\bf{Q}}\Delta{\bf{Q}}}\right)}+\frac{1}{2}\|\varepsilon\odot\Delta\|_{F}^{2}+\rho_{MCP}\left({{\bf{L}}^{(t)}+\Delta;\lambda}\right)$
(32) s.t. $\displaystyle\Delta_{i,j}\leq-{\bf{L}}^{(t)}_{i,j},\forall i\neq j$
$\displaystyle\Delta{\bf{1}}={\bf{0}}.$

Next, we further simplify the problem by replacing the symmetric unknown
matrix $\Delta\in\mathcal{S}^{p}$ with a vector
$\delta\in\mathbb{R}^{p(p-1)/2}$ by using Definition 1 of the linear operator
$\mathcal{P}:\mathbb{R}^{p(p-1)/2}$. According to this definition, given any
${\bf{w}}\geq 0$, $\mathcal{P}\left({{{\bf{w}}}}\right)$ satisfies
$\sum_{j}\mathcal{P}\left({{{\bf{w}}}}\right)_{i,j}=0$ for any
$i\in\left\\{{1,\dots,p}\right\\}$, and
$\mathcal{P}\left({{{\bf{w}}}}\right)_{i,j}=\mathcal{P}\left({{{\bf{w}}}}\right)_{j,i}$
for any $i\neq j$. In other words, $\mathcal{P}\left({{{\bf{w}}}}\right)$ is a
Laplacian matrix for any ${\bf{w}}$ which has nonnegative values. The simple
example below illustrates the definition of $\mathcal{P}$. Let
${\bf{w}}=[w_{1},w_{2},w_{3}]\in\mathbb{R}^{3}$ be a nonnegative vector. Then

$\displaystyle\mathcal{P}\left({{{\bf{w}}}}\right)=\begin{bmatrix}\sum_{i=1,2}w_{i}&-w_{1}&-w_{2}\\\
-w_{1}&\sum_{i=1,3}w_{i}&-w_{3}\\\ -w_{2}&-w_{3}&\sum_{i=2,3}w_{i}\\\
\end{bmatrix}.$ (33)

Since the Laplacian is a symmetric matrix, we can rewrite the sparse promoting
term in (32) as:

$\displaystyle\rho_{MCP}\left({\mathcal{P}\left({{{\bf{w}}}}\right);\lambda}\right)$
$\displaystyle\overset{(a)}{=}\sum_{i\neq
j}MCP\left({{\bf{L}}_{ij};\gamma,\lambda}\right)$ (34)
$\displaystyle\overset{(b)}{=}2\sum_{i>j}MCP\left({{\bf{L}}_{ij};\gamma,\lambda}\right)$
$\displaystyle\overset{(c)}{=}2\sum_{k=1}^{p\left({p-1}\right)/2}MCP\left({-w_{k};\gamma,\lambda}\right)$
$\displaystyle\overset{(d)}{=}2\tilde{\rho}_{mcp}\left({{{\bf{w}}};\lambda}\right),$

where ${\bf{L}}=\mathcal{P}\left({{{\bf{w}}}}\right)$. In (34), $(a)$ is
obtained by substituting (11), $(b)$ is obtained from the symmetry of
${\bf{L}}$ (i.e., since ${\bf{L}}\in\mathcal{S}^{p}$), and $(c)$ is obtained
by substituting (23). Finally, $(d)$ is obtained since
$\tilde{\rho}_{MCP}\left({{{\bf{w}}};\lambda}\right)$ defined in (25) is an
even function. Noting that for ${\bf{w}},\delta\in\mathbb{R}^{p(p-1)/2}$, such
that ${\bf{L}}=\mathcal{P}\left({{{\bf{w}}}}\right)$ and
$\Delta=\mathcal{P}\left({{\delta}}\right)$, we can write:

  * •

$\rho_{mcp}\left({{\bf{L}}+\Delta;\lambda}\right)\overset{\eqref{r mcp for
vec}}{=}2\tilde{\rho}_{mcp}\left({{\bf{w}}+\delta;\lambda}\right)$.

  * •

$\Delta\in\left\\{{\Delta\in\mathcal{S}^{p}|\Delta{\bf{1}}={\bf{0}},\Delta_{i,j}\leq-{\bf{L}}^{(k)}_{i,j},\forall
i\neq j}\right\\}$ $\overset{\eqref{Lx
matrix}}{\iff}\delta\in\left\\{{\delta\in\mathbb{R}^{p(p-1)/2}|\quad\delta\geq-{\bf{w}}}\right\\}$.

  * •

$\frac{1}{2}\|\varepsilon\odot\Delta\|_{F}^{2}=\|\tilde{\varepsilon}\odot\delta\|_{2}^{2}$,
which holds trivially because the diagonal of $\varepsilon$ is 0 and it is
symmetric.

Hence, (32), which is equivalent to (18), is also equivalent to (24), which is
the problem we solve in practice in the inner Newton problem.

##  Appendix C NONLINEAR PRECONDITIONED CONJUGATE GRADIENTS

Generally, iterates of projection methods are obtained by

$\delta^{(t+1)}=\Pi_{U}\left\\{{\delta^{(t)}+\zeta{\bf{d}}^{(t)}}\right\\},$
(35)

where $\delta^{(t)}$ is the current iterate, ${\bf{d}}^{(t)}$ is the descent
(search) direction, $\zeta>0$ is a step-size parameter, $U$ is some feasible
set, and $\Pi_{U}\left\\{{{\bf{v}}}\right\\}$ is a projection of ${\bf{v}}$
onto $U$. That is, using the $\ell_{2}$ norm we have
$\Pi\\{{\bf{v}}\\}=\arg\min_{\bf{x}}\|{\bf{x}}-{\bf{v}}\|_{2}$. Since there is
only a nonnegativity constraint in (24), the projection here is simple to
calculate and is given by:

$\Pi_{{\bf{v}}+{\bf{w}}\geq{\bf{0}}}\left\\{{{\bf{v}}}\right\\}=\textrm{max}\left\\{{{\bf{v}},-{\bf{w}}}\right\\}.$
(36)

The step-size $\zeta>0$ is computed by a linesearch procedure, and the descent
direction ${\bf{d}}^{(t)}$ is given by

$\displaystyle{\bf{d}}^{(t)}=\left\\{\begin{array}[]{lr}-\nabla_{\delta}f_{N}\left({\delta^{(t)}}\right)&\textrm{if}\quad
t=0\\\
-\nabla_{\delta}f_{N}\left({\delta^{(t)}}\right)+\beta_{t}{\bf{d}}^{(t-1)}&\textrm{if}\quad
t>0\end{array},\right.$ (39)

where $\beta_{t}$ is a scalar and ${f_{N}(\delta^{(t)})}$ is the objective
defined in (24) at the point $\delta^{(t)}$. Since 1952, there have been many
formulas for the scalar $\beta_{t}$, and in this paper, we use the DY method
Dai and Yuan, (1999), given by

$\beta^{DY}_{t}=\textrm{max}\left\\{{0,\frac{\|\nabla_{\delta}f^{(t)}_{N}\|^{2}_{2}}{\left({\nabla_{\delta}f^{(t)}_{N}-\nabla_{\delta}f^{(t-1)}_{N}}\right)^{T}{\bf{d}}^{(t-1)}}}\right\\}.$
(40)

One remarkable property of the DY method is that it produces a descent
direction at each iteration and converges globally for convex problems in the
case where linesearch is used Dai and Yuan, (1999).

To apply the projected NLCG, the gradient of $f_{N}(\delta)$ (defined in (24))
is needed and is given by

$\displaystyle\begin{aligned}
\nabla_{\delta}f_{N}=\mathcal{P}^{*}\left({\nabla_{\Delta}f_{N}}\right)=\quad&\mathcal{P}^{*}\left({{\bf{S}}-{\bf{Q}}+{\bf{Q}}\mathcal{P}\left({{\delta}}\right){\bf{Q}}}\right)+2\tilde{\varepsilon}^{2}\odot\delta+2\nabla_{\delta}\tilde{\rho}_{mcp}\left({{\bf{w}}+\delta;\lambda}\right),\end{aligned}$
(41)

where $\mathcal{P}\left({{\delta}}\right)=\Delta$, and $\mathcal{P}^{*}$ is
the adjoint of $\mathcal{P}$, that is

$\left[{\mathcal{P}^{*}\left({{\bf{Y}}}\right)}\right]_{k}\stackrel{{\scriptstyle\triangle}}{{=}}Y_{ii}+Y_{jj}-Y_{ij}-Y_{ji}$
(42)

for $k$ as in Definition 1. Here we see the advantage of the proximal Newton
approach. In (41) the gradient of $\tilde{\rho}$ is computed at
${\bf{w}}+\delta$, instead of ${\bf{w}}$ alone if we were to naively majorize
$\tilde{\rho}$ in the outer Newton problem. So, the inner solver is able to
adapt to the curvature of $\tilde{\rho}$ and to the different sections in its
definition.

###  C.1 Diagonal Preconditioning

The NLCG method that we use in this work is essentially an accelerated
gradient-based method. To further accelerate it with rather minimal
computational effort, a diagonal approximation to the Hessian can be used,
acting as a preconditioner in the NLCG method. The aim of the diagonal
preconditioner is to properly scale the gradient and speed up the rate of
convergence of the iterative inner method for the Newton direction. In this
paper, we suggest using the diagonal of the Hessian with respect to $\delta$
as a preconditioner with NLCG. Below we show the derivation of this
preconditioner.

Our diagonal preconditioner matrix, ${\bf{D}}$, is defined as the diagonal of
the Hessian of $f_{N}$. To ease the derivation, we ignore the $\varepsilon$
term of (24) in this section. Following the Hessian of $f$ at (17) and
considering the linear operator $\mathcal{P}$ as a matrix in
$\mathbb{R}^{p^{2}\times p(p-1)/2}$ (i.e., ignoring the vector-to-matrix
reshaping in its definition), the Hessian of $f_{N}$ is

${\bf{H}}_{N}=\mathcal{P}^{\top}\left({\bf{Q}}\otimes{\bf{Q}}\right)\mathcal{P}\in\mathbb{R}^{p(p-1)/2\times
p(p-1)/2}.$ (43)

To get its $(k,k)$ diagonal entry for our preconditioner, we multiply
${\bf{H}}_{N}$ by a unit vector from both sides:

$[{\bf{D}}]_{k,k}=[{\bf{H}}_{N}]_{k,k}=\mathbf{e}_{k}^{T}{\bf{H}}_{N}\mathbf{e}_{k},$
(44)

where $\mathbf{e}_{k}\in\mathbb{R}^{p(p-1)/2}$ is the unit vector of zeros,
except the $k$-th entry which equals 1. Equivalently to (44), we can write

$[{\bf{H}}_{N}]_{k,k}={\rm
Tr}{\left({\mathcal{P}(\mathbf{e}_{k}){\bf{Q}}\mathcal{P}(\mathbf{e}_{k}){\bf{Q}}}\right)},$
(45)

where the matrix $\mathcal{P}(\mathbf{e}_{k})\in\mathbb{R}^{p\times p}$ is a
matrix that is essentially the reshape of the $k$-th column of the operator
$\mathcal{P}$, if we consider it as a matrix. If we look at the example of
(33), then taking $k=1$ corresponds to

$\mathcal{P}(\mathbf{e}_{1})=\begin{bmatrix}\;\;1&-1&0\\\ -1&\;\;1&0\\\
\;\;0&\;\;0&0\\\ \end{bmatrix}.$

Note that $k=1$ in (33) influences only the $2\times 2$ sub-matrix of the
indices $(1,2)$ in $\mathcal{P}(\mathbf{e}_{1})$. More generally,
$\mathcal{P}(\mathbf{e}_{k})$ results in a matrix of zeros, except a $2\times
2$ nonzero submatrix of the indices $(i,j)$ where $k$ corresponds to the pair
$(i,j)$ according to Definition 1. Hence, multiplying
$\mathcal{P}(\mathbf{e}_{k}){\bf{Q}}$, will result in a matrix of zeros except
the $i$-th and $j$-th _rows_ , which equal to $\textbf{q}_{i}-\textbf{q}_{j}$
and $\textbf{q}_{j}-\textbf{q}_{i}$, respectively, where $\textbf{q}_{i}$ is
the $i$-th column/row of (the symmetric) ${\bf{Q}}$. According to (45) we now
have

$\displaystyle[{\bf{H}}_{N}]_{k,k}$ $\displaystyle=$ $\displaystyle{\rm
Tr}{\left({(\mathcal{P}(\mathbf{e}_{k}){\bf{Q}})^{2}}\right)}$ (46)
$\displaystyle=$ $\displaystyle{\rm
Tr}{\left({\begin{pmatrix}[{\bf{Q}}]_{i,i}-[{\bf{Q}}]_{j,i}&[{\bf{Q}}]_{i,j}-[{\bf{Q}}]_{j,j}\\\
[{\bf{Q}}]_{j,i}-[{\bf{Q}}]_{i,i}&[{\bf{Q}}]_{j,j}-[{\bf{Q}}]_{j,i}\end{pmatrix}^{2}}\right)}$
$\displaystyle=$
$\displaystyle\left({\left[{{\bf{Q}}}\right]_{i,i}+\left[{{\bf{Q}}}\right]_{j,j}-2\left[{{\bf{Q}}}\right]_{i,j}}\right)^{2}=[{\bf{D}}]_{k,k},$

where $\left({i,j}\right)$ corresponds to the index $k$ according to
Definition 1. In the last equality we used the symmetry of ${\bf{Q}}$, i.e.,
$[{\bf{Q}}]_{j,i}=[{\bf{Q}}]_{i,j}$. Below, Algorithm 2 summarizes the
projected and preconditioned NLCG method for solving (18). Note that the
linesearch in the last bullet of the algorithm is applied before the
projection, hence taking a smaller steplength $\xi$ means we effectively
enlarge the Hessian approximation at each step.

Input:
${\bf{L}}^{(t)}=\mathcal{P}\left({{{\bf{w}}}}\right),{\bf{Q}},\text{Free}\left({{\bf{L}}^{(t)}}\right),\lambda$

Result: Newton direction $\Delta^{(t)}$

Initialization: Set $\delta^{(0)}={\bf{0}}$

while  _Stopping criterion is not achieved_ do

  * •

Compute the preconditioner ${\bf{D}}$ as in (46).

* •

Compute the gradient $\nabla_{\delta}f_{N}$ as in (41).

* •

Apply the preconditioner
${\bf{g}}=\left({{\bf{D}}^{-1}\cdot\nabla_{\delta}f_{N}}\right)$. • Zero non-
free elements: $g_{k}\leftarrow 0\;\forall
k\in\text{Free}\left({{\bf{L}}^{(t)}}\right)$. • Compute ${\bf{d}}$ as in
(39), with ${\bf{g}}$ in the role of $\nabla_{\delta}f_{N}$.

* •

Update:
$\delta^{(t+1)}=\textbf{max}\left({\delta^{(t)}+\zeta{\bf{d}},-{\bf{w}}}\right)$,
where $\zeta$ achieved by linesearch.

end while

Algorithm 2 Projected & preconditioned NLCG

##  Appendix D CONVERGENCE GUARANTEES: PROOF OF THEOREM 2

In this section, we provide theoretical proof for the convergence of our
method to a stationary point (Theorem 2). We note that the MCP penalty is non-
convex; hence one typically gets to a local minimum using a monotone series of
iterations. Some results in our convergence proof follow Hsieh et al., (2014)
(e.g., Lemma 1 and Lemma 3). However, in the work (Hsieh et al.,, 2014), the
$\ell_{1}$ norm is used as a penalty, and the proof heavily relies on its
convexity. Here, we had to adapt the proofs for these Lemmas to the concave
MCP penalty, which is far from trivial and has not been done before, to the
best of our knowledge.

This section is organized as follows. In subsection D.1 we define the
considered setup of proximal methods for MCP regularized penalties. In
subsection D.2 we provide some spectral bounds on our problem, which are
needed for the proofs that follow. In subsection D.3 we provide a lemma
proving that our method is a fixed point iteration. In subsection D.4 we prove
that the objective value sufficiently decreases between the iterations, which
is the main point needed for the convergence proof. Lastly, in subsection D.5
we finalize the convergence proof given the previous auxiliary lemmas.

###  D.1 The Proximal Methods Setup

Let us write the problem (15) as follows

$\displaystyle\underset{{\bf{L}}\in\mathcal{F}}{\text{minimize
}}F({\bf{L}})=f({\bf{L}})+\rho_{MCP}\left({{\bf{L}};\lambda}\right),$ (47)

where $f$ is the smooth part of the objective, defined in (16), $\mathcal{F}$
and $\rho_{MCP}\left({{\bf{L}};\lambda}\right)$ are defined in (10) and (11),
respectively. In the setup of general proximal methods, at each iteration $t$,
we minimize the penalized and constrained quadratic objective, as follows:

$\displaystyle\Delta^{(t)}=\underset{\Delta\in\mathcal{S}^{p}}{\text{argmin}}$
$\displaystyle G(\Delta)=f({\bf{L}}^{(t)})+\big{\langle}\nabla
f({\bf{L}}^{(t)}),\Delta\rangle+\frac{1}{2}\big{\langle}\Delta,\Delta\rangle_{{\bf{M}}^{(t)}}+\rho_{MCP}\left({{\bf{L}}^{(t)}+\Delta;\lambda}\right)$
(48) s.t. $\displaystyle({\bf{L}}^{(t)}+\Delta)\in\mathcal{F},$

where ${\bf{M}}^{(t)}$ is some positive definite matrix that is an
approximation of the true Hessian, or the true Hessian itself. In particular,
Proximal Gradient Descent uses ${\bf{M}}={\bf{I}}$, while proximal Newton will
have ${\bf{M}}=\nabla^{2}f({\bf{L}})+\mbox{diag}(\mbox{vec}(\varepsilon))$, as
in (18) and Algorithm 1. Once the quadratic problem (48) is solved, we apply
linesearch as in (19), finding a step size $\alpha^{(t)}$ that ensures a
reduction in the objective $F$. We now prove the convergence of this framework
to our Laplacian estimation problem.

#### Optimality Conditions

We can state that by the definition of $\mathcal{F}$ our optimality condition
is given by

$\left\\{\begin{array}[]{ll}0\in\frac{\partial
F({\bf{L}}^{*})}{\partial{\bf{L}}^{*}_{i,j}}+\mu_{i,j},&i>j\\\
{\bf{L}}^{*}_{j,i}={\bf{L}}^{*}_{i,j}\geq 0,&i\neq j\\\
{\bf{L}}^{*}_{i,i}=\sum_{j\neq
i}-{\bf{L}}^{*}_{i,j}&i=1,...,p\end{array}\right.,$ (49)

where $\frac{\partial F({\bf{L}}^{*})}{\partial{\bf{L}}_{i,j}}=[\nabla
f]_{i,j}+\frac{\partial
MCP({\bf{L}}^{*}_{ij};\gamma,\lambda)}{\partial{\bf{L}}_{ij}}$ is the sub-
differential of $F$ w.r.t ${\bf{L}}_{i,j}$, $\mu_{i,j}\geq 0$ is an inequality
Lagrange multiplier for the constraint ${\bf{L}}_{i,j}\leq 0$, and the sub-
differential of the MCP function in (14) is given by

$\frac{\partial MCP(x;\gamma,\lambda)}{\partial
x}=\left\\{\begin{array}[]{lc}\in[-\lambda,\lambda]&x=0\\\
\lambda\cdot\text{sign}\left({x}\right)-\frac{x}{\gamma}&0<|x|\leq\gamma\lambda\\\
0&|x|>\gamma\lambda\end{array}\right..$ (50)

If the conditions in (49) are fulfilled, then ${\bf{L}}^{*}$ is a stationary
point of (15). Since the MCP penalty is non-convex, we expect (15) to have
multiple minima, and we show that Algorithm 1 converges to one of them. We
initialize our algorithm with ${\bf{L}}^{(0)}\in\mathcal{F}$ and keep
${\bf{L}}^{(t)}\in\mathcal{F}$ in all iterations.

###  D.2 Spectral Bounds

Let us first note that the $-\log\det({\bf{L}}+{\bf{J}})$ term in (15) acts as
a barrier function. That is, if ${\bf{L}}+{\bf{J}}$ gets closer to being
singular, then $\det({\bf{L}}+{\bf{J}})\to 0$ and therefore,
$F({\bf{L}})\to\infty$ following the $\log$ term. Therefore, we can say that
as long as $F({\bf{L}}^{(t)})$ is bounded, there exists some $\beta_{1}>0$ for
which ${\bf{L}}^{(t)}+{\bf{J}}\succeq\beta_{1}{\bf{I}}$ for all $t$.
Furthermore, by Ying et al., (2021), we have that for every Laplacian matrix
${\bf{L}}$

$\langle{\bf{S}},{\bf{L}}\rangle=\langle\mathbf{v},\mathbf{w}\rangle,$ (51)

where ${\bf{L}}=\mathcal{P}\mathbf{w},\;\;\mathbf{w}\geq 0$, and

$\mathbf{v}_{k}=[\mathcal{P}^{*}{\bf{S}}]_{k}=\left[\mathcal{P}^{*}\frac{1}{n}\sum_{q=1}^{n}\mathbf{x}_{q}\mathbf{x}_{q}^{T}\right]_{k}=\frac{1}{n}\sum_{q=1}^{n}\left([\mathbf{x}_{q}]_{i}-[\mathbf{x}_{q}]_{i}\right)^{2}>0$
(52)

holds with a very high probability, and $\mathcal{P}^{*}$ is the adjoint of
$\mathcal{P}$ in (42). Following that, we have

$\langle\mathbf{v},\mathbf{w}\rangle\geq\min_{j}\\{\mathbf{v}_{j}\\}\sum_{q}{\mathbf{w}_{q}}=\frac{1}{2}\min_{j}\\{\mathbf{v}_{j}\\}{\rm
Tr}{\left({{\bf{L}}}\right)}.$ (53)

Now, for any ${\bf{L}}\in\mathcal{F}$ such that $F({\bf{L}})\leq
F({\bf{L}}^{(0)})$ we have

${\rm
Tr}{\left({{\bf{L}}{\bf{S}}}\right)}-\log|{\bf{L}}|_{+}+\rho\left({{\bf{L}};\lambda}\right)\leq
F({\bf{L}}^{(0)}),$

and since $\rho\left({{\bf{L}};\lambda}\right)\geq 0$ we can write:

${\rm Tr}{\left({{\bf{L}}{\bf{S}}}\right)}\leq
F({\bf{L}}^{(0)})+\log|{\bf{L}}|_{+}.$

Furthermore, using (53) and
$\log|{\bf{L}}|_{+}\leq(p-1)\log(\lambda_{max}({\bf{L}}))$, we get

$\min_{k}\\{\mathbf{v}_{k}\\}\frac{1}{2}\lambda_{max}({\bf{L}})\leq\min_{k}\\{\mathbf{v}_{k}\\}\frac{1}{2}{\rm
Tr}{\left({{\bf{L}}}\right)}\leq
F({\bf{L}}^{(0)})+(p-1)\log(\lambda_{max}({\bf{L}})).$ (54)

Now, since in (54) the left-had-side grows much faster than the right-hand-
side, then there is and upper-bound on $\lambda_{max}({\bf{L}})$, depending on
${\bf{S}}$ and ${\bf{L}}^{(0)}$, assuming (52). Hence, we can conclude that if
the iterations $\\{{\bf{L}}^{(t)}\\}$ are monotonically decreasing starting
from ${\bf{L}}^{(0)}\in\mathcal{F}$, there exist upper and lower bounds

$\beta_{1}{\bf{I}}\preceq{\bf{L}}^{(t)}+{\bf{J}}\preceq\beta_{2}{\bf{I}}.$
(55)

This also means that the level set
$\mathcal{R}=\\{{\bf{L}}\in\mathcal{F}:F({\bf{L}})\leq F({\bf{L}}^{(0)})\\}$
is compact, and the Hessian is bounded, i.e.,
$\nabla^{2}f({\bf{L}})=({\bf{L}}+{\bf{J}})^{-1}\otimes({\bf{L}}+{\bf{J}})^{-1}\preceq\frac{1}{\beta_{1}^{2}}{\bf{I}}=\theta{\bf{I}}$,
where $\theta>0$. The last upper bound comes from the properties of
eigenvalues of Kronecker products.

We summarize the conclusions above in the following lemma:

######  Lemma 1.

Assume that the iterations ${\bf{L}}^{(t)}$ are monotonically deceasing with
respect to the objective in (47), starting from a feasible ${\bf{L}}^{(0)}$.
Then there exist constants $\beta_{1},\beta_{2}$ such that
$\beta_{1}{\bf{I}}\preceq{\bf{L}}^{(t)}+{\bf{J}}\preceq\beta_{2}{\bf{I}}$, and
the level-set $\mathcal{R}=\\{{\bf{L}}\in\mathcal{F}:F({\bf{L}})\leq
F({\bf{L}}^{(0)})\\}$ is compact. Also, the Hessian is bounded:
$\nabla^{2}f({\bf{L}})\preceq\frac{1}{\beta_{1}^{2}}{\bf{I}}=\theta{\bf{I}}$

###  D.3 Fixed Point Iteration

The following lemma shows that if ${\bf{L}}^{(t)}$ is feasible, and the
quadratic minimization (48) ends with $\Delta^{(t)}=0$, then ${\bf{L}}^{(t)}$
is a minimum point of (47).

######  Lemma 2.

Assume that ${\bf{L}}^{(t)}\in\mathcal{F}$. If $\Delta^{(t)}=0$ is a minimizer
of $G(\Delta)$ in (48), then ${\bf{L}}^{(t+1)}={\bf{L}}^{(t)}$ is a minimizer
of $F({\bf{L}})$ in (47).

###### Proof.

Since $\Delta^{(t)}$ is a minimizer of (48), the following optimality
conditions are held

$\left\\{\begin{array}[]{ll}0\in\frac{\partial
G(\Delta^{(t)})}{\partial\Delta_{i,j}}+\mu_{i,j},&i>j\\\
{\bf{L}}^{(t)}_{j,i}+\Delta^{(t)}_{j,i}={\bf{L}}^{(t)}_{i,j}+\Delta^{(t)}_{i,j}\geq
0,&i\neq j\\\ {\bf{L}}^{(t)}_{i,i}+\Delta^{(t)}_{i,i}=\sum_{j\neq
i}-({\bf{L}}^{(t)}_{i,j}+\Delta^{(t)}_{i,j})&i=1,...,p\end{array}\right.,$
(56)

where $\frac{\partial G(\Delta)}{\partial\Delta_{i,j}}=[\nabla
f]_{i,j}+[\langle{\bf{M}}\Delta\rangle]_{i,j}+\frac{\partial
MCP({\bf{L}}^{(t)}_{i,j}+\Delta_{i,j};\gamma,\lambda)}{\partial\Delta_{ij}}$,
and $\langle{\bf{M}}\Delta\rangle$ denotes a matrix vector multiplication of
${\bf{M}}$ with the vectorized $\Delta$, i.e.,
$\mbox{mat}({\bf{M}}\cdot\mbox{vec}(\Delta))$. Placing $\Delta^{(t)}=0$ in
(56) yields the same conditions as in (49) regardless of the specific choice
of the matrix ${\bf{M}}\succ 0$, hence ${\bf{L}}^{(t)}$ is a stationary point
of (47). ∎

###  D.4 Sufficient Decrease in Objective

We now show that we have a decrease in the objective $F$ following the
minimization in (48) and the linesearch in (19). This lemma is the key to our
convergence proof.

######  Lemma 3.

Assume that the Hessian is bounded, $\nabla^{2}{f}\preceq\theta{\bf{I}}$ for
$\theta>0$. Also assume that the iterates $\\{{\bf{L}}^{(t)}\\}$ are defined
by (48) followed by a linesearch (19), starting from
${\bf{L}}^{(0)}\in\mathcal{F}$, with ${\bf{M}}^{(t)}$ satisfying
${\bf{M}}^{(t)}\succeq\lambda^{{\bf{M}}}_{min}{\bf{I}}\succ\gamma^{-1}{\bf{I}}$,
where $\lambda^{{\bf{M}}}_{min}$ is a constant and $\gamma$ is the MCP
parameter. Then

$F({\bf{L}}^{(t+1)})-F({\bf{L}}^{(t)})\leq-c\cdot\|\Delta^{(t)}\|_{F}^{2},$
(57)

where $c>0$. Furthermore, the linesearch parameter $\alpha^{(t)}$ in (19) can
be chosen to be bounded away from zero, i.e.,
$\alpha^{(t)}\geq\alpha_{min}>0$.

###### Proof.

For $\Delta^{(t)}$, the solution of (48), and some linesearch parameter
$0<\alpha<1$ we have:

$\begin{array}[]{ll}\big{\langle}\nabla
f({\bf{L}}^{(t)}),\Delta\big{\rangle}+\frac{1}{2}\big{\langle}\Delta^{(t)},\Delta^{(t)}\rangle_{{\bf{M}}^{(t)}}+\rho_{MCP}\left({{\bf{L}}^{(t)}+\Delta^{(t)};\lambda}\right)\\\
\\\ \quad\quad\leq\langle\nabla
f({\bf{L}}^{(t)}),\alpha\Delta^{(t)}\rangle+\frac{1}{2}\langle\alpha\Delta^{(t)},\alpha\Delta^{(t)}\rangle_{{\bf{M}}^{(t)}}+\rho_{MCP}\left({{\bf{L}}^{(t)}+\alpha\Delta^{(t)};\lambda}\right).\end{array}$
(58)

That is because the objective in (48) is lower for $\Delta^{(t)}$ than for
$\alpha\Delta^{(t)}$ (because $\Delta^{(t)}$ is assumed to be the minimum).
Also, since $\mathcal{F}$ is convex, then
${\bf{L}}^{(t)}+\alpha\Delta\in\mathcal{F}$. Thus, we have

$(1-\alpha)\big{\langle}\nabla
f({\bf{L}}^{(t)}),\Delta\big{\rangle}+(1-\alpha^{2})\frac{1}{2}\big{\langle}\Delta,\Delta\rangle_{{\bf{M}}^{(t)}}+\rho_{MCP}\left({{\bf{L}}^{(t)}+\Delta;\lambda}\right)-\rho_{MCP}\left({{\bf{L}}^{(t)}+\alpha\Delta;\lambda}\right)\leq
0.$ (59)

After dividing by $(1-\alpha)>0$ we get:

$\big{\langle}\nabla
f({\bf{L}}^{(t)}),\Delta\big{\rangle}\leq-(1+\alpha)\frac{1}{2}\big{\langle}\Delta,\Delta\rangle_{{\bf{M}}^{(t)}}-\frac{1}{1-\alpha}\left(\rho_{MCP}\left({{\bf{L}}^{(t)}+\Delta;\lambda}\right)-\rho_{MCP}\left({{\bf{L}}^{(t)}+\alpha\Delta;\lambda}\right)\right).$
(60)

Now, by the definition of $MCP$ in (14), for every $x,y$ such that
$sign(x)=sign(x+y)$ the function is twice differentiable in the interval
$(x,x+y)$, and we can write

$MCP(x+y)=MCP(x)+MCP^{\prime}(x)y-\frac{1}{2\gamma}(\omega y)^{2},$ (61)

where $0\leq\omega\leq 1$. Given that the feasible $\mathcal{F}$ states that
${\bf{L}}^{(t)}\leq 0$ and ${\bf{L}}^{(t)}+\Delta\leq 0$ at all times, we can
write

$\rho_{MCP}\left({{\bf{L}}^{(t)}+\Delta;\lambda}\right)=\rho_{MCP}\left({{\bf{L}}^{(t)};\lambda}\right)+\langle\nabla\rho_{MCP},\Delta\rangle-\frac{1}{2\gamma}\langle\Omega\odot\Delta,\Omega\odot\Delta\rangle$
(62)

where $0\leq\Omega\leq 1$ is a matrix of the $\omega$ values as in (61), and
$\odot$ is an element-wise product. Similarly, let

$\rho_{MCP}\left({{\bf{L}}^{(t)}+\alpha\Delta;\lambda}\right)=\rho_{MCP}\left({{\bf{L}}^{(t)};\lambda}\right)+\langle\nabla\rho_{MCP},\alpha\Delta\rangle-\frac{1}{2\gamma}\langle\Omega_{\alpha}\odot\alpha\Delta,\Omega_{\alpha}\odot\alpha\Delta\rangle,$
(63)

where $0\leq\Omega_{\alpha}\leq 1$. To show a decrease in the function $F$
following the solution of (48) and a linesearch, we write

$\displaystyle F({\bf{L}}^{(t)}+\alpha\Delta^{(t)})-F({\bf{L}}^{(t)})=$
$\displaystyle
f({\bf{L}}^{(t)}+\alpha\Delta^{(t)})+\rho_{MCP}\left({{\bf{L}}^{(t)}+\alpha\Delta;\lambda}\right)-f({\bf{L}}^{(t)})-\rho_{MCP}\left({{\bf{L}}^{(t)};\lambda}\right)$
(64) $\displaystyle=$ $\displaystyle\langle\nabla
f({\bf{L}}^{(t)}),\alpha\Delta^{(t)}\rangle+\frac{1}{2}\langle\alpha\Delta^{(t)},\alpha\Delta^{(t)}\rangle_{\nabla^{2}f({\bf{L}}^{(t)})}+O(\|\Delta\|^{3})$
(65)
$\displaystyle+\rho_{MCP}\left({{\bf{L}}^{(t)}+\alpha\Delta;\lambda}\right)-\rho_{MCP}\left({{\bf{L}}^{(t)};\lambda}\right)$
$\displaystyle\leq$
$\displaystyle-\frac{\alpha(1+\alpha)}{2}\big{\langle}\Delta,\Delta\rangle_{{\bf{M}}^{(t)}}+\frac{1}{2}\langle\alpha\Delta^{(t)},\alpha\Delta^{(t)}\rangle_{\nabla^{2}f({\bf{L}}^{(t)})}+O(\|\Delta\|^{3})$
(66)
$\displaystyle-\frac{\alpha}{1-\alpha}\left(\rho_{MCP}\left({{\bf{L}}^{(t)}+\Delta;\lambda}\right)-\rho_{MCP}\left({{\bf{L}}^{(t)}+\alpha\Delta;\lambda}\right)\right)$
$\displaystyle+\rho_{MCP}\left({{\bf{L}}^{(t)}+\alpha\Delta;\lambda}\right)-\rho_{MCP}\left({{\bf{L}}^{(t)};\lambda}\right).$

Equation (65) is a simple Taylor expansion, and the inequality in (66) stems
from (60). At this point, we will focus on the last two lines of (66):

$\displaystyle-\frac{\alpha}{1-\alpha}\left(\rho_{MCP}\left({{\bf{L}}^{(t)}+\Delta;\lambda}\right)-\rho_{MCP}\left({{\bf{L}}^{(t)}+\alpha\Delta;\lambda}\right)\right)+\rho_{MCP}\left({{\bf{L}}^{(t)}+\alpha\Delta;\lambda}\right)-\rho_{MCP}\left({{\bf{L}}^{(t)};\lambda}\right)$
(67)
$\displaystyle=\frac{1}{1-\alpha}\left((1-\alpha+\alpha)\rho_{MCP}\left({{\bf{L}}^{(t)}+\alpha\Delta;\lambda}\right)-\alpha\rho_{MCP}\left({{\bf{L}}^{(t)}+\Delta;\lambda}\right)\right)-\rho_{MCP}\left({{\bf{L}}^{(t)};\lambda}\right)$
(68)
$\displaystyle=\frac{1}{1-\alpha}\left((1-\alpha)\rho_{MCP}\left({{\bf{L}}^{(t)};\lambda}\right)-\frac{1}{2\gamma}\langle\alpha^{2}\Omega_{\alpha}^{2}-\alpha\Omega^{2}\odot\Delta,\Delta\rangle\right)-\rho_{MCP}\left({{\bf{L}}^{(t)};\lambda}\right)$
(69)
$\displaystyle=\frac{1}{1-\alpha}\left(-\frac{1}{2\gamma}\langle\alpha^{2}\Omega_{\alpha}^{2}-\alpha\Omega^{2}\odot\Delta,\Delta\rangle\right)\leq\frac{\alpha}{2\gamma}\langle\max\\{\Omega_{\alpha}^{2},\Omega^{2}\\}\odot\Delta,\Delta\rangle\leq\frac{\alpha}{2\gamma}\langle\Delta,\Delta\rangle,$
(70)

where $\Omega_{\alpha}$ and $\Omega$ are defined in (62)-(63). Going back to
the bound on the decrease in $F$ in (66) we get:

$\displaystyle F({\bf{L}}^{(t)}+\alpha\Delta^{(t)})-F({\bf{L}}^{(t)})\leq$
$\displaystyle-\frac{\alpha(1+\alpha)}{2}\big{\langle}\Delta,\Delta\rangle_{{\bf{M}}^{(t)}}+\frac{1}{2}\langle\alpha\Delta^{(t)},\alpha\Delta^{(t)}\rangle_{\nabla^{2}f({\bf{L}}^{(t)})}+\frac{\alpha}{2\gamma}\langle\Delta,\Delta\rangle+O(\|\Delta\|^{3})$
(71) $\displaystyle=$
$\displaystyle-\frac{\alpha}{2}\langle\Delta,\Delta\rangle_{({\bf{M}}^{(t)}-\gamma^{-1}{\bf{I}})}-\frac{\alpha^{2}}{2}\langle\Delta,\Delta\rangle_{({\bf{M}}^{(t)}-\nabla^{2}f({\bf{L}}^{(t)}))}+O(\|\Delta\|^{3}).$
(72)

Two conditions are needed to keep the above bound negative and guarantee a
decrease in the objective. One is by having
${\bf{M}}^{(t)}\succ\gamma^{-1}{\bf{I}}$, which is easy to guarantee, e.g.,
via the choice of $\varepsilon$ in our case. The other condition, assuming
that we have ${\bf{M}}^{(t)}\succ\gamma^{-1}{\bf{I}}$, can be written as

$\alpha({\bf{M}}-\gamma^{-1})+\alpha^{2}\nabla^{2}f({\bf{L}}^{(t)})\succ 0,$

and then we have

$F({\bf{L}}^{(t)}+\alpha\Delta^{(t)})-F({\bf{L}}^{(t)})\leq-\alpha(\lambda_{min}^{{\bf{M}}}-\frac{1}{\gamma}-\alpha\theta)\|\Delta^{(t)}\|_{F}^{2}.$
(73)

where $\theta$ is the upper bound on the Hessian, and
$\lambda_{min}^{{\bf{M}}}>\gamma^{-1}$ is the smallest eigenvalue of
${\bf{M}}$. This value is negative for any
$0\leq\alpha_{min}<\frac{\lambda_{min}^{{\bf{M}}}-\gamma^{-1}}{\theta}$. ∎

The result above is relevant to Algorithm 1, where we have
${\bf{M}}^{(t)}=\nabla^{2}f({\bf{L}}^{(t)})+\mbox{diag}(\epsilon^{2})$. Hence,
the second term in (72) almost vanishes if $\epsilon$ is small. On the other
hand, $\epsilon$ allows us to make sure that the first term in (72) remains
negative. We note that one can also monitor which entries are in the concave
region of the MCP (smaller than $\gamma\lambda$ in magnitude) and set
$\varepsilon_{i,j}=\gamma^{-1}$ only for them, following (70).

The analysis above essentially is also suitable for Algorithm 2 if we set
$\beta_{t}=0$, i.e., use a projected quasi-Newton method with ${\bf{D}}$ as a
diagonal approximation to the Hessian. This is because in Algorithm 2, the
Newton problem is also a smooth objective with an MCP prior and positivity
constraints, and the change of variables to $\delta$ is only for convenience
of symmetry and does not change the problem.

###  D.5 Convergence

Before stating the main convergence theorem, we prove the following auxiliary
lemma:

######  Lemma 4.

Let $\\{{\bf{L}}^{(t)}\\}$ be series of points produced by (48) followed by
the linesearch (19), with
$\lambda^{{\bf{M}}}_{max}{\bf{I}}\preceq{\bf{M}}^{(t)}\preceq\lambda^{{\bf{M}}}_{min}{\bf{I}}\preceq
0$. Let $\\{{\bf{L}}^{(t_{j})}\\}$ be any converging subseries of
$\\{{\bf{L}}^{(t)}\\}$ and let $\bar{\bf{L}}$ be its limit. Then
$\bar{\bf{L}}$ is a stationary point of $F$ in (15).

###### Proof.

Since the sub-series $\\{{\bf{L}}^{(t_{j})}\\}$ convergences to
$\bar{{\bf{L}}}$, then $\\{F({\bf{L}}^{(t_{j})})\\}$ convergences to
$F(\bar{{\bf{L}}})$. Therefore,
$\\{F({\bf{L}}^{(t_{j}+1)})-F({\bf{L}}^{(t_{j})})\\}\rightarrow 0$, which
implies following (73) that
$||{\bf{L}}^{(t_{j})}-{\bf{L}}^{(t_{j}+1)}||^{2}_{F}\rightarrow 0$ and
$\lim_{j\rightarrow\infty}{\bf{L}}^{(t_{j}+1)}=\bar{{\bf{L}}}$. According to
(48) and the linesearch (19),
${\bf{L}}^{(t+1)}={\bf{L}}^{(t)}+\alpha^{(t)}\Delta^{(t)}$ and because,
$||{\bf{L}}^{(t_{j}+1)}-{\bf{L}}^{(t_{j})}||^{2}_{F}\rightarrow 0$ and
$\alpha\geq\alpha_{min}>0$, this implies $||\Delta^{(t)}||^{2}_{F}\rightarrow
0$. We know that $\Delta^{(t)}$ satisfies (56), and by taking the limit as
$j\rightarrow\infty$ and using Lemma 2 we get that $\bar{{\bf{L}}}$ is a
stationary point of (15). ∎

Now, after providing the auxiliary lemmas above, we are ready to state the
convergence theorem.

Theorem 2. Let $\\{{\bf{L}}^{(t)}\\}$ be a series of points produced by a
proximal quadratic minimization, constrained by $\mathcal{F}$ with some SPD
matrices
$0\prec\lambda^{{\bf{M}}}_{min}{\bf{I}}\preceq{\bf{M}}^{(t)}\preceq\lambda^{{\bf{M}}}_{max}{\bf{I}}$
as the Hessian (as in (48)), followed by the linesearch (19), starting from
${\bf{L}}^{(0)}\in\mathcal{F}$. Then any limit point $\bar{\bf{L}}$ of
$\\{{\bf{L}}^{(t)}\\}$ is a stationary point of (15).

###### Proof.

By Lemma 3, the series $F({\bf{L}}^{(t)})$ is monotonically decreasing, hence
it is also convergent as $F$ is bounded from below. Since the level set
$\mathcal{R}$ is compact by Lemma 1, and ${\bf{L}}^{(t)}$ is bounded in
$\mathcal{R}$, there is a converging subseries $\\{{\bf{L}}^{(t_{j})}\\}$,
which let $\bar{\bf{L}}$ be its limit point. By Lemma 4, $\bar{\bf{L}}$ is a
stationary point of $F$ in (15). Since $F$ is continuous then
$F({\bf{L}}^{(t_{j})})\longrightarrow F(\bar{\bf{L}})$. Since the limit of
$\\{{\bf{L}}^{(t)}\\}$ and $F({\bf{L}}^{(t)})$ equal to that of any of their
corresponding subseries, then
$\\{{\bf{L}}^{(t)}\\}\longrightarrow\bar{\bf{L}}$ and
$F({\bf{L}}^{(t)})\longrightarrow F(\bar{\bf{L}})$. ∎