considers the fundamental problem of finding the complex
solutions of a system
of polynomial equations in complex variables
. The system is triangular in the sense
that for ,
As in , we assume
that the system is regular: this means that
for each , if is a zero of
and is the leading coefficient
of in , then .
Thus is a -dimensional system.
But unlike , we
do not assume that the system is square-free:
indeed the goal of this paper is to demonstrate
new techniques that can properly
determine the multiplicity of the root clusters of ,
up to any resolution.
this paper, we use boldface symbols to denote vectors and tuples; for
We are interested in finding clusters of solutions of triangular systems
and in counting the total multiplicity of solutions in clusters.
Solving triangular systems is a fundamental task in polynomial equations
solving, since there are many algebraic techniques to
decompose the original system into triangular systems.
The problem of isolating the complex solutions of a polynomial system
in an initial region-of-interest (ROI) is defined as follows: let
denote the set of solutions of in
, regarded as
Local Isolation Problem (LIP):
Given: a polynomial map , a polybox ,
Output: a set of
pairwise disjoint polydiscs of radius where
- each is a
This is “local” because we restrict attention to roots in a ROI
. There are two issues with (LIP) as formulated above:
deciding if is a singleton, and deciding if
such a singleton lies in , are two “zero problems” that
require exact computation. Generally, this can only be decided if
is algebraic. Even in the algebraic case, this may be very
expensive. In [27, 4] these
two issues are side-stepped by defining the local clustering problem
which is described next.
Before proceeding, we fix some general notations for this paper. A
polydisc is a vector of
complex discs. The center of is the vector of the
centers of its components and the radius
is the vector of the radii
of its components. If is any positive real number, we denote
by the polydisc
that has the same center
as and radius .
We also say if each component of
is . A (square complex) box is
a complex interval where
and ; the width
of is and the center of is
is the set which is represented by the
vector of boxes.
The center of is the vector of the centers of its
components; the width of is the max of
the widths of its components.
If is any positive real number, we denote by
the polybox that has the same
center than and width .
It is also convenient to identify as the subset
of ; a similar remark applies
We introduce three notions to define the local solution clustering
problem. Let be a solution of
The multiplicity of in , also called the
intersection multiplicity of in is classically
defined by localization of rings as in [28, Def. 1,
p. 61], we denote it by . An
equivalent definition uses dual spaces, see [12, Def. 1,
For any set , we denote by the
multiset of zeros (i.e., solutions) of in , and
the total multiplicity of .
If is a polydisc, we
call a cluster if it is non-empty, and is
an isolator of the cluster. If in addition, we have that
, we call
a natural cluster and call a natural
In the context of numerical algorithm, the notion of cluster of
solutions is more meaningful than that of solution with multiplicity
since the perturbation of a multiple solution generates a cluster. We
thus “soften” the problem of isolating the solutions of a triangular
system of polynomial equations while counting their multiplicities by
translating it into the local solution clustering problem defined as
Local Clustering Problem (LCP):
Given: a polynomial map , a polybox ,
Output: a set of pairs
- the s are pairwise disjoint
polydiscs of radius ,
In this (LCP) reformulation of (LIP), we have removed the two “zero
problems” noted above: we output clusters to avoid the first problem,
and we allow the output to contain zeroes outside the ROI to
avoid the second one. We choose for simplicity; it is easy
to replace the factor of by for any desired .
In the remaining of this section we explain our
contribution, summarize previous work and the local univariate
clustering method of . In
Sec. 2, we define the notion of tower
of clusters together with a recursive method to compute the sum of
multiplicities of the solutions it contains.
Sec. 3 analyzes the loss of precision induced by
Our algorithm for solving
the local clustering problem for triangular systems is
introduced in Sec. 4. The implementation and
experimental results are presented in Sec. 5.
1.1 Our contributions
We propose an algorithm for solving the complex clustering
problem for a triangular system with a
zero-dimensional solution set. To this end, we propose a formula to
count the sum of multiplicities of solutions in a cluster.
Our formula is derived from a result of  that
links the intersection multiplicity of a solution of a triangular
system to multiplicities in fibers.
We define towers of clusters to encode clusters of solutions
of a triangular system in stacks (or towers) of clusters of roots of
Our algorithm exploits the triangular form of :
the standard idea is to recursively find roots
of the form of
, then substituting them into
to obtain a univariate polynomial
If is a root of , then we
have have extended the solution to
of the original .
The challenge is to extend this idea to compute clusters
of zeros of from clusters of
zeros of .
Moreover, we want to allow
the coefficients of each to be oracle numbers.
The use of oracle numbers allows us to treat polynomial
systems whose coefficients are algebraic numbers and beyond.
To compute clusters of roots of a univariate polynomial given as an
oracle, we rely on the recent algorithm described in
based on a predicate introduced in  that combines
Pellet’s theorem and Graeffe iterations to determine the number of
roots counted with multiplicities in a complex disc; this predicate is
called soft because it only requires the polynomial to be known
as approximations. It is used in a subdivision framework combined with
Newton iterations to achieve a near optimal complexity.
We implemented our algorithm and made it available
as the Julia
package Ccluster.jl .
Our experiments show that it advantageously compares
to major homotopy solvers
for solving random dense triangular systems
in terms of solving times and reliability
(i.e. getting the correct number of solutions
and the correct multiplicity structures).
Homotopy solving is more general because
it deals with any polynomial system.
We also propose experiments with
triangular systems obtained with elimination
1.2 Related work
There is a vast literature on solving polynomial systems and we can only refer
to book surveys and references therein, see for instance
[13, 25]. On the algebraic side,
symbolic tools like Groebner basis, resultant, rational univariate
parametrization or triangularization,
find an equivalent triangular system or set
thus reducing the problem to the univariate case.
these methods handle all input, in particular with solutions with
multiplicities, and are certified but at the price of a high complexity that
limits their use in practice.
Implementations of hybrid symbolic-numeric solvers are available
for instance in Singular
or in Maple via RootFinding[Isolate].
On the numerical side, one can find subdivision and homotopy methods. The main
advantage of subdivision methods is their locality: the practical complexity
depends on the size of the solving domain and the number of solutions in this
domain. Their main drawback is that they are only practical for low dimensional
systems. On the other hand, homotopy methods are efficient for high dimensional
systems, they are not local but solutions are computed independently from one
Numerical methods only work for restricted classes of systems and the
certification of the output remains a challenge. Multiprecision arithmetic, interval
analysis, deflation and -theory are now classical tools to address this
certification issue [15, 22, 6, 26].
In the univariate case, practical certified algorithms are now available for
real and complex solving that match the best known complexity bounds together
with efficient implementations [19, 17].
For the bivariate case, the problem of solving a triangular system can be seen
as a univariate isolation in an extension field. The most recent contributions
in this direction presenting algorithms together with complexity analysis are
Only a few work address the specific problem of solving triangular polynomial
systems. The solving can then be performed coordinate by coordinate by
specialization and univariate solving in fibers. When the systems only have
regular solutions, extensions of classical univariate isolation algorithms to
polynomial with interval coefficients have been proposed
[9, 14, 7]. In the presence of multiple
solutions, one approach is to use a symbolic preprocessing to further decompose
the system in regular sub-systems. Another approach is the sleeve method with
separation bounds . The authors of 
propose a formula to compute the multiplicity of a solution of a triangular
system: the latter multiplicity is the product of the multiplicities of the
components of a solution in the fibers. Then, by using square free factorization
of univariate polynomials specialized in fibers, they describe an algorithm to
retrieve the real solutions of a triangular system with their multiplicities.
In , the method of Local Generic Position is adapted to the
special case of triangular systems with the advantage of only using resultant
computations (instead of Goebner basis), multiplicities are also computed.
1.3 Definitions and Notation
Convention for Vectors.
We introduce some general conventions for vectors that will
simplify the following development.
Vectors are indicated by bold fonts. If
is an -vector, and , then the -th component
the -vector is denoted .
Thus , and
“” is an idiomatic way
of saying that is an -vector.
Because of the subscript convention, we
will superscripts such as , etc, to distinguish
among a set of related -vectors.
Normed Vector Spaces.
In order to do error analysis, we need to treat
and as normed vector spaces:
for and ,
and denote the infinity norm
on polynomials and vectors, respectively.
We use the following perturbation convention:
let . Then we will write
to denote some polynomial
that satisfies .
Similarly, denotes some vector
that satisfies .
If then and
are called -bit approximations
of and , respectively.
We define the degree sequence of
to be where is the degree of in .
If (), let
denote the polynomial that results
from the substitution
The result is a polynomial
specialization of by .
Note that is a polynomial in at most variables.
In particular, when , then is a constant
(called the evaluation of at ).
For instance, suppose ,
then is a polynomial in
If is a box
with center and width ,
we denote by the disc with center and
radius . Note that contains .
If is a polybox,
let be the polydisc
Oracle Computational Model.
We use two kinds of numbers in our algorithms:
an explicit kind which is standard in computing,
and an implicit kind which we call “oracles”.
Our explicit numbers are dyadic numbers (i.e., bigFloats),
. A pair
of integers represents the nominal value of .
However, we also want this pair to represent the interval
. To distinguish between them,
we write for the nominal value, and
for the interval of width .
Call an -bit dyadic interval
if (so the interval has width at most ).
Note that and are different despite having the
same nominal value. As another example, note that
is the interval .
When we say a box, disc, polybox, etc, is dyadic, it means that
all its parameters are given by dyadic numbers.
The set of closed intervals with dyadic endpoints is
denoted . Also, let denote
-dimensional dyadic boxes.
The implicit numbers in our algorithms are functions: for any
real number , an oracle for is a function
such that is an -bit
dyadic interval containing . There is normally no confusion in
identifying the real number with any
oracle function for . Moreover, we write
instead of .
E.g., if is a real algebraic number with defining polynomial
and isolating interval , we may define
an oracle for in a fairly standard way.
Next, an oracle for a complex number
is a function such that
are oracles for and . Again, we may
identify with any oracle , and write
instead of . For polynomials
in variables, we assume a sparse representation,
with fixed support ,
with coefficients , and
are power products.
An oracle for amounts to having oracles for
each coefficient of . Moreover
may be written
is the interval polynomial whose coefficients are .
Call a dyadic interval polynomial.
1.4 Oracles for Root Cluster of Univariate Polynomials
The starting point for this paper is the fundamental result
that the Local Clustering Problem (LCP) has been
solved in the univariate setting:
Proposition 1 (See [4, 5])
There is an algorithm
that solves the Local Clustering Problem
when is a univariate oracle polynomial.
In other words, the output of
is a set
such that each is a natural -isolator, and
To make this result the basis of our multivariate
clustering algorithm, we need to generalize this result.
In particular, we need to be able to further
refine each output of this algorithm.
If represents the cluster of roots,
we want to get better approximation of , i.e.,
we want to treat like number oracles.
Fortunately, the algorithm in [4, 5]
already contains the tools to do this. What is lacking is
a conceptual framework to capture this.
Our goal is to extend the concept of number oracles
to “cluster oracles”. To support the several modifications
which are needed, we revise our previous view of “oracles as
functions”. We now think of an oracle
as a computational object with state information,
and which can transform itself in order to
update its state information.
For any , recall that is
a dyadic object that is at least -bit accurate.
E.g., if is the oracle for ,
is an interval containing of width .
But now, we say that oracle is transformed to a
new oracle which we shall denote by “”
whose state information is . In general,
let denote the state information in .
for a cluster of roots of a univariable polynomial
, its oracle has state
that is a pair where is a dyadic disc
and is the total multiplicity of .
Thus is automatically a natural cluster.
We say is -bit accurate if the radius
of is at most .
Intuitively, we expect to be
an oracle for that is -bit accurate.
Unfortunately, this may be impossible unless is a singleton
cluster. In general, we may have to split into two
or more clusters. We therefore need one more extension:
the map returns a set
of cluster oracles with the property that
(union of multisets), and
each is -bit accurate.
We generalize Proposition 1 so that it outputs
a collection of cluster oracles:
Proposition 2 (See [4, 5, 17])
Let be an oracle for a univariate polynomial
There is an algorithm
that returns a set
of cluster oracles such that
and each is -bit accurate.
3 Error Analysis of Approximate Specializations
The proofs for this section is found in the Appendix.
Given and ,
our basic goal is to bound the evaluation error
in terms of
This will be done by induction on .
Our analysis aims not just
to produce some error bound, but to express this error
in terms that are easily understood, and which reveals
the underlying inductive structure.
Towards this end, we introduce the following -bound
function: if is a positive integer and ,
i.e., for each .
The support of is where
where . Here,
We assume that .
Our induction variable is .
For , let
E.g., if then .
With this notation, we can write
where each .
E.g., if then and so .
Assume that we are given
and . Also the degree sequences satisfies
, that is the inequality holds
componentwise. Then we may define these quantities