Improving ADMMs for Solving Doubly Nonnegative Programs through Dual Factorization

Augmented Lagrangian methods are among the most popular first-order approaches to handle large scale semidefinite programs. In particular, alternating direction methods of multipliers (ADMMs), which are a variant of augmented Lagrangian methods, gained attention during the past decade. In this paper, we focus on solving doubly nonnegative programs (DNN), which are semidefinite programs where the elements of the matrix variable are constrained to be nonnegative. Starting from two algorithms already proposed in the literature on conic programming, we introduce two new ADMMs by employing a factorization of the dual variable. It is well known that first order methods are not suitable to compute high precision optimal solutions, however an optimal solution of moderate precision often suffices to get high quality lower bounds on the primal optimal objective function value. We present methods to obtain such bounds by either perturbing the dual objective function value or by constructing a dual feasible solution from a dual approximate optimal solution. Both procedures can be used as a post-processing phase in our ADMMs. Numerical results for DNNs that are relaxations of the stable set problem are presented. They show the impact of using the factorization of the dual variable in order to improve the progress towards the optimal solution within an iteration of the ADMM. This decreases the number of iterations as well as the CPU time to solve the DNN to a given precision. The experiments also demonstrate that within a computationally cheap post-processing, we can compute bounds that are close to the optimal value even if the DNN was solved to moderate precision only. This makes ADMMs applicable also within a branch-and-bound algorithm.


Introduction
In a semidefinite program (SDP) one wants to find a positive semidefinite (and hence symmetric) matrix such that linear -in the entries of the matrix -constraints are fulfilled and a linear objective function is minimized.If the matrix is also required to be entrywise nonnegative, the problem is called doubly nonnegative program (DNN).Since interior point methods fail (in terms of time and memory required) when the scale of the SDP is big, augmented Lagrangian approaches became more and more popular to solve this class of programs.Wen, Goldfarb and Yin [15] as well as Malick, Povh, Rendl and Wiegele [10] and De Santis, Rendl and Wiegele [3] considered alternating direction methods of multipliers (ADMMs) to solve SDPs.One can directly apply these ADMMs to solve DNNs, too, by introducing nonnegative slack variables for the nonnegativity constraints in order to obtain equality constraints only.However, this increases the size of the problem significantly.
In this paper, we first present two ADMMs already proposed in the literature (namely ConicADMM3c by Sun, Toh and Yang [14] and ADAL+ [15]) to specifically solve DNNs.Then we introduce two new methods: DADMM3c, which is convergent and employs a factorization of the dual matrix to avoid spectral decompositions, and DADAL+ taking advantage of the practical benefits of DADAL [3].Note that there are examples for which a 3-block ADMM (like DADAL+) diverges.However, the question of convergence of 3-block ADMMs for SDP relaxations arising from combinatorial optimization problems is still open.
In case the DNN is used as relaxation of some combinatorial optimization problem, one is interested in dual bounds, i.e. bounds that are the dual objective function value of a dual feasible solution.In case of a minimization problem this is a lower bound, in case of a maximization problem an upper bound.Having bounds is in particular important if one intends to use the relaxation within a branch-and-bound algorithm.This, however, means that one needs to solve the DNN to high precision such that the dual solution is feasible and hence the dual objective function value is a reliable bound.Typically, first order methods can compute solutions of moderate precision in reasonable time, whereas progressing to higher precision can become expensive.To overcome this drawback, we present two methods to compute a dual bound from a solution obtained by the ADMMs within a post-processing phase.
In the following section we state our notations and introduce the formulation of standard primal-dual SDPs and DNNs.In Section 2 we go through the two existing ADMMs for DNNs we mentioned before, and in Section 3 we introduce the tool of dual matrix factorization used in the new ADMMs DADAL+ and DADMM3c presented later in the same section.In Section 4 we present two methods for obtaining dual bounds from a solution of a DNN that satisfies the optimality criteria to moderate precision only.Section 5 shows numerical results for instances of DNN relaxations of the stable set problem.We evaluate the impact of the dual factorization within the methods as well as the two post-processing schemes for obtaining dual bounds.Section 6 concludes the paper.

Problem Formulation and Notations
Let S n be the set of n-by-n symmetric matrices, S + n ⊂ S n be the set of positive semidefinite matrices and S − n ⊂ S n be the set of negative semidefinite matrices.Denoting by X, Y = trace(XY ) the standard inner product in S n , we write the standard primal-dual pair of SDPs where C ∈ S n , b ∈ R m , A : S n → R m is the linear operator (AX) i = A i , X with A i ∈ S n , i = 1, . . ., m and A ⊤ : R m → S n is its adjoint operator, so A ⊤ y = i y i A i for y ∈ R m .When in the primal SDP (1) the elements of X are constrained to be nonnegative, then the SDP is called a doubly nonnegative program (DNN).To be more precise the primal DNN is given as min C, X Introducing S as the dual variable related to the nonnegativity constraint X ≥ 0, we write the dual of the DNN (3) as We assume that both the primal DNN (3) and the dual DNN (4) have strictly feasible points (i.e.Slater's condition is satisfied), so strong duality holds.Under this assumption, (y, S, Z, X) is optimal for (3) and (4) if and only if hold.We further assume that the constraints formed through the operator A are linearly independent.
Let v ∈ R n and M ∈ R m×n .In the following, M (i, :) is defined as the i-th row of M and M (:, j) as the j-th column of M .Further we denote by Diag(v) the diagonal matrix having v on the main diagonal.The vector e i is defined as the i-th vector of the standard basis in R n .Whenever a norm is used, we consider the Frobenius norm in case of matrices and the Euclidean norm in case of vectors.Let S ∈ S n .We denote the projection of S onto the positive semidefinite and negative semidefinite cone by (S) + and (S) − , respectively.The projection of S onto the nonnegative orthant is denoted by (S) ≥0 .Moreover we denote by λ(S) the vector of the eigenvalues of S and by λ min (S) and λ max (S) the smallest and largest eigenvalue of S, respectively.

ADMMs for Doubly Nonnegative Programs
In this section, we present two different ADMMs for solving DNNs.Let X ∈ S n be the Lagrange multiplier for the dual equation A ⊤ y + Z + S − C = 0 and σ > 0 be fixed.Then the augmented Lagrangian of the dual DNN ( 4) is defined as In the classical augmented Lagrangian method applied to the dual DNN (4) the problem max L σ (y, S, Z; X) where X is fixed and σ > 0 is a penalty parameter is addressed at every iteration.Once Problem ( 6) is (approximately) solved, the multiplier X is updated by the first order rule and the process is iterated until convergence, i.e., until the optimality conditions ( 5) are satisfied within a certain tolerance (see [1,Chapter 2] for further details).
If the augmented Lagrangian L σ (y, S, Z; X) is maximized with respect to y, S and Z not simultaneously but one after the other, this yields the well known alternating direction method of multipliers (ADMM).The number of blocks of an ADMM is the number of blocks of variables for which Problem ( 6) is maximized separately, so we consider a 3-block ADMM.Such an ADMM has been specialized and used by Wen, Goldfarb and Yin [15] to address DNNs and in the following we refer to this method as ADAL+ as it is an alternating direction augmented Lagrangian method for DNNs.Details will be given in Section 2.1.Even though in all our numerical tests this algorithm reaches the desired precision of our stopping criteria, it has been recently shown in [2] that an ADMM with more than two blocks may diverge.
In order to overcome this theoretical issue, Sun, Toh and Yang [14] proposed to update the third block twice per iteration, or, in other words, to maximize L σ (y, S, Z; X) with respect to the variable y two times in one iteration.Their algorithm, named ConicADMM3c and detailed in Section 2.2, is the first theoretically convergent 3-block ADMM proposed in the context of conic programming.

ADAL+
In the following, we refer to the ADMM presented by Wen, Goldfarb and Yin in [15] and applied to the dual DNN (4) as ADAL+.As already mentioned, ADAL+ iterates the maximization of the augmented Lagrangian with respect to each block of dual variables.To be more precise the new point (y k+1 , S k+1 , Z k+1 , X k+1 ) is computed by the following steps: The update of y in ( 8) is derived from the first-order optimality condition of the problem on the right-hand side of (8), so y k+1 is the unique solution of As shown in [15], the update of S according to ( 9) is equivalent to min where Hence, S k+1 is obtained as the projection of U k+1 onto the nonnegative orthant, namely Then, the update of Z in ( 10) is conducted by considering the equivalent problem min with ), or, in other words, by projecting W k+1 ∈ S n onto the (closed convex) cone S − n and taking its additive inverse (see Algorithm 1).Such a projection is computed via the spectral decomposition of the matrix W k+1 .
Finally, it is easy to see that the update of X in (11) can be performed considering the projection of W k+1 ∈ S n onto S + n multiplied by σ k , namely We report in Algorithm 1 the scheme of ADAL+.
The other optimality conditions (namely X ∈ S + n , Z ∈ S + n , S ∈ S n , S ≥ 0, ZX = 0) are satisfied up to machine accuracy throughout the algorithm thanks to the projections employed in ADAL+.

ConicADMM3c
A major drawback of ADAL+ is that it is not necessarily convergent.By considering two updates of the variable y within one iteration, Sun, Toh and Yang are able to prove that the algorithm ConicADMM3c proposed in [14] and detailed in Algorithm 2 is a 3-block convergent ADMM: Under certain assumptions on the penalty parameter σ, they show that the sequence {(y k , S k , Z k ; X k )} produced by ConicADMM3c converges to a KKT point of the primal DNN (3) and the dual DNN (4).Note that also the order of the updates on the blocks of variables is different with respect to ADAL+.The convergence analysis is based on the fact that ConicADMM3c is equivalent to a specific convergent 2-block ADMM.
With respect to ADAL+, ConicADMM3c has the drawback that fewer optimality conditions are satisfied up to machine accuracy throughout the algorithm.Additionally to r P , r D , r P P and r CS , the stopping criterion of ConicADMM3c has to take into account the errors related to the primal feasibility X ∈ S + n and the complementarity condition ZX = 0.In fact, as the second update of y is performed after the update of Z, the spectral decomposition of W k+1 cannot be used to update X as in ADAL+ and both the complementarity condition ZX = 0 and the positive semidefiniteness of X are not satisfied by construction.(We will give a summary on the conditions satisfied throughout the algorithms in Table 1 in a subsequent section.)From a computational point of view this slows down the convergence of the scheme, which will be confirmed in our computational evaluation in Section 5.

Dual Matrix Factorization
In this section, we present our new variants of ADAL+ and ConicADMM3c, namely DADAL+ and DADMM3c, where a factorization of the dual variable Z is employed.We adapt the method introduced by De Santis, Rendl and Wiegele in [3].In particular, we look at the augmented Algorithm 2 Scheme of ConicADMM3c from [14] 1: Choose σ > 0, ε > 0, X ∈ S + n , Z ∈ S + n , S ∈ S n with S ≥ 0 2: δ = max{r P , r D , r P P , r P D , r CS , r CZ } 3: while δ < ε do 4: δ = max{r P , r D , r P P , r P D , r CS , r CZ } 10: Update σ 11: end while Lagrangian problem where the positive semidefinite constraint on the dual matrix Z is eliminated by considering the factorization Z = V V ⊤ .To be more precise, in each iteration of the ADMMs for fixed X, we focus on the problem max L σ (y, S, V ; X) where Compared to (6) n is fulfilled automatically.Note that the number of columns r of the matrix V represents the rank of Z.
The use of the factorization of the dual variable in ADAL+ should improve the numerical performance of the algorithm when dealing with structured DNNs, as it happens in the comparison of the algorithm DADAL with ADAL when dealing with structured SDPs [3].For what concerns ConicADMM3c, we will see in Section 3.2 that using the factorization of the dual variable allows to avoid any spectral decomposition along the iterations of the algorithm, without compromising the theoretical convergence of the method.
Note that Problem ( 13) is unconstrained with respect to the variables y and V .In particular, the following holds.
Proposition 1 implies that fulfilling the necessary optimality conditions with respect to y is equivalent to solve one system of linear equations.
As in [3], we consider Algorithm 3 in order to update y and V (and hence Z) for fixed S and X.In particular in Algorithm 3, starting from (y, S, V ; X), we move V along an ascent direction D V ∈ R n×r with a stepsize α.While doing this, we update y in such a way that we keep its optimality conditions of (13) satisfied, so ∇ y L σ (y, S, V + αD V ; X) = 0 holds for the updated y (see [3,Proposition 2]).We stop as soon as the necessary optimality conditions with respect to V (see Proposition 1) are fulfilled to a certain precision.
As in the algorithm DADAL presented in [3], in our implementation we set D V either to the gradient of L σ (y, S, V ; X) or to the gradient scaled with the inverse of the diagonal of the Hessian of L σ (y, S, V ; X).In order to determine a stepsize α, at Step 4 in Algorithm 3 we could perform an exact linesearch to maximize L σ (y(V + αD V ), S, V + αD V ; X) with respect to α.This is a polynomial of degree 4 in α, so we can interpolate it from five different points in order to get its analytical expression and by this determining the maximizer explicitly.In practice we evaluate L σ (y(V + αD V ), S, V + αD V ; X) for 1000 different values of α ∈ (0, 10) and take the α corresponding to the maximum value of L σ .
Compute stepsize α 5: As output of Algorithm 3, we get y and V (and therefore also Z = V V ⊤ ) that have been updated through the maximization of the augmented Lagrangian (13) with respect to V .This leads to a new point (y, S, V ; X).
This update can be used within ADAL+ and ConicADMM3c as detailed in the following.

DADAL+
First we consider DADAL+, our version of ADAL+ where the use of the factorization of the dual variable Z leads to a double update of Z.As a further enhancement of the algorithm ADAL+ devised in [15], we propose to perform also a double update of the dual variable y.
To be more precise, we replace the first update of y in ADAL+ with a update of y and V with Algorithm 3 in DADAL+.Furthermore in DADAL+ we update y not only before, but also a second time after the computation of S. This second update is performed by applying the closed formula solution of the maximization problem in (8).Note that the second update of y is performed before the update of Z so that by computing the spectral decomposition of W = X/σ − C + A ⊤ y + S, we can simultaneously update Z and X and both the complementarity condition ZX = 0 and the positive semidefiniteness of X are satisfied up to machine accuracy throughout the algorithm in the same way it is the case in ADAL+.The scheme of DADAL+ is detailed in Algorithm 4.

DADMM3c
We now investigate the use of the dual factorization within the algorithm ConicADMM3c and call the modified algorithm DADMM3c.In ConicADMM3c, the effort spent to compute the spectral decomposition of W = X/σ − C + A ⊤ y + S is not that well exploited as it is used to update only the dual matrix Z but not the primal matrix X. Hence in DADMM3c we update Z and Algorithm 4 Scheme of DADAL+ Update (y, V ) with Algorithm 3 δ = max{r P , r D , r P P , r CS } 13: Update σ 14: end while y by employing the factorization Z = V V ⊤ and performing Algorithm 3 instead of updating them by a spectral decomposition and a closed formula as it is done in ConicADMM3c.If we assume that the update of y and V is done such that Problem ( 13) is solved to optimality, the theoretical convergence of the method is maintained.Note that the computation of any spectral decomposition is avoided.The scheme of the algorithm DADMM3c is detailed in Algorithm 5.
Algorithm 5 Scheme of DADMM3c Update (y, V ) with Algorithm 3 Update σ 12: end while A limit of DADMM3c is that the rank of Z is not updated throughout the iterations.This means that the maximization of L(y, S, V ; X) with respect to V is performed keeping r fixed to the initial value that, in our implementations, is given by the Pataki bound [11] It is still an open question to update the rank of Z in a beneficial way.
On the other hand, note that in DADAL+ the rank of Z is determined at every iteration through the eigenvalue decomposition in the second update of Z.
As already mentioned in Section 2, some of the optimality conditions are satisfied throughout the algorithms ADAL+/DADAL+ and ConicADMM3c/DADMM3c.A summary is presented in Table 1.
Table 1: Optimality conditions (5) satisfied through the algorithms by construction are indicated by a checkmark, all others by an x-mark.

Computation of Dual Bounds
When solving combinatorial optimization problems, DNN relaxations very often yield high quality bounds.These bounds can then be used within a branch-and-bound framework in order to get an exact solution method.In this section we want to discuss how we can obtain lower bounds on the optimal objective function value of the primal DNN (3) from a dual solution of moderate precision only.
Thanks to weak and strong duality results, the objective function value of every feasible solution of the dual DNN ( 4) is a lower bound on the optimal objective function value of the primal DNN (3) and the optimal values of the primal and the dual DNN coincide.Therefore every dual feasible solution and in particular the optimal dual solution give rise to a dual bound.
Note that the dual objective function value serves as a dual bound only if the DNN relaxation is solved to high precision.If the DNN is solved to moderate precision, the dual objective function value might not be a bound as the dual solution might be infeasible.However, solving the DNN to high precision comes with enormous computational costs.
So unfortunately ADAL+, DADAL+, ConicADMM3c and DADMM3c are not suitable to produce a bound fast.Running an ADMM typically gives approximate optimal solutions rather quickly, while going to optimal solutions with high precision can be very time consuming.As the dual constraint A ⊤ y + Z + S − C = 0 does not necessarily hold in every iteration of the four algorithms (see Table 1), obtaining a dual feasible solution with sufficiently high precision with ADMMs may take extremely long.
To save time, but still ensure that we obtain a dual bound, we will stop the four methods at a certain precision.After that we will use one of two procedures in a post-processing phase in order to obtain a bound.In Section 4.1 we will describe how to obtain a bound with a method already presented in the literature.In Section 4.2 we present a new procedure for obtaining a dual feasible solution and hence a bound from an approximate optimal solution.

Dual Bounds through Error Bounds
In this section we present the method to obtain lower bounds on the primal optimal value of an SDP of the form (1) introduced by Jansson, Chaykin and Keil [6].We adapt this method for DNNs in order to use it in a post-processing phase of the four ADMMs presented above.We start with the following lemma from [6, Lemma 3.1].
Lemma 1.Let Z, X be symmetric matrices of dimension n that satisfy for some z, x ∈ R. Then the inequality Proof.Let Z = QΛQ ⊤ be an eigenvalue decomposition of Z with QQ ⊤ = I for some Q ∈ R n×n and Λ = Diag(λ(Z)).Then Because of (15), we have 0 At this point we can present the following theorem of [6, Theorem 3.2] adapted for DNNs.
Theorem 1.Consider the primal DNN (3), let X * be an optimal solution and let p * be its optimal value.Given y ∈ R m and S ∈ S n with S ≥ 0, set holds.
Proof.Let X * be optimal for the primal DNN (3).Then Since S ≥ 0 and X * ≥ 0, the inequality is satisfied and Lemma 4.1 implies which proves (17).
Theorem 1 justifies to compute dual bounds via Algorithm 6.If the matrix Z defined in ( 16) is positive semidefinite, then (y, Z, S), is a dual feasible solution and b ⊤ y is already a bound.Otherwise, we decrease the dual objective function value b ⊤ y of the infeasible point (y, Z, S) by adding the negative term x λ k ( Z) to it.In this way, we obtain a bound (EB in Algorithm 6) as proved by Theorem 1.
Note that for the computation of the bound of Theorem 1 it is not necessary to have a primal optimal solution X * at hand, only an upper bound on the maximum eigenvalue of an optimal solution is needed.Such an upper bound is known for example if there is an upper bound on the maximum eigenvalue of any feasible solution.

Dual Bounds through the Nightjet Procedure
Next we will present a new procedure to obtain bounds.In contrast to the procedure described in the previous section, this approach will also provide a dual feasible solution.The key ingredient to obtain such dual feasible solutions will be the following lemma.
has an optimal solution ỹ, let S = C − Z − A ⊤ ỹ.Then (ỹ, S, Z) is a dual feasible solution.If (18) is unbounded, then also (4) is unbounded.If (18) is infeasible, then there is no dual feasible solution with Z.
Proof.If (18) has an optimal solution ỹ, then it is easy to see that S ≥ 0 by construction.Furthermore S ∈ S n because C, Z, A ⊤ y ∈ S n .Therefore (ỹ, S, Z) is a dual feasible solution.
If (18) is unbounded, then the same values of y that make the objective function value of (18) arbitrarily large can be used to make the objective function value of (4) arbitrary large, hence also (4) is unbounded.Furthermore it is easy to see that (18) is feasible if there is a dual feasible solution with Z. Hence if (18) is infeasible, then there is no dual feasible solution with Z.
Let (y, S, Z, X) be any solution (not necessarily feasible) to the primal DNN (3) and the dual DNN (4).In the back of our minds we think of them as the solutions we obtained by ADAL+, DADAL+, ConicADMM3c or DADMM3c, so they are close to optimal solutions but not necessarily dual or primal feasible.We want to obtain ỹ, S and Z satisfying dual feasibility We use Lemma 2 within the Nightjet procedure for obtaining such solutions in the following way.From the given Z we obtain the new positive semidefinite matrix Z by projecting Z onto the positive semidefinite cone.Then we solve the linear program (18).
If ( 18) is infeasible, then we are neither able to construct a feasible dual solution nor to construct a dual bound.If (18) is unbounded, then also the dual DNN ( 4) is unbounded and hence the primal DNN (3) is not feasible.If (18) has an optimal solution ỹ, then we obtain a dual feasible solution (ỹ, S, Z) with the help of Lemma 2. Furthermore the dual objective function value b ⊤ ỹ is a bound in this case, so we can return a dual feasible solution and a bound.The Nightjet procedure is detailed in Algorithm 7.

Algorithm 7 Scheme of the Nightjet Procedure
return "No dual feasible solution and no bound found" 6: end if To summarize, we have presented two different approaches to determine dual bounds for the primal DNN (3) from given y, S and Z.
Note that the approaches are in the following sense complementary to each other: In the first approach from Jansson, Chaykin and Keil we fix y and S and obtain the bound from a newly computed Z, but we do not obtain a dual feasible solution.In our second approach, the Nightjet procedure, we fix Z to be the projection of Z onto the positive semidefinite cone and then construct a feasible ỹ and S from that.
Furthermore note that in the approach of Jansson, Chaykin and Keil the obtained bound is always less or equal to the dual objective function value of y, because a negative term is added to b ⊤ y, the dual objective function value using y.In contrast to that, it can happen in the Nightjet procedure that the bound is larger and hence better than b ⊤ y.However, the Nightjet procedure comes with the drawback that it might be unable to produce a feasible solution.In this case one should continue running the ADMM to a higher precision and apply the procedure to the improved point.

Numerical Experiments
In this section we present a comparison of the four ADMMs using the two procedures presented in Section 4 as post-processing phase.Towards that end we consider instances of one fundamental problem from combinatorial optimization, the stable set problem.

The Stable Set Problem and an SDP Relaxation
Given a graph G, let V (G) be its set of vertices and E(G) its set of edges.A subset of V (G) is called stable, if no two vertices are adjacent.The stability number α(G) is the largest possible cardinality of a stable set.It is NP-hard to compute the stability number [8] and it is even hard to approximate it [5], therefore upper bounds on the stability number are of interest.One possible upper bound is the Lovász theta function ϑ(G), see for example [12].The Lovász theta function is defined as the optimal value of the SDP X ∈ S + n , where J is the n-by-n matrix of all ones.Note that ϑ(G) -as SDP of polynomial size -can be computed to arbitrary precision in polynomial time.Hence ϑ(G) is a polynomial computable upper bound on α(G).
Several attempts of improving ϑ(G) towards α(G) have been done.One of the most recent ones is including the so called exact subgraph constraints into the SDP of computing ϑ(G), which make sure that for small subgraphs the solution is in the respective squared stable set polytope [4].This approach is a generalization of one of the first approaches to improve ϑ(G) in [13], which consisted of adding the constraint X ≥ 0. Compared to ϑ(G) this leads to an even stronger bound on α(G) as the copositive cone is better approximated.We denote by ϑ + (G) the optimal objective function value of the DNN Note that in the DNN (20) the matrix AA ⊤ is a diagonal matrix, which leads to an inexpensive update of y in the methods discussed.

Dual Bounds for ϑ + (G)
As already discussed in Section 4, for a combinatorial optimization problem like the stable set problem, bounds on the objective function value are of huge importance.
The bound according to Jansson, Chaykin and Keil [6] can be used for computing bounds on ϑ + (G) very easily: We can set x = 1, as for every feasible solution X of (20) we have trace(X) = 1 and X ∈ S + n and hence λ max (X) ≤ 1.The computation of the dual bound with the Nightjet procedure simplifies drastically.In particular there is no need to solve the linear program (18), since the solution can be computed explicitly.To be more precise, the following holds.Lemma 3. We consider the primal DNN (20) to compute ϑ + (G) and the dual of it.Let y t be the dual variable for the constraint trace(X) = 1 and y e be the dual variable for the constraint X ij = 0 for every edge e = {i, j} ∈ E(G).Furthermore let Z ∈ S + n and let Proof.We first consider the dual of (20) in more detail.To be consistent with our notation we replace the objective function max J, X of (20) with the equivalent objective function − min −J, X in order to consider a primal minimization problem as in the primal DNN (3).We introduce one dual variable y t for the constraint trace(X) = 1 and one dual variable y e for the constraint X ij = 0 for every edge e = {i, j} ∈ E(G).Then the dual of (20) is given as Now we apply Lemma 2 for (20).Thus we replace the dual variable Z with the fixed Z ∈ S + n and the linear program (18) becomes Clearly this linear program is bounded and detecting infeasibility or constructing an optimal solution is straightforward.Indeed, let M = max Zij | {i, j} ∈ E(G) , then it is easy to see that ( 22) is infeasible if M > −1.If −1 < M < 0 holds, then we can redefine Z as Z = 1 −M Z, and obtain a new Z for which M = −1.If M ≥ 0, then we can not update Z in a straightforward way.If M ≤ −1, then ( 22) is feasible and we can construct the optimal solution as Then we let S = C − Z − A ⊤ ỹ and due to Lemma 2 this yields a feasible dual solution (ỹ, S, Z).
Hence, for computing a dual bound for ϑ + (G) it is not necessary to solve the linear program ( 18), but the solution of it can be written down explicitly.This explicit solution is used by the Nightjet procedure for ϑ + (G) to obtain ỹ.The computation of Z and S is the same as in the original Nighjet procedure.The pseudocode of the Nightjet procedure applied to the computation of ϑ + (G) can be found in Algorithm 8.
Algorithm 8 Scheme of the Nightjet Procedure for ϑ + (G) return "No dual feasible solution and no bound found" 5: else if 0 > M > −1 then

Comparison of the Evolution of the Dual Bounds
In the following, we give a numerical comparison of the two procedures for the computation of bounds for ϑ + (G) on one instance from the second DIMACS implementation challenge [7], namely johnson8 2 4. For this instance the stability number α(G) and ϑ + (G) coincide and both are equal to 4.
In Figure 1, we show the evolution of the bounds along the iterations for ADAL+, DADAL+, ConicADMM3c and DADMM3c.For each algorithm we report the dual objective function value (dualOfv), the bound computed according to Jansson, Chaykin and Keil [6] (EB) and the bound computed by the Nightjet procedure (N B) at every iteration.
Note that in some iterations the dual objective function value is not a bound on ϑ + (G) = 4 and hence also not on α(G).This is due to the fact that the solution considered is not dual feasible.(The criteria are satisfied only to moderate precision.) We observe that for ADAL+, DADAL+ and ConicADMM3c the Nightjet bound is always less or equal than the error bound and in several iterations it is significantly better, in particular at the iterations in the beginning.Hence our Nightjet procedure is an effective tool to obtain dual bounds.Note that every ADMM keeps Z positive semidefinite along the iterations (see Table 1) and this may be in favor of the Nightjet procedure.

Computational Setup
In our numerical experiments we compare the performance of ADAL+, DADAL+, ConicADMM3c and DADMM3c on 66 instances of the DNN (20) to compute ϑ + (G).The graphs are taken from the second DIMACS implementation challenge [7].Note that in that challenge the task was to find a maximum clique of several graphs, so we consider the complement graphs of the graphs in [7].In Table 2, for each instance on a graph G, we report its name (Problem) and its dimension (the number of vertices n and the number of edges m of G).
We implemented the four algorithms detailed in Sections 2 and 3 in MATLAB R2019a.In all computations, we set the accuracy level ε to 10 −5 and we set a time limit of 3600 seconds CPU time.In both DADAL+ and DADMM3c we perform two iterations of Algorithm 3 in order to update (y, V ).
It is known that the performance of ADMMs strongly depends on the update of the penalty parameter σ.In all implementations, we use the strategy described by Lorenz and Tran-Dinh [9], so in iteration k we set The experiments were carried out on an Intel Core i7 processor running at 3.1 GHz under Linux.

Comparison between ADAL+ and DADAL+
In Table 3 we report the results obtained with ADAL+ and DADAL+ on the 66 instances of computing ϑ + (G) detailed in Table 2.We include the following data for the comparison: For each instance, we report its name (Problem) and its stability number (α) and for each of the two algorithms, we report the dual objective function value obtained (d ofv), the bound obtained by computing the error bound described in Section 4.1 (EB), the bound obtained by applying the Nightjet procedure described in Section 4.2 (N B), the number of iterations (it) and the CPU time needed to satisfy the stopping criterion (time).
As a further comparison, we report in Figure 2 the performance profiles of ADAL+ and DADAL+ with respect to the number of iterations and the CPU time.These performance profiles are obtained in the following way.Given our set of solvers S and a set of problems P, we compare the performance of a solver s ∈ S on problem p ∈ P against the best performance obtained by any solver in S on the same problem.To this end we define the performance ratio r p,s = t p,s / min{t p,s ′ | s ′ ∈ S}, where t p,s is the measure we want to compare, and we consider a cumulative distribution function ρ s (τ ) = |{p ∈ P | r p,s ≤ τ }|/|P|.The performance profile for s ∈ S is the plot of the function ρ s .
Note that both ADAL+ and DADAL+ stopped on 7 instances because of the time limit.In the performance profiles, we exclude those instances where at least one of the solvers exceeds the time limit.
It is clear from the results on Table 3 and from the performance profiles that DADAL+ performs much less iterations than ADAL+.However, this does not always correspond to an improvement in terms of computational time as the double update of y is an expensive operation.
With respect to the CPU time, Figure 2 shows that the performance of the two algorithms is similar, even if DADAL+ slightly outperforms ADAL+ as its curve is always above the other one.
If we consider the dual objective function value in Table 3 we see that in fact the dual objective function value obtained by ADAL+ and DADAL+ is often not a bound, for example on the instances hamming6 4, c fat200 1, san200 0 7 1, san400 0 9 1, c fat500 1 and c fat500 5.This shows that a procedure for obtaining a bound from the approximate solution is indeed of major importance.
Regarding the quality of the bounds, the Nightjet procedure is able to obtain better bounds with respect to the error bounds, both when applied as post-processing phase for ADAL+ and for DADAL+, for the vast majority of the instances.The improvement is particularly impressive when looking at those instances where the time limit is exceeded.We want to further highlight that the bound obtained from the Nightjet procedure comes from a newly computed feasible dual solution.This means that applying the Nightjet procedure as post-processing does not only guarantee a bound generally better than the one obtained by the error bounds, but it also provides a dual feasible solution.

Comparison between ConicADMM3c and DADMM3c
In Table 4 we report the results obtained with ConicADMM3c and DADMM3c on the 66 instances of computing ϑ + (G) detailed in Table 2.
As before, we report the name of the instances, the stability number and, for each algorithm, the dual objective function value obtained, the bounds obtained by computing the error bound and by applying the Nightjet procedure, the number of iterations and the CPU time needed to satisfy the stopping criterion.

Conclusions
In this paper we propose to use a factorization of the dual matrix within two ADMMs for conic programming proposed in the literature.In particular we use a first order update of the dual variables in order to improve the performance of the ADMMs considered.
Termini on November 12, 2019 for being one hour delayed; this helped Elisabeth Gaar to focus and set the ground stone for the Nightjet procedure.
. A ⊤ y + Z = C Z ∈ S + n , r D , r P P , r P D , r CS , r CZ } 4: while δ < ε do 5: then it is not possible to construct a dual feasible solution with this Z.If 0 > M > −1, then we can redefine Z as Z = 1 −M Z, and obtain a new Z for which M = −1.If M ≤ −1, then we obtain a dual feasible solution with ỹt

Figure 1 :
Figure 1: Evolution of the computed bounds on the instance johnson8 2 4.