Chordal decomposition in operator-splitting methods for sparse semidefinite programs

We employ chordal decomposition to reformulate a large and sparse semidefinite program (SDP), either in primal or dual standard form, into an equivalent SDP with smaller positive semidefinite (PSD) constraints. In contrast to previous approaches, the decomposed SDP is suitable for the application of first-order operator-splitting methods, enabling the development of efficient and scalable algorithms. In particular, we apply the alternating direction method of multipliers (ADMM) to solve decomposed primal- and dual-standard-form SDPs. Each iteration of such ADMM algorithms requires a projection onto an affine subspace, and a set of projections onto small PSD cones that can be computed in parallel. We also formulate the homogeneous self-dual embedding (HSDE) of a primal-dual pair of decomposed SDPs, and extend a recent ADMM-based algorithm to exploit the structure of our HSDE. The resulting HSDE algorithm has the same leading-order computational cost as those for the primal or dual problems only, with the advantage of being able to identify infeasible problems and produce an infeasibility certificate. All algorithms are implemented in the open-source MATLAB solver CDCS. Numerical experiments on a range of large-scale SDPs demonstrate the computational advantages of the proposed methods compared to common state-of-the-art solvers.


Introduction
Semidefinite programs (SDPs) are convex optimization problems over the cone of positive semidefinite (PSD) matrices.Given b ∈ R m , C ∈ S n , and matrices A 1 , . . ., Am ∈ S n , the standard primal form of an SDP is (2) In the above and throughout this work, R m is the usual m-dimensional Euclidean space, S n is the space of n × n symmetric matrices, S n + is the cone of PSD matrices, and •, • denotes the inner product in the appropriate space, i.e., x, y = x T y for x, y ∈ R m and X, Y = trace(XY ) for X, Y ∈ S n .SDPs have found applications in a wide range of fields, such as control theory, machine learning, combinatorics, and operations research [8].Semidefinite programming encompasses other common types of optimization problems, including linear, quadratic, and second-order cone programs [10].Furthermore, many nonlinear convex constraints admit SDP relaxations that work well in practice [39].
It is well-known that small and medium-sized SDPs can be solved up to any arbitrary precision in polynomial time [39] using efficient second-order interior-point methods (IPMs) [2,22].However, many problems of practical interest are too large to be addressed by the current state-of-the-art interior-point algorithms, largely due to the need to compute, store, and factorize an m × m matrix at each iteration.
A common strategy to address this shortcoming is to abandon IPMs in favour of simpler first-order methods (FOMs), at the expense of reducing the accuracy of the solution.For instance, Malick et al. introduced regularization methods to solve SDPs based on a dual augmented Lagrangian [28].Wen et al. proposed an alternating direction augmented Lagrangian method for large-scale SDPs in the dual standard form [40]. Zhao et al. presented an augmented Lagrangian dual approach combined with the conjugate gradient method to solve large-scale SDPs [45].More recently, O'Donoghue et al. developed a first-order operatorsplitting method to solve the homogeneous self-dual embedding (HSDE) of a primal-dual pair of conic programs [29].The algorithm, implemented in the C package SCS [30], has the advantage of providing certificates of primal or dual infeasibility.
A second major approach to resolve the aforementioned scalability issues is based on the observation that the large-scale SDPs encountered in applications are often structured and/or sparse [8].Exploiting sparsity in SDPs is an active and challenging area of research [3], with one main difficulty being that the optimal (primal) solution is typically dense even when the problem data are sparse.Nonetheless, if the aggregate sparsity pattern of the data is chordal (or has sparse chordal extensions), Grone's [21] and Agler's theorems [1] allow one to replace the original, large PSD constraint with a set of PSD constraints on smaller matrices, coupled by additional equality constraints.Having reduced the size of the semidefinite variables, the converted SDP can in some cases be solved more efficiently than the original problem using standard IPMs.These ideas underly the domain-space and the range-space conversion techniques in [16,24], implemented in the MATLAB package SparseCoLO [15].
The problem with such decomposition techniques, however, is that the addition of equality constraints to an SDP often offsets the benefit of working with smaller semidefinite cones.One possible solution is to exploit the properties of chordal sparsity patterns directly in the IPMs: Fukuda et al. used Grone's positive definite completion theorem [21] to develop a primal-dual path-following method [16]; Burer proposed a nonsymmetric primal-dual IPM using Cholesky factors of the dual variable Z and maximum determinant completion of the primal variable X [11]; and Andersen et al. developed fast recursive algorithms to evaluate the function values and derivatives of the barrier functions for SDPs with chordal sparsity [4].Another attractive option is to solve the sparse SDP using FOMs: Sun et al. proposed a first-order splitting algorithm for partially decomposable conic programs, including SDPs with chordal sparsity [35]; Kalbat & Lavaei applied a first-order operator-splitting method to solve a special class of SDPs with fully decomposable constraints [23]; Madani et al. developed a highly-parallelizable first-order algorithm for sparse SDPs with inequality constraints, with applications to optimal power flow problems [27].
In this work, we embrace the spirit of [23,27,29,35] and exploit sparsity in SDPs using a first-order operator-splitting method known as the alternating direction method of multipliers (ADMM).Introduced in the mid-1970s [17,19], ADMM is related to other FOMs such as dual decomposition and the method of multipliers, and it has recently found applications in many areas, including covariance selection, signal processing, resource allocation, and classification; see [9] for a review.More precisely, our contributions are: 1. Using Grone's theorem [21] and Agler's theorem [1], we formulate domain-space and range-space conversion frameworks for primal-and dual-standard-form sparse SDPs with chordal sparsity, respectively.These resemble the conversion methods developed in [16,24] for IPMs, but are more suitable for the application of FOMs.One major difference with [16,24] is that we introduce two sets of slack variables, so that the conic and the affine constraints can be separated when using operator-splitting algorithms.2. We apply ADMM to solve the domain-and range-space converted SDPs, and show that the resulting iterates of the ADMM algorithms are the same up to scaling.The iterations are cheap: the positive semidefinite (PSD) constraint is enforced via parallel projections onto small PSD cones -a computationally cheaper strategy than that in [35], while imposing the affine constraints requires solving a linear system with constant coefficient matrix, the factorization/inverse of which can be cached before iterating the algorithm.3. We formulate the HSDE of a converted primal-dual pair of sparse SDPs.In contrast to [23,27,35], this allows us to compute either primal and dual optimal points, or a certificate of infeasibility.In particular, we extend the algorithm proposed in [29] to exploit the structure of our HSDE, reducing its computational complexity.The resulting algorithm is more efficient than a direct application of the method of [29] to either the original primal-dual pair (i.e.before chordal sparsity is taken into account), or the converted problems: in the former case, the chordal decomposition reduces the cost of the conic projections; in the latter case, we speed up the projection onto the affine constraints using a series of block-eliminations.4. We present the MATLAB solver CDCS (Cone Decomposition Conic Solver), which implements our ADMM algorithms.CDCS is the first open-source first-order solver that exploits chordal decomposition and can detect infeasible problems.We test our implementation on large-scale sparse problems in SDPLIB [7], selected sparse SDPs with nonchordal sparsity pattern [4], and randomly generated SDPs with block-arrow sparsity patterns [35].The results demonstrate the efficiency of our algorithms compared to the interior-point solvers SeDuMi [34] and the first-order solver SCS [30].
The rest of the paper is organized as follows.Section 2 reviews chordal decomposition and the basic ADMM algorithm.Section 3 introduces our conversion framework for sparse SDPs based on chordal decomposition.We show how to apply the ADMM to exploit domain-space and range-space sparsity in primal and dual SDPs in Section 4. Section 5 discusses the ADMM algorithm for the HSDE of SDPs with chordal sparsity.CDCS and our numerical experiments are presented in Section 6. Section 7 concludes the paper.

A review of graph theoretic notions
We start by briefly reviewing some key graph theoretic concepts (see [6,20] for more details).A graph G(V, E) is defined by a set of vertices V = {1, 2, . . ., n} and a set of edges E ⊆ V × V.A graph G(V, E) is called complete if any two nodes are connected by an edge.A subset of vertices C ⊆ V such that (i, j) ∈ E for any distinct vertices i, j ∈ C, i.e., such that the subgraph induced by C is complete, is called a clique.The number of vertices in C is denoted by |C|.If C is not a subset of any other clique, then it is referred to as a maximal clique.
A chord is an edge joining two non-adjacent nodes in a cycle.
An undirected graph G is called chordal (or triangulated, or a rigid circuit [38]) if every cycle of length greater than or equal to four has at least one chord.Chordal graphs include several other classes of graphs, such as acyclic undirected graphs (including trees) and complete graphs.Algorithms such as the maximum cardinality search [36] can test chordality and identify the maximal cliques of a chordal graph efficiently, i.e., in linear time in terms of the number of nodes and edges.Non-chordal graphs can always be chordal extended, i.e., extended to a chordal graph, by adding additional edges to the original graph.Computing the chordal extension with the minimum number of additional edges is an NP-complete problem [42], but several heuristics exist to find a good chordal extensions efficiently [38].

Sparse matrix cones and chordal decomposition
The sparsity pattern of a symmetric matrix X ∈ S n can be represented by an undirected graph G(V, E), and vice-versa.For example, the graphs in Fig. 2 correspond to the sparsity patterns illustrated in Fig. 3.With a slight abuse of terminology, we refer to the graph G as the sparsity pattern of X.Given a clique C k of G, we define a matrix where C k (i) is the i-th vertex in C k , sorted in the natural ordering.Given X ∈ S n , the matrix E C k can be used to select a principal submatrix defined by the clique For example, the chordal graph in Fig. 1(b) has a maximal clique C 1 = {1, 2, 4}, and for Y ∈ S 3 we have Given an undirected graph G(V, E), let E * = E ∪ {(i, i), ∀ i ∈ V} be a set of edges that includes all self-loops.We define the space of sparse symmetric matrices represented by G as S n (E, 0) : and the cone of sparse PSD matrices as Moreover, we consider the cone given by the projection of the PSD cone onto the space of sparse matrices S n (E, 0) with respect to the usual Frobenius matrix norm (this is the norm induced by the usual trace inner product on the space of symmetric matrices).In is not difficult to see that X ∈ S n + (E, ?) if and only if it has a positive semidefinite completion, i.e., if there exists M 0 such that For any undirected graph G(V, E), the cones S n + (E, ?) and S n + (E, 0) are dual to each other with respect to the trace inner product in the space of sparse matrices S n (E, 0) [38].In other words, If G is chordal, then S n + (E, ?) and S n + (E, 0) can be equivalently decomposed into a set of smaller but coupled convex cones according to the following theorems: Theorem 1 (Grone's theorem [21]) Let G(V, E) be a chordal graph, and let {C 1 , C 2 , . . ., Cp} be the set of its maximal cliques.Then, X ∈ S n + (E, ?) if and only if Theorem 2 (Agler's theorem [1]) Let G(V, E) be a chordal graph, and let {C 1 , C 2 , . . ., Cp} be the set of its maximal cliques.Then, Z ∈ S n + (E, 0) if and only if there exist matrices Note that these results can be proven individually, but can also be derived from each other using the duality of the cones S n + (E, ?) and S n + (E, 0) [24].In this paper, the terminology chordal (or clique) decomposition of a sparse matrix cone will refer to the application of Theorem 1 or Theorem 2 to replace a large sparse PSD cone with a set of smaller but coupled PSD cones.Chordal decomposition of sparse matrix cones underpins much of the recent research on sparse SDPs [4,16,24,27,35,38], most of which relies on the conversion framework for IPMs proposed in [16,24].
To illustrate the concept, consider the chordal graph in Fig. 1(b).According to Grone's theorem, Similarly, Agler's theorem guarantees that (after eliminating some of the variables) Note that the PSD contraints obtained after the chordal decomposition of X (resp.Z) are coupled via the elements X 22 , X 44 , and X 24 = X 42 (resp.Z 22 , Z 44 , and Z 24 = Z 42 ).

The Alternating Direction Method of Multipliers
The computational "engine" employed in this work is the alternating direction method of multipliers (ADMM).ADMM is an operator-splitting method developed in the 1970s, and it is known to be equivalent to other operator-splitting methods such as Douglas-Rachford splitting and Spingarn's method of partial inverses; see [9] for a review.The ADMM algorithm solves the optimization problem where f and g are convex functions, x ∈ R nx , y ∈ R ny , A ∈ R nc×nx , B ∈ R nc×ny and c ∈ R nc .Given a penalty parameter ρ > 0 and a dual multiplier z ∈ R nc , the ADMM algorithm finds a saddle point of the augmented Lagrangian Lρ(x, y, z) := f (x) + g(y) by minimizing L with respect to the primal variables x and y separately, followed by a dual variable update: The superscript (n) indicates that a variable is fixed to its value at the n-th iteration.Note that since z is fixed in (4a) and (4b), one may equivalently minimize the modified Lagrangian Under very mild conditions, the ADMM converges to a solution of (3) with a rate O( 1 n ) [9, Section 3.2].ADMM is particularly suitable when (4a) and (4b) have closed-form expressions, or can be solved efficiently.Moreover, splitting the minimization over x and y often allows distributed and/or parallel implementations of steps (4a)-(4c).

Chordal decomposition of sparse SDPs
The sparsity pattern of the problem data for the primal-dual pair of standard-form SDPs ( 1)-( 2) can be described using the so-called aggregate sparsity pattern.We say that the pair of SDPs ( 1)-( 2) has aggregate sparsity pattern G(V, E) if In other words, the aggregate sparsity pattern G is the union of the individual sparsity patterns of the data matrices C, A 1 , . . ., Am.Throughout the rest of this paper, we assume that the aggregate sparsity pattern G is chordal (or that a suitable chordal extension has been found), and that it has p maximal cliques C 1 , . . ., Cp.In addition, we assume that the matrices A 1 , . .., Am are linearly independent.It is not difficult to see that the aggregate sparsity pattern defines the sparsity pattern of any feasible dual variable Z in (2), i.e. any dual feasible Z must have sparsity pattern G. Similarly, while the primal variable X in ( 1) is usually dense, the value of the cost function and the equality constraints depend only on the entries X ij with (i, j) ∈ E, and the remaining entries simply guarantee that X 0. Recalling the definition of the sparse matrix cones S n + (E, ?) and S n + (E, 0), we can therefore recast the primal-form SDP (1) as and the dual-form SDP (2) as This nonsymmetric formulation was first proposed by Fukuda et al. [16], and was later discussed in [4,24,35].Note that ( 6) and ( 7) are a primal-dual pair of linear conic problems because the cones S n + (E, ?) and S n + (E, 0) are dual to each other.

Domain-space decomposition
As we have seen in Section 2, Grone's theorem allows us to decompose the sparse matrix cone constraint X ∈ S n + (E, ?) into p standard PSD constraints on the submatrices of X defined by the cliques C 1 , . . ., Cp.In other words, These p constraints are implicitly coupled since The primal optimization problem ( 6) is then equivalent to the SDP Adopting the same terminology used in [16], we refer to (9) as the domain-space decomposition of the primal-standard-form SDP (1).
Remark 1 In the domain-space decomposition of [16,24], the primal matrix X is eliminated by replacing the constraints with Redundant constraints in ( 11) can be eliminated using the running intersection property [6] of the cliques [16], and the decomposed SDP can be solved efficiently by IPMs in certain cases [16,24].However, effectively applying FOMs to (9) after eliminating X is not straightforward.In [35] an SDP with a quadratic objective had to be solved at each iteration to impose the PSD constraints, requiring an additional iterative solver.Even when this problem is resolved, e.g. by using the algorithm of [29], the size of the KKT system enforcing the affine constraints is increased dramatically by the consensus conditions (11), sometimes so much that memory requirements are prohibitive on desktop computing platforms [16].In contrast, we show in Section 4 that if a set of slack variables X k are introduced in ( 8) and X is retained in (9), then the PSD constraint can be imposed via projections onto small PSD cones.At the same time, the affine constraints require the solution of an m × m linear system of equations, as if no consensus constraints were introduced.This makes our conversion framework more suitable to FOMs than that of [16,24].

Range-space decomposition
A range-space decomposition of the dual-standard-form SDP (2) can be formulated by applying Agler's theorem to the sparse matrix cone constraint Z ∈ S n + (E, 0) in ( 7): We then introduce slack variables V k , k = 1, . . ., p and rewrite  Similar comments as in Remark 1 hold, and the slack variables V 1 , . . ., Vp are essential to formulate a decomposition framework suitable for the application of FOMs.The rangespace decomposition of ( 2) is then given by max y,Z,V1,...,Vp b, y Remark 2 Although the domain-and range-space decompositions ( 9) and ( 12) have been derived individually, they are in fact a primal-dual pair of SDPs.The duality between the original SDPs ( 1) and ( 2) is inherited by the decomposed SDPs ( 9) and ( 12) by virtue of the duality between Grone's and Agler's theorems.This elegant picture is illustrated in Fig. 4.

ADMM for domain-and range-space decompositions of sparse SDPs
In this section, we demonstrate how ADMM can be applied to solve the domain-space decomposition (9) and the range-space decomposition (12) efficiently.Furthermore, we show that the resulting domain-and range-space algorithms are equivalent, in the sense that one is just a scaled version of the other.Throughout this section, δ K (x) will denote the indicator function of a set K, i.e.
To ease the exposition further, we consider the usual vectorized forms of ( 9) and ( 12).Specifically, we let vec : S n → R n 2 be the usual operator mapping a matrix to the stack of its column and define the vectorized data Note that the assumption that A 1 , . .., Am are linearly independent matrices means that A has full row rank.For all k = 1, . . ., p, we also introduce the vectorized variables and define "entry-selector" matrices H k := E C k ⊗ E C k for k = 1, . . ., p that project x onto the subvectors x 1 , . . ., xp, i.e. such that Note that for each k = 1, . . ., p, the rows of H k are orthonormal, and that the matrix while (12) becomes

ADMM for the domain-space decomposition
We start by moving the constraints Ax = b and x k ∈ S k in (13) to the objective using the indicator functions δ 0 (•) and δ S k (•), respectively, i.e., we write This problem is in the standard form for the application of ADMM.Given a penalty parameter ρ > 0 and a Lagrange multiplier λ k for each constraint x k = H k x, k = 1, . . ., p, we consider the (modified) augmented Lagrangian and group the variables as X := {x}, Y := {x 1 , . . ., xp}, and Z := {λ 1 , . . ., λp}.According to (4), each iteration of the ADMM requires the minimization of the Lagrangian in (16) with respect to the X -and Y-blocks separately, and follows by an update of the multipliers Z.At each step, the variables not being optimized over are fixed to their most current value.Note that splitting the primal variables x, x 1 , . . ., xp in the two blocks X and Y defined above is essential to solving the X and Y minimization subproblems (4a) and (4b); more details will be given in Remark 3 after describing the Y-minimization step in Section 4.1.2.

Minimization over X
Minimizing the augmented Lagrangian ( 16) over X is equivalent to the equality-constrained quadratic program Letting ρy be the multiplier for the equality constraint (we scale the multiplier by ρ for convenience), and defining the optimality conditions for ( 17) can be written as the KKT system Recalling that the product H T k H k is a diagonal matrix for all k = 1, . . ., p we conclude that so is D, and since A has full row rank by assumption ( 19) can be solved efficiently, for instance by block elimination.In particular, eliminating x shows that the only matrix to be inverted/factorized is Incidentally, we note that the first-order algorithms of [29,40] require the factorization of a similar matrix with the same dimension.Since this matrix is the same at every iteration, its Cholesky factorization (or any other factorization of choice) can be computed and cached before starting the ADMM iterations.For some families of SDPs, such as the SDP relaxation of MaxCut problems and sum-of-squares (SOS) feasibility problems [46], the matrix AD −1 A T is diagonal, so solving (19) is inexpensive even when the SDPs are very large.If factorizing AD −1 A T is too expensive, the linear system (19) can alternatively be solved by an iterative method, such as the conjugate gradient method [33].

Minimization over Y
Minimizing the augmented Lagrangian (16) over Y is equivalent to solving p independent conic problems of the form min In terms of the original matrix variables X 1 , . . ., Xp, each of these p sub-problems amounts to a projection on a PSD cone.More precisely, if P S k denotes the projection onto the PSD cone S k and mat(•) = vec −1 (•), we have Since the projection P S k can be computed with an eigenvalue decomposition, and since the size of each cone S |C k | + is small for typical sparse SDPs (such as SDP relaxations of MaxCut problems), the variables x 1 , . . ., xp can be updated efficiently.Moreover, the computation can be carried out in parallel.In contrast, the algorithms for generic SDPs developed in [28,29,40] require projections onto the original large PSD cone S n + .
Remark 3 As anticipated in Remark 1, retaining the global variable x in the domain-space decomposed SDP to enforce the consensus constraints between the entries of the subvectors x 1 , . . . ,xp (i.e., x k = H k x) is fundamental.In fact, it allowed us to separate the conic constraints from the affine constraints in (13) when applying the splitting strategy of ADMM, making the minimization over Y easy to compute and parallelizable.In contrast, when x is eliminated as in the conversion method of [16,24], the conic constraints and the affine constraints cannot be easily decoupled when applying the first-order splitting method: in [35] a quadratic SDP had to be solved at each iteration, impeding the scalability of the algorithm.

Updating the multipliers Z
The final step in the n-th ADMM iteration is to update the multipliers λ 1 . . ., λp with the usual gradient ascent rule: for each k = 1, . . ., p, λ This computation is cheap and easily parallelized.

Summary & Stopping conditions
The ADMM algorithm is stopped after the n-th iteration if the relative primal/dual error measures are smaller than a specified tolerance, tol .The reader is referred to [9] for a detailed discussion of stopping conditions for ADMM algorithms.In conclusion, a primal-form SDP with domain-space decomposition (13) can be solved using the steps summarized in Algorithm 1.

ADMM for the range-space decomposition
An ADMM algorithm similar to Algorithm 1 can be developed for the range-space decomposition ( 14) of a dual-standard-form sparse SDP.As in Section 4.1, we start by moving Algorithm 1 ADMM for the domain-space decomposition of sparse SDPs 1: Set ρ > 0, tol > 0, a maximum number of iterations nmax, and initial guesses x (0) , x end if 15: end for all but the consensus equality constraints z k = v k , k = 1, . . ., p, to the objective using indicator functions.This leads to Given a penalty parameter ρ > 0 and a Lagrange multiplier λ k for each of the constraints z k = v k , k = 1, . . ., p, we consider the (modified) augmented Lagrangian L(y, v 1 , . . ., vp, z 1 , . . ., zp, λ 1 , . . ., λp) := − b, y and consider three groups of variables, X := {y, v 1 , . . ., vp}, Y := {z 1 , . . ., zp}, and Z := {λ 1 , . . ., λp}.Similar to Section 4.1, each iteration of the ADMM algorithm for (14) consists of minimizations over X and Y, and an update of the multipliers Z.Each of these steps admits an inexpensive closed-form solution, as we demonstrate next.

Minimization over X
Minimizing (26) over block X is equivalent to solving the equality-constrained quadratic program min y,v1,...,vp Let ρx be the multiplier for the equality constraint.After some algebra, the optimality conditions for ( 27) can be written as the KKT system plus a set of p uncoupled equations for the variables v k , The KKT system (28) is the same as (19) after rescaling x → −x, y → −y, c → ρ −1 c and b → ρb.Consequently, the numerical cost of these operations is the same as in Section 4.1.1,plus the cost of (29), which is cheap and can be parallelized.Moreover, as in Section 4.1.1,the factors of the coefficient matrix required to solve the KKT system (28) can be pre-computed and cached, before iterating the ADMM algorithm.

Minimization over Y
As in Section 4.1.2,the variables z 1 , . . ., zp are updated with p independent projections, where P S k denotes projection on the PSD cone S |C k | + .Again, these projections can be computed efficiently and in parallel.
Remark 4 As anticipated in Section 3.2, introducing the set of slack variables v k and the consensus constraints z k = v k , k = 1, . . ., p is essential to obtain an efficient algorithm for range-space decomposed SDPs.The reason is that the splitting strategy of the ADMM decouples the conic and affine constraints, and the conic variables can be updated using the simple conic projection (30).

Updating the multipliers Z
The multipliers λ k , k = 1, . . ., p, are updated (possibly in parallel) with the cheap gradient ascent rule

Summary & Stopping conditions
Similarly to Section 4.1.4,we stop our ADMM algorithm after the n-th iteration if the relative primal/dual error measures Algorithm 2 ADMM for dual form SDPs with range-space decomposition 1: Set ρ > 0, tol > 0, a maximum number of iterations nmax and initial guesses y (0) , z end if 16: end for are smaller than a specified tolerance, tol .The ADMM algorithm to solve the range-space decomposition ( 14) of a dual-form sparse SDP is summarized in Algorithm 2.

Equivalence between the primal and dual ADMM algorithms
Since the computational cost of ( 29) is the same as ( 23), all ADMM iterations for the dualform SDP with range-space decomposition ( 14) have the same cost as those for the primalform SDP with domain-space decomposition (13), plus the cost of (31).However, if one minimizes the dual augmented Lagrangian (26) over z 1 , . . ., zp before minimizing it over y, v 1 , . . ., vp, then (29) can be used to simplify the multiplier update equations to Given that the products H 1 x, . . ., Hpx have already been computed to update v 1 , . . ., vp in (29), updating the multipliers λ 1 , . . ., λp requires only a scaling operation.Recalling that ( 19) and ( 28) are scaled versions of the same KKT system, after swapping the order of the minimization, the ADMM algorithms for the primal and dual standard form SDPs can be considered as scaled versions of each other; see Fig. 4 for an illustration.In fact, the equivalence between ADMM algorithms for the original (i.e., before chordal decomposition) primal and dual SDPs was already noted in [41].
Remark 5 Although the iterates of Algorithm 1 and Algorithm 2 are the same up to scaling, the convergence performance of these two algorithms differ in practice because first-order methods are sensitive to the scaling of the problem data and of the iterates.

Homogeneous self-dual embedding of domain-and range-space decomposed SDPs
Algorithms 1 and 2, as well as other first-order algorithms that exploit chordal sparsity [23,27,35], can solve feasible problems, but cannot detect infeasibility in their current formulation.Although some recent ADMM methods reslove this issue [5,25], an elegant way to deal with an infeasible primal-dual pair of SDPs-which we pursue here-is to solve their homogeneous self-dual embedding (HSDE) [44].
The essence of the HSDE method is to search for a non-zero point in the intersection of a convex cone and a linear space; this is non-empty because it always contains the origin, meaning that the problem is always feasible.Given such a non-zero point, one can either recover optimal primal and dual solutions of the original pair of optimization problems, or construct a certificate of primal or dual infeasibility.HSDEs have been widely used to develop IPMs for SDPs [34,43], and more recently O'Donoghue et al. have proposed an operator-splitting method to solve the HSDE of general conic programs [29].
In this section, we formulate the HSDE of the domain-and range-space decomposed SDPs ( 13) and ( 14), which is a primal-dual pair of SDPs.We also apply ADMM to solve this HSDE; in particular, we extend the algorithm of [29] to exploit chordal sparsity without increasing its computational cost (at least to leading order) compared to Algorithms 1 and 2.

Homogeneous self-dual embedding
To simplify the formulation of the HSDE of the decomposed (vectorized) SDPs ( 13) and ( 14), we let S := S 1 × • • • × Sp be the direct product of all semidefinite cones and define When strong duality holds, the tuple (x * , s * , y * , v * , z * ) is optimal if and only if all of the following conditions hold: 1. (x * , s * ) is primal feasible, i.e., Ax * = b, s * = Hx * , and s * ∈ S. For reasons that will become apparent below, we introduce slack variables r * = 0 and w * = 0 of appropriate dimensions and rewrite these conditions as 2. (y * , v * , z * ) is dual feasible, i.e., A T y * + H T v * = c, z * = v * , and z * ∈ S. Again, it is convenient to introduce a slack variable h * = 0 of appropriate size and write 3. The duality gap is zero, i.e.
The idea behind the HSDE [44] is to introduce two non-negative and complementary variables τ and κ and embed the optimality conditions ( 34), ( 35) and ( 36) into the linear system v = Qu with u, v and Q defined as Any nonzero solution of this embedding can be used to recover an optimal solution for ( 9) and ( 12), or provide a certificate for primal or dual infeasibility, depending on the values of τ and κ; details are omitted for brevity, and the interested reader is referred to [29].
The decomposed primal-dual pair of (vectorized) SDPs ( 13)-( 14) can therefore be recast as the self-dual conic feasibility problem find (u, v) where, writing

A simplified ADMM algorithm
The feasibility problem ( 38) is in a form suitable for the application of ADMM, and moreover steps (4a)-(4c) can be greatly simplified by virtue of its self-dual character [29].Specifically, the n-th iteration of the simplified ADMM algorithm for (38) proposed in [29] consists of the following three steps, where P K denotes projection onto the cone K: Note that (39b) is inexpensive, since K is the cartesian product of simple cones (zero, free and non-negative cones) and small PSD cones, and can be efficiently carried out in parallel.The third step is also computationally inexpensive and parallelizable.On the contrary, even when the preferred factorization of I + Q (or its inverse) is cached before starting the iterations a direct implementation of (39a) may require substantial computational effort because is a very large matrix (e.g., n 2 + 2n d + m + 1 = 2360900 for the instance of rs365 in Section 6.3, which would take over 10 4 GB to store Q as a dense double-precision matrix).Yet, as we can see in (37), Q is highly structured and sparse, and these properties can be exploited to speed up step (39a) using a series of block-eliminations and the matrix inversion lemma [10,Section C.4.3].

Solving the "outer" linear system
The affine projection step (39a) requires the solution of a linear system (which we refer to as the "outer" system for reasons that will become clear below) of the form where and we have split Note that û2 and ω 2 are scalars.After one step of block elimination in (40) we obtain Moreover, applying the matrix inversion lemma [10,Section C.4.3] to (43) shows that Note that the vector M −1 ζ and the scalar 1 + ζ T (M −1 ζ) depend only on the problem data, and can be computed before starting the ADMM iterations (since M is quasi-definite it can be inverted, and any symmetric matrix obtained as a permutation of M admits an LDL T factorization).Instead, recalling from ( 42) that ω 1 − ω 2 ζ changes at each iteration because it depends on the iterates u (n) and v (n) , the vector M −1 (ω 1 − ω 2 ζ) must be computed at each iteration.Consequently, computing û1 and û2 requires the solution of an "inner" linear system for the vector M −1 (ω 1 − ω 2 ζ), followed by inexpensive vector inner products and scalar-vector operations in ( 45) and (44).

Solving the "inner" linear system
Recalling the definition of M from (41), the "inner" linear system to calculate û1 in (45) has the form where σ 1 and σ 2 are the unknowns and represent suitable partitions of the vector (45) (which is to be calculated), and where we have split Applying block elimination to remove σ 1 from the second equation in (46), we obtain Recalling the definition of Â and recognizing that is a diagonal matrix, as already noted in Section 4.1.1,we also have Block elimination can therefore be used once again to solve (48), and simple algebraic manipulations show that the only matrix to be factorized (or inverted) is Note that this matrix depends only on the problem data and the chordal decomposition, so it can be factorized/inverted before starting the ADMM iterations.In addition, it is of the "diagonal plus low rank" form because A ∈ R m×n 2 with m < n 2 (in fact, often m n 2 ).This means that the matrix inversion lemma can be used to reduce the size of the matrix to factorize/invert even further: letting P = I + 1 2 D be the diagonal part of (50), we have In summary, after a series of block eliminations and applications of the matrix inversion lemma, step (39a) of the ADMM algorithm for (38) only requires the solution of an m × m linear system of equations with coefficient matrix Algorithm 3 ADMM for the HSDE of sparse SDPs with chordal decomposition 1: Set tol > 0, a maximum number of iterations nmax and initial guesses for the variables û(0) , u (0) , v (0) 2: Data preprocessing: chordal extension, chordal decomposition and factorization of the matrix in (51).3: for n = 1, . . ., nmax do 4: Compute û(n+1) using the sequence of block eliminations ( 40)-(51).5: Compute u (n+1) using (39b).6: Compute v (n+1) using (39c).Certificates of primal or dual infeasibility (with tolerance tol ) are then given, respectively, by the points These stopping criteria are identical to those used by many other conic solvers, e.g., SCS [30].The complete ADMM algorithm to solve the HSDE of the primal-dual pair of domain-and range-space decomposed SDPs is summarized in Algorithm 3.

Summary of computational gains
Algorithm 3 is clearly more efficient than a direct application of the ADMM algorithm of [29] to the decomposed primal-dual pair of (vectorized) SDPs ( 13)- (14).In fact, the cost of the conic projection (39b) is the same for both algorithms, but the sequence of block eliminations and applications of the matrix inversion lemma we have described greatly reduces the cost of the affine projection step: we only need to invert/factorize an m × m matrix, instead of the (n 2 + 2n d + m + 1) × (n 2 + 2n d + m + 1) matrix Q (as we noted before, n 2 + 2n d + m + 1 is usually very large).
Furthermore, it can be checked that when we exploit the special structure of the matrix I + Q, the overall computational cost of (39a) coincides (to leading order) with the cost of the affine projection step when the algorithm of [29] is applied to the original primal-dual pair (1)-(2), i.e., before chordal decomposition is applied.This means that our algorithm should also outperform the algorithm of [29] when it is applied to the original primal-dual pair of SDPs ( 1)-( 2): the cost of the affine projection is the same, but the conic projection in Algorithm 3 is cheaper because we work with smaller PSD cones.
Finally, note from ( 20) and (51) that the matrices to be inverted/factorized in Algorithms 1-3 have the same dimensions.Moreover, since matrix D is diagonal, AD −1 A T and I + A I + 1 2 D −1 A T have the same sparsity pattern, and hence the number of flops required to factorize them is the same.In addition, the computational cost of the conic projection step in Algorithms 1-3 is dominated by the projection on the PSD cones S |C k | , k = 1, . . ., p, which are the same in all three algorithms (they are determined by the chordal decomposition of the aggregate sparsity pattern of the original problem).Consequently, each iteration of our ADMM algorithm for the HSDE formulation has the same leading order cost as applying the ADMM to the primal and dual problem alone, with the major advantage of being able to detect infeasibility.

Implementation and numerical experiments
We implemented Algorithms 1-3 in an open-source MATLAB solver which we call CDCS (Cone Decomposition Conic Solver).We refer to our implementation of Algorithms 1-3 as CDCS-primal, CDCS-dual and CDCS-hsde, respectively.This section briefly describes CDCS, and presents numerical results of some sparse SDPs from SDPLIB [7], some large and sparse SDPs with nonchordal sparsity patterns from [4], and randomly generated SDPs with block-arrow sparsity pattern.Such problems have also been used as benchmarks in [4,35].
In order to highlight the advantages of chordal decomposition, first-order algorithms, and their combination, the three algorithms in CDCS are compared to the interior-point solver SeDuMi [34], and to the single-threaded direct implementation of the first-order algorithm of [29] provided by the conic solver SCS [30].All solvers were called with termination tolerance tol = 10 −3 , number of iterations limited to 2 000, and their default remaining parameters.The purpose of comparing CDCS to a low-accuracy IPM is to demonstrate the advantages of combining FOMs with chordal decomposition, while a comparison to the high-performance first-order conic solver SCS highlights the advantages of chordal decomposition alone.Accurate solutions ( tol = 10 −8 ) were also computed using SeDuMi; these can be considered "exact", and used to assess how far the solution returned by CDCS is from optimality.All experiments were carried out on a PC with a 2.8 GHz Intel Core i7 CPU and 8GB of RAM.

CDCS
To the best of our knowledge, CDCS is the first open-source first-order conic solver that exploits chordal decomposition for the PSD cones and is able to handle infeasible problems.Cartesian products of the following cones are supported: the cone of free variables R n , the non-negative orthant R n + , second-order cones, and PSD cones.The current implementation is written in MATLAB and can be downloaded from https://github.com/oxfordcontrol/cdcs .Note that although many steps of Algorithms 1-3 can be carried out in parallel, our implementation is sequential.Interfaces with the optimization toolboxes YALMIP [26] and SOSTOOLS [31] are also available.

Implementation details
CDCS applies chordal decomposition to all PSD cones.Following [38], the sparsity pattern of each PSD cone is chordal extended using the MATLAB function chol to compute a symbolic Cholesky factorization of the approximate minimum-degree permutation of the cone's adjacency matrix, returned by the MATLAB function symamd.The maximal cliques of the chordal extension are then computed using a .mexfunction from SparseCoLO [15].
As far as the steps of our ADMM algorithms are concerned, projections onto the PSD cone are performed using the MATLAB routine eig, while projections onto other supported cones only use vector operations.The Cholesky factors of the m × m linear system coefficient matrix (permuted using symamd) are cached before starting the ADMM iterations.The permuted linear system is solved at each iteration using the routines cs lsolve and cs ltsolve from the CSparse library [12].

Adaptive penalty strategy
While the ADMM algorithms proposed in the previous sections converge independently of the choice of penalty parameter ρ, in practice its value strongly influences the number of iterations required for convergence.Unfortunately, analytic results for the optimal choice of ρ are not available except for very special problems [18,32].Consequently, in order to improve the convergence rate and make performance less dependent on the choice of ρ, CDCS employs the dynamic adaptive rule , Here, k p and k d are the primal and dual residuals at the k-th iteration, while µ incr , µ decr and ν are parameters no smaller than 1.Note that since ρ does not enter any of the matrices being factorized/inverted, updating its value is computationally cheap.
The idea of the rule above is to adapt ρ to balance the convergence of the primal and dual residuals to zero; more details can be found in [9,Section 3.4.1].Typical choices for the parameters (the default in CDCS) are µ incr = µ decr = 2 and ν = 10 [9].

Scaling the problem data
The relative scaling of the problem data also affects the convergence rate of ADMM algorithms.CDCS scales the problem data after the chordal decomposition step using a strategy similar to [29].In particular, the decomposed SDPs ( 13) and ( 14) can be rewritten as: CDCS solves the scaled problems  Section 5]), an optimal point for (53) can be recovered from the solution of (54) according to

Sparse SDPs from SDPLIB
Our first experiment is based on benchmark problems from SDPLIB [7]: two Lovász ϑ number SDPs (theta1 and theta2); two infeasible SDPs (infd1 and infd2); four largescale sparse SDPs, two MaxCut problems (maxG11 and maxG32) and two SDP relaxations of box-constrained quadratic programs (qpG11 and qpG51).Table 1 reports the dimensions of these problems, as well as chordal decomposition details.Problems theta1 and theta2 are dense, so have only one maximal clique; all other problems are sparse and have many maximal cliques of size much smaller than the original cone.
The numerical results are summarized in Tables 2-6.Table 2 shows that the small dense SDPs theta1 and theta2, were solved in approximately the same CPU time by all solvers.Note that since these problems only have one maximal clique, SCS and CDCS-hsde use similar algorithms, and performance differences are mainly due to the implementation (most notably, SCS is written in C).Table 3 confirms that CDCS-hsde successfully detects infeasible problems, while CDCS-primal and CDCS-dual do not have this ability.
The results for the four large-scale sparse SDPs are shown in Tables 4 and 5.All algorithms in CDCS were faster than either SeDuMi or SCS, in particular-as one would  expect-for problems with smaller maximum clique size.Notably, CDCS-dual and CDCShsde solved maxG11, maxG32, and qpG11 in less than 100 s, a speedup of approximately 11×, 48×, and 64× over SCS.Table 6 reports the average CPU time per iteration for CDCS and SCS.This metric gives a fairer comparison of the performance of the algorithms because, in contrast to the total CPU time, it does not depend on the exact stopping conditions.All algorithms in CDCS are faster than SCS for the large-scale sparse SDPs (maxG11, maxG32, qpG11 and qpG51), and in particular CDCS-hsde improves on SCS by approximately 1.5×, 7.9×, 7.9×, and 2.7× for each problem, respectively.This is to be expected since the conic projection step in CDCS is more efficient due to smaller semidefinite cones, but the results are remarkable considering that CDCS is written in MATLAB, while SCS is implemented in C. In fact, the performance of CDCS could be improved even further with a parallel implementation of the projections onto small PSD cones.Finally, note that although FOMs are only meant to provide moderately accurate solutions, the objective value returned by CDCS-hsde was always within 0.2% of the highaccuracy optimal value computed using SeDuMi.This is an acceptable difference in many practical applications.Fig. 5 Aggregate sparsity patterns of the nonchordal SDPs in [4]; see Table 7 for the matrix dimensions.

Nonchordal SDPs
In our second experiment, we solved six large-scale SDPs with nonchordal sparsity patterns form [4]: rs35, rs200, rs228, rs365, rs1555, and rs1907.The aggregate sparsity patterns of these problems, illustrated in Fig. 5, come from the University of Florida Sparse Matrix Collection [13].Table 7 demonstrates that all six sparsity patterns admit chordal extensions with maximum cliques that are much smaller than the original cone.
The numerical results are presented in Table 8 and Table 9.For all problems, the algorithms in CDCS (primal, dual and hsde) are all much faster than either SCS or SeDuMi.For the largest instance, rs1555, CDCS-hsde returned successfully winthin 20 minutes, 100 times faster than SCS, which stopped after 38 hours having reached the maximum number of iterations without meeting the convergence tolerance.In fact, SCS never terminates succesfully, while the objective value returned by CDCS is always within 2% of the high-accuracy solutions returned by SeDuMi (when this could be computed).The average CPU time per iteration in CDCS-hsde is also 20×, 21×, 26×, and 75× faster than SCS, respectively, for problems rs200, rs365, rs1907, and rs1555.In addition, the results show that the average CPU time per iteration for CDCS (primal, dual, and hsde) is independent of the original problem size and, perhaps not unexpectedly, it seems to depend mainly on the dimension of the largest clique.In fact, in all of our algorithms the complexity of the conic projection-which dictates the overall complexity when m is fixed to a moderate value like in the examples presented here-is determined by the size of the largest maximal clique, not the size of the cone in the original problem.

Random SDPs with block-arrow patterns
In our last experiment, we consider randomly generated SDPs with the block-arrow aggregate sparsity pattern illiustrated in Figure 6.Such a sparsity pattern, used as a benchmark case in [4,35], is chordal.Its parameters are: the number of blocks, l; the block size,  The CPU times for the different solvers considered in this work are shown in Figure 7.In all three test scenarios, CDCS is more than 10 times faster than SeDuMi even when a low tolerance is used, and on average it is also faster than SCS.In addition, we observed that the optimal value returned by all algorithms in CDCS (primal, dual and hsde) was always within 0.1% of the high-accuracy value returned by SeDuMi, which can be considered to be negligible in practice.

Conclusion
In this paper, we have presented a conversion framework for large-scale SDPs characterized by chordal sparsity.This framework is analogous to the conversion techniques for IPMs of [16,24], but is more suitable for the application of FOMs.We have then developed efficient ADMM algorithms for sparse SDPs in either primal or dual standard form, and for their homogeneous self-dual embedding.In all cases, a single iteration of our ADMM algorithms only requires parallel projections onto small PSD cones and a projection onto an affine subspace, both of which can be carried out efficiently.In particular, when the number of constraints m is moderate the complexity of each iteration is determined by the size of the largest maximal clique, not the size of the original problem.This enables us to solve large, sparse conic problems that are beyond the reach of standard interior-point and/or other first-order methods.
All our algorithms have been made available in the open-source MATLAB solver CDCS.Numerical simulations on benchmark problems, including selected sparse problems from SDPLIB, large and sparse SDPs with a nonchordal sparsity pattern, and SDPs with a blockarrow sparsity pattern, demonstrate that our methods can significantly reduce the total CPU time requirement compared to the state-of-the-art interior-point solver SeDuMi [34] and the efficient first-order solver SCS [30].We remark that the current implementation of our algorithms is sequential, but many steps can be carried out in parallel, so further computational gains may be achieved by taking full advantage of distributed computing architectures.Besides, it would be interesting to integrate some acceleration techniques (e.g., [14,37]) that are promising to improve the convergence performance of ADMM in practice.
Finally, we note that the conversion framework we have proposed relies on chordal sparsity, but there exist large SDPs which do not have this property.An example with applications in many areas are SDPs from sum-of-squares relaxations of polynomial optimization problems.Future work should therefore explore whether, and how, first order methods can be used to take advantage other types of sparsity and structure.
while the standard dual form ismax y, Z b, y s.t.Z + m i=1 A i y i = C, Z ∈ S n + .

Fig. 4
Fig. 4 Duality between the original primal and dual SDPs, and the decomposed primal and dual SDPs.

Fig. 6
Fig. 6 Block-arrow sparsity pattern (dots indicate repeating diagonal blocks).The parameters are: the number of blocks, l; block size, d; the width of the arrow head, h.

Fig. 7
Fig. 7 CPU time for SDPs with block-arrow patterns.Left to right: varying number of constraints; varying number of blocks; varying block size. p.

Table 1
Details of the SDPLIB problems considered in this work.

Table 2
[29]lts for two small SDPs, theta1 and theta2, in SDPLIB.vectorsbandĉ by positive scalars ρ and σ, and the primal and dual equality constraints by positive definite, diagonal matrices D and E. Note that such a rescaling does not change the sparsity pattern of the problem.As already observed in[29], a good choice for E, D, σ and ρ is such that the rows of Ā and b have Euclidean norm close to one, and the columns of Ā and c have similar norms.If D and D −1 are chosen to preserve membership to the cone R n 2 × K and its dual, respectively (how this can be done is explained in[29,

Table 3
Results for two infeasible SDPs in SDPLIB.An objective value of +Inf denotes infeasiblity.Results for the primal-only and dual-only algorithms in CDCS are not reported since they cannot detect infeasibility.

Table 4
Results for two large-scale sparse SDPs from MaxCut problems in SDPLIB, maxG11 and maxG32.
†: maximum number of iterations reached.

Table 5
Results for two large-scale sparse SDPs from box-constrained QPs in SDPLIB, qpG11 and qpG51.
***: the problem could not be solved due to memory limitations.†:maximum number of iterations reached.

Table 6
Average CPU time per iteration (in seconds) for the SDPs from SDPLIB tested in this work.

Table 7
[4]mary of chordal decomposition for the chordal extensions of the nonchordal SDPs form[4].
***: the problem could not be solved due to memory limitations.†:maximum number of iterations reached.