Multivariate Bernstein Subdivision solver

{{J.-P. Pavone}}

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 $M_j (f ; x_j) = 0$ in $[a_j, b_j]$ or $a_j$ (resp. $b_j$) if $m_j (f ; x_j) = 0$ (resp. $M_j (f ; x_j) = 0$) has no root on $[a_j, b_j]$,
  • $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

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

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:

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 $M$ behaves like Newton iteration near a simple root, that is we have a quadratic convergence.

can be considered.

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

Implementation details

See synaps/msolve/sbdslv.h.



#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]).