# Resultant-based methods for solving bivariate polynomial equations

H. Khalil

This section describes an algorithm, which uses generalized eigenvalues and eigenvectors, for solving a system of two bivariate polynomials , with . It treats multiple roots, which are usually considered as obstacles for numerical linear algebra techniques. The numerical difficulties are handled with the help of singular value decomposition (SVD) of the matrix of eigenvectors associated to each eigenvalue. The extension to equations could be treated similarly.

## Generalized eigenvalues and eigenvectors

Let be two square matrices of size . A {generalized eigenvalue} of and is a value in the set

A vector is called a {generalized eigenvector} associated to the eigenvalue if . The matrices and have generalized eigenvalues if and only if . If , then can be finite, empty, or infinite. Note that if then . Moreover, if is invertible then , which is the ordinary spectrum of .\

Recall that an matrix with polynomial coefficients can be equivalently written as a polynomial with matrix coefficients. If , we obtain , where are matrices.

 Definition: The companion matrices of are the two matrices and given by:

And we have the following interesting property :

 Proposition: With the above notation, the following equivalence holds:

## Intersection of planar curves

From now on we will assume that is an algebraically closed field and suppose given two bivariate polynomials and in . Our aim will be to study and compute their common roots. This problem can be interpreted more geometrically. Polynomials and actually define two algebraic curves in the affine plane (having coordinates ) and we would like to know their intersection points.

Hereafter we will consider systems having only a {finite} number of roots. This condition is not quite restrictive since it is sufficient (and necessary) to require that polynomials and are coprime in (otherwise one can divide them by their ).

In the case where one of the curve is a line, i.e. one of the polynomial has degree 1, the problem is reduced to solving a univariate polynomial. Indeed, one may assume e.g. that and thus we are looking for the roots of . Let us denote them by . Then we know that

where the 's are assumed to be distinct and is a non-zero constant in . The integer , for all , is called the {multiplicity} of the root , and it turns out that if do not vanish at infinity in the direction of the -axis (in other words, if the homogeneous part of of highest degree do not vanish when ). This later condition can be avoid in the projective setting: let be the homogeneous polynomial obtained by homogenizing with the new variable , then we have

where is an integer corresponding to the multiplicity of the root at infinity, and . Moreover, it turns out that the roots and their corresponding multiplicities can be computed by eigenvalues and eigenvectors computation.

In the sequel we will generalize this approach to the case where and are bivariate polynomials of arbitrary degree. For this we first need to recall the notion of multiplicity in this context. Then we will show how to recover the roots from multiplication maps.

### Intersection multiplicity

Let be two coprime polynomials in , and the corresponding algebraic plane curves, the ideal they generate in and the associated quotient ring. We denote by the distinct intersection points in of and (i.e. the distinct roots of the system ).

A modern definition of the intersection multiplicity of and at a point is (see [Ful84])

where the point is here abusively (but usually) identified with its corresponding prime ideal in , denoting the ordinary localization of the ring by this prime ideal. As a result, the finite -algebra (which is actually finite if and only if and are coprime in ) can be decomposed as the direct sum

and consequently

The intersection multiplicities can be computed using a resultant. The main idea is to project algebraically'' the intersection points on the -axis. To do this let us see both polynomials and in where , that is to say as univariate polynomials in with coefficients in the ring which is a domain; we can rewrite

Their resultant (with respect to ) is an element in which is non-zero, since and are assumed to be coprime in , and which can be factorized, assuming that and do not vanish simultaneously, as

where the 's are distinct elements in and (as sets). For instance, if all the are distinct then we have and for all . Moreover, we have (see e.g. [Ful84] (S)):

 Proposition: For all , the integer equals the sum of all the intersection multiplicity of and at the points such that .

As a corollary, if all the 's are distinct (this can be easily obtained by a linear change of coordinates ) then is nothing but the valuation of at . Another corollary of this result is the well-know Bezout theorem for algebraic plane curves, which is better stated in the projective context. Let us denote by and the homogeneous polynomials in obtained from and by homogenization with the new variable .

 Proposition: [Bezout theorem] If and are coprime then the algebraic projective plane curves associated to and intersect in points in , counted with multiplicities.

{proof} It follows straightforwardly from the previous proposition and the well-known fact that is a homogeneous polynomial in of degree (see e.g. [Lang02], [VWe48II]). {proof}

### Intersection points

We take again the notation of A[x]} and [*] . We denote by the Sylvester matrix of and seen in (assuming that both as degree at least one in ); thus . From [*] , we deduce immediately that vanishes at a point if and only if either there exists such that , or either (in which case the intersection point is at infinity). Therefore one may ask the following question : being given a point such that , how to recover all the possible point such that ?

Suppose given such a point and assume that and are not both zero. If is of dimension one, then it is easy to check that there is only one point such that and moreover any element in is a multiple of the vector . Therefore to compute , we only need to compute a basis of and compute the quotient of its second coordinate by its first one. In case of multiple roots, this construction can be generalized as follows:

 Proposition: Let be any basis of the kernel , be the matrix whose row is the vector , be the -submatrix of corresponding to the first columns and be the -submatrix of corresponding to the columns number . Then is the set of roots of the equations (i.e. the set of -coordinates of the intersection points of and above .

In this theorem, observe that the construction of and is always possible because we assumed that both polynomials and depend on the variable , which implies that the number of rows in the matrix is always strictly greater than the .

The previous theorem shows how to recover the -coordinates of the roots of the system from the Sylvester matrix. However, instead of the Sylvester matrix we can use the Bezout matrix to compute (see section [*] ). This matrix have all the required properties to replace the Sylvester matrix in all the previous results, except in proposition [*] . Indeed, the Bezout matrix being smaller than the Sylvester matrix, the condition given in remark [*] is not always fulfilled; using Bezout matrix in proposition [*] requires that for all roots of ,

This condition may failed: take for instance in the system

Note that the use of the Bezout matrix gives, in practice, a faster method because it is a smaller matrix than the Sylvester matrix, even if its computation takes more time (the save in the eigen-computations is greater).

## The algorithm

We are now ready to describe our resultant-based solver. According to proposition [*] , we replace the computation of the resultant, its zeroes, and the kernel of (or the Bezout matrix) for each root by the computation of generalized eigenvalues and eigenvectors of the associated companion matrices and . This computation can be achieved with the QZ algorithm [GLVL96]. Here are some numerical problems which might occur: {enumeratealpha}
• In order to compute the -coordinates of intersection points, it is necessary to compute the real generalized eigenvalues of . Numerically, some of these values can contain a nonzero imaginary part. Generally such a problem is encountered for multiple eigenvalues.
• What does it mean by a numerical multiple -coordinate?
• How do we choose the linearly independent vectors among the computed eigenvectors? {enumeratealpha} For difficulty a'' we choose , and any eigenvalue whose imaginary part is of absolute value less than will be considered real.
For difficulty b'', we choose possibly a different and gather the points which lay in an interval of size . According to the following proposition, the algorithm remains effective.

 Proposition: Assume that with corresponding generalized eigenvectors with . We define the matrix whose its rows are Then, if is the first left block of and the second block, then the generalized eigenvalues of and give the set of ordinates of the intersections above ; i.e. this algorithm can be used to gather clusters of -coordinates and regarding them as the projection of only one point.

For difficulties c'', it is necessary to use the singular value decomposition: If a generalized eigenvalue of and , and the matrix of the associated eigenvectors, then the singular value decomposition of is

and is the rank of . Let

The matrix is also a matrix of eigenvectors, because the changes are done on the columns, and V is invertible. Hence and have the same rank , and is a matrix of linearly independent eigenvectors. Consequently can be used to build and .

This leads to the following complete algorithm:

{Algorithm { {Numerical approximation of the root of a square bivariate system} {{-0.3cm}}

{Input:} Two bivariate polynomials

• {{Compute the Bezout matrix of and , and compute the matrices and ;}}
• {{Compute the generalized eigenvalues and eigenvectors of to get the coordinates of the roots;}}
• {{Eliminate the imaginary eigenvalues and the values at infinity, and gather the close points (with the input );}}
• For each of {{all the close points represented by , take the matrix of the associated eigenvectors. If their number then gives a basis of ;}}
• Compute a singular value decomposition of ;
• {{Compute the rank of by testing the until is found (the rank is thus determined to be );}}
• {{Let be the first columns of . Define as the first block in , and as the second block in ;}}
• {{Compute the generalized eigenvalues of to get the coordinates;}}
{Output:} an approximation of the real roots of , . }}{}

## Implementation

The function which implement this algorithm is
template<class C> Seq<VectDse<C> > solve(p,q,EigenRes<C>());

synaps/resultant/bivariate_res.h.

## Example

#include <synaps/msolve/projection/EigenRes.h>

#ifdef SYNAPS_WITH_LAPACK
using std::cout;
using std::cin;
using std::endl;

typedef MPol<double> mpol_t;

int main(int argc, char **argv)
{

mpol_t p1("16*x0^4-160*x0^2*x1^2+400*x1^4-32*x0^3+160*x0*x1^2+6*x0^2-50*x1^2+10*x0+1.5625"),
p2("x0^2-x0*x1+x1^2-1.2*x0-0.0625");
cout<<"p1="<<p1<<endl;
cout<<"p2="<<p2<<endl;

Seq<VectDse<double> > res;
Seq<int> mux,muy;
res= solve(p1,p2, EigenRes<double>(1e-3));
res= solve(p1,p2,mux, EigenRes<double>(1e-6));
res= solve(p1,p2,mux,muy, EigenRes<double>(1e-7));

cout<<"Intersection points:"<<endl;
for(unsigned i=0;i<res.size();++i)
std::cout<<"("<<res[i][0]<<","<<res[i][1]<<") of multiplicity "<<mux[i]<<" "<<muy[i]<<std::endl;
}
#else
#include <iostream>
int main()
{
std::cout <<"  LAPACK is not connected"<<std::endl;
return(77);
}
#endif //SYNAPS_WITH_LAPACK

SYNAPS DOCUMENTATION