Partial Correlation Network Estimation by Semismooth Newton Methods
Authors: DongWon Kim
Sungdong Lee
Joong-Ho Won
Venue: NeurIPS 2025 poster ()
Keywords: graphical model selection
semismooth Newton methods
parallel algorithm
global convergence
$\ell_1$-regularization
convex optimization
Primary area: optimization
TL;DR: We develop a scalable second-order algorithm for pseudolikelihood-based partial correlation network estimation.

Abstract: We develop a scalable second-order algorithm for a recently proposed  $\ell_1$-regularized pseudolikelihood-based partial correlation network estimation framework. While the latter method admits statistical guarantees and is inherently scalable compared to likelihood-based methods such as graphical lasso, the currently available implementations rely only on first-order information and require thousands of iterations to obtain reliable estimates even on high-performance supercomputers. In this paper, we further investigate the inherent scalability of the framework and propose locally and globally convergent semismooth Newton methods. Despite the nonsmoothness of the problem, these second-order  algorithms converge at a locally quadratic rate, and require only a few tens of iterations in practice.  Each iteration reduces to solving linear systems of small dimensions or linear complementary problems of smaller dimensions, making the computation also suitable for less powerful computing environments. Experiments on both simulated and real-world genomic datasets demonstrate the superior convergence behavior and computational efficiency of the proposed algorithm, which position our method as a promising tool for massive-scale network analysis sought for in, e.g., modern multi-omics research.

Partial Correlation Network Estimation by Semismooth Newton Methods
Dongwon Kim Department of Statistics Scoul National University dongvonida@snu ac kr
Sungdong Lee" Department of Medicine National University of Singapore sdlee23@nus wwtp:/mmuuuuuuu......ee. edu..g
Joong-Ho Won Department of Statistics Scoul National University won. j@snu. ac kr
Abstract
We develop scalable second-order algorithm for a recently proposed 4 regularized pseudolikelihodd-be partial correlation network estimation frame- work. While the latter method admits statistical guarantees and inherently scalable compared to likelihood-basee methods such as graphical lasso, the cur- rently available implementations rely only on first-order information and require thousands of iterations to obtain reliable estimates percomputers. In this paper, we further investigate the inherent scalability of the framework and propose locally and globally convergent semismooth Newton meth- ods. Despite the nonsmoothness of the problem, these second-order algorithms converge at locally quadratic rate, and require only few tens of iterations in practice. Each iteration reduces to solving linear systems small dimensions or linear complementary problems of smaller dimensions, making the computation also suitable for less powerful computing environments. Experiments on both simulated and real-world genomic datasets demonstrate the superior convergence behavior and computational efficiency of the proposed algorithm, which position our method as a promising tool for massive-scale network analysis sought for e.g., modern multi-omics research.
-eerformmme
Introduction
partial correlation network refers to graph in which each node (vertex) corresponds to a variable and an edge between two nodes represents the partial corrrelation between the corresponding pair of variables. Partial correlation is defined as the correlation coefficient between the residuals resulting from the linear regression of each of the two variables with the rest of the variables, thereby measuring the degree association between the two adjusting for the effect of the remaining variables. Learning partial correlation networks form data is particularly important in multi-omics studies of biological systems, where spurious edges due to shared regulatory factors should be avoided to detect direct
biological interaction of causal implications 227 24 22
Partial correlation network estimation closely connected covariance selection (6, which estimates nonzero entries of the inverse covariance precision matrix of the variables. This is because nonzero entry of the precision matrix implies a nonzero partial correlation between the corresponding variables. Moreover, if normality of the data distribution is assumed, then the zero entries of the precision matrix imply conditional independence. For this reason, when the number of variables large and it is reasonable to assume that the number of edges is small, the graphical lasso [ 34 that maximizes the Gaussian likelihood under & penalty on the entries of the precision matrix or solve the following optimization problem
min log det O (r((@) + x||O-pll1...
(1)
'Equal contribution
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

is by far the most popular choice. Here, S" denotes the space of p-dimensional symmetric matrices, and S = (1/n)XX* X is the sample covariance matrix of the data X = [rI....,Za] E RExP, The ① D denotes matrix that takes the off-diagonal elements of θ with zero diagonal elements. The & norm || sums up the absolute values of the elements of its input matrix. The solution ê to estimates the target precision matrix O¹, such that zi
+id (00 (0")-1).
However, many methods that directly solve ① suffer from computational bottleneck when p reaches few thousands, which is not desirable for analyzing modern massive-scaled data such as arising from multi-omics. The root cause is that the Keachh-Kuh--Tukkke (KKT) optimality condition of (1)
o 1 + XZ
= 0, Z E 0
where 01(@- denotes the subdifferential of the & norm at @ involves matrix inversion.
Direct computation of the inverse θ requires O(r²) arithmetic operations, Since most existing
can be understood as solving (2) in an iterative fashion, such a costly computation essentially occurs every iteration. Worse, matrix inversion is difficult to parallelize or distribute, limiting scalability.
algorithms for solving []] (e.g.. 5 8 20 13)
A recently proposed framework called ACCORD [9] overcomes this limitation by employing a computation-friendly loss function in a one-to-one transformation of the precision matrix variable O. It solves a convex optimization problem
min - log det LD + (1/2) tr(S2" Al|S>
(3)
where Ltp denotes the diagonal matrix that takes the diagonal elements of 52. Precision matrix o and the optimization variable Ω are related by Ω = o 0.1/2 and o = Rp Note Ω may not be symmetric. However, the sparsity patterns of 22 and o coincide with each other, so the & penalty on f2_D is legitimate. Under some regulatory conditions, the estimate 52 that uniquely 1/ minimizes (3) is a consistent estimator of = 1 and its sparsity pattern. The ACCORD framework lies in line with fi-penalized pseudolikelihoodd-asee methods for partial correlation network estimation. These methods are computionally more scalable than the graphical lasso owing to the use of pseudolikelihood, which in turn more robust to non-normality. SPACE [2] is the first of its kind, which directly targets partial correlations. While SPLICE (32 and SYMLASSO p follow this line, these methods are not convex, and convergence of the corresponding optimization algorithms are not established. CONCORD 17 convexifies SPACE by changing the target variable to precision matrix O. Nevertheless, consistency of the estimated precision matrix requires quite strong conditions. ACCORD relaxes those by asymmetrizing the target variable and establishes non-asymptotic bounds on estimation error and sample complexity [19]
The key computational attraction of ACCORD is its scalability. In 3 the log-determinant barrier is applied only to the diagonal elements of Ω2 hence reduces to merely a sum of logarithms. This makes the KKT conditions for (3) much more tractable than (2). In [19] proximal gradient algorithm proposed to solve 3 at a linear rate. The major computational bottleneck of this algorithm is sparse-dense matrix-matrix multiplication, which easily allows for distributed computation [18). This contrasts ACCORD with the graphical lasso. For massive-scale data from, e-g.. multi-omic studies, [19] implements the algorithm on distributed-memook systems equipped with the message passing inferface (MPI). As a showcase, the liver hepatocellular carcinoma (LIHC) cohort from The Cancer Genome Atlas [1] with p s 300, 000 molecular features, for which the most scalable graphical lasso algorithm [14 fails, are analyzed. Despite the success, this computation takes 200 node-hours on a powerful HPC system in a supercomputing center, accessibility to which is generally limited to public. For wider uses of ACCORD among a large number of scientists, computational algorithm that can operate in smaller but more accessible environment, such as a GPU workstation, and possesses a convergence rate faster than the linear rate of the proximal gradient algorithm (to compensate for the limited computing power) is on demand.
high-performance computing (HPC)
This paper addresses such demand by proposing a semismooth Newton method for solving (3). Because of its superlinear convergence property, semismooth Newton methods [3] have been widely We develop
employed for learning problems with non-smooth objective functions [21 23 33 35
the semismooth Newton method for ACCORD through two steps. First, we show that problem
(3) is rowwise separable, which allows it to decompose into separate p-dimensional minimization problems for each row of f2 Second, we show that the Karush- Kuhn-Tucke (KKT) condition for

---

each minimization problem can be formulated as finding a root of a nonlinear multivariate function F. which nondifferentiable but strongly semismooth. This observation ensures that the iteration converges at locally quadratic rate. The key to implementation is to find suitable element of Clarke's generalized Jacobian (which is a set) of F. We achieve this goal by deriving the Bouligand derivative of F. The latter allows us to globalize the semismooth Newton algorithm, for which the initial point can be chosen with liberty. A connection to the primal-dual active set (PDAS) strategy [12 is also made. We then examine the performance of our algorithm in both synthetic and the LIHC data, with concluding remarks provided thereafter. Proofs the results are provided in the Appendix.
2
Locally convergent semismooth Newton method
2.1 KKT condition as root-finding
In this subsection, we show that the KKT condition for (3) can be formulated finding root of a multivariate nonlinear function. We begin with observing that the objective function in (3) is
separable in rows. Let St be the i-th row of f. Then, (3) can be rewritten finding minimizer of
Σ (-log det 22, (1/2) tr 225) + X11
P
where i and - denote the diagonal and the off- di-gooal parts of I4 respectively. Equation (4) clearly shows that the objective function is separable in Z2 s. Henceforth, will derive the optimality condition for i = without loss generality. Let w = (wi.w",
where wi = -2 and w-1 -1 Then, the objective function minimize (w) = g(w) h(w)), where
g(w) = (1/2)w" Sw,
h(w) log + A|w-11l1,
which are both convex functions. By introducing dual variable d (d1, d_,)T, the optimality condition 0f(w*) can be expressed as
Sw + 0, d' @h(w")
(5)
where 0ʃ and Ōh denote the subdifferentials of f and h, respectively. Moreover, the condition d' Öh(w") can be equivalently expressed using the proximal operator of h [15
w proxx d²),
(6)
which can be explicitly written as
prox, (w) arg min (M(v) (1/2)]lu w [1]] = - ww++ Tx(w1)7))
where (x = (sign(w,) max( ( - X, 0), #1 is an element-wise thesshhresn operator.
Combining and (6), the KKT condition for can be equivalently written
F(w", d3)) 0,
where F: 22 2p is defined
Sw + d
x
F(w,d)
W1 prox log(.) (wi + d1)) W-1 (2(w-1 + d-1)
prox log() (og))()))))))))))))))))))
(8)
2
Therefore, solving 3 amounts to finding a root of nonlinear equation 7 Map F has a unique root
almost surely, long the sample is drawn from continuous distribution [19 Theorem 3.1].
2.2 Newton's method for semismooth functions
The formulation 7 8 of the KKT condition naturally suggests Newton's method for finding the
root. However, because of the presence of the soft thresholding operator in (8). it not differentiable everywhere, which means that the plain Newton iteration cannot be applied. In this section, we show that, despite the nondifferentiablity, function in K is strongly semismooth, hence semismooth version of Newton's 's method can be employed. We begin by defining semismoothness [26] 31).
3

---

Definition 2.1 (Semismoothness), Function F:; R" + R' is said to be semismooth at x E R" if it is locally Lipschitz, all directional derivatives F (x; v) == lim ((a F(x) exist, and for any G E OF + v), the following holds:
G (x;)) ([[v]))
Furthermore, F is strongly semismooth at r if it is semismooth at a and satisfies the condition for any G OF +
Ge The F is said be (strongly) semismooth on U E m if it (strongly) semismooth at any x E U.
Ge - Fa; v) = ((1l=
Theorem 2.2. The map (w, d) + F(w, d) defined in 8 is strongly semismooth.
The semismooth Newton method 31 relies on the notion of the generalized Jacobian 3
Definition 2.3 (Generalized Jacobian). For function F : I >> Rd that is locally Lipschitz continuous around z E Rm", Clarke's 's generalized Jacobian at is defined as
lim VF(x") I² +x, ) k
0F(z) conv
where (x(x) = {y: F is differentiable at y). conv denotes the convex hull operation, and V F(v) represents the Jacobian F at y. We call element of OF(x) generalized Jacobian of F at r.
Thus for locally Lipschitz functions, version of Newton's 's method for solving F(x) = 0 can be defined as
x A+1 = x² G(a G(a*) E OFF(*")
provided that G(z*) is nonsingular. If F semismoothness, then the Newton iteration 9 exhibits superlinear convergence for suitable initialization BI] It is strongly semismooth, then it converges at locally quadratic rate:
Theorem 2.4 (31). Suppose that x² is root of F Rm RR. and that F is strongly semismooth an open neighborhood U containing x², with a generalized Jacobian G(x) E aF(x). If nonsingular for all x E U and the set {[(G(x)) : E U} is bounded, then the semismooth Newton method ( converges quadratically to x", provided that ||x00 - xx" || is sufficiently small.
To verify the assumptions of Theorem for the target problem 7 we provide sufficient condition.
Proposition 2.5. Suppose (w", d²) satisfies equation (7) for F defined in (8. Let A =
(2. W + dī > A). If A., the submatrix of X taking the columns indexed by A has full column rank, then there exists nonsingular G(w, d) E OF(w, d) for any (w,d) E R in neighborhood of (w- (d))
In particular, we may choose (see Appendix D for the derivation)
G(w, d) (1 S
1-J -1).
(10)
where J diag(J:) is diagonal matrix with J11 1 +dl +d; and JH 1i w d;1 > X.
Ja 0 if (w;; + d.l < À for i = .,...,P. the primal-dual active algorithm M 15).
p. In this case, semismooth Newton iteration is equivalent
3 Globalizing semismooth Newton method
3.1 B-semismooth Newton method
In this section, we propose globalization strategy for the semismooth Newton method (9). Although the quadratic convergence result (Theorem 2.4) is attractive, this result is only local and is silent on how to choose the initial point の = (w',, ; dd). An algorithm that converges from any initial point while maintaining locally quadratic convergence near the solution is more practical and desirable. To this end, it turns out that globalizing iteration (9) per se difficult, while globalizing the semismooth Newton based on Bouligand derivatives or B-derivatives (8 is manageable.

---

Definition 3.1 (B-diffrretiabillit). A function F I' + R' is said Bouligand differentiable (B-differentiable) at E I is locally Lipschitz, the directional derivative F' v) limego F(x+tu)-F(az exists for all and fulfill the approximation property [[ + v) - F(z)
(x; I| = ((Ix]]) as v + O. If is B-differentiable for any 2 E R", then it is said B- differentiable. The directional derivative that satisfies these conditions is called the Bouligand derivative (B-derivative) of at x in the direction v.
Using the directional derivative F"(x; v), Newton method for finding the root of B-differentiable function F can defined as
k+ x + + F'(x::e*)- -F(=").
(11)
Our target function is B-differentiable:
Proposition 3.2. Function F as defined in in the direction of v (Uwi Ua) given by
8 is B-differentiable. Moreover, its B-derivative at (w, d)
Sen + Udd
(z 1/2 + 1/2) Ud,i
i A
(12)
iE I
min (Weies Udd i I
max Uww.s - -dd
i I
where: 2 = 2(w1, d2) = twl 2244
td and A = A(w, d), I = Iw, d), I+ = II, d). and
== (w:(w,)) d) such that
A = E (2..... PP) |ui + dil {}.
T = E (2...., |ز + dd < {};
(13)
I e (2. wi + d }}.
I {i E (2..... p} w, d, -A}.
3.2 Solvability of the B-semismooth Newton method via linear complementarity problem
Applying
12 to
finding the Newton direction such that F' F'(w, d; v) = - -F d) reduces to
Su + vd SSw d (a + 1/2 + (2 - 1/2)ua, = -wi + prox log (w] + d1)
(14)
(15)
min(tw.i, ) wi; i I+, max cai) = Wi. E IT,
(16)
where z is the same as in 12) Determining the existence of solution to 14 16 and finding one
a quadratic programming problem in u = (ww va) E R However, under the condition similar to that of Proposition 3.2, only need to solve much smaller linear complementarity problem (LCP)
[4]. Let A = A(w, d) == (1) A and Á A(w, d) == AUI+ UI- from (3) and X A X A be
the submatrix of X obtained in a similar way to X A- in Proposition 2.5
Lemma 3.3 (cf. 10), Assume X ā has full column rank. Given (w, d) E 22, v = (wa, va) E 2
solves 14 16 if and only if
y x ((-Vd,I) + wII) ,V WI y y= + #WI
(- .
WI )T)
solves the LCP
x x y 0, (x,y) = O, and y Mx +
(17)
where M MMSS w. d) and qSS w. d) are given in Appendix E.2 The M is symmetric positive definite and (17 admits a unique solution.
A LCP can be solved efficiently using standard techniques such as Lemke's pivoting methods 4. Let SR.C be the submatrix of S with rows indexed by and columns indexed by C. With the solution to (T7) available, the Newton direction b in (TT) for the F in (8) can be found in closed form as stated in the following main theorem.
5

---

Theorem 3.4. In 13 let A = Alw", ,d"), ) A(we, and I= = (w" dd)). Assume X A and
has full column rank. Let x, be the unique solution to for an iterate (w ). ) In addition to (3). let A = fi E (2....,P) w + d > X) and ==
= {i E f2..... pt wt + d < }
Then, the B-semismooth Newton iteration (1) is well defl-ee with
-TT = A-=Wwj wzo, UAA 1(1 -1-)1 d²A,
(
Xx+
) +
XI-
w
((vE,) ( (V5,)) (16,))) )T, , -u
,". (16220) ==== =NN== (p0-(pp-00000
,0 ,0,, (1/,zz)
- )'0))
vd,1 = b² bu - c², N (-Or
(-CNT' 1
NT
where bb, ch E R. Pk R,, C CER'R R(P-IAIXIA and N1 RllÍIIÁÍ,, given in Appendix E.2 are determined by wh dh and S. In particular, N1 is symmetric positive definite.
In particular, given that I UI- = for I+ = E (2..... w/ di = A}. I-' = {i (2....,P) p) w di =},, the convergence of (1I becomes locally quadratie:
Corollary 3.5 (30) Corollary 3.3]). Let (w" ,d)) be the unique solution to Flw, d) = O. If I+
I- = 0 and X- from Proposition 2.5 has full column rank, then the B-semismooth Newton iteration I is well-defined and converges quadratically (w" , .d") in its neighborhood.
Remark 3.6. Multiplication with matrix N only requires a linear system solve involving the symmetric positive definite submatrix N1 of size [Á| |AA. In the sparse estimation context, the size of the active set Á is small, hence solving (T] is not a bottleneck for scalability.
Remark 3.7. Given that F' in (8) is locally Lipschitz and directionally differentiable, there exists generalized Jacobian G such that F" d: v) G(w,, d)v where satisfies (IT) B1] More specifically, the B-semismooth Newton method can be regarded as the semismooth Newton method that employs the generalized Jacobian Gn (10) such that Jn for i AUfi UI :x;-
and Ju = for i E I U {i E I+ UI- x; > 0}. where x is the solution to the LCP Furthermore, since it is numerically unlikely that (ws + de precisely equals to X, in most scenarios I+ -I- and the B-semismooth Newton method coincides with the PDAS algorithm.
[7
3.3 Globalizing B-semismooth Newton with search
Globalization of the B-semismooth Newton method (TT) is based on the observation that, under the (unique) existence of the solution to 7 solving this nonlinear equation is equivalent to minimizing 8(w, d) = 1 F(w, d)112- Given the Newton direction x² obtained by solving (14 (T)) for (w, d) (w*, ) we can apply line search to construct a damped Newton algorithm that monotonically decreases the objective sequence {d(w>, dd Following 25 that finds an integer m, such that
we consider the line search algorithm
0aa + p - 0(x*) S -2oppm0000 0(a)),
(18)
for algorithm parameters a, p E (0,
Also, Theorem 3.4 suggests that the full column rank of the submatrix X A is essential in computing the direction xx for A A(w=, dd)). If this condition is not met, then we can resort to the proximal gradient algorithm of [19] which converges linearly to the solution. Note that this algorithm does not necessarily descent (w,, d) in one iteration. Thus, we need repeat the step for ((w,, d) to decrease.
The resulting damped B B-semismooth Newton algorithm is summarized in Algorithm I The main question is whether the iterate sequence {(w*, dd²) will converge the solution (w dd and whether the Newton steps will eventually take place; these are all addressed in the following theorem.
Theorem 3.8. Suppose I+ UI- 0 and X á- from Proposition 2.5 has full column rank. Let {Ww*, ) } be the sequence of iterates generated by Algorithm| / Then, the proximal gradient step occurs only a finite number of times, and if we let {p denote the corresponding step sizes,
1. if lim sup nu 0, then (w* dd) + (W", d'), ) the unique solution to F(w, d) 0;

---

Algorithm 1: Damped B-semismooth Newton method with line search
Input:
Initial guess = (w": .d") for all rows, po: a E (0, 1), and tolerance tol
foreach row in St simultaneously do
l 2-4-0 0; while ((=")== > tol do
Compute sets Av A(w", d. [w d²)), and I I' (w",, (w", (0"
Set A A
if X A.. has fill column rank then
Find the search direction v according to Theorem 3.4 m + D; repeat
m + 1;
until 0(x Pk v0 S (1 20 (0(=)) 1+1 x² + Pk v
13
$ + x; repeat
16
Update 3 with proximal gradient algorithm
until (()) < A(*");
19
k 1;1;
2. if lim s0 pm O and (wi, (d) is accumulation point of {w ) then (wt, dtt) = (w" d and (wh d²) - (w",d"))
Moreover. if (w* ) + (w", d, then the convergence is locally quadratic.
Comparison to the SSNAL framework There is line of work that falls into the semismooth Newton-based augmented Lagrangian method (SSNAL) framework, developed for solving the lasso- related regression problems [21] 23 35 As the name suggests, this framework solves the problem via the inexact augmented Lagrangian method (ALM). ALM incurs an optimization subproblem. and SSNAL solves this subproblem via semismooth Newton. While the subproblem is solved at locally superlinear (quadratic it is strongly semismooth) rate, the convergence rate of the the outer ALM iteration is locally linear, even if the subproblem be solved at a locally quadratic rate. Global convergence is also established.
In contrast, Algorithm formulates the KKT condition as strongly semismooth equation and apply semismooth Newton. Globalization is achieved via damping instead of ALM. So our algorithmic design is fundamentally different from SSNAL: there is no subproblem to solve inexactly, and the convergence is locally quadratic. (If the KKT equation were merely semismooth, then the local convergence would be superlinear.)
Another distinction is that SSNAL typically constructs augmented Lagrangian via the dual of the original (primal) problem, while we focus on the primal. It is worth pointing out that (33 applies the SSNAL framework directly to the primal addition the dual, but design and analysis focus on low-dimensional settings where the number of observations significantly exceeds the number of features (i.e., n > p). As such, it not readily applicable to high-dimensional settings, which are the primary focus of our work.
Experiments
4.1 Convergence behavior Damped B-semismooth Newton method
We present numerical experiments to empirically validate the convergence and descent property of the proposed damped B-semismooth Newton method. Specifically, we compare Algorithm againstt the proximal gradient algorithm proposed in [19), by tracking the following three metrics over the

---

iteration:
IIF(w" d") )²2 (objective function of Algorithm
f(w*) f()), (objective function ACCORD)
(19)
1w"-w"ll 11 (squared distance to the optimum)
The @ estimates the optimal solution w" by the output of Algorithm
D convergence. Convergence The rank test for line 6 in
declared if all rows satisfies the criterion [| dI2
10
Algorithm is based on the sufficiency check [A& S n. while more stringent tests (e.g.. condition number singular value thresholds) are also possible. For the line search parameter (line 11), Pk: = max 0.003 k. 0.4) and a = 0.001. The algorithm is implemented in PyTorch 1.13.1. Matrix operations are parallelized six NVIDIA RTX 6000 Ada generation GPUs using CUDA 11.7.
To assess the convergence behavior, generated dimensional zero-mean multivariate Gaussian distribution with precision matrix O that contained 3% non-zero entries, locations of which were sampled uniformly at random. The regularization parameter was set to 0.1, 0.15, and 0.2 to ensure that the estimated precision matrix exhibits sparsity level comparable to O*.
500 observations from a p 1, 000-
[[F(ut ghiP
1.8
16-
20-4
Iteration
iteration
Itoration
method
B-semi
Prox
0.1
0.15
0.2
Figure 1: Convergence comparison: damped B B-semismooth Newton proximal gradient method
The results are summarized in Figure where the vertical axis is plotted on logarithmic scale. The proximal gradient method exhibits an approximately linear convergence rate with respect to both the ACCORD objective function f and the squared distance to the optimum, consistent with the theory 19. In contrast, the damped B- B-semismooth Newton method demonstrates superlinear convergence for both metrics, achieving approximately order of magnitude faster convergence. Furthermore, for [l> dd*) ][2 22 which Algorithm has descent property, the convergence appears locally quadratic (as predicted by theory) and occurs within a few tens of iterations.
As Algorithm inclldess two components addition the semismooth Newton steps, namely the proximal gradient steps and the LCP solve step, we also conducted an ablation study to assess the contribution of each component. See Appendix G.I for details.
4.2 Graph density-convergence trade-off
In Algorithm [ the B-semismooth Newton update cannot be applied to the row in which rank(X A,) min(n, p). Since the size of the active set is inherently linked to the sparsity of the true precision matrix @* (hence that of I2"), we vary the regularization parameter À to study the robustness of the algorithm (in terms of the occurrence of the Newton step) to the density of the underlying partial correlation graph.
To conduct this investigation, we define the graph density as Density [F) where [F] denotes the number of undirected edges and p is the number of nodes in the graph implied by O'. Based on this definition, we constructed graph structure sampled uniformly at random with densities of 0.03, 0.05, and 0. Note that the graph density in Section 4.1 is roughly 0.03. The partial correlation matrix corresponding to the graph was constructed to have the minimum eigenvalues of least 0.2, with all non-zero entries having absolute values no less than 0.1. Based on this matrix, we generate 100 observations from a 1, 000-dimensional zero-mean multivariate Gaussian distribution with precision matrix **

---

The hyperparameter À was chosen so that the estimated solution exhibits a sparsity level comparable to that of the true precision matrix O while the line search parameters were set to p 0.5, 0.001. We then assesses whether convergence is achieved within 300 iterations for each case. This limit of iterations was chosen because from Section 4.1 we see B-semismooth Newton method typically converges within few tens of iterations. The results are presented in Table !
Density \ A
0.2
0.15
0.1
0.09
0.08
0.07
Density \ 0.07
0.06
0.05
0.03
Y
Y
0,05
Y
Y
0.10
9.1%
97.7%
100%
0.10
Y
N
Table 1: Convergence within 300 iterations for varying Table 2: Proportion of rows updated by graph density and À (Yes / No)
proximal method after 300 iterations
With one exception, I converged within 300 iterations. The exception occurred when the underlying graph was the densest and the regularization parameter was the smallest. Usually the denser the underlying graph becomes, the smaller regularization parameter required to recover its structure. To further investigate this phenomenon, we summarize in Table the proportion of rows updated by the proximal gradient method after 300 iterations for smaller values of X.
This result shows that there is trade-off between the density of the underlying graph to estimate and the convergence speed of Algorithm [ Nevertheless, for consistent estimation of the support of the underlying graph, sparsity assumption is essential, which is tied to the identifiability condition of Proposition 3.2 77. Under this condition, the Newton step occurs eventually. Adaptive selection of the regularization parameter À under the sparsity regime, e.g., continuation strategy 7 [5]. is an interesting future direction of study.
Further experiments on other types of random graphs were also conducted, the results of which can be found in Appendix G.3
4.3 Application to large-scale genomic data
In this subsection, we assess the performance of Algorithm | in estimating the partial correlation graph from large-scale multi-omic dataset. Specifically, we apply Algorithm I and proximal gradient method (ACCORD-FBS) to the LIHC dataset from The Cancer Genome Atlas 1. consisting of 305, 471 features including RNA transcription levels and the DNA methylation status of human genomes from n = 365 samples. Estimating partial correlations such high-dimensional biological data facilitate the identification and validation of epigenomic and transcriptomic regulatory mechanisms human cells. This dataset was analyzed in [19] in which an HPC version of ACCORD-FBS was we employed our own implementation of ACCORD-FBS for the GPU workstation of the previous subsection, utilizing the row-wise separability.
national supercomputing centers for several hours. We emphasize that
331
104
231
101
337
234
10.1
104
40
3.0
340
120
140
150
Iteration
Iteration
method
Algorithm 1 .... ACCORD-FDS
-
0.2
0.45
-de- 0.6
Figure 2: Convergence comparison on a large-scale multi-omic cancer dataset
The regularization parameter was set to À = 0.45, selected via the extended mation criterion 19: detailed sensitivity analysis of estimation performance around this value of X provided Appendix G.2 test robustness of Algorithm T we further examined the optimization performance for À E {0.20, 0.60), corresponding to relatively dense and sparse estimates. Figure
infor-

---

illustrates the convergence behavior of Algorithm I and ACCORD-FBS, evaluated using the first two metries defined (9). While ACCORD-FBS typically requires over thousand iterations to converge [99, our algorithm reached the solution within only a few dozen iterations. Despite the reduced number of iterations, the computational cost per iteration remains nearly unchanged, as the dominant operation in each step involves only the inversion of a relatively small submatrix of dimension | x [Ag].. Due to the rapid convergence of the B-semismooth Newton method, each iteration involves significantly fewer row updates compared to ACCORD-FBS.
Method \ A
0.2
0.45
0.6
Time (s)
|F|l
Time (s)
|||
Time (s)
||]
Algorithm ACCORD-FBS
5341
1.260 x 10*
3352
0
2083
6778
4.709 x 10
5150
4.035
102
3452
2.671 x 10?
Table Wall-clock time (in seconds) and ||]'| until convergence 100 iterations
As shown Table Algorithm| consistently achieved shorter wall-clock time and smaller [|] than ACCORD-FBS for all tested values of À within 100 iterations at convergence, whichever came first.
Algorithi
Rows Updated Methoz
300002
B-semi
ACCORD-FBS
250002
Praximal
Nonzero Entries
200002
250002
100002
30²
20002
Iberation
Iteratioo
Figure 3: Interaction between sparsity and update strategy in Algorithm I
As expected, more iterations were required to converge as the regularization parameter À decreases, and more likely the proximal gradient step occurs (rather than the Newton step) in Algorithm [] As means of dissecting this phenomenon, we fix at 0.45 and monitor the number of nonzero elements in the iterate sequence 5 and the number of rows updated by either the two methods; see Figure 3 In the early stage, updates the number of nonzero entries in s2 decreased and the proportion of rows updated by the Newton step increased. This transition aligns with intuition and enables Algorithm to take advantage of the fast convergence behavior of the B-semismooth Newton method in later iterations.
performed primarily by the proximal gradient method. Over time,
Conclusion
We explored the applicability of semismooth Newton methods in partial correlation network esti- mation to derive a scalable algorithm with superlinear convergence. These methods exhibit locally quadratic convergence to the solution, where the global convergence is also achievable by combin- ing with the proposed damping strategy. They also enjoy parallelism by leveraging the rowwise separability of the ACCORD objective function. Simulations and genomic data applications GPUs validate the superior computational performance of the semismooth methods compared to the proximal method. In conclusion, semismooth Newton methods suggest promising outlook for fast and efficient partial correlation analysis large-scale data with modest computing resources.
6 Acknowledgements
Lee was partially supported by the Sejong Science Fellowship through the National Research Foundation of Korea (NRF) funded by the Ministry of Science and ICT (MSIT, No. RS-2024- 00347617), Singapore Ministry of Education (MOF-000244-00,, MOR-000617-00), and National Medical Research Council of Singapore (MOH-000986-00), D. Kim and J.-H. Won were supported
10

---

by the AI-Bio Research Grant through Seoul National University (No. 0413-20230050), Artificial Intelligence Graduate School Program (Seoul National University), and by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT, No. RS-2024-000336611)
References
[1] Adrian Ally, Miruna Balasundaram, Rebecca Carlsen, Eric Chuah, Amanda Clarke, Noreen Dhalla, Robert A Holt, Steven JM Jones, Darlene Lee, Yussanne Ma, et al. Comprehensive and integrative genomic characterization of hepatocellular carcinoma. Cell, 169(7)11327-3411, 2017.
[2] Dejun Chu, Changshui Zhang, Shiliang Sun, and Qing Tao. Semismooth Newton algorithm for efficient projections onto norm ball. In Proceedings of the 37th International Conference on Machine Learning, volume 119, pages 1974-1983. PMLR, 2020.
[3] Frank H Clarke. Optimization and nonsmooth analysis. SIAM, 1990.
[4] Richard W Cottle, Jong-Shi Pang, and Richard Stone. The linear complementarity problem. SIAM, 2009.
[5] Alexandre d' Aspremont, Onureena Banerjee, and Laurent El Ghaoui. First-order methods for sparse covariance selection. SIAM Journal on Matrix Analysis and Applications, 30(1):56-66, 2008.
[6] Arthur P Dempster. Covariance selection. Biometrics, 28(1):157-175, 1972.
[7] Qibin Fan, Yuling Jiao, and Xiliang Lu. primal dual active set algorithm with continuation for compressed sensing. IEEE Transactions on Signal Processing. (2(23):6276-6285 2014.
[8] Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Sparse inverse covariance estimation with the graphical Lasso. Biostatistics, 9(3):432 441, 2008.
[9] Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Applications of the Lasso and grouped Lasso to the estimation sparse graphical models. Technical report, Department of Statisties, Stanford University, 2010.
[10] Esther Hans and Thorsten Raasch. Global convergence of damped semismooth Newton methods for & tikhonov regularization. Inverse Problems, (1(2):25005, 2015.
[11] Michael Hintermüller. Semismooth Newton methods and applications. Technical report, Department of Mathematics, Sumbboldt-Universiyy of Berlin, 2010.
[12] Michael Hintermüller, Kazufumi Ito, and Karl Kunisch. The primal-dual active set strategy as a semismooth Newton method. SIAM Journal Optimization, 13(3):865-88) 888, 2002.
[13] Cho-Jui Hsieh, Mátyás A Sustik, Inderjit Dhillon, and Pradeep Ravikumar. QUIC: quadratic approximation for sparse inverse covariance estimation. Journal of Machine Learning Research, 15(1):2911 1-477 2014.
[14] Cho-Jui Hsieh, Mátýás A Sustik, Inderjit S Dhillon, Pradeep K Ravikumar, and Russell Poldrack. BIG & QUIC: Sparse inverse covariance estimation for million variables. Advances in Neural Information Processing Systems, 26:3 (6:3165-3173. 2013.
[15] Jian Huang, Yuling Jiao, Bangti Jin, Jin Liu, Xiliang Lu, and Can Yang. A unified primal dual active set algorithm for nonconvex sparse recovery. Statistical Science, 36(2):215-23) 2021.
[16] Kaifeng Jiang, Defeng Sun, and Kim-Chuan Toh. A partial proximal point algorithm for nuclear norm regularized matrix least squares problems. Mathematical Programming Computation, 6(3):281-325, 2014.
[17] Kshitij Khare, Sang- Yun Oh, and Bala Rajaratnam. A convex pseudolikelihood framework for high dimensional partial correlation estimation with convergence guarantees. Journal of the Royal Statistical Society Series B: Statistical Methodology, 77(4):803-82) 825, 2015.

---

[18] Penpom Koanantakool, Alnur Ali, Ariful Azad, Aydin Buluc, Dmitriy Morozov, Leonid Oliker, Katherine Yelick. and Sang-Yuu Yun Oh. Commumication-avoiii distributed massive-scale sparse inverse covariance estimation. In Proceedings of the Twenty- First International Conference on Artificial Intelligence and Statistics, volume 84, pages 1376- 1386. PMLR, 2018.
optimization methods for
[19] Sungdong Lee, Joshua Bang, Youngrae Kim, Hyungwon Choi, Sang- Yun Oh, and Joong- Ho Won. Learning massive-scale partial correlation networks in clinical multi-omics studies with HP-ACCORD. The Annals of Applied Statistics, 2025. to appear.
[20] Lu Li and Kim-Chuan Toh. An inexact interior point method for L -eeguularree sparse covari- ance selection. Mathematical Programming Computation, 2(3):291-314 2010.
[21] Xudong Li, Defeng Sun, and Kim-Chuan Toh. A highly efficient semismooth Newton augmented Lagrangian method for solving Lasso problems. SIAM Journal on Optimization, 28(1):433 458, 2018.
[22] Faming Liang and Bochao Jia. Sparse graphical modeling for high dimensional data: paradigm of conditional independence tests. Chapman and Hal//CRC, 2023.
[23] Meixia Lin, Yong-Jin Liu, Defeng Sun, and Kim-Chuan Toh. Efficient sparse semismooth Newton methods for the clustered Lasso problem. SIAM Journal Optimization, 29(3): 2052, 2019.
[24] Chenjin Ma, Mengyun Wu, and Shuangge Ma. Analysis of cancer omics data: a selective review of statistical techniques. Briefings in Bioinformatics, 23(2):bbab585, 2022.
[25] José Mario Martínez and Liqun Qi. Inexact Newton methods for solving nonsmooth equations. Journal of Computational and Applied Mathematics, 60(1-2): 1995.
[26] Robert Mifflin. Semismooth and semiconvex functions in constrained optimization. SIAM Journal on Control and Optimization, 15(6)995-9972, 1977.
[27] Felicia SL Ng. David Ruau, Lorenz Wernisch, and Berthold Göttgens. A graphical model approach visualizes regulatory relationships between genome-wide transcription factor binding profiles. Briefings in Bioinformatics, 19(1):162-173; 2018.
[28] Jong-Shi Pang. Newton's 's method for B-diffcrentiable equations. Mathematics of Operations Research, 15(2):311-341 1990.
[29] Jie Peng, Pei Wang. Nengfeng Zhou, and Ji Zhu. Partial correlation estimation by joint sparse regression models. Journal of the American Statistical Association, 104(486):735-746 2009.
[30] Liqun Qi. Convergence analysis of some algorithms for solving nonsmooth equations. Mathe- matics of Operations Research. 18(1):227-244 1993.
[31] Liqun Qi and Jie Sun. A nonsmooth version of Newton's 's method. Mathematical Programming, 58(1):353 367, 1993.
[32] Guilherme V Rocha, Peng Zhao, and Bin Yu. path following algorithm for sparse pseudo- likelihood inverse covariance estimation (SPLICE). arXiv preprint arXiv:0807.3734 2008
[33] Joong-Ho Won. Proximity operator of the matrix perspective function and its applications. Advances Neural Information Processing Systems, 33:6305 10300566314-- 631,, 2020.
[34] Ming Yuan and Yi Lin. Model selection and estimation in the gaussian graphical model. Biometrika, 94(1):19-35, 2007.
[35] Yangjing Zhang, Ning Zhang, Defeng Sun, and Kim-Chuan Toh. An efficient Hessian based algorithm for solving large-scale sparse group Lasso problems. Mathematical Programming, (79(1):22--26, 2020.
12

---

Supplementary
Material for 'Partial Correlation Network Estimation by Semismooth Newton Methods"
A Proof of Theorem 2.2
Recall that F is given in 8
Sw d
F(w,d) ()
W1 prox - log(-) (w) + d1) W-1 Tx(w-1 + d-1)
where S is the sample covariance matrix, Ty is the soft-thresholding operator, and
x 2x²+ +4
prox log(-)(z)
2
A.1 F is locally Lipchitz
To prove the local Lipschitz continuity of F, it suffices to show that each component of F's Lipschitz continuous. Note that each row of Sw + d is differentiable, with gradient (S,, 1)", which norm is bounded by v S +I and ||||] is largest eigenvalue of S. wi prox- log())(w) + di) is
also differentiable with gradient (1 2 (d)
Wi+d
wi 2 tdi where the absolute
value of each component are within [0, 1]. Finally, the soft-thresholding operator (w4 + d,) is Lipschitz continuous with a constant of 1, so we can conclude that each component of F is Lipschitz
continuous.
A.2 Directional differentiability of F
In here, only establish the directional differentiability of F, but also derive its directional derivative explicitly. A vector-valued function F R - R' is directionally differentiable at E ITT in the direction v if and only each of its components F; (for i I.....J I) is directionally differentiable at in the direction v. Note that it is not necessary to consider the differentiable components, as the directional derivative of differentiable function is given by F'(a;;)) w) = F=(x) .v. Accordingly, the directional derivative of the linear term Sw + d is simply +u Similarly, the directional derivative of the composite term w1 prox logt-) (w1+ +d1)
a point x in the direction e
is (= + 1/2) + (z 1/2)va,1 where
w 4:
2 (w) )244
Now, we compute directional derivative of the nonfdifferentiable component. We use the definition of active sets and inactive sets in (13). Then, the direction derivatives falls into one of the following cases:
i) E I+ (w, d) and twws + 0a, 0. Then i E A+ + tw d tua). In this case, F(w; + -Ww d; + tva.s) F(ws, di) = - tua,i + x wi- d, = -td The same argument works for i E I (w, d) and w, + Va,, < O.
ii E I (w, d) and uwst + Va, 0. Then, E I (w tow, d tua). In this case, F(w; + tuw di tua,)) - F(w,, di) = wi + to - Wi = tW The same argument works for I (wi d) and
Va,i Ud, 0.
i E A(w, d). Then, we can find to such that A(w, d) c A(w+ to d t+v for all < <ts to with
to S min inf (t > 0: |w5 + towss + d3 + wass > A}.
jEA(u,d)
For 1 5 to and E A(w, d) Alw Www, d + tva).
F(w; + tu de tva,s) F(ws, di) tUd,i
iv) i E (w, d). Then, we can find t 0 such that (w,,)) d) C. I + tu d + tua) for all
0 t 5 t with
t S min inf {t 0: jk + town + d, + toas < X}.
(
13

---

For 0 < <tsts t and E (w, d) C (w + tw d d tea),
F(w; + tow,ss d2 + tva,s) F F(wi, di) = t0w
Note that to and 1 are independent of Hence, directional derivative of F is (2).
A.3 Strong Semismoothness of
Now, prove that F is strongly semismooth. From BI Corollary 4.4], it is sufficient check the component-wise strong semismoothness of F. First, since Sw + d is an affine function, strong semismoothness holds trivially.
Next, verify the strong semismoothness of ((w) d1) v ( di)² 4)/2. It suffices to consider only the non-linear term H(w,, d1) == V(wi + d1) 4. The derivatives of (w1, d1) with respect to and dı are ((w++d)444 The directional derivative of H (wi, d1) in the direction of vector v = (wwr va) is given by
ti
0H Owr
01
WI + d1
(w dw Uws va)
Uw
Vww + d1)2 (v- + va).
dd
and the generalized Jacobian of Hw1 ow d1 va) is a matrix given by the partial derivatives with respect to wi and d1
V (H + hw,, d1 + ha) = (Viorttutdtttaa
the
th +d:
(with1+ +b22)))++
Then,
- Wi + how + di + ha
W1 + [Nw + ha).
Vh (w1, hi; bw: ha) = ((wi
how + di haa)
v Vwi di )2 +
Using a Taylor expansion, we obtain
w1 +ww tdl the
with
H(w + hw d1
/wi++h_ +di +hu)242
ha) - (w1, d1) = (l])))) Thus, Vh H ('1,,, di; hw: ha) = (114))))
Lastly, we have the facts that the soft thresholding operator is strongly semismooth everywhere [6 and w - w + d does not introduce any discontinuity or nonlinear behavior Tx(w), Therefore, w-1 (w + d-1) is strongly semismooth.
Since all components of are strongly semismooth, it follows that F is strongly semismooth.
B Proof of Theorem 2.4
Before proving Theorem 2.4 we present a well-known condition that are equivalent to semismooth- ness from [IT] Theorem 2.9, Theorem 2.1
Proposition B.1. Let F U m be defined on the open set U C R". Then, for E U, the following statements are equivalent:
(a) F semismooth æ.
(b) F' locally Lipschitz continuous at x, all directional derivatives F F(;x; exist, and for any
G E (F(x+
F + +) Gull o(([=]) as O.
Furthermore, the following statements are equivalent:
(c) F is strongly semismooth at x.
(d) F is locally Lipschitz continuous at v. all directional derivatives F'(x;;) exist, and for any
G OF (x vv,
[F + - F(x) -Gull O(1ee|1) 0.
14

---

Note that the Newton iterates satisfy
||**+| x²11 S
F(t")
(20)
provided that 5² E U. Let rr) denote ball of radius centered x² contained in U, and let m be such that || I m for all x E B(x",r). Assume that F is semismooth at x" and let
n E (0, be arbitrary. Then, there exists p (0,r) such that
- 2 ||||
1 ||||
|F(x +u) - F(x') - G( + e)
(21)
for all |||||| p. Consequently, if we choose r0 such that || x²1 || p. then by induction from
(20 and (21) for xº x², have || - x=' < p and z++1 B(x", p). It follows that the
iterates are well-defined.
Assume further that F is strongly semismooth x², meaning there exists constant such that
[[(x + - F((x)) G(x + b)aol s c||ul|² c. ||l
(22)
for all ||a|| p.
Then, if we choose x such that ||~0~x"1 then by induction from 20 and 22 for v = x x²,
we have
|10 x¹1l ml| ml|F((*) F(x*) G(at"(*** x*)l S m m||a* xx x011?
where | x² p. Therefore, the iterates remain well-defined and satisfy
|| x 11/11> x ==11²==0(1))
indicating quadratic convergence of the sequence {a*} to x=
Thus, under the assumption of strong semismoothness of F,, the Newton iterates converge quadrati- cally to x*.
C Proof of Proposition 2.5
Let
(.....,PP) A = A(w, d) = E (2. ................................ .p......p bن + de > }} I[w = ii (2. .p] |M, + d, S A}.
Since F is continuous, there exists neighborhood N of (w", ,d") such that for all (w, d) E N, have A = {1} UA = (1) U = AA- Since X X has full column rank, X A also has full column rank for all (w, d) sufficiently close to (w",d'
One of the generalized Jacobian Giw,, d) E .FF(,, d) can be expressed as
X XT AuTXI
X
0
XAUI
I
G(w, d)
2 2+2
II
-TA
where =
4i 4dl d:
Note that we can permute the rows and columns of as
4XTXa T x X
LXTX,
L
I Tx
II
G(w,d)=
((2 + !)1 + 0,4) 0
((2
2
1
TA)
17
0
0,
where (xr + XA)=
I1
xx
0
M
M2
Let G
where Mi =
M2 M2 - L1 0).
M3
Ma
15

---

1 XI
xx
0
II
M3 ((2=
1 + +0A)
and M4 =
0
((=
2
TA)
0
1I
0
O0
Note that ((=
- 1 La) -1 is ((1((z
2)1
LA
and
Xx
XI
0
17
00
1I
0
((2
) - IA)
0
((1((z
- ()) LA)
=
I.
I
00
I
1x7Xx
X
XI
We utilize Schur complement of G to formulate the inverse of G. Note that
G/M
M1 MM M3
x - TA
0)
0z
0 (0/(= - })); - TA)
II
LXTX2
0
0
((2 + (() 0w)
17
0
X
0
(ct)))
n
Given that (xxxa) nonsingular, (x] -n ( ): + +Oa) is also nonsingular since |z| < 1
Thus, can construct the inverse of G
(z) (-M7 (MM Ma(G/MMA) MT + +M 3 M
(G/M))-1
-(G/Ma)-
D Derivation of Generalized Jacobian for F
Let F 22 + R.2p be
F(w,d) W1 - prox log(-) ) (wi + di) ( prox 10()))) I x+ V2² +
Sw
w-1 w-1 - T (w + d_1)
where Tx(-) is the componentwise soft ooftthheesholding operator.
prox lox()) (w) + +d11) Ty T2(w-1 +d-))
Let Fi (w, d)- Sw + d and F2 d)
Then, it is straightforward to see that 0f:
S,
0F 3d
I and
Ou
i Orrx
log(() (w1 + d1)
pprox_ log(.) (w) + d1)
W1 dl
+
Ows
dd1
2 (wi dd +
Note that for i= = 2.... .P.
|w; + d:l >
0T, (wi + di)
0Tx(wi + di)
[0
14i + dg = X.
Ow.
Od;
(wi + dil <>
By adopting the convention
8T, (T,(w;+d))
HTs (w)ed,))
0 the kink {lux dil A}. the Jacobians
aF2 Owi
Od;
simplify to aFz as
- J and
-J where
= diag( (Jx) with
3d
1
W1 + d1
wi + dil X Wi di X, i= 2...... P.
J11
+
Jst
2
(w/ + + +4'
Thus,
aF
S - J -3)
OF(w, d)
Orl
which is equivalent to the example of the generalize Jacobian (00).
16

---

E Proofs of Theorem 3.4
E.1 Proof of Proposition 3.2
In Appendix A.2 we established more than the existence of the directional derivative; we demon- strated stronger property that F(x te) F(x) (F: (2;; v) is exactly zero for sufficiently small
t > O. Specifically, the main result in can summarized as follows.
Let (w, d) E RR There exists > 0 such that, for all 0 < t 8,
Su + Rd
(2 1/2) + (2
1/2)vd,1
F(w + tw d tva) F F(w, d)
-Vd
i i A i T i i I i E
Uuii min (wodd -d,
mmxx imaxfeu,i, Ud,,s
The right-hand side defines the component-wise directional derivative F" d; v) and F' d: 2v) = (w,, d: v) for À 2 O. Consequently, F satisfies the approximation property | +v) F(x)
(' w| = of [[l]]] as v - O. Since F is locally Lipschitz (see Appendix A. and the directional
derivative exists at every v, we conclude that F is B-differentiable.
E.2
Proofs of Lemma 3.3 and Theorem 3.4
Proof of Lemma
3.3
For the "only if" part, suppose (Bw, Ud) R satisfies
14
16
We
then observe that T6 is equivalent to
(Ww + wi > 00 Udd + wi = 0) U (Wws + wi = 0n - + wi 2 0),
(ww + wy 00 Ud, + wi 0) U (Wwir +w 0 - w S 0),
which in turn is equivalent to
W.I= + wr+ 2 0, - -d + + 2 0, (tw.I+ + Wi+ -d + + ) = 0,
(23)
- + WI- $ 0, - w/- S (W
w/-) 0.
Rewrite 14
as
Su + vi = SSw d,
1 (WI di) -(ww d1)2
(2 1/2)ow.1 + 1/2) 12,
where z(w1, d1)
w1 +d
2 2.
The second equation can be equivalently written
Ud11 blw.1 +
(24)
where
-1
1 1--2)// ((w) di) (wn +,)
b =
2 + 2
> 0, c= -
2
2
(25)
Let E = Sww - and SR.c denote the submatrix of S with rows indexed by R and columns indexed by C. Then, for the first equation,
S1, + S1,2 + ... + Swpp + Ud,1 22,1 + S2,2 + +... S1 Nw. 0d,2
E
F2
So + + =
Sp + Sp.2 + +... + Uw vd.p
Ep
(Si,1 + b) b) + S1,2 + ... + Su. S2,10 + S2.20
Uwep Ud.2
Eg
:
U, + Sp.20w2 + .... Sp + vd
Ep
17

---

Invoking the definition of Z in 112 we can rearrange this equation
(S1,1 + b) Uw + -w,-1
Ei
S,.1% ST
-4,1- + Vd,A
EA
1w,1 - Rw,-1 U
Ext
S1-,1
ST
-- u1 + Ud
EI
S20,1 +
52.
200 -0 ww,-1 + Ud,2o
Ero
In the above system, transpose the known V,,I0 and Va,A (from of the equation to yield
and
to the right-hand side
(S1,1 + b) + ST. +UT AUTTUU
E, St. Uw -
SA,1V + ST A.AUUIUZ W.IAU
EA S
A. Ud,A
ST,, + ST. AAZ+UZ AUTTUUU + VATT.
Ext I ,0 20
SI S1-,1 + 51- -UZ--U------ AUTZUU- VdII-
Er- ST Wi
S29,1 STO,AUUUUU AUI
EzO So I 1
Now, rearranging the above equation gives the following.
S1,1
S
S1,1+
-1
0 11
0w,1
SA.1
SA.A
SAI+
-AAI-
A,2o
SIt,1
SIt,A
SI+,It
-St,II
σg ,20
-S
-SI-AA
-SI-II+
SI-I-
Oz
Vw.I
$20,1
S2a,4
ST.It
-S
Iz0.ro
=N
Ej ST
1,
(26)
E ST Uu,70 VdiA
0
V
Ex +
S2 5t 10
-Ud
F-- S I0 MI
Ez 5'1
To ro WwiI0
Since X a is full-rank,
S1,1
SI.A
S
S1
SA.1
SA.A
SATI
SAT
= =11nn)
XIUAuI+UZ-)
(27)
ST
Sx
ST
II1
ST II
-T-,11
SI-,A
SI- II+
SI I--
is nonsingular (in fact, positive definite). Denote the 4 x principal block submatrix of N by
S1,1 +
S1.A
St S SSAI
SA11
SA.
SA.T
ERAAÍIIĀ
N₂ =
(28)
Sz+,1
SI+.A
SI+
-St+I I--
SI-,1
-SI-AA
-SI-It
SI-22 I-
N
so that N -
Cc
with
ITe To
C= (Szo,1 Szo S. I+, - -
(29)
Then N₁ is also nonsingular since O, and thus
Ni1
0
N-1
-CN-1
Izo
Therefore,
1w,1
0
VerAA
0
ST + Uy = E +
Vw,I+
N -1
1
P
U.II-
(30)
UdZ
Ud.20
0
18

---

S1 +
SI.AA SA.A)
let the 2 2 principal diagonl block submatrix of N1 be N₂
then the
SA.1
bottom right 2 x block Ni 1--1ip is the Schur complement of N2 in N1:
S 1I+ (
-1t,I-) SI-1-1
S -S1,11
Sx-,d -SL-Ad
[S1zz -51,I-
I
M =
SI-I+ --S--x
N2¹
(31)
SAI+
-SAII
Thus,
St Ud
I UwwI
M
Vd,I+ UdII
INT
(32)
where [N[ la:4 refers to the submatrix of N1
spanning the third and fourth row blocks.
From the property of the Schur complement, M is positive definite if and only if N1 is positive definite. To see this, write
S 111+b +6
SI.A
S1.1
S1I
SA.1
SA.A
SAII
SAI
N1 To
To:
ST
SI
STI I1
SI SI-II-
ST-,1
SI-,A A
ST- I+
where
11,1
IA.A
Io
II+
is symmetric and orthogonal. The middle matrix is
27 that is positive definite.
If we let
-Ua,I + WI+
Ud WT- ) y
U.IIT WTT
(33)
-U,,I- - WT
then we see
(x, =0 and y Mx + 9
where
M
WI+ 1----
+ [N 1sa
Wh+ --wz-
(34)
from 23 and 32 In other words, the x and in (33) jointly solve the LCP (77 with the M and
q given in (31) and (34). fact, (17) has unique solution, since dimensiimnniional LCP of the form (17) has unique solution if M is P- P-matrix, meaning all of its principal minors are positive. symmetric positive definite matrix a P-matrix.
For the "if" part, suppose yx solves the LCP (17) with the M and given in 31 and yx, y) is unique. Let
Then
-x+
+
WH+ WI-
for
XI-
(35)
Cd.II
XI-
X
and
0
Uw
N-1
Va.I
-UdII
3:1
Then
UW,IT -W,I--
N-1 P
Xx+ + Wht
-X
WI-
3:4
lanpP -M
(wIt]
Mx = Mx
(-wz-)
=y- - (
-(-wr-
(wwy)
19

---

which implies 14
16
Proof of Theorem
3.4 From the uniqueness of the solution to I.P 17) the wh in 11 is well-define::
We know
-Ho,
1(1
- 1 d
w
from 8
and 12 That
XI+
d
follows from 350 the M and for solving (x, in LCP 17 are given in (31 and (34). Also, from (4
bmu
ch,
d
Finally, with the above quantities computed, we have where & and e* can be deduced from from 30 and (26).
25
N-1
pk
uk
I
u d,z9
with
E
ST .ookk
ST
A,zo uh wi70 ck
E
4,I0 uk
N-1
9) and pk
d.A
N-1
-CN,
E
ST
It I,,zo th
w.20
- + ST
.70 vit
Ex. S10 Io u 2,
WiI0
where N1 and C' are given in
28 and 29 and E -Sw dt.
F Proofs of Theorem 3.8
F1 Proof of Theroem 3.8
To verify the global convergence of Algorithm ! we utilize some regularity conditions of F in 28 25 Consider following lemmas:
Lemma F1. The proximal gradient step in Algorithm Joccurs only a finite number of times.
Lemma F22 There exists a constant >> O and a neighborhood N of x² such that for all æ E N.
[I; w)ll 2 L|
(36)
holds for all vectors u P.Pp. Lemma F.3. For every sequence f converging to xt, {ur} converging vf, and positive scalar sequence {1x) converging to zero,
0( tht*) 0(x*)
lim sup 1-000
lim
Proof of Theorem 3.8
We exclude the trivial case where x exactly hits the optimal solution during
the iteration, in other words, 0(**) > O. First, we show that my in (18) always exists for all Newton iterations. To this assume that
0xx* + p A(x²) ~2op"0(x*),
(37)
holds for all E {m;} where m, is an infinite increasing integer sequence. Then, dividing both sides by nm and taking the limit as j + xc, we have
0x + , 0xx²)
-200(I & lim +0000
F' ('( = -20(2*).
ms
20

---

as we assumed that a E (0, 1) and 8(x²) > O, we have a contradiction. So, there always exists my such that (18) holds. Since (0(x*)) } is strictly decreasing sequence bounded above zero, its limit exists and therefore,
lim {0(x²) - 0(z 11 2 20 lim siimsupp===((a"" O. k
1-+00
If lim supk- n O, then limg->o 0(x²) O. Since we assume that I+ UI- 0 and X has full column rank, every (in particular, the unique) generalized jacobian OF(x*) is nonsingular. Then, by Clarke's is inverse function theorem [], F admits locally Lipschitz continuous inverse in a neighborhood of x². Consequently, for sufficiently large k, we have x² = F ==
F-1 (0) x²
If lim sup p = 0, then we have pm 0=)((x²) + 0 and
0(. + 4pm p=)=-1ok) 0(x²) > -2ap"M
Consider a subsequence of fo* k E }} that converges to xt such that At = A For all k E K large enough, æ* E N. where N is the neighborhood of x specified Lemma F.2 By Lemma F.2 it follows that ||e* S +||=(>>))|, which implies that the corresponding subsequence {u²) keEK is bounded. Thus, we can find a subsequence k1, l2 such that x²t and S converges into some x! and vi, respectively. Then, by Lemma
E.3 we have
0(r + pmne (mme-1ph)) ((aht) 2 -20((x))
lim F( F' ; 1-000
lim
pm
On the other hand, since tt is Newton direction for 11
lim F(x"t)? 1-000
(acku)TFFF F (x*tt = lim -20(x=)) -20(xx))
10
Therefore, 0(x²) = 0 and since the solution is unique, the accumulation point xf coincides with x".
Finally, we establish the convergence of the full sequence {-*}. Since x is the unique solution to
F(x:) 0, there exists 8 > such that x* is the unique minimizer of [[)(lll for | ==11| S 6.
Recall that for all & large enough, | < +=(((=)] from Lemma F22 Since [|F(x*)]I + 0 as
k oo, it follows that ||u*|| + 0. Therefore, there exists constant 8 (0,8) such that, have
||a x²1| S | 8 8'
whenever ||* x=11 8'. Let μ denote the minimum value of [[F(z)] over the region 8' < || x11| S 6. Suppose that to x= < 5' and | )l| < μ. Then, since the sequence
[==()))] I is strictly decreasing, it follows that [2 ||| < But, since ||z =0 < 8 8'.
we also have that
||o - xx'1l S ||a x || x < 6.
Hence, by the definition of we have ||2:001 -=1 s 8'. By induction, it follows that - --=1= 8' and || [[ < for all k 2 kn- Since x" is the unique accumulation point of {x*} within this neighborhood, we conclude that the entire sequence {a*} converges to x".
We now characterize the local convergence behavior of the iterates. By [30 Theorem 4.3], since x -> x², the unique solution F(x) = O, and F"(::)) :) in (12) semicontinuous degree 2 30
Lemma 2.3], the sequence (a*) converges quadratically x" in its neighborhood and pd eventually becomes 1.
F.2 Proof of Lemma M
The proximal gradient algorithm of 199 generates iterate sequence to converge linearly to the unique solution x². We distinguish between two contexts in which the proximal gradient step is invoked: the inner loop (lines 16-16) and the outer loop (lines 12-17).
First consider the inner loop. Let (x denote the sequence of iterates computed by the proximal gradient update. This update is applied iteratively from the initial point 2. x* Assume
(2 0. (Otherwise we done.) Since xprox x as j 0o, ((ppiox) + 0 j 30, and
21

---

[F d: w) inf B by linear convergence, there exists finite index J such that 0(x S ((x Thus, the inner loop terminates finitely many iterations.
Next, consider the proximal gradient phase of the outer loop. This phase is potentially invoked until the iterate enters neighborhood, N of x= (w" dd) where A = A holds (such neighborhood exists due to the assumption of Theorem 3.8p. If the iterate does not enter N after a finite number of invocations, then (xxptxx) A O, contradiction.
Once the (outer) sequence {a*] enters N, within this neighborhood the rank condition for local B-semismooth Newton convergence is satisfied, so the algorithm remains in the Newton phase. Consequently, the proximal gradient step is invoked only a finite number of times throughout the entire procedure.
Proof Lemma
E2
Let
A E (2. .......77 . w0 + d X}, I= E (2. ... .P) wm + d S A}.
Recall that the Bouligand derivative F' d; v) can be expressed with the generalized Jacobian Ga(w, d) as
X XT XuuXI 2+
X XAUZ
1
0
X7 XAuz
I
F' ("(w,d:0) GA(w,
V
0
0
II
0
TA
for the differentiable points (lw, + del # A) where 2 = 2 (d: Since we know that the
wy dd
value of the (p + i) -th element of (w, d: v) for i E I= is Dd,t or owii we can construct the following inequality:
d)ull
Thus, suffices to show that the generalized Jacobian Ga(w, d) is nonsingular. Then, it can be deduced that
| (w, (w, d) d)ull 2 ||=||.
and thus, L in (36 holds for the minimum value of ||G (w, d)ll. Since F is continuous and It -U----- 0. there exists neighborhood A of (w" dd) such that for all (w, d) N, we have A= Ar Since A. has full column rank, x also has full column rank for all (w, d) sufficiently close to (w",d") d. We already proved that such G has an inverse in Appendix
F.4
Proof of Lemma 3
Note that
(x + txo") - (x(
lim 8-000
< lim F(a F' (xtt:
(38)
21x
k
holds as the equality within its domain of x for 1,... PP + 1, and the equality also holds for
i p ..... 2p where |uل + ddi # X (j= 2.... .P). Ft' are differentiable and Ft's continuous at these points. Hence, it suffices to prove that (38) also holds for non-differential points w++ d) = x. Suppose that (w* dd²) converges to (wt, dt) and (u o/)) converges to (ui, vd).. Without loss of generality, let wt + di = A. Let k1, k2.-.. be any increasing sequence of a positive integer. Then, for any {p- + 1,... .2p). we can find some subquence kim, where
Fmm + Ekmn Wis
Km dleme + Ex V kmm)) F, Fm (dkms)))
lim
(39)
tRim
exists and falls into at least one of the following cases:
22

---

1. mms dlimy > A and whi + dfiins Uierss
(1/w
Rome
Amm
2 X. Then,
kmn
Komz
Renis
39
lim 1-y00
Uw
t tomg
(Uw
+ vd
2/ lim F' (xrMin emms))
tkmm
2. wms climy < A and whis + dfime Wieess (1w
Rome
Remm 2 A. Then,
Komz
(wmtmm + cless + tkms (Uw (owm
Kims A)
t &mmy
Uw
+ Od
39
lim 1-500
tkmy
S
lim (tkma Uw 1-500
Amms tkkis
lim F" (Itims olim)).
10000
3. mmi + d > X and wime + drime + thmmy (wi
S Then, Roms
vd
limm
kmy Tw
+ (wmtime + dlims - 2)
lim 10000
tking
Kmy Ow
Rims
Boms
lkmms
(Uw
+ 1 -- 1+20
lim 1+000
lim ;
lkms
4. wmi + k < A and wm; + drime + Uyomy (w
Rome
Romn Ud
S A. Then,
lim (tkme Uw
Kmy ) t/km U/
lim F (xtttme whtmss ). 1-000
5. wm. + dk X. Then,
lim (tkm, Fiw
kmp tkm max km + vd
Kmy 0})} tt)mm
11000
min (ut, -nd) S lim F/
1+00
Hence, conclude that
Fa + (100) - Fi(x*))
lim sup
lim (0 k5500
k+200
holds for Then, using the continuity of F it follows that
0xx txx*)) 0(x*)
lim sup
lim F(a*)(F'((ztto" 14000
k-000
with some simple arithmetic.
G Additional Experiments
G.1
Ablation study: roles of the proximal update and LCP solving
Algorithm Srrrithhiinnlludsss inncludes two components other than the semismooth Newton steps, that are the proximal update steps and the LCP solve step. To disentangle the roles of these two components, we replicated the setup of Section 1.1 and ran each configuration for 50 iterations under three ablations: (i) disabling only the proximal update steps, (ii) disabling only the LCP solve step, and (iii) disabling both. The LCP solve step is invoked only when, for given row, I+ UI + 0. In this ablation study, "disabling the LCP solve step" means reclassifying all indices in I+ UI- as inactive by assigning them to T thereby bypassing the LCP solve subroutine. Recall that with both components enabled, Algorithm terminated within 30 iterations.
At X 0.10. disabling the proximal steps prevented the algorithm from converging within 50 iterations, regardless of whether the I.CP solve component was enabled. By contrast, convergence was achieved whenever the proximal gradient component retained, even without the LCP solve
23

---

comeponent. In the original run at this value of X, the proximal update was invoked in 229 out of 1000 rows in the first iteration, indicating its importance in the denser regime. For À 0.15, 0.2. the algorithm converged under all ablation settings. In the original runs, at most one proximal update (out of 1000 rows) occurred in the first iteration, after which the semismooth Newton updates predominated. Thus, while the proximal gradient steps can be dispensable in these sparser regimes, it remains a practical safeguard unless the sparsity level is known a priori.
The LCP solve step was activated only on the kink set {i: w + dil = }. Since exact equality is numerically rare, such activations almost never occurred. Consistent with this observation, Table shows that Algorithm [ is largely insensitive to this event, and Table S reports that the LCP solve occurs for negligible fraction of rows. Even with 10,000 Newton row updates per iteration, the LCP solve occurred in about one or fewer row on average, as reflected by the small average and maximum per-iteration fractions. Thus, the practical impact of the LCP module is limited, although remains theoretically indispensable for correctness and robustness, especially on unseen data.
Full method
Proximal
No LCP
No Proximal No LCP
0.10
Converged (30)
Diverged
Converged (30)
Diverged Converged (26)
0.15
Converged (22)
Converged (26)
Converged (22)
0.20
Converged (14)
Converged (14)
Converged (14)
Converged (14)
Table 4: Convergence within 50 iterations under ablations of the proximal and LCP modules.
A
Avg. frac. LCP / iter.
Max frac. LCP/ iter.
Last iter. LCP occurred
Total iters.
0.20
0.00005
0.0011
62
159
0.45
0.00004
0.0012
29
69
0.60
0,00012
0.0045
31
63
Table Incidence of LCP updates in the original runs. "Frac." denotes the fraction of rows updated by LCP in an iteration.
G.2 Sensitivity to the choice of regularization parameter
We assessed the convergence behavior of Algorithm across range of regularization parameters under the same experimental setup as described Section [.3] focusing on real-world datasets. The baseline choice À = 0.450 was selected by the extended Bayesian information criterion in [19]. To examine sensitivity, we varied À from 0.325 to 0.575 with an increment of 0.025 and recorded the iteration count, the number of nonzeros in the estimated precision matrix, and the resulting fraction of nonzero elements (density). The result summarized in Table
The change of the density remained modest within the range À E (0.4, indicating that the
sparsity pattern relatively insensitive around the "optimal" À = 0.450. Across the tested values, A Algorithm convergee A within a few tens of iterations; while the iteration count tended to increase mildly as λ decreases (i.e., as the solution becomes denser), the overall convergence behavior remained stable.
G.3 Convergence across densities and alternate random graph topologies
In this section, we further investigate the trade-off between graph density and convergence using alternative random graph structures. In Section 4.2 the ground truth precision matrix O was constructed based an Erdös-Rényt graph, where edges are sampled uniformly random. To capture a broader variety of graph structures, we additionally consider two alternative random graph topologies: hub networks and scale-free graphs.
hub network was generated using the following procedure: (1) for each cluster, (300 x density) hub nodes were selected, each of which was connected to 66 randomly chosen nodes; (2) the remaining (100- 300 x density) nodes within each cluster were connected via an Aadös-Rénys graph, using the remaining intra-cluster edges.
24

---

Iterations
Nonzeros
Density (%)
0.325
103
5,494,644
5.89
10-3
0.350
69
4,806,952
5.15 x 10
0.375
78
4.214,966
4.52 10
0,400
67
3,699,259
3.96 x 10
0,425
94
3,248.215
3,48 10-3
0.450
2,852,968
3.06 10
0.475
67
2,505,706
2.68 x 10
0,500
2,200,810
2.36 10
0,525
56
1,930,194
2,07 X 10
0.550
66
1,692,240
1.81 10
0.575
65
1,481,495
1.59 x 10
Table 6: Iteration count, number of nonzero elements, and the fraction of nonzero elements in the estimated precision matrix for the LIHC dataset with varying regularization paremeter A.
A scale-free graph was generated so that the degree k of the nodes has probability P(k) ~ k-2.3
Each graph comprised 10 clusters with 100 nodes per cluster. The total number of edges was determined by the graph density defined in Section 4.2 with 80% of the edges allocated to intra- cluster connections and the remaining 20% inter-cluster connections. The resulting graphs are summarized in T
Density
Hub Nodes
Intra-cluster Edges per Cluster
0.03
1198
0.05
15
1998
0.10
30
3996
Table 7: The number of hub nodes and intra-cluster edges per cluster according to graph density
chosen Following the experimental setup described in Section 4.2 the regularization parameter so that the estimated precision matrix exhibits sparsity level similar to that of the true matrix O We then assessed whether Algorithm converges within 300 iterations across various graph types and density levels (0.03, 0.05, and 0.1). The results are summarized in Tables 8
We observe that the convergence behavior in different graph density levels is primarily determined by the overall density of the true precision matrix (i.e., graph structure), rather than by the specific graph structures. Intuitively, although the graph structures considered in this experiment may exhibit diverse local structural patterns, their overall density was controlled to remain constant across all settings. Since the applicability of the B-semismooth Newton method depends more on overall graph density than on local structural variations, the observed results are consistent with this interpretation.
Hub Network
Density \ A
0.2
0.15
0.1
0.09
0.08
0.07
Density \
0.07
0.06
0.05
0.03
Y
Y
0.05
Y
Y
0.1
3%
83.3%
100%
0.1
Y
N
Table 8:
Convergence within 300 iterations for varying graph density and À
Table 9: Proportion of rows updated by proximal method after 300 iterations
25

---

Scale-free graph
Density
À
0.15
0.1
0.09
0.08
0.07
Density
A
0.07
0.06
0.05
0.03
0.05
Y
0.1
1.6%
85.5%
99.8%
0.1
Y
N
Table 10: Convergence within 300 iterations for varying graph density and À
Table 11: Proportion rows updated by proximal method after 300 iterations
NeurIPS Paper Checklist
1. Claims
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Question:
Answer: [Yes]
Justification: All claims made in the abstract and introduction are justified with theorems and proofs.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [No]
Justification: We do not provide explicit discussions on the limitations of our method as they align with the original method in
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested few datasets with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss.
26

---

While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer:
[Yes]
Justification:
The theoretical results of the paper provided with assumptions and proofs.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proof's provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Answer: [Yes]
Justification:
The code for the main experiments is provided.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
27

---

(c) If the contribution is a new model (e.g.. a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described supplemental material?
Answer: [Yes]
Justification:
Sufficient access to the data and the code is provided.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips. ec/
public/guide Satia//guuidos/oooeuubibsssssss wwssonooiicyyy ollicyy for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/ htips.cccpublic//eee Seeguidddde/odSSubmmm /wtp///ooeeSeseeeeeeeeeeee wetpimiissioooooiiiiccoc
for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification:
The experiment details are specified.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [No]
Justification: Reporting error bars on the experiments is not relevant with our research.
28

---

Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by enror bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the eπror bar is the standard deviation or the standard error of the mean.
It is OK to report I-sigma error bars, but should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful figures symmetric error bars that would yield results that are out of range (c.g. negative error rates).
show in tables
error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
Question:
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Answer: [Yes]
Justification: Computer resources for the experiments are specified.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't 't make it into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics httppp///weerrppp...ooo/////eeea
Answer: [Yes]
Justification:
We have reviewed the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. there special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive socictal impacts and negative societal impacts of the work performed?
Answer: [NA]

---

Justification:
This research discusses an optimization algorithm with no societal impacts.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negatíve societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of MI.).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (c.g., pretrained language models, image generators, scraped datasets)?
Question:
Answer: [NA]
Justification:
The release of data or models does not have such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification:
The data used in the experiment is properly credited.
Guidelines:
The answer NA means that the paper does not existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
30

---

The name of the license (e.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, palasets.papersstt paperswi thcode. com con/datases has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is available online, the authors are encouraged reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: This paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
This research does not involve crowdsourcing nor human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
This research does not involve human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
31

---

Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if' applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification:
The core method of this research does not involve L.Ms..
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy for what should or should not be described.
32


Reviews:
[{"review_id": "QPD5FQhuEv", "reviewer": "Reviewer_r5ah", "summary": "This paper proposes a new semismooth Newton method for efficiently estimating partial correlation networks using an ℓ₁-regularized pseudolikelihood-based framework. Building on the existing ACCORD framework (which itself is scalable but relies on slower first-order methods), the authors derive a second-order optimization algorithm that significantly improves convergence speed and reduces reliance on expensive hardware like supercomputers.\nThe method is applied to both synthetic and real-world datasets, showing improved performance over the state-of-the-art proximal gradient method (ACCORD-FBS) in terms of convergence rate, iteration count, and wall-clock time.", "questions": "1.  How sensitive is the method to the choice of regularization parameter λ, particularly in real-world datasets?\n2. Have you tested the algorithm on non-Gaussian data distributions? How does performance compare in those settings? Include experiments on non-Gaussian synthetic data (e.g., t-distributions or contaminated normals) to validate robustness.", "limitations": "yes", "rating": 3, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n1. The authors develop a semismooth Newton approach with a strong theoretical foundation.\n2. The algorithm is shown to have locally quadratic convergence, which is a major improvement over linear convergence of previous methods.\n3. The paper provides clear derivations, conditions for global convergence, and strategies for damping and line search.\n\nWeekness\n\n1. The method assumes submatrices like X_A have full column rank. While reasonable in many high-dimensional biological datasets, this may not hold in low-sample or noisy regimes. Include a robustness analysis or empirical study on how violations of this assumption affect convergence and solution accuracy would be helpful\n2. One motivation for pseudolikelihood methods is robustness to non-Gaussian data, but this is not empirically evaluated in the paper.", "was_revised": true, "final_justification": "I cannot access the contribution of an efficient second-order algorithm for ACCORD, given that ACCORD [19] has not shown its wide applicability and correctness yet, which is arXiv paper beyond the reviewing scope of this paper.\n\nTo clarify, I do acknowledge the technical difficulty of developing a second-order solver that is more efficient and with provable performance for ACCORD, that's why I initially gave 4. But I find that difficulty does not necessarily mean a significant contribution, given that the significance of the original algorithm ACCORD has not been shown yet, which is what I did not realize before the rebuttal. That's why I lowered my score to 3. However, my confidence score remains low as before.", "rebuttal": "### Answer to Weakness 1\n\nIn fact the desired robustness study has been already presented in Section 4.2. We conducted systematic experiments by varying both the graph density (ratio of the number of edges in the true graph to that of the complete graph) and the regularization parameter $\\lambda$; see Tables 1 and 2. These experiments allow us to examine how the size of the active set—and by implication, the likelihood of rank-deficiency in $X_{\\bar{\\bar{\\mathcal{A}}}_k}$— which in turn affects the applicability of the semismooth Newton step and overall convergence.\n\nThe results reveal a clear pattern:\n\n- Our algorithm converges reliably within 300 iterations in the majority of settings;\n- Failure occurs when the true graph is dense (e.g., density of 0.10) and the regularization parameter $\\lambda$ is too small to sufficiently regularize the solution. This leads to a large active set, increasing the chance of rank-deficiency of $X_{\\bar{\\bar{\\mathcal{A}}}_k}$  and preventing semismooth Newton steps from being applied.\n\nTo further assess robustness under more complex structural variations, Appendix F presents additional experiments on hub networks and scale-free graphs, which exhibit diverse local connectivity patterns (Section 4.2 is about uniformly random graphs). Importantly, we observe that convergence behavior remains consistent across these graph types when the overall graph density is held fixed. This strongly suggests that the convergence of our method—and the applicability of the semismooth Newton step—is primarily governed by global sparsity (active set size) rather than local graph topology.\n\nTaken together, these findings provide empirical supports that our globalized algorithm remains robust in realistic settings as long as the model is appropriately regularized. \n\nObtaining a statistical guarantee of continuation strategies *a la* [7] under noisy regime is challenging even in the regression setting [15]. Failure to converge within, say, 300 iterations, and/or dominance of the proximal gradient steps in our globalized algorithm as examined in Section 4.2 and Appendix F may suggest a stopping criterion that adapts to the sparsity level of noisy data, while further study is warranted. Thank you for the insightful comment!\n\n### Answer to Question 1\n\nPer your suggestion, we have investigated the convergence behavior across a range of regularization parameters $\\lambda$ under the same experimental settings described in the paper, focusing on real-world datasets.\n\nSpecifically, for the LIHC dataset, we used the value $\\lambda = 0.450$, selected via the extended Bayesian Information Criterion proposed in [19]. To further assess the sensitivity of our algorithm with respect to the choice of $\\lambda$, we systematically varied $\\lambda$ from 0.325 to 0.575 with increment of 0.025.\n\nThe table below presents the number of iterations, the number of nonzero elements in the estimated precision matrix, and its corresponding density. \n\n| λ | Iterations | Nonzeros | Density (%) |\n| --- | --- | --- | --- |\n| 0.325 | 103 | 5,494,644 | 5.89 × 10⁻³ |\n| 0.350 | 69 | 4,806,952 | 5.15 × 10⁻³ |\n| 0.375 | 78 | 4,214,966 | 4.52 × 10⁻³ |\n| 0.400 | 67 | 3,699,259 | 3.96 × 10⁻³ |\n| 0.425 | 94 | 3,248,215 | 3.48 × 10⁻³ |\n| 0.450 | 64 | 2,852,968 | 3.06 × 10⁻³ |\n| 0.475 | 67 | 2,505,706 | 2.68 × 10⁻³ |\n| 0.500 | 71 | 2,200,810 | 2.36 × 10⁻³ |\n| 0.525 | 56 | 1,930,194 | 2.07 × 10⁻³ |\n| 0.550 | 66 | 1,692,240 | 1.81 × 10⁻³ |\n| 0.575 | 65 | 1,481,495 | 1.59 × 10⁻³ |\n\nBased on the leading significant digits of the estimated sparsity levels (density), the variation remains modest within the range $\\lambda \\in [0.4, 0.5]$, suggesting that this interval corresponds to a region where the final sparsity pattern is relatively insensitive to the choice of the “optimal” $\\lambda$.\n\nWe also observed that our algorithm consistently converges within a few tens of iterations across all tested values of $\\lambda$. Although the number of iterations tends to slightly increase as $\\lambda$ decreases, the overall convergence behavior remains stable. \n\n### Answer to Weakness 2 and Question 2\n\nWe investigated the performance of ACCORD on simulated data with $n=500$ and $p=1{,}000$, generated from a multivariate $t$-distribution with degrees of freedom $\\nu = 5$, where the graph structure of the precision matrix $\\Theta^*$ was defined by an Erdős–Rényi graph with 1,000 randomly selected edges. We present the overall result from 50 replications following [19], including the mean (standard deviation) of the Area Under Precision-Recall Curve (AUPRC) obtained from multiple values of $\\lambda$ and number of True positive (TP) / False Positive (FP) edges detected for the selected $\\lambda$. The $\\lambda$ was selected with the extended psuedo-Bayesian Information Criterion (ep-BIC) for ACCORD and the plain BIC for GLASSO, following [19]. The result is as follows.\n\n| Method | AUPRC | # TP edges | # FP edges |\n| --- | --- | --- | --- |\n| ACCORD | 0.722 (0.082) | 704 (28.9) | 264 (277.9) |\n| GLASSO | 0.695 (0.098) | 698 (26.6) | 399 (712.9) |\n\nWe could observe that by comparing AUPRC, ACCORD showed better overall performance on detecting the precision graph. Comparing the graphs selected with (ep)-BIC, ACCORD detected significantly less false positive edges than GLASSO while both method detected similar number of true positive edges. \n\nWe note that the statistical performance of the ACCORD model is not our primary focus of the present manuscript, as the aim of the study is to develop a second-order algorithm that solves the ACCORD optimization problem (3) order of magnitude faster than the proximal gradient method proposed in [19]."}, {"review_id": "L6J2VgkB6b", "reviewer": "Reviewer_fHCF", "summary": "This paper studies partial correlation network estimation by applying graph lasso within the recently proposed ACCORD framework (Lee et al., 2019). While the method in Lee et al. (2019) is computationally scalable, it faces a bottleneck for graph lasso due to sparse-dense matrix-matrix multiplication. The authors further explore the inherent scalability of this framework by employing semismooth Newton methods and propose a scalable second-order algorithm that achieves superlinear convergence. Theoretically, they demonstrate local quadratic convergence. Simulations and real data applications comparing their method with the proximal gradient method in Lee et al. (2019) are also provided.", "questions": "1.  While the proposed method and algorithm does significantly improve the cur- rent algorithms based on glasso, it is known that glasso has its own problem in partial correlation network estimation. For instance, it relies on assumption of equal regularization for all nodes but this is not appropriate for real network, and the performance is not good as it may retrieves very small false edges and sensitive to outliers. See some references here:\n\n[1] Huang. S, Jin. J, & Yao.Z. (2016). Partial correlation screening for estimating large precision matrices, with applications to classification.” Ann. Statist. 44 (5) 2018 - 2057.\n\n[2] Cheng, C., Ke, Y., & Zhang, W. (2025). Large Precision Matrix Estimation with Unknown Group Structure. Journal of the American Statistical Association, 1–24.\n\nGiven this, how can you justify such improvement really contribute to the study of partial correlation network estimation? People are more likely to use other methods that has a fair computational cost but more accurate performance.\n\n2. The locally quadratic convergence result is quite attractive. How can we heuristic under stand it? What is the intuition behind this quadratic convergence?", "limitations": "Yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strength: The paper is well-organized with some theory to support their method and demonstrate the phenomenon appears in simulation. The improvement compared to Lee et al., 2019 is significant.\n\nWeakness:\n\n1. The writing of theoretical part is hard to understand. More explanations on the lemma and theorem and clearly state the connection of the lemmas could help improve the readability given the abundant notations crowded in lemmas or theorems. \n\n2.  The numerical part only compare their method with the one in Lee et al., 2019. How about other algorithms and also the algorithms beyond glasso framework? In addition, only convergence rate or computational cost is compared, nothing is\nreported about accuracy. \n\n3. The motivation of the paper is this: Glasso is a great method, therefore, we wish to improve it. However, despite the popularity of the glasso, study on the glasso is largely theoretical, with limited numerical success for high dimensional data (simulated and real). In fact, it was reported that the glasso estimation of the precision matrix tends to be overly sparse, and alternative approaches may have much better performance.    My question is,  if the glasso is not an effective algorithm, how to justify the work that tries to improve its (and only its) computing side? It seems that the more important problem is to improve its accuracy.  In fact, the paper only investigate the computing speed of the proposed method, without investigation on the estimation accuracy.", "was_revised": true, "final_justification": "I would like to thank the authors for prepare a carefully-written rebuttal and a carefully-written response to my most recent questions. \nI wish to clarify that, by accuracy, I meant numerical experiments. The authors pointed to some papers on GLASSO on the literature, \nbut those papers focus on theory. My experience is that, numerically, GLASSO is quite poor, compared to other methods. Overall, I think the paper focus on computation side, but does not adequate address the accuracy the procedure. I would encourage an extensive numerical study, where the authors can compare with representative methods in the literature. I will keep my rating.", "rebuttal": "### Answer to W1\n\nWe appreciate your comments on the clarity of our theoretical development. We acknowledge that the presentation of Lemma 3.3, Theorem 3.4, and Theorem 3.8 may have been hard to follow due to dense notation and limited space. We welcome the opportunity to clarify the logical structure and agree that more detailed explanations would improve readability.\n\nThe theoretical foundation of our method is presented in two parts:\n- In Section 2, we study the local convergence of the semismooth Newton method under the standard framework based on generalized Jacobians. In particular, Theorem 2.4 shows that if a generalized Jacobian is nonsingular near a root of a (strongly) semismooth function, then the Newton iteration converges locally quadratically. To ensure this condition in our problem, Proposition 2.5 provides a sufficient condition based on the rank of a submatrix of the design matrix, which is linked to the model identifiability.\n- Section 3 extends the analysis to ensure global convergence of the semismooth Newton-type method when applied to the  ACCORD objective. While the local theory guarantees fast convergence near the solution, it does not inform us how to choose a suitable initial point. To address this, we propose a damping strategy and work with the Bouligand (B-) derivative, a specific choice of the generalized Jacobian. This leads to the development of the B-semismooth Newton method, which is amenable to globalization.\n\nWe summarize below how the key lemmas and theorems connect within this framework:\n\n- Proposition 2.5 provides a sufficient condition for the nonsingularity of generalized Jacobians, thereby ensuring local quadratic convergence per Theorem 2.4.\n- Lemma 3.3 reformulates the B-semismooth Newton direction as a solution to a linear complementarity problem (LCP). This result is key to formulating the Newton step as a tractable numerical subproblem: solving the B-derivative equation is only a descriptive step in general B-semismooth Newton algorithm. Here we provide a constructive way to solve it - solve a well-defined LCP in reduced dimension.\n- Theorem 3.4 shows that the well-defined LCP admits a unique solution, so the B-semismooth Newton iteration is indeed well-defined.\n- Corollary 3.5 confirms local quadratic convergence of the B-semismooth Newton iteration under the same rank condition as Proposition 2.5.\n- Theorem 3.8 establishes global convergence of Algorithm 1. Even when the Newton direction is undefined (e.g., due to local rank deficiency), the algorithm uses a proximal gradient algorithm, which occurs only finitely many times. Eventually, the method enters the region where local quadratic convergence (guaranteed by Theorem 3.4 and Corollary 3.5) takes over.\n\nIn the revised version, we will add brief guiding remarks before theorems and propositions to explain their roles and improve the readability of the theoretical development.\n\n### Answer to W2\n\nWe would like to clarify that the scope of the present manuscript is enhancing the algorithm for solving the ACCORD optimization problem (3) proposed in [19]. ACCORD is a pseudolikelihood method designed for scalable estimation of partial correlation networks, distinct from the the graphical lasso (glasso) framework in learning partial correlation graphs, and their comparison with respect to convergence property and model performance is already presented in [19] along with other preceding pseudolikelihood-based models. Regarding the improvement of the solution algorithm for ACCORD, however, no alternative beyond the proximal gradient method in [19] has been proposed, to the best of our knowledge. As such, we  focus on comparing the computational performance between this baseline algorithm and our semi-smooth Newton approach. If there is no time and memory limits, then both algorithms produce the same optimal solution for the same data, which makes the present comparison meaningful.\n\n### Answer to W3\n\nAs we pointed out above, glasso and ACCORD are different frameworks/methods for learning partial correlation graphs from data. As we mentioned in Introduction, ACCORD is closer to pseudolikelihood-based methods such as SPACE [27] or CONCORD [17] than glasso, which are known to be more robust to non-normality. Also, It is *not* that the glasso estimates are overly sparse and the statistical performance of glasso is bad, but that it is not *computationally* tractable in modern, ultra-high dimensional applications, say when $p > 100,000$. ACCORD [19] is a recently proposed, viable alternative to glasso, yielding a  tractable and scalable optimization problem in computing the estimate with a similar statistical performance to glasso. While the statistical properties are studied in depth in [19], in our view, the optimization algorithm in [19], which is a first-order proximal gradient algorirhtm, is not sufficiently efficient. So the present manuscript focuses on the computational aspect of the ACCORD framework, and exploits the semismooth nature of the optimization problem to produce faster and more scalable algorithm. We demonstrate that our algorithm reaches the same solution as [19] within substantially fewer iterations and lower runtime—an important practical advantage for large-scale partial correlation graph estimation. For the statistical performance of ACCORD estimator, including the accuracy (estimation error) measured in various norms, please refer to [19].\n\n### Answer to Q1\n\nAgain, we emphasize that our contribution lies in improving the computational performance of the estimation algorithm, rather than in the methodology itself. Nonetheless, the “equal regularization for all nodes” problem can be easily addressed by introducing an individual regularization coefficient (hyperparameter) $\\lambda_{ij}$ for each *edge*. The theory of [19] applies with essentially no modification as long as $\\lambda_{ij}$’s diminish with $n$ and $p$ at the same rate. Refinements like SLOPE (Bogdan et al., 2015) to control false discovery rates can also be made easily. \n\nThe method of Cheng, Ke, and Zhang (2025) is essentially the symmetric lasso loss (Friedman et al., 2010) with sparse-group lasso penalty (Simon et al, 2013), and the focus of the paper is how to find groups from data. As pointed out in [17], symmetric lasso is a nonconvex loss function and its convergence is problematic. In contrast, the ACCORD loss is convex and the sparse-group lasso penalty can be easily imposed to make the overall optimization problem convex. In many multi-omics studies, groups are known a priori in biological basis. Statistical analysis of the resulting estimator, however, requires a distinct set of assumptions, which is common in the sparse learning literature (see Wainwright, 2019). \n\nThe method of Huang et al. (2016) requires frequent access to *all* the principal submatrices of the sample covariance matrix. When $p > 100,000$ as we focus, constructing the $p\\times p$ sample covariance matrix itself is troublesome, which often requires a distributed storage [19]. In this setting, access to principal submatrices faces a communication bottleneck. Our algorithm only requires storing the $n \\times p$ data matrix $X$ and $n \\ll p$. Further, each row of $X$ is independently fed to the algorithm to run in parallel. Methodologically, the method of Huang et al. (2016) depends on the ordering of the variables. Different ordering may result in quite different outcomes. \n\n**REFERENCES**\n\nBogdan, M., Van Den Berg, E., Sabatti, C., Su, W., & Candès, E. J. (2015). SLOPE—adaptive variable selection via convex optimization. *Ann. Appl. Stat.* 9 (3), 1103.\n\nFriedman, J., Hastie, T. and Tibshirani, R. (2010) *Applications of the lasso and grouped lasso to the estimation of sparse graphical models*. Technical Report. Stanford University, Stanford.\n\nSimon, N., Friedman, J., Hastie, T. and Tibshirani, R., 2013. A sparse-group lasso. *J.  Comput. Graph. Stat.*, *22*(2), pp.231-245.\n\nWainwright, M. J. (2019). *High-dimensional statistics: A non-asymptotic viewpoint* (Vol. 48). Cambridge University Press.\n\n### Answer to Q2\n\nThank you for your insightful question. The locally quadratic convergence behavior is indeed one of the key advantages of our proposed second-order method.\n\nTo build intuition, let us contrast this with first-order methods such as gradient descent. In first-order algorithms, each iteration at best reduces the error by a constant multiplicative factor, i.e., convergence is sublinear and becomes linear when the objective is strongly convex. Practically, this means that each iteration improves the solution by approximately one significant digit. As a result, achieving high-precision solutions requires a large number of iterations, especially near the optimum.\n\nIn contrast, second-order methods such as the semismooth Newton method exploit curvature information via (generalized) Hessians or generalized Jacobians. Once sufficiently close to the solution, these methods enjoy quadratic convergence, meaning that the number of accurate digits in the solution approximately doubles for each iteration.\n\nThis dramatic acceleration near the solution is particularly beneficial for solving non-smooth problems like ours, where naive first-order approaches often suffer from slow convergence or stagnation. The value of the first-order method is that it is global — it converges (after many iterations) to the solution wherever it begins, while for second-order methods to quadratically converge, they should start near the solution. So we combined the two in order to globalize the semismooth Newton algorithm so that it converges from everywhere. In our experiments, this theoretical advantage translates into practice: the algorithm consistently converges to high-precision solutions within a few tens of iterations, even for high-dimensional instances.\n\nWe hope this clarifies the intuition behind the observed quadratic convergence."}, {"review_id": "ngXK4oG7rr", "reviewer": "Reviewer_7HYy", "summary": "This paper introduces a semismooth Newton method for partial correlation network estimation, establishing a framework that ensures both local and global convergence while mitigating the high iteration costs of first-order methods. The method reformulates the KKT conditions as a nonlinear equation, employs Bouligand derivatives for global convergence, and addresses linear complementarity problems to optimize iteration direction. Experiments on synthetic and liver cancer genomic data demonstrate the method's quadratic convergence and computational efficiency, achieving a significant speedup over proximal gradient methods.", "questions": "Please see the detailed comments on Weaknesses and Questions.", "limitations": "No. Please see the detailed comments on Weaknesses and Questions.", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n1. The background analysis is comprehensive, focusing on computational bottlenecks in high-dimensional network estimation, such as matrix inversion complexity, with a clearly defined objective.\n2. The innovative framework combines second-order information from semismooth Newton methods with B-derivative globalization, theoretically guaranteeing local quadratic convergence and global stability.\n3. The experimental design is compelling, demonstrating effectiveness and scalability through comparisons with proximal gradient methods, tests across different graph densities, and applications to large-scale genomic data.\n\nWeakness and Questions：\nMajor problems:\n1. Lack of systematic analysis for hyperparameter selection and algorithm switching criteria: The paper does not clarify how regularization parameters (e.g., λ) or line search step sizes are chosen (e.g., via cross-validation). The switching conditions between the two-stage algorithm (local Newton vs. global B-derivative strategy), such as error thresholds or active set rank criteria, lack theoretical justification, potentially affecting convergence stability.\n2. Incomplete ablation studies: Core components (e.g., semismooth Newton iteration, linear complementarity problem solving) are not ablated, precluding quantification of each module’s contribution to convergence speed (e.g., how removing B-derivative globalization affects convergence rate).\n3. Unexplored applicability of quasi-Newton methods: While semismooth Newton methods avoid matrix diagonalization, the paper does not compare quasi-Newton methods (e.g., L-BFGS) for optimizing semismooth functions, missing a cross-validation of computational complexity and convergence properties.\n4. Limited generalizability validation: Experiments are confined to partial correlation network estimation, without verifying the method’s applicability to unstructured optimization problems (e.g., general convex/non-convex functions). Dependencies on problem structures (e.g., row-wise separability) need further clarification.\n\nMinor Problems:\nIntuitive explanations for theoretical derivations (e.g., generalized Jacobian construction) are insufficient, potentially hindering reader comprehension.", "was_revised": true, "final_justification": "Many thanks to the authors for the rebuttal. Most of my concerns were addressed. So I will raise my score to 5. Good luck.", "rebuttal": "## Answer to Major Problems\n### Problem 1.\nWe address each component of your concerns in turn:\n\n(1) Regularization parameter $\\lambda$\n\nFor the LIHC dataset, we used $\\lambda = 0.45$, selected via the extended pseudo-Bayesian Information Criterion (epBIC) as in Section 6 of [19], which employed the first-order ACCORD-FBS algorithm. To assess robustness across sparsity levels, we also tested $\\lambda = 0.2$ (less sparse) and $\\lambda = 0.6$ (more sparse). In both cases, our method converged stably and efficiently, demonstrating reliable performance. We will add how the value $\\lambda=0.45$ was chosen in the revision.\n\n(2) Line search step size\n\nOur algorithm uses a backtracking line search satisfying a sufficient descent condition. In practice, convergence is stable over a wide range of initial step sizes. For visualization purposes (e.g., Figure 3), we selected the step size schedule ($\\rho_k=0.7-0.003k$ and $\\sigma=0.001$) that best illustrates the rapid local convergence of the B-semismooth Newton method. This choice is purely for presentation and does not affect correctness or convergence guarantees.\n\n(3) Switching condition\n\nFirst of all, note that the switching is between the proximal gradient step and the B-semismooth Newton step, not between the local semismooth Netwon step and the B-semismooth Newton step. The former is absorbed into the latter; see L199.\n\nLet $|\\cdot|$ be the size of a set. As noted in the paper, the Newton step becomes applicable when the submatrix $X_{\\bar{\\bar{\\mathcal{A}}}\\_k}$ is of full column rank, i.e., $\\text{rank}(X\\_{\\bar{\\bar{\\mathcal{A}}}\\_k}) = |\\bar{\\bar{\\mathcal{A}}}\\_k|$ . Since $\\text{rank}(X_{\\bar{\\bar{\\mathcal{A}}}_k})\\leq\\min(n,|\\bar{\\bar{\\mathcal{A}}}_k|)$ always holds, a sufficient condition for full column rank is $|\\bar{\\bar{\\mathcal{A}}}_k|\\leq n$.\n\nMotivated by this, we adopt a simplified switching criterion: we switch to the Newton step when $|\\bar{\\bar{\\mathcal{A}}}_k|\\leq n$. Despite its simplicity, this rule has consistently performed well in our experiments without causing convergence issues. While more refined checks (e.g., based on condition numbers) may offer stronger theoretical guarantees, we find our criterion strikes a good balance between computational simplicity and practical robustness. We will clarify this in the revision.\n\n### Problem 2.\n\nTo assess the individual contribution of the proximal step and LCP solving modules, we re-ran the experiment of Section 4.1 for 50 iterations under three ablated configurations: disabling (i) only the proximal step, (ii) only the LCP solving step, and (iii) both. Recall that with both modules enabled, Algorithm 1 terminated within 30 iterations. The result is shown in the following table.\n\nTable : Convergence within 50 iterations under ablation of proximal and LCP updates\n|λ|Full method|No Proximal|No LCP|No Proximal + No LCP|\n|---|---|---|---|---|\n|0.1|✔  (30)|✘|✔   (30)|✘|\n|0.15|✔  (22)|✔  (26)|✔  (22)|✔  (26)|\n|0.2|✔  (14)|✔  (14)|✔  (14)|✔  (14)|\n\nAt $\\lambda = 0.10$, the algorithm failed to converge within 50 iterations whenever the proximal step was disabled, regardless of solving LCP. In contrast, convergence was successfully achieved in settings where the proximal step remained enabled even without LCP solving. In the original run at this $\\lambda$ value, 229 out of 1000 rows were updated via the proximal step during the first iteration alone, indicating substantial reliance on this mechanism.\n\nAt $\\lambda = 0.15$ and $\\lambda = 0.20$, the algorithm converged under all ablation settings. In these cases, at most one proximal update out of 1000 occurred in the first iteration of the original run, followed on by Newton updates thereafter. This suggests that in these sparse settings, the proximal steps may not be necessary. The catch here is, however, that we do not know the sparsity of the solution a priori, and the proximal module is needed as a safeguard.\n\nImportantly, in all configurations where convergence was achieved, the number of iterations to convergence remained comparable to that of the original results. \n\nOn the other hand, the LCP step kicks in only when $|\\omega_i+d_i|=\\lambda$, and it is numerically subtle to declare equality. The above table suggests that Algorithm 1 is largely insensitive to this equality test. The following table summarizes the occurrences of the LCP solving step in the original run, which provides a partial reason:\n|λ|Avg. fraction of LCP per iteration|Max. fraction of LCP per Iteration|Last iteration LCP occurred|Total # iterations|\n|---|---|---|---|---|\n|0.2|0.00005|0.0011|62|159|\n|0.45|0.00004|0.0012|29|69|\n|0.6|0.00012|0.0045|31|63|\n\nThis result demonstrates that the LCP solving step arises only in a negligible fraction of rows. For example, even when the Newton step is applied to as many as 10,000 rows per iteration, the average number of rows updated by LCP solving remains around 1 or fewer, as reflected by the low average/max fractions. Thus its contribution to the search direction is contained.\n\nDespite how rarely LCP occurs, it remains theoretically indispensable to ensure correctness and robustness of the algorithm before we run it on unseen data.\n\nSummary:\n- The proximal step is essential for achieving convergence.\n- Solving the LCP plays a secondary role.\n- The empirical convergence rate is largely unaffected.\n\nWe will incorporate these findings into the revised version to clarify the specific roles and necessity of each module.\n\n### Problem 3.\n\nThe success of (Quasi-)Newton methods such as L-BFGS heavily depends on the assumption that the objective function is twice continuously differentiable everywhere. The ACCORD problem involves non-differentiable points due the $\\ell_1$ penalty; hence, methods that do not account this non-differentiability may not be suitable. This is why we are employing semismooth Newton methods.\n\nThere are a few semismooth quasi-Newton algorithms, but they are highly problem-specific (Mannel and Rund, 2021; Muoi et al., 2013). It would be great if we could develop a quasi-Semismooth Newton algorithm for ACCORD, but we believe that requires nontrivial efforts and is a subject of a separate research paper. \n\nREFERENCES\n\nMannel, F., & Rund, A. (2021). A hybrid semismooth quasi-Newton method for nonsmooth optimal control with PDEs. *Optim. Eng.*, *22*(4), 2087-2125. \n\nMuoi, P. Q., Hào, D. N., Maass, P., & Pidcock, M. (2013). Semismooth Newton and quasi-Newton methods in weighted ℓ1-regularization. *J. Inverse Ill-Posed Probl.*, *21*(5), 665-693.\n\n### Problem 4.\n\nThe focus of the present manuscript is a semismooth Newton method for solving the ACCORD problem [19] and its globalization. Optimization of *general* convex/nonconvex functions is not the scope. We also believe an effective optimization algorithm must exploit the structure of the problem. Otherwise we may have had to stop at Theorem 2.4, leaving much room for practical issues.  As such, row-wise separability contributes to the parallelism in the algorithm and facilitates scalability without compromising theoretical rigor. Note that the semismooth Newton method remains applicable even in the absence of row-wise separability, as long as the optimality map $F$ defined by the KKT conditions is (strongly) semismooth. For example, if the $\\ell_1$ penalty is replaced by the Frobenius norm, then the problem stops being separable but remains semismooth.\n\n## Answer to Minor Problem\n\nDue to space constraints in the original submission, we could not include a detailed derivation of certain technical components, which may have hindered understanding of the algorithm’s mechanics.\n\nHere we provide a more detailed and intuitive explanation of how the generalized Jacobian is constructed, aiming to clarify the mathematical structure and logic, especially for readers less familiar with semismooth Newton methods in nonsmooth optimization.\n\nWe want to solve $F(\\omega^\\*,d^*)=0$ where  $F:\\mathbb{R}^{2p}\\to\\mathbb{R}^{2p}$ is defined as a nonlinear operator \n$\nF(\\omega,d)=\\begin{pmatrix}\\mathbf{S}\\omega+d\\\\\\\\\\omega_1-\\text{prox}\\_{-\\log (\\cdot)}(\\omega\\_1+d\\_1)\\\\\\\\\\omega_{-1}-T_{\\lambda}(\\omega_{-1}+d_{-1})\\end{pmatrix},\\quad\\text{prox}\\_{-\\log (\\cdot)}(x) = \\frac{x + \\sqrt{x^2+4}}{2}.\n$\n\nNote that if $F$ is real-valued and convex, then the generalized Jacobian $\\partial F$ reduces to the usual convex subdifferential. \n\nLet $F_1(\\omega,d)=\\mathbf{S}\\omega+d$ and $F_2(\\omega,d)=\\begin{pmatrix}\\omega_1-\\text{prox}\\_{-\\log (\\cdot)}(\\omega_1+d_1)\\\\\\\\\\omega_{-1}-T_{\\lambda}(\\omega_{-1}+d_{-1})\\end{pmatrix}$. \n\nThen, it is straightforward to see that $\\frac{\\partial F_1}{\\partial\\omega}=\\mathbf{S}$, $\\frac{\\partial F_1}{\\partial d} = \\mathbf{I}$ and $\\frac{\\partial\\text{prox}\\_{-\\log (\\cdot)}(\\omega_1+d_1)}{\\partial\\omega_1} = \\frac{\\partial\\text{prox}_{-\\log (\\cdot)}(\\omega_1+d_1)}{\\partial d_1} =\\frac{1}{2} + \\frac{\\omega_1+d_1}{2\\sqrt{(\\omega_1+d_1)^2+4}}$.\n\nNote that for $i=2,\\cdots,p$,\n$\n\\frac{\\partial T_\\lambda(\\omega_{i}+d_{i})}{\\partial\\omega_{i}}=\\frac{\\partial T_\\lambda(\\omega_{i}+d_{i})}{\\partial d_{i}}=\\begin{cases}0&|\\omega_i+d_i|>\\lambda\\\\\\\\ [0,1]&|\\omega_i + d_i|=\\lambda\\\\\\\\ 1&|\\omega_i + d_i|<\\lambda\\end{cases}.\n$\n\nIf we set $\\frac{\\partial T_\\lambda(\\omega_{i}+d_{i})}{\\partial\\omega_{i}}=0$ when $|\\omega_i+d_i|=\\lambda$, then $\\frac{\\partial F_2}{\\partial\\omega}=\\mathbf{I}-\\mathbf{J}$ where $\\mathbf{J}=\\text{diag}(J_{ii})$ is a diagonal matrix with $J_{11}=\\frac{1}{2}+\\frac{\\omega_1+d_1}{2\\sqrt{(\\omega_1+d_1)^2+4}}$ and  $J_{ii}=1$ if $|\\omega_i+d_i|>\\lambda$, $J_{ii}=0$ if $|\\omega_i+d_i|\\leq\\lambda$ for $i=2,\\dots,p$.\n\nThus, \n$\n\\partial F(\\omega,d)=\\begin{pmatrix}\\frac{\\partial F_1}{\\partial\\omega}&\\frac{\\partial F_1}{\\partial d}\\\\\\\\\\frac{\\partial F_2}{\\partial \\omega}&\\frac{\\partial F_2}{\\partial d}\\end{pmatrix}=\\begin{pmatrix}\\mathbf{S}&\\mathbf{I}\\\\\\\\\\mathbf{I}-\\mathbf{J}&-\\mathbf{J}\\end{pmatrix},\n$\nwhich is equivalent to the example of the generalize Jacobian (10) in the paper."}, {"review_id": "jEOQmB8Ihw", "reviewer": "Reviewer_jUps", "summary": "The paper focuses on the problem of partial correlation network estimation. The classical approach to this is the graphical lasso problem.\nUnfortunately, as the authors note, the graphical lasso struggles to scale to datasets where p>1000, owing to cubic scaling that results from the KKT condition. This is problematic, as in areas like genetics one routinely encounters problems where p>100000, for which graphical lasso would be prohibitively expensive. \n\nRecently, a framework called ACCORD proposes solving an equivalent problem that is much easier to solve for large p. Existing approaches for solving this problem are based on the proximal gradient method, which while possessing scalable iterations, is slow to converge to an acceptable solution. \nIndeed, a problem with p~300000, required more than 200 node hours of compute.\n\nTo address this the problem, the authors develop a bespoke semismooth-Newton method for solving the equivalent problem. Unlike the proximal gradient approach, the proposed method exhibits local quadratic convergence which ensures rapid convergence to a high accuracy solution once the iterates are close enough to the solution. In addition, global convergence from any initialization is ensured. \nDespite being a Newton type method, the authors show that thanks to the structure of the problem, the resulting subproblems reduce to solving small linear systems or linear complementarity problems.\nNumerical experiments validate the theory in the paper and demonstrates superior performance to the current state-of-art proximal gradient algorithm for the problem.", "questions": "- Can you clarify the relation of your use of the semismooth Newton method in relation to references [1,3] in the above? \nThis would solidify the novelty of the approach taken by the paper and would make connections with prior uses of semsimooth Newton for sparse learning problems.", "limitations": "Yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n\n- The presentation in the paper is excellent. The authors thoroughly explain the limitations of existing methodology for solving the partial correlation network estimation problem, and carefully explain each component of their algorithm and how it addresses the problem. Similarly, the proofs of the main results are presented clearly and with sufficient detail.\n\n- The algorithm developed in the paper is efficient and comes with strong theoretical guarantees that address the limitations of existing algorithms.\n\n- Numerical experiments show that the proposed algorithm outperforms the existing baseline on a real world large-scale genetics task, which shows its potential to yield better performance on real world tasks that are of interest to practitioners. \n\n**Weaknesses**\n\n- My main issue with the paper is that it misses some important related work. Semismooth Newton methods have been used for solving sparse learning problems for sometime, see for example [1-3].\nThus, the authors should cite these works and clarify how their approach differs from this prior work.\n\n**Overall**\n\nThe paper makes a good contribution by developing a better solver for an important problem that frequently arises in large-scale genomics. I believe the method will be of interest to members of the community working in this area. \nThe main issue I would like to see the authors address is my question about how their method compares prior semismooth Newton methods for large-scale sparse learning.\n\n**References**\n\n[1] Li, X., Sun, D. and Toh, K.C., 2018. A highly efficient semismooth Newton augmented Lagrangian method for solving Lasso problems. SIAM Journal on Optimization, 28(1), pp.433-458.\n\n[2] Lin, M., Liu, Y.J., Sun, D. and Toh, K.C., 2019. Efficient sparse semismooth Newton methods for the clustered Lasso problem. SIAM Journal on Optimization, 29(3), pp.2026-2052.\n\n[3] Zhang, Y., Zhang, N., Sun, D. and Toh, K.C., 2020. An efficient Hessian based algorithm for solving large-scale sparse group Lasso problems. Mathematical Programming, 179, pp.223-263.", "was_revised": true, "final_justification": "The issues regarding the novelty of the method and contributions relative to prior methods utilizing Semismooth Newton for sparse learning was addressed. As I already thought the paper was worth accepting, I'll be maintaining my initial recommendation.", "rebuttal": "### Answer to Weakness 1 and Question 1 \n\nThe line of work you referred to, such as those solving Lasso, Clustered Lasso, and Sparse Group Lasso problems ([1]–[3]), falls into the semismooth Newton-based augmented Lagrangian method (SSNAL) framework. As the name suggests, this framework solves the problem via the inexact augmented Lagrangian method (ALM). ALM incurs an optimization subproblem, and SSNAL solves this subproblem via semismooth Newton. While the subproblem is solved at a locally superlinear (quadratic if it is strongly semismooth) rate, the convergence rate of the the outer ALM iteration is locally linear, even if the subproblem can be solved at a locally quadratic rate. Global convergence is also established.\n\nIn contrast, our algorithm formulates the main problem (KKT equation) as a strongly semismooth equation and apply semismooth Newton. Globalization is achieved via damping instead of ALM. So our algorithmic design is fundamentally different from SSNAL: there is no subproblem to solve, and the convergence is locally quadratic. (If the KKT equation was merely semi smooth, then the local convergence is superlinear.) \n\nAnother distinction is that SSNAL typically constructs augmented Lagrangian via the dual of the original (primal) problem, while we focus on the primal. It is worth pointing out that [2] applies the SSNAL framework directly to the primal in addition to the dual, but its design and analysis focus on low-dimensional settings where the number of observations significantly exceeds the number of features (i.e.,$n \\gg p$) As such, it is not readily applicable to high-dimensional settings, which are the primary focus of our work.\n\nWe appreciate the opportunity to clarify these differences and will explicitly cite and discuss the relevant SSNAL literature in the revised version to better position our contribution in the context of prior work."}]

Author rebuttal:

Meta-review: The paper considers partial correlation network estimation. The classical approach, graphical lasso, scales cubically and hence struggles to handled datasets where p>1000. Instead, the recently proposed ACCORD framework solves an equivalent problem that is much easier to solve for large p. Prior work based on ACCORD solve the equivalent problem using proximal gradient, which enjoys scalable iterations but converge too slowly. This paper, therefore, develops a bespoke semismooth-Newton method for solving the equivalent problem. The method enjoys quadratic local convergence and can be implemented to be scalable. Experiments on synthetic and large liver-cancer genomics show large speedups from far fewer iterations compared to proximal-gradient baselines. The method also reduces hardware demands compared with prior MPI/cluster settings.

Concerns mostly focus on scope and presentation: clarifying assumptions (e.g., rank conditions), documenting hyperparameter choices and switch criteria, broadening comparisons (e.g., quasi-Newton), and adding robustness studies (non-Gaussian data). The rebuttal added sensitivity analyses, ablations, and non-Gaussian results, and clarified how ACCORD addresses glasso bottlenecks at extreme dimension.
Decision: Accept (poster)