Soodeh Habibi

Michal Kočvara

Michael Stingl

Loraine-An interior-point solver for low-rank semidefinite programming

Keywords: Semidefinite optimization, interior-point methods, preconditioned conjugate gradients, truss topology optimization, MAXCUT problem, Lasserre relaxations, sensor network localization AMS CLASSIFICATION 90C22, 90C51, 65F08, 74P05

The aim of this paper is to introduce a new code for the solution of large-and-sparse linear Semidefinite Programs (SDPs) with low-rank solutions and/or low-rank data. We propose to use a preconditioned conjugate gradient method within an interior-point SDP algorithm and an efficient preconditioner fully utilizing the low-rank information. The efficiency is demonstrated by numerical experiments using the truss topology optimization problems, Lasserre relaxations of the MAXCUT problems and the sensor network localization problems.

Introduction

The first efficient solvers for semidefinite optimization emerged more than twenty years ago. All of the general purpose solvers have been using second-order algorithms, predominantly the interior-point (IP) method. In every iteration of such an algorithm, one has to solve a system of linear equations Hx = g of size n × n, where n is the dimension of the unknown vector x in the linear semidefinite programming (SDP) problem

min x∈R n c x subject to n ∑ i=1 x i A (k) i -B (k) 0 , k = 1, . . . , p with A (k) i , B (k) ∈ R m×m .
The assembly and solution of this linear system is a well-known bottleneck of these solvers. This is the case even for problems with sparse data and sparse linear systems when the assembly of the system requires O pnm 3 flops, and the solution by sparse Cholesky factorization O (n α ) flops with some α ∈ [START_REF] Achtziger | On the maximization of the fundamental eigenvalue in topology optimization[END_REF][START_REF] Bellavia | A relaxed interior point method for low-rank semidefinite programming problems with applications to matrix completion[END_REF].

CONTACT M. Kočvara. Email: m.kocvara@bham.ac.uk On the other hand, many real-world applications lead to large-scale SDP, often unsolvable by current general-purpose software. There are, essentially, three ways how to approach this conundrum:

(1) Reformulating the problem to make it more suitable for general-purpose solvers. This can be done, for instance, by facial reduction [START_REF] Krislock | Explicit sensor network localization using semidefinite representations and facial reductions[END_REF][START_REF] Permenter | Partial facial reduction: simplified, equivalent SDPs via approximations of the PSD cone[END_REF][START_REF] Zhu | Sieve-SDP: a simple facial reduction algorithm to preprocess semidefinite programs[END_REF] or by decomposition of large matrix inequalities (and thus reducing m) into several smaller ones [START_REF] Kim | Exploiting sparsity in linear and nonlinear matrix inequalities via positive semidefinite matrix completion[END_REF][START_REF] Kočvara | Decomposition of arrow type positive semidefinite matrices with application to topology optimization[END_REF]. (2) Using a different algorithm, such as spectral bundle [START_REF] Helmberg | A spectral bundle method for semidefinite programming[END_REF], ADMM [START_REF] Oliveira | ADMM for the SDP relaxation of the QAP[END_REF], augmented Lagrangian [START_REF] Malick | Regularization methods for semidefinite programming[END_REF][START_REF] Zhao | A Newton-CG augmented lagrangian method for semidefinite programming[END_REF], optimization on manifolds [START_REF] Journée | Low-rank optimization on the cone of positive semidefinite matrices[END_REF], randomized algorithms [START_REF] Yurtsever | Scalable semidefinite programming[END_REF] or techniques of nonlinear programming [START_REF] Burer | A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization[END_REF][START_REF] Burer | Local minima and convergence in low-rank semidefinite programming[END_REF]. (3) Using one of the second-order algorithms with Cholesky factorization replaced by an iterative method for the solution of the linear systems Hx = g; see, e.g., [START_REF] Kočvara | On the solution of large-scale SDP problems by the modified barrier method using iterative solvers[END_REF][START_REF] Toh | Solving large scale semidefinite programs via an iterative solver on the augmented systems[END_REF][START_REF] Zhang | Modified interior-point method for large-and-sparse low-rank semidefinite programs[END_REF].

The last approach is particularly attractive whenever n m, and addresses both bottlenecks of a second-order SDP solver.

• The assembling of the system matrix: an iterative solver only needs matrix-vector multiplications, and so the matrix does not have to be explicitly assembled and stored. • The solution of the linear system itself: an iterative solver can handle very large systems, as compared to a sparse Cholesky solver, under the assumption of good conditioning of the matrix or existence of a good preconditioner.

When the data matrices A (k)

i are of very low rank (such as rank one), the complexity of the linear system assembling can be substantially reduced by using this fact. This is particularly true for rank-one dense matrices. This fact is, however, rarely utilized in standard SDP software, due to the complications related to a different data input.

Our goal in this paper is to introduce Loraine, a new general-purpose interior-point SDP solver targeted to problems with low-rank data and low-rank solutions. It employs special treatment of low-rank data matrices and, in particular, an option to use an iterative Krylov type method for the solution of the linear system. As the choice of an efficient preconditioner in a Krylov type method is problem dependent, in particular, in the context of optimization algorithms. We will focus on problems with expected very low-rank solutions.

SDP problems with low-rank solutions are common in relaxations of optimization problems in different areas such as combinatorial optimization [START_REF] Benson | Solving large-scale sparse semidefinite programs for combinatorial optimization[END_REF], approximation theory [START_REF] Candès | The power of convex relaxation: Near-optimal matrix completion[END_REF][START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF], control theory [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF], structural optimization [START_REF] Stingl | A sequential convex semidefinite programming algorithm with an application to multiple-load free material optimization[END_REF], and power systems [START_REF] Madani | ADMM for sparse semidefinite programming with applications to optimal power flow problem[END_REF].

It is not our goal to find a low-rank solution in case of non-unique solutions; indeed, it is well-known that an interior-point method will converge to a maximal complementary solution. However, if the method will converge to a low-rank solution, we will use this fact to improve its computational complexity and memory demands. The preconditioner is not only limited to problems with (theoretically exactly) low-rank solutions, it can be efficiently applied when the solution has a few outlying eigenvalues.

There is one major difference to other algorithms for SDP problems with low-rank solutions, such as SDPLR [START_REF] Burer | A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization[END_REF][START_REF] Burer | Local minima and convergence in low-rank semidefinite programming[END_REF], optimization on manifolds [START_REF] Journée | Low-rank optimization on the cone of positive semidefinite matrices[END_REF] and other approaches [START_REF] Bellavia | A relaxed interior point method for low-rank semidefinite programming problems with applications to matrix completion[END_REF][START_REF] Yurtsever | Scalable semidefinite programming[END_REF]. All these methods require the knowledge of a guaranteed (possibly tight) upper bound on the rank of the solution. If the estimate of the rank was lower than the actual rank, the algorithm would not find a solution. In our approach, the rank information is merely used to speed up a standard, well-established algorithm. If our estimated rank is smaller than the actual one, we will only see it in possibly more iterations of the Krylov solver; the convergence behaviour of the optimization algorithm will be unchanged.

The actual optimization algorithm used in Loraine is a "standard" primal-dual interiorpoint method. The algorithm mimics the reliable primal-dual algorithm with Nesterov-Todd direction, as used, e.g., in the SDPT3 solver [START_REF] Todd | On the Nesterov-Todd direction in semidefinite programming[END_REF][START_REF] Toh | Solving large scale semidefinite programs via an iterative solver on the augmented systems[END_REF].

Our preconditioner is motivated by the work of Zhang and Lavaei [START_REF] Zhang | Modified interior-point method for large-and-sparse low-rank semidefinite programs[END_REF]. They present a preconditioner for the CG method within a standard IP method that makes it more efficient for large-and-sparse low-rank SDPs. We partly follow their approach though there are some major differences: we abandon a crucial assumption on the sparsity of the data matrices (Assumption 2 in [START_REF] Zhang | Modified interior-point method for large-and-sparse low-rank semidefinite programs[END_REF]); our method allows for SDPs with more than one LMI and with (simple) linear constraints; the complexity of computing our preconditioner is substantially lower; we can also efficiently solve problems with a few outlying eigenvalues in the optimal solution.

To test Loraine and the underlying algorithms, we have generated a library of problems arising from truss topology optimization, sensor network localization and Lasserre relaxations of the MAXCUT problem. These problems are characterized by very low-rank solutions and scalability, in the sense that one can generate a range of problems of various dimension but of the same type. By comparing Loraine with other SDP solvers, we will demonstrate the efficiency of our approach.

The paper is organized as follows. Section 2 introduces the SDP problems we want to solve and some basic assumptions. In Sections 3 an IP method is briefly presented. Next, in Section 4, we introduce the preconditioners used with the iterative solvers. Section 5 presents the solver Loraine. Sections 6-8 describe our applications, truss topology optimization problems, the sensor network localization problem and Lasserre hierarchies of the MAXCUT problem. Finally, in Section 9, we present the results of our numerical experiments. We also give comparisons with other existing SDP software and demonstrate the high efficiency of our solver.

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 eigenvalues of a given matrix X ∈ S m ++ are ordered as

λ 1 (X) ≤ • • • ≤ λ m (X)
, and its condition number is defined as κ(X) = λ 1 (X)/λ m (X). The notation "vec" and "⊗" refer to the (non-symmetrized) vectorization and Kronecker product, respectively, while "svec" and "smat" refer to symmetrized vectorization and its inverse operation. With this notation the following identity holds true: vec(AXB ) = (A ⊗ B)vec(X). The symbol • denotes the Frobenius inner product of two matrices, A • B = tr(A B). Finally, e i ∈ R n is a vector with one on the i-th position and zeros otherwise.

Low-rank semidefinite optimization

While Loraine is a general-purpose SDP solver, its main feature is the ability to efficiently solve SDP problems with low-rank solutions. The first part of the paper is thus focused on this feature.

We consider a (primal) semidefinite optimization problem with matrix variables X i ∈ S m i , i = 1, . . . , p, and with linear constraints min

X 1 ,...,X p , x lin p ∑ i=1 C i • X i + d x lin (1) subject to p ∑ i=1 A (i) j • X i + (D x lin ) j = b j , j = 1, . . . ,
j + S i = C i , i = 1, . . . , p Dy + s lin = d S i 0, i = 1, . . . , p s lin ≥ 0 .
Here

A (i) j ∈ S m i , j = 1, . . . , n, i = 1, . . . , p, C i ∈ S m i , i = 1, . . . , p, b ∈ R n , D ∈ R ν×n , d ∈ R ν are data of the problem.

Basic assumptions

We make the following assumptions.

Assumption 2.1. There exist strictly feasible

X i ∈ S m i ++ , x lin ∈ R ν ++ , y ∈ R n , S i ∈ S m i ++ , s lin ∈ R ν
++ , i = 1, . . . , p satisfying the equality constraints in ( 1) and ( 2) (Slater's condition).

Assumption 2.2. Define the matrices

A i = [vecA (i) 1 , . . . , vecA (i)
n ], i = 1, . . . , p. We assume that matrix-vector products with A i and A i may each be computed in O(n) flops and memory.

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

Assumption 2.4. The dimensions of X i , i = 1, . . . , p, are much smaller than the number of constraints, i.e., m i n.

Remark 1. Assumption 2.3 is satisfied, for instance, in case of box constraints on variable y in the dual formulation. Assumption 2.4 is not really needed as a "mathematical" assumption, the iterative methods presented in the paper would work without it. However, the complexity of these methods would then be comparable with standard software using Cholesky factorization. The presented methods will be only superior to the standard software under this assumption.

Low-rank assumption

As mentioned in the Introduction, our approach is focused on linear systems with matrices having few large outlying eigenvalues. This, of course, includes matrices with very low rank but also matrices with "approximate" low rank. These matrices frequently appear in optimization algorithms, such as interior-point methods, when the "exact" low rank is only attained at the (unreached) optimum. Moreover, the approach allows us to solve problems with matrices of high rank but with a few large outlying eigenvalues, such as the problems in Section 7.

Because our approach covers SDP problems with (genuinely) low-rank solutions and for the lack of better terminology for the larger class of matrices with outlying eigenvalues, we still call our main assumption "low-rank" and the problems we are targeting "SDP problems with low-rank solutions".

Our main assumption used in the design of an efficient preconditioner concerns the eigen-value distribution of the solution X * = (X * 1 , . . . , X * p ) of (1). We assume that X * has a very small number of large outlying eigenvalues. This includes problems where X * i are of a very low rank. Assumption 2.5. Let X * be the solution of [START_REF] Achtziger | On the maximization of the fundamental eigenvalue in topology optimization[END_REF]. We assume that, for i = 1, . . . , p, X * i has k outlying eigenvalues, i.e., that

(0 ≤) λ 1 (X * i ) ≤ • • • ≤ λ (X * i ) m i -k λ (X * i ) m i -k+1 ≤ • • • ≤ λ (X * i ) m i ,
where k is very small, typically smaller than 10 and, often, equal to 1.

Remark 2. In order to emphasize the difference between our algorithm and other "low-rank algorithms", such as SDPLR [START_REF] Burer | A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization[END_REF] and the low-rank preconditioner introduced in [START_REF] Zhang | Modified interior-point method for large-and-sparse low-rank semidefinite programs[END_REF]-that otherwise served as our motivation-we include the following remarks.

(1) We do not assume anything about the distribution of the first (m i -k) eigenvalues of X * i ; they may not be equal to zero, they may not be clustered.

(2) We make no assumption about the "rank" of the solution x * lin of (1), i.e., about the number of active constraints in Dy ≤ d in (2). This can be arbitrarily large or small.

(3) We do not assume that the matrix A A is easily invertible. Indeed, it is not in our application in Section 6, contrary to examples in [START_REF] Zhang | Modified interior-point method for large-and-sparse low-rank semidefinite programs[END_REF] where this matrix is diagonal.

Interior-point method

In this section, we will describe a primal-dual predictor-corrector interior-point method for SDP which is using the Nesterov-Todd (NT) direction, as implemented in Loraine. We will closely follow the articles by Todd, Toh and Tütüncü [START_REF] Todd | On the Nesterov-Todd direction in semidefinite programming[END_REF] and by Toh and Kojima [START_REF] Toh | Solving some large scale semidefinite programs via the conjugate residual method[END_REF] and only repeat formulas needed to present the structure of the matrix that is required to develop the preconditioners in Section 4. The basic framework of our interior-point solver thus, more or less, mimics that of the software SDPT3. We use the notation of problems (1) and (2); however, for the sake of simplicity, we ignore the linear constraints in the following development and consider only one matrix variable, i.e., p = 1. Linear constraints and several matrix variables can, of course, be formally included into a single linear matrix inequality.

Basic framework

Let A be an m 2 × n matrix defined by

A := [vecA 1 , . . . , vecA n ].
Our algorithm follows the standard steps of a primal-dual interior-point method. We write down optimality conditions for (1) with relaxed complementarity condition:

A vec(X) = b, (3)
Ay -S = C, XS = σ µI ,
where µ = X•S m and σ is a centering parameter to be specified below. This system of nonlinear equations is solved repeatedly until convergence.

Newton direction

The system (3) is solved approximately by Newton's method, where in every iteration of the method we solve the following system of linear equations in variables (∆X, ∆y, ∆S):

A vec(∆X) = r p , (4a) A∆y + vec(∆S) = vec(R d ), (4b)
H W (X∆S + ∆XS) = σ µI -H W (XS) . (4c)
Here

r p := b -A X, R d := C -S - n ∑ i=1 y i A i ,
and H W is a linear transformation guaranteeing symmetry of the resulting matrix, in particular

H W (M) = 1 2 W MW -1 +W -M W
with some invertible matrix W ∈ S m ++ known as the scaling matrix [START_REF] Zhang | On extending some primal-dual interior-point algorithms from linear programming to semidefinite programming[END_REF]. By assuming (X, y, S) to be the current iterate, we define the so-called Nesterov-Todd (NT) scaling matrix

W = S -1 2 S 1 2 XS 1 2 1 2 S -1 2 = X 1 2 X -1 2 SX -1 2 1 2 X 1 2
satisfying X = W SW and S = W X -1 W . This is the scaling matrix used in Loraine. Instead of solving the linear system of 2m 2 + n equations (4a)-(4c) directly, we can solve a Schur complement equation (SCE) involving only ∆y. The general bottleneck of any interiorpoint method for SDP is assembling and solving this SCE

H∆y = r. ( 5a
)
Here H is the Schur complement matrix with elements

H i j = A i •WA j W, i, j = 1, . . . , n, (5b)
and

r = r p + A vec (W R d W +W SW ) . ( 5c
)
For details of computing W , r, ∆X and ∆S efficiently for the NT scaling, see [START_REF] Toh | Solving some large scale semidefinite programs via the conjugate residual method[END_REF]. By considering (5a) and (5b), we have the following linear system

(H∆y) i = A i • W n ∑ j=1 ∆y i A j W = r i , i = 1, . . . , n. (6)
Finally, vectorizing the matrix variables allows (6) to be written as

A (W ⊗W )A ∆y = r . (7)

Preconditioners

As we mentioned before, the general bottleneck of the algorithm defined above is the assembling and solving of the Schur complement equation [START_REF] Borchers | SDPLIB 1.2, a library of semidefinite programming test problems[END_REF]. One way to solve this equation is using an iterative method, in particular, a Krylov-type method such as the method of conjugate gradients (CG) or the minimum residual method (MINRES). The Schur complement matrix H is large and the assembling of it may be expensive. Besides, and more importantly, it becomes increasingly ill-conditioned as the IP method makes progress toward the solution. So a successful CG-based solution of our linear systems must rely on an efficient preconditioner. In this section, we introduce two preconditioners implemented in Loraine and targeted pro-SDP problems with low-rank solutions.

Rank of H

The Schur complement matrix H can be written as a sum

H = p ∑ i=1 H i lmi + H lin ,
where H lin is associated with linear constraints. Recall from ( 7) that each H i lmi is computed as

H i lmi = A i (W i ⊗W i )A i , (8)
where W i is the NT scaling matrix.

In the next three Lemmata 4.1-4.3, we will omit the subscript, for simplicity, before returning to the full notation in Theorem 4.4. The next lemma shows that W has the same rank as X. So, as X is low-rank, W will be low-rank. Later in this section, we will use this low-rank property of W to present the preconditioners.

Let (X, S) be a pair of solutions to [START_REF] Achtziger | On the maximization of the fundamental eigenvalue in topology optimization[END_REF]. By complementarity, we know that XS = 0, and rank(X) + rank(S) = k + σ ≤ m. Assume k + σ = m, i.e., strict complementarity. From XS = 0, X 0, S 0, we know that X and S are simultaneously diagonalizable, i.e., there exists U such that X = UΛ X U , S = UΛ S U , where Λ X and Λ S are the diagonal matrices of eigenvalues and U U = I. Lemma 4.1. Let X, S be as above. Let W ∈ S m be any matrix such that X = W SW . Then rank(W ) = rank(X).

Proof. Assume, without loss of generality, that eigenvalues of X and S are sorted such that Λ X = Diag (λ 1 , . . . , λ k , 0, . . . , 0) and Λ S = Diag (0, . . . , 0, µ 1 , . . . , µ m-k ). Then X = W SW is equivalent to UΛ X U = WUΛ S U W , and so, as

U U = I, Λ X = U WUΛ S U WU. Define Z = U WU, so that Λ X = Z Λ S Z. Because (Λ X ) i,i = 0 for i > k, it must hold that z :,i z :,i = 0, i > k, where z :,i is the i-th column of Z.
As Z is positive semidefinite, it means that Z i, j = 0 for i, j > k, i.e., the leading k × k submatrix of Z is a rank-k non-zero matrix and the rest of the matrix is zero. Hence, rank(Z) = k and, as W = UZU , so is the rank of W .

The proof of the next lemma is straightforward.

Lemma 4.2. Let X ∈ S m such that rank X = k, k ≤ m. Then rank(X ⊗ X) = k 2 . Lemma 4.3. Let Y ∈ S m such that rankY = k, k ≤ m, and A ∈ R n×m , n < m. Then rank(AYA ) ≤ k. Proof. Because rankY = k, we have Y = ∑ k i=1 y (i) (y (i) ) with some y (i) ∈ R m , i = 1, . . . , k. Hence AYA = A ∑ k i=1 y (i) (y (i) ) A = ∑ k i=1 z (i) (z (i) ) with z (i) = Ay (i) . Therefore rank(AYA ) = rank ∑ k i=1 z (i) (z (i) ) ≤ k.
A strict inequality occurs, trivially, when n < k. It can also occur when z (i) , i = 1, . . . , k, are linearly dependent.

By combining the above results we get the following theorem.

Theorem 4.4. Let, for some i ∈ {1, . . . , p}, W i be the scaling matrix from Section 3 and let rankW i = k. Let further H i lmi be defined as in [START_REF] Burer | A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization[END_REF]. Then rank H i lmi ≤ k 2 . Remark 3. All results in this section have been presented under the assumption of the exact rank of the involved matrices. This is, of course, only the limit case in our algorithms. The actual matrices have outlying eigenvalues and converge to the low-rank matrices. The message of Theorem 4.4 remains the same, though: when W i has k outlying eigenvalues, then H i lmi will have at most k 2 outlying eigenvalues. This is the fact on which we build the preconditioner.

H α preconditioner

The preconditioner introduced in this section has been motivated by the work of Zhang and Lavaei [START_REF] Zhang | Modified interior-point method for large-and-sparse low-rank semidefinite programs[END_REF] and, in its derivation, we partly follow their paper; however, the new preconditioner differs in a substantial detail, as explained below. Also, our assumption about the matrix W is weaker.

In Section 3, we introduced the scaling matrix W . This matrix becomes progressively illconditioned as the IP method approaches the solution. This ill-conditioning of W is a result of Lemma 4.1. Assuming that the solution matrices X * i , i = 1, . . . , p, have ranks k i , the rank of matrices W i will also tend to k i .

The main idea of the preconditioner is to utilize Assumption 2.5 and decompose matrices W i accordingly to their expected rank as follows:

W i = V s i V l i Λ s i 0 0 τ i I V s i V l i W 0 i +V l i (Λ l i -τI)(V l i ) U i U i (9)
with

τ i satisfying λ 1 (W i ) ≤ τ i < λ m i -k i (W i ).
Here U i are m i × k i matrices of full column rank.

Recall now the form of the Schur complement matrix for problem (1):

H = p ∑ i=1 A i (W i ⊗W i )A i + D X lin S -1 lin D, (10)
in which X lin = diag(x lin ) and S lin = diag(s lin ).

Remark 4. Formally, we could treat the linear constraints, if present, as matrix constraints with diagonal matrices and perform the above decomposition for these constraints, too. How-ever, it is rather unlikely that the corresponding part of the solution, vector x lin would be "low rank", i.e., that only a very few of the linear constraints would be active. This fact could then ruin the whole idea. Here we propose to treat the linear constraints as "full rank" and add the matrix D X lin S -1 lin D to the preconditioner as is. In two of our applications, truss topology optimization and sensor network localization, the linear constraints only consist of upper and lower bounds (in the dual formulation), hence the matrix D X lin S -1 lin D is diagonal and satisfies Assumption 2.3. In the third application, relaxation of the MAXCUT problem, the matrix D X lin S -1 lin D is chordal and sparse ([13, Lemma 2.4]), so sparse Cholesky factorization leads to zero fill-in and is thus very efficient.

By substituting the splittings (9) for i = 1, . . . , p into the matrix [START_REF] Candès | The power of convex relaxation: Near-optimal matrix completion[END_REF], we get

H = p ∑ i=1 A i W 0 i ⊗W 0 i +U i U i ⊗W 0 i +W 0 i ⊗U i U i +U i U i ⊗U i U i A i + D X lin S -1 lin D.
Using the identity A (Φ⊗Ξ)A = A (Ξ⊗Φ)A for any Φ, Ξ ∈ R m i ×m i with A defined as above ([39, Lemma 6] and [15, Chap. 4.2, Problem 25]), we obtain

H = p ∑ i=1 A i (W 0 i ⊗W 0 i )A i + p ∑ i=1 A i (U i ⊗ Γ i )(U i ⊗ Γ i ) A i + D X lin S -1 lin D,
where Γ i is any matrix satisfying

Γ i Γ i = 2W 0 i +U i U i . Next we define V i = A i (U i ⊗ Γ i ) for i = 1, . . . , p. Then, H = p ∑ i=1 A i (W 0 i ⊗W 0 i )A i + D X lin S -1 lin D + Ṽ Ṽ , (11)
where [START_REF] Fawzi | Sparse sums of squares on finite abelian groups and improved semidefinite lifts[END_REF], to define the H α preconditioner, which is

Ṽ = [V 1 , . . . ,V p ]. We now approximate A i (W 0 i ⊗ W 0 i )A i by τ 2 i I in
H α = p ∑ i=1 τ 2 i I + D X lin S -1 lin D A α + Ṽ Ṽ . (12)
By using the Sherman-Morrison-Woodbury (SMW) formula, its inverse will be

H -1 α = A -1 α (I -Ṽ Θ -1 Ṽ A -1 α ) , (13)
where

Θ = I + Ṽ A -1 α Ṽ is the Schur complement.
Notice that, by Assumption 2.3, the inverse of A α is inexpensive.

Complexity of the PCG

• We compute the Schur complement

Θ = I + [A 1 (U 1 ⊗ Γ 1 ) . . . A p (U 1 ⊗ Γ p )] A -1 α [A 1 (U 1 ⊗ Γ 1 ) . . . A p (U p ⊗ Γ p )],
Algorithm 1 Solution of the linear system H∆y = r by PCG Given data k > 0 (solution rank), W i ∈ S m ++ , i = 1, . . . , p, r ∈ R n (right-hand side) and an initial iterate ∆y 0 . for i = 1, . . . , p do computing the decomposition [START_REF] Fawzi | Sparse sums of squares on finite abelian groups and improved semidefinite lifts[END_REF] 3:

Compute eigenvalue decomposition W i = V i Λ i V i and set τ i = λ min (W i ).

4:

Form the matrices W 0 i and U i by (9).

5:

Compute the Cholesky factorization

Γ i Γ i = 2W 0 i +U i U i . 6:
end for Compute the matrix-vector product with H using (6).

11:

Compute the matrix-vector product with H -1 α by means of the SMW in ( 13), using Θ -1 = L -L -1 . 12: end procedure block-wise, with the (i, j)-block Θ i j = B i B j , where

B i = R -1 A i (U i ⊗ Γ i ), i = 1, . . . , p,
and R is the Cholesky factor of A α . Then, we multiply

(C i =)R -1 A i , which is O(n) flops. The sparsity of A i is maintained in C i , and so C i d is still O(n) flops for an arbi- trary vector d. Hence, as (U i ⊗ Γ i ) has m i k i columns, the product C i (U i ⊗ Γ i ) requires O(nm i k i ) flops. Finally, we compute B i B j , which requires O(m i m j k i k j ) flops, if B i
and B j are sparse as well. If this is not the case, the complexity would be O(m i m j k i k j n), however the additional factor n in the complexity formula can be avoided by reordering calculations appropriately. Define m = max • The number of PCG iterations is a little unpredictable, in particular, when approaching the solution of the problem. However, our experience shows that this number is usually well below 10 and does not exceed 100 in the worst problems. Dropping the lowerorder terms and assuming n ≈ m 2 , we will need O m 3 k3 flops in the preconditioner and O m 3 k flops in one PCG iteration. We come to the conclusion that the cost of the preconditioner is about the same as the cost of the solution of the linear system by the PCG method.

H β and hybrid preconditioner

Recall the definition of the matrix H:

H = p ∑ i=1 A i (W 0 i ⊗W 0 i )A i + Ṽ Ṽ + D X lin S -1 lin D .
It turns out that, in our numerical examples, the last term is dominating in the first iterations of the IP algorithm, before the low-rank structure of W is clearly recognized. This is demonstrated in Figure 1 that shows distribution of eigenvalues of matrices D X lin S -1 lin D (in red) and defined as follows

∑ p i=1 A i (W 0 i ⊗ W 0 i )A i + Ṽ Ṽ (in
H β = p ∑ i=1 τ 2 i I + D X lin S -1 lin D, (14)
in which τ i is defined as in the previous section. This matrix is easy to invert by Assumption 2.3; in fact, the matrix is diagonal in many problems. It is therefore an extremely "cheap" preconditioner that is efficient in the first iterations of the IP algorithm.

For relevant problems, we therefore recommend to use a hybrid preconditioner: we start the IP iterations with H β and, once it becomes inefficient, switch to the more expensive but more efficient H α . See more details in Section 9.2 and an example in Section 9.4.2.

The code Loraine

Algorithms presented in the previous sections have been implemented in a code Loraine (for LOw-RAnk INtErior point). The implementation was done in MATLAB and Julia programming environments and is available as open source 1 .

Special features

Loraine is a general-purpose solver for any linear SDP with linear equality and inequality constraints. Compared to other general-purpose SDP software, it particularly targets two classes of problems.

Problems with low-rank solutions

Loraine is using the preconditioned conjugate gradient method, as described in detail in Section 4. In particular, the user can choose between the direct and iterative solver, the type of preconditioner and the expected rank of the solution. The direct solver relies on the implementation of the (sparse or dense) Cholesky factorization provided by MATLAB or Julia.

Problems with low-rank data input

This feature is only useful when a direct solver is used. In this case, the bottleneck of an interior-point algorithm is the computation of the Schur complement matrix in (5b) and [START_REF] Borchers | SDPLIB 1.2, a library of semidefinite programming test problems[END_REF]. When the data matrices A i are of low rank (typically of rank one), the complexity of (5b) can be drastically reduced. The user has a choice -to provide vectors (a j i ) k , k = 1, . . . , r, defining matrices

A (i) j = ∑ r k=1 (a j i ) k (a j i ) k (field Avec in the Data input section below); -to indicate that the input matrices A (i)
j are all expected to be of rank one. Then their decomposition to a vector-vector product can be automatically computed by Loraine.

Data input

The notation in the input file is related to the dual formulation of the problem (2) without slack variables. The input for the MATLAB version of Loraine is a MATLAB structure with the following fields: for inequality constraints; if this vector is not present and ν > 0, inequality type constraints are assumed. A . . . matrices A (i) j stored in MATLAB sparse format as cells A{i,j} Avec . . . vectors (a j i ) k , k = 1, . . . , r, defining matrices

A (i) j = ∑ r k=1 (a j i ) k (a j i ) k C . . . matrix C of linear constraints d . . .

vector d of linear constraints

We also provide several converters, e.g., from SDPA input files. For the Julia version of Loraine, the MATLAB structure can be exported and read by the Julia code.

Loraine options

The following parameters can be changed by the user: kit = 1; % 0..direct solver; 1..iterative solver tol_cg = 1e-2; % initial tolerance for iterative solver tol_cg_up = 0.5; % tolerance update tol_cg_min = 1e-6; % minimal tolerance for CG solver cg_type = 'minres';% 'minres' or 'cg' implemented eDIMACS = 1e-6; % epsilon for DIMACS error stopping criterion prec = 1; % 0..no; 1..H_alpha; 2..H_beta; 4..hybrid erank = 1; % estimated rank of the solution mup = 1000; % initial penalty parameter mu for l1-penalization verb = 2; % 2..full output, 1..short output, 0..no output datarank = 0; % -1..A_i expected rank one; 0..full rank; 1..A_i by vectors initpoint = 0;% 0..Loraine heuristics, 1..SDPT3-like heuristics 6. Application 1: Truss topology optimization

Truss notation

The notation used throughout Section 6 is specific to this application and unrelated to the rest of the paper.

By truss, we understand a mechanical structure, an assemblage of pin-jointed uniform straight bars made of elastic material, such as steel or aluminium. The bars can only carry axial tension and compression. We denote by m the number of bars and by N the number of joints. The positions of the joints are collected in a vector y of dimension n := dim • N where dim is the spatial dimension of the truss. The material properties of bars are characterized by their Young's moduli E i , the bar lengths are denoted by i and bar volumes by t i , i = 1, . . . , m.

Let f ∈ R n be a load vector of nodal forces. The response of the truss to the load f is measured by nodal displacements collected in a displacement vector u ∈ R n . Some of the displacement components may be restricted: a node can be fixed in a wall, then the corresponding displacements are prescribed to be zero. The number of free nodes multiplied by the spatial dimension will be denoted by n and we will assume that f ∈ R n and u ∈ R n .

We introduce the bar stiffness matrices K i and assemble them in the global stiffness matrix of the truss

K(t) = m ∑ i=1 t i K i = m ∑ i=1 t i E i 2 i δ i δ i , i = 1, . . . , m (15)
with "position vectors" δ i ∈ R n , i = 1, . . . , m, and Young's moduli E i ∈ R, i = 1, . . . , m, characterizing the stiffness of the material. The truss must satisfy the equilibrium equation

K(t)u = f . (16)
Assumption 6.1. K(1) with 1 = (1, 1, . . . , 1) ∈ R m is positive definite and the load vector f is in the range space of K(1).

Truss topology problem, rank of the solution

Let 0 ≤ t i ≤ t i , i = 1, . . . , m, and γ be a positive constant. The basic "minimum volume" truss topology optimization (TTO) problem reads as follows:

min t∈R m m ∑ i=1 t i (17)
subject to

γ -f -f K(t) 0 t i ≤ t i ≤ t i , i = 1, . . . , m .
The dual to ( 17) can be written as max

X∈S n+1 ρ∈R m ,ρ∈R m -γ f f 0 • X - m ∑ i=1 ρ i t i + m ∑ i=1 ρ i t i (18)
subject to Using SDP complementarity, it is straightforward to prove the following result.

0 0 0 K i • X -ρ i + ρ i = 1, i = 1, . . . , m X 0 ρ i ≥ 0, ρ i ≥ 0, i = 1, . . . ,
Theorem 6.2. There exists a solution (X * , ρ * , ρ * ) ∈ S n+1 ×R m ×R m of the dual SDP problem [START_REF] Kim | Exploiting sparsity in linear and nonlinear matrix inequalities via positive semidefinite matrix completion[END_REF] such that the rank of X * is one. For t i > 0 this solution is unique.

Truss topology problem, vibration constraints

Formulation of the basic truss topology problem as an SDP is rather academic. It was shown, e.g., in [START_REF] Jarre | Optimal truss design by interior-point methods[END_REF] that it is equivalent to a convex nonlinear programming problem that can be solved very efficiently by interior point methods. The SDP formulation gains significance once we add more, important and practical, constraints to the problem. In particular, it was shown in [START_REF] Achtziger | On the maximization of the fundamental eigenvalue in topology optimization[END_REF][START_REF] Achtziger | Structural topology optimization with eigenvalues[END_REF] that a constraint on natural (free) vibrations of the optimal structure leads to a linear matrix inequality and that this is, arguably, the best way how to treat this constraint. The free vibrations are the squares of the eigenvalues of the generalized eigenvalue problem

K(t)w = λ (M(t) + M 0 ) w , (19)
where

M(t) = m ∑ i=1
t i M i is the so-called mass matrix that collects information about the mass distribution in the truss. The matrices M i are positive semidefinite and have the same sparsity structure as K i . The non-structural mass matrix M 0 is a constant, typically diagonal matrix with very few nonzero elements.

Low vibrations are dangerous and may lead to structural collapse, hence a typical vibration constraint is λ min ≥ λ for a given λ > 0 where λ min is the smallest eigenvalue of problem [START_REF] Kočvara | Decomposition of arrow type positive semidefinite matrices with application to topology optimization[END_REF]. This constraint can be equivalently written as a linear matrix inequality K(t)λ (M(t) + M 0 ) 0. We will thus get the following SDP formulation of the truss topology design with a vibration constraint:

min t∈R m m ∑ i=1 t i (20) subject to γ -f -f K(t) 0
K(t) -λ (M(t) + M 0 ) 0 t i ≤ t i ≤ t i , i = 1, . . . , m .
It has been shown, e.g., in [START_REF] Stingl | A sequential convex semidefinite programming algorithm with an application to multiple-load free material optimization[END_REF] that, when t i > 0, the rank of the optimal dual variable to the second matrix inequality is equal to the multiplicity of the smallest eigenvalue of ( 19). This multiplicity depends on the geometry of the optimal structure but is, typically, very low, usually not bigger than 1 or 2. Problem [START_REF] Kočvara | On the solution of large-scale SDP problems by the modified barrier method using iterative solvers[END_REF] thus fits in our framework of SDPs with very low-rank solutions.

Application 2: Sensor network localization with noisy data

The problem

Assume there are m distinct (sensor) points in x i ∈ R d , whose locations are to be determined, and ν other fixed (anchor) points, whose locations a 1 , . . . , a ν are known. The Euclidean distance d i j between the ith and jth sensor point is known if (i, j) ∈ I x = {(i, j) | x ix j = δ i j ≤ ρ}, where ρ is a fixed parameter called the radio range. Similarly, the Euclidean distance d k j between the kth anchor and jth sensor point is known if

(k, j) ∈ I a = {(k, j) | a k -x j = δ i j ≤ ρ}. The sensor network localization (SNL) problem is to Find x i ∈ R d , i = 1, . . . , m, for which (21) x i -x j 2 = δ 2 i j , (i, j) ∈ I x a k -x j 2 = δ 2 k j , (k, j) ∈ I a .
Typical networks of this type consist of a large number of densely deployed sensor nodes which gather local data and communicate with other nearby nodes. The sensor data from the nodes are relevant only if we know to what location they refer. The problem arises in applications as different as the habitat monitoring system in the Great Duck Island, detecting volcano eruptions, industrial control in semiconductor manufacturing plants, structural health monitoring, battlefield surveillance, moving object tracking, asset location and the problem of determining protein molecule structure [START_REF] Biswas | Semidefinite programming approaches for sensor network localization with noisy distance measurements[END_REF]. Figure 2-left shows an example of a network with 9 anchors and 512 sensors. In many applications the data δ i j are noisy. Assume perturbation of δ i j and δ ik by random noises ε i j and ε ik , respectively:

δ i j = δi j |1 + ε i j |, (i, j) ∈ I x δ ik = δ ik |1 + ε ik |, (i, j) ∈ I a
where δi j are the true distances. Obviously, the problem ( 21) may now become infeasible.

SDP relaxation

In a series of papers, Biswas, Ye and their co-authors [START_REF] Biswas | Semidefinite programming approaches for sensor network localization with noisy distance measurements[END_REF][START_REF] Biswas | Semidefinite programming for ad hoc wireless sensor network localization[END_REF][START_REF] Wang | Further relaxations of the semidefinite programming approach to sensor network localization[END_REF] proposed a relaxation of the problem based on the semidefinite programming formulation. The next lines follow the development from Biswas and Ye [START_REF] Biswas | Semidefinite programming for ad hoc wireless sensor network localization[END_REF]. Let X = [x 1 x 2 . . . x m ] be a d × m unknown matrix. Then

x i -x j 2 = (e i -e j ) T X T X(e i -e j ) a k -x j 2 = (a k ; -e j ) T I d X T [I d X](a k ; -e j )
and the original problem ( 21) can be equivalently written as

Find x i ∈ R d , i = 1, . . . , m, for which ( 22
)
(e i -e j ) T X T X(e i -e j ) = δ 2 i j (a k ; -e j ) T I d X X T Y (a k ; -e j ) = δ 2 k j Y = X T X .
We now relax the equality Y = X T X to Y X T X, which is equivalent to the linear matrix inequality

Z = I d X X T Y 0 .
The semidefinite optimization problem min

Z,u,v ∑ (i, j)∈I x (u i j + v i j ) + ∑ (k, j)∈I a (u k j + v k j ) (23)
subject to

Z 1:d,1:d = I d (0; e i -e j ) T Z(0; e i -e j ) -u i j + v i j = δ 2 i j , (i, j) ∈ I x (a k ; -e j ) T Z(a k ; -e j ) -u k j + v k j = δ 2 k j , (k, j) ∈ I a Z 0 u i j , v i j ≥ 0, (i, j) ∈ I x , u k j , v k j ≥ 0, (k, j) ∈ I a
is called the full SDP relaxation (FSDP) of ( 21). It has been shown by So and Ye [START_REF] So | Theory of semidefinite programming for sensor network localization[END_REF] that for certain "good" problems with no data noise this relaxation is exact.

For larger sensor network localization problems, the FSDP relaxation is difficult to solve numerically. The SDP problem has a large number of variables and the constraint matrix is thus large and full. In order to solve the problem, one may exploit the sparsity of I x and I a at the relaxation modelling level. Several approaches have been published on this subject. A simple, yet powerful approach was proposed by Wang et al. [START_REF] Wang | Further relaxations of the semidefinite programming approach to sensor network localization[END_REF] and called edge-based relaxation (ESDP). However, ESDP does not maintain the localizability of FSDP: the solution can have high rank even for uniquely localizable problems. For problems with no data noise, Krislock and Wolkowicz [START_REF] Krislock | Explicit sensor network localization using semidefinite representations and facial reductions[END_REF] proposed a very efficient algorithm based on facial reductions. This approach, however, cannot be extended to the problems with noisy data.

As mentioned above, for many noiseless problems the relaxation is exact, i.e., the rank of the matrix Z is two. Then, for problems with small noise, Z is expected to have two large eigenvalues and several, perhaps many, small nonzero ones. The problem thus satisfies our low-rank Assumption 2.5 and can be efficiently solved by Loraine. This will be demonstrated in Section 9.

Application 3: Lasserre relaxations of the MAXCUT problem

The problem

Let Γ be an undirected n-node graph and let the arcs (i, j) be associated with nonnegative weights a i j . The task is to find a cut of the largest weight, that is, to partition the set of nodes into two parts, S and S such that the total weight of all arcs linking S and S is as large as possible.

The MAXCUT problem is formulated as a quadratic optimization problem max

x 1 4 n ∑ i, j=1 a i j (1 -x i x j ) | x 2 i = 1, i = 1, . . . , n (24)
or, equivalently, with X = xx T , as a linear SDP with a rank-one constraint:

max X 1 4 n ∑ i, j=1 a i j (1 -X i j ) | X ii = 1, i = 1, . . . , n; rank X = 1; X 0 . (25)
It is well-known that problem [START_REF] Madani | ADMM for sparse semidefinite programming with applications to optimal power flow problem[END_REF] is NP-hard. The celebrated result by Goemans and Williamson [START_REF] Goemans | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF] shows that the (Shor's) relaxation of ( 25) by ignoring the rank constraints will deliver a solution with objective function not bigger than 1 0.87865 times the global optimum. The relaxation is, however, almost never exact which means that the solution X * of the relaxed (25) has, typically, an unknown rank bigger than one.

Tighter approximations with low-rank solutions can be obtained by higher order relaxations, introduced by Lasserre [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]. The original problem ( 24) is a polynomial optimization problem and we can just use the Lasserre machinery to obtain tight approximations of its global solution. It was shown by Fawzi et al. [START_REF] Fawzi | Sparse sums of squares on finite abelian groups and improved semidefinite lifts[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 ( 24) 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 unweighted complete graphs.

The SDP relaxations can be written in the following form:

min y∈R s y q (26) subject to s ∑ i=1 M i y i + I 0 .
Here, for the relaxation of order one, Q is the Laplacian matrix of the graph, q = svec(Q), y = svec(X), s = n(n + 1)/2, and M i ∈ R n+1×n+1 are the (pseudo-)Hankel matrices; in this case ( 26) is just re-writing of the relaxed problem (25) (without the rank constraint). For higher-order relaxations, the dimensions s and m get bigger by adding higher-order terms in M i and y; for details, see [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]. For an exact relaxation (i.e., relaxation delivering the global solution of ( 24)), the rank of the optimal moment matrix M * := ∑ s i=1 M i y * i + I is equal to the number of global solutions; in particular, if the global solution is unique, the rank of M * is one. So, while problem [START_REF] Mittelmann | An independent benchmarking of SDP and SOCP solvers[END_REF] has a low-rank solution, it is in the "wrong" form for our IP solver. While the solution M * of ( 26) is of low rank, our solver expects the solution of an SDP problem in the primal form (1) to be of low rank. In [START_REF] Habibi | An interior-point method for Lasserre relaxations of unconstrained binary quadratic optimization problems[END_REF] we show that the dual to (26) can be re-written in the dual form

max z∈R n (svec(I)) z (27)
subject to smat(z) 0

Mz = q ,
where n = m(m + 1)/2 and M = (svec(M 1 ), . . . , svec(M s )) T , M ∈ R s× n . Now, due to complementarity, we know that the dual solution to this problem (the Lagrangian multiplier to the constraint smat(z) 0 and the solution to [START_REF] Mittelmann | An independent benchmarking of SDP and SOCP solvers[END_REF]) is of low rank, assuming high enough relaxation order and unique solution of [START_REF] Madani | ADMM for sparse semidefinite programming with applications to optimal power flow problem[END_REF]. Also in [START_REF] Habibi | An interior-point method for Lasserre relaxations of unconstrained binary quadratic optimization problems[END_REF] we show how to treat the equality constraints in [START_REF] Mosek Aps | The MOSEK optimization toolbox for MATLAB manual[END_REF] in the interior point algorithm; in particular, we use the 1 -penalty method for this purpose. Problem [START_REF] Mosek Aps | The MOSEK optimization toolbox for MATLAB manual[END_REF] now satisfies the assumptions needed for the low-rank preconditioner.

Numerical experiments

Problem database

To test Loraine and the underlying algorithms, we generated a library of problems arising from truss topology optimization, sensor network localization and Lasserre relaxations of the MAXCUT problem. These problems are characterized by very low-rank solutions and scalability, in the sense that one can generate a range of problems of various dimension but of the same type.

TTO problems

We generated a set of truss topology optimization problems of various sizes, both with and without the vibration constraint. All problems have the same geometry, boundary conditions and loads and only differ in the number of potential nodes and bars. The geometry and loading for TTO problems without vibration (formulation ( 17)) are as shown in Figure 3-left, while the data for problems with the vibration constraint (formulation ( 20)) differ only in the horizontal orientation of the load vector. We consider two groups of problems with the following naming convention:

tru<n> standard TTO problem [START_REF] Journée | Low-rank optimization on the cone of positive semidefinite matrices[END_REF] with t = 0, discretized by n × n nodes; vib<n> same as above but with the vibration constraint [START_REF] Kočvara | On the solution of large-scale SDP problems by the modified barrier method using iterative solvers[END_REF].

In each group, all nodes in the initial ground structure are connected by potential bars. Table 1 presents the dimensions of the generated problems. For illustration, Figure 3 shows optimal results for problems tru25 and vib19.

SNL problems

We consider SNL problems of growing dimension with nine anchors and n sensors randomly distributed around the letters O and B; see Figure 2. The radio range ρ (see the beginning of Section 7) was set to 1 5 of the maximal vertical distance of the anchors. The known distances between the sensors were perturbed by 1% Gaussian noise. Table 2 (left) shows the dimensions of the generated problems.

MAXCUT problems

Recall that the solution of the MAXCUT problem is of rank one only if the problem has a unique global solution. Obviously, a complete unweighted graph may have many "symmetric" solutions. 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. The dimensions of the corresponding SDP problems (27) (for relaxation order two) are shown in Table 2 (right).

Loraine setting

The algorithms presented in this article were implemented in MATLAB and Julia codes. The codes use DIMACS stopping criteria [START_REF] Mittelmann | An independent benchmarking of SDP and SOCP solvers[END_REF] with ε DIMACS = 10 -5 . The stopping parameter for the CG method ε CG in Algorithm 1 has been set and updated as follows: we start with ε CG = 0.01 and multiply it after every major iteration by 0.5, until it reaches the value 10 -6 ; then we continue with this value till convergence of the optimization algorithm.

Unless stated otherwise, we use the following hybrid preconditioner called H hyb : we start with H β [START_REF] Helmberg | A spectral bundle method for semidefinite programming[END_REF] until the criterion below is satisfied and then switch to H α [START_REF] Goemans | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF]:

N cg_iter > kp √ n/10 & N iter > √ n/60 .
Here N iter is the current IP iteration and N cg_iter the number of CG steps needed to solve the corrector equation in this iteration. This criterion is, obviously, purely heuristic. We further choose the value of τ in H α and H β as

τ = λ 1 (W i ) + 0.5 mean(λ 1 (W i ), . . . , λ m-k i (W i )) .
To solve the problems, we used the default setting of Loraine parameters as shown in Section 5, with the following exceptions, as explained in the respective paragraphs below: SNL: tol_cg = 1e-1, eDIMACS = 1e-4, erank = 12, initpoint = 1 MAXCUT: eDIMACS = 1e-5

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.

Other software

Later in this section, we present results obtained by Loraine, together with (sometimes brief) comparison with other SDP software. For this comparison, we use the solver MOSEK [START_REF] Mosek Aps | The MOSEK optimization toolbox for MATLAB manual[END_REF] which is also based on an interior-point algorithm. While this comparison may seem unfair, as MOSEK uses a direct solver to solve the linear systems, we believe that it gives a good perspective of the clear advantage of iterative solvers with good preconditioners, when available and when applicable. We further present a comparison with our MATLAB implementation of the ADMM method, based on [START_REF] Wen | Alternating direction augmented lagrangian methods for semidefinite programming[END_REF].

For the truss problems, we also offer a comparison with SDPNAL+ [START_REF] Yang | SDPNAL+: a majorized semismooth Newton-CG augmented Lagrangian method for semidefinite programming with nonnegative constraints[END_REF], an implementation of the semi-smooth Newton-CG augmented Lagrangian method; and SDPLR [START_REF] Burer | A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization[END_REF] using a nonlinear programming algorithm based on low-rank factorization of the variables.

Numerical results

Truss topology problem

We now present results for all tru and vib problems using Loraine. These are shown in Table 3. The table presents the number of iterations of the respective algorithm and the total number of CG iterations using the hybrid preconditioner H hyb . This is followed by CPU time spent in the optimization solver and time per iteration. Notice that Loraine solves two linear systems per iteration (predictor and corrector). For all problems we set the expected rank of the dual solutions k = 1. The next Table 4 presents results of Loraine with direct solver and rank-one input data, MOSEK, SDPNAL+ and SDPLR; we again give the number of iterations of the solver, CPU time and (for Loraine and MOSEK) CPU time per one iteration. The data matrices for the truss problems are of rank one; see [START_REF] Horn | Topics in matrix analysis[END_REF]. While Loraine is much more efficient for these problems with an iterative solver, the goal of this comparison is to show that when handling the rank-one data efficiently, we get a code that is (at least per one iteration) faster than other software even with a direct solver. We only present results for problems tru3-tru17; for the larger problems MOSEK exceeded the available 64GB RAM. SDPNAL+ and SDPLR were used with default stopping tolerance. None of these codes can solve problems bigger than tru11 before reaching the maximum number of iterations, either internal or global. While Loraine solutions satisfy DI-MACS criteria with 10 -5 , other codes often finish with lower precision. The ADMM is not included in this comparison as, even for the smallest problems, the method stagnated at a point far from the solution.

In Figure 4-left we give a comparison of the CPU time per iteration for Loraine (iterative) and MOSEK. While MOSEK actually exceeds the theoretical complexity with coefficient 2.2, the complexity of Loraine is almost linear for these problems.

Sensor network localization problems

Recall from Section 7 that problems with noisy data lead to optimal solution Z * of [START_REF] Laurent | Lower bound for the number of iterations in semidefinite hierarchies for the cut polytope[END_REF] with two outlying and a number of small, nonzero eigenvalues. Selecting k = 2 therefore did not lead to an efficient preconditioner and we had to increase the estimate of the rank of Z * . For the problems we solved, k = 12 turned out to be a good compromise between the quality of the preconditioner and the CPU time needed to apply it. Due to ill-conditioning of the problems, we relaxed the stopping criterion to eDIMACS=1e-4; the same criterion was used for the ADMM code.

Table 5 presents the results for Loraine, MOSEK and ADMM. While ADMM turned out to be more successful than for the truss problems, it is still much slower than the two IP codes. In this table, the number of CG iterations sums up the iterations with preconditioners H β and H α . For instance, in problem SNL-64, the average number of CG iterations with H α was 10 and the maximal number (in the last IP iteration) was 23. In Figure 4-right we give a comparison of the CPU time per iteration for Loraine and MOSEK. For these problems, we also demonstrate the usefulness of the hybrid preconditioner (see Sections 4 and 9.2). Figure 5 shows the iteration log for problem SNL-80; the first IP iterations are solved using H β ; once it requires too many iterations, we switch to H α .

MAXCUT problems

The results obtained for the MAXCUT problems are summarized in Table 6. These problems are relatively "simple": the ADMM method is rather efficient now (this is well-known from other studies) and MOSEK converges in a very small number of iterations. However, the sheer size prevents MOSEK with a direct solver to solve larger problems, while ADMM is still about three times slower than Loraine. Moreover, the ADMM behaviour is rather unpredictable, regarding the number of iterations.

Conclusions

We have demonstrated the efficiency of Loraine for three classes of well-known SDP problems. While we tested our software for other problems from standard libraries, we do not report it here for one of the following reasons:

(1) The available collections, such as SDPLIB [START_REF] Borchers | SDPLIB 1.2, a library of semidefinite programming test problems[END_REF], are not representative enough in the sense that they do not include enough problems of the same type. We did not want to pick up single problems. (2) Loraine mainly benefits from the use of an iterative solver for the linear system using a specific preconditioner targeted to problems with low-rank solutions. For many such problems, however, no preconditioner is needed for the convergence of the iterative method; these "easy" problems were not included in our testing. This is the case of matrix completion problems [START_REF] Bellavia | A relaxed interior point method for low-rank semidefinite programming problems with applications to matrix completion[END_REF][START_REF] Zhang | Modified interior-point method for large-and-sparse low-rank semidefinite programs[END_REF] or some problems arising in the relaxation of combinatorial optimization. It is not the case of problems from Section 9; for those, the preconditioner is essential for convergence. (3) The problems do not satisfy our dimensional assumption n m. For these problems, the preconditioner is still very efficient (assuming low-rank solutions) but the use of iterative methods does not bring any significant advantage to direct solvers. This is the case of many problems found in SDPLIB.

7 :

 7 Compute the Schur complement Θ = I + Ṽ A -1 α Ṽ and its Cholesky factorization Θ = LL . 8: end procedure 9: procedure PCG Use standard PCG algorithm to solve H∆y = r until H∆yr / r ≤ ε CG .At each PCG iteration:10:

  i=1,...,p m i and k = max i=1,...,p k i . Forming and factorizing the Schur complement (i.e., forming the preconditioner) requires O p 2 m 2 k2 + pn m k + m 3 k3 flops; the last term comes from the Cholesky factorization of Θ. • Each iteration of the PCG method is dominated by the matrix-matrix product in (6) requiring O ∑ p i=1 m 3 i flops and application of the preconditioner (13) requiring O ∑ p i=1 m i k i 2 + O n ∑ p i=1 m i k i flops.

Figure 1 .

 1 Figure 1. Problem tru7e; eigenvalues of D X lin S -1 lin D (red) and ∑ p i=1 A i (W 0 i ⊗W 0 i )A i + Ṽ Ṽ (blue) in iteration 2, 17 and 30 (left to right).

  type . . . 'sdp' [field used for future extensions] name . . . 'name of the problem' nvar . . . number of variables n c . . . objective vector b nlmi . . . number of linear matrix inequalities p msizes . . . dimensions of the LMIs m 1 , . . . , m p nlin . . . number of linear constraints ν lsi_op . . . (optional) binary vector of length nlin, set to zero for equality constraints and one

Figure 2 .

 2 Figure 2. Sensor network localization problem: initial data (left) and optimal solution of the SDP relaxation (right). Black: anchors; magenta: exact positions of sensors; red (left): known distances between anchors and sensors; blue (left): known mutual distances of sensors; blue (right): computed positions of sensors.

Figure 3 .

 3 Figure 3. Initial design for tru* problems (left); optimal results for problems tru25 (middle) and vib19 (right).

Figure 4 .

 4 Figure 4. Loraine (red) and MOSEK (blue) in tru * problems (left) and SNL * problems (right). CPU time per iteration vs number of variables in log-log scale.*** Loraine v0.1 *** Number of variables: 40839 Preconditioner : 4 Matrix size(s) : 642 Expected rank : 12 Linear constraints : 81674 it obj error cg_iter time/iter 1 4.37357010e+02 3.32e+05 4 0.34 2 4.79694156e+02 7.79e+04 4 0.40 3 4.75844194e+02 8.44e+03 38 0.81 ... 9 1.64854663e+00 7.24e+02 45 0.71 10 -4.55049809e+00 1.73e+02 69 0.91 11 -1.14083700e+01 2.25e+01 112 1.23 Switching to preconditioner 1 12 -1.95288379e+01 8.02e+00 14 2.74 13 -2.43383931e+01 4.02e+00 13 2.65 14 -2.97196245e+01 2.32e+00 14 2.70 ... 25 -4.17128514e+01 1.71e-03 40 4.61 26 -4.17274123e+01 6.84e-04 58 5.95 27 -4.17345508e+01 2.33e-04 81 7.67 28 -4.17372049e+01 7.21e-05 111 9.71 *** Total CG iterations: 1064 *** Total CPU time: 80.70 seconds

Figure 5 .

 5 Figure 5. Problem SNL-80 solved using the hybrid preconditioner.

1 :

 1 procedure SETTING UP PRECONDITIONER H α

	2:

Table 1 .

 1 Problems tru<n> and problems vib<n>, number of variables n, size of the LMI constraint m and number of linear constraints.

	problem	n	m	lin. constr.	problem	n	m	lin. constr.
	tru3	36	13	72	vib3	36	(13, 12)	72
	tru5	300	41	600	vib5	300	(41, 40)	600
	tru7	1176	85	2352	vib7	1176	(85, 84)	2352
	tru9	3240	145	6480	vib9	3240	(145, 144)	6480
	tru11	7260	221	14520	vib11	7260	(221, 220)	14520
	tru13	14196	313	28392	vib13	14196	(313, 312)	28392
	tru15	25200	421	50400	vib15	25200	(421, 420)	50400
	tru17	41616	545	83232	vib17	41616	(545, 544)	83232
	tru19	64980	685	129960	vib19	64980	(685, 684)	129960
	tru21	97020	841	194040	vib21	97020	(841, 840)	194040
	tru23	139656 1013	279312	vib23	139656 (1013, 1012)	279312
	tru25	195000 1201	390000	vib25	195000 (1201, 1200)	390000

Table 2 .

 2 Problems SNL-<n> and MAXCUT-<n>, number of variables n, size of the LMI constraint m and number of linear constraints.

	problem	n	m	lin. constr.	problem	n	m	lin. constr.
	SNL-16 SNL-32 SNL-48 SNL-64 SNL-80 SNL-96 SNL-112 SNL-128 109464 1026 1677 130 6818 258 14886 386 26767 514 40840 642 58956 770 82137 898	3348 13630 29766 53528 81674 117906 164268 218922	MAXCUT-20 MAXCUT-25 MAXCUT-30 108811 22366 53301 MAXCUT-35 199396 MAXCUT-40 337431 MAXCUT-45 537166 1036 211 326 466 631 821 MAXCUT-50 814726 1276	12390 30550 63860 119070 204180 328440 502350

Table 3 .

 3 Loraine in tru<n> and vib<n> problems

	problem iter CG iter time time/iter	problem	iter CG iter	time time/iter
	tru3	16	122	0.01	0.00	vib3	20	209	0.04	0.00
	tru5	21	190	0.08	0.00	vib5	31	411	0.21	0.01
	tru7	27	236	0.24	0.01	vib7	39	501	0.67	0.02
	tru9	31	333	0.65	0.02	vib9	47	663	2.2	0.05
	tru11	36	370	1.6	0.04	vib11	59	995	6.1	0.10
	tru13	45	500	4.5	0.10	vib13	69	1153	16	0.23
	tru15	52	882	11	0.22	vib15	80	1480	37	0.46
	tru17	53	980	24	0.45	vib17	94	1781	92	0.98
	tru19	64	1310	51	0.80	vib19	108	2116	189	1.75
	tru21	65	1325	84	1.29	vib21	120	2844	401	3.34
	tru23	72	2450	189	2.63	vib23	130	2720	718	5.52
	tru25	87	2148	317	3.64	vib25	143	3752 1321	9.24

Table 4 .

 4 Loraine direct and other solvers in tru<n> problems

			Loraine direct		MOSEK		SDPNAL+	SDPLR
	problem iter	time time/iter iter	time time/iter	iter time iter	time
	tru3	16	0.01	0.001	14	0.22	0.02	440	0.5	15	0.01
	tru5	21	0.09	0.004	15	0.23	0.02	9674	29	18	1.7
	tru7	27	0.46	0.02	17	0.78	0.05	6324	55	19	21
	tru9	31	4.3	0.14	20	5.6	0.28 47579	468	20	179
	tru11	36	31	0.86	26	44	1.69	maxit		22 2469
	tru13	45	239	5.31	29	235	8.10			maxit
	tru15	52 1216	23.38	32 1079	33.72
	tru17	51 4583	89.86	37 4671	126.24

Table 5 .

 5 Loraine, MOSEK and ADMM in SNL-<n> problems

			Loraine		MOSEK	ADMM
	problem iter CG iter time iter	time	iter	time
	SNL-16	21	395	1.4	13	1.0	6239	12
	SNL-32	25	532	7.4	15	21 19343	362
	SNL-48	27	676	21	17	150 14328	1050
	SNL-64	27	705	40	16	611 81938	17554
	SNL-80	28	1064	84	21 2052	>50000
	SNL-96	29	969	120	27 6579	time
	SNL-112	30	1293	224	memory
	SNL-128	30	1316	299

Table 6 .

 6 Loraine, MOSEK and ADMM in MAXCUT-<n> problems

			Loraine		MOSEK	ADMM
	problem	iter CG iter time iter	time	iter	time
	MAXCUT-20	18	559	3.8	6	9 3486	15
	MAXCUT-25	20	728	11	7	78 4314	45
	MAXCUT-30	21	1032	28	9	607 4837	105
	MAXCUT-35	23	2183	96	9 2911 3030	126
	MAXCUT-40	27	2275	186	memory	1280	92
	MAXCUT-45	25	2521	335			7976 1034
	MAXCUT-50	24	2540	528			8783 1832

github.com/kocvara/Loraine.m

Funding

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.