Nathan Krislock
email: nathan.krislock@inria.fr

Jérôme Malick
email: jerome.malick@inria.fr

Frédéric Roupin
email: roupin@lipn.univ-paris13.fr

Improved semidefinite bounding procedure for solving Max-Cut problems to optimality

Keywords: combinatorial optimization, semidefinite programming, quasi-Newton algorithm, triangle inequalities

We present an improved algorithm for finding exact solutions to Max-Cut and the related binary quadratic programming problem, both classic problems of combinatorial optimization. The algorithm uses a branch-(and-cut-)and-bound paradigm, using standard valid inequalities and nonstandard semidefinite bounds. More specifically, we add a quadratic regularization term to the strengthened semidefinite relaxation in order to use a quasi-Newton method to compute the bounds. The ratio of the tightness of the bounds to the time required to compute them can be controlled by two real parameters; we show how adjusting these parameters and the set of strengthening inequalities gives us a very efficient bounding procedure. Embedding our bounding procedure in a generic branch-and-bound platform, we get a competitive algorithm: extensive experiments show that our algorithm dominates the best existing method.

Introduction

Maximizing a quadratic function over the vertices of an hypercube is an important problem of discrete optimization. As far as formulation is concerned, it is the simplest problem of nonlinear (mixed-)integer programming. However, this problem is NP-hard [START_REF] Pardalos | Quadratic programming with one negative eigenvalue is NP-hard[END_REF] and it is considered a computational challenge to be solved to optimality, even for instances of moderate size.

The Max-Cut problem is one of the most famous NP-hard problems of this type, as is evidenced by its consideration in such seminal works as [START_REF] Karp | Reducibility among combinatorial problems[END_REF] and [START_REF] Goemans | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF], and in the survey chapter [START_REF] Palagi | Computational approaches to Max-Cut[END_REF] of the recent handbook [START_REF]Handbook on Semidefinite, Conic and Polynomial Optimization[END_REF]. Given a graph G = (V, E) with edge weights w ij for ij ∈ E and w ij = 0 for ij ∈ E, Max-Cut is the problem of finding a bipartition of the nodes V such that the sum of the weights of the edges across the bipartition is maximized. Let n = |V | be the cardinality of V ; we can state Max-Cut as

maximize i≤j w ij 1 -x i x j 2 subject to x ∈ {-1, 1} n . (1)
We can rewrite the problem of Max-Cut as

(MC) maximize x T Qx subject to x ∈ {-1, 1} n (2)
where the matrix Q is defined as Q := 1 4 L, and L is the Laplacian matrix of the weighted graph G; see, e.g., [START_REF] Anjos | Introduction to semidefinite, conic and polynomial optimization[END_REF].

There have been many methods for finding exact solutions of Max-Cut using semidefinite programming (SDP), including early efforts of [START_REF] Helmberg | Solving quadratic (0,1)-problems by semidefinite programs and cutting planes[END_REF], the QCR method [START_REF] Billionnet | Using a mixed integer quadratic programming solver for the unconstrained quadratic 0-1 problem[END_REF], and most recently the state-of-the-art Biq Mac method [START_REF] Rendl | Solving Max-Cut to optimality by intersecting semidefinite and polyhedral relaxations[END_REF]. This paper builds on this line of research: we present an improved branch-and-bound algorithm inspired by the above mentioned approaches and introducing new techniques. As shown by the extensive numerical experiments of this paper, our new algorithm dominates the best existing methods.

The algorithm uses a branch-and-(cut-and)-bound paradigm, using the standard valid triangle inequalities together with nonstandard semidefinite bounds. Bounds of the same type have been used recently in several papers (for general binary quadratic problems [START_REF] Malick | Spherical constraint in Boolean quadratic programming[END_REF][START_REF] Malick | On the bridge between combinatorial optimization and nonlinear optimization: a family of semidefinite bounds leading to Newton-like methods[END_REF] and specifically for the kcluster problem [START_REF] Malick | Solving k-cluster problems to optimality using adjustable semidefinite programming bounds[END_REF]). More precisely, we apply here the bounds to the Max-Cut problem that have been introduced in [START_REF] Malick | On the bridge between combinatorial optimization and nonlinear optimization: a family of semidefinite bounds leading to Newton-like methods[END_REF] for general quadratic optimization problems with quadratic and linear constraints; however, we provide a different motivation for these bounds, looking at them from a quadratic penalty point of view (see Section 3). This reveals the key role of two positive parameters, α and ε, controlling the tightness of the bound.

Our main contribution is an improved bounding procedure obtained by reducing the α and ε parameters to zero and iteratively adding triangle inequality cuts (see Section 4). We show that our bounding procedure converges to the classic reinforced semidefinite bound for the Max-Cut problem, in an efficient and robust way.

We embed our bounding procedure within a generic branch-and-bound platform for solving Max-Cut to optimality (see Section 5). We compare the overall algorithm with the best existing exact resolution method, the Biq Mac solver [START_REF] Rendl | Solving Max-Cut to optimality by intersecting semidefinite and polyhedral relaxations[END_REF] on the 328 instances in the Biq Mac Library [START_REF] Wiegele | Biq Mac Library-A collection of Max-Cut and quadratic 0-1 programming instances of medium size[END_REF] with size n < 500. To prove unambiguously that the new bounding procedure is of high practical value, we make a fair and complete computational study between the two codes: we have compiled and run both our code and the Biq Mac code (kindly provided by the authors of Biq Mac) on the same machine, and have used the same libraries and compilation flags for both codes.

We finish this section with a remark about our choice to focus on Max-Cut for presenting our approach. With the help of classic reformulation techniques, more general binary quadratic optimization problems can be recast as Max-Cut. We could, for example, add a linear term to the objective in problem [START_REF]Handbook on Semidefinite, Conic and Polynomial Optimization[END_REF], or use variables in {0, 1}. For instance, the binary homogeneous quadratic optimization problem minimize

x T Qx subject to x ∈ {0, 1} n

can be reformulated as Max-Cut by considering a change of variable and by increasing the size of the problem by one; see, e.g., [START_REF] Helmberg | Solving quadratic (0,1)-problems by semidefinite programs and cutting planes[END_REF]Section 2]. Therefore, to keep the presentation as simple as possible, we only present our developments for the Max-Cut problem [START_REF] Anderson | LAPACK Users' Guide, third edn[END_REF]. On the other hand, we consider both Max-Cut problems and binary quadratic optimization problems in our numerical experiments.

Semidefinite formulation and relaxation. The inner product of two matrices X and Y is defined and denoted by X, Y = trace(X T Y ). The notation X 0 means that X is symmetric positive semidefinite. Introducing the symmetric positive semidefinite rank-one matrix X = xx T , we observe that we can write the Max-Cut problem as (MC) maximize Q, X subject to diag(X) = e, X 0,

X = xx T , (4)
where diag(X) is the vector of the diagonal entries of the matrix X, and the vector of all ones is represented by e-for simplicity, we allow the size of the vector e to be determined from the context. The classic semidefinite relaxation is obtained by dropping the nonconvex rank-one constraint in problem ( 4):

(SDP) maximize Q, X subject to diag(X) = e, X 0.

(

Since the feasible set of problem [START_REF] Bonnans | Numerical Optimization[END_REF] is strictly larger than that of problem (4), we can easily see that the optimal value of problem [START_REF] Bonnans | Numerical Optimization[END_REF] provides an upper bound on the weight of the optimal cut:

(MC) ≤ (SDP).

There exist some theoretical results quantifying the tightness of the above inequality (see, in particular, the seminal work [START_REF] Goemans | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF]). However, in practice, the ratio of the tightness of the bound to the time needed to compute the bound is too low to allow this bound be used efficiently to solve the Max-Cut problem to optimality (see, e.g., [START_REF] Helmberg | Solving quadratic (0,1)-problems by semidefinite programs and cutting planes[END_REF][START_REF] Rendl | Solving Max-Cut to optimality by intersecting semidefinite and polyhedral relaxations[END_REF][START_REF] Wiegele | Nonlinear optimization techniques applied to combinatorial optimization problems[END_REF]).

Strengthening inequalities. In order to improve the performance of the bound, it is necessary to tighten the bounds by adding valid inequalities to the SDP relaxation [START_REF] Bonnans | Numerical Optimization[END_REF]. The most popular class of inequalities are the triangle inequalities, defined for 1 ≤ i < j < k ≤ n by

X ij + X ik + X jk ≥ -1, X ij -X ik -X jk ≥ -1, -X ij + X ik -X jk ≥ -1, -X ij -X ik + X jk ≥ -1,
(see [START_REF] Helmberg | Solving quadratic (0,1)-problems by semidefinite programs and cutting planes[END_REF]Section 4] and [START_REF] Rendl | Solving Max-Cut to optimality by intersecting semidefinite and polyhedral relaxations[END_REF]Section 3.2]). They represent the fact that for any 3-cycle of vertices in the graph, we can only have an even number of edges cut (i.e., either no edges are cut, or exactly two edges are cut). There are

4 n 3 = 4 n(n -1)(n -2) 6
triangle inequalities in total. Adding any subset of those valid inequalities can strengthen the SDP bound [START_REF] Bonnans | Numerical Optimization[END_REF]. In particular, let I be a subset of the triangle inequalities, and A I : S n → R |I| be the corresponding linear operator describing the inequalities in this subset. We define the (SDP I ) problem as

(SDP I ) maximize Q, X subject to diag(X) = e, X 0, A I (X) + e ≥ 0. (6)
Adding such valid constraints gives an improved bound on the value of the maximum cut:

(MC) ≤ (SDP I ) ≤ (SDP). (7)
Note that the dual of problem ( 6) is given by minimize e T y + e T z subject to Q -Diag(y) + A * I (z) 0, z ≥ 0, [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF] where Diag(y) is the diagonal matrix with the vector y along its diagonal, and A * I is the adjoint linear operator of A I . There is a very large number of triangle inequalities: adding all of them would make the (SDP I ) problem intractable even for moderate values of n. However, to get the best possible (SDP I ) bound, we only need to add the inequalities that are active at the optimal solution of the problem having all triangle inequalities, which we call (SDP I all ). Such a set of active inequalities is obviously unknown, but we will see in Section 4 how we approximate it (see Algorithm 1 and Theorem 1). For the moment, we just note that in practice we have to select a (small) number of promising inequalities.

In addition to the triangle inequalities, other cutting planes could be used to tighten the relaxation: hypermetric inequalities [START_REF] Helmberg | Solving quadratic (0,1)-problems by semidefinite programs and cutting planes[END_REF], general gap inequalities [START_REF] Laurent | Gap inequalities for the cut polytope[END_REF], or even a sophisticated choice based on the objective function. The nature of the valid inequalities has no impact on the theoretical developments of the next section; we will deal with a general set of inequalities I.

However, in the computational experiments, as is done in [START_REF] Rendl | Solving Max-Cut to optimality by intersecting semidefinite and polyhedral relaxations[END_REF], we only use triangle inequalities as it greatly simplifies the selection of inequalities and leads to satisfactory results.

Best existing method. Biq Mac [START_REF] Rendl | Solving Max-Cut to optimality by intersecting semidefinite and polyhedral relaxations[END_REF][START_REF] Wiegele | Nonlinear optimization techniques applied to combinatorial optimization problems[END_REF] is currently the best solver for solving Max-Cut problems to optimality. This method uses a branchand-bound approach and solves the strengthened SDP relaxation (SDP I ) in equation ( 6) with a nonsmooth optimization algorithm. By dualizing the triangle inequality constraints in problem [START_REF] Borwein | Convex Analysis and Nonlinear Optimization : Theory and Examples[END_REF], the authors of the Biq Mac solver obtain the nonsmooth convex dual function

θ I (z) = e T z + maximize Q + A * I (z), X subject to diag(X) = e, X 0, (9)
where z ∈ R |I| + is the dual variable. The value θ I (z) provides an upper bound on the value of the maximum cut for each z ∈ R |I| + ; the goal is to minimize θ I (z) to obtain the tightest such upper bound. The two main ingredients of the bounding procedure of Biq Mac are as follows:

1. the function value θ I (z) and a subgradient g ∈ ∂θ I (z) are evaluated by solving the semidefinite programming problem in equation ( 9) by a customized primal-dual interior-point method;

2. the nonsmooth convex function θ I is minimized by a bundle method [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms[END_REF].

The extensive numerical experiments of [START_REF] Rendl | Solving Max-Cut to optimality by intersecting semidefinite and polyhedral relaxations[END_REF][START_REF] Wiegele | Nonlinear optimization techniques applied to combinatorial optimization problems[END_REF] show that Biq Mac dominates all other approaches (see the six tested approaches in [START_REF] Rendl | Solving Max-Cut to optimality by intersecting semidefinite and polyhedral relaxations[END_REF]Section 7]).

New family of semidefinite bounds

We present the family of semidefinite bounds that we will use in our branchand-bound method for solving Max-Cut to optimality. By construction, these bounds are less tight than corresponding usual SDP bounds (see Section 3.1). On the other hand, Section 3.2 presents their useful properties that will lead us to the bounding procedure of the next section (Algorithm 1) whose bounds converge to the usual SDP bounds in the limit. First we need some more notation. For a real number a, we denote its nonnegative part by a + = max{a, 0}. We extend this definition to vectors as follows: if x ∈ R n , then (x + ) i = (x i ) + , for i = 1, . . . , n. For a given symmetric matrix A, we denote its positive semidefinite part by A + ; more specifically, if the eigendecomposition of A is given by A = U Diag(λ)U T , with the eigenvalues λ ∈ R n , and orthogonal matrix U ∈ R n×n , then we have

A + = U Diag(λ + )U T .
We denote similarly a -, x -, and A -.

The new semidefinite bounds

We begin with the following simple fact.

Lemma 1. Let 0 ≤ ε ≤ 1 and X 0 such that all the diagonal entries X ii lie in [1 -ε, 1 + ε]. Then (1 -ε) 2 n ≤ X 2 F ≤ (1 + ε) 2 n 2 .
Proof. Let X 0 as defined above. We have

X 2 F = ij X 2 ij = i X 2 ii + i =j X 2 ij ≥ i X 2 ii ≥ n(1 -ε) 2 ,
Now take (i, j) and extract the submatrix whose determinant is nonnegative det

X ii X ij X ij X jj = X ii X jj -X 2 ij ≥ 0. Thus X 2 ij ≤ X ii X jj ≤ (1 + ε) 2 . Also, for each i, we have X 2 ii ≤ (1 + ε) 2 , so that X 2 F = ij X 2 ij ≤ (1 + ε) 2 n 2 .
A particular case of the above result is when X is feasible in (SDP I ): we have diag(X) = e, so we take ε = 0 in Lemma 1 to get

n 2 -X 2 F ≥ 0.
Adding a multiple of this nonnegative term to the objective function, we obtain the regularized problem

(SDP α I ) maximize Q, X + α 2 n 2 -X 2 F subject to diag(X) = e, X 0, A I (X) ≥ -e. ( 10
)
Proposition 1. The following holds: for all I and for all α ≥ 0,

(MC) ≤ (SDP I ) ≤ (SDP α I ). ( 11
)
If α = 0, the second inequality is an equality. If α > 0, the two inequalities are equalities if and only if there exists a rank-one optimal solution of (SDP I ) or (SDP α I ). Moreover, we have

α ′ ≤ α ⇒ (SDP α ′ I ) ≤ (SDP α I ), (12)
and

I ′ ⊆ I ⇒ (SDP α I ′ ) ≥ (SDP α I ). ( 13
)
Proof. Inequality [START_REF] Helmberg | Solving quadratic (0,1)-problems by semidefinite programs and cutting planes[END_REF] clearly holds since the feasible set of (SDP α I ) is included in that of (SDP α I ′ ). We now prove [START_REF] Golub | Matrix Computations, 3rd edn[END_REF], from which we can deduce the remaining easily. For any feasible matrix X in problem [START_REF] Dolan | Benchmarking optimization software with performance profiles[END_REF], we have

n 2 -X 2 F ≥ 0. Therefore α ′ ≤ α yields Q, X + α ′ (n 2 -X 2 F ) ≤ Q, X + α(n 2 -X 2 F )
which gives (SDP α ′ I ) ≤ (SDP α I ). Taking now α ′ = 0 gives the second inequality in [START_REF] Goemans | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF]; the first inequality comes from [START_REF] Boyd | Convex Optimization[END_REF]. For the case of equality, note first from problem (4) that (MC) = (SDP I ) if and only if there is a rank-one optimal solution of (SDP I ). Finally, we use a corollary of [18, Theorem 1] which states that if X 0 and diag(X) = e, then X F = n if and only if rank(X) = 1, completing the proof. This proposition says, first, that the bound (SDP α I ) is strictly weaker than (SDP I ) (except in the very special situation that (SDP I ) has a rankone optimal solution), and, second, that making α smaller reduces the difference. We will come back to this second point when describing our algorithm where we gradually decrease α.

Let us now fix α and I, and let us dualize the affine constraints of (SDP α I ). We define the Lagrangian function with respect to the primal variable X 0 and the dual variables y ∈ R n and z

∈ R |I| + as L(X; y, z) := Q, X + α 2 n 2 -X 2 F + y, e -diag(X) + z, e + A I (X) = Q -Diag(y) + A * I (z), X - α 2 X 2 F + e T y + e T z + α 2 n 2 .
The associated dual function is then defined by

F α I (y, z) := max X 0 L(X; y, z). ( 14
)
By weak duality, we have an upper bound on (MC) for given (I, α, y, z) with α ≥ 0 and z ≥ 0. More precisely, we have

(MC) ≤ (SDP I ) ≤ (SDP α I ) ≤ F α I (y, z). (15)
These F α I (y, z) for given (I, α, y, z) are the bounds that we are going to use to solve Max-Cut to optimality; we study them in the next section.

Mathematical study of the bounds

F α I (y, z)
This first result gives us a useful explicit expression of the bound F α I (y, z) for α > 0. Let us consider the positive semidefinite matrix

X I (y, z) := [Q -Diag(y) + A * I (z)] + . ( 16
)
Proposition 2. Let α > 0 and let I be a set of inequalities. Then, for F α I and X I defined in equations ( 14) and ( 16), respectively, we have

F α I (y, z) = 1 2α X I (y, z) 2 F + e T y + e T z + α 2 n 2 , for all y, z. (17)
Proof.

Let M = Q -Diag(y) + A * I (z) and f (X) = M, X -α 2 X 2 F .
Since α > 0, we have that f is a concave function. Moreover, max X 0 f (X) is a strictly feasible convex optimization problem, so strong duality holds (see, e.g., [START_REF] Boyd | Convex Optimization[END_REF]); that is, max

X 0 f (X) = min S 0 max X f (X) + S, X .
Then, from the observation that min

S 0 max X f (X) + S, X = min S 0 1 2α M + S 2 F = 1 2α [M ] + 2 F ,
the remainder of the proof follows easily.

This proposition makes it clear that the bounds that we consider here correspond (up to a change of sign and notation) to the ones introduced in [START_REF] Malick | On the bridge between combinatorial optimization and nonlinear optimization: a family of semidefinite bounds leading to Newton-like methods[END_REF] for general 0-1 quadratic optimization problems (see Θ(λ, µ, α) in Theorem 2 of [START_REF] Malick | On the bridge between combinatorial optimization and nonlinear optimization: a family of semidefinite bounds leading to Newton-like methods[END_REF]). The introduction of the bounds we give here is different though: it is more direct, as it does not rely on the reformulation of the initial problem with the so-called spherical constraint (see [START_REF] Malick | Spherical constraint in Boolean quadratic programming[END_REF]).

Fixing α and I, we obtain the best bound F α I (y, z) when (y, z) is an optimal solution of the problem

(DSDP α I ) minimize 1 2α X I (y, z) 2 F + e T y + e T z + α 2 n 2 subject to y free, z ≥ 0, (18)
which is the dual problem of (SDP α I ).

Observe that when α = 0 the dual function has the usual expression

F 0 I (y, z) = e T y + e T z, if Q -Diag(y) + A * I (z) 0, +∞, otherwise,
which leads us back to the dual problem [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF]. We see that the condition Q-Diag(y)+A * I (z) 0 is exactly X I (y, z) = 0, which opens another view on the approach we consider. Indeed, from equation ( 17) the function F α I (y, z) can be viewed as a penalization of the standard dual function F 0 I (y, z). The parameter α is interpreted as a penalization parameter controlling the loss in the tightness of the bound.

In view of solving problem [START_REF] Malick | Spherical constraint in Boolean quadratic programming[END_REF], the next proposition is crucial.

Proposition 3. Let α > 0. Then the dual function F α I is convex and differentiable, with partial gradients given by

∇ y F α I (y, z) = e - 1 α diag (X I (y, z)) , ∇ z F α I (y, z) = e + 1 α A I (X I (y, z)) .
Proof. Let f : R n → R be defined by f

(x) = 1 2 x + 2 2 .
Then f is a convex and differentiable function invariant under permutations of the entries of x. The gradient of f is ∇f (x) = x + . From [6, Sect. 5.2], we have that the function g : S n → R defined by g(X) = f (λ(X)) is convex and differentiable with gradient ∇g(X) = U Diag(∇f (λ(X)))U T , where U is any n × n orthogonal matrix satisfying X = U Diag(λ(X))U T ; thus ∇g(X) = X + . Now note that

F α I (y, z) = 1 α g(Q -Diag(y) + A * I (z)) + e T y + e T z + α 2 n 2 ,
so we immediately have that F α I is convex and differentiable. Now we simply apply the chain rule, ∇ x [g(Mx)] = M * ∇g(Mx), where M : R p → S n is any linear operator, to get the desired result.

As expected, this proposition is in the same vein as [START_REF] Malick | On the bridge between combinatorial optimization and nonlinear optimization: a family of semidefinite bounds leading to Newton-like methods[END_REF]Theorem 2] establishing differentiability results for the corresponding bounds. The proof given here is original though, using the differentiability of spectral functions as a direct argument.

The final proposition considers what could happen to the value of the dual function F α I when reducing the value of α with fixed variables y and z.

Proposition 4. Let 0 ≤ ε ≤ 1 and 0 < α ′ < α. Consider (y, z) such that ∇ y F α I (y, z) ∞ ≤ ε. Then F α ′ I (y, z) ≤ F α I (y, z) + n 2 2 (α -α ′ ) α ′ α(1 + ε) 2 -α ′ and F α ′ I (y, z) ≥ F α I (y, z) + n 2 2 (α -α ′ ) α ′ α n (1 -ε) 2 -α ′ .
Proof. We have

F α ′ I (y, z) -F α I (y, z) = 1 2 1 α ′ - 1 α X I (y, z) 2 F - n 2 2 (α -α ′ ).
We apply Lemma 1 to X = 1 α X I (y, z), which satisfies max i |X ii -1| ≤ ε by Proposition 3 and the assumption

∇ y F α I (y, z) ∞ ≤ ε. This gives α 2 (1 -ε) 2 n ≤ X I (y, z) 2 F ≤ α 2 (1 + ε) 2 n 2 .
We deduce the first desired bound, as

F α ′ I (y, z) -F α I (y, z) ≤ 1 2 1 α ′ - 1 α α 2 (1 + ε) 2 n 2 - n 2 2 (α -α ′ ) = n 2 2 (α -α ′ ) α ′ α(1 + ε) 2 -α ′ .
Using the other inequality, we get in a similar way the other bound.

This result tells us that F α ′ I (y, z) may increase over F α I (y, z), but not much if αα ′ is small enough and α ′ is not too small. This gives the intuition that we should not decrease α too much or too quickly. In our numerical experiments, we take α ′ = α 2 and never take α smaller than 10 -5 (see Section 5.2). Later we will present Lemma 3 which will give us a more precise picture.

Computing the bounds

For a given set of inequalities I and a tightness level α > 0, the computation of the value of the bounding function F α I and its gradient essentially corresponds to the computation of X I (y, z), which, in turn, reduces to computing the positive eigenvalues and corresponding eigenvectors of the symmetric matrix Q -Diag(y) + A * I (z) (see equation ( 16)). To compute this partial eigenvalue decomposition, we use the LAPACK [START_REF] Anderson | LAPACK Users' Guide, third edn[END_REF] routine DSYEVR-for improved performance we use the Intel Math Kernel Library (MKL) rather than the Automatically Tuned Linear Algebra Software (ATLAS) package.

When requesting a partial eigenvalue decomposition of a matrix, DSYEVR first tridiagonalizes the matrix, then finds the positive eigenvalues using bisection, and finally uses the inverse iteration method to find the corresponding eigenvectors; see [START_REF] Golub | Matrix Computations, 3rd edn[END_REF] for more information.

In view of the differentiability result (Proposition 3), we can solve problem [START_REF] Malick | Spherical constraint in Boolean quadratic programming[END_REF] with any classic nonlinear optimization algorithm that can handle nonnegativity constraints. Among all possible algorithms and software, quasi-Newton methods are known to be efficient [START_REF] Bonnans | Numerical Optimization[END_REF][START_REF] Nocedal | Numerical Optimization[END_REF]. In our numerical experiments, we use the FORTRAN code L-BFGS-B [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF][START_REF] Zhu | Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization[END_REF] which has been recently upgraded to version 3.0 [START_REF] Morales | Remark on "Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound constrained optimization[END_REF]. We use the default parameters of the code. Note, though, that we do an initial scaling of the problem by normalizing the constraints we dualize.

In summary, the two main ingredients needed to use our bounds in practice are the following:

1. the function value F α I (y, z) and gradient ∇F α I (y, z) are evaluated by computing a single partial eigenvalue decomposition;

the smooth convex function F α

I is minimized by a quasi-Newton method. Thus, our method is built using reliable and efficient numerical codes: the eigensolver DSYEVR and the quasi-Newton solver L-BFGS-B.

Improved semidefinite bounding procedure

This section presents the bounding procedure using the bounds F α I (y, z) introduced in the previous section, provides a theoretical analysis, and finally a numerical illustration.

Description of the bounding procedure

We have three ways to improve the tightness of the bound F α I (y, z): adding triangle inequalities, reducing the parameter α, or reducing the tolerance ε in the stopping criteria of the quasi-Newton method:

1. Adding inequalities. In view of (13), we can improve the bound by adding violated (triangle) inequalities to I; that is,

find i such that A i (X) + 1 < 0 where X = 1 α X I (y, z).
Recall that we only need to add the triangle inequalities that are active at the optimal solution that satisfies all triangle inequalities. We do this by adding a predefined number of the most violated inequalities each time to improve the bound as quickly as possible. We add just a few inequalities at the beginning when we are far from the optimal solution, then we add more inequalities when we are closer to this optimal solution. Moreover, we borrow an idea from the description of Biq Mac in [START_REF] Wiegele | Nonlinear optimization techniques applied to combinatorial optimization problems[END_REF] where, instead of using the current iterate X k to find violated triangle inequalities, a convex combination of the current iterate with the previous iterate is used:

X test = λX k-1 + (1 -λ)X k .
However, while Biq Mac uses λ = 0.9, we take the opposite strategy and use λ = 0.2 since we found this works well in our code.

2. Reducing α. We can reduce α > 0 and solve problem (18) again, warm-starting with the current (y, z). In practice, we reduce α if the number of violated triangle inequalities is small for the current value of α. In other words, our strategy consists of changing the target when we can no longer make good progress by adding inequalities.

3. Reducing ε. We can request more accuracy from the quasi-Newton method by reducing the tolerance ε > 0 in the stopping test:

max ∇ y F α I (y, z) ∞ , [∇ z F α I (y, z)] -∞ < ε.
Note that the expressions of the gradients of

F α I (Proposition 3) gives max diag 1 α X I (y, z) -e ∞ , A I 1 α X I (y, z) + e -∞ < ε.
How we control the decrease of ε and α, and how we add inequalities, is important for the overall efficiency. Our algorithm, described formally in Algorithm 1, can be viewed as a succession of cutting plane methods (alternating adding cuts and reducing α and ε). We give details in Section 5.1 on the parameters we chose to efficiently interlace the decrease of α and ε with the management of the set of enforced inequalities.

Algorithm 1 Improved semidefinite bounding algorithm Input: Scalars α 1 > 0, ε 1 > 0, and 0 < ScaleAlpha, ScaleEps < 1.

y 0 ← 0 ∈ R n {initial y variables} and z 0 ← ∅ {initial z variables} I 0 ← ∅ {initial set of triangle inequalities I} for k = 1, 2, . . . do Starting from (y k-1 , z k-1
), use a descent method (e.g., a quasi-Newton method) to compute (y k , ẑk ) such that

max diag(X k ) -e ∞ , A I k-1 (X k ) + e -∞ < ε k , where X k ← 1 α k X I k-1 (y k , ẑk ). Update the bound: F k ← F α k I k-1 (y k , ẑk ).
Remove triangle inequalities that are not active:

I - k-1 ← {i ∈ I k-1 : (ẑ k ) i = 0 and A i (X k ) + 1 > ε k } .
Add triangle inequalities that are violated:

Let X test k ← λX k-1 + (1 -λ)X k (taking X 0 = X 1 ) and let i 1 , . . . , i ℓ be the indices i ∈ I k-1 such that A i (X test k ) + 1 ≤ ViolTol k , where ViolTol k ≤ 0, ordered such that A i 1 (X test k ) ≤ • • • ≤ A i ℓ (X test k ). Let I + k-1 ← {i 1 , . . . , i K } , where K = min{ℓ, MaxIneqAdded k }.
Update set of inequalities:

I k ← I k-1 \ I - k-1 ∪ I + k-1 .
Initialize multipliers for added inequalities to zero:

for each i ∈ I k , (z k ) i ← (ẑ k ) i , if i ∈ I k-1 , 0, if i ∈ I k-1 .
if the number of inequalities added, K, is less than MinIneqViol then

α k+1 ← ScaleAlpha • α k , ε k+1 ← ScaleEps • ε k else α k+1 ← α k , ε k+1 ← ε k end if end for

Convergence of the bounding procedure

We now analyze Algorithm 1, starting with the case when the set of triangle inequalities I is fixed. Lemma 2. Let I be a (fixed) set of valid triangle inequalities. Consider five sequences indexed by k:

α k > 0, ε k > 0, y k ∈ R n , z k ∈ R |I| + , X k 0 such that max ∇ y F α k I (y k , z k ) ∞ , [∇ z F α k I (y k , z k )] -∞ < ε k , (19)
and

X k := X I (y k , z k )/α k . If both (α k
) k and (ε k ) k tend to 0, and if the sequence (X k , y k , z k ) k converges, then its limit is a primal-dual solution of (SDP I ), and

lim k→+∞ F α k I (y k , z k ) = (SDP I ). ( 20
)
Proof. First note that the feasible set of problem ( 6) is closed and bounded (by n), so there exists a primal optimal solution. Moreover, problem (6) satisfies Slater's constraint qualification: the identity matrix satisfies the affine constraints (equalities and inequalities) and strictly satisfies the positive semidefinite constraint. Therefore (see, e.g., [START_REF] Boyd | Convex Optimization[END_REF]), there is no duality gap between between primal-dual problems ( 6) and ( 8), and there exists a dual optimal solution. Thus the following conditions are necessary and sufficient for optimality of (X, S, y, z):

Q -Diag(y) + A * I (z) = S (21)
X, S = 0, X 0, S 0 [START_REF] Nocedal | Numerical Optimization[END_REF] diag(X) = e, A I (X) + e ≥ 0 (23)

z T (A I (X) + e) = 0, z ≥ 0 (24)
Assume that (X k , y k , z k ) k converges to ( X, ȳ, z). Let us show that the triple satisfies the above optimality conditions. By construction of X k ,

Q -Diag(y k ) + A * I (z k ) = α k X k + S k , for all k,
where

S k := [Q -Diag(y k ) + A * I (z k )] -. The sequence (S k ) k converges to S := [Q -Diag(ȳ) + A * I (z)] -,
by continuity of the operator [•] -. Passing to the limit shows that ( X, S, ȳ, z) satisfies ( 21) since α k vanishes. By construction we also have S k , X k = 0, X k 0, and S k 0. Since the cones of positive semidefinite and negative semidefinite matrices are closed, passing to the limit, we get that X and S satisfy [START_REF] Nocedal | Numerical Optimization[END_REF].

By Proposition 3, condition [START_REF] Malick | On the bridge between combinatorial optimization and nonlinear optimization: a family of semidefinite bounds leading to Newton-like methods[END_REF] can be written as

max { diag(X k ) -e ∞ , [A I (X k ) + e] -∞ } < ε k .
Passing to the limit we have that X satisfies [START_REF] Palagi | Computational approaches to Max-Cut[END_REF]. By assumption, we have that z k ≥ 0 for all k, so we also have that z ≥ 0. Finally, we write

|z T k (A I (X k ) + e)| ≤ z k 1 A I (X k ) + e ∞ < z k 1 ε k .
Since ε k vanishes, at the limit, we have that z and X satisfy [START_REF] Pardalos | Computational aspects of a branch and bound algorithm for quadratic zero-one programming[END_REF]. Altogether, this shows that ( X, ȳ, z) is indeed a primal-dual solution.

Finally [START_REF] Lemaréchal | Semidefinite relaxations and Lagrangian duality with application to combinatorial optimization[END_REF] yields

F α k I (y k , z k ) = e T y k + e T z k + α k ( X k 2 + n 2 )/2.
Therefore we have lim k→+∞ F α k I (y k , z k ) = e T ȳ + e T z. Since there is no duality gap for the primal-dual problems ( 6) and ( 8), we conclude that (20) holds.

The parameters α and ε control the level of tightness of the bound in that smaller values give tighter upper bounds. Furthermore, we can reach the usual SDP bound, (SDP I ), in the limit as α and ε both approach zero. In practice, we interlace the decrease of these control parameters with the process of adding violated inequalities; the overall convergence follows from the next lemma. Lemma 3. Let the sequence (α k , ε k , X k , y k , ẑk , z k , I k , F k ) k be generated by Algorithm 1. Then the following holds: (i) For all k = 1, 2, . . ., we have:

X k = 1 α k X I k (y k , z k ) and F k = F α k I k (y k , z k ).
(ii) There exists γ > 0 such that:

F k+ℓ ≤ F k + γα k , for all k, ℓ. (iii) If α k → 0, then (F k ) k is a convergent sequence.
Proof. The first item follows from the observation that, based on the definition of I k and z k in Algorithm 1, we have

A * I k (z k ) = A * I k-1 (ẑ k
) and e T z k = e T ẑk ; this follows from the fact that inequalities are removed only if they have a zero multiplier (ẑ k ) i , and multipliers for added inequalities are initialized with (z k ) i = 0.

For the second item, first let r := ScaleAlpha -1 > 1. Since

∇ y F α k I k-1 (y k , ẑk ) ∞ < ε k , by Proposition 4 we have F α k+1 I k-1 (y k , ẑk ) ≤ F α k I k-1 (y k , ẑk ) + n 2 2 (α k -α k+1 ) α k+1 α k (1 + ε k ) 2 -α k+1 . ( 25
)
From the previous item, we have

F k+1 = F α k+1 I k (y k+1 , ẑk+1 ) ≤ F α k+1 I k (y k , z k ) = F α k+1 I k-1 (y k , ẑk ),
where the inequality follows from the fact that (y k+1 , ẑk+1 ) is generated using a descent method starting from (y k , z k ). Therefore, if [START_REF] Pardalos | Quadratic programming with one negative eigenvalue is NP-hard[END_REF] we have

α k+1 = ScaleAlpha • α k = r -1 α k , by inequality
F k+1 ≤ F k + n 2 2 (r -1) (1 + ε k ) 2 -r -1 α k . Since 0 < ScaleEps < 1, we have ε k ≤ ε 1 , for all k. Letting δ := n 2 (r -1) (1 + ε 1 ) 2 -r -1 /2 > 0,
we have that

F k+1 ≤ F k + δα k , if α k+1 = ScaleAlpha • α k , F k+1 ≤ F k , if α k+1 = α k , (26)
for all k. Let ℓ and k be given, and let k 1 , . . . , k p be the

p indices k ≤ k i < k + ℓ such that α k i +1 = ScaleAlpha • α k i .
Then, from repeated application of the inequalities (26), we obtain

F k+ℓ ≤ F k + δ α k 1 + α k 2 + • • • + α kp = F k + δ α k + r -1 α k + • • • + r -(p-1) α k = F k + δ 1 + 1 r + • • • + 1 r p-1 α k ≤ F k + δ r r -1 α k .
Letting γ := δr r-1 > 0, we obtain the desired result. For the third item, suppose α k → 0. Then, by the previous item, we can conclude that lim sup k→+∞

F k ≤ lim inf k→+∞ F k , so the sequence (F k ) k converges. Theorem 1. Let the sequence (α k , ε k , X k , y k , ẑk , z k , I k , F k ) k be generated by Algorithm 1. If (α k
) k and (ε k ) k both converge to zero, and ( X, ȳ, z, Ī) is an accumulation point of the sequence (X k , y k , z k , I k ) k , then ( X, ȳ, z) is a primal-dual solution of (SDPĪ), and the sequence F k converges to the classic semidefinite bound:

lim k→+∞ F α k I k (y k , z k ) = (SDPĪ). ( 27
)
Proof. Take a subsequence of (X k , y k , z k , I k ) k that converges to ( X, ȳ, z, Ī).

Since the I k are finite sets, the convergence of a subsequence to Ī means that there are infinitely many indexes k i of this subsequence such that

I k i = Ī.
Observe now that the sequence (X k i , y k i , z k i ) i satisfies the assumptions of Lemma 2 with I = Ī. Therefore, lim i→+∞

F k i = lim i→+∞ F α k i I k i (y k i , z k i ) = lim i→+∞ F α k i Ī (y k i , z k i ) = (SDPĪ).
Since α k → 0, from Lemma 3 we also have that the entire sequence (F k ) k converges, so we conclude that lim k→+∞ F k = (SDPĪ).

We now have a key observation about our bounding procedure Algorithm 1. Although, by Proposition 1, the bounds F α I (y, z), for α > 0, are weaker than the usual SDP bound (SDP I ), we have, by Theorem 1, that we can get as close as we like using Algorithm 1.

At the limit, the quality of the bounds is only governed by the selection of promising inequalities in I k . Looking at a special case, we give some intuition to the fact that the selection of Algorithm 1 approximates an "optimal" set of inequalities I: enforcing the inequalities in I would in fact enforce all the inequalities. Corollary 1. Let the assumptions of Theorem 1 hold. If the sequence (X k , y k , z k , I k ) k converges to ( X, ȳ, z, Ī) and if ViolTol k → 0, then (SDPĪ) is equal to the semidefinite bound (SDP I all ) with all the inequalities, and ( X, ȳ, z) is a primal-dual solution of (SDP I all ), where z ∈ R 4( n 3 ) is obtained from z by expanding with zeros the entries that are not indexed by Ī.

Proof. We use the assumption of the convergence of the whole sequence. For k large enough, the set I k is equal to Ī. Thus no more inequalities are added at iteration k, which means

A i (X test k ) + 1 ≥ ViolTol k , for all i ∈ Ī.
Since (X test k ) k also converges to X, we pass to the limit to get A i ( X) + 1 ≥ 0 for all i ∈ Ī. We also know by Theorem 1 that X is optimal, hence feasible, for (SDPĪ). Therefore we have that

A i ( X) + 1 ≥ 0, for any inequality i, (28)
and that X is in fact feasible for (SDP I all ). We observe now ( X, ȳ, z) satisfies the optimality conditions of (SDP I all ) (i.e., conditions [START_REF] Morales | Remark on "Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound constrained optimization[END_REF][START_REF] Nocedal | Numerical Optimization[END_REF][START_REF] Palagi | Computational approaches to Max-Cut[END_REF][START_REF] Pardalos | Computational aspects of a branch and bound algorithm for quadratic zero-one programming[END_REF] from the proof of Lemma 2). Indeed [START_REF] Wiegele | Nonlinear optimization techniques applied to combinatorial optimization problems[END_REF] implies that ( 23) is satisfied for X = X and I = I all , and so are ( 21) and ( 24) by construction of z. Thus we have (SDP I ) = (SDP I all ) and the rest of the proof follows from Theorem 1.

Numerical illustration of the bounding procedure

This section presents a numerical comparison of our bounding procedure Algorithm 1 with the one of Biq Mac [START_REF] Rendl | Solving Max-Cut to optimality by intersecting semidefinite and polyhedral relaxations[END_REF]. We take one problem of the Biq Mac Library [START_REF] Wiegele | Biq Mac Library-A collection of Max-Cut and quadratic 0-1 programming instances of medium size[END_REF]-specifically Beasley bqp250.6, the problem highlighted in [27, Figure 1]-and we plot the convergence curve for the two bounding procedures in Figure 1. The code for each solver was compiled on the same machine using the same numerical libraries. This plot is typical of the behaviour of the two solvers: our bounding procedure has a slower start since we take a large value of α at the beginning, then the convergence curves intersect, after which we see that our bounding procedure is able to compute good bounds in less time than that of Biq Mac. Therefore, the computing time needed for our solver to get to a zone of "useful bounds" is smaller than the one used by Biq Mac-this is crucial for the relative efficiency of our branch-and-bound method and is the reason for the good numerical results of the next section.

In Figure 1, one can see large decreases in the value of the dual function on the curve for our bounding procedure: this corresponds to the iterations when α is decreased. We also note that our bounding procedure is able to attain a tighter bound than the one of Biq Mac. This may be due to a combination of: a different selection of inequalities, and the slow convergence of the bundle method near the end of the solving process.

Solving Max-Cut to optimality

In this section, we describe our implementation of a branch-and-bound method for solving Max-Cut to optimality. The scheme of our algorithm is simple and follows the usual branch-and(-cut-and)-bound paradigm-the novelty of our approach is essentially our bounding procedure that we described in Algorithm 1. Finally, we provide a complete numerical comparison to the leading Biq Mac solver [START_REF] Rendl | Solving Max-Cut to optimality by intersecting semidefinite and polyhedral relaxations[END_REF].

Branch-and-bound implementation

While Biq Mac uses a dedicated code, we embed our bounding procedure in a generic code: we use the BOB Branch & Bound platform [START_REF] Cun | The PNN Team: BOB: A unified platform for implementing branch-and-bound like algorithms[END_REF], which provides an easy and flexible way to implement a branch-and-bound algorithm. The BOB platform only requires the user to implement the following: (1) a bounding procedure, (2) a method for generating a feasible solution, and (3) a method for generating subproblems. Here we give some details about these points.

We have already described our bounding procedure in Section 4, but we would like to mention here how we decide when to stop our bounding procedure within the branch-and-bound method. Suppose the value of the best-known solution for the current subproblem is given by β. If we are able to compute a bound that is below β, we can stop our bounding procedure and prune this subproblem from the branch-and-bound tree. Indeed, in our method, stopping the run of the quasi-Newton solver at any time gives the valid bound F α k I k (y k , z k ). Similarly, Biq Mac can stop their feasible interiorpoint method at any time and get a valid bound; this is because Biq Mac uses the following form of their dual function θ I (z) = e T z + minimize e T y subject to Diag(y) Q + A * I (z), implying that (M C) ≤ θ I (z) ≤ e T z + e T y, for any feasible vector y.

We can also stop the bounding procedure when we detect that we will not be able to reach the value of β in a reasonable time. Biq Mac estimates the time to reach β using a linear approximation to the time-bound curve [START_REF] Wiegele | Nonlinear optimization techniques applied to combinatorial optimization problems[END_REF]. In our method, we fit a hyperbolic function to the time-bound curve, given its hyperbolic nature that can be observed in Figure 1. If this hyperbola has a horizontal asymptote that is far above the value of β, we terminate the bounding procedure and add the current subproblem, together with its computed bound and fractional solution, to the list of subproblems on which we will branch. We also terminate the bounding procedure if the difference between consecutive bounds, F k-1 -F k , is less than one, we are still at least two away from the lower bound β, and α k is very small (on the order of 10 -4 ).

In all subproblems in the branch-and-bound tree, except the root problem, we wait until the termination of the bounding procedure to compute a feasible solution (upper-bound) using the Goemans-Williamson (GW) random hyperplane algorithm [START_REF] Goemans | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF] then local swapping, as is done in Biq Mac. However, we do not use the technique used in Biq Mac of repeatedly running the GW algorithm (while progress is still made) on a matrix obtained by shifting the current matrix X in the direction xx T , where x ∈ {-1, 1} n is the best cut found so far. In the root problem, to get an upper-bound β, we run the GW procedure up to two times during the bounding procedure -this value of β is then used to inform us when we should terminate the bounding procedure in the root problem.

For generating subproblems, Biq Mac considers two branching strategies: "easy first" (R2), and "difficult first" (R3); for details, see [START_REF] Rendl | Solving Max-Cut to optimality by intersecting semidefinite and polyhedral relaxations[END_REF]. We use two different but similar "easy first" and "difficult first" branching strategies in the {0, 1}-model: we branch on the variable in the first row/column of X with fractional value that is furthest from, or closest to, 1 2 , respectively. For the sake of simplicity, we also refer to our "easy first" and "difficult first" branching strategies as (R2) and (R3), respectively. In our numerical results, we compare our method using (R2) branching and our method using (R3) branching to Biq Mac using (R2) branching and Biq Mac using (R3) branching-in total, we compare four methods.

Finally, we note that we have used a best-first exploration strategy in that we always take the subproblem from the current list of subproblems that has the greatest computed bound.

Description of the numerical tests

Comparison. There exist several efficient methods for solving Max-Cut and binary quadratic programming problems, such as [START_REF] Pardalos | Computational aspects of a branch and bound algorithm for quadratic zero-one programming[END_REF] and [START_REF] Billionnet | Using a mixed integer quadratic programming solver for the unconstrained quadratic 0-1 problem[END_REF]. However, the numerical tests of [START_REF] Rendl | Solving Max-Cut to optimality by intersecting semidefinite and polyhedral relaxations[END_REF] show that the results of the other methods are dominated by Biq Mac, so we only compare our solution method to Biq Mac.

Test problems. We use the test problems in the Biq Mac Library [START_REF] Wiegele | Biq Mac Library-A collection of Max-Cut and quadratic 0-1 programming instances of medium size[END_REF], which consists of both Max-Cut problems and binary quadratic optimization problems. Some are randomly generated instances, others come from a statistical physics application. We refer to [START_REF] Rendl | Solving Max-Cut to optimality by intersecting semidefinite and polyhedral relaxations[END_REF]Section 6] for a description of the data set. Since solving the instances in the Biq Mac Library with n = 500 is beyond the reach of current solvers (including Biq Mac and our method) we restrict our tests to the 328 instances in the Biq Mac Library with n < 500.

Machine. In our tests we used a Dell T-7500 (using a single core) with 4 GB of memory running the Linux operating system. We implemented our algorithm in C / FORTRAN and have used the Intel Math Kernel Library (MKL) for the eigenvalue computations. We have compiled and ran both our code and the Biq Mac code (kindly provided by the authors) on the same machine, and have used the same libraries (i.e., MKL) and compilation flags for both codes.

Parameters. We used Biq Mac with its default parameters, except for using both the non-default (R2) and the default (R3) branching strategies; additionally, we changed the time limit to 100 hours. We used the following values for Algorithm 1 in our tests:

1. Reducing α. We use ScaleAlpha = 0.5 and start with α 1 = 10; however, if no inequalities are added after the first iteration in the root problem (i.e., |I 1 | = 0 or |I 2 | = 0), we start with a smaller α 1 .

We continue to reduce α k until the minimum value of α k = 10 -5 .

2. Reducing ε. With (R3) branching, we use ε 1 = 0.08 and ScaleEps = 0.93, with a minimum value of ε k = 0.05. When using (R2) branching, we found that it is better to be less aggressive in decreasing ε, so we take ε 1 = 0.2 and ScaleEps = 0.95, with a minimum value of ε k = 0.1.

3.

Inequalities. We have used MaxIneqAdded k = 20k (resp. 30k) and MinIneqViol = 30 (resp. 60) for n < 150 (resp. for n ≥ 150). This means that we usually add less inequalities than Biq Mac.

Performance profiles. We plot the results using performance profiles [START_REF] Dolan | Benchmarking optimization software with performance profiles[END_REF], that we briefly describe here. Let P be the set of problems used to benchmark the solvers. For each problem p ∈ P, we define t min p as the minimum time required to solve p over all the solvers. Then, for each solver, we consider the performance profile function θ, which is defined as

θ(τ ) = 1 |P| p ∈ P : t p ≤ τ t min p , for τ ≥ 1, ( 29
)
where t p is the time required for the solver to solve problem p. The function θ is therefore a cumulative distribution function, and θ(τ ) represents the probability of the solver to solve a problem from P within a multiple τ of the minimum time required by all solvers considered.

Computational results

We report aggregated results from the 328 test-problems and almost 1600 hours of computing time: Figure 2 gives the performance profile, Table 1 counts the number of times a solver attains the minimum solution time, and Table 2 summarizes the computing times. For the full tables of numerical results, see Appendix A.

Compared to Biq Mac, our algorithm performs well for a large part of the problems: 241 out of the 328 problems are solved strictly faster by using our solver, which is around 75% of the test-problems. When considering Table 1: The number of times a solver attains the minimum solution time for each set of problems from the Biq Mac Library [START_REF] Wiegele | Biq Mac Library-A collection of Max-Cut and quadratic 0-1 programming instances of medium size[END_REF] having n < 500.

Problem sets

Biq Mac Our method (R2) (R3) (R2) (R3) bqp50 (Table 3 only the "hard problems" (those that Biq Mac does not solve at the root node), this percentage increases to 85% (226 out of 269).

Table 2 shows that, on some problems (e.g., g05 or w100) our method is roughly equivalent to Biq Mac, and on some others it is two times quicker (e.g., gk) or even 6 times quicker (e.g., be250.d10). In particular, we obtain substantial improvements, in both time and number of nodes, for large-sized problems, as can be seen in the detailed results in Appendix A; see, e.g., bqp250 (Table 5), be250.d10 (Table 18), and ising300 (Table 36). Additionally, we note that both our methods are able to solve all test problems within the 100 hour time limit (Biq Mac (R3) fails to solve 6 instances within this time limit).

The performance profile in Figure 2 shows that both of our methods dominate both of the Biq Mac methods in terms of speed, and also robustness, since the curves for our methods are constantly above the ones for Biq Mac.

Conclusions

In this paper, we presented an improved semidefinite bounding procedure to efficiently solve Max-Cut and binary quadratic programming problems to optimality. We precisely analyzed its theoretical convergence properties, and we conducted a complete numerical comparison with the leading Biq Mac method. Our algorithm is shown to be often faster and more robust than Biq Mac; in particular, our solver was able to solve around 75% of the problems of the Biq Mac Library (with n < 500) in less time than Biq Mac.

Moreover, there is still room for improvement in our current implementation. In particular, we would like to investigate how we can make the eigenvalue decomposition more efficient by somehow using prior information. We could also try decoupling the optimization of the y and z dual variables, as is done in Biq Mac, to possibly make the bound computation more efficient.

By focusing on Max-Cut and binary quadratic programming problems, we have shown unambiguously the interest of our semidefinite-based method to solve these classes of problems to optimality. Because our approach is based on the general bounds of [START_REF] Malick | On the bridge between combinatorial optimization and nonlinear optimization: a family of semidefinite bounds leading to Newton-like methods[END_REF], it can also be applied to binary quadratic problems with additional linear or quadratic constraints. In our future work, we will consider this generalization. Due to the inherent flexibility of our method, we believe that our method has a strong potential to handle other classes of problems, even those for which semidefinite-based methods are not yet competitive.

A Tables of numerical results

Table 3: CPU time (s) and nodes to solve the bqp50 problems.

Figure 1 :

 1 Figure 1: Comparison of the bounding procedure of Biq Mac [27] and our bounding procedure Algorithm 1 on the Beasley bqp250.6 problem. The behaviour shown here is typical.

Figure 2 :

 2 Figure 2: Performance profile curves (29) for the four solvers on the 328 problems.

Table 2 :

 2 Minimum / mean / maximum CPU time (s) for Biq Mac and our method to solve problems in the Biq Mac Library[START_REF] Wiegele | Biq Mac Library-A collection of Max-Cut and quadratic 0-1 programming instances of medium size[END_REF].

	Problem sets	Min	Biq Mac Mean	Max	Min	Our method Mean	Max
	bqp50	0.11	0.70	5.03	0.17	0.22	0.26
	bqp100	0.92	3.77	22.29	1.40	2.58	11.12
	bqp250	526.18	7052.28 54277.92	47.13 1736.11 14714.12
	gkaa	0.06	0.55	1.37	0.06	0.65	2.11
	gkab	1.11	1650.88 10719.44	0.10	787.00	5104.21
	gkac	0.16	0.62	1.02	0.13	0.85	1.95
	gkad	1.53	32.98	79.33	1.63	18.06	55.72
	gkae	126.98	6116.31 25721.41	21.77 3897.82 16654.15
	be100.d100	16.70	42.88	108.69	3.79	33.51	114.88
	be120.d30	4.52	32.64	113.27	5.01	17.68	76.30
	be120.d80	39.12	146.09	211.55	32.60	142.54	268.07
	be150.d30	39.21	367.92	844.33	7.74	285.44	767.08
	be150.d80	441.62	506.15	561.98	363.55	500.04	608.20
	be200.d30	1079.71	8223.51 43850.81	102.61 4852.01 32612.82
	be200.d80	1210.52 10692.92 35247.12	596.22 6759.69 25518.98
	be250.d10	328.64	2437.95	3754.20	27.98	409.71	905.98
	g05.n60	0.26	5.74	13.27	0.67	7.42	24.84
	g05.n80	6.39	63.09	274.50	2.80	63.10	296.10
	g05.n100	93.45	803.88	4197.29	96.05	721.25	3382.27
	pm1s80.d090	0.58	4.76	23.96	0.96	3.92	13.40
	pm1s100.d010	1.15	49.87	130.28	1.85	34.70	88.48
	pm1d80.d090	24.56	71.90	212.47	13.74	56.31	138.19
	pm1d100.d090	106.24	945.40	2868.06	56.94	620.76	1800.63
	w100.d010	1.37	23.15	131.45	1.71	23.47	152.09
	w100.d050	109.99	563.27	1152.49	81.02	475.32	1061.46
	w100.d090	38.67	836.25 2945.89	26.50	997.06	4675.78
	pw100.d010	1.64	65.34	228.70	2.06	34.88	113.86
	pw100.d050	122.67	715.95 1798.61	81.34	732.48	2506.69
	pw100.d090	209.02	606.82	1297.98	201.61	509.48	1061.18
	ising100	7.74	12.76	19.86	2.62	3.04	3.87
	ising150	42.55	60.31	85.74	7.74	10.70	12.89
	ising200	150.39	233.99	403.19	19.72	27.26	37.59
	ising250	165.79	981.56	2161.49	47.88	138.27	406.84
	ising300	991.28	4268.83 11858.56	97.43	455.77	1628.66
	t2g10	1.63	2.22	3.04	3.50	4.14	5.30
	t2g15	56.55	66.36	73.34	38.87	42.45	46.03
	t2g20	1014.38	4676.45 11132.18	351.37 1159.24	2748.55
	t3g5	3.22	4.42	5.49	7.18	7.39	7.71
	t3g6	52.34	159.36	367.12	41.85	218.22	566.66
	t3g7	81.37	5931.61 17601.56	123.46 3838.06 11234.02

Table 4 :

 4 CPU time (s) and nodes to solve the bqp100 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	1	2.40	2.40	1	1	1.40	1.72	1	1
	2	2.80	2.80	1	1	9.54	2.06	1	1
	3	1.14	1.14	1	1	7.39	1.60	3	1
	4	0.92	0.92	1	1	3.08	1.42	1	1
	5	1.43	1.43	1	1	3.99	1.69	1	1
	6	139.09	22.29	169	7	40.10	11.12	13	3
	7	1.87	1.87	1	1	5.60	1.80	1	1
	8	1.62	1.62	1	1	1.83	1.63	1	1
	9	1.11	1.11	1	1	3.10	1.46	1	1
	10	2.12	2.12	1	1	5.84	1.67	1	1
	Total:	5/10 5/10			1/10	4/10

Table 5 :

 5 CPU time (s) and nodes to solve the bqp250 problems.

	Prob.	Biq Mac time (R2)	(R3)	Biq Mac nodes (R2) (R3)	Our time (R2)	(R3)	Our nodes (R2) (R3)
	1	3709.00	1790.82	291	41	660.04	171.08	37	3
	2	3094.08	1063.38	173	25	502.86	166.32	31	3
	3	4220.53	526.18	309	11	65.49		47.13	1	1
	4	3425.14	1615.88	243	37	53.04	185.95	1	3
	5	3590.11	727.07	293	17	2122.92		55.17	101	1
	6	3937.78	13014.86	309	323	2953.55	1213.87	331	19
	7	3872.44	1119.27	297	25	917.86	173.05	55	3
	8	54277.92	>360000.00	4341	>15387	40039.83	14714.12	5571	307
	9	3398.68	3011.68	341	73	3241.21	373.52	263	7
	10	3363.66	2452.87	233	59	2651.67	393.81	189	7
	Total:	0/10		0/10			1/10		9/10

Table 6 :

 6 CPU time (s) and nodes to solve the gkaa problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	1	1.03	1.03	1	1	0.25	0.33	1	1
	2	0.19	0.19	1	1	0.58	0.42	1	1
	3	1.37	1.37	1	1	1.74	0.93	1	1
	4	0.82	0.82	1	1	1.45	1.07	1	1
	5	0.18	0.18	1	1	0.29	0.29	1	1
	6	0.08	0.08	1	1	0.07	0.06	1	1
	7	0.06	0.06	1	1	0.06	0.07	1	1
	8	0.68	0.68	1	1	3.54	2.11	1	1
	Total:	5/8	5/8			2/8	2/8

Table 7 :

 7 CPU time (s) and nodes to solve the gkab problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2)	(R3)	Our nodes (R2) (R3)
	1	1.11	1.11	1	1	1.83		0.10	13	1
	2	5.97	9.54	43	33	4.76		2.18	29	3
	3	12.69	33.03	47	83	11.95		9.44	53	13
	4	29.59	108.85	89	173	24.36	18.42	79	27
	5	66.27	300.14	161	335	39.20		42.15	85	43
	6	131.01	1044.74	175	853	86.56		94.96	163	67
	7	545.12	2509.01	551	1459	212.73	266.60	331	113
	8	1630.28	11068.75	1623	6009	851.14	789.10	1379	331
	9	3367.32	30173.75	3609	14173	2008.11	1608.04	3083	703
	10	10719.44 178991.20	5881	43609	5104.21	7968.24	6309	2107
	Total:	0/10	0/10			4/10		6/10

Table 8 :

 8 CPU time (s) and nodes to solve the gkac problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	1	0.16	0.16	1	1	0.13	0.15	1	1
	2	0.17	0.17	1	1	0.27	0.25	1	1
	3	0.38	0.38	1	1	0.44	0.43	1	1
	4	0.97	0.97	1	1	0.73	0.86	1	1
	5	0.72	0.72	1	1	1.36	0.93	1	1
	6	1.02	1.02	1	1	1.97	1.54	1	1
	7	0.94	0.94	1	1	2.99	1.95	1	1
	Total:	5/7	5/7			2/7	0/7

Table 9 :

 9 CPU time (s) and nodes to solve the gkad problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	1	1.53	1.53	1	1	1.96	2.10	1	1
	2	15.93	9.59	5	3	1.63	2.65	1	1
	3	34.73 11.38	13	3	1.85	2.79	1	1
	4	124.40 35.21	77	9	3.25	4.23	1	1
	5	112.72 74.91	129	21	86.70	55.72	125	11
	6	73.86 10.09	31	3	2.14	3.03	1	1
	7	118.68 79.33	123	25	84.20	46.04	111	9
	8	57.30 11.26	21	3	3.08	2.94	1	1
	9	127.02 50.54	139	15	100.76	34.08	103	7
	10	124.43 45.92	117	13	114.69	30.97	93	5
	Total:	1/10	1/10			4/10	5/10

Table 10 :

 10 CPU time (s) and nodes to solve the gkae problems.

	Prob.	Biq Mac time (R2)	(R3)	Biq Mac nodes (R2) (R3)	Our time (R2)	(R3)	Our nodes (R2) (R3)
	1	1391.44	126.98	129	5	285.89		21.77	37	1
	2	1902.51	13839.63	413	665	1047.08	2581.48	373	87
	3	1436.94	25243.00	241	289	896.16	1159.87	237	43
	4	1393.70	3876.28	311	177	928.96	869.92	293	31
	5	25721.41	>360000.00	4717	>19270	16654.15 41480.79	4977	1705
	Total:	0/5		0/5			3/5		2/5

Table 11 :

 11 CPU time (s) and nodes to solve the be100.d100 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	1	145.63	19.84	107	5	8.07	11.06	1	1
	2	133.79	18.53	83	5	3.79	6.27	1	1
	3	113.25	16.70	61	5	4.83	5.73	1	1
	4	134.86	66.40	153	17	214.01	40.08	101	5
	5	108.61	37.96	89	11	206.99	25.98	55	3
	6	123.99	17.93	69	5	8.13	6.82	1	1
	7	122.58	30.69	109	9	203.55	22.37	79	3
	8	115.35	108.69	141	31	157.96 114.88	119	15
	9	120.99	89.84	153	27	127.67	84.72	107	11
	10	90.89	22.22	53	7	154.37	23.57	57	3
	Total:	0/10	2/10			3/10	5/10

Table 12 :

 12 CPU time (s) and nodes to solve the be120.d30 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	1	223.63 113.27	161	19	432.05	76.30	125	5
	2	92.97	15.33	23	3	8.99	7.18	1	1
	3	120.97	16.76	31	3	11.50	7.76	1	1
	4	248.71	40.07	137	7	14.49	12.88	1	1
	5	147.27	25.69	43	5	6.25	8.12	1	1
	6	40.13	14.98	9	3	7.25	6.57	1	1
	7	6.10	6.10	1	1	9.42	6.77	1	1
	8	4.52	4.52	1	1	5.01	6.62	1	1
	9	168.79	64.48	67	11	293.80	38.91	61	3
	10	164.08	25.17	55	5	9.13	12.15	1	1
	Total:	2/10	2/10			2/10	6/10

Table 13 :

 13 CPU time (s) and nodes to solve the be120.d80 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	1	194.02 605.39	169	121	176.44	387.06	155	35
	2	211.55 211.99	151	43	210.19	155.44	127	17
	3	221.43 117.90	135	21	265.80	114.50	131	9
	4	246.84	39.12	163	7	284.60	34.34	37	3
	5	207.09	43.70	91	7	120.71	32.60	19	3
	6	218.58 120.64	137	23	284.45	97.73	109	7
	7	185.74 402.85	171	75	239.83	290.84	163	27
	8	194.87 896.84	173	191	268.07	500.34	159	47
	9	234.79 182.72	171	35	278.10	177.83	137	15
	10	241.88 170.65	183	31	323.77	128.65	157	9
	Total:	2/10	0/10			1/10	7/10

Table 14 :

 14 CPU time (s) and nodes to solve the be150.d30 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	1	552.72	170.68	145	15	784.90	111.72	101	7
	2	487.97	162.74	125	15	750.87	90.88	125	7
	3	229.30	39.21	29	3	11.00	7.74	1	1
	4	656.06	111.98	181	9	27.47	32.53	1	3
	5	560.44	344.08	193	31	769.36	101.60	181	7
	6	529.45 2214.69	215	219	511.70	487.36	219	35
	7	572.02	340.56	163	29	703.18	128.91	121	9
	8	607.62 3461.01	243	359	665.17	909.82	243	63
	9	844.33 5863.95	231	671	767.08 1072.17	227	87
	10	528.52 4004.32	253	427	466.50 1004.20	225	65
	Total:	1/10	0/10			3/10	6/10

Table 15 :

 15 CPU time (s) and nodes to solve the be150.d80 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	1	499.99 1207.49	225	129	483.87	593.74	195	49
	2	536.48 7004.88	243	715	475.02 1872.56	233	155
	3	441.62 3404.46	257	329	443.02 1019.10	239	71
	4	526.56	591.11	167	59	468.46	363.55	151	27
	5	482.92 1604.64	255	137	542.07	572.02	237	39
	6	501.99 2033.05	223	201	536.32	802.06	205	59
	7	525.18 5880.48	273	565	608.20 1211.89	217	95
	8	507.49 1963.59	209	189	522.82	760.43	205	53
	9	477.27 3903.36	253	413	544.90 1615.86	235	121
	10	561.98 1634.39	217	157	480.67	568.80	207	43
	Total:	6/10	0/10			3/10	1/10

Table 16 :

 16 CPU time (s) and nodes to solve the be200.d30 problems.

	Prob.	Biq Mac time (R2)	(R3)	Biq Mac nodes (R2) (R3)	Our time (R2)	(R3)	Our nodes (R2) (R3)
	1	43850.81	>360000.00	7927	>25405	32612.82 56902.36	8985	2399
	2	2553.97	87629.96	461	2627	1479.82 10624.76	503	355
	3	1570.50	9464.23	305	351	798.46	1345.49	285	43
	4	2597.34	90241.20	527	2251	844.91	7671.43	369	253
	5	1630.94	52812.92	339	1121	808.48	4523.10	343	155
	6	2054.03	34443.12	395	801	778.60	3057.63	293	119
	7	1571.54	1079.71	195	25	642.16	102.61	69	3
	8	4503.36	83747.75	685	2097	1135.38	5076.93	373	209
	9	7737.93	74078.57	1217	2639	2551.76	6329.74	693	251
	10	14656.52	319573.82	2359	11431	7407.25 22016.54	1935	861
	Total:	0/10		0/10			9/10		1/10

Table 17 :

 17 CPU time (s) and nodes to solve the be200.d80 problems.

	Prob.	Biq Mac time (R2)	(R3)	Biq Mac nodes (R2) (R3)	Our time (R2)	(R3)	Our nodes (R2) (R3)
	1	1210.52	25458.48	307	951	1572.38	4581.21	321	137
	2	35247.12	>360000.00	5907	>34105	25518.98 81821.55	6581	3597
	3	26318.86	>360000.00	6069	>23691	14474.43 64771.02	6067	2789
	4	1866.05	77647.69	403	1785	1062.46	7627.42	389	305
	5	9189.31	241577.21	1671	6189	4949.74 17366.65	1417	743
	6	1580.43	3643.09	355	87	1051.23	596.22	239	21
	7	1344.09	77552.67	345	1963	712.30	6756.64	333	265
	8	25950.40	>360000.00	4941	>14429	16204.96 44663.86	5053	1895
	9	1499.53	11130.84	259	211	930.82	935.36	253	37
	10	2722.89	111291.76	509	3035	1574.58	8805.80	441	363
	Total:	1/10		0/10			8/10		1/10

Table 18 :

 18 CPU time (s) and nodes to solve the be250.d10 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	1	3703.05	2098.62	319	49	2473.27	182.31	255	3
	2	3611.26	1516.95	277	35	3618.13	176.91	367	3
	3	3996.90	896.11	187	19	34.85	51.04	1	1
	4	4157.19	1522.45	353	35	1518.58	151.37	153	3
	5	3754.20	4516.16	305	113	2731.70	903.37	361	15
	6	3471.50	3285.36	225	83	2440.30	432.71	215	7
	7	2394.79	328.64	81	7	27.98	39.76	1	1
	8	3675.89	3826.01	295	87	2980.22	529.83	307	9
	9	3636.68	3992.81	311	97	3243.97	905.98	411	13
	10	3664.60	3708.73	451	85	3920.43	751.78	443	13
	Total:	0/10	0/10			2/10	8/10

Table 19 :

 19 CPU time (s) and nodes to solve the g05.n60 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	0	6.91	3.34	15	5	8.74	2.15	19	3
	1	1.98	2.02	3	3	0.89	1.00	1	1
	2	9.08	6.87	17	11	7.08	7.83	23	11
	3	0.26	0.26	1	1	0.67	0.99	1	1
	4	13.27	21.45	43	41	24.84	28.31	51	31
	5	0.34	0.34	1	1	5.89	5.00	5	3
	6	11.04	8.42	29	15	12.08	9.31	37	11
	7	13.36	4.48	33	7	7.52	12.83	15	7
	8	7.95	5.23	17	9	6.98	4.95	19	7
	9	13.26	14.39	43	23	15.35	11.80	39	25
	Total:	3/10	5/10			1/10	3/10

Table 20 :

 20 CPU time (s) and nodes to solve the g05.n80 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	0	49.16	66.56	77	61	55.22	43.71	89	47
	1	14.25	6.39	13	5	2.80	3.66	1	1
	2	34.09	20.54	43	17	18.88	22.24	45	15
	3	274.50 512.58	385	533	296.10 386.78	571	355
	4	40.71	52.76	69	45	43.36	36.54	79	29
	5	52.10	53.65	71	49	67.16	75.93	109	77
	6	36.44	40.81	59	33	42.73	36.54	71	31
	7	24.45	18.48	29	15	28.41	18.26	43	17
	8	51.92	64.20	71	61	73.96	37.41	131	47
	9	80.64 141.93	113	141	119.53	73.64	199	87
	Total:	3/10	0/10			2/10	5/10

Table 21 :

 21 CPU time (s) and nodes to solve the g05.n100 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2)	(R3)	Our nodes (R2) (R3)
	0	665.41	2194.05	547	1301	1060.48	663.41	1105	417
	1	4197.29	14856.75	3533	9565	4527.49	3382.27	5345	2293
	2	151.21	448.26	125	237	291.61	209.21	399	113
	3	1540.31	4409.16	1157	2803	2143.00	1304.29	2533	859
	4	93.45	133.74	75	67	104.66	96.05	105	51
	5	126.56	339.26	117	171	325.38	182.40	313	97
	6	197.72	422.34	181	227	346.61	237.33	355	131
	7	351.58	1176.15	297	633	479.94	540.45	567	283
	8	344.37	904.24	261	509	643.44	331.47	719	257
	9	370.94	852.80	279	503	599.54	326.16	671	275
	Total:	5/10	0/10			0/10	5/10

Table 22 :

 22 CPU time (s) and nodes to solve the pm1s80.d090 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	0	0.65	0.65	1	1	1.64	1.19	1	1
	1	0.58	0.58	1	1	1.63	1.50	1	1
	2	26.71	23.96	27	17	26.74	13.40	41	13
	3	1.19	1.19	1	1	1.73	1.90		1
	4	0.65	0.65	1	1	0.96	2.33	1	1
	5	0.90	0.90	1	1	1.05	1.22	1	1
	6	14.74	12.53	13	9	7.88	6.61	13	3
	7	1.33	1.33	1	1	7.41	4.41	11	3
	8	17.59	4.64	15	3	3.22	8.08	3	7
	9	1.13	1.13	1	1	5.10	6.20	1	1
	Total:	7/10 7/10			1/10	2/10

Table 23 :

 23 CPU time (s) and nodes to solve the pm1s100.d010 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	0	26.44	27.62	15	15	55.87	16.78	49	13
	1	130.28 143.84	93	93	190.01	88.48	193	47
	2	18.80	20.78	9	9	24.77	12.27	29	9
	3	100.00 110.81	59	59	286.00	67.01	305	45
	4	80.33	88.18	51	51	63.50	39.82	75	25
	5	27.45	30.03	15	15	73.63	33.63	89	23
	6	80.07	88.12	51	51	137.16	71.86	171	47
	7	1.15	1.15	1	1	1.85	2.55	1	1
	8	10.88	12.16	5	5	34.03	5.95	29	3
	9	23.26	25.69	11	11	9.39	13.08	5	5
	Total:	2/10	1/10			1/10	7/10

Table 24 :

 24 CPU time (s) and nodes to solve the pm1d80.d090 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	0	33.17	27.54	51	23	37.43	15.07	59	17
	1	111.02 293.45	173	275	114.29 171.07	243	131
	2	64.88 131.77	99	113	43.10	55.73	113	57
	3	30.36	31.32	39	27	30.85	18.44	45	19
	4	37.74	46.06	65	37	31.43	44.47	75	45
	5	212.47 419.49	307	405	138.19 177.21	365	185
	6	105.25 114.51	147	109	67.88	58.97	149	65
	7	64.04 211.43	111	183	109.78	98.84	253	111
	8	41.18	44.06	59	37	30.98	33.19	65	37
	9	29.88	24.56	35	21	28.31	13.74	47	15
	Total:	2/10	0/10			4/10	4/10

Table 25 :

 25 CPU time (s) and nodes to solve the pm1d100.d090 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2)	(R3)	Our nodes (R2) (R3)
	0	1229.18	8419.09	911	2281	1414.77	1026.52	1507	667
	1	2868.06	17301.82	2157	4717	3236.25	1800.63	3575	1207
	2	1863.13	11310.46	1375	3523	2047.23	1092.72	2443	823
	3	376.93	1955.87	281	553	379.67	386.88	587	275
	4	1651.17	4804.95	1203	2241	1138.05	948.92	1771	729
	5	381.03	1900.61	271	475	335.83	276.72	485	189
	6	350.15	1289.96	243	305	246.31	207.73	389	139
	7	106.24	499.16	109	109	80.74	114.22	141	65
	8	134.02	358.58	99	81	116.56	56.94	133	39
	9	494.13	2500.91	355	637	430.41	336.99	529	241
	Total:	1/10	0/10			1/10	8/10

Table 26 :

 26 CPU time (s) and nodes to solve the w100.d010 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	0	27.70	7.09	13	3	5.28	12.13	1	1
	1	57.66	17.95	35	3	43.34	21.59	23	3
	2	131.45 296.92	77	137	595.77 152.09	491	47
	3	43.14	13.46	23	5	35.67	13.81	33	5
	4	2.12	2.12	1	1	1.85	2.47	1	1
	5	18.31	16.32	9	7	24.89	9.76	17	3
	6	2.35	2.35	1	1	2.22	2.64	1	1
	7	50.86	31.53	27	13	24.06	26.43	23	11
	8	1.37	1.37	1	1	1.71	2.08	1	1
	9	12.45	7.85	5	3	2.30	2.69	1	1
	Total:	2/10	3/10			5/10	1/10

Table 27 :

 27 CPU time (s) and nodes to solve the w100.d050 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	0	945.25 2164.97	737	1217	1335.05	1061.46	969	459
	1	211.42	515.41	191	263	290.43	311.70	229	109
	2	175.46	357.04	143	177	186.81	202.08	175	63
	3	612.54 2087.57	489	1159	825.00	563.16	949	319
	4	1040.41	2704.91	835	1499	964.91	720.19	1125	357
	5	1027.07	2436.14	789	1331	1137.77	802.39	1349	413
	6	220.65	587.59	171	305	188.11	189.44	223	93
	7	137.39	196.16	105	93	84.50	86.61	93	39
	8	1152.49	2879.08	875	1637	1069.53	775.14	1255	465
	9	109.99	135.18	91	65	81.02	83.06	87	33
	Total:	3/10	0/10			3/10	4/10

Table 28 :

 28 CPU time (s) and nodes to solve the w100.d090 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2)	(R3)	Our nodes (R2) (R3)
	0	539.51	1163.23	411	629	1038.44	770.67	653	271
	1	2945.89 14097.89	2525	8485	5897.16	4675.78	4413	1821
	2	1364.49	4619.19	1087	2711	1167.26 1519.29	1179	765
	3	1487.39	6598.30	1199	3939	1637.23	1630.03	1789	897
	4	717.48	1943.66	589	1065	800.85	608.70	845	293
	5	79.69	38.67	57	17	28.48	26.50	23	9
	6	89.86	190.11	103	87	131.40	138.67	131	57
	7	377.18	949.78	271	509	374.73	299.03	391	157
	8	205.76	439.04	213	209	187.90	225.81	209	105
	9	596.26	1637.01	459	903	571.00	473.28	637	267
	Total:	4/10	0/10			2/10	4/10

Table 29 :

 29 CPU time (s) and nodes to solve the pw100.d010 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	0	32.32	28.47	15	11	94.80	27.52	63	13
	1	144.42 112.86	107	47	82.52	49.65	83	15
	2	53.49	57.39	33	25	174.22	49.84	123	21
	3	95.28 120.84	61	55	109.82	35.61	99	15
	4	15.11	8.02	7	3	2.64	3.10	1	1
	5	20.20	8.33	9	3	2.50	3.20	1	1
	6	56.65	51.35	29	23	37.24	23.02	35	11
	7	89.73	65.23	57	27	60.77	42.09	57	15
	8	1.64	1.64	1	1	2.06	2.15	1	1
	9	228.70 244.33	157	111	172.00	113.86	191	43
	Total:	1/10	1/10			2/10	7/10

Table 30 :

 30 CPU time (s) and nodes to solve the pw100.d050 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2)	(R3)	Our nodes (R2) (R3)
	0	1798.61 7433.26	1447	4679	3732.88	2506.69	2949	1161
	1	529.82 1787.33	413	1033	692.52	628.58	745	357
	2	831.04 2315.65	697	1381	1136.32	710.04	1095	373
	3	145.66	371.98	119	197	164.25	157.69	185	85
	4	911.92 2763.95	765	1629	1188.62	970.65	1305	549
	5	232.10	374.56	169	201	232.79	218.71	259	117
	6	1726.73	5507.96	1401	3327	1751.65	1359.81	2031	737
	7	454.48	913.14	333	485	384.06	300.39	429	153
	8	122.67	196.11	121	91	103.26	81.34	113	37
	9	406.52 1093.26	315	611	390.92	397.39	445	219
	Total:	4/10	0/10			1/10	5/10

Table 31 :

 31 CPU time (s) and nodes to solve the pw100.d090 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2)	(R3)	Our nodes (R2) (R3)
	0	753.88 4508.30	559	1175	1122.32	677.50	979	321
	1	971.02 6800.50	723	1827	1113.34	883.57	1237	485
	2	362.03 1957.67	267	501	354.55	334.49	405	199
	3	209.02 1261.49	155	301	201.61	255.36	215	131
	4	544.85 2324.07	401	575	396.38	363.92	467	195
	5	647.55 3797.79	533	919	605.58	658.84	733	337
	6	577.00 2826.12	391	717	448.18	405.48	541	261
	7	1297.98	8382.66	981	2253	1112.62	1061.18	1375	621
	8	381.88 2154.22	303	515	347.05	313.59	393	175
	9	323.01 1536.21	247	369	247.90	288.67	289	145
	Total:	0/10	0/10			3/10	7/10

Table 32 :

 32 CPU time (s) and nodes to solve the ising100 problems.

	Prob.	Biq Mac time (R2) (R3)	Biq Mac nodes (R2) (R3)	Our time (R2) (R3)	Our nodes (R2) (R3)
	251005555	165.65 19.86	113	7	2.86	3.15	1	1
	251006666	141.74	9.41	65	3	2.62	2.96	1	1
	251007777	145.33 10.01	75	3	2.96	2.76	1	1
	301005555	142.05 18.78	99	7	3.87	4.51	1	1
	301006666	69.03	7.74	115	3	4.36	2.98	1	1
	301007777	33.13 10.76	11	3	3.51	3.14	1	1
	Total:	0/6	0/6			3/6	3/6

Acknowledgments We are grateful to Angelika Wiegele for the discussions we have had and for providing us with a copy of the Biq Mac solver for our numerical comparison.


Olivier Collet

Master equation of proteins in interaction with implicit or explicit solvent

Keywords: numbers: 87.15.hm, 87.15.kr, 82.20.Yn

Theoretical studies of protein folding on lattice models relie on the assumption that water close to amino-acids is always in thermal equilibrium all along the folding pathway. Within this framework, it has always been considered that out-of-equilibrium properties, such as folding time, could be evaluated equivalently from an averaging over a collection of trajectories of the protein with water described either explicitly or through a mean-field approach. To critically assess this hypothesis, we built a two-dimensional lattice model of a protein in interaction with water molecules that can adopt a wide range of conformations. This microscopic description of the solvent has been used further to derive an effective model by averaging over all the degrees of freedom of the solvent. At thermal equilbrium, the two descriptions are rigourously equivalent, predicting the same folded conformation of the protein. The model allows exact calculations of some relaxation properties using the master equations associated to both solvent descriptions. The kinetic patterns associated to the folding pathways are remarkably different. In this work we demonstrate, that an effective description of the solvent can not described properly the folding pathway of a protein. The microscopic solvent model, that describes correctly the microscopic routes, appears to be the only candidate to study folding kinetics.

Proteins folding is a hot topic of the biophysics field and the question of the mechanisms which governs its kinetics is still in debate. The two ingredients guiding a protein towards its native structure are the distribution of the intrachain interactions and the solvent effect [1]. Thus, the solvation of hydrophobic compounds at thermal equilibrium has been widely studying to extract the key role of the water in protein folding [2][3][4][5][6][7].

However, in most of the works using lattice models for protein, the effect of the solvent is usually taken into account by a temperature independent, structure-less, parameter which simply increases the strength of the some intrachain contact [8,9]. The lattice model using such couplings also provided numerous kinetics works using Monte Carlo simulations [START_REF] Gutin | Proc. Natl. Acad. Sci[END_REF][START_REF] Gutin | [END_REF][12][13] or evolutions of the master equation [14][START_REF] Zwanzig | Proc. Natl. Acad. Sci[END_REF][START_REF] Cieplak | [END_REF][17][18][19][20]. In these works, as the potential associated to each protein structure results from an average over the degrees of freedom, the kinetics is guided by transition rates between chain conformations in interaction with an effective solvent whose mean energy does not depend on the temperature. However, numerous experimental results showing the importance of the relaxation of the first shell solvent on the folding kinetics illuminate the importance of the degrees of freedom of the first shell solvent for the fast kinetics of folding [21][22][23][24] In a few recent works, the contribution of the solvent effect on the configurational Hamiltonian becomes temperature dependent because they result from an average over the water configurations. In these models, the solvent around the proteins is modeled by its energy spectra which takes into account of the formation or breakage of the water hydrogen bonds. Such an approach gave an explanation of cold denaturation [25,26].

Recently, we calculated the kinetics of the folding of a protein model in interaction with implicit water model [27] where the role of the hydrogen bonds of the first shell solvent was taken into account.

Here, we compare the results obtained for the kinetics of folding of protein where each chain structure is in interaction with a highly degenerated microscopic solvent in one hand and with the equivalent effective solvent in an other hand. A micro-state of the system "proteinsolvent" is given by the conformation of the protein and the position of the atoms of water (the solvent configuration). We calculated the evolution in time of the probability of occurrence of each protein-solvent micro-state toward the native structure of the protein, using a master equation approach and starting far from equilibrium. The waiting times to observe the native structure with a probability p, noted t mic;p are calculated as functions of the temperature for p varying from 0.10 to 0.30 by steps of 0.05.

On an other hand, the effective solvent model is introduced by integrating out the water degrees of freedom of the same solvent model and by computing the free energy of hydration of each protein conformation. Simulations start from the equivalent initial condition to that stated for the microscopic model calculation. The evolution in time of the probabilities of the protein conformations is also computed using a master equation and the waiting times, t eff;p , to observe the native structure with a probability p are also calculated under the same conditions.

The curve of t mic;p (T ) and t eff;p (T ), shown in fig. 1, present clearly different shapes and orders of magnitude.

Model. The protein is modeled by a self avoiding walk chain whose the twelve beads are positioned on the nods of a two-dimensional lattice. The number of intrachain contacts of the chain conformation m is

C m = i≥j+3 ∆ (m) ij
where ∆ (m) ij = 1, if the monomers i and j are first neighbors on the lattice, and 0 otherwise. The accessible surface area to the solvent is

A m = 2N + 2 -2C m [26].
The bulk contribution is taken as a mean effect which increases with respect of the number of intrachain contacts and the first shell contribution increases with the accessible surface area to the solvent. One configuration (β = 0) of the first shell is well ordered and the other are disorganized.

The Hamiltonian of the micro-state where the protein is in conformation m, the first shell in configuration β and the bulk in structure α is:

H mic mαβ = i>j B ij ∆ (m) ij + 2C m ε bk + σ(β)A m ε sh for 1 ≤ α ≤ g 2Cm bk and σ(0) = 0 σ(β) = 1 if 1 ≤ β ≤ g Am sh
The intrachain couplings B ij between monomers i and j are drawn at random from a Gaussian distribution centered on B 0 = -2 with standard deviation equals 1 [8].

The values of the solvent parameters are ranked as follow ε bk = 0.4, ε sh = 0.8, g bk = 3.3 and g sh = 3.5 [27]. The Hamiltonian of each protein structure m takes two values following that of β. The ground state, noted (-), is associated to a value of σ = 0 (for β = 0) and the exited state, (+), to σ = 1 (for β > 0). The Hamiltonian and the degeneracy of the macro-states (mσ) are :

H mac mσ = i>j B ij ∆ (m) ij + 2C m ε bk + σA m ε sh g mσ = g 2Cm bk α=1 g Am sh β=0 δ(σ -σ(β)) = g σAm sh g 2Cm bk
where δ(x) = 1 if x = 0 and 0 otherwise. Then, the partition function may be written as :

Z(T ) = m g 2Cm bk α=1 g Am sh β=0 exp - H mic mαβ T (micro -states) = m 1 σ=0 g mσ exp - H mac mσ T (macro -states) = m exp - H eff m (T ) T (effective)
The temperature dependent effective Hamiltonian of conformation m, where the solvent degrees of freedom have been integrating out, is given by :

H eff m = -T ln g 2Cm bk α=1 g Am sh β=0
exp -H mic mαβ T that is to say,

H eff m = i>j B ij ∆ (m) ij + 2C m ε bk -T ln g 2Cm bk 1 + g Am sh exp - A m ε sh T
The effective Hamiltonian has been written to satisfy to the partition function which determines the properties of the system at equilibrium but out-of-equilibrium, this function have no more a clear physical meaning.

The probability of occurrence of the conformation m at equilibrium is : P eq m = exp -H eff m /T /Z(T ). The native conformation is the structure of the largest probability determined by a full enumeration of the conformational space.

Dynamics of the microscopic solvent model. The out-of-equilibrium probability of occurrence of the micro-state (mαβ) at time t is denoted P mic mαβ (t). It evolves following the master equation [27]:

dP mic mαβ (t) dt = m ′ α ′ β ′ X mαβ,m ′ α ′ β ′ P mic m ′ α ′ β ′ (t)
where

X mαβ,m ′ α ′ β ′ = V (0) mm ′ τ mic m,m ′ a T (H mic mαβ , H mic m ′ α ′ β ′ )
is the microscopic transition rate between the configurations (m ′ α ′ β ′ to (mαβ and the diagonal terms are:

X mαβ,mαβ = - m ′ α ′ β ′ X m ′ α ′ β ′ ,mαβ with V (0)
mm ′ = 1 and τ mic m,m ′ = τ c if the chain structures m and m ′ are connected by a one-monomer move and V (0) mm = 1 and τ mic mm = τ s for the solvent modifications which keep the chain structure unchanged. The acceptance function, a T (x, x ′ ) = [1 + exp((x -x ′ )/T )] -1 , satisfies to the microscopic detailed balance condition. These probabilities evolves using an Euler algorithm [START_REF] Cieplak | [END_REF] :

P mic mαβ (t+δt) = P mic mαβ (t)+δt m ′ α ′ β ′ X mαβ,m ′ α ′ β ′ P mic m ′ α ′ β ′ (t)
The out-of-equilibrium probability of occurrence of the macro-state (mσ), at time t is :

P mac mσ (t) = αβ P mic mαβ (t) δ(σ -σ(β))
Then, at time t + δt, it becomes:

P mac mσ (t + δt) = αβ P mic mαβ (t + δt) δ(σ -σ(β)) = αβ P mic mαβ (t)δ(σ -σ(β)) + α,β δt m ′ α ′ β ′ X mαβ,m ′ α ′ β ′ P mic m ′ α ′ β ′ (t) δ(σ -σ(β))
Using the following equality

σ ′ α ′ β ′ δ(σ ′ - σ(β ′ ))f β ′ = α ′ β ′ f β ′ , it comes: P mac mσ (t + δt) = αβ P mic mαβ (t)δ(σ -σ(β)) +δt m ′ σ ′ α ′ β ′ α,β δ(σ ′ -σ(β ′ ))δ(σ -σ(β)) V (0) mm ′ τ mic m,m ′ a T (H mic mαβ , H mic m ′ α ′ β ′ )P mic m ′ α ′ β ′ (t) = P mac mσ (t) + δt m ′ σ ′ Y mσ,m ′ σ ′ P mac m ′ σ ′ with Y mσ,m ′ σ ′ = g mσ V (0) mm ′ τ mic m,m ′ a T (H mac mσ , H mac m ′ σ ′ ) (1)
The increment of time is chosen as δt = 1/ max mσ {Y mσ,mσ } ≪ 1 in order to maintain the sum of the probabilities equals to 1. The form of the transition rates Y mσ,m ′ σ ′ implies that the probability distribution converges to P mac mσ (t) → g mσ exp(-H mac mσ /T )/Z(T ). The simulations starts with the initial condition : P mac mσ (0) = g mσ / mσ g mσ which set the same initial weight to any micro-state.

The time needed to observe the native structure, with a probability p is noted t mic;p or in other words :

1 σ=0 P mac Nat,σ (t mic;p ) = p.
Dynamics of the effective solvent model. Consider now the evolution of the probability of occurrence, P eff m (t) of the chain structure m in interaction with an effective solvent starting with the initial effective probability P eff m (0) = σ P mac mσ (0). at time t = 0. The effective probabilities evolve following the master equation :

dP eff m (t) dt = m ′ V mm ′ P eff m ′ (t)
where

V mm ′ = V (m ′ → m)
is the transition rate from conformations m ′ to m. In order to satisfy to the condition of the convergence towards the equilibrium probability distribution, a solution for the rate is :

V mm ′ = V (0) mm ′ τ eff a T (H eff m , H eff m ′ ) (2)
where τ eff is the effective time associated to a chain move. Defining V mm =m ′ =m V m ′ m and using the Euler algorithm the evolution equation reads:

P eff m (t + δt eff ) = P eff m (t) + δt eff m ′ V mm ′ P eff m ′ (t)
with δt eff = 1/ max m {V mm }. Obviously, the probability distribution tends towards P eff m (t) → P eq m . The probability of the native structure reaches p at a time, denoted by t eff;p .

In contrast with the previous approach, the solvent degrees of freedom are integrated first, here, and afterwards the transition rates are calculated using the effective potential. This is the procedure usually applied in lattice model of protein where the attractive term between monomers results indeed of an average of the solvent degrees of freedom.

Discussion. We address now, the question of a possible equivalence between both descriptions after rewriting the effective characteristic time of the transition between two chain structures would be a, time independent function of the parameters associated to the two connected chain conformations. At this purpose, we require to satisfy the following equality for every protein conformation:

P eff m (t) = P mac m-(t) + P mac m+ (t) ⇒ dP eff m (t) dt = dP mac m-(t) dt + dP mac m+ (t) dt ⇒ m ′ V (0) mm ′ τ eff a T (H eff m , H eff m ′ )P eff m ′ (t) = m ′ σ ′ σ g mσ V (0) mm ′ τ mic m,m ′ a T (H mac mσ , H mac m ′ σ ′ )P mac m ′ σ ′ (t)
As only chain moves may be considered, the above equation leads to :

τ eff = τ c a T (H eff m , H eff m ′ )P eff m ′ (t) σ ′ σ g mσ a T (H mac mσ , H mac m ′ σ ′ )P mac m ′ σ ′ (t)
where it appears clearly that τ eff depends on the time in this equation.

This definitively proves that the kinetics of folding of the proteins can not be understood using protein-solvent models where the degrees of freedom of water have been integrating out in a conformational free energy of solvation.

However, we mention that some relations may be found for extreme temperatures. At very low temperature, as H eff m → H mac m--T ln g m-, and as we may assume that the exited states always have a nil non-equilibrium probability, the above equation leads to an effective characteristic time, only depending on the ground states of the two connected chain conformations:

τ mm ′ eff = 1 + exp(∆H - mm ′ /T ) g m-+ g m ′ -exp(∆H - mm ′ /T ) τ c with ∆H - mm ′ = H mac m--H mac m ′ -.
Putting this results into eq.2, it comes V mm ′ = Y m-,m ′ -and assuming that only the ground states of the protein structures are visited, the two kinetics becomes equivalent. In a similar way, at very high temperature (under which the protein is unfolded), as H eff m → H mac m+ -T ln g m+ and as we may assume that only the excited states have a non nil probability to occur, an equivalent relation between both characteristic time (where the -are replaced by some +) may be found which leads to V mm ′ = Y m+,m ′ + . In both extreme temperature cases, it is possible to rewrite effective kinetics equations equivalent to the microscopic ones but that is not feasible at medium temperature.

Conclusion.

Waiting times to observe some proportions of folded proteins have been calculated using a microscopic description of the solvent and the equivalent mean effect on the chain conformation weights. In both cases, the evolution of the system depends on the ratio of the difference of (free) energies, induced by the attempted moves, over the temperature.

In the first simulations, the (huge) configurational space is composed of all the protein and solvent microstates. The result of the acceptance function of move depends on the energy associated to the microscopic configurations. In other words, the protein and the solvent evolves by performing structural changes between microscopic realizations of the system but the calculations converge slowly towards the equilibrium distribution of the protein conformation, as the value of δt is very small.

In the latter, the conformational space is smaller but the folding takes places in an "free energy" landscape. However, for not too small temperature, the solvation entropy contribution to the values of the effective Hamiltonian leads to free energy values smaller than the ground state energy for each protein structures. As a consequence, the folding takes place in a conformational space in which the values of the effective Hamiltonian are not associated to a physical realization. Here, the simulations converge very fast towards equilibrium (as δt eff ≪ δt) but... by following non-physical routes.

As a consequence, an microscopic solvent model is the only good candidate to study the out-of-equilibrium folding of proteins and the effective solvent may only be restricted to the study of equilibrium properties.

FIG. 1 :

 1 FIG.1: Time after which the probability of occurrence of the native state equals some values p as function of the temperature starting from an equiprobability of each protein-solvent configurations. Calculations are performed with the explicit (left) and the implicit (right) solvent models.