arXiv CSopen access

Second-Order KKT Guarantees for Bregman ADMM in Nonconvex and Non-Lipschitz Optimization

Second-Order KKT Guarantees for Bregman ADMM in Nonconvex and Non-Lipschitz Optimization Shuang Li, Zhihui Zhu, and Qiuwei Li*

arXiv:2606.28307v1 [math.OC] 26 Jun 2026

June 29, 2026

Abstract We analyze Bregman ADMM for nonconvex linearly constrained problems under two-sided relative smoothness, a condition that replaces the standard Lipschitz gradient assumption with a Hessian comparison relative to a Bregman kernel. This setting covers polynomial objectives arising in matrix and tensor models for which a global Lipschitzgradient constant need not exist. We show that on an invariant open state-space domain, one iteration of Bregman ADMM defines a smooth primal–dual fixed-point map whose strict-saddle KKT points are unstable fixed points; consequently, from random initialization the iterates converge to a strict saddle with probability zero. Combined with existing first-order convergence results, this yields almost-sure second-order stationarity of limiting KKT points. We extend the analysis to a multi-block star consensus formulation for distributed optimization. The technical novelty lies in a determinant reduction with a Bregman-specific symmetrization and scaling step in the two block spectral argument, together with a null space cancellation exploiting the star graph structure in the consensus case. Numerical experiments on distributed matrix factorization illustrate the theory, and a symmetric tensor factorization example demonstrates the broader Bregman proximal splitting idea beyond the separable consensus setting.

1

Introduction

We consider the following optimization problem . minimize f (x, y) = f1 (x) + f2 (y) s.t. Ax + By = b,

x∈Rn ,y∈Rm

(1)

where A ∈ Rp×n , B ∈ Rp×m , and b ∈ Rp . The functions f1 : Rn → R and f2 : Rm → R can be nonconvex and their gradients need not be globally Lipschitz continuous. The augmented Lagrangian is ρ L(x, y, λ) =f1 (x) + f2 (y) + ⟨λ, Ax + By − b⟩ + ∥Ax + By − b∥22 , (2) 2 where λ ∈ Rp is the Lagrangian dual variable and ρ > 0 is the augmented penalty parameter. Throughout the paper, we treat ρ as a fixed constant (i.e., it does not vary with the iteration index). A standard approach for (1) is the alternating direction method of multipliers (ADMM; see Algorithm 1) [1, 2]. Each iteration alternates between minimizing over x and over y, followed by a dual ascent step for λ. Algorithm 1 ADMM 0

1: Initialization: (y0 , λ ) k 2: For k = 1, 2, . . . , generate (xk , yk , λ ) by

xk = argmin L(x, yk−1 , λk−1 ), x

y = argmin L(xk , y, λk−1 ), k

λk =λ

y k−1

(3)

+ ρ(Axk + Byk − b).

* SL is with the Department of Electrical and Computer Engineering, Iowa State University, Ames, Iowa 50014, USA. ZZ is with the Department of Computer Science and Engineering, Ohio State University, Columbus, Ohio 43210, USA. QL is with DAMO Academy, Alibaba Group US, Bellevue, Washington 98004, USA.

1

Common variants include proximal ADMM [3, 4], which adds proximal regularization in the primal updates, and linearized ADMM [5], which replaces the primal minimizations by alternating gradient-type steps. These splitting schemes are well suited to structured models where separability can be exploited, with applications spanning low-rank matrix recovery [6], tensor factorization, consensus optimization, and beyond; see [7] for an overview. First-order convergence and complexity for nonconvex ADMM. For nonconvex linearly constrained problems and related consensus formulations, convergence of ADMM to first order stationary points has been studied under various regularity conditions. [8] establishes convergence for certain nonconvex consensus and sharing formulations under regularity assumptions on the objective and subproblems and sufficiently large augmented Lagrangian penalty parameters, while [9] analyzes a Bregman modification of ADMM for nonconvex composite problems and proves convergence to stationary points of the associated augmented Lagrangian under suitable assumptions. [10] further proves global convergence of ADMM for nonconvex nonsmooth problems by combining conditions that ensure lower boundedness, sufficient descent and a subgradient bound, together with the Kurdyka Łojasiewicz (KŁ) inequality. More recently, [11] develops finite iteration convergence guarantees under restricted strong convexity without requiring smoothness or differentiability, and [12] proposes communication efficient decentralized primal dual methods for smooth nonconvex consensus optimization. Second-order stationarity and strict-saddle avoidance. Moving beyond first-order stationarity is useful because many structured nonconvex problems have benign landscapes in which all second-order stationary points are globally optimal; see, e.g., low-rank matrix recovery [13] and tensor decomposition [14]. For ADMM-type methods, [5] uses a dynamical-systems viewpoint and the stable manifold theorem [15, Theorem III.7] to show that linearized ADMM avoids strict saddles almost surely under a Lipschitz-gradient assumption. Beyond ADMM, strict-saddle avoidance has been extended to broader classes of first-order and proximal-type algorithms. In particular, [16] shows that a stable-manifold argument continues to apply under vanishing stepsizes (covering, among others, mirror-descenttype updates), and [17] establishes that proximal methods avoid (active) strict saddles for weakly convex nonsmooth objectives. We note that [17] works in the nonsmooth weakly convex regime and introduces the notion of active strict saddle, which accounts for the combinatorial structure of the subdifferential; their setting and ours are largely complementary: we treat smooth (but non-Lipschitz-smooth) objectives within ADMM splitting, whereas they handle nonsmooth objectives via proximal point iterations. Bregman geometry and relative smoothness. A key limitation of Lipschitz-gradient analyses is that ∥∇2 f ∥ may be unbounded in many applications of interest (e.g., polynomial objectives arising in matrix/tensor models), so the standard smoothness constant does not exist globally. Relative smoothness and Bregman first-order methods [18, 19] provide a principled alternative by replacing Euclidean geometry with a problem-adapted distance-generating function h (see also the Bregman divergence [20, 21]). Bregman-divergence-based ADMM methods have been studied for firstorder convergence, e.g., [9, Theorem III.8] and [22, Theorem 3.9]; we discuss specific Bregman ADMM variants next. Bregman ADMM variants. Early Bregman variants of ADMM were proposed to better match the problem geometry and to simplify subproblems by replacing Euclidean proximal terms with Bregman divergences; see, e.g., the Bregman ADMM of [23]. More recent work has developed first-order convergence theory in nonconvex and multi-block settings. For instance, [24] studies an improved Bregman ADMM with convergence-rate guarantees, [25] considers a Bregman proximal-linearized scheme for separable sums coupled by a difference of functions term, and [26] analyzes a distributed multi-block variant for nonconvex and nonsmooth sharing problems; see also [22, Theorem 3.9] for a representative first-order result. What is missing from this line of work is a mechanism excluding convergence to strict saddles. Our contribution fills this gap: under relative (bi-)smoothness, strict saddles are unstable fixed points of the induced smooth map, and random initialization therefore rules out convergence to them almost surely. Combined with standard first-order convergence assumptions, this yields second-order stationarity of limiting KKT points. At the same time, the extension from the Euclidean and Lipschitz setting of [5] to Bregman ADMM is not a trivial substitution of proximal terms. The Euclidean analysis in [5] also reduces its instability proof, after a suitable transformation, to a real symmetric spectral argument. In the Bregman setting, however, replacing the Euclidean proximal term by a Bregman divergence introduces metric blocks from the kernel Hessians, so the determinant reduction and sign change argument must be rederived in the Bregman geometry. In the two block case, this requires a diagonal

2

rescaling that symmetrizes the reduced spectral problem and a correspondingly rescaled strict saddle direction. In the consensus case, the proof is again not a direct lift of the Euclidean argument: the star consensus structure creates a hub and peripheral coupling pattern, and the key sign change relies on a null space cancellation of the consensus penalty, since the equal x test direction lies in the kernel of the star graph Laplacian. These are the main new technical ingredients behind the strict saddle result in the Bregman setting. Contributions. Our main contributions are as follows. • Strict-saddle avoidance for Bregman ADMM beyond Lipschitz smoothness. We analyze Bregman ADMM (Algorithm 2) under relative (bi-)smoothness with respect to suitable Bregman kernels. By viewing one iteration as a smooth primal–dual fixed-point map and showing that strict saddles are unstable fixed points, we prove that with random initialization the probability that the iterates converge to a strict saddle is zero, without requiring a Lipschitz gradient assumption. The two block proof requires a determinant reduction, symmetrization, and scaling step specific to the Bregman setting in the spectral argument; see the discussion around Section 5. • Second-order stationarity of limit points as a corollary of strict-saddle avoidance. Under additional assumptions guaranteeing first-order convergence to a KKT point (e.g., the boundedness and KŁ and subanalyticity conditions in [9, Theorem III.8] and [22, Theorem 3.9]), our instability result implies that limiting KKT points are almost surely second-order stationary. In the Euclidean special case when the two Bregman kernels are chosen as h1 (x) = 21 ∥x∥22 and h2 (y) = 12 ∥y∥22 , Algorithm 2 reduces to proximal ADMM. A recent work of [27] establishes strict saddle avoidance and almost sure second order convergence for Euclidean proximal ADMM. Relative to this Euclidean proximal ADMM result, our contribution is to extend the strict saddle avoidance mechanism to the Bregman setting under relative smoothness, which covers polynomial objectives without a global Lipschitz constant, and to develop the corresponding multi block consensus extension. • Multi-block consensus extensions with the same second-order stationarity guarantee. We extend the framework to a global consensus formulation and prove the same strict-saddle avoidance guarantees for Consensus Bregman ADMM; see Section 4. Here the key new ingredient is a star-consensus structural reduction, together with a hub/nonhub null-space argument that isolates the negative-curvature direction despite the consensus penalty. • Practical implications for polynomial objectives and distributed factorization problems. Using known relative(bi-)smoothness constructions for polynomial objectives, our theory applies to a broad class of non-Lipschitz models arising in matrix/tensor optimization. We report numerical experiments on distributed matrix factorization and symmetric tensor factorization in Section 6. The paper is organized as follows. Section 2 reviews stationarity notions and tools from Bregman geometry and relative smoothness. Section 3 presents Bregman ADMM and the strict-saddle avoidance result. Section 4 extends the algorithm and guarantees to a global consensus setting. Section 5 first presents the single-block core mechanism through the proof for the Bregman augmented Lagrangian method (Bregman ALM), and then gives proof sketches explaining how the argument extends to the two-block Bregman ADMM and consensus settings. Section 6 reports numerical experiments. The appendix contains the full two-block Bregman ADMM and consensus proofs. Notation. We write bold lowercase letters (e.g., x, y, z) for vectors, bold uppercase letters (e.g., A, B, H) for matrices, and T for tensors. Rn denotes n-dimensional Euclidean space. ∥ · ∥2 is the Euclidean (vector) norm and ∥ · ∥F the Frobenius (matrix/tensor) norm. ⟨·, ·⟩ denotes the standard inner product. null(M) denotes the null space of M, and I denotes an identity matrix of compatible dimension. ∇f and ∇2 f denote the gradient and Hessian of f ; for functions of two variables we write ∇2xx f, ∇2yy f for the diagonal Hessian blocks and ∇2xy f, ∇2yx f for the cross Hessian blocks. For a matrix M, λi (M) denotes an eigenvalue indexed by i; when M is real symmetric, λmin (M) denotes its smallest eigenvalue. M ≻ 0 (resp. M ⪰ 0) means M is positive definite (resp. positive semidefinite). For matrices M1 , . . . , Mk , diag(M1 , . . . , Mk ) denotes the block-diagonal matrix with diagonal blocks M1 , . . . , Mk . Dg denotes the Jacobian matrix (total derivative) of a differentiable map g, and g k its k-fold composition; ∇x f denotes the partial . . gradient with respect to x. We write [J] = {1, . . . , J} for the index set, and use = throughout to mean “is defined as”.

2

Preliminaries

Definition 2.1. Consider the following equality-constrained program with A ∈ Rp×n , b ∈ Rp , and twice continuously differentiable f :

3

minimize f (x) s.t. Ax = b. n x∈R

1. x is a (first-order) stationary point if Ax = b and ∇f (x) + A⊤ λ = 0 for some multiplier λ ∈ Rp ; 2. x is a second-order stationary point if it is a first-order stationary point and d⊤ ∇2 f (x)d ≥ 0 for all feasible directions d ∈ null(A); 3. x is a strict saddle if it is a first-order stationary point but not a second-order stationary point, i.e., d⊤ ∇2 f (x)d < 0 for some d ∈ null(A). Definition 2.2 (Multi-block stationarity). The notions in Definition 2.1 extend to any linearly constrained problem minz F (z) subject to Cz = b with twice continuously differentiable F by viewing all decision variables as a single block z and replacing A with C. In particular, for Problem (1), set z = (x, y), F (z) = f1 (x) + f2 (y), and C = [A B]. 1. (x⋆ , y⋆ ) is a first order stationary point if and only if Ax⋆ + By⋆ = b, ∇f1 (x⋆ ) + A⊤ λ⋆ = 0, and ∇f2 (y⋆ ) + B⊤ λ⋆ = 0 for some λ⋆ . 2 ⋆ 2. (x⋆ , y⋆ ) is a second order stationary point if and only if it is a first order stationary point and d⊤ x ∇ f1 (x )dx + ⊤ 2 ⋆ dy ∇ f2 (y )dy ≥ 0 for all (dx , dy ) satisfying Adx + Bdy = 0. ⋆ 2 ⋆ ⊤ 2 3. (x⋆ , y⋆ ) is a strict saddle if and only if it is a first order stationary point and d⊤ x ∇ f1 (x )dx +dy ∇ f2 (y )dy < 0 for some (dx , dy ) satisfying Adx + Bdy = 0; cf. (46) and (47) in Section A.

The same convention applies to the consensus formulation after stacking all local variables and writing the agreement constraints as a linear system. We refer to any primal dual point satisfying the feasibility and first order conditions in the relevant stacked formulation as a KKT point. Bregman divergence and relative smoothness. kernel h [20, 21].

We recall the Bregman divergence induced by a strongly convex

Definition 2.3. The Bregman divergence induced by a differentiable, strongly convex function h is defined as Dh (x, y) = h(x) − h(y) − ⟨∇h(y), x − y⟩

(4)

Geometrically, Dh (x, y) is the gap between h(x) and the first-order Taylor approximation of h at y, so it measures discrepancy in the geometry induced by the kernel rather than by the Euclidean norm. When h(x) = 21 ∥x∥22 , the Bregman divergence reduces to the standard squared Euclidean distance Dh (x, y) = 12 ∥x − y∥22 . More generally, Bregman divergences underpin mirror-descent and Bregman-proximal schemes: [18] introduces the descent lemma under relative smoothness, [19] extends the framework to nonconvex objectives, and [28] constructs polynomial kernels ensuring relative smoothness for broad classes of problems. These ideas also lead to Bregman variants of ADMM [9, 22, 23]. One key analytic advantage is that the smoothness of f can be measured relative to h. Examples used later.

We record several kernels and the corresponding Bregman divergences used below.

• Quadratic kernel (Euclidean). If h(x) = 12 ∥x∥22 , then Dh (x, y) = 12 ∥x − y∥22 . • Polynomial kernels used for relative smoothness. For the kernel in Proposition 2.1, where d ≥ 2 is the polynomial degree, 1 1 h(x) = ∥x∥d2 + ∥x∥22 + 1, ∇h(x) = ∥x∥d−2 x + x, 2 d 2 so Dh (x, y) can be computed by substituting these expressions into (4).

4

• Bi-variable/product kernels and bi-Bregman divergence. For the product-type bi-convex kernel used in Proposition 2.2, h(x, y) = hx (x)hy (y), the bi-Bregman divergences simplify to scalar multiples of the corresponding one-variable Bregman divergences: Dhx (x1 , x2 ; y) = hy (y) Dhx (x1 , x2 ), Dhy (y1 , y2 ; x) = hx (x) Dhy (y1 , y2 ). This structure is used in the multi-block extensions of Section 4. Definition 2.4. A C 2 function f (x) is relatively smooth with constant Lf with respect to a C 2 strongly convex function h(x) if Lf ∇2 h(x) ± ∇2 f (x) ⪰ 0,

∀x ∈ Rn .

(5) 2

2

Our convention in (5) is a two-sided Hessian comparison: besides the standard upper bound ∇ f (x) ⪯ Lf ∇ h(x), we also assume the lower bound ∇2 f (x) ⪰ −Lf ∇2 h(x). The lower bound is what guarantees strong convexity of the Bregman-proximal subproblems used later in the paper. When h(x) = 21 ∥x∥22 , (5) reduces to the standard Lipschitzsmoothness condition ∥∇2 f (x)∥ ≤ Lf . Under suitable standard assumptions in addition to relative smoothness, Bregman gradient and proximal gradient methods admit descent and convergence guarantees, including convergence to first order stationary points in the nonconvex setting [18,19], thereby avoiding the need for a global Lipschitz constant in the gradient. This framework is particularly useful for polynomial models, which are central to the applications below. In particular, any finite degree polynomial objective is relatively smooth under a suitable polynomial kernel h [28]; the corresponding two variable construction yields relative bi-smoothness, as stated next. Proposition 2.1 ( [28]). We call f (x) a dth-degree polynomial if the maximum total degree among its monomials of the form xa1 1 · · · xann is d. Let f (x) be a dth-degree polynomial with d ≥ 2 and h(x) = d1 ∥x∥d2 + 12 ∥x∥22 + 1. Then f is Lf -relatively smooth to h for some Lf > 0 depending on the coefficients of f . Remark 2.1 (Two-sided condition). The relative smoothness convention in [28] is already the two sided Hessian comparison Lf ∇2 h(x) ± ∇2 f (x) ⪰ 0. Thus Proposition 2.1 is a specialization of the polynomial kernel construction in [28]. If one starts instead from the standard one sided upper bound convention, the lower bound can also be obtained by applying the same polynomial kernel argument to −f and taking Lf ← max(Lf , L−f ). The same reasoning applies to the bi variable case (Proposition 2.2). Since (1) involves two primal blocks, we also use the following bi-smoothness conditions. Definition 2.5. A twice continuously differentiable kernel h(x, y) is called strongly bi-convex on its domain if there exists mh > 0 such that ∇2xx h(x, y) ⪰ mh I, ∇2yy h(x, y) ⪰ mh I, ∀x, y. Definition 2.6. The bi-Bregman divergence induced by a strongly bi-convex function h is defined as Dhx (x1 , x2 ; y) =h(x1 , y) − h(x2 , y) − ⟨∇x h(x2 , y), x1 − x2 ⟩,

(6)

Dhy (y1 , y2 ; x) =h(x, y1 ) − h(x, y2 ) − ⟨∇y h(x, y2 ), y1 − y2 ⟩. (7) y x 2 2 2 2 When h(x, y) = (∥x∥2 + ∥y∥2 )/2, we have Dh (x1 , x2 ; y) = ∥x1 − x2 ∥2 /2 and Dh (y1 , y2 ; x) = ∥y1 − y2 ∥2 /2, respectively. Definition 2.7. A twice continuously differentiable function f (x, y) is relatively bi-smooth with constants (Lxf , Lyf ) with respect to a strongly bi-convex function h(x, y) if Lxf ∇2xx h(x, y) ± ∇2xx f (x, y) ⪰ 0, Lyf ∇2yy h(x, y) ± ∇2yy f (x, y) ⪰ 0, ∀x, y.

(8)

As in the one-variable case, (8) is a two-sided assumption: it controls each diagonal Hessian block of f both above and below by the corresponding diagonal Hessian block of h. Throughout the paper, “relative smoothness” and “relative bi-smoothness” refer to these two-sided conditions unless stated otherwise. Proposition 2.2 ( [28]). We call f (x, y) a (d1 , d2 )th-degree polynomial if it is a d1 th-degree polynomial in x and a d2 th-degree polynomial in y. Let f (x, y) be a (d1 , d2 )th-degree polynomial with d1 , d2 ≥ 2 and h(x, y) = ( d11 ∥x∥d21 + y y d2 1 1 1 2 2 x x 2 ∥x∥2 + 1)( d2 ∥y∥2 + 2 ∥y∥2 + 1). Then f is (Lf , Lf )-relatively bi-smooth with respect to h for some Lf , Lf > 0. 5

3

Bregman ADMM

We consider Bregman ADMM for the two-block linearly constrained problem (1). Under the one variable relative smoothness condition of Definition 2.4, applied separately to each block, together with the kernel strong convexity, stepsize, and invariant domain assumptions stated in Theorem 3.1, the probability that the iterates converge to a strict saddle is zero. Combined with first order convergence guarantees (e.g., [9, Theorem III.8]), this yields almost sure convergence to second order stationary points. Bregman ADMM alternates between two Bregman-proximal primal updates and a dual ascent step. The divergences Dh1 and Dh2 act as geometry-adapted regularizers that stabilize each block subproblem and allow the analysis to rest on the relative smoothness condition of Definition 2.4 rather than global Lipschitz continuity of ∇f . When h1 (x) = 12 ∥x∥22 and h2 (y) = 21 ∥y∥22 , Bregman ADMM reduces to proximal ADMM. The algorithm is stated in Algorithm 2. Algorithm 2 Bregman ADMM 1: Input: Bregman divergence kernels h1 , h2 and parameters η, ρ. 0 2: Initialization: (x0 , y0 , λ ) k 3: For k = 1, 2, . . . , recursively generate (xk , yk , λ ) by 1 xk = argmin L(x, yk−1 , λk−1 ) + Dh1 (x, xk−1 ), η x 1 k y = argmin L(xk , y, λk−1 ) + Dh2 (y, yk−1 ), η y

(9)

λk =λk−1 + ρ(Axk + Byk − b).

First-order convergence of Bregman ADMM has been analyzed for the two-block scheme in [9, Theorem III.8] and for the multi-block scheme in [22, Theorem 3.9]. These works P use Bregman distance regularization, with constraints written as Ax = By, Ax + By = 0, or more generally i Ai xi = 0. Up to block order, sign, and affine shift conventions, our updates coincide with these Bregman ADMM forms by using the scaled Bregman regularizers η1 Dh1 and η1 Dh2 and replacing the homogeneous constraint by Ax+By = b. Under Assumption 3.1 below, we complement these results by showing that the iteration defines a smooth fixed-point map g whose strict saddles are unstable fixed points; consequently, the set of initializations that converge to strict saddles has Lebesgue measure zero. Here and throughout, random initialization means that the initial algorithmic state, for instance (x0 , y0 , λ0 ) for Algorithm 2, is drawn from a distribution on an open state space domain Ω that is absolutely continuous with respect to Lebesgue measure. Assumption 3.1 (Open-domain self-map condition). For the algorithm under consideration, the objective and Bregman kernels are C 2 on an open primal domain, and the induced fixed-point map g is well defined on an open state-space domain Ω satisfying g(Ω) ⊆ Ω. Remark 3.1 (Sufficient conditions for Assumption 3.1). Assumption 3.1 is the technical condition needed for the implicit function theorem and the dynamical systems argument used below to apply on the whole state space. Two common sufficient scenarios are: • Global full domain case. If fj , hj ∈ C 2 on the whole primal space and the Bregman proximal subproblems are globally well defined, for instance under the strong convexity, two sided relative smoothness, and stepsize assumptions of Theorem 3.1, then one may take Ω to be the full primal dual state space; the self map and domain invariance parts of Assumption 3.1 are automatic. • Invariant open domains. In problems with natural open primal domains (e.g., when the kernels are Legendre-type so the primal variables stay in the interior of dom hj ), the Bregman-proximal subproblems typically return minimizers in the same interior; taking Ω to be the corresponding product with the unrestricted dual space then yields g(Ω) ⊆ Ω. If one also wants to conclude almost sure convergence to second order stationary points from first order convergence, then one sufficient route is to assume boundedness of the generated sequence, possibly obtained from coercivity in 6

applications, together with KŁ and subanalyticity type assumptions and the other regularity assumptions required by first order Bregman ADMM convergence results such as [9, Theorem III.8]. Theorem 3.1. Assume Assumption 3.1 holds for the fixed-point map induced by Algorithm 2, and assume also that each kernel hj is strongly convex. Let fj be Lj -relatively smooth w.r.t. hj in the sense of Definition 2.4 (applied to the j-th block), j = 1, 2, and set η < minj 1/Lj . With random initialization, the probability that the iterates of Algorithm 2 converge to a strict saddle is zero. We stress that Theorem 3.1 does not assert convergence of the iterates; it states that the set of initializations leading to convergence to a strict saddle has Lebesgue measure zero. Convergence itself requires separate conditions; once convergence to a KKT point, i.e., a primal-dual point satisfying the feasibility and first-order conditions in Definition 2.2, is established, the theorem excludes strict saddles from the set of possible limits. The proof of Theorem 3.1 is given in Section A. We emphasize that C 2 regularity alone does not guarantee either invariance of the state-space domain or boundedness of the iterates; Remark 3.1 gives representative sufficient conditions. Under additional assumptions guaranteeing first order convergence, for example assumptions analogous to [9, Theorem III.8], which include boundedness of the generated sequence and the regularity and KŁ assumptions needed for convergence to a stationary point of the augmented Lagrangian, Theorem 3.1 implies that the limiting KKT points are almost surely second order stationary. By Proposition 2.1, the strict saddle avoidance conclusion of Theorem 3.1 applies to finite degree polynomial objectives f1 and f2 with corresponding polynomial kernels. If, in addition, a first order convergence guarantee is available for the resulting Bregman ADMM iterates, then any limiting KKT point is almost surely second order stationary. For the polynomial kernel h(x) = d1 ∥x∥d2 + 21 ∥x∥22 + 1, the global full domain scenario in Remark 3.1 applies because the objectives and kernels are C 2 and the subproblems are globally strongly convex under the stepsize condition; see Remark 3.2. Remark 3.2 (Global strong convexity of Bregman-proximal subproblems). For the polynomial kernel h(x) = d1 ∥x∥d2 + 1 2 2 n d 2 ∥x∥2 +1 of Proposition 2.1, we have ∇ h(x) ⪰ In for all x ∈ R (since the ∥x∥2 term contributes a positive semidefinite Hessian and the quadratic term contributes In ). By two-sided relative smoothness, ∇2 f (x) ⪰ −Lf ∇2 h(x) for all x, so the Bregman-proximal Hessian satisfies 1  1  1 − Lf ∇2 h(x) ⪰ − L f In ≻ 0 ∇2 f (x) + ∇2 h(x) ⪰ η η η uniformly over Rn whenever η < 1/Lf . Thus each subproblem is globally strongly convex with a unique minimizer, and the fixed-point map g is well defined on the entire primal-dual space. For non-polynomial objectives satisfying relative smoothness, Assumption 3.1 must be verified on a case-by-case basis; the Legendre-type scenario in Remark 3.1 provides a useful template. Since Algorithm 2 reduces to proximal ADMM when h1 (x) = 21 ∥x∥22 and h2 (y) = 12 ∥y∥22 , Theorem 3.1 covers proximal ADMM as a special case. Remark 3.3 (On single-variable vs. bi-variable relative smoothness). Since the objective in (1) is separable, f (x, y) = f1 (x)+f2 (y), the cross-Hessian ∇2xy f vanishes and relative smoothness of each fj in its own variable (Definition 2.4) suffices for Theorem 3.1. For the consensus formulation in Section 4, each local objective fj (xj , yj ) may couple xj and yj , so the stronger bi-smoothness condition (Definition 2.7) is used in Theorem 4.2.

4

Extension to Multi-block Case

The two-block analysis already contains the essential mechanism. We now show how the same approach extends to a global consensus model, which covers the distributed settings most relevant for our applications.

4.1

Global consensus problem

Consider the following unconstrained centralized problem minimize x,{yj }J j=1

J X

fj (x, yj ),

(10)

j=1

where fj are possibly nonconvex functions, which do not necessarily have a Lipschitz continuous gradient. One can distribute this problem across a network of J nodes in a “star topology”, where J − 1 agents are connected to and 7

can communicate with a central node (which we denote by 1), but do not communicate directly among themselves [5]. This results in the following global consensus problem minimize

{xj ,yj }J j=1

J X

fj (xj , yj ) s.t. xj = x1 , ∀j ≥ 2.

(11)

j=1

Problem (11) is a generalized consensus formulation: only the shared variables xj are constrained to agree, while the local variables yj remain agent-specific. The constraints xj = x1 encode a star-topology (global consensus) model in which agent 1 plays the role of a hub and all other agents enforce agreement with x1 . More general consensus topologies (e.g., edge-wise constraints xi = xj over a connected graph) can be written as linear consensus constraints and reformulated into a global consensus form by introducing an auxiliary global variable, but extending the proof to a specific reformulation requires checking the induced fixed-point-map structure in addition to strong convexity of the resulting Bregman-proximal subproblems. For clarity, we focus on the star/global consensus case.

4.2

Landscape of global consensus problem

Theorem 4.1 below connects the geometric landscapes of the unconstrained centralized problem (10) and its equalityconstrained distributed counterpart (11). Theorem 4.1 (Corollary 1 in [29]). [x, y1 , · · · , yJ ] is a first-order/second-order stationary point of problem (10) iff [x, x, · · · , x, y1 , · · · , yJ ] is a first-order/second-order stationary point of problem (11). Moreover, if problem (10) satisfies the strict saddle property (every first order stationary point is either a local minimizer or a strict saddle) and has no spurious local minima, then for every second-order stationary point [x, x, · · · , x, y1 , · · · , yJ ] of (11), [x, y1 , · · · , yJ ] is a global minimizer of (10). Theorem 4.1 links the landscapes of (10) and (11). In particular, obtaining a second-order stationary point of the distributed consensus problem (11) yields a second-order stationary point of the centralized objective (10); in many structured matrix/tensor models, such points are globally optimal under a strict-saddle or no-spurious-minima landscape [13, 14]. Most available algorithms for global consensus problems, however, provide guarantees only for first-order stationarity. Under a Lipschitz gradient assumption, [5] proves that linearized ADMM avoids strict saddles almost surely from random initialization; this assumption can fail for polynomial objectives arising in matrix/tensor factorizations, including distributed matrix factorization [30].

4.3

Convergence of Consensus Bregman ADMM

To address these issues, we extend Bregman ADMM (see Algorithm 2) to the multi-block case. We write the augmented Lagrangian for (11) as L(x1 , y1 , · · · , xJ , yJ , λ2 , · · · , λJ ) =

J X

fj (xj , yj ) +

j=1

J X j=2

J

⟨λj , xj − x1 ⟩ +

ρX ∥xj − x1 ∥22 . 2 j=2

(12)

Applying Bregman-proximal updates to the augmented Lagrangian (12) yields Algorithm 3 below. The agents are updated sequentially (Gauss-Seidel ordering): for j ≥ 2, the xj -update uses the already-computed xk1 . More precisely, the ellipsis in L(xk1 , y1k , . . . , xj , yjk−1 , . . .) means that all blocks with indices ℓ < j take their already-updated iteration-k values, while blocks with ℓ > j retain their iteration-(k−1) values; in particular, for j = 1 no preceding blocks exist, so the Lagrangian is evaluated at (x1 , y1k−1 , xk−1 , . . .). This ordering is also reflected in the proof in 2 Section B, whose Jacobian factorization follows the same update sequence. Theorem 4.2. Assume Assumption 3.1 holds for the fixed-point map induced by Algorithm 3, and assume also that each kernel hj is strongly bi-convex. Let fj be (Lxj , Lyj )-relatively bi-smooth w.r.t. hj , j ∈ [J], and set η < minj∈[J] min(1/Lxj , 1/Lyj ). With random initialization, the probability that the iterates of Algorithm 3 converge to a strict saddle is zero. The proof of Theorem 4.2 is given in Section B. As with Theorem 3.1, the theorem does not assert convergence but rather excludes strict saddles from the set of possible limits; convergence to a KKT point must be established separately. Combined with Theorem 4.1 and additional assumptions guaranteeing first order convergence, for example 8

Algorithm 3 Consensus Bregman ADMM 1: Input: Bregman divergence kernels hj , j ∈ [J] and parameters η, ρ. 0 2: Initialization: (x01 , y10 , · · · , x0J , yJ , λ02 , · · · , λ0J ) 3: For k = 1, 2, . . . 4: For j = 1, 2, . . . , J 1 ; yjk−1 ) + L(xk1 , y1k , · · · , xj , yjk−1 , · · · , xk−1 , yJk−1 , λk−1 xkj = argmin Dhxj (xj , xk−1 , · · · , λk−1 ), 2 J j J η xj 1 yjk = argmin Dhy j (yj , yjk−1 ; xkj ) + L(xk1 , y1k , · · · , xkj , yj , · · · , xk−1 , yJk−1 , λk−1 , · · · , λk−1 ). 2 J J η yj 5:  End For k 6: λ2 · · · 7: End For

  λkJ = λk−1 2

···

  + ρ xk2 − xk1 λk−1 J

···

(13)

 xkJ − xk1 .

assumptions analogous to [22, Theorem 3.9], which include boundedness or coercivity, regularity, and KŁ and subanalyticity type conditions for convergence to a stationary point of the augmented Lagrangian, Theorem 4.2 implies that any limiting KKT point of Algorithm 3 corresponds almost surely to a second-order stationary point of the centralized problem (10). Numerical illustrations are given in Section 6.

5

Strict Saddle Avoidance via Fixed Point Instability

This section provides the high-level ideas underlying the proofs of Theorem 3.1 and Theorem 4.2; the detailed proofs are presented in Section A and Section B, respectively. Specifically, Section 5.2 and Section 5.3 develop the corresponding proof sketches. To this end, we first present a fully expanded proof of the single block Bregman augmented Lagrangian method (Bregman ALM) prototype in Theorem 5.2, whose role is to make the common fixed point template transparent. The argument highlights four ingredients used again in the two block and consensus analyses: construction of a smooth fixed point map, nonsingularity of its Jacobian, identification of KKT points with fixed points, and instability of strict saddle fixed points. We begin with the following single block equality constrained problem: minimize f (x) s.t. Ax = b n x∈R

where A ∈ Rp×n and f is a nonconvex function. The augmented Lagrangian form of (14) is ρ L(x, λ) = f (x) + ⟨λ, Ax − b⟩ + ∥Ax − b∥22 , 2

(14)

(15)

where λ ∈ Rp is the Lagrangian dual variable and ρ > 0 is the augmented penalty parameter. In the single-block case, Bregman ADMM reduces to the Bregman augmented Lagrangian method (Bregman ALM; see Algorithm 4). Algorithm 4 Bregman Augmented Lagrangian Method 1: Input: Bregman divergence kernel h and parameters η, ρ. 0 2: Initialization: (x0 , λ ) 3: For k = 1, 2, . . . 1 xk = argmin L(x, λk−1 ) + Dh (x, xk−1 ), η x k

λ =λ

k−1

(16)

k

+ ρ(Ax − b).

Our analysis follows the dynamical systems viewpoint of Lee et al. [31]: an iterative method defines a smooth fixed point map on the algorithmic state. Once strict saddles are shown to be unstable fixed points of this map, local invertibility of the fixed point map allows the stable manifold argument of [31] to imply that the set of initializations converging to such points has Lebesgue measure zero; see also the stable manifold theorem [15]. 9

Definition 5.1. Let g be a C 1 mapping from Ω to Ω. Its set of unstable fixed points is Ag = {z ∈ Ω : g(z) = z, max |λi (Dg(z))| > 1}. i

(17)

Theorem 5.1 ( [31]). Let g be a C 1 mapping from Ω to Ω and suppose det(Dg(z)) ̸= 0 for all z ∈ Ω (equivalently, g is a local diffeomorphism on Ω). Then the set of initial points that converge to unstable fixed points has zero Lebesgue measure:   Leb z0 ∈ Ω : lim g k (z0 ) ∈ Ag = 0. k→∞

Under the random initialization convention stated in Section 3, Theorem 5.1 shows that the set of initial points that converge to an unstable fixed point has probability zero. In our applications, we verify that strict saddles correspond to such unstable fixed points, so strict-saddle avoidance follows.

5.1

Prototype: strict saddle avoidance for Bregman ALM

Theorem 5.2. Assume Assumption 3.1 holds for the fixed-point map induced by Algorithm 4, and assume also that the kernel h is strongly convex. Let f be Lf -relatively smooth w.r.t. h, and set η < 1/Lf . With random initialization, the probability that the iterates of Algorithm 4 converge to a strict saddle is zero. Note that Theorem 5.2 (and likewise Theorems 3.1 and 4.2) does not assert that the iterates converge; it states that the set of initializations leading to convergence to a strict saddle has measure zero. Convergence itself requires separate conditions (e.g., the KŁ/boundedness assumptions cited earlier); once convergence to a KKT point is established, the theorem excludes strict saddles from the set of possible limits. We apply Theorem 5.1 to the fixed-point map induced by Algorithm 4. The proof proceeds in four steps: 1. Construct the fixed-point map g (so that iterates satisfy zk = g(zk−1 )). 2. Compute the Jacobian Dg (typically by differentiating optimality conditions via the implicit function theorem). 3. Show det(Dg) ̸= 0 on the domain, i.e., the fixed-point map is locally invertible everywhere. 4. Show that any strict saddle corresponds to an unstable fixed point of g. 5.1.1

Constructing the fixed-point map g

To simplify notation, we rewrite (16) as (xk , λk−1 ) = g1 (xk−1 , λk−1 ), (xk , λk ) = g2 (xk , λk−1 ), which implies that

. (xk , λk ) = (g2 ◦ g1 )(xk−1 , λk−1 ) = g(xk−1 , λk−1 ).

(18)

Fixed points and KKT points. For (14), the algorithmic state is z = (x, λ). A point (x⋆ , λ⋆ ) ∈ Ω is a fixed point of the fixed-point map g if and only if it satisfies the KKT conditions Ax⋆ = b,

∇f (x⋆ ) + A⊤ λ⋆ = 0.

At a fixed point, the g1 update returns x+ = x, where x+ denotes the updated primal variable, so the Bregman gradient term ∇h(x+ ) − ∇h(x) vanishes; the x update optimality condition and the dual update equation then reduce to the KKT system. Conversely, if (x⋆ , λ⋆ ) satisfies the KKT conditions, then it satisfies the optimality conditions of (16); since the x-subproblem is strongly convex on the primal domain, its minimizer is unique and therefore x+ = x⋆ , i.e., (x⋆ , λ⋆ ) is a fixed point.

10

5.1.2

Computing the Jacobian matrix Dg

Using the chain rule, we have Dg = Dg2 Dg1 . We first compute Dg2 . For notational convenience, denote (x, λ+ ) = g2 (x, λ). Since g2 is determined by λ+ = λ + ρ(Ax − b),

(19)

we have  Dg2 (x, λ) =

In ∇ x λ+

 ∇ λ λ+

 =

In ρA

 .

Ip

(20)

We now compute Dg1 . To simplify notation, denote (x+ , λ) = g1 (x, λ). Then   ∇x x+ ∇λ x+ Dg1 (x, λ) = Ip

(21)

Relative smoothness implies the positive-semidefinite inequality (Definition 2.4) Lf ∇2 h(x) ± ∇2 f (x) ⪰ 0; in particular, ∇2 f (x+ ) ⪰ −Lf ∇2 h(x+ ). Since h is strongly convex, ∇2 h(x+ ) ≻ 0. Hence for η < 1/Lf , 1  1 − Lf ∇2 h(x+ ) ≻ 0. ∇2 f (x+ ) + ρA⊤ A + ∇2 h(x+ ) ⪰ η η By Assumption 3.1, the x update is well defined. The positive definiteness above makes the x subproblem strongly convex, hence its minimizer is unique, and the implicit function theorem applies to the optimality condition. We use the following standard form of the implicit function theorem [32, Theorem 1.3.1]: if F (u, v) = 0 and ∂v F (u0 , v0 ) is nonsingular, then locally v is a C 1 function of u, and its derivative is given by (22). Dv(u) = − ∂v F (u, v(u))

−1

∂u F (u, v(u)).

(22)

∇h(x+ ) − ∇h(x) = 0. η

(23)

The optimality condition of the first block of (16) is ∇f (x+ ) + A⊤ λ + ρA⊤ (Ax+ − b) +

We take u = (x, λ) and v = x+ , with F (u, v) equal to the left-hand side of (23). The Jacobian ∂v F equals ∇2 f (x+ )+ ρA⊤ A + η1 ∇2 h(x+ ), which is nonsingular by the inequality above. Applying the implicit function theorem to (23), −1   1 1 2 ∇ h(x) ∇2 f (x+ ) + ρA⊤ A + ∇2 h(x+ ) η η  −1 1 ∇λ x+ = − ∇2 f (x+ ) + ρA⊤ A + ∇2 h(x+ ) A⊤ . η

∇ x x+ =



Therefore, plugging (24) into (21), we get " −1 1 2 2 + ⊤ + ∇ f (x ) + ρA A + ∇ h(x ) η Dg (x, λ) =

#

1

Ip

1 2 η ∇ h(x)

(24)

−A⊤ Ip



Finally, by the chain rule, we get that Dg(x, λ) = Dg2 (g1 (x, λ))Dg1 (x, λ), which further implies # −1   "  1 2 ⊤ 1 2 2 + ⊤ + In ∇ f (x ) + ρA A + ∇ h(x ) η ∇ h(x) −A η Dg(x, λ) = . ρA Ip Ip Ip 11

(25)

(26)

5.1.3

Showing that det(Dg) is nonzero globally

Because Dg(x, λ) = Dg2 (g1 (x, λ))Dg1 (x, λ) and Dg1 , Dg2 are square matrices, it suffices to show the global nonsingularity of both Dg1 , Dg2 . Dg1 is nonsingular for all (x, λ) ∈ Ω because it is a product of two nonsingular matrices in (25): the first factor (∇2 f (x+ ) + ρA⊤ A + η1 ∇2 h(x+ ))−1 is nonsingular by the positive definiteness i h1 2 ∇ h(x) −A⊤
= det( η1 ∇2 h(x)) ̸= 0 since h is strongly convex. shown above, and the second factor has det η 0 Ip  In 0  Moreover, Dg2 (x, λ) = ρA Ip by (20), so det(Dg2 ) = 1. 5.1.4

Showing any strict saddle lies in the unstable set

Lemma 5.1. For any stationary point x⋆ of Problem (14), there exists λ⋆ ∈ Rp such that (x⋆ , λ⋆ ) is a fixed point of . the mapping g = g2 ◦ g1 . Proof of Lemma 5.1. By Definition 2.1, for any stationary point x⋆ of (14), there exists λ⋆ such that Ax⋆ = b and ∇f (x⋆ ) + A⊤ λ⋆ = 0. This implies that (x+ , λ+ ) = (x⋆ , λ⋆ ) and (x, λ) = (x⋆ , λ⋆ ) satisfy (23) (defining g1 ) and (19) (defining g2 ) and therefore (x⋆ , λ⋆ ) is a fixed point of g2 ◦ g1 . Lemma 5.2. Let x⋆ be a strict saddle of Problem (14), and let λ⋆ be a corresponding multiplier such that (x⋆ , λ⋆ ) ∈ Ω. Then the Jacobian matrix Dg(x⋆ , λ⋆ ) has an eigenvalue with magnitude greater than 1. Proof of Lemma 5.2. To simplify notation, we denote H = ∇2 h(x⋆ ), F = ∇2 f (x⋆ ). We compute the Jacobian matrix Dg(x⋆ , λ⋆ ) by plugging (x, λ) = (x+ , λ+ ) = (x⋆ , λ⋆ ) into (26): 1   ⊤ (F + ρA⊤ A + η1 H)−1 In η H −A ρA Ip Ip Ip  −1  1  1 ⊤ ⊤ F + ρA A + η H η H −A = −ρA Ip Ip     −1 F + ρA⊤ A + η1 H F + ρA⊤ A A⊤ =I− −ρA −ρA Ip . = I − Φ.

Dg(x⋆ , λ⋆ ) =



It suffices to show that Φ has a real negative eigenvalue (equivalently, det(Φ + µI) = 0 for some µ > 0). For µ > 0, taking the Schur complement with respect to the bottom-right block µIp gives   (1 + µ)F + (1 + µ)ρA⊤ A + µη H A⊤ =0 0 = det(Φ + µI) ⇐⇒ det −(1 + µ)ρA µIp   µ (1 + µ)2 p ⊤ ⇐⇒ µ det (1 + µ)F + H + ρA A = 0 η µ   µ (1 + µ)2 ⊤ ⇐⇒ det (1 + µ)F + H + ρA A = 0. η µ Note that (1 + µ)2 µ . ρA⊤ A J(µ) = (1 + µ)F + H + η µ

(27)

is real symmetric for every µ > 0, and its entries depend continuously on µ. Remark 5.1 (Intuition for the test direction). The strict-saddle condition (Definition 2.1, item 3) provides z ∈ null(A) with z⊤ ∇2 f (x⋆ )z < 0. Since Az = 0, we have z⊤ A⊤ Az = ∥Az∥2 = 0. Thus the penalty term in z⊤ J(µ)z vanishes for every µ > 0, so along this feasible direction only the Hessian and Bregman kernel terms remain. At µ → 0+ the Bregman kernel term µη H vanishes as well, so the sign of z⊤ J(µ)z is governed by the negative Hessian 2

alone, giving λmin (J(µ)) < 0. At µ → ∞ both µη H ≻ 0 and (1+µ) ρA⊤ A ⪰ 0 dominate, forcing J(µ) ≻ 0. µ Continuity of eigenvalues then yields a zero crossing of λmin (J(µ)) at some intermediate µ⋆ > 0. 12

First, we show that for µ → 0+ , λmin (J(µ)) < 0. Using the third item of Definition 2.1 and that F = ∇2 f (x⋆ ), there exists z ∈ null(A) such that z⊤ Fz < −σ∥z∥22 for some σ > 0. Now   µ ⊤ 1 ⊤ µ ⊤ 2 2 2 ⊤ ⊤ z J(µ)z = (1 + µ)z Fz + z Hz < −σ(1 + µ)∥z∥2 + z Hz = −σ∥z∥2 − µ σ∥z∥2 − z Hz . η η η The first term −σ∥z∥22 is strictly negative and independent of µ, while the second term is O(µ); hence the right-hand side is negative for all sufficiently small µ > 0. Since J(µ) is a symmetric matrix, this implies that λmin (J(µ)) < 0 when µ is a sufficiently small positive number. Second, note that lim

µ→∞

H J(µ) =F+ + ρA⊤ A ≻ 0, µ η

where the positive definiteness follows from F ⪰ −Lf H (relative smoothness) and η1 H ≻ Lf H (since η < 1/Lf and H ≻ 0), giving F + η1 H ≻ 0, plus ρA⊤ A ⪰ 0. Since J(µ) is symmetric, its eigenvalues are real and continuous in µ [33, Theorem 5.1]; the intermediate value theorem then gives λmin (J(µ)) = 0 for some µ > 0, i.e., det(J(µ)) = 0. By Theorem 5.1, the set of initial points from which Algorithm 4 converges to a strict saddle has Lebesgue measure zero. Since random initialization is absolutely continuous with respect to Lebesgue measure, the probability that the iterates converge to a strict saddle is zero.

5.2

Proof sketch for the two-block case (Theorem 3.1)

The full proof is in Section A; here we outline the key steps and highlight where the argument differs from the singleblock case above. Step 1: Fixed-point map and Jacobian. One iteration of Algorithm 2 is decomposed as g = g3 ◦ g2 ◦ g1 (primal xupdate, primal y-update, dual ascent). By the implicit function theorem applied to the two primal optimality conditions, g1 and g2 are C 1 on Ω; their Jacobians are nonsingular because the corresponding subproblem Hessians with respect to the updated primal variables are positive definite. The dual update g3 has a nonsingular block triangular Jacobian, and hence Dg = Dg3 Dg2 Dg1 is well defined and nonsingular on Ω. Step 2: Jacobian at a fixed point. with

At a strict-saddle fixed point (x⋆ , y⋆ , λ⋆ ), the Jacobian evaluates to Dg = I − Φ Φ = M−1 N,

where M collects the Bregman-regularized subproblem Hessian blocks, while N collects the remaining curvature and primal dual coupling blocks arising from the augmented Lagrangian and the Gauss Seidel ordering (the explicit expression is given in (48) of Section A). It suffices to find µ > 0 with det(Φ + µI) = 0. Step 3: Schur reduction and symmetrization. After Schur-complementing the dual block (as in the single-block case), the determinant condition reduces to the singularity of a 2 × 2 block matrix C(µ) ∈ R(n+m)×(n+m) . Unlike the single-block case, C(µ) is not symmetric: its (1, 2) and (2, 1) blocks carry different scalar coefficients depending p on µ. The remedy is a diagonal similarity D = diag(In , α(µ)Im ) with α(µ) = (2 + µ + 1/µ)/(2 + 1/µ), which symmetrizes C(µ) to a real-symmetric matrix J(µ) = D−1 C(µ)D without changing its determinant. Step 4: Sign change via a µ-dependent test direction. The symmetrization changes the natural test direction (dx , dy ) to (dx , s(µ)dy ) with s(µ) = α(µ) → 1 as µ → 0+ . Using the KKT constraint Adx + Bdy = 0 and the 2 ⊤ 2 strict-saddle negativity d⊤ x ∇ f1 dx + dy ∇ f2 dy < 0, one shows λmin (J(µ)) < 0 for small µ > 0, while J(µ)/µ → diag(F1 + η1 H1 + ρA⊤ A, F2 + η1 H2 + ρB⊤ B) ≻ 0 gives λmin (J(µ)) > 0 for large µ. The intermediate value theorem then yields the desired zero crossing. Equivalently, det(Φ + µI) = 0 for some µ > 0, so Dg = I − Φ has an eigenvalue with magnitude greater than one.

13

5.3

Proof sketch for the consensus case (Theorem 4.2)

The full proof is in Section B. Here we indicate how the single block mechanism and the two block spectral argument extend to the star consensus formulation of Section 4. Step 1: Fixed point map and nonsingularity. One iteration of Algorithm 3 defines a smooth fixed point map on the state (x1 , y1 , . . . , xJ , yJ , λ2 , . . . , λJ ). The update decomposes into 2J + 1 elementary maps, one xj update and one yj update for each agent, followed by the dual update. Relative bi-smoothness, strong bi-convexity of the kernels, and the stepsize condition make each local Bregman proximal Hessian positive definite, so the Jacobian of each xj and yj update is nonsingular by the implicit function theorem. The dual update has a nonsingular block triangular Jacobian. Hence every elementary Jacobian is nonsingular, and the chain rule gives det(Dg) ̸= 0 on Ω. Step 2: Jacobian at a strict saddle fixed point. At a strict saddle KKT point of the consensus problem, the Jacobian again has the form Dg = I − Φ. As in the two block case, it is enough to find µ > 0 such that det(Φ + µI) = 0. After a positive row column scaling that symmetrizes the determinant condition and after eliminating the dual blocks by Schur complementation, the condition reduces to the singularity of a real symmetric matrix J̃(µ) built from the local Hessian blocks and the star consensus coupling. Step 3: Star consensus cancellation. The key new feature is the star graph structure. A feasible strict saddle direction has equal x components across all agents, say (dx , d1y , . . . , dx , dJy ), because the constraints impose xj = x1 . This equal x direction lies in the kernel of the star graph Laplacian. Consequently, the consensus penalty and the Schur complement correction vanish on the test direction, leaving only the block Hessian contribution. As µ → 0+ , the corresponding quadratic form has a well-defined limit, and this limit is exactly the strict saddle negative curvature. Step 4: Large µ positivity and instability. For large µ, the scaled matrix J̃(µ)/µ converges to a block diagonal matrix whose diagonal blocks are the positive definite Bregman proximal Hessian blocks. Thus λmin (J̃(µ)) > 0 for large µ, while the previous step gives λmin (J̃(µ)) < 0 for small positive µ. Continuity of the eigenvalues gives a zero crossing, hence an eigenvalue of Dg with magnitude greater than one. Applying Theorem 5.1 then excludes convergence to strict saddle KKT points almost surely.

6

Numerical Experiments

We consider two polynomial models: distributed matrix factorization and symmetric tensor factorization. Both fall outside the globally Lipschitz-smooth regime. The former fits the separable consensus framework of Theorem 4.2 and admits an explicit relative-bi-smoothness certificate; the latter is included because the same splitting idea leads to explicit Bregman-proximal updates even though the resulting objective is not separable across agents. For distributed matrix factorization, the corresponding centralized problem has a benign strict saddle landscape under standard assumptions, so second order stationary points are globally optimal; the tensor example is used as an algorithmic illustration of the same Bregman proximal splitting idea beyond the separable setting.

6.1

Distributed matrix factorization

  Given a data matrix Z ∈ Rn×m with column partitioning Z = Z1 Z2 · · · ZJ , where Zj ∈ Rn×mj , suppose rank(Z) = r and let X⋆ ∈ Rn×r span the column space of Z. The centralized and distributed matrix factorization problems are Centralized : minimize X,{Yj }

Distributed : minimize

{Xj },{Yj }

J X

∥XYj⊤ − Zj ∥2F ;

(28)

∥Xj Yj⊤ − Zj ∥2F s.t. Xj = X1 ∀j ≥ 2.

(29)

j=1 J X j=1

The distributed formulation enforces exact consensus on the shared factor while allowing local right factors, a standard model in distributed low-rank matrix recovery [30]. 14

. Applicability of Theorem 4.2. Each local loss fj (Xj , Yj ) = ∥Xj Yj⊤ − Zj ∥2F is a polynomial of bidegree (2, 2) in the entries of (Xj , Yj ), hence relatively bi-smooth under a suitable kernel by Proposition 2.2. In contrast, many standard consensus algorithms (e.g., linearized ADMM [5], decentralized gradient descent [34], and EXTRA [35]) require global Lipschitz gradient assumptions, which do not hold for such polynomial objectives. The centralized formulation (28) has a strict-saddle landscape with no spurious local minima under the rank condition above and standard identifiability conditions [30, 36], so second-order stationary points are global minimizers. With Xj ∈ Rn×r and Yj ∈ Rmj ×r , we use the product-type bi-kernel    . hj (Xj , Yj ) = 12 ∥Xj ∥2F + 1 12 ∥Yj ∥2F + 1 , (30)

Kernel choice and block subproblems.

. . which is strongly bi-convex. Writing hj = hx,j · hy,j with hx,j (Xj ) = 12 ∥Xj ∥2F + 1 and hy,j (Yj ) = 12 ∥Yj ∥2F + 1, the product-kernel identity in Section 2 gives   Dhxj (Xj , X′j ; Yj′ ) = hy,j (Yj′ ) Dhx,j (Xj , X′j ) = 21 12 ∥Yj′ ∥2F + 1 ∥Xj − X′j ∥2F ,   Dhy j (Yj , Yj′ ; X′j ) = hx,j (X′j ) Dhy,j (Yj , Yj′ ) = 12 21 ∥X′j ∥2F + 1 ∥Yj − Yj′ ∥2F . (31) αk

; Yjk−1 ) contributes a quadratic proximal regularizer 2j ∥Xj − In Algorithm 3, the Bregman term η1 Dhxj (Xj , Xk−1 j βk . Xk−1 ∥2F with iteration-dependent coefficient αjk = η1 ( 21 ∥Yjk−1 ∥2F + 1); the subsequent Yj -update uses 2j ∥Yj − j . Yjk−1 ∥2F with βjk = η1 ( 12 ∥Xkj ∥2F + 1). For a non-hub agent (j ≥ 2), setting the gradient of the Xj -subproblem to zero yields the r × r normal equation
Xkj 2Yjk−1⊤ Yjk−1 + γjk Ir = 2Zj Yjk−1 − λk−1 + ρXk1 + αjk Xk−1 , j j . where γjk = ρ + αjk (the hub-node update has the same structure with (J−1)ρ replacing ρ). Since γjk > 0, the coefficient matrix 2Yjk−1⊤ Yjk−1 + γjk Ir ∈ Rr×r is positive definite, so each block update reduces to an r × r linear k k ⊤ k solve. The Yj -update carries no consensus constraint; its normal equation is Yjk (2Xk⊤ j Xj + βj Ir ) = 2Zj Xj + k k−1 βj Y j . Results. We generate a noiseless rank-r instance Z = AB⊤ with A ∈ Rn×r and B ∈ Rm×r having independent standard normal entries, and partition Z into J blocks. The primal variables {X0j }, {Yj0 } and the multipliers are initialized randomly. In Figure 1, we use a larger instance with n = 50, mj = 50 for all j, r = 4, J = 100, and (η, ρ) = (1, 1000).1 The objective decreases rapidly to approximately 6.43 × 10−25 and the consensus residual to approximately 7.31 × 10−28 after 600 iterations, showing that the algorithm both fits the data and enforces agreement among the distributed factors at machine precision.

6.2

Symmetric tensor factorization

As a second application, we consider the symmetric tensor factorization problem minimize ∥T − n×r U∈R

r X

ui ⊗ ui ⊗ ui ∥2F ,

(32)

i=1

where U = [u1 , · · · , ur ], u ⊗ u ⊗ u denotes the symmetric rank-one tensor with (i, j, k)-th entry ui uj uk , and T ∈ Rn×n×n is a given symmetric tensor. Symmetric tensor decomposition arises widely in signal processing and statistical estimation, including method of moments estimators [37–40]. Problem (32) is structurally more coupled than a standard CP model because the same factor matrix U enters all three modes. Consequently, the usual mode-wise alternating least-squares strategy for nonsymmetric CP, whose 1 The guarantee of Theorem 4.2 applies when η satisfies the corresponding relative-bi-smoothness stepsize bound. For the kernel above, one may y take Lx j = Lj = 4, so η < 1/4 is sufficient. We use the larger practical stepsize η = 1 for speed; the normal equations above show that the block subproblems remain strongly convex because γjk > 0 and βjk > 0.

15

Objective:

102

Objective

j

Consensus:

10 2

XjkYjk j

Zj 2F

Xjk X1k 2F

103 10 1 10 5

Consensus residual

106

10 6

10 9

10 10

10 13

10 14

10 17

10 18

10 21

10 22

10 25 0

100

200

300

Iteration

400

500

600

Figure 1: Consensus Bregman ADMM on distributed matrix factorization. block updates reduce to linear least-squares problems, is not directly available in the symmetric formulation [37]. A common remedy is to introduce separate mode factors and enforce symmetry through exact coupling constraints, or to encourage symmetry through penalties in related penalized formulations. The exact constraint formulation restores blockwise multilinearity while preserving equivalence with the original symmetric model. Related formulations appear in numerical optimization approaches to symmetric CP decomposition and, more broadly, in structured/coupled factorization frameworks [41, 42]. Motivated by this observation, we introduce auxiliary factors V, W and impose symmetry via consensus constraints, thereby rewriting (32) as minimize ∥T −

U,V,W∈Rn×r

r X

ui ⊗ vi ⊗ wi ∥2F s.t. U = V = W,

(33)

i=1

where V = [v1 , · · · , vr ], W = [w1 , · · · , wr ], and u ⊗ v ⊗ w is a rank-one tensor with (i, j, k)-th entry ui vj wk . Motivated by the consensus update in Algorithm 3, we apply the same Bregman proximal update pattern to (33), treating (U, V, W) as three consensus coupled copies of the factor. Remark P 6.1 (Scope of the theoretical guarantee). The consensus model (11) assumes a separable objective of the form j fj (xj , yj ), with each term depending only on agent j’s own variables. The symmetric tensor splitting (33) P does not have this structure, because the single loss ∥T − i ui ⊗ vi ⊗ wi ∥2F still couples the three factor copies. For this reason, Theorem 4.2 does not apply to (33) as stated. Even so, the Bregman-proximal subproblems remain strongly convex, as the normal equations below show, so the algorithm is perfectly well defined. We include this example primarily as an algorithmic illustration: it demonstrates the practical value of the Bregman regularization in stabilizing subproblems for polynomial objectives beyond the Lipschitz-smooth regime. We note that the key obstacle to extending Theorem 4.2 to the non-separable case is that the fixed-point Jacobian acquires nonzero off-diagonal blocks between agents (through the coupled loss), so the block-sparsity exploited in Section B no longer holds. A natural conjecture is that the same instability conclusion persists whenever the coupling is “sufficiently sparse” relative to the consensus graph—formalizing this is an open direction discussed in the conclusion. With U, V, W ∈ Rn×r , we use the three-factor product kernel     . h(U, V, W) = 21 ∥U∥2F + 1 12 ∥V∥2F + 1 12 ∥W∥2F + 1 .

Kernel choice and block subproblems.

(34)

Fixing two factors and applying the product-kernel identity from Section 2 to the third gives a scaled Frobenius proximal term. (The identity in Section 2 is stated for two-factor product kernels; for the three-factor case one simply treats the product of the two fixed factors as a single scalar, reducing to the same formula.) For the U-variable:    DhU (U, U′ ; V′ , W′ ) = 12 12 ∥V′ ∥2F + 1 21 ∥W′ ∥2F + 1 ∥U − U′ ∥2F , (35) 16

and analogously for V and W by cyclic permutation of the three factors. In Algorithm 3, η1 DhU (U, Uk−1 ; Vk−1 , Wk−1 ) pk

contributes a proximal regularizer 2U ∥U − Uk−1 ∥2F with coefficient . 1 pkU = η



k−1 2 1 ∥F + 1 2 ∥V



k−1 2 1 ∥F + 1 2 ∥W



, 2

and similarly pkV , pkW by cyclic permutation. The loss can be written via the mode-1 matricization T(1) ∈ Rn×n (the matrix whose rows index the first mode and whose columns enumerate all index pairs of the remaining two modes) as 2 . ∥T(1) − UQ⊤ ∥2F , where Q = W ⊛ V ∈ Rn ×r is the Khatri–Rao (column-wise Kronecker) product [43]. Since U is the hub node with J−1 = 2 consensus constraints, setting the gradient of the U-subproblem to zero gives the r × r normal equation
k−1 k + Wk−1 ) + pkU Uk−1 , Uk 2Qk−1⊤ Qk−1 + γU + λk−1 Ir = 2T(1) Qk−1 + λk−1 W + ρ(V V . k . where Qk−1 = Wk−1 ⊛ Vk−1 and γU = 2ρ + pkU > 0. The V- and W-updates (non-hub agents) have the same structure with mode-2 and mode-3 matricizations T(2) , T(3) and the corresponding Khatri–Rao products Wk−1 ⊛ Uk and Vk ⊛ Uk . Pr Results. We form T = i=1 u⋆i ⊗ u⋆i ⊗ u⋆i from a random U⋆ ∈ Rn×r and initialize (U0 , V0 , W0 ) as Gaussian perturbations of U⋆ with noise level 0.01 and zero multipliers. In Figure 2, we use a larger instance with n = 64, r = 8, and (η, ρ) = (8 × 103 , 100). The reconstruction error decreases to approximately 2.55 × 10−17 and the symmetry residual ∥U − V∥2F + ∥U − W∥2F to approximately 3.18 × 10−22 after 1000000 iterations. In other words, the split formulation does what it is supposed to do: it drives the factors toward agreement while continuing to reduce the tensor-fitting objective. The run is slower than in the matrix-factorization example, which is unsurprising here because the admissible Bregman regularization is much stronger and the underlying polynomial has higher degree.

10 3

ui vi wi 2F 10 3 Consensus: U V 2F + U W 2F 10 6

10 6

10 9

Reconstruction:

i

Consensus residual

Reconstruction error

100

10 12

10 9

10 15

10 12

10 18

10 15

10 21 0.0

0.2

0.4

0.6

Iteration

0.8

1.0 1e6

Figure 2: Consensus Bregman ADMM on symmetric tensor factorization.

7

Conclusion

The main point of the paper is that the strict-saddle-avoidance mechanism familiar from smooth first-order methods survives in the Bregman ADMM setting, even when global Lipschitz smoothness is unavailable. What is nontrivial is not the dynamical-systems principle itself, but the verification that it continues to apply after replacing Euclidean proximal terms by Bregman regularization. In the two block case this requires a determinant reduction and symmetrization and scaling step specific to the Bregman setting in the spectral argument, while in the consensus case it relies on the special hub and peripheral structure of the star graph and the associated null space cancellation of the 17

consensus penalty. Once these ingredients are in place, the iteration can again be viewed as a smooth primal–dual map on an invariant open domain, so strict saddle KKT points correspond to unstable fixed points and random initialization rules out convergence to them almost surely. Combined with standard first-order convergence results, this gives second-order stationarity of limiting KKT points. Several natural questions remain. In the two-block setting, our analysis assumes a separable objective f (x, y) = f1 (x) + f2 (y), allowing nonzero cross-Hessian terms ∇2xy f would widen the scope of the theory. In the consensus setting, the current result still relies on sum-separability across agents, whereas the symmetric tensor splitting in Remark 6.1 suggests that interesting non-separable consensus formulations should also be accessible. The key technical barrier is that non-separability introduces off-diagonal Jacobian blocks between agents, breaking the star-shaped block sparsity used in Section B. A promising approach is to treat these inter-agent couplings as perturbations of the separable Jacobian and show that the negative-curvature direction still dominates for small enough coupling, we leave a rigorous treatment to future work. It would also be worthwhile to understand whether the same dynamical picture persists for stochastic or mini-batch variants of Bregman ADMM.

Acknowledgment This work was supported in part by NSF grants ECCS-2409701 and ECCS-2409702.

References [1] R. Glowinski and A. Marroco. Sur l’approximation, par éléments finis d’ordre un, et la résolution, par pénalisation-dualité d’une classe de problèmes de Dirichlet non linéaires. Revue française d’automatique, informatique, recherche opérationnelle. Analyse numérique, 9(R2):41–76, 1975. [2] Daniel Gabay and Bertrand Mercier. A dual algorithm for the solution of nonlinear variational problems via finite element approximation. Computers & Mathematics with Applications, 2(1):17–40, 1976. [3] MH Xu. Proximal alternating directions method for structured variational inequalities. Journal of Optimization Theory and Applications, 134(1):107–117, 2007. [4] Max L. N. Goncalves, Jefferson G. Melo, and Renato D. C. Monteiro. Convergence rate bounds for a proximal ADMM with over-relaxation stepsize parameter for solving nonconvex linearly constrained problems, 2017. [5] Mingyi Hong, Meisam Razaviyayn, and Jason Lee. Gradient primal-dual algorithm converges to second-order stationary solution for nonconvex distributed optimization over networks. In Jennifer Dy and Andreas Krause, editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 2009–2018. PMLR, 2018. [6] Yangyang Xu, Wotao Yin, Zaiwen Wen, and Yin Zhang. An alternating direction algorithm for matrix completion with nonnegative factors. Frontiers of Mathematics in China, 7(2):365–384, 2012. [7] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, Jonathan Eckstein, et al. Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends® in Machine learning, 3(1):1–122, 2011. [8] Mingyi Hong, Zhi-Quan Luo, and Meisam Razaviyayn. Convergence analysis of alternating direction method of multipliers for a family of nonconvex problems. SIAM Journal on Optimization, 26(1):337–364, 2016. [9] Fenghui Wang, Zongben Xu, and Hong-Kun Xu. Convergence of Bregman alternating direction method with multipliers for nonconvex composite problems, 2014. [10] Yu Wang, Wotao Yin, and Jinshan Zeng. Global convergence of ADMM in nonconvex nonsmooth optimization. Journal of Scientific Computing, 78:29–63, 2019. [11] Rina Foygel Barber and Emil Y. Sidky. Convergence for nonconvex ADMM, with applications to CT imaging. Journal of Machine Learning Research, 25(38):1–46, 2024.

18

[12] Gabriel Mancino-Ball, Yangyang Xu, and Jie Chen. A decentralized primal-dual framework for non-convex smooth consensus optimization. IEEE Transactions on Signal Processing, 71:525–538, 2023. [13] Rong Ge, Jason D Lee, and Tengyu Ma. Matrix completion has no spurious local minimum. In Advances in Neural Information Processing Systems 29, pages 2973–2981, 2016. [14] Rong Ge, Furong Huang, Chi Jin, and Yang Yuan. Escaping from saddle points—online stochastic gradient for tensor decomposition. In Conference on Learning Theory, pages 797–842, 2015. [15] Michael Shub. Global stability of dynamical systems. Springer Science & Business Media, 2013. [16] Ioannis Panageas, Georgios Piliouras, and Xiao Wang. First-order methods almost always avoid saddle points: The case of vanishing step sizes. In Advances in Neural Information Processing Systems, volume 32, 2019. [17] Damek Davis and Dmitriy Drusvyatskiy. Proximal methods avoid active strict saddles of weakly convex functions. Foundations of Computational Mathematics, 22(2):561–606, 2022. [18] Heinz H Bauschke, Jérôme Bolte, and Marc Teboulle. A descent lemma beyond Lipschitz gradient continuity: First-order methods revisited and applications. Mathematics of Operations Research, 42(2):330–348, 2016. [19] Jérôme Bolte, Shoham Sabach, Marc Teboulle, and Yakov Vaisbourd. First order methods beyond convexity and Lipschitz gradient continuity with applications to quadratic inverse problems. SIAM Journal on Optimization, 28(3):2131–2151, 2018. [20] Yair Censor and Arnold Lent. An iterative row-action method for interval convex programming. Journal of Optimization theory and Applications, 34(3):321–353, 1981. [21] Lev M Bregman. The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming. USSR computational mathematics and mathematical physics, 7(3):200–217, 1967. [22] Fenghui Wang, Wenfei Cao, and Zongben Xu. Convergence of multi-block Bregman ADMM for nonconvex composite problems. Science China Information Sciences, 61(12):122101, 2018. [23] Huahua Wang and Arindam Banerjee. Bregman alternating direction method of multipliers. In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K.Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 27. Curran Associates, Inc., 2014. [24] Peng-Jie Liu, Jin-Bao Jian, Hu Shao, Xiao-Quan Wang, Jia-Wei Xu, and Xiao-Yu Wu. A Bregman-style improved ADMM and its linearized version in the nonconvex setting: Convergence and rate analyses. Journal of the Operations Research Society of China, 12(2):298–340, 2024. [25] Tan Nhat Pham, Minh N. Dao, Andrew Eberhard, and Nargiz Sultanova. Bregman proximal linearized ADMM for minimizing separable sums coupled by a difference of functions. Journal of Optimization Theory and Applications, 203(2):1622–1658, 2024. [26] Tian-Tian Cui, Ya-Zheng Dang, and Yan Gao. Distributed multi-block partially symmetric Bregman ADMM for nonconvex and nonsmooth sharing problem. Journal of the Operations Research Society of China, 2025. [27] Jingyu Gao and Xiurui Geng. A note on the second-order convergence of proximal ADMM for non-convex optimization. Journal of Optimization Theory and Applications, 209(1):1–17, April 2026. [28] Qiuwei Li, Zhihui Zhu, Gongguo Tang, and Michael B. Wakin. Provable Bregman-divergence based methods for nonconvex and non-Lipschitz problems, 2019. arXiv preprint arXiv:1904.09712. [29] Qiuwei Li, Zhihui Zhu, Gongguo Tang, and Michael B Wakin. The geometry of equality-constrained global consensus problems. In ICASSP, pages 7928–7932, 2019. [30] Zhihui Zhu, Qiuwei Li, Xinshuo Yang, Gongguo Tang, and Michael B Wakin. Distributed low-rank matrix factorization with exact consensus. Advances in Neural Information Processing Systems, 32, 2019.

19

[31] Jason D Lee, Ioannis Panageas, Georgios Piliouras, Max Simchowitz, Michael I Jordan, and Benjamin Recht. First-order methods almost always avoid strict saddle points. Mathematical Programming, 176(1–2):311–337, 2019. [32] Steven G Krantz and Harold R Parks. Introduction to the implicit function theorem. In The Implicit Function Theorem, pages 1–12. Springer, 2003. [33] Tosio Kato. Perturbation theory for linear operators, volume 132. Springer Science & Business Media, 2013. [34] Jinshan Zeng and Wotao Yin. On nonconvex decentralized gradient descent. IEEE Signal Process., 66(11):2834– 2848, 2018. [35] Wei Shi, Qing Ling, Gang Wu, and Wotao Yin. EXTRA: An exact first-order algorithm for decentralized consensus optimization. SIAM Journal on Optimization, 25(2):944–966, 2015. [36] Maher Nouiehed and Meisam Razaviyayn. Learning deep models: Critical points and local openness. INFORMS Journal on Optimization, 4(2):148–173, 2022. [37] Tamara G Kolda and Brett W Bader. Tensor decompositions and applications. SIAM review, 51(3):455–500, 2009. [38] Pierre Comon, Gene Golub, Lek-Heng Lim, and Bernard Mourrain. Symmetric tensors and symmetric tensor rank. SIAM Journal on Matrix Analysis and Applications, 30(3):1254–1279, 2008. [39] Animashree Anandkumar, Rong Ge, Daniel Hsu, Sham M Kakade, and Matus Telgarsky. Tensor decompositions for learning latent variable models. Journal of Machine Learning Research, 15(1):2773–2832, 2014. [40] Shuang Li and Qiuwei Li. Local and global convergence of general Burer-Monteiro tensor optimizations. Proceedings of the AAAI Conference on Artificial Intelligence, 36(9):10266–10274, 2022. [41] Tamara G Kolda. Numerical optimization for symmetric tensor decomposition. Mathematical Programming, 151(1):225–248, 2015. [42] Laurent Sorber, Marc Van Barel, and Lieven De Lathauwer. Structured data fusion. IEEE Journal of Selected Topics in Signal Processing, 9(4):586–600, 2015. [43] Shuangzhe Liu and Goetz Trenkler. Hadamard, khatri-rao, kronecker and other matrix products. International Journal of Information and Systems Sciences, 4(1):160–177, 2008.

20

A PPENDIX A

Proof of Theorem 3.1

Proof of Theorem 3.1. We verify the two conditions of Theorem 5.1: det(Dg) ̸= 0 everywhere, and every strict-saddle KKT point in Ω is an unstable fixed point of g. One iteration of Algorithm 2 decomposes into three elementary maps g1 , g2 , g3 (defined below). By Assumption 3.1, the full update, including the intermediate elementary updates g1 , g2 , g3 , is well defined on the relevant open state space domains. The relative smoothness and stepsize assumptions make the Bregman proximal Hessians positive definite, so any minimizer is unique and the implicit function theorem applies to the optimality conditions. The nonsingularity of Dg is verified below. For instability, the Bregman proximal term contributes a positive definite shift that cannot cancel the negative curvature at a strict saddle, so the linearized map has an eigenvalue of magnitude greater than one. Constructing the fixed-point map g.

Represent (9) as (x+ , y, λ) = g1 (x, y, λ), (x, y+ , λ) = g2 (x, y, λ),

(36)

+

(x, y, λ ) = g3 (x, y, λ). The elementary updates are well defined on the corresponding open domains by Assumption 3.1; the positive definiteness of the Bregman proximal Hessians gives local uniqueness and smooth dependence of the primal solution maps. Then Algorithm 2 iterates the composite map (x+ , y+ , λ+ ) = g(x, y, λ),

(37)

. g = g3 ◦ g2 ◦ g1 .

(38)

with

Computing the Jacobian matrix Dg of the fixed-point map for Algorithm 2.

By the chain rule, we have

Dg = Dg3 Dg2 Dg1 . • Computing Dg3 . Since (x, y, λ+ ) = g3 (x, y, λ) with g3 defined by the optimality condition of the third-block of (9) λ+ = λ + ρ(Ax + By − b),

(39)

then we get 



In

Dg3 (x, y, λ) =  ρA

Im ρB

• Computing Dg2 . Since (x, y+ , λ) = g2 (x, y, λ), then we get  In Dg2 (x, y, λ) = ∇x y+ ∇y y+

.

(40)

Ip

 ∇λ y +  , Ip

Before invoking the implicit function theorem, we note that relative smoothness means Li ∇2 hi (z) ± ∇2 fi (z) ⪰ 0 for i = 1, 2, hence ∇2 fi (z) ⪰ −Li ∇2 hi (z). Since each hi is strongly convex, ∇2 hi (z) ≻ 0 on the domain. Hence for η < 1/Li , 1  1 − Li ∇2 hi (z) ≻ 0, i = 1, 2. ∇2 fi (z) + ∇2 hi (z) ⪰ η η 21

Adding the positive semidefinite quadratic terms from the augmented Lagrangian (e.g., ρA⊤ A or ρB⊤ B) preserves positive definiteness. Therefore each primal subproblem in (9) is strongly convex with a unique minimizer, and the implicit function theorem applies to the optimality conditions. Specifically, ∇x y+ , ∇y y+ , ∇λ y+ are obtained by applying the implicit function theorem to the optimality condition of the second block of (9) ∇f2 (y+ ) + B⊤ λ + ρB⊤ (Ax + By+ − b) +

∇h2 (y+ ) − ∇h2 (y) = 0. η

(41)

Applying the implicit function theorem to (41) gives 1 ∇2 f2 (y+ )∇x y+ + ρB⊤ A + ρB⊤ B∇x y+ + ∇2 h2 (y+ )∇x y+ = 0 η −1 
1 ρB⊤ A ; =⇒∇x y+ = − ∇2 f2 (y+ ) + ρB⊤ B + ∇2 h2 (y+ ) η 1 1 ∇2 f2 (y+ )∇y y+ + ρB⊤ B∇y y+ + ∇2 h2 (y+ )∇y y+ − ∇2 h2 (y) = 0 η η  −1   1 1 2 =⇒∇y y+ = ∇2 f2 (y+ ) + ρB⊤ B + ∇2 h2 (y+ ) ∇ h2 (y) ; η η 1 ∇2 f2 (y+ )∇λ y+ + B⊤ + ρB⊤ B∇λ y+ + ∇2 h2 (y+ )∇λ y+ = 0 η  −1 1 =⇒∇λ y+ = − ∇2 f2 (y+ ) + ρB⊤ B + ∇2 h2 (y+ ) B⊤ . η Therefore, Dg2 (x, y, λ)  In = ∇x y+ ∇y y+  In  =

 ∇λ y +  Ip 



In  −ρB⊤ A

−1  2 ∇ f2 (y+ ) + ρB⊤ B + η1 ∇2 h2 (y+ ) Ip

 In  =



In  −ρB⊤ A

 2 −1 ∇ f2 (y+ ) + ρB⊤ B + η1 ∇2 h2 (y+ )

−B⊤  Ip  In 1 2 ⊤  −B η ∇ h2 (y) Ip

1 2 η ∇ h2 (y)

Im

Ip

  (42) Ip

where ∇2 f2 (y+ ) + ρB⊤ B + η1 ∇2 h2 (y+ ) is positive definite by relative smoothness, strong convexity of h2 , and η < 1/L2 (and adding the positive semidefinite term ρB⊤ B preserves positive definiteness). • Computing Dg1 . Since (x+ , y, λ) = g1 (x, y, λ), then we get  ∇x x + ∇y x +  Im Dg1 (x, y, λ) =

∇λ x +

 ,

Ip where ∇x x+ , ∇y x+ , ∇λ x+ can be obtained by applying the implicit function theorem to the optimality condition of the first block of (9) ∇f1 (x+ ) + A⊤ λ + ρA⊤ (Ax+ + By − b) +

∇h1 (x+ ) − ∇h1 (x) = 0. η

Applying the implicit function theorem to (43) gives 1 1 ∇2 f1 (x+ )∇x x+ + ρA⊤ A∇x x+ + ∇2 h1 (x+ )∇x x+ − ∇2 h1 (x) = 0 η η 22

(43)

 −1   1 2 1 2 2 + ⊤ + =⇒∇x x = ∇ f1 (x ) + ρA A + ∇ h1 (x ) ∇ h1 (x) ; η η 1 ∇2 f1 (x+ )∇y x+ + ρA⊤ B + ρA⊤ A∇y x+ + ∇2 h1 (x+ )∇y x+ = 0 η  −1
1 2 + 2 + ⊤ + =⇒∇y x = − ∇ f1 (x ) + ρA A + ∇ h1 (x ) ρA⊤ B ; η 1 ∇2 f1 (x+ )∇λ x+ + A⊤ + ρA⊤ A∇λ x+ + ∇2 h1 (x+ )∇λ x+ = 0 η  −1 1 2 + 2 + ⊤ + =⇒∇λ x = − ∇ f1 (x ) + ρA A + ∇ h1 (x ) A⊤ . η +

Therefore,  ∇x x +  Dg1 (x, y, λ) =

∇y x + Im

∇λ x +

 

Ip   =

∇2 f1 (x+ ) + ρA⊤ A + η1 ∇2 h1 (x+ )

−1



1 2 η ∇ h1 (x)

 

Im

−ρA⊤ B Im

Ip

−A⊤

  (44)

Ip

where ∇2 f1 (x+ ) + ρA⊤ A + η1 ∇2 h1 (x+ ) is positive definite by relative smoothness, strong convexity of h1 , and η < 1/L1 . Finally, by the chain rule, we get Dg(x, y, λ) = Dg3 (g2 (x, y, λ))Dg2 (g1 (x, y, λ))Dg1 (x, y, λ). Showing that det(Dg) is nonzero globally.

(45)

Because

Dg(x, y, λ) = Dg3 (g2 (x, y, λ))Dg2 (g1 (x, y, λ))Dg1 (x, y, λ), with Dg1 , Dg2 and Dg3 being square matrices, it suffices to show the global nonsingularity of all Dg1 , Dg2 and Dg3 . We begin with Dg1 . Fix (x, y, λ) ∈ Ω and suppose (α, β, γ) ∈ Rn ×Rm ×Rp satisfies Dg1 (x, y, λ)[α⊤ β ⊤ γ ⊤ ]⊤ = 0. We show α = β = γ = 0. Observe that        α ∇ x x + ∇y x + ∇λ x + α 0   β  = 0 =⇒ β = 0, γ = 0, ∇x x+ α = 0 =⇒ α = 0 Im Dg1 (x, y, λ)  β  = 0 ⇐⇒  γ Ip γ 0 where the last step uses the nonsingularity, for every state in Ω, of ∇x x + =

 −1   1 2 1 ∇ h1 (x) , ∇2 f1 (x+ ) + ρA⊤ A + ∇2 h1 (x+ ) η η

which follows from relative smoothness: the first factor is nonsingular by the positive definiteness established above, and η1 ∇2 h1 (x) ≻ 0 by strong convexity of h1 , so the product is nonsingular. (Rows 2–3 give β = 0 and γ = 0; substituting into row 1 yields ∇x x+ α = 0, hence α = 0.) Next consider Dg2 . Fix (x, y, λ) ∈ Ω and suppose Dg2 (x, y, λ)[α⊤ β ⊤ γ ⊤ ]⊤ = 0. We show again that α = β = γ = 0. Observe that        α In α 0 Dg2 (x, y, λ)  β  = 0 ⇐⇒ ∇x y+ ∇y y+ ∇λ y+   β  = 0 =⇒ α = 0, γ = 0, ∇y y+ β = 0 =⇒ β = 0 γ Ip γ 0

23

where the last step uses the nonsingularity, for every state in Ω, of −1 ∇y y = ∇2 f2 (y+ ) + ρB⊤ B + η1 ∇2 h2 (y+ ) +





 1 2 ∇ h2 (y) , η

which follows from relative smoothness: the first factor is nonsingular by positive definiteness, and η1 ∇2 h2 (y) ≻ 0 by strong convexity of h2 , so the product is nonsingular. (Rows 1 and 3 give α = 0 and γ = 0; row 2 then gives ∇y y+ β = 0, hence β = 0.) Finally, for Dg3 , fix (x, y, λ) ∈ Ω and suppose Dg3 (x, y, λ)[α⊤ β ⊤ γ ⊤ ]⊤ = 0. Then            In α α α 0 0   β  = 0 =⇒  β  = 0 Im Dg3 (x, y, λ)  β  = 0 ⇐⇒  γ γ 0 0 ∇ x λ + ∇ y λ+ ∇ λ λ+ ∇ λ λ+ γ Finally, rows 1–2 give α = 0 and β = 0, and the third row gives ∇λ λ+ γ = 0; since ∇λ λ+ = Ip by (39), we obtain γ = 0. Showing any strict-saddle KKT point in Ω lies in the unstable set. Let (x⋆ , y⋆ , λ⋆ ) ∈ Ω be a strict-saddle KKT point of Problem (1). We first show that it is a fixed point of g = g3 ◦ g2 ◦ g1 , i.e., g(x⋆ , y⋆ , λ⋆ ) = (x⋆ , y⋆ , λ⋆ ). By the KKT conditions and strict-saddle condition for this point,  ⋆ ⋆  = b, Ax + By ⋆ ∇f1 (x ) + A⊤ λ⋆ = 0,   ∇f2 (y⋆ ) + B⊤ λ⋆ = 0, 2 ⋆ ⊤ 2 ⋆ d⊤ x ∇ f1 (x )dx + dy ∇ f2 (y )dy < 0

(46)

for some (dx , dy ) : Adx + Bdy = 0.

(47)

Note that the Hessian ∇2 F (z⋆ ) is block-diagonal because f (x, y) = f1 (x) + f2 (y) is separable, so the strict-saddle condition involves no cross terms ∇2xy f . Second, (x+ , y+ , λ+ ) = (x⋆ , y⋆ , λ⋆ ) and (x, y, λ) = (x⋆ , y⋆ , λ⋆ ) satisfy the optimality conditions (43), (41) and (39) (which define the mappings g1 , g2 , g3 , respectively). Therefore gj (x⋆ , y⋆ , λ⋆ ) = (x⋆ , y⋆ , λ⋆ ) for j = 1, 2, 3, which follows from the well-posedness (unique minimizers) of the subproblems on the primal domain. It remains to show that the Jacobian matrix Dg(x⋆ , y⋆ , λ⋆ ) has an eigenvalue with magnitude greater than 1. To simplify notation, we denote H1 = ∇2 h1 (x⋆ ),

H2 = ∇2 h2 (y⋆ ),

F1 = ∇2 f1 (x⋆ ),

F2 = ∇2 f2 (y⋆ ).

We compute the Jacobian matrix Dg(x⋆ , y⋆ , λ⋆ ) by plugging (x, y, λ) = (x+ , y+ , λ+ ) = (x⋆ , y⋆ , λ⋆ ) into (45): Dg(x⋆ , y⋆ , λ⋆ ) =Dg3 (x⋆ , y⋆ , λ⋆ )Dg2 (x⋆ , y⋆ , λ⋆ )Dg1 (x⋆ , y⋆ , λ⋆ )    In In In  −ρB⊤ A Im  (F2 + ρB⊤ B + η1 H2 )−1 Im = ρA ρB Ip Ip  1  1 ⊤ −1 ⊤ ⊤ (F1 + ρA A + η H1 ) η H1 −ρA B −A    Im Im Ip Ip    In In In   −ρB⊤ A Im (F2 + ρB⊤ B + η1 H2 )−1 Im = ρA ρB Ip Ip 24

 In −B⊤   Ip

 −B⊤  Ip

 1 η H2

 Ip

 (F1 + ρA⊤ A + η1 H1 )−1 

 In 

Im

1 1 η H2



−ρA⊤ B Im

−A⊤

 In 

 

Ip 

Ip 

η H1

Ip 

In In  −ρB⊤ A Im −B⊤  (F2 + ρB⊤ B + η1 H2 )−1 Im = ρA ρB Ip Ip Ip 1  ⊤ ⊤ 1 −1 ⊤ (F1 + ρA A + η H1 ) η H1 −ρA B −A 1     Im η H2 Ip Ip  −1  1 1 ⊤ ⊤ ⊤ F1 + η H1 + ρA A η H1 −ρA B −A 1  F2 + η1 H2 B⊤   = η H2 Ip −ρA −ρB Ip     −1 1 F1 + η H1 + ρA⊤ A F1 + ρA⊤ A ρA⊤ B A⊤ F2 + η1 H2 B⊤   F2 B⊤  =I −  −ρA −ρB 0 −ρA −ρB Ip . =I − Φ, where the third equality from the end follows from Lemma A.1, and in the last line we define  F1 + η1 H1 + ρA⊤ A .  Φ= −ρA

−1  F1 + ρA⊤ A ⊤  B −ρA Ip

F2 + η1 H2 −ρB

ρA⊤ B F2 −ρB

 A⊤ B⊤  0

(48)

Lemma A.1. If F2 + ρB⊤ B + η1 H2 is positive definite (which holds on the primal domain by relative smoothness, 



In

strong convexity of h2 , and η < 1/L2 ), then  −ρA

F2 + η1 H2 −ρB

B⊤  is nonsingular on the primal domain, and its Ip

inverse is given by 

−1

In

 −ρA

F2 + η1 H2 −ρB

B⊤  Ip



 In 

In

= ρA

Im ρB



In  −ρB⊤ A

(F2 + ρB⊤ B + η1 H2 )−1

Ip

Ip

 Im

−B⊤  Ip

Proof of Lemma A.1. Multiplying the four factors below shows that their product equals the identity matrix:      In In In In    −ρB⊤ A Im −B⊤  = I. F2 + η1 H2 B⊤   (F2 + ρB⊤ B + η1 H2 )−1 Im ρA ρB Ip Ip −ρA −ρB Ip Ip 



In

Therefore 

F2 + η1 H2

B⊤ 

−ρB

Ip

−ρA

is nonsingular as long as F2 + ρB⊤ B + η1 H2 is nonsingular.

Second, we reduce the problem of showing that Dg(x⋆ , y⋆ , λ⋆ ) has an eigenvalue of magnitude greater than 1 to showing that Φ (see (48)) has a real negative eigenvalue, or equivalently, that det(Φ + µI) = 0 for some µ > 0. The following chain of equivalences reduces this singularity condition to that of a real-symmetric matrix J(µ); the key steps are Schur complementation (to eliminate the dual block) and a diagonal similarity transform (to symmetrize the off-diagonal blocks). Using det(UV) = det(U) det(V) and det(U−1 ) = det(U)−1 , and multiplying by the nonzero  F1 + η1 H1 + ρA⊤ A factor det   −ρA

−1 F2 + η1 H2 −ρB

 B⊤  Ip

, we obtain the following equivalent determinant conditions:

det(Φ + µI) = 0 25

  −1   F1 + η1 H1 + ρA⊤ A F1 + ρA⊤ A ρA⊤ B A⊤   F2 + η1 H2 B⊤   F2 B⊤  + µI = 0 ⇔ det  −ρA −ρB 0 −ρA −ρB Ip     F1 + η1 H1 + ρA⊤ A F1 + ρA⊤ A ρA⊤ B A⊤ F2 B⊤  + µ  F2 + η1 H2 B⊤  = 0 ⇔ det  −ρA −ρB 0 −ρA −ρB Ip   µ ⊤ ⊤ ⊤ (1 + µ)F1 + (1 + µ)ρA A + η H1 ρA B A (1 + µ)F2 + µη H2 (1 + µ)B⊤  = 0 ⇔ det  −(1 + µ)ρA −(1 + µ)ρB µIp   µ ⊤ ⊤ (1 + µ)F1 + (1 + µ)ρA A + η H1 ρA B A⊤ (1 + µ)ρB⊤ A (1 + µ)F2 + (1 + µ)ρB⊤ B + µη H2 B⊤  = 0 ⇔ det  −(1 + µ)ρA −(1 + µ)ρB µIp    ⊤ µ ⊤ ⊤ (1 + µ)F1 + (1 + µ)ρA A + η H1 ρA B 1 A A ⇔ det + (1 + )ρ (1 + µ)ρB⊤ A (1 + µ)F2 + (1 + µ)ρB⊤ B + µη H2 B⊤ A µ   (2 + µ1 )ρA⊤ B (1 + µ)F1 + (2 + µ + µ1 )ρA⊤ A + µη H1 ⇔ det =0 (2 + µ + µ1 )ρB⊤ A (1 + µ)F2 + (2 + µ + µ1 )ρB⊤ B + µη H2   q (2 + µ1 )(2 + µ + µ1 )ρA⊤ B (1 + µ)F1 + (2 + µ + µ1 )ρA⊤ A + µη H1  = 0, q ⇔ det  (1 + µ)F2 + (2 + µ + µ1 )ρB⊤ B + µη H2 (2 + µ1 )(2 + µ + µ1 )ρB⊤ A

 A⊤ B =0 B⊤ B

The first two equivalences use Φ + µI = M−1 (N + µM) (from Φ = M−1 N), so det(Φ + µI) = det(M−1 ) det(N + µM). Here M is nonsingular: taking the Schur complement with respect to its bottom-right block Ip yields a blocklower-triangular matrix with positive definite diagonal blocks F1 + η1 H1 + ρA⊤ A and F2 + η1 H2 + ρB⊤ B, so det(M−1 ) ̸= 0. The third adds N and µM blockwise. The fourth applies the block row operation R2 ← R2 − B⊤ R3 (where Ri denotes the i-th block row), which preserves the determinant. Since R3 = [−(1+µ)ρA, −(1+µ)ρB, µIp ], subtracting B⊤ R3 from R2 = [0, (1 + µ)F2 + µη H2 , (1 + µ)B⊤ ] gives [(1 + µ)ρB⊤ A, (1 + µ)F2 + (1 + µ)ρB⊤ B + µη H2 , (1 + µ)B⊤ − µB⊤ ] =[(1 + µ)ρB⊤ A, (1 + µ)F2 + (1 + µ)ρB⊤ B + µη H2 , B⊤ ]. The fifth applies the Schur complement with respect to µIp . Writing the 3 × 3 matrix after the fourth step as h ⊤i with Q = A and R⊤ = −(1 + µ)ρ[A B], the Schur formula gives B⊤  P det ⊤ R

 P

Q  R⊤ µIp

   Q = µp det P − µ1 QR⊤ . µIp

Since µp ̸= 0, the condition reduces to det(P − µ1 QR⊤ ) = 0, where   A⊤ A − µ1 QR⊤ = 1 + µ1 ρ B⊤ A

 A⊤ B . B⊤ B

The sixth collects like terms: (1 + µ) + (1 + µ1 ) = 2 + µ + µ1 on the diagonal ρA⊤ A, ρB⊤ B blocks and 1 + (1 + µ1 ) = 2 + µ1 on the off-diagonal ρA⊤ B block; call the resulting 2 × 2 block matrix C(µ). Its (1, 2) and (2, 1) blocks carry different coefficients 2 + µ1 and 2 + µ + µ1 , so C(µ) is not yet symmetric. The last step symmetrizes C(µ): Setting q D = diag(In , αIm ) with α = (2 + µ + µ1 )/(2 + µ1 ), the similarity D−1 C(µ)D scales the (1, 2) block by α and q the (2, 1) block by α−1 , equalizing both to (2 + µ1 )(2 + µ + µ1 ) ρA⊤ B, yielding the symmetric matrix J(µ); and det(C(µ)) = 0 ⇔ det(J(µ)) = 0 since det(D) ̸= 0.

26

Therefore, the problem is now reduced to showing J(µ) is singular for some µ > 0, where  q  µ 1 ⊤ (2 + µ1 )(2 + µ + µ1 )ρA⊤ B 1 + (2 + µ + µ )ρA A + η H1 . (1 + µ)F  q J(µ) = (2 + µ1 )(2 + µ + µ1 )ρB⊤ A (1 + µ)F2 + (2 + µ + µ1 )ρB⊤ B + µη H2 The matrix J(µ) is real-symmetric and varies continuously with µ [33, Theorem 5.1], so its eigenvalues are real and continuous. It suffices to prove that λmin (J(µ)) crosses zero for some µ > 0, since λmin (J(µ⋆ )) = 0 implies det(J(µ⋆ )) = 0. Remark A.1 (Intuition for the µ-dependent test direction). In the single-block ALM (Lemma 5.2) one can take a fixed z ∈ null(A): the penalty ρA⊤ A vanishes on null(A) independently of µ. With two blocks the cross-coupling ρA⊤ B complicates matters. Before symmetrization, the reduced matrix C(µ) has asymmetric off-diagonal coefficients 2 + µ1 and 2 + µ + µ1 . The diagonal similarity transform equalizes these coefficients in the symmetric matrix . J(µ), but it also changes the natural strict-saddle test direction to the scaled direction (dx , s(µ)dy ), where s(µ) = q (2 + µ + µ1 )/(2 + µ1 ). This scaled direction lets the KKT constraint Adx + Bdy = 0 eliminate the leading cross

term; after expansion, the remaining penalty contribution is a single non-negative term ρ · (positive factor) · µ∥Bdy ∥22 that vanishes as µ → 0+ . Since s(µ) → 1 as µ → 0+ , the test direction reduces to the strict-saddle pair (dx , dy ) and ⊤ the quadratic form tends to d⊤ x F1 dx + dy F2 dy < 0. . Lemma A.2. Let y(µ) = (dx , s(µ)dy ) with dx , dy defined in KKT condition (47) and s(µ) =

q

1 2+µ+ µ q . 1 2+ µ

Then

λmin (J(µ⋆ )) = 0 for some µ⋆ > 0, hence det(J(µ⋆ )) = 0. . Proof of Lemma A.2. Define ϕ(µ) = y(µ)⊤ J(µ)y(µ). Expanding the quadratic form using the definition of J(µ) and y(µ) = (dx , s(µ)dy ), and using the KKT constraint Adx + Bdy = 0 to eliminate the cross term, yields 1 µ ⊤ ⊤ 2 ⊤ (d H1 dx + s(µ)2 d⊤ )(s(µ)2 − 1)∥Bdy ∥22 y H2 dy ) + (1 + µ)(dx F1 dx + s(µ) dy F2 dy ) + ρ(2 + µ + η x µ 2 + µ + µ1 µ 2 ⊤ ⊤ 2 ⊤ = (d⊤ µ∥Bdy ∥22 x H1 dx + s(µ) dy H2 dy ) + (1 + µ)(dx F1 dx + s(µ) dy F2 dy ) + ρ η 2 + µ1

ϕ(µ) =

=

2µ + µ2 + 1 µ ⊤ ⊤ 2 ⊤ (dx H1 dx + s(µ)2 d⊤ µ∥Bdy ∥22 y H2 dy ) + (1 + µ)(dx F1 dx + s(µ) dy F2 dy ) + ρ η 2µ + 1

where in the second line we used s(µ)2 − 1 = 2+µ 1 (from the definition of s(µ)), and in the last line we rewrote µ

1 2+µ+ µ 1 2+ µ

=

2µ+µ2 +1 2µ+1 .

We first use ϕ(µ) to show that λmin (J(µ)) < 0 for sufficiently small µ > 0. We then prove

separately that J(µ) ≻ 0 for sufficiently large µ. The intermediate value theorem applied to λmin (J(µ)) gives the desired zero crossing. • First, since limµ→0+ s(µ) = 1, we have ⊤ lim ϕ(µ) = d⊤ x F1 dx + dy F2 dy < 0,

µ→0+

by KKT condition (47). 2

• Second, we show that ϕ(N ) > 0 for some sufficiently large N . If dy ̸= 0, then since limµ→∞ s(µ) = 21 , we have µ ϕ(µ) 1 ⊤ 1 ρ 1 ⊤ 2 = dy H2 dy + d⊤ y F2 dy + ∥Bdy ∥2 = dy 2 µ→∞ µ 2η 2 2 2 lim



 1 H2 + F2 + ρB⊤ B dy > 0, η

where the last inequality uses η1 H2 + F2 ≻ 0 (from relative smoothness, strong convexity of h2 , and η < 1/L2 ) and B⊤ B ⪰ 0. Therefore ϕ(N ) > 0 for some sufficiently large N . If dy = 0, then dx ̸= 0 (since (dx , dy ) ̸= (0, 0)),

27

2 ⊤ 2 and the KKT condition (47) gives Adx = 0. In this case s(µ)2 d⊤ y H2 dy = 0, s(µ) dy F2 dy = 0, and ∥Bdy ∥2 = 0, so   ϕ(µ) ⊤ 1 = dx H1 + F1 dx > 0, lim µ→∞ µ η

by the positive definiteness of η1 H1 +F1 +ρA⊤ A (from relative smoothness, strong convexity of h1 , and η < 1/L1 ) and Adx = 0. Therefore ϕ(N ) > 0 for some sufficiently large N . • Finally, we combine the two cases to conclude that J(µ) is singular for some µ⋆ > 0. From the first item we have, by continuity of ϕ on (0, ∞), that ϕ(µ) < 0 for all sufficiently small µ > 0, which gives y(µ)⊤ J(µ)y(µ) < 0 and therefore λmin (J(µ)) < 0. From the second item we have ϕ(N ) > 0 for some sufficiently large N > 0. In addition, we show λmin (J(N )) > 0 for large N via a large-µ positive-definiteness argument: dividing J(µ) by √ µ, the off-diagonal coefficient is O( µ) and therefore vanishes after scaling by 1/µ, while the diagonal blocks converge. Define   F1 + ρA⊤ A + η1 H1 0 J(µ) . J̄ = lim = ≻ 0, 0 F2 + ρB⊤ B + η1 H2 µ→∞ µ where the positive definiteness follows from relative smoothness (5), strong convexity of hi (i = 1, 2), and η < 1/Li . By continuity of eigenvalues for symmetric matrices, λmin (J(µ)/µ) → λmin (J̄) > 0, so λmin (J(N )) > 0 for some sufficiently large N . Since J(µ) is a real-symmetric matrix that varies continuously in µ, its eigenvalues are continuous functions of µ (see, e.g., [33, Theorem 5.1]). Therefore, by the intermediate value theorem applied to λmin (J(µ)), there exists µ⋆ ∈ (0, N ) such that λmin (J(µ⋆ )) = 0, i.e., J(µ⋆ ) is singular.

This completes the proof that strict-saddle KKT points in Ω are unstable fixed points of the fixed-point map of Algorithm 2. By Theorem 5.1, the set of initial points in Ω from which the iteration converges to such a strict-saddle KKT point has Lebesgue measure zero; since random initialization is absolutely continuous with respect to Lebesgue measure, the probability that the iterates converge to such a strict-saddle KKT point is zero.

B

Proof of Theorem 4.2

Proof of Theorem 4.2. We apply Theorem 5.1 to the fixed-point map g of Algorithm 3, verifying det(Dg) ̸= 0 everywhere and that every strict-saddle KKT point in Ω is an unstable fixed point. The consensus formulation couples local primal blocks through agreement constraints, yielding a Jacobian with large but structured block form. The argument parallels the two-block case: positive definiteness of the local Bregmanproximal Hessians gives nonsingularity, and the consensus coupling cannot cancel the negative curvature at a strict saddle, so the linearized map retains an expanding direction. Constructing the fixed-point map g.

Define the fixed-point map of Algorithm 3 as

+ + + + + (x+ 1 , y1 , · · · , xJ , yJ , λ2 , · · · , λJ ) = g(x1 , y1 , · · · , xJ , yJ , λ2 , · · · , λJ ).

(49)

We first write the Jacobian and instability calculation explicitly for the representative case J = 3, which is the smallest case showing the star coupling between the hub and multiple non-hub agents. The same block argument extends to arbitrary J ≥ 2; the case J = 2 is the simpler one-non-hub specialization. For a general number of agents, write z = (x1 , y1 , x2 , y2 , . . . , xJ , yJ , λ2 , . . . , λJ ). The star-consensus constraints are xj − x1 = 0 for j ≥ 2, with multipliers λj . Each constraint couples only the pair (x1 , xj ) through λj , so different non-hub agents do not interact directly: the block (xj , yj , λj ) interacts with (x1 , y1 ), but there is no coupling between (xj , yj , λj ) and (xℓ , yℓ , λℓ ) for j ̸= ℓ, j, ℓ ≥ 2. Accordingly, Dg(z) has the same star-shaped block sparsity as in the J = 3 case. All nonsingularity calculations in the J = 3 case detailed next are agent-wise and use only strong convexity and the implicit function theorem for each local (xj , yj )-subproblem; the same calculation is repeated for j = 2, . . . , J. The instability construction is also unchanged: the Schur-reduced 28

matrix has the same replicated j-block structure, and the general-J test direction is obtained by repeating the same x-component across agents (see Lemma B.2 below for the explicit coupling-matrix form used in that step). + + + + + + + (x+ 1 , y1 , x2 , y2 , x3 , y3 , λ2 , λ3 ) = g(x1 , y1 , x2 , y2 , x3 , y3 , λ2 , λ3 ).

Computing the Jacobian Dg at a fixed point of Algorithm 3. inequality (Definition 2.7)

(50)

Relative bi-smoothness is the positive-semidefinite

Lxj ∇2xx hj (x, y) ± ∇2xx fj (x, y) ⪰ 0,

Lyj ∇2yy hj (x, y) ± ∇2yy fj (x, y) ⪰ 0,

∇2xx fj (x, y) ⪰ −Lxj ∇2xx hj (x, y),

∇2yy fj (x, y) ⪰ −Lyj ∇2yy hj (x, y)

in particular, By assumption, each bi-Bregman kernel hj is strongly bi-convex, so ∇2xx hj (x, y) ≻ 0 and ∇2yy hj (x, y) ≻ 0 on the primal domain. Hence for η < min(1/Lxj , 1/Lyj ), 1 ∇2yy fj (x, y) + ∇2yy hj (x, y) ≻ 0. η

1 ∇2xx fj (x, y) + ∇2xx hj (x, y) ≻ 0, η

Adding the consensus penalty terms preserves positive definiteness, so each block subproblem in Algorithm 3 is strongly convex with a unique minimizer. For the global nonsingularity argument below, it is therefore enough to differentiate each elementary update map only with respect to its updated variable. For the instability argument, we only need the Jacobian at a fixed point. At such a point, differentiating the bi-Bregman terms produces no leftover + mixed-anchor terms, because the anchor and updated variables coincide (x+ j = xj and yj = yj ). Thus the fixed-point Jacobian of g satisfies  + + + + + + + + ∂x1

∂x1  ∂y  1+  ∂x1  ∂x+  2  ∂x+1  ∂y2  ∂x 1 Dg(x1 , y1 , x2 , y2 , x3 , y3 , λ2 , λ3 ) =   ∂x+  ∂x31  +  ∂y3  ∂x1  ∂λ+  2  ∂x+1 ∂λ3 ∂x1

∂x1 ∂y1 ∂y1+ ∂y1 ∂x+ 2 ∂y1 ∂y2+ ∂y1 ∂x+ 3 ∂y1 ∂y3+ ∂y1 ∂λ+ 2 ∂y1 ∂λ+ 3 ∂y1

∂x1 ∂x2 ∂y1+ ∂x2 ∂x+ 2 ∂x2 ∂y2+ ∂x2 ∂x+ 3 ∂x2 ∂y3+ ∂x2 ∂λ+ 2 ∂x2 ∂λ+ 3 ∂x2

∂x1 ∂y2 ∂y1+ ∂y2 ∂x+ 2 ∂y2 ∂y2+ ∂y2 ∂x+ 3 ∂y2 ∂y3+ ∂y2 ∂λ+ 2 ∂y2 ∂λ+ 3 ∂y2

∂x1 ∂x3 ∂y1+ ∂x3 ∂x+ 2 ∂x3 ∂y2+ ∂x3 ∂x+ 3 ∂x3 ∂y3+ ∂x3 ∂λ+ 2 ∂x3 ∂λ+ 3 ∂x3

∂x1 ∂y3 ∂y1+ ∂y3 ∂x+ 2 ∂y3 ∂y2+ ∂y3 ∂x+ 3 ∂y3 ∂y3+ ∂y3 ∂λ+ 2 ∂y3 ∂λ+ 3 ∂y3

∂x1 ∂λ2 ∂y1+ ∂λ2 ∂x+ 2 ∂λ2 ∂y2+ ∂λ2 ∂x+ 3 ∂λ2 ∂y3+ ∂λ2 ∂λ+ 2 ∂λ2 ∂λ+ 3 ∂λ2

∂x1 ∂λ3  ∂y1+  ∂λ3   ∂x+ 2  ∂λ3  ∂y2+   ∂λ3  . ∂x+ 3 ∂λ3   ∂y3+  ∂λ3   ∂λ+ 2  ∂λ3  ∂λ+ 3 ∂λ3

(51)

At a fixed point, the optimality conditions of (13) and L(·) in (12) read: 1 + + ∇x f1 (x+ 1 , y1 ) − λ2 − λ3 + ρ(2x1 − x2 − x3 ) + (∇x h1 (x1 , y1 ) − ∇x h1 (x1 , y1 )) = 0 η 1 + + + + ∇y f1 (x+ 1 , y1 ) + (∇y h1 (x1 , y1 ) − ∇y h1 (x1 , y1 )) = 0 η 1 + + + ∇x f2 (x+ 2 , y2 ) + λ2 + ρ(x2 − x1 ) + (∇x h2 (x2 , y2 ) − ∇x h2 (x2 , y2 )) = 0 η 1 + + + + ∇y f2 (x+ 2 , y2 ) + (∇y h2 (x2 , y2 ) − ∇y h2 (x2 , y2 )) = 0 η 1 + + + ∇x f3 (x+ 3 , y3 ) + λ3 + ρ(x3 − x1 ) + (∇x h3 (x3 , y3 ) − ∇x h3 (x3 , y3 )) = 0 η 1 + + + + ∇y f3 (x+ 3 , y3 ) + (∇y h3 (x3 , y3 ) − ∇y h3 (x3 , y3 )) = 0 η + + −λ+ 2 + λ2 + ρ(x2 − x1 ) = 0

(58)

+ + −λ+ 3 + λ3 + ρ(x3 − x1 ) = 0

(59)

29

(52) (53) (54) (55) (56) (57)

Before applying the implicit function theorem, to simplify notations, we define 1 2 . + Γ1xx = ∇2xx f1 (x+ 1 , y1 ) + 2ρI + ∇xx h1 (x1 , y1 ), η 1 2 . + + + Γ1yy = ∇2yy f1 (x+ 1 , y1 ) + ∇yy h1 (x1 , y1 ), η 1 2 . + Γ2xx = ∇2xx f2 (x+ 2 , y2 ) + ρI + ∇xx h2 (x2 , y2 ), η 1 2 . + + + Γ2yy = ∇2yy f2 (x+ 2 , y2 ) + ∇yy h2 (x2 , y2 ), η 1 2 . + Γ3xx = ∇2xx f3 (x+ 3 , y3 ) + ρI + ∇xx h3 (x3 , y3 ), η 1 2 . + + + Γ3yy = ∇2yy f3 (x+ 3 , y3 ) + ∇yy h3 (x3 , y3 ), η . Fjxx = ∇2xx fj (x+ j , yj ),

. F1xy = ∇2xy f1 (x+ 1 , y1 ),

. H1xx = ∇2xx h1 (x1 , y1 ),

. + F1yx = ∇2yx f1 (x+ 1 , y1 ),

. H1yy = ∇2yy h1 (x+ 1 , y1 ),

. F2xy = ∇2xy f2 (x+ 2 , y2 ),

. H2xx = ∇2xx h2 (x2 , y2 ),

. + F2yx = ∇2yx f2 (x+ 2 , y2 ),

. H2yy = ∇2yy h2 (x+ 2 , y2 ),

. F3xy = ∇2xy f3 (x+ 3 , y3 ),

. H3xx = ∇2xx h3 (x3 , y3 ),

. + F3yx = ∇2yx f3 (x+ 3 , y3 ),

. H3yy = ∇2yy h3 (x+ 3 , y3 ).

. + Fjyy = ∇2yy fj (x+ j , yj ),

∀j = 1, 2, 3.

(60)

Note that Hjxx = ∇2xx hj (xj , yj ) is evaluated at the old iterate (xj , yj ), whereas Hjyy = ∇2yy hj (x+ j , yj ) is evaluated + + at (xj , yj ), i.e., with the already-updated xj . This asymmetry reflects the Gauss–Seidel ordering in Algorithm 3: the xj -subproblem uses the bi-Bregman divergence Dhxj (xj , xk−1 ; yjk−1 ) whose anchor point is (xk−1 , yjk−1 ) = (xj , yj ), j j while the subsequent yj -subproblem uses Dhy j (yj , yjk−1 ; xkj ) whose anchor point is (xkj , yjk−1 ) = (x+ j , yj ). Apply implicit function theorem to (52): Γ1xx

∂x+ 1 ∂x1

=

1 1 H η xx

Apply implicit function theorem to (53):

∂x+ Γ1xx 1 = −F1xy ∂y1

Apply implicit function theorem to (54): −ρ

∂x+ ∂x+ 1 + Γ2xx 2 = 0 ∂x1 ∂x1

H1yy ∂x+ ∂y+ F1yx 1 + Γ1yy 1 = ∂y1 ∂y1 η

−ρ

∂x+ ∂x+ 1 + Γ2xx 2 = 0 ∂y1 ∂y1

F1yx

∂x+ 1 ∂x1

+ Γ1yy

∂y1+ ∂x1

=0

Γ1xx

∂x+ 1 = ρI ∂x2

F1yx

∂x+ ∂y+ 1 + Γ1yy 1 = 0 ∂x2 ∂x2

−ρ

∂x+ ∂x+ H2 1 + Γ2xx 2 = xx ∂x2 ∂x2 η

Γ1xx

∂x+ 1 =0 ∂y2

F1yx

∂x+ ∂y+ 1 + Γ1yy 1 = 0 ∂y2 ∂y2

−ρ

∂x+ ∂x+ 1 + Γ2xx 2 = −F2xy ∂y2 ∂y2

Γ1xx

∂x+ 1 = ρI ∂x3

F1yx

∂x+ ∂y+ 1 + Γ1yy 1 = 0 ∂x3 ∂x3

−ρ

∂x+ ∂x+ 1 + Γ2xx 2 = 0 ∂x3 ∂x3

Γ1xx

∂x+ 1 =0 ∂y3

F1yx

∂y+ ∂x+ 1 + Γ1yy 1 = 0 ∂y3 ∂y3

−ρ

∂x+ ∂x+ 1 + Γ2xx 2 = 0 ∂y3 ∂y3

Γ1xx

∂x+ 1 =I ∂λ2

F1yx

∂x+ ∂y+ 1 + Γ1yy 1 = 0 ∂λ2 ∂λ2

−ρ

∂x+ ∂x+ 1 + Γ2xx 2 = −I ∂λ2 ∂λ2

Γ1xx

∂x+ 1 =I ∂λ3

F1yx

∂y+ ∂x+ 1 + Γ1yy 1 = 0 ∂λ3 ∂λ3

−ρ

∂x+ ∂x+ 1 + Γ2xx 2 = 0. ∂λ3 ∂λ3

(61)

30

(62)

(63)

Apply implicit function theorem to (56):

Apply implicit function theorem to (55): F2yx

∂x+ 2 ∂x1

+ Γ2yy

∂y2+ ∂x1

F3yx

∂y+ ∂x+ 3 + Γ3yy 3 = 0 ∂x1 ∂x1

∂x+ ∂y+ F2yx 2 + Γ2yy 2 = 0 ∂y1 ∂y1

∂x+ ∂x+ −ρ 1 + Γ3xx 3 = 0 ∂y1 ∂y1

F3yx

∂x+ ∂y+ 3 + Γ3yy 3 = 0 ∂y1 ∂y1

∂x+ ∂y+ 2 + Γ2yy 2 = 0 ∂x2 ∂x2 + + H2yy ∂y ∂x F2yx 2 + Γ2yy 2 = ∂y2 ∂y2 η F2yx

−ρ

∂x+ ∂x+ 1 + Γ3xx 3 = 0 ∂x2 ∂x2

F3yx

∂x+ ∂y+ 3 + Γ3yy 3 = 0 ∂x2 ∂x2

−ρ

∂x+ ∂x+ 1 + Γ3xx 3 = 0 ∂y2 ∂y2

F3yx

∂y+ ∂x+ 3 + Γ3yy 3 = 0 ∂y2 ∂y2

F3yx

∂x1

+ Γ3xx

∂x+ 3 ∂x1

F2yx

∂x+ ∂y+ 2 + Γ2yy 2 = 0 ∂x3 ∂x3

−ρ

∂x+ ∂x+ H3 1 + Γ3xx 3 = xx ∂x3 ∂x3 η

F2yx

∂y+ ∂x+ 2 + Γ2yy 2 = 0 ∂y3 ∂y3

−ρ

∂x+ ∂x+ 1 + Γ3xx 3 = −F3xy ∂y3 ∂y3

F2yx

∂x+ ∂y+ 2 + Γ2yy 2 = 0 ∂λ2 ∂λ2

−ρ

∂x+ ∂x+ 1 + Γ3xx 3 = 0 ∂λ2 ∂λ2

F2yx

∂x+ ∂y+ 2 + Γ2yy 2 = 0. ∂λ3 ∂λ3

−ρ

∂x+ ∂x+ 1 + Γ3xx 3 = −I. ∂λ3 ∂λ3

(64)

∂x+ ∂y+ 3 + Γ3yy 3 = 0 ∂x3 ∂x3 + H3yy ∂y+ ∂x F3yx 3 + Γ3yy 3 = ∂y3 ∂y3 η

(65)

∂x+ ∂y+ 3 + Γ3yy 3 = 0 ∂λ2 ∂λ2

F3yx

∂x+ ∂y+ 3 + Γ3yy 3 = 0. ∂λ3 ∂λ3

(66)

=0

−

∂λ+ ∂x+ ∂x+ 3 +ρ 3 −ρ 1 =0 ∂x1 ∂x1 ∂x1

∂λ+ ∂x+ ∂x+ − 2 +ρ 2 −ρ 1 =0 ∂y1 ∂y1 ∂y1

−

∂λ+ ∂x+ ∂x+ 3 +ρ 3 −ρ 1 =0 ∂y1 ∂y1 ∂y1

−

∂λ+ ∂x+ ∂x+ 2 +ρ 2 −ρ 1 =0 ∂x2 ∂x2 ∂x2

−

∂λ+ ∂x+ ∂x+ 3 +ρ 3 −ρ 1 =0 ∂x2 ∂x2 ∂x2

−

∂λ+ ∂x+ ∂x+ 2 +ρ 2 −ρ 1 =0 ∂y2 ∂y2 ∂y2

−

∂λ+ ∂x+ ∂x+ 3 +ρ 3 −ρ 1 =0 ∂y2 ∂y2 ∂y2

−

∂x+ ∂x+ ∂λ+ 2 +ρ 2 −ρ 1 =0 ∂x3 ∂x3 ∂x3

−

∂x+ ∂x+ ∂λ+ 3 +ρ 3 −ρ 1 =0 ∂x3 ∂x3 ∂x3

−

∂λ+ ∂x+ ∂x+ 2 +ρ 2 −ρ 1 =0 ∂y3 ∂y3 ∂y3

−

∂λ+ ∂x+ ∂x+ 3 +ρ 3 −ρ 1 =0 ∂y3 ∂y3 ∂y3

−

∂λ+ ∂x+ ∂x+ 2 + ρ 2 − ρ 1 = −I ∂λ2 ∂λ2 ∂λ2

−

∂λ+ ∂x+ ∂x+ 3 +ρ 3 −ρ 1 =0 ∂λ2 ∂λ2 ∂λ2

−

∂λ+ ∂x+ ∂x+ 2 + ρ 2 − ρ 1 = 0. ∂λ3 ∂λ3 ∂λ3

−

∂λ+ ∂x+ ∂x+ 3 + ρ 3 − ρ 1 = −I. (68) ∂λ3 ∂λ3 ∂λ3

−

∂λ+ 2 ∂x1

+ρ

∂x1

−ρ

∂x+ 1

F3yx

Apply implicit function theorem to (59):

Apply implicit function theorem to (58): ∂x+ 2

Apply implicit function theorem to (57):

=0

−ρ

=0

∂x+ 1

∂x1

(67)

Representing (61)-(68) in matrix form, we further get ∂x+ 1 ∂x1   ∂y + 1 0   ∂x1  + 0   ∂x2 ∂x1  0   ∂y2+   0   ∂x1  ∂x+ 0   ∂x31  0   ∂y3+  ∂x1 0   ∂λ +  −I  ∂x21 ∂λ+ 3 ∂x1



Γ1xx F1yx  −ρI   0  −ρI   0  −ρI −ρI 

0 Γ1yy 0 0 0 0 0 0

0 0 Γ2xx F2yx 0 0 ρI 0

0 0 0 Γ2yy 0 0 0 0

0 0 0 0 Γ3xx F3yx 0 ρI

0 0 0 0 0 Γ3yy 0 0

0 0 0 0 0 0 −I 0

31

∂x+ 1 ∂y1 ∂y1+ ∂y1 ∂x+ 2 ∂y1 ∂y2+ ∂y1 ∂x+ 3 ∂y1 ∂y3+ ∂y1 ∂λ+ 2 ∂y1 ∂λ+ 3 ∂y1

∂x+ 1 ∂x2 ∂y1+ ∂x2 ∂x+ 2 ∂x2 ∂y2+ ∂x2 ∂x+ 3 ∂x2 ∂y3+ ∂x2 ∂λ+ 2 ∂x2 ∂λ+ 3 ∂x2

∂x+ 1 ∂y2 ∂y1+ ∂y2 ∂x+ 2 ∂y2 ∂y2+ ∂y2 ∂x+ 3 ∂y2 ∂y3+ ∂y2 ∂λ+ 2 ∂y2 ∂λ+ 3 ∂y2

∂x+ 1 ∂x3 ∂y1+ ∂x3 ∂x+ 2 ∂x3 ∂y2+ ∂x3 ∂x+ 3 ∂x3 ∂y3+ ∂x3 ∂λ+ 2 ∂x3 ∂λ+ 3 ∂x3

∂x+ 1 ∂y3 ∂y1+ ∂y3 ∂x+ 2 ∂y3 ∂y2+ ∂y3 ∂x+ 3 ∂y3 ∂y3+ ∂y3 ∂λ+ 2 ∂y3 ∂λ+ 3 ∂y3

∂x+ 1 ∂λ2 ∂y1+ ∂λ2 ∂x+ 2 ∂λ2 ∂y2+ ∂λ2 ∂x+ 3 ∂λ2 ∂y3+ ∂λ2 ∂λ+ 2 ∂λ2 ∂λ+ 3 ∂λ2

∂x+ 1 ∂λ3  ∂y1+  ∂λ3   ∂x+ 2  ∂λ3  ∂y2+   ∂λ3   ∂x+ 3 ∂λ3  + ∂y3  ∂λ3   ∂λ+ 2  ∂λ3  ∂λ+ 3 ∂λ3



 H1

xx

η

  0    0   0 =   0    0   0 0



−F1xy

ρI

0

ρI

0

I

H1yy η

0

0

0

0

H2xx η

−F2xy

0

0

 0   −I 0    0 0  0 −I   0 0  −I 0  0 −I

0 0

0

H2yy η

0

0 −F3xy 0 0

0

0

0

H3xx η

0 0 0

0 0 0

0 0 0

0 0 0

H3yy η

I

0

(69)

where Γjxx , Γjyy , Fjxy , Fjyx , Fjxx , Fjyy , Hjxx , Hjyy for j = 1, 2, 3 are defined in (60). By the relative smoothness assumption and stepsize choice of η, we have from (60) that 1 . Γ1xx = 2ρI + H1xx + F1xx ≻ 0, η

1 . Γjxx = ρI + Hjxx + Fjxx ≻ 0, ∀j ≥ 2, η

. 1 Γjyy = Hjyy + Fjyy ≻ 0, ∀j ≥ 1. η

+ Therefore, at any fixed point (where x+ j = xj and yj = yj ), we can compute the fixed-point Jacobian Dg by solving the invertible linear system (69).

Lemma B.1. For the case J = 3, suppose g is the mapping function of Algorithm 3. At any fixed point + (x⋆1 , y1⋆ , x⋆2 , y2⋆ , x⋆3 , y3⋆ , λ⋆2 , λ⋆3 ) where x+ j = xj and yj = yj for j = 1, 2, 3, the Jacobian Dg is given by  H1

xx

Γ1xx F1yx  

−ρI   0 Dg =  −ρI   0  −ρI −ρI

0 Γ1yy 0 0 0 0 0 0

0 0 Γ2xx F2yx 0 0 ρI 0

0 0 0 Γ2yy 0 0 0 0

0 0 0 0 Γ3xx F3yx 0 ρI

0 0 0 0 0 Γ3yy 0 0

−1

0 0 0 0  0 0  0 0  0 0  0 0  −I 0  0 −I

η

  0    0    0   0    0   0

ρI

0

ρI

0

I

H1yy η

0

0

0

0

H2xx η

−F2xy

0

0

 0   −I 0    0 0 .  0 −I   0 0  −I 0  0 −I (70)

0 0

0

0

H2yy η

0

0 −F3xy 0 0

0

0

0

H3xx η

0 0 0

0 0 0

0 0 0

0 0 0

H3yy η

0

Showing that det(Dg) is nonzero globally. We decompose one full iteration of Algorithm Q3 into 2J + 1 elementary maps (for J = 3: seven maps) and show each has a nonsingular Jacobian; then det(Dg) = k det(Dgk ) ̸= 0 follows by the chain rule. Define the elementary maps acting on the full state (x1 , y1 , x2 , y2 , x3 , y3 , λ2 , λ3 ): • hxj : updates xj via the Bregman-proximal step (13), holding all other state variables fixed. • hyj : updates yj via the bi-Bregman-proximal step in (13), using the already-updated x+ j , holding all other variables fixed. • hλ : dual ascent, updating λ2 , λ3 using the already-updated primal variables. Then g = hλ ◦ hy3 ◦ hx3 ◦ hy2 ◦ hx2 ◦ hy1 ◦ hx1 Nonsingularity of Dhxj . The map hxj only changes the xj component  of the state, so Dhxj is a block matrix with  identity in every row except the xj -row. Hence det(Dhxj ) = det

∂x+ j ∂xj

. From the implicit function theorem applied

to the xj optimality condition (with all other variables treated as fixed parameters), differentiating with respect to xj gives
−1 1 j  ∂x+ j = Γjxx η Hxx , ∂xj 32



−F1xy

I

where Γjxx ≻ 0 by the relative bi-smoothness assumption and the stepsize condition, and Hjxx = ∇2xx hj (xj , yj ) ≻ ∂x+

0 by the strong convexity of hj . Since ∂xjj is a product of two positive definite (hence nonsingular) matrices, det(Dhxj ) ̸= 0. Nonsingularity of Dhyj . By the same argument applied to the yj optimality condition, differentiating with respect to yj gives
−1 1 j  ∂yj+ = Γjyy η Hyy , ∂yj  + ∂y j + j 2 where Γyy ≻ 0 and Hyy = ∇yy hj (xj , yj ) ≻ 0. Hence det(Dhyj ) = det ∂yjj ̸= 0. + + Nonsingularity of Dhλ . The dual update is λ+ j = λj + ρ(xj − x1 ). Its full Jacobian is block triangular: the + primal variables are left unchanged, and the derivative of each λj with respect to its own multiplier λj is I, while derivatives with respect to the updated primal variables appear only in lower-left blocks. Hence det(Dhλ ) = 1, so Dhλ is nonsingular. Q Since all seven elementary Jacobians are nonsingular, we conclude det(Dg) = k det(Dhk ) ̸= 0 everywhere in the domain.

Strict-saddle KKT points. First, by Definition 2.2, let z⋆ ∈ Ω be a strict-saddle KKT point of the consensus problem for J = 3, with primal components (x⋆1 , y1⋆ , . . . , x⋆3 , y3⋆ ) and multipliers λ⋆2 , λ⋆3 . Then there exist dx , d1y , d2y , d3y such that   x⋆1 = x⋆2 = x⋆3 ,    ⋆ ⋆ ⋆ ⋆   ∇x f1 (x1 , y1 ) − λ2 − λ3 = 0, (71) ∇x f2 (x⋆2 , y2⋆ ) + λ⋆2 = 0,   ∇x f3 (x⋆3 , y3⋆ ) + λ⋆3 = 0,    ∇ f (x⋆ , y⋆ ) = 0, ∀j ∈ {1, 2, 3}, y j j j   3 X dx i⊤ 2 ⋆ ⋆ [d⊤ d ]∇ f (x , y ) < 0. (72) i i x y i diy

i=1 ⋆ ⋆ ⋆ ⋆ ⋆ Now, (x1 , y1 , · · · , x3 , y3 , λ2 , λ⋆3 ) is a fixed point of g since (71) satisfies the fixed point equations (52)-(59). It remains to show that the spectral radius of Dg(x⋆1 , y1⋆ , x⋆2 , y2⋆ , x⋆3 , y3⋆ , λ⋆2 , λ⋆3 ) is larger than 1.

Set ⋆ ⋆ + + + + + + + ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ (x+ 1 , y1 , x2 , y2 , x3 , y3 , λ2 , λ3 ) =(x1 , y1 , x2 , y2 , x3 , y3 , λ2 , λ3 ) = (x1 , y1 , x2 , y2 , x3 , y3 , λ2 , λ3 ),

(73)

and evaluate Dg (70) at (73) H1 xx −1  η



Γ1xx F1yx  −ρI   0 Dg =  −ρI   0  −ρI −ρI 

Γ1xx F1yx  −ρI   0 =I −  −ρI   0  −ρI −ρI 

0 Γ1yy 0 0 0 0 0 0 0 Γ1yy 0 0 0 0 0 0

0 0 Γ2xx F2yx 0 0 ρI 0 0 0 Γ2xx F2yx 0 0 ρI 0

0 0 0 Γ2yy 0 0 0 0 0 0 0 Γ2yy 0 0 0 0

0 0 0 0 Γ3xx F3yx 0 ρI 0 0 0 0 Γ3xx F3yx 0 ρI

0 0 0 0 0 Γ3yy 0 0 0 0 0 0 0 Γ3yy 0 0

0 0 0 0 0 0 −I 0 0 0 0 0 0 0 −I 0

0 0  0  0  0  0  0 −I

−F1xy

ρI

0

ρI

0

I

I

H1 yy η

0

0

0

0

0

H2 xx η

−F2xy

0

0

−I

0

0

0

−F3xy H3 yy η

0

  0   0    0    0    0   0 0

−1  1 Fxx + 2ρI 0 1  0   Fyx  −ρI 0    0 0    0   −ρI  0 0   0   −ρI −I −ρI

33

0 0

0

H2 yy η

0

0

0

H3 xx η

0 0 0

0 0 0

0 0 0

0 0 0

F1xy F1yy 0 0 0 0 0 0

−ρI 0 F2xx + ρI F2yx 0 0 ρI 0

0 0 F2xy F2yy 0 0 0 0



0 0

0 −I 0

  0  0   0   −I   0  0 −I

−ρI 0 0 0 F3xx + ρI F3yx 0 ρI

0 0 0 0 F3xy F3yy 0 0

−I 0 I 0 0 0 0 0

 −I 0  0  0  I  0  0 0

. =I − Φ.

Here, in view of (60), the second equality follows from that 1 . Γ1xx = 2ρI + H1xx + F1xx , η

1 . Γjxx = ρI + Hjxx + Fjxx , ∀j ≥ 2, η

. 1 Γjyy = Hjyy + Fjyy , ∀j ≥ 1. η

Therefore, to show that the spectral radius of Dg is larger than 1, it suffices to show that Φ has a real negative eigenvalue, which is further equivalent to det(Φ + µI) = 0

for some µ > 0.

Write Dg = I − Φ = I − M−1 N, where M is the block-lower-triangular matrix appearing with the inverse in the computation above and N the remaining factor. The first ⇔ below uses Φ + µI = M−1 N + µI = M−1 (N + µM), so det(Φ+µI) = det(M−1 ) det(N+µM). Since M is block-lower-triangular with diagonal blocks built from Γjxx and Γjyy (positive definite by relative bi-smoothness and the stepsize condition) and identity blocks from the dual update, M is nonsingular and det(M−1 ) ̸= 0, reducing the condition to det(N + µM) = 0. The second ⇔ symmetrizes N + µM via a two-sided row-column scaling: blocks corresponding to the yj -variables and the non-hub pxj -blocks √ (j ≥ 2) are scaled by s2 = 1 + µ, and those corresponding to the Lagrangian dual variables λj by s1 = (1 + µ)ρ; since all factors are positive, the zero-determinant condition is preserved. That is, det(Φ + µI) = 0 ⇐⇒  1 Fxx + 2ρI + µΓ1xx  (1 + µ)F1yx   −(1 + µ)ρI   0  det  −(1 + µ)ρI   0   −(1 + µ)ρI −(1 + µ)ρI ⇐⇒  1 Fxx + 2ρI + µΓ1xx  s2 F1yx   −s2 ρI   0  det  −s2 ρI   0   −s1 I −s1 I |

F1xy 1 Fyy + µΓ1yy 0 0 0 0 0 0 s2 F1xy F1yy + µΓ1yy 0 0 0 0 0 0

−ρI 0 F2xx + ρI + µΓ2xx (1 + µ)F2yx 0 0 (1 + µ)ρI 0

0 0 F2xy F2yy + µΓ2yy 0 0 0 0

−ρI 0 0 0 F3xx + ρI + µΓ3xx (1 + µ)F3yx 0 (1 + µ)ρI

0 0 0 0 F3xy F3yy + µΓ3yy 0 0

−I 0 I 0 0 0 −µI 0

−s2 ρI 0 F2xx + ρI + µΓ2xx s2 F2yx 0 0 s1 I 0

0 0 s2 F2xy F2yy + µΓ2yy 0 0 0 0 {z

−s2 ρI 0 0 0 F3xx + ρI + µΓ3xx s2 F3yx 0 s1 I

0 0 0 0 s2 F3xy F3yy + µΓ3yy 0 0

−s1 I 0 s1 I 0 0 0 −µI 0

J(µ)

 −I 0   0   0  =0 I   0   0  −µI  −s1 I 0   0   0  =0 s1 I   0   0  −µI }

. p . √ where s1 = (1 + µ)ρ and s2 = 1 + µ. Therefore, the problem reduces to showing J(µ) is singular for some µ > 0. The last two diagonal blocks of J(µ) are −µI, so λmin (J(µ)) → −∞ as µ → ∞ and a direct sign-change argument on λmin (J(µ)) fails; instead we Schur-reduce the dual block. Since the bottom-right diagonal blocks of J(µ) are −µI (invertible for µ > 0), partition   S(µ) T(µ) J(µ) = , T(µ)⊤ −µI where S(µ) is the top-left 6-primal-block matrix and T(µ) contains the two dual coupling columns of J(µ). Then for any µ > 0: det(J(µ)) = det(−µI) · det(J̃(µ)),

. J̃(µ) = S(µ) + µ1 T(µ)T(µ)⊤ ,

(74)

so J(µ) is singular ⇐⇒ J̃(µ) is singular. The matrix J̃(µ) is real-symmetric and varies continuously in µ for µ > 0.

34

Ordering the primal rows as (x1 , y1 , x2 , y2 , x3 , y3 ) and the dual columns as (λ2 , λ3 ), the coupling matrix T(µ) reads   −s1 I −s1 I  0 0     s1 I 0  . p   , s1 = (1 + µ)ρ, T(µ) =   0 0    0 s1 I  0 0 where column j − 1 encodes the constraint xj − x1 = 0 via −s1 I in the x1 -rows and +s1 I in the xj -rows (j = 2, 3). Remark B.1 (Intuition for the consensus test direction). For the J = 3 calculation, the test direction d = (dx , d1y , dx , d2y , dx , d3y , 0, 0) must satisfy two conditions at once: it should lie in the null space of the consensus coupling so that the ρ-penalty cancels, and it should expose negative curvature of the block-diagonal Hessian diag(∇2 f1 , ∇2 f2 , ∇2 f3 ). Both are possible because a strict-saddle KKT point of (11) is feasible and therefore satisfies x⋆1 = x⋆2 = x⋆3 , so the negative-curvature direction from (72) may be chosen with a common x-component dx across the three agents. This equal-component structure places the x-stack (dx , dx , dx ) in the kernel of the three-node star-graph Laplacian, so T(µ)⊤ dp = 0 for every µ, where dp denotes the primal part of d. Equivalently, with dx,stack = (dx , dx , dx ), the consensus penalty contributes ρ d⊤ x,stack (Lstar ⊗ In )dx,stack = 0 on 2 2 2 this equal-x direction. At µ = 0 only the three Hessian blocks ∇ f1 , ∇ f2 , ∇ f3 survive, recovering the strict-saddle negativity from (72). In contrast to the two-block ADMM case (Remark A.1), no µ-dependent rescaling of the y-components is needed: the off-diagonal ρ-terms couple only x-blocks, and setting the dual components to zero decouples the primal form from the Schur complement correction µ1 TT⊤ . By feasibility of the consensus constraints in (11), every stationary point, hence every strict saddle, satisfies x⋆1 = = x⋆3 . Therefore the negative-curvature direction from (72) can be chosen with all x-components equal to a common dx . Indeed, the constraints xj = x1 for j = 2, 3 force every feasible perturbation to satisfy δxj = δx1 , so any direction witnessing the strict-saddle condition (72) must already have this equal-x form. Set

x⋆2

d = (dx , d1y , dx , d2y , dx , d3y , 0, 0) ̸= 0, where dx , d1y , d2y , d3y satisfy (72) and the dual components are zero. Write dp = (dx , d1y , dx , d2y , dx , d3y ) for the primal part of d. Since each column of T(µ) has a −s1 I in the x1 -rows and a +s1 I in exactly one other x-block, and all x-components of d equal dx , we get   (−s1 I)dx + (s1 I)dx T(µ)⊤ dp = = 0 for all µ ≥ 0. (−s1 I)dx + (s1 I)dx . Setting ϕ(µ) = d⊤ J(µ)d and writing d = (dp , 0), we expand using the block structure of J(µ): ⊤ 2 ⊤ ϕ(µ) = d⊤ p S(µ)dp + 2 dp T(µ) · 0 − µ∥0∥ = dp S(µ)dp .

Since T(µ)⊤ dp = 0, the Schur complement formula (74) gives ⊤ 2 ⊤ ⊤ 1 d⊤ p J̃(µ)dp = dp S(µ)dp + µ ∥T(µ) dp ∥ = dp S(µ)dp ,

so ⊤ ϕ(µ) = d⊤ J(µ)d = d⊤ p S(µ)dp = dp J̃(µ)dp

for all µ > 0,

and by continuity of S(µ) also at µ = 0. We start by showing that ϕ(0) < 0. Substituting µ = 0 into J(µ) and expanding the quadratic form ϕ(0) = d⊤ J(0)d (using that the dual components of d are zero and the x-components repeat across agents) gives ϕ(0) = d⊤ J(0)d

35

⊤  1 Fxx + 2ρI dx d1y   F1yx    dx   −ρI  2  dy   0   = dx   −ρI  3  dy   0     0   −√ρI √ 0 − ρI  ⊤  1 Fxx + 2ρI dx d1y   F1yx    dx   −ρI   = d2y   0    dx   −ρI d3y 0  ⊤  1 Fxx F1xy dx d1y  F1yx F1yy    dx   0 0   = d2y   0 0    dx   0 0 d3y 0 0 

F1xy F1yy 0 0 0 0 0 0

−ρI 0 F2xx + ρI F2yx 0 0 √ ρI 0

0 0 F2xy F2yy 0 0 0 0

F1xy F1yy 0 0 0 0

−ρI 0 F2xx + ρI F2yx 0 0

0 0 F2xy F2yy 0 0

−ρI 0 0 0 F3xx + ρI F3yx 0 √ ρI

−ρI 0 0 0 F3xx + ρI F3yx   0 dx d1y  0      0   d2x    0  dy  3  dx  Fxy 3 d3y Fyy

√ − ρI 0 √ ρI 0 0 0 0 0   0 dx d1y  0      0   d2x    0  dy  3  dx  Fxy 3 d3y Fyy

0 0 0 0 F3xy F3yy 0 0

√   − ρI dx d1y  0      0   d2x  dy  0    √    ρI   d3x    0   dy  0  0  0 0

0 0 0 0 0 0 0 F2xx F2xy F2yx F2yy 0 0 0 F3xx 0 0 F3yx   3 X dx i⊤ 2 ⋆ ⋆ [d⊤ = d ]∇ f (x , y ) < 0 (by (72)). i i x y i diy i=1

Second, we will show that ϕ(µ) > 0. µ p √ To compute the limit, divide J(µ) by µ and use that s1 = (1 + µ)ρ and s2 = 1 + µ, so the off-diagonal coupling √ blocks are O( µ) and therefore vanish after scaling by 1/µ, while each diagonal primal block is of the form Fj∗∗ + const · ρI + µΓj∗∗ and hence converges to Γj∗∗ . This yields lim

µ→∞

⊤  1 Γxx dx d1y      dx    2    dy   ϕ(µ) J(µ) ⊤   lim =d lim d= dx   µ→∞ µ µ→∞ µ  3  dy      0  0  ⊤  1 Γxx dx d1y      dx     = d2y      dx   d3y 



Γ1yy

Γ1yy

Γ2xx

Γ2xx

Γ2yy

Γ2yy

Γ3xx

Γ3xx

 dx  d1y     dx   2  dy     dx   3 3  dy  Γyy    0  −I 0 −I   dx  d1y     dx    2  > 0.  dy     dx  3 d3y Γyy

Here, in view of (60), the last line is because 1 . Γ1xx = 2ρI + H1xx + F1xx ≻ 0, η

1 . Γjxx = ρI + Hjxx + Fjxx ≻ 0, ∀j ≥ 2, η

. 1 Γjyy = Hjyy + Fjyy ≻ 0, ∀j ≥ 1, η

which follows from relative bi-smoothness, strong bi-convexity of the kernels, and η < minj (1/Lxj , 1/Lyj ). Since ϕ(µ) is a continuous real function of µ, the above limit implies ϕ(N ) > 0 for some sufficiently large N . 36

⊤ Although J̃(µ) is defined only for µ > 0, the identity above gives d⊤ p J̃(µ)dp = dp S(µ)dp along this test direction, and the right-hand side has a continuous extension to µ = 0. Therefore the two cases above translate to statements about λmin (J̃(µ)): 2 • Small µ: ϕ(0) < 0 and continuity give ϕ(µ) < 0 for small µ > 0, so λmin (J̃(µ)) ≤ d⊤ p J̃(µ)dp /∥dp ∥ < 0. p √ √ • Large µ: Since s1 = (1 + µ)ρ and s2 = 1 + µ, both coefficients are O( µ). The entries of T(µ) are multiples ⊤ 2 ⊤ of s1 , so T(µ)T(µ) = O(µ) and the term (1/µ )T(µ)T(µ) = O(1/µ) → 0. For S(µ): its diagonal blocks are √ Fj∗∗ + const · ρI + µΓj∗∗ , so S(µ)/µ has diagonal blocks → Γj∗∗ , while each off-diagonal block is O( µ), giving √ S(µ)/µ off-diagonal terms of order O(1/ µ) → 0. Therefore

S(µ) 1 J̃(µ) µ→∞ = + 2 T(µ)T(µ)⊤ −−−−→ diag Γ1xx , Γ1yy , Γ2xx , Γ2yy , Γ3xx , Γ3yy ≻ 0, µ µ µ where the positive definiteness is from (60) and the stepsize condition. Since the limit is ≻ 0, there exist ϵ > 0 and N > 0 such that J̃(µ)/µ ≻ ϵI for all µ > N , hence J̃(µ) ≻ ϵµI ≻ 0, and in particular λmin (J̃(N )) > 0. Since J̃(µ) is a continuous real-symmetric matrix for µ > 0, its eigenvalues vary continuously with µ (see [33, Theorem 5.1]). By the intermediate value theorem applied to λmin (J̃(µ)), there exists µ⋆ ∈ (0, N ) such that λmin (J̃(µ⋆ )) = 0, i.e., J̃(µ⋆ ) is singular, and by (74) so is J(µ⋆ ). The preceding paragraphs complete the detailed J = 3 instability calculation. We now record why the same argument extends to arbitrary J ≥ 2. Lemma B.2 (Extension to general J). Under the hypotheses of Theorem 4.2, for any J ≥ 2: 1. det(Dg(z)) ̸= 0 for all z ∈ Ω; 2. every strict-saddle KKT point of (11) in Ω is an unstable fixed point of g. Proof. The proof uses exactly the same two ingredients as the J = 3 case. Nonsingularity. One full iteration decomposes into 2J + 1 elementary maps: one xj -update and one yj -update for each agent, followed by the dual update. Each Jacobian is block triangular with diagonal blocks given by the same implicit-function-theorem matrices as in the displayed J = 3 computation, namely Γjxx , Γjyy , and identity blocks for untouched variables and dual updates. Since every Γ block is positive definite by relative bi-smoothness and the stepsize condition, every elementary Jacobian is nonsingular; hence det(Dg) ̸= 0 for arbitrary J by the chain rule. Instability. The coupling matrix T(µ) ∈ R2Jn×(J−1)n has column j − 1 (j = 2, . . . , J) with −s1 I in the x1 -block and +s1 I in the xj -block (all other blocks zero). With the test direction d = (dx , d1y , . . . , dx , dJy , 0, . . . , 0) (all x-components equal to dx ), each column of T(µ) contributes (−s1 I)dx + (s1 I)dx = 0, hence T(µ)⊤ dp = 0. The general-J structure of the Schur-reduced matrix J̃(µ) = S(µ) + µ1 T(µ)T(µ)⊤ ∈ R2Jn×2Jn is as follows. Since T(µ) has nonzero entries only in the x-rows, the additive correction µ1 T(µ)T(µ)⊤ affects only x–x blocks, with all y-rows and y-columns equal to zero. Restricting to the x-coordinates (x1 , x2 , . . . , xJ ), we obtain   (J−1)I −I −I · · · −I  −I I 0 ··· 0   s2 1  s21   ⊤ −I 0 I ··· 0  T(µ)T(µ) =   = 1 Lstar ⊗ In , µ x-blocks,x-blocks µ  .. µ .. .. . . .. ..   .  . . −I 0 0 ··· I where Lstar ∈ RJ×J is the graph Laplacian of the star with hub node 1. The constant vector (dx , . . . , dx ) lies in ker(Lstar ⊗ In ), so µ1 T(µ)T(µ)⊤ contributes nothing to d⊤ p J̃(µ)dp when all x-components of d are equal. Moreover, √ since s21 = (1 + µ)ρ, the entries of T(µ) are O( µ), so T(µ)T(µ)⊤ = O(µ) and µ1 T(µ)T(µ)⊤ = O(1). The √ cross-agent entries of J̃(µ) are therefore at most O( µ) (from S(µ)) plus O(1) (from µ1 TT⊤ ), both of which are o(µ). Therefore, dividing by µ and sending µ → ∞ yields the block-diagonal limit
J̃(µ) µ→∞ −−−−→ diag Γ1xx , Γ1yy , . . . , ΓJxx , ΓJyy ≻ 0, µ 37

where each Γ block is positive definite by relative bi-smoothness and the stepsize condition, exactly as in the J = 3 case. For the ρ-cancellation at µ = 0, write dx,stack = (dx , . . . , dx ). The consensus penalty contribution is PJ 2 ρ d⊤ x,stack (Lstar ⊗In )dx,stack = ρ j=2 ∥dx −dx ∥2 = 0, because dx,stack is in the kernel of the star-graph Laplacian. Combined with the strict-saddle negativity of the corresponding quadratic form at µ = 0, this gives the same sign change and hence the same intermediate-value argument for arbitrary J ≥ 2. This completes the proof that strict-saddle KKT points in Ω are unstable fixed points of the fixed-point map of Algorithm 3. By Theorem 5.1, the set of initial points in Ω from which the iteration converges to such a strict-saddle KKT point has Lebesgue measure zero; since random initialization is absolutely continuous with respect to Lebesgue measure, the probability that the iterates converge to such a strict-saddle KKT point is zero.

38