11footnotetext: Institute of Science and Technology Austria
(ISTA)22footnotetext: Université de Limoges, CNRS, XLIM, UMR 7252, F-87000
Limoges, France.33footnotetext: Sorbonne Université, CNRS, LIP6, F-75005
Paris, France.

# Verifying feasibility of degenerate semidefinite programs

Vladimir Kolmogorov1
<EMAIL_ADDRESS>Simone Naldi2,3
<EMAIL_ADDRESS>Jeferson Zapata1
<EMAIL_ADDRESS>
###### Abstract

This paper deals with the algorithmic aspects of solving feasibility problems
of semidefinite programming (SDP), aka linear matrix inequalities (LMI). Since
in some SDP instances all feasible solutions have irrational entries,
numerical solvers that work with rational numbers can only find an approximate
solution. We study the following question: is it possible to certify
feasibility of a given SDP using an approximate solution that is sufficiently
close to some exact solution? Existing approaches make the assumption that
there exist rational feasible solutions (and use techniques such as rounding
and lattice reduction algorithms).

We propose an alternative approach that does not need this assumption. More
specifically, we show how to construct a system of polynomial equations whose
set of real solutions is guaranteed to have an isolated correct solution
(assuming that the target exact solution is maximum-rank). This allows, in
particular, to use algorithms from real algebraic geometry for solving systems
of polynomial equations, yielding a hybrid (or symbolic-numerical) method for
SDPs. We experimentally compare it with a pure symbolic method in [22]; the
hybrid method was able to certify feasibility of many SDP instances on which
[22] failed. We argue that our approach may have other uses, such as refining
an approximate solution using methods of numerical algebraic geometry for
systems of polynomial equations.

Keywords. Symbolic-numeric algorithms; linear matrix inequalities; weak
feasibility; facial reduction; spectrahedra.

##  1 Introduction

Semidefinite programming (SDP) is a nontrivial generalization of linear
programming (LP) asking to minimize linear functions defined over the space of
real symmetric matrices and restricted to the PSD cone (the convex cone of
positive semidefinite matrices). Despite its simple definition the complexity
analysis of SDP and of its feasibility problem (LMI, linear matrix
inequalities) is still a fairly open question, and likewise the development of
efficient but reliable algorithms, which is at the heart of research in convex
optimization.

Beyond its interest as conic optimization problem, the importance of SDP is
enhanced by its numerous applications, the most classical of which is perhaps
the semidefinite relaxation of MAX-CUT by Goemans and Williamson [13]. More
generally, SDP is used as numerical tool for solving hard non-convex
optimization problems through the so-called moment-SOS hierarchy [20],
consisting in relaxing polynomial optimization to a sequence of SDP problems
of increasing size and with good convergence properties [28]. Specific
applications include checking stability of systems in control theory [36] and
analyzing the convergence rate of numerical algorithms for convex optimization
[10]. In both cases one needs to find a certain Lyapunov function, and the
search for such function is cast as a small-scale SDP.

Unfortunately, naively relying on numerical SDP solvers may lead to incorrect
conclusions; examples for small dynamical systems can be found in [42]. This
motivates the need for algorithms that can certify feasibility of a given SDP.
Below we discuss two classes of such algorithms: symbolic methods that compute
an exact solution in a Rational Univariate Representation using techniques
from real algebraic geometry, and symbolic-numerical (or hybrid) methods that
certify feasibility using a numerical SDP solver as a subroutine.

Symbolic methods  Because the constraints in SDP are nonlinear in the entries
of the matrix, but still algebraic (defined by polynomial inequalities), the
solution of a SDP program is in general not rational but defined over some
algebraic extension of the base field, which we usually assume to be
$\mathbb{Q}$. The degree of this extension, known as the algebraic degree of
SDP, represents a measure of the inner complexity of the program, for which
exact formulas are known [34].

The underlying algebraic structure of SDP has motivated in the last years the
development of computer algebra algorithms whose arithmetic complexity is
essentially quadratic in some multilinear bound on the mentioned algebraic
degree [22, 31]. The strategy of the series of work [21, 22, 24, 32] relies on
reducing the problem of finding a matrix in a section of the PSD cone to the
one of finding low rank elements in some affine space of matrices. This can be
cast as a real-root finding problem, asking to compute at least one point in
every real connected component of an algebraic set, modeled using systems of
polynomial equations with finitely-many solutions and hence solved using
Gröbner-bases-based algorithms [5]. By its purely symbolic nature, these
algorithms are not comparable, in terms of target problem size, with numerical
methods for SDP that are mainly based on variants of interior-point method
(IPM), and their correctness depend on genericity assumptions on input data.

Concerning algebraic approaches to SDP, let us mention the method in [33]
based on homogenization and projective geometry for identifying the
feasibility type of a LMI and the certificate in [27] based on the theory of
sums of squares polynomial for checking SDP infeasibility.

Hybrid methods  These methods can potentially handle much larger SDP
instances. Most of the existing work has focused on the case when the set of
feasible solutions contains rational points [17, 26, 30, 40, 41, 42, 9]. Some
of these papers [30, 9] employ a lattice reduction technique, namely the
Lenstra-Lenstra-Lovász (LLL) algorithm. The paper [9] discusses an extension
to quadratic fields $\mathbb{Q}[\sqrt{\ell}]$, but the method seems to assume
that $\ell$ is known.

Let us finally mention the method by numerical algebraic geometry in [19]
which does not strictly fall into the two classes above. It is based on SDP
duality and applies homotopy continuation for numerically tracking the central
path.

Our contributions  In this paper we develop an alternative hybrid method which
is not restricted to the case when the feasible region contains rational
points. The method essentially reduces the problem of certifying feasibility
of an SDP to that of solving a system of polynomial equations whose set of
real solutions has a zero-dimensional component with a desired solution. Our
main result can be stated informally as follows (a precise formulation is
given later in Theorems 6 and 7).

######  Theorem 1 (informal).

There exists an algorithm that, given a system $\mathscr{A}(X)=b,\,X\succeq 0$
and an approximate solution $\tilde{X}$, constructs a system of polynomial
equations over symmetric matrices $X\in\mathbb{R}^{n\times n}$ and auxiliary
variables $Y\in\mathbb{R}^{n\times(n-r)}$ of the following form:

$\displaystyle\mathscr{A}(X)$
$\displaystyle\;\;=\;\;b\;,\;\;\;X^{\top}\;\;=\;\;X$ (1a)
$\displaystyle\Pi^{\top}X\Pi\left(\begin{matrix}Y\\\
I_{n-r}\end{matrix}\right)$ $\displaystyle\;\;=\;\;0^{n\times(n-r)}$ (1b)
$\displaystyle X_{ij}$
$\displaystyle\;\;=\;\;\tilde{X}_{ij}\qquad\forall(i,j)\in K$ (1c)

where $\Pi$ is a permutation matrix, $r\in[n]$, and $K$ is a subset of
$\\{(i,j)\in[n]\times[n]\>:\>i\leq j\\}$. (Quantities $\Pi$, $r$ and $K$ are
computed from $\mathscr{A}$, $b$ and $\tilde{X}$). If $\tilde{X}$ is
sufficiently close to some maximum-rank solution $X^{\ast}\in
P{\stackrel{{\scriptstyle\mbox{\tiny\tt~{}def~{}}}}{{=}}}\\{X\succeq
0\>:\>\mathscr{A}(X)=b\\}$ then the set of real solutions of system (1) has a
zero-dimensional component $\\{(X,Y)\\}$ with $X\in P$.

There may be several ways to exploit this result. One possibility could be to
refine the approximate solution $\tilde{X}$ and the associated matrix
$\tilde{Y}$ (which is easily computable from $\tilde{X}$) by applying the
Newton method for system (1). For some polynomial systems this method
converges very slowly near the optimum or even diverges; in that case one
could try to solve (1) using techniques from numerical algebraic geometry such
as those implemented in the software Bertini [4]. The approximate solution
$(\tilde{X},\tilde{Y})$ could allow the method to focus on the desired zero-
dimensional component, and use the regime known as the endgame.

In this paper we explore a different direction: we solve system (1) using
exact methods from real algebraic geometry, similar to [22]. Accordingly, our
experimental results focus on a comparison with [22]. Note that the latter
paper also reduces the problem to a system of polynomial equations of the form
(1a)-(1b), but with the following differences: (1) the algorithm in [22] has
no a priori information about a feasible rank and thus needs to exhaustively
search over all (exponentially many) kernel profiles until a solution is
found; (2) it provides a guarantee that for one of these kernel profiles the
set of real solutions contains a (possibly positive-dimensional) connected
component whose points $(X,Y)$ satisfy $X\in P$; (3) it searches for a
minimum-rank solution, whereas we search for a maximum-rank solution.

As in [22], we use a critical point method to find one real solution per
component. This method works under certain assumptions (such as complete
intersection of the corresponding variety or finiteness of the number of
solutions of some critical equations); on the instances that we tested these
assumptions were often violated for the method in [22] but not for our hybrid
method, and accordingly the hybrid method was able to certify feasibility of
many instances on which [22] failed.

Outline  The paper is organized as follows. Section 2 contains the main
theoretical ingredients and notations. In Section 3 we make a link between
determinantal varieties and facial reduction for SDP. Section 4 describes a
hybrid algorithm, which is applied to a collection of instances of degenerate
SDPs of the literature in Section 5.

##  2 Preliminaries

Notation for matrices  All norms $||\cdot||$ for vectors and matrices used in
this paper are 2-norms. The Frobenius norm for matrices is denoted as
$||\cdot||_{F}$.

For a matrix $A\in\mathbb{R}^{p\times q}$ we denote by $\sigma_{i}(A)\geq 0$
the $i$-th singular value of $A$; if $i>\min\\{p,q\\}$ then $\sigma_{i}(A)=0$
by definition. We denote by $\mathbb{S}^{p}_{k}$ the set of $p\times p$
symmetric matrices with coefficients in a field $k$ and for $k=\mathbb{R}$ we
simply write $\mathbb{S}^{p}$. If $A$ is symmetric of size $p$
($A\in\mathbb{S}^{p}$) then $\lambda_{i}(A)$ denotes the $i$-th largest
eigenvalue of $A$. A matrix $A\in\mathbb{S}^{p}$ is positive semidefinite
($A\succeq 0$) if $\lambda_{p}(A)\geq 0$ and positive definite ($A\succ 0$) if
$\lambda_{p}(A)>0$. Recall that if $A\succeq 0$ then
$\lambda_{i}(A)=\sigma_{i}(A)$ for all $i$.

The following inequalities are well known (see e.g. [14, Section 8]):

$\displaystyle|\sigma_{i}(A)-\sigma_{i}(B)|$ $\displaystyle\leq$
$\displaystyle||A-B||\qquad\forall A,B\in\mathbb{R}^{p\times q},\;\forall i$
(2a) $\displaystyle|\lambda_{i}(A)-\lambda_{i}(B)|$ $\displaystyle\leq$
$\displaystyle||A-B||\qquad\forall A,B\in\mathbb{S}^{p},\;\forall i$ (2b)

We will also denote $\rho(A)=\sigma_{r}(A)$ where $r=\text{rank}(A)$. Clearly,
for every non-zero matrix $A$ we have $\rho(A)>0$.

For a set of rows $I\subseteq[p]$ let $A_{I}\in\mathbb{R}^{|I|\times q}$ be
the submatrix of $A$ formed by the rows in $I$. Similarly, for a set of
columns $J\subseteq[q]$ let $A^{J}\in\mathbb{R}^{p\times|J|}$ be the submatrix
of $A$ formed by the columns in $J$.

Semidefinite programming  Let us equip the vector space $\mathbb{S}^{n}$ of
real symmetric matrices of size $n$ with the Frobenius inner product:

$\left\langle
M,N\right\rangle_{\mathbb{S}}:=\text{Trace}(MN)=\sum_{i,j}M_{ij}N_{ij}.$

for $M=(M_{ij}),\,N=(N_{ij})\in\mathbb{S}^{n}$. The cone of positive
semidefinite matrices (PSD cone) is
$\mathbb{S}^{n}_{+}=\\{M\in\mathbb{S}^{n}:M\succeq 0\\}$ is a convex closed
basic semialgebraic set in $\mathbb{S}^{n}$, and its interior is the open
convex cone of positive definite matrices:
$\mathbb{S}^{n}_{++}:=\text{int}(\mathbb{S}^{n}_{+})\subseteq\mathbb{S}^{n}$.
We denote by $\left\langle u,v\right\rangle_{\mathbb{R}}=\sum u_{i}v_{i}$ the
Euclidean inner product on $\mathbb{R}^{m}$.

A semidefinite program in standard primal form is given by

$\begin{array}[]{rcll}p^{*}&:=&\inf_{X\in\mathbb{S}^{n}}&\left\langle
C,X\right\rangle_{\mathbb{S}}\\\ &&\text{s.t.}&\mathscr{A}(X)=b\\\
&&&X\in\mathbb{S}^{n}_{+}\end{array}$ (3)

where $\mathscr{A}\,\,:\,\,\mathbb{S}^{n}\to\mathbb{R}^{m}$ is a linear map
defined by $\mathscr{A}(X)=(\left\langle
A_{1},X\right\rangle,\ldots,\left\langle A_{n},X\right\rangle)$ for some real
symmetric matrices $A_{1},\ldots,A_{n}\in\mathbb{S}^{n},$
$b\in\mathbb{R}^{m}$, and $C\in\mathbb{S}^{n}$. The feasible sets of the
program (3) is denoted by

${{P}}:=\\{X\in\mathbb{S}^{n}\,\,:\,\,\mathscr{A}(X)=b,X\succeq 0\\}.$ (4)

Problem (3) with $P\neq\varnothing$ is called strongly feasible if $P$
contains a matrix $X$ with $X\succ 0$, and weakly feasibly otherwise. We will
be mainly interested in weakly feasible SDPs. (Note that if (3) is strongly
feasible then $P$ contains rational-valued matrices.)

Algebraic and semialgebraic geometry  We refer to [8] for basics of
(computational) algebraic geometry, and we recall here the main definitions
needed in this work. A (complex) algebraic variety is a set of the form

$V=\\{x\in\mathbb{C}^{n}:f_{1}(x)=0,\ldots,f_{s}(x)=0\\}$ (5)

for polynomials $f_{1},\ldots,f_{s}\in K[x]$ with coefficients in some
subfield $K\subseteq\mathbb{C}$ (we usually fix $K=\mathbb{Q}$ or
$K=\mathbb{R}$). Of course defining $I=\langle
f_{1},\ldots,f_{s}\rangle\subseteq K[x]$ the ideal generated by
$f_{1},\ldots,f_{s}$, remark that every $f\in I$ vanishes over $V$, in other
words $V=V(I):=\\{x\in\mathbb{C}^{n}:f(x)=0,\forall f\in I\\}$. The vanishing
ideal of a set $W\subseteq\mathbb{C}^{n}$ is defined as

$I(W)=\\{f\in\mathbb{C}[x]:f(x)=0,\,\forall x\in W\\}$

and by Hilbert’s Nullstellensatz [8, §4.1], one has the correspondence
$I(V(J))=\sqrt{J}$, where $\sqrt{J}:=\\{f\in\mathbb{C}[x]:f^{m}\in J,\exists
m\in\mathbb{N}\\}$ is the radical ideal of $J$. An ideal $I$ is called radical
if $I=\sqrt{I}$. The smallest algebraic variety containing a set $W$ is
denoted by $\overline{W}$ and called its Zariski-closure.

A variety $V\subseteq\mathbb{C}^{n}$ is called irreducible if it is not the
union of two proper subvarieties, and every variety is finite union of
irreducible varieties called its irreducible components. The dimension of an
algebraic variety $V\subseteq\mathbb{C}^{n}$ is the integer $d=\dim(V)$ such
that the intersection $V\cap H_{1}\cap\cdots\cap H_{d}$ of $V$ with $d$
generic hyperplanes $H_{i}$ is finite. The degree of an irreducible variety
$V\subseteq\mathbb{C}^{n}$ is the cardinality of the intersection of $V$ with
$\dim(V)$ generic hyperplanes, and the degree of any variety is the sum of the
degrees of its irreducible components. A variety whose irreducible components
have the same dimension is called equidimensional. An ideal $I$ with
$V(I)\subseteq\mathbb{C}^{n}$ is in complete intersection if it can be
generated by $\text{codim}(V(I))=n-\dim(V(I))$ many polynomials, and with
abuse of notation $V$ is said in complete intersection if this is the case for
$I(V)$.

Let $V\subseteq\mathbb{C}^{n}$ be a variety with $I(V)=\langle
f_{1},\ldots,f_{s}\rangle$. For a points $p\in V$ the tangent space $T_{p}V$
at $p$ to $V$ is the kernel of the Jacobian matrix
$J(f_{1},...,f_{s}):=({\partial f_{i}}/{\partial x_{j}})_{ij}$ evaluated at
$p$. A point $p\in V$ is said to be a regular point if the dimension of the
tangent space $T_{p}V$ at $p$ coincides with the local dimension of $V$ at $p$
(the maximum of the dimensions of the irreducible components containing $p$).
When $V$ is equidimensional of dimension $d$, the Jacobian
$J(f_{1},...,f_{s})$ has rank $\leq n-d$ and equality holds exactly at regular
points. We denote the set of regular points as $\text{Reg}(V)$, and its
complement, $\text{Sing}(V):=V\setminus\text{Reg}(V)$, is referred to as the
set of singular points of $V$.

For $K=\mathbb{R}$ and $V$ as in (5), the set
$\mathbb{R}V:=V\cap\mathbb{R}^{n}$ is called a real algebraic variety, and it
is defined by the vanishing of finitely many real polynomials. Given
$g_{1},\ldots,g_{t}\in\mathbb{R}[x]$, the set

$S=\\{x\in\mathbb{R}^{n}:g_{1}(x)\geq 0,\ldots,g_{t}(x)\geq 0\\}$

is called a basic closed semialgebraic set.

Determinantal and incidence varieties  We refer to [12, 16] for the general
theory of determinantal varieties and we recall below the main definitions in
the context of LMI.

Let $A_{1},\ldots,A_{m}\in\mathbb{S}^{n}$, $b\in\mathbb{R}^{m}$ be the
matrices and the vector defining the map $\mathscr{A}$ in Equation 3. The set

${\mathscr{D}}_{r}=\\{X\in\mathbb{S}^{n}_{\mathbb{C}}\,\,:\,\,\text{rank}(X)\leq
r,\mathscr{A}(X)=b\\}$

is called the determinantal variety associated with $\mathscr{A}$. It is the
algebraic variety defined by the vanishing of the $(r+1)\times(r+1)$ minors of
the matrix $X$ (plus linear constraints $\mathscr{A}(X)=b$). For some
$r\in\\{0,1,\ldots,m\\}$, the real variety
$\mathbb{R}{\mathscr{D}}_{r}:={\mathscr{D}}_{r}\cap\mathbb{S}^{n}$ intersects
the feasible set $P$.

From a computational point of view, minors are not easy to handle. Instead,
the rank constraint is equivalent to the existence of a matrix of high rank
whose columns generate the kernel of the original matrix. The following set

${\mathscr{V}}_{r}=\\{(X,Y)\in\mathbb{S}^{n}_{\mathbb{C}}\times\mathbb{C}^{n\times(n-r)}\,\,:\,\,\mathscr{A}(X)=b,XY=0,\text{rank}(Y)=n-r\\}$

is associated with the variety ${\mathscr{D}}_{r}$; indeed,
${\mathscr{D}}_{r}$ is the projection of ${\mathscr{V}}_{r}$ to
$\mathbb{S}^{n}_{\mathbb{C}}$. Remark that ${\mathscr{V}}_{r}$ is not an
algebraic variety, since the constraint on the rank of $Y$ is not closed in
the Zariski topology, but instead a constructible set (an open subset of a
Zariski closed set). We use the following view of ${\mathscr{V}}_{r}$ from the
coordinate charts (see [22]):

${\mathscr{V}}_{r,\iota}=\\{(X,Y)\in\mathbb{S}^{n}_{\mathbb{C}}\times\mathbb{C}^{n\times(n-r)}\,\,:\,\,\mathscr{A}(X)=b,XY=0,Y_{[n]-\iota}=I_{n-r}\\}$
(6)

for $\iota\subseteq[n]=\\{1,\ldots,n\\}$, of cardinality $|\iota|=r$, and
$Y_{\iota^{\prime}}$ is the matrix obtained by selecting the rows of $Y$ in
$\iota^{\prime}=[n]-\iota$. We define
${\mathscr{W}}_{r,\iota}=\\{X\in\mathbb{S}^{n}_{\mathbb{C}}\,\,:\,\,\exists
Y\in\mathbb{C}^{n\times(n-r)}\;\mbox{s.t.}\;(X,Y)\in{\mathscr{V}}_{r,\iota}\\}$
to be the projection of ${\mathscr{V}}_{r,\iota}$ to
$\mathbb{S}^{n}_{\mathbb{C}}$. Note that ${\mathscr{V}}_{r,\iota}$ is an
algebraic set. Set $\mathbb{R}{\mathscr{W}}_{r,\iota}$ may not be algebraic;
it is a constructible set. Let
$\mathbb{R}{\mathscr{V}}_{r,\iota}:=(\mathbb{S}^{n}\times\mathbb{R}^{n\times(n-r)})\cap{\mathscr{V}}_{r,\iota}$
and
$\mathbb{R}{\mathscr{W}}_{r,\iota}:=\mathbb{S}^{n}\cap{\mathscr{W}}_{r,\iota}$
be the real traces of these two complex varieties. It can be seen that
$\mathbb{R}{\mathscr{W}}_{r,\iota}$ is the projection of
$\mathbb{R}{\mathscr{V}}_{r,\iota}$ to $\mathbb{S}^{n}$. (This follows from
the following observation: if $(X,Y)\in{\mathscr{V}}_{r,\iota}$ with
$X\in\mathbb{S}^{n}$ then
$(X,{\mathscr{R}e}\,Y)\in\mathbb{R}{\mathscr{V}}_{r,\iota}$).

Incidence varieties are used in [22] to solve feasibility problems of generic
instances of LMI, and the following Theorem is a refinement of results in
[22].

######  Theorem 2.

Let $X^{\ast}\in\mathbb{S}^{n}_{+}$ be a minimum-rank solution of (4), and let
$\iota\subseteq[n]$ be a maximal set of linearly independent columns of
$X^{\ast}$, with $r=|\iota|$.

  * (a)

Real variety $\mathbb{R}{\mathscr{W}}_{r,\iota}$ has a connected component $C$
with $X^{\ast}\in C\subseteq{{P}}$.

  * (b)

Real variety $\mathbb{R}{\mathscr{V}}_{r,\iota}$ has a connected component
$C^{\prime}$ with $X^{\ast}\in\pi_{\mathbb{S}^{n}}(C^{\prime})\subseteq{{P}}$
where $\pi_{\mathbb{S}^{n}}(C^{\prime})=\\{X\>:\>\exists Y\mbox{ s.t.\
}(X,Y)\in C^{\prime}\\}$ is the projection of $C^{\prime}$ to
$\mathbb{S}^{n}$.

###### Proof.

The definition of $\iota$ implies that there exists
$Y^{\ast}\in\mathbb{R}^{n\times(n-r)}$ with $X^{\ast}Y^{\ast}=0$ and
$Y^{\ast}_{[n]-\iota}=I_{n-r}$. Clearly, we have
$(X^{\ast},Y^{\ast})\in\mathbb{R}{\mathscr{V}}_{r,\iota}$ and hence
$X^{\ast}\in{\mathscr{W}}_{r,\iota}$. Let $C$ be the component of
${\mathscr{W}}_{r,\iota}$ containing $X^{\ast}$.

(a)  We need to show that $C\subseteq P$. If $r=0$ then all points
$(X,Y)\in\mathbb{R}{\mathscr{V}}_{r,\varnothing}$ satisfy $X=0$, hence
$C=\\{X^{\ast}\\}=\\{0\\}$ and the claim is straightforward. Suppose that
$r>0$ and the claim is false, then there exists continuous curve
$\\{X^{(t)}\\}_{t\in[0,1]}\subseteq C\subseteq{\mathscr{W}}_{r,\iota}$ with
$X^{((0)}=X^{\ast}$ and $X^{(1)}=X$. For each $t\in[0,1]$ there must exist
matrix $Y^{(t)}$ with $(X^{(t)},Y^{(t)})\in\mathbb{R}{\mathscr{V}}_{r,\iota}$.
Denote $\lambda^{(t)}_{i}=\lambda_{i}(X^{(t)})$ for $t\in[0,1]$, so that
$\lambda^{(t)}_{1}\geq\ldots\geq\lambda^{(t)}_{n}$. These values satisfy the
following properties for each $t\in[0,1]$:

  * (i)

At least $n-r$ values in $\lambda^{(t)}_{1},\ldots,\lambda^{(t)}_{n}$ are
zeros, or equivalently $\text{rank}(X^{(t)})\leq r$. This follows from
conditions $X^{(t)}Y^{(t)}=0$ and $\text{rank}(Y^{(t)})=n-r$, which imply that
$\text{rank}(X^{(t)})\leq n-\text{rank}(Y^{(t)})\leq n-(n-r)=r$.

  * (ii)

If $\lambda^{(t)}_{n}\geq 0$ then
$\lambda^{(t)}_{r+1}=\ldots=\lambda^{(t)}_{n}=0$ and $\lambda^{(t)}_{r}>0$.
The first claim follows from (ii); if $\lambda^{(t)}_{r}=0$ then
$\text{rank}(X^{(t)})\leq r-1$ contradicting the minimality of $r$ (note that
$X^{(t)}\succeq 0$ and hence $X^{(t)}\in P$).

  * (iii)

For each $i\in[n]$, $\lambda^{(t)}_{i}$ is a continuous function of $t$. This
holds by continuity of the curve and by (2b).

Conditions $X^{(0)}\in P$ and $X^{(1)}\notin P$ give $\lambda_{n}^{(0)}\geq 0$
and $\lambda_{n}^{(1)}<0$. Denote
$s=\sup\\{t\in[0,1]\>:\>\lambda_{n}^{(t)}\geq 0\\}$, then
$\lambda_{n}^{(s)}\geq 0$ by continuity and hence $s<1$. By (ii) we have
$\mu{\stackrel{{\scriptstyle\mbox{\tiny\tt~{}def~{}}}}{{=}}}\lambda_{r}^{(s)}>0$.
By continuity, there exists $u\in(s,1]$ such that $\lambda_{r}^{(t)}>\mu/2$
for all $t\in[s,u]$. Condition (i) thus gives
$\lambda_{r+1}^{(t)}=\ldots=\lambda_{n}^{(t)}=0$ for all $t\in[s,u]$, and thus
$\sup\\{t\in[0,1]\>:\>\lambda_{n}^{(t)}\geq 0\\}\geq u>s$, which is a
contradiction.

(b)  Let $C^{\prime}$ be the connected component of ${\mathscr{V}}_{r,\iota}$
containing $(X^{\ast},Y^{\ast})$. Clearly, we have
$\pi_{\mathbb{S}^{n}}(C^{\prime})\subseteq C$, and so the claim holds by part
(a). ∎

##  3 Incidence varieties and facial reduction

We consider the feasibility problem of the primal SDP in (3), that is the
system:

$\mathscr{A}(X)=b,\,\,\,\,\,\,\,\,X\succeq 0.$ (7)

Let $P=\\{X\in\mathbb{S}^{n}\,\,:\,\,\mathscr{A}(X)=b,X\succeq 0\\}$ be the
set of feasible solutions of (4), and let $r$ be the maximum rank of a matrix
in $P$. From now on, we assume that $P\neq\varnothing$ and $r>0$. A matrix
$X\in P$ of rank $r$ will be called a maximum-rank solution.

Below we recall some basic facts about facial reduction for SDPs, see e.g.
[11]. Let $\mathcal{F}$ be the minimal face of the SDP cone
$\mathbb{S}^{n}_{+}$ containing $P$. It can be represented by a full-rank
matrix $U\in\mathbb{R}^{n\times r}$ as follows:

$\mathcal{F}=\\{X\in\mathbb{S}^{n}_{+}\>:\>\text{range}(X)\subseteq\text{range}(U)\\}=\\{UZU^{\top}\>:\>Z\in\mathbb{S}^{r}_{+}\\}.$
(8)

In other words, $\mathcal{F}$ is linearly isomorphic (through $UZU^{\top}\to
Z$) to a copy of the cone $\mathbb{S}^{r}_{+}$, and in particular it has
dimension $\binom{r+1}{2}$. Any matrix in the relative interior of $P$ is in
the relative interior of $\mathcal{F}$ and hence has rank $r$; every other
matrix in $\mathcal{F}$ has rank $\leq r$.

Let $\iota\subseteq[n]$ be a maximal set of linearly independent rows of $U$,
with $|\iota|=r$. We assume for notational convenience that $\iota=[r]$ (this
can be achieved by renaming rows and columns of $X$). Since $U$ is full-rank,
we can apply Gaussian elimination to compute the column reduced echelon form
of $U$, without changing its range. Accordingly, we can assume w.l.o.g. that

$U=\begin{pmatrix}I_{r}\\\ U_{0}\end{pmatrix},\qquad\text{with
}U_{0}\in\mathbb{R}^{n-r\times r}.$ (9)

Then for the matrix $Y_{U}:=-{U_{0}}^{\top}$ we have

$\displaystyle\mathcal{F}=\\{X\in\mathbb{S}^{n}_{+}\,\,:\,\,\text{range}(X)\subseteq\text{range}(U)\\}$
$\displaystyle=\\{X\in\mathbb{S}^{n}_{+}\,\,:\,\,\left(\begin{smallmatrix}Y_{U}^{\top}&I_{n-r}\end{smallmatrix}\right)X=0^{(n-r)\times
n}\\}$ (10a)
$\displaystyle=\\{X\in\mathbb{S}^{n}_{+}\,\,:\,\,X\left(\begin{smallmatrix}Y_{U}\\\
I_{n-r}\end{smallmatrix}\right)=0^{n\times(n-r)}\\}$ (10b)

Now consider the following system of equations (in variables
$X\in\mathbb{S}^{n}_{\mathbb{C}}$ and $Y\in\mathbb{C}^{n\times(n-r)}$):

$\displaystyle\mathscr{A}(X)$ $\displaystyle=b$ (11a) $\displaystyle
X\left(\begin{matrix}Y\\\ I_{n-r}\end{matrix}\right)$
$\displaystyle=0^{n\times(n-r)}$ (11b)

Let
${\mathscr{V}}_{r,\iota}\subseteq\mathbb{S}^{n}_{\mathbb{C}}\times\mathbb{C}^{n\times(n-r)}$
be the set of complex solutions of (11). Note that it is essentially the same
incidence variety defined in the previous section (for $\iota=[r]$), modulo
notation: matrix $Y$ in eq. (6) has a larger size, but all additional entries
are fixed to constants. We define ${\mathscr{W}}_{r,\iota}$,
$\mathbb{R}{\mathscr{V}}_{r,\iota}$ and $\mathbb{R}{\mathscr{W}}_{r,\iota}$ in
the same way as before:

$\displaystyle{\mathscr{W}}_{r,\iota}=\pi_{\mathbb{S}^{n}_{\mathbb{C}}}({\mathscr{V}}_{r,\iota})$
$\displaystyle=\\{X\in\mathbb{S}^{n}_{\mathbb{C}}\,\,:\,\,\exists
Y\in\mathbb{C}^{r\times(n-r)}\;\mbox{s.t.}\;(X,Y)\in{\mathscr{V}}_{r,\iota}\\}$
(12a) $\displaystyle=\\{X\in\mathbb{S}^{n}_{\mathbb{C}}\,\,:\,\,\exists
Y\in\mathbb{C}^{r\times(n-r)}\;\mbox{s.t.}\;\mathscr{A}(X)=b,X\left(\begin{smallmatrix}Y\\\
I_{n-r}\end{smallmatrix}\right)=0^{n\times(n-r)}\\}$ (12b)
$\mathbb{R}{\mathscr{V}}_{r,\iota}=(\mathbb{S}^{n}\times\mathbb{R}^{n\times(n-r)})\cap{\mathscr{V}}_{r,\iota}\hskip
40.0pt\mathbb{R}{\mathscr{W}}_{r,\iota}=\mathbb{S}^{n}\cap{\mathscr{W}}_{r,\iota}=\pi_{\mathbb{S}^{n}}(\mathbb{R}{\mathscr{V}}_{r,\iota})$
(13)

By construction,
$P\subseteq\mathcal{F}\cap\\{X\in\mathbb{S}^{n}\,\,:\,\,\mathscr{A}(X)=b\\}\subseteq\mathbb{R}{\mathscr{W}}_{r,\iota}$.

######  Lemma 3.

Let $X^{*}=UZ^{*}U^{\top}\in P$ be a maximum-rank solution, with $U$ as in (9)
and $Z^{*}\in\mathbb{S}^{r}_{++}$ positive definite. Define open neighborhood
$\Omega^{*}$ of $X^{\ast}$ via
$\Omega^{\ast}=\\{X\in\mathbb{S}^{n}\>:\>||X-X^{\ast}||<\rho(X^{\ast})\\}$.
The following holds:

  * (a)

$P\cap\Omega^{*}=\mathbb{R}{\mathscr{W}}_{r,\iota}\cap\Omega^{*}$.

  * (b)

For every $X\in\mathbb{R}{\mathscr{W}}_{r,\iota}\cap\Omega^{*}$ there exists a
unique matrix $Y$ such that $(X,Y)\in{\mathscr{V}}_{r,\iota}$. This matrix is
$Y=Y_{U}$.

###### Proof.

Part (a). If $X\in P$ then $\mathscr{A}(X)=b$ and
$(X,Y_{U})\in{\mathscr{V}}_{r,\iota}$ (since $P\subseteq{\mathcal{F}}$), and
thus $X\in\mathbb{R}{\mathscr{W}}_{r,\iota}$. This shows that
$P\subseteq\mathbb{R}{\mathscr{W}}_{r,\iota}$ and hence
$P\cap\Omega^{*}\subseteq\mathbb{R}{\mathscr{W}}_{r,\iota}\cap\Omega^{*}$.

Now suppose that $X\in\mathbb{R}{\mathscr{W}}_{r,\iota}\cap\Omega^{*}$. Then
there exists $Y$ such that $(X,Y)\in{\mathscr{V}}_{r,\iota}$, i.e. $(X,Y)$
satisfies system (11). Eq. (11b) implies that $X$ has at least $n-r$ zero
eigenvalues. Since $\text{rank}(X^{\ast})=r$, we have
$\lambda_{r}(X^{\ast})=\rho(X^{\ast})>0$. By (2b), for all $X\in\Omega^{\ast}$
we have
$|\lambda_{r}(X)-\lambda_{r}(X^{\ast})|\leq||X-X^{\ast}||<\rho(X^{\ast})=\lambda_{r}(X^{\ast})$
and thus $\lambda_{r}(X)>0$, i.e. $X$ has at least $r$ strictly positive
eigenvalues. We conclude that $X$ has exactly $n-r$ zero eigenvalues and
exactly $r$ strictly positive eigenvalues, therefore $X\succeq 0$ and
$\text{rank}(X)=r$. Combined with the fact that $\mathscr{A}(X)=b$, this
implies that $X\in P$. We showed that converse inclusion
$\mathbb{R}{\mathscr{W}}_{r,\iota}\cap\Omega^{*}\subseteq P\cap\Omega^{*}$ is
also true.

Part (b). As shown in (a), we have $X\in P\subseteq\mathcal{F}$. Therefore,
$(X,Y_{U})\in{\mathscr{V}}_{r,\iota}$, i.e. it satisfies (11). Suppose that
there exists $Y\neq Y^{\prime}$ such that $(X,Y)$ and $(X,Y^{\prime})$ satisfy
(11). Denote $C=\left(\begin{smallmatrix}Y&Y^{\prime}\\\
I_{n-r}&I_{n-r}\end{smallmatrix}\right)$, then $XC=0^{n\times 2(n-r)}$.
Clearly, we have $\text{rank}(C)\geq n-r+1$: indeed, the first $n-r$ columns
of $C$ are linearly independent, and there exists a column among the remaining
ones which is not in the range of the first $n-r$ columns, by the assumption
$Y\neq Y^{\prime}$. This implies that $\text{rank}(X)\leq n-\text{rank}(C)\leq
n-(n-r+1)=r-1$, which contradicts condition $\text{rank}(X)=r$ shown in part
(a).

∎

##  4 Description of the method

Motivated by Lemma 3, we propose the following hybrid method for certifying
the feasibility of system (7) (see Algorithm 1).

1 Solve system (7) numerically to get solution $\tilde{X}\in\mathbb{S}^{n}$
which is close to some maximum-rank solution $X^{\ast}$.

2 Find numerically a maximal set $\iota\subseteq[n]$ of linearly independent
columns of $\tilde{X}$, let $r=|\iota|$. Assume for notational convenience
that $\iota=[r]$ (this can be achieved by renaming variables). Write
$\tilde{X}=\left(\begin{smallmatrix}\tilde{S}&\tilde{R}^{\top}\\\
\tilde{R}&\tilde{W}\end{smallmatrix}\right)$ where
$\tilde{S}=\tilde{X}^{\iota}_{\iota}\in\mathbb{S}^{r}$. If $\tilde{S}$ is
singular then terminate with failure, otherwise define
$\tilde{Y}=-\tilde{S}^{-1}\tilde{R}^{\top}\in\mathbb{R}^{r\times(n-r)}$. As we
show later, $\tilde{X}\left(\begin{smallmatrix}\tilde{Y}\\\
I_{n-r}\end{smallmatrix}\right)\approx 0^{n\times(n-r)}$.

3 Write system (11) as  $\mathcal{Q}(Y)X_{\text{hvec}}=q$ (14) where
$X_{\text{hvec}}$ is a vector of dimension $k=\frac{{n(n+1)}}{2}$ obtained by
vectorizing the lower triangular part of $X$, $\mathcal{Q}(Y)$ is a matrix
that depends linearly on $Y$, and $q$ is fixed vector.

4 Find numerically a maximal set $J\subseteq[k]$ of linearly independent
columns of $\mathcal{Q}(\tilde{Y})$.

Augment system (11) by fixing variables in $J^{\prime}:=[n]-J$:
$\displaystyle\mathscr{A}(X)$ $\displaystyle=b$ (15a) $\displaystyle
X\left(\begin{matrix}Y\\\ I_{n-r}\end{matrix}\right)$
$\displaystyle=0^{n\times(n-r)}$ (15b) $\displaystyle(X_{\text{hvec}})_{j}$
$\displaystyle=(\tilde{X}_{\text{hvec}})_{j}\qquad\forall j\in J^{\prime}$
(15c) Compute one point per connected component of the real algebraic variety
defined by Equation 15 as in [22] (cf. Section 4.2). By [22, Lemma 3.2], the
equations in (15b) at positions $i,j$ with $i-j>r$ are implied by other
equations and thus can be omitted.

Algorithm 1 Hybrid algorithm for system $\mathscr{A}(X)=b,\,\,\,X\succeq 0$.

The next observation justifies Step 2.

######  Lemma 4.

Consider a PSD matrix $X=\left(\begin{smallmatrix}S&R^{\top}\\\
R&W\end{smallmatrix}\right)\in\mathbb{S}^{n}_{+}$ with $S\in\mathbb{S}^{r}$,
and suppose that $\iota=[r]$ is a maximal set of linearly independent columns
(so that the rank of $X$ is $r$). Then $S$ is positive definite, and matrix
$Y=-S^{-1}R^{\top}$ satisfies $X\left(\begin{smallmatrix}Y\\\
I_{n-r}\end{smallmatrix}\right)=0^{n\times(n-r)}$.

###### Proof.

The first claim follows from Lemma 8 given later. Since the last $n-r$ columns
of $X$ linearly depend on the first $r$ columns, we must have

$\begin{pmatrix}R^{\top}\\\ W\end{pmatrix}=\begin{pmatrix}S\\\
R\end{pmatrix}Z$

for some matrix $Z\in\mathbb{R}^{r\times(n-r)}$, i.e. $R^{\top}=SZ$ and
$W=RZ$. Hence $Z=S^{-1}R^{\top}$ and $W=RS^{-1}R^{\top}$. Using these
equations, one can now check that $SY+R^{\top}=0^{r\times(n-r)}$ and
$RY+W=0^{(n-r)\times(n-r)}$, which yields the second claim. ∎

######  Definition 5.

We say that sets $\iota\subseteq[n]$ and $J\subseteq[k]$ are valid for a
maximum-rank solution $X^{\ast}$ if they can be outputs of steps 2 and 4,
respectively, assuming that step 1 computes the vector $X^{\ast}$, and the
computations in steps 2 and 4 are performed exactly (i.e. a maximal set of
linearly independent columns is computed exactly rather than numerically).

######  Theorem 6.

Each maximum-rank solution $X^{\ast}$ has an open neighborhood
$\tilde{\Omega}\subset\mathbb{S}^{n}$ with the following property. Suppose
that step 1 produces point $\tilde{X}\in\tilde{\Omega}$, and lines 2 and 4
output sets $\iota$ and $J$, respectively, which are valid for $X^{\ast}$.
Then the set of real solutions of system (15) has a zero-dimensional component
with a solution $(X,Y)$ where $X\in\mathbb{S}^{n}_{+}$.

###### Proof.

By assumption, set $\iota$ is valid for $X^{\ast}$, i.e. $\iota$ is a maximal
set of linearly independent columns in $X^{\ast}$ (and so $|\iota|=r$). Assume
for notational convenience that $\iota=[r]$. By Lemma 4, $\iota$ is also a
maximal set of linearly independent rows in $X^{\ast}$. Choose matrix
$U\in\mathbb{R}^{n\times r}$ characterizing face ${\mathcal{F}}$ as in (8).
Since $\text{range}(X^{\ast})\subseteq\text{range}(U)$ and
$\text{rank}(X^{\ast})=r$, set $\iota$ is a maximal set of linearly
independent rows in $U$. Thus, we can choose $U$ as in (9).

Denote
$\Omega^{\ast}=\\{X\in\mathbb{S}^{n}\>:\>||X-X^{\ast}||<\rho(X^{\ast})\\}$. By
Lemma 3(a,b), for any $(X,Y)$ we have

$\begin{array}[]{c}\mathcal{Q}(Y)X_{\text{hvec}}=q\\\
X\in\Omega^{\ast}\end{array}\;\;\Longleftrightarrow\;\;\begin{array}[]{c}\mathcal{Q}(Y_{U})X_{\text{hvec}}=q\\\
X\in\Omega^{\ast}\\\
Y=Y_{U}\end{array}\;\;\Longrightarrow\;\;X\in\mathbb{S}^{n}_{+}$ (16)

To simplify notation, let us write
$X_{\text{hvec}}=:\left(\begin{smallmatrix}X^{J}\\\
X^{J^{\prime}}\end{smallmatrix}\right)$ and
$\mathcal{Q}(Y_{U})=:\mathcal{Q}=:\begin{pmatrix}Q^{J}&Q^{J^{\prime}}\end{pmatrix}$.
By assumption, $J$ is a maximal set of linearly independent columns for matrix
$\mathcal{Q}$. Let $I$ be a maximal set of linearly independent rows of
$\mathcal{Q}$, then $|I|=|J|=\text{rank}(\mathcal{Q})$. The linear system
$\mathcal{Q}X_{\text{hvec}}=q$ has at least one solution, namely
$X_{\text{hvec}}=X^{\ast}_{\text{hvec}}$. Thus, removing equations
corresponding to rows $i\notin I$ will not affect the set of feasible
solutions. This implies that system
$\mathcal{Q}X_{\text{hvec}}=q=\mathcal{Q}X^{\ast}_{\text{hvec}}$ is equivalent
to the system

$\begin{pmatrix}Q^{J}_{I}&Q^{J^{\prime}}_{I}\end{pmatrix}\begin{pmatrix}X^{J}-(X^{\ast})^{J}\\\
X^{J^{\prime}}-(X^{\ast})^{J^{\prime}}\end{pmatrix}=0,$

which is in turn equivalent to
$X^{J}-(X^{\ast})^{J}=-(\mathcal{Q}^{J}_{I})^{-1}\mathcal{Q}^{J^{\prime}}_{I}(X^{J^{\prime}}-(X^{\ast})^{J^{\prime}})$
since matrix $\mathcal{Q}^{J}_{I}$ is nonsingular. Let us denote
$\mathcal{R}^{J}_{I}=-(\mathcal{Q}^{J}_{I})^{-1}\mathcal{Q}^{J^{\prime}}_{I}$.

By adding constraint $X^{J^{\prime}}=\tilde{X}^{J^{\prime}}$ to (16) we can
now conclude the following for any $(X,Y)$:

$\begin{array}[]{c}(X,Y)\mbox{ satisfies~{}\eqref{eq:fixed-variables}}\\\
X\in\Omega^{\ast}\end{array}\;\;\Longleftrightarrow\;\;\begin{array}[]{c}X^{J}-(X^{\ast})^{J}=\mathcal{R}^{J}_{I}(\tilde{X}^{J^{\prime}}-(X^{\ast})^{J^{\prime}})\\\
X^{J^{\prime}}=\tilde{X}^{J^{\prime}}\\\ X\in\Omega^{\ast}\\\
Y=Y_{U}\end{array}\;\;\Longrightarrow\;\;X\in\mathbb{S}^{n}_{+}$ (17)

Let us define
$\delta_{\iota,J,I}=\frac{\rho(X^{\ast})}{4\max\\{||\mathcal{R}^{J}_{I}||,1\\}}$.
We claim that for any $\tilde{X}\in\mathbb{S}^{n}$ with
$||\tilde{X}-X^{\ast}||_{F}<\delta_{\iota,J,I}$ the system in the middle of
(17) has exactly one feasible solution $(X,Y)$. Indeed, we only need to verify
that $X$ defined by the equalities in this system satisfies
$X\in\Omega^{\ast}$. This holds since
$||X^{J}-(X^{\ast})^{J}||\leq||\mathcal{R}^{J}_{I}||\cdot||\tilde{X}^{J^{\prime}}-(X^{\ast})^{J^{\prime}}||<\frac{1}{4}\rho(X^{\ast})$
and $||X^{J^{\prime}}-(X^{\ast})^{J^{\prime}}||<\frac{1}{4}\rho(X^{\ast})$,
implying $||X-X^{\ast}||\leq||X-X^{\ast}||_{F}\leq
2||X^{J}-(X^{\ast})^{J}||+2||X^{J^{\prime}}-(X^{\ast})^{J^{\prime}}||<\rho(X^{\ast})$.

We can now define set $\tilde{\Omega}$ in Theorem 6 as follows:
$\tilde{\Omega}=\\{\tilde{X}\in\mathbb{S}^{n}\>:\>||\tilde{X}-X^{\ast}||_{F}<\delta\\}$
where $\delta=\min_{\iota,J,I}\delta_{\iota,J,I}$ and the minimum is taken
over (finitely many) valid choices of $\iota,J,I$. ∎

###  4.1 Numerical algorithms: implementing steps 2 and 4

To implement Algorithm 1, we need a procedure for numerically computing a
maximal set of linearly independent columns of a given matrix. This can be
done, for example, by a rank-revealing Gaussian elimination algorithm [35,
44]. Such algorithm takes matrix $A\in\mathbb{R}^{p\times q}$ and tolerance
value $\varepsilon$ as an input, and produces an integer $r$ and subset
$J\subseteq[q]$ of size $r$ with the following guarantees: 111The algorithm
actually produces subsets $I\subseteq[p]$ and $J\subseteq[q]$ of size $r$ such
that $\sigma_{r}(A_{I}^{J})\geq\sigma_{r}(A)/c_{pq}$. (This follows by
combining Lemma 3.1 and Theorem 4.1 in [44].) This implies (18c) since
$\sigma_{r}(A^{J})\geq\sigma_{r}(A^{J}_{I})$ by a well-known property of
singular values, see e.g. [14, Corollary 8.6.3].

$\displaystyle\sigma_{r}(A)$ $\displaystyle\geq$ $\displaystyle\varepsilon$
(18a) $\displaystyle\sigma_{r+1}(A)$ $\displaystyle\leq$ $\displaystyle
c_{pq}\varepsilon$ (18b) $\displaystyle\sigma_{r}(A^{J})$ $\displaystyle\geq$
$\displaystyle\sigma_{r}(A)/c_{pq}$ (18c)

where constant $c_{pq}>1$ depends polynomially on dimensions $p,q$. Below we
analyze Algorithm 1 assuming that steps 2 and 4 use the method above with
tolerance values $\varepsilon_{1}$ and $\varepsilon_{2}$, respectively. The
next result shows that if $||\tilde{X}-X^{\ast}||$ is sufficiently small then
from $\tilde{X}$ it is possible to compute sets $\iota,J$ which are valid for
$X^{\ast}$.

######  Theorem 7.

(a) Suppose that $\delta<\varepsilon_{1}<(\rho^{\ast}-\delta)/c_{nn}$ where
$\delta=||\tilde{X}-X^{\ast}||$ and $\rho^{\ast}=\rho(X^{\ast})$. Then set
$\iota$ computed in step 2 is a maximal set of linearly independent columns
for $X^{\ast}$. Furthermore,

$\lambda_{r}(\tilde{S})\;\;\geq\;\;\phi(\delta)\;\;:=\;\;\frac{((\rho^{\ast}-\delta)/c_{nn}-\delta)^{2}}{n||X^{\ast}||}-\delta$
(19)

(b) Suppose the precondition in (a) holds and $\phi(\delta)>\delta$. Write
$X^{\ast}=\begin{pmatrix}S^{\ast}&(R^{\ast})^{\top}\\\
R^{\ast}&W^{\ast}\end{pmatrix}$ where $S^{\ast}=(X^{\ast})^{\iota}_{\iota}$,
and let $Y^{\ast}=-(S^{\ast})^{-1}(R^{\ast})^{\top}$. Let $||\mathcal{Q}||$ be
the 2-norm of linear operator
$\mathcal{Q}:\mathbb{R}^{r\times(n-r)}\rightarrow\mathbb{R}^{p\times q}$, and
denote

$\psi(\delta)=||\mathcal{Q}||\cdot\left[\frac{||X^{\ast}||}{\phi(\delta)-\delta}+1\right]\cdot\frac{\delta}{\phi(\delta)}$
(20)

If
$\psi(\delta)<\varepsilon_{2}<(\rho(\mathcal{Q}(Y^{\ast}))-\psi(\delta))/c_{pq}$
then set $J$ computed in step 4 is a maximal set of linearly independent
columns for $\mathcal{Q}(Y^{\ast})$.

Note that $\lim_{\delta\rightarrow 0}\phi(\delta)=const>0$ and
$\lim_{\delta\rightarrow 0}\psi(\delta)=0$. Thus, the preconditions in (a,b)
can be satisfied for sufficiently small $\delta>0$ as long as
$\rho(\mathcal{Q}(Y^{\ast}))>0$. (The case when
$\rho(\mathcal{Q}(Y^{\ast}))=0$ is not interesting: we then have
$\mathcal{Q}(Y^{\ast})=0$ and thus ${\mathcal{A}}(X)=0$ for all $X$). We
remark that linear operator $\mathcal{Q}$, function $\psi(\cdot)$ and value
$\rho(\mathcal{Q}(Y^{\ast}))-\psi(\delta)$ depend on $\iota$, but there are
only finitely many subsets $\iota\subseteq[n]$ which give maximal linearly
independent columns for $X^{\ast}$, so we could take the minimum value over
such $\iota$’s when formulating the final condition.

In the remainder of this section we prove Theorem 7. We start with part (a).
Let $r^{\ast}=\text{rank}(X^{\ast})$ be the true rank and $r$ be the rank
computed by the numerical procedure in step 2. Condition (18c) gives
$\sigma_{r}(\tilde{X}^{\iota})\geq\sigma_{r}(\tilde{X})/c_{nn}\geq(\sigma_{r}(X^{\ast})-||\tilde{X}-X^{\ast}||)/c_{nn}=(\rho^{\ast}-\delta)/c_{nn}$.
Therefore,
$\sigma_{r}((X^{\ast})^{\iota})\geq\sigma_{r}(\tilde{X}^{\iota})-||X^{\ast}-\tilde{X}||\geq(\rho^{\ast}-\delta)/c_{nn}-\delta>0$,
implying that columns in $\iota$ are linearly independent in $X^{\ast}$ (and
so $r\leq r^{\ast}$). If $r<r^{\ast}$ then $\sigma_{r+1}(\tilde{X})\leq
c_{nn}\varepsilon_{1}$ by condition (18b) and
$\sigma_{r+1}(\tilde{X})\geq\sigma_{r^{\ast}}(\tilde{X})\geq\sigma_{r^{\ast}}(X^{\ast})-||\tilde{X}-X^{\ast}||=\rho^{\ast}-\delta$
by (2a), which is a contradiction. This shows that $r=r^{\ast}$.

The last claim (19) will follow from the result below for matrix $X=\tilde{X}$
and value
$\tau=\sigma_{r}(\tilde{X}^{\iota})\geq(\rho^{\ast}-\delta)/c_{nn}>\delta$.

######  Lemma 8.

Consider symmetric matrices $X^{\ast}\in\mathbb{S}^{n}_{+}$,
$X\in\mathbb{S}^{n}$ with $||X^{\ast}-X||=\delta$. Let $\iota\subseteq[n]$ be
a subset of size $r$ with $\sigma_{r}(X^{\iota})=\tau\geq\delta$. Then
$\lambda_{r}(X^{\iota}_{\iota})\geq\frac{(\tau-\delta)^{2}}{n||X^{\ast}||}-\delta$.

Note that this lemma also implies the first part of Lemma 4 (by setting
$X=X^{\ast}$).

###### Proof.

We will use a well-known fact that for any matrix $A\in\mathbb{R}^{p\times q}$
with $p\geq q$ we have
$\sigma_{q}(A)=\min\\{||Au||\>:\>u\in\mathbb{R}^{q},||u||=1\\}$.

Consider vector $u\in\mathbb{R}^{\iota}$ with $||u||=1$, then
$||X^{\iota}u||\geq\tau$. Let $v\in\mathbb{R}^{n}$ be the vector with
$v_{i}=u_{i}$ for $i\in\iota$ and $v_{i}=0$ for $i\in[n]-\iota$. We have
$||Xv||=||X^{\iota}u||\geq\tau$, and hence

$||X^{\ast}v||\geq||Xv||-||(X^{\ast}-X)v||\geq\tau-||X^{\ast}-X||\cdot
1=\tau-\delta$

Let us write
$X^{\ast}=\sum_{i=1}^{n}\lambda_{i}v_{i}v_{i}^{\top}=\sum_{i=1}^{n}w_{i}w_{i}^{\top}$
where $||X^{\ast}||=\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{n}\geq
0$, $v_{1},\ldots,v_{n}$ are orthonormal vectors, and
$w_{i}=\sqrt{\lambda_{i}}v_{i}$ (then $||w_{i}||\leq\sqrt{\lambda_{1}}$).
Denote $\alpha_{i}=|v^{\top}w_{i}|$, then
$v^{\top}X^{\ast}v=\sum_{i=1}^{n}\alpha_{i}^{2}$. We can write

$\displaystyle\tau-\delta$
$\displaystyle\leq||X^{\ast}v||=||\sum_{i=1}^{n}w_{i}(v^{\top}w_{i})||\leq\sum_{i=1}^{n}\alpha_{i}||w_{i}||$
$\displaystyle\leq\sqrt{\lambda_{1}}\cdot\sum_{i=1}^{n}\alpha_{i}\leq\sqrt{\lambda_{1}}\cdot\sqrt{n\sum_{i=1}^{n}\alpha_{i}^{2}}=\sqrt{\lambda_{1}n(v^{\top}X^{\ast}v)}$

where the last inequality is tight if and only if
$\alpha_{1}=\ldots=\alpha_{n}$. Therefore,

$u^{\top}X^{\iota}_{\iota}u=v^{\top}Xv=v^{\top}X^{\ast}v+v^{\top}(X-X^{\ast})v\geq\frac{(\tau-\delta)^{2}}{n\lambda_{1}}-||v||^{2}\cdot||X-X^{\ast}||=\frac{(\tau-\delta)^{2}}{n\lambda_{1}}-\delta$

∎

Next, we proof part (b) of Theorem 7. We will need the following result.

######  Theorem 9 ([1, Theorem 7.12]).

Let $A,B$ be square matrices of the same size such that $A$ is nonsingular and
$||A-B||\leq\frac{1}{||A^{-1}||}$. Then $B$ is also nonsingular and

$||B^{-1}||\leq\frac{||A^{-1}||}{1-||A^{-1}||\cdot||A-B||}$
$||A^{-1}-B^{-1}||\leq\frac{||A^{-1}||^{2}\cdot||A-B||}{1-||A^{-1}||\cdot||A-B||}$

For brevity, let us denote $S=S^{\ast}$ and $R=R^{\ast}$. We can write

$\tilde{Y}-Y=S^{-1}R^{\top}-\tilde{S}^{-1}\tilde{R}^{\top}=(S^{-1}-\tilde{S}^{-1})R^{\top}+\tilde{S}^{-1}(R^{\top}-\tilde{R}^{\top})$

From part (a) we know that $||\tilde{S}^{-1}||\leq\frac{1}{\phi(\delta)}$, and
so
$||\tilde{S}^{-1}-S||\leq||\tilde{X}-X^{\ast}||=\delta<\phi(\delta)\leq\frac{1}{||\tilde{S}^{-1}||}$.
Using Theorem 9 with $A=\tilde{S}$, $B=S$ we get

$\displaystyle||\tilde{Y}-Y^{\ast}||$
$\displaystyle\leq||S^{-1}-\tilde{S}^{-1}||\cdot||R^{\top}||+||\tilde{S}^{-1}||\cdot||R^{\top}-\tilde{R}^{\top}||$
$\displaystyle\leq\frac{(\frac{1}{\phi(\delta)})^{2}\cdot\delta}{1-\frac{1}{\phi(\delta)}\cdot\delta}\cdot||X^{\ast}||+\frac{1}{\phi(\delta)}\cdot\delta=\left[\frac{||X^{\ast}||}{\phi(\delta)-\delta}+1\right]\cdot\frac{\delta}{\phi(\delta)}$

Denote $\mathcal{Q}^{\ast}=\mathcal{Q}(X^{\ast})$ and
$\tilde{\mathcal{Q}}=\mathcal{Q}(\tilde{X})$, then
$||\tilde{\mathcal{Q}}-\mathcal{Q}^{\ast}||\leq||\mathcal{Q}||\cdot||\tilde{Y}-Y^{\ast}||\leq\psi(\delta)$.
We now conclude that $J$ is a maximal set of linearly independent columns of
$\mathcal{Q}^{\ast}$ by the same argument as in part (a).

###  4.2 Symbolic algorithms: implementing step 5

Step 5 of Algorithm 1 and Theorem 6, similarly to [22], reduce the original
problem of validating the feasibility of a LMI to a classical question in real
algebraic geometry, the computation of one point per connected component of
the real trace $\mathbb{R}V=V\cap\mathbb{R}^{n}$ of a complex variety

$V=\\{x\in\mathbb{C}^{n}:f_{1}(x)=0,\ldots,f_{c}(x)=0\\},$

see e.g. [3, § 12.6]. A standard technique consists in computing the critical
points of a well-chosen polynomial map $\varphi:\mathbb{C}^{n}\to\mathbb{C}$
(usually of low degree) restricted to $V$, often called critical point method
(CPM). The expected output of a CPM is a finite set $T\subset\mathbb{R}^{n}$
such that $T$ intersects every connected component of $\mathbb{R}V$.

Usual choices for $\varphi$ are linear functions: when $\varphi$ is a generic
linear function, and since $\mathbb{R}V$ has finitely many connected
components, there are finitely many critical points of its restriction to
$\mathbb{R}V$. In case $\mathbb{R}V$ is non-compact, $\varphi$ might not admit
critical points on the (unbounded) connected components $C\subseteq V$
satistfying $\varphi(C)=\mathbb{R}$. In this case, the recursive method in
[43] allows to compute such components. However, our situation is special
inasmuch as our target component is an isolated real point, so that the
mentioned recursion is not needed on Step 5 of Algorithm 1.

The CPM method for polynomials $f:C^{n}\to\mathbb{C}^{c}$ and polynomial map
$\varphi:\mathbb{C}^{n}\to\mathbb{C}$ can be described as follows. Recall that
$J(f)=({\partial f_{i}}/{\partial x_{j}})_{ij}$ is the Jacobian matrix of $f$.
Define the extended Jacobian of $(\varphi,f)$ via

$J(\varphi,f):=\begin{pmatrix}\frac{\partial\varphi}{\partial
x_{1}}&\cdots&\frac{\partial\varphi}{\partial x_{n}}\\\ &J(f)&\end{pmatrix}.$

Now consider the following Lagrange system associated with $\varphi$ and $V$:

$\mathscr{L}_{u}:\begin{cases}f_{1}=\cdots=f_{c}=0\\\ z^{T}J(\varphi,f)=0\\\
z^{T}u-1=0\end{cases}$ (21)

where $z=(z_{0},z_{1},\ldots,z_{c})\in\mathbb{C}^{c+1}$ are new variables and
$u\in C^{c+1}$ is a fixed vector. Intuitively, constraint $z^{T}u=1$ for a
generic $u$ is essentially equivalent to the constraint $z\neq 0$.

Let $V(\mathscr{L}_{u})\subseteq\mathbb{C}^{n+c+1}$ be the set of solutions of
system (21), and let $\Pi_{x}(x,z):\mathbb{C}^{n+c+1}\to\mathbb{C}^{n}$ be the
projection to the first $n$ variables (with $\Pi_{x}(x,z)=x$). Under certain
conditions (e.g. those described below) the set $\Pi_{x}(V(\mathscr{L}_{u}))$
will be finite. In this case this set can be computed, for example, via the
software msolve based on Gröbner basis computations, and represented through a
so-called rational univariate representation. This is given by a sequence

$(q,q_{0},q_{1},\ldots,q_{n})\in\mathbb{Q}[t]^{n+2},$

with $q_{0},q$ coprime, and such that

$T=\left\\{\left(\frac{q_{1}(t)}{q_{0}(t)},\ldots,\frac{q_{n}(t)}{q_{0}(t)}\right)\in\mathbb{R}^{n}:q(t)=0\right\\}.$

In other words, the coordinates of vectors $x\in\Pi_{x}(V(\mathscr{L}_{u}))$
are represented by the evaluation of $n$ univariate rational functions to the
roots of a univariate polynomial $q$.

Properties of CPM  With some abuse of notation, let us denote

$\displaystyle\mathrm{Reg}(f)=\\{x\in
V\>:\>\text{rank}\,J(f)=c\\}\qquad\mathrm{Sing}(f)=\\{x\in
V\>:\>\text{rank}\,J(f)<c\\}$

Note that $\mathrm{Reg}(f)=\mathrm{Reg}(V)$ and
$\mathrm{Sing}(f)=\mathrm{Sing}(V)$ (respectively the set of regular and
singular points of $V=V(f)$), assuming that the following condition holds.

######  Assumption 1.

$V$ is equidimensional of codimension $c$ and $I(V)=\langle
f_{1},\ldots,f_{c}\rangle$. (Then $V$ is a complete intersection of dimension
$d=n-c$, as defined in Section 2).

A critical point of the restriction of $\varphi$ to $\mathrm{Reg}(f)$ is a
point $x\in\mathrm{Reg}(f)$ such that the differential of the restriction of
$\varphi$ to $\mathrm{Reg}(f)$ at $x$ is not surjective, that is an element of
the constructible set

$\displaystyle\mathrm{Crit}(\varphi,f)$
$\displaystyle={\\{x\in\mathrm{Reg}(f):\,\text{rank}\,J(\varphi,f)<c+1\\}}$
$\displaystyle={\\{x\in\mathrm{Reg}(f):\,\text{rank}\,J(\varphi,f)=c\\}}$

We also denote
$\mathrm{Crit}(\varphi,V)=\\{x\in\mathrm{Reg}(V):\,\text{rank}\,J(\varphi,f)=c\\}$
under 1.

######  Proposition 10.

Define

$\displaystyle\mathrm{Crit}_{u}(\varphi,f)$
$\displaystyle=\\{x\in\mathrm{Crit}(\varphi,f)\>:\>\exists
z\in\mathbb{C}^{c+1}\mbox{ s.t.\ }z^{T}J(\varphi,f)=0,z^{T}u=1\\}$ (22)
$\displaystyle\mathrm{Sing}_{u}(\varphi,f)$
$\displaystyle=\\{x\in\mathrm{Sing}(f)\>:\>\hskip 9.0pt\exists
z\in\mathbb{C}^{c+1}\mbox{ s.t.\ }z^{T}J(\varphi,f)=0,z^{T}u=1\\}$

  * (a)

There holds
$\Pi_{x}(V(\mathscr{L}_{u}))=\mathrm{Crit}_{u}(\varphi,f)\cup\mathrm{Sing}_{u}(\varphi,f)$.

  * (b)

If Assumption 1 holds then for a generic linear form
$\varphi\in\mathbb{C}[x]_{1}$, the set $\mathrm{Crit}(\varphi,V)$ is finite.

  * (c)

If $\mathrm{Crit}(\varphi,f)$ is finite then
$\mathrm{Crit}_{u}(\varphi,f)=\mathrm{Crit}(\varphi,f)$ for a generic $u$.

  * (d)

If $\mathrm{Sing}(f)$ is finite then
$\mathrm{Sing}_{u}(\varphi,f)=\mathrm{Sing}(f)$ for a generic $u$.

Consequently, if Assumption 1 holds and $\mathrm{Sing}(V)$ is finite then for
generic $\varphi,u$ set $\Pi_{x}(V(\mathscr{L}_{u}))$ is finite and satisfies
$\Pi_{x}(V(\mathscr{L}_{u}))=\mathrm{Crit}(\varphi,V)\cup\mathrm{Sing}(V)$.

###### Proof.

(a)  The following can be checked using simple linear algebra:

  * •

Suppose that $x\in\mathrm{Reg}(f)$. Then $x\in\mathrm{Crit}_{u}(\varphi,f)$ if
and only if there exists $z\in\mathbb{C}^{c+1}$ s.t. $z^{T}J(\varphi,f)=0$ and
$z^{T}u=1$, or equivalently s.t. $(x,z)$ satisfies (21).

  * •

Suppose that $x\in\mathrm{Sing}(f)$. Then $x\in\mathrm{Sing}_{u}(\varphi,f)$
if and only if there exists $z\in\mathbb{C}^{c+1}$ s.t. $z^{T}J(\varphi,f)=0$
and $z^{T}u=1$, or equivalently s.t. $(x,z)$ satisfies (21).

These facts imply that
$\Pi_{x}(V(\mathscr{L}_{u}))\cap\mathrm{Reg}(f)=\mathrm{Crit}_{u}(\varphi,f)$
and
$\Pi_{x}(V(\mathscr{L}_{u}))\cap\mathrm{Sing}(f)=\mathrm{Sing}_{u}(\varphi,f)$,
which yields the desired claim.

(b)  The claim follows from [2, Lem.7].

(c, d)  For $\mathcal{X}\in\\{\mathrm{Crit}(\varphi,f),\,\mathrm{Sing}(f)\\}$
denote by $\mathcal{X}_{u}$ the corresponding set in (22). We claim that if
$\mathcal{X}$ is finite then $\mathcal{X}_{u}=\mathcal{X}$ for a generic $u$.
Indeed, assume $\mathcal{X}=\\{x_{1},\ldots,x_{k}\\}$. By definition of
$\mathcal{X}$, for each $x_{i}\in\mathcal{X}$ there exists
$z_{i}\in\mathbb{C}^{c+1}-\\{0\\}$ with $z_{i}^{T}\cdot
J(\varphi,f)|_{x_{i}}=0$. Define
$\mathscr{U}_{i}=\\{u\in\mathbb{C}^{c+1}\>:\>z_{i}^{T}u\neq 0\\}$, and let
$\mathscr{U}=\mathscr{U}_{1}\cap\ldots\cap\mathscr{U}_{k}$. Clearly,
$\mathscr{U}$ is a non-empty Zariski-open set, and for each $u\in\mathscr{U}$
we have $\mathcal{X}_{u}=\mathcal{X}$. ∎

######  Remark 1.

The precondition of the last statement in Proposition 10 can be shown to hold
for input system $f$ in certain circumstances. For example, this is the case
for the method of [22] (described in Section 2) applied to the LMI problem
$\mathscr{A}(X)=b,X\in\mathbb{S}^{n}_{+}$ with generic input data
$(\mathscr{A},b)$: if $r_{\min}$ is the minimum rank then for at least one
subset $\iota\subseteq[n]$ of size $|\iota|=r_{\min}$ system (11) will satisfy
Assumption 1, and the variety $V$ will be smooth (i.e.
$\mathrm{Sing}(V)=\varnothing$), by [22, Thm. 2.1]. Accordingly, in this case
the CPM method will solve the problem with generic $\varphi,u$, i.e. it will
find a rational univariate representation of at least one minimum-rank
solution.

However, for non-generic instances $(\mathscr{A},b)$ there are no such
guarantees. It may happen, for example, that $\mathrm{Sing}(f)$ is positive-
dimensional, in which case the Zariski closure of
$\Pi_{x}(V(\mathscr{L}_{u}))$ is expected to be also positive-dimensional
causing the CPM method to fail.

Radical ideals  Instead of directly applying CPM to $f$, one can first compute
the radical ideal of $V(f)$ and attempt to derive a minimal set of polynomials
$g$ such that $\langle g\rangle=I(V(f))=\sqrt{\langle g\rangle}$.
Subsequently, CPM is applied to $g$. This approach is referred to as CPMrad.

In our experiments, we utilized Macaulay2 software [15] to compute the radical
and then employed the command mingens that attempts to find a smaller set of
its generators. Note that it is not guaranteed to produce a minimal generating
set (unless the radical ideal is homogeneous). We observed that CPMrad was
capable of solving some systems that CPM could not resolve. In the latter
case, the ideal $\langle f\rangle$ was not radical, while the former was able
to construct a polynomial system whose ideal was radical and complete
intersection.

##  5 Numerical Results

This section compares our method, denoted as “Hybrid”, with the method in [22]
that we call “HNS”.222The method of [22] is implemented in software [23], but
some details are different, e.g. [23] relies on FGb (deprecated) instead of
the more recent softsware msolve for solving 0-dimensional polynomial systems.
For a more fair comparison we used msolve both for the new method and for the
algorithm in [22]. Note that the method in [22] has been improved and adapted
to the case when the ${{P}}$ is non-generic in [25], but this variant is not
implemented and it can be considered computationally more demanding, inasmuch
as it would involve computation of bivariate rational representations of
algebraic curves, instead of univariate representations.  Their details are
described below.

Hybrid For the implementation Algorithm 1 requires a numerical solution
$\tilde{X}$ in the neighborhood
$\tilde{\Omega}=\\{\tilde{X}\in\mathbb{S}^{n}\>:\>||\tilde{X}-X^{\ast}||_{F}<\delta\\}$
of an exact solution $X^{*}$ as in theorem 6. Since we are interested in
weakly feasible SDPs, we cannot use standard interior-point methods that are
designed to work for instances in which both primal and dual problems are
strongly feasible. Fortunately, weakly feasible SDPs can be tackled via facial
reduction approaches [6, 7, 37, 45, 38, 39, 46, 29]. We chose to use the
method in [18] which combines facial reduction with a numerical algebraic
geometry algorithm based on the Bertini software [4]. The resulting polynomial
system (15) was solved with a critical point method (CPM or CPMrad) as
described in Section 4.2.

HNS This method constructs system (11) for every subset $\iota\subseteq[n]$ of
increasing cardinality $r=|\iota|$, ranging from $0$ to $n-1$, and applies the
critical point method to each such variety. It stops if it finds a PSD matrix
(with entries in a rational univariate representation). By Theorem 2, we can
expect to obtain a minimal-rank solution in $P$ for at least one $\iota$ of
size $|\iota|=r_{\min}$. Note that in [22] only CPM was used, which was shown
to be sufficient for generic instances $(\mathscr{A},b)$. Since we deal with
non-generic LMI instances, we tested both CPM and CPMrad.

The Lagrange system (21) was solved using library msolve [5]. We say that this
computation succeeds if the projection to variables $X$ is zero-dimensional,
and at least one of the solutions (represented in a rational univariate
representation) is a PSD matrix. Otherwise the computation fails.

Instances  We applied the methods to various instances of weakly feasible SDPs
extracted from the existing literature; their description can be found in
https://git.ista.ac.at/jzapata/hybrid-method. Note that these instances are
quite sparse; to test how the methods cope with non-sparse instances, for each
SDP linear map $\mathscr{A}(X)=(\left\langle
A_{1},X\right\rangle,\ldots,\left\langle A_{n},X\right\rangle)$ we created two
input instances denoted as clean (I) and rotated (T). The latter uses the map

$\mathscr{A}_{T}(X)=\left(\left\langle
T^{\top}A_{1}T,X\right\rangle,\ldots,\left\langle
T^{\top}A_{n}T,X\right\rangle\right)=\mathscr{A}(TXT^{\top})$

for some random matrix $T\in GL(n,\mathbb{R})$. This transformation does not
affect feasibility of the system, but it allows us to generate more weakly
feasible SDPs that are less sparse and may become more challenging for
numerical methods.

Results  The outcomes of our experiments are shown in Table 1. Values
$r_{\min}$ / $r_{\max}$ given in the second column are the ranks of solutions
found by HNS / Hybrid (recall that these methods search for minimum- and
maximum-rank solutions, respectively). The results are given in the format
“$\,t\,/\,t_{\tt rad}$” where $t$ is the runtime of CPM and $t_{\tt rad}$ is
the runtime of CPMrad (in seconds). For Hybrid we report additionally the time
for computing approximate solution; it is given in square brackets. Note that
this time is included in both $t$ and $t_{\tt rad}$. If the method fails then
the corresponding runtime is crossed; for example, entry “
$\xcancel{5.8}\,/\,3.1$” in the first row for HNS means that CPM failed for
all considered polynomial systems (corresponding to different subsets
$\iota\subseteq[n])$ while CPMrad succeeded for at least one system.

In cases where the method reaches the predefined time limit of 10 minutes, the
computation is terminated, and the respective output is labeled with $\infty$.

Discussion of results  Hybrid was able to solve most of our instances both
with CPM and CPMrad. In contrast, HNS failed more often, especially when CPM
was used. Interestingly, in several cases HNS with CPMrad solved the clean
version of the problem but failed on the rotated version. We conjecture that
in these cases Macaulay2 did not find minimal generators of the radical ideal.

We can conclude that in many cases polynomial systems constructed by Hybrid
appear to be easier to solve compared to systems in HNS, and also Hybrid is
usually faster than HNS. (Note that the latter needs to solve many more
systems). However, there are several exceptions: HNS with CPMrad could solve
clean versions of the last four rows in Table 1, whereas Hybrid failed or did
not terminate on these instances.

In Appendix A we expand the first example from Table 1 (the clean version of
DruWo2017). We show, in particular, that in the failed cases varieties
$\mathrm{Sing}(f)$ for given polynomial systems $f$ are positive-dimensional.

We also confirm that system (15) after fixing variables in Hybrid for this
example indeed has a zero-dimensional real component, as predicted by Theorem
6. We observe that the system also has two other real components, which are
positive-dimensional. After applying CPM, the system becomes 0-dimensional,
and so Hybrid with CPM succeeds.

SDP | $n/r_{\min}/r_{\max}$ | clean instances (I) | rotated instances (T)
---|---|---|---
HNS | Hybrid | HNS | Hybrid
DruWo2017-2.3.2P | 3/1/2 |  ~~5.8~~ / 3.1 | 1.1/1.6 [0.5] |  ~~4.6~~ / ~~12.7~~ | 1.2 / 1.6 [0.6]
Gupta2013-12.3P | 3/1/2 |  ~~5.9~~ / 3.1 | 1.7/1.9 [1.0] |  ~~4.8~~ / ~~12.4~~ | 1.7 / 2.0 [1.1]
Hauenstein2.6P | 3/1/2 |  ~~5.8~~ / 3.1 | 1.5/1.9 [0.8] |  ~~4.4~~ / ~~12.8~~ | 2.0 / 1.9 [1.3]
Helmberg2000-2.2.1P | 3/1/2 | 2.6 / 2.8 | 1.7/1.9 [1.2] | 2.1 / 3.4 | 1.8 / 1.8 [1.3]
LauVall2020-2.5.1P | 2/1/1 | 1.6 / 1.8 | 1.1/1.3 [0.7] | 1.2 / 2.0 | 1.1 / 1.5 [0.7]
LauVall2020-2.5.2P | 3/1/2 |  ~~6.2~~ / 3.1 | 1.6/1.9 [1.0] |  ~~4.4~~ / ~~12.5~~ | 1.7 / 2.4 [1.0]
Pataki2017-4P | 3/1/2 |  ~~5.9~~ / 3.0 | 1.5/1.9 [0.9] |  ~~4.5~~ / ~~12.6~~ | 1.6 / 2.0 [0.9]
deKlerk2002-2.1P | 2/1/1 | 1.4 / 1.8 | 1.1/1.6 [0.7] | 1.5 / 1.9 | 1.1 / 1.3 [0.7]
DruWo2017-2.3.2D | 3/1/2 | 2.4 / 2.8 | 1.7/2.0 [1.2] | 2.6 / 6.2 | 2.0 / 2.1 [1.3]
Gupta2013-12.3D | 3/1/2 | 2.4 / 2.8 | 1.6/2.0 [1.1] | 2.2 / 3.2 | 1.8 / 2.0 [1.3]
HNS2020-4.1D | 4/2/2 | 8.7 / 9.7 | 2.4/2.9 [1.7] | 7.0 / 17.6 | 2.9 / 3.1 [2.0]
Hauenstein2.6D | 3/1/2 | 2.5 / 2.8 | 1.8/1.5 [1.2] | 2.0 / 6.2 | 1.9 / 2.2 [1.3]
Helmberg2000-2.2.1D | 3/1/2 |  ~~5.9~~ / 3.0 | 1.6/1.8 [1.0] |  ~~4.1~~ / 12.3 | 1.8 / 2.1 [1.1]
Pataki2017-4D | 3/1/2 | 2.4 / 2.9 | 1.7/1.9 [1.2] | 2.0 / 3.4 | 2.5 / 1.8 [1.6]
Permenter2018-4.3.1D | 5/0/1 | 11.5 / 1.0 | 5.5/8.1 [4.5] | 1.0 / 1.5 | 6.6 / 6.7 [5.6]
Permenter2018-4.3.2D | 4/2/2 | 1.0 / 1.5 | 2.3/2.7 [1.6] | 9.1 / 129.6 | 2.6 / 3.1 [1.5]
PatakiCleanDim3P | 3/1/2 |  ~~4.7~~ / 3.6 | 1.6/1.5 [1.0] |  ~~4.6~~ / ~~12.5~~ | 1.7 / 1.9 [1.0]
PatakiCleanDim4P | 4/1/- |  ~~49.5~~ / 4.4 |  ~~2.3~~ /~~3.3~~ [1.3] |  ~~77.3~~ / ~~133.4~~ |  ~~27.2~~ / $\infty$ [2.6]
PatakiCleanDim5P | 5/1/- |  $\infty$ / 6.9 |  ~~154.6~~ /~~519.2~~ [3.0] |  $\infty$ / $\infty$ |  $\infty$ / $\infty$ [9.5]
PatakiCleanDim6P | 6/1/- |  $\infty$ / 12.4 |  $\infty$/$\infty$ [11.5] |  $\infty$ / $\infty$ |  $\infty$ / $\infty$ [68.4]
HeNaSa2016-6.2P | 6/2/- | 33.1 / 36.3 |  $\infty$/$\infty$ [9.3] |  $\infty$ / $\infty$ |  $\infty$ / $\infty$ [10.5]
Table 1: Comparison between Hybrid and HNS for clean and rotated instances.

##  6 Conclusions and future work

We presented a hybrid method for certifying weakly feasible SDP problems that
uses an approximate numerical SDP solver and an exact solver for polynomial
systems over reals. Our numerical results indicate that the hybrid method can
outperform a pure exact algorithm from [22] when given a good approximate
solution.

In our current experiments scalability was limited both by the numerical SDP
solver that we employed (which was the facial reduction algorithm [18] based
on Bertini [4]), and by the exact solver for polynomial systems (msolve). We
believe that our approach can handle larger instances if the method in [18] is
replaced with an alternative facial reduction algorithm, and an exact solver
for polynomial systems is replaced with a numerical solver; this is left as a
future work. Potential candidates for the latter could be algorithms from
numerical algebraic geometry such as Bertini. Using such solver could address
the issue discussed in the previous section: if a system of polynomial
equations has a zero-dimensional component and a positive-dimensional
component, a numerical solver would focus on the desired zero-dimensional
component if the previous step found a good approximation.

Acknowledgments  This work has been partially supported by the ANR JCJC
project ANR-21-CE48-0006-01 “HYPERSPACE”.

##  Appendix A Example: DruWo2017

In this section we illustrate how Hybrid and HNS work on (the clean version
of) DruWo2017. (Recall that on this example HNS with CPM fails, and all other
methods succeed.) Our conclusions are as follows:

  * •

Hybrid. The resulting polynomial system $F$ has
$\mathrm{Sing}(F)=\\{(X^{\ast},Y^{\ast})\\}$ where $(X^{\ast},Y^{\ast})$ is
the desired solution, and CPM succeeds. We write down the rational univariate
representation of the solution obtained with msolve.

  * •

HNS. In this case there are two polynomial systems to consider, $F_{1}$ and
$F_{2}$. Sets $\mathrm{Sing}(F_{1})$ and $\mathrm{Sing}(F_{2})$ are both
positive-dimensional, which causes CPM to fail for $F_{1}$ and $F_{2}$.
Furthermore, the ideals $\langle F_{1}\rangle$ and $\langle F_{2}\rangle$ are
not radical. We then compute generators of $\sqrt{\langle F_{1}\rangle}$ and
$\sqrt{\langle F_{2}\rangle}$ with Macaulay2. These generators turn out to be
non-minimal in the first case, and minimal in the second case. CPMrad fails on
$F_{1}$ but succeeds on $F_{2}$.

The example comes from [11, Example 2.3.2], and can be written as follows
(over matrices $X\in\mathbb{S}^{4}$):

$\displaystyle X_{33}$ $\displaystyle=0$ (23a) $\displaystyle 2X_{13}+X_{22}$
$\displaystyle=1$ (23b) $\displaystyle X$ $\displaystyle\succeq 0$ (23c)

Its solution set $P$ is

$P=\left\\{\begin{pmatrix}x&y&0\\\ y&1&0\\\
0&0&0\end{pmatrix}\>:\>\begin{pmatrix}x&y\\\ y&1\end{pmatrix}\succeq
0\right\\}\hskip 30.0pt$

###  A.1 Hybrid method

The minimal face ${\mathcal{F}}$ containing $P$ is given by

${\mathcal{F}}=\left\\{\begin{pmatrix}x&y&0\\\ y&z&0\\\
0&0&0\end{pmatrix}\>:\>\begin{pmatrix}x&y\\\ y&z\end{pmatrix}\succeq
0\right\\}$

This face is described by matrix

$U=\begin{pmatrix}1&0\\\ 0&1\\\ 0&0\end{pmatrix}$

as in (8),(9) (with $\iota=\\{1,2\\}$), which corresponds to matrix
$Y_{U}=\left(\begin{smallmatrix}0\\\ 0\end{smallmatrix}\right)$. Suppose that
in step 1 we found approximate solution $\tilde{X}$ with entries
$\tilde{X}_{11}=1$ and $\tilde{X}_{12}=0$ (and other entries are close to the
unique true values). Also suppose that step 2 correctly finds
$\iota=\\{1,2\\}$, which would give
$\tilde{Y}\approx\left(\begin{smallmatrix}0\\\ 0\end{smallmatrix}\right)$. If
$\tilde{X}_{ij}\neq X^{\ast}_{ij}$ for all $i,j$ then the only way to get a
system that contains a 0-dimensional real component with a correct solution is
to fix variables $X_{11}=\tilde{X}_{11}$, $X_{12}=\tilde{X}_{12}$. Suppose
that step 4 correctly identifies these variables, then in step 5 we obtain the
system

$\begin{array}[]{rcl}X_{33}&=&0\\\ 2X_{13}+X_{22}&=&1\\\ X_{11}&=&1\\\
X_{12}&=&0\end{array}\qquad\begin{pmatrix}X_{11}&X_{12}&X_{13}\\\
X_{12}&X_{22}&X_{23}\\\
X_{13}&X_{23}&X_{33}\end{pmatrix}\begin{pmatrix}Y_{1}\\\ Y_{2}\\\
1\end{pmatrix}=0^{3\times 2}$ (24)

Making substitutions $X_{11}=1$, $X_{12}=0$, $X_{33}=0$,
$X_{13}=\frac{1}{2}-\frac{1}{2}X_{22}$ and eliminating these variables yields
the system

$\displaystyle
F:\left\\{\begin{array}[]{rcl}Y_{1}-\frac{X_{22}}{2}+\frac{1}{2}&=0\\\
X_{23}+X_{22}Y_{2}&=0\\\
X_{23}Y_{2}-Y_{1}\left(\frac{X_{22}}{2}-\frac{1}{2}\right)&=0\end{array}\right.$
(28)

Its set of solutions is

$V(F)=\left\\{(X_{22},X_{33},Y_{1},Y_{2})=\left(2t+1,-2\sqrt{-\tfrac{1}{2t+1}}t^{2}\pm\sqrt{-\tfrac{1}{2t+1}}t,t,\sqrt{-\tfrac{1}{2t+1}}t\right)\>:\>t\in\mathbb{C}\setminus\left\\{-\tfrac{1}{2}\right\\}\right\\}$

The solutions are real for $t\in\\{0\\}\cup(-\infty,-\tfrac{1}{2})$.
Accordingly, the real variety $\mathbb{R}V(F)$ contains a zero-dimensional
component with the correct solution $(1,0,0,0)$ (as predicted by Theorem 6)
but also two one-dimensional components corresponding to
$t\in(-\infty,-\tfrac{1}{2})$ with “plus” and “minus” signs.

To compute $\mathrm{Sing}(F)$, we need to evaluate the rank of the Jacobian
matrix

$J(F)=\begin{pmatrix}-\frac{1}{2}&0&1&0\\\ Y_{2}&1&0&X_{22}\\\
-\frac{Y_{1}}{2}&Y_{2}&\frac{1}{2}-\frac{X_{22}}{2}&X_{23}\end{pmatrix}$

at points in $V(F)$. This rank equals 2 at the point corresponding to $t=0$,
and is 3 at all other points. Hence, we have

$\mathrm{Sing}(F)=\\{\left\\{(X_{22},X_{23},Y_{1},Y_{2})=(1,0,0,0)\right\\}$

Next, we formed the Lagrangian system $\mathscr{L}_{u}$ as in (21) for a
randomly chosen linear map $\varphi$ and vector $u\in\mathbb{Z}^{4}$, and
solved it using msolve. It returned the following rational parametrization
encoding its zero-dimensional solution set:

$\Pi_{{X_{22},X_{23}}}(V(\mathscr{L}_{u}))=\left\\{\left(\frac{q_{1}(t)}{q^{\prime}(t)},\frac{q_{2}(t)}{q^{\prime}(t)}\right):q(t)=0\right\\}$

where

$\displaystyle q_{1}(t)$ $\displaystyle=-129t^{4}-240t^{3}-76t^{2}-16t-4,$
$\displaystyle q_{2}(t)$
$\displaystyle=t\left(-135t^{4}+880t^{3}+60t^{2}+16t+4\right),$ $\displaystyle
q(t)$ $\displaystyle=27t^{5}-220t^{4}-20t^{3}-8t^{2}-4t.$

The floating-point representation (up to 5 decimal places) is

$\Pi_{{X_{22},X_{23}}}(V(L_{u}))\approx\left\\{\begin{array}[]{l}(-10.18376,17.84481),\\\
(-0.08682,-0.16012),\\\ (1,0),\\\ (0.04640-0.13399i,0.08358+0.16089i),\\\
(0.04640+0.13399i,0.08358-0.16089i)\end{array}\right\\}$

###  A.2 HNS method

Since HNS searches for a minimum rank solution, we need to use $|\iota|=1$.
Consequently, three possibilities arise for the kernel matrix:

$\begin{pmatrix}Y_{11}&Y_{12}\\\ 1&0\\\ 0&1\\\
\end{pmatrix},\begin{pmatrix}1&0\\\ Y_{21}&Y_{22}\\\ 0&1\\\
\end{pmatrix},\begin{pmatrix}1&0\\\ 0&1\\\ Y_{31}&Y_{32}\\\ \end{pmatrix}$

It can be checked that for the third matrix $Y$ system $XY=0$ does not have
solutions with $X\in P$. Below we analyze the two remaining systems.

  1. (1)

Making substitutions $X_{33}=0$, $X_{13}=\frac{1}{2}-\frac{1}{2}X_{22}$ and
eliminating these variables yields

$\begin{pmatrix}X_{11}&X_{12}&\frac{1}{2}-\frac{X_{22}}{2}\\\
X_{12}&X_{22}&X_{23}\\\ \frac{1}{2}-\frac{X_{22}}{2}&X_{23}&0\\\
\end{pmatrix}\begin{pmatrix}Y_{11}&Y_{12}\\\ 1&0\\\ 0&1\\\
\end{pmatrix}=0^{3\times 2}\quad\Longleftrightarrow\quad
F_{1}:\left\\{\begin{array}[]{rcl}X_{12}+X_{11}Y_{11}&=&0\\\
X_{22}+X_{12}Y_{11}&=&0\\\ X_{11}Y_{12}-\frac{X_{22}}{2}+\frac{1}{2}&=&0\\\
X_{23}+X_{12}Y_{12}&=&0\\\
-Y_{12}\left(\frac{X_{22}}{2}-\frac{1}{2}\right)&=&0\end{array}\right.$

Recall that the equation at position (3,1) on the LHS is implied by other
equations (see step 5 of Algorithm 1), and accordingly we removed it on the
RHS. This polynomial system has the following parametric description for its
set of solutions:

$V(F_{1})=\left\\{(X_{11},X_{12},X_{22},X_{23},Y_{11},Y_{12})=\left(\tfrac{1}{t^{2}},-\tfrac{1}{t},1,0,t,0\right)\>:\>t\in\mathbb{C}\setminus\left\\{0\right\\}\right\\}$

If we compute the Jacobian and replace the parametric description we can see
the matrix has rank 4:

$J(F_{1})=\begin{pmatrix}Y_{11}&1&0&0&X_{11}&0\\\ 0&Y_{11}&1&0&X_{12}&0\\\
Y_{12}&0&-\tfrac{1}{2}&0&0&X_{11}\\\ 0&Y_{12}&0&1&0&X_{12}\\\
0&0&-\frac{Y_{12}}{2}&0&0&\tfrac{1}{2}-\tfrac{X_{22}}{2}\end{pmatrix}=\begin{pmatrix}t&1&0&0&\frac{1}{t^{2}}&0\\\
0&t&1&0&-\frac{1}{t}&0\\\ 0&0&-\frac{1}{2}&0&0&\frac{1}{t^{2}}\\\
0&0&0&1&0&-\tfrac{1}{t}\\\ 0&0&0&0&0&0\end{pmatrix}$

Moreover, we can see the ideal $I=\langle F_{1}\rangle$ is not radical as
$(X_{22}-1)^{2}\in I$ but $X_{22}-1\not\in I$. Using Macaulay2 we can compute
its radical ideal and obtain

$\sqrt{I}=\langle
X_{23},Y_{12},X_{22}-1,X_{12}^{2}-X_{11},X_{12}Y_{11}+1,X_{12}+X_{11}Y_{11}\rangle$

but in this occasion the set of generators is not minimal since

$X_{12}+X_{11}Y_{11}=X_{12}(X_{12}Y_{11}+1)-Y_{11}(X_{12}^{2}-X_{11})$

  1. (2)

Similarly, for the second system we have

$\begin{pmatrix}X_{11}&X_{12}&\frac{1}{2}-\frac{X_{22}}{2}\\\
X_{12}&X_{22}&X_{23}\\\ \frac{1}{2}-\frac{X_{22}}{2}&X_{23}&0\\\
\end{pmatrix}\begin{pmatrix}1&0\\\ Y_{21}&Y_{22}\\\ 0&1\\\
\end{pmatrix}=0^{3\times 2}\quad\Longleftrightarrow\quad
F_{2}:\left\\{\begin{array}[]{rcl}X_{11}+X_{12}Y_{21}&=&0\\\
X_{12}+X_{22}Y_{21}&=&0\\\ X_{23}Y_{21}-\frac{X_{22}}{2}+\frac{1}{2}&=&0\\\
X_{23}+X_{22}Y_{22}&=&0\\\ X_{23}Y_{22}&=&0\end{array}\right.$

Again, the ideal $I=\langle F_{2}\rangle$ is not radical but when we obtain
its radical in Macaulay2 we get a minimal set of generators:

$\sqrt{I}=\langle
X_{23},Y_{22},X_{12}+Y_{21},X_{22}-1,X_{11}-Y_{21}^{2}\rangle$

To summarize, CPM fails for $F_{1},F_{2}$ while CPMrad fails for $F_{1}$ but
succeeds for $F_{2}$.

######  Remark 2.

As Table 1 shows, HNS with both CPM and CPMrad fails for the rotated version
of DruWo2017. The resulting polynomial systems were too big to analyze by
hand; we conjecture that in these cases Macaulay2 also failed to find minimal
generators of radical ideals, as for system $F_{1}$.

## References

  * [1] K. Atkinson.  An introduction to numerical analysis.  Wiley, 2nd edition, 1991.
  * [2] B. Bank, M. Giusti, J. Heintz, and L.-M. Pardo.  Generalized polar varieties: geometry and algorithms.  Journal of Complexity, 21(4):377–412, 2005.
  * [3] S. Basu, R. Pollack, and M.-F. Roy.  Algorithms in real algebraic geometry, volume 10 of Algorithms and Computation in Mathematics.  Springer-Verlag, second edition, 2006.
  * [4] D. Bates, J. Hauenstein, A. Sommese, and C. Wampler.  Numerically Solving Polynomial Systems with Bertini.  SIAM, Philadelphia, PA, USA, 2013.
  * [5] J. Berthomieu, C. Eder, and M. Safey El Din.  msolve: A library for solving polynomial systems.  In Proceedings of the 2021 International Symposium on Symbolic and Algebraic Computation, pages 51–58, 2021.
  * [6] J. Borwein and H. Wolkowicz.  Facial reduction for a cone-convex programming problem.  J. Aust. Math. Soc. (Ser. A), 30:369–380, 1981.
  * [7] J. Borwein and H. Wolkowicz.  Regularizing the abstract convex program.  J. Math. Anal. Appl, 83:495–530, 1981.
  * [8] D. Cox, J. Little, and D. O’Shea.  Ideals, varieties, and algorithms: an introduction to computational algebraic geometry and commutative algebra.  Springer, 2007.
  * [9] M. Dostert, D. de Laat, and P. Moustrou.  Exact semidefinite programming bounds for packing problems.  SIAM J. Optim., 31(2):1433–1458, 2021.
  * [10] Y. Drori and M. Teboulle.  Performance of first-order methods for smooth convex minimization: A novel approach.  Mathematical Programming, 145(1):451–482, 2014.
  * [11] D. Drusvyatskiy and H. Wolkowicz.  The many faces of degeneracy in conic optimization.  Foundations and Trends in Optimization, 3(2):77–170, 2017.
  * [12] D. Eisenbud.  Linear sections of determinantal varieties.  American J. Math., 110(3):541–575, 1988.
  * [13] M. Goemans and D. Williamson.  Improved approximation algorithms for maximum cuts and satisfiability problems using semidefinite programming.  Journal of the ACM, 42:1115–1145, 1995.
  * [14] G. Golub and C. Van Loan.  Matrix Computations.  Johns Hopkins Studies in the Mathematical Sciences. Johns Hopkins University Press, 2013.
  * [15] D. Grayson and M. Stillman.  Macaulay2, a software system for research in algebraic geometry.  Available at http://www.math.uiuc.edu/Macaulay2/.
  * [16] J. Harris and L. W. Tu.  On symmetric and skew-symmetric determinantal varieties.  Topology, 23(1):71–84, 1984.
  * [17] J. Harrison.  Verifying nonlinear real formulas via sums of squares.  In Theorem Proving in Higher Order Logics, pages 102–118. Springer, 2007.
  * [18] J. D. Hauenstein, A. C. L. Jr., S. McPherson, and Y. Zhang.  Numerical algebraic geometry and semidefinite programming.  Results in Applied Mathematics, 11, 2021.
  * [19] J. D. Hauenstein, A. C. Liddell Jr, S. McPherson, and Y. Zhang.  Numerical algebraic geometry and semidefinite programming.  Results in Applied Mathematics, 11:100166, 2021.
  * [20] D. Henrion, M. Korda, and J. B. Lasserre.  The Moment-SOS Hierarchy. Lectures In Probability, Statistics, Computational Geometry, Control And Nonlinear Pdes, volume 4.  World Scientific, 2020.
  * [21] D. Henrion, S. Naldi, and M. Safey El Din.  Real root finding for determinants of linear matrices.  J. of Symbolic Computation, 74:205–238, 2015.
  * [22] D. Henrion, S. Naldi, and M. Safey El Din.  Exact algorithms for linear matrix inequalities.  SIAM J. Optim., 26(4):2512–2539, 2016.
  * [23] D. Henrion, S. Naldi, and M. Safey El Din.  SPECTRA: a Maple library for solving linear matrix inequalities in exact arithmetic.  Optimization Methods and Software, 34(1):62–78, 2019.
  * [24] D. Henrion, S. Naldi, and M. Safey El Din.  Real root finding for low rank linear matrices.  Appl. Algebra Eng. Comm. Comput., 31(2):101–133, 2020.
  * [25] D. Henrion, S. Naldi, and M. Safey El Din.  Exact algorithms for semidefinite programs with degenerate feasible set.  J. Symb. Comput., 104:942–959, 2021.
  * [26] E. L. Kaltofen, B. Li, Z. Yang, and L. Zhi.  Exact certification in global polynomial optimization via sums-of-squares of rational functions with rational coefficients.  J. Symbolic Comput., 47(1):1–15, 2012.
  * [27] I. Klep and M. Schweighofer.  An exact duality theory for semidefinite programming based on sums of squares.  Mathematics of Operations Research, 38(3):569–590, 2013.
  * [28] J.-B. Lasserre.  Global optimization with polynomials and the problem of moments.  SIAM J. Optim., 11(3):796–817, 2001.
  * [29] B. Lourenço, M. Muramatsu, and T. Tsuchiya.  Solving SDP completely with an interior point oracle.  Optimization Methods and Software, pages 1–47, 2021.
  * [30] D. Monniaux and P. Corbineau.  On the generation of Positivstellensatz witnesses in degenerate cases.  In International Conference on Interactive Theorem Proving, pages 249–264. Springer, 2011.
  * [31] S. Naldi.  Exact algorithms for determinantal varieties and semidefinite programming.  Theses, Université de Toulouse, Sept. 2015.
  * [32] S. Naldi.  Solving rank-constrained semidefinite programs in exact arithmetic.  J. Symb. Comput., 85(C):206–223, 2018.
  * [33] S. Naldi and R. Sinn.  Conic programming: infeasibility certificates and projective geometry.  J. Pure Appl. Algebra, page 106605, 2020.
  * [34] J. Nie, K. Ranestad, and B. Sturmfels.  The algebraic degree of semidefinite programming.  Mathematical Programming, 122(2):379–405, 2010.
  * [35] C.-T. Pan.  On the existence and computation of rank-revealing LU factorizations.  Linear Algebra Appl., 316(1-3):199–222, 2000.
  * [36] A. Papachristodoulou and S. Prajna.  A tutorial on sum of squares techniques for systems analysis.  In American Control Conference, pages 2686–2700, 2005.
  * [37] G. Pataki.  Strong duality in conic linear programming: facial reduction and extended duals.  In Computational and Analytical Mathematics, volume 50, pages 613–634. Springer, New York, 2013.
  * [38] F. Permenter, H. Friberg, and E. Andersen.  Solving conic optimization problems via self-dual embedding and facial reduction: a unified approach.  SIAM. J. Optim., 27:1257–1282, 2017.
  * [39] F. Permenter and P. Parrilo.  Partial facial reduction: simplified, equivalent SDPs via approximations of the PSD cone.  Mathematical Programming (Ser. A), 171:1–54, 2018.
  * [40] H. Peyrl and P. Parrilo.  Computing sum of squares decompositions with rational coefficients.  Theoretical Computer Science, 409(2):269 – 281, 2008.  Symbolic-Numerical Computations.
  * [41] A. Platzer, J.-D. Quesel, and P. Rümmer.  Real world verification.  In CADE, pages 485–501, 2009.
  * [42] P. Roux, Y.-L. Voronin, and S. Sankaranarayanan.  Validating numerical semidefinite programming solvers for polynomial invariants.  Formal Methods in System Design, 53(2):286–312, 2018.
  * [43] M. Safey El Din and E. Schost.  Polar varieties and computation of one point in each connected component of a smooth real algebraic set.  In ISSAC’03, pages 224–231. ACM, 2003.
  * [44] L. Schork and J. Gondzio.  Rank revealing gaussian elimination by the maximum volume concept.  Linear Algebra and its Applications, 592:1–19, 2020.
  * [45] H. Waki and M. Muramatsu.  Facial reduction algorithms for conic optimization problems.  J. Optim. Theory. Appl., 158:188–215, 2013.
  * [46] Y. Zhu, G. Pataki, and Q. Tran-Dinh.  Sieve-sdp: a simple facial reduction algorithm to preprocess semidefinite programs.  Math. Prog. Comput., 11(3):503–586, 2019.