Soodeh Habibi
email: s.habibi@bham.ac.uk

Michal Kočvara
email: m.kocvara@bham.ac.uk

Michael Stingl
email: michael.stingl@fau.de

Solving Lasserre relaxations of unconstrained binary quadratic optimization problems by an interior-point method

Keywords: Binary quadratic optimization, Lasserre hierarchy, semidefinite optimization, interior-point methods, preconditioned conjugate gradients, MAXCUT problem. MSC Classification: 90C22, 90C51, 65F08, 74P05

come

Solving Lasserre relaxations of unconstrained binary quadratic optimization problems by an interior-point method

Introduction

Unconstrained binary quadratic optimization problems (UBQP) represent a surprisingly wide class of important optimization problems; see, e.g., the comprehensive overview [START_REF] Kochenberger | The unconstrained binary quadratic programming problem: a survey[END_REF]. The famous MAXCUT problem is a typical representative of this class. It is thus not surprising that they attract a great deal of attention among algorithm and software developers. Some of the most efficient algorithms for finding a global solution of UBQP combine branch-and-bound techniques with semidefinite programming (SDP) relaxations to obtain good lower bounds. These are represented by software like Biq Mac [START_REF] Rendl | A branch and bound algorithm for max-cut based on combining semidefinite and polyhedral relaxations[END_REF], BiqCrunch [START_REF] Krislock | Biqcrunch: A semidefinite branch-andbound method for solving binary quadratic problems[END_REF] and BiqBin [START_REF] Gusmeroli | BiqBin: a parallel branch-and-bound solver for binary quadratic problems with linear constraints[END_REF].

The SDP relaxations are typically based on Shor's relaxation [START_REF] Ben-Tal | Lectures on Modern Convex Optimization: Analysis, Algorithms, and Engineering Applications[END_REF] which is equivalent to the first-order Lasserre relaxation [START_REF] Lasserre | An Introduction to Polynomial and Semi-algebraic Optimization[END_REF]. While higher-order relaxations would deliver much tighter lower bounds (if not exact solutions), the dimensions of the arising SDP problems are considered prohibitively large already for medium-sized UBQPs. For this reason, many authors proposed various techniques to strengthen the first-order relaxations; see, e.g., [START_REF] Gvozdenović | The operator ψ for the chromatic number of a graph[END_REF][START_REF] Ghaddar | A dynamic inequality generation scheme for polynomial programming[END_REF][START_REF] Campos | Partial lasserre relaxation for sparse Max-Cut[END_REF] and the references therein.

Our goal is to show that (at least) the second-order Lasserre relaxations are solvable by a specialized SDP software. We will also demonstrate that the second-order relaxations are, indeed, superior to the first-order ones and, in many cases, already deliver global solutions of the UBQP. This may allow one to avoid the branch-and-cut machinery, whose behaviour is known to be sometimes unpredictable, and replace it entirely by robust algorithms of known complexity.

Our aim is to solve UBQP min x∈R s x Qx subject to x i ∈ B , i = 1, . . . , s

with a symmetric matrix Q ∈ R s×s , where B is either the set {0, 1} or the set {-1, 1}. We do not assume any sparsity in Q, it is a generally dense matrix. In order to find a global optimum, we use Lasserre hierarchy of SDP problems-relaxationsof growing dimension [START_REF] Lasserre | An Introduction to Polynomial and Semi-algebraic Optimization[END_REF]. The SDP relaxations have the form

min y∈R n q y (2)
subject to M (y) := n i=1

y i M i -M 0 0 .
Here M is a so-called moment matrix, a (generally) dense matrix of a very specific form. In particular, if the solution of ( 1) is unique and the order of the relaxation is big enough, then rank M (y * ) ∈ {1, 2}, where y * is a solution of (2) (see Section 3 below). These problems are known to be difficult to solve due to the quickly growing dimension of the problem with the order of the relaxation; see, e.g., [START_REF] Kim | Binary quadratic optimization problems that are difficult to solve by conic relaxations[END_REF]. Since we do not assume any sparsity in Q, we cannot use sparse techniques, such as those proposed in [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF][START_REF] Wang | TSSOS: A moment-SOS hierarchy that exploits term sparsity[END_REF][START_REF] Wang | CS-TSSOS: Correlative and term sparsity for large-scale polynomial optimization[END_REF]. Instead, we propose to use our recently developed software Loraine [START_REF] Habibi | Loraine-an interior-point solver for lowrank semidefinite programming[END_REF]. Loraine uses a primal-dual predictor-corrector interior-point (IP) method together with an iterative algorithm for the solution of the resulting linear systems. The iterative solver is a preconditioned Krylov-type method with a preconditioner utilizing low rank of the solution.

Several authors observed and confirmed by numerical experiments that SDP reformulation of UBQP (the first-order relaxation) can be efficiently solved by an SDP variant of the ADMM method, see [START_REF] Wen | Alternating direction augmented Lagrangian methods for semidefinite programming[END_REF]. We will show that not only this observation can be extended to higher-order relaxations but that inexact ADMM results can be efficiently used as a warm start for the IP algorithm in Loraine. This is due to the choice of the preconditioner used within the iterative solver in Loraine. Based on this observation, we will propose a new, hybrid ADMM-Loraine algorithm that, for a certain class of problems, will be superior to the single algorithms.

The paper is organized as follows. Section 2 introduces the low-rank SDP solver Loraine and the assumptions needed for its efficiency. In Section 3 we discuss the various forms of SDP relaxations of UBQP. Then, in Section 4 we briefly describe the ADMM algorithm for SDP and introduce the new hybrid algorithm, using an inexact ADMM result as a warm start for Loraine. The last Section 5 is devoted to numerical experiments using instances of the MAXCUT problems and randomly generated UBQPs.

Notation

We denote by S m , S m + and S m ++ , respectively, the space of m×m symmetric matrices, positive semidefinite and positive definite matrices. The notation "svec" and "smat" refer to the symmetrized vectorization and its inverse operation, respectively. The symbol • denotes the Frobenius inner product of two matrices, A • B = trace(A B). Finally, the notation e n (or just e) is used for the vector of all ones.

The solver Loraine

Loraine1 is a general-purpose solver for any linear SDP developed by the authors and implemented in MATLAB and Julia. Compared to other general-purpose SDP software, it particularly targets at problems with low-rank solutions. To solve the arising systems of linear equations, it uses the preconditioned conjugate gradient method, as described in detail in [START_REF] Habibi | Loraine-an interior-point solver for lowrank semidefinite programming[END_REF]. The preconditioner, introduced in [START_REF] Zhang | Modified interior-point method for large-and-sparse low-rank semidefinite programs[END_REF], is based on the assumption that the solution matrix has a small number of outlying eigenvalues. The user can choose between the direct and iterative solver, the type of preconditioner and the expected rank of the solution.

Loraine was developed for problems of the type

min y∈R n , S∈S m , s lin ∈R ν b y (3) subject to n i=1 y i A i + S = C Dy + s lin = d S 0, s lin ≥ 0 with the Lagrangian dual max X∈S m , x lin ∈R ν C • X + d x lin (4) subject to A i • X + (D x lin ) i = b i , i = 1, . . . , n X 0, x lin ≥ 0 .
In the following, we call (3) a problem of primal form and (4) a problem of dual form.

Loraine is efficient under the following assumptions:

Assumption 1 Problems (3),(4) are strictly feasible, i.e., there exist

X ∈ S m ++ , x lin ∈ R ν ++ , y ∈ R n , S ∈ S m ++ , s lin ∈ R ν ++ , such that A i • X + (D x lin ) i = b i , n i=1 y i A i + S = C and Dy + s lin = d (Slater's condition).
Assumption 2 Define the matrix A = [svec A 1 , . . . , svec An]. We assume that any matrixvector products with A and A may each be computed in O(n) flops and memory.

Assumption 3

The inverse (D D) -1 exists and (D D) -1 together with the matrix-vector product with (D D) -1 may each be computed in O(n) flops and memory.

Assumption 4

The dimension of X is much smaller than the number of constraints in (4), i.e., m n.

Assumption 5 Let X * be the solution of (4). We assume that X * has k outlying eigenvalues, i.e., that

(0 ≤) λ 1 (X * ) ≤ • • • ≤ λ(X * ) m-k λ(X * ) m-k+1 ≤ • • • ≤ λ(X * )m ,
where k is very small, typically smaller than 10 and, often, equal to 1. This includes the particular case when the rank of X * is very small.

The last Assumption 5 is not satisfied by problem (2) so, in the next section, we will:

(i) re-write SDP relaxation (2) in the dual form (4) by introducing auxiliary variables and additional linear equality constraints; (ii) treat the new linear equality constraints by 1 -penalty approach, in order to replace equalities by inequalities, as required by the interior-point algorithm; (iii) show that the matrix associated with the new linear inequality constraints is block diagonal with blocks of very small size and thus satisfies Assumption 2; (iv) show that Assumption 4 is naturally satisfied for the re-written problem.

3 Forms of SDP relaxations

Lasserre hierarchy of moment problems

Let B = {-1, 1}. By introducing a matrix variable X = xx , problem (1) can be equivalently written as the following SDP problem with a rank constraint:

min X∈S s Q • X (5) subject to X i,i = 1, i = 1, . . . , s X 0 rank X = 1 .
Remark 1 For B = {0, 1}, we first use the substitution x = 2x -1 to get a problem with objective function x Qx + 2e Qx and constraints xi ∈ {-1, 1}. Then we set

Q = 0 e Q Qe Q ∈ S s+1 , introduce a variable X ∈ S s+1 , X = 1 x 1 x
and solve problem (5) with Q, X replaced by Q, X.

In view of the above remark, in the rest of the paper we will only consider the case B = {-1, 1}. Now, to each unique component of the matrix X = xx we assign a variable y k . Considering the constraints X i,i = 1 and the symmetry of X, we thus have y ∈ R r with r = s(s -1)/2. Replacing the variable X by y, the constraint X 0 can then be written as r i=1 M i y i + I 0, where M i ∈ S s are suitable matrices; for details, see [START_REF] Lasserre | An Introduction to Polynomial and Semi-algebraic Optimization[END_REF].

By ignoring the rank constraint, we can now define the first-order Lasserre relaxation of (5) as the following moment problem:

min y∈R r q y (6) subject to r i=1 M i y i + I 0 ,
where q = svec Q.

Example 1 For instance, for s = 3 we use the vector (monomial basis) x 1 , x 2 , x 3 ; then r = 3 and

M 1 =   0 1 0 1 0 0 0 0 0   , M 2 =   0 0 1 0 0 0 1 0 0   , M 3 =   0 0 0 0 0 1 0 1 0  
and the moment matrix r i=1 M i y i + I is associated with the original variables as

(M (y) =) r i=1 M i y i + I =   1 y 1 y 2 y 1 1 y 3 y 2 y 3 1   ←→   x 1 x 2 x 3     x 1 x 2 x 3   =   x 2 1 x 1 x 2 x 1 x 3 x 1 x 2 x 2 2 x 2 x 3 x 1 x 3 x 2 x 3 x 2 3   .
The Lasserre hierarchy consists in adding higher-order monomials in the monomial basis and, therefore, the higher-order relaxation problems will all be of the form (6) only of growing dimension.

Example 2

The second-order relaxation in Example 1 will lead to a (6 × 6) moment matrix associated with the dyadic product of the vector of monomials of order up to two

(x 1 , x 2 , x 3 , x 2 1 , x 2 2 , x 2 3 , x 1 x 2 , x 1 x 3 , x 2 x 3 ) with itself.
As some of the terms in the higher-order moment matrix are repeated and some can be simplified using the constraint x 2 i = 1 (e.g., x 1 x 2 2 is reduced to x 1 ) and thus even more terms are repeated, the number of variables in the higher-order relaxations is smaller than the number of the elements of the upper triangle of the matrix (see also Example 3 below).

Denote the order of the Lasserre relaxation by ω. The dimensions n, m grow quickly with the order of the relaxation. In particular, we have

s 4 24 n ≤ 2 s -1 and m = ω i=1 s i . (7)
(The authors are not aware of any exact formula for n or the lower bound on n; the lower bound in ( 7) is our estimate based on numerical experiments. Obviously, the higher ω the bigger n until it reaches the "saturation point" 2 s -1.)

The dimension of the problem makes it very challenging for standard SDP solvers based on second-order methods, such as interior-point methods; see [START_REF] Kim | Binary quadratic optimization problems that are difficult to solve by conic relaxations[END_REF] and the numerical examples in the last section of this paper.

It was shown by Fawzi et al. [START_REF] Fawzi | Sparse sums of squares on finite abelian groups and improved semidefinite lifts[END_REF] and Sakaue et al. [START_REF] Sakaue | Exact SDP relaxations with truncated moment matrix for binary polynomial optimization problems[END_REF], and numerically confirmed in [START_REF] Kim | Binary quadratic optimization problems that are difficult to solve by conic relaxations[END_REF], that for this type of problems, the sequence of these approximations is finite and the upper bound on the order of the relaxation to obtain exact solution of (1) is n/2 . Laurent [START_REF] Laurent | Lower bound for the number of iterations in semidefinite hierarchies for the cut polytope[END_REF] showed that this is also a lower bound for MAXCUT problems with unweighted complete graphs.

To determine whether a relaxation of order ω is exact, i.e., whether its solution is the solution of (1), we use the rank of the moment matrix; see [START_REF] Lasserre | An Introduction to Polynomial and Semi-algebraic Optimization[END_REF]Thm. 6.19]. In particular, assuming that (1) has p global minimizers, the rank of the moment matrix will be less than or equal to p. Therefore, if we assume that (1) has a unique solution, then the moment matrix of the exact relaxation will be of rank one.

Remark 2 Notice that, when B = {-1, 1}, every global minimizer x will have a "symmetric" counterpart -x, as there are only quadratic terms in the problem. In this case, by unique minimizer, we understand "unique up to the multiple by -1" and the rank of the exact relaxation will be two.

Above, we spoke about an "exact relaxation" and a "rank of the moment matrix". When solving the problem by an interior point method, we can only speak about a numerical rank defined below.

Definition 2 Consider A ∈ S m
+ with eigenvalues λ i (A) and with the largest eigenvalue λmax(A). Let ε > 0. The numerical rank of A is defined as

rε(A) = |I|, I = i ∈ {1, . . . , m} | λ i (A) λmax(A) ≥ ε .
Definition 3 For problems with a unique solution and B = {-1, 1}, the relaxation of order ω is called exact if the numerical rank of the corresponding moment matrix is less than or equal to two.

Re-writing the problem

As mentioned in the Introduction, problem ( 6) is in the "wrong" form for our solver Loraine. While the solution M * := M (y * ) of ( 6) has low rank (in particular, rank M * = 2 when the solution of ( 1) is unique), Loraine expects the solution of an SDP problem in the dual form (4) to have low rank. Our next goal is thus to rewrite (6) in the dual form; in other words, to rewrite the dual to [START_REF] Lasserre | An Introduction to Polynomial and Semi-algebraic Optimization[END_REF] in the primal form (3). The dual to ( 6) is the problem

max Z∈S m -I • Z subject to M i • Z = q i , i = 1, . . . , n Z 0 ,
and, by vectorizing the matrices as z = svec Z and M = (svec M 1 , . . . , svec M n ) , M ∈ R n× n , it can be further written as

min z∈R n (svec I) z (8)
subject to smat(z) 0 Mz = q

with n = m(m + 1)/2. Now, from the fact that the problem ( 8) is the dual problem to (6), we have the following lemma:

Lemma 1 The vector y in (6) is the Lagrangian multiplier to Mz = q in (8).

Problem ( 8) is now almost in the right form for our solver: it is in the primal form and the dual solution is a low-rank matrix. Notice that, despite Q being possibly a dense matrix, the data in [START_REF] Ghaddar | A dynamic inequality generation scheme for polynomial programming[END_REF] are always very sparse. One obstacle remains-the linear equality constraints in [START_REF] Ghaddar | A dynamic inequality generation scheme for polynomial programming[END_REF]. We are addressing it in the next section.

Equality constraints by 1 -penalty

Interior-point methods have been designed for optimization problems with inequality constraints and cannot directly handle equality constraints. There are various ways how to overcome this, starting with writing equality as two inequalities; see [START_REF] Anjos | On handling free variables in interior-point methods for conic linear optimization[END_REF] for an overview focused on SDP interior-point algorithms. In our algorithm, we treat the equality constraints by the following 1 -penalty approach.

Introducing 1 -penalty for the linear equality constraints in [START_REF] Ghaddar | A dynamic inequality generation scheme for polynomial programming[END_REF], we obtain the following problem

min z∈R n (svec I) z + µ Mz -q 1 (9)
subject to smat(z) 0 , with a penalty parameter µ > 0. Recall that this penalty is exact, in the sense that there exists µ 0 > 0 such that the solution of ( 9) with µ ≥ µ 0 is equivalent to the solution of (8) (see, e.g., [START_REF] Di Pillo | Exact penalty functions in constrained optimization[END_REF]). We now introduce two new variables, r ∈ R n , s ∈ R n , satisfying

Mz -q = r -s, r ≥ 0, s ≥ 0 .
After substitution, (9) reads as

min z∈R n , r∈R n , s∈R n (svec I) z + µ n i=1 (r i + s i ) (10)
subject to smat(z) 0 Mz -q = r -s r ≥ 0, s ≥ 0 .

Using the identity r = Mz -q + s to eliminate variable r, we arrive at our final problem

min z∈R n , s∈R n (svec I) z + µ n i=1 ((Mz -q) i + 2s i ) (11)
subject to smat(z) 0

Mz -q + s ≥ 0 s ≥ 0 .
Problem [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF] is now in the primal form [START_REF] Krislock | Biqcrunch: A semidefinite branch-andbound method for solving binary quadratic problems[END_REF]. By Lemma 1, we know that y in ( 6) is the Lagrangian multiplier to the equality constraint in [START_REF] Ghaddar | A dynamic inequality generation scheme for polynomial programming[END_REF]. The question is how to obtain y from [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF], the 1 -penalty reformulation of (8).

Lemma 2 A solution y of (6) can be obtained from the solution of (11) as y i = µ + λ i , i = 1, . . . , n, in which µ is the penalty parameter at the solution and λ is the Lagrangian multiplier to the constraint Mz -q + s ≥ 0.

Proof We start with the optimality conditions of the problem [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF]. Using suitable matrices E i ∈ S m , we can write smat(z) = n i=1 E i z i . The Lagrangian of ( 11) is then

L 1 (z, s, Γ, λ, η) = -f (z) + µ n i=1 ((Mz -q) i + 2s i ) -Γ • n i=1 E i z i - n i=1 λ i   ñ j=1 M ij z j -q i + s i   - n i=1 η i s i ,
where Γ 0 and λ, η ≥ 0 are the corresponding Lagrangian multipliers and f (z) = (svec I) z. By using the KKT conditions at the solution we have

∇z L 1 (z * , s * , Γ * , λ * , η * ) = 0 : -(∇f (z * )) j -µ n i=1 M ij - n i=1 M ij λ * i -Γ * • E j = 0, j = 1, . . . , n
and so

(∇f (z * )) j = - n i=1 M ij (µ + λ * i ) -Γ * • E j , j = 1, . . . , n . (12)
Now, the Lagrangian of the problem (8) is

L 2 (z, Γ, y) = -f (z) -Γ • n i=1 E i z i - n i=1 y i   ñ j=1 M ij z j -q i   ,
where, as above, Γ 0 and y ≥ 0 are the corresponding Lagrangian multipliers. From the KKT conditions, we have ∇z L 2 (z * , γ * , y * ) = 0 :

-(∇f (z * )) j -Γ * • E j - n i=1 M ij y * i = 0, j = 1, . . . , n so (∇f (z * )) j = - n i=1 M ij y * i -Γ * • E j , j = 1, . . . , n . (13)
Thus, from ( 12) and ( 13), y can be written as

y i = µ + λ i , i = 1, . . . , n (14)
and we are done.

The question remains how to initialize and update the penalty parameter µ in order to get an exact solution of problem (8) using our interior-point solver. There are essentially two options.

(A) For the current value of µ, fully solve the penalized problem [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF] by Loraine, check the residuum of the equality constraint in [START_REF] Ghaddar | A dynamic inequality generation scheme for polynomial programming[END_REF] for the optimal solution of (11) and, if bigger than a prescribed threshold, increase µ and repeat this process. (B) Incorporate penalty update within the interior-point solver, following, e.g., [START_REF] Gould | An interior-point 1-penalty method for nonlinear optimization[END_REF] and [START_REF] Curtis | A penalty-interior-point algorithm for nonlinear constrained optimization[END_REF]. That is, run the IP algorithm with an initial µ and, as soon as the residuum of the equality constraint [START_REF] Ghaddar | A dynamic inequality generation scheme for polynomial programming[END_REF] gets bigger than a value dependent on the IP barrier parameter, increase µ and continue in IP iterations. Option (A) is simple and robust but may be expensive: we may need to solve several SDP problems to full optimality. Option (B) is more attractive from the efficiency point of view but is also more sensitive to the choice of parameters. We implemented (B) but had mixed experience with it so we returned to (A). The good news is that a class of BQP problems usually requires a single choice of µ in order to guarantee convergence to an optimal point with the norm of the residuum of the equality constraint smaller than a prescribed value; see the numerical experiments in the last section. Needless to say that one cannot start with a too big value of µ due to the potential ill-conditioning of the resulting SDP problem and resulting difficulties of the IP algorithm to solve the problem at all.

Loraine and problem (11)

Recall the definition of the matrix M = (svec M 1 , . . . , svec M n ) T with matrices M i introduced at the beginning of this section. We have the following result.

Lemma 3 There exists a permutation matrix P ∈ R n × n such that P M MP is a block diagonal matrix with small full blocks. In particular, M M is a sparse chordal matrix.

Proof Every matrix M i localizes the corresponding variables in the moment matrix n i=1 M i y i +I in [START_REF] Lasserre | An Introduction to Polynomial and Semi-algebraic Optimization[END_REF], whereas the elements of the moment matrix are either ones or the variables y i ; see also Example 1. As every M i is associated with exactly one variable y i , the nonzeros in this matrix are unique to the matrix. Therefore every column of matrix M will either contain all zeros or a single nonzero number one at the i-th position, i = 1, . . . , n. We can now re-order the columns of M such that the first columns will only contain zeros, the next columns element 1 at the first position, the next ones element 1 at the second position, etc. This re-ordering will define the permutation matrix P . The claims follow.

Example 3 Consider a problem with s = 3 and relaxation order ω = 2. Utilizing the constraints x 2 i = 1, we only consider variables associated with the monomial basis b(x) = (x 1 , x 2 , x 3 , x 1 x 2 , x 2 x 3 , x 1 x 3 ):

X = b(x)b (x) =          x 2 1 x 1 x 2 x 1 x 3 x 2 1 x 2 x 1 x 2 x 3 x 2 1 x 3 x 1 x 2 x 2 2 x 2 x 3 x 1 x 2 2 x 2 2 x 3 x 1 x 2 x 3 x 1 x 3 x 2 x 3 x 2 3 x 1 x 2 x 3 x 2 x 2 3 x 1 x 2 3 x 2 1 x 2 x 1 x 2 2 x 1 x 2 x 3 x 2 1 x 2 2 x 1 x 2 2 x 3 x 2 1 x 2 x 3 x 1 x 2 x 3 x 2 2 x 3 x 2 x 2 3 x 1 x 2 2 x 3 x 2 2 x 2 3 x 1 x 2 x 2 3 x 2 1 x 3 x 1 x 2 x 3 x 1 x 2 3 x 2 1 x 2 x 3 x 1 x 2 x 2 3 x 2 1 x 2 3          =         1 x 1 x 2 x 1 x 3 x 2 x 1 x 2 x 3 x 3 x 1 x 2 1 x 2 x 3 x 1 x 3 x 1 x 2 x 3 x 1 x 3 x 2 x 3 1 x 1 x 2 x 3 x 2 x 1 x 2 x 1 x 1 x 2 x 3 1 x 1 x 3 x 2 x 3 x 1 x 2 x 3 x 3 x 2 x 1 x 3 1 x 1 x 2 x 3 x 1 x 2 x 3 x 1 x 2 x 3 x 1 x 2 1         i.e., variables y ∈ R 7 corresponding to (x 1 , x 2 , x 3 , x 1 x 2 , x 2 x 3 , x 1 x 3 , x 1 x 2 x 3
), the unique elements of X. The matrices M i simplify to

M 1 =        
0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0

        , M 2 =        
0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0

        , . . . , M 7 =         0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0        
and thus

M =          
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 1 0 0 0 0

          and MP =          
0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1

         
with a suitable permutation matrix P . Then P M MP is block diagonal with block sizes 2 and 3.

Remark 3

The size of the blocks in the matrix P M MP is given by the count of the corresponding variables in the upper triangle of the moment matrix. This number grows with the relaxation order.

Theorem 4 Assume that (1) has a unique solution and that Lasserre hierarchy is exact for order ω ≥ 1. Then problem [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF] corresponding to ω satisfies Assumptions 1-5.

Proof A strictly feasible point for (11) can be found, e.g., by choosing z as a unit vector and s positive, arbitrarily large. Similarly, the moment problem ( 6) is trivially strictly feasible by choosing y = 0. Hence Assumption 1 holds. Every matrix E i contains at most two nonzero elements, hence Assumption 2 is trivially satisfied. Assumption 3 holds by Lemma 3: because the matrix M M is chordal and sparse, sparse Cholesky factorization leads to zero fill-in and is thus very efficient. Assumption 4 is satisfied by construction of the problem, as the number of variables n is proportional to the square of the size of the matrix inequality m. Finally, under the assumption of unique solution to (1) and exact relaxation, the dual variable associated with the matrix inequality in [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF] has rank at most two. Hence Assumption 5 is satisfied.

ADMM for SDP relaxations

The Alternating Direction Method of Multipliers (ADMM) is a popular alternative to the interior-point methods for convex optimization problems. Its extension to semidefinite optimization was introduced in [START_REF] Wen | Alternating direction augmented Lagrangian methods for semidefinite programming[END_REF]. It proved to be efficient for certain SDP problems, among others many UBQPs ([15, Sec 4.3]). We thus briefly describe this algorithm and offer a numerical comparison with Loraine. Furthermore, we will introduce a new "hybrid" approach combining ADMM with Loraine.

The ADMM-SDP algorithm

The method consists of minimizing the augmented Lagrangian of the problem

L ρ (y, S, X) = -b y + X • n i=1 y i A i + S -C + 1 2ρ n i=1 y i A i + S -C 2 F
alternately, first with respect to y, then with respect to S, and then updating X by

X (k+1) = X (k) + n i=1 y (k+1) i A i + S (k+1) -C ρ .
While the minimization with respect to y leads to an unconstrained problem and the minimum is obtained by solving the first-order optimality system of (linear) equations, the minimization of L ρ with respect to S can be formulated as a projection of the matrix C -n i=1 y i A i -ρX onto the semidefinite cone S m + . For Z ∈ S m , let A(Z) be a vector with elements A i • Z, i = 1, . . . , n and [Z] + the projection of Z on the semidefinite cone. With a penalty parameter ρ > 0 and a relaxation parameter σ ∈ (0, 1+ √ 5

2 ), the update of the primal-dual point (y, S, X) in the k-th iteration of ADMM is defined in Algorithm 1.

For details, in particular the choice and update of ρ, see [START_REF] Wen | Alternating direction augmented Lagrangian methods for semidefinite programming[END_REF]. In the numerical experiments below, we use our MATLAB implementation2 of the algorithm.

Algorithm 1 ADMM for SDP

Given an initial point (y (0) , S (0) , X (0) ), parameters ρ and σ, and a stopping parameter ε ADMM .

1: for k = 0, 1, 2, . . . do

2:

Update y by y

(k+1) = -(AA ) -1 ρ(A(X (k) ) -b) + A(S (k) -C) 3: Update S by S = C - n i=1 y (k+1) i A i -ρX (k) , S (k+1) = S + 4: Update X by X = 1 ρ S (k+1) -S , X (k+1) = (1 -σ)X (k) + σ X 5: Update ρ 6:
Check convergence 7: end for

ADMM as a warm starter for Loraine

When testing the ADMM algorithm with the MAXCUT problems considered in this paper, we have made the following two observations:

(i) The method is more efficient when solving the reformulated (bigger) problem ( 11) than the smaller one [START_REF] Lasserre | An Introduction to Polynomial and Semi-algebraic Optimization[END_REF]. Moreover, for problem [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF], the convergence is very fast in the first iterations, before it slows down: the algorithm gets quickly a "good" approximation of the solution, both primal and dual. (ii) Due to the nature of the ADMM method (projections on the semidefinite cone), the rank of the approximate solution is always "numerically exact"; in particular, the rank of X (k) is low when we are close enough to the solution. The fact that a low-rank approximation of the solution can be obtained relatively quickly, together with the fact that ADMM is a primal-dual algorithm leads to the following idea: use a low-precision ADMM solution as a warm start for Loraine. Presuming that the ADMM approximation of X is already of the expected low rank (or very close to it), the preconditioner H α will be extremely efficient during the remaining iterations of Loraine. We therefore propose the following hybrid ADMM-Loraine Algorithm 2.

Algorithm 2 ADMM-Loraine 1: Solve (11) by ADMM with stopping tolerance ε ADMM . Save the primal-dual approximate solution (y, S, X) ADMM . 2: Solve ( 11) by Loraine with -initial point (y, S, X) ADMM -initial value of the stopping criterion for the CG method reduced to 10 -6 .

Notice that, in order to get a primal-dual approximation of the solution, we have to solve the same problem formulation by ADMM and Loraine, i.e., formulation [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF] with the 1 -penalty and the same value of the penalty parameter µ, despite the fact that ADMM could directly handle the linear equality constraints.

Numerical experiments

In this section we report the results of our numerical experiments with relaxations of MAXCUT problems and randomly generated UBQPs. The presented algorithms will be compared with the state-of-the-art SDP solver MOSEK [START_REF] Mosek Aps | The MOSEK Optimization Toolbox for MATLAB Manual[END_REF] that is relying on the interior-point method and uses a direct linear solver for the Schur complement equation.

The stopping criterion for our implementations of Loraine and ADMM was based on the DIMACS errors [START_REF] Mittelmann | An independent benchmarking of SDP and SOCP solvers[END_REF] that measure the (normalized) primal and dual feasibility, duality gap and complementary slackness. The algorithms were stopped when all these errors were below 10 -6 . MOSEK was used with default settings.

All problems were solved on an iMac desktop computer with 3.6 GHz 8-Core Intel Core i9 and 64 GB 2667 MHz DDR4 using MATLAB R2022b.

MAXCUT problems

For the first set of numerical experiments, we use MAXCUT problems, the standard test problems for UBQP algorithms. Let Γ be an undirected n-node graph and let the arcs (i, j) be associated with nonnegative weights a ij . The MAXCUT problem is formulated as the following UBQP:

max x 1 4 n i,j=1 a ij (1 -x i x j ) | x 2 i = 1, i = 1, . . . , n (15)
and can be thus solved by the techniques presented above. The solution of the SDP relaxation (2) is of rank two whenever the relaxation is exact and the solution to [START_REF] Wen | Alternating direction augmented Lagrangian methods for semidefinite programming[END_REF] is unique. (Recall that uniqueness of x means uniqueness up to the multiple by -1.) Obviously, a complete unweighted graph (with a ij ∈ {0, 1}) may have many "symmetric" solutions. Therefore we generated two sets of problems of increasing size. Firstly, to avoid the non-uniqueness, we generated undirected, weighted, generally complete graphs with weights randomly distributed between 0 and 12 with 20-50 nodes, using the MATLAB command graph. It turned out that for all these problems the relaxation order two in the Lasserre hierarchy is already high enough to deliver the optimal solution of the MAXCUT problem [START_REF] Wen | Alternating direction augmented Lagrangian methods for semidefinite programming[END_REF]. Secondly, we generated a set of problems of the same sizes but for unweighted graphs. The dimensions of the corresponding SDP problems for relaxation order two are shown in Table 1.

We first present the results for the weighted problems, see Table 2. These problems are relatively "simple": the ADMM method is rather efficient, and MOSEK converges in a very small number of iterations, six to eight. However, the sheer size prevents MOSEK with a direct solver to solve larger problems. We present results for ADMM applied to both problem formulations, the original one (6) and the re-written one of larger dimension and with 1 -penalty [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF]. Perhaps surprisingly, ADMM is more efficient for the latter problem, despite the bigger dimension. Moreover, unlike for ADMM applied to [START_REF] Lasserre | An Introduction to Polynomial and Semi-algebraic Optimization[END_REF], the number of iterations of ADMM applied to [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF] is almost independent of the problem size. The last four columns of Table 2 show Table 1 Problems MAXCUT-<n>, relaxation order ω = 2: dimensions for formulations [START_REF] Lasserre | An Introduction to Polynomial and Semi-algebraic Optimization[END_REF] and [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF] problem [START_REF] Lasserre | An Introduction to Polynomial and Semi-algebraic Optimization[END_REF] problem ( 11 the results for the hybrid ADMM-Loraine algorithm: we present the number of iterations of ADMM plus Loraine ('iter'), the time of ADMM ('timeA'), time of Loraine ('timeL') and the total time of the hybrid method. For these experiments, we have used the stopping tolerance of ADMM ε ADMM = 5 • 10 -3 (problems 20-25) and ε ADMM = 5 • 10 -4 (problems 30-50). To better demonstrate the behaviour of the hybrid method, in Figure 1 we show the output of both codes. We can see that the ADMM method, indeed, gets very quickly a good approximation of objective value, while Loraine uses the warm start very efficiently.

Next we will try to solve the unweighted MAXCUT problems. Recall that these problems often have nonunique solutions, in particular problems with almost dense graphs. Moreover, relaxation order ω = 2 may not be high enough to obtain an exact solution. We thus cannot expect the preconditioner to be efficient in those cases. This is, indeed, demonstrated in Table 3. This table, in addition to columns identical to Table 2, shows also the rank of the solution of the order-2 relaxation. As expected, Loraine is less efficient for problems with higher solution rank, in particular problems MAXCUT-30 and MAXCUT-50. Still, it can solve these problems reliably. Again, the hybrid method is superior for these problems. The convergence behaviour of the algorithms and their estimated complexity are further illustrated in Figure 2. ---------------------------------------------------------------------- ---------------------------------------------------------------------- -----------------------------------------------------------------------

Randomly generated problems

We further conducted numerical experiments with matrices Q randomly generated. We are aware of the fact that problems with random data may not always be representative and may sometimes lead to false conclusions regarding algorithm behaviour. However, we believe that these results still well demonstrate the efficiency of our approach. It has been observed in [START_REF] Kim | Binary quadratic optimization problems that are difficult to solve by conic relaxations[END_REF] that problems with rank-one matrix Q may require relaxation order of up to ω = s/2 to reach the exact solution; a typical example is Q = ee . On the other hand, for problems with matrix Q of rank 3, ω = 2 was always sufficient in experiments performed in ([10, Fig. 4]). We have thus considered two main classes of problems: (A) problems with B = {-1, 1} and with rank Q = 1 generated by the following MATLAB code rng(0); q = randn(s,1); Q = q * q';

(B) problems with B = {-1, 1} and with a full-rank indefinite Q generated by the following MATLAB code rng(0); q = randn(s,1); Q = q * q'; for k=1:s-1 if ceil(k/2) * 2 == k q = randn(s,1); Q = Q -q * q'; else q = randn(s,1); Q = Q + q * q'; end end Apart from random numbers generated with normal distribution, we also performed tests with uniform distribution (function rand) and lognormal distribution (function logncdf); in both cases the results and conclusions were rather similar to the above choice and are thus not reported here.

The exactness of the relaxation was measured by the numerical rank of the dual solution to the matrix inequality in problem (11)-when the rank was equal to 1 or 2 (depending on the set B) the relaxation order was considered sufficient; see Definition 3.

Remark 4 For Q constructed as in (B), relaxation order ω = 2 was sufficient to get an exact solution of (1). This observation, though, cannot be extended to any full-rank matrix Q. For instance, for Q = ee +Diag(d), d ∈ R s , the lowest relaxation order will be ω = s/2 . This is because x 2 i = 1 and thus the diagonal elements of Q will be irrelevant in the optimization process and the resulting problem will be equivalent to that with a rank-one matrix Q = ee .

Full-rank Q and relaxation order ω = 2

Using the MATLAB code from point (B) above, we generated problems of growing dimension s = 10 . . . 50 and solved the corresponding SDP relaxations for order ω = 2. The dimensions of the generated problems are reported in Table 4; the table shows problem sizes for the original SDP relaxation [START_REF] Lasserre | An Introduction to Polynomial and Semi-algebraic Optimization[END_REF] and for the re-written problem [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF]. The computational results are presented in Table 5 and clearly demonstrate the efficiency of Loraine. The computational complexity of Loraine (applied to [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF]) and MOSEK and ADMM (applied to ( 6)) is further illustrated in Figure 3. We do not report on the hybrid ADMM-Loraine algorithm. That is because ADMM applied to problem [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF] (as required by the hybrid method) appears to be much less efficient than for the MAXCUT problems and so the hybrid algorithm is not even competitive to ADMM applied to [START_REF] Lasserre | An Introduction to Polynomial and Semi-algebraic Optimization[END_REF].

Rank-one Q and higher relaxation order

Let us recall from [START_REF] Gvozdenović | The operator ψ for the chromatic number of a graph[END_REF] the dependence of problem sizes on the relaxation order: for the original moment problem (6), we have On the other hand, in the re-written problem [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF] with 1 -penalty term, we have n +n variables with n = m 2 , LMI of size m and 2n linear inequality constraints; see Table 6. This table, in particular, shows that the number of variables grows much more quickly for problem [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF] than for the original problem [START_REF] Lasserre | An Introduction to Polynomial and Semi-algebraic Optimization[END_REF] where it, eventually, reaches the finite limit 2 s -1. Therefore we cannot expect Loraine with iterative solver applied to [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF] to be as efficient for higher-order relaxations as it is for ω = 2. This is clearly demonstrated in Table 7 comparing Loraine with a direct solver applied to problem (6) (this code was slightly faster than MOSEK for these problems, due to direct handling of the rank-one data matrices) with Loraine with the iterative solver applied to problem [START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF].

Conclusions

Our numerical experiments demonstrate the ability of an interior point method with a specially designed solver of the linear system to solve higher-order Lasserre relaxations of UBQP. The approach is particularly efficient for relaxation order two which, for most of the tested problems, was high enough to deliver either the exact solution of the UBQP or a good approximation of it.

We have also introduced a new, hybrid algorithm that uses an approximate solution obtained by ADMM as a warm start for the used interior point method. This algorithm is rather efficient for problems for which ADMM itself is relatively efficient, such as the MAXCUT problems. On the other hand, it may be inefficient once ADMM fails to converge quickly to an approximate solution.

  ================================= A D M M for S D P ================== Number of LMI Constraints: 1 Number of Variables: 1065901 Maximal Constraint Size: 1276 -

Figure 1

 1 Figure1The printout of low-precision ADMM and warm-started Loraine for the weighted MAXCUT-50.

Figure 2

 2 Figure 2 CPU times in log-log scale for weighted (left) and unweighted (right) MAXCUT problems. Loraine (blue), MOSEK (brown), ADMM (green) and hybrid ADMM-Loraine (red).

Figure 3

 3 Figure 3 CPU times in log-log scale for random BQP problems. Loraine (blue), MOSEK (brown) and ADMM (green).

  Definition 1 The Lasserre relaxation of order ω of the problem (1) is given by problem[START_REF] Lasserre | An Introduction to Polynomial and Semi-algebraic Optimization[END_REF] with general dimensions n, m, replacing r, s, respectively; i.e., y ∈ R n and M i ∈ S m , i = 1, . . . , n.

Table 2

 2 Loraine, MOSEK, ADMM and ADMM-Loraine in weighted MAXCUT-<n> problems, relaxation order ω = 2

	MAXCUT	Loraine for (11)	MOSEK (6)	ADMM for (6)	ADMM for (11)	ADMM-Loraine for (11)
	problem	iter	CG it	time	time	iter	time	iter	time	iter	timeA	timeL	time
	20	18	559	3.8	9	3615	10	3042	13	628+4	3.9	0.8	4.7
	25	20	728	11	78	4732	33	2735	28	181+5	2.5	2.9	5.4
	30	21	1032	28	607	6770	99	3537	80	795+5	20	5.8	26
	35	23	2183	96	2911	5255	164	3030	126	863+4	42	11	53
	40	27	2275	186	mem	9611	500	1280	92	914+7	73	132	205
	45	25	2521	335		16901	1400	2639	358	755+4	104	33	137
	50	24	2540	528		19521	2951	3296	745	727+5	162	58	220

Table 3

 3 Loraine, MOSEK, ADMM and ADMM-Loraine in unweighted MAXCUT-<n> problems, relaxation order ω = 2

	MAXCUT		Loraine for (11)	MOSEK (6)	ADMM for (11)	ADMM-Loraine for (11)
	problem	iter	CG it	time	time	iter	time	iter	timeA	timeL	time	rank
	20	17	735	4.3	9	1981	8	738+7	3.4	1.5	4.9	2
	25	18	858	14	82	805	9	328+4	3.5	2.8	6.3	4
	30	23	6482	145	635	2218	53	566+8	13	135	148	57
	35	21	2604	128	3357	5251	229	1222+8	48	15	63	2
	40	22	3538	305	mem	7137	575	1855+7	142	35	177	2
	45	22	3225	536		6794	845	1081+4	139	49	188	4
	50	25	12599	2593		5723	1162	1677+5	361	365	726	13

Table 5

 5 Randomly generated UBQP problems, relaxation order ω = 2: Loraine, MOSEK and ADMM

			Loraine for (11)	MOSEK for (6)	ADMM for (6)
	UBQP size	iter	CG iter	time	iter	time	iter	time
	10	10	256	0.1	6	0.2	1006	0.4
	15	10	725	1.3	5	1.1	1928	2.2
	20	11	423	2.4	6	9.3	2888	9.7
	25	13	294	5.5	7	81	5622	48
	30	13	326	14	7	496	8498	147
	35	15	530	43	memory	11672	396
	40	14	730	106			14362	880
	45	16	896	230			14320	1369
	50	16	1114	431			21701	3251

Table 6

 6 Problem dimensions for formulations[START_REF] Lasserre | An Introduction to Polynomial and Semi-algebraic Optimization[END_REF] and[START_REF] Waki | Algorithm 883: SparsePOP-a sparse semidefinite programming relaxation of polynomial optimization problems[END_REF] as functions of the relaxation order; here s is the number of variables in BQP[START_REF] Kochenberger | The unconstrained binary quadratic programming problem: a survey[END_REF] and ω the relaxation order and s4 24 n ≤ 2 s -1

	problem (6)						problem (11)
	variables	matrix size	variables		matrix size	lin. constraints
	n	ω i=1	s i	ω i=1	s i	2	+ n	ω i=1	s i	2n

Table 7

 7 Randomly generated UBQP problem, n = 9, relaxation order ω = 2, . . . , 5: problem dimensions and CPU times for Loraine-direct and Loraine-iterative

		Loraine-direct for (6)		Loraine-iterative for (11)
	ω	vars	matrix size	time	vars	matrix size	lin constr	CG iter	time
	2	255	46	0.09	1081	46	510	99	0.08
	3	465	130	0.77	8515	130	930	203	0.97
	4	510	256	2.87	32896	256	1020	302	15.7
	5	511	382	7.04	73153	382	1022	701	129

github.com/kocvara/Loraine.m

https://github.com/kocvara/ADMM for SDP

Solving relaxations of unconstrained BQP by interior-point method

Acknowledgments. This work has been supported by European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Actions, grant agreement 813211 (POEMA). The second author has been partly supported by the Czech Science Foundation through project No. 22-15524S.

Solving relaxations of unconstrained BQP by interior-point method