Global Optimality in Bivariate Gradient-based DAG Learning
Authors: Chang Deng
Kevin Bello
Pradeep Kumar Ravikumar
Bryon Aragam
Venue: NeurIPS 2023 poster (2023)
Keywords: global optimization
nonconvex optimization
graphical models
directed acyclic graphs
structure learning
TL;DR: We prove that a simple path-following optimization scheme globally converges to the global minimum of the population loss in the bivariate setting of DAG Learning.

Abstract: Recently, a new class of non-convex optimization problems motivated by the statistical problem of learning an acyclic directed graphical model from data has attracted significant interest. While existing work uses standard first-order optimization schemes to solve this problem, proving the global optimality of such approaches has proven elusive. The difficulty lies in the fact that unlike other non-convex problems in the literature, this problem is not "benign", and possesses multiple spurious solutions that standard approaches can easily get trapped in. In this paper, we prove that a simple path-following optimization scheme globally converges to the global minimum of the population loss in the bivariate setting.

Global Optimality in Bivariate Gradient-based DAG Learning
Bryon Aragam *Booth School of Business, University of Chicago, Chicago, IL 60637 Machine Learning Department, Carnegie Mellon University, Pittsburgh, PA 15213
Chang Dengi* Kevin Bellofd
Pradeep Ravikumar!
Abstract
Recently, new class of non-convex optimization problems motivated by the statistical problem of learning an acyclic directed graphical model from data has attracted significant interest. While existing work uses standard first-order optimization schemes to solve this problem, proving the global optimality of such approaches has proven elusive. The difficulty lies in the fact that unlike other non-convex problems the literature, this problem is not "benign", and possesses multiple spurious solutions that standard approaches can easily get trapped in. In this paper, we prove that a simple path-following optimization scheme globally converges to the global minimum of the population loss in the bivariate setting.
1 Introduction
Over the past decade, non-convex optimization has become major topic of research within the machine learning community, in part due to the successes of training large-scale models with simple first-order methods such as gradient descent- with their stochastic and accelerated variants- in spite of the non-convexity of the loss function. A large part of this research has focused on characterizing which problems have benign loss landscapes that are amenable to the use of gradient- based methods, i.c., there are no spurious local minima, or they can be easily avoided. By now, several theoretical results have shown this property for different non-convex problems such as: learning two hidden unit ReLU network [4]] learning (deep) over- pee-parameteriiedd quadratic neural networks 4 27], low-rank matrix recovery 19 3), learning a two-layer ReLU network with
single non-overlapping convolutional filter [6), semidefinite matrix completion 4 20 learning
neural networks for binary classification with the addition of single special neuron and learning deep networks with independent ReL.U activations [26 11, to name a few.
Recently, new class of non-convex optimization problems due Zheng al. 51 have emerged in the context of learning the underlying structure of a structural equation model (SEM) or Bayesian network. This underlying structure typically represented by a directed acyclic graph (DAG), which makes the learning task highly complex due to its combinatorial nature. In general, learning DAGs is well-known be NP-complete g 10 The key innovation in Zheng al. [ST] was the introduction of differentiable function h, whose level set at zero exactly characterizes DAGs. Thus, replacing the challenges of combinatorial optimization by those of non-convex optimization. Mathematically, this class of non-convex problems take the following general form:
min f(O) subject to h(W(()) = 0, e
(1)
where θ R represents the model parameters, f R - R is (possibly non-convex) smooth loss function (sometimes called a score function) that measures the fitness of θ, and Rdxd - [0, 20)
changleng@@hicgggoocc...e
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

is smooth non-convex function that takes the value of zero if and only if the induced weighted adjacency matrix of d nodes, W(), corresponds to a DAG.
Given the smoothness of f and h. problem ①1 can be solved using off-the-shelf nonlinear solvers, which has driven series of remarkable developments in structure learning for DAGG. Multiple empirical studies have demonstrated that global or near-global minimizers for ① can often be found in
variety of settings, such as linear models with Gaussian and non-Gaussian noises [c.g.. 0g....111 34 1] and
non-linear models, represented by neural networks, with additive Gaussian noises [e.g.. 29 52 49
The empirical success for learning DAGs via (II which started with the NOTEARS method of Zheng
et al. [SI], bears a resemblance to the advancements in deep learning, where breakthroughs like AlexNet significantly boosted the field's recognition, even though there were notable successes before it.
Importantly, the reader should note that the majority of applications in MI consist of solving single unconstrained non-convex problem. In contrast, the class of problems (1) contains a non-convex constraint. Thus, researchers have considered some type of penalty method such as the augmented Lagrangian [51 (2), quadratic penalty (35), and a log-barrier (Л. In all cases, the penalty approach consists in solving a sequence of unconstrained non-convex problems, where the constraint enforced progressively [see e.g. 2 for background]. In this work, we will consider the following form of penalty:
min 9uz (+) = (0) + h(W())).
was shown by Bello et al. 1 that due to the invexity property of h solutions to (2) will converge to DAG as μl - 0. However, no guarantees on local/global optimality were given.
With the above considerations in hand, one is inevitably led to ask the following questions:
(i) Are the loss landscapes 9u. (0) benign different ps?
(ii) Is there a (tractable) solution path [0x] that converges to a global minimum of ([)
Due to the NP-completeness of learning DAGs, one would expect the answer to (i) to be negative in its most general form. Moreover, it is known from the classical theory of constrained optimization [e.g. 2 that if we can exactly and globally optimize ( for each (4g, then the answer to (īi) is affirmative. This is not practical algorithm, however, since the problem ( is nonconvex. Thus we seek solution path that can be tractably computed in practice, e.g. by gradient descent.
In this work, we focus on perhaps the simplest setting where interesting phenomena take place. That is, a linear SEM with two nodes (i.e., d = 2), f is the population least squared loss (i.e., f is convex), and θ₂ is defined via gradient flow with warm starts. More specifically, we consider the case where O is obtained by following the gradient flow of 9ps. with initial condition Ox-1-
Under this setting, to answer (i), it is easy to see that for a large enough /4s, the convex function dominates and
expect a benign landscape, i.e., a (almost) convex landscape. Similarly, when μı approaches zero, the invexity of f kicks in and we could expect that all stationary points are (near) global minimizers | That is, the extremes μ4g- - Do and μ4 + O, the landscapes seem well-behaved, and the reader might wonder if it follows that for any μr E [0, oo) the landscape is well-behaved. We answer the latter in the negative and show that there always exists a 7 > 0 where the landscape of 9ps. is non- benibn for any 14 < T, namely, there exist three stationary points: i) saddle point, ii) spurious local minimum, and iii) The global minimum. In addition, each of these stationary points have wide basins of attractions, thus making the initialization of the gradient flow for 9M, crucial. Finally, we answer (ii) in the affirmative and provide an explicit scheduling for 14 that guarantees the asymptotic convergence of θ₂ to the global minimum of (I] Moreover, we show that this scheduling cannot be arbitrary as there exists a sequence of (14) that leads [O₂} to a spurious local minimum.
Overall, we establish the first set of results that study the optimization landscape and global optimality for the class of problems ( We believe that this comprehensive analysis in the bivariate case provides a valuable starting point for future research in more complex settings.
An invex function is any function where all its stationary points are global minima. It worth noting that the composite objective in (2) is not necessarily invex, even when f convex. 'This transition or path. from an optimizer of simple function to an optimizer of a function that closely resembles the original constrained formulation, is also known homotopy.

---

65
Sodtle Foint
Global
64
10
60
(a) Contour plot
(b) Stationary points
Figure I: Visualizing the nonconvex landscape. (a) A contour plot of 9μ. for a = 0.5 and μ = 0.005 (see Section 2 for definitions). We only show a section of the landscape for better visualization. The solid lines represent the contours, while the dashed lines represent the vector field - -yg- (b) Stationary points of gμ. rly: (1)
and r(x; p)) 0 (see Section 4 for definitions).
Remark We emphasize that solving ① in the bivariate case is not an inherently difficult problem. Indeed, when there are only two nodes, there only two DAGs to distinguish and one can simply fit f under the only two possible DAGs, and select the model with the lowest value for f. However, evaluating f for each possible DAG structure clearly cannot scale beyond 10 or 20 nodes, and is not a standard algorithm for solving (T). Instead, here our focus is on studying how (T) is actually being solved in practice namely, by solving unconstrained non-convex problems in the form of (2). Previous work suggests that such gradient- based approaches indeed scale well to hundreds and even thousands of nodes le.g.
1.1 Our Contributions
More specifically, we make the following contributions:
1. We present homotopy-based optimization scheme (Algorithm 2 to find global minimizers of the program (T) by iteratively decreasing the penalty coefficentt according given schedule. Gradient flow is used to find the stationary points of (2) at each step, starting from the previous solution.
2. We prove that Algorithm 2 converges globally (i.e. regardless of initialization for W) to the global minimum (Theorem D.
3. We show that the non-convex program | is indeed non-benign, and naive implementation of black-box solvers are likely to get trapped in bad local minimum. See Figure] I (a).
4. Experimental results verify our theory, consistently recovering the global minimum of (L). regardless of initialization or initial penalty value. We show that our algorithm converges to the global minimum while naïve approaches can get stuck.
The analysis consists of three main parts: First, explicitly characterize the trajectory of the stationary points of (2). Second, we classify the number and type of all stationary points (Lemma and use this isolate the desired global minimum. Finally, we apply Lyapunov analysis identify the basin of attraction for each stationary point, which suggests a schedule for the penalty coefficient that ensures that the gradient flow is initialized within that basin the previous solution.
1.2 Related Work
The class of problems ① falls under the umbrella of score-based methods, where given a score function f. the goal is to identify the DAG structure with the lowest score possible g 22.. We shall note that learning DAGs is very popular structure model in a wide range domains such as biology 40), genetics [50], and causal inference 44 39 to name few.
3

---

Score-based methods that consider the combinatorial constraint. Given the ample set of score- based methods the literature, we briefly mention some classical works that attempt to optimize by considering the combinatorial DAG constraint. In particular, we have approximate algorithms such as the greedy search method of Chickering et al. [10), order search methods 45
41 [8], the
LP-relaxation method of Jaakkola et al. [24), and the dynamic programming approach of Loh and
Bühlmann [3I]. There are also exact methods such as GOBNILP [12] and Bene [42], however, these
algorithms only scale up to R 30 nodes.
Score-based methods that consider the continuous non-convex constraint h. The following works are the closest to ours since they attempt to solve problem in the form of (T. Most of these developments either consider optimizing different score functions such as ordinary least squares [51 (2), the log-likelihood [29] 34), the evidence lower bound [49], regret function [53]; or consider different differentiable characterizations of acyclicity h [9 However, none of the aforementioned works provide any type of optimality guarantee. Few studies have examined the optimization intricacies of problem (T). Wei et al. [47] investigated the optimality issues and provided local optimality guarantees under the assumption of convexity in the score f and linear models. On the other hand, Ng et al. (35) analyzed the convergence to (local) DAGs of generic methods for solving nonlinear constrained problems, such as the augmented Lagrangian and quadratic penalty methods. In contrast to both, our work is the first to study global optimality and the loss landscapes of actual methods used in practice for solving
Bivariate causal discovery. Even though in a two-node model the discrete DAG constraint does not pose major challenge, the bivariate setting has been subject major research in the area of and references therein.
causal discovery. See for instance 36 16 32
Penalty and homotopy methods. There exist classical global optimality guarantees for the penalty method if f and h were convex functions, see for instance 12. 3 37.. However, to our knowledge, there are no global optimality guarantees for general classes of non-convex constrained problems, let alone for the specific type of non-convex functions h considered in this work. On the other hand, homotopy methods (also referred as continuation or embedding methods) are in many cases capable of finding better solutions than standard first-order methods for non-convex problems, albeit they typically do not come with global optimality guarantees either. When homotopy methods come with global optimality guarantees, they are commonly computationally more intensive as it involves discarding solutions, thus, closely resembling simulated annealing methods, see for instance [15] Authors in [1 characterize family of non-convex functions where a homotopy algorithm provably converges to a global optimum. However, the conditions for such family of non-convex functions are difficult to verify and are very restrictive; moreover, their homotopy algorithm involves Gaussian smoothing, making it also computationally more intensive than the procedure we study here. Other examples of homotopy methods in machine learning include 7 18 global optimality guarantees are given.
46
7
in all these cases, no 2
2 Preliminaries
The objective we consider can be easily written down as follows:
S(W) = Ex
[|X - W X
2
where X E I² is random vector and W E 22x2 Although not strictly necessary for the developments that follow, we begin by introducing the necessary background on linear SEM that leads to this objective and the resulting optimization problem of interest.
The bivariate model. Let = (2), E I² denote the random variables in the model, and let N (N1, N2) E R² denote vector of independent errors. Then linear SEM over X is defined as X = W. X + N, where W. E R2x2 is a weighted adjacency matrix encoding the coefficients in the linear model. In order to represent valid Bayesian network for X [see e.g. 39 44 for details), the matrix W. must be acyelic: More formally, the weighted graph induced by the adjacency matrix W. must be a DAG. This (non-convex) acyclicity constraint represents the major computational hurdle that must overcome in practice (cf. Remark I

---

The goal is to recover the matrix W. from the random vector X. Since W. is acyclic, we can assume the diagonal of W. zero (i.e. no self-loops). Thus, under the bivariate linear model, it then suffices to consider two parameters x and y that define the matrix of parameters
W = ((x,v) = 0
(4)
For notational simplicity, we will use f(W) and faz y) interchangeably, similarly for h(W) and h(2, y). Without loss of generality, we write the underlying parameter as
W. = 0 a
(5)
which implies
[X 1x2=ax-- N2
X = W
In general, we only require N; to have finite mean and variance, hence we do not assume Gaussianity, We assume that [ar[N]] Var [N2]]. and for simplicity, we consider E[[]] = 0 and Cov[N] = I, where / denotes the identity matrix. Finally, in the sequel we assume w.l.o.g. that a > 0.
The population least squares. In this work, consider the population squared loss defined by 3 If we equivalently write / in terms of x and y, then we have: f(W)- = (a)--yy)++++(--0))) In fact, the population loss can be substituted with empirical loss. In such a case, our algorithm can still attain the global minimum, WG. of problem (6). However, the output WG will serve as an empirical estimation of W., An in-depth discussion on this topic can be found in Appendix B
The non-convex function h. We use the continuous acyclicity characterization of Yu et al. 149 i.c. h(W)= Tr((I + W o W)d) d, where denotes the Hadamard product. Then, for the bivariate case, we have h(W) this version of h. that is, we can use the same techniques used throughout this work for other existing formulations of h, such as the trace of the matrix exponential (51), and the log-det formulation (I]. Nonetheless, here we consider that the polynomial formulation of Yu et al. [4 is more amenable for the analysis.
y² /2. We note that the analysis presented in this work is not tailored to
Our restriction to the bivariate case highlights the simplest setting in which this problem exhibits nontrivial behaviour. Extending our analysis higher dimensions remains challenging future direction, however, we emphasize that even in two-dimensions this problem is nontrivial. Our approach is similar to that taken in other parts of the literature that started with simple cases (e.g. single-neurd neuron models in deep learning).
Remark 2.
Remark 3. It is worth noting that our choice of the population least squares is not arbitrary. Indeed, for linear models with identity error covariance, such as the model considered in this work, it is known that the global minimizer of the population squared loss is unique and corresponds to the underlying matrix W.. See Theorem 7 in (31%.
Gluing all the pieces together, we arrive to the following version of (1) for the bivariate case: 1 min (zx y) = _ ((1 az)² + x² + (a x²²² +1) subject to h(r,,) v) = x²y² z.y
0.
(6)
2
Moreover, for any 2 0, we have the corresponding version of 2 expressed as:
y2y² min g, = uz v) + M(z, v) 2/2 ((1 - a1)= + + (a x)² + 11) 2,8
2
To conclude this section, we present a visualization of the landscape of gp (xx, y) in Figure (a), for a 0.5 and 0.005. We can clearly observe the non-benign landscape of 9μ., i.e., there exists a spurious local minimum, a saddle point, and the global minimum. In particular, we can see that the basin of attraction of the spurious local minimum is comparable to that of the global minimum, which is problematic for a local algorithm such as the gradient flow (or gradient descent) as it can easily get trapped in a local minimum if initialized in the wrong basin.
'Following the notation in 1 for the bivariate model we simply have θ = (x, y) and W(O)= = ()a). Although is used to represent edge weights and not data in 3 this distinction should not lead to confusion
5

---

Algorithm 1: GradientFlow(S,.2. 20)
set =0) 2g
2: _ =t) (z(())
3: return lim 2t)
Algorithm 2: Homotopy algorithm for solving (.).
a² 4(p2+1j3 4
Input: Initial Wo = W (Iv, Mo), Mo E
Output: { 1,0 k==
W na t GradientFlow( Wo)
2 for k = 1, ......d.
Let μı= = (2/a)a/3
2/3 4/3
pk-1
W < GradiontFlov(quu W w
5 end
3 A Homotopy-Based Approach and Its Convergence to the Global Optimum
To fix notation, let us write W == W == (
Vug 0 ). and let Wc denote the global minimizer (
In this section, we present our main result, which provides conditions under which solving series of unconstrained problems 7 with first-order methods will converge to the global optimum We of (6). in spite of facing non- benign landscapes. Recall that from Remark 3 we have that W6 = -00) Since we use gradient flow path to connect WHy, and Wps we specify this path Procedure [ for clarity. Although the theory here assumes continuous-time gradient flow with /. - co, see Section for an iteration complexity analysis for (discrete-time) gradient descent, which a straightforward consequence of the continuous-tin theory.
In Algorithm 2 we provide an explicit regime of initialization for the homotopy parameter 140 and a specific scheduling for μt such that the solution path found by Algorithm 2 will converge to the global optimum of (6). This is formally stated in Theorem whose proof is given in Section
Theorem 1. For any initialization Wo and a E TR, the solution path provided in Algorithm converges the global optimum of 66. ie.,
lim W WG. k+000
few observations regarding Algorithm 2 Observe that when the underlying model parameter >>>>>>>>>> O, the regime of initialization for μ4 is wider; on the other hand, if is closer to zero then the interval for No is much narrower. As a concrete example, if a 2 then suffices to have Mo E [0.008, 1); whereas if 0.1 then the regime is about (e E [0.0089, 0.01). This matches the intuition that for a "stronger" value of a it should be easier to detect the right direction of the underlying model. Second, although in Line 3 we set μv in specific manner, actually suffices to have
Mak [ Mk (
a./ (4px-1)4a)) -.--1)
We simply chose particular expression from this interval for clarity of presentation; see the proof in Section 5 for details.
As presented, Algorithm is of theoretical nature in the sense that the initialization for μg and the decay rate for μl in Line 3 depend on the underlying parameter a, which in practice is unknown. In Algorithm 3 we present a modification that is independent of and W.. By assuming instead lower bound on a, which is standard assumption in the literature, we can prove that Algorithm also converges to the global minimum:
Corollary 1. Initialize ps z !fa / 55227 then for any initialization Wo. Algorithm oulputs the global optimal solution to (6), i.e.
lim Wie k
WG.
For more details on this modification, see Appendix

---

Algorithm 3: Practical (i.e. independent of a and W,) homotopy algorithm for solving (I).
Input: Initial Wo W (To.30)
Output: (W 1
H4 f 1/27
W nu f- GradientFlow(gion Wo)
3 for k = 1,2,........ .22......
Let jık (2/ 5pm)
2/3 4/3
Mk-1
Wn, & GradientFlore W Pz .
end
1164
02
635
81
635
62
825
03
626
21
045
625
03
033
03
033
05
035
04
04
820
31
035
63
826
83
035
94
66
(8)
(b)
(c)
7
Figure 2: The behavior of r(y μ1) for differrntt μ. Here, for there exists single solution to
r(yi jJ) = 0, which implies there is one stationary point in Equation (7). When μ = T, two solutions found for r(y: (4) O, suggesting that there
stationary points in Equation 7.. Conversely, when μ < T, we observe three solutions for r(y: μ)) = O, indicating that there are three stationary points in Equation 7 local optimum, saddle point, and global optimum.
4 A Detailed Analysis of the Evolution of the Stationary Points
The homotopy approach in Algorithmm relies heavily on how the stationary points of 7 behave with respect /4k. In this section, dive deep into the properties of these critical points.
By analyzing the first-order conditions for Ju. we first narrow our attention to the region fo < S a, < S - }. By solving the resulting equations, we obtain an equation that only involves the variable
r(y: p1)
+ +)2 (a² 1).
(8)
y
(y²
Likewise, we can find an equation only involving the variable æ:
aa?
((:; (4) (a(a² + 1) x²)²
(9)
To understand the behavior of the stationary points of 9 (W), t( ()) in the range x E [0, a} and the properties r μt) in the interval y E [0, c
examine the characteristics of
In Figures 2 and 3 we show the behavior of r(y; (4) and ((;; μ)) for Theorems 5
appendix establish the existence of a > 0 with the following useful property:
Corollary 2. There exists μ < T such that the equation Vgn(W) has three different solutions, denoted as W W** W*** Then,
0
[
0 0
lim W
lim
W P
lim W #
0
0
0
Note that the interesting regime takes place when μ < v. Then, we characterize the stationary points either local minima or saddle points:
7

---

14.
1.
(a)
(b)
(c)
Figure 3: The behavior of t(z; p)) for different μ. Here, for μ > T, there exists single solution to
t(z; p1) = 0. which implies there is one stationary point in Equation 77). When μ т, two solutions found for t(z; (4) O, suggesting that there are two stationary points in Equation 7.. Conversely,
when μ < т, we observe three solutions for t(x; μ)) = O, indicating that there are three stationary points in Equation 7 local optimum, saddle point, and global optimum.
Lemma 1. Let < T, then 9p (W) has two local minima at W W and saddle point at W
With the above results, it has been established that W converges the global minimum WG as + O. In the following section for the proof of Theorem [] we perform a thorough analysis on how track Wi and avoid the local minimum W by carefully designing the scheduling for μig
5 Convergence Analysis: From continuous to discrete
We now discuss the iteration complexity of our method when gradient descent used in place of gradient flow. We begin with some preliminaries regarding the continuous-time analysis.
5.1 Continuous case: Gradient flow
The key to ensuring the convergence of gradient flow to W is to accurately identify the basin of attraction of W". The following lemma provides region that lies within such basin of attraction.
Lemma 2. Define By = {(x,y) 1 < < 0, < y }. Run Algorithm / with imput f = (r, y), 20 = (z()), (() where (I(0), v(0) E B then Vt 2 O. we have that v(t) E B, and lim (I)) = (1h., 9.)
In Figure (b), the lower-right rectangle corresponds to B,. Lemma 2 limpliss that the gradient flow with any initialization inside B will converge to W Me+z at last. Then, by utilizing the previous solution W ne the initial point, long as it lies within region B Bu the gradient flow can converge to W pR+1" thereby achieving the goal of tracking W A4:1- Following the scheduling for /4 prescribed
in Algorithm//provddees provides a sufficient condition to ensure that will happen.
The following lemma, with proof in the appendix, is used for the Proof of Theorem It provides lower bound for yj," and upper bound for y
Lemma 3. If < T. then Uy > VTa and (3 10/33Va212 (44))/33) > y/-
Proof of Theorem 1 Consider that we
iteration & + of Algorithm 2 then 14411 < Pg- If
[4 > and p4+1 > then there is only one stationary point for 9us (x,, y) and gus (x, y) thus, I will converge such stationary point. Hence, let us assume 1k+1 S т. From Theorem 6 in the appendix, we know that x0 22ats
< Fmgg Then, the following relations hold:
1/3
143
(at/3 - V22/s - (44x
Muues V 2 2
Vis
4a
2
Here (1) and (3) are due Lemma implics that (I yi is in the region {(x) y) I < a, S y y J. By Lemma 2 the
3 and (2) follows from / x for S S 1. Then

---

I procedure will converge to (2) i Mig ). Finally, from Theorems 5 and 6 if limg->o P4: = 0,
Then lime- Sp = a, 1om Vya 0, thus, converging to the global optimum, i.c..
lim W WG. k+000
5.2 Discrete case: Gradient Descent
In Algorithms 2 and 4 gradient flow is employed to locate the next stationary points, which is not practically feasible. A viable alternative is to execute Algorithm replacing the gradient flow with gradient descent. Now, at every iteration k, Algorithm uses gradient descent to output W Me c stationary point 9u4. initialized at W lift 1 and a step size of mp,- = 1/(2e((a + + 3a²). The tolerance parameter ég: can significantly influence the behavior of the algorithm and must be controlled for different iterations. A convergence guarantee is established via a simplified theorem presented here. A more formal version of the theorem and comprehensive description of the algorithm (i.e., Algorithm 6 can be found in Appendix C
Theorem 2 (Informal). For any Edist O, set μo satisfy mild condition, and use updating rule EN: 222))2/3((ate/// (atcc////)/ min Baap P4 3/2 }. pk+11 = (2pz)2/3 (a (a-ce/1u)477) and let K = K K(H,, a, Edist) E O In Po
afdiet
Then, for any initialization Wo, following the updated procedure above for k 0,..., K, we have: I|W me.ex
uk.ex - WGlle S Edist
that is, W pn,en is Edist close in Frobenius norm global optimum WG. Moreover, the total number of gradient descent steps is upper bounded by O (muaa² + a² + No) na + =dme
Experiments
We conducted experiments to verify that Algorithms 2 and both converge to the global minimum of
77) Our purpose to illustrate two main points: First, we compare our updating scheme as given in Line 3 of Algorithm 2 against faster-decreasing updating scheme for μıg- In Figure wee illustrate how naive faster decrease of μ can lead spurious local minimum. Second, in Fig we show that regardless of the initialization, Algorithms and [andaalaays Jlways return the global minimum. In the supplementary material, we provide additional experiments where the gradient flow is replaced with gradient descent. For more details, please refer to Appendix F
0.0
3278
231
Figure 4: Trajectory of the gradient flow path for two different update rules for μ with same initialization and Noo Here, "good scheduling" uses Line of Algorithm [ while "bad scheduling" uses a faster decreasing scheme for μ4 which leads the path to a spurious local minimum.
7 Acknowledgments and Disclosure of Funding
K. B. was supported by NSF under Grant # 2127309 to the Computing Research Association for the CIFellows 2021 Project. B.A. was supported by NSF IIS- ISS-195330,, NIH RO1GM1 and the Robert H. Topel Faculty Research Fund at the University of Chicago Booth School of Business. This work was done in part while B.A. was visiting the Simons Institute for the Theory of Computing. P.R. was supported by ONR via N0001418 110011111661, and NSF via IIS- 16, IIS- 195553,,,,, IIS-221
We are also grateful for the support of the University Chicago Research Computing Center for assistance with the calculations carried out in this work.

---

125
033
022
210
115
41
4.5
(a) Random Initialization
(b) Random Initialization
(c) Random Initialization
Figure 5: Trajectory of the gradient flow path with the different initializations. We observe that under a proper scheduling for pas they all converge to the global minimum.
References
[1] Bello, K., Aragam, B. and Ravikumar, P. [2022], DAGMA: Learning dags via M-matrices and a log-determinant acyclicity characterization, 'Advances in Neural Information Processing Systems'.
[2] Bertsekas, D. P. [1997], 'Nonlincar programming' Journal of the Operational Research Society 48(3), 334-334.
[3] Bhojanapalli, S., Neyshabur, B. and Srebro, N. [2016), 'Global optimality of local search for low rank matrix recovery', Advances in Neural Information Processing Systems 29.
[4] Boumal, N., Voroninski, V. and Bandeira, A. [2016]. 'The non-convex burer-monteiro approach works smooth semidefinite programs' Advances in Neural Information Processing Systems 29.
[5] Boyd, S., Boyd, S. P. and Vandenberghe, L. [2004], Convex optimization, Cambridge university press.
[6] Brutzkus, A. and Globerson, A. [2017], Globally optimal gradient descent for convnet with gaussian inputs, in "International conference on machine learning', PMLR, pp. 605-614.
[7] Chen, W. Drton, M. and Wang, Y. S. [2019], 'On causal discovery with equal-variance assumption', Biometrika 106(4), 973-980.
[8] Chickering, D. M. [1996]. Learning Bayesian networks is NP-complete, in 'Learning from data', Springer, pp. 121-130.
[9] Chickering, D. M. [2003], 'Optimal structure identification with greedy search', JMLR 3, 507554.
[10] Chickering, D. M., Heckerman, D. and Meek, C. [2004), 'Large-sample learning of Bayesian networks is NP-hard', Journal of Machine Learning Research 5, 1287-1330.
[11] Choromanska, A., Henaff, M., Mathieu, M., Arous, G. B. and LeCun, Y. [2015], The loss surfaces of multilayer networks, in Artificial intelligence and statistics', PMLR, pp. 192-204.
[12] Cussens, J. [2012], 'Bayesian network learning with cutting planes", arXiv preprint arXiv: 1202.3713.
[13] De Sa, C., Re, C. and Olukotun, K. [2015], Global convergence of stochastic gradient descent for some non-convex matrix problems, in 'International conference on machine learning', PMI., pp. 2332-2341.
[14] Dontchev, L.., Rockafellar, T. and Rockafellar, T. [2009], Implicit functions and solution mappings: A view from variational analysis, Vol. 11, Springer.
10

---

[15] Dunlavy, D. M. and O'Leary, P. [2005]. Homotopy optimization methods for global optimiza- tion, Technical report, Sandia National Laboratories (SNL), Albuquerque, NM, and Livermore, CA.
[16] Duong, and Nguyen. T. [2022], Bivariate causal discovery via conditional divergence, in "Conference on Causal Learning and Reasoning'.
[17] Gargiani, M., Zanelli, A., Tran-Dinh, Q., Diehl, M. and Hutter, F. [2020], 'Convergence analysis of homotopy-sgd for non-convex optimization', aaXXiv preprint arXiv:201 1.10998.
[18] Garrigues, P. and Ghaoui, L. [2008], 'An homotopy algorithm for the lasso with online observa- tions", Advances in neural information processing systems 21.
[19] Ge, R., Jin, C. and Zheng, Y. [2017], No spurious local minima in nonconvex low rank problems: A unified geometric analysis, in 'International Conference on Machine Learning', PMLR, pp. 1233- 1242.
[20] Ge, R., Lee, D. and Ma, T. [2016], 'Matrix completion has no spurious local minimum', Advances in neural information processing systems 29.
[21] Hazan, E., Levy, K. Y. and Shalev-Shwartz, S. [2016], graduated optimization for stochastic non-convex problems, in 'International conference on machine learning', PMLR, pp. 1833-1841.
[22] Heckerman, D., Geiger, D. and Chickering, D. M. [1995], 'Learning bayesian networks: The combination of knowledge and statistical data', Machine learning 20(3), 197-243.
[23] Iwakiri, H., Wang, Y., Ito, S. and Takeda, A. [2022], 'Single loop gaussian homotopy method for non-convex optimization', arXiv preprint 2035:7200.
[24] Jaakkola, T., Sontag, D., Globerson, A. and Meila, M. [2010], Learning bayesian network structure using Ip relaxations, 'Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics', Vol. 9 of Proceedings of Machine Learning Research, pp. 358-365.
[25] Jiao, R., Lin, N., Hu, Z., Bennett, D. A., Jin, L. and Xiong, M. [2018]. 'Bivariate causal discovery and its applications gene expression and imaging data analysis', Frontiers Genetics 9, 347.
[26] Kawaguchi, K. [2016], 'Deep learning without poor local minima', Advances in neural informa- tion processing systems 29.
Advances
[27] Kazemipour, A., Larsen, B. W. and Druckmann, S. [2019], 'Avoiding spurious local minima in deep quadratic networks', anXiv preprint 2tt:2000000009.
[28] Khalil, H. K. [2002]. 'Nonlinear systems third edition', Patience Hall 115.
[29] Lachapelle, S., Brouillard, P., Deleu, T. and Lacoste-Julien, S. [2020], Gradient-bart based neural dag learning, in 'International Conference on Learning Representations'
[30] Liang, S., Sun, R., Lee, J. D. and Srikant, R. [2018], 'Adding one neuron can eliminate all bad local minima', Advances in Neural Information Processing Systems 31.
[31] Loh, P.L. and Bühhlann, P. [2014], 'High-dimensionall learning of linear causal networks via inverse covariance estimation', The Journal of Machine Learning Research 15(1), 3065-3105.
[32] Mooij, J. M., Peters, J., Janzing, D., Zscheischler, J. and Schölkopf, B. [2014], 'Distinguishing cause from effect using observational data: methods and benchmarks', Arxiv.
[33] Nesterov, Y. et al. [2018], Lectures convex optimization, Vol. 137, Springer.
[34] Ng, I., Ghassami, A. and Zhang, K. [2020], On the role of sparsity and dag constraints for learning linear DAGs, in 'Avvances in Neural Information Processing Systems'.

---

[35] Ng, I., Lachapelle, S., Ke, N. R., Lacoste-Julient S. and Zhang, K. [2022], On the convergence of continuous constrained optimization for structure learning, in 'International Conference on Artificial Intelligence and Statistics', PMLR, pp. 8176-8198.
[36] Ni, Y. [2022], 'Bivariate causal discovery for categorical data via classification with optimal label permutation', Arxiv.
[37] Nocedal, J. and Wright, S. [1999), Numerical optimization, Springer.
[38] Park, Y. W. and Klabjan, D. [2017], 'Bayesian network learning via topological order', The Journal of Machine Learning Research 18(1). 3451-3482.
[39] Pearl, [2009]. Causality: Models, Reasoning. and Inference, 2nd edn, Cambridge University Press.
[40] Sachs, K., Perez, O., Pe'er, D., Lauffenburger, D. A. and Nolan, G. P. [2005], 'Causal protein- signaling networks derived from multiparameter single-cell data', Science 308(5721), 523-529.
[41] Scanagatta, M., de Campos, C. P., Corani, G. and Zaffalon, M. [2015], Learning bayesian networks with thousands of variables., in 'NIPS', pp. 1864-1872.
[42] Silander, T. and Myllymaki, P. [2006], A simple approach for finding the globally optimal bayesian network structure, in 'Proceedings of the 22nd Conference on Uncertainty in Artificial Intelligence'.
[43] Soltanolkotabi, M., Javanmard, A. and Lee, J. [2018], 'Thcoretical insights into the opti- mization landscape of over-parameterized shallow neural networks', IEEE Transactions on Information Theory 65(2), 742-769.
[44] Spirtes, P. Glymour, C. N., Scheines, R. and Heckerman, D. [2000], Causation, prediction, and search, MIT press.
[45] Teyssier, M. and Koller, D. [2005], Ordering-basee search: A simple and effective algorithm for learning bayesian networks, in 'Proceedings of the Twenty-First Conference on Uncertainty in Artificial Intelligence'.
[46] Wauthier, F. and Donnelly, P. [2015], greedy homotopy method for regression with nonconvex constraints, in 'Artificial Intelligence and Statistics', PMLR, pp. 1051-1060.
[47] Wei, D., Gao, T. and Yu, Y. [2020), DAGs with no fears: A closer look at continuous optimization for learning bayesian networks, in 'Advances Neural Information Processing Systems'.
[48] Wu, C., Luo, J. and Lee, J. D. [2018], No spurious local minima in two hidden unit relu network, in 'International Conference Learning Representations'
[49] Yu, Y., Chen, J.. Gao, T. and Yu, M. [2019], Dag-gnn: Dag structure learning with graph neural networks, in "International Conference on Machine Learning'. PMLR, pp. 7154-7163.
[50] Zhang, B., Gaiteri, C., Bodea, L.-G., Wang, Z., McElwee, J., Podtelezhmikov, A. A., Zhang, C., Xie, T., Tran, L., Dobrin, R. et al. [2013], 'Integrated systems approach identifies genetic nodes and networks in late-onset alzheimers disease', Cell 153(3), 707-720.
[51] Zheng, X., Aragam, B., Ravikumar, P. K. and Xing, E. P. [2018], DAGs with NO TEARS: Continuous optimization for structure learning, in 'Advances in Neural Information Processing Systems'.
[52] Zheng, X., Dan, C., Aragam. B.. Ravikumar, P. and Xing, E. [2020], Learning sparse nonpara- metric DAGs, in 'International Conference on Artificial Intelligence and Statistics', PMLR, pp. 3425-3425.
[53] Zhu, S., Ng. I. and Chen, Z. [2020), Causal discovery with reinforcement learning, in 'Interna- tional Conference on Learning Representations'
12

---

Algorithm 4: Find path {WAL pa.
via particular scheduling for when is unknown.
Input: po 41a
8²
a²
Output: { Me Jk=
ā 416 +
1/ VE s.t. ā ã<o
2 W Hg f GradientFlow 0) for
k 1,2,...00 do
Let uK+1 E (2/ā)²
2/3 4/3
uk MM
WPe4 + GradientFlou( (g W pp.)
6 end
refurn (W Wpp 1k=
A Practical Implementation of Algorithml 2
We present a practical implementation our homotopy algorithm in Algorithm 4 The updating scheme for μa is now independent of the parameter a, but presented, the initialization for μo still depends on a. This is for the following reason: It is possible to make the updating scheme independent of without imposing any additional assumptions on a, evidenced by Lemma below. The initialization for po, however, is trickier, and we must consider two separate cases:
1. No assumpttion the initial choice of Hg matters.
In this case, if a is too small, then the problem becomes harder and
2. Lower bound on a. If we are willing to accept a lower bound on a, then there is an initialization for μu that does not depend on
In Corollary we illustrate this last point with the additional condition that a > 55/27. This essentially amounts to assumption on the minimum signal, and is quite standard in the literature on learning SEM.
Lemma 4. Under the assumption 11 140 solution to (6. i.e.
a²
a. the Algorithm 4 outputs the global optimal
lim W k-000
WG.
turns out that the assumption in Lemma is not overly restrictive, as there exist pre-determined sequences of {tu ) that can ensure the effectiveness of Algorithm for any values of greater than certain threshold.
B From Population Loss to Empirical Loss
The transformation from population loss to empirical can be thought from two components. First, with given empirical loss, Algorithms and still achieve the global minimum, WG, of problem 6, but now the output from the Algorithm is an empirical estimator à, rather than ground truth Theorem and Corollary [ would continue to be valid. Second, the global optimum, We, of the empirical Toss possess the same DAG structure as the underlying W.. The finite-sample findings in Section (specifically, Lemmas 18 and 19) of Loh and Bühlmann (31, which offer sufficient conditions on the sample size to ensure that the DAG structures of WG and W, are identical.
C From Continuous to Discrete: Gradient Descent
Previously, gradient flow was employed address the intermediate problem 7 method that poses implementation challenges in computational setting. In this section, we introduce Algorithm 6 that leverages gradient descent to solve (7) each iteration. This adjustment serves practical considerations. We start with the convergence results of Gradient Descent.
Definition 1. f L-smooth, if f differentiable and Vx, y dom(l) such that [V](()) ff(y))2 L x vl12.
13

---

Algorithm 5: Gradient Descent(f. "7, Eo,e)
Input: function f, step size 1/, initial point Wo, tolerance
Output: We
t-0
2 while Vf(W)))l2 edo
W W411 Wt Ww-n)
1+t+1
5 end
Algorithm 6: Homotopy algorithm using gradient descent for solving
Input: Initial W-1 (x- 8-1), 1Ao
a²
(1+,8)*
a² (1-8)" 8)
aa2+1133 (1-3)a
(1+0)²
7/ Pola +1)+3aa. 60 min Bupto Mo
3/2
Output: {W pr 100
W po <- Gradient Descent(9po: 70. W 1-1, co)
2 for 1, ..
Let μı (2 22/3 11/t
(a-ek 1/M4 134/3
Let ?) px(02-1)+334
Let EN = min min(Bau (aaux u
W "iiir + Gradient Descent(gpe: Me, Way 1 Ex)
end
Theorem 3 (Nesterov et al. 33). If function is T-smooth, then Gradient Descent (Algorithm 5 with step size = L,T, finds an c-first-oddrr stationary point (Lee. [[V S e) in (x+) - ///
iterations.
One of the pivotal factors influencing the convergence of gradient descent is the selection of the step size. Theorem 3 select a step size 77 L. Therefore, our initial step is to determine the smoothness of gu (W) within our region of interest, A {0 < < a, 0 S y S -
Lemma 5. Consider the function 9μ (W) as defined in Equation 7 7within the region A = (0 < S a. 0.0 S S ). It follows that for all μ 2 0, the function gu (W) is u(a² +1) + 3a- -mooth..
Since gradient descent is limited to identifying the stationary point of the function. Thus, we study the gradient of (W) = (W) + h(W), i.c. (gn(W) has the following form
(pt p(x - a) - y²x (p(a² + 1)y ap + yx²
(g)
As gradient descent is limited to identifying the stationary point of the function, we, therefore, focus on (W) l| & e. This can be expressed in the subsequent manner:
(V) se - S u(ı - a) x <e and - E S p(a² + 1)y - ap + yx² S
As result,
{(.,)) (W) < C {(x
na -
< x <
μa +
pa - E
< 5
μa +
# 3²
2 3وو
x² +
pla2² 22 +1)
12 ula² 1)
Here we denote such region An,
pa
aa +
μa -
pa +
An,e = ((xx)) y)
< r I
(10)
μ y2
u 1²
x2 + u(a² u +1)
z² + p +
Figure 6 and 7 illustrate the region A Ay,s.
14

---

=0.6. 0.009, c
Org XG2
Og (
VuI 11+
VR-g -1
ym sp
nyme
0.2
0.1
0.0
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Figure 6: An example of Ap., is depicted for a = 0.6, 0.009, and = 0.00055. The yellow
region signifies stationary points, denoted A,,, and defined by Equation (TOp. Al,,e the disjoint union of A and A Al... which are defined by Equations 2D and 22 respectively.
Piot -0.6, 0.009, -0.00055
0.10
1-11-
9-11-
0.04
0.02
Figure Here is localized illustration of A,,c that includes the point (z 81). This region, referred to as A1 is defined in Equation (2T).
Given that the gradient descent can only locate stationary points within the region A,.r during each iteration, the boundary of Ay.c becomes critical component of our analysis. To facilitate clear presentation, it is essential establish some pertinent notations.
15

---

ja
(Ila)
1 y²²
pa
y (1a2 + + x
(11b)
If the system of equations yields only a single solution, we denote this solution as (2) Vp)) If it yields two solutions, these solutions are denoted as (# fu)) (x U2*),, with an < Fy. In the event that there three distinct solutions to the system of equations, these solutions are denoted as (o)s Vi),, (177 V7"), (z "p"", where 1 a.* < xx < x
pa
(12a)
# ++1 y?
pa +
y
(12b)
maa² + +)+
If the system of equations yields only a single solution, we denote this solution as (stt,e: Mj,ee) it yields two solutions, these solutions are denoted as (2/ Vy.e). (aftee Mwee), with oy < Ty.e. In the event that there are three distinct solutions to the system of equations, these solutions are denoted as (#j Uj,c), (# Wü,c), (1 a Vyle), where Fne < Tne < Ffy
T pa +e E
(13a)
11 + y2
μa -C
y pla² + +)++00 x
(13b)
If the system of equations yields only a single solution, we denote this solu- tion as as (= Wi..e (. (2 My with Tne three distinct solutions to system equations, these solutions are denoted as (it Vi (ztiee ) ( where xme < <x <xi.e.
If it yields two solutions, these solutions are denoted Tpee_ In the event that there are
Remark 4. There always exists at least one solution to the above system of equations. When is sufficiently small, the above system of equations always yields three solutions, as demonstrated in Theorem 5 and Theorem 9
The parameter can substantially influence the behavior of the systems equations 12a 12b and ( (1b)). A crucial consideration is to ensure that remains adequately small. To facilitate this, we introduce a new parameter, 3, whose specific value will be determined later. At this stage, we merely require that 3 should lie within the interval (0, 1). We further impose a constraint on to satisfy the following inequality:
S Bap
(14)
Following the same procedure when we deal with e = O. Let us substitute obtain an equation that only involves the variable y
12a into 12b.. then we
c/p
(na--)2/1 f)2/uu
re(y: p1)
(a² (12
(15)
y
(y² +1)
Let us substitute
12b into 12a
then obtain equation that only involves the variable
e/U
(wa + +c)²/u
Le (a (1)
(16)
(u(a² + 11 + x22²
T
Proceed similarly for equations
(3a) and (3b
- e/F
(μa + +c/2/1
Te (v:
(a² +1)
(17)
11
(y² + 11)2
16

---

L (x: a + a///
(na -))2/p
(18)
(u(n² 11
+ x²)²
Given the substantial role that the system of equations 12a and 12b play in our analysis, the existence
o in these equations complicates the analysis. this can be avoided by considering the worst-case
scenario, i.e., when Sap. With this particular choice of e, we reformulate (15) and (16) as
follows, denoting them as r <) and ro( t) respectively.
()( a(1 + 9)
(a² +1 paa - $)
(19)
(2² + 11)²
y
a(1 33 x
pa 2 (1 + (32
to(x; pb)
(20)
(u(a² + +1) + 22)2
The functions r, (y: 14), Tr_ (y: 1r), and ( (u) possess similar properties to rly; (4) as defined in Equation (8), with more details available in Theorem 7 and g Additionally, the functions (x(:;; μt),
to (x; and to (z; (i) share similar characteristics with f(;; (1) as defined in Equation with more
details provided in Theorem |
As illustrated in Figure 6 the e-stationary point region Ap.e can be partitioned into two distinct areas, of which only the lower-right one contains (17 Ff.) and it is of interest to our analysis. Moreover,
(z e) and (2 Vf..e) are extremal point of two distinct regions. The upcoming corollary substantiates this intuition.
Corollary 3. p < (T is defined in Theorem 5(v)), assume e satisfies 14
G satisfies
1+0
S
1-3
a² + L, systems of equations
(1bb) at least have two solutions. Moreover, An, A U
A A.c
p,c - An, n {(x,
x 2 xeee 1/ S Wp,.e)
(21)
A = An, n { y)
x R.:: y 2 yptel
(22)
Corollary suggests that Ap.c can be partitioned into two distinct regions, namely ALe and Ah..- Furthermore, for every (z, y) belonging to Al.., it follows that x 2 xyee and y S ye.r- Similarly, for every (z, y) that lies within Ab.c. the condition S x" and y 2 yoe holds. The region A represents the "correct" region that gradient descent should identify. In this context, identifying the region equates to pinpointing the extremal points of the region. As result, our focus should be on the extremal points of A, and e specifically at (3%,.e V).e) and (1 V/ Furthermore, the key to ensuring the convergence of the gradient descent to the Al. is to accurately identify the "basin of attraction" of the region Al.e- The following lemma provides region within which, regardless of the initialization point of the gradient descent, it converges inside A! Al..c
Lemma 6. Assume < T (T is defined in Theorem 5(v).
14.9 I-3
< ca²+1. Define By = { y)
a.e <x a. S < V. J. Run Algorithmm 5 with input f = g (x. p(a 33 Wo
(z(0), y(0),, where (c()), E By.e. then after at most 2(w(a" +1)+3aa (on (=(0),((0)) (x iterations, (zt, y1) E A/,.c
Lemma can be considered the gradient descent analogue of Lemma 2 It plays pivotal role in the proof of Theorem 4 In Figure/6) the lower-right rectangle corresponds to B Lemma 6implies that the gradient descent with any initialization inside Bp will converge A A +1,001 at last. Then, by utilizing the previous solution W WeLLee as the initial point, as long as it lies within region Bpei1ee the gradient descent can converge to Al p which is stationary points region that contains W* W thereby achieving the goal of tracking W W Following the scheduling for pa prescribed in Algorithm provides a sufficient condition to ensure that will happen.
We now proceed to present the theorem which guarantees the global convergence of Algorithm 6
2
Theorem 4. Ifō E (0,
E (0, 1),
< (1 - 5) + 11) and μr satisfies
1-3
a²
a²
(1 3)
c2(1 (1115)3(1 B) a²
5 no
4(a2 + 1)3 4(a² + 1)3
8)²
4
(1+ (33)²
17

---

Set the updating rule
< = min (BauB 1kk
3/23
2/3 (a + ck/ )2/3 (a E ph 34)
ph+1 =22
Then pk+1 S (1- (1. Moreover, for any Editt > 0, running Algorithm Gatter K (Po, a, ó, Edist) outer iteration
|| mu.ex - W S Edist
(23)
where
K (HHo a, 8, Eami) 2 In(1/(1 5))
max In 8oo
nd0
In
a² - 144) In(
72p0
3(1 6)mb) ), 2 a Edtin
In(
16656på 1 In( 46656p0
Edist
a Ediat.
The total gradient descent steps are
a. 2( (a² + 1) 30 19 (Wpr.)))
gux+1 (W (
3(4
6)
216
216
2/3
72
<2(uo(a² +1++ 3a²)
max
(1-(1/2)//)aa²)
guo (WI8)
Edist
OEdist
oEdist
82
(1-
50 (baa? a² po)
1
3oa8
Ediat
d
Q
Proof. Upon substituting gradient flow with gradient descent, it becomes possible to only identify an e-stationary point for (W)) This modification necessitates specifying the stepsize 7) for gradient descent, as well an updating rule for μ. The adjustment procedure used can substantially influence the result of Algorithm 6 In this proof, we will impose limitations on the update scheme (ss the stepsize Thk and the tolerance EK to ensure their effective operation within Algorithm 6 The approach employed for this proof closely mirrors that of the proof for Theorem
[ albeit with more careful scrutiny. In this proof, will work out all the requirements for μ4, e, r.. Subsequently, will verify that our selection Theorem conforms to these requirements.
In the proof, we occasionally use μr, pky When we employ inequality or equality holds for any p, Conversely, when we use Pk: Ck. it indicates we are examining how to set these parameters for distinct iterations.
it signifies that the given
Establish the Bound yice 2 V First, let us consider re(v/: S i.e.
re(VTaH1) a + t/1u
a(a - e/m)²
(a² + +1
< 0
V/
4μ²
This is always true when μ > 4/a?, and we require
4 when μ S
S 2/2 + ap 2222//22 1²
7
Now we name it condition 1
Condition 1.
S 24 3/ + aj 2 22,55/2 2 u²a² when μ <
Under the assumption that Condition is satisfied. Since r. (r μ) is increasing function with interval y E [JAb,e, )u and we know y/b. S V I yub.r and based on Theorem 7 we have Bbb Vu. Yub.e. Te Te(T;; μ) S (e 0. Therefore, v vp..
18

---

Ensuring the Correct Solution Path via Gradient Descent Following the argument when we prove Theorem I we strive to ensure that the gradient descent, when initiated at (I (Zpe.rr Bon.e)) ) will
converge within the "correct" 6++1 -stttnattyy point region (namely, [[V (W)l22 61+1) which includes (m (. For this to occur, we necessitate that:
(1)
(a + <x/Mx)) 1/3 (1)
(5)
Mk++1 (2p
) (a-Ex/uu))
(a -
Dpu.ch"
ypu.e (24)
Here (I), (5) are due Corollary B} (2) comes from the boundary we established earlier; (3) is based on the constraints we have placed on P and pk+1, which we will present as Condition 2 subsequently; (4) from the Theorem 7(n) and relationship yaa.es < gib.pa.c. Also, from the Lemma 9 maxpsr 2/e S minp>0 me.e- Hence, by invoking Lemma 6 we can affirm that our
gradient descent consistently traces the correct stationary point. Now we state condition make happen,
Condition 2.
(a ck/14)213 (a E 24/3
- (1-8)ux 2 p4+1 2 (2p2)2/3
In this context, our requirement extends beyond merely ensuring that 14 decreases. We further stipulate that it should decrease by a factor of 5. Next, we impose another important constraint Condition 3.
3/2
E/ S P3
Updating Rules
Now we are ready to check our updating rules satisfy the conditions above
cs = min (aiu 1k 3/22
Mk+1 = =22 2/3 (a + es/ 1A )/33
(a ck/ Mk 4
Check for Conditions First, we check the condition 2 condition 2 requires
(1-8)pu 2 (2 (a + cx/14x)2/8
Ml (a + cx/
(1 8)3
(a Ex/MU)) 2k/Mxx)4/3 4//3
(a ek p1))
4
Note that S 3aapp < ope
(a + c 1xx)
(1 3)
Fk (a c// 1z S 1* (1 B) a²
Therefore, once the following incquality is true, Condition 2 is satisfied.
(1 + 1)2
513 U4 4
a² () 3)
F1 (1 - 3) a
(1 +1+1 B)²
Because S. 140 S a- (1- (1-8))(11)* from the condition we impose for μa- Consequently, Condition is satisfied under our choice of cg.
7 (1+3)²
Now we focus on the Condition I Because c a.5pu, if we can ensure a,344 S 2/1 3/2
+ apk
2 2ap 5/2 uda² holds, then we can show Condition is always satisfied.
a33uk 22 3/2 + app 2/ 2apu 12
5/2
2 2app
n 22 3/2 + - -3)apk
4 4(2a/1
5/2
<4mm + (1 1-3))322 2 + 4(1 3)ap 5/2
0 544 + +1)Ø + -B))a"² 4(1 B)apk 6/2
0 0 + 1)m4 - 4(1 + 3)auk 1/2 + (1 - 3)20² when
0 m S 4/a²
(1 3)a 1/2 11-8)
5/4
(a2 + 14 4(a² 11)
19

---

We also notice that
(1 + 2)
(1 (1- 8)
(1+
2
5
(i-
S s +1
(a² (a²+13 +1)
4(a² +1
[++8] 2 Because 1-3 1 (1- - 5)(a2 + 1). the inequality above always holds and this incquality implies
that for any Pa 2
(1 B)a 1/2 3)
F4
(a² 1) 1k
4(a² +1
Therefore, Condition 2 holds. Condition also holds because of the choice of 6/.-
Bound the Distance
Let 72/a², and assume that μ satisfies the following
J S mint (1- (11122) ,B2a²}
(25)
Note that when j satisfies
25 then 1/3/2 S Bap, s0 = 1/3/2
1 (1 (1/2) c
a²
1 - (1/2)2/0)
a²
s
72
4
c/p= = vp s 2
(26)
Then
1
1la +
(aa - (/1)(1//)11c cx);
cuu
(M(0² + 1) + (a - 1(a p(a+c/11))
-//)=/)/()))0 (2 cq1)322
cp
cH
(u(a² + + (a -1x)((L-cc))
(a e/1)/
20h 11
(1 cu)
(a + a/2)²
20h
(a
a/2)a(1 cpi)
36
=p
a²
a²
q41)4
72
36
=p
0
a²
a²
cui)
Then we know (a - e/p) - cp) < x/c Now we can bound the distance [ Well. it is important to note that
Wmecc WGll
= (Emmoca a)2 + (Yuupca))
S <max VV - aa22 + (yrr.ce) v (oterra a02 (2
We use the fact that Teeees < Fppiica < a, Fpnnen < x on and Upeice < Vjueies Next, we can
separately establish bounds for these two terms. Due to (24 Mpsiie < (2/2) 1/3 (ates/pa) 1/3 1 and (a- - <x/ cpis) peice
(a-cx/14)373
V (Irunee a)2 + (Bii S V / + (a- - (ct/ee))(11-0 cpx))²
Given that if xn furt S a, then / (1fherr - a)2 + (U
)2 2 v (1imee - a)2 + (Viferr)) )a. There-
fore, if ar 2 a, we can use the fact that ot S +
en
In this case,
us
V (stere - aa) + (U )2 2 VeEtt + (en/ua)²= / /uat1 MA S /(2 - p8
20

---

As result, we have
W Weice - WGll S mmxxf M6+1 + (a (a Ex/Mss(1 cum)) v(2 8)
M8+1 + (a - (a (1 cux))² <(1 8)pk (acHu + 1/4 cut 3//
<(1 6)pz 3 + pk + e (11)))
=(4 8)14 + 22 + +32
[ Me.Ek - WGll S maxt VuN+I + (a (a ex/ sx/Ma)] - (1cuu)))) V(2 8)
<max{ /(1 8)/44 + c + 3c2 v(2 - (4-
-V (4-4) + 30 11 + +3c²
Just let
(4-8) $ (4 11-8)(1 5)5Mk 4u
Edist
k 2
In(3(4 - -8)10/800001 D)/p//Ecistt
(27)
3
In(1/(1 8))
In(46656pz/ (a" 2))
3a 3a²c2 14 S 2 5)2k pg
Edist
k
(28)
21n(1/(1 (())-))) 8))
In In(46656pB/ 2))
3e 3e²p1 1 s 3c2 8)3h p
Edist
(29)
31n(1/(1 8))
We use the fact that μ E (1 111-84 830$1b Ng. In order to satisfy 25
a²
In a11-(11/2)
72pm
S po(1 - 8) - 72 (1-((11) (1/2)1114 = ソス
(30)
In 1-1
In (po/ (B2a))
N& - 1o - 5) & f8² = トト
(31)
In 1
Consequently, running Algorithm after K (No, a, 8, Edist) outer iteration
Wuneex Wellz S Edist
where
K(po, a, 8, Edts)) 2 In(1/(1 - 0))
1
max In
Pd z2a²
In
7240
In(
3(4 - 87p0
1 In
46656/6
1
In
46656p
a² - (1/2)1/4)
22
a²
3
a
By Lemma 6. l iteration of Algorithm need the following step of gradient descent
2(ua(a² 11 + 3mm) (Wmnraa
Irs (Wiy, 19821))
2
21

---

Let (Po, a, 6, Edis)) satisfy l & (.s S B8a < (R(o ( Hence, the total number of gradient steps required by Algorithm can be expressed follows:
a. dadee 20u + + 3o Kouge me SMR11 (W 1.6111))
%
(ome (Winer) Ripo.odisai decait)
(W (Wma11ick1))
Ripo dinian)
(Bok (Wygle)) ) 90u11 (WML)1
Sa(po(a" + 1) 30²)
Σ
Σ
o
of
Ripoio Siedua)-1
(Bok (Wuges) ) - Dok (W 1R
Loex (Wprier) ) BuN (Was 11:5811
=2(no(a" ? 3a²)
K(po Siediaa diddiot!-
(Buets (W (ggee)) Bor+1 Weet 1)
King Rediia)
(Bee ( (Waries) ) 9ux-1 (W pett ck41
3 +) 3
89
k=0
K Klppo.a..eeaue
lang- (Wpa.cg) Bua+i (W Pe41
loness (Wezice) gman (W (WPs41 28411
$2(po("" 4 1) )a²)
Σ
Σ
8og
=2(po(a" P 1) )
10k (W (1
goat
Pi Ailpaian
S2(mp(a" + 1) ) 66a6
(1en (W 8uki1 (eT
=2(mo(aa + 1) 3
(Ano (Wug.sa) w. H (W KKpo aisedia
P (klegreriaam))
S2(Ho(a" + )
Sup (oa
Note from and (30. the following should holds
(:s minf Edist
dedist
2/3
a² ,,800²
(11/2)////
3(4 144-81 8) 216 216
72
Therefore,
2(He(a² 1) + ) (W (Wpn.ca)
-p (W
H4+1 cc-11 ))
2
A=
3(4 8)
216
216
2/3
72
<2(1u(aa² +1) + 3a²)
+ max
Ino (W
36a5
Edist
2
AEdist
aEdist
8 (1 (1/2) )))aa?
22

---

D Additional Theorems and Lemmas
Theorem (Detailed Property of r(ys p)). For r(yi p)) in ( then
① For 0, limy-ot r(y; v)) r( (( ()) < 0
iii) For p > 0, r(vA, (4) < O.
(ii) For p >
dr(y; 1s)
< 0
dy
For < p S o²
dr(y: p)) dy
Vlb < y Yub
(32a)
dr dr(y; 11) dy
0
Otherwise
(32b)
where
(4p)
(4 (a + v0213 (4p))//3)
(a 1/2/3- (4p)2/3) Yub
Mbb 2
Moreover,
Vlb S Vu S yub
(iv) For 0 < u <
let p(p)) r(yub, 1)). then p' (μ) 0 and there exist unique solution to p(u) denoted as Additionally, T < 2.
(v) There exists > 0 such that, Vy > the equation r(y; p.) 0 has only one solution. At # = T, the equation r(y; μ4) = 0 has two solutions, and Vu < т, the equation r p1) = has three solutions. Moreover, μ < o
(vi) Vp < t, the equation r(y: p) has three solution, i.e. U < Vp < Vpa
dy du
dy
dy and
lim 2 0, lim 1 n=0 #-0
0, lim U/
dp
dp
+
Moreover,
5/ yhh V yp yub Uu
Theorem (Detailed Property of ((2; For (x;; p)) in 9 then
ii For > 0, limg-+0+ ((:; 11) = xi, l(a, 14) <0
(ii If < (o(va=TTaa) 2(a2+1)
(a(vT4I++a)
then w 1), < 0.
Or
(iii) For a²
+11
dt(x:; x p) dx
< 0
For 0 < 4 S +11
dt(x: p1)
Rib < Fub
(33a)
dt(x; p)
Otherwise
(33b)
23

---

where
(4pa)2/8(1
V1 2
(4x)1/8(a241)
(4pa) + v1
(4n)1/3(a2+1)
a473
a4/3
21b
Lub
2
Moreover,
Ibb S VM(0² + +1) zub
(iv) For < < F1y and let g(p) = t( p.), then d'(p) > 0 and there exist a unique
solution to q()) 0, denoted and < Ta=+IJ S 27
(v) There exists T > such that, Vp > the equation (x; n) = 0 has only one solution. At
μ = 5, the equation t(x; p) O has two solutions, and Vp < т, the equation t(x; μ1) has three solutions. Moreover, < T[ S 27
(vi) Vp < T, Mx; p)
has three stationary points, i.e. x < x < x
dx
dx
0
> 0 and lim z = a, lim 4 0, lim 2/ 0 1000
du
dp
Besides,
a(Va² + +1-0) 0)
and alva2+1+0 + a)
max PST
# S
S min Pn0 P
2Va² +1
2Va² +1
It also implies that t( ol voa
2/0*+I ;1) 2 and max 2xr < min x7
Lemma 7. Algorithm with input f (x, y). 20 (z(0), y()) where (z(0). Cus in (1). then vt 2 O. (z(t), v()) E Cpa- Moreover, lim, (x(t),((t) = (2), V/
Lemma 8. For any (z, v) E C,s in 441) and (I, v) # (2f, v7,)
gp(2,,) g vfi))
Theorem (Detailed Property of rely: p)). For (y):: (1)
then
i) For > 0,6 > O, y Te (e 11) = 0o, y( ग O
(ii) For p > (a-e/p)² Matt/aazz then dr.(w:) dy
(a-e/a)" O. For 0 < f [()=/11)²
dr. dre p)
Mbbue < y < Dub.uue
(34a)
dy
dr dre 1)
< 0
Otherwise
(34b)
dy
where
1/3
(a - e/1)2 (aa=c/u)²)
1/3
(a e/H222
2/3
ylb,uu,
2
E/13
V
e/11
244))//3
24
(a - e/
1/3
(a ( e/)22
2/3
yub,u.
+
144
2
a - e/μ
a - e/μ4
Also,
(a 1/3
2/3
Mb.
(a 2/3 2/3
ybb S V yub
Theorem 8 (Detailed Property of (o(w: p)). For Toly p) in 19 then
24

---

ii For > 0, > O, lim, - re (4) 00
i For (1(10)) 4(1+8}2 then dra(wi))
a1(1-0)
< 0. For 0 u 4(143)2
dy
dr drøly: p))
0
nlb < < Bub.u.8
(35a)
dy
dr dr p dy
Otherwise
(35b)
where
(4
[al1 aal1-33* 833
1/3
(44)2/3
1 B
2/3
Dbb
1
1 1 1
a2/3
(1-19)2
(4/)1/3
(a(1 (8)² 1/3
(4/)1/3
1 B (1-
2/3
1 V 1 -
1 1+
a2/3
Also,
(4u)2/3 (1 + 3/1/3
u. S 2/3 (1 /3)233
ybb S V S yub.p..
Theorem 9 (Detailed Property of ta 14)). For tz(2; n)) in (00). then
i) For > O, lim (o(x; 11) = 0o, (p(a; u) < 0
a² (0+1)*
(ii) For p
4(a2 (B-1)2
dl d(
< 0
dx
a² (9+1)*
For < KS +1a (3-1)²
d p)) dr
Fib. 2x < Fub
(36a)
dt ( (x: p1) dr
Otherwise
(36b)
where
1
4azu(1 + -8)
1/3
(3 +1
1 8
2/3
=
2
1 1- 3
1
a2/3
(1+ 3)3
4aju(1 + -3)
1/3
(4p) + +1)
1- B
2/3
Dub.uu,8
2
1- 8
1 1++
a2/3
(1 3y²
(iii 1f0 < 3 < v(a2+1)-1
ra(1;1) v(a+=1j+1' then there exists a Tg O such that, Vp TB, the equation
has only one solution. At p TB. the equation r3(2;;) O has two solutions, and Vi < Tg. the equation ra (2; (1) = 0 has three solutions. Moreover, <
2 1)+
4/a +11+ 11-1)
(iv) 0 < 3 <
+1++11+11 then Vp < Tg, ts(z: p)) = 0 has three stationary points, i.e. FB < FB Fy.g. Besides,
a((1 3) +1 () (2² 1 (B
max - I
MSTa
2 +1
+1+ + (1 (a² +1)=
min p25
2/a²
25

---

It implies that
max PS
min F
Lemma 9. Under the same setting as Corollary
maxa MST
n.e
< min THe 1SO
E Technical Proofs
E.1 Proof of Theorem 3
Proof. For the sake of completeness, we have included the proof here. Please note that this proof can also be found in 33
Proof. We use the fact that f is L-smooth function if and only if for any W,Y Y dom(s)
f(W) s f(Y) + (V w) + / ll W 1
Let W W and Wt. then using the updating rule W W (V
f(W++1) sf(W)) (Vf(W") w - W3) + 2
| wº1
=f(W')- - 7.
112 + 2L
/Vf(W))
1
-s(W")- 2T. VVS(W))
Therefore,
f 2 - -
2L(((W9)- - f(Ww))
min 0st5n-1
[V(
2 - S(W*))
2 n ([([(W)-S(W')
min 0st5n-1
||//(W
E.2
Proof of Theorem 5
Proof.
For any y >
a²
lim
r(ys 11)
lim
(a² + 1)
2
r(
< 0.
a
1
(ii)
a²
r(TT.1))
(a² +1)
Vμ
4μ
a
2 2 02 <
Vp
a
26

---

(iii)
dr(y:1) p1)
a
+ (2/² + 1)3
dy
y2
4a pas a(y² + 11)
(32 + 12)
a((4aj))213, 22 + (Aap)) y(y + (1) + (y² + ((aau)11 - 11)
212 + 11)3
For μ 2 4: ((4ap)) y²² 11) <0 dr(yan)
0.
dy
For < 4. lib < < Yub. ((4ap) 1/8y y- y² p1 > 0 +
dr(w:)) dy
> 0. For <
< ybb or Yub < (Aaa))113)))) 1/² - 14) S 0 + dr(w
< 0.
Note that
dr(y: p) + ((4ap)) y- 3² 1)) + (4ajz) 1/3 y /
du
y
The intersection between line (4a and function y + & are exactly yob and Yub. and Vib < VP < Bub.
(iv) Note that for 0 < u < a²
Ör
-a? 11² p
and Yib < vH < Yub
0p
(u+ + 33)3
dr then
< Let p(1) = r ryyub p.), because
Or By w=
= 0, then
8y
yub
dp(w)
dr dr(bubi 14
dr
dyub
dr
Ör
du
dp
By
djp
0y
Oy
Bub
B=
y= Bub
Also note that when / 4
a²
Uub
//h, p(1))
r(uub, 11) r(vu,, 14) < 0, and also
a then
(44)
2a 1/3 = (4pa))13
Yub
2
Thus,
aaa 2
r((44a)1133 (4) = (4 ((4ua)2233 + u22 - -(a² +1)
(4a)/3(4 n1/372 - (a² + 1)
(4ua)2/3
a²
(4a)1/3) (a² + +1)
1//3 (4a)1/3
(10)473 it easy when μ - 0, ((Aua)2/3 11) co. We know (4pa)) u) oo as μ 0 because of the monotonicity of r(y: p) in Theorem ( Combining all of these, i.e.
Because (2a)175
2
r(yub, p))
dp(p) O,
lim p(u) = oo, n-0
pl T <
dp
There exists
4 such that p(T) =
(v) From Theorem 5(iv) for P > then p(n) = r(yub: > 0. and for μ = r, then
p(1)
r[uu: O. For μ < r, then p(u)) r(yub, 1)) < 0, combining Theorem get the conclusions.
5()
5
27

---

(vi) By Theorem (v))
4p < т, there exists three stationary points such that 0 < U/ < Ylb <
vp
f Yub
Up Because
dr(wa)]
dr(3:w)
0, then
dy
dy
ym
Bub
dr(y; p)) dy
dr(y: 11)
dr ni)
O,
0,
#0
dy
dy
y=
y=
By implicit function theorem [14], for solution to equation r(y ()) = O, there exists unique continuously differentiable function such that y y(u) and satisfies r(v(u), p) Therefore,
Ōr
2 D (1r + y2)33
Ör
a
4a py
dy(p) du
0r/0p
dy
y²
(v2
Ör
Therefore by Theorem 5(11)
dy
dy
dy
0
0
<0
du yー
du
dp 1-
y=
Because lim Vlb lim Yub 0, then lim 00 D/ lim 87 0. Let us
consider r( a (1 - cp),, 1)) where e = 32 La 1a2++0 +1
and / < 20
r( a² + m))
uaaa 2
cu) + 11)2 - (a² + 1)
-cu)
a²
a++I
(a2+1)2
=(a² + 11(
qp
uua
cp
la 41'117 (1 cpu) + 11)
ua?
20(a² + +1)p
11=+112 22)
16(a² + 1)4
=c(a² + +1)p=
a²
16(a² + +11
#1
a
By Theorem 5(ii)
then . 1
cui
< 4/ then
lim
a² + (1 cp), S lim U/ S -
a² +1+1
+++
a² +1
Consequently,
lim Vp a n=
a a² +1
E.3
Proof of Theorem 6
Proof.
For μ u>000000000 0,
a²
lim t(2;p) lim x++
1 =0
--
u(a² + 1)
2a?
t(a,u) = (u(a² + 1) + 22)² <0
(ii)
d² 4p(a² + 1)2
1( (Vuu(a² + 1),1)
a
Va² + Vp
28

---

Ift(y/ + 1), O, then
1 (a² + 1/3//2
+ 22(a² - = m= a( a(va2 + F F
2
Vp
2(a² +1
(1120==110)) or 14 (alva2+1t then t( (Vu(a² + 11, <
when /
2((2211)
(iii)
dt 11 dæ
4pa²x
x²
(n(a² + 1) +233
4uaa2 3 a(u(aa² +1)
x² (wla² + +)
a((p(a² + +) + (a(a² + + (4 1/3 x + (4pa) 2/332 x²) ((4 1/3 x- -u(a² +1) 11) - 1))
x² <2((22 +1) + 22
For p >
2 +1p then (4ua) (/ 11a² + - x² < 0 +
d((a..p))
< 0. For < a²
dz
and Øib
Tub, then 1/3 (4 +122+1)22"
dt(z.ja
For ua-
dx
I < Tbb or æ Tubs (4ua) - uia² + 1) - -x² < 0
d((x,)))
< 0.
dx
We use the same argument as before to show that
Fbb < VMla² + 1) < Tub
(iv) Note that for 0 < μ
a² 1(a2-1)3
01 2² x² uıaa² 1)
and Tb < Vi(a² +1)
Fub
0μ (m(a² + 1) + 2²²33
then 0t
> O. Let ((p)) = +t(iib: p.), because
at on -=zm
= 0, then
dq(p)
dt(2mn, 1)
Ot
dayb
Ot
0t
> 0
dp
dp
dr
dp
z=
r=
2
(4paa)*
(4pa) 1/3
Note that μ I(a Lub Fbb
4(
+11z
1 >>> 0.
When j (a(va==T=o) then (( u((22 + 1), < <0 Theorem 6(i) It implies that +1a
2
g(u)) O when f - 0+. By Theorem 6(ii)) ((u) ((01,11)) < ((V p(a2 1), < 0. Combining all of the theses, i.e.
dq(w)
0,
lim ((4) < 0, q( 4(a² 113 > 0 #
du
There exists 7
4(a 1)s. ((7) 0. Such 7 is the same as in Theorem 5(iv)
a²
(v) We follow the same proof from the proof of Theorem 5(v)
(vi) By Theorem 6(v) Vp Po. there exists three stationary points such that 2/ < 21b
z < Fub < x < a. Because At(2gp)
((( dx
0, then
dx
2=
dt p1) dx
dt(2;; n
dt(z;; p))
* 0,
# O,
#
dx
dx
xx=z
x=
x==z
29

---

By implicit function theorem [14]. for solutions to equation t(2; O, there exists unique contiuously differentiable function such that x(p) and satisfies t(z(1), μ1) 0. Therefore,
1 +11=2+11) (p(a +
dx
ōt/Op a²
dp 01/0r
4ua (p(aa +1)
Therefore, by Theorem 6(11))
dx
dr
0
>0
dp
du 2=2
x=
Because sxx Jib < z < Fub and lim
lim=-+0+ 4-+0 Fib lim Fub 0.
lim x #
... = 0 lim 1-+ #
Let us consider t(a(1 cu), 14) where c= 32 and < 2c
t(a(1 - cc));:
u
a(1 - -41)
(n(a² + 1) + a² cy1)2)2
c/
pa
- cfi
(u(a² + 1) + a2 c1))312
ua?
2c14
(a²(1 - cc))222
16 2cu
a / > 0
By Theorem 6(11) It implies
a(1 cu) S x/
taking μ + 0+ on both side,
a lim a(1 - cp) S u0
lim x S a
Hence, li 5/1 a.
When j r, because t( = 0 and Fub > Vu(a² + > Fbb. t(x; (1) is increas- ing function between Fubb] then t( t(vula² + 1):1) > t(rib: p) 0. Moreover,
1( u(a² + 1), aih and x are continuous function w.r. 14, 38 > which is really small, such that and xt Zbb. hence da -=
δ and ((yu(a² + +1,,1) > 0, ((II,, 1) < 0 (by Theorem 6(iv)
< 0. It implies when μ decreases, then 2r increases. This
relation holds until x = V Vl(a2 +1)
16 14) = 1( ((Vii(('' 11a2 + +1)p)= 0
2
*Ø = (1
+1 a)
+11
and Vu(a² + 1)
a(vo=+I-a) 2va*+1
Note that when P
<
alva==1 2
2(a
t( + < 0, it implies that xn
Vula² 1) and
da
> 0, thus de-
4p
x=zz
creasing j leads to decreasing ot
We can conclude
max 1 < a( +1 a) UST
+1
30

---

alvo=_I-a)
al alva++I-0 2(a=+1)
Note that V.s.t.
2(a
50
t({ [(v=++Iaa)
11)
2 O.
2(a++1)
Note that when / a=+1 i.e. (x
a²
uuaa² + 1) then
dr
du
-=
It implies that when j decreases, zz also decreases. It holds true until x vu(a² + 1). The same analysis can be applied to zn like above, we can conclude that
min a(va² +1 a)
2v +1
Hence
alv a(va² +1- a)
alva² +
max 1 MST
min
2v +1 2v +1
p20
E.4 Proof of Theorem 78 and 9
Proof.
The proof similar to the proof of Theorem and Theorem 6
E.5 Proof of Lemma
Proof.
[ + y/²
2zy p(a² + 1) + x²
( v)
2ny
Let \1 (z,)),,))), A2 V 9p xr y) be the eigenvalue of matrix V² 9u (2. y). then
( \ (2,8))
Tr 9/ (z, u)) - = + x² + u(a² + 1) + x² >
Now we calculate the product of eigenvalue
\v A (2 v))
>2 9
det((V'g 9m, (W))
-Ip + y² (p((2 +1)+ 22) 4x
μa pao 4x² 2 > 0
x
4//13/3 > zy
2
am ,//,2/3 >
ap
<>
2
+y + 11 > (44 1/3
Note that for (st y (x 3/ ", they satisfy 11a and (1b), this fact is used in third equality and
second " By 32b). we know A1 ( 9r (z. v)) >2 22 (x, v)) > O for (2i Wi). (1771 w
and a1 2 gr (z. 5) Ae ( (z. v))
for (z v/ ). then
\ gr, v) 0. >2 V g.,
0
for (z Vi)), (2rt""
Aa 9 or (v yn (z, v)) < O for (1 u12))
and
(x, = 0
Then (2 w/), (x7" ")*" are locally minima, (I vx²) is saddle point for g₁ (W).
31

---

E.6 Proof of Lemma
2
9L
Og A
Figure 8: Stationary points when μ < 7
Proof. Let define the functions below
(2) Va(
0 << S
(37a)
pa
U)u((x) ula² + 1) x
0 x <
(37b)
x1)(3)) = pa
0< < क+
(38a)
y² 4
xn2(z)) v 11( - - -aa + 1)) 0 0<g< < 7.4
(38b)
with simple calculations,
Yul 2 Uu2 + 1(1; u)) 2 2 0 20 20,
[x xxp
and
Pp1 2 F2 + r(y;u) 0 y E 11 3
1y
0² +
Here we divide
B, into three parts, Cal, C22, C33
Cμ1 =((x, , < < t Vp < マ < // u <r < <ソ< ","""
(39)
Cμ2 =x y) < 2/, 0 $ < y
< yn
(40)
Cus (00x((x)))))))) y>
< S t/ Uu
$ y و y
(z
x
a. $ yu2)
(41)
Also note that
0g, v)
Og (وا v) Öy
v(r.y) E Cp1
0,
Ox ag v)
Og (2, v)
v(a,y) E Cpa
0,
0
Ox
Oy
The gradient flow follows
8g.. (e()
2'(1) v((1)
2ig 0g,, (1(1)).9(()) dy
(z)()(t).
32

---

then
(2'(()
v(x,y) Cμ1 7
< 0, V || > 0 (v(11)
(42)
('((t) v(())
V(, Cμ2
> 0, VV > 0
(43)
Note that || V is not diminishing and bounded away from O. Let us consider the (z(0), (0)) E C,i, since Vg (x, y) 0, (g,(x, v) < 0 in 42 and boundness of C,a, it implies there exists finite to > such that
(e(to). pa (z((t,, v(t) Cşa for 0 St <tt
where αC₁n is defined as
a = [(r, wa < x u y wnn) u (() < S 4, y = mea) ) S Cμs
For the same reason, if (z(0), (()) E C22, there exists a finite time t > 0,
(c(ta). 2C (z(t), v()) C22 for 0 S t
where OCpa is defined as
aC = {1, | < S x yuut < x = yau } E Cua
then by lemma
7 oo (z(1), v(()) = (1 2(1)
E.7 Proof of Lemma
3
Proof.
This just result the Theorem
E.8 Proof of Lemma 5
Proof. Note that
H y² 2ry
2xy ula²
0 ula² +1)
y²
2ry x²
V
(W)
1) x²
2xy
Let . llop is the spectral norm, and it satisfies triangle inequality
(W) s =((t ulaa + +1)
y²²
2xy x
+
2xy
op
y²2
20y
=m(a² + 1) +
2xry
The spectral norm of the second term area A is bounded by
(x + 3²) + / (x + 2² y2x
2aa + v 4a + 12
max (x.
-a
2
2
We use x² S a 12 < in the inequality. Therefore,
|| gp (W)ll s 3a² + u(a² +1)
Also, according to 15 33 for any f. if 2 exists, then f is L smooth if and only if 2 - L.
With this, we conclude the proof.
E.9
Proof of Lemma
Proof. First we prove Vt 2 0, (2((t,)((t)) Cus, because if (z(t), vt)) & Cps, then there exists finite such that
(z(1), E 0Cps
where 8023 is the boundary of Cus- defined as
aCp3 (z, y yus (z) y y)2(x), <x 5 5
33

---

W.L.O.G, let us assume (z(0), (0) E aCps and (z(0),3(0) # (z 8/2) Here are four different cases,
if = yn (x(00) < x(0) < z
if y(0) = ynΛ (2(00) z00) <
(g,(c(().2())
if y Una (,,( < z(0) < <
if y yas yus (2()) < x00) S
This indicates that - (x((t,,), v()) are pointing to the interior of C4a. then (z(t), v()) can not escape Cya- Here we can focus our attention in C23, because Vt z 0, (z(t), v(t)) C3s. For Algorithm D
df di
V[(Zz) = IV
In our setting, v(r,y) E Cps
{
V (2, y) # Vgu (z. y) =0
(z.0) # (o yin) (z. = ( vf())
so
dgu(a((t).
V < 0 (a, y) # (a Vir)
di
9
(x, (z y
Plus, (zt Vn) is the unique stationary point of (W) in Cps- By lemma 8
ga(2, v) > (x( U/h (x, v) o (2)) Mn)
By Lyapunov asymptotic stability theorem [28]. and applying it to gradient flow for Ip (1, y) in Cus- we can conclude (x)))) (1 VT)
E.10 Proof of Lemma g
Proof. For any (r,y) E Cps in and (z, y) * (z Fii) in Algorithm W.L.O.G, we can assume E (x :27), the analysis details can also be applied to x E (2 a). It is obvious that ă, < 2z+1 and Mi+1 < Uj. Also, lim (Tj, D5) (1%,02) WH). Otherwise either ī; * x or ūj # V/ hold, Algorithm 7 continues until lim (Tj, Mn) limy + (1ne(T)), an (ī,)), i.e. (ūj, Mz) converges to (Irje V())
Moreover, note that for any j = 0, 1....
(.(3j-1.0j-11 > 9p (. Uz) > 9x(
Because
0g ( (0g-1 (z)) where 8 E Ūy1,
(=($))1. (x(ā5-1,00)
Oy
Note that
( Oy
+ gn ( Ø5-1) ga (.(
By the same reason,
9 > 9 (ā Fz))
By Lemma 192 is local minima, and there exists Tp 0 and any (x,,)) I [1 y) (arp l22 Tw Iu (z, > ( wD) Since lim (ā, Dj) = (I w1), there exists J > 0 such that Vij > J. (1 [2) (4 w I Tup combining them all
9x (z. y) 9p (3., 0z) g (1 v()
34

---

Algorithm 7: Path goes to (7))
Input: (x. Cus. Fiu (t), (2) as (38a), (37b)
Output: {āz, My 100
1j=0
(5o, ūo) + (I, y)
for j 1,2,... do
D1 + Mp2 (1)
2j n1 (1--1)
end
E.11 Proof of Lemma 4
Proof. From the proof of Theorem T any any scheduling for P4 satisfies following will do the job
(2/a)2/814
Pk-1 S P < Pk-1
4/3
Note that in Algorithm we have ā
V4(uo + 8) a, then it is obvious
(2/a)2/$112) 4/3 1- < <2//a>/3/ 4/3 M4-1
The same analysis for Theorem can be applied here.
E.12 Proof of Lemma 6
Proof. By the Theorem 3 and Lemma 5 and the fact that A1 is p-stationary point region, we use the same argument as proof of Lemma to demonstrate the gradient descent will never go A2 p.e D
E.13 Proof of Lemma 9
Proof. By Theorem (iv
mAx 2h, I min 21, 120
uSTa
We also know from the proof of Corollary
x+ Hit
< pB
and Fn xhe Consequently,
max *t USTy
p.c S min xt
#>0
Because T3 >>
max n.e S max the I min 14.6
PSTS pS0
E.14
Proof of Corollary
Proof. Note that
a²
4 2 + 1)3
when a >
27 then
4
140 #
1| it satisfies condition in Lemma 4 we obtain the
same result.
E.15
Proof of Corollary
Proof. Use Theorem ( and Theorem ((i)
35

---

E.16 Proof of Corollary 3
Proof. It easy know that
To(un p)) > re(vi)) r((n
and
tØ(x: < (z: t(z:p)
and when /a < т, there are three solutions to r(y: p) by Theorem Also, we know from
Theorem 7
lim (3(
lim (3s μ)) w-0+
2
Note that when
T-B
I+f
S
a +1
a(1 3)
a2 (11-38 -)
ro(vA:m)
(a² +
I
4p >
vp
4p
Therefore,
0 0 r>(V/c))))))))) re(Vrr))))))) r(vu:)
Also, we know that for yub defined in Theorem 5(ii) we know ((Bubi from Theorem 5(v Therefore,
re(Lubip) (yub: p4) ()ubi 0
Besides, v < Yub. By monotonicity of re p)) and (2(3: p)) from the Theorem 7 and Theorem 8(i) it implies that there are at least two solutions (3 μ)) and (3 u). From the geometry of (y p). re (yr rly; μ)) and (z(x; te (x: μ), t( p). it is trivial to know that xp.c S z/.
Vi Vp Tp. E4 Vn,e I y/""
Finally, for every point (z, E A there exists pair c1, 62 each satisfying |c1 S c cndd |eal such that (z, y) is the solution to
I pa + # + y²
pa + E2
x + ula2 + 1)
We can repeat the same analysis above show that x S Vu. 2 Applying the same logic to V(x, u) A we find whe x. Vi y. Thus, (2 v57 the extreme point of A and
(xe)) y*") is the extreme point of A2 uee. we get the results.
F Experiments Details
In this section, we present experiments to validate the global convergence of Algorithm 6 Our
goal is twofold: First, we aim to demonstrate that irrespective of the starting point, Algorithm 6 using gradient descent consistently returns the global minimum. Second, we contrast our updating scheme for & as prescribed in Algorithm 6 with an arbitrary updating scheme for μ4 This comparison illustrates how inappropriate setting of parameters in gradient descent could lead to incorrect solutions.
36

---

F1 Random Initialization Converges to Global Optimum
Algorktim a =2 3--00.00
0.4
033
0.2
-01
-02
-03
0.0
-04
-30
-25
05
2.0
20
-20
-05
60
3.5
(a) Random Initialization
(b) Random Initialization
Algorthmm 22 5.00 6-10 go
32
11
93
(c) Random Initialization
(d) Random Initialization 4
Figure 9: Trajectory of the gradient descent path with the different initializations for 2. We observe that regardless of the initialization, Algorithm 6 alwaay converges to the global minimum.
Initial No a² (1-8)
(1+3)²
37

---

Algottiin d =2.5 gi 3.21 bu1 ug- 20429
gerite efar +0.5 C.0k 6m0.1 6.0425
0.
00
0.2
0.3
-04
-02
0.9
94
05
0.5
0.6
(a) Random Initialization
(b) Random Initialization 2
g- 0429
Aoowithm faria v0.5 0.01 a 1.1 you 3.049
0.-
--2
0.2
06
0.
3.2
14
07
21
06
(c) Random Initialization
(d) Random Initialization
Figure 10: Trajectory the gradient descent path with the different initializations for a = 0.5. We observe that regardless of the initialization, Algorithm Gawaas converges to the global minimum.
Initial Nu 4 (1-8)
(1+3)*
38

---

F22
Wrong Specification of 8 Leads to Spurious Local Optimial
Aigurtirs 6 fora w0. p- 0.01 5-0.4 0.0222
Ago sorithrr 6 fara -0. p.00. 5.0.1. 0,0422
14-
0.3
(a) 0.4
(b) 5 0.1
Figure 11: Trajectory of the gradient descent path for two difference ó. Left: 8 violates requirement
1-8 $ (1- (a² 1) Theorem 4 leading to spurious local minimum. Right: 8 follows
1+8
requirement
1+8
$ (1 () + 1) in Theorem 4 leading to global minimum. Initial No
1-B
a² (1-8)
4
(1+8)²
F3
Wrong Specification of 3 Leads to Incorrect Solution
Algertion 6 Etr +2 5-050- 6. p.000000022
Agoithm 5 for =2 $w0.01 6 0020
24-
64-
0.-
(a) 8 0.6
(b) 8 = 0.01
Figure 12: Trajectory of the gradient descent path for two difference 3. Left: 3 violates requirement
1+
(a + 1) in Theorem 4 leading to incorrect solution. Right: 3 follows requirement
1+8
< (1 5)( 1) in Theorem leading to global minimum. Initial μe 4 (1+D)*
a² (1-))
1-7
39

---

E4
Faster decrease of 14 Leads to Incorrect Solution
Aggriiin for w0. 2-001 0-00.5 u. -0.020)
Algurkten for ei 5 0.00 d 0:15 u -0.0001
44-
02-
0.1
62
00 0.2
0.2
(a) Bad scheduling
(b) Good scheduling
Figure 13: Trajectory of the gradient descent path for two difference update rules for μi with the same initialization. Left: "Bad scheduling" uses faster-decreasing scheme for μ4g, leading to incorrect solution, even a non-local optimal solution. Right: "Good scheduling" follows updating rule for P4 in Algorithm leading the global minimum. Initial (oD
a² (1- (1-6)(118)*
(1+37
40


Reviews:
[{"review_id": "8VABV1Tu6W", "reviewer": "Reviewer_NK3q", "summary": "the paper provides a theoretical study on the loss landscape and convergence in gradient-based DAG learning framework. By focusing on linear functions with the number of variable d = 2, They provide a homotopy-based optimization scheme to guarantee the global optimality. Some numerical validations are provided.", "questions": "- \"...bears a resemblance to the success of training deep models, which started with AlexNet for image classification\": First of all, modern deep models have success before alexnet on imagenet, so this is not accurate. Second, such a statement is over the top and should be removed.\n- L82: indeed only two DAGs are in the entire search space yet the analysis is quite complex. There is no indication on how such analysis and approach would scale or behave with an increasing number of variables, hence it is hard to judge the usefulness of these analysis yet. In comparison, although the combinatorial search is not used, its expected performance is the same with larger d. \n- Eq 4: x is already used to represent data, and should avoid using it as parameter (lack of notation clarity).\n- Eq 5: while re-ordering variable to obtain such a upper triangular structure is not an issue here yet, in datasets W* would generally not be so. How does a non-upper triangular structure impact analysis, for example a's usage afterward?\n- Some derivation are not given fully. For example, f(W) and h(W) skips many steps in between.\n- L199: what is k here?\n- the regime of initialization for \\mu: very interesting observation.\n- optimization landscape with nonlinear functions: current analysis focuses on linear function. How would it change for nonlinear functions (such neural networks), even for d=2? Given the nonlinear nature, one can imagine it would make the analysis even hard, given two layers of neural networks are the most complex cases that has some theoretical understanding. \n- can the homotopy-based optimization scheme be used to improve convergence of NOTEARS related algorithms?", "limitations": "not discussed but not needed.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "the paper focuses on the theoretical study on the global optimality and convergence. This is an important but difficult questions. The theoretical results mainly shows the initialization regime as requirements for convergence, which makes sense given the nonconvex nature, although I did not check the proof in details (given its length).\n\nIllustrative examples and results are provided. ", "weaknesses": "Unfortunately, the current results are only applicable to d=2 and linear functions (understandably, of course). No discussions on how nonlinear or  a large number of variables would affect the loss landscape, and/or how the homotopy algorithm would be affected.\n\nWould be more interesting if the insights of the theoretical results can be used in practical algorithms.  \n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to express our gratitude to the reviewer for their time, effort, and valuable suggestions.\n\n> No discussions on how nonlinear or a large number of variables would affect the loss landscape, and/or how the homotopy algorithm would be affected.\n> \n\nThanks for this insightful question! See **Common Concern I**.\n\n> Would be more interesting … be used in practical algorithms.\n> \n\nIt’s excellent point!  We go deeper into this in the discussion below.\n\n> can the homotopy-based optimization scheme be used to improve convergence of NOTEARS related algorithms?\n> \n\nIndeed, our homotopy-based scheme is proposed precisely because it is what is used in practice [1][4]! Thus, our theory provides guarantee that homotopy methods are indeed a viable approach in practice. Concretely, our results show that tuning the penalty schedule is important, and that global initialization is possible if the signal strength is sufficiently large. When done properly, convergence is very fast (exponential).\n\nIn more detail: When considering both Algorithm 2 and Algorithm 3, any given initialization can converge towards the global optimum, consequently, there's no need to prioritize specific initializations. Furthermore, note that any initial penalty parameter, $\\mu_0$, fall into a specific range can achieve this, negating the need to set $\\mu_0$ to an exceedingly large value to ensure a favorable landscape for $g_{\\mu}(\\Theta)$(i.e., one that only possesses a singular stationary point), as a result, reducing the number of outer iterations in Algorithms 2 and 3. As for the decay rate of $\\mu_k$, it can be quite rapid -  in specific terms, an exponential decay. Our analysis offers clarity on the explicit dependence between convergence rates and various parameters, e.g. Theorem 2 in the paper and Theorem 4 in appendix for details. In principle, this knowledge can be used to accelerate NOTEARS-related algorithms, which is an important direction for future work. To illustrate, a greater weight in $W_*$ can allow for a more progressive decay of $\\mu_k$ and accommodate a smaller initialization for $\\mu_0$, as a result, faster convergence. Empirically, these theoretical conclusions are aligned with our experimental findings (arbitrary initialization, and convergence in a few iterations.) The experiments detailed in [1] confirm this – they begin with a zero initialization point and achieve their objectives in a mere four iterations, utilizing $\\mu$ values set at $\\{1,0.1,0.001,0\\}$, even for hundreds and thousands of nodes.\n\n> \"...bears a resemblance … over the top and should be removed.\n> \n\nThanks for flagging this. in hindsight we see how this is a bit over the top. We will take out this comment in the final version.\n\n> L82: indeed only two DAGs … the same with larger d.\n> \n\nWe appreciate the reviewer's insightful question. See **Common Concern II** for more details.\n\n> Eq 5: while re-ordering … for example a's usage afterward?\n> \n\nThank you for the insightful observation! Absolutely, in a broader context, $W_*$ could be neither an upper nor a lower triangular matrix, posing a considerable challenge for analysis. Determining how this affects our analysis is an interesting future direction. Nonetheless, there is another finding that could be useful. The underlying structure decides the nature of the analysis. For instance, for a three-node system with variables *$X_1, X_2, X_3$,* the analysis for a chain structure like $X_1\\rightarrow X_2\\rightarrow X_3$ should mirror that of $X_2\\rightarrow X_1\\rightarrow X_3$, given their similar structure in essence. Conversely, the analysis would be different if we compare the former chain with a collider structure like $X_1\\rightarrow X_2\\leftarrow X_3$. Yet, when we constrain the dimensions to $d=2$, the two choices are $X_1\\rightarrow X_2$ or $X_2\\rightarrow X_1$, thus, $W_*$ will either be an upper or lower triangular matrix. Consequently, the scenarios become equivalent, suggesting that our analysis for one suffices for both. \n\n> the regime of initialization for $\\mu$: very interesting observation.\n> \n\nThanks for acknowledging this interesting observation. This is based on our analysis of $g_{\\mu}(\\Theta)$, and such regime yields several compelling implications. First, as $a$ increases, the regime typically expands, suggesting the problem get easier when the underlying structure has stronger signal (e.g. $a$ in Equation (5)). Second, any initial penalty parameter, $\\mu_0$, fall into a specific range can achieve this, negating the need to set $\\mu_0$ to an exceedingly large value to ensure a favorable landscape for $g_{\\mu}(\\Theta)$. Finally, it is essential to notice that such regime is the key for us to figure out an universal initialization $\\mu_0 = \\frac{1}{27}$ for Algorithm 3, since $\\frac{a^2}{4(a^2+1)^3}$ is upper bounded by $\\frac{1}{27}$.\n\n> optimization landscape with nonlinear functions… theoretical understanding.\n> \n\nThis is a very good question! When neural networks are introduced into the mix, the loss function $g_{\\mu}(\\Theta) = \\mu f(\\Theta )+h(W(\\Theta))$ inherits a dual-layer complexity. The first layer arises from the inherent nonconvexity of the neural network as seen in the loss $f(\\Theta)$, while the second stems from the nonconvexity of the DAG constraint present in $h(W(\\Theta))$. Such compounded intricacies can be expected to substantially alter the landscape of $g_{\\mu}(\\Theta)$, potentially introducing numerous local optima and saddle points. Although we currently are uncertain on this matter, it's an interesting direction to explore how advancements in the study of two-layer neural networks can be integrated with the techniques we employed, such as the Lyapunov asymptotic stability theorem and the implicit function theorem, to better understand the theoretical facets of DAG learning using neural networks.\n\n[4] Ng, Lachapelle, Ke, Julien, and Zhang. \"On the convergence of continuous constrained optimization for structure learning.\" 2022."}, {"review_id": "K7aRaZfQNE", "reviewer": "Reviewer_xJpx", "summary": "This paper presents a novel approach to the non-convex optimization problems associated with learning the structure of a structural equation model (SEM) or Bayesian network. Considering the equivalent penalty form, the authors propose a homotopy-based optimization scheme that finds global minimizers of the problem by iteratively decreasing the penalty coefficient according to a given schedule. They prove that this algorithm converges globally to the global minimum, regardless of the initialization for W. The authors also demonstrate that the non-convex program is non-benign, meaning that naïve implementation of black-box solvers are likely to get trapped in a bad local minimum.", "questions": "\n\nHow does the landscape look like in high dimensional setting?  Can you provide any experiments to more than two variables?\n\n", "limitations": "Does not apply", "rating": 5, "confidence": 1, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "\nI am not an expert in this field. But the paper's findings seems to have significant implications in the context of learning the structure of SEMs or Bayesian networks.  The paper is well-structured and the authors clearly explain their methodology and findings. They also provide a clear visualization of the non-convex landscape and the solution trajectory.\n\n", "weaknesses": "The authors' approach is primarily focused on the bivariate case, which may limit its applicability in more complex settings. Some numerical study on large-scale problems is desirable. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely thank the reviewer for the insightful critiques and comprehensive understanding of our work, and for providing such useful feedback. We will try our best to address the reviewer’s concern.\n\n> The authors' approach is primarily focused on the bivariate case, which may limit its applicability in more complex settings. Some numerical study on large-scale problems is desirable.\n> \n\nWe appreciate the reviewer's insightful question about the scalability of the Homotopy method on large-scale problems. Indeed, this is exactly what is done in practice on large-scale problems (see discussion in Line 84-87). Since this point is closely related to your follow-up questions, please see below for more details.\n\n> Can you provide any experiments to more than two variables?\n> \n\nMany existing papers show that the Homotopy method can perform well empirically, even with hundreds or thousands of nodes.  For example, the study presented in [1] employs the Homotopy method to solve Equation (1) and has achieved state-of-the-art results, as demonstrated in Figures 4, 5, and 6 in [1], and further elaborated in their appendix. The experiments, encompassing both linear and nonlinear models, offer compelling evidence in support of the homotopy method's effectiveness in practice. Furthermore, while [2] and [3] don't adopt the Homotopy algorithm in its exact form, their approaches share a similar spirit of the homotopy algorithm, solving Equation (2) repeatedly with previous solutions. For more empirical results, see Figure 1 and Appendix H in [2], and Figure 3,7,8 in [3]. More discussions can be found in **Common Concern II**.\n\n> How does the landscape look like in high dimensional setting?\n> \n\nThanks for this insightful question. In higher-dimensional settings, the landscape remains non-benign (e.g. even for $d = 3$, experiments indicate that multiple local optima and saddle points persist), however, existing work [1,3] has shown that the homotopy method is still very effective at finding good (in some cases, near optimal) solutions. This is in fact part of the motivation for our study. Extending our results to higher dimensions is an important direction for future work, and it is worth pointing out that our techniques are indeed generalizable in principle: The main tools we use (the implicit function theorem and the Lyapunov asymptotic stability theorem) possess broad applicability, extending beyond two dimensions. While the direct translation of our findings to other contexts remains challenging, they undeniably pave the way for future explorations. More discussions can be found in **Common Concern I**."}, {"review_id": "bIrup6CSIV", "reviewer": "Reviewer_Y11K", "summary": "This paper studies the problem of learning the correct Directed Acyclic Graph (DAG) that describes the data using continuous optimization. The connection of this problem with continuous methods comes from a prior work of Zheng et al., where they introduce a differentiable function $h$ whose level set at 0 exactly characterizes DAGs. The problem then becomes to optimize a score function of the data subject to the constraint that h is 0. This is a non-convex problem, since h is non-convex. The standard way of solving it is by converting it to an unconstrained problem and penalizing the solutions with high values of h. This paper provides a way of solving a sequence of these problems with varying penalty parameters, so that the final output provably converges to the true solution, when we have two nodes in the DAG.  The model is analyzed in the population case, i.e. when we have access to infinite data.  The authors first show that for sufficiently small values of the penalty parameter, the loss landscape is not benign, as there exists a saddle point and a spurious local minimum in addition to the global minimum. Using this observation, they show a way to pick a sequence of penalty parameters, which provides a solution path that avoids the local minima and converges provably to the true model asymptotically. Their main result is phrased in terms of gradient flow, but they also provide a finite step analysis for gradient descent. ", "questions": "-It would be interesting to see whether similar results hold for higher values of $d$. Have the authors tried to run simulations for $d=3$ or $4$ and observe similar solution landscapes? What about other choices of loss function or penalty?\n-Is it always true that there is a unique DAG that describes the data? For example, for two nodes, suppose we have $X_2 = X_1 + \\epsilon$, where $\\epsilon\\sim \\mathcal{N}(0,1)$. Then, clearly it is also true that $X_1 = X_2 + \\epsilon'$, where $\\epsilon' \\sim \\mathcal{N}(0,1)$. In that case, it seems that we have two solutions $x=1,y=0$ and $x=0,y=1$. I'm curious where would this come up in the analysis. ", "limitations": "Yes", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "Consider the optimization problem \n\\begin{equation}\\label{eq1}\n\\min_\\Theta g_\\mu(\\Theta) := \\mu f(\\Theta) + h(W(\\Theta))\n\\end{equation}\nwhere $f$ is the least squares loss, $\\mu$ is the penalty parameter, $\\Theta$ represents the model parameters and $W(\\Theta)$ represents the adjacency matrix of the DAG. $f$ is a function that was introduced in prior work, which is $0$ if and only if $W(\\Theta)$ corresponds to a DAG. \nThis paper is the first that provides theoretical guarantees about the optimization landscape of this unconstrained problem for various values of $\\mu$, in the case of $2$ nodes. It is worth noting that this is actually the method that is used in practice to learn DAGs, since enumeration methods quickly become intractable when the number of nodes in the DAG increases.\nThe authors establish interesting properties of the solutions space, namely that if $\\mu$ is sufficiently small, the loss is not benign. \nAlso, the idea of successively solving different optimization problems and using the solution of the last as a starting point for the next is an intriguing one and could have applications beyond this work. The proof strategy is also very intuitive and elegant, by choosing $\\mu$ so that the optimizer of $g_\\mu$ lies in the basin of attraction of the global minimum of the next problem. \nAll arguments are very clearly explained, which makes this paper enjoyable to read. ", "weaknesses": "-The result is only proven for $d=2$ nodes in the DAG. This makes most calculations tractable, since there is a closed form for the loss $g_\\mu$ (a quadratic polynomial in two variables). Hence, while the proof strategy is elegant and novel, the technical arguments required to establish the claims do not have significant innovation. \n  ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the time and effort into carefully reviewing our paper and provide such valuable feedback. \n\n> The result is only proven for $d = 2$ nodes in the DAG. This makes most calculations tractable, since there is a closed form for the loss  (a quadratic polynomial in two variables). Hence, while the proof strategy is elegant and novel, the technical arguments required to establish the claims do not have significant innovation.\n> \n\nWe appreciate reviewer for acknowledging the analysis is elegant and novel! See discussion in **Common Concern II**.\n\n> It would be interesting to see whether similar results hold for higher values of $d$. Have the authors tried to run simulations for $d = 3$ or $d = 4$ and observe similar solution landscapes?\n> \n\nThanks for the keen question! Please refer to **Common Concern I** for more details.\n\n> What about other choices of loss function or penalty?\n> \n\nWe appreciate this interesting question! \n\nCrucially, the proof techniques in our work are notably general. Specifically, the implicit function theorem offers a means to trace the trajectory of stationary points, and the Lyapunov asymptotic stability theorem consistently can be used to identify the basin of attraction. These techniques clearly generalize to more general loss functions, which is an intriguing feature of our results. Extending our analysis in this way is an important future direction. Furthermore, it is worth emphasizing that under our model, the least squares loss is the most appropriate loss function: This is because its global minimizer is unique and equals the underlying matrix $W_*$.\n\nIn relation to alternative penalties, we explored the matrix exponential $h_{\\text{expm}}(W) = \\operatorname{tr}(e^{W\\circ W}) - d$  as described in [3] and the log-det formulation $h_{\\text{ldet}}(W) = -\\operatorname{logdet}(sI-W\\circ W) +d\\log s$ where $\\rho(W\\circ W)<s$(with $\\rho$  representing the spectral radius) from [1]. Firstly, all these penalties are designed to enforce the DAG constraint on $W$. They are largely equivalent in application, implying that understanding one offers insights into the others. This equivalence is evident in experimental results; for instance, the loss landscapes are similar to each other. When $\\mu$ falls below a certain threshold, the presence of local optima and saddle points for  $g_{\\mu}(\\Theta)$ remains consistent irrespective of whether $h_{\\text{expm}}(W)$ or $h_{\\text{ldet}}(W)$ is employed. However, incorporating $h_{\\text{expm}}(W)$ or $h_{\\text{ldet}}(W)$ introduces exponential or logarithmic terms to $g_{\\mu}(W)$, adding layers of complexity to the analysis, as well as additional constraints in the case of $h_{\\text{ldet}}$ since its validity rests on the condition $\\rho(W\\circ W)< s$. \n\n> Is it always true that there is a unique DAG that describes the data?\n> \n\nWe appreciate your important question.  This is the identification problem in the DAG learning literature. In essence, the model is identifiable if different DAGs cannot generate the same distribution. In our setting, a linear model with equal (or known) error variances is identifiable (this is well-known, e.g. [5,6]). Therefore, solving Equation (6) will yield a unique DAG that accurately generates the data distribution. Although identifiability is not universal for all models, many have been confirmed as identifiable, e.g. [4,5,6,7,8]. Extending our techniques to these models is an important direction for future work. We hope this clarifies your question.\n\n> For example, for two nodes, suppose we have $X_2 = X_1 +\\epsilon$, where $\\epsilon \\sim N(0,1)$. Then, clearly it is also true that $X_1 = X_2 +\\epsilon'$ , where $\\epsilon'\\sim N(0,1)$. In that case, it seems that we have two solutions $x = 1,y = 0$  and $x = 0,y = 1$. I'm curious where would this come up in the analysis.\n> \n\nThis is a great question! The issue boils down to independence of the noise terms (see L154), which is violated by your example (see below for details). Another way to interpret this assumption is that $N_2$ is independent of $X_1$, or more generally, that each noise term is independent of the parents in the structural equation model. See [5] or [6] for more discussion.\n\nDetails: Your two models are (a)  $X_2 = X_1 +\\epsilon$ and (b) $X_1 = X_2 +\\epsilon'$ with $\\epsilon'=-\\epsilon$. Recalling our model assumptions that $(N_1,N_2)$ are independent and after transcribing the notation, in (a) this is equivalent to $X_1$ independent of $\\epsilon$. Clearly, $X_2$ is not independent of $\\epsilon$. But then in (b), $\\epsilon'$ cannot be independent of $X_2$, which violates our model assumptions. So, there is only one model with independent noise terms. In fact, this argument generalizes quite substantially to LinGAM and additive noise models. In the Gaussian case there are some subtleties, but see [5-7] for details on this case.\n\nThis fact does not explicitly come up in our analysis, but is implicit in the fact that only (a) will minimize the LS loss, a fact which follows from known results such as [5] or [6].\n\n[4] Peters, Jonas, Joris M. Mooij, Dominik Janzing, and Bernhard Schölkopf. \"Causal discovery with continuous additive noise models.\" (2014).\n\n[5] Loh, Po-Ling, and Peter Bühlmann. \"High-dimensional learning of linear causal networks via inverse covariance estimation.\" *The Journal of Machine Learning Research* 15, no. 1 (2014): 3065-3105.\n\n[6] Peters, Jonas, and Peter Bühlmann. \"Identifiability of Gaussian structural equation models with equal error variances.\" *Biometrika* 101, no. 1 (2014): 219-228.\n\n[7] Aragam, Bryon, and Qing Zhou. \"Concave penalized estimation of sparse Gaussian Bayesian networks.\" *The Journal of Machine Learning Research* 16, no. 1 (2015): 2273-2328.\n\n[8] Shimizu, Shohei, Patrik O. Hoyer, Aapo Hyvärinen, Antti Kerminen, and Michael Jordan. \"A linear non-Gaussian acyclic model for causal discovery.\" *Journal of Machine Learning Research* 7, no. 10 (2006)."}, {"review_id": "UIuIl6S2ap", "reviewer": "Reviewer_u4D3", "summary": "The authors give a simple optimization algorithm for DAG-learning-inspired optimization problems that avoids the limitations of known techniques.", "questions": "Discussion around Equation (1) is very clean.\n\nLine 49:\nCan you please explain Equation (2) a bit more?\n\nLine 64:\nIs it really easy to see? :)\n\nLine 90:\nPlease explain \"homotopy.\"\n\nRemark 2:\nCan you please further explain the difficulties that you mention here?\n\nCan you please elaborate on Equation (7)?\n\nLine 217:\nCan you please elaborate on bounding $a$?\n\nPlease add more details in the caption of Figure 2.\n\nLine 235:\nWhy is the interesting regime the one where $\\mu < \\tau$?\n\nLine 258:\nNot clear!\n\nLine 261:\n\"we know\" instead of \"we known.\"", "limitations": "Yes.", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "Originality:\nWork is original.\nI particularly liked the reduction from a combinatorial problem to a non-convex optimization one.\n\nQuality:\nSimple and strong paper.\n\nClarity:\nClear writing, which can be improved; see questions.\n\nSignificance:\nSignificant topic/contributions.", "weaknesses": "No significant weaknesses found.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to express our gratitude to the reviewer for acknowledging the value of our contributions and our clear presentation.\n\n> Line 49: Can you please explain Equation (2) a bit more?\n> \n\nEquation (2) is a penalized version of Equation (1), where $h(W(\\Theta))$ acts as a penalty. This change turns the initial constrained problem into a set of simpler unconstrained problems. The role of $\\mu_k$ in Equation (2) is to control the magnitude of constraint violation by penalty. When $\\mu_k$ is large, $f(\\Theta)$ takes priority, so the solution minimizes this but may not meet the constraint $h(W(\\Theta)) = 0$. By decreasing $\\mu_k$, we give more weight to the penalty. As we approach $\\mu_k = 0$, the solution focuses more on satisfying $h(W(\\Theta))=0$.\n\n> Line 64: Is it really easy to see? :)\n> \n\nWe will be sure to expand on this in the final version. Intuitively, let's recall that $f(\\Theta)$ is presumed to be a convex function in our study. When $\\mu_k$ is large, $f(\\Theta)$ dominates in Equation (2). Consequently, $g_{\\mu_k}(\\Theta)$ behaves similarly to $f(\\Theta)$ as a \"convex\" function, implying a benign loss landscape for $g_{\\mu_k}(\\Theta)$ when $\\mu_k$ is large. This can be formally established as long as $f(\\Theta)$ is a strongly convex function, which our population LS loss is.\n\n> Line 90: Please explain \"homotopy.\"\n> \n\n“homotopy optimization\", also known as \"continuation optimization\" [1, 2], is a strategy for finding solutions to complex optimization problems. The term \"homotopy\" is borrowed from topology (a branch of mathematics), which is a continuous transformation from one function to another. This concept is applied to optimization to form a bridge from a problem with a known or easy-to-find solution to a more complex problem. In our content, resolving $g_{\\mu_k}(\\Theta)$ becomes straightforward with larger values of  $\\mu_k,$ while it becomes challenging when $\\mu_k$ is small.\n\n> Remark 2: Can you please further explain the difficulties that you mention here?\n> \n\nThe main challenge lies in the increasing complexity of the loss landscape of $g_{\\mu}(\\Theta)$ as the model dimensions increase. As a result, analyzing the basin of attraction for the global minimum of $g_{\\mu}(\\Theta)$ becomes increasingly difficult. Moreover, the penalty is an order-$d$ matrix polynomial, which becomes more complex as $d$ increases.\n\n> Can you please elaborate on Equation (7)?\n> \n\nEquation (7) is a detailed version of Equation (2). In the context of our study, we focus on the bivariate case, where the loss function is represented by the expected least square, given by $f(x,y) = \\frac{1}{2}((1-ay)^2+y^2+(a-x)^2+1)$. The constraint function is denoted by $h(x,y) = \\frac{x^2y^2}{2}.$ \n\nTherefore, we can express $g_{\\mu}(x,y) = \\mu f(x,y)+h(x,y)=\\frac{\\mu}{2}((1-ay)^2+y^2+(a-x)^2+1)+\\frac{x^2y^2}{2}$. The function $g_{\\mu}(x,y)$ is used recurrently in our Algorithm 2 or 3, so our principal interest lies in understanding its properties.\n\n> Line 217: Can you please elaborate on bounding $a$?\n> \n\nThanks for bringing it up! The lower bound on $a$ is indeed a standard and sensible assumption. One can view the magnitude of $a$ as a measure of problem difficulty, with $a$ essentially behaving like the “signal strength\" of the underlying structure, larger $a$ indicates the structure is easier to learn. For a more direct insight, consider Equation (6)'s local optimal solution, namely $(x_1^*,y_1^*) = (0,\\frac{a}{a^2+1})$. The corresponding loss is  $\\frac{1}{2}(a^2+1+\\frac{1}{a^2+1})$. Now, the global optimal solution is $(x_0^*,y_0^*) = (a,0)$, with loss $f(x^*_0,y^*_0) = 1$. The loss difference between the local and global optima is therefore $\\frac{1}{2}(a^2+1+\\frac{1}{a^2+1}) - 1$, an increasing function of $a$. As $a$ enlarges, both the loss difference and the basin of attraction for the global optimum increase in general, making the optimization easier.\n\n> Please add more details in the caption of Figure 2.\n> \n\nThanks for the suggestion. Here, for $\\mu>\\tau,$ there exists a single solution to $r(y;\\mu) = 0$, which implies there is one stationary point in Equation (7). When $\\mu=\\tau,$ two solutions are found for $r(y;\\mu) = 0$, suggesting that there are two stationary points in Equation (7). Conversely, when $\\mu<\\tau,$ we observe three solutions for $r(y;\\mu) = 0$, indicating that there are three stationary points in Equation (7)-  a local optimum, a saddle point, and a global optimum.\n\n> Line 235: Why is the interesting regime the one where $\\mu<\\tau$?\n> \n\nThanks for bringing up this point! In this regime where multiple stationary points exist in $g_{\\mu}(\\Theta)$ - a local optimum, a saddle point, and a global optimum, as described in Lemma 1, the usage of gradient descent for $g_{\\mu}(\\Theta)$ may result in trap within the local optimal solution. This scenario is precisely what our methodology aims to avoid, our updating scheme of $\\mu_k$ in Algorithm 2 or 3 can do the job. Conversely, when $\\mu>\\tau,$ only a single stationary point is observed. As demonstrated in Theorem 3 (Appendix), gradient descent invariably converges to this point, which consistently represents the global optimum. Hence, in such conditions, there are no concerns to be addressed.\n\n> Line 258: Not clear!\n> \n\nWe are happy to explain this in a more accessible way. The fundamental idea is to identify the basin of attraction for the global optimum solution$(x^*_{\\mu_{k+1}},y^*_{\\mu_{k+1}})$ of $g_{\\mu_{k+1}}(\\Theta)$, and make sure previous initialization point $(x^*_{\\mu_{k}},y^*_{\\mu_{k}})$ fall into such region. We then use gradient descent or flow to converge towards $(x^*_{\\mu_{k+1}},y^*_{\\mu_{k+1}})$. This proof aligns with this reasoning and affirms that our updating scheme for $\\mu_k$ can make it happen!\n\n[1] Lin, Yang, Zhang, and Zhang. \"Continuation Path Learning for Homotopy Optimization.\" 2023.\n\n[2] Hazan, Levy, and Shalev-Shwartz. \"On graduated optimization for stochastic non-convex problems.\" 2016."}]

Author rebuttal:
To all reviewers,

We thank all reviewers for their time put into reading our work and their valuable comments. We appreciate the consensus that our paper is well-written and theoretical contributions are delivered clearly. Finally, we appreciate Reviewer Y11K acknowledging our analysis elegant and novel. We next respond to the common concerns.

*Due to character limits, any minor concerns omitted will certainly be addressed based on the reviewers' suggestions in the final version.*

*Regarding the other major concerns raised, we will respond to each reviewer individually.* 


**Common Concern I**

> How landscape is affected for high dimensional setting and how the homotopy algorithm would be affected.
> 

This is a very insightful question. In higher-dimensional settings, the landscape remains challenging, however, previous work has convincingly demonstrated the utility and applicability of homotopy-based approaches in practice. Indeed, the original NOTEARS implementation as well as more recent advances use some form of homotopy algorithm. See [1,2,3] for details. Although homotopy methods generally faces challenges with getting trapped in local optima, there are certain scenarios where this can be potentially circumvented. Based on our own experiments, global convergence remains feasible in these instances:

- Increasing the signal strength (e.g. $a$ in Equation (5)) appears to prevent the Homotopy method from being trapped into local optima. This weight can be seen as the signal strength of the underlying structure, and enhancing it could also alleviate the "trapping into local optimum" issue. This aligns with our findings in the case of $d = 2$; when $a > \sqrt{5/27}$, there exists a universal updating scheme for $\mu_k$ that functions effectively, and global convergence is guaranteed as stated in Corollary 1.
- Enhancing the sparsity of the underlying structure indeed seems to mitigate the risk of the Homotopy method getting trapped into local optima. This intriguing observation coincides with prior experimental results in [1][2][3], which indicate that sparse graphs are generally easier to learn. As such, this important discovery further reinforces the value of investigating sparsity as a potential direction for future work.

In higher-dimensional settings, the landscape remains non-benign. For instance, even when $d = 3$, experiments indicate the presence of multiple local optima and saddle points. Determining how to incorporate sparsity structure and signal strength into the analysis to ensure global convergence of the homotopy method is an important direction for future work.

**Common Concern II**

> How such analysis and approach would scale or behave with an increasing number of variables?
> 

This is really insightful point!

- Firstly, many existing papers suggest that the Homotopy method can perform well empirically, even with hundreds or thousands of nodes. This is described in our paper (specifically, Remark 1, L84-87). The study presented in [1] employs the Homotopy method to solve Equation (1) and has achieved state-of-the-art results, as demonstrated in Figures 4, 5, and 6 in [1], and further elaborated in their appendix. Furthermore, while [2] and [3] don't adopt the Homotopy algorithm in its exact form, their approaches share the similar spirit of homotopy algorithm, solving Equation (2) repeatedly with previous solutions. It is essential to note, however, that a theoretical understanding of the convergence dynamics and the loss landscape of $g_{\mu}(\Theta)$ remains largely unexplored. Our investigation into the bivariate scenario aspires to address this gap, offering a theoretical framework to explain and improve these known empirical findings.
- Secondly, we would like to emphasize that the current analysis is intuitive but far from trivial and has the potential to inform more complex settings. Initially, we delineate the conditions leading to multiple stationary points of $g_{\mu}(\Theta)$ as evidenced by Theorems 5, 6(v) – a process that relies heavily on a clever analysis of two related objectives $r(y_{\text{ub}};\mu)$ and $t(x_{\text{lb}};\mu)$ , especially given the absence of explicit formula for the solutions to $r(y;\mu) = 0$  and $t(x;\mu) = 0$ (see proof of Theorems 5,6 in appendix). Subsequently, we leverage the implicit function theorem to identify the evolving trajectories of these stationary points as a function of $\mu$ (see Theorem 5, 6(vi)), pinpointing the precise trajectory to pursue. Finally, our use of the Lyapunov asymptotic stability theorem allows us to define the basin of attraction for the correct stationary point, detailed in the proof of Lemma 7. Each of these components not only presents an intuitive understanding but can also serve as foundational blocks for more intricate cases. It's noteworthy that tools like the implicit function theorem and the Lyapunov asymptotic stability theorem possess broad applicability, extending beyond the confines of linear models. While the direct translation of our findings to other contexts remains challenging, they undeniably pave the way for future explorations.

[1] Bello, Kevin, Bryon Aragam, and Pradeep Ravikumar. "Dagma: Learning dags via m-matrices and a log-determinant acyclicity characterization." *Advances in Neural Information Processing Systems* 35 (2022): 8226-8239.

[2] Ng, Ignavier, AmirEmad Ghassami, and Kun Zhang. "On the role of sparsity and dag constraints for learning linear dags." *Advances in Neural Information Processing Systems* 33 (2020): 17943-17954.

[3] Zheng, Xun, Bryon Aragam, Pradeep K. Ravikumar, and Eric P. Xing. "Dags with no tears: Continuous optimization for structure learning." *Advances in neural information processing systems* 31 (2018).

Meta-review: This work addresses the DAG learning problem. Despite being a nonconvex optimization problem with the possibility of multiple spurious solutions, the authors claim that their proposed algorithm can find the global optimal solution for the population loss but only in the bivariate setting. However, it's important to note that the theoretical result is limited to cases where $d=2$, which is not practical in most real-world scenarios.

All reviewers have generally provided positive feedback. While the work is interesting and may have some scientific impact, it currently lacks practical usability.
Decision: Accept (poster)