Multivariate Bernstein Subdivision solver |
We consider here the problem of computing the solutions of a polynomial system
in a box 
. 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 
of degree 
in the variable
, can be decomposed as:
where 
is the tensor product Bernstein basis on
the domain and 
are the
control coefficients of 
on 
.
The polynomial is represented in this basis by
the 
order tensor of control coefficients 
. The size of , denoted by
, is 
.
De Casteljau algorithm also applies in each of the direction , , 
so that we can split this
representation in these directions. We use the following properties to
isolate the roots:
and 
, let
,
and any , we have
As a direct consequence, we obtain the following corollary:
of the equation 
in the domain 
, we have 
where
(resp. 
) is
either a root of 
or 
in 
or 
(resp.
) if (resp.
) has no root on ,
on 
.
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.
Preconditioners.
Solving the system
is equivalent to solving the system 
,
where 
is an 
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
, we transform locally the system 
into a system 
, where 
is the Jacobian matrix of
at a point 
of the domain ,
where it is invertible. A direct computation shows that locally
(in a neighborhood of 
the level-set of
are orthogonal to
the -axes. We can prove that the reduction
based on the polynomial bounds 
and behaves like Newton iteration near a simple
root, that is we have a quadratic convergence.
Reduction strategy
Here are several reduction strategies,, that can be considered.-
Convex hull reduction: A method called Interval Projected Polyhedron (or IPP) can be used 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 
defined in
corollary ?.
The guarantee that the computed intervals contain the root of
, 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.
Subdivision strategy
Here some simple rules that can be used to 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 
(or the difference between the
control coefficients) is maximal. Then, a value 
for splitting the domain in the direction 
,
is chosen
-
either where
has a local maximum,
-
or where
(1) The right-hand side of this equation can be easily computed, from the sum of all the control coefficients of
.
-