Fast Projected Newton-like Method for Precision Matrix Estimation under Total Positivity
Authors: Jian-Feng CAI
José Vinícius De Miranda Cardoso
Daniel P. Palomar
Jiaxi Ying
Venue: NeurIPS 2023 poster (2023)
Keywords: MTP2
Total Positivity
Generalized graph Laplacian
Precision matrix estimation
Nonnegative partial correlations

Abstract: We study the problem of estimating precision matrices in Gaussian distributions that are multivariate totally positive of order two ($\mathrm{MTP}_2$). The precision matrix in such a distribution is an M-matrix. This problem can be formulated as a sign-constrained log-determinant program. Current algorithms are designed using the block coordinate descent method or the proximal point algorithm, which becomes computationally challenging in high-dimensional cases due to the requirement to solve numerous nonnegative quadratic programs or large-scale linear systems. To address this issue, we propose a novel algorithm based on the two-metric projection method, incorporating a carefully designed search direction and variable partitioning scheme. Our algorithm substantially reduces computational complexity, and its theoretical convergence is established. Experimental results on synthetic and real-world datasets demonstrate that our proposed algorithm provides a significant improvement in computational efficiency compared to the state-of-the-art methods.

Fast Projected Newton-like Method for Precision Matrix Estimation under Total Positivity
Jian-Feng Cail 1,2 José Vinícius de M. Cardoso', Daniel P. Palomar', Jiaxi Ying 1,2+
Hong Kong University of Science and Technology'
HKUST Shenzhen-Hong Kong Collaborative Innovation Research Institute?
Abstract
We study the problem of estimating precision matrices Gaussian distributions that multivariate totally positive of order two (MTP2). The precision matrix in such a distribution an M-matrix. This problem can be formulated as a sign- constrained log-determinant program. Current algorithms are designed using the block coordinate descent method the proximal point algorithm, which becomes computationally challenging in high-dimensional cases due to the requirement to solve numerous nonnegative quadratic programs or large-scale linear systems. To address this issue, we propose novel algorithm based on the two-metric projection method, incorporating a carefully designed search direction and variable partitioning scheme. Our algorithm substantially reduces computational complexity, and its theoretical convergence is established. Experimental results on synthetic and real-world datasets demonstrate that our proposed algorithm provides significant improvement in computational efficiency compared to the state-of setee-of-theeeee the-art methods.
Introduction
We consider the problem of estimating the precision matrix (i.e., inverse covariance matrix) in multivariate Gaussian distribution, where all the off off-diagonal elements of the precision matrix are nonpositive. The resulting precision matrix a symmetric M-matrix. Such property also known as multivariate totally positive of order (MTP2) [I] For ease of presentation, nonpositivity constraints This model arises in a variety of applications such as taxonomic reasoning 2], graph signal processing 3 factor analysis in psychometrics 4 and financial markets S
call the the off-diagonal elements of the precision matrix as MTP₂ constraints.
Estimating precision matrices under MTP, constraints is an active research topic. Recent results in 2 show that MTP₂ constraints lead to a drastic reduction on the number of observations required for the maximum likelihood estimator (MLE) to exist in Gaussian distributions and Ising models. This advantage is crucial in high-dimensional regimes with limited observations. Growing interest in estimating precision matrices under MTP; constraints is seen in graph signal processing 77.. A precision matrix fulfilling MTP₂ constraints can be regarded as generalized graph Laplacian, where eigenvalues and eigenvectors are interpreted as spectral frequencies and Fourier bases, facilitating the computation of graph Fourier transform [10] The MTP2 property has also been studied portfolio allocation (5) and structure recovery 12
Estimating precision matrices under MTP, constraints can be formulated as a sign-constrained log-determinant program. Existing algorithms for solving this problem, such as the block coordinate descent (BCD) 3 and proximal point algorithm (PPA) 99 are efficient for low low-dimensional
"Corresponding author.
37th Conference on Neural Information Processing Systems (NeurIPS 2023).

---

problems. However, they become time-consuming in high- eigh-dimensional scenarios due to the necessity of solving numerous nonnegative quadratic programs or large-scale linear systems. An alternative is the gradient projection method (13), which offers computational efficiency per iteration. Nevertheless, it often grapples with slow convergence rates in high-dimensional cases. Thus, there is demand for efficient and scalable algorithms for precision matrix estimation under MTP2 constraints.
1.1 Contributions
In this paper, we propose a fast projected Newton-like algorithm for estimating precision matrices under MTP₂ constraints. While second-order algorithms generally require fewer iterations than first-order methods, they often encounter computational challenges due to the necessity of computing a large (approximate) Hessian matrix inverse or equivalently solving the linear system. Our main contributions in this paper are threefold:
Our proposed algorithm, rooted in the two-metric projection method, stands apart from current for solving the target problem. Utilizing a well-designed
BCD and PPA approaches 2 search direction and variable partitioning scheme, our algorithm avoids the need to solve nonnegative quadratic programs or linear systems, yielding significant computational reduction compared to BCD and PPA algorithms. As second-order method, our algorithm maintains the same per-iteration computational complexity as the gradient projection method.
3
g g
2. We establish that our algorithm converges the minimizer of the target problem. Furthermore, under a mild assumption, we prove the convergence of the set of free variables to the support of the minimizer within finite iterations and provide the convergence rate of our algorithm.
3. Numerical experiments both synthetic and real-world datasets provide compelling evidence that our algorithm converges to the minimizer considerably faster than the compared methods. We apply the proposed method to financial time-series data and observe significant performance terms of modularity value on the learned financial networks.
Notation: Lower case bold letters denote vectors and upper case bold letters denote matrices. Both X1, and [X]s Jay denote the (i, entry of X. [p] denotes the set {.,....P}, and [p] denotes Il....,p] Let @ be the Kronecker product and @ be the entry-wise product. supp X) = {i, E Ix, + 0}. || |m = maxi,j XXay and || Il ming, Xxyl- S and Sº. denote the sets of symmetric positive semi-definite and positive definite matrices with dimension p p. vec(X) and X denote the vectorized version and transpose of X.
the set {1, pp) x
pp}
Problem Formulation and Related Work
In this section, we first introduce the problem formulation, then present related works.
2.1 Problem formulation
Let y = (B1...., Bp) be a p-dimensional random vector following y ~ *N(, where Σ is the covariance matrix. We focus on the problem of estimating the precision matrix O = --1 given n i.i.d. observations y ",.,..,,(()) Let S Σ (+) (y( be the sample covariance matrix. Throughout the paper, the sample covariance matrix is assumed to have strictly positive diagonal elements, which holds with probability one. This is because some diagonal element SJj is zero if and
(")),
only if the j-th element of y()) must be zero for every i E [p], which holds with probability zero.
We consider solving the following sign-constrained log determinant program:
X = arg min XEMP
- -log det (x) + tr (XS) - Eizy Ais [Xil,
subject to X. 0, v (i.j) E,
where Aij is the regularization parameter, E is the disconnectivity set with each node pair forced to disconnect, and MP is the set of all p-dimensional, symmetric, non-singular M-matrices defined by
M X E I S 0, Vi H))
The disconnectivity set in I can be obtained in several ways: (i) it is often the case that some edges between nodes must not exist due prior knowledge; (ii) it can be estimated from initial estimators; (iii) it can be obtained in some tasks of learning structured graphs such as bipartite graph [6].

---

2.2 Related work
Estimating precision matrices under Gaussian graphical models has been extensively studied in the literature. One well-known method is graphical lasso (17-20), which minimizes the f--regularized Gaussian negative log-likeling likelihood. Various algorithms were proposed to solve this problem including first-order methods 17
18 21 and second-order methods 30 The graphical lasso optimization problem is unconstrained and nonsmooth, while Problem ① is smooth and constrained. The difficulties in solving the two problems are inherently different, and the algorithms mentioned above cannot be directly extended to solve Problem IT.
Recent studies 2 3 D employed BCD and PPA-type algorithms estimate precision matrices under MTP₂ constraints. In (2) the primal variable is updated one column/row at a time by solving a nonnegative quadratic program, cycling until convergence. The work [ follows a similar approach but addresses the dual problem, improving memory efficiency. Both works target Problem (T) without the disconnectivity constraint. A BCD-type algorithm was proposed in []] to accommodate disconnectivity constraints. However, the computational complexity of these algorithms, at O(p*) operations per cycle, becomes prohibitive for high-dimensional problems. Alternatively, recent work O introduced inexact PPA algorithm to solve transformed problem, derived using the noft-thresholding technique. However, this algorithm demands the computation of an inexact Newton direction from a p² xp linear system at every iteration within the inner loop, presenting computational difficulties in high-dimensional scenarios.
The proposed algorithm adopts the two-metric projection framework [3]. incorporating distinct metrics for search direction and projection. A representative method in this framework the projected Newton algorithm [35], originally designed for nonnegativity constrained problems. However, it is unsuitable for Problem ① due to its O(p") operations needed to compute the inverse of the Hessian at each iteration. To mitigate computation and memory costs, the projected quasi-Newton algorithm with limited limited-memory memory
(L-BFGS) was introduced in (36), requiring O((m + p)p?) operations per iteration, with m the number of iterations stored for Hessian approximation. By leveraging the structure of Problem (T), this paper carefully designs the search direction and variable partitioning scheme to substantially reduce computation and memory costs, achieving the same orders per iteration the gradient projection method [3]
3 Proposed Algorithm
In this section, we propose a fast projected Newton-like algorithm to solve Problem ①.. The constraints in Problem (I) can be rewritten X E Q + where Qt is defined as
Ω = {X [RxxP| X.j = O, (i.j) E X., S 0, Vi Vi
The set Ω is convex and closed, thus this constraint can be handled by projection Ps defined by (i,5) E
(Po(A))
Aij min(A:j. 0) if (i,j) & & and #
if i
(3)
The positive definite set S is not closed and cannot be managed by a projection, which will be handled using the line search method in Section 3.3 Let f denote the objective function of Problem (.] To address Problem II. start with the gradient projection method, expressed as:
X = Po (x - 7M V (Xx)),
where Yo is the step size. To accelerate convergence, one may consider
X Pa Xx 7 Px),
(5)
in which Pa E IRPxP is a search direction defined by vec (Pk) = M,-" 'vee (VS(X))),
'vvee
(6)
where M, E TRP² xp is positive definite symmetric matrix, incorporating second-order derivative information. If we adopt Mg the Hessian matrix, then Ps. becomes the Newton direction as shown in Proposition[..1 with the proof provided in Appendix B. and iterate (5) can be viewed as a natural adaptation of the unconstrained Newton's 's method. Regrettably, the convergence of such an iterate to the minimizer cannot be guaranteed, as Pş may not be descent direction here, which is supported by numerical results in Section 5.1.I Similar observations have also been reported in [35 37
3

---

Proposition 3.1. If M2 is constructed as the Hessian matrix H₂ of Problem 1 then the search direction Pk defined in (6) can be written as Ps V
Remark 3.2. We refer iterate (5) as the two-metric projection method [4 it adopts two distinct metrics: the search direction Pk induced by the quadratic norm |l-II defined by || =
(vec (A), Ve (A)) and the projection Pa with respect to the Frobenius norm ||-1lp-
3.1 Identifying the sets of restricted and free variables In order guarantee iterate (5) converge the minimizer, we partition the variables into two groups, Le., restricted and free variables, and update the two groups separately. We first define set T(X,c) with respect to X E 8 and c E R++ T(X,c) = ((i.j) (alp)))ppPl---e - S O, [Vf())l < 0}. Then at the k-th iteration, we identify the set of restricted variables based on X as follows,
Is : T(x:; UE),
(8)
where E is the disconnectivity set from Problem ① and é is small positive scalar. For any (i,5) E (T (x), X3 in the next iterate is likely to be outside the feasible set (i.e., X > 0) if we remove the projection Pa. as it near zero and moves towards the positive direction using the negative of the gradient as the search direction. Therefore, we set all restricted variables to zero.
To establish the theoretical convergence of the algorithm, the positive scalar ek in is specified as
E/ = min (2(1 - amm? [ Ilmin 0).
where is a positive scalar (See Lemma B. in Appendix), a E (0, 1) is parameter in the line search condition, and To represents the set T(X&, 5). In the rare event that Ta is empty, particularly in sparse settings, we define E = implying an empty 7(xx, f1.) in (8) The parameter satisfies 0 8 < _espp.))) *kyl where X is the minimizer of Problem ( Such condition can be ensured by setting sufficiently small positive 5. Then in (9) is nearly equal to 6. From an implementation view, we can directly set a small positive 6g. resulting in the algorithm performing well in practice. The set of free variables, denoted by Ig. the complement of Ig.
3.2 Computing approximate Newton direction
While the (approximate) Newton direction usually demands a considerably higher computational cost than the gradient, our designed direction maintains the same computational order as the gradient.
At the k-th iteration, we first partition X into two groups, [Kl]z, and [Xu]zę> where [Xu]I, E [R and [X&]z RIZE denote two vectors containing all elements of Xx in the sets In and Ig. respectively. Then we can rewrite the search direction Pk in (6) as follows,
pvec, (P.) = Q pvec; (V fXxx))
(10)
where pvec, stacks P& into vector, similar to veo (x)). but places elements from Iī. first, followed by those from I.- Qx is obtained by permuting the rows and columns of M;" in (6. To enhance computational efficiency, we propose constructing Qx and M, as follows:
Q = [M -1 I
[M Ig
(11)
[M N,
[M TTT.
O
D
where D₁ II is a positive definite diagonal matrix, and " Izzzz E RIZEEXIZEI is principal submatrix of H.)), preserving rows and columns indexed by ZE. Here, H₃ is the Hessian matrix at Xg. The construction of M 11 in (TT) is crucial for defining the search direction, enabling computation and memory costs comparable to the gradient projection method while effectively incorporating second-order derivative information.
Next, compute the approximate Newton direction Pa over the IE and present the iterate X (X+41lze- We define projection Pzz (A) follows,
{ A if E
[Pig (A),,, = {o otherwise.
(12)

---

Leveraging the well-designed gradient scaling matrix Qs in (TI. we can efficiently compute the approximate Newton direction, as demonstrated in Proposition 3.3
with proof in Appendix B.2
Proposition 3.3, If Mr is constructed by (T)). then the search direction P& defined in (6) over I can be written as [Pxlz PI (Vf()))) Xelzg Zg
Using the search direction from Proposition 3.3 update X over I as follows,
1 = Psi [Xx - P (x)
(13)
For the restricted variables in the set In, we directly set them to zero, i.e., [x++1I1a = 0.
3.3 Computing step size
We adopt an Armijo-like rule for step size selection, ensuring the global convergence of our algorithm.
Based on the iterate proposed Section 3.2 we define (7())) with [Xx(m))z,, = 0 and
7k X ( (Vf [X*(T))TzE Pa (Xxl
(14)
We test step sizes E 130,80 B,, ,?.,.... with 3 E (0, 1). until we find the smallest t E N such that (xu),,, with 7h 3r, satisfies (x(7k) E $4 and the line search condition:
f(Xa(T)) f(Xx) an (V (Vx)) [Polzz) a (VVVf(X))) (Klz,),
(15)
where E 1) is scalar. We then set X 1 (x())). Positive definiteness of X can be verified during Cholesky factorization for objective function evaluation. is worth mentioning that working with the positive semi- -emindefiniseness constraint on X instead of positive definiteness would not change anything in the algorithm if we keep the line search, as the positive definiteness is automatically enforced due to the form of the objective function.
The line search condition (15) is variant of the Armijo rule. Condition (5) can be always satisfied for small enough step size shown in Proposition 3.4 Define the feasible set of Problem as
ur X E RPKP| X MM, X.j 0, v (i,j) E}.
(16)
For any given X uP, define the lower level set of the objective function f for Problem !
Lj = (X uP f(X) S f(X))
(17)
Proposition 3.4. For any Xx E Li, there exists a F such that (x( S.+ and the line search condition (15) holds for any 7 E (0, 7x).
The proof of Proposition 3s4 available in Appendix B.3
We demonstrate that Xx(Tx) ensures
decrease of the objective function value Proposition 3.5 proved in Appendix B.4
Proposition 3.5. For any X E Ly. if (T()) satisfies the line search condition
15 then we have
f(Xx(T))) f (X+) arkm a (
where m is a positive scalar (See Lemma B. in Appendix).
3.4 Computation and memory costs
In each iteration. our algorithm calculates the gradient, performs two matrix multiplications, and conducts two projections, with respective computational costs of O O(p²), and O(p²). In current implementation of the line search method, we first conduct the Cholesky factorization LL using MATLAB's s "chol" function. This function can simultaneously verify the positive definiteness of X. Next, we calculate the log- determminnnn function as log (e)(X)) = 2, log(L.)). The Cholesky factorization the most computationally demanding step. generally requiring (pp)) costs for p X p matrix.
5

---

Algorithm Fast Projected Newton-like (FPN) algorithm
1: Input: Regularization parameter Aej. a, 3, S, and Xo:
2: for k 0,..2..... do
3: Identify the restricted set Iş and free set Ig according to 8
4: Compute the search direction over Ig. [Pxlag = [Xx (Xx)) Xxlzg'
5: +
6: repeat
7:
Update X (Xe+ilz, and [Xa- Prr((xaIEE 3t (1(Pal)));
8: 1 ++;
9:
until X X8+1 > 0, and
((x+t) -x
(Pllzy)))))))))) - (1x
(xalz.)
10: end for
To mitigate the computational burden associated with Cholesky factorization, we suggest a more efficient method for evaluating the log-determinant function and verifying positive definiteness, as presented in 38) This method, which leverages Schur complements and sparse linear system solving, can tackle problems of up to 10 dimensions. Furthermore, it worthwhile to investigate more efficient strategies for computing an approximate log determiint function. In this context, the approach proposed in (39] offers a nearly linear scaling of execution time with the number of non-zero entries, while maintaining a high level of accuracy.
In summary, algorithm has overall complexity of O(p²) per iteration. BCD-type algorithms 2 g need (pp") operations per cycle, while projected quasi-Newton with L-BFGS [36 requires O(mm+ + +))p") operations per iteration, with m as the stored iteration count for Hessian approximation. The PPA algorithm [9] requires computing an inexact Newton direction from a p² x p² linear system during each inner loop iteration, with the exact complexity not established. In addition, our algorithm, gradient projection and BCD-type methods
2 3 g need ((p²) memory costs, while projected quasi-Newton with L-BFGS 36 requires (mmp² and PPA demands Opp*).
Convergence Analysis
Prior to delving into the convergence analysis, first establish the uniqueness of the minimizer for Problem ( and determine the sufficient and necessary conditions for point to be the minimizer.
Theorem 4.1. The minimizer of Problem and only if it satisfies
is unique, and point X E M is the minimizer if
[x = (i.j) E, [Vf(X"))vE s 0, and Vf(X = 0,
(18)
where V ((.,)) E [X|] = 0}.
The proof of Theorem is available in Appendix B.5 The following theorem shows that algorithm converges to the minimizer of Problem (T.
Theorem 4.2. The sequence {Xa} generated by Algorithm 7 with Xo E L, converges to the minimizer X* of Problem [ with {f((x)) monotonically decreasing.
The proof of Theorem 4.2 is available in Appendix B.6 It is worth noting that constructing an initial
point X, E Ly, as defined in (I7). is straightforward. Please refer to the proof of Theorem 4.2 for more details on this. The theoretical analysis on support set convergence and sequence convergence rate relies on the following assumption.
Assumption 4.3, For the minimizer X of Problem Tb. we assume that the gradient of the objective function f at X satisfies
[f((') <0, (i.j) VIE,
where V= = ((.5) E [p]] x
0}. and E is the disconnectivity set.

---

Theorem 4.4. Under Assumption 4.3 to the support of the minimizer of Problem (T) in finite iterations. In other words, there exists some Ko E N such that I (oopp(X")for any k 2 Ko-
the set of free variables I& generated by Algorithm 7 converges
The proof of Theorem 4.4 is provided in Appendix [B77 Theorem 4.4 demonstrates that the set free variables constructed by our variable partitioning scheme exactly identify the support of X in finite iterations. Now we establish the convergence rate of our algorithm. Define
R
H
IgI I
[HA] [Hx] [
IVIN
III
(19)
Theorem 4.5. Under Assumption
4.3
the sequence {Xx) generated by Algorithm satisfies
X
x Mr.
2(1 a) 2
limsup 8-+00
S lim
1 1--min/
min
mm
Xx
x+1
10
Ma
Mr
where my and Mg as the smallest and largest eigenvalues R
H
R2 respectively. T
The proof of Theorem 4.3 is provided in Appendix B.8 Theorem 4.5 reveals that the convergence rate of our algorithm depends the condition number / & of R2 Zzzz R2 Replacing Rg with an identify matrix, (i.e., using the projected gradient method) results in a rate dependent on
the condition number of H Izzz-- The condition number of R H Jz R could be larger, as Rg could approximate [H Lepgg well. Thus, the gradient scaling matrix Rg i.e., [Hg lzzzg in leads our algorithm to converge faster than the projected gradient method.
It is important note that our algorithm generally does not achieve superlinear convergence, despite the incorporation of second-order information. Superlinear convergence necessitates that the inverse gradient scaling matrix progressively approximate the Hessian at the minimizer [40]. However, this is a condition that our constructed scaling matrix does not meet. Despite this, we should note that constructing scarch direction to achieve superlinear convergence proves significantly more computationally demanding than our approach, as it cannot leverage the special structure of the Hessian decrease the computational load.
Experimental Results
We conduct experiments on synthetic and real-world data verify the performance of our algorithm. All experiments were conducted on 2.1OGHZ Xeon Gold 6152 machines and Linux OS, and all methods were implemented in MATLAB. State-of-the-art methods for comparisons include:
BCD
2
Updates each column/row of primal variable using nonnegative quadratic program.
optGL
Similar to BCD but solves nonnegative quadratic programs on the dual variable.
GGL
Similar to BCD but handles disconnectvvity constraints, while BCD and optGL cannot.
PGD
3
Projected gradient descent method with backtracking line search.
APGD
Accelerated projected gradient algorithm with extrapolation step.
PPA
Inexact proximal point algorithm with Newton-CG method.
33 Projected quasi-Newton method using limited- meetrry BFGS.
Note that all state-of-the-arr methods listed above converge to the minimizer of Problem (1). and we focus on the comparison of computational time for those methods. To that end, we report the relative error of the objective function value as a function of the run time, which is calculated by
If(())) f(X*) (X
(20)
where f is the objective function of Problem and X is its minimizer. The X is computed by running the state-of -ia---ffthhheeeee
the-an art method GGL
3 until it converges to point X E MP satisfying
Xxl, 0 V (i, E
[V f x)) I 0,
S 1
(21)
where A := (i. )) E [p| || < 10-8). Through the comparison with the sufficient and
necessary conditions the unique minimizer of Problem (1) presented in Theorem we can see that any point X₃ satisfying the conditions in (1 is sufficiently close to the minimizer.

---

Newton direction
10
Proposed
Iterations
(a) BA graph degree one
(b) graph of degree two
(c) Convergence comparison
Figure 1: Illustration of BA graphs of degree (a) and degree two (b) with 500 nodes. (c) Convergence comparison between our algorithm and the Newton direction from Proposition 3.1
5.1 Synthetic data
We generate independent samples y ..... y(()) E I2 from a multivariate Gaussian distribution with zero mean and precision matrix θ where O E Mr is the underlying precision matrix associated with a graph consisting
nodes. We use the Barabási-Albert (BA) model [42 to generate the support of the underlying precision matrix. To help readers to know well the BA graphs, present two examples in Figure I
More details about experimental setting are provided in Appendix A
5.1.1 Comparisons of search directions
We evaluate the convergence of algorithms with different search directions. Figure (c) demonstrates
that algorithm, using the direction from Proposition 3.3 converges the minimizer, aligning with our theoretical convergence results in Theorem 4.2 In contrast, the algorithm using iterate (5)
and the Newton direction from Proposition stops decreasing the objective function value after few iterations, indicating that this direction cannot be consistently regarded descent direction. The algorithm selects the step size through the Armijo rule, i.e., it is continually reduced until decrease the objective function achieved.
5.1.2 Comparisons computational time
We evaluate the computational time of our algorithm and state-of-the-art methods on synthetic datasets, averaging results over 10 realizations. We plot markers every 10 iterations for PGD, APGD, PQN-LBFGS, and FPN, while marking each cycle of updating all columns/rows for BCD-type algorithms (BCD, GGL, and optGL) and each outer iteration for PPA.
Figure 2 compares the computational time of various methods for solving Problem 4 on BA graphs degree one. Our proposed FPN significantly outperforms all state-of-the-art methods in convergence time for node counts ranging from 1000 to 5000. BCD and GGL efficient at 1000 nodes, being faster than PGD and APGD, and competitive with PQN-LBFGS and PPA. However, at 5000 nodes, BCD and GGL become slower due to the O(p*) operations per cycle required to solve nonnegative quadratic programs, leading to rapidly increasing computational costs.
RED
ACSI
PFA
TDA
DPA
APOD
Proporal
(sacedd
(a)p= 1000
(b) p 3000
(c)
5000
Figure 2: Relative errors of the objective function values versus time on BA graphs degree one.

---

COL
12
PLD PON-LAOCS
POD
PURI
POJ LECOS
21
130
20011
240
CPU (wcoed))
CPU stane (xcoesd
CPL tinz (scoonda)
(a) p 1000
(b) p 3000
(c) 5000
Figure 3: Relative errors of the objective function value versus time BA graphs of degree two.
OBR PCD PON- DNPOS
Cet
PGD
NCTI
1QK-11FFO
PGD LXN-S
Tan EIDRR CoUImdo
18:11
21 5-11
2.
CEU tima (sacedd
Number ofmoorrreeeeeeeeeeee
(a)
(b)
(c)
Figure 4: Relative enrors of objective function values versus time data sets: (a) BA graph of degree one and (b) BA graph of degree two, imposing disconnectivity constraints. (c) Run time versus numbers of nonzero elements in precision matrices across varying regularization parameter values, with the underlying precision matrix having 2998 nonzero elements (indicated by vertical red line).
Figure presents the computational time of different methods on BA graphs of degree two. As with degree graphs, FPN outperforms state-of oftte-offthe-art methods in computational time for varying node counts. PQN-LBFGS and FPN require fewer iterations to converge to the minimizer than PGD and APGD, particularly in high dimensions (e.g., p = 5000), indicating faster convergence. This because both PON-LBFGS and FPN utilize second-order information and approximate Newton direction, overcoming low convergence rates of first-order methods in high-dimensional cases.
Figure 4 (a) and (b) compare the computational time of various algorithms solving Problem (T) with disconnectivity constraints, where FPN consistently converges fastest. (c) evaluates the impact of the estimated precision matrix's 's sparsity level on run time. BCD, GGL, PPA, and FPN exhibit stable run time across varying sparsity levels, highlighting their robustness regarding regularization parameter settings, while other methods display increased run time as sparsity decreases.
5.2 Real-world data
We perform experiments on two real-world datasets: the concepts dataset and a financial time-series dataset. For the concepts dataset, we compare computational time of different algorithms solving Problem (T. The experimental results the financial time-series dataset are provided in Appendix where we examine the performance of our method in graph edge recovery.
The concepts dataset 4 from Intel Labs, comprises 1000 nodes and 218 semantic features, with p 1000 and n 218. Nodes represent concepts like "house," "coat," and "whale," while semantic features are questions like "Can fly?", "Is alive?", and "Can you use it?". Responses, collected via Amazon Mechanical Turk, range from "definitely no" to "definitely yes" on a five-point scale.
Figure 5 (a) compares the run time of various algorithms solving Problem ①I on the concepts dataset. Our proposed algorithm converges to the minimizer considerably faster than state-of-the of-the------rt hetart algorithms, which is consistent with the observations synthetic experiments. Note that all compared algorithms can reach the minimizer of Problem (I. and thus learn the same graph.

---

10
DGL.
BCD
8-
SpAGT:
PPA
L8FOS
Proposed
10
1071
1377
00
20
30
CPU time (seconda)
(b)
Figure 5: (a) Relative errors of the objective function values versus time for different algorithms in solving Problem ① on the concepts data set, consisting of 1000 nodes. (b) The connected subgraph illustrated by the minimizer on the concepts data set, which consists of 132 nodes.
Figure (b) displays connected subgraph illustrated by the minimizer of Problem Interestingly, it is observed that the learned graph forms semantic network, where related concepts closely connected. For instance, insect concepts such as "bee", "butterfly". "flea", "mosquito", and "spider" are grouped together, while human- melane concepts like "baby", "husband", "child", "girls", and "man" form another group. Moreover, the network connects "penguin" closely birds like "owl' and "crow" and sea animals like "goldfish" and "seal", highlighting its aquatic bird nature. Overall, the learned network effectively captures concept relationships.
Conclusions and Discussions
In this paper, have introduced fast projected Newton-like method for estimating precision matrices under MTP2 constraints. Our algorithm, leveraging the two-metric projection method, stands out from existing BCD and PPA-type approaches for addressing the target problem. The proposed algorithm
only straightforward to implement but also efficient terms of computation and memory usage. We have provided theoretical convergence analysis and conducted extensive experiments, which clearly demonstrate the superior efficiency of our algorithm in computational time, outperforming state-of-ther methods. Moreover, we have observed significant performance of our method in terms of modularity value the learned financial time-series graphs
Finally, we discuss the limitations our paper. Our algorithm is proven to converge to the minimizer without any assumptions; however, require Assumption 4.3 establish support set convergence in finite iterations and to determine the convergence rate. This assumption is relatively mild, as Theorem 4.1 shows that the minimizer X must satisfy [VS(X*))s; & for each (i,j) E V.EE The only additional requirement in Assumption 4 is the strictness of this inequality. However, the conditions for ensuring this strictness remain unclear. As this assumption is equivalent to the strict complementary slackness condition in optimization theory, exploring verifiable conditions to guarantee Assumption 4.3 could enrich our algorithm's insights.
7
Acknowledgements
This work supported by the Hong Kong Research Grants Council GRF 16207820, 163 16310620, and 16306821, the Hong Kong Innovation and Technology Fund (ITF) MHP/009/20, and the Project of Hetao Shenzhen-Howe Kong Science and Technology Innovation Cooperation Zone under Grant BZZB-KCCZB-202003: We would also like thank the anonymous reviewers for their valuable feedback the manuscript.
10

---

References
[1] S. Fallat, S. Lauritzen, K. Sadeghi, C. Uhler, N. Wermuth, and P. Zwiernik, "Total positivity in
Markov structures," " The Annals of Statistics, vol. 45, no. 3, pp. 1152-1184, 2017.
[2] M. Slawski and M. Hein, "Estimation of positive definite M M-matrices and structure learning for attractive Gaussian Markov random fields," Linear Algebra and its Applications, vol. 473, pp. 145-179, 2015.
[3] H. E. Egilmez, E. Pavez, and A. Ortega, "Graph learning from data under Laplacian and structural constrints," IEEE Journal of Selected Topics in Signal Processing, vol. 11, no. 6, pp. 2017.
[4] S. Lauritzen, C. Uhler, and P. Zwiernik, "Maximum likelihood estimation in Gaussian models under total positivity," The Annals of Statistics, vol. 47, no. 4, pp. 1835-1863, 2019.
[5] R. Agrawal, U. Roy, and C. Uhler, "Covariance Matrix Estimation under Total Positivity for Portfolio Selection," Journal of Financial Econometrics, vol. 20, no. 2, pp. 367-389, 09 2020.
[6] S. Lauritzen, C. Uhler, and P. Zwiernik, "Total positivity in exponential families with application to binary variables," The Annals of Statistics, vol. 49, no. 3, pp. 1436-1459, 2021.
[7] E. Pavez, H. E. Egilmez, and A. Ortega, "Learning graphs with monotone topology properties and multiple connected components," IEEE Transactions on Signal Processing, vol. 66, no. 9, pp. 209223992222 2413, 2018.
[8] E. Pavez and A. Ortega, "Generalized Laplacian precision matrix estimation for graph signal processing," " in 2016 IEEE International Conference on Acoustics, Speech and Signal Processing,
2016, pp. 6350-6354.
[9] 7. Deng and A. M.-C. So, "A fast proximal point algorithm for generalized graph laplacian learning," in ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2020, pp. 5425-54555555
[10] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, "The emerging field of signal processing on graphs: Extending high high-dimensional data analysis to networks and other irregular domains," IEEE Signal Processing Magazine, vol. 30, no. pp. 83-98, 2013.
[11] J. Kelner, F. Kochler, R. Meka, and A. Moitra, "Learning some popular Gaussian graphical models without condition number bounds," in Advances in Neural Information Processing Systems, vol. 33, 2020, pp. 10 10986-10998.
[12] Y. Wang, U. Roy, and C. Uhler, "Learning high-dimensional Gaussian graphical models under total positivity without adjustment of tuning parameters," in International Conference on Artificial Intelligence and Statistics, vol. 108, 2020, pp. 2698-2708.
[13] J. Ying. J. V. d. M. Cardoso, and D. P. Palomar, "Adaptive estimation of graphical models under total positivity," in International Conference on Machine Learning, 2023, pp. 40 054-400074
[14] S. Kumar, J. Ying, J. V. d. M. Cardoso, and D. P. Palomar, "A unified framework for structured graph learning via spectral constraints," Journal of Machine Learning Research, vol. 21, no. 22, pp. 1 1-60, 2020.
[15] V. d. M. Cardoso, Ying, and D. Palomar, "Learning bipartite graphs: Heavy tails and multiple components," in Advances in Neural Information Processing Systems, vol. 35, 2022,
pp. 14 14044-14057.
[16] S. Kumar, J. Ying, J. V. d. M. Cardoso, and D. P. Palomar, "Structured graph learning via Laplacian spectral constraints," in Advances in Neural Information Processing Systems, 2019, pp. 11 647-11 658.
[17] O. Banerjee, L. E. Ghaoui, and A. d' Aspremont, "Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data," Journal of Machine Learning Research, vol. 9, no. Mar, pp. 485-516, 2008
[18] A. d' A'spremoot, O. Banerjee, and L. El Ghaoui, "First-order methods for sparse covariance selection," SIAM Journal on Matrix Analysis and Applications, vol. 30, no. 1, pp. 56-66, 2008
11

---

[19] P. Ravikumar, M. J. Wainwright, G. Raskutti, and B. Yu, FHigh-dimensionn covariance estimation by minimizing & ---penalieed log-determinant divergence," Electronic Journal of Statistics, vol. 5. pp. 935-980, 2011.
[20] M. Yuan and Y. Lin, "Model selection and estimation in the Gaussian graphical model," Biometrika, vol. 94, no. 1, pp. 19-35, 2007.
[21] J. Duchi, S. Gould, and D. Koller, "Projected subgradient methods for learning sparse Gaussians," in Conference on Uncertainty in Artificial Intelligence, 2008, p. 153-160.
[22] J. Friedman, T. Hastie, and R. Tibshirani, "Sparse inverse covariance estimation with the graphical lasso," Biostatistics, vol. 9, no. 3, pp. 432-441, 2008.
[23] Z. Lu, "Smooth optimization approach for sparse covariance selection," SIAM Journal on Optimization, vol. 19, no. 4, pp. 1807-1827, 2009.
[24] Q. Sun, K. M. Tan, H. Liu, and T. Zhang, "Graphical nonconvex optimization via an adaptive convex relaxation," in Proceedings of the 35th International Conference on Machine Learning, vol. 80, 2018, pp. 4810-4817.
[25] K. Scheinberg, S. Ma, and D. Goldfarb, "Sparse inverse covariance selection via alternating linearization methods," in Advances in Neural Information Processing Systems, vol. 23, 2010, pp. 2101-2109
[26] A. Eftekhari, D. Pasadakis, M. Bollhöfer, S. Scheidegger, and O. Schenk, "Block-enhanced precision matrix estimation for large-scale datasets," Journal of computational science, vol. 53,
101389, 2021.
[27] J. Chen, P. Xu, L. Wang. J. Ma, and Q. Gu, "Covariate adjusted precision matrix estimation via nonconvex optimization," " in Proceedings of the 35th International Conference on Machine Learning, vol. 80, 2018, pp. 921-930.
[28) P. Xu, J. Ma, and Q. Gu, "Speediigg up latent variable Gaussian graphical model estimation via nonconvex optimization," in Advances in Neural Information Processing Systems, vol. 2017, pp. 1933-1944.
[29] E. Treister and S. Turek, "A block-coordinate descent approach for large-scale sparse inverse covariance estimation," Advances in neural information processing systems, vol. 27, 2014.
[30] Q. T. Dinh, A. Kyrillidis, and V. Cevher, "A proximal Newton framework for composite minimization: Graph learning without Cholesky decompositions and matrix inversions," in International Conference Machine Learning, vol. 28, 2013, pp. 271-279.
[31] C.-J. Hsieh, M. A. Sustik, I. S. Dhillon, and P. Ravikumar, "QUIC: quadratic approximation for sparse inverse covariance estimation," " The Journal of Machine Learning Research, vol. 15, no. pp. 2911 1-2947, 2014.
[32] L. Li and K.-C. Toh, "An inexact interior point method for L1 Le-reullarzeed sparse covariance selection," Mathematical Programming Computation, vol. 2, no. 3, pp. 291-315, 2010.
[33] C. Wang, D. Sun, and K.-C Toh, "Solving log-determinant optimization problems by a Newton- CG primal proximal point algorithm," " SIAM Journal on Optimization, vol. 20, no. 6, pp. 2994-3013, 2010.
[34] E. M. Gafni and D. P. Bertsekas, "Two-metric projection methods for constrained optimization."
SIAM Journal on Control and Optimization, vol. 22, no. 6, pp. 936-964, 1984.
[35] D. P. Bertsekas, "Projected Newton methods for optimization problems with simple constraints."
SIAM Journal on Control and Optimization, vol. 20, no. 2, pp. 221-246, 1982.
[36] D. Kim, S. Sra, and I. S. Dhillon, "Tackling box-constrained optimization via a new projected Quasi-Newton approach," SIAM Journal on Scientific Computing, vol. 32, no. 6, pp. 3548-3563, 2010.
[37] "Fast Newton-type methods for the least squares nonnegative matrix approximation problem," in SIAM International Conference on Data Mining, 2007, pp. 343-354.
[38] C.-J. Hsieh, M. A. Sustik, I. S. Dhillon, K. Ravikumar, and R. Poldrack, "BIG & QUIC: Sparse inverse covariance estimation for million variables." in Advances in Neural Information Processing Systems, vol. 26, 2013, pp. 3165-3173.
12

---

[39] I. Han, D. Malioutov, and J. Shin, "Large-scale log- degeretrmmnant computation through stochastic Chebyshev expansions," in International Conference on Machine Learning. PMLR, 2015, pp. 908-917.
[40] Nocedal and S. Wright, Numerical optimization. Springer Science & Business Media, 2006.
[41] Y. Nesterov, Introductory lectures on convex optimization: basic course. Springer Science
& Business Media, 2003, vol. 87.
[42] A.-L. Barabási and R. Albert, "Emergence of scaling in random networks," Science, vol. 286, no. 5439, pp. 509-512, 1999.
[43] B. Lake and J. Tenenbaum, "Discovering structure by learning sparse graphs," in Proceedings of the 33rd Annual Cognitive Science Conference. 2010, pp. 778-.83.
[44] Ying. J. d. M. Cardoso, and D. Palomar, "Nonconvex sparse graph learning under Laplacian constrained graphical model," in Advances in Neural Information Processing Systems, vol. 33, 2020, pp. 7101-7113
[45] E. J. Candes, M. B. Wakin, and S. P. Boyd, "Enhancing sparsity by reweighted & minimization," Journal of Fourier analysis and applications, vol. 14, no. 5, pp. 877-905, 2008.
[46] Ying. J. M. Cardoso, and D. Palomar, "Minimax estimation of Laplacian constrained precision matrices," in International Conference Artificial Intelligence and Statistics, vol. 130, 2021, pp. 3736-3744.
[47]
"A fast algorithm for graph learning under attractive Gaussian Markov random fields," in 202/ I 55th Asilomar Conference on Signals, Systems, and Computers, 2021, pp. 1520-1524.
[48] M. E. Newman, "Modularity and community structure in networks," " Proceedings of the National Academy of Sciences, vol. 103, no. 23, pp. 8577-8582, 2006.
[49] R. J. Plemmons, "M-matrix characterizations. I-nonsingular M-matrices," Linear Algebra and its applications, vol. 18, no. 2, pp. 175-188, 1977.
[50] P. Drábek and J. Milota, Methods of nonlinear analysis: applications to differential equations. Springer Science Business Media, 2007.
[51] D. Bertsekas, Nonlinear programming, 3rd ed. Athena Scientific, Belmont, MA, 2016.
[52] "On the Goldstein-Levviin-PPllaa gradient projection method," IEEE Transactions on Automatic Control, vol. 21, no. 2, pp. 174-184, 1976.
13

---

In the following sections, we provide more details on experimental settings and additional experimental results for financial time-series data in Appendix (n]a the proofs for Propositions 3.1
3.3
and 3.5 as well as Theorems 4.1 4.4 a in Appendix
A Additional Experiments
We first provide in-depth description of the experimental settings, then present numerical experiments carried out on the financial time-series data.
A.1 Experimental settings
We examine Barabasi-Albert (BA) graphs 42 as a model for the support structure of the underlying precision matrix. BA models hold significant position in network science due their ability to generate random scale-free networks via a preferential attachment mechanism. This mechanism ensures that newly introduced nodes are more likely to connect with nodes that possess a higher degree during the network's evolution. Scale-free networks serve as suitable models for various systems, including the Internet, protein interaction networks, citation networks, as well as the majority of social and online networks
In BA graph with degree r, each new node connects to pre-existing nodes, with the probability of connection being proportional to the number of edges the existing nodes currently have. In this paper, we consider values of and
Following the procedures detailed in 44 we assign a positive weight to each edge of a graph and set a zero weight for disconnected nodes. Positive weights are uniformly sampled from ((2 5). This process results weighted adjacency matrix A containing all the graph weights. Then we adopt the procedures outlined in [ for generating the underlying precision matrix + We first set
e = 81 A, with 1.055mmax (A),
(22)
where Amax (A) represents the largest eigenvalue of A. In this context, the weight is zero two nodes are disconnected, and follows uniform distribution U 5) when nodes are connected. Lastly, we define θ DÔD, where D is a diagonal matrix selected such that the covariance matrix O-1 has unit diagonal elements.
We set the regularization parameter Xıj in Problem as follows
A
Vi
(23)
ii +
where X represents an estimator, c is set to 10-3, and > a parameter that adjusts the sparsity. The function in (23) is closely related to the re- re-weighted -ormr regularization [14 45 (6).. which effectively enhances the sparsity of the solution and reduces the estimation bias resulting from the f-normm 477.. We employ the maximum likelihood estimator as X, i.e. the minimizer of Problem (T) without the sparsity regularization. Note that one can solve the maximum likelihood estimator with relatively large tolerance obtain coarse estimator, and alternative monotonically decreasing functions may be explored in 23 For PQN-LBFGS, we utilize the previous 50 updates to compute the search direction. For FPN, we set 6 = 10 15 in K for identifying the set of restricted variables.
For calculating the relative error as defined in 20 there is no additional computational cost for PGD, APGD, PON-LBFGS and FPN, since the objective function value is already evaluated within the backtracking line search during each iteration. In contrast, BCD, optGI., and GGL. not require evaluating the objective function, leading extra computational cost when computing the relative eтror in (20) for comparison purposes. However, this additional cost can be considered negligible, as these methods only compute the objective function value after completing a cycle. Moreover, the number of cycles needed for BCD, optGL, and GGL is substantially fewer than the number of iterations required by other methods.
14

---

A.2 Financial time-series data
We carry out numerical experiments on a financial time-series dataset to evaluate the performance of method in recovering graph edges. The applicability of MTP2 models financial time-series data is well-established,
market factors lead to positive dependencies among stocks [
The dataset consists of 201 stocks composing the S&P 500 index, spanning the period from January 2017, to January 1, 2020, yielding 753 observations per stock, i.e., 201 and 753. We construct the log-returns data matrix as Xi,j = log P.j - log P.-1.j- where Pt, represents the closing price of the j-th stock on the i-th day. The stocks are categorized into five sectors based on the Global Industry Classification Standard (GICS) system: Consumer Staples, Utilities, Industrials, Information Technology, and Energy.
Directly assessing the correctness the learned graph edges is not feasible in financial time-series data, the underlying graph structure remains unknown. However, expect stocks from the same sector to have interconnected edges. To measure the performance of edge recovery, we employ the modularity metric [48]. Given graph G (V, E), where V represents the vertex and E denotes the edge set, the modularity is defined as:
d,d3 2|
Σ JEE
Au
b(c,,))
(24)
2
where Ai = if (i.j) E E,, and O otherwise. d, represents the number edges connected to node i, c; indicates the type node and 6(:,)) refers the Kronecker delta function, with 5(a, b) = if a and otherwise. A stock graph with high modularity exhibits dense connections among stocks within the same sector and sparse connections between stocks in distinct sectors. A higher modularity value implies a more faithful representation of the underlying stock network.
(a) Glasso, Q = 0.47
(b) FPN, Q = 0.65
(c) FPN(E), Q = 0.67
Figure Financial time-series graphs learned via (a) Glasso, (b) FPN, and (c) FPN(E).
Figure demonstrates that the performances of our proposed FPN and FPN(E) are better than that of Glasso 22 since the majority of connections in graphs learned through FPN and FPN(E) occur between nodes within the same sector. In contrast, only a few connections (depicted as gray-colored edges) exist between nodes from different sectors. Both FPN and FPN(E) achieve higher modularity values compared Glasso, indicating that the former have higher degree of interpretability than the latter. Furthermore, we observe that FPN(E) moderately enhances the performance of FPN.
We fine-tune the sparsity regularization parameter for each method based on the modularity value, allowing only limited number of isolated nodes. Note that increasing the regularization parameter for Glasso would result in numerous isolated nodes that cannot be grouped. FPN(E) refers to the application of FPN for solving Problem T with disconnectivity set E. This set is obtained through hard thresholding the MLE, which is also used in computing regularization weights in (23
15

---

B
Proofs
this section, we provide the proof' for Propositions 3.1 3.3 and 4.5
3.4
and 3.5 as well as Theorems 4.1
B.1
Proof of Proposition 3.1
Proof. The Hessian matrix at X has the form:
H₂ Xx Xx
Then we obtain
= X Xk.
Following from the property of the Kronecker product that vee (ABC) (CT & A vec (), can obtain
vec (Px) = H vec (Vf(X))) = vec (X (Xx)))
As result, we have
Ps X+
completing the proof.
B.2 Proof of Proposition 3.3
Proof. Following from 10 and we obtain
P
T H
(Xx)) [H 'vee (Prg (VI(X))))
where projection Pz; is defined in 12 [A]Tg E I2 and [vec (A)]zz E
For ease of presentation, by slight abuse of notation, both ! represent a vector containing all elements of the set I
Following from the fact that vec (ABC) = (CT & A) vee B) and H X 2 Xg. we have vec (Pzg (V](())))
(X&PIE (VS (Xx)) Xx)-
By collecting the elements in the set IG, we obtain
P PL = XxP (Vf X xz
completing the proof.
B.3
Proof of Propositio 3.4
To prove Propositio 3.4 we first establish Lemma B.1 below to show that the lower level set of the objective function is compact. We note that Lemma B. depends on the condition that the sample covariance matrix has strictly positive diagonal elements, which is assumed throughout the paper and holds with probability one.
Lemma B.1. The lower level set L, defined in 17 is nonempty and compact, and for any X E L,, we have
mI
X MI,
where and M two positive scalars.
Proof. We first show that the largest eigenvalue of any X E Ly can be upper bounded by M. The objective function of [ can be written as
S(X)
log (eeddee((xi)) x)+ (XG),
where G = S 5-2 A with À defined by
{
Aas if
[Als
(25)
0 othrewise.
16

---

Define two constants
# = min Gas and
v =
max Gigg
ig (3))EE
Note that μ > 0 because Siz= > 0 for each E [p]. Then one has
tr (xG) Σ Ga + Eij 2 ptr (X) + Σ Xij.
(26)
(#)
where the inequality follows from the fact that Xy for # j, and X2z 0 for (i.j) E. We denote the largest eigenvalue of X by Amax (X). Then we have
1 (m((G)) + p
Amax () S (t
(27)
where the second inequality follows from 26 In what follows, we bound the terms (XG) and Σizj Xiy. For any X E Ly, has f(X*)
2 log det( (X) tr (XG). Therefore, the term
(XG) can be bounded by
tr (xG) S f(x)) + log det x)= S f(xx) + plog (max(())
(28)
Let X' == min (1))E Aij. Then one has
Ei E A
1 x (f(x)) + plog (Amax (x)) ),
(29)
x
where the last inequality follows from Σ i Aigl S tr (XG). because (XS) > since
X E SP and S E S.. Together with (27).
(28) and 29 we obtain
Amax x) x) S 11 (1+ ()(X)) plog (Amax (X)))
Since log (Amax (X)) grows much slower than Amax Amax (X) can be upper bounded by a constant M. which depends on f(X"). S, and A.
We denote the smallest eigenvalue of X by Amin (X). For any X E Ly, one has
f(X²) 2 f(X) log det( 2 log Amia (X) (p 1) log M.
As a result, we have Amin (X) 2 - 1(x)) M (p-1), which shows that Amin (X) can be lower bounded by a positive constant m = e ((x")) M M-(--1) Finally, we show that the lower level set L is compact. First, L, is closed because is a continuous function. Second, b is also bounded because it is subset of {X RPxP mI = X = MT).
Proof of Propositio 3.4 We first show that (*())) E S holds for a small enough step size. The (x(w)) can be equivalently written in the form (1())) = Zy - Gr where Z₃ E RP*P defined by [Za &ali if (i, Ig. and 0 otherwise. G3 is some search direction only over Ig. which is symmetric matrix. It is known that a matrix A is nonsingular M-matrix if and only if A is Z-matrix and there exists a vector æ > 0 with Ax > [49]. Following from the fact that X₃ is nonsingular M-matrix, there exists > 0 such that
Z 2 X 0.
Therefore, Zμ also a nonsingular M-matrix, implying that Z₁ E S.+. As a result, if 7k < Amin Zz)) where ((Gx) is the spectral radius of Gk, then Xx(T)) E S.+. The definition of X Xe(%k) further guarantees that X((1k) ) E it. We can verify that f(X) + for any X e cl(U)) U", where X is positive semidefinite and singular. Thus, we consider an X E UP with f(X)) sufficiently large such that fxxx (7k)) < fX)).. Therefore, X²(*k) E Ly.
Next we prove that the line search condition 15 holds for a small enough step size. Recall that Zţ is the complement of the set Ig defined in For any X₂ E T.g, the set Ie can represented as
5
= T 3=
U B
1=1
17

---

where 0, . 1,....5,. _,, are defined
8/ = {(1,5) E E- - < 1 [x] 50, [vf(X))) 1, 20, [Pal, <0}.
[ - {(.5) E c S [Xn]ys S [vf(X [)) 20, [Palyy 2 0}.
8/ {(.3) 89 [ Ekk [P < <0},
_ = {(.,)) e 80 [Xx]es <kk
[Pt] 2 0),
B, ((i,)) 501 [X*],,] 0}.
Each subset 5) disjoint with each other. of Xg. whiles the other subsets include the indexes of the off-diagonal elements. Then one has, for
B%() only contains the indexes of the diagonal elements
any (i h and N O,
0 s [x*(T)]) [xa]is [Pn
7/ Pk) [Xx].] I Y>k [Paly
(30)
For any (i E if the step size satisfies
2 and 7ik > O, one has [xx(>>)) - -xxx);; -> [Polej- aorany (i,j) E 13())
Ek
0 9h -
min (1J)
(31)
[
then one has [Xx(Tx)] [xa]e] 7 Ph is' Similarly, for any (i,5) E B and th
[x (mx)) - [xx] ~ [x)]i Finally, for any (i.j) ( [Xx], iy must be on the diagonal
of X₁. Therefore, we can directly remove the projection Pa. and obtain [Xx(1e)] [xa]
-7h [p Therefore, if (31 holds, one has
([Vf(x))] z [X*(T>)] [Xa] (xx]
- (VF((x))) [Px]zz)
(32)
where the inequality follows from 30) Recall that [X*(1u)]z, = O, which can be equivalently expressed as the form of projected gradient descent:
[x((1)) T = -Pa (1xx)) - Tk De @ VS(X
(33)
where we write Pa for [Po(A)] with slight abuse of notation, and D₁ is defined as
[Dxlig = [x([V)(X))ls))
Ck
(i.j)) Is \E
Then one has
X (x()n) -xx 1
<:x PPx Z
[Ph] zzz
7/6 D @ vf(X+) [Xa(T)) -[Xu]za) I
(34)
< [Px] Tz2
[Px]
+ as [V](()))] I [xx] )m))
where ag = f[xx)( TeXe min
Ck
and To represents T(.. Note that In = Ts U and
[Kaly = (1x)ls = 0 for any (i.j) E
It is worth mentioning that we primarily focus cases with non-empty Tr. In the rare event that Th is empty, such as when To is empty (in this situation we define e = 6. leading to an empty Tt), we can simply ignore the terms related to I₁ in (34) and the following proof.
Furthermore, one has
< Ph ze.
Ph xz)
- [|
12211 M
[V((x))]
IE 11?
zEZE
(35)
= >maxx
[M,')
(zşzz!
([V
((X)()))
[Pk Iz).
18

---

The largest eigenvalue of [Ma' Izzzz can be bounded by
Amax (LM
Amax H
Amax (X Xx) M²,
(36)
where the first inequality follows from the Eigenvalue Interlacing Theorem, and the second inequality follows from Lemma
B.1 Combining
34
35
and (36) yields
[x(Tk) - Xx l ( [po]Tz) + ax (V)(()) Tx [xx]
Since (x(7)) E Lyg, we obtain
I ((xx(7)))
xx(*k) Xx) + 2m?
(x((Tk) Xx l
(37)
where the inequality follows from
mmax V f(x)) Amax
(x
1) X 2max (x (x-1) S m²
AX
L,,
where the inequality follows from Lemma B.I
Let 7n = min 2( 2(1-a)mi a)m²
Ek
Amin (Zs)
min (
Note that 7 is bounded away from
M2
[Plla
p(Gx)
zero, because each [Piley is bounded following from Lemma
B.
To this end, for any E (0,
(0))..)
((Xx(T))) - (07 z
Pk
x a [<(V]((x)))
Iu [ I
(38)
where the inequality follows from [[v
32 the definition of <s in
2 the inequality
(xxx)(X))
Ta \ min" because Tk C Ts. and the inequality [M
min
([Vf(X)) z
Ph r
2 0. because
is positive definite and Izzg
(xxx xz) [Polag)
([V
[Mr I4
[Vf(x))
Ig)
(39)
completing the proof.
B.4 Proof of Proposition 3.5
Proof. The line search condition 15 leads to
f (Xx(7a)) f (X < - on [V[(x))] (
[M,'] [M''ZzTT ĀØ
(40)
where the inequality follows from
([V](())) [X&]
[xx]
Is )
200
We further have
Amin
[M,
Amin (H 1 Amin (X & Xx)
2 m
(41)
where the first and second inequalities follow from Eigenvalue Interlacing Theorem and Lemma B. respectively. We complete the proof by combining 40 and (41)
B.5 Proof of Theorem 4.1
Proof. We first prove that Problem T has at least one minimizer. It is known by the Weierstrass' extreme value theorem (50] that the of minima is nonempty for any lower semicontinuous function with nonempty compact lower level set. Therefore, the existence of the minimizers for Problem (I) can be guaranteed by Lemma B.I
On the other hand, we show that Problem I has at most one minimizer by its strict convexity. We have f(x X @ X- Thus 2 f2f f(X) > 0 for any X in the feasible region of Problem
19

---

① defined in (T).. Therefore, Problem 1 Together with the existence of the minimizers, we conclude that Problem (T) has an unique minimizer.
is strictly convex, and thus has at most one minimizer.
The Lagrangian of Problem | is
C (X, = logdet tr + (Y X),,
(42)
where Y is KKT multiplier with Yú, = for i E [p]. and À is defined in (25).
For convex optimization problem with Slater's condition holding, pair is primal and dual optimal if and only if the KKT conditions hold. Thus, (x*, YY) is primal and dual optimal if and only if it
satisfies the KKT conditions of (T) as below,
(X - - 0;
(43)
x1-= 0, V
(44)
x = 0, X = 0, G 2 0, Vi * and (i.j) $ E:
(45)
Y = 0, ViE [p)
(46)
We first prove that the minimizer X must satisfy all conditions in 18 Note that the KKT conditions
=0. First, X = 0 for any (i,5) & following from 44 Second, for any (i.j) with * j and (i.j) E V we have X + 0
(43)
46) must hold for the minimizer X Let V = { (i. E [p]
X
and (i,j) & E. Following from (45) we further obtain Yū = O. Together with (46). we conclude that
Y13 0 for any (i V Following from 443
V f Ve + j. thus obtain Yi 0 according
O. Since X is positive definite,
(i, i) E v for any i [p]. Then, for any (i, E VIE, we have
to 445 Following from (43) we get [Vf(x*) ViE I O.
Now we prove that any point X E MP satisfying the conditions in (8) must be the minimizer, i.e., the KKT conditions (43) 46) hold for X We construct Y by Y =VV[((())) First, is
straightforward to check that the conditions (43 and (44) hold. Second, we have
[Y']y- = VV }]v O.
(47)
We know that (i, E vo for any i E [p], since X E M Together with (47). we obtain that the condition (46) holds.
Finally, following from the fact that (i, i) E V for any i E [p] and E C V, we obtain
{(i E 1 (i.5) $ E) II
where I = {(,,5 E |p] (i.j) E V" i f j}. and Is = {i.,) E Wp| E. & E}. For
any (i,j) E Ii, we have Y = 0 according to 447), and X;; < <0 since X E MP. Thus the condition (45) holds for any (i.j) E Z, For any (i,j) Iz, we have X; = 0. and Y3 --V )[s 2
according to the second condition in (18) Thus the condition 455 also holds for any (i,j) E Iz.
Totally, the condition (45) holds. To sum up, all KKT conditions (43 (46) hold for any X E MP
satisfying the conditions
and thus we conclude that X * is the minimizer of Problem TT.
B.6 Proof of Theorem 4.2
Proof. Following from Proposition 3.4, for any iterate X₃ E Ly, a small enough step size ensures that the line search condition (15) holds, which leads a sufficient decrease of the objective function as shown in Theorem 3.5, and the next iterate X6+1 E Ly- When the sequence starts with Xo E Ly. each point of the sequence (Xx) admits X E T.y. Note that it is easy to construct an initial point Xo E Ly. because we could consider an X E LP in (I7), which is close to cl(UP) \ UP such that f(X") sufficiently large, following from the fact that f(() = +oo for any X E cl(UP) APP Lemma B.1 shows that the lower level set L.y is compact, thus the sequence (Xx) has at least one limit point. For every limit X of the sequence { {X&), we have X E MP. Define I' = T(X* e UE, where T(X*, (X*,e") ,,)) is equal to
T(X* (()==')= = ((.5) E 1p| - c < [X SO, [VI(''), <0}.
20

---

and c" is defined as e* := min (2(1 - a ammm [V TyVe min' 6), where Tg denotes T(X*; (X*:6) .6.
Since f is continuous and {f(X))) keeps decreasing and is bounded, it follows that {
converges and
lim f(Xx) f(Xx+1)=0 =
Following from Proposition 3.4 the line search condition a strictly positive "Yk it holds
15 ensures that each X E S4 and, for
f(Xx) - f(X++1) 2 on
[Pllzy + (V x (Xalza) [M E
(48)
an [Vf
Vf(Xx)) xz + a (VV [Xxlz,]
Recall that Ix = -.. We note that the following facts hold: a E (0, 1) is constant, [ zçig is positive definite, lis = for any (i.j) E E. [x Σ 0 for any (i.j) E Is (E, and [ < O for any (i.j) E Za E due to the definition of T(X,e) Thus, the two terms on the right-hand side of (48) are both nonnegative. Moreover, the right-hand side of (48) approaches zero and only both [V](x)) and [Xa]z, simultaneously go to zero. Therefore, we deduce that every limit point X must satisfy:
[V( = 0, and X z- = 0.
(49)
We show that every limit point X + is the minimizer of Problem ① according Theorem 4.1 Let
V = (i.j) E [r] [x] ig = 0}. First, for any (1,5) E &, we have xx] = 0, because of the
projection Pn in each iteration. Note that [x+] = 0 for any (i.j) E I according to (49 For any
(i,j) VV, ie., [x+] * O, we must have (i, E {T* Together with [V((x)]
holds. For any (i.j) E VİE, we must have
E ( (')
Recall that [V((*)) < 0 for any (1,5) E T(X*, (*), and VV(X)), = 0 for any (t,j) E {I*]S. Overall, we obtain
VV Vig 2 O.
To sum up, all the conditions in Theorem 4.1 hold for every limit point X and thus every limit point is the minimizer of Problem (I]
Since the minimizer of Problem (1 is unique, we obtain that the limit point of the sequence {Xx] is also unique, and thus } is convergent. Therefore, we conclude that the sequence {Xa) converges the unique minimizer of Problem (T). The monotone decreasing of the sequence (f(X))) can be
established by Proposition 3.5
B.7 Proof of Theorem 4.4
Proof. We prove that the support of X is consistent with the set Zg for sufficiently large k. Without loss of generality, we specify the constant 5 (9) as
min (i.j) (su))sppp((')) (
k
(50)
where w. E (0, 1) is a constant. Note that sequence {Xx} converges to X Under Theorem 4.3 that [V[(X))] is strictly negative for any (i,j) (X+) E, there must exist some and Kj E N4 such that
[F((x)),,
(i,j) e supp (x')
(51)
holds for any le 2 K1, where a (0, 1) is constant, and m is a positive constant defined in Lemma B. We consider a neighbourhood of X defined by
<( = (X E RPP X - X ST).
(52)
21

---

where r is a positive constant defined
min (c
|[ 18b a, 8). ((p))sm)p(())
min
(53)
where 1 w is positive constant. There must exist K2 E N+ such that X E <X(xx; holds for any l: 2 K2. Take K, max max(Ki, K2). For any h 2 K. and (i.j) supp(X )*), one has
[[x*] - [ (g) S [[x]
Xx - X Sr. [x
Thus one can obtain, for any & 2 Ko.
X
min (k.)))) x))
| - - 8> O, (i,j) supp(X"),
(54)
where the last incquality follows from and 53 [xx]a]] sy4 S
Recall that for any (i.j) E T(Xx: 8).
Then one has
T(X 6) S supp"(X').
(55)
Then the E in p can be bounded by
2 min (2 a0llaamml|[V[i supp (x o) min(a, 8),
(56)
where the first and second inequalities follow from 55 and (51). respectively. For any Ko and (i,j) E TXx,, U E, one has
[ S c 6,
where the second inequality follows from the definition of Gk- Thus we obtain
T(X k, ()) C supp"(X')
(57)
On the other hand, for any &: 2 K2 and (i,j) E supp (Xp) one has
[ k/ [X+]
[
Xx x 1| S S Ek:
(58)
where the last inequality follows from (53) and (56). Therefore, one has
supp (upp"(X")=e ( Ex) UE U B (
Note that B supp supp"(Xx)) diagonal which must be nonzero. Moreover, following from 51 one has
because any (i,j) E B corresponds the element the
(supp" (x*) (X*) (E) n B(" Ø, and (supp" (oupp"("") (x*) (B( Ø.
Therefore, we can obtain supp
C TT(X (x) U. Together with we obtain
supp "X*) = T(X k. 6n -E= Iz.
(59)
Equivalently, we have
supp X = Ií 4 2 Ko. We can see that the sets I, and Ig are fixed for any k 2 Kn, Therefore, following from the iterate of X
(60)
X I--, = [X+tiz2, O, vk 2 Ko.
Moreover, together with 54 and (60). we obtain that for any A 2 Ko, EXxtIIze,, # 0. Take
k. = K, + 1. We obtain
(xp) Ig, vk ko.
completing the proof.
22

---

B.8 Proof of Theorem 4.5
Proof. Theorem 4.5 is a direct extension of the result in 51 Let g be twice continuously
differentiable and consider the following iterate
xx+1 2g Yk D Volak))
where D is positive definite and symmetric. Then the sequence { admits the following convergence result [5I]
FK++1
xx Dr lim sup max (11-
lim sup 10-000
(11-70m1 7a (111-1-M111
(61)
| x² 11 D.
R-000
where x² is the limit of the sequence {æx), which satisfies that Vg(x=) and g(==)) is positive definite, and myg and M2 are the smallest and largest eigenvalues of (Dz) g(xx) (Dx) respectively. This conclusion is extended from the convergence result for quadratic objective function. Conceptually, this makes function is very close to positive definite quadratic function in the neighborhood of non-singular local minimum.
because twice continuously differentiable objective
Following from Theorem 4.4 for any lo Kos iterate (3) be written
Xx+1 R++
Xx
Th R
(62)
z
which reduces to an iterate of an unconstrained optimization algorithm some subspace. Following from results (TT, we obtain
k+1 X
X
lim sup k-100
I
I
R.
lim sup max (11 (1- (1-
(2 12: -xkMas))
(63)
X
[X*
R k
where and Mg are the smallest and largest eigenvalues of R
H z R
respectively.
Following from Theorem 4.2 and R H
obtain I
T R₁ X
X
1 M
Vk
2 Ko.
(64)
When & 2 ko. the line search condition reduces to
/ ( an
[Px Px]TT)
Similar to the unconstrained case in [52], the step size must satisfy the line search condition if
Th 2 min (1, 2(1 3)3/Mr) as + co. Then together with (63) and (64), we complete the proof.
23


Reviews:
[{"review_id": "PVZzpVdMs9", "reviewer": "Reviewer_AisR", "summary": "The authors propose and analyze a variable metric projected gradient method for a variant of the graphical Lasso problem where the off-diagonal elements are additionally constraint to be non-positive, and in which specific entries of the precision matrix are known to be zero.\nThe proposed metric which is used to precondition the gradients comes from the curvature of the function $X \\mapsto -\\operatorname{log det}(X)$ , which is part of the objective, and thus incorporates second-order information of the objective function, unlike more standard approaches for the problem.\n\nThe authors put their result in the context of the recent literature, and show in numerical experiments on covariances created from multivariate Gaussian whose precision matrix sparsity comes from a Barabási-Albert random graph that the proposed method outperforms other methods used to optimize the same objective function, such as block coordinate descent or the proximal point algorithm in terms of clock-time convergence speed.\n\nFrom a theoretical point of view, the authors show that the proposed algorithm indeed converges to the unique minimizer of the (convex) objective, that the proposed stepsize rule indeed works, and that the algorithm identifies the true support in a finite number of iterations under an additional reasonable assumption on the minimizer.", "questions": "- When discussing the positive definiteness constraint in line 112, you mention that it is not closed and furthermore, cannot be managed by a projection. However, by looking at (1), it seems that due to the objective, it is clear that the minimizer will not be singular due to the existence of the term $- \\operatorname{log det}(X)$  in the minimization: If at least one of the eigenvalues of X is 0, $\\operatorname{log det}(X)$ is finite. Lemma B.1 touches upon this, but couldn't the positive definiteness constraint be relaxed to a psd constraint?\n- In the considered model (1), the role of the disconnectivity constraint remains slightly unclear. Appendix A.1. and also the discussion of the financial time series example leaves the impression that the disconnectivity set is found by a similar method as the weights of the $\\ell_1$-term $\\lambda_{ij}$. Can you make the distinction clearer?\n - In Theorem 4.2, it will be convenient if one the authors can refer definition equation for $L_f$.", "limitations": "The this discussion of limitations is adequate in general. However, it would have been beneficial if the \"hidden - $\\lambda$-dependence\" on the sparsity of the estimated covariance matrix had been explained better in the paper.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "I consider the proposed efficient estimator for the sparse precision matrix with off-diagonal non-negativity constraint of interest for the machine learning community due to wide potential applicability of the setting to engineering problems. \nThe specific gradient preconditioning through the weighted metric induced by the second derivatives of the $\\ell_2$-error is original for the problem, and the computational advantages compared to state-of-the-art methods are well explained ($O(p^3)$ time complexity) and showcased in sufficiently timing experiments. The choice of the comparison methods is adequate and the experimental setups are informative.\n\nThe paper is clearly and well written. A strength of the paper is the suitable co-design that heavily takes into account the problem structure, which makes the Hessian approximation used in the variable metric preconditioning simple enough to compute in a reasonable amount of time, while at the same time delivering a benefit in per-iteration improvement.", "weaknesses": "While the problem as formulated in (1) is convex as stated, the underlying estimator pipeline is somewhat inherently non-convex due to the fact that the $\\lambda_{ij}$ correspond to weights of a reweighted $\\ell_1$-term whose weights are chosen based on the MLE. This detail could be mentioned more prominently in the paper. \n\nFrom a computational point of view, the line search method deployed has the computational disadvantage that a lot of function evaluations, which might be expensive due to their need for spectral decompositions, might be needed to make the method work.\n\nIt remains unclear whether the convergence rate result Theorem 4.5 is optimal. Due to the incorporation of second-order information, a discussion of why superlinear convergence cannot be shown or observed would be further informative.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> 1. From a computational point of view, the line search method deployed has the computational disadvantage that a lot of function evaluations, which might be expensive due to their need for spectral decompositions, might be needed to make the method work.\n\n**Reply:** Thank you for your perceptive observations. We acknowledge the potential computational intensity of the line search procedure, due to the need to compute the log-determinant function and verify positive definiteness. We genuinely appreciate your attention to this critical aspect.\n\nIn our current implementation, we first conduct the Cholesky decomposition $\\mathit{X} = \\mathit{L} \\mathit{L}^\\top$ using MATLAB's \"chol\" function. This function can simultaneously verify the positive definiteness of $\\mathit{X}$. Next, we calculate the log-determinant function as $\\log \\det (\\mathit{X}) = 2 \\sum_i \\log ( L_{ii} )$. The Cholesky decomposition is the most computationally demanding step, generally requiring $O(p^3)$ costs for a $p \\times p$ matrix. However, the sparsity of the matrix can substantially reduce the computational cost in MATLAB, as the \"chol\" function utilizes CHOLMOD. This function can process a sparse positive definite matrix of 10,000 dimensions within approximately 2 seconds. Nevertheless, the computational load could rise significantly as we scale to higher dimensions.\n\nTo address the computational challenge, we draw attention to a more efficient method for evaluating the log-determinant function and verifying positive definiteness, as presented in [R1]. This method, leveraging Schur complements and sparse linear system solving, can tackle problems of up to 1,000,000 dimensions. Furthermore, it's worthwhile to investigate more efficient strategies for computing an approximate log-determinant function. In this context, the approach proposed in [R2] ensures a nearly linear scaling of execution time with the number of non-zero entries, while maintaining a high level of accuracy.\n\nThis discussion will be integrated into our revised manuscript.\n\nReferences:\n\n[R1] C.-J. Hsieh, M. A. Sustik, I. S. Dhillon, P. K. Ravikumar, and R. Poldrack. \"BIG & QUIC: Sparse inverse covariance estimation for a million variables.\" Advances in neural information processing systems, 2013.\n\n[R2] I. Han, D. Malioutov, and J. Shin. \"Large-scale log-determinant computation through stochastic Chebyshev expansions.\" International Conference on Machine Learning, 37:908-917, 2015.\n\n> 2. It remains unclear whether the convergence rate result Theorem 4.5 is optimal. A discussion of why superlinear convergence cannot be shown or observed would be further informative.\n\n**Reply:** We appreciate your insightful comment. We would like to clarify that our convergence rate result is not optimal, as our algorithm generally does not attain superlinear convergence, despite the incorporation of second-order information. Superlinear convergence requires the inverse gradient scaling matrix to progressively approximate the Hessian at the minimizer, a condition that our algorithm does not satisfy. \n\nHowever, constructing a search direction to achieve superlinear convergence is significantly more computationally demanding than our method, as it cannot leverage the special structure of the Hessian to decrease the computational burden.\n\nWe will include this discussion in the revised manuscript for further clarity.\n\n> 3. Couldn't the positive definiteness constraint be relaxed to a PSD constraint?\n\n**Reply:** Thank you for your insightful comment. We acknowledge your suggestion that the positive definiteness constraint can be relaxed to a positive semi-definiteness (PSD) constraint. This would indeed create a closed set that could be managed through projection.\n\nHowever, we chose not to follow this route due to the computational burden associated with projecting onto the intersection of the two sets - corresponding to the PSD and non-positive constraints. Implementing such a procedure, for example via Dykstra's projection algorithm, would involve alternating projection onto each set until convergence is achieved. This would significantly increase the computational load compared to our current method, which ensures positive definiteness through a line search method and manages the non-positive constraint via projection.\n\nWe will enhance the revised manuscript with a more detailed explanation of this aspect. \n\n> 4. In the considered model (1), the role of the disconnectivity constraint remains slightly unclear. Can you make the distinction clearer?\n\n**Reply:** Thank you for your valuable feedback. Both the disconnectivity set and the $\\ell_1$-term weights are derived based on the Maximum Likelihood Estimation (MLE) in our experiments, but they are procured through distinct techniques and serve different roles. The disconnectivity set is identified by applying hard thresholding to the MLE, and it is used to incorporate prior insights about the graph structure. On the other hand, the $\\ell_1$-term weights, which are derived from a monotonically decreasing function of the MLE, are designed to promote sparsity in the solution.\n\nIt's worth noting that the disconnectivity set and the $\\ell_1$-term weights can be determined through various methodologies. In addition to the MLE, other estimators can also be effectively employed. For instance, one could use the thresholded sample covariance matrix supported by various thresholding techniques. Moreover, the disconnectivity set can also be obtained in the context of learning structured graphs, such as bipartite graphs, by incorporating inherent disconnection in graphs.\n\nWe will provide a more detailed distinction between these two aspects in the revised manuscript.\n\n> 5. In Theorem 4.2, it will be convenient if the authors can refer definition equation for $L_f$.\n\n**Reply:** Thank you for your helpful suggestion. We will refer definition equation for $L_f$ in the revised manuscript to make it more convenient for readers.\n"}, {"review_id": "t1AXysYSNZ", "reviewer": "Reviewer_7MEo", "summary": "(1) This paper investigates the problem of estimating precision matrices in Gaussian distributions  that are multivariate totally positive of order two (MTP2), and propose an algorithm based on the two-metric projection method.\n(2) Empirical tests are conducted. ", "questions": "Please address these questions I posed to in the previous message.", "limitations": "Yes, they have.", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "- **Originality:** The techniques used in this paper reduces computational complexity and its theoretical convergence is established.\n \n-  **Quality:** The notations, problem statements, and mathematical proofs appear rigorous, although I have not scrutinized them word by word.\n \n- **Clarity:** The readability of this article is good.\n \n- **Significance:** The techniques (and tricks) may inspire and guide future researchers.", "weaknesses": "1. The main issue with this paper is that it does not cite the highly relevant paper [A]. The bound provided by theoretical analysis of this paper is a supremum of a limit point, which generally can not guarantee the convergence, while the theoretical result of [A] is much more attractive. \n\n2. Also, in the empirical study, the algorithm presented in the paper has not been compared to [A].  \n\n[A] Ying, Jiaxi, José Vinícius De Miranda Cardoso, and Daniel P. Palomar. \"Adaptive Estimation of Graphical Models under Total Positivity.\" (2023). ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> 1a. The main issue with this paper is that it does not cite the highly relevant paper [A].\n\n> [A] J. Ying, J. V. D. M. Cardoso, and D. P. Palomar. \"Adaptive Estimation of Graphical Models under Total Positivity.\" ICML, 2023.\n\n**Reply:** We sincerely appreciate your insightful feedback. We apologize for our initial oversight in not citing this relevant work—this was definitely not our intention to overlook its contributions to the field.\n\n**It's worth noting that paper [A] appeared after an earlier version of our work had been posted on arXiv and it cited this earlier version.**\n\nWe would like to clarify that our work and paper [A] have different emphases and distinct methodologies, although both aim to advance the field of MTP2 graphical models.\n\n**Paper [A] primarily focuses on the enhancement of precision matrix estimation and graph edge recovery**, prioritizing these aspects over computational efficiency. The authors present a multiple-stage estimation method that progressively refines the estimate at each stage by solving a constrained log-determinant program, using a projected gradient descent (PGD) algorithm. However, as a first-order algorithm, PGD is often hindered by computational inefficiency due to its slow convergence rate.\n\nConversely, **our work emphasizes enhancing computational efficiency, recognizing this as a crucial aspect when dealing with large-scale problems.** Our paper is dedicated to the development of an efficient and scalable algorithm tailored for estimating large-scale MTP2 graphical models. In comparison with the PGD algorithm used in paper [A], our proposed second-order algorithm achieves a significantly faster convergence rate, while maintaining the same per-iteration computational complexity. Our extensive experiments demonstrate that our algorithm significantly boosts computational efficiency compared to state-of-the-art methods, including PGD.\n\n**We provide a more detailed comparison of the two papers in our “global” response.** In the revised manuscript, we will ensure to appropriately cite paper [A] and clearly articulate its relevance to our work. We appreciate your constructive feedback, which helps enhance the clarity and thoroughness of our paper.\n\n> 1b. The bound provided by theoretical analysis of this paper is a supremum of a limit point, which generally can not guarantee the convergence, while the theoretical result of [A] is much more attractive.\n\n**Reply:** Thank you for your insightful feedback. In terms of the theoretical results of the two papers, we would like to provide the following clarifications. **Our theoretical analysis shows that the proposed algorithm is guaranteed to converge to the global minimizer**, with the set of free variables converging to the support of the global minimizer within finite iterations, as declared in Theorem 4.2 and Theorem 4.4, respectively. We also establish the sequence convergence rate in Theorem 4.5. \n\nIn contrast, paper [A] emphasizes enhancing estimation performance. Their Theorem 4.4 characterizes the estimation error between their proposed estimator and the underlying precision matrix. However, **they fail to provide a guarantee for the algorithm's convergence to a global minimizer and do not establish the sequence convergence rate.** In fact, the formulation in paper [A] is nonconvex, and their algorithm only guarantees convergence to a stationary point. We will include a more detailed comparison in the revised manuscript.\n\n> 2. Also, in the empirical study, the algorithm presented in the paper has not been compared to [A].\n\n**Reply:** Thank you for your insightful comment. Our empirical studies primarily focus on comparing the runtime of various algorithms in terms of their convergence to the global minimizer of a common convex log-determinant program, as illustrated in Figures 2-5 of our paper. **We clarify that a direct comparison with paper [A] is not feasible**, as its algorithm reaches a stationary point of its nonconvex formulation problem, which is fundamentally different from our approach.\n\nHowever, **we have made an indirect comparison with paper [A] in our first submission through our comparison with the PGD method**, which paper [A] employs to solve a log-determinant program at each stage. Our extensive experiments demonstrate that our algorithm significantly improves computational efficiency compared to PGD. Moreover, since paper [A]'s method requires the repeated use of PGD across different stages, presents a considerable advantage in terms of computational efficiency.\n\nTo further address your concern, **we provide below an additional comparison with paper [A] in the experiment on financial time-series data** (which is presented in Appendix A). We compared both the graph edge recovery performance and runtime, using the modularity metric for edge recovery, where a higher value indicates a more accurate graph representation. \n\nThe table below provides the modularity values and runtime for the method in paper [A], FPN, and FPN($\\mathcal{E}$), where the latter two are proposed in our paper. The method in paper [A] leads to a higher modularity value than FPN due to its iterative refinement strategy. FPN($\\mathcal{E}$), which imposes a disconnected set, yields a modularity value comparable to that of the method in [A]. In terms of computational time, both FPN and FPN($\\mathcal{E}$) significantly outperform the method in paper [A], due to the relative inefficiency of the PGD used in [A] and the repeated need to solve the log-determinant program.\n\n| Methods |  paper [A]   |  FPN  |  FPN($\\mathcal{E}$)  |\n|-------|------|------|------|\n| **Modularity** | 0.667 | 0.646 | 0.672 |\n| **Runtime** | 215.57 | 6.80 |  10.17 |  \n\n**The graphs learned via the three methods are displayed in the attached PDF.**\n\n We will include this detailed explanation and additional experimental results in our revised manuscript. We appreciate your feedback, which will undoubtedly enhance the quality of our paper.\n"}, {"review_id": "mHvcSHb3eL", "reviewer": "Reviewer_KQq4", "summary": "This paper presents a novel projected Newton-like method for computing the inverse covariance matrix of a multivariate Gaussian distribution, for which the off-diagonal elements are non-positive. The method is shown to be an order of magnitude faster than existing state-of-the-art methods for this task, while the authors also prove convergence to the global optimum. ", "questions": "Please see above.", "limitations": "-", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "Well-written.\nFast algorithm. \nComplete theory. \nConvincing experiments. ", "weaknesses": "I would like the authors to spend a little more time motivating and introducing the problem.\nIn particular, discussing the following questions on the level of introduction, might make the paper more appealing.\n1. What features of the data lead to such an inverse covariance matrix structure? \n2. What is the exact \"reduction on the number of observations\" for a model with the MTP_2 structure? (line 23). \n3. What is the role, in terms of data / applications, of the disconnectivity set Epsilon? \n4. Why is problem (1) preferred to other alternatives? \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> I would like the authors to spend a little more time motivating and introducing the problem. In particular, discussing the following questions on the level of introduction, might make the paper more appealing.\n\n**Reply:** Thank you for your valuable feedback and constructive suggestions. We appreciate your recommendation to provide a more detailed introduction to motivate the problem and provide context. Accordingly, we commit to expanding our introduction in the revised manuscript to address each of your suggested questions, ultimately making the paper more appealing and comprehensive.\n\n> 1. What features of the data lead to such an inverse covariance matrix structure?\n\n**Reply:** Thank you for your insightful query. The structure of the inverse covariance matrix is primarily determined by the inherent relationships within the data. In particular, conditional independence among variables is represented as zeroes in this matrix. This structure can be seen in data with inherent sparsity, where only a small subset of variables directly influences a given variable. In the revised manuscript, we will elaborate on how the nature of the data leads to such a structure. \n\nWe will provide a more detailed explanation of this in the revised manuscript.\n\n> 2. What is the exact \"reduction on the number of observations\" for a model with the $\\mathrm{MTP}_2$ structure? (line 23).\n\n**Reply:** Thank you for your pertinent query. In the context of a Gaussian distribution, the $\\mathrm{MTP}_2$ structure leads to a significant reduction in the number of observations required for the existence of a maximum likelihood estimator (MLE). Specifically, the requirement reduces from $n \\geq p$ to $n \\geq 2$ [R1, R2], where $n$ and $p$ represent the number of observations and the problem dimension, respectively. \n\nIn the case of Ising models, the MLE can exist with only $n = p$ observations under $\\mathrm{MTP}_2$ constraints [R3], which is a substantial reduction from the $2^p$ observations required without these constraints. \n\nWe will provide a more detailed explanation of this aspect in the revised manuscript.\n\nReferences:\n\n[R1] M. Slawski, and M. Hein. \"Estimation of positive definite M-matrices and structure learning for attractive Gaussian Markov random fields.\" Linear Algebra and its Applications, 473: 145-179, 2015.\n\n[R2] S. Lauritzen, C. Uhler, and P. Zwiernik. \"Maximum likelihood estimation in Gaussian models under total positivity.\" Annals of Statistics, 47(4): 1835-1863, 2019.\n\n[R3] S. Lauritzen, C. Uhler, and P. Zwiernik. \"Total positivity in exponential families with application to binary variables.\" Annals of Statistics, 49(3): 1436-1459, 2021.\n\n> 3. What is the role, in terms of data / applications, of the disconnectivity set Epsilon?\n\n**Reply:** Thank you for your insightful question. The disconnectivity set plays a crucial role in integrating prior insights about the structure of the model. It enables us to incorporate pre-existing knowledge about inherent disconnections in the data, especially in cases where certain variables are known to be conditionally independent. \n\nFurthermore, it is particularly significant in learning structured graphs, such as bipartite graphs, where nodes are partitioned into two disjoint sets, leading to inherent disconnection within each set.\n\nWe will include expanded discussions on the role of the disconnectivity set in our revised manuscript.\n\n> 4. Why is problem (1) preferred to other alternatives?\n\n**Reply:** Thank you for your pertinent question. Problem (1) is preferred due to its unique characteristics and benefits. It is designed to minimize the Gaussian negative log-likelihood, regularized by the weighted $\\ell_1$-norm of the off-diagonal entries of the precision matrix. As a result, we obtain a log-determinant program that can be solved in polynomial time using various methods. This estimator has been well-studied regarding its behavior in high-dimensional scenarios, particularly when dealing with multivariate Gaussian data. In this context, the non-zero structure of the precision matrix encapsulates conditional dependence. Interestingly, this estimator remains robust even for non-Gaussian data, as it corresponds to minimizing a regularized log-determinant Bregman divergence.\n\nMoreover, Problem (1) stands out due to its ability to incorporate the disconnectivity constraint and sign constraint. The disconnectivity constraint provides a mechanism for integrating prior knowledge about the graph structure and facilitates learning structured graphs. On the other hand, the sign constraint contributes to reducing the sample size necessary for estimation, hence enhancing the effectiveness of the method.\n\nIn the revised manuscript, we will delve into these aspects further to enhance the clarity and readability of our work. We appreciate your insightful question, which will undoubtedly improve the overall quality of our paper.\n\n"}, {"review_id": "zHQIjhKgOU", "reviewer": "Reviewer_oo93", "summary": "This paper proposed a novel algorithm for solving the problem of estimating precision matrices in Gaussian distributions with MTP2 constraint. The problem formulation the paper considered is a sign-constrained log-determinant program with disconnectivity constraints. The algorithm is based on a carefully designed search direction and variable partitioning scheme. The author provide the choices of approximate Newton direction and the computing step size with the theoretical convergence analysis, where the author finally shows the convergence in finite iterations. The experimental results indicates better computational efficiency of the proposed algorithm compared to the state-of-the-art methods.", "questions": "See weakness.", "limitations": "The authors have adequately addressed the limitations.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well written and clear, the approach is well supported by the theoretical analysis. The experiment and application seem strong and nice. The authors also made comparisons to other state-of-the-art methods that makes the paper more complete.", "weaknesses": "- Some quantities could be elaborated in the main paper. e.g., $D_k$ in line 153. I am not sure if the form of the diagonal matrix is important or a minor issue that will only be useful for theoretical analysis.\n\n- It might also be good to mention the special case of the algorithm without the disconnectivity constraint and how it compares with other existing methods.\n\n- It seems to me that the choice of regularization parameter $\\lambda_{ij}$ doesn't affect the results in Theorem 4.4 & 4.5. Could you elaborate more on this?\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> Strengths:\n> The paper is well written and clear, the approach is well supported by the theoretical analysis. The experiment and application seem strong and nice. The authors also made comparisons to other state-of-the-art methods that makes the paper more complete.\n\n**Reply:** Thank you for your kind words and for acknowledging the clarity, theoretical support, experimental strength, and comparisons made in our paper. We are grateful for your time and efforts in reviewing our paper and providing us with your valuable feedback.\n\n> Weaknesses:\n> 1. Some quantities could be elaborated in the main paper. e.g., $D_k$ in line 153. I am not sure if the form of the diagonal matrix is important or a minor issue that will only be useful for theoretical analysis.\n\n**Reply:** Thank you for your insightful comment.  The primary purpose of constructing this diagonal matrix is to ensure the theoretical convergence of our algorithm. In practical implementation, however, there is no need to explicitly compute this diagonal matrix. This crucial distinction between the theoretical role and practical implementation of the diagonal matrix is indeed worth emphasizing. \n\nWe appreciate your suggestion and will clarify this point in our revised manuscript.\n\n> 2. It might also be good to mention the special case of the algorithm without the disconnectivity constraint and how it compares with other existing methods.\n\n**Reply:** We appreciate your suggestion to emphasize the performance of our algorithm when the disconnectivity constraint is not applied. Indeed, such experimental comparisons have been conducted and are displayed in Figures 2 and 3, demonstrating the superior computational efficiency of our algorithm relative to existing methods.\n\nIn light of your valuable feedback, we will refine our manuscript to better highlight this aspect, emphasizing the robustness and efficiency of our algorithm under various conditions.\n\n> 3. It seems to me that the choice of regularization parameter $\\lambda_{ij}$ doesn't affect the results in Theorem 4.4 \\& 4.5. Could you elaborate more on this?\n\n**Reply:** Thank you for your insightful observation regarding the role of the regularization parameter in Theorems 4.4 and 4.5.\n\nYou're correct in noting that the choice of the regularization parameter does not impact the conclusion of these theorems. This is primarily because it acts as a linear term in our optimization problem. Consequently, it does not influence the convergence behavior of our algorithm.\n\nWe appreciate your constructive comment. In our revised manuscript, we will provide a more detailed explanation on this matter.\n\n\n"}]

Author rebuttal:
Dear Reviewers,

We sincerely thank you for dedicating your time to review our paper and for your insightful comments. These have significantly contributed to improving the clarity and overall quality of our manuscript.

In this global response, we embrace the opportunity to clarify the distinction between our paper and paper [A], and our unique contributions relative to [A]. This issue was raised by one of the reviewers. It is worth mentioning that Paper [A] was published after an earlier version of our work had already been posted on arXiv and it cites this earlier version.

While both papers aim to advance the field of MTP2 graphical models, they do so with different emphases and distinct methodologies, thus offering different contributions.

### 1. Divergent Emphases and Methodologies Between Our Paper and Paper [A]: 

Paper [A] primarily targets the design of an estimator that can concurrently achieve superior performance in precision matrix estimation and graph edge recovery. **The authors' emphasis in paper [A] leans more towards enhancing estimation performance, even at the expense of computational efficiency.** They present a multiple-stage estimation method, founded on the framework of multiple convex relaxation, which progressively refines estimates across stages. This requires the repeated solving of the constrained log-determinant program across stages. Their theoretical analysis offers an upper bound of the estimation error between their proposed estimator and the underlying precision matrix, but does not establish a guarantee for the algorithm's convergence or the convergence rate. Their experimental results primarily demonstrate the performance of estimating precision matrices and recovering graph edges using metrics like relative estimation error, true positive rate, false positive rate, and F-score. However, this method is computationally challenging in high-dimensional scenarios due to the relatively slow convergence of the PGD and its repeated need to use PGD to solve the log-determinant program across stages.

In contrast, **our paper emphasizes the enhancement of computational efficiency, recognizing it as a critical requirement when dealing with large-scale problems.** Our focus is on developing an efficient and scalable algorithm tailored for estimating large-scale MTP2 graphical models. Both per-iteration complexity and convergence rate are of paramount importance in this endeavor. We propose a novel second-order algorithm, grounded in the two-metric projection method. This algorithm substantially reduces per-iteration complexity by leaving a special structure of the Hessian. It also achieves a rapid convergence rate by incorporating second-order information into the search direction. Our empirical study primarily compares the runtime of various algorithms in terms of their convergence to the global minimizer of a common convex log-determinant program, as depicted in Figures 2-5. It's important to note that all algorithms share the same estimation performance in these experiments, given that they all solve the same problem.

### 2. Unique Contributions of Our Paper:

Existing algorithms designed to solve the sign-constrained log-determinant program, such as Projected Gradient Descent (PGD) and Block Coordinate Descent (BCD), encounter computational challenges in high-dimensional scenarios due to their slow convergence rate or high per-iteration computational complexity. To mitigate these issues, **we propose a novel projected Newton-like algorithm, which achieves a significantly faster convergence rate than first-order algorithms while maintaining the same per-iteration complexity.** We provide solid evidence of its guaranteed convergence to the global minimizer and ascertain its theoretical convergence rate. Our extensive experiments demonstrate that our proposed algorithm significantly enhances computational efficiency compared to state-of-the-art methods.

It is worth noting that paper [A] employs the PGD algorithm to tackle the constrained log-determinant program at each stage. However, being a first-order method, it is often hindered by computational inefficiency due to its slow convergence rate. In contrast, **our proposed algorithm can boost the computational efficiency of the method outlined in paper [A]** by applying it to solve the log-determinant program across stages, thereby enabling its application to higher-dimensional problems. However, it's important to acknowledge that PGD, as demonstrated in paper [A], has a distinct advantage in handling complex constraints, such as the diagonally dominant M-matrices, a capability currently not supported by our algorithm.

In the revised manuscript, we will ensure to appropriately cite paper [A] and clearly articulate its relevance to our work. Once again, we express our gratitude to all reviewers for your valuable feedback, which greatly contributes to improving the clarity and thoroughness of our manuscript.

Reference:

[A] J. Ying, J. V. D. M. Cardoso, and D. P. Palomar. "Adaptive Estimation of Graphical Models under Total Positivity." ICML, 2023.




Meta-review: In this paper, the authors considered the problem of sparse precision matrix estimation under total positivity and developed a fast projected Newton-type algorithm for efficiently solving the formulated sign-constrained log-determinant program. After discussions, the reviewers reached a majority consensus that the proposed algorithm is theoretically sound and numerically convincing. One of the reviewers expressed a major concern regarding the missing comparison to a closely relevant work, which has been thoroughly addressed by the authors in their responses. Based on the reviews, authors responses and post-rebuttal discussions, it was assessed that this submission makes a sufficiently novel and solid contribution to the addressed problem and thus is worthy of publication.
Decision: Accept (poster)