Y De Castro

F Gamboa

D Henrion

J.-B Lasserre

EXACT SOLUTIONS TO SUPER RESOLUTION ON SEMI-ALGEBRAIC DOMAINS IN HIGHER DIMENSIONS

Keywords: super resolution, signed measure, semidefinite programming, total variation, semialgebraic domain

We investigate the multi-dimensional Super Resolution problem on closed semi-algebraic domains for various sampling schemes such as Fourier or moments. We present a new semidefinite programming (SDP) formulation of the 1 -minimization in the space of Radon measures in the multi-dimensional frame on semi-algebraic sets.

While standard approaches have focused on SDP relaxations of the dual program (a popular approach is based on Gram matrix representations), this paper introduces an exact formulation of the primal 1 -minimization exact recovery problem of Super Resolution that unleashes standard techniques (such as moment-sum-of-squares hierarchies) to overcome intrinsic limitations of previous works in the literature. Notably, we show that one can exactly solve the Super Resolution problem in dimension greater than 2 and for a large family of domains described by semi-algebraic sets.

1. Introduction 1.1. Super Resolution. The early formulation of the Super Resolution problem can be identified as the ability of faithfully reconstruct a high-dimensional sparse vector from the observation of a low-pass filter. This situation models important applications in imaging spectroscopy [START_REF] Harris | Super-resolution imaging spectroscopy[END_REF], image processing [START_REF] Park | Super-resolution image reconstruction: a technical overview[END_REF], radar imaging [START_REF] Odendaal | Two-dimensional superresolution radar imaging using the music algorithm, Antennas and Propagation[END_REF], or astronomy [START_REF] Makovoz | Point-source extraction with mopex[END_REF]. As a theoretical baseline, suppose that one wants to reconstruct a vector x by solving a system of linear equations: If the number s of non-vanishing components of x is small and the matrix A enjoys some geometric property, namely the Null Space Property [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF], that depends only on its kernel and the sparsity s, then one can exactly reconstruct x by minimizing the 1 -norm within the affine subspace of all solutions of the linear system. Conditions on m, N , s and properties have been extensively studied, see for example [START_REF]For most large underdetermined systems of linear equations the minimal l 1norm solution is also the sparsest solution[END_REF], [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF] and reference therein. Less than ten year ago, Super Resolution has seeded the ideas of compressed sensing theory [START_REF]Compressed sensing[END_REF], [START_REF] Candes | Near-optimal signal recovery from random projections: universal encoding strategies?[END_REF], [START_REF]The Dantzig selector: statistical estimation when p is much larger than n[END_REF]. In this theory, the matrix A is randomized and one is interested both in the construction of probability distributions allowing to show relevant properties, such as the Restricted Isometry Property [START_REF] Candes | Near-optimal signal recovery from random projections: universal encoding strategies?[END_REF], and the stability of the reconstruction process. This area of research is very fruitful and leads to many practical applications in signal and image processing, see for example [START_REF] Herman | High-resolution radar via compressed sensing[END_REF], [START_REF] Tian | Compressed sensing for wideband cognitive radios[END_REF], or [START_REF] Chen | Prior image constrained compressed sensing (piccs): a method to accurately reconstruct dynamic ct images from highly undersampled projection data sets[END_REF]. To the best of our knowledge, the first mathematical works on Super Resolution are due to Donoho et al in the early ninety, see [START_REF] Donoho | Uncertainty principles and signal recovery[END_REF] and [START_REF] Donoho | Superresolution via sparsity constraints[END_REF]. In these papers, the term Super Resolution appeared because the matrix A is related to a discretization of some Fourier transform. As a matter of fact, when inverting a discrete Fourier transform, the separation of two close spikes in a sparse signal is made possible by minimizing the 1 -norm while a linear inversion method is not able to do so. Beside, in these years, many applied researchers were performing Fourier inversion of non negative sparse signals using some entropy regularization, see for example [START_REF] Skilling | Prior knowledge must be used, Maximum-entropy and Bayesian spectral analysis and estimation problems[END_REF] and [START_REF] Burg | Maximum entropy spectral analysis[END_REF]. At this time, the respective roles of sparsity and non negativity in regard of spikes separation from non linear Fourier inversion methods was not completely clear. An important step for understanding these roles has been taken by lifting the linear equation (1.1) up to the more abstract measure set up, see [START_REF] Gassiat | Problèmes des moments et concentration de mesure[END_REF], [START_REF]Sets of superresolution and the maximum entropy method on the mean[END_REF] and [START_REF] Doukhan | Superresolution rates in Prokhorov metric[END_REF]:

(1.2) a i , µ = X a i (x)µ(dx) = b i where µ ∈ R(X) + and i = 1, . . . , m .

Here, (a i (x)) is a vector of continuous function defined on X, a given compact subset of R n , b = (b i ) ∈ R m and R(X) + is the set of all nonnegative Radon measures on X. In this frame, there exists very special points b such that the set of all members of R(X) + satisfying (1.2) for b = b reduces to a singleton {µ b }. Furthermore, µ b is a discrete measure concentrated on very few points. Hence, if one deals with a b close to a point b the set of all solutions of (1.2) (or (1.1) with positivity constraint and a design matrix A discretizing the vectorial function (a i (x)) m i=1 ) is very small. So that, two different methods for selecting a member of this set will lead to similar solutions. We refer to [START_REF] Anastassiou | Weak convergence and the Prokhorov radius[END_REF], [START_REF]Sets of superresolution and the maximum entropy method on the mean[END_REF] and [START_REF] Doukhan | Superresolution rates in Prokhorov metric[END_REF] where quantitative evaluations on the size of the set of solutions is performed and to [START_REF] Gamboa | The maximum entropy method on the mean: applications to linear programming and superresolution[END_REF] and [START_REF] Lewis | Superresolution in the Markov moment problem[END_REF] for related evaluations in another context. Notice that the structure of points b can be completely described in the case where the family of functions (a i ) is a Chebyshev system, T -system for short, this includes the case of discrete Fourier transform and moments, see [START_REF] Borwein | Polynomials and polynomial inequalities[END_REF] or [START_REF] Karlin | Tchebycheff systems: With applications in analysis and statistics[END_REF] for an exhaustive overview on these systems of functions. A more involved situation is when (1.2) does not enjoy the non negativity assumption on the measure. This means that one wishes to solve the linear equation:

(1.3) a i , µ = X a i (x)µ(dx) = b i where µ ∈ R(X) and i = 1, . . . , m .
Here, R(X) is the set of signed Radon measures on X. This is the frame of the present paper. Surprisingly, as shown by authors of this paper [START_REF] Castro | Exact reconstruction using beurling minimal extrapolation[END_REF] and under some assumptions on the family of functions (a i ), there exists pairs (b The study of the solutions to (1.1) and (1.3) uncovers that 1 -minimization faithfully reconstruct objects concentrated on very few points. However, the analysis in Super Resolution differs dramatically from Compressed Sensing. For instance, it is well known that sparse 1 -minimization cannot be successful in the ultrahighdimensional setting [START_REF] Verzelen | Minimax risks for sparse regressions: Ultra-high dimensional phenomenons[END_REF] where N exp(m). Indeed it was shown in [START_REF] Chafaï | Interaction between compressed sensing, random matrices and high dimensional geometry[END_REF] that one needs at least m ≥ (cst)s log(N/s) measures to faithfully uncover s sparse vectors by 1 -minimization. Hence, the analysis of Compressed Sensing in terms of high-dimensional random geometry [START_REF] Chafaï | Interaction between compressed sensing, random matrices and high dimensional geometry[END_REF] cannot be extended to the space of measure. Moreover, observe that Compressed Sensing aims at recovering a sparse signal from random projection while, in Super Resolution, the sampling scheme is deterministic.

, µ b ) ∈ R m × R(X) such that b i = a i , µ b , for i = 1, • • • , m,
Admittedly their analysis differ but we can bridge the gap between (1.1) and (1.3) by considering their dual formulations. From the point of view of convex analysis, we see that the dual form of these linear programs aims at reconstructing a dual certificate [START_REF] Candes | Near-optimal signal recovery from random projections: universal encoding strategies?[END_REF][START_REF] Castro | Exact reconstruction using beurling minimal extrapolation[END_REF], i.e. an ∞ -constraint linear combination of (a i ), where (a i ) are the lines of A in (1.1) and a family of continuous functions in (1.3). As pointed out by authors of this paper [START_REF] Castro | Exact reconstruction using beurling minimal extrapolation[END_REF], a parallel between Compressed Sensing and Super Resolution exists where the lines of A are the evaluation of a vector of continuous function at some prescribed points. In this frame, Super Resolution can be seen as a Compressed Sensing problem where the dimension N goes to infinity. This analogy persists with the notion of dual certificate, i.e. a solution to the following dual program (1.4). Indeed, the dual programs given by the constraints (1.1) and (1.3) (while minimizing the 1 -norm and the total variation norm respectively) share the same expression:

(1.4) sup u∈R m b u s.t. a u ∞ ≤ 1 ,
where a = A in Compressed Sensing (1.1) and a(x) = (a i (x)) m i=1 is a vector of continuous function in Super Resolution (1.3). Define a dual certificate P as:

(1.5)

P = a u ,
where u is a solution to the dual program (1.4). From the duality properties, note that P is a sub-gradient of the 1 -norm at a solution to the primal program (1.3). Hence, we can ensure that a target measure µ is a solution to (1.3) if we are able to construct a dual polynomial (1.5) that interpolates the phases of the weights of µ at its support points. From a theoretical point of view, one of the main issue in Super Resolution consists in exhibiting such a dual certificate P. In the Fourier frame, notice that an important construction, for target discrete measures whose support satisfy a separable condition, is given in the fundamental paper [START_REF] Candes | Towards a mathematical theory of superresolution[END_REF] where a huge step has been taken. Indeed, the authors are the first to give a sharp condition on the support points of the target measure in order to warrant the existence of a dual certificate. Moreover, their proof is based on interpolating, by a Jackson kernel, the phases of the weights of the target measure at its support points and, henceforth, explicitly construct a dual certificate (1.5). In the Compressed Sensing frame, observe that the same method has been investigated in [START_REF] Kahane | Idempotents et échantillonnage parcimonieux[END_REF] using a Dirichlet kernel. In the present paper, we will not deal with this issue but rather with the practical resolution of the convex program (1.4).

In the Super Resolution frame, remark that the program (1.4) has finitely many variables but infinitely many constraints. This last point can be a severe limitation in practice. As a matter of fact, a difficult task is to construct a tractable program that deals with the ∞ -constraint of (1.4). Standard formulations [START_REF] Candes | Towards a mathematical theory of superresolution[END_REF] are based on Gram matrix representations, see below. Incidentally, these procedures cannot be extended to dimension greater than 2 or to semi-algebraic domains X. To cope with this issue, we consider a new parametrization of the primal program based on works of the authors [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]. Note that our method relies on infinitely many parameters but relaxations involving only a finite number of parameters are proved to lead to the exact solution of the primal program. 1.2. Previous works. During the last years, theoretical guarantees for exact recovery [START_REF] Bredies | Inverse problems in spaces of measures[END_REF][START_REF] Castro | Exact reconstruction using beurling minimal extrapolation[END_REF][START_REF] Candes | Towards a mathematical theory of superresolution[END_REF], bounds on the support recovery from inaccurate samplings [START_REF] Azais | Spike detection from inaccurate samplings[END_REF][START_REF] Fernandez-Granda | Support detection in super-resolution[END_REF], prediction of the Fourrier coefficient from noisy observations [START_REF] Tang | Near minimax line spectral estimation[END_REF], and noise robustness [START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF] have been showed. These works prove that discrete measures can be recovered, in a robust manner, from few samples using an 1 -method.

From a numerical point of view, a solution to 1 -minimization is often computed using the dual program described by (1.4). Then, the ∞ -norm constraint of the dual program (1.4) is equivalently formulated as a nonnegative constraint on (trigonometric) polynomials. This point of view unleashes Gram matrix representations [START_REF] Candes | Towards a mathematical theory of superresolution[END_REF] or Toeplitz matrix representations [START_REF] Tang | Near minimax line spectral estimation[END_REF] to handle the constraint of non-negativity of (trigonometric) polynomials on domains. However, these formulations are limited to the frame of the real line and the torus in dimension one (see [START_REF] Dumitrescu | Positive trigonometric polynomials and signal processing applications[END_REF] for instance) since they rely on the Fejér-Riesz theorem. As a matter of fact, the literature of Super Resolution has been focused on the fact that a semidefinite programming (SDP) formulation for the dual problem can be given as long as there exists a spectral factorization for a globally nonnegative trigonometric polynomial [START_REF] Candes | Towards a mathematical theory of superresolution[END_REF] (Bounded Real Lemma using Fejér-Riesz theorem) or a spectral decomposition for semi-definite Toeplitz matrices [START_REF] Tang | Near minimax line spectral estimation[END_REF] (Caratheodory-Toeplitz theorem). Hence, except on the real line and the torus in dimension one, there is no exact SDP formulation of the Super Resolution problem for the dual form. However, relaxed SDP versions of the dual form in dimension greater than 2 are discussed in [XCM + 14] and they have been used on the 2-sphere in [START_REF] Bendory | Exact recovery of dirac ensembles from the projection onto spaces of spherical harmonics[END_REF][START_REF]Super-resolution on the sphere using convex optimization[END_REF].

Contribution.

To the best of our knowledge, the present paper is the first to overcome this limitation and expand the scope of Super Resolution implementation to the multi-dimensional frame in general basic semi-algebraic domains. Indeed, we operate a smart method to tackle numerically the solution of equation (1.3) with minimal 1 -norm. This method uses both a re-parametrization in terms of moment sequences [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF] and the so-called sum-of-squares (SOS) decompositions of nonnegative multidimensional polynomials, used widely in systems control theory during the last decade, see for example [START_REF] Henrion | Positive polynomials in control[END_REF]. Notice that, in the scope of Super Resolution, this technique is new and, contrary to other approaches, focuses on the primal program through a truncation of the moment sequences.

More precisely, given the real numbers b i , i = 1, . . . , m, consider the infinitedimensional optimization problem:

(1.6) inf µ T V s.t. a i , µ = b i , i = 1, . . . , m µ ∈ R(X),
where . T V is the total variation norm of measures (to be defined later). Notice that, under standard assumptions, problem (1.6) is feasible. That is

(1.7) ∃µ ∈ R(X) such that for i = 1, . . . , m , a i , µ = b i .
Our main contribution concerns the numerical resolution of the total variation minimization problem (1.6). We extend the univariate (n = 1) trigonometric SDP formulation of [START_REF] Candes | Towards a mathematical theory of superresolution[END_REF] to a much more general SDP formulation in dimension n ≥ 2, for measures supported on basic semialgebraic sets.

To this end, we use the Jordan decomposition of the signed measure µ = µ + -µ - as a difference of two nonnegative measures supported on X and we follow [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF] to define a hierarchy of finite-dimensional primal-dual SDP problems:

• the primal problems correspond to SDP relaxations of the conditions that must satisfy finitely many moments of the two nonnegative measures on X; • the dual problems correspond to SDP strengthenings using SOS multipliers of the conditions that two distinguished polynomials are nonnegative on X. The moment-SOS hierarchy is indexed by an integer k, called relaxation order, which is the (half of the) number of moments used to represent the measures in the primal problem, or equivalently, the (half of the) degree of the SOS representations of the polynomials in the dual problem. The larger is the relaxation order k, the larger is the size of the SDP problems, the number of variables and constraints growing polynomially in O(k n ).

The primal SDP problem features the matrices of moments of the two nonnegative measures. If the rank of each moment matrix, as a function of k, stabilizes to a certain constant value, then the corresponding measure is atomic, with the number of atoms equal to the rank. Therefore, the total variation minimization problem (1.6) has been solved succesfully, and this is certified by the polynomials solving the SDP problem. Numerical linear algebra can then be used to retrieve the support of the optimal measure.

In the sequel, we present some examples for which our method is the first to give an SDP formulation of the Super Resolution phenomena. As a matter of fact, our procedure encompasses a larger class of measurements than the class of standard moments discussed previously. Our numerical experiments are carried out with the Matlab interface GloptiPoly 3 which is designed to generate semidefinite relaxations of measure LP problems with polynomial data. So we assume that the functions a i (x) in LP problem (2.1) are multivariate polynomials, and for notational simplicity, we let a i (x) :=

x αi = x α 1 i 1 • • • x α n i
n where α i ∈ N n are given. Note that the choice of monomials is only motivated for notational simplicity, and that other choices of polynomials (e.g. Chebyshev polynomials) are typically preferable numerically1 . SDP relaxations are then solved with SeDuMi or MOSEK, implementations of a primal-dual interior-point algorithm. For reproducibility purposes, our Matlab codes (using the public-domain interface GloptiPoly and the SDP solver SeDuMi) of the numerical examples presented next are available for download at homepages.laas.fr/henrion/software/tvsdp.tar.gz 1.3.1. Disconnected domain. We want to recover the measure:

µ := δ -3/4 + δ 1/2 -δ 1/8
on the disconnected set X := [-1, -1/2] ∪ [0, 1] which can be modeled as the polynomial superlevel set X = {x ∈ R : g 1 (x) ≥ 0} for the choice:

g 1 (x) := -(x + 1)(x + 1/2)x(x -1).
In LP (2.1) we let a i (x) := x i and b i := (-3/4) i +(1/2) i -(1/8) i for i = 0, 1, 2 . . . , 9.

Solving the SDP relaxation of order k = 5 on our standard PC takes 0.2 seconds, and optimality is certified from the solution of the primal moment problem with a rank 2 moment matrix for µ + and a rank 1 moment matrix for µ -, from which the 3 points can be extracted using numerical linear algebra. On Figure 1 we represent the degree 9 polynomial 9 i=0 u i x i certifying optimality, constructed from the solution of the dual SOS problem. Indeed we can check that the polynomial attains the value +1 at the points x = -3/4, and x = 1/2 (in red), it attains the value -1 at the point x = 1/4 (in blue), while taking values between -1 and +1 on X. Notice in particular that the polynomial is larger than +1 around x = -1/4, but this point is not in X. 1.3.2. Low-pass filters in dimension greater than 3. In the Fourier frame, the recent SDP formulations of 1 -minimization in the space of complex valued measures are based on the Fejér-Riesz theorem. As a consequence, they cannot handle dimensions greater than 3. Observe that our procedure can bypass this limitation. For sake of readability, we present an example in dimension 2 although it can be extended to any dimension. We want to recover the measure

µ := δ (-1/2,1/2) + δ (1/2,-1/2) + δ (1/2,1/2) + δ (0,0) -δ (0,-1/2) -δ (1/2,0)
on the box X := [-1, 1] 2 , from the knowledge of moments of degree up to 12, i.e. a i (x) = x i for i = 0, 1, . . . , 12. Solving the SDP relaxation of order k = 6 on our standard PC takes less than 3 seconds, and optimality is certified from the solution of the primal moment problem with a rank 4 moment matrix for µ + for and a rank 2 moment matrix for µ -from which the 6 points of the support of the optimal measure µ can be extracted using numerical linear algebra with a relative accuracy around 10 -6 . On Figure 2 we represent the degree 12 polynomial certifying optimality, constructed from the solution of the dual SOS problem. Indeed we can check that the polynomial attains the value +1 at the 3 points x ∈ {(-1/2, 1/2), (1/2, -1/2), (0, 0)}, it attains the value -1 at the 2 points x ∈ {(0, -1/2), (1/2, 0)} (in blue), while taking values between -1 and +1 on X.

Figure 3. Degree 6 polynomial certificate for the sphere example, with 3 points (red) in the support of the positive part, and 3 points (blue) in the support of the negative part of the optimal measure. 1.3.3. Localization of points on the sphere. Recent extensions of Super Resolution to spike deconvolution on the 2-sphere from spherical harmonic measurements has been investigated in [START_REF] Bendory | Exact recovery of dirac ensembles from the projection onto spaces of spherical harmonics[END_REF][START_REF]Super-resolution on the sphere using convex optimization[END_REF]. In these paper, the authors give a sufficient condition for exact recovery using 1 -minimisation and they investigate spikes localization when the measurements are perturbed by additive noise.

From a numerical point of view, they used a relaxed version of the dual program (bounded real lemma in dimension d = 3 and a Gram representation of the ∞constraint appearing in the dual). Our work naturally extends to this frame and provides an exact formulation of the primal form.

For sake of numerical code simplicity, we have considered polynomials on R 3 restricted to the domain X given by the 2-sphere (note that one could have used homogenous spherical harmonics instead as in [START_REF] Bendory | Exact recovery of dirac ensembles from the projection onto spaces of spherical harmonics[END_REF]). We want to recover the measure

µ := δ (1,0,0) + δ (0,1,0) + δ (0,0,1) -δ ( √ 2 2 , √ 2 2 ,0) -δ ( √ 2 2 ,0, √ 2 2 ) -δ (0, √ 2 2 , √ 2
2 )

that is supported on the positive orthant just for better visualization purposes. In LP (2.1), the a i consist of 3-variate monomials of degree up to 5, i.e. m = 56. Solving the SDP relaxation of order k =6 on our standard PC takes less than 20 seconds, and optimality is certified from the solution of the primal moment problem with a rank 3 moment matrix for µ + and a rank 3 moment matrix for µ -, from which the 6 points can be extracted using numerical linear algebra. On Figure 3 we represent the degree 6 polynomial certifying optimality on the 2-sphere, constructed from the solution of the dual SOS problem. Indeed we can check that the polynomial attains the value +1 at the 3 prescribed points, it attains the value -1 (in red), at the 3 others prescribed points (in blue), while taking values between -1 and +1 on X.

From a theoretical point of view, the minimal separation condition appearing in [START_REF] Bendory | Exact recovery of dirac ensembles from the projection onto spaces of spherical harmonics[END_REF] requires a degree 2 polynomial. So our example satisfy the sufficient aforementioned condition.

Primal and dual LP formulation

General model and notation. Let n be a positive integer. Denote by R[x]

the set of all polynomials on R n , and for d ∈ N, R d [x] the set of all polynomials on R n with degree not greater than d. Further, we use the following notation:

• X ⊂ R n , is a given closed basic semi-algebraic set:

X := {x ∈ R n : g j (x) ≥ 0, j = 1, . . . , n X }
where g j ∈ R[x], j = 1, . . . , n X , are given polynomials whose degrees are denoted by d j , j = 1, . . . , n X . It is assumed that X is compact with an algebraic certificate of compactness. For example, one of the polynomial inequalities g j (x) ≥ 0 should be of the form:

R 2 - n i=1
x 2 i ≥ 0 , for R a sufficiently large constant. Let g X := (g j ) j=1,...,n X . • Let a = (a i ) m i=1 be a linearly independent family of polynomials of degree at most d on X. Notice that m ≤ (1 + d) n .

• For monomials we use the multi-index notation

x α := n j=1 x αj j for every x = (x 1 , . . . , x n ) ∈ R n and α = (α 1 , . . . , α n ) ∈ N n .
• C (X), the space of continuous functions on X, a Banach space when equipped with the sup-norm:

f = sup x∈X |f (x)| .
• R(X), the space of signed Radon measures on X, a Banach space isometrically isomorphic to the topological dual C (X) * when equipped with the total variation norm:

µ T V = sup P E∈P |µ(E)| ,
where the supremum is taken over all partitions P of X into a finite number of disjoint measurable subsets. • C (X) + ⊂ C (X) and R(X) + ⊂ R(X) the respective positive cones of nonnegative continuous functions on X and (nonnegative) Radon measures on X. We use the standard notation f ≥ 0 and µ ≥ 0 for membership in C (X) + and R(X) + , respectively. • To denote the integration of a function against a measure, we use the duality bracket:

f, µ = X f dµ ,
for all f ∈ C (X), and µ ∈ R(X). With the usual Jordan decomposition:

µ = µ + -µ -
into a sum of two nonnegative Borel measures µ + , µ -, the optimization problem (1.6) can be rewritten equivalently as a linear programming (LP) problem in the convex cone R(X) + , namely:

(2.1)

p * = inf 1, µ + + 1, µ - s.t. a i , µ + -a i , µ -= b i , i = 1, . . . , m µ + ∈ R(X) + µ -∈ R(X) + . If a = (a i ) i=1,...,m ∈ C (X) m and b = (b i ) i=1,...,m ∈ R m , problem (2.1
) is the dual of the following LP problem in the convex cone C (X) + :

(2.2)

d * = sup b u s.t. z + (x) := 1 + a (x)u ∈ C (X) + z -(x) := 1 -a (x)u ∈ C (X) +
where the maximization is w.r.t. u = (u i ) i=1,...,m ∈ R m . Remark that LP problem (2.2) can be also written as:

d * = sup b u s.t. a (x)u ∞ ≤ 1.
Lemma 1. There is no duality gap between primal LP (2.1) and dual LP (2.2), i.e. p * = d * .

Proof: Define the vector r(µ + , µ -) ∈ R m+1 by:

r(µ + , µ -) := ( 1, µ + + 1, µ -, a 1 , µ + -a 1 , µ -, . . . a m , µ + -a m , µ -)
and the set

R := {r(µ + , µ -) : (µ + , µ -) ∈ R(X) + × R(X) + } ⊂ R m+1 .
By [Bar02, Theorem 7.2], p * = d * provided that p * is finite and R is closed. Finiteness of p * follows from Assumption 1.7 and nonnegativity of the objective function 1, µ + + 1, µ -. To prove that R is closed we have to show that for any sequence

(µ n + , µ n -) n∈N ∈ R(X) + × R(X) + such that r(µ n + , µ n -) → s ∈ R m+1
as n → ∞, one has s = r(µ, ν) for some finite measures µ, ν ∈ R(X) + . Since the supports of all the measures are contained in a compact set, and since 1, µ n + + 1, µ n -→ s 0 all measures µ n + , µ n -are uniformly bounded. Therefore, from the weak-* compactness (and weak-* sequential compactness) of the unit ball (Banach-Alaoglu's Theorem), there is a subsequence (µ n k + , µ n k -) k∈N that converges weakly-* to an element (µ, ν) ∈ R(X) + × R(X) + . In particular, as all a i are continuous,

lim k→∞ r(µ n k + , µ n k -) = r(µ, ν),
which proves that R is closed.

Lemma 2. For the dual LP problem (2.2) the supremum is attained.

Proof:

The feasibility set

U := {u ∈ R m : a (x)u ∞ ≤ 1}
of the LP problem (2.2) is a closed convex subset of a finite-dimensional Euclidean space, and it contains the origin. Since the objective function in LP (2.2) is continuous on U , the optimum is attained if U is bounded. Suppose that U is not bounded. Then there exists a sequence (u n

) n∈N ⊂ U such that u n → ∞ as n → ∞. Write u n = λ n v n , with v n = 1. Notice that 0 < λ n → ∞ and v n ∈ U because 0 ∈ U and U is convex. Then a (x)u n ∞ = a (x)λ n v n ∞ = = λ n a (x)v n ∞ ≤ 1, so that a (x)v n ∞ ≤ λ -1 n → 0 as n → ∞. Since v n = 1, there exists a subsequence n k and v with v = 1 such that v n k → v as k → ∞ and a (x)v ∞ = 0.
By linear independence of a, this implies that v = 0, a contradiction.

As a consequence of strong duality of Lemma 1, for any optimal primal-dual pair (µ, u) we have the complementarity conditions

z + , µ + = 0, z -, µ -= 0 implying jointly with Lemma 2 that spt µ + ⊂ { x ∈ X : a (x) u = 1 } and spt µ -⊂ { x ∈ X : a (x) u = -1 } for some continuous function x → a (x)u,
and where spt µ denotes the support of µ, that is, the smallest closed set S ⊂ R n such that µ(R n \ S) = 0. Lemma 3. Problem (1.6) has an optimal atomic measure supported on at most 2(m + 1) points.

Proof: Let µ + be a nonnegative measure solving problem (2.1) and let b +i := a i , µ + , with a 0 := 1, i = 0, 1, . . . , m. If b + = 0 then µ + = 0 is trivially atomic (with no atoms), so assume b +0 = 0, and consider the probability measure μ+ := µ + /b +0 which satisfies the m equality constraints a i , μ+ = b+i := b +i /b +0 , i = 1, . . . , m. From [Bar02, Proposition 9.4] there exists a probability measure μ+ satisfying the same equality constraints a i , μ+ = b+i and which is supported on (at most) m+1 points of X. The same reasoning can be applied to any nonnegative measure µ -solving problem (2.1), which has a discrete counterpart μ-supported on (at most) m + 1 points of X. The result follows by considering the union of these two discrete supports, which consists of (at most) 2(m + 1) points of X.

Primal and dual SDP formulation

Problem (2.1) is an instance of a generalized moment problem. As such it can be solved by a converging hierarchy of finite-dimensional primal-dual semidefinite programming (SDP) problems, as described comprehensively in [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]. In the sequel, we extract the key instrumental ingredients to the construction of the hierarchy.

3.1. Primal moment SDP. Recall from paragraph 2.1 that

X := {x ∈ R n : g j (x) ≥ 0, j = 1, . . . , n X }
is a basic semi-algebraic set with a straighforward certificate of compactness, and let g X := (g j ) j=1,...,n X denote its defining polynomials. Given a measure µ ∈ R + (X), the real number (3.1)

y α := x α , µ
is called its moment of order α ∈ N n . Conversely, given a real valued sequence y := (y α ) α∈N n , if identity (3.1) holds for all α ∈ N n , we say that y has a representing measure µ ∈ R + (X). Equivalently, sequence y belongs to the infinite-dimensional moment cone

M (X) := {(y α ) α∈N n : y α = x α , µ , µ ∈ R(X) + }.
In the sequel we describe a procedure to approximate this convex cone. Given k ∈ N, let R[x] k denote the space of real polynomials of degree at most k. Let us identify a polynomial p(x) = α p α x α ∈ R[x] k with its vector p of coefficients in the monomial basis. Define the Riesz functional y as the linear functional acting on polynomials as follows:

p ∈ R[x] k → y (p) = α p α y α = p y ∈ R.
Note that if sequence y has a representing measure µ, then y (p) = p, µ . Define the moment matrix of order k as the Gram matrix of the quadratic form p ∈ R[x] k → y (p 2 ) ∈ R, i.e. the matrix M k (y) such that y (p 2 ) = p M k (y)p. By construction this matrix is symmetric and linear in y. Given a polynomial g ∈ R[x], define its localizing matrix of order k as the Gram matrix of the quadratic form p ∈ R[x] k → y (gp 2 ) ∈ R, i.e. the matrix M k (g y) such that y (gp 2 ) = p M k (g y)p. By construction this matrix is symmetric and linear in y. For j = 1, . . . , n X , let k j denote the smallest integer not less than half the degree of polynomial g j , and let k X := max{1, k 1 , . . . , k n X }. With these notations, and for k ≥ k X , define the finite-dimensional moment cone

M k (g X ) := {(y α ) |α|≤2k : M k (y) 0, M k-kj (g j y) 0, j = 1, . . . , n X }
where 0 means positive semidefinite.

Let y + resp. y -denote the sequence of moments

y +α := x α µ + (dx), y -α := x α µ - (dx)
of µ + (resp. µ -), indexed by α ∈ N n . Primal measure LP (2.1) can be written as a primal moment LP:

p * = min y +0 + y -0 s.t. A(y + , y -) = b y + ∈ M (X) y -∈ M (X)
where the linear system of equations A(y , y -) = b models the linear moment constraints. The moment relaxation of order k ≥ max{k X , d} of the primal moment LP then reads:

(3.2)

p * k = min y +0 + y -0 s.t. A(y + , y -) = b y + ∈ M k (g X ) y -∈ M k (g X )
where the minimization is w.r. ) is an optimal solution of (2.1) and p k = p * . Given a moment matrix M k (y * + ) satisfying the rank constraint of Theorem 1, there is a numerical linear algebra algorithm that extracts the r + points of the support of the corresponding atomic measure µ * + , and similarly for µ * -. The algorithm is described e.g. in [Las10, Section 4.3] and it is implemented in the Matlab toolbox GloptiPoly 3.

A certificate of optimality can be obtained by solving the dual problem to primal SDP problem (3.2), and this is described next.

3.2. Dual SOS SDP. For a given integer k, let Σ[x] k ⊂ R[x] 2k denote the space of SOS (sums of squares) polynomials of degree at most 2k. If p ∈ Σ[x] k this means that there exists q j ∈ R[x] k , j = 1, . . . , n p , such that p = np j=1 q 2 j . Let P(X) := {p ∈ R[x] : p(x) ≥ 0, ∀x ∈ X} denote the infinite-dimensional cone of nonnegative polynomials on X, and for k ≥ k X , define the finite-dimensional SOS cone, also called quadratic module

P k (g X ) := {p 0 + n X j=1 g j p j , p j ∈ Σ[x] k-kj , j = 0, 1, . . . , n X } ⊂ R[x] 2k .
Under the above assumptions on the polynomial family g X defining X, from Putinar's theorem, see e.g. [Las10, Theorems 2.14 and 3.8], it holds that P(X) ∩ R[x] κ is the closure of (∪ k≥k X P k (g X )) ∩ R[x] κ for all κ ≥ k X . Observe also that M (X) (resp. M k (g X )) is the dual cone to P(X) (resp. P k (g X )). Whereas testing whether a given polynomial belongs to P(X) is a difficult task, testing whether a given polynomial belongs to P k (g X ), for a fixed k, amounts to solving an SDP problem.

Dual continuous function LP (2.2) can be written as a positive polynomial LP

d * = max b u s.t. 1 + a (x) u ∈ P(X) 1 -a (x) u ∈ P(X)
and its SOS strengthening of order k ≥ k X reads:

(3.4)

d * k = max b u s.t. 1 + a (x) u ∈ P k (g X ) 1 -a (x) u ∈ P k (g X )
where the maximization is w.r.t. a vector u ∈ R m . It turns out that this is SOS problem (3.4) is an SDP problem dual to the moment problem (3.2):

Lemma 4. There is no duality gap between SDP problems (3.2) and (3.4), i.e. p * k = d * k , and both (3.2) and (3.4) have an optimal solution. Proof: We first show that (3.2) has an optimal solution. Recall that one of constraints g j (x) ≥ 0 that define X states that M -x 2 ≥ 0 for some M > 1.

From the constraint M k-kj (g j y + ) 0 one deduces that y+ (M -x2 i ) ≥ 0, and y+ (M x t i -x t+2 i ) ≥ 0 for every t = 1, . . . , 2k -2. Hence y + (x 2k i ) ≤ M k y +0 for every i = 1, . . . , n. With similar arguments, y -(x 2k i ) ≤ M k y -0 for every i = 1, . . . , n. By [Las10, Proposition 3.6] |y +α | ≤ M k y +0 and |y -α | ≤ M k y -0 for all α ∈ N n 2k . Next, differential equation of order n, it was proved in [START_REF] Neustadt | Optimization, a moment problem, and nonlinear programming[END_REF] that there is always an n-atomic optimal measure solving problem (2.1).

In practice, Theorem 1 should be used as follows:

1 Let k = max{d, k X }.

2 Solve SDP problem (3.2) and its dual (3.4) with a primal-dual algorithm.

3 If the rank condition (3.3) of Theorem 1 is satisfied, then extract the measure from the solution of (3.2) and the polynomial certificate from the solution of (3.4). Otherwise, let k = k + 1, and go to 1.

We conjecture that if the data a, b in problem (2.1) are generic, then there is a finite value of k for which the rank condition of Theorem 1 is satisfied. The rationale behind this assertion follows from a result by Nie [START_REF]Optimality conditions and finite convergence of lasserre's hierarchy[END_REF] on generic finite convergence for the moment-SOS SDP hierarchy for polynomial optimization over compact basic semi-algebraic sets. Translated in the present context for a fixed family of data a, results in [START_REF]Optimality conditions and finite convergence of lasserre's hierarchy[END_REF] yield that there is a set of polynomials {h 1 , . . . , h L } ⊂ R[u], such that, given a feasible solution u of (2.2), if h (u) = 0 for all = 1, . . . , L, then indeed

1 + a T (x)u = p 1 0 (x) + n X j=1 p 1 j (x) g j (x), x ∈ R n , and
1 -a T (x)u = p 2 0 (x) + n X j=1 p 2 j (x) g j (x), x ∈ R n ,
for some SOS polynomials p k j , k = 1, 2 and j = 1, . . . , n X . So if the optimal solution u * of (2.2) satisfies h (u * ) = 0, = 1, . . . , L, then d * = d * k for some index k (i.e. finite convergence takes place). Similarly, by [START_REF] Nie | Certifying convergence of lasserre's hierarchy via flat truncation[END_REF] the rank-condition (3.3) of Theorem 1 also holds generically for polynomial optimization (which however is a context different from the present context). Put differently, finite convergence would not hold only if every optimal solution u of (2.2) would be a zero of some polynomial of the family {h 1 , . . . , h L } ⊂ R[u]. But so far we have not proved that at least one optimal solution u * of (2.2) is not a zero of some of the polynomials h , at least for generic b.

Of course, finite convergence occurs for trigonometric polynomials on X = [0, 2π], which follows from the Fejér-Riesz theorem and this was exploited in the landmark paper [START_REF] Candes | Towards a mathematical theory of superresolution[END_REF]. Similarly, but apparently not so well-known, the Fejér-Riesz theorem also holds in dimension n = 2. Indeed it follows from Corollary 3.4 in [START_REF] Scheiderer | Sums of squares on real algebraic surfaces[END_REF] that every non-negative bivariate trigonometric polynomial can be written as a sum of squares of trigonometric polynomials 2 . So again for trigonometric polynomials on X = [0, 2π] 2 , finite convergence of the hierarchy (3.4) takes place, i.e., d * k = d * . Note however that in contrast to the one-dimensional case, there is no explicit upper bound on the degrees of the sum of squares which are required, so that even in the two-dimensional Fourier case we do not have an a priori estimates on the smallest value of k for which d * k = d * an for which we can guarantee that the rank condition of Theorem 1 is satisfied.

Generally speaking, even if our genericity conjecture is true, we do not have a priori estimates on the smallest value of k for which Theorem 1 holds. As mentioned above, this also true even in the two-dimensional case on [0, 2π] 2 where finite convergence is guaranteed in all cases.

Appendix

We first recall some standard results of convex analysis. x α A jα , j = 1, . . . , n X , for some appropriate real symmetric matrices A jα .

Consider a sequence (q t ) t∈N ⊂ P k (g X ) such that q t → q ∈ R[x] 2k as t → ∞. That is

q t (x) = p 0t (x) + n X j=1 p jt (x) g j (x), ∀x ∈ R n ,
for some p jt ∈ Σ[x] k-vj , j = 0, . . . , n X , for all t ∈ N. More precisely, coefficient-wise (5.1)

q tα = Q 0t , A 0α + n X j=1 Q jt , A jα , ∀ α ∈ N n 2k ,
for some appropriate matrices Q jt ∈ S k-vj +

. Let y = (y α ) α∈N n 2k be the moments y α := X x α dx of the measure uniformly supported on X. Observe that since X has nonempty interior, Put differently M k (y) 0 and M k-vj (g j y) 0, j = 1, . . . , n X . The convergence q t → q implies q t , y → q, y as t → ∞. Hence there is some η such that η ≥ q t , y for all t ∈ N. This in turn implies η ≥ q t , y = p 0t , y + ) * , j = 1, . . . , n X , one may invoke Lemma 7 and conclude that the sequences (Q 0t ) t∈N ⊂ S k + and (Q jt ) t∈N ⊂ S k-vj + are norm-bounded. Therefore there is a subsequence (t ) ∈N and matrices Q 0 ∈ S k + and Q j ∈ S k-vj + , j = 1, . . . , n X , such that

Q 0t → Q 0 , Q jt → Q j , j = 1, . . . , n X
as → ∞. Taking the limit for the subsequences (q t α ) ∈N and (Q jt ) ∈N in (5.1) yields coefficient-wise

q α = Q 0 , A 0α + n X j=1 Q j , A jα , ∀ α ∈ N n 2k ,
which proves that q ∈ P k (g X ), the desired result.

(1. 1 )

 1 Ax = b where x ∈ R N , b := Ax , b ∈ R m and m N .

  and µ b is the unique solution of (1.3) minimizing the total variation norm with b = b . Such µ b are sparse in the sense that they are measures with finite support.

Figure 1 .

 1 Figure 1. Degree 9 polynomial certificate for the univariate example, with 2 points (red) in the support of the positive part, and 1 point (blue) in the support of the negative part of the optimal measure.

Figure 2 .

 2 Figure 2. Degree 12 polynomial certificate for the bivariate example, with 4 points (red) in the support of the positive part, and 2 points (blue) in the support of the negative part of the optimal measure.

  Lemma 6. ([FK94, Corollary I.1.3]) Let C ⊂ R n be a closed convex cone with dual C * = {y : x, y ≥ 0, ∀x ∈ C}. Then int C * = ∅ ⇔ C ∩ (-C) = {0}. Lemma 7. ([FK94, Corollary I.1.6]) Let C ⊂ R n be a closed convex cone whose dual C * has nonempty interior. Then for all y ∈ int C * , the set {x ∈ C : x, y ≤ 1} is compact. Lemma 8. The convex cone P k (g X ) is closed. Proof: Let S n + be the convex cone of real symmetric matrices of size n that are positive semidefinite. Let N n k be the set of n-dimensional integer vectors α such that n i=1 α i ≤ k and let v k (x) := (x α ) α∈N n k be a vector of monomials of degree up to k. Next let v k (x) v k (x) T = α∈N n 2k x α A 0α and v k-vj (x) v k-vj (x) T g j (x) = α∈N n 2k

X

  p(x) dx > 0 ∀0 = p ∈ Σ[x] k , and X p(x) g j (x) dx > 0 ∀0 = p ∈ Σ[x] k-vj , j = 1, . . . , n X .

p

  jt g j , y= Q 0t , M k (y) + n X j=1 Q jt , M k-vj (g j y) . (5.2) Therefore sup t Q 0t , M k (y) ≤ η, sup t Q jt , M k-vj (g j y) ≤ η, j = 1, . . . , n X .As 0 ≺ M k (y) ∈ int(S k + ) * , and 0 ≺ M k-vj (g j y) ∈ int(S k-vj +

  t. a vector (y + , y -) of moments of degree at most 2k. For fixed k, problem (3.2) is a finite-dimensional linear programming problem in the convex cone of positive semidefinite matrices, i.e. an SDP problem. When k varies, the number of moments, as well as the size of the moment and localizing matrices in problem (3.2) are binomial coefficients growing in O(k n ).Theorem 1. For a given relaxation order k ≥ max{d, k X }, let (y * := rank M k (y * + ), (resp. r -:= rank M k (y * -)) points.

	then p * k = p * and LP (2.1) has an optimal solution (µ * + , µ * -) with µ * + (resp. µ * -)
	atomic supported at r Proof: By [Las10, Theorem 3.11], y * + (resp. y * -) is the vector of moments up to order 2k, of a measure µ * + (resp. µ * -) supported on rank M k (y * + ) points (resp. rank M k (y * -)) points of X. Therefore (µ * + , µ * -) is a feasible solution of (2.1) with value p k ≤ p * , which proves that (µ * + , µ * -
	It can be shown that (p * k ) is a monotonically nondecreasing converging sequence of lower bounds on p * , i.e. p * k+1 ≥ p * k and lim k→∞ p * k = p * . However, in the context
	of solving LP (2.1), a more relevant result is the following:
	+ , y * -) denote the
	solution of the moment SDP (3.2). If
	(3.3) rank M k-k X (y * + ) = rank M k (y * + ) and rank M k-k X (y * -) = rank M k (y * -)

+

A numerical analysis of the impact of the basis is however out of the scope of our work.

We are grateful to Markus Schweighofer for providing this reference.

in a minimizing sequence (y s + , y s -), s ∈ N, of (3.2) one has y s +0 +y s -0 ≤ y 1 +0 +y 1 -0 =: ρ for all s, and so |y s +α | ≤ M k ρ and |y s -α | ≤ M k ρ for all α ∈ N n 2k , and all s = 1, . . .. From this we deduce that there is a subsequence (y st + , y st -), t ∈ N, that converges to some (y * + , y * -) as t → ∞, with value y * +0 + y * -0 = p * k . In addition by a simple continuity argument, M k (y * + ) 0 and M k-kj (g j y * + ) 0, j = 1, . . . , n X . Similarly M k (y * -) 0 and M k-kj (g j y * -) 0, j = 1, . . . , n X , which proves that (y * + , y * -) is an optimal solution of (3.2).

Next, the set of optimal solutions y * := {(y * + , y * -)} of (3.2) is compact. This follows from

And so every sequence in y * has a converging subsequence. From [Bar02, Chapter IV. Theorem 7.2] one also deduces that there is no duality gap between (3.2) and (3.4).

It remains to prove that (3.4) has an optimal solution. Consider a maximizing sequence (u t ) t∈N , with b T u t → p * k = p * as t → ∞. By feasibility in (3.4), one has a(x) T u t ∞ ≤ 1 for all t and therefore (u t ) ⊂ U := {u ∈ R n : a(x) T u ∞ ≤ 1} and U is compact (see the proof of Lemma 2). Therefore there exists u * ∈ U and a subsequence (t ) ∈N such that u t → u * ∈ U as → ∞. In particular b T u * = p * . Moreover, since by Lemma 8 in the Appendix the convex cone

, which proves that u * is an optimal solution of (3.4).

Assume that the rank conditions of Theorem 1 is satisfied at some relaxation order k, and let (µ * + , µ * -) denote the atomic measures optimal for problem (2.1), obtained from the solution of the primal SDP problem (3.2). Let u * denote an optimal solution of the dual SDP problem (3.4). The duality result of Lemma 4 implies that Supp µ * + ⊂ {x ∈ X : a (x)u * = 1} and Supp µ * -⊂ {x ∈ X : a (x)u * = -1} so that the polynomial a (x)u * can be used as a certificate of optimality. We formulate this in the following dual to Theorem 1.

Lemma 5. Assume that the rank conditions (3.3) of Theorem 1 hold. Let us denote by u * the optimal solution of SOS SDP (3.4). Then the polynomial z * + (x) := 1 + a (x)u * vanishes at the r + points of the support of µ * + , and the polynomial z * -(x) := 1 -a (x)u * vanishes at the r -points of the support of µ * -.

Proof: Let us denote by {x k + } k=1,...,r+ ⊂ X the points of the support of the optimal measure µ * + , computed from the moments y * + solving optimally moment SDP (3.2). By complementarity of the solutions of primal-dual SDP (3.2) and (3.4), it holds z * + , µ * -= 0 and hence z * + , δ x k

The proof is similar for z * -and µ * -.

Discussion

We would like to point out that the developments in this paper were inspired by a previous work on optimal control for linear systems formulated as a primal LP (2.1) on measures and a dual LP on continuous functions (2.2), and solved numerically with primal-dual moment-SOS SDP hierarchies [START_REF] Claeys | Moment lmi approach to ltv impulsive control, Decision and Control (CDC)[END_REF][START_REF]Measures and lmis for impulsive nonlinear optimal control[END_REF]. Formulating optimal control problems as moment problems was a classical research topic in the 1960s, where optimal control laws were sought in measures spaces (completions of Lebesgue spaces) to allow for oscillations and concentrations, see e.g. [START_REF] Krasovskii | Theory of motion control[END_REF] or the overview in [Fat99, Section III]. In the case of linear optimal control of an ordinary


Didier Henrion

Simone Naldi

Mohab Safey

El Din

Keywords: Computer algebra, symbolic computation, linear matrix inequalities, semidefinite programming, low rank matrices, real algebraic geometry

This document describes our freely distributed Maple library spectra, for Semidefinite Programming solved Exactly with Computational Tools of Real Algebra. It solves linear matrix inequalities with symbolic computation in exact arithmetic and it is targeted to small-size, possibly degenerate problems for which symbolic infeasibility or feasibility certificates are required.

Introduction

Given symmetric matrices A 0 , A 1 , . . . , A n of size m with rational coefficients, let

S := {x ∈ R n : A(x) := A 0 + A 1 x 1 + • • • + A n x n 0}
denote the corresponding convex spectrahedron, defined by the linear matrix inequality (LMI) enforcing that A is positive semidefinite, or equivalently that the eigenvalues of A, as functions of x, are all nonnegative. Spectrahedra are a broad generalization of polyhedra [START_REF] Vinzant | What is... a spectrahedron ? Notices of the[END_REF]. Like polyhedra, spectrahedra have facets, edges and vertices. However, while the facets of a polyhedron are necessarily flat, the facets of a spectrahedron can be curved outwards or inflated, see Figure 1 for an example of a spectrahedron of dimension n = 3 defined by an LMI of size m = 5. Optimization of a linear function on a spectrahedron is called semidefinite programming (SDP), a broad generalization of linear programming (LP) with many applications in control engineering, signal processing, combinatorial optimization, mechanical structure design, etc, see [START_REF] Vandenberghe | Semidefinite programming[END_REF][START_REF] Todd | Semidefinite optimization[END_REF]. The algebra and geometry of spectrahedra is an active area of study in real algebraic geometry, especially in connection with the problem of moments and the decomposition of real multivariate polynomials as sums-of-squares (SOS), see [START_REF] Lasserre | Moments, Positive Polynomials and Their Applications[END_REF][START_REF] Blekherman | Semidefinite Optimization and Convex Algebraic Geometry[END_REF] and references therein.

Our software spectra aims at either proving that S is empty, or finding at least one point in S , using exact arithmetic. Contrary to numerical algorithms which are based on approximate computations and floating point arithmetic -such as the projection and rounding heuristics of e.g. [START_REF] Peyrl | Computing sum of squares decompositions with rational coefficients[END_REF] or [START_REF] Kaltofen | Exact certification in global polynomial optimization via sums-of-squares of rational functions with rational coefficients[END_REF] spectra is exclusively based on computations with exact arithmetic. Since exact computations are potentially expensive, spectra should be used when the number n of variables or the size m of the LMI are small. It should not be considered as a competitor to numerical algorithms such as interior-point methods for SDP.

It should be primarily used either in potentially degenerate situations, for example when it is expected that S has empty interior, or when a rigorous certificate of infeasibility or feasibility is required.

The input provided to spectra is the set of matrices A 0 , A 1 , . . . , A n with rational coefficients describing the pencil A and hence the spectrahedron S . If S is empty, spectra returns the empty list. Otherwise, the output generated by spectra is a finite set

Z := q 1 (z) q 0 (z) , q 2 (z) q 0 (z) , • • • , q n (z) q 0 (z) : q(z) = 0, z ∈ C . (1)
described by a collection of univariate polynomials with integer coefficients q, q 0 , q 1 , . . . , q n ∈ Z[z] and such that Z meets S in at least one real point x * . Such a description is called a rational parametrization. It allows to isolate the (generally irrational) coordinates of x * in rational intervals of length given a priori, as small as required.

If S is not empty, spectra is guaranteed to compute a point x * minimizing the rank of A in S . It solves exactly the (non-convex) optimization problem

r(A) := min rank A(x) s.t. x ∈ S . (2)
This is in sharp contrast with interior-point methods which are designed to compute a point of maximal rank.

The outline of the remainder of the paper is as follows. In Section 2 we survey some background material and extract the essential theoretical results of [START_REF] Henrion | Exact algorithms for linear matrix inequalities[END_REF] on which spectra relies.

In Section 3 we provide instructions to download and install spectra, and we illustrate its use on two elementary examples. More advanced examples are described in Section 4. The performance of spectra on larger examples is reported in Section 5. Finally, in Section 6 we describe formally the exact input and output syntaxes of SolveLMI, the main function of spectra.

Background material and main theoretical results

The algorithm implemented in spectra computes points in the determinantal varieties

D r := {x ∈ C n : rank A(x) ≤ r} for r = 0, 1, . . . , m -1. By construction it holds D 0 ⊂ D 1 ⊂ • • • D m-1 .
Since the determinant of A vanishes on the boundary ∂S of S , it holds

∂S ⊂ D m-1 ∩ R n .
When S is not empty, the value r(A) of the optimization problem (2) is the minimum integer r such that D r ∩ R n intersects S , namely the smallest rank on the spectrahedron. Our main geometrical result [4, Theorem 2] states that the spectrahedron contains at least one of the connected components of the real part of the determinantal variety of smallest rank:

Theorem 1 (Smallest rank on a spectrahedron) Assume that S is not empty. Let C be a connected component of D r(A) ∩ R n such that the intersection C ∩ S is not empty. Then C ⊂ S and hence

C ⊂ (D r(A) \D r(A)-1 ) ∩ R n .
As a consequence of this result, an algorithm computing at least one point in each connected component of D r(A) ∩ R n will compute at least one point in the spectrahedron S . Since the value of r(A) is not known beforehand in general, spectra proceeds iteratively by computing at least one point in the real determinantal variety D r ∩ R n for increasing values of the expected rank r = 0, 1, . . . , m -1.

More specifically, spectra computes points in the determinantal varieties D r by projecting onto the subspace of x variables the incidence varieties

V r := {(x, y) ∈ C n × C m(m-r) : A(x)Y (y) = 0, rank Y (y) = m -r}
for r = 0, 1, . . . , m -1. The reader familiar with SDP duality will recognize the classical complementarity conditions, see e.g. [START_REF] Vandenberghe | Semidefinite programming[END_REF][START_REF] Todd | Semidefinite optimization[END_REF]. The dual matrix

Y (y) =    y 1,1 • • • y 1,m-r . . . . . . y m,1 • • • y m,m-r   
depends linearly on the dual variables y, and some normalization constraint should be added to ensure that rank Y (y) = m -r. Unlike D r , the incidence variety V r , up to genericity conditions on the pencil A , turns to be smooth and equidimensional. This crucial geometric property allows for the application of a recursive method which is guaranteed to find at least one point in each connected component of the incidence variety. This leads to the main algorithmic result [4, Theorem 3] on which spectra relies:

Theorem 2 (Exact algorithm for finding a point in a spectrahedron) Suppose that for each r = 0, 1, . . . , m -1, the incidence variety V r is smooth and equidimensional and that its defining polynomial system generates a radical ideal. Suppose also that the determinantal variety D r has the expected dimension n -m-r+1

2

. Then, there is a probabilistic algorithm that takes A as input and returns:

1. either the empty list, if and only if S is empty, or 2. a vector x * such that A(x * ) = 0, if and only if the linear system A(x) = 0 has a solution, or 3. a rational parametrization q, q 0 , q 1 , . . . , q n ∈ Z[z] such that there exists z * ∈ R with q(z * ) = 0 and:

• A(q 1 (z * )/q 0 (z * ), . . . , q n (z * )/q 0 (z * )) 0 and

• rank A(q 1 (z * )/q 0 (z * ), . . . , q n (z * )/q 0 (z * )) = r(A).

The probabilistic nature of the algorithm comes from random changes of variables performed during the procedure, allowing to put the sets D r in generic position.

Recall that the incidence varieties V r are defined by enforcing a full column rank constraint on the dual matrix Y (y). In spectra this is achieved as follows [4, Section 3.1]: given a subset of m -r dinstinct integers between 1 and r, we enforce the submatrix of Y (y) whose rows are indexed by these integers to be equal to the identity matrix of size m -r. For a given value of r, there are m r distinct choices of row indices and hence the same number of normalized incidence varieties. For each value of r, the algorithm in spectra processes iteratively these normalized incidence varieties.

Finally, let us explain briefly how spectra is able to certify the correctness of the output. This explanation was not included in our paper [START_REF] Henrion | Exact algorithms for linear matrix inequalities[END_REF], but we believe it is useful for readers interested in the implementation details. For each computed solution (x * , y * ) belonging to a connected component of an incidence variety, spectra uses exact arithmetic to decide whether A(x * ) is positive semidefinite and to evaluate the rank of A(x * ). If A(x * ) is not positive semidefinite, then the point x * is discarded. From Theorem 1 we know that at least one computed point x * lies on the spectrahedron S , and this point is of minimal rank, i.e. it solves problem (2).

We first build the following characteristic polynomial:

s → p(s, x) = det(sI m + A(x)) = s m + p 1 (x)s m-1 + • • • + p m-1 (x)s + p m (x),
where I m is the identity matrix of size m. The coefficient p k (x) ∈ Q[x] has degree k in x and it is the k-th elementary symmetric polynomial of the eigenvalues of A(x), for k = 1, . . . , m. For example, p 1 (x) is the trace of A(x) and p m (x) is the determinant of A(x). This computation is done only once. Let x * ∈ R n be given. The rank defect of A(x * ) is equal to the number of consecutive zeros in the sequence p m (x * ), p m-1 (x * ), . . . Moreover, by Descartes' rule of signs, A(x * ) 0 if and only if p k (x * ) ≥ 0 for all k = 1, . . . , m. Hence, computing exactly the rank of A(x * ) and checking its positive semidefiniteness amounts to determining the signs of p k (x * ) for k = 1, . . . , m.

Whereas this sign determination is a delicate issue when using floating arithmetic and approximate computation, it can be done exactly with spectra, since we represent the point x * with a rational univariate parametrization with coefficients in Z. Indeed, suppose that x * belongs to the finite set Z described as in (1) by the integer coefficient polynomials q, q 0 , q 1 , . . . , q n . Together with the rational intervals isolating each entry of x * , spectra computes rational intervals isolating each coefficient p k (x * ). Each isolating interval is gradually reduced, until it is so small that at the interval bounds the coefficient takes 1) distinct signs, in which case it vanishes somewhere in the interval, or 2) the same sign, in which case it does not vanish in the whole interval.

3 Getting started spectra is freely available as a library for Maple version 16 and above. It can be downloaded in the form of single binary file SPECTRA.mla from the following page homepages.laas.fr/henrion/software/spectra spectra relies on FGb, a library for fast computation of Gröbner bases, whose Maple interface must be installed, see [START_REF] Faugère | FGb: a library for computing Gröbner bases[END_REF]. spectra does not work without FGb.

In a Maple worksheet, from the directory containing the file SPECTRA.mla, please type the command > with(SPECTRA); to activate the main function SolveLMI.

Half disk

Let A(x) =   1 + x 1 x 2 0 x 2 1 -x 1 0 0 0 x 1  
with n = 2 and m = 3. The corresponding spectrahedron

S = {x ∈ R 2 : A(x) 0} = {x ∈ R 2 : 1 -x 2 1 -x 2 2 ≥ 0, x 1 ≥ 0}
is a half disk. To find a point in S , we use spectra as follows:

> A := Matrix([[1+x1, x2, 0], [x2, 1-x1, 0], [0, 0, x1]]): > SolveLMI(A); [[x1 = [0, 0], x2 = [1, 1]]]
This returns the point x = [0, 1] ∈ S in interval arithmetic notation, i.e.

x 1 ∈ [0, 0], x 2 ∈ [1, 1]
and for each component in x we obtain rational (exact) lower and upper bounds. Here the bounds coincide as the point has rational coordinates.

At this point, matrix A(x) is guaranteed to have minimal rank over all points in S . This rank can be obtained as follows:

> SolveLMI(A,{rnk}); [[x1 = [0, 0], x2 = [1, 1], rnk = 1]]

Degenerate spectrahedra

Let us modify the bottom right entry in the matrix of the previous section, so that now

A(x) =   1 + x 1 x 2 0 x 2 1 -x 1 0 0 0 x 1 -1  
and the corresponding spectrahedron S = {x ∈ R 2 : A(x) 0} = {[1, 0]} reduces to a point in the plane. spectra can easily deal with such a degenerate case:

> A := Matrix([[1+x1, x2, 0], [x2, 1-x1, 0], [0, 0, x1-1]]): > SolveLMI(A); [[x1 = [1, 1], x2 = [0, 0]]]
Now let us modify further the bottom right entry, letting

A(x) =   1 + x 1 x 2 0 x 2 1 -x 1 0 0 0 x 1 -2  
so that the corresponding spectrahedron is empty. spectra returns the empty list, and this is a certificate of emptiness:

> A := Matrix([[1+x1, x2, 0], [x2, 1-x1, 0], [0, 0, x1-2]]): > SolveLMI(A); []
Since spectra is based on exact arithmetic, it is not sensitive to numerical roundoff errors or small parameter changes: Displayed with 10 significant digits, the latter point reads:

x 1 ∈ [ 3689348814741899533 36893488147419103232 , 4611686018427401391 4611686018427387904 ] ≈ 1.000000000, x 2 ∈ [ -350142318592414077 2475880078570760549798248448 , - 2801138548739304423
19807040628566084398385987584 ] ≈ -0.1414213562 • 10 -9 . The above point is an irrational solution, and the rational intervals are provided so that their floating point approximations are correct up to the number of digits specified in the Maple environment variable Digits, which is by default equal to 10. Use the command > Digits:=100: prior to calling SolveLMI if you want an approximation correct to 100 digits. At the price of increased computational burden, spectra then provides larger integer numerators and denominators in the coordinate intervals.

Examples

Irrational spectrahedron

In general, each coordinate of a point computed by spectra is an algebraic number, i.e. the root of a univariate polynomial with integer coefficients.

For the classical univariate matrix

A(x 1 ) =     1 x 1 0 0 x 1 2 0 0 0 0 2x 1 2 0 0 2 x 1    
the spectrahedron reduces to the irrational point

x 1 = √ 2. The simple call > A:=Matrix([[1, x1, 0, 0], [x1, 2, 0, 0], [0, 0, 2*x1, 2], [0, 0, 2, x1]]): > SolveLMI(A); [[x1 = [26087635650665550353 / 18446744073709551616, 13043817825332807945 / 9223372036854775808]]]
returns an interval enclosure valid to 10 digits. We can however obtain an exact representation of this point via a rational parametrization:

> SolveLMI(A, {par}); [[x1 = [..], par = [_Z^2-2,_Z,[2]]]]
The output parameter par contains three univariate polynomials q, q 0 , q 1 such that the computed point is contained in the finite set

Z = {q 1 (z)/q 0 (z) : q(z) = 0} = {2/z : z 2 -2 = 0} = {± √ 2}
as in [START_REF] Blekherman | Semidefinite Optimization and Convex Algebraic Geometry[END_REF]. Here obviously the rational interval isolates the irrational point

x 1 = √ 2.

Algebraic degree

The algebraic degree of semidefinite programming was studied in [START_REF] Nie | The algebraic degree of semidefinite programming[END_REF]. Let us consider the spectrahedron of Example 4 in this reference, for which

A(x) =     1 + x 3 x 1 + x 2 x 2 x 2 + x 3 x 1 + x 2 1 -x 1 x 2 -x 3 x 2 x 2 x 2 -x 3 1 + x 2 x 1 + x 3 x 2 + x 3 x 2 x 1 + x 3 1 -x 3    
The following point can be easily found with spectra, and it has rank 2, which is guaranteed to be the minimal rank achieved in the spectrahedron: With the following instruction we can indeed certify that there is no point of rank 1 or less:

> A:=Matrix([[1+x3,
> SolveLMI(A, {}, [1]); []
The command > SolveLMI(A, {par});

returns the following rational univariate parametrization (1) of the above rank 2 point: q(z) = 16144z 10 + 35160z 9 + 14536z 8 -17690z 7 -16278z 6 -2001z 5 + 1556z 4 + 454z 3 + 23z 2 -4z -1 q 0 (z) = 161440z 9 + 316440z 8 + 116288z 7 -123830z 6 -97668z 5 -10005z 4 + 6224z 3 + 1362z 2 + 46z -4 q 1 (z) = 97248z 9 + 146144z 8 -18192z 7 -134826z 6 -63302z 5 + 4048z 4 + 6758z 3 + 846z 2 -49z -14 q 2 (z) = 34456z 9 + 37516z 8 -8734z 7 -22150z 6 -8223z 5 -3978z 4 -1324z 3 + 104z 2 + 103z + 13 q 3 (z) = -35160z 9 -29072z 8 + 53070z 7 + 65112z 6 + 10005z 5 -9336z 4 -3178z 3 -184z 2 + 36z + 10

The degree of the polynomial q in this parametrization can be obtained with the command > SolveLMI(A, {deg});

We can obtain more points in the spectrahedron as follows:

> SolveLMI(A, {all, rnk, deg}, [START_REF] Chua | Gram spectrahedra[END_REF]);

This returns 4 feasible solutions of rank r = 2, all parametrized by the above degree 10 polynomial. Notice that this degree matches with the algebraic degree of a generic semidefinite programming problem with parameters (m, n, r) = (4, 3, 2), which is 10 according to [8, Table 2].

Reproducibility

Consider the matrix

A(x) = 1 + x 1 x 2 x 2 1 -x 1
modeling the unit disk. Two consecutive calls to SolveLMI return two distinct points: After another call, or on your own computer, these intervals should still differ as spectra makes random changes of coordinates to ensure that the geometric objects computed are in general position. This kind of behavior is expected when there are infinitely many points of minimal rank in the spectrahedron.

Convex quartic

Let

A(x) =     1 + x 1 x 2 0 0 x 2 1 -x 1 x 2 0 0 x 2 2 + x 1 x 2 0 0 x 2 2 -x 1     .
The spectrahedron S = {x ∈ R 2 : A(x) 0} is the orange region whose boundary is the internal oval of the smooth quartic determinantal curve {x ∈ R 2 : det A(x) = 0} represented in black on Figure 2. With the following instructions

> A:=Matrix([[1+x1,x2,0,0],[x2,1-x1,x2,0],[0,x2,2+x1,x2],[0,0,x2,2-x1]]): > SolveLMI(A,{},[3]); > SolveLMI(A,{},[3]); > ...
we compute several points on the boundary of S , they are plotted in red on Figure 2. Note the third input argument [START_REF] Faugère | FGb: a library for computing Gröbner bases[END_REF] which specifies to SolveLMI the expected rank of the computed point. Since the determinantal curve is smooth, we know that the rank of A(x) equals 3 on the whole curve, and in particular on the boundary of S . Since the rank is specified, spectra does not have to process iteratively the incidence varieties corresponding to points of smaller ranks, thereby reducing the computational burden to find at least one point in the spectrahedron.

Each of these points is represented by a rational univariate parametrization of degree 12, obtained with the instruction For example, for the point (x 1 , x 2 ) ≈ (-0.9689884394, -0.2434013983) the polynomial q in the rational parametrization (1) is q(z) = 5506034827600 z 12 -4608031295324 z 10 -192908794368 z 9 + 25693318717857 z 8 + + 4774492660608 z 7 -17188212283956 z 6 -23438418515712 z 5 + 64967482316484 z 4 --11285164470528 z 3 -11887630039728 z 2 + 296990121024.

Recall that the algebraic degree of a point x * in S is the degree of the minimal algebraic extension of the ground field (here the rational numbers) required to represent x * . The algebraic degree depends on the size of the pencil A but also on the rank r of A(x * ). With (m, n, r) = (4, 2, 3) and generic data, the algebraic degree is 12, cf. [8, Table 2], which indeed coincides with the degree of the exact representation of x * computed by spectra.

Polynomial sums of squares

Deciding whether a multivariate real polynomial is non-negative is difficult in general. A sufficient condition, or certificate for non-negativity, is that the polynomial can be expressed as a sum of squares (SOS) of other polynomials. Finding a polynomial SOS decomposition amounts to finding a point in a specific spectrahedron called Gram spectrahedron, see e.g.

[2] and references therein.

As an example, consider the homogeneous ternary quartic

f (u) = u 4 1 + u 1 u 3 2 + u 4 2 -3u 2 1 u 2 u 3 -4u 1 u 2 2 u 3 + 2u 2 1 u 2 3 + u 1 u 3 3 + u 2 u 3 3 + u 4 3 .
The polynomial f belongs to a series of examples provided by C. Scheiderer in [START_REF] Scheiderer | Sums of squares of polynomials with rational coefficients[END_REF] to answer (in the negative) the following question by B. Sturmfels: let f be a polynomial with rational coefficients which is an SOS of polynomials with real coefficients; is it an SOS of polynomials with rational coefficients? Scheiderer's counterexamples prove that, generally speaking, there is no hope in obtaining nonnegativity certificates over the rationals. However, certificates exist in some algebraic extension of the field of rational numbers.

In the graded reverse lexicographic ordered monomial basis, the Gram matrix of f is the matrix

A(x) =         1 0 x 1 0 -3/2 -x 2 x 3 0 -2x 1 1/2 x 2 -2 -x 4 -x 5 x 1 1/2 1 x 4 0 x 6 0 x 2 x 4 -2x 3 + 2 x 5 1/2 -3/2 -x 2 -2 -x 4 0 x 5 -2x 6 1/2 x 3 -x 5 x 6 1/2 1/2 1        
depending linearly on 6 real parameters. The Gram spectrahedron S = {x ∈ R 6 : A(x) 0} parametrizes the set of all SOS decompositions of f . We deduce by the discussion above that S does not contain rational points. In particular, its interior is empty. Let us use spectra to compute points in S and hence to get positivity certificates for f : We obtain an irrational point x ∈ S whose coordinates are algebraic numbers of degree Let us compute more non-negativity certificates of rank 2:

> A := Matrix([[1,0,x1,0,-3/2-x2,x3], [0,-2*x1,1/2,x2,-2-x4,-x5], [x1,1/2,1,x4,0,x6], [0,x2,x4,-2*x3+2,x5,1/2], [-3/2-x2,-2-x4,0,x5,-2*x6,1/2], [x3,-x5,x6,1/2,1/2,1]]): > SolveLMI(A,
> SolveLMI(A,{rnk,deg,par,all}, [START_REF] Chua | Gram spectrahedra[END_REF]);

In addition to the point already obtained above, we get another point. The user can compare this output with [10, Ex. 2.8]: it turns out that these are the only 2 points of rank 2. Other points in the Gram spectrahedron have rank 4 and they are convex combinations of these 2.

Performance

Exponential bit-size spectrahedron

For a given n ∈ N, consider the spectrahedron

S n = x ∈ R n : 1 2 2 x 1 0, 1 x 1 x 1 x 2 0, • • • , 1 x n-1 x n-1 x n 0 . For every x * ∈ S n it holds x * n ≥ (x * n-1 ) 2 ≥ • • • ≥ (x * 1 ) 2 n-1 ≥ 2 2 n
, which shows that exponentially many bits are required to represent a point. It is elementary to check that each of the above n matrices of size 2 can have rank 1, and hence that we can compute a point of rank n as follows:

> with(LinearAlgebra): Recall from Section 2 that spectra examines iteratively a family of m r = 2n n incidence varieties, a number growing exponentially with n. For example there are 12870 = 16 8 incidence varieties to test to solve our problem for n = 8. Hence we could expect spectra to perform poorly on this example. However, on our standard desktop PC equipped with Intel i7 processor at 2.5GHz and 16GB RAM, we were able to handle spectrahedra of size 2n = 10 in 29 seconds, and of size 2n = 12 in 505 seconds. For m = 2, r = 1 and values of n ranging from 30 to 100, we obtain the timings reported on Figure 3. This corresponds to spectrahedra whose boundaries belong to determinantal varieties of increasing dimension. Moreover, the singularity locus of the determinant has positive dimension, it is a linear subspace of co-dimension 3. We observe a polynomial dependence of the computational time as a function of the number of variables, with exponent around 3.

Random spectrahedra

When m = 3, r = m -1 = 2 and values of n ranging from 30 to 80, we obtain the timings reported on Figure 4, depending polynomially on n with an exponent around 4. 6 Input and output syntax

Input

The calling sequence of function SolveLMI is as follows:

> SolveLMI(A, options, ranks); where • A is a symmetric matrix of size m with rational coefficients, depending affinely on n variables;

• options (optional) is a set that can contain the following keywords:

all : compute as many solutions as possible, which can be computationally demanding; when this option is not specified, the algorithm is stopped as soon as one solution is computed, which is typically much faster;

rnk : return the rank of A at every computed solution;

par : return the rational univariate parametrization of every computed solution;

deg : return the algebraic degree of every computed solution;

• ranks (optional) is a list of nonnegative integers corresponding to expected ranks of computed solutions. The default value is [0, 1, . . . , m -1]. The algorithm is run for each value r in ranks by solving the quadratic system of equations A(x)Y (y) = 0 for a vector x and a matrix Y (y) with m rows and m -r columns whose entries are stored in a vector y. It may happen that the rank of A(x) at a computed solution x is strictly less than r.

Output

Let us denote by x 1 , x 2 , . . . , x n the variables on which matrix A depends affinely. They are gathered in a vector x ∈ R n . When the input argument options is empty, the output returned by SolveLMI is where a i , b i are rational numbers, displayed as ratios of integers. This means that each coordinate x i belongs to the interval [a i , b i ] ensuring a floating point approximation of x valid to a number of digits equal to the Maple environment variable Digits. When a i = b i this implies that x i is a rational number.

When options contains the keyword all, more points can be returned, in the form of a list such that nops(out) is the number of computed points, out [START_REF] Blekherman | Semidefinite Optimization and Convex Algebraic Geometry[END_REF] is the first point, out [START_REF] Chua | Gram spectrahedra[END_REF] is the second point, etc.

When options contains the keyword rnk, the rank of A at x is returned:

Figure 1 :

 1 Figure 1: A spectrahedron.

>

  A := Matrix([[1+x1, x2, 0], [x2, 1-x1, 0], [0, 0, x1-1-10^(-20)]]): > SolveLMI(A); [] > A:=Matrix([[1+x1, x2, 0], [x2, 1-x1, 0], [0, 0, x1-1+10^(-20)]]): > SolveLMI(A); [[x1 = [36893488147418995335 / 36893488147419103232, 4611686018427401391 / 4611686018427387904], x2 = [-350142318592414077 / 2475880078570760549798248448, -2801138548739304423 / 19807040628566084398385987584]]

>

  A:=Matrix([[1+x1,x2],[x2,1-x1]]): > SolveLMI(A); [[x1 = [-21201056044062027875/36893488147419103232, -662533001376936933/1152921504606846976], x2 = [-7548363607018988253/9223372036854775808, -1887090901754742967/2305843009213693952]]] > SolveLMI(A); [[x1 = [-10862500438565607907/590295810358705651712, -21725000877131177215/1180591620717411303424], x2 = [-576363141759828805/576460752303423488, -9221810268157244495/9223372036854775808]]]

Figure 2 :

 2 Figure 2: Quartic curve (black) with sample points (red) on the boundary of its spectrahedron (orange).

  {rnk, deg, par}); [[[x1 = [..], x2 = [..], x3 = [..], x4 = [..], x5 = [..], x6 = [..]], rnk = 2, deg = 3, par = [8*z^3+8*z+1, 24*z^2-8, [16*z+3, -24*z^2+8, 8*z^2+6*z+8, -16*z^2+6*z+16, -16*z-3, 16*z+3]]]

  > A:=DiagonalMatrix([<<1,2>|<2,x1>>,<<1,x1>|<x1,x2>>,<<1,x2>|<x2,x3>>,..]): > SolveLMI(A,{},[n]);

Finally, we report

  on randomly generated examples. The rational entries of A are generated as quotients of integers drawn uniformly in the interval [-100, 100]. Here is the script we used to generate a random symmetric pencil given the number n of variables and the size m: > var:=[seq(cat('x',i),i=1..n)]: > A:=Matrix(m,m): > for i from 1 to m do for j from i to m do A[i,j]:=randpoly(var, degree=1, dense, coeffs=rand(-100..100)): A[j,i]:=A[i,j]; od: od: For each instance, given the expected rank r, we execute the command > SolveLMI(A,{},[r]);

Figure 3 :Figure 4 :

 34 Figure 3: Timings for random instances of size m = 2 and rank r = 1, as a function of the number of variables n.

•

  eithter the empty list [] in which case S is empty, or • a rational enclosure of a single point x ∈ S , in the form > SolveLMI(A) [[x1 = [a1, b1], x2 = [a2, b2], ..., xn = [an, bn]]]

>

  out := SolveLMI(A, {all}) [[x1 = [a11, b11], x2 = [a12, b12], ..., xn = [a1n, b1n]], [x1 = [a21, b21], x2 = [a22, b22], ..., xn = [a2n, b2n]], ...]

  To generate reproducible outputs, the instruction randomize can be used to seed the random number generator used by Maple:

	> randomize(31415926):
	> SolveLMI(A);
	[[x1 = [-35204733513421104993/36893488147419103232, -35204733513421000447/36893488147419103232],
	x2 = [-2758579864857623899/9223372036854775808, -5517159729715231413/18446744073709551616]]]
	> randomize(31415926):
	> SolveLMI(A);
	[[x1 = [-35204733513421104993/36893488147419103232, -35204733513421000447/36893488147419103232],
	x2 = [-2758579864857623899/9223372036854775808, -5517159729715231413/18446744073709551616]]]

  At this point, the Gram matrix A has rank 2, and hence f is an SOS of 2 polynomials.

												3,
	and which belongs to the finite set
	16z + 3 24z 2 -8	,	-24z 2 + 8 24z 2 -8	,	8z 2 + 6z + 8 24z 2 -8	,	-16z 2 + 6z + 16 24z 2 -8	,	-16z -3 24z 2 -8	,	16z + 3 24z 2 -8	: 8z 3 -8z -1 = 0

* This work was partly funded by the ERC Advanced Grant Taming.

When options contains the keyword par, a rational univariate parametrization of x is returned:

.,qn]]]

This parametrization is such that q, q 0 , q 1 , q 2 , . . . , q n are univariate polynomials with integer coefficients such that x belongs to the finite set

The intervals [a i , b i ] are provide to isolate the computed point from this set of points.

When options contains the keyword deg, the degree of the polynomial q in the rational univariate parametrization of each computed point x is also returned: