SYNAPS (SYmbolic Numeric ApplicationS) » Multivariate Solvers

Multivariate Bernstein Subdivision solver

We consider here the problem of computing the solutions of a polynomial system

$\left\{ \begin{array}{l} f_1 (x_1, \ldots, x_n) = 0\\ \vdots\\ f_s (x_1, \ldots, x_n) = 0 \end{array} \right.$

in a box $I : = [a_1, b_1] \times \cdots \times [a_n, b_n] \subset \mathbbm{R}^n$ . This solver uses the representation of multivariate polynomials in the Bernstein basis, analysis of sign variations and univariate solvers to localise the real roots of a polynomial system. The output is a set of small-enough boxes, which may contain these roots.

By a direct extension to the multivariate case, any polynomial $f (x_1, \ldots, x_n) \in \mathbbm{R} [x_1, \ldots, x_n]$ of degree $d_i$ in the variable $x_i$ , can be decomposed as:

$f (x_1, \ldots, x_n) = \sum_{i_1 = 0}^{d_1} \cdots \sum_{i_n = 0}^{d_n} b_{i_1, \ldots, i_n} \hspace{0.25em} B_{d_1}^{i_1} (x_1 ; a_1, b_1) \cdots B_{d_n}^{i_n} x (x_n ; a_n, b_n) .$

where $( \hspace{0.25em} B_{d_1}^{i_1} (x_1 ; a_1, b_1) \cdots B_{d_n}^{i_n} (x_n ; a_n, b_n))_{0 \le i_1 \le d_1, \ldots, 0 \le i_n \le d_n}$ is the tensor product Bernstein basis on the domain $I : = [a_1, b_1] \times \cdots \times [a_n, b_n] \subset \mathbbm{R}^n$ and $\mathbf{b} = (b_{i_1, \ldots, i_n})_{0 \le i_1 \le d_1, \ldots, 0 \le i_n \le d_n}$ are the control coefficients of $f$ on $I$ . The polynomial $f$ is represented in this basis by the $n^{\text{\tmtextrm{th}}}$ order tensor of control coefficients $\mathbf{b} = (b_{i_1, \ldots, i_n})_{0 \le i \le d_1, 0 \le j \le d_2, 0 \le k \le d_3}$ . The size of $I$ , denoted by $|I|$ , is $|I| = \max \{|b_i - a_i | ; i = 1, \ldots, n\}$ .

De Casteljau algorithm also applies in each of the direction $x_i$ , , $i = 1, \ldots, n$ so that we can split this representation in these directions. We use the following properties to isolate the roots:

Definition: For any $f \in \mathbbm{R} [ \mathbf{x}]$ and $j = 1, \ldots, n$ , let
$\begin{eqnarray*} m_j (f ; x_j) & = & \sum_{i_j = 0}^{d_j} \min_{\{0 \le i_k \le d_k, k \neq j\}} b_{i_1, \ldots, i_n} \hspace{0.25em} B^{i_j}_{d_j} (x_j ; a_j, b_j)\\ M_j (f ; x_j) & = & \sum_{i_j = 0}^{d_j} \max_{\{0 \le i_k \le d_k, k \neq j\}} b_{i_1, \ldots, i_n} x \hspace{0.25em} B^{i_j}_{d_j} (x_j ; a_j, b_j) . \end{eqnarray*}$

Theorem: [Projection Lemma] For any $\mathbf{u} = (u_1, \ldots, u_n) \in I$ , and any $j = 1, \ldots, n$ , we have
$m (f ; u_j) \le f ( \mathbf{u}) \le M (f ; u_j) .$

As a direct consequence, we obtain the following corollary:

Corollary: For any root $\mathbf{u} = (u_1, \ldots, u_n)$ of the equation $f ( \mathbf{x}) = 0$ in the domain $I$ , we have $\underline{\mu}_j \le u_j \le \overline{\mu}_j$ where

$\underline{\mu}_j$ (resp. $\overline{\mu}_j$ ) is either a root of $m_j (f ; x_j) = 0$ or in $[a_j, b_j]$ or $a_j$ (resp. $b_j$ ) if (resp. ) has no root on ,

$m_j (f ; u) \le 0 \le M_j (f ; u)$ on $[ \underline{\mu}_j, \overline{\mu}_j]$ .

The solver implementation contains the following main steps. It consists in

[1] applying a preconditioning step on the equations;

[2] in reducing the domain;

[3] if the reduction ratio is too small, to split the domain

until the size of the domain is smaller than a given epsilon.

As we are going to see, we have several options for each of these steps, leading to different algorithms with different behaviors, as we will see in the next sections. Indeed the solvers that we will consider are parameterized by the

{Preconditioner}: A transformation of the initial system into a system, which has a better numerical behavior.

{Reduction strategy}: The technique used to reduce the initial domain, for searching the roots of the system.

{Subdivision strategy}: The technique used to subdivide the domain, in order to simplify the forthcoming steps, for searching of roots of the system.

equivalent to solving the system $M \hspace{0.25em} \mathbf{f} = 0$ , where $M$ is an $s \times s$ invertible matrix

As such a transformation may increase the degree of some equations, with respect to some variables, it has a cost, which might not be negligible in some cases.

Moreover, if for each polynomials of the system not all the variables are involved, that is if the systems is sparse with respect to the variables, such a preconditioner may transform it into a system which is not sparse anymore. In this case, we would prefer a partial preconditioner on a subsets of the equations sharing a subset of variables.

The following preconditioners are curently avialable:

{Global transformation:} We minimise the distance between the equations, considered as vectors in an affine space of polynomials of a given degree.

{Local straightening:} In this part we consider square systems, for which .

If we are ``near'' a simple solution of the system $\mathbf{f} = 0$ , we transform locally the system $\mathbf{f}$ into a system $J^{- 1} \mathbf{f}$ , where $J = (\partial_{x_i} f_j ( \mathbf{u}_0)_{1 \le i, j \le s}$ is the Jacobian matrix of $\mathbf{f}$ at a point $\mathbf{u}$ of the domain $I$

, where it is invertible. A direct computation shows that locally (in a neighborhood of $\mathbf{u}_0$ the level-set of $f_i$

$i = 1, \ldots, s$ are orthogonal to the $x_i$

-axes: This transformation was also discussed in {{gs-isbaspe-01}}. We can prove that the reduction based on the polynomial bounds $m$

and

behaves like Newton iteration near a simple root, that is we have a quadratic convergence.

can be considered.

{Convex hull reduction}: In {{sp-csnps-93}}, a method called Interval Projected Polyhedron (or IPP) is described, in order to reduce the domain of search. It is based on the property that the convex hull property. A new reduced domain is computed, by intersecting the convex hull of the projected set of control points. With our notations, in considering the control polygons defining and instead of these polynomials.

{Extreme root reduction}: A direct improvement of the convex hull reduction consist in computing the first (resp. last) root of the polynomial , (resp. ), in the interval . The current implementation of this reduction steps allows us to consider the convex hull reduction, as one iteration step of this reduction process.

{Univariate solver reduction} Here, we compute all the roots of the polynomials and and keep the intervals $[ \underline{\mu}, \overline{\mu}]$ defined in corollary [*] .

The guarantee that the computed intervals contain the root of $f$

, is achieved by controlling the rounding mode of the operations during the de Casteljau computation.

This method will be particularly interesting in the cases where more than one interval have to be considered. This may happen at the beginning of the search but is less expected locally near a single root.

subdivide a domain. We will show in the next section their impact on the performance of the solver

{Parameter bisection}: The domain is then split in half in a direction for which is maximal.

{Image bisection}: Instead of choosing the size of the interval as a criterion for the direction in which we split, we may choose a and a such that $\delta_j (f_k ; u_j) = M_j (f_k ; u_j) - m_j (f_k ; u_k)$ (or the difference between the control coefficients) is maximal. Then, a value for splitting the domain in the direction , is chosen

either where $\delta_j (f_k ; u_j)$ has a local maximum,

or where
$\begin{equation} \int_{v = 0}^u \delta_j (f_k ; v) dv = \frac{1}{2} \int_{v = 0}^1 \delta_j (f_k ; v) dv \end{equation}$
The right-hand side of this equation can be easily computed, from the sum of all the control coefficients of $\delta_j (f_k ; v)$ .

Implementation details

See synaps/msolve/sbdslv.h.

Example

#include <synaps/init.h>
#include <synaps/mpol.h>
#include <synaps/base/Seq.h>
#include <synaps/msolve/sbdslv.h>

typedef MPol<double> mpol_t;

int main (int argc, char **argv)
{
  using std::cout; using std::endl;
  
  mpol_t p("x0^2*x1-x0*x1-1"), q("x0^2-x1^3-2");

  std::list<mpol_t> I; I.push_back(p); I.push_back(q);

  VectDse<double> dmn(4,"-2 2 -2 2");
  //  cout << solve(I, SBDSLV<double,SBDSLV_RDL>())<<endl;
  cout << solve(I, SBDSLV<double,SBDSLV_RDL>(1e-3),dmn)<<endl;

  //| [[1.66068,0.911523], [-1.4228,0.290084]
  //
  // The result is a sequence of vectors with coefficients of double type.
  // The size of the vectors is the number of variables.
  // Here there are no solutions in the default [0,1]x[0,1] and
  // 2 real solutions in the box [-2,2]x[-2,2]).

}

SYNAPS DOCUMENTATION