# Determinantal sets, singularities and application to optimal control in medical imagery

Bernard Bonnard, Jean-Charles Faugère, Alain Jacquemard, Mohab Safey El Din, Thibaut Verron

► To cite this version:

Bernard Bonnard, Jean-Charles Faugère, Alain Jacquemard, Mohab Safey El Din, Thibaut Verron. Determinantal sets, singularities and application to optimal control in medical imagery. International symposium on symbolic and algebraic computations, Jul 2016, Waterloo, Canada. pp.103-110, 10.1145/2930889.2930916. hal-01307073v3

HAL Id: hal-01307073

https://inria.hal.science/hal-01307073v3

Submitted on 5 Jul 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Determinantal Sets, Singularities and Application to Optimal Control in Medical Imagery

Bernard Bonnard<sup>be</sup>
Bernard.Bonnard@u-
bourgogne.fr

Jean-Charles Faugère<sup>dac</sup>
Jean-Charles.
Faugere@inria.fr

Alain Jacquemard<sup>dbac</sup>
Alain.Jacquemard@u-
bourgogne.fr

Mohab Safey El Din<sup>acd</sup>
Mohab.Safey@lip6.fr

Thibaut Verron<sup>acd</sup>
Thibaut.Verron@lip6.fr

<sup>a</sup>Sorbonne Universités, UPMC Univ Paris 06, UMR 7606, LIP6, F-75005, Paris, France

<sup>b</sup>Institut de Mathématiques de Bourgogne, UMR CNRS 5584, Dijon F-21078, France

<sup>c</sup>CNRS, UMR 7606, LIP6, F-75005, Paris, France

<sup>d</sup>Inria, Paris Center, PolSys Project

<sup>e</sup>Inria, Sophia Antipolis, McTao Project

## ABSTRACT

Control theory has recently been involved in the field of nuclear magnetic resonance imagery. The goal is to control the magnetic field optimally in order to improve the contrast between two biological matters on the pictures.

Geometric optimal control leads us here to analyze meromorphic vector fields depending upon physical parameters, and having their singularities defined by a determinantal variety. The involved matrix has polynomial entries with respect to both the state variables and the parameters. Taking into account the physical constraints of the problem, one needs to classify, with respect to the parameters, the number of real singularities lying in some prescribed semi-algebraic set.

We develop a dedicated algorithm for real root classification of the singularities of the rank defects of a polynomial matrix, cut with a given semi-algebraic set. The algorithm works under some genericity assumptions which are easy to check. These assumptions are not so restrictive and are satisfied in the aforementioned application. As more general strategies for real root classification do, our algorithm needs to compute the critical loci of some maps, intersections with the boundary of the semi-algebraic domain, etc. In order to compute these objects, the determinantal structure is exploited through a stratification by the rank of the polynomial matrix. This speeds up the computations by a factor 100. Furthermore, our implementation is able to solve the application in medical imagery, which was out of reach of more general algorithms for real root classification. For instance, computational results show that the contrast problem where one of the matters is water is partitioned into three distinct classes.

## Keywords

Polynomial system solving; Real geometry; Applications

Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from permissions@acm.org.

ISSAC '16, July 19-22, 2016, Waterloo, ON, Canada

© 2016 ACM. ISBN 978-1-4503-4380-0/16/07...$15.00

DOI: http://dx.doi.org/10.1145/2930889.2930916

## 1. INTRODUCTION

**Motivations and problem description.** Nuclear Magnetic Resonance (NMR) is a powerful tool in medical imagery. In order to distinguish two biological matters on a picture, it is required to optimize the contrast between the two matters. Because of its importance in medical sciences, this contrast imaging problem has received a lot of attention. The pioneering work of [3] has established geometric optimal control techniques as a major tool for designing optimal control strategies for the problem of improving the contrast.

These strategies depend on the biological matters under study. In NMR imagery the main physical parameters involved are the longitudinal and transversal relaxation times of each matter. This approach is formalized in [2]. It requires us to solve the following real root classification problem.

Consider a $k \times k$ matrix $M$ whose coefficients are polynomials in $Q[X_1, ..., X_n, G_1, ..., G_t]$, and assume that $n = (k-r+1)^2$ with $r \in \{1, ..., k-1\}$. Let $\pi: \mathbb{C}^n \times \mathbb{C}^t \rightarrow \mathbb{C}^t$ be the canonical projection. Let $V_r \subset \mathbb{C}^n \times \mathbb{C}^t$ be the set of points at which $M$ has rank $r$, and let $V$ be the union of the singular locus of $V_r$ and of the critical points of $\pi$ restricted to $V_r$. Generically, this variety has dimension $n+t - (k-r+1)^2 = t$. Also consider a semi-algebraic set $B$ in $\mathbb{R}^n \times \mathbb{R}^t$. Assume that $B$ has non-empty interior and that there exists a Zariski-open set $O \subset \mathbb{C}^t$ such that $V \cap \pi^{-1}(g)$ is a non-empty finite set for $g \in O$. Further assume that the projection $\pi$ restricted to $V \cap B$ is proper ([10, Def. 2.10.1]). We aim at describing connected open sets $C_1, ..., C_\ell \subset \mathbb{R}^t$ such that $\bigcup_{i=1}^\ell C_i$ is dense in $\mathbb{R}^t$ (for the Euclidean topology) and the cardinality of $V \cap \pi^{-1}(g)$ inside $B \cap \pi^{-1}(g)$ is invariant when $g$ ranges over $C_i$.

For our application, the size of the matrix is 4, and the target rank is $k-1=3$. The number of variables and parameters are $n=4$ and $t=4$ for the general case. Since the system is homogeneous in the parameters, we may set one of the parameters to 1, reducing the problem to $t=3$. An important particular case is when one of the matters is water; then the corresponding relaxation times are 1, and the number of parameters is $t=2$.

Such a determinantal structure is general enough to design dedicated algorithms. We also mention that the optimal control problems in [6] lead to algebraic classification problems with similar structures.
**State-of-the-art.** The modeling through an optimal control problem is introduced in [2]. The so-called Bloch modeling and saturation method for tackling this problem is developed therein.

In [2], four experimental important cases are studied (all parameters of the classification problem are fixed). Among other properties, it has been observed there that the number of singularities is constant when water is involved. This led to the following questions:

1. Is this number of singularities preserved for any choice of a second matter, the first one being water?

2. If not, how many different classes of pairs of matters can we distinguish through the analysis of those singularities?

Answering these questions leads to the real root classification problem described above. Symbolic computation techniques are good candidates to solve them.

### Properties of Cylindrical Algebraic Decomposition
(CAD) adapted to a given polynomial family allow to solve real root classification problems. Hence the CAD algorithm [7] can be used in our context. However, the complexity of computing a CAD is doubly exponential in the number of variables ([4, 9]); its implementations are usually limited to non-trivial problems involving 4 variables and cannot tackle our application.

The complexity of computing a CAD can be much improved when taking into account equational constraints (see e.g. [22]). In the context of real root classification problems, this leads us to take advantage of the presence of equations to compute closed sets in the parameter space ($\mathbb{R}^f$ using our notation) containing the boundaries of the regions $C_1, \dots, C_f$, hence substituting the recursive (doubly exponential) projection steps of CAD with more involved projection techniques. In the past ten years, several works have focused on this problem [20, 26] using various computer algebra tools such as Gröbner bases, regular chains, etc. We also mention [24] which uses evaluation/interpolation techniques to compute those closed sets in the parameters space.

While the implementation of [26] is able to solve our classification problem for the case of water, none of the implementations were able to classify the singular locus of $V$ in the general case (the number of parameters is 3).

Our strategy is not as general as the ones in [26] or [20]. It exploits properties of sets defined by minors of matrices with polynomial entries. Such structures have been used for computing sample points in each connected component of the real trace of determinantal varieties [16, 17, 18] or for solving linear matrix inequalities [15]. These works are based on dedicated strategies for computing critical loci of some projections restricted to determinantal varieties. Such computations are naturally related to real root classification problems and real quantifier elimination (see e.g. [19]). Finally, our computations rely on Gröbner bases. Several works have shown some connection between Gröbner bases and determinantal ideals [14] and critical point computation [13, 25].

**Main results.** Our main results are twofold:

* an algorithm solving the special real root classification problem described above and which exploits the determinantal structure of the input data arising in contrast imaging problems;

* its successful use for solving the challenging application to the contrast problem in the general case;

* answers to the questions raised by the experimental data involving water: the answer to question 1. is no, and the answer to question 2. is that there are 3 classes of second matters that we can separate, depending on whether there are 1, 2 or 3 singularities.

We start by describing our algorithmic contribution. Recall that we are given a matrix, denoted by $M$, with polynomial entries. As in [15, 16, 17, 18], it is based on splitting computations according to the rank of $M$.

More precisely, in order to solve our real root classification problem, we need to identify where the number of real solutions inside $B$ of the determinantal system describing $V$ changes depending on the values of the parameters. The real root classification problem we want to solve involves inequalities defining a semi-algebraic set with non-empty interior. In this context, we use standard tools from real geometry, such as Thom's first isotopy lemma, which reduce our classification problem to computing the singular points of $V$, the critical points of the projection of the parameter space restricted to $V$, and the intersection of $V$ with the boundary of the semi-algebraic set $B$.

This computation may be difficult because generically, the variety $V_r$ has singularities corresponding to points where $rank(M) < r$: this is proved using Bertini's theorem and [5, Prop. 1.1], as in [16, Prop. 2]). Hence, observe that the variety $V$ is naturally split according to the rank of $M$. This is the very basic idea on which our algorithm relies: we compute critical loci of the projection on the parameter space restricted to the variety $V$ by distinguishing those points at which $M$ has rank less than $r$ from those at which $M$ has rank exactly $r$. Incidentally, it raises the question of how these higher rank deficiencies should be interpreted from the application viewpoint. To the best of our knowledge, this question was never raised in the optimal control community.

Our algorithms need to compute projection of algebraic sets, which is done using elimination algorithms such as Gröbner bases or triangular sets for example. We have performed experiments for both these tools, using the package FGb [12] in Maple and an implementation of F5 [11] for Gröbner bases, and using the package RegularChains [21] in Maple for triangular sets.

Regarding the contrast imaging problem, we illustrate the behaviour of our algorithm in the case of water, giving the whole classification. Using Gröbner bases to perform the eliminations, the computation takes 10 s on an 2 GHz Intel Xeon CPU. The RealRootClassification command of the Maple RegularChain library needs 1600 s to find this classification.

We also ran our algorithm on the general case. While none of the available implementation is able to tackle this classification problem directly, ours can find the polynomials separating the open sets $C_i$ within 4 h using FGb, or 2 min using F5, and the projection step of the CAD can be done in 4 h. We see similar speed-ups when using triangular sets to perform the elimination.

This illustrates how our dedicated algorithms take advantage of the special structure of the problem, to achieve speed-ups when compared with more general techniques.

**Conclusion.** We propose an algorithmic strategy refining general roots classification strategies for the case of singularities of a determinantal variety, under genericity hypotheses. We give an overview of the results for the application at the end of the paper<fnref n="1" />. We were able to give a full classification in the case of water, answering the questions raised by the experiments. For the general case, the separating polynomials were found, it remains to perform the analysis of the subdivisions in order to obtain the full classification. This work also raised questions concerning the interpretation of higher rank deficiencies from the viewpoint of control theory.

**Structure of the paper.** In Section 2, we present the mathematical background around NMR imagery and the contrast problem. Sec-

<footnote n="1">¹The full results and the source code which produced them are available at http://mercurey.gforge.inria.fr/</footnote>
tion 3 deals with the dedicated classification algorithm. Finally, in Section 4, we report on experimental results obtained when solving the application to the contrast imaging problem.

We acknowledge the support from the ANR-DFG research program Explosys (Grant No.ANR-14-CE35-0013-01; GL203/9-1). Mohab Safey El Din is member of and supported by Institut Universitaire de France.

## 2. MODELING THE DYNAMICS

The model we describe below has been introduced in [3] in order to apply techniques from geometric optimal control theory to the control of the spin dynamics by NMR. Up to some normalization, each spin 1/2 particle is governed by the Bloch equation

$$ egin{cases} \dot{x} = -\Gamma_x u_y z \\ \dot{y} = -\Gamma_y - u_x z \\ \dot{z} = \gamma(1-z) + u_x y - u_y x, \end{cases} $$

where the state variable $q = (x, y, z)$ represents the magnetization vector which must lie in the Bloch ball defined by $|q| \le 1$, and the parameters $(\Gamma, \gamma)$ are related to the physical relaxation times. The parameters must also satisfy $2\Gamma \ge \gamma > 0$. The control $u = (u_x, u_y)$ represents the magnetic field whose magnitude is bounded by a maximum value $\mu$.

In the context of the contrast imaging problem, this leads to the simultaneous control of two non-interacting spins with different relaxation time parameters. The contrast by saturation method consists in bringing the magnetization vector of the first spin toward the center of the Bloch ball while maximizing the modulus of the magnetization vector of the other matter. The matter with a zero magnetization is black on the picture, while the other matter with a maximum modulus of the magnetization vector is bright.

Using the symmetry of revolution [2] which allows to eliminate one state variable for each matter, we obtain the system

$$ egin{cases} \dot{y}_1 = -\Gamma_1 y_1 - u_x z_1 \\ \dot{z}_1 = \gamma (1-z_1) + u_x y_1 \\ \dot{y}_2 = -\Gamma_2 y_2 - u_x z_2 \\ \dot{z}_2 = \gamma_2 (1-z_2) + u_x y_2, \end{cases} \quad |u| \le \mu \quad (1) $$

and the optimal control problem is: starting from the equilibrium point $N = ((0,1), (0,1))$, saturate the first spin, that is $q_1(T) = 0$, where $T$ is the transfer time while maximizing $|q_2(T)|^2$, where $|q_2(T)|$ represents the final contrast. It is a standard Mayer problem in optimal control, studied in [3] through the analysis of the Hamiltonian dynamics given by the Pontryagin Maximum Principle [23]. We summarize this analysis below.

Writing (1) as $\dot{q} = F(q) + uG(q)$, $|u| \le \mu$, the optimality conditions associated with the Maximum Principle lead us to construct the optimal solution as a concatenation of bang-arcs where the control is $u = \pm \mu$, and singular arcs solutions of $X_e = F + u_s G$ where the control $u_s$ is the rational fraction $-D'/D$ with

$$ D = \det(F, G, [G, F], [[G, F], G]) \\ D' = \det(F, G, [G, F], [[G, F], F]), $$

where $[\cdot, \cdot]$ denotes the Lie bracket of vector fields. Explicitly, with $d_i = \gamma_i - \Gamma_i (i \in \{1, 2\})$,

$$ D = \det \begin{bmatrix} -\Gamma_1 y_1 & -z_1 - 1 & d_1 z_1 - \Gamma_1 & 2 d_1 y_1 \\ -\gamma_1 z_1 & y_1 & d_1 y_1 & -2 d_1 z_1 + \Gamma_1 - d_1 \\ -\Gamma_2 y_2 & -z_2 - 1 & d_2 z_2 - \Gamma_2 & 2 d_2 y_2 \\ -\gamma_2 z_2 & y_2 & d_2 y_2 & -2 d_2 z_2 + \Gamma_2 - d_2 \end{bmatrix}. $$

The localization of the singularities of $\{D=0\}$ inside the Bloch ball is important to understand the geometry of the hypersurface, as well as the dynamics of the vector field $X_e$ which is closely linked to the presence of such singularities. Indeed, generically when approaching the surface $D=0$ along a singular arc, the con-

trol $u = -D'/D$ is such that $|u| > \mu$. Hence the control policy switches from singular to bang. Therefore this surface is related to the complexity of the optimal law, as a concatenation of singular and bang arcs. See [3] for numerical simulations related to this phenomenon.

# 3. ALGORITHM

## 3.1 Classification strategy

We consider the polynomial algebra $\mathbb{Q}[\mathbf{X}, \mathbf{G}]$ with variables $\mathbf{X} = (X_1, \dots, X_n)$ and parameters $\mathbf{G} = (G_1, \dots, G_t)$. Let $F$ and $H$ be families of polynomials in $\mathbb{Q}[\mathbf{X}, \mathbf{G}]$. Let $V_{\mathbb{R}} = V_{\mathbb{R}}(F)$, $V = V_{\mathbb{C}}(F)$ be respectively the set of zeroes of $F$ in $\mathbb{R}^{n+t}$ and in $\mathbb{C}^{n+t}$. Let $B$ be the closed semi-algebraic set defined by $H$:

$$ B = \{ (\mathbf{x}, \mathbf{g}) \in \mathbb{R}^{n+t} \mid \forall h \in H, h(\mathbf{x}, \mathbf{g}) \le 0 \}, $$

and let $B_0 = \cup_{h \in H} V_{\mathbb{C}}(h)$. Let $\pi: \mathbb{C}^{n+t} \to \mathbb{C}^t$ be the projection onto the affine space with coordinates $\mathbf{G}$. Let sing($V$) be the singular locus of $V$, crit($\pi, V$) be the set of critical points of $\pi$ restricted to $V$, and $K(\pi, V) = \pi(\text{sing}(V) \cup \text{crit}(\pi, V)) \cap \mathbb{R}^t$.

Given a subset $A$ of a real or complex affine space, $\bar{A}$ and $\partial A$ are used to denote the closure and the boundary of $A$ for the Euclidean topology respectively.

Assume that the following hypotheses are satisfied:

H1 There exists a nonempty Zariski-open subset $O_1$ of $\mathbb{C}^t$ such that for all $\mathbf{g} \in O_1$, the fiber $V \cap \pi^{-1}(\mathbf{g})$ is a nonempty finite subset of $\mathbb{C}^{n+t}$;

H2 The restriction of the projection $\pi$ to $B$ is proper ([10, Def. 2.10.1]);

H3 The intersection $V \cap B_0$ has dimension at most $t-1$ in $\mathbb{C}^{n+t}$;

H4 The variety $V$ is equidimensional of dimension $t$.

We want to find a non-zero polynomial $P \in \mathbb{Q}[\mathbf{G}]$ such that, on each connected component $U$ of $\mathbb{R}^t \setminus V_{\mathbb{R}}(P)$, for $\mathbf{g} \in U$, the cardinality of $V \cap B \cap \pi^{-1}(\mathbf{g})$ does not depend on $\mathbf{g}$.

In Lemma 1, we describe a well-known strategy for computing these objects (see for example [20, 26]).

**LEMMA 1.** Let $F$ and $H$ be polynomial systems satisfying hypotheses H1, H2, H3 and H4. Let $C_B = \pi(V \cap B_0)$, $U$ a non-empty connected open subset of $\mathbb{R}^t$ which does not meet $C_B \cup K(\pi, V)$, and $\mathbf{g} \in U$. Then $V \cap \pi^{-1}(\mathbf{g})$ is finite, and $\forall \mathbf{g}' \in U$, $\#(V \cap \pi^{-1}(\mathbf{g}')) = \#(V \cap \pi^{-1}(\mathbf{g}))$.

**PROOF.** We will construct a Whitney stratification of $V \cap B([1, Def. 9.7.1])$ with certain properties. First note that since $V$ is $t$-equidimensional by H4, $V$ has real dimension at most $t$ ([1, Prop. 2.8.2]). Let $S_{=t}$ be the intersection of the points where $V_{\mathbb{R}}$ has local dimension $t$ and of the interior of $B$. There exists a Whitney stratification $(S_i)$ of the semi-algebraic set $V \cap B$ such that $S_{=t}$ is the union of strata of dimension $t$ ([1, Th. 9.7.11]). Let $S_{\le t}$ be the union of the other strata, they all have real dimension less than $t$. By construction, this is a semi-algebraic set which is the union of $(V \cap \partial B) \subset (V \cap B_0)$ and of the singular locus of $V \cap B$, and it has dimension less than $t$ (by H3 for $V \cap B_0$). Its image through $\pi$ has dimension less than $t$, and so it has codimension at least 1.

Now consider $S_{=t}$. By Hyp. H1, for any $\mathbf{g} \in O_1$, $\pi^{-1}(\mathbf{g}) \cap V$ is non-empty, hence $\pi(V)$ contains the non-empty Zariski-open set $O_1 \subset \mathbb{C}^t$. The intersection $O_1 \cap \mathbb{R}^t$ is a non-empty Zariski-open set of $\mathbb{R}^t$, contained in $\pi(V)$, hence $\pi(V) \cap \mathbb{R}^t$ has real dimension $t$. Let $U_0$ be its interior. The subset $S_{=t} \cap \pi^{-1}(U_0)$ is a locally closed semi-algebraic set. If it is empty, then there is nothing to prove. Otherwise, by construction it has dimension $t$; and the projection $\pi$ restricted to this subspace is proper, by Hyp. H2. Thom's
isotopy lemma ([8]) states that for any nonempty connected open set $U$ of $\mathbb{R}^t$ not meeting $K(\pi, V)$, and for any $g \in U$, there exists a semi-algebraic diffeomorphism $h = (h_0, \pi) : V \cap B \cap \pi^{-1}(U) \to \pi^{-1}(g) \times U$. By Hyp. $H1$, if $U$ is nonempty, $\pi^{-1}(g)$ is finite, and the cardinality of the fibers is constant on $U$. $\square$

So computing the wanted decomposition of the parameter space is equivalent to computing a polynomial $P \in \mathbb{Q}[\mathbf{G}]$ such that $V(P)$ covers $\pi(V \cap B_0)$ and $K(\pi, V)$.

## 3.2 The determinantal problem

Let $k$ be an integer greater than 1, $r_0 \in \{1, \dots, k-1\}$, and $n = (k-r_0+1)^2$. Let $t \in \mathbb{N}$, and let $M(\mathbf{X}, \mathbf{G})$ be a $k \times k$ matrix with polynomial entries in $n$ variables $\mathbf{X} = (X_1, \dots, X_n)$ and $t$ parameters $\mathbf{G} = (G_1, \dots, G_t)$. As before, let $\pi: \mathbb{C}^{n+t} \to \mathbb{C}^t$ be the projection onto the affine space with coordinates $\mathbf{G}$.

Let $\{h(\mathbf{X}, \mathbf{G}) \le 0 \mid h \in H\}$ be a system of inequalities, with $H \subset \mathbb{Q}[\mathbf{X}, \mathbf{G}]$. Let $V_{-1} = \emptyset$ by convention, and for any $r \in \{0, \dots, k\}$, we define the variety

$$V_r = \{ (\mathbf{x}, \mathbf{g}) \in \mathbb{C}^{n+t} \mid \mathrm{rank}(M(\mathbf{x}, \mathbf{g})) \le r \}$$

and the constructible set $V_r = V_r \setminus V_{r-1}$, that is the set of points at which the matrix $M$ has rank exactly $r$.

Let $V$ be the union of the singular locus of $V_{r_0}$ and of the set of critical points of $\pi$ restricted to $V_{r_0}$. We want to classify the cardinality of the real fibers by $\pi$ of the semi-algebraic set $V \cap \{(\mathbf{x}, \mathbf{g}) \mid \forall h \in H, h(\mathbf{x}, \mathbf{g}) \le 0\}$.

Assume that $V$ and $H$ satisfy hypotheses $H1$, $H2$, $H3$ and $H4$. Further assume that:

H5 There exists a non-empty Zariski-open subset $O_2 \subset \mathbb{C}^t$ such that $V \cap \pi^{-1}(O_2) = V_{r_0-1} \cap \pi^{-1}(O_2)$;

H6 For any $r \in \{0, \dots, k-1\}$, the ideal defined by the $(r+1)$-minors of $M$ is radical;

H7 For any $r \in \{0, \dots, k-1\}$, the variety $V_r$ is equidimensional with dimension $n+t-(k-r+1)^2$.

These properties are generic ([16, Prop. 2]).

**LEMMA 2.** *Assuming Hyp. H6, for any $r \in \{1, \dots, k\}$, $V_{r-1} \subset \mathrm{sing}(V_r)$.*

**PROOF.** We will prove the following stronger statement: let $(\mathbf{x}, \mathbf{g}) \in V_{r-1}$, then all partial derivatives of all $(r+1)$-minors of $M$ vanish at $(\mathbf{x}, \mathbf{g})$. This will prove that the Jacobian of $V_r$ at $(\mathbf{x}, \mathbf{g})$ is the zero matrix, and in particular has rank 0. By Hyp. H6 the ideal of all $(r+1)$-minors of $M$ is radical, so we can use the Jacobian criterion to characterize $\mathrm{sing}(V)$, so $(\mathbf{x}, \mathbf{g}) \in \mathrm{sing}(V_r)$.

Let $(\mathbf{x}, \mathbf{g}) \in V_{r-1}$. For any $(r+1) \times (r+1)$-submatrix of $M$, the result we want to prove depends only on the coefficients of the sub-matrix. So w.l.o.g., we may assume that $r=k-1$. Let $(\mathbf{x}, \mathbf{g}) \in V_{k-2}$, this means that all $(k-1)$-minors of $M$ vanish at $(\mathbf{x}, \mathbf{g})$. Consider a matrix of polynomial indeterminates $U$:

$$\tilde{M} = \begin{bmatrix} U_{1,1} & \dots & U_{1,k} \\ \vdots & \ddots & \vdots \\ \dot{U}_{k,1} & \dots & \dot{U}_{k,k} \end{bmatrix}$$

and let $\tilde{D}$ be its determinant. Then for any $i, j \in \{1, \dots, k\}$, $\frac{\partial \tilde{D}}{\partial U_{i,j}} = (-1)^{i+j} \tilde{M}_{i,j}(U)$ where $\tilde{M}_{i,j}$ is the $(k-1)$-minor of $\tilde{M}$ obtained by removing row $i$ and column $j$.

For any $i, j \in \{1, \dots, k\}$, let $m_{i,j}$ (resp. $M_{i,j}$) be the coefficient at row $i$ and column $j$ of the matrix $M$ (resp. the minor obtained by removing row $i$ and column $j$ from $M$). By the derivation chain rule, for any

$v \in \{X_1, \dots, X_n, G_1, \dots, G_t\}$, $\frac{\partial D}{\partial v} = \sum_{i,j=1}^k (-1)^{i+j} \frac{\partial m_{i,j}}{\partial v} \cdot \tilde{M}_{i,j}(\mathbf{m})$, which equals $\sum_{i,j=1}^k (-1)^{i+j} \frac{\partial m_{i,j}}{\partial v} \cdot M_{i,j}(\mathbf{X})$. Since by hypothesis all $(k-1)$ minors of $M$ vanish at $(\mathbf{x}, \mathbf{g})$, all partial derivatives $\frac{\partial D}{\partial v}$ vanish at $(\mathbf{x}, \mathbf{g})$. $\square$

In the following subsections, we will describe two algorithms DeterminantCritVals and DeterminantBoundary, which, given such a matrix $M$, a target rank $r_0$ and inequalities $H$, compute a polynomial whose zeroes cover $K(\pi, V)$, and a polynomial whose zeroes cover $\pi(V \cap B_0)$ respectively. By Lemma 1, the zeroes of the product of these polynomials will subdivide the parameter space into connected components where the cardinality of real fibers is constant. These algorithms are probabilistic, because they will rely on the choice of generic linear forms to ensure linear independence. However, the algorithms could be made deterministic by testing that these linear forms are generic enough for our purpose, and repeating the random draw otherwise.

The algorithms will also need to compute the projection of algebraic sets onto coordinate subspaces. For this purpose, we assume that we are given a routine Elimination, which, given a system of polynomials

$F \subset \mathbb{Q}[V_1, \dots, V_N]$ and a set of variables $V' \subset \{V_1, \dots, V_N\}$, computes a system of generators of $\langle F \rangle \cap \mathbb{Q}[V']$. Such a routine can be implemented using Gröbner bases or regular chains, for example.

## 3.3 Incidence varieties

We decompose the problem depending on the rank of the matrix. The classical technique that we use to model properties on the rank relies on incidence varieties.

**DEFINITION 3.** Let $r \in \{0, \dots, k-1\}$. The incidence variety of rank $r$ associated with $M$ is the variety $V_r \subset \mathbb{C}^{n+t} \times (\mathbb{P}^{k-1}(\mathbb{C}))^{k-r}$ defined by:

$$M = \begin{bmatrix} y_{1,1} & \dots & y_{1,k-r} \\ \vdots & \ddots & \vdots \\ y_{k,1} & \dots & y_{k,k-r} \end{bmatrix} = \begin{bmatrix} 0 & \dots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \dots & 0 \end{bmatrix} \quad (2)$$

with the additional condition that the matrix $(y_{i,j})$ has rank $k-r$.

The projection of $V_r$ onto the affine space with coordinates $(\mathbf{X}, \mathbf{G})$ is $V_r$. Let $(u_{1,1}, \dots, u_{k-r,k}) \in \mathbb{C}^{k(k-r)}$, we define the variety $V'_{r,u}$ as the intersection of $V_r$ and the complex solutions of

$$\begin{bmatrix} u_{1,1} & \dots & u_{1,k} \\ \vdots & \ddots & \vdots \\ u_{k-r,1} & \dots & u_{k-r,k} \end{bmatrix} \cdot \begin{bmatrix} y_{1,1} & \dots & y_{1,k-r} \\ \vdots & \ddots & \vdots \\ y_{k,1} & \dots & y_{k,k-r} \end{bmatrix} = \mathrm{Id}_{k-r} \quad (3)$$

In the rest of Section 3, $\mathcal{F}_{r,u}$ denotes the union of Equations (2) and (3).

**LEMMA 4.** *For any $r \in \{0, \dots, k-1\}$, the varieties $V_r$ and $V'_{r,u}$ are birational ([1, Sec. 3.4]).*

**PROOF.** We need to define a morphism $f: W \to V'_{r,u}$ with $W$ a non-empty Zariski-open subset of $V_r$, and such that $f$ is inverse to the projection $V'_{r,u} \to V_r$ onto the affine space with coordinates $(\mathbf{X}, \mathbf{G})$. Let $W$ be the open subset of $V_r$ defined as the non-vanishing locus of the top-left $r$-minor of $M$. Consider the block decompositions, where $A, Y_{(1)}$ and $U_{(1)}$ are $r \times r$ matrices:

$$M = \begin{bmatrix} A & B \\ C & D \end{bmatrix} \quad Y = \begin{bmatrix} Y_{(1)} \\ Y_{(2)} \end{bmatrix} \quad U = [U_{(1)} \quad U_{(2)}]$$

* Over $W$, $A$ is invertible, let $\Delta = \det(A)$. Let $M/A$ be the Schur complement of $A$ in $M$, Equations (2) and (3) can be rewritten

$$\begin{bmatrix} \Delta \mathrm{Id}_r & A^{-1}B \\ 0 & M/A \\ U_{(1)} & U_{(2)} \end{bmatrix} \cdot \begin{bmatrix} Y_{(1)} \\ Y_{(2)} \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ \mathrm{Id}_r \end{bmatrix}$$
We may restrict to the open subset of $V'_{r,u}$ where $Y_{(2)}$ is invertible, then eliminating $Y_{(1)}$ yields that

$$ egin{cases} Y_{(2)} = (U_{(2)} - U_{(1)}A^{-1}B)^{-1} \\ Y_{(1)} = \frac{-1}{\Delta}A^{-1}BY_{(2)} \end{cases} $$

which defines the wanted morphism $W \to V'_{r,u}$. $\square$

**PROPOSITION 5.** Let $r_1 = r_0 - 1$. Let $\phi: V'_{r_1,u} \to \mathbb{C}^t$ be the projection onto the affine space with coordinates $\mathbf{G}$. Assuming that hypotheses H1 to H7 hold, there exists a Zariski-open subset $\mathcal{U} \subset \mathbb{C}^{k(k-r_1)}$ such that if $\mathbf{u} \in \mathcal{U} \cap \mathbb{Q}^{k(k-r_1)}$, $K(\pi, V_{r_1}) = K(\phi, V'_{r_1,u})$.

**PROOF.** Let $P = (\mathbf{x}, \mathbf{g}, \mathbf{y}) \in V'_{r_1,u}$.

If $M(\mathbf{x}, \mathbf{g})$ has rank less than $r_1$, then by Lemma 2, $(\mathbf{x}, \mathbf{g}) \in \text{sing}(V_{r_1})$, hence $\mathbf{g} \in K(\pi, V_{r_1})$.

Since $M(\mathbf{x}, \mathbf{g})$ has rank less than $r_1$, its kernel $L_1$ has dimension at least $k - r_1 + 1$. Equations (3) encode that the vectors $y_i$ given by the columns of matrix $(y_{i,j})$ generate a $r_1$-dimensional vector space $L_2$. So there exists $\mathbf{y}_0 \in L_1 \cap L_2$, and for all $a \in \mathbb{C}$, $(\mathbf{x}, \mathbf{g}, \mathbf{y}_1 + \mathbf{a}\mathbf{y}_0, \mathbf{y}_2)$ belongs to the fiber above $(\mathbf{x}, \mathbf{g})$ in $V_{r_1}$. So this fiber has dimension at least 1, while the generic fiber has dimension 0 by hypothesis H1. So $(\mathbf{x}, \mathbf{g})$ is a critical value of the projection of $V_{r_1}$ onto $\mathbb{R}^{n+t}$, hence $(\mathbf{g}) \in K(\phi, V'_{r_1,u})$.

So we may assume that $M(\mathbf{x}, \mathbf{g})$ has rank exactly $r_1$. There is a $r_1 \times r_1$ submatrix $A$ of $M(\mathbf{x}, \mathbf{g})$ which is invertible, without loss of generality we may assume that it is the top-left $r_1 \times r_1$ submatrix. In an open neighborhood of $(\mathbf{x}, \mathbf{g})$, $V=r_1$ is described by the vanishing of the entries of $M/A$, that is the determinants of the $(r_1+1) \times (r_1+1)$ submatrices containing $A$. The same computations as in the proof of Lemma 4 give the following equations describing $V'_{r_1,u}$ in the open neighborhood of $(\mathbf{x}, \mathbf{g}, \mathbf{y})$ where $\Delta = \det(A)$ does not vanish:

$$ egin{cases} M/A = 0 \\ Y_{(2)} = (U_{(2)} - U_{(1)}A^{-1}B)^{-1} \\ Y_{(1)} = \frac{-1}{\Delta}A^{-1}BY_{(2)} \end{cases} \quad (4) $$

and the truncated Jacobian matrix in $(\mathbf{X}, \mathbf{Y})$ of this system can be written

$$ \begin{bmatrix} \text{Jac}_{\mathbf{X}}(M/A) & 0 & 0 \\ \star & \text{Id}_{r_1(k-r_1)} & \star \\ \star & 0 & \text{Id}_{r_1(k-r_1)} \end{bmatrix} $$

where $\text{Jac}_{\mathbf{X}}(M/A)$ is the truncated Jacobian matrix in $\mathbf{X}$ of the $(k-r_1)^2$ entries of $M/A$, which define $V_{r_1} \setminus \{(\mathbf{x}, \mathbf{g} | \Delta = 0)\} \in \mathbb{C}^{n+t}$. By hypothesis H6, the ideal defined by the entries of $M/A$, which is a subideal of the ideal of all $(r_1+1)$-minors of $M$, is radical. Since the Schur complement appears by multiplication with invertible matrices with entries in the localized ring $\mathbb{Q}[\mathbf{X}, \mathbf{g}]_{\Delta}$ (using the same notations as in the proof of Lemma 4):

$$ \begin{bmatrix} \text{Id}_{r_1} & 0 \\ -C & \text{Id}_{k-r_1} \end{bmatrix} \cdot \begin{bmatrix} A^{-1} & 0 \\ 0 & \text{Id}_{k-r_1} \end{bmatrix} \cdot M = \begin{bmatrix} \text{Id}_{r_1} & A^{-1}B \\ 0 & M/A \end{bmatrix}, $$

Equations (4) describe the localization of $\langle \mathcal{F}_{r_1,u} \rangle$ in $\mathbb{Q}[\mathbf{X}, \mathbf{g}]_{\Delta}$, so this ideal is radical as well. So we can use the Jacobian criterion on $V_{r_1}$ near $(\mathbf{x}, \mathbf{g})$ and on $V'_{r_1,u}$ near $(\mathbf{x}, \mathbf{g}, \mathbf{y})$. Both Jacobians matrices have the same rank and both varieties have the same local codimension $(k-r_1)^2$ (by Lemma 4 and hypothesis H7), so

$$ K(\pi, V_{r_1}) \cap \phi(V'_{r_1,u}) = K(\phi, V'_{r_1,u}) \cap \pi(V_{r_1}) $$

The image $\phi(V'_{r_1,u}) \cap \pi(V_{r_1})$ is a Zariski-open subset $\mathcal{O}_{\mathbf{u}}$ of $\pi(V_{r_1})$. It remains to prove that if $\mathbf{u}$ is sufficiently generic, then all irreducible components of $K(\pi, V_{r_1})$ meet this open subset.

Let $C_1, \dots, C_a$ be these irreducible components, and let $(\mathbf{x}_1, \mathbf{g}_1) \in \pi^{-1}(C_1), \dots, (\mathbf{x}_a, \mathbf{g}_a) \in \pi^{-1}(C_a)$. For any $(\mathbf{x}, \mathbf{g}) \in V_{r_1}$, the proof of [16, Prop. 4, Sec. 6] shows that there exists a non-empty Zariski-open subset $\mathcal{U}_{(\mathbf{x}, \mathbf{g})} \subset \mathbb{C}^{k(k-r_1)}$ such that if $\mathbf{u} \in \mathcal{U}_{(\mathbf{x}, \mathbf{g})} \cap \mathbb{Q}^{k(k-r_1)}$, then

($\mathbf{x}, \mathbf{g}) \in \mathcal{O}_{\mathbf{u}}$; namely, $\mathcal{U}$ is the set of $\mathbf{u}$ such that rank $\begin{bmatrix} M(\mathbf{x}, \mathbf{g}) \\ (u_{i,j}) \end{bmatrix} = k$. Taking the finite intersection of the non-empty Zariski-open subsets $\mathcal{U}_{(\mathbf{x}_i, \mathbf{g}_i)}$ for $i \in \{1, \dots, a\}$ yields the wanted subset $\mathcal{U}$. $\square$

## 3.4 Locus of rank exactly $r_0$

Recall that by $\mathcal{H}_5$, $\pi(V \cap V_{r_0})$ has codimension at least 1, and that we want to compute a polynomial whose zeroes cover $K(\pi, V)$ and $\pi(V \cap B_0)$. So we may multiply the result by the equation of one hypersurface covering $\pi(V \cap V_{r_0})$, it will naturally cover $\pi(V \cap V_{r_0}) \cap K(\pi, V)$ and $\pi(V \cap V_{r_0}) \cap \pi(V \cap B_0)$.

**Algorithm RankExactly:**

Input $M \subset \mathbb{Q}[\mathbf{X}, \mathbf{G}]^{k \times k}$, $r_0 \in \{1, \dots, k-1\}$

Output $P_1 \in \mathbb{Q}[\mathbf{G}] \setminus \{0\}$ s.t. $\pi(V \cap V_{r_0}) \subset V(P_1)$

**Procedure**

1. res $\leftarrow 1$

2. $F_{V,0} \leftarrow \{(r_0+1)-\text{minors of } M\}$

3. $J \leftarrow \text{Jac}_{\mathbf{X}}(F_{V,0})$

4. $F_{V,1} \leftarrow F_{V,0} \cup \{(k-r_0)^2-\text{minors of } J\}$

5. Pick at random $u_1, \dots, u_{k(k-r_0)} = \mathbf{u} \in \mathbb{Q}^{k(k-r_0)}$

6. $F_0 \leftarrow \mathcal{F}_{k-r_0, \mathbf{u}}$

7. $\{M_1, \dots, M_N\} \leftarrow \{r_0-\text{minors of } M\}$

8. For $i \in \{1, \dots, N\}$ do

9. $F_1 \leftarrow F_0 \cup \mathcal{F}_{V,1} \cup \{M_1, \dots, M_{i-1}, u_i \cdot M_i - 1\}$

10. $G \leftarrow \text{Elimination}(F_1, \{u, \mathbf{X}, \mathbf{Y}\})$

11. Multiply res by 1 polynomial from $G$

12. End for

13. Return res

## 3.5 Singularities

**Algorithm DeterminantCritVals:**

Input $M \subset \mathbb{Q}[\mathbf{X}, \mathbf{G}]^{k \times k}$, $r_0 \in \{1, \dots, k-1\}$

Output $P_c \in \mathbb{Q}[\mathbf{G}] \setminus \{0\}$ s.t. $K(\pi, V) \subset V(P_c)$

**Procedure**

1. res $\leftarrow \text{RankExactly}(M, r_0)$

2. Pick at random $u_1, \dots, u_{k(k-r_0+1)} = \mathbf{u} \in \mathbb{Q}^{k(k-r_0+1)}$

3. $F_0 \leftarrow \mathcal{F}_{r_0-1, \mathbf{u}}$

4. $J \leftarrow \text{Jac}_{\mathbf{X}, \mathbf{Y}}(F_0)$

5. $F_1 \leftarrow F_0 \cup \{k(k-r_0+1) + (k-r_0+1)^2-\text{minors of } J\}$

6. $G \leftarrow \text{Elimination}(F_1, \{\mathbf{X}, \mathbf{Y}\})$

7. Multiply res by 1 polynomial from $G$

8. Return res

**PROPOSITION 6.** *Algorithm DeterminantCritVals is correct.*

**PROOF.** By definition, $V \subset V_{r_0}$. Using the decomposition $V_{r_0} = V_{r_0} \cup V_{r_0-1}$, we decompose the variety $V$ as $V = (V \cap V_{r_0}) \cup (V \cap V_{r_0-1})$.

The subspace $\pi(V \cap V_{r_0})$ is covered by the output of RankExactly, so we may restrict to $V \cap V_{r_0-1}$, which is the whole variety $V_{r_0-1}$ by Lemma 2. Recall that by hypothesis $\mathcal{H}7$, $V_{r_0-1}$ is t-equidimensional.

By Prop. 5, in order to compute $K(\pi, V_{r_0-1})$, we can compute polynomials whose zeroes cover $K(\phi, V'_{r_0-1, \mathbf{u}})$ with $\mathbf{u}$ sufficiently generic instead.

By hypotheses $\mathcal{H}6, \mathcal{H}7$ and the proof of Prop. 5, $V'_{r_0-1, \mathbf{u}}$ is t-equidimensional and $\mathcal{F}_{r_0-1, \mathbf{u}}$ is a set of generators of its ideal, so we can use the Jacobian criterion to compute equations defining $K(\phi, V'_{r_0-1, \mathbf{u}})$. $\square$

## 3.6 Boundary

**Algorithm DeterminantBoundary:**

**Input**

• $M \subset \mathbb{Q}[\mathbf{X}, \mathbf{G}]^{k \times k}$
• $r_0 \in \{1, \dots, k-1\}$

• $H \subset \mathbb{Q}[\mathbf{X}, \mathbf{G}]$ (set of constraints on the variables)

**Output** $P_b \in \mathbb{Q}[\mathbf{G}] \setminus \{0\}$ s.t. $\pi(V \cap B_0) \subset V(P_b)$

**Procedure**

1. res $\leftarrow$ RankExactly($M$, $r_0$)

2. Pick at random $u_1, \dots, u_{k(k-r_0+1)} = \mathbf{u} \in \mathbb{Q}^{k(k-r_0+1)}$

3. $F_0 \leftarrow \mathcal{F}_{r_0-1, \mathbf{u}}$

4. For $h$ in $H$ do

5. $F_1 \leftarrow F_0 \cup \{h\}$

6. $G \leftarrow$ Elimination($F$, $\{\mathbf{X}, \mathbf{Y}\}$)

7. Multiply res by 1 polynomial from $G$

8. End for

9. Return res

**PROPOSITION 7.** *Algorithm DeterminantBoundary is correct.*

**PROOF.** As in Section 3.5, we write: $V = (V \cap V_{=r_0}) \cup (V \cap V_{r_0-1})$. Since $B_0 = \cup_{h \in H} V(h)$, the intersection $V \cap \partial B$ is contained in the union of the varieties $V(\langle F \rangle + \langle h \rangle)$ for $h$ ranging over $H$, and the equation of the projections can be obtained with polynomial elimination. $\square$

**REMARK 8.** *For the real root classification problem, the subdivision is given by the product of the outputs of DeterminantCritical and DeterminantBoundary. In order to avoid repeating computations, we may skip the call to RankExactly in either subroutine (but not both), and initialize res to 1 instead.*

# 4. THE CONTRAST PROBLEM

## 4.1 The case of water

With the notations of Section 2, the variety $V$ is the complex algebraic variety defined by

$$D = \frac{\partial D}{\partial y_1} = \frac{\partial D}{\partial y_2} = \frac{\partial D}{\partial z_1} = \frac{\partial D}{\partial z_2} = 0.$$

With the notations of Section 3, we want to classify the singularities of the set of points where $M$ has rank at most $r_0 = 3$. Our semi-algebraic constraints are that the solutions are within the Bloch ball, that is

$$\mathcal{B}: \begin{cases} h_1 = y_1^2 + (z_1 + 1)^2 \le 1 \\ h_2 = y_2^2 + (z_2 + 1)^2 \le 1. \end{cases}$$

Since the equations are homogeneous in $\Gamma_1, \Gamma_2, \gamma_1, \gamma_2$, and the parameters are supposed to be non-zero, we may normalize by setting $\gamma_1 = 1$. In the case where the first matter is water, we further simplify by setting $\Gamma_1 = \gamma_1 = 1$, leaving free the two parameters $\Gamma_2, \gamma_2$ corresponding to the second matter. We recall that we also assume that $2\Gamma_2 \ge \gamma_2$ and that $(\gamma_2, \Gamma_2) \ne (1, 1) = (\gamma_1, \Gamma_1)$ (that is, the second matter is not water).

**THEOREM 9.** *Consider the 9 polynomials:*

$$f_1 = \Gamma_2 - 1, \quad f_2 = 3\Gamma_2 - 2\gamma_2 - 1,$$

$$f_3 = 3\Gamma_2^2 - 5\Gamma_2\gamma_2 + \gamma_2^2 + 2\Gamma_2 - 2\gamma_2 + 1,$$

$$f_4 = 2\Gamma_2^2 - 5\Gamma_2\gamma_2 + 2\gamma_2^2 - 2\Gamma_2 + 3\gamma_2,$$

$$f_5 = 2\gamma_2^2 - (3\Gamma_2 + 11)\gamma_2^2 + (9\Gamma_2 + 6 - 3\Gamma_2^2)\gamma_2 + 2\Gamma_2(\Gamma_2 + 2)(\Gamma_2 - 1),$$

$$f_6 = \Gamma_2 - 2\gamma_2 + 1, \quad f_7 = 2\Gamma_2 - \gamma_2 - 1,$$

$$f_8 = \gamma_2 - 2 + \Gamma_2, \quad f_9 = 2\Gamma_2^2 - 5\Gamma_2\gamma_2 + 2\gamma_2^2 + 1.$$

The zeroes of their product divide the subset of $\mathbb{R}^2$ defined by $2\Gamma_2 > \gamma_2 > 0$ into connected components where the cardinality of $V_{\mathbb{R}} \cap \pi^{-1}(\gamma_2, \Gamma_2)$ is constant.

<figure><img src="image_2.png" /><figcaption>**Figure 1:** Curves involved in the definition of the semi-algebraic set $G$. The blue (resp. green) sample points correspond to points in $\mathcal{G}_1^- \cup \mathcal{G}_1^+$ (resp. $\mathcal{G}_2^- \cup \mathcal{G}_2^+$). The circled numbers in each area correspond to the number of singularities in $B$ for parameters in the area. Parameters in the red area are physically irrelevant.</figcaption></figure>

Let $\psi: (y_1, z_1, y_2, z_2) \mapsto (-y_1, z_1, -y_2, z_2)$ be the symmetry fixing $\Pi = \{y_1 = y_2 = 0\}$, and let us consider the semi-algebraic sets (see Fig. 1):

$$\begin{align*} \mathcal{G}_1^- &= \{0 < \gamma_2 < 2\Gamma_2, \Gamma_2 < 1, f_2 > 0\}, \\ \mathcal{G}_1^+ &= \{0 < \gamma_2 < 2\Gamma_2, \Gamma_2 > 1, f_2 < 0, f_4 > 0\}, \\ \mathcal{G}_2^- &= \{0 < \Gamma_2 < 1, f_6 > 0, f_3 < 0\}, \\ \mathcal{G}_2^+ &= \{\Gamma_2 > 1, f_6 < 0, f_5 > 0\}, \\ \mathcal{G} &= \mathcal{G}_1^- \cup \mathcal{G}_1^+ \cup \mathcal{G}_2^- \cup \mathcal{G}_2^+. \end{align*}$$

**THEOREM 10.** *For all $(\gamma_2, \Gamma_2)$ such that $2\Gamma_2 > \gamma_2 > 0$, the center of the Bloch ball $\mathcal{B}$ is a singularity of $\{D=0\}$. And provided $(\gamma_2, \Gamma_2) \in \mathcal{G}$, there exist at most two other singularities:*

1. provided $(\gamma_2, \Gamma_2) \in \mathcal{G}_1^- \cup \mathcal{G}_1^+$ there is one other singularity lying on $\Pi \cap \mathcal{B}$;

2. provided $(\gamma_2, \Gamma_2) \in \mathcal{G}_2^- \cup \mathcal{G}_2^+$ there are two other singularities in $\mathcal{B}$, $\psi$-symmetric, outside $\Pi$.

The configuration is illustrated in Figs. 1 and 2. Observe that the number of singularities inside $\mathcal{B}$ is an invariant of the contrast problem. Two of the pairs of biological matters studied in [2], water-cerebrospinal fluid (normalized parameters $[\gamma_2 = \frac{5}{4}, \Gamma_2 = \frac{25}{3}]$) and water-fat (normalized parameters $[\gamma_2 = \frac{25}{2}, \Gamma_2 = 25]$) correspond to points outside $\mathcal{G}$, and their invariant is 1 in both cases (see Fig. 2). But our results give answers to our guiding questions: there exist pairs of matters for which this algebraic invariant can differ; and any pair (water, matter) belongs to one of 3 classes, depending on whether the number of singularities inside $\mathcal{B}$ is 1, 2 or 3.

**PROOF OF THEOREM 9.** Let $V_3 = \{p \in \mathbb{C}^4 \times \mathbb{R}^2 | \operatorname{rank}(M) = 3\}$ and $V_2 = \{p \in \mathbb{C}^4 \times \mathbb{R}^2 | \operatorname{rank}(M) < 3\}$, where $p = (y_1, y_2, z_1, z_2, \gamma_2, \Gamma_2)$. We apply the strategies described in Section 3.

We study the generic case $V_2 \cap V$ first. This set does cover a dense subset of $\mathbb{R}^2$. Its intersection with the boundary of $\mathcal{B}$ is given by the vanishing of either $h_1$ or $h_2$. The projection on $(\Gamma_2, \gamma_2)$ of the set of points of $V_2 \cap V$ such that $h_1 = 0$ is described by $0 = \gamma_2 f_1^2 f_2 f_3$ which gives us polynomials $f_1, f_2$ and $f_3$.

The projection on $(\Gamma_2, \gamma_2)$ of the set of points of $V_2 \cap V$ such that $h_2 = 0$ is described by $0 = (2\Gamma_2 - \gamma_2) f_1^2 f_4 f_5$ which gives us new polynomials $f_4$ and $f_5$.

Next, we consider the incidence variety $\mathcal{V}_2$ associated with the matrix $M$:
<table><caption>Table 1: Timings</caption><thead><tr><th>System</th><th>RegularChain (direct)</th><th>Gröbner (direct)</th><th>RegularChain (new algo.)</th><th>FGb (new algo.)</th><th>F5 (new algo.)</th><th>CAD</th></tr></thead><tbody><tr><td>Water</td><td>1600 s</td><td>100 s</td><td></td><td>10 s</td><td>1 s</td><td>50 s</td></tr><tr><td>General</td><td>&gt;24 h</td><td>&gt;24 h</td><td>90 × 200 s</td><td>46 × 200 s</td><td>110 s</td><td>4 h (projection step)</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2: Positions of the parameters corresponding to the pairs water-cerebrospinal fluid (red circle) and water-fat (red square) and the set $\mathcal{G}$ (with the same conventions as in Fig. 1). For both these pairs, there is only 1 singularity in $\mathcal{B}$.</figcaption></figure>

$$M \cdot \begin{pmatrix} \lambda_{1,1} & \lambda_{1,2} \\ \lambda_{2,1} & \lambda_{2,2} \\ \lambda_{3,1} & \lambda_{3,2} \\ \lambda_{4,1} & \lambda_{4,2} \end{pmatrix} = \begin{pmatrix} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \end{pmatrix}$$

with random linear equations ensuring that the matrix $(\lambda_{i,j})$ has rank 2.

Out of the surface $\gamma_2 = 0$, this affine variety is a complete intersection (it has dimension 2 and it is given by 9 equations in 11 variables, including the saturation by $\gamma_2$). The set of critical values of $\pi$ is described by $0 = (2\Gamma_2 - \gamma_2)(\Gamma_2 + 1)f_1^2 f_6^2 f_7^2$ which gives us new polynomials $f_6$ and $f_7$ ($\Gamma_2 + 1$ has no solutions within our constraint range).

This completes the study of $V \cap V_2$. We now move on to the study of $V \cap V_3$. As described in the algorithm, we define the incidence variety of rank 3 of $M$, and we saturate successively by the 3-minors of $M$. Only the first of these subcases is nonempty, and it is described by $0 = (2\Gamma_2 - \gamma_2)f_8f_9$ which gives us $f_8$ and $f_9$. $\square$

**PROOF OF THEOREM 10.** Observe first by means of a trivial evaluation that $O$ is a singularity of $\{D=0\}$. We now focus on singularities in $B^* = B \setminus \{O\}$. Theorem 9 provides a list of 9 polynomials to which we add our constraints $2\Gamma_2 \ge \gamma_2 > 0$. Let $\xi = \gamma_2\Gamma_2(\gamma_2 - 2\Gamma_2)\prod_{i=1}^9 f_i$. The complementary of $\{\xi = 0\}$ is the union of a sequence of connected open semi-algebraic sets where the number of singularities is constant. The routine CylindricalAlgebraDecompose of the Maple package RegularChains[Semi-AlgebraicSetTools] provides 1533 sample points. Excluding those at which $\xi$ vanishes and those outside our physical constraints domain, remains a set $K_c$ of 570 points. At each point of $K_c$ we locate the singularities by computing a Gröbner basis.

We get 187 points of $K_c$ such that there exists at least one singularity in $B^*$. We have a set $K_s$ of 31 points, each of them corresponding to a couple of $\psi$-symmetric singularities outside the symmetry plane $\Pi$, and a set $K_p$ of 156 points corresponding to a unique singularity on $\Pi \cap B^*$. For parameters at which $\xi$ does not vanish, the number of singularities in $B^*$ is at most two.

<figure><img src="image_4.png" /><figcaption>Figure 3: The curves involved in the decomposition of the parameter space (with the same conventions as in Fig. 1).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Magnification of Fig. 3 near (1, 1).</figcaption></figure>

Points of $K_s$ (resp. $K_p$) are represented in green (resp. blue) in Figs. 1, 3 and 4. Let us evaluate on $K_c$ the condition $(\Gamma_2 < 1, f_2 > 0, f_4 < 0)$ or $(\Gamma_2 > 1, f_2 < 0, f_4 > 0)$. Indeed the set of points of $K_c$ satisfying this condition coincides with $K_p$. This proves item 1). The proof of item 2) is similar. $\square$

## 4.2 The general case

The variety $V$ and the semi-algebraic set $B$ are defined as in the previous section. We normalize again by $\gamma_1 = 1$, we assume that $2\Gamma_1 \ge 1, 2\Gamma_2 \ge \gamma_2 > 0, (\gamma_2, \Gamma_2) \ne (1, \Gamma_1)$, and that $\Gamma_1 \ne 1, \Gamma_2 \ne \gamma_2$ (case of water).

**THEOREM 11.** Splitting the subset of $\mathbb{R}^3$ defined by $2\Gamma_2 > \gamma_2 > 0$ and $2\Gamma_1 > 1$ into open subsets where the number of real singularities of $V$ in the fibers is constant, can be done by cutting out 12 irreducible surfaces, consisting of 5 planes, 3 quadrics, two surfaces of degree 9 and one of degree 14.

These polynomials were obtained by applying the algorithms from Section 3 to our system. The elimination steps were done using both Gröbner bases with FGb or with F5, and with triangular sets with RegularChains. Table 1 presents some timings for these
methods (for computations done with interpolation, we give the results as $a \times b$ where a is the interpolation degree and b the time taken for each specialized computation). It shows that our algorithm, implemented with either Gröbner bases or regular chains, is more efficient than the direct (general) approach, and it allows to deal with the previously untractable general case. Furthermore, the implementation using Gröbner bases appears to be faster.

## 5. REFERENCES

[1] J. Bochnak, M. Coste, and M.-F. Roy. *Real algebraic geometry*, volume 36 of *Ergebnisse der Mathematik und ihrer Grenzgebiete (3)* [*Results in Mathematics and Related Areas (3)*]. Springer-Verlag, Berlin, 1998. Translated from the 1987 French original, Revised by the authors.

[2] B. Bonnard, M. Chyba, A. Jacquemard, and J. Marriott. Algebraic geometric classification of the singular flow in the contrast imaging problem in nuclear magnetic resonance. *Math. Control Relat. Fields*, 3(4):397–432, 2013.

[3] B. Bonnard, O. Cots, S.J. Glaser, M. Lapert, D. Sugny, and Y. Zhang. Geometric optimal control of the contrast imaging problem in nuclear magnetic resonance. *IEEE Trans. Automat. Control*, 57(8):1957–1969, 2012.

[4] C. W. Brown and J. H. Davenport. The complexity of quantifier elimination and cylindrical algebraic decomposition. In *ISSAC 2007*, pages 54–60. ACM, New York, 2007.

[5] W. Bruns and U. Vetter. *Determinantal Rings*. Springer, 2014.

[6] M. Chyba, J-M. Coron, P. Gabriel, A. Jacquemard, G. Patterson, G. Picot, and P. Shang. Optimal geometric control applied to the protein misfolding cyclic amplification process. *Acta Applicandae Mathematicae*, 135(1):145–173, 2015.

[7] G. E. Collins. Quantifier elimination for real closed fields by cylindrical algebraic decomposition. In *Automata theory and formal languages (Second GI Conf., Kaiserslautern, 1975)*, pages 134–183. Lecture Notes in Comput. Sci., Vol. 33. Springer, Berlin, 1975.

[8] M. Coste and M. Shiota. Thom's first isotopy lemma: a semialgebraic version, with uniform bound. In *Real analytic and algebraic geometry (Trento, 1992)*, pages 83–101. de Gruyter, Berlin, 1995.

[9] J. H. Davenport and J. Heintz. Real quantifier elimination is doubly exponential. *J. Symbolic Comput.*, 5(1-2):29–35, 1988.

[10] M. Demazure. *Bifurcations and catastrophes*. Universitext. Springer-Verlag, Berlin, 2000. Geometry of solutions to nonlinear problems, Translated from the 1989 French original by D. Chillingworth.

[11] J.-C. Faugère. A new efficient algorithm for computing Gröbner bases without reduction to zero ($F_5$). In *Proceedings of the 2002 International Symposium on Symbolic and Algebraic Computation*, pages 75–83 (electronic), New York, 2002. ACM.

[12] J.-C. Faugère. FGb: A Library for Computing Gröbner Bases. In Komei Fukuda, Joris Hoeven, Michael Joswig, and Nobuki Takayama, editors, *Mathematical Software - ICMS 2010*, volume 6327 of *Lecture Notes in Computer Science*, pages 84–87, Berlin, Heidelberg, September 2010. Springer Berlin / Heidelberg.

[13] J.-C. Faugère, M. Safey El Din, and P.-J. Spaenlehauer. Critical points and Gröbner bases: the unmixed case. In

*Proceedings of the 2012 International Symposium on Symbolic and Algebraic Computation (ISSAC 2012)*, pages 162–169, 2012.

[14] J.-C. Faugère, M. Safey El Din, and P.-J. Spaenlehauer. On the complexity of the Generalized MinRank Problem. *Journal of Symbolic Computation*, 55:30–58, 2013.

[15] D. Henrion, S. Naldi, and M. Safey El Din. Exact algorithms for linear matrix inequalities. Submitted., August 2015.

[16] D. Henrion, S. Naldi, and M. Safey El Din. Real root finding for low rank linear matrices. Submitted., June 2015.

[17] D. Henrion, S. Naldi, and M. Safey El Din. Real root finding for rank defects in linear Hankel matrices. In *Proceedings of the 2015 ACM on International Symposium on Symbolic and Algebraic Computation*, page 26p., Bath, United Kingdom, July 2015.

[18] D. Henrion, S. Naldi, and M. Safey El Din. Real root finding for determinants of linear matrices. *Journal of Symbolic Computation*, 74:205–238, 2016.

[19] H. Hong and M. Safey El Din. Variant quantifier elimination. *Journal of Symbolic Computation*, 47(7):883–901, 2012. International Symposium on Symbolic and Algebraic Computation (ISSAC 2009).

[20] D. Lazard and F. Rouillier. Solving parametric polynomial systems. *Journal of Symbolic Computation*, 42(6):636–667, 2007.

[21] F. Lemaire, M. Moreno Maza, and Y. Xie. The regularchains library. In *Maple conference*, volume 5, pages 355–368, 2005.

[22] S. McCallum. On projection in CAD-based quantifier elimination with equational constraint. In *Proceedings of the 1999 International Symposium on Symbolic and Algebraic Computation (Vancouver, BC)*, pages 145–149 (electronic). ACM, New York, 1999.

[23] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko. *The mathematical theory of optimal processes*. Translated from the Russian by K. N. Trirogoff; edited by L. W. Neustadt. Interscience Publishers John Wiley & Sons, Inc. New York-London, 1962.

[24] J. I. Rodriguez and X. Tang. Data-discriminants of likelihood equations. In *Proceedings of the 2015 ACM on International Symposium on Symbolic and Algebraic Computation*, ISSAC '15, pages 307–314, New York, NY, USA, 2015. ACM.

[25] P.-J. Spaenlehauer. On the complexity of computing critical points with Gröbner bases. *SIAM Journal on Optimization*, 24(3):1382–1401, 2014.

[26] L. Yang, X. Hou, and B. Xia. A complete algorithm for automated discovering of a class of inequality-type theorems. *Sci. China Ser. F*, 44(1):33–49, 2001.