Multivariate Polynomials, Duality and Structured Matrices

Bernard Mourrain^{∗}  Victor Y. Pan ^{†}

December 14, 1998

Abstract

We first review the basic properties of the well known classes of Toeplitz, Hankel, Vandermonde, and other related structured matrices and re-examine their correlations to operations with univariate polynomials. Then we define some natural extensions of such classes of matrices based on their correlations to multivariate polynomials. We describe these correlations in terms of the associated operators of multiplication in the polynomial ring and its dual space, which allows us to generalize these structures to the multivariate case. Multivariate Toeplitz, Hankel, and Vandermonde matrices, Bezoutians, algebraic residues and relations among them are studied. Finally, we show some applications of structured matrices, duality and residues to rootfinding problems for a system of multivariate polynomial equations, where the dual space of linear forms, algebraic residues, Bezoutians and other structured matrices play an important role. The developed techniques enable us to obtain a better insight into the major problems of multivariate polynomial computations and to improve substantially the known techniques of the study of these major problems. From the algorithmic point, we yield acceleration by order of magnitude of the known methods for some fundamental problems of solving multivariate polynomial systems of equations.

Keywords: Structured matrices, polynomial equations, ideals, roots, quotient algebra, dual algebra, residues, basis of idempotents, Jacobians

1991 AMS Subject classification: 12E12, 65F15, 68Q40, 47B35, 14Q20, 68Q25

∗INRIA, SAGA, 2004 route des Lucioles, B.P. 93, 06902 Sophia Antipolis, mourrain@sophia.inria.fr, (partially supported by European ESPRIT project FRISCO, LTR 21.024)

†Department of Mathematics and Computer Science Lehman College, City University of New York, Bronx, NY 10468, USA, VPAN@LCVAX.LEHMAN.CUNY.EDU (Supported by NSF Grants CCR 9625344 and 9732206 and PSC CUNY Awards 668365 and 669363.)

1 Introduction

The main goal of this paper is to summarize and to develop various techniques in the areas of algebraic residues, dual spaces of linear maps and structured matrices and to demonstrate the power of application of these techniques to algorithmic study of polynomial systems of equations; in particular we accelerate the known solution algorithms by order of magnitude. Let us comment on the structure of our presentation and on some specific new results of this paper.

It is well known that the important classes of Toeplitz, Hankel, Vandermonde, and some other structured matrices have a natural characterization in terms of the associated linear operators of scaling and displacements. We will study some extensions of the classes of such matrices, based on their correlations to the fundamental operations with polynomials, such as polynomial multiplication, multipoint evaluation, interpolation, and rootfinding. We will start with a review of the simpler and well known correlations to operations with univariate polynomials and then will use the patterns of this study as basic samples for our extended study where we involve multivariate polynomials. This will enable us to give a natural introduction to some other large topics, intensively studied and having important applications, such as duality, algebraic residues, and Bezoutians.

The correlations between structured matrices and univariate polynomials have been effectively used for the acceleration of structured matrix computation. We will extend these results to the structured matrices associated with multivariate polynomials and will exploit matrix structure to improve substantially the known methods and algorithms for polynomial systems of equations.

In particular, in section 3.2.3, using these tools enables us to simplify substantially the known derivations of the fundamental upper bounds by Bézout and Bernstein on the number D of the roots of a given polynomial system. In section 3.3, we specify our iterative algorithm outlined in the conference paper [24]. The algorithm quadratically converges right from the start to a selected root of a polynomial system of equations that has D distinct and simple roots, and we approximate such a root by using order of D^{2} arithmetic operations (up to a polylogarithmic factor in D). (Hereafter, we will use the abbreviation “ops” for “arithmetic operations”. We say ”ops” rather than ”flops” to cover also rational computations with finite precision.) The algorithm can be applied recursively to compute several roots using order of D^{2}δ ops (up to a polylogarithmic factor) to approximate δ roots for δ ≤ D. In section 3.4, we devise algorithms, also running in D^{2} time (up to a polylog factor), that compute the numbers of distinct roots and distinct real roots of a given polynomial system of equations with real input coefficients. This improves by order of magnitude the known algorithms (not involving structured matrices and algebraic residues), which all require at least order of D^{3} time to solve any of the cited computational problems.

Thus, we reached our main technical goal of developing the basic techniques for the improvement of computations with multivariate polynomials by using the associated structured matrices, dual spaces and algebraic residues. We were able to demonstrate the power of such techniques already in the present paper; in our subsequent works we will show how to accentuate this power further (in particular, by removing the assumption that the residue associated with a given polynomial system is known or readily available) and to elaborate and ameliorate the resulting algorithms from numerical and algebraic points of view. Our progress in these directions has been reported in our recent conference papers [4], [25]. In our present paper we have not touched these aspects and only provided an illustrative example for our approach. Some of the presented techniques appeared earlier in less developed form. In particular, some extensions of the well-known classes of structured matrices associated with univariate polynomials were presented in [33], but they only worked in much more restricted cases, and the restrictions do not allow to apply them to solving polynomial systems.

We will use the following order of presentation. Section 1 deals with structured matrices associated with univariate polynomials. In section 2, we present a natural generalization to the multivariate case, and in section 3, we show some applications to the polynomial root-finding problem in the multivariate case. Various advanced techniques that we use have not been collected together so far, so we detail our presentation and give many comments and some illustrative examples, in particular, one basic example is used thoughout sections 2 and 3.

In this section, we will recall the basic results on the classical matrix structure, presenting them from a univariate polynomial point of view. This will give us a sample pattern, which we will use as a springboard for developing similar techniques in the multivariate case. The reader is referred to appendix A, for the summary of the basic definitions, and to appendix B, for the summary of the estimates for the complexity of some fundamental polynomial and matrix computations.

1.1 Toeplitz operators and matrices

Consider a polynomial t = t_{0}+t_{1} x+⋯+t_{2d} x^{2d} and the map of multiplication by this polynomial t in the ring R = ℂ[x] of polynomials in the variable x with coefficients from the complex field ℂ:

M_{t} : R → R
p ↦ t p.

The matrix M of this map in the monomial basis (obtained by computing the polynomials M_{t}(1), M_{t}(x), M_{t}(x^{2}), …) has the form

1 ; ⋮ ; x^{d} ; ⋮ ; x^{2d} ; ⋮ [ \begin{matrix}t_{0} & & & 0 \\ ⋮ & ... \\ t_{d} & & t_{0} & · \\ ⋮ & ... & ⋮ & · \\ t_{2d} & & t_{d} & · \\ & ... & ⋮ & · \\ 0 & & t_{2d} & · \end{matrix} ] \} T  (1)

The matrix M infinitely continues rightward and downward. Its rows and columns are indexed by the monomials (x^{i}), and its (i, j)-th entry is the coefficient of x^{i} in the polynomial x^{j} t(x). The entries of M are invariant in their shift along the diagonal direction. This property characterizes the class of Toeplitz matrices:

Definition 1.1.1 A matrix T = (t_{i,j}) is a Toeplitz matrix if for all i, j the entry t_{i,j} depends only on i − j, that is, if t_{i,j} = t_{i+1,j+1} for all pairs of (i, j) and (i + 1, j + 1) for which the entries t_{i,j} and t_{i+1,j+1} are defined.

It is immediately observed that any h×k Toeplitz matrix T where max{h, k} ≤ d+ 1 can be obtained as a submatrix of the matrix M defined in (1). Let E = {1, …, x^{d}} and F = {x^{d}, …, x^{2d}} be two linear subspaces of R and let π_{E} (resp. π_{F}) be the projection of R into the vector space generated by E (resp. F). Then the matrix T is just the matrix of the map

T_{t} = π_{F} ∘ M_{t} ∘ π_{E}.

The projections π_{E} and π_{F} select the first columns and the middle rows of M, respectively.

Proposition 1.1.2 A Toeplitz operator (associated with a Toeplitz matrix) is the projection of the multiplication of a fixed polynomial by a polynomial. This is a map from R to R.

Problem 1.1.1 Compute the product of an n × n Toeplitz matrix by a vector.

We observe that such a product is a subvector of the coefficient vector of the product of two polynomials of R. Therefore, by theorem B.1.1 of appendix B, we may solve problem 1.1.1 in O(n log(n)) ops.

Hereafter, we will use the abbreviation f.p.s. for formal power series and will write W^{t} to denote the transpose of an operator, matrix or vector W. Similarly to the map M_{t}, we define the map

M_{t}^{t} : S → S
q(∂) ↦ t(x) ⋆ q(∂) = π_{+}(t(∂^{−1})q(∂)),

where S is the ring of f.p.s. in the variable ∂, ∂^{i} is the differential form: p ↦ \frac{1}{i!}p^{(i)}(0), and π_{+} is the projection of an f.p.s. in ∂ and ∂^{−1} into an f.p.s. in S obtained by deleting all the monomials in ∂^{−1}, that is, π_{+} is the projection on the monomials of positive degree in ∂. The matrix of this map is the transpose of the matrix of M_{t}, where we can extract the transpose of the matrix T:

[ \begin{matrix}t_{0} & ⋯ & t_{d} & ⋯ & t_{2d} & 0 \\ & ... & ⋮ & ... & ⋮ & ... \\ & & t_{0} & ⋯ & t_{d} \\ & & & ... & ⋮ & ... \\ 0 & & & & t_{0} \end{matrix} ] .

1.2 Hankel operators and matrices

Next, consider the multiplication map defined by a f.p.s. h(∂) = h_{0} + h_{1}∂ + ⋯ + h_{2d}∂^{2d} + ⋯ in ∂ as follows: for any polynomial p ∈ ℂ[x] we compute the product, p(∂^{−1})h(∂), and project it into the space of monomials of positive degree. (Then again, the reader may think of ∂ as a variable and of ∂^{−1} as its reciprocal, and we interpret ∂^{i} as the linear map p ↦ \frac{1}{i!}p^{(i)}(0).) Here is the matrix M representing such maps:

1 ; ⋮ ; ∂^{d} ; ⋮ ; ∂^{2d} ; · [ \begin{matrix}h_{0} & & h_{d} & · \\ ⋮ & … & ⋮ & · \\ h_{d} & & h_{2d} & · \\ ⋮ & … & · & · \\ h_{2d} & & · & · \\ · & · & · & · \end{matrix} ] \} H  (2)

The matrix M infinitely continues rightward and downward in this case. Its columns are indexed by monomials in x and its rows by monomials in ∂. The (i, j)-th entry of this matrix is the coefficient of ∂^{i} in ∂^{−j}h(∂), which explains why its entries are invariant in their shifts into the antidiagonal direction. This property characterizes the class of Hankel matrices.

Definition 1.2.1 A matrix H = (h_{i,j}) is a Hankel matrix if its entry h_{i,j} depends only on i + j, that is, if h_{i+1,j−1} = h_{i,j} for all pairs (i, j) for which the entries are defined.

Definition 1.2.2 The space of linear forms (maps) from R to ℂ, that is, the dual space of the ring of polynomials R, is denoted by R̂. The ring of such maps from ℂ[x] to ℂ[∂] is generalized to the ring of maps from the ring of f.p.s. in x to the ring of f.p.s. in ∂. We denote the latter ring by both S and ℂ[[∂]].

The matrix M is the matrix of the map

H_{h} : ℂ[x] → ℂ[[∂]]  (3)
p(x) ↦ p(x) ⋆ h(∂) = π_{+}(p(∂^{−1})h(∂)).

where π_{+} is the projection into the space of monomials of positive degree in ∂.

We immediately observe that any h×k Hankel matrix H where max{h, k} ≤ d+1 is a submatrix of the above matrix M, defined in (2), (for some h ∈ ℂ[∂]). Let E = {1, x, …, x^{d}}, F = {1, ∂, …, ∂^{d}} be the two monomial sets in x and ∂, respectively, and let π_{E} and π_{F} be the corresponding projections on the vector spaces generated by these sets. Then the matrix H is the matrix of the following map:

H_{h} = π_{F} ∘ H_{h} ∘ π_{E}.

The projections π_{E} and π_{F} select the first columns and the middle rows of the matrix M, respectively.

Proposition 1.2.3 A Hankel operator (associated with a Hankel matrix) can be defined as the projection of the multiplication of a (projected) polynomial by a fixed Laurent polynomial.

Problem 1.2.1 Compute the product of a (d+1)×(d+1) Hankel matrix by a vector.

We observe that such a product is a subvector of the coefficient vector of the product of a fixed polynomial h(∂) by a polynomial in ∂^{−1}. Therefore, by theorem B.1.1 of appendix B, we may solve problem 1.2.1 in O(d log(d)) ops.

1.3 Bezoutians

Next, let us study linear maps from ℂ[[∂]] to ℂ[x]. First, consider a polynomial in two variables x and y:

Θ(x, y) = ∑_{i=0,j=0}^{d−1} θ_{i,j} x^{i} y^{j}.

To any element Λ(∂) ∈ ℂ[[∂]], we associate the constant coefficient in ∂ (that is, the ∂-free term) of the product

Θ(x, ∂^{−1}) Λ(∂).

This defines a map B from ℂ[[∂]] to ℂ[x]. We immediately verify that the matrix of this map (which can be obtained by computing the constant coefficients in ∂ of Θ(x, ∂^{−1}) ∂^{j} : B(1) = ∑_{i=0}^{d−1} θ_{i,0}x^{i}, B(∂) = ∑_{i=0}^{d−1} θ_{i,1}x^{i},…) is precisely the coefficient matrix [θ_{i,j}]_{0≤i,j≤d−1} of Θ(x, y).

A fundamental example of such a polynomial is the Bezoutian defined as follows:

Definition 1.3.1 Let p and q be two polynomials of ℂ[x]. The term Bezoutian of p and q is used for both the bivariate polynomial

Θ_{q,p}(x, y) = \frac{p(x) q(y) − p(y) q(x)}{x − y} = ∑_{0≤i,j≤d−1} θ_{i,j}^{q,p}x^{i}y^{j}

and the matrix

B_{q,p} = [ \begin{matrix}θ_{0,0}^{q,p} & ⋯ & θ_{0,d−1}^{q,p} \\ ⋮ & & ⋮ \\ θ_{d−1,0}^{q,p} & ⋯ & θ_{d−1,d−1}^{q,p} \end{matrix} ] .

B_{q,p} : ℂ[[∂]] → ℂ[x] denotes the associated map, B_{q,p}(Λ) ↦ ∂-free term of Θ_{q,p}(x, ∂^{−1})Λ(∂). The image of this map can be expressed as the product

[1, x, …, x^{d−1}] B_{q,p}[λ_{0}, …, λ_{d−1}]^{t},

where Λ(∂) = ∑_{i=0}^{∞}λ_{i}∂^{i}.

In particular, if p = p_{0} + p_{1} x + ⋯ + p_{d} x^{d}, then the polynomial Θ_{1,p} is of the form

Θ_{1,p}(x, y) = ∑_{i=0}^{d−1} x^{i} Θ_{i}^{p}(y),

where Θ_{i}^{p}(y) = p_{i+1} + p_{i+2} y + ⋯ + p_{d} y^{d−i−1}. This polynomial is also called the i-th Horner polynomial, for it corresponds to the i-th polynomial, appearing in the so-called Horner rule for polynomial evaluation. It can be also written as

Θ_{i}^{p}(y) = π_{+}(y^{−i−1}p(y)),  (4)

where π_{+} is the projection on the set of polynomials in y. We immediately observe that the matrix B_{1,p} associated with Θ_{1,p} is a triangular Hankel matrix of the form

[ \begin{matrix}p_{1} & ⋯ & p_{d} \\ ⋮ & … \\ p_{d} & & 0 \end{matrix} ] .  (5)

More generally, we have the decomposition

Θ_{q,p}(x, y) = \frac{p(x) q(y) − p(y) q(x)}{x − y}

= \frac{p(x) − p(y)}{x − y} q(y) − \frac{q(x) − q(y)}{x − y} p(y) = Θ_{1,p}(x, y) q(y) − Θ_{1,q}(x, y) p(y).

This implies

B_{q,p}(Λ) = B_{1,p}(q ⋆ Λ) − B_{1,q}(p ⋆ Λ)

for any Λ(∂) ∈ ℂ[[∂]] or, in terms of operators,

B_{q,p} = B_{1,p} ∘ M_{q}^{t} − B_{1,q} ∘ M_{p}^{t}.  (6)

In term of matrices, this yields the Barnett formula,

B_{q,p} = [ \begin{matrix}p_{1} & ⋯ & p_{d} \\ ⋮ & … \\ p_{d} & & 0 \end{matrix} ] [ \begin{matrix}q_{0} & ⋯ & q_{d−1} \\ & ... & ⋮ \\ 0 & & q_{0} \end{matrix} ] − [ \begin{matrix}q_{1} & ⋯ & q_{d} \\ ⋮ & … \\ q_{d} & & 0 \end{matrix} ] [ \begin{matrix}p_{0} & ⋯ & p_{d−1} \\ & ... & ⋮ \\ 0 & & p_{0} \end{matrix} ] .

1.4 Vandermonde operators and matrices

Consider the linear space R_{d} of polynomials of degree at most d and d + 1 distinct points in ℂ : Ξ = {ξ_{0}, …, ξ_{d}}. Also consider the next two bases of R_{d}:

• the basis of monomials ⟨1, x, …, x^{d}⟩

• and the basis of Lagrange interpolation polynomials

⟨L_{i} = L_{i}(x) = ∏_{j≠i} \frac{x − ξ_{l}}{ξ_{i} − ξ_{j}} , i = 0, …, d⟩.

Any polynomial p ∈ R_{d} can be decomposed in the latter basis as follows:

p(x) = ∑_{i=0}^{d} p(ξ_{i}) L_{i}(x).  (7)

We deduce from this decomposition that the (d + 1) × (d + 1) matrix of the basis transformation from (x^{i})_{i=0,...,d} to (L_{i}(x))_{i=0,1,...,d} is the Vandermonde matrix,

V (Ξ) = [ \begin{matrix}1 & ξ_{0} & ⋯ & ξ_{0}^{d} \\ 1 & ξ_{1} & ⋯ & ξ_{1}^{d} \\ ⋮ & & & ⋮ \\ 1 & ξ_{d} & ⋯ & ξ_{d}^{d} \end{matrix} ] .

Remark 1 Many authors use the name “Vandermonde matrix” for V^{t}(Ξ), the transpose of V (Ξ).

Problem 1.4.1 Multiply the matrix V (Ξ) by a vector p = (p_{0}, …, p_{d})^{t} or, equivalently, evaluate a polynomial p(x) = ∑_{i=0}^{d} p_{i}x^{i} on the set of points Ξ = {ξ_{0}, …, ξ_{d}}.

Clearly, the multiplication of the row vector (1, ξ_{i}, …, ξ_{i}^{d})^{t} by the vector p = (p_{0}, …, p_{d}) amounts to evaluating the polynomial p(x) = p_{0} + ⋯ + p_{d} x^{d} at the point ξ_{i}. Equivalently, the coefficients p(ξ_{i}) of p = p(x) in the Lagrange basis can be obtained by means of the evaluation of p = p(x) at the points ξ_{i}.

Problem 1.4.2 Solve linear system V (Ξ)v = w or, equivalently, interpolate to a polynomial p(x) from its values w_{0}, …, w_{d} on a set Ξ = {ξ_{0}, …, ξ_{d}}.

The known algorithms solve problems 1.4.1 and 1.4.2 in O(d log^{2}d) ops (see [3], pp. 25-26).

Evaluation at a point is an example of a linear form (map), and equation (7) shows that the dual basis of (L_{i})_{i=0,...,d} (that is, the linear forms (maps) that compute the coefficients of a polynomial p in this basis) is the set of linear forms (1_{ξ_{i}})_{i=0,...,d} of the evaluation at ξ_{i}: 1_{ξ_{i}}(p) = p(ξ_{i}). Such an evaluation will play important role in the following, so we will next define it formally:

Definition 1.4.1 For any point ξ ∈ ℂ, let 1_{ξ} ∈ R̂ ⊂ R̂_{d−1} denote the linear form that corresponds to the evaluation at ξ:

1_{ξ} : R → ℂ
p ↦ p(ξ).

Note that R̂ is the subset of the dual space R̂_{d} made by the linear forms on the vector space of polynomials of degree at most d and that the coordinates of the evaluation 1_{ξ} ∈ R̂_{d} in the dual basis ⟨1, ∂, …, ∂^{d}⟩ of R̂_{d} are obtained by computing 1_{ξ}(x^{i})_{i=0,...,d}. This yields the vector (1, ξ, ξ^{2}, …, ξ^{d}). In terms of polynomials in ∂, we have

1_{ξ} = 1 + ξ ∂ + ⋯ + (ξ ∂)^{d} = \frac{1 − (ξ ∂)^{d+1}}{1 − ξ y} .

Thus, the matrix of the basis transformation from the basis (1_{ξ_{i}})_{i=0,...,d} to the dual basis ⟨1, ∂, …, ∂^{d}⟩ of ⟨1, x, …, x^{d}⟩ is given by

V^{t}(Ξ) = [ \begin{matrix}1 & 1 & ⋯ & 1 \\ ξ_{0} & ξ_{1} & ⋯ & ξ_{d} \\ ⋮ & & & ⋮ \\ ξ_{0}^{d} & ξ_{1}^{d} & ⋯ & ξ_{d}^{d} \end{matrix} ] .

Problem 1.4.3 Multiply V^{t}(Ξ) by a vector.

Problem 1.4.4 Solve the linear system V^{t}(Ξ)v = w.

Problems 1.4.3 and 1.4.4 can be solved in O(d log^{2} d) ops (see [3], pp.141-144). Problem 1.4.3 can be also solved at this cost by reduction to problem 1.4.1 (see theorem B.2.1 of appendix B.2). A slower but technically interesting approach relies on the observation that the multiplication of the matrix V^{t}(Ξ). by a vector Λ = [λ_{0}, …, λ_{d}] amounts to the computation, in the monomial basis, of the polynomial

∑_{i=0}^{d} λ_{i} \frac{1 − (ξ_{i} ∂)^{d+1}}{1 − ξ_{i} ∂} .

If the interpolation points are the d-th roots of unity, we arrive at a special Vandermonde matrix, sometimes called the Fourier matrix. In this special case, problems 1.4.1-1.4.4 represent forward and inverse discrete Fourier transforms (DFTs) and can be solved by using O(d log d) ops. The inverse of the Fourier matrix is the transpose of its conjugate (up to the factor d).

1.5 Relations between Bezoutians and Hankel matrices

The Hankel operators correspond to some maps from ℂ[x] to ℂ[[∂]], whereas the Bezoutians define some maps from ℂ[[∂]] to ℂ[x]. It is natural to ask if there is a relationship between the maps of these two classes. This is what we are going to examine next. We will use the basic concept of the ideal I = (p), generated by p ∈ R, that is, the set of polynomials {pq, q ∈ R}.

In order to relate these two classes of operators to each other, we will next describe the elements h(∂) ∈ ℂ[[∂]] such that H_{h} of (3) vanishes on all multiples of a fixed polynomial p(x) = p_{0} + p_{1} x + ⋯ + p_{d}x^{d} of degree exactly d (that is, on the ideal generated by p): H_{h}(p v) = 0 for all elements v ∈ R.

Proposition 1.5.1 The class of f.p.s. h ∈ ℂ[[∂]] such that H_{h} vanishes on the ideal (p) generated by a polynomial p = p_{0} +p_{1} x+⋯+p_{d} x^{d} of degree d (p_{d} ≠ 0) coincides with the class of rational functions

h(∂) = \frac{∂^{−1}r(∂)}{p(∂^{−1})} = h_{0} + h_{1} ∂ + ⋯ + h_{d−1} ∂^{d−1} + ⋯,  (8)

where r(x) = ∑_{i=0}^{d−1} r_{i}x^{i} is a polynomial in R_{d−1}.

Proof. First, note that the rational fraction h(∂) = \frac{r_{0}∂^{d−1}+r_{1}∂^{d−2}+···+r_{d−1}}{p_{d}+p_{d−1} ∂+···+p_{0} ∂^{d}} is an f.p.s. in ∂, having no terms ∂^{−i} for i > 0, since p_{d} ≠ 0.

To show that H_{h} vanishes on the ideal (p) for h(∂) of (8), observe that

h(∂)p(∂^{−1})v(∂^{−1}) = ∂^{−1} r(∂^{−1})v(∂^{−1}),

for v ∈ R, has only terms with negative powers of ∂ since r(x) and v(x) are polynomials. Therefore, p(x)v(x) ⋆ h(∂) = 0 for any polynomial v(x) ∈ R.

Finally, let H_{h} vanish on (p), for an f.p.s. h = h(∂). This means that

π_{+}(p(∂^{−1}) h(∂^{−1})) = 0,

that is, p(∂^{−1})h(∂) is an f.p.s. in ∂^{−1}, with no constant term: p(∂^{−1})h(∂) = ∂^{−1} r(∂^{−1}), where r(∂) is an f.p.s. from ℂ[[∂]]. Furthermore, by replacing ∂^{−1} by x, we obtain that r(x) = x^{−1}p(x)h(x^{−1}) = π_{+}(x^{−1}p(x)h(x^{−1})), so that r(x) is a polynomial of degree less than deg(p(x)) = d. □

The proposition implies that the class of the f.p.s. h ∈ ℂ[[∂]] such that H_{h} vanishes on (p) is the class of all multiples of the f.p.s. τ = τ(∂) = \frac{∂^{−1}}{p(∂^{−1})} \frac{∂^{d−1}}{p_{d}+p_{d−1} ∂+···+p_{0} ∂^{d}}, called the (algebraic) residue of p. (This concept extends the concept of the residue of an analytic function.) We will next give a characterization of this residue that can be easily generalized to the multivariate case.

Proposition 1.5.2 Let p = p_{0} + p_{1}x + ⋯ + p_{d}x^{d} be a fixed polynomial of degree exactly d. Then the residue τ = τ_{p}(∂) is the unique element of ℂ[[∂]] that satisfies:

1. τ vanishes on the multiples of p,

2. B_{1,p}(τ) = 1,

where B_{1,p} is the map defined in definition 1.3.1.

Proof. Property 1. of τ follows from the definition of τ and proposition 1.5.1. Now, by the definition of τ = τ_{p}(∂), the element τ_{p}(∂) = ∑_{i=0}^{∞} τ_{i}∂^{i} = ∑_{i=0}^{∞} τ(x^{i})∂^{i} of ℂ[[∂]] has the form

\frac{1}{p_{d}} ∂^{d−1} + τ_{d}∂^{d} + ⋯,

that is, τ_{0} = ⋯ = τ_{d−2} = 0, τ_{d−1} = \frac{1}{p_{d}}, which means that the linear form (map) associated with τ vanishes on 1, x, …, x^{d−2} and equals \frac{1}{p_{d}} on x^{d−1}.

Now we obtain from definition 1.3.1 that

B_{1,p}(τ) = [1, x, ..., x^{d−1}]B_{1,p}[0, ...0, 1/p_{d}]^{t}.

As B_{1,p} is of the form (5), we immediately check that

B_{1,p} [0, …, 0, \frac{1}{p_{d}}]^{t} = [1, 0, …, 0]^{t},

which implies property 2. of τ, that is B_{1,p}(τ) = 1. On the other hand, since B_{1,p} is nonsingular, the latter vector equation implies that τ_{0} = ⋯ = τ_{d−2} = 0, τ_{d−1} = \frac{1}{p_{d}}. Furthermore, as any polynomial q can be decomposed uniquely as q = p u + r, where r ∈ R_{d−1} = ⟨1, x, …, x^{d−1}⟩, the linear form τ is entirely determined by the values τ_{0}, …, τ_{d−1} and the fact that it vanishes on the multiples of p. This proves the uniqueness of the element of ℂ[[∂]] satisfying properties 1. and 2. and thus completes the proof of the proposition. □

Proposition 1.5.3 The set (Θ_{i}^{p})_{i=0,...,d−1} is the dual basis of the monomial basis (x^{i})_{i=0,...,d−1} for the inner product associated to τ:

τ(x^{i} Θ_{j}^{p}(x)) = \{ \begin{matrix}1 if i = j \\ 0 otherwise. \end{matrix}  (9)

Proof. For 0 ≤ i, j ≤ d − 1, we have (see (4))

τ(x^{i}Θ_{j}^{p}(x)) = τ(x^{i}π_{+}(x^{−j−1}p(x))) = τ(x^{i−j−1}p(x)).

The last equation holds because x^{i} (x^{−j−1}p(x) − π_{+}(x^{−j−1}p(x))) is in the vector space R_{−d,d−2} and τ vanishes on this vector space. If i > j, then x^{i−j−1}p(x) is in the ideal (p) generated by p in R, and τ vanishes on this ideal. On the other hand, if i < j, then x^{i−j−1}p(x) is in the vector space R_{−d,d−2}, and τ vanishes on this vector space too. For i = j, we obtain τ(x^{−1}p(x)) = τ(p_{d} x^{d−1}) = 1, which proves relations (9). □

We immediately deduce from this result the following corollary.

Corollary 1.5.4 Let B_{1} = B_{1,p} and let H_{1} = H_{τ} be the Hankel matrix of the map H_{τ} of (3) for h = τ. Then

B_{1} H_{1} = H_{1} B_{1} = 𝕀_{d},

where 𝕀_{d} is the d × d identity matrix.

Proof. From (9), we deduce that

∑_{j=0}^{d} x^{j} τ(x^{i} Θ_{j}^{p}(x)) = x^{i}.

On the other hand, the left-hand side of this equation equals B_{1,p}(x^{i} ⋆ τ). Thus, if we compose the two maps H_{τ} : R_{d−1} → ℂ[[∂]] and B_{1,p} : ℂ[[∂]] → R_{d−1}, we obtain that

B_{1,p} ∘ H_{τ}(x^{i}) = B_{1,p}(x^{i} ⋆ τ) = x^{i},

for i = 0, …, d − 1. In other words,

B_{1,p} ∘ H_{τ} = 𝕀_{R_{d−1}}

or, equivalently, B_{1} H_{1} = 𝕀_{d}, which shows that the inverse of the Bezoutian B_{1} is the Hankel matrix H_{1} and vice versa. □

In view of corollary 1.5.4, the Barnett formula of section 1.3 (that is, the matrix version of (6)) gives us an expression for the inverse of a Hankel matrix via the 2d+2 scalars p_{0}, q_{0}, …, p_{d}, q_{d}. Such an expression extends the Gohberg-Semencul celebrated formula, known for the inverses of Toeplitz matrices (cf. [3]).

2 Generalization to the multivariate case

Our next goal, is the extension of the approach and the results of the previous section to the multivariate case. We will start with recalling some definitions and techniques used in [1], [8], [11], [21]-[24], [32]. Then, in sections 2.8, 2.10–2.12, we will develop some new techniques to be used in section 3.

2.1 Polynomial ring

The definitions of the previous section and appendix A can be immediately extended to the n-variate case, for any natural n. In this case, R = ℂ[x] is replaced by the ring ℂ[x_{1}, …, x_{n}] of multivariate polynomials in x_{1}, …, x_{n}; x and ∂ are assumed to be vectors, rather than scalars, x = (x_{1}, …, x_{n}) and ∂ = (∂_{1}, …, ∂_{n}). We keep denoting by R_{d} the subspace of all polynomials of degree at most d. Instead of working in the complex space ℂ, we could have allowed the vector spaces over any algebraically closed field 𝕂, and then R would denote the space of multivariate polynomials in x, with coefficients from 𝕂. Our results of this section would be easily extended, but, to simplify our presentation, we will state them for 𝕂 = ℂ. We will let L = ℂ[x_{1}^{±1}, …, x_{n}^{±1}] denote the ring of Laurent’s polynomials in the variables x_{1}, …, x_{n}. For any element p of R, let

M_{p} : R → R  (10)
r ↦ p r

denote the operator of multiplication by p in R.

Hereafter, I = (p_{1}, …, p_{n}) denotes the ideal of R = ℂ[x] generated by the elements p_{1}, …, p_{n}, that is, the set of polynomial combinations ∑_{i} p_{i}q_{i} of these elements. A = R/I denotes the quotient ring defined in R by I, and ≡ denotes the equality in A. We assume that the set of the common zeros of the n polynomials p_{1}, …, p_{n} (that is, the set of the roots of the polynomial system p_{1} = ⋯ = p_{n} = 0) is finite and denote it Z = Z(I) = {ζ_{1}, …, ζ_{d}}. This implies that the vector space A has a finite dimension D, D ≥ d. (D is the number of the roots counted with their multiplicities.)

2.2 The quotient algebra

Our main objective is the analysis of the structure of A, in particular in order to devise efficient algorithms for computing the zeros in Z(I).

The first operator that comes naturally in this study is the operator of multiplication by an element of A, based on (10). For any element a ∈ A, we define the map

\begin{matrix}\overline{M}_{a} : A & → & A \\ b & ↦ & a b. \end{matrix}

An important property of this operator is given in the next theorem (see [1], [32], [22]):

Theorem 2.2.1 The set of the eigenvalues of the linear operator \overline{M}_{a} is exactly {a(ζ_{1}), …, a(ζ_{d})}.

Proof. Let p(x) = ∏_{ζ∈Z(I)}(a(x) − a(ζ)). This polynomial vanishes on Z(I), so that (according to the Nullstellensatz, see [9]) there exists d = d_{p} ∈ ℕ such that p(x)^{d} ∈ I.

Consequently, we have

∏_{ζ∈Z(I)}(\overline{M}_{a} − a(ζ)𝕀)^{d} = 0,

where 𝕀 is the identity map : b → b, and the minimal polynomial of \overline{M}_{a} divides ∏_{ζ∈Z(I)}(T − a(ζ))^{d}, for indeterminate T. This implies that an eigenvalue of \overline{M}_{a} is necessarily in the set {a(ζ_{1}), …, a(ζ_{d})}. On the other hand, we will show in theorem 2.4.1, by using duality, that for any ζ ∈ Z(I), a(ζ) is an eigenvalue of the transpose of \overline{M}_{a}. □

Example Let n = 2,

p_{1} = x_{1}^{2} + 2 x_{2}x_{1} − x_{1} − 1, p_{2} = x_{1}^{2} + x_{2}^{2} − 8 x_{1}.

We check (by hand computation) that a basis of A = ℂ[x_{1}, x_{2}]/(p_{1}, p_{2}) is (1, x_{1}, x_{2}, x_{1} x_{2}) and that the matrix of multiplication by x_{1} in this basis is

M_{x_{1}} = [ \begin{matrix}0 & 1 & 0 & −\frac{14}{5} \\ 1 & 1 & 0 & −\frac{12}{5} \\ 0 & 0 & 0 & \frac{1}{5} \\ 0 & −2 & 1 & \frac{29}{5} \end{matrix} ] .

The eigenvalues of M_{x_{1}} are the first coordinates of the roots, that is,

6.8200982, −0.19395427 + 0.20520688 i, −0.19395427 − 0.20520688 i, 0.36781361.

Theorem 2.2.1 reduces the nonlinear problem of solving a polynomial system of equations to a well known problem of linear algebra. The reduction, however, involves the analysis of the structure of A and the properties of the operators of multiplication, and this leads us to the study of the dual space, the multivariate Bezoutians, and structured matrices associated with multivariate polynomials. This is needed, in particular, in order to express explicitly the matrices of multiplication associated with the operator \overline{M}_{a}. (Such matrices are called multiplication tables.) The main difficulties stem from the requirement to work modulo the ideal I, and the dual space, Bezoutians, and structured matrices are effective tools for handling this nontrivial problem.

Definition 2.2.2 Hereafter, ℕ denotes the set of nonnegative integers, and we fix a subset E ⊂ ℕ^{n}, such that (x^{α})_{α∈E} is a basis of A.

2.3 Dual space

Let R̂ denote the dual of the ℂ-vector space R, that is, the space of linear forms

λ : R → ℂ
p ↦ λ(p).

(R will be the primal space for R̂.) The evaluation at a fixed point ζ is a well-known example of such a linear form:

1_{ζ} : R → ℂ
p ↦ p(ζ).

Another class of linear forms is obtained by using differential operators. Namely, for any a = (a_{1}, …, a_{n}) ∈ ℕ^{n}, consider the map

∂^{a} : R → ℂ
p ↦ \frac{1}{∏_{i=1}^{n} a_{i}!} (d_{x_{1}})^{a_{1}} ⋯ (d_{x_{n}})^{a_{n}} (p)(0),  (11)

where d_{x_{i}} is the derivative with respect to the variable x_{i}. We denote this linear form by ∂^{a} = (∂_{1})^{a_{1}} ⋯ (∂_{n})^{a_{n}} and for any (a_{1}, …, a_{n}) ∈ ℕ^{n}, (b_{1}, …, b_{n}) ∈ ℕ^{n} observe that

\frac{1}{∏_{i=1}^{n} a_{i}!}∂^{a} (∏_{i=1}^{n} x_{i}^{b_{i}}) (0) = \{ \begin{matrix}1 if ∀i, a_{i} = b_{i}, \\ 0 otherwise. \end{matrix}

It immediately follows that (∂^{a})_{a∈ℕ^{n}} is the dual basis of the primal monomial basis. By applying Taylor’s expansion formula at 0, we decompose any linear form Λ ∈ R̂ as

Λ = ∑_{a∈ℕ^{n}} Λ(x^{a}) ∂^{a}.

The map Λ → ∑_{a∈ℕ^{n}} Λ(x^{a}) ∂^{a} defines a one-to-one correspondence between the set of linear forms Λ and the set ℂ[[∂_{1}, … ∂_{n}]] = ℂ[[∂]] = {∑_{a∈ℕ^{n}} λ_{a}∂_{1}^{a_{1}} ⋯ ∂_{n}^{a_{n}}} of f.p.s. (formal power series) in the variables ∂_{1}, …, ∂_{n}.

As in the univariate case, we will identify R̂ with ℂ[[∂_{1}, …, ∂_{n}]]. The evaluation at 0 corresponds to the constant 1, under this definition. It will also be denoted 1_{0} = ∂^{0}.

Example

(1 + ∂_{1}^{2} ∂_{2})(1 + 2 x_{1} x_{2} + 10 x_{1}^{2} x_{2}) = 11.

We will continue using this example throughout sections 2 and 3.

Let us next examine the structure of the dual space. We can multiply a linear form by a polynomial (we say that R̂ is an R-module) as follows. For any p ∈ R and λ ∈ R̂, we define p ⋆ Λ as

p ⋆ Λ : R → ℂ
q ↦ Λ(p q).

What kind of operation does this multiplication induces on the formal power series representation? For any pair of elements p ∈ R and d ∈ ℕ, d > 1, we have

(d_{x_{i}})^{d} (x_{i} p)(0) = (d_{x_{i}})^{d−1} (p + x_{i} d_{x_{i}}p) (0)
= (d_{x_{i}})^{d−2} (2 d_{x_{i}}(p) + x_{i} (d_{x_{i}})^{2} (p)) (0)
= d (d_{x_{i}})^{d−1} (p)(0) + x_{i} (d_{x_{i}})^{d} (p)(0)
= d (d_{x_{i}})^{d−1} p(0).

Also we surely have d_{x_{i}}(x_{i} p)(0) = d p(0). Consequently, for any pair of elements p ∈ R, d = (d_{1}, …, d_{n}) ∈ ℕ^{n}, where d_{i} ≠ 0 for a fixed i, we obtain that

x_{i} ⋆ ∂^{d}(p) = ∂^{d}(x_{i} p)
= ∂_{1}^{d_{1}} ⋯ ∂_{i−1}^{d_{i−1}}∂_{i}^{d_{i}−1}∂_{i+1}^{d_{i+1}} ⋯ ∂_{n}^{d_{n}}(p),

that is, x_{i} acts as the inverse of ∂_{i} in ℂ[[∂]]. This is the reason why in the literature such a representation is referred to as the inverse systems (see, for instance, [21]). If d_{i} = 0, then x_{i} ⋆ ∂^{d}(p) = 0, which allows us to redefine the product p ⋆ Λ as follows:

Proposition 2.3.1 For any p, q ∈ R and any Λ(∂) ∈ ℂ[[∂]], we have

p ⋆ Λ(q) = Λ(p q) = π_{+}(p(∂^{−1}) Λ(∂))(q).

Example (continued) (x_{1} ⋆ (1 + ∂_{1}^{2} ∂_{2})) (1 + 2 x_{1} x_{2} + 10 x_{1}^{2} x_{2})

= (1 + ∂_{1}^{2} ∂_{2})(x_{1} + 2 x_{1}^{2}x_{2} + 10x_{1}^{3} x_{2}) (12) = ∂_{1} ∂_{2}(1 + 2 x_{1} x_{2} + 10 x_{1}^{2} x_{2}) = 2. (13)

For any linear form Λ ∈ R̂, let

H_{Λ} : R → R̂
r ↦ r ⋆ Λ

denote the operator of multiplication by Λ, from R to R̂.

2.4 The dual of the quotient algebra

Now, let Â be the dual space of A. It is possible to identify the set Â with the elements of R̂ that vanish on I. Thus, the set Â will be also denoted I^{⊥}. Now, for any element a ∈ A, we can describe the transposed operator \overline{M}_{a}^{t}:

\begin{matrix}\overline{M}_{a}^{t} : Â & → & Â \\ Λ & ↦ & a ⋆ Λ = Λ ∘ \overline{M}_{a}. \end{matrix}

The matrix associated to this operator is the transpose of the matrix associated to the matrix \overline{M}_{a}.

We have already described the eigenvalues of this operator in theorem 2.2.1 and will give now a description of its eigenvectors (see [22], [32]):

Theorem 2.4.1 The common eigenvectors of the operators \overline{M}_{a}^{t}, for a ∈ A, are (up to a scalar factor) the evaluations 1_{ζ_{1}}, …, 1_{ζ_{d}}, where 1_{ζ} : p → p(ζ).

Proof. For any pair of polynomials a, b ∈ R and any ζ_{i} ∈ Z(I), we have

\overline{M}_{a}^{t}(1_{ζ_{i}})(b) = 1_{ζ_{i}}(a b) = a(ζ_{i}) 1_{ζ_{i}}(b),

that is, \overline{M}_{a}^{t}(1_{ζ_{i}}) = a(ζ_{i}) 1_{ζ_{i}}. Moreover, 1_{ζ_{i}} is in Â, because ζ_{i} is a common root of the polynomials in I. Therefore, for any a ∈ R, 1_{ζ_{i}} is an eigenvector of \overline{M}_{a}^{t} associated with the eigenvalue a(ζ_{i}).

Conversely, let us prove that the common eigenvectors of (\overline{M}_{x_{i}}^{t})_{i=1,...,n} are (up to a scalar factor) exactly 1_{ζ_{1}}, …, 1_{ζ_{d}}. Let Λ ∈ Â be a non-zero common eigenvector of (\overline{M}_{x_{i}}^{t})_{i=1,...,n} for the eigenvalues (γ_{i})_{i=1,...,n}: x_{i} ⋆Λ−γ_{i} Λ = 0. Then, for any monomial x^{α} of R, we have

x_{i} ⋆ Λ(x^{α}) = Λ(x_{i} x^{α}) = γ_{i}Λ(x^{α}).

By induction, this implies that Λ(x^{α}) = γ^{α} Λ(1). In other words, Λ = Λ(1) 1_{γ}, where γ = (γ_{1}, …, γ_{n}) ∈ ℂ^{n} and 1_{γ} ∈ R̂ is the evaluation at γ. As Λ ∈ Â ≡ I^{⊤}, we have Λ(p) = Λ(1) 1_{γ}(p) = Λ(1) p(γ) = 0, for any p ∈ I, which implies that γ ∈ Z(I). □

Both theorems 2.2.1 and 2.4.1 reduce the solution of a polynomial system to matrix eigenproblem, but theorem 2.4.1 has an advantage compared to theorem 2.2.1: Each eigenvector of an operator \overline{M}_{a}^{t} defines all the coordinates of a root (whereas each eigenvalue of M_{a} defines only one coordinate or the inner product of the vector of a root by a fixed vector defined by a ∈ A). Indeed, the evaluations 1_{ζ_{i}} at the roots ζ_{i} ∈ Z(I) are eigenvectors of \overline{M}_{a}^{t}. From these evaluations 1_{ζ_{i}}, we can recover the coordinates ζ_{i,j} = 1_{ζ_{i}}(x_{j}) of the root 1_{ζ_{i}}. We will specify this computation in section 3.1.

2.5 Quasi-Toeplitz and quasi-Hankel matrices

Definition 2.5.1 Let E and F be two finite subsets of ℕ^{n} and let M = (m_{α,β})_{α∈E,β∈F} be a matrix whose rows are indexed by the elements of E and columns by the elements of F. To the end of this definition, let i denote the i-th coordinate vector of ℕ^{n}, i = 1, …, n.

• M is an (E, F) quasi-Toeplitz matrix iff, for all α ∈ E, β ∈ F, the entries m_{α,β} = t_{α−β} depend only on α − β, that is, if for every i = 1, …, n, we have m_{α+i,β+i} = m_{α,β}, provided that α, α + i ∈ E; β, β + i ∈ F; such a matrix M is associated with the polynomial T_{M}(x) = ∑_{u∈E+F} t_{u} x^{u}.

• M is an (E, F) quasi-Hankel matrix iff, for all α ∈ E, β ∈ F, the entries m_{α,β} = h_{α+β} depend only on α + β, that is, if for every i = 1, …, n, we have m_{α−i,β+i} = m_{α,β} provided that α, α − i ∈ E; β, β + i ∈ F; such a matrix M is associated with the Laurent polynomial H_{M}(∂) = ∑_{u∈E−F} h_{u}∂^{u}.

If we work with the Laurent polynomials, these definitions can be immediately extended to subsets E, F of ℤ^{n}, ℤ denoting the set of all integers.

For E = [0, …, h − 1] and F = [0, …, k − 1], definition 2.5.1 turns into the usual definition of h × k Hankel (resp. Toeplitz) matrices (see sections 1.1 and 1.2). For E and F forming rectangles in ℕ^{n}, the quasi-Toeplitz matrices appeared in [33] under the name of multilevel Toeplitz matrices. For our study of the multivariate polynomial systems the latter class is not sufficiently general, and we need our definition 2.5.1 due to [23] (cf. also [24]).

The two classes of quasi-Toeplitz and quasi-Hankel matrices can be immediately transformed into each other by their pre- or post-multiplicaiton by the reflection matrix, having its antidiagonal filled with ones and having all other entries zero.

Definition 2.5.2 Let π_{E} : L → L be the projection map such that

π_{E}(x^{α}) = x^{α}

if α ∈ E and π_{E}(x^{α}) = 0 otherwise. We also let π_{E} : ℂ[[∂]] → ℂ[[∂]] denote the projection map such that π_{E}(∂^{α}) = ∂^{α} is α ∈ E and π_{E}(∂^{α}) = 0 otherwise.

We can describe the quasi-Toeplitz and quasi-Hankel operators in terms of polynomial multiplication (see [24], [23]).

Proposition 2.5.3 The matrix M is an (E, F) quasi-Toeplitz (resp. an (E, F) quasi-Hankel) matrix, if and only if it is the matrix of the operator π_{E} ∘ M_{T_{M}} ∘ π_{F} (resp. π_{E} ∘ H_{H_{M}} ∘ π_{F}).

Proof. We will give a proof only for an (E, F) quasi-Toeplitz matrix M = (M_{α,β})_{α∈E,β∈F}. (The proof is similar for a quasi-Hankel matrix.) The associated polynomial is T_{M}(x) = ∑_{u∈E+F} t_{u}x^{u}. For any vector v = [v_{β}] ∈ ℂ^{F}, let v(x) denote the polynomial v(x) = ∑_{β∈F} v_{β}x^{β}. Then

T_{M}(x) v(x) = ∑_{u∈E+F,β∈F} x^{u+β} t_{u} v_{β}

= ∑_{α=u+β∈E+2 F} x^{α} (∑_{β∈F} t_{α−β} v_{β}) ,

where we assume that v_{β} = 0 if u ∉ E + F, t_{u} = 0 if u ∉ E + F. Therefore, for α ∈ E, the coefficient of x^{α} equals

∑_{β∈F} t_{α−β} v_{β} = ∑_{β∈F} M_{α,β} v_{β},

which is precisely the coefficient α of Mv. □

Due to proposition 2.5.3, multiplication of an (E, F) quasi-Toeplitz (resp. quasi-Hankel) matrix by a vector v = [v_{β}] ∈ ℂ^{F} reduces to (Laurent’s) polynomial multiplication.

Algorithm 2.5.1 (Multiplication of the (E, F) quasi-Toeplitz (resp. quasi-Hankel) matrix M = (M_{α,β})_{α∈E,β∈F} by a vector v = [v_{β}] ∈ ℂ^{F}). Multiply the two polynomials, T_{M} = ∑_{u∈E+F} t_{u} x^{u} (resp. H_{M}(∂) = ∑_{u∈E−F} h_{u}∂^{u}) by v(x) = ∑_{β∈F} v_{β}x^{β} (resp. v(∂^{−1}) = ∑_{β∈F} v_{β}∂^{−β}), and output the projection of the product on x^{E} (resp. ∂^{E}).

Hereafter, C_{P olMult}(E, F) denotes the number of ops required to multiply a polynomial with a support in E by a polynomial with a support in F. (We will estimate C_{P olMult}(E, F) in appendix B.) Algorithm 2.5.1 can be performed by using C_{P olMult}(E + F, F), resp. C_{P olMult}(E − F, −F), ops. According to the estimates of the appendix B, this means O(N log^{2} N + C_{M,N}) ops, where N = ⌊E − 2 F ⌉ (resp. ⌊E + 2F ⌉) and where C_{M,N} bounds the cost of evaluating the polynomial H_{M} (resp. T_{M}) on a fixed set of N points.

The displacement rank analysis (cf. [3]) can be also generalized to the multivariate case. Instead of the well-known displacement matrices

Z = ( \begin{matrix}0 · · · · · · · · · · · · & 0 \\ 1 ... & ⋮ \\ 0 ... ... & ⋮ \\ ⋮ ... ... ... & ⋮ \\ 0 · · · · · · 0 1 & 0 \end{matrix} )

and Z^{t}, we use the following operators (one per variable):

Z_{i}^{E} = π_{E} M_{x_{i}}π_{E}  (12)

and

Z_{−i}^{E} = π_{E} M_{x_{i}^{−1}}π_{E},  (13)

respectively. The displacement rank (that is, the rank of the matrix obtained by applying the displacement operator of the matrix) is bounded by the sum in i of the sizes of the boundary of E and F in the direction i (see [24], [23]).

Example Let the sets E and F correspond to the set of the monomials in x_{1}, x_{2} graphically represented as follows:

\begin{matrix}∘ & ∘ \\ ∘ & ∘ & ∘ & ∘ & ∘ \end{matrix}

Then the displacement rank is less than 2 × 2 = 4 in the direction x_{1} and is less than 2 × 5 = 10 in the direction x_{2}.

In other words, the flatter the sets E and F in a fixed direction, the smaller the displacement rank in this direction.

If E = F = {(α_{1}, …, α_{n}) ∈ ℕ^{n} ; 0 ≤ α_{i} ≤ d_{i} − 1}, the displacement rank for the operator associated to Z_{i} is ≤ 2 ∏_{j≠i} d_{j}.

2.6 Multivariate Bezoutians

In this section and in the next one, we will recall some basic definitions from the theories of Bezoutians and algebraic residues (compare the special univariate cases of sections 1.3 and 1.5), referring the reader to [8], [11] for further details and to section 3 for some applications.

In addition to the vector of variables x, consider another vector of varialbes, y = (y_{1}, …, y_{n}), and write x^{(0)} = x, x^{(1)} = (y_{1}, x_{2}, …, x_{n}), …, x^{(n)} = y. For a polynomial q ∈ R, define θ_{i}(q) = \frac{q(x^{(i)})−q(x^{(i−1)})}{y_{i}−x_{i}}, the discrete differentiation of q. For a sequence of n + 1 polynomials q, p_{1}, …, p_{n} ∈ R, construct the following polynomial in x and y:

Θ_{p}(q) = Θ_{q,p} = det ( \begin{matrix}q(x) & θ_{1}(q) & ⋯ & θ_{n}(q) \\ ⋮ & ⋮ & & ⋮ \\ p_{n}(x) & θ_{1}(p_{n}) & ⋯ & θ_{n}(p_{n}) \end{matrix} ) = ∑_{α,β} θ_{α,β}^{q,p} x^{α} y^{β},  (14)

where det(A) denotes the determinant of a matrix A, p = (p_{1}, …, p_{n}), and α and β vary in fixed ranges. This polynomial of ℂ[x, y] is called the Bezoutian of q, p_{1}, …, p_{n}. It defines a map B_{q,p}:

B_{q,p} : R̂ → R
Λ ↦ ∑_{α,β} θ_{α,β}^{q,p} x^{α} Λ(y^{β}).

Using the representation of Λ as a formal power series in ∂_{1}, …, ∂_{n}, we obtain the value of B_{q,p}(Λ(∂)) by computing the term free of ∂_{1}, …, ∂_{n} in the product

Θ_{q,p}(x, ∂^{−1})Λ(∂).

This construction extends the construction of section 1.3 to the multivariate case. The matrix of the map B_{q,p} in the monomial basis is the matrix of the coefficients [θ_{α,β}^{q}].

If (x^{α})_{α∈E} is a basis of A, then for any q in R, the polynomial Θ_{p}(q) can be rewritten as

Θ_{p}(q) ≡ ∑_{α,β∈E} B_{α,β}^{q,p} x^{α}y^{β}.  (15)

This polynomial is obtained from (14) by reducing Θ_{q,p} modulo I.

To simplify the notation, we will write B_{α,β}^{q}, dropping the superscript p for a fixed ideal (p).

Example (continued) We have

Θ_{p}(1) = x_{1}x_{2} + 2 x_{2}^{2} + (−2 y_{1} + y_{2}) x_{1} + (y_{1} + 2 y_{2} − 1) x_{2}
−2 y_{1}^{2} + y_{1}y_{2} + 16 y_{1} − y_{2}
≡ 5 x_{1}x_{2} + (y_{2} − 2 y_{1} + 14) x_{1} + (2 y_{2} + y_{1} − 1) x_{2}  (16)
+5 y_{1}y_{2} − y_{2} + 14 y_{1} − 4.

Definition 2.6.1 The matrix

B_{q,p} = [B_{α,β}^{q,p}]_{α,β∈E},  (17)

associated to the polynomial Θ_{p}(q) of (15), is called the Bezoutian matrix or the Bezoutian of q, p. This is the matrix of the map

\overline{B}_{q,p} : Â → A
Λ ↦ ∑_{α,β∈E} B_{α,β}^{q} x^{α}Λ(y^{β})

in the monomial basis (x^{α})_{α∈E} and its dual. When p is fixed, we will write B_{q} and \overline{B}_{q} instead of B_{q,p} and \overline{B}_{q,p}.

Example (continued) In the basis (1, x_{1}, x_{2}, x_{1} x_{2}) of A = ℂ[x_{1}, x_{2}]/(p_{1}, p_{2}) we have

B_{1} = [ \begin{matrix}−4 & 14 & −1 & 5 \\ 14 & −2 & 1 & 0 \\ −1 & 1 & 2 & 0 \\ 5 & 0 & 0 & 0 \end{matrix} ] .

The rows of this matrix are filled with the coefficients of the monomials in x_{1}, x_{2} in (16). It is a symmetric matrix, which is a property of the Bezoutians.

2.7 Bezoutians and algebraic residues

We will next define the residue and recall some fundamental properties of the multivariate Bezoutians and residues, to end with some correlations between primal and dual multiplication tables in the next section.

Definition 2.7.1 The residue of p = (p_{1}, …, p_{n}) is the unique linear form τ in the set of linear forms on R such that

1. τ vanishes on (p),

2. B_{1,p}(τ) − 1 ∈ (p).

This definition extends the characterization of the residue of proposition 1.5.2, given in the univariate case; we now consider all polynomials modulo the ideal (p), in particular, B_{p}(q) is modulo (p). This is not a constructive definition; we prove the existence of τ but give no general recipe for computing τ yet.

Consider the decomposition Θ_{1,p} ≡ ∑_{α,β∈E} B_{α,β}^{1} x^{α} y^{β} and let us write w_{α}(y) = ∑_{β∈E} B_{α,β}^{1}y^{β}, so that

Θ_{1,p} ≡ ∑_{α∈E} x^{α} w_{α}(y).

Then we have the following property:

Proposition 2.7.2 The set (w_{α})_{α∈E} is the dual basis of (x^{α}) for τ:

τ(x^{α} w_{β}) = \{ \begin{matrix}1 if α = β \\ 0 otherwise. \end{matrix}

Example (continued) The residue is defined on (1, x_{1}, x_{2}, x_{1} x_{2}) by

τ(1) = τ(x_{1}) = τ(x_{2}) = 0, τ(x_{1} x_{2}) = \frac{1}{5}.

It vanishes on all multiples of p_{1}, p_{2}. According to (16), the dual basis of (1, x_{1}, x_{2}, x_{1} x_{2}) is given by

w_{1} = 5 y_{1}y_{2} − y_{2} + 14 y_{1} − 4, w_{x_{1}} = y_{2} − 2 y_{1} + 14, w_{x_{2}} = 2 y_{2} + y_{1} − 1, w_{x_{1}x_{2}} = 5.

Again, we are going to study the properties of the dual basis but do not give yet any algorithm for actually computing this basis. According to proposition 2.7.2, for any a ∈ A, we have the relations

a ≡ ∑_{α∈E} τ(a x^{α}) w_{α} ≡ ∑_{α∈E} τ(a w_{α}) x^{α}.  (18)

We also have the following simple but fundamental property ([8], [11]):

Θ_{1,p} ≡ ∑_{α∈E} x^{α} w_{α}(y) ≡ ∑_{α∈E} w_{α}(x) y^{α} mod (p(x), p(y)),  (19)

which shows that B_{1} is a symmetric matrix.

Moreover, we recall from [8], [11] that for any polynomial q ∈ R we have

Θ_{q,p} = Θ_{1,p}(x, y) q(x) ≡ Θ_{1,p}(x, y) q(y) mod (p(x), p(y)).  (20)

In particular, we set q(x) = x_{i} for i = 1, …, n, and for any pair, ζ and η, of distinct roots of the polynomial system p = 0, we set x = ζ, y = η ∈ Z(I) and deduce that

Θ_{1,p}(ζ, η) = 0.  (21)

If ζ = η, then Θ_{1,p}(ζ, η) = J(ζ), where J = J_{p} = (∂p_{i}/∂x_{j}) is the Jacobian of p.

2.8 Bezoutians and multiplication tables in primal and dual bases

The notion of dual basis (for τ), defined in the previous section, should not be confused with the following notion of dual basis in the dual space Â:

Definition 2.8.1 Given a basis (b_{i})_{i=1,...,D} of A, let (b̂_{i})_{i=1,...,D} denote the dual basis of (b_{i}), that is, the basis set of linear forms in R̂ that compute the coefficients of any a ∈ A in the primal basis.

The next proposition relates the matrix B_{a} = B_{a,p} of the map \overline{B}_{a} of definition 2.6.1 with q = a (cf. (17)), to the transformations between the primal bases (x^{α}) and (w_{α}) and their dual bases (x̂^{α}) and (ŵ_{α}), respectively.

Proposition 2.8.2 The matrix of the map \overline{B}_{a} of definition 2.6.1,

1. from the basis (x̂^{α}) of Â to the basis (x^{α}) of A is B_{a} = (τ(a w_{α} w_{β})),

2. from the basis (ŵ_{α}) to the basis (w_{α}) is H_{a} = (τ(a x^{α} x^{β})).

Proof. According to proposition 2.7.2, the coordinates of \overline{B}_{a}(x̂^{β}) in the basis (x^{α})_{α∈E} are given by

τ(\overline{B}_{a}(x̂^{β}) w_{α}).

The identities (20) and (19) imply that Θ_{p}(a) ≡ a(x)Θ_{p}(1), and \overline{B}_{a}(x̂^{β}) ≡ a \overline{B}_{1}(x̂^{β}) ≡ a w_{β}. Therefore,

τ(\overline{B}_{a}(x̂^{β}) w_{α}) = τ(a \overline{B}_{1}(x̂^{β}) w_{α}) = τ(a w_{α} w_{β}).

In other words, we have B_{α,β}^{a} = τ(\overline{B}_{a}(x̂^{β}) w_{α}). This proves the first part of the proposition.

The coordinates of \overline{B}_{a}(ŵ_{β}) in the basis (w_{α})_{α∈E} are given by

τ(\overline{B}_{a}(ŵ_{β}) x^{α}).

According to identities (20) and (19), we also have

τ(\overline{B}_{a}(ŵ_{β}) x^{α}) = τ(a \overline{B}_{1}(ŵ_{β}) x^{α}) = τ(a x^{α} x^{β}),

which proves the second part of the proposition. □

Now, we deduce some simple correlations between multiplication tables in the bases (x^{α}) and (w_{α}).

Definition 2.8.3 For any a in A, let M_{a} = (M_{α,β}^{a}) denote the matrix of the map \overline{M}_{a} in the basis (x^{α}) and let N_{a} = (N_{α,β}^{a})_{α,β∈E} denote its matrix in the basis (w_{α}).

Proposition 2.8.4 The matrix N_{a} of multiplication by a in A, in the basis (w_{α}), is the transpose M_{a}^{t} of the matrix M_{a} of multiplication by a in A, in the basis (x^{α}).

Proof. For any α ∈ E, we have

b x^{β} ≡ ∑_{γ∈E} M_{γ,β}^{a} x^{γ}, b w_{β} ≡ ∑_{γ∈E} N_{γ,β}^{a} w_{γ},

and

\begin{matrix}M_{α,β}^{a} = τ(b x^{β} w_{α}), \\ N_{α,β}^{a} = τ(a x^{α} w_{β}). \end{matrix}

Therefore, N_{a} = M_{a}^{t}. □

The proposition also implies that the matrix of the transposed map \overline{M}_{a}^{t} in the dual basis (x̂^{α}) of (w_{α}) is M_{a}.

2.9 Multivariate (stipped) Vandermonde matrices

Vandermonde matrices can be immediately generalized to the multivariate case, in the following way.

Definition 2.9.1 For a set (x^{α})_{α∈E} of D monomials and a set ξ = (ξ_{1}, …, ξ_{D}) of D points of ℂ^{n}, define the Vandermonde matrix of ξ on E by

V_{E}(ξ) = [ξ_{i}^{α}]_{i=1,...,D; α∈E}.

The rows of this matrix are the vectors [x^{α}]_{α∈E} of monomials evaluated at points ξ_{i} (for i = 1, …, D).

V_{E}(ξ) is the matrix of the coefficients (of (∂^{α})_{α∈E}) in the f.p.s. representing the evaluations 1_{ξ_{i}} at the points ξ_{i}. Multivariate Vandermonde matrices are also called stripped because the powers ξ_{i}^{α} do not generally fill the entire D-dimensional box of integers.

Algorithm 2.9.1 (Multiplication of a multivariate Vandermonde matrix V_{E}(ξ) by a vector v and the solution in v of a linear system V_{E}(ξ)v = w, for given ξ, E and w). Perform multipoint evaluation at the node-points ξ_{i}, i = 1, …, D, of the associated multivariate polynomial with the coefficient vector w (resp. perform the converse operation of multivariate polynomial interpolation).

See [6] and [12], for a record (asymptotic) bounds on the number of ops involved in algorithm 2.9.1.

2.10 Relations between Hankel and Bezoutian matrices

Motivated by applications of matrix computations to the solution of polynomial systems, we are particularly interested in studying multiplication tables (see theorems 2.2.1, 2.4.1).

Definition 2.10.1 For any Λ in Â, let H_{Λ} denote the quasi-Hankel matrix of residues,

H_{Λ} = (Λ(x^{α+β}))_{α,β∈E}.

For any element a in A, we will also write H_{a} = H_{a⋆τ}, where τ is the residue of p.

Let us extend corollary 1.5.4, by relating the multivariate Bezoutian B_{1} with the quasi-Hankel matrix of residues H_{1}.

Theorem 2.10.2 The inverse of H_{1} is B_{1}.

Proof. By definition, w_{α}(x) = ∑_{β∈E} B_{α,β}^{1} x^{β}. Therefore, by using proposition 2.7.2, we obtain that

τ(w_{α} x^{β}) = ∑_{γ∈E} B_{α,γ}^{1} τ(x^{γ+β})

equals 1 if α = β and is 0 otherwise. This is precisely the coefficient (α, β) of the matrix B_{1} H_{1}. Thus, we have

B_{1} H_{1} = 𝕀_{D},

where 𝕀_{D} is the D × D identity matrix. □

Example (continued) We have

τ(1) = τ(x_{1}) = τ(x_{2}) = 0,

τ(x_{1}x_{2}) = \frac{1}{5}, τ(x_{1}^{2}) = −\frac{2}{5}, τ(x_{2}^{2}) = \frac{2}{5}, τ(x_{1}^{2}x_{2}) = \frac{29}{25}, τ(x_{1}^{2}x_{2}^{2}) = −\frac{12}{25}, τ(x_{1}^{2}x_{2}^{2}) = −\frac{398}{125},

and

H_{1} = [ \begin{matrix}0 & 0 & 0 & \frac{1}{5} \\ 0 & −\frac{2}{5} & \frac{1}{5} & \frac{29}{25} \\ 0 & \frac{1}{5} & \frac{2}{5} & −\frac{12}{25} \\ \frac{1}{5} & \frac{29}{25} & −\frac{12}{25} & −\frac{398}{125} \end{matrix} ] .

Here is the polynomial associated to this quasi-Hankel matrix:

P = \frac{1}{5}∂_{1}∂_{2} − \frac{2}{5}∂_{1}^{2} + \frac{2}{5}∂_{2}^{2} + \frac{29}{25}∂_{1}^{2}∂_{2} − \frac{12}{25}∂_{1}∂_{2}^{2} − \frac{398}{125}∂_{1}^{2}∂_{2}^{2}.

To illustrate the multiplication property, we compute for instance the coordinates of the vector [1, 0, −1, 0]^{T} H_{1}, which are the coefficients of 1, ∂_{1}, ∂_{2}, ∂_{1}∂_{2} in the next product:

(1 − ∂_{2}^{−1}) P =

2 ∂_{1}^{2}∂_{2}^{−1} − ∂_{1} − 2 ∂_{2} − \frac{39}{5} ∂_{1}^{2} + \frac{17}{5} ∂_{2}∂_{1} + 2 ∂_{2}^{2} + \frac{543}{25} ∂_{2}∂_{1}^{2} − \frac{12}{5} ∂_{2}^{2}∂_{1} − \frac{398}{25} ∂_{2}^{2}∂_{1}^{2},

which yields the vector [0, −1, −2, \frac{17}{5}]. We may verify that H_{1} is the inverse of the Bezoutian B_{1} of the example of section 2.6:

[ \begin{matrix}−4 & 14 & −1 & 5 \\ 14 & −2 & 1 & 0 \\ −1 & 1 & 2 & 0 \\ 5 & 0 & 0 & 0 \end{matrix} ] [ \begin{matrix}0 & 0 & 0 & \frac{1}{5} \\ 0 & −\frac{2}{5} & \frac{1}{5} & \frac{29}{25} \\ 0 & \frac{1}{5} & \frac{2}{5} & −\frac{12}{25} \\ \frac{1}{5} & \frac{29}{25} & −\frac{12}{25} & −\frac{398}{125} \end{matrix} ] = [ \begin{matrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{matrix} ]

The matrices B_{1} and H_{1} express the transformation from the basis (x^{α}) to the dual basis (w_{α})_{α∈E}:

Proposition 2.10.3 For any a ∈ A, if v is the coordinate vector of a in the monomial basis (x^{α})_{α∈E} and w is the coordinate vector of a in the dual basis (w_{α})_{α∈E}, then we have

v = B_{1} w, w = H_{1} v.

Let us relate the matrices above to multiplication tables (compare section 2.8).

Proposition 2.10.4 For any linear form Λ ∈ Â and any a ∈ A, we have

H_{a⋆Λ} = M_{a}^{t}H_{Λ} = H_{Λ} M_{a},  (22)

where M_{a} is the matrix of definition 2.8.3. In particular, we have

H_{a} = H_{1} M_{a} = M_{a}^{t} H_{1}.  (23)

Proof. For any pair a, b ∈ R, we define the operator

H_{a⋆Λ}(b) = b ⋆ (a ⋆ Λ) = a b ⋆ Λ = H_{Λ}(a b)
= a ⋆ (b ⋆ Λ) = a ⋆ H_{Λ}(b).

Therefore, the operator H_{a⋆Λ} can be decomposed as follows:

H_{a⋆Λ} = H_{Λ} ∘ M_{a} = M_{a}^{t} ∘ H_{Λ}.

In terms of matrices, this yields equations (22). □

A similar relation is also valid for the Bezoutian matrices (see definition 2.6.1):

Theorem 2.10.5 For any a ∈ A, we have

B_{a} = B_{1} M_{a}^{t} = M_{a} B_{1}.  (24)

Proof. According to (20), in terms of operators (see definition 2.6.1 with a = q) we have ∀Λ ∈ Â that

\overline{B}_{a}(Λ) = ∑_{α,β∈E} B_{α,β}^{a} x^{α} Λ(y^{β})
= a(x) ∑_{α,β∈E} B_{α,β}^{1}x^{α} Λ(y^{β}) = a(x) \overline{B}_{1}(Λ)
= ∑_{α,β∈E} B_{α,β}^{1} x^{α} Λ(a(y) y^{β}) = \overline{B}_{1}(a ⋆ Λ).

Thus, we can decompose the map \overline{B}_{a} as follows:

\overline{B}_{a} = \overline{M}_{a} ∘ \overline{B}_{1} = \overline{B}_{1} ∘ \overline{M}_{a}^{t}.

In terms of matrices, this implies (24). □

According to proposition 2.10.3, the latter theorem can be also reformulated as follows: For any pair of elements a and b of A, let v be the coordinate vector of b in (w_{α})_{α∈E}. Then the coordinate vector of a b in the monomial basis (x^{α})_{α∈E} is B_{a}v.

We will use relations (24) and (23) in section 3, in order to transform the eigenproblem of sections 2.2 and 2.3 into a generalized structured eigenproblem (see in particular our demonstration in section 3.1.2).

Proposition 2.10.6 If a b ≡ 1 in A, then

B_{a} H_{b} = B_{b} H_{a} = 𝕀_{D}.

Proof. According to (24), we have

B_{a}H_{b} = B_{1} M_{a}^{t} M_{b}^{t} H_{1}B_{1}H_{1} = 𝕀_{D},

for M_{a} M_{b} = M_{a b} = 𝕀_{D}. Similarly, we deduce that B_{b} H_{a} = 𝕀_{D}. □

Proposition 2.10.7 For any a ∈ A, we have the relations

• B_{a} = B_{1} H_{a} B_{1},

• H_{a} = H_{1} B_{a} H_{1}.

Proof. According to (24) and (23) and proposition 2.10.2, we have

B_{a} = B_{1}M_{a}^{t} and M_{a}^{t} = H_{a}H_{1}^{−1} = H_{a}B_{1},

which implies the first relation of this proposition. The other relation is obtained by inverting the first one and applying proposition 2.10.6. □

2.11 Some relations between the values of the Jacobian at the roots and structured matrices, in the case of simple roots

Let us assume that the roots ζ ∈ Z are simple. Then J(ζ_{i}) ≠ 0, where J = J_{p} = (\frac{∂p_{i}}{∂x_{j}} ) is the Jacobian of p = (p_{1}, …, p_{n}).

Let V_{E}(Z) be the Vandermonde matrix, defined in section 2.9. We recall that for any vector v = [v_{α}]_{α∈E}, the product V_{E}(Z) v is the vector [v(ζ_{1}), …, v(ζ_{D})] of the evaluations of the polynomial v(x) = ∑_{α} v_{α} x^{α} at the roots ζ_{i} ∈ Z(I).

Proposition 2.11.1 For any polynomial a ∈ R, we have

B_{a} = V_{E}(Z)^{−1} diag (a(ζ_{1}) J(ζ_{1}), …, a(ζ_{D}) J(ζ_{D})) V_{E}(Z)^{−t},

where diag (l_{1}, …, l_{D}) represents the D×D diagonal matrix, with the diagonal entries l_{1}, …, l_{D}.

Proof. As the rows of V_{E}(Z) are given by the values of the monomial vector [x^{α}] at the roots ζ_{i} ∈ Z(I), the matrix V_{E}(Z) B_{a} V_{E}^{t}(Z) is the matrix

[Θ_{a,p}(ζ_{i}, ζ_{j})]_{i,j=1,...,D}.

According to equations (20), (21), we have Θ_{a,p}(ζ, η) = Θ_{a,p}(ζ, η) = 0 if ζ ≠ η.

If η = ζ, then, by construction, Θ_{a,p}(a)(ζ, ζ) = a(ζ) J(ζ). Consequently, (Θ_{1,p}(ζ_{i}, ζ_{j})) is the diagonal matrix

diag (a(ζ_{1}) J(ζ_{1}), …, a(ζ_{D}) J(ζ_{D})) .
□

Corollary 2.11.2 If the roots of the system p = 0 are simple, then

H_{1} = V_{E}(Z)^{t} diag (\frac{1}{J(ζ_{1})}, …, \frac{1}{J(ζ_{D})}) V_{E}(Z).

Proof. We have B_{1} = V_{E}(Z)^{−1} diag (J(ζ_{1}), …, J(ζ_{D})) V_{E}(Z)^{−t}, according to proposition 2.11.1, and we deduce from theorem 2.10.2 that

H_{1} = B_{1}^{−1} = V_{E}(Z)^{t}diag (\frac{1}{J(ζ_{1})}, …, \frac{1}{J(ζ_{D})}) V_{E}(Z).

□

If we substitute these relations into (24), we obtain the following property:

Corollary 2.11.3 If the roots of the system p = 0 are simple, then

M_{a} = V_{E}^{−1}(Z)diag (a(ζ_{1}), …, a(ζ_{d})) V_{E}(Z).  (25)

According to theorem 2.10.5, we have H_{a} = H_{1} M_{a}, which yields:

Corollary 2.11.4 If the roots of the system p = 0 are simple, then

H_{a} = V_{E}(Z)^{t}diag (\frac{a(ζ_{1})}{J(ζ_{1})}, …, \frac{a(ζ_{D})}{J(ζ_{D})}) V_{E}(Z).  (26)

2.12 Relations between Bezoutians and idempotents

As in section 2.11, we assume that the roots ζ ∈ Z are simple and denote by J the Jacobian of p. For any ζ ∈ Z, J(ζ) ≠ 0.

Proposition 2.12.1 If the roots of the system p = 0 are simple, then the vectors

e_{ζ} = \frac{1}{J(ζ)}Θ_{1,p}(x, ζ), ζ ∈ Z,

form a basis, consisting of orthogonal idempotents of A, whose sum equals 1, that is, e_{ζ}^{2} ≡ e_{ζ}, e_{ζ}e_{η} ≡ 0 if ζ ≠ η, and ∑_{ζ∈Z(I)} e_{ζ} ≡ 1.

Proof. According to equation (20), for any q ∈ R and any ζ ∈ Z(I), we have

Θ_{1,p}(x, ζ) q(x) ≡ Θ_{1,p}(x, ζ) q(ζ)

in the quotient ring B. Therefore,

Θ_{1,p}(x, ζ) Θ_{1,p}(x, ζ) ≡ J(ζ)Θ_{1,p}(x, ζ),

and e_{ζ} = \frac{1}{J(ζ)} Θ_{1,p}(x, ζ) ≡ e_{ζ}^{2} is an idempotent (J(ζ) ≠ 0 assuming all roots of the system p = 0 are simple). Moreover, according to (21), we have

Θ_{1,p}(x, ζ) Θ_{1,p}(x, η) ≡ Θ_{1,p}(x, ζ)Θ_{1,p}(ζ, η) ≡ 0,

for any pair of distinct roots ζ, η ∈ Z(I), which shows that e_{ζ} e_{η} ≡ 0 unless ζ = η. We obtain from the definition of the residue τ and from the Euler-Jacobi identity (cf. [11]) that

Θ_{1,p}(τ) ≡ 1 (by definition 2.7.1)
≡ ∑_{ζ∈Z} \frac{1}{J(ζ)} Θ_{1,p}(x, ζ) ≡ ∑_{ζ∈Z} e_{ζ} (by the Euler−Jacobi identity).

This shows that the sum of the idempotents equals 1 in A, and thus they form a basis of A (which is of dimension D). □

Now let us recover the root ζ from the idempotent e_{ζ}. By definition, we have

e_{ζ} = \frac{1}{J(ζ)}Θ_{1,p}(x, ζ) = \frac{1}{J(ζ)} ∑_{α∈E} x^{α} (∑_{β} B_{α,β}^{1} ζ^{β}) ,

so that the coordinate vector [e_{ζ}] of e_{ζ} in the basis (x^{α})_{α∈E} is

[e_{ζ}] = \frac{1}{J(ζ)}B_{1} [ζ^{α}]_{α∈E}.

Equivalently, we have

[ζ^{α}]_{α∈E} = J(ζ) H_{1} [e_{ζ}].  (27)

Corollary 2.12.2 The coordinates of e_{ζ} in the dual basis (w_{α}) are \frac{1}{J(ζ)}[ζ^{α}].

Algorithm 2.12.1 (Recovery of the root vector ζ from the idempotent vector e_{ζ}). Multiply e_{ζ} by the quasi-Hankel matrix H_{1} and output the ratios of the coordinates of the resulting vector.

Let us estimate the computational cost of performing the algorithm. If v = H_{1}[e_{ζ}] = \frac{1}{J(ζ)}[ζ^{α}]_{α∈E} = [v_{1}, v_{x_{1}}, …, v_{x_{n}}, v_{x_{1}^{2}}, …], then the i-th coordinate of ζ is

ζ_{i} = \frac{v_{x_{i}}}{v_{1}} .

Thus, the roots can be computed from the idempotent e_{ζ} in at most C_{P olMult}(E, 2E) ops, by using algorithm 2.5.1 applied for F = 2E.

3 Applications

In this section, we exploit the properties of and the relations among multivariate polynomial roots, dual space of linear maps, residue and structured matrices in order to devise algorithms for solving polynomial systems of equations. First we focus on structured generalized eigenproblem, involving quasi-Hankel and Bezoutian matrices. Then we consider quasi-Toeplitz matrices that generalize the Sylvester matrices. They are used for computing a basis (x^{α})_{α∈E} of A, the multiplication tables, and the first coefficients of the dual basis of (x^{α})_{α∈E}, for generic input. Using the machinery of the previous section enables us to yield better insight into the subject and simplify substantially the proofs of some known fundamental results. Finally, we focus on iterative methods converging to idempotents (and based on using quasi-Hankel matrices) and on application of structured matrices to counting distinct roots and real roots of a polynomial system. In this part, we improve dramatically the known computational complexity estimates, though the algorithms are proposed in preliminary form and require further elaboration for their implementation.

3.1 Reduction of solving a polynomial system to matrix eigenproblems

Let us restate theorem 2.2.1 and 2.4.1 in terms of matrices rather than their associated operators. For a fixed element a ∈ A, we consider the operator of multiplication by a :

\begin{matrix}\overline{M}_{a} : A & → & A \\ b & ↦ & a b, \end{matrix}

whose matrix in the monomial basis (x^{α})_{α∈E} is denoted M_{a}. The transposed operator from Â to Â is defined by the map:

\begin{matrix}\overline{M}_{a}^{t} : Â & → & Â \\ Λ & ↦ & a ⋆ Λ = Λ ∘ \overline{M}_{a}, \end{matrix}

and its matrix in the dual basis is M_{a}^{t}. We have the following theorem, whose first two parts restate theorems 2.2.1 and 2.4.1 in terms of matrices (see [1], [22]):

Theorem 3.1.1

1. The eigenvalues of the matrices of the linear operators M_{a} and M_{a}^{t} are {a(ζ_{1}), …, a(ζ_{d})}.

2. The common eigenvectors of the matrices (M_{x_{i}}^{t})_{i=1,...,n} are (up to a scalar) [ζ_{i}^{α}]_{α∈E}.

3. If n = m, then the common eigenvectors of the matrices (M_{x_{i}})_{i=1,...,n} are (up to a scalar factor) J(x) e_{1}, …, J(x) e_{d}, where J(x) is the Jacobian of p_{1}, …, p_{n}, and e_{i} are the idempotents associated with the roots.

Part 1 amounts to theorem 2.2.1. Part 2 is deduced from theorem 2.4.1: the coordinates of the evaluation 1_{ζ_{i}} at the root ζ_{i} in the dual basis of (x^{α})_{α∈E} are precisely [ζ_{i}^{α}]_{α∈E}. The third part is proved in [11].

As a consequence of theorem 3.1.1, we may compute easily the roots ζ_{i} from the eigenvectors of M_{x_{i}}^{t}, as in algorithm 2.12.1:

Proposition 3.1.2 If (x^{α})_{α∈E} = (1, x_{1}, …, x_{n}, …) contains the monomials 1, x_{1}, …, x_{n} and if v = [v_{α}]_{α∈E} = (v_{1}, v_{x_{1}}, …, v_{x_{n}}, …,) is a common vector of (M_{x_{i}})_{i=1,...,n}, then

ζ = (\frac{v_{x_{1}}}{v_{1}} , …, \frac{v_{x_{n}}}{v_{1}} )

is a root of p = 0.

Algorithm 3.1.1 (Computation of the roots of the polynomial system p = 0). Compute and output the roots as the scaled common eigenvectors of the matrices M_{a}^{t} for a ∈ R.

Example (continued) Here is the normalized matrix V of eigenvectors (with eight digit accuracy) of the matrix M_{x_{1}}^{t}:

[ \begin{matrix}1.0 & 1.0 & 1.0 & 1.0 \\ 6.8200982 & −0.19395427 + 0.20520688 i & −0.19395427 − 0.20520688 i & 0.36781361 \\ −2.8367388 & −0.61937124 − 1.3895199 i & −0.61937124 + 1.3895199 i & 1.6754769 \\ −19.346814 & 0.40526841 + 0.14240419 i & 0.40526841 − 0.14240419 i & 0.61626304 \end{matrix} ] .

The columns of this matrix are the vectors [ζ_{i}^{α}]_{α∈E} for ζ_{i} ∈ Z(I). Thus, we immediately deduce that the four roots of p_{1}(x_{1}, x_{2}) = p_{2}(x_{1}, x_{2}) = 0 are given by the next table:

\begin{matrix}ζ_{1} & ζ_{2} & ζ_{3} & ζ_{4} \\ 6.8200982 & −0.19395427 + 0.20520688 i & −0.19395427 − 0.20520688 i & 0.36781361 \\ −2.8367388 & −0.61937124 − 1.3895199 i & −0.61937124 + 1.3895199 i & 1.6754769 \end{matrix} .

We immediately check that V_{2,i}V_{3,i} = V_{4,i} for i = 1, 2, 3, 4. The structure of these multiplication matrices is not easy to analyze directly, however. Thus, we will multiply the matrices M_{a}^{t} by two fixed invertible matrices A and B in order to transform the problem into an equivalent generalized eigenproblem, (AM_{a}^{t}B − λ A B)v = 0, where the structure can be explicitly exploited. We will give some examples of such a reduction involving structured matrices.

3.1.1 Reduction by using Hankel matrices

According to (22), for any Λ ∈ Â and any a ∈ R, we have

H_{a⋆Λ} = M_{a}^{t}H_{Λ},

so that solving the eigenproblem (H_{a⋆Λ} − λH_{Λ})v = 0 yields the eigenvector H_{Λ} v of M_{a}^{t}. Let us next exploit this matrix equation assuming that we have a normal form algorithm Nf that projects R onto ⟨x^{α}⟩_{α∈E} along I, that is, computes the unique element of ⟨x^{α}⟩_{α∈E} in the same class modulo I.

Algorithm 3.1.2 (Solution of a polynomial system via the solution of a generalized eigenvector problem defined by using Hankel matrices). Fix two exponents α_{0}, α_{1} ∈ E. Then proceed as follows:

1. For all monomials x^{α+β} with α, β ∈ E, compute in the normal form Nf(x^{α+β}) of x^{α+β}:

• the coefficient of x^{α_{0}}, which we denote σ_{0}(x^{α+β}),

• the coefficient of x^{α_{1}}, which we denote σ_{1}(x^{α+β}).

2. Construct the two quasi-Hankel matrices:

• H_{σ_{0}} = (σ_{0}(x^{α+β}))_{α,β∈E},

• H_{σ_{1}} = (σ_{1}(x^{α+β}))_{α,β∈E}.

3. Solve the polynomial system p=0 via the solution of the generalized eigenvector

problem:

(H_{σ_{1}} − λH_{σ_{0}}) v = 0.  (28)

Let us specify stage 3. The coefficient of x^{α} in A (for any α ∈ E) is output by the linear form p → τ(w_{α} p) = w_{α} ⋆ τ(p). Thus, we have

H_{σ_{i}} = M_{w_{α_{i}}}^{t}H_{1},

for i = 0, 1. Therefore, if v is a generalized eigenvector of (28), then ṽ = H_{1}v is a generalized eigenvector of (M_{w_{α_{1}}}^{t} − λM_{w_{α_{0}}}^{t}) ṽ = 0, and its associated eigenvalue is

equal to \frac{w_{α_{1}}(ζ)}{w_{α_{0}}(ζ)} (if w_{α_{0}}(ζ) ≠ 0) for one of the roots ζ ∈ Z(I).

According to theorem 3.1.1, the common eigenvectors of M_{w_{α_{1}}}^{t} − λM_{w_{α_{0}}}^{t} for all pairs α_{0}, α_{1} ∈ E are the multiples of the vectors [ζ^{α}]_{α∈E} for ζ ∈ Z(I). The roots ζ are easily computed from these vectors, by using algorithm 3.1.1.

Example (continued) Suppose that we have computed the following normal forms in the basis (1, x_{1}, x_{2}, x_{1} x_{2}) of A = ℂ[x_{1}, x_{2}]/(p_{1}, p_{2}):

Nf(1) = 1, Nf(x_{1}) = x_{1}, Nf(x_{2}) = x_{2}, Nf(x_{1}x_{2}) = x_{1}x_{2},

Nf(x_{1}^{2}) = 1 + x_{1} − 2 x_{1}x_{2}, Nf(x_{2}^{2}) = −1 + 7 x_{1} + 2 x_{1}x_{2},

Nf(x_{2}x_{1}^{2}) = −\frac{14}{5} − \frac{12 x_{1}}{5} + \frac{x_{2}}{5} + \frac{29 x_{1}x_{2}}{5} ,

Nf(x_{1}x_{2}^{2}) = \frac{7}{5} + \frac{6 x_{1}}{5} + \frac{2 x_{2}}{5} − \frac{12 x_{1}x_{2}}{5} ,

Nf(x_{1}^{2}x_{2}^{2}) = \frac{198}{25} + \frac{209 x_{1}}{25} − \frac{12 x_{2}}{25} − \frac{398 x_{1}x_{2}}{25} .

We choose the monomial x^{α_{0}} = x_{1} x_{2} and x^{α_{1}} = x_{2}, which yields the following matrices:

H_{σ_{0}} = [ \begin{matrix}0 & 0 & 0 & 1 \\ 0 & −2 & 1 & \frac{29}{5} \\ 0 & 1 & 2 & −\frac{12}{5} \\ 1 & \frac{29}{5} & −\frac{12}{5} & −\frac{398}{25} \end{matrix} ] , H_{σ_{1}} = [ \begin{matrix}0 & 0 & 1 & 0 \\ 0 & 0 & 0 & \frac{1}{5} \\ 1 & 0 & 0 & \frac{2}{5} \\ 0 & \frac{1}{5} & \frac{2}{5} & −\frac{12}{25} \end{matrix} ] ,

and we obtain

H_{σ_{1}} H_{σ_{0}}^{−1} = [ \begin{matrix}−\frac{1}{5} & \frac{1}{5} & \frac{2}{5} & 0 \\ \frac{1}{5} & 0 & 0 & 0 \\ −\frac{2}{5} & \frac{14}{5} & −\frac{1}{5} & 1 \\ 0 & 0 & \frac{1}{5} & 0 \end{matrix} ] .

We have σ_{0} = w_{α_{0}} ⋆ τ = (2 x_{2} + x_{1} − 1) ⋆ τ and σ_{α_{1}} = w_{α_{1}} ⋆ τ = 5 τ. Therefore, H_{σ_{0}} = 5 H_{1}, and H_{σ_{1}} = H_{2 x_{2}+x_{1}−1}, so that

H_{σ_{1}} H_{σ_{0}}^{−1} = M_{\frac{1}{5} (2 x_{2}+x_{1}−1)}^{t}.

Indeed, the first row of the latter matrix represents the polynomial \frac{1}{5}(2 x_{2} + x_{1} − 1), the second row is x_{1} × \frac{1}{5}(2 x_{2} + x_{1} − 1), which is reduced to \frac{1}{5} in A. This implies that x_{1}^{−1} ≡ (2 x_{2} + x_{1} − 1).

3.1.2 Reduction by using Bezoutian matrices

The relations (23) on Bezoutians imply that

B_{a} = B_{1} M_{a}^{t}.

As in algorithm 3.1.2, assume that we have a normal form algorithm that computes an element in A reduced modulo I.

Algorithm 3.1.3 (Solution of a polynomial system via the solution of a generalized eigenvector problem defined by using Bezoutian matrices).

1. Compute the polynomials Θ_{1,p} and Θ_{x_{1},p} and their normal forms in x and y.

2. Compute the matrices B_{1} and B_{x_{1}} associated with these normal forms.

3. Solve the polynomial system p=0 via the solution of the generalized eigenvector

problem

(B_{x_{1}} − λ B_{1}) v = 0.

The generalized eigenvector of the pencil (B_{x_{1}}, B_{1}) (computed at stage 3) yields immediately the eigenvectors [ζ_{i}^{α}]_{α∈E}, and then scaling immediately gives us the coordinates of the roots ζ_{i} (cf. algorithm 3.1.1).

Example (continued) The Bezoutian of x_{1} is

B_{x_{1}} = [ \begin{matrix}0 & −2 & 1 & 0 \\ −2 & 12 & 0 & 5 \\ 1 & 0 & 0 & 0 \\ 0 & 5 & 0 & 0 \end{matrix} ] and B_{1}^{−1} B_{x_{1}} = M_{x_{1}}^{t} = [ \begin{matrix}0 & 1 & 0 & 0 \\ 1 & 1 & 0 & −2 \\ 0 & 0 & 0 & 1 \\ −\frac{14}{5} & −\frac{12}{5} & \frac{1}{5} & \frac{29}{5} \end{matrix} ] .

The first row of this matrix represents multiplication by x_{1} in A, the second represents multiplication by x_{1}^{2}, and so on.

3.2 Computation of multiplication matrices and the dual space

3.2.1 Sylvester’s matrices

As a basic pattern, we will first revisit the construction of the well-known Sylvester matrix in the univariate case.

Given two univariate polynomials, p_{0} = p_{0,0} + ⋯ + p_{0,d_{0}} x^{d_{0}} of degree d_{0} and p_{1} = p_{1,0} + … + p_{1,d_{1}} x^{d_{1}} of degree d_{1}, we will define the multiplication by p_{0} modulo p_{1} by the map:

\overline{M}_{p_{0}} : A → A
a ↦ a p_{0},

in the basis ⟨1, …, x^{d_{1}−1}⟩ of A = ℂ[x]/(p_{1}). The matrix of this map is defined via the Sylvester matrix S of p_{0} and p_{1}, that is, the matrix of the coefficients of the polynomials

p_{0}, x p_{0}, …, x^{d_{1}−1} p_{0}, p_{1}, x p_{1}, …, x^{d_{0}−1} p_{1}

in the monomial basis. The matrix S takes the following form:

\overbrace{p_{0} ⋯ x^{d_{1}−1}p_{0} p_{1} ⋯ x^{d_{0}−1}p_{1}}^{d_{0}+d_{1}} ; [ \begin{matrix}p_{0,0} & & 0 & p_{1,0} & & p_{1,1−d_{0}} \\ ⋮ & ... & & ⋮ & ... \\ p_{0,d_{1}−1} & ⋯ & p_{0,0} & p_{1,d_{1}−1} & ⋯ & p_{1,d_{1}−d_{0}} \\ p_{0,d_{1}} & ⋯ & p_{0,1} & p_{1,d_{1}} & ⋯ & p_{1,d_{1}−d_{0}+1} \\ ⋮ & & ⋮ & & ... & ⋮ \\ p_{0,d_{0}+d_{1}−1} & ⋯ & p_{0,d_{0}} & 0 & & p_{1,d_{1}} \end{matrix} ] \begin{matrix}1 \\ x \\ · \\ x^{d_{1}−1} \\ · \\ · \\ x^{d_{0}+d_{1}−1} \end{matrix} \} d_{0} + d_{1} (29)

under the convention that p_{0,i} = 0 if i > d_{0}, p_{1,j} = 0 if j < 0. Let V_{0}, V_{1}, and V denote the vector spaces generated by the monomials {1, …, x^{d_{1}−1}}, {1, …, x^{d_{0}−1}}, and {1, …, x^{d_{0}+d_{1}−1}}, respectively. Then the Sylvester matrix is the matrix of the map

S : V_{0} × V_{1} → V
(q_{0}, q_{1}) ↦ p_{0} q_{0} + p_{1} q_{1},

in the corresponding monomial basis. The determinant of this (d_{0} + d_{1}) × (d_{0} + d_{1}) matrix is the resultant of p_{0} and p_{1}.

For computing the matrix M_{p_{0}} of the multiplication by p_{0} modulo p_{1}, we have to reduce the polynomials p_{0}, x p_{0}, …, x^{d_{1}−1} p_{0} modulo p_{1}, by subtracting some multiples of p_{1}, to obtain linear combinations of the monomial basis (1, …, x^{d_{1}−1}) of A. The partition of the Sylvester matrix into four blocks as in (29),

S = [ \begin{matrix}U & V \\ Z & W \end{matrix} ] ,

enables us to interpret these operations in terms of matrix operations and thus to

analyze the structure of the matrix of multiplication. The block P_{0} = [ \begin{matrix}U \\ Z \end{matrix} ] rep-

resents the multiples of p_{0}, and the block P_{1} = [ \begin{matrix}V \\ W \end{matrix} ] represents the multiples of

p_{1}. Therefore, reducing the multiples of p_{0} by p_{1} consists in subtracting some linear combinations of the columns of P_{1} from the columns of P_{0} so that Z is replaced by a zero block. These operations on the columns of the Sylvester matrix are given explicitly by the following formula:

[ \begin{matrix}U & V \\ Z & W \end{matrix} ] [ \begin{matrix}𝕀_{d_{1}} & 0 \\ −W^{−1}Z & 𝕀_{d_{0}} \end{matrix} ] = [ \begin{matrix}U − V W^{−1}Z & V \\ 0 & W \end{matrix} ] ,

and we have the following property:

Proposition 3.2.1 The matrix M_{p_{0}} of multiplication by p_{0} modulo p_{1} in the monomial basis ⟨1, x, …, x^{d_{1}−1}⟩ is the Schur complement of W in S:

M_{p_{0}} = U − V W^{−1} Z.

Note that the blocks U, V, W, and Z have Toeplitz structure, U and W are triangular, and if d_{0} ≤ d_{1} (resp. d_{0} ≥ d_{1}), then so is Z (resp. V) also. Thus, we have the following algorithm:

Algorithm 3.2.1 (Multiplication by a polynomial modulo a polynomial, in the univariate case). Given three polynomials p_{0}, p_{1} and a of degrees d_{0}, d_{1} and less than d_{1}, respectively, compute the coefficient vector of the polynomial ap_{0} mod p_{1} as the matrix-by-vector product:

M_{p_{0}}a = (U − V W^{−1}Z)a,

where a is the coefficient vector of the polynomial a and U, V, W and Z are the above blocks of the Sylvester matrix S of (29).

The computation reduces to multiplication of the Toeplitz matrices Z of size d_{0} × d_{1} and U of size d_{1} × d_{1} by the vector a, solving the triangular Toeplitz system

W q = Za

of d_{0} equations, multiplying the Toeplitz matrix V by the solution q of this system, and subtracting the vectors V q from Ua.

With application of the algorithms of section B.1, one may perform algorithm 3.2.1 in O(d log d) ops, where d = max(d_{0}, d_{1}).

We are going to extend this approach to the multivariate case. Let us mention some of the main difficulties that are peculiar to the multivariate case but do not occur in the univariate case:

• We lose the notion of the leading monomial of the highest degree.

• We have no natural monomial basis for representing the quotient modulo a set of polynomials.

• When we homogenize the polynomials, we may introduce spurious solutions (at

infinity) to a polynomial system of equations. For the latter reasons and many others, we need to restrict our study to the cases where we may describe easily the structure of the matrices. These are the generic cases of two types that we are going to describe.

3.2.2 The generic multivariate case

In order to generalize the Sylvester matrix construction to the multivariate case, we consider n + 1 polynomials p_{0}, …, p_{n} and n + 1 vector spaces V_{0}, …, V_{n} generated by the monomials x^{F_{i}} = {x^{α}, α ∈ F_{i}}, where F_{i} is the set of the exponents,

F_{i} = {β_{i,1}, β_{i,2}, …}.

Let V be a vector space containing all the monomials of the polynomials p_{i} x^{β_{i}}, for β_{i} ∈ F_{i}, so that we can define the following map:

S : V_{0} × ⋯ × V_{n} → V  (30)
(q_{0}, …, q_{n}) ↦ ∑_{i=0}^{n} p_{i} q_{i}.

Let the matrix of the map S in the monomial basis of V_{0} × ⋯ × V_{n} and V be also denoted by S and take the following form:

V \{ \begin{matrix}x^{α_{1}} \\ · \\ · \\ · \\ x^{α_{N}} \end{matrix}  [ \begin{matrix}· & & · & & & · \\ · & & · & & & · \\ x^{β_{0,1}}p_{0} & ⋯ & x^{β_{1,1}}p_{1} & ⋯ & · · · · · · · · · & x^{β_{n,1}}p_{n} & ⋯ \\ · & & · & & & · \\ · & & · & & & · \end{matrix} ] .^{^{V_{0}}  ^{V_{1}}  ^{V_{n}}} (31)

Let us decompose such a matrix S into blocks S = [S_{0}, …, S_{n}], where S_{i} involves only the coefficients of p_{i}. The matrix S_{i} is a submatrix of the matrix of multiplication by p_{i}, defined in section 2.5. More precisely, S_{i} is the matrix of the map

π_{F} ∘ M_{p_{i}} ∘ π_{F_{i}}.

Thus, it is a quasi-Toeplitz matrix.

Algorithm 3.2.2 ( Multiplication of the matrix S of (31) by a vector). For every j, compute the products x^{β_{i,j}}p_{i} q_{i} for all i and sum them together in i. Output the sum for every j.

The complexity of this algorithm is bounded by C_{P olMult}(F_{0}, F) + ⋯ + C_{P olMult}(F_{n}, F) (due to using the algorithms of section B.1).

It is possible to consider the global matrix S as a quasi-Toeplitz matrix by adding a new variable x_{0}. The sum ∑_{i=0}^{n} p_{i} q_{i} can be computed from the product of p = ∑_{i} p_{i} x_{0}^{i} by ∑_{i=0}^{n} q_{i} x_{0}^{n−i}. Indeed, this sum is the coefficient of x_{0}^{n} in the product. Let F′ and F′′ be the sets of the exponents of the monomials in x_{0}^{n} x^{F} and ∪_{i=0}^{n}x_{0}^{n−i}x^{E_{i}}, respectively. Then the matrix S is the matrix of the operator

π_{F′} ∘ M_{p} ∘ π_{F′′}.

Remark 2 We can extend easily the construction of the map S to the case where the number of polynomials p_{0}, …, p_{m} is greater than n + 1 (m ≥ n).

Operators of this type have been extensively used in the literature, in order, for instance, to define resultants (see [20], [34], [15]). Let us recall that the vanishing of the resultant is the necessary and sufficient condition on the coefficients of the polynomials p_{0}, …, p_{n}, under which these polynomials have a common root (in a projective variety X). Two main examples appear in the literature:

• The classical case corresponds to X = ℙ^{n}, the projective space of dimension n. In this case, the polynomials p_{0}, …, p_{n} of degree d_{0}, …, d_{n} are homogenized, and the vanishing of the resultant gives a necessary and sufficient condition on their coefficients under which the homogenized polynomials have a common zero in ℙ^{n}. This case is referred to as Macaulay case (see [20]).

• In the second case, the variety X = T is a toric variety, and the map S is used to define the toric resultant of the polynomials p_{0}, …, p_{n}. The polynomials can also be homogenized in a toric sense, and the resultant gives a necessary and sufficient condition on their coefficients under which the toric-homogenized polynomials have a common zero in the toric variety T (see [15]). We refer to this case as the toric case.

Let us describe more carefully the monomials with exponents in F_{i} used in the construction of S.

The Macaulay case Let us fix integers d_{0}, …, d_{n}, and ν = d_{0} + ⋯ + d_{n} − n. For any d ∈ ℕ, let R_{d} denote the set of polynomials of degree not greater than d. Let p_{0}, …, p_{n} be polynomials of degrees d_{0}, …, d_{n}, respectively. To construct the map S that yields the resultant of these polynomials, we follow Macaulay’s work and choose V_{i} = R_{ν−d_{i}}, V = R_{ν}, so that we define the map

S : R_{ν−d_{0}} × ⋯ × R_{ν−d_{n}} → R_{ν}
(q_{0}, …, q_{n}) ↦ ∑_{i=0}^{n} p_{i} q_{i}.

The toric case In this case, we replace the constraints on the degree of the polynomials by the constraints on the support of the polynomials p_{i} (that is, the set of the exponents of the monomials with non-zero coefficients in p_{i}). Let C_{0}, …, C_{n} be polytopes in ℤ^{N} and let p_{0}, …, p_{n} ∈ L be Laurent’s polynomials, whose supports are in C_{0}, …, C_{n}, respectively. In order to construct the map S that yields the toric resultant, we fix (at random) a direction δ ∈ ℚ^{n}. For any polytope C, let C^{δ} denote the polytope obtained from C, by removing its facets whose normals have positive inner products with δ (see [5],[27]). For F_{i} = (∑_{j≠i} C_{j})^{δ} and F = (∑_{j} C_{j})^{δ}, we define the map

S : ⟨x^{F_{0}}⟩ × ⋯ × ⟨x^{F_{n}}⟩ → ⟨x^{F} ⟩
(q_{0}, …, q_{n}) ↦ ∑_{i=0}^{n} p_{i} q_{i}.

Many other examples of this type can be obtained by means of convenient choices of the vector spaces V_{0}, …, V_{n}, and V. We are going to examine the properties of these maps in the generic cases.

Definition 3.2.2 A property is generically true in the Macaulay case (or in the toric case), if this property is true for an (algebraically) open subset of the set of all possible values of the coefficients, according to the constraints on the degree (or on the support) of the polynomials.

Given polynomials p_{1}, …, p_{n}, we will compute from the matrix S:

• a basis (x^{α})_{α∈E} of the quotient A = R/(p_{1}, …, p_{n}),

• the table of multiplication by a polynomial p_{0} in A, from the matrix S (note that the matrix S of S is not a square matrix anymore),

• the dual basis of the monomial basis (x^{α})_{α∈E} of A.

These constructions will be valid for generic values of the coefficients of p_{1}, …, p_{n} but may fail for specific values. A more sophisticated method, described in [22], circumvents this difficulty by the compression of pencils of matrices.

3.2.3 A basis of A

First, we will define a subset E_{0} of exponents such that x^{E_{0}} is generically a basis of A = R/(p_{1}, …, p_{n}). For that purpose, we choose p_{0} = u_{0} + u_{1} x_{1} + ⋯ + u_{n} x_{n} (or p_{0} = u_{0} + u_{1} x_{1} + ⋯ + u_{n} x_{n} + u_{−1} x_{1}^{−1} + ⋯ u_{−n} x_{n}^{−1} in the toric case), where u_{i} are parameters. We also choose subsets E_{i} ⊂ F_{i} for i = 0, …, n, such that

(a) |E_{0}| + ⋯ + |E_{n}| = |F |

(b) and the matrix of the map

S̃: ⟨x^{E_{0}}⟩ × ⋯ × ⟨x^{E_{n}}⟩ → ⟨x^{F} ⟩
(q_{0}, …, q_{n}) ↦ ∑_{i=0}^{n} p_{i} q_{i}

takes the form

S̃ = \begin{matrix}\overbrace{E_{0}} & \overbrace{E_{1} ⋯ E_{n}} \\ \{ \begin{matrix}E_{0} \\ F′ \end{matrix} & [ \begin{matrix}U & V \\ Z & W \end{matrix} ] , \end{matrix}  (32)

where W is generically invertible.

In order to prove this generic property, it is sufficient to specify the coefficients of polynomials p_{i}, for which the property is satisfied.

Theorem 3.2.3 If condition (a) and (b) are satisfied, then for generic values of the coefficients of p_{1}, …, p_{n}, (x^{α})_{α∈E} is a generating set of A, and we have

dim_{ℂ}(A) ≤ |E_{0}|.

Proof. As W is generically invertible, the same process as in section 3.2.1 enables us to reduce modulo (p) the elements x^{α} p_{0} for α ∈ E_{0}, in ⟨x^{α}⟩_{α∈E_{0}}. As this is valid for any value of the parameter u_{i}, we can reduce in ⟨x^{α}⟩_{α∈E_{0}} modulo (p) the monomial x^{α} x_{i} (resp. x^{α} x_{i}^{−1} in the toric case), for any variable x_{i} and any α ∈ E_{0}. By induction, for any polynomial p in R (or for any Laurent polynomial in L in the toric case) and any α ∈ E_{0}, we can reduce modulo (p) the polynomial x^{α} p in ⟨x^{α}⟩_{α∈E_{0}}.

Therefore, as 1 ∈ (x^{α})_{α∈E_{0}} in the Macaulay case (or because any Laurent’s monomial p ∈ L is of the form p = p′x^{α} with α ∈ E_{0} and p′ ∈ L in the toric case), we can reduce modulo (p) any polynomial or Laurent polynomial p in ⟨x^{α}⟩_{α∈E_{0}} (where p ∈ R, in the Macaulay case, or p ∈ L, in the toric case). This proves that (x^{α})_{α∈E_{0}} is a generating set of A = R/(p_{1}, …, p_{n}) (A = L/(p_{1}, …, p_{n}) in the toric case). □

Let us show more details on how we choose the subsets E_{i} in the Macaulay case and in the toric case.

Macaulay case Let us choose E_{i} such that the matrix S̃ becomes the identity matrix (see [20]), when we replace the polynomial p_{i} by x_{i}^{d_{i}}. We can choose, for instance,

E_{0} = {(α_{1}, …, α_{n}); 0 ≤ α_{i} ≤ d_{i} − 1, i = 1, …, n},
E_{1} = {α = (α_{1}, …, α_{n}); |α| ≤ ν − d_{1}; 0 ≤ α_{i} ≤ d_{i} − 1, i = 2, …, n},
⋮
E_{n} = {α = (α_{1}, …, α_{n}); |α| ≤ ν − d_{n}},

where |α| = |α_{1}| + ⋯ + |α_{n}|.

According to theorem 3.2.3, (x^{α})_{α∈E_{0}} is generically a generating set of A, and

dim_{ℂ}(A) ≤ |E_{0}| = ∏_{i=1}^{n} d_{i},

which gives Bezout theorem.

Toric case In the toric case, the polynomial p_{i} is replaced by p_{i}^{t} = ∑_{α} a_{i,α} t^{w_{α}}x^{α} (where t is a new variable and w_{α} ∈ ℚ_{+}). The subsets of the exponents E_{i} are chosen so that the corresponding matrix S(t) = (s_{i,j}(t)) satisfies

deg_{t}(s_{i,i}(t)) < deg_{t}(s_{i,j}(t)) for i ≠ j

(see [15],[5] for more details). The set E_{0} is the set of the exponents in the mixed cells of a regular triangulation of C_{1} ⊕ ⋯ ⊕ C_{n}, so that, by construction, |E_{0}| is the mixed volume of C_{1}, …, C_{n}. This yields Bernstein theorem (part 1) (see [2], [19]).

Part 2 of Bernstein theorem shows that generically the number of common zeros of the system p_{1} = … = p_{n} = 0 is at least |E_{0}|. Thus, we deduce that dim_{𝕂}(A) ≥ |E_{0}|, and we have the following theorem:

Theorem 3.2.4 For generic values of the coefficients of p_{1}, …, p_{n}, (x^{α})_{α∈E_{0}} is a basis of A, in both Macaulay and toric case.

Note that we gave simpler proofs than in the articles [13], [28].

3.2.4 Matrices of multiplication in A

In this section, we still let S denote the map (30), constructed with using the fixed polynomials p_{1}, …, p_{n} and vector spaces V_{1}, …, V_{n}, V and with various choices of polynomial p_{0} and vector space V_{0} = ⟨x^{α}⟩_{α∈E_{0}}. The set of monomials (x^{α})_{α∈E_{0}}, defined in the previous section, is a basis of A.

For any polynomial p_{0}, we can also construct the table of multiplication by p_{0}, starting with a submatrix of S. Namely, we choose any subset E′_{i} ⊂ F_{i}, i = 1, …, n, such that simultaneously (a′) |E′_{1}| + ⋯ + |E′_{n}| = |F | − |E_{0}|,

(b′) and the corresponding columns in the matrix of S are linearly independent.

Generically, this is always possible, which we can show by giving a specific example. Decomposing again the matrix of the map

S̃: ⟨x^{E_{0}}⟩ × ⟨x^{E′_{1}}⟩ × ⋯ × ⟨x^{E′_{n}}⟩ → ⟨x^{F} ⟩
(q_{0}, …, q_{n}) ↦ ∑_{i=0}^{n} p_{i} q_{i}

in the form (32), we obtain the following property:

Theorem 3.2.5 For generic values of the coefficients of p_{1}, …, p_{n}, the matrix of multiplication by p_{0} in A is given by

M_{p_{0}} = U − V W^{−1} Z.

Proof. First, we will show that W is invertible. Otherwise, there exists a vector v ≠ 0 in the kernel of W. Then we have

[ \begin{matrix}V \\ W \end{matrix} ] v = [ \begin{matrix}w \\ 0 \end{matrix} ] ,

and w is not 0, because the columns [ \begin{matrix}V \\ W \end{matrix} ] of the matrix S are linearly indepen-

dent (condition (b′)). This implies that there is a non-zero polynomial of the form w(x) = ∑_{i=1}^{n} p_{i}q_{i} in ⟨x^{α}⟩_{α∈E_{0}}, which contradicts the fact that (x^{α})_{α∈E_{0}} is a basis of A. Consequently, W is invertible and, by the same argument as in section 3.2.1, U − V W^{−1}Z is the matrix M_{p_{0}} of multiplication by p_{0} in the basis (x^{α})_{α} of A. □

Example (continued) Let p_{0} = x_{1}, x^{E_{0}} = (1, x_{1}, x_{2}, x_{1} x_{2}), x^{E_{1}} = x^{E_{2}} = (1, x_{1}, x_{2}), and

x^{F} = (1, x_{1}, x_{2}, x_{1}x_{2}, x_{1}^{2}, x_{2}^{2}, x_{1}^{3}, x_{2}^{3}, x_{1}^{2}x_{2}, x_{1}x_{2}^{2}).

Then we have

S̃ = [ \begin{matrix}0 & 0 & 0 & 0 & −1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & −1 & −1 & 0 & −8 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & −1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 2 & 0 & −1 & 0 & 0 & −8 \\ 0 & 1 & 0 & 0 & 1 & −1 & 0 & 1 & −8 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 & 2 & 1 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 2 & 0 & 1 & 0 \end{matrix} ] .

We may verify that

U − V W^{−1} Z = [ \begin{matrix}0 & 1 & 0 & −\frac{14}{5} \\ 1 & 1 & 0 & −\frac{12}{5} \\ 0 & 0 & 0 & \frac{1}{5} \\ 0 & −2 & 1 & \frac{29}{5} \end{matrix} ]

is the matrix of multiplication M_{x_{1}} of the example of section 2.8.

3.2.5 The dual basis

It is possible to construct the dual basis (σ_{α})_{α∈E_{0}} of (x^{α})_{α∈E}, from the matrix S. Let

σ_{α} = ∑_{β∈ℕ^{n}} σ_{α,β} ∂^{β}

be the f.p.s. representing σ_{α} in ℂ[[∂]]. Then we have the following property:

Proposition 3.2.6 The coefficients [σ_{α,β}]_{α∈E_{0},β∈F} of (∂^{β})_{β∈F} in the dual basis (σ_{α})_{α∈E_{0}} are given by the matrix

[ 𝕀_{D} −V W^{−1} ] .

Proof. Let [σ_{α}] = [σ_{α.β}]_{β∈F} denote the vector of the first coordinates of σ_{α} and let Σ denote the matrix Σ = [σ_{α,β}]_{α∈E_{0},β∈F}. As E_{0} ⊂ F, we represent this matrix as a 1 × 2 block matrix Σ = [Σ′|Σ′′], where Σ′ = [σ_{α,β}]_{α,β∈E_{0}} and Σ′′ = [σ_{α,β}]_{α∈E_{0},β∈F −E_{0}}. The linear forms σ_{α} vanish on the multiples of p_{1}, …, p_{n}, which implies that

[Σ′|Σ′′] [ \begin{matrix}V \\ W \end{matrix} ] = 0

or, equivalently,

Σ′ V + Σ′′ W = 0.  (33)

The set (σ_{α})_{α∈E_{0}} is the dual basis of (x^{α})_{α∈E_{0}}, which means that, for any α, β ∈ E_{0}, σ_{α}(x^{β}) = σ_{α,β} equals 1 if α = β and 0 otherwise. In other words, Σ′ = 𝕀_{D} is the identity matrix and, according to (33) we have

Σ′′ = −V W^{−1}.
□

Algorithm 3.2.3 (Computation of the normal form of a multivariate polynomial). For any polynomial p ∈ ⟨x^{β}⟩_{β∈F}, compute its normal form by multiplying the matrix [𝕀_{D}| − V W^{−1}] by the coordinate vector of p.

Proposition 3.2.7 Algorithm 3.2.3 can be performed by using C_{LinSolve}(W)+C_{P olMult}(E_{0}, F)+ D ops, where C_{LinSolve}(W) denotes the arithmetic complexity of solving a linear system of equations with the coefficient matrix W.

Proof. The normal form of a polynomial p = ∑_{β∈F} p_{β} x^{β} is by definition

∑_{α∈E_{0}} σ_{α}(p) x^{α}.

The coefficients σ_{α}(p) = ∑_{β∈F} σ_{α,β} ∂^{α}(p) = ∑_{β∈F} σ_{α,β} p_{β} are obtained by multiplication of Σ = [𝕀_{D}| − V W^{−1}] by the vector [p_{β}]_{β∈F}. □

Example (continued) Let us be given the matrix

[𝕀_{4}| − V W^{−1}] = [ \begin{matrix}1 & 0 & 0 & 0 & 1 & −1 & \frac{33}{5} & \frac{14}{5} & −\frac{14}{5} & \frac{7}{5} \\ 0 & 1 & 0 & 0 & 1 & 7 & \frac{34}{5} & \frac{12}{5} & −\frac{12}{5} & \frac{6}{5} \\ 0 & 0 & 1 & 0 & 0 & 0 & −\frac{2}{5} & −\frac{1}{5} & \frac{1}{5} & \frac{2}{5} \\ 0 & 0 & 0 & 1 & −2 & 2 & −\frac{68}{5} & \frac{11}{5} & \frac{29}{5} & −\frac{12}{5} \end{matrix} ] .

The normal form of x_{1}x_{2}^{2} is defined by the last column of this matrix:

Nf(x_{1}x_{2}^{2}) = \frac{7}{5} + \frac{6}{5} x_{1} + \frac{2}{5} x_{2} − \frac{12}{5} x_{1}x_{2},

as found in the example of section 3.1.1. The linear form σ_{x_{1}x_{2}} (in the last row of this matrix) turns into

σ_{x_{1}x_{2}} = ∂_{1}∂_{2} − 2 ∂_{1}^{2} + 2 ∂_{2}^{2} − \frac{68}{5} ∂_{1}^{3} + \frac{11}{5} ∂_{2}^{3} + \frac{29}{5} ∂_{1}^{2}∂_{2} − \frac{12}{5} ∂_{1}∂_{2}^{2} + · · · .

3.3 Iterative methods in A

Efficient iterative methods for solving the system p = 0 rely on fast multiplication in A. We will next present some examples where the structure of the matrices can be exploited.

3.3.1 Fast multiplication in A

For a large class of polynomial ideals, specified, for instance, in [24], we may effectively compute the residue. Hereafter, we assume that the residue τ is known, and we will use it for computing efficiently the product of two elements in A.

For any element f ∈ A, let [f] denote the coordinate vector of f in the basis (x^{α})_{α∈E}. Let us write w_{α}(x) = ∑_{β∈E} B_{α,β}^{1}x^{α} to denote the dual basis of (x^{α})_{α∈E} and B_{1} = (B_{α,β}^{1})_{α,β∈E} to denote the Bezoutian of 1.

We want to compute the product [f g] in A where

\begin{matrix}f := ∑_{α∈E} f_{α} x^{α}, \\ g := ∑_{α∈E} g_{α} x^{α}. \end{matrix}

We may first compute the polynomial f g and then reduce it to a linear combination of the elements of the monomial basis ⟨x^{α}⟩ in order to obtain [f g]. We may also proceed directly by using the projection formula:

f g = ∑_{α∈E} τ(f g x^{α})w_{α}
= ∑_{α∈E} f g ⋆ τ(x^{α}) w_{α}.

In this case, we have to compute the coefficients of the linear form f g ⋆ τ and then shift from the basis (w_{α})_{α∈E} to the monomial basis (x^{α})_{α∈E}. By using relations (23), we may also proceed in an equivalent way, based on the formula

[f g] = M_{g}[f] = H_{1}^{−1} H_{f⋆τ}[g] = B_{1} H_{f⋆τ}[g].  (34)

As we want to compute the coefficients f g ⋆ τ(x^{α}) = τ(f g x^{α}) for α ∈ E, we need to know the value of τ for the monomials x^{α+β+γ} for α, β, γ ∈ E. Let τ̃ := ∑_{u∈3 E} τ(x^{u}) ∂^{u} denote the leading part of the series τ associated with the residue τ. We first compute

g ⋆ τ̃ = π_{+}(g(∂^{−1}) τ̃(∂))
= π_{+}((∑_{α∈E} g_{α} ∂^{−α})( ∑_{u∈3 E} τ_{u}∂^{u}))

and then

f g ⋆ τ̃ = f ⋆ (g ⋆ τ̃) = π_{+}(f(∂^{−1}) g(∂^{−1}) τ̃(∂))
= π_{+} ((∑_{α∈E} f_{α}∂^{−α}) (∑_{β∈E} g_{β}∂^{−β}) ( ∑_{u∈3 E} τ_{u} ∂^{u})) .

The coefficients λ_{α} of ∂^{α} in f g ⋆ τ̃ for α ∈ E are precisely the coefficients of f g in the dual basis (w_{α})_{α∈E} of A. Summarizing, we obtain the following algorithm:

Algorithm 3.3.1 (Multiplication by a polynomial modulo an ideal in a monomial basis). To obtain the coefficients of f g in the basis (x^{α}):

• Compute the coefficient vector Λ = [λ_{α}]_{α∈E} of ∂^{α} for α ∈ E, by multiplying the Laurent polynomial f(∂^{−1})g(∂^{−1}) by τ̃(∂).

• Multiply the vector Λ = [λ_{α}]_{α∈E} by the matrix B_{1} = H_{1}^{−1}, that is, solve the linear system of equations H_{1} v = Λ. Output the vector v.

3.3.2 Fast inversion in A

Similar techniques can be used to compute the inverse (reciprocal) of an invertible element f ∈ A. By relation (34), for g = f^{−1}, we have

[1] = H_{1}^{−1} H_{f⋆τ}[f^{−1}] or, equivalently, H_{1}[1] = H_{f⋆τ}[f^{−1}],

where H_{1}[1] is the coefficient vector of (∂^{α})_{α∈E} in τ̃. This yields the following algorithm:

Algorithm 3.3.2 (Inverse of a polynomial modulo an ideal in a monomial basis). To obtain the coefficients of f^{−1} in the basis (x^{α}):

• Compute the coefficient vector u = [λ_{α}]_{α∈E} of ∂^{α} for α ∈ E, in τ̃.

• Compute the coefficients of ∂^{α+β} for α, β ∈ E in the Laurent polynomial f(∂^{−1}) τ̃(∂), to dedine the matrix H_{f⋆τ}.

• Solve the linear system H_{f⋆τ} v = u. Output the vector v.

3.3.3 Computing selected roots of a polynomial system

As before, let Z denote the set of all common roots of the system p = 0. Then again, we assume here that the roots are simple.

By decomposing any element h of A in the basis of idempotents e_{ζ} (see section 2.12), we obtain that

h(x) = ∑_{ζ∈Z} h(x) e_{ζ} ≡ ∑_{ζ∈Z} h(ζ) e_{ζ}.

The second equation follows since e_{ζ}h(x) ≡ e_{ζ}h(ζ). Squaring h in the quotient ring A gives us that

h^{2} ≡ ∑_{ζ∈Z} h(ζ)^{2} e_{ζ}.

Here and hereafter, for any element a ∈ A, [a] denotes the vector of the coefficients of a in the basis (x^{α})_{α∈E}. In particular, [1] = (1, 0, ⋯, 0)^{t} if the basis starts with the monomial 1. Let || · || denote a norm in ℂ^{D}, say, the Euclidean or Hermitian norm,

||v|| = (v, v) = (∑_{i=1}^{D} |v_{i}|^{2})^{1/2}, v = (v_{i}), i = 1, …, D.

By minor abuse of notation, for any element a ∈ A, we will let ||a|| denote ||[a]||. Let h ∈ R and assume that there is a unique root ζ ∈ Z, for which the absolute value of h(ζ) is maximum, so that

|h(ζ)/h(η)| − 1 ≥ ρ,  (35)

for some fixed positive ρ and for any η ∈ Z distinct from ζ. (Since all the roots in Z are assumed to be distinct, we may, in principle, ensure the latter relation with a high probability, by means of a random linear substitution of the vector of the variables x. For the uniform distribution of all roots in a box in ℂ^{n}, we easily yield 1/ρ = O(D).) Then, by iteratively computing and normalizing the squares, we obtain

h_{0} = h, h_{i+1} ≡ h_{i}^{2}/||h_{i}^{2}||, i = 0, 1, …, k − 1,

and arrive at the bounds

ϵ_{k} := ∥\frac{h_{k}}{||h_{k}||} − \frac{e_{ζ}}{||e_{ζ}||}∥ ≤ \frac{c}{(1 + ρ)^{2^{k}}} .

Therefore, we ensure the bound ϵ_{k} ≤ 2^{−b} in k = k(ρ, b) = O(log(b/ρ)) recursive steps for any positive b, which shows that, right from the start, the process quadratically converges to a multiple of the idempotent e_{ζ}.

Proposition 3.3.1 In the case of simple roots and for h ∈ R such that |h(ζ)| > |h(η)| for any η ≠ ζ ∈ Z(I), the latter process of repeated squaring and normalization in A, always converges quadratically right from the start to a multiple of the idempotent e_{ζ}.

We refer the reader to [31] and [7] for some preceding works on a similar approach in the univariate case.

Algorithm 2.12.1 computes the root vector ζ as the product of H_{1} by the idempotent vector e_{ζ}. Therefore, the transition from e_{ζ} to the root ζ of the system p = 0 can be performed in C_{LinSolve}(H_{1}) ops.

Thus, we have the following algorithm:

Algorithm 3.3.3 (Computation of the simple root that maximizes the modulus of a fixed polynomial). Assume that all the roots in Z(I) are simple and h ∈ R is such that there exists ζ ∈ Z(I), with |h(ζ)| > |h(η)| for any η ≠ ζ ∈ Z(I).

• Set u_{0} := h and fix a positive tolerance value ϵ = 2^{−b}.

• Recursively, for k = 0, 1, ⋯, N −1, apply algorithm 3.3.1 to compute v_{k+1} ≡ u_{k}^{2} and u_{k+1} = \frac{v_{k}}{||v_{k}||} in A, until the norm ||u_{k+1} − u_{k}|| becomes smaller than ϵ,

• Compute and output the product of the matrix H_{1} by u_{N}.

The algorithm outputs a multiple of the vector [ζ^{α}]_{α∈E}, whose scaling gives us the root ζ for which |h(ζ)| is maximal (compare algorithm 3.1.1). The overall cost of the approximation of the root within an error norm 2^{−b} is bounded by O(D^{2} log(D) log(b/ρ)) ops.

Algorithm 3.3.3 can be applied recursively to approximate several roots of the polynomial system p = 0. In each application, the basic element h(x) of A should be updated; we propose the following policy:

h(x) = ∏_{ζ} h_{ζ}(x)

where the product is over all previously computed approximations ζ to the roots of the system p = 0 and h_{ζ}(ζ) = 0 for all ζ.

One may choose h_{ζ}(x) to be linear functions in x, h_{ζ}(x) = ∑_{i=1}^{n} h_{i,ζ}x_{i} − h_{0,ζ} with fixed or random coefficients h_{1,ζ}, …, h_{n,ζ} and with h_{0,ζ} = ∑_{i=1}^{n} h_{i,ζ}ζ_{i}. Under this choice, assumption (35) is satisfied for all computed approximations ζ to the roots, in the case of generic input, and we yield the following result.

Theorem 3.3.2 In δ recursive applications of algorithm 3.3.1, δ distinct roots of the generic polynomial system p = 0 can be approximated at the cost of performing O(aδD^{2} log(D)) ops, where a is the average number of iterations per root, a = O(log(\frac{b}{ρ})) for the approximation within the error 2^{−b} of roots ζ satisfying (35). In particular, a = O(log D) if \frac{b}{ρ} = D^{O(1)}.

A distinct extension to approximation of several roots is shown in [25].

3.3.4 Computing the closest root

Suppose that we seek a root of the system p = 0 whose coordinate x_{1} is the closest to a given value u ∈ ℂ. Let us assume that u is not a projection of any root of the system p = 0, so that x_{1} − u has reciprocal in A. Let ρ_{1}(x) denote such a reciprocal. We have ρ_{1} (x)(x_{1} − u) ≡ 1 and ρ_{1}(ζ) = \frac{1}{ζ_{1}−u}. Therefore, a root whose coordinate x_{1} is the closest to u_{1} is a root for which |ρ_{1}(ζ)| is the largest. Consequently, iterative squaring of ρ_{1} = ρ_{1}(x) shall converge to this root.

The polynomial ρ_{1} can be computed in the following way. Let \overline{M}_{x_{1}−u} denote the multiplication by x_{1} −u in A. Then ρ_{1} = (\overline{M}_{x_{1}−u})^{−1}(1), and according to the matrix equation (24), we have

[ρ_{1}] = H_{1} (H_{x_{1}} − u H_{1})^{−1}[1].

[ρ_{1}] defined by the latter equation can be computed within C_{LinSolve}(H_{x_{1}−u})+C_{P olMult}(−2E, E) ops, by using the black box algorithms of the appendix.

One may compute several roots of the polynomial system by applying the latter computation (successively or concurrently) to several initial values u.

Example (continued) We illustrate this approach by first computing the root for which x_{1} is maximal. We start with u_{0} = x_{1}. After 4 iterations, we obtain

u_{4} = 7.6055995 + 7.7975926x_{1} − 0.46159096x_{2} − 15.740471x_{1}x_{2}.

By multiplying the coefficient vector of this polynomial by H_{1} and dividing by the first coordinate, we obtain

[ζ_{1}^{α}]_{α∈E} = [1., 6.820095, −2.836734, −19.34680],

where ζ_{1} = (6.820095, −2.836734).

If we start with

u_{0} ≡ (x_{1} − \frac{1}{2})^{−1} ≡ −\frac{78}{35} − \frac{228}{35} x_{1} − \frac{32}{35} x_{2} − \frac{16}{7} x_{1}x_{2},

the algorithm should converge to the root with the coordinate x_{1} closest to \frac{1}{2}. Indeed, after 4 iterations, we obtain

u_{4} = 0.15292071 + 0.89409187x_{1} + 0.16270766x_{2} + 0.29923055x_{1}x_{2},

and after multiplication by H_{1} and normalization, we arrive at

[ζ_{4}^{α}]_{α∈E} = [1., 0.3678148, 1.675476, 0.6162664],

where ζ_{4} = (0.3678148, 1.675476) is the root with the coordinate x_{1} closest to \frac{1}{2}.

3.3.5 Splitting the set of roots

In addition to the repeated squaring iteration of algorithm 3.3.3, we will also consider iteration associated to a slight modification of the so-called Joukovski map (see [17],[7]): z ↦ \frac{1}{2}(z + \frac{1}{z}) and its variant z ↦ \frac{1}{2}(z − \frac{1}{z}).

The two attractive fixed points of this map are 1 and −1; for its variant, they turn into i and −i, i denoting here \sqrt{−1}.

Algorithm 3.3.4 (Sign iteration). Fix u_{0} = h ∈ ⟨x^{α}⟩_{α∈E}. Then recusively compute the polynomials u_{n+1} ≡ \frac{1}{2}(u_{n} − \frac{1}{u_{n}}) ∈ A, n = 0, 1, ….

Each iteration step of the algorithm can be performed within C_{LinSolve}(H_{u_{n}}) + C_{P olMult}(−3E, E) ops, by using algorithm 3.3.2 and the black box algorithms of the appendix.

Hereafter, ℜ(h) and ℑ(h) denote the real and imaginary parts of a complex h, respectively.

Proposition 3.3.3 Assume that for any root ζ ∈ Z, ℜ(h(ζ)) ≠ 0. Then the sequence (u_{n}) converges quadratically to σ = ∑_{ℑ(h(ζ))>0} e_{ζ} − ∑_{ℑ(h(ζ))<0} e_{ζ}, and we have

∥u_{n} − σ∥ ≤ Kρ^{2^{n}}

(for some constant K), where

ρ^{+} = max_{ℑ(h(ζ))>0,ζ∈Z(I)} |\frac{h(ζ) − i}{h(ζ) + i}| ,

ρ^{−} = max_{ℑ(h(ζ))<0,ζ∈Z(I)} |\frac{h(ζ) + i}{h(ζ) − i}| ,

and ρ = max{ρ^{+}, ρ^{−}}.

Proof. We apply the classical convergence analysis of the Joukovski map (see [17]) to the matrices of multiplication by u_{n} in A, whose eigenvalues are {u_{n}(ζ), ζ ∈ Z(I)}. □ This iteration can be applied to count the number of roots in a half-space, based on the following proposition:

Proposition 3.3.4 The rank of the matrix H_{σ^{+}} is the number of roots such that ℑ(h(ζ)) > 0 (where the roots are counted with their multiplicity).

Proof. As H_{1} is invertible, the rank of H_{σ^{+}} = H_{1}M_{σ^{+}} is the rank of M_{σ^{+}}, that is, the dimension of σ^{+}A, which equals ∑_{ℑ(h(ζ))>0} e_{ζ}A. The dimension of A_{ζ} = e_{ζ}A is the mutiplicity of ζ, which yields the proposition. □ The ranks can be computed efficiently by the algorithm supporting theorem B.3.2 (of appendix B).

Note that we may replace the iteration of algorithm 3.3.4 by the iteration u_{n+1} ≡ \frac{1}{2}(u_{n} + \frac{1}{u_{n}}) ∈ A and then replace ℑ by ℜ in the previous two propositions.

By 4n − 4m successive applications of the above splitting procedures, we can compute efficiently the number of roots in a fixed box in ℂ^{n} bounded by 2n − 2m pairs of parallel hyperplanes, 0 ≤ m < n. See [25] for more advanced applications of these techniques, which in particular enable us to restrict the approximation of the roots by algorithm 3.3.3 to a fixed box in ℂ^{n}.

3.4 Traces and real roots

In this section, we will keep assuming that the residue τ is known, will suppose that the coefficients of the polynomials p_{i} are real, and will study the problem of computing the numbers of distinct roots and of real roots. We will next define a special element of R̂, called the trace.

Definition 3.4.1 The linear form Tr is defined over any fixed field 𝕂 by

Tr : R → 𝕂
p ↦ trace(\overline{M}_{p}),

where trace(\overline{M}_{p}) is the usual trace of the linear operator \overline{M}_{p}.

By using this linear form, we define the quasi-Hankel matrix

H_{Tr} = [Tr(x^{α+β})]_{α,β∈E}.

In order to compute H_{Tr}, assuming that we know the table of the multiplications by x_{i} in A (i = 1, …, n), we may compute the values of x^{γ} (for γ = α +β and α, β ∈ E) by induction, for we have x^{γ} = x_{i} x^{γ′} with |γ′| < |γ| and Tr(1) = D = dim_{ℝ}(A). By using the linearity of the trace, we compute all the coefficients of H_{Tr} (see, for instance, [30]). Alternatively, we may apply the following theorem (see [11]):

Theorem 3.4.2 Let J ∈ R be the Jacobian of the polynomials p_{1}, …, p_{n}. Then

Tr = J ⋆ τ.

Example (continued) According to the example of section 2.8, we have

Tr(x_{1}) = 1 + \frac{29}{5} = \frac{34}{5}

and also

τ(x_{1} J) = τ(−16 − 16 x_{1} + 4 x_{2} + 34 x_{1}x_{2}) = \frac{34}{5}.

Algorithm 3.4.1 (Application of the trace to a monomial set). Compute and output H_{Tr} = [Tr(x^{α+β})]_{α,β∈E} as the product of the linear form

τ̃ = ∑_{α∈3 E} τ_{α} ∂^{α}

by J(∂^{−1}).

The number of ops involved in this algorithm is bounded by C_{P olMult}(3 E, −E). Once the matrix H_{Tr} is computed, we apply the following theorem, due to Hermite

(see [18], [26], [10]):

Proposition 3.4.3 (Hermite). Let J be the Jacobian of p = (p_{1}, …, p_{n}) and let B_{J} be the Bezoutian of J. Then

• the rank of H_{Tr} or B_{J} is the number of distinct roots,

• the signature of H_{Tr} or B_{J} is the number of real roots.

Algorithm 3.4.2 (Computation of the numbers of distinct roots and real roots). For a polynomial system p_{1} = ⋯ = p_{n} = 0 , define the matrix H_{Tr}, then compute the numbers of the distinct roots and the real roots of the system by applying proposition 3.4.3 and the algorithm supporting theorem B.3.2 (of appendix B).

Example (continued) The normal form of the Jacobian J is

J = −8 + 40 x_{1} − 2 x_{2} + 20 x_{1}x_{2}.

Note that τ(J) = \frac{1}{5} × 20 = 4 is the dimension of A. The matrix H_{Tr} is given by

H_{Tr} = [ \begin{matrix}4 & \frac{34}{5} & −\frac{12}{5} & −\frac{448}{25} \\ \frac{34}{5} & \frac{1166}{25} & −\frac{448}{25} & −\frac{16492}{125} \\ −\frac{12}{5} & −\frac{448}{25} & \frac{194}{25} & \frac{6976}{125} \\ −\frac{448}{25} & −\frac{16492}{125} & \frac{6976}{125} & \frac{234354}{625} \end{matrix} ] .

The matrix of the Bezoutian B_{J} is

B_{J} = [ \begin{matrix}−4 & −50 & 52 & −40 \\ −50 & 602 & −36 & 200 \\ 52 & −36 & 6 & −10 \\ −40 & 200 & −10 & 100 \end{matrix} ] .

The rank and the signature of both matrices are 4 and 2, respectively. The number of distinct roots is 4, and the number of distinct real roots is 2.

The overall randomized cost of computing the numbers of distinct roots and of real roots is O(D^{2}) up to a polylogarithmic factor.

4 Conclusion

We have deduced the results of sections 3.3 and 3.4 assuming that the residue τ associated with the ideal I = (p_{1}, ⋯, p_{n}) is known (or readily available). This somewhat restricts the class of polynomial systems to which the results apply. A major research challenge is an extension of these results to a more general class of polynomial systems of equations having a finite number of solutions. Substantial progress in this direction based on further extension of the techniques of this paper combined with some other new techniques has been reported in [25]. In [4] some further elaboration of the presented approach towards some practical problems of multivariate polynomial rootfinding and optimization was shown. Another research challenge (still widely open) is to extend the results of section 3.3 to approximating all the D roots of the system at the cost O(D^{2}) (up to a polylogarithmic factor).

Our goal, throughout this paper, was to demonstrate the power of the application of dual spaces, algebraic residues and the generalization of the structure of Toeplitz and Hankel matrices to the solution of a polynomial system in the multivariate case. In order to be able to yield the latter generalization, we reinterpreted the associated operators in terms of operations in the polynomial ring and in its dual. Multivariate Bezoutians and residues come naturally under these studies, and the algebraic interpretation of the associated operators yielded the relations between these matrices. Section 3 was devoted to applications of the developed techniques to yield order of magnitude improvement of the known algorithms for some fundamental problems of multivariate polynomial rootfinding. We hope that this work and our subsequent progress will motivate new interest in this recently open and challenging area.

References

[1] W. Auzinger and H. Stetter, An elimination algorithm for the computation of all zeros of a system of multivariate polynomial equations, in Proc Intern. Conf. on Numerical Math., vol. 86 of Int. Series of Numerical Math., Birkhäuser, 1988, pp. 11–30.

[2] D. Bernstein, The number of roots of a system of equations, Funct. Anal. and Appl., 9 (1975), pp. 183–185.

[3] D. Bini and V.Y. Pan, Polynomial and Matrix Computations, Vol. 1: Fundamental Algorithms, Birkhäuser, Boston, 1994.

[4] D. Bondifalat, B. Mourrain, and V. Y. Pan, Controlled iterative methods for solving polynomial systems, Proc. Annual ACM-SIGSAM Intern. Symp. on Symb. and Algebr. Comp. (ISSAC’98), ACM Press, New York, 1998, pp. 252-259.

[5] J. Canny and I. Emiris, An efficient algorithm for the sparse mixed resultant, in Proc. Intern. Symp. Applied Algebra, Algebraic Algor. and Error-Corr. Codes (Puerto Rico), G. Cohen, T. Mora, and O. Moreno, eds., vol. 263 of Lect. Notes in Comp. Science, Springer Verlag, 1993, pp. 89–104.

[6] J. Canny, E. Kaltofen, and Y. Lakshman, Solving systems of non-linear polynomial equations faster, in Proc. of the Annual ACM-SIGSAM Int. Symp. on Symb. and Alg. Comp. (ISSAC’89), ACM Press, New York, 1989, pp. 121–128.

[7] J. Cardinal, On two iterative methods for approximating the roots of a polynomial, in Proc. AMS-SIAM Summer Seminar on Math. of Numerical Analysis, (Park City, Utah, 1995), J. Renegar, M. Shub, and S. Smale, eds., vol. 32 of Lectures in Applied Math., American Math. Soc. Press, Providence, Rhode Island, 1996, pp. 165–188.

[8] J. Cardinal and B. Mourrain, Algebraic approach of residues and applications, in Proc. AMS-SIAM Summer Seminar on Math. of Numerical Analysis, (Park City, Utah, 1995), J. Renegar, M. Shub, and S. Smale, eds., vol. 32 of Lectures in Applied Math., American Math. Soc. Press, Providence, Rhode Island, 1996, pp. 189–210.

[9] D. Cox, J. Little, and D. O’Shea, Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra, Undergraduate Texts in Mathematics, Springer, 1992.

[10] M. Demazure, Charles Hermite, déjà …, Notes informelles de calcul formel 8, Centre de Math., Ecole Polytechnique (France), 1987.

[11] M. Elkadi and B. Mourrain, Approche effective des résidues algébriques, Rapport de Recherche 2884, INRIA, Sophia Antipolis, 1996.

[12] I. Z. Emiris and V. Y. Pan, The Structure of Sparse Resultant Matrices, in Proc. Annual ACM-SIGSAM Intern. Symp. on Symbolic and Algebraic Computation (ISSAC’97), ACM Press, New York, 1997, pp. 189-196.

[13] I. Emiris and A. Rege, Monomial bases and polynomial system solving, in Proc. Annual ACM-SIGSAM Intern. Symp. on Symbolic and Algebraic Computation (ISSAC’94), ACM Press, New York, 1994, pp. 114–122.

[14] P. Fuhrmann, A Polynomial Approach to Linear Algebra, Springer-Verlag, 1996.

[15] I. Gelfand, M. Kapranov, and A. Zelevinsky, Discriminants, Resultants and Multidimensional Determinants, Birkhäuser, Boston-Basel-Berlin, 1994.

[16] G. H. Golub and C. F. Van Loan, Matrix Computations (third edition), John Hopkins Univ. Press, Baltimore, Maryland, 1996.

[17] P. Henrici. Applied and Computational Complex Analysis, volume I. Wiley, 1988.

[18] C. Hermite, Remarques sur le théorème de Sturm, C. R. Acad. Sci. de Paris, 36 (1853), pp. 52–54.

[19] A. Kushnirenko, A Newton polyhedron and the number of solutions of a system of k equations in k unknowns, Uspekhi Matem. Nauk, 30 (1975), pp. 266–267.

[20] F. Macaulay, Some formulae in elimination, Proc. London Math. Soc., 1 (1902), pp. 3–27.

[21] F. Macaulay, The Algebraic Theory of Modular Systems, Cambridge Univ. Press, 1916.

[22] B. Mourrain, Computing isolated polynomial roots by matrix methods, J. Symb. Comp., December, 1998.

[23] B. Mourrain and V. Y. Pan, Multidimensional structured matrices and polynomial systems, Calcolo, 33 (1996), pp. 389-401.

[24] , Solving special polynomial systems by using structured matrices and algebraic residues, in Proc. of the workshop on Foundations of Computational Mathematics (Rio de Janeiro, 1997), F. Cucker and M. Shub, eds., Springer, 1997, pp. 287–304.

[25] B. Mourrain and V. Y. Pan, Asymptotic acceleration of solving multivariate polynomial systems of equations, in Proc. 30th Ann. ACM Symp. on Theory of Computing, ACM Press, New York, 1998, pp. 488-496.

[26] P. S. Pedersen, M.-F. Roy, and A. Szpirglas, Counting real zeros in the multivariate case, in Effective Methods in Algebraic Geometry (MEGA’92), A. Galligo and F. Eyssette, eds., Progress in Math., Nice (France), Birkhäuser, 1993, pp. 203–223.

[27] P. S. Pedersen and B. Sturmfels, Product formulas for resultants and Chow forms, Math. Zeitschrift, 214 (1993), pp. 377–396.

[28] , Mixed monomial basis, in Effective Methods in Algebraic Geometry (MEGA’94), vol. 143 of Progress in Math., Santander (Spain), Birkhäuser, 1994, pp. 285–306.

[29] P. Penfield Jr., R. Spencer, and S. Duinker, Tellegen’s Theorem and Electrical Networks, M.I.T. Press, Cambridge, Massachussetts, 1970.

[30] F. Rouillier, Systèmes polynomiaux, PhD thesis, Univ. de Rennes, 1996.

[31] J. Sebastiao e Silva, Sur une méthode d’approximation semblable à celle de Graeffe, Portugal Math. J., 2 (1941), pp. 271–279.

[32] H. J. Stetter, Eigenproblems are at the heart of polynomial system solving, SIGSAM Bulletin, 30, 4, ACM Press, New York (1996), pp. 22–25.

[33] E. E. Tyrtyshnikov, A unifed approach to old and new theorems on distribution and clustering, Linear Algebra and Its Applications, 232 (1998) pp. 1-43.

[34] B. Van der Waerden, Modern Algebra, Vol. II, Frederick Ungar Publishing Co, 1948.

A Polynomials, Laurent’s polynomials, and dual spaces (univariate case). Basic Definitions

Consider univariate polynomials p = p(x) = ∑_{i=0}^{d} p_{i}x^{i} ∈ R = ℂ[x], represented by vectors of their complex coefficients p = (p_{0}, …, p_{d}). Let the subspace R_{d} denote the vector space (of dimension d + 1) of polynomials in R of degree at most d.

A fixed polynomial p(x) of R generates the ideal I = (p(x)) in R, formed by all polynomial multiples q(x) of p(x). Let A = R/I denote the quotient ring of polynomials reduced modulo p(x) (that is, modulo the ideal I). If p(x) is of degree d, then A is isomorphic to R_{d−1}, as a vector space.

By introducing the reciprocal x^{−1}, we arrive at the ring of Laurent’s polynomials ℂ[x, x^{−1}] = L and let L_{−c,d} denote the subspace of Laurent’s polynomials of the form ∑_{i=−c}^{d} λ_{i} x^{i}.

A polynomial p ∈ R_{d} can be represented by the vector of its d + 1 coefficients or, equivalently, by the values p(0), p′(0), …, \frac{1}{d !}p^{(d)}(0). In other words, a primal basis of R_{d} is ⟨1, x, …, x^{d}⟩, and its dual basis (that is, the set of linear forms (maps) that compute the coefficients of p in the primal basis) is the set of linear forms

⟨p ↦ \frac{1}{i !}p^{(i)}(0)⟩_{i=0,...,d}.

We introduce a new variable ∂ and let ∂^{i} denote the i^{th} element, p ↦ \frac{1}{i !}p^{(i)}(0), of this dual basis. Thus, a linear form on R_{d}, that is, an element Λ of the dual space R̂_{d} of R_{d}, is represented by a polynomial

Λ = ∑_{i=0}^{d} λ_{i} ∂^{i}.

For any p ∈ R_{d}, we have Λ(p) = ∑_{i=0}^{d} λ_{i} \frac{1}{i !} \frac{d^{(i)}}{dx^{i}} (p)(0) and λ_{i} = Λ(x^{i}).

Next, consider linear forms Λ ∈ R̂ on the primal space R. The restrictions of the linear forms to R_{d} ⊂ R are the elements of R̂_{d}, which can be represented by polynomials in ∂ of degree at most d. This is valid for any d; therefore, an element Λ ∈ R̂ is a formal power series (f.p.s.) in ∂:

Λ = ∑_{i=0}^{∞} Λ(x^{i})∂^{i}.

Such a ring of f.p.s. in the variable ∂ is denoted S = ℂ[[∂]].

For any p(x) ∈ ℂ[x] and Λ(∂) ∈ ℂ[[∂]], we define

p(x) ⋆ Λ(∂) = π_{+}(p(∂^{−1})Λ(∂)),

where π_{+} : ℂ[∂^{−1}][[∂]] → ℂ[[∂]] is the projection on the monomials having nonnegative exponents in ∂.

Example

(1 + x^{2}) ⋆ (∂^{3} + 3 ∂ − 2) = ∂^{3} + 4 ∂ − 2.

Contrary to [14], we introduce a new variable ∂ for the “inverse” of x, which we consider as an element of the dual space.

B Algorithms and complexity for some linear algebra computations

We will recall the known estimates for the computational cost of performing some basic algorithms used in this paper.

B.1 Polynomial multiplication

In sections 1 and 2, we reduced multiplication of various structured matrices by vectors to polynomial multiplication. Now, let us recall the known arithmetic complexity bounds for the latter operation. As before, let C_{P olMult}(E, F) denote the number of arithmetic operations required for the multiplication of a polynomial with support in E by a polynomial with support in F.

Theorem B.1.1 ([3], p. 13) Let E_{d} = [0, …, d] ⊂ ℕ. Then

C_{P olMult}(E_{d}, E_{d}) = O(d log(d)).

Theorem B.1.2 ([3], p. 64) Let E_{d} = {(α_{1}, …, α_{n}) ; 0 ≤ α_{i} ≤ d_{i} − 1}. Then we have

C_{P olMult}(E_{d}, E_{d}) = O(M log(M)),

where d = max(d_{1}, …, d_{n}), M = c^{n}, and c = 2 d + 1.

Theorem B.1.3 [6], [12]. Let E_{d,n} be a fixed set of exponents having total degree at most d in n variables and let S be a given set containing the set E_{d,n} + E_{d,n}. Then

C_{P olMult}(E_{d,n}, E_{d,n}) = O(C_{P olMult}(E_{T} , E_{T}) log T),

where T = (_{n}^{n+d}) is the maximum possible number of distinct monomials of degree at most d in n variables. Furthermore, C_{P olMult}(E_{d,n}, E_{d,n}) = O(s(log^{2}(s) + n) + n\sqrt{d}) where s is the cardinality of the set S.

Remark Theorems B.1.1 and B.1.2 can be extended to the computations over any ring of constants (rather than over the complex field) at the expense of increasing their complexity bounds by factors at most log log(d) and log log(c), respectively. Theorem B.1.3 applies over any field of constants having characteristic 0.

B.2 Tellegen’s theorem on duality of multiplication of a matrix and its transpose by a vector

Theorem B.2.1 [29]. Let W be a square n×n matrix with no zero rows or columns. Let C_{W} ops suffice to compute the product W v for a vector v. Then C_{W} ops also suffice to compute the product W^{t}v = (v^{t}W)^{t}.

B.3 Solving a linear system of equations and the computation of the rank and the signature of a matrix

Application of the conjugate gradient algorithm ([16], sect. 10.2) gives us the following result:

Theorem B.3.1 Let W be a nonsingular N × N matrix. Performing 2N multiplications of W and W^{t} by vectors and O(N^{2}) other arithmetic operations suffice to compute the solution v to a linear system W v = w.

In section 3.4, we needed an algorithm for computing the rank and the signature of an N × N real symmetric (and quasi-Hankel) matrix W. Such an algorithm may start with tridiagonalizing the matrix. In exact arithmetic, this can be done by means of the Lanczos algorithm, whose randomized computational cost is dominated by the cost of performing O(N) multiplications of W by vectors and O(N^{2}) other ops (cf. [3], p.118, or [16], ch. 9]). Then both the rank and the signature of W can be computed at the cost O(N), by using the three-term recurrence relations for the characteristic polynomials of the leading principal submatrices of W (cf. [16], p. 439). We arrive at the following result.

Theorem B.3.2 Let W be an N-by-N real symmetric matrix. Then application of Lanczos randomized algorithm (which uses O(N) multiplications of W by vectors and O(N^{2}) other ops) and performing O(N) additional ops suffice to compute the rank and the signature of W.

If W is a structured (resp. and real symmetric) matrix, whose multiplication by a vector is expressed in terms of polynomial multiplication, one may combine theorems B.1.1-B.1.3 and B.3.2 in order to express the arithmetic cost of the solution of the linear system W v = w and the randomized arithmetic cost of computing the rank (resp. and signature) of W in terms of the size of W.