Ioannis Z Emiris
email: emiris@di.uoa.gr

Christos Konaxis
email: ckonaxis@di.uoa.gr

Clément Laroche
email: claroche@di.uoa.gr

Implicit representations of high-codimension varieties

Keywords: Implicitization, space curve, Chow form, resultant, conical hypersurface, randomized algorithm, Maple implementation

Implicitization usually focuses on plane curves and (hyper)surfaces, in other words, varieties of codimension 1. In this paper we shift the focus on space curves and, more generally, on varieties of codimension larger than 1, and discuss approaches that are not sensitive to base points.

Our first contribution is a direct generalization of an implicitization method based on interpolation matrices for objects of high codimension given parametrically or as point clouds. Our result shows the completeness of this approach which, furthermore, reduces geometric operations and predicates to linear algebra computations.

Our second, and main contribution is an implicitization method of parametric space curves and varieties of codimension > 1, which exploits the theory of Chow forms to obtain the equations of conical (hyper)surfaces intersecting precisely at the given object. We design a new, practical, randomized algorithm that always produces correct output but possibly with a non-minimal number of surfaces. For space curves, which is the most common case, our algorithm returns 3 surfaces whose polynomials are of near-optimal degree; moreover, computation reduces to a Sylvester resultant. We illustrate our algorithm through a series of examples and compare our Maple code with other methods implemented in Maple. Our prototype is not faster but yields fewer equations and is more robust than Maple's implicitize. Although not optimized, it is comparable with Gröbner bases and matrix representations derived from syzygies, for degrees up to 6.

Introduction

In manipulating geometric objects, it is essential to possess robust algorithms for changing representation. This paper considers two fundamental representa-tions and, in particular, switching to implicit representation from the parametric representation, including the case of base points. One of our methods also treats point cloud representations as input. While many algorithms allow to switch representation of hypersurfaces of codimension 1, our algorithms construct implicit equations expressed as matrix determinants for parametric curves and (hyper)surfaces of codimension > 1. More precisely, for varieties of higher codimension, one cannot expect a single implicit equation. Instead, we compute a number of implicit (hyper)surfaces whose set-theoretic intersection equals the given variety. For instance, rational space curves are defined as the intersection of two or, at most, three implicit surfaces.

Implicit representations of geometric objects allow to answer to the membership and intersection problems. With such a representation, one can check easily whether given points are lying on the object or not. They are used for applications requiring geometric boolean operations such as intersections of several objects, differences and, in particular, use ray shooting and ray tracing methods. We introduce two implicitization methods for varieties of codimension higher than 1. Our first method is a direct generalization of interpolation matrices for objects given parametrically, but also as point clouds. Interpolation matrix methods usually are very simple but, since they can accept inputs represented as point clouds or can be used for approximate implicitization (as in [START_REF] Dokken | Approximate implicitization[END_REF]), they provide a very relevant, alternative representation especially for industrial applications today. In that context, interpolation matrices are sometimes called collocation matrices. The construction of the interpolation matrix representing the object reduces geometric operations to linear algebra.

The second implicitization method, and our main contribution, is the following. Given a variety in parametric form, and starting with the powerful and classic theory of Chow forms, which generalizes resultant theory, we design an original, practical algorithm that uses randomization in order to compute implicit hypersurfaces containing the variety. Computing the Chow form is expensive in the general case; high complexity is one reason why very few implementations exist.

Our algorithm cannot guarantee to yield the optimal implicit equations but is reasonably fast and usually its output is very close to being optimal. For space curves, our method yields 3 implicit surfaces defining the curve set-theoretically, which is the minimum number of equations in general. The degree of the defining polynomials is close to being optimal. The output given by our method possesses strong geometric structure: all computed implicit equations correspond to conical (hyper)surfaces. The question is eventually reduced to the standard question of computing a (sparse) resultant; it is the complexity bottleneck of the algorithm.

The method generalizes for any variety of codimension higher than 1, even in the presence of base points or self-intersections. In particular, self-intersections and base points do not affect the degree of the output equations. We produce a number of hypersurfaces each containing the variety set-theoretically and whose degree is bounded in terms of mixed volume. The number of such hypersurfaces that define the variety set-theoretically is unknown in the general case but experiments indicate that n + 1 are sufficient, where n is the ambient dimension. This complies with a result by Kronecker [START_REF] Kronecker | Grundzüge einer arithmetischen Theorie der algebraischen Grössen[END_REF] but uses one more hypersurface than the optimal number n [16, Ch.V]. The specific output of our algorithm makes it difficult to exploit this result.

We implement our algorithm in Maple and Sage, illustrate it through a series of examples and compare our code with other methods implemented in Maple. We provide a comparison between various existing methods and ours in Table 1, specifying advantages and disadvantages of each method. Our prototype is not faster but yields fewer equations and is more robust than Maple's native command implicitize. Although not optimized, for curves in any ambient dimension its speed is comparable to Gröbner bases and matrix representations derived from syzygies, for degrees up to 6. This paper extends certain results and ideas presented in preliminary form in the second part (Section 4) of [START_REF] Emiris | Matrix representations by means of interpolation[END_REF].

The rest of the paper is organized as follows: The next section sketches the most relevant existing work. Section 3 generalizes interpolation matrices to codimension > 1. Sections 4 and 5 employ the Chow form to handle highercodimension parametric varieties in 3 and higher dimension. Section 6 examines experimentally the performance of our methods and their implementation, and compares them to other methods of implicitization. We conclude with future work and open questions.

Preliminaries and previous work

This section introduces basic concepts and sketches the most relevant existing work.

Unless specified otherwise, we work in the projective or affine spaces of dimension n over the complex field. A homogeneous parameterization of a projective variety V of dimension d is given as

P d → P n t = (t 0 : • • • : t d ) → (F 0 (t) : • • • : F n (t))
where F j , j = 0, . . . , n is a homogeneous polynomial of some unique degree δ. A rational parameterization is obtained by dehomogenization of a homogeneous parameterization:

f j (t 1 , . . . , t d ) = F j (1, t 1 , . . . , t d ) F 0 (1, t 1 , . . . , t d ) , j = 1, . . . , n.
There is no difficulty in switching between rational and homogeneous parameterizations but the latter deals with the points at infinity. Most of the time, we work in the projective space with homogeneous parameterization. For simplicity, we assume that t 0 is not a common factor of the polynomials F j hence δ = max{deg(F j (t) |t0=1 ), j = 0, . . . , n}. We call δ the degree of the parameterization (both rational and homogeneous).

In general, approaches to implicitization include resultants, Gröbner bases, moving lines and surfaces (based on the theory of syzygies), and interpolation techniques. Resultants, and their matrix formulae, have led to practical methods to express the implicit surface equation, e.g., in [START_REF] Manocha | The implicit representation of rational parametric surfaces[END_REF][START_REF] Manocha | Algorithms for implicitizing rational parametric surfaces[END_REF], but require the assumption that no base points exist.

Here, the resultant (resp. homogeneous resultant), without any more precision, should be understood as in the following: Definition 1. The resultant is the abstract map R of d+1 polynomials h 0 , . . . , h d in d variables with symbolic coefficients (resp. homogeneous polynomials in d+1 variables) satisfying:

1. R is a homogeneous polynomial of degree i =i deg h i in the symbolic coefficients of h i and 2. R vanishes if and only if the polynomials h i have a common root (resp. a non-trivial common root).

Such a map is known to exist and is unique up to a non-zero constant factor (see e.g. [START_REF] Gelfand | Discriminants, Resultants and Multidimensional Determinants[END_REF]). Its computation may vary: if d = 1, it is the determinant of the Sylvester, Bézout or hybrid matrices. In general it can be computed as the quotient of the determinant of the Macaulay or of the sparse resultant matrix divided by the determinant of one of its minors (cf. [START_REF] Dickenstein | Explicit formulas for the multivariate resultant[END_REF]).

Implicit matrix representations are quite robust, since they do not require developing the implicit equation; instead, they reduce geometric operations on the object to matrix algebra. This is the case of interpolation matrices below, but also of the method based on the theory of syzygies, e.g., [START_REF] Sederberg | Implicitization using moving curves and surfaces[END_REF][START_REF] Busé | Implicit matrix representations of rational Bézier curves and surfaces[END_REF]. It handles base points and yields matrices whose entries are polynomials in the implicit variables and they indirectly represent implicit objects: their rank drops exactly on the curve or surface. Most of the time, the entries of matrices constructed with such techniques are linear in the implicit variables, though there exists methods for constructing implicit matrices with quadratic entries instead. They allow for geometric operations, such as surface-surface intersection [START_REF] Busé | The surface/surface intersection problem by means of matrix based representations[END_REF] and, more recently, ray shooting [START_REF] Shen | A line/trimmed NURBS surface intersection algorithm using matrix representations[END_REF], to be executed by linear algebra. The advantage of these matrices is that they are much smaller than interpolation matrices, and allow for inversion by solving an eigenproblem on them.

Chow forms

Chow forms have been studied in computer algebra, in particular for varieties of codimension 2, since they provide a method to describe the variety by a single polynomial [START_REF] Dalbec | Introduction to chow forms[END_REF][START_REF] Gelfand | Discriminants, Resultants and Multidimensional Determinants[END_REF]. The Chow form of a variety V is basically a polynomial R V which indicates when linear subspaces intersect V . For example, the Chow form of a space curve in projective 3-dimensional space is a polynomial in the indeterminates u ij that vanishes whenever the planes

H 0 = u 00 x 0 + u 01 x 1 + u 02 x 2 + u 03 x 3 = 0, (1)
H 1 = u 10 x 0 + u 11 x 1 + u 12 x 2 + u 13 x 3 = 0,
intersect on the curve. If the space curve is given parametrically, the Chow form represents the variety in terms of R V . It can be computed by a symbolic resultant of the system of linear equations [START_REF] Bajaj | Hermite interpolation of rational space curves using real algebraic surfaces[END_REF] where the set of variables X = (x i ) i is substituted with the parametric equations; the resultant eliminates the parameters and yields a polynomial in the variables U = (u i ) i . The implicit hypersurfaces in X containing the variety have to be extracted through rewriting rules. These make implicitization algorithms that rely on the computation of R V impractical for varieties of high degree and/or dimension.

Due to their complexity, very few implementations exist for computing the Chow forms themselves. Among them, [START_REF] Staglianò | Macaulay2 package "resultants[END_REF] is an implementation in Macaulay2, based on the formula of the Chow form in the Grassmannian space using the Plücker coordinates, and [15, Subroutine 7] is an algorithm using polynomial ring tools and based on a Poisson-like formula of the Chow form.

To formally define the Chow form let Gr(k + 1, n + 1) denote the Grassmannian space of k-dimensional linear projective subspaces of P n . For a variety V ⊂ P n of codimension c, let B(V ) ⊂ P n × Gr(c, n + 1) be the set of (p, L) such that p belongs both to V and to the projective linear subspace L of dimension c -1. Then we obtain V by forgetting the second component in B(V ) and we obtain an hypersurface Z(V ) := {L ∈ Gr(c, n + 1) | L ∩ V = ∅} of the Grassmannian space Gr(c, n + 1) by forgetting the first component in B(V ).

V Z(V ) B(V )
Z(V ) is called the Chow variety of V and has the advantage of being an hypersurface in the Grassmannian space, so it is determined by a unique implicit equation up to a constant factor: the Chow form R V . Despite being determined by a unique equation, Z(V ) describes the variety V of unconstrained (co)dimension, see Proposition 1. Note that when V is a variety of codimension 1, we have Z(V ) V ; this explains why the theory of Chow form is effective only for codimension c > 1. On the other hand, the Chow form of a zero-dimensional variety V = {v 1 , . . . , v k } is also known as the u-resultant. Definition 2. Let V ⊂ P n be a d-dimensional irreducible variety and H 0 , . . . , H d be linear forms where

H i = u i0 x 0 + • • • + u in x n , i = 0, . . . , d (2)
and u ij are new variables, 0 i d, 0 j n. The Chow form R V of V is a polynomial in the variables u ij such that

R V (u ij ) = 0 ⇔ V ∩ {H 0 = 0, . . . , H d = 0} = ∅.
The intersection of the d + 1 hyperplanes H i defined in equation ( 2) is generically a (n -d -1)-dimensional linear subspace L of P n , i.e., an element of the Grassmannian Gr(n -d, n + 1) = Gr(c, n + 1), where c is the codimension of V . Proposition 1. [14, Prop.2.5,p.102] A d-dimensional irreducible subvariety V ⊂ P n is uniquely determined by its Chow form. More precisely, a point ξ ∈ P n lies in V if and only if any (n -d -1)-dimensional plane containing ξ belongs to the Chow variety Z(V ) defined by R V .

One standard way to compute the Chow form would be to proceed as follows. Consider a variety V as in Proposition 1, parameterized as

x j = F j (t), j = 0, . . . , n, t = (t 0 : • • • : t d ), and d + 1 hyperplanes H 0 = • • • = H d = 0
, where H i is defined as in equation [START_REF] Busé | Implicit matrix representations of rational Bézier curves and surfaces[END_REF]. Substituting x j = F j (t) in every equation H i = 0, we would obtain an overdetermined system of equations in the parameters t. We would then saturate these equations by the parameterization polynomials F j (t), j = 0, . . . , n, that is, removing the common solutions of the parameterization (base points) e.g. by using a Gröbner basis method. This step can be ignored if V has no base point (e.g. V is a curve and gcd(F 0 , . . . , F n ) = 1). Then we could eliminate the parameters t by using resultants. This reduces the computation of any Chow form to the computation of a resultant: Corollary 2. Consider any V ⊂ P n of dimension d, with parameterization x j = F j (t), j = 0, . . . , n. Then, the Chow form R V is the resultant of the hyperplane equations H i (0 i d), where one eliminates t.

Here the resultant should be understood, in the sense of Definition 1, as a polynomial that eliminates the variables of the input equations. In other words, while the variables of R V are u ij , its coefficients are polynomials in the coefficients of the parametric homogeneous functions F j .

The coordinates used to represent points in the Grassmannian, and, hence, to describe R V , are most commonly defined as the maximal minors of the (d + 1)×(n+1) matrix whose rows are the normals to the hyperplanes H i . These are known as Plücker coordinates or brackets and are the variables of R V . Brackets are denoted as [j 0 , j 1 , . . . , j d ], where indices correspond to columns of the matrix. Equivalently, the dual Plücker coordinates or dual brackets can be used; these are the maximal minors of a (n -d) × (n + 1) matrix whose rows are n -d points that span the intersection L of the hyperplanes H i . Dual brackets are denoted as [[j 0 , j 1 , . . . , j n-d-1 ]], where indices correspond to columns of the matrix. Brackets and dual brackets with complementary index sets are equal up to sign. There are algorithms to recover the implicit or affine equations of hypersurfaces intersecting on V from R V [START_REF] Gelfand | Discriminants, Resultants and Multidimensional Determinants[END_REF], [START_REF] Sturmfels | Algorithms in Invariant Theory[END_REF].

Assuming that R V is a polynomial in the Plücker coordinates, to obtain a representation for V as intersection of implicit hypersurfaces from its Chow form, one may apply a rewriting method. There are two such methods, namely [14, Cor.2.6,p.102], and [6, Prop.3.1], see also [START_REF] Sturmfels | Algorithms in Invariant Theory[END_REF]. They are not straightforward procedures and, in the case of implicitization, typically yield more implicit polynomials than necessary. However, all implicit polynomials in X are of optimal degree, equal to the degree of V .

To illustrate the approach in [START_REF] Dalbec | Introduction to chow forms[END_REF], let us focus on varieties of codimension 2 in P 3 , i.e., space curves. Consider the planes H 0 , H 1 in equation ( 1). The Chow form R V is a polynomial in the brackets [j 0 , j 1 ], where [j 0 , j 1 ] denotes the maximal minor indexed by the columns 0 j 0 , j 1 3, of the matrix

U := u 00 u 01 u 02 u 03 u 10 u 11 u 12 u 13 .
R V is then rewritten as a polynomial in the dual brackets [[j 0 , j 1 ]], using the relations:

[0, 1] = [[2, 3]], [0, 2] = -[[1, 3]], [0, 3] = [[1, 2]], [1, 2] = [[0, 3]], [1, 3] = -[[0, 2]], [2, 3] = [[0, 1]].
The dual brackets are then substituted by the determinant u 0j0 u 1j1 -u 0j1 u 1j0 of the corresponding minor of U . Finally, the result is expanded as a polynomial whose variables are polynomials in the u 10 , u 11 , u 12 , u 13 and its coefficients are polynomials in the u 00 , u 01 , u 02 , u 03 . The latter polynomials are all of degree equal to the degree of V and form a system of implicit equations of V .

Example 1. As an example, the Chow form of the twisted cubic curve (see also Example 2) with parameterization

(x 0 : x 1 : x 2 : x 3 ) = (s 3 : s 2 t : st 2 : t 3 ), (s, t) ∈ P 1 , (3)
is given by the following determinant in the primal brackets:

det   [0, 1] [0, 2] [0, 3] [0, 2] [0, 3] + [1, 2] [1, 3] [0, 3] [1, 3] [2, 3]   ,
which is also known as the Bézout resultant of the system of equations (1) where we have substituted the X variables with the parameterization in (3). Rewriting this determinant in the dual brackets we obtain:

det   [[2, 3]] -[[1, 3]] [[1, 2]] -[[1, 3]] [[1, 2]] + [[0, 3]] -[[0, 2]] [[1, 2]] -[[0, 2]] [[0, 1]]  
Substituting the dual brackets by the determinant of the corresponding minors of U and collecting the terms in the variables {u 0j0 } and {u 1j1 }, we obtain the Chow form of the twisted cubic:

(u 2 02 u 03 -u 01 u 2
+(2 u 2 01 u 03 -2 u 00 u 02 u 03 )u 10 u 2 12 + (u 00 u 2 01 -u 2 00 u 02 )u 2 12 u 13 ,
where the polynomial coefficients in the {u 0j } variables are the defining implicit equations of the twisted cubic.

Here we adopt a practical method that avoids conversion assuming we have a parametric representation of V . We compute a polynomial whose vanishing is a necessary but not always sufficient condition for a point to lie on the variety. The algorithm we propose in Sections 4 and 5 avoids the computation of the Chow form polynomial and the need for rewriting techniques. For space curves, we shall achieve the result of Proposition 1 by computing only a few selected subsets of Z(V ).

Interpolation matrices

This section describes a direct method to reduce implicitization to linear algebra by constructing an interpolation matrix M , given a plane curve or a (hyper)surface in parametric form or as a point cloud. The matrix is indexed by all possible monomials in the implicit equation (columns) and different values (rows) at which all monomials get evaluated. The vector of coefficients of the implicit equation is in the kernel of M , even in the presence of base points. This idea has been extensively used, e.g. [START_REF] Bajaj | Hermite interpolation of rational space curves using real algebraic surfaces[END_REF][START_REF] Dokken | Approximate implicitization[END_REF]. The matrix is somewhat different in [START_REF] Corless | Numerical implicitization of parametric hypersurfaces with linear algebra[END_REF], which is the method implemented in Maple for implicitization via the implicitize command. The latter method consists of expressing the implicit equations as the kernel of a carefully chosen integral form. It accepts nonalgebraic parameterization and can still return formal implicit formulae provided that the integral form can be expressed with a formal formula (the integrand is polynomial in the parametric input equations). Alternatively, this method can perform floating-point computations and return approximate implicit formulae.

In [START_REF] Emiris | Sparse implicitization by interpolation: Characterizing non-exactness, and an application to computing discriminants[END_REF], sparse elimination theory is employed to predict the implicit monomials and build the interpolation matrix. The monomial set is determined quite tightly for parametric models, by means of the sparse resultant of the parametric polynomials, thus exploiting the sparseness of the parametric and implicit polynomials.

More specifically, if the input object is affinely parameterized by rational functions

x j = f j (t), j = 1, . . . , n, t = (t 1 , . . . , t n-1 ) (4)
then it is possible to predict the implicit monomials. This set is included in the predicted (implicit) polytope computed by software ResPol1 [START_REF] Emiris | An oraclebased, output-sensitive algorithm for projections of resultant polytopes[END_REF]. If the input is a point cloud, we consider a coarse estimation of the monomial set by guessing the total degree of the variety and taking all the monomials of that degree or lower. Let S be the predicted set of implicit monomials and |S| its cardinality.

The set S is used to construct a numerical matrix M , expressing a linear system whose unknowns are the coefficients c i (i = 1, . . . , |S|) of the monomials S in the implicit polynomial, as discussed above. If the input object is a parameterization, we obtain the linear system in the c i by substituting each x j by its rational parametric expression x j = f j (t) in the equation

|S| i=1 c i x ai = 0, x ai := x ai1 1 • • • x ain n .
We then evaluate the parameters t at generic points (randomized in practice)

τ k ∈ C n-1 , k = 1, . . . , µ, µ
|S|, avoiding values that make the denominators of the parametric expressions close to 0. Each evaluation gives a linear equation in the coefficients c i .

Letting m i = m i (t) denote the monomial x ai after substituting each x j by its parametric expression in (4), and m i | t=τ k its evaluation at t = τ k , we end up with a matrix M of the form:

M =     m 1 | t=τ1 • • • m |S| | t=τ1 . . . • • • . . . m 1 | t=τµ • • • m |S| | t=τµ     . (5)
Typically µ = |S| for performing exact kernel computation, and µ = 2|S| for approximate numeric computation.

If the input object is given as a point cloud, we take µ random points out of it (µ |S|) and use them instead of the points evaluated at the parameters τ k , k = 1, . . . , µ.

Let M be the (|S| -1) × |S| numeric matrix obtained by evaluating the monomials S at |S| -1 points τ k , k = 1 . . . , |S| -1. We obtain the |S| × |S| matrix M (x), which is numeric except for its last row, by appending the row of monomials S to matrix M :

M (x) = M S(x) , (6)
where we use the notation S(x) to emphasize that this is the only symbolic row of M (x). Notice that matrices M , M and M (q), for a point q lying on the hypersurface, have the same kernel. Matrix M (x) has an important property:

Lemma 3. [11, Lem. 7]
Assuming M is of full rank, then det M (x) equals the implicit polynomial up to a constant.

Interpolation matrices in higher codimension

This section generalizes the construction of interpolation matrices described in Section 2.2 to the case of varieties of codimension greater than 1.

Let V ⊂ C n be a variety of any codimension given parametrically or as a point cloud. Given a set of monomials S, we randomly pick µ (µ |S|) points x k ∈ V, k = 1, . . . , µ on V either from the input point cloud, or as evaluations

x k = (f 1 (τ k ), . . . , f n (τ k )) of the input parametric equations at random points τ k .
Then we construct the interpolation matrix

M (x) = M S(x) as in Section 2.2,
where M is a numeric submatrix, x = (x 1 , . . . , x n ) are symbolic coordinates and S(x) is a row of symbolic monomials in x.

Recall that the support of a polynomial is the set of monomials appearing with nonzero coefficient. We define the following set of polynomials:

P := {p ∈ C[x] | support(p) ⊂ S and ∀ ξ ∈ V, p(ξ) = 0}.
P is a C-vector space. We assume that S contains all the monomials of a set of generators of the ideal

I(V ), i.e., V = {ξ ∈ C n | ∀p ∈ P, p(ξ) = 0}. This construction is summarized in Algorithm 1.
(2) Pick µ random points (x k ) 1 k µ on V with |S| µ 2|S| (3) Construct M := (m i (x k )) 1 i |S|,1 k µ (4) Construct M (x) := M S(x)
Suppose also that the points x k , k = 1, . . . , µ, defining the rows of M are chosen generically enough, so that for all h ∈ C[x] with monomials in S, we have h(x k ) = 0, ∀k ∈ {1, . . . , µ} ⇐⇒ h ∈ P.

Then, the matrix M (x) has a drop-of-rank property but weaker than that of Lemma 3 in the sense that M is not of full rank.

Lemma 4. Assume we built M (x) using a set S that contains all the monomials of a set of generators of the ideal I(V ). Then, for ξ ∈ C n , ξ belongs to V if and only if Rank(M (ξ)) < Rank(M (x)), where it holds Rank(M (x)) = Rank(M ) + 1.

Proof. Using the basis S of monomials, we consider the canonical complex space of dimension |S|, C |S| . It is isomorphic to the space of polynomials

{p ∈ C[x] | support(p) ⊂ S}.
By abuse of notation, the image of P under this isomorphism will also be called P. By the hypothesis on genericity of x k , k = 1, . . . , µ, we have that Ker(M ) = P. So, for ξ ∈ C n , we have:

ξ ∈ V ⇐⇒ ∀p ∈ P, p(ξ) = 0 (by the hypothesis that P characterizes V ) ⇐⇒ Ker(M (ξ)) = P = Ker(M ) ⇐⇒ Rank(M (ξ)) = Rank(M ) < Rank(M (x)).
In practice, taking µ = |S| and random points x k , k = 1, . . . , µ, is enough to satisfy the hypothesis of Lemma 4. The set of monomials S is hard to determine optimally. One can estimate bounds on the degree of V and take all monomials of degree up to that bound, which usually leads to more monomials and a larger matrix M (x) than needed. If the input is a rational parameterization (f 1 , . . . , f n ) of V , we have an upper bound of deg(V ) given by deg(V )

d i=1 max j (deg ti f j )
where deg ti f j is the maximum of the degrees of the numerator and denominator of f j in the i-th parameter.

The drop-of-rank property readily leads to a computation of the implicit equations representing V set-theoretically, either by computing all the maximal nonzero minors of M (x) containing the last line, which is inefficient in practice, or by computing the nullspace of M .

The matrix representation given in [START_REF] Busé | The surface/surface intersection problem by means of matrix based representations[END_REF] has the same drop-of-rank property but is of smaller size. The construction of their matrix, unlike M (x), relies on syzygy computations and is thus slower. However, that method is overall more efficient because of the faster rank computation at each point evaluation.

Space curves

This section derives implicit representations of parametric space curves by Chow forms. Our methods avoid complex computations, such as the rewriting algorithm. In particular, we avoid the explicit computation of the Chow form and instead focus on a proper subset of the Chow variety that is enough to describe the space curve. Indeed, the Chow form of a space curve vanishes on a space line L if and only if L intersects the space curve. The method presented here provides sets of such lines, not all of them, but enough to be able to retrieve the space curve from them. This is how we proceed:

Suppose that we have a space curve V parameterized as

x j = F j (t), j = 0, . . . , 3, t = (t 0 : t 1 ).

Let the line L be defined by a symbolic point ξ = (ξ 0 : • • • : ξ 3 ) and a sufficiently generic point G ∈ V . Define two planes Aff(G, ξ, P 0 ) and Aff(G, ξ, P 1 )

that intersect along L, by choosing two random points P 0 and P 1 and let H 0 (x 0 : • • • : x 3 ) and H 1 (x 0 : • • • : x 3 ) be their respective implicit equations, as in [START_REF] Bajaj | Hermite interpolation of rational space curves using real algebraic surfaces[END_REF]. The coefficients of H 0 and H 1 are now linear polynomials in ξ. The (homogeneous) Sylvester resultant (see [START_REF] Cox | Using Algebraic Geometry[END_REF]Ch.3, Prop.1.7]) of this system, where we set x j = F j (t), eliminates t and returns a polynomial in ξ which vanishes on V (but not only on V ), thus offering a necessary but not sufficient condition, see Algorithm 2.

Lemma 5. Let δ = deg F j (t), j = 0, . . . , 3 and R G be the Sylvester resultant of

H 0 (F 0 (t) : . . . : F 3 (t)), H 1 (F 0 (t) : . . . : F 3 (t)), (7)
where H 0 , H 1 are defined as above. Then R G is of degree 2δ and factors into a degree δ polynomial defining a surface S G V ⊃ V , and a polynomial E δ L , where E L is a linear polynomial defining the plane passing through points G, P 0 , P 1 .

Proof. The degree of the Sylvester resultant in the coefficients of each of the H 0 , H 1 , is δ. ξ is involved linearly in the coefficients of both H 0 and H 1 , since it is taken to lie in the intersection of the two planes. Hence the degree of the sought polynomial in ξ is 2δ.

It vanishes only in two cases: if ξ belongs to the plane defined by G, P 0 and P 1 , or if ξ belongs to a line passing by G and intersecting V . Hence we can divide (possibly several times) the sought polynomial in ξ by the equation E L of the plane defined by G, P 0 and P 1 , thus obtaining an equation of the conical surface S G V of vertex G and directrix V . Since such a conical surface is of degree δ, its equation is R G /E δ L .

Theorem 6. Let F : P 1 → P 3 , be a homogeneous parameterization of a space curve V and S G k V , k = 1, 2, 3 be three conical surfaces obtained by the method above with 3 different random points G k / ∈ V . We distinguish two cases.

1. If V is not planar and the points G k are not collinear, then V is the only

1-dimensional component of S G1 V ∩ S G2 V ∩ S G3 V . 2. If V is contained in a plane P and if G 1 is not in P, then V = P ∩ S G1 V .
Proof. Case [START_REF] Bajaj | Hermite interpolation of rational space curves using real algebraic surfaces[END_REF]. Since S G k V are three different cones -or cylinders when the vertices G k are at infinity -, they have no 2-dimensional component in common. We first reduce the problem to the case where the vertices are P x := (0 : 1 : 0 : 0), P y := (0 : 0 : 1 : 0) and P z := (0 : 0 : 0 : 1). We shall prove that an algebraic space curve is the intersection of the 3 cylinders spanned by the curve itself and of directions x 1 , x 2 and x 3 respectively.

Let C be an irreducible component of

S G1 V ∩ S G2 V ∩ S G3 V . Since G 1 , G 2 , G 3
are not collinear, there exists a map φ ∈ PGL(4, C) that sends G 1 to P x , G 2 to P y and G 3 to P z . By linearity of φ (in particular, φ preserves alignment), we have φ S G1 V = S Px φ(V ) and similar equalities for φ S G2 V and φ S G3 V . Also,

φ(V ) is a homogeneous variety parameterized by φ • F . Thus, if C ⊂ S G3 V , then φ(C) ⊂ S Pz φ(V )
. By this argument, we only have to prove that there is no homogeneous space curve φ(V ) for which S Px φ(V ) ∩ S Py φ(V ) ∩ S Pz φ(V ) contains a different curve than φ(V ). For convenience, φ(V ) is denoted by V and φ(C) is denoted by C in what follows.

Remark 1. Since V does not lie in the plane at infinity (x 0 = 0), we can switch to the affine setting.

The parameterization in this affine setting is

f : t 1 ∈ C → F 1 (1 : t 1 ) F 0 (1 : t 1 ) , F 2 (1 : t 1 ) F 0 (1 : t 1 ) , F 3 (1 : t 1 ) F 0 (1 : t 1 ) ∈ C 3 .
For convenience, we use the same notations for both the homogeneous parameterization and its restriction to this affine setting t 0 = x 0 = 1. We now have simpler expressions for our surfaces:

• S Px V = {(x 1 , f 2 (t 1 ), f 3 (t 1 )) | x 1 , t 1 ∈ C}, • S Py V = {(f 1 (t 1 ), x 2 , f 3 (t 1 )) | x 2 , t 1 ∈ C}, • S Pz V = {(f 1 (t 1 ), f 2 (t 1 ), x 3 ) | x 3 , t 1 ∈ C}.
Since C ⊂ S Px V , there is a (not necessarily rational) parameterization of C given by q : t 1 ∈ C → (q 1 (t 1 ), f 2 (ϕ(t 1 )), f 3 (ϕ(t 1 ))), where q 1 and ϕ are continuous piecewise smooth maps.

Remark 2. We see that ϕ (resp. q 1 ) is not locally constant: otherwise, a part of C would be included in a straight line (resp. a plane), which contradicts the fact that V is not planar.

We can thus pick a small open disc I ⊂ C such that q |I is injective and so, without loss of generality, we assume that ϕ |I = Id |I . Also, since V is not planar and f is rational, the sets of singular points of the three maps π x1 : t 1 → (f 2 (t 1 ), f 3 (t 1 )), π x2 : t 1 → (f 1 (t 1 ), f 3 (t 1 )) and π x3 : t 1 → (f 1 (t 1 ), f 2 (t 1 )) are finite. Shrinking I if necessary, we assume there is no such singular point in q(I). Now, we have a local curve q(I) = {(q

1 (t 1 ), f 2 (t 1 ), f 3 (t 1 )) | t 1 ∈ I} included in S Px V ∩ S Py V ∩ S Pz V .
Using the fact that it lies on S Py V , we have another parameterization of q(I) given by r = (r 1 , f 2 , r 3 ) with r i injective on I and r i (I) ⊂ f i (C). Similarly, q(I) ⊂ S Pz V gives a third parameterization s = (s 1 , s 2 , f 3 ) with s i injective on I and s i (I) ⊂ f i (C).

Comparing s with q and r, we have s = (s 1 , f 2 , f 3 ). Lastly, since π x1 is regular on I, there is only one branch in π x1 (s(I)) and so s -1 1 (s(I)) = f -1 1 (s(I)). The curve C is thus locally contained in V ; it follows that C = V .

Case (2). Since G 1 ∈ P, the conical surface S G1 V consists only of the union of lines transversal to P: S G1 V = ∪ x∈V Line(G 1 , x). Each of these lines intersects P only in one point x ∈ V so the curve V is exactly P ∩ S G1 V .

Note that Theorem 6 is also valid over the reals, the key argument being the existence of local smooth maps around the (dense) set of regular points.

Varieties of arbitrary codimension

In this section we generalize the construction of Section 4 to varieties of arbitrary codimension. Let V ⊂ P n be a d-dimensional variety parameterized as

x j = F j (t), j = 0, . . . , n, t = (t 0 :

• • • : t d ),
and G = {G 1 , . . . , G n-d-1 } be a set of n -d -1 sufficiently generic points not in V . Choose d + 1 sets of d random points P i = {P i1 , . . . , P id }, none of these points lying in V , such that for i = 0, . . . , d the points in G and in P i form an affinely independent set. Let H i , i = 0, . . . , d be the hyperplane defined as the span of the points ξ, G, P i . Substitute x j = F j (t) in each H i to obtain the system of equations

H 0 (F 0 (t) : • • • : F n (t)) = • • • = H d (F 0 (t) : • • • : F n (t)) = 0. ( 8
)
The resultant of the polynomial system (8) eliminates t and returns a polynomial R G in ξ which vanishes on V (but not only on V ), thus offering a necessary but not sufficient condition, see Algorithm 2.

There are three issues we have to examine when generalizing the algorithm. We do so in the following three subsections.

Computing the resultant in several variables

The resultant computation is the bottleneck of the algorithm. To compute the resultant for d = 1, we can use Sylvester determinantal formula. For arbitrary d there exist rational formulae yielding the resultant as the ratio of two determinants, namely the Macaulay determinant [START_REF] Cox | Using Algebraic Geometry[END_REF] or the sparse resultant matrix and one of its minors [START_REF] White | Macaulay-style formulas for sparse resultants[END_REF]. These formulae are optimal for generic coefficients. For arbitrary coefficients, an infinitesimal perturbation may be applied.

Another option is to use interpolation in conjunction with information on the resultant support. This might be obtained from degree bounds on R G , as explained below, or by the computation of the monomials of the polynomials in (8) using software ResPol from [START_REF] Emiris | An oraclebased, output-sensitive algorithm for projections of resultant polytopes[END_REF]. The latter is analogous to the basic approach for defining the interpolation matrix by the support obtained from the resultant polytope, see Section 2.2.

If we choose to interpolate the resultant, an issue arises at sampling: all generated points ξ lie on V , whereas we are trying to compute a hypersurface containing V . But the kernel of M should have large dimension and, among the kernel vectors, we may choose one or more "small" independent vectors to define the implicit equations. We use independent vectors so as to obtain distinct surfaces. Indeed, with independent vectors, the tangent spaces of the surfaces differ and thus the surfaces intersect transversally. Here "small" may refer to the number of non-zero vector entries, or to the total degree of the monomials corresponding to its non-zero entries.

By continuing developing d times the exterior product, we obtain the relation

∩ d i=0 H i = E(ξ)(ξ ∧ (ΛG)),
where E is a polynomial of degree

d in ξ. So the resultant vanishes if dim(∩ i H i ) > n -d, that is if E(ξ) = 0.
Example 8 shows how this polynomial E is computed for a concrete variety. Since R G (ξ) = 0 if and only if S G V (ξ) = 0 or E(ξ) = 0, the factor E is present and raised to a power p δ d , in the expression of R G .

How many hypersurfaces are sufficient

Computing the resultant and factoring out the extraneous factor given by Lemma 7, yields one implicit polynomial defining a hypersurface that contains the given variety. To achieve the hypothesis of Proposition 1, we must iterate for a few distinct pointsets G thus obtaining implicit polynomials S G V = 0 whose intersection is V .

Unfortunately, we do not yet have an a priori bound ρ for the number of equations needed to describe the variety set-theoretically as in Theorem 6. Experimental results indicate that for curves in C n , n + 1 hypersurfaces of the type S G V = 0 are required. This bound also extends to surfaces in C 4 , where 5 such hypersurfaces are sufficient. The theoretical result in [16, Ch.V] indicates that n hypersurfaces suffice for any variety in C n . It is not clear how to apply this result to the specific type of hypersurfaces (conical) obtained by our method.

These equations are obtained using random pointsets G: the hypothesis of genericity is important. Indeed, each pointset G k (k = 1, . . . , ρ) must obviously consist of affinely independent points not in V in order to define S G k V properly. It is also possible that more implicit polynomials are required for specific (bad) choices of pointsets G k . In particular, if there is a common affine subspace

L ⊂ (∩ k Aff(G k )) and if L ∩ V = ∅, then L ⊂ ∩ k S G k V
and the equations of these conical hypersurfaces are not defining V set-theoretically. To avoid these degenerate cases, it may be interesting to choose the random pointsets G k such that any n + 1 points picked from those pointsets are always affinely independent.

Degree bounds

Before stating the implicitization algorithm, we first examine the degrees of the factors of the resultant polynomial R G .

We showed that E p appears as a factor of R G , where p is possibly very high. For curves (d = 1), p indeed achieves the upper bound δ d . However, in our tests with 1 < d < n and in presence of base points, the factor defining S G V also appears to a power q > 1 in the expression of R G :

R G degree δ d +dδ d = (S G V ) q degree δ d ×q × E p degree d×p
Note that when V is a properly parameterized curve, the inequalities become equalities and we have q = 1 and p = δ. The algorithm works correctly on non-properly parameterized varieties. However, a non-proper parameterization decreases the degree of S G V by some factor (the generic number of preimages) and increases the power degree q by that same factor. In practice, the extraneous factor E seems to always appear with some power p close to its upper bound δ d .

A tighter bound on the degree of R G can be obtained by considering sparse resultants and mixed volumes [START_REF] Cox | Using Algebraic Geometry[END_REF]Ch.7]. Let

M V -i = M V (H 0 , . . . , H i-1 , H i+1 , . . . , H d ), 0 i d,
be the mixed volume of all polynomials excluding H i . The degree of the sparse resultant in the coefficients of H i is known to equal M V -i , therefore its total degree equals

d i=0 M V -i .
For curves in P n , i.e., when d = 1, Algorithm 2 utilizes the Sylvester determinant for computing the resultant instead of the Macaulay or sparse resultant determinant in the general case. This fact, in conjunction with Lemma 7 for determining the extraneous factors of the resultant, make the algorithm much more efficient for d = 1 than in the general case of arbitrary d.

Algorithm 2: Implicit representations of a d-dimensional variety V ⊂ P n Input : -V , parameterized by x j = F j (t), j = 0, . . . , n, -Number of iterations: ρ (= 3 when n = 3) Output: ρ polynomials vanishing on V Repeat ρ times:

(1) Define a (n -d -1)-dimensional linear subspace by affinely independent random points

G = {G 1 , . . . , G n-d-1 } none G k ∈ V (1 j n -d -1)
, and consider symbolic point ξ = (ξ 0 , . . . , ξ n ). (2) Define d + 1 hyperplanes H i through G, ξ, P i , for random point sets P i affinely independent from points in G.

(3) Set x j = F j (t) in the equations of H i : their resultant R G , where we eliminate t, is the sought polynomial in ξ. Compute it by formula, or by interpolation as in Section 2.2. (4) Compute the extraneous factor E using the exterior algebra recursive formula. Then divide R G by E as many times as possible to obtain S G V .

Experiments and performance

This section presents the features of Algorithm 2 and its practical performance through a number of examples. In all these examples we use our Maple implementation of Algorithm 22 on a laptop with a 2GHz Intel Celeron processor running Maple 18. Our algorithm was also implemented in Sage with similar runtimes as in Maple. In some of these examples we compare different existing approaches for implicitization with ours: (Impl) Maple's command implicitize from the package algcurves which is an implementation of the algorithm presented in [START_REF] Corless | Numerical implicitization of parametric hypersurfaces with linear algebra[END_REF].

For all these algorithms the input is the parametric equations. In Figure 2, and whenever we need to consider generic parameterizations, Maple's command randpoly is used to generate these parametric equations. (MI) also needs a superset of the support of the implicit equations. For this, in the examples we use the exact set of monomials. Even though the algorithm's runtime is very competitive, its usage is slower than the usage of the other algorithms: (MI) becomes inefficient as soon as the user needs to check membership of more than ∼ 20 points. This is due to the rank computations performed by (MI) vs the polynomial evaluations of all other algorithms except (Mrep). The latter's usage also requires rank computations but as its matrices are much more compact, these are done more efficiently than (MI). On the upside, (MI) can also accept a point cloud representation as input. Algorithm (Impl) also requires as input the target total degree of the implicit polynomial. Unless specified otherwise, we use the lowest degree that is sufficient for having a generating set of the variety's ideal. Note that Algorithm 2, (CF), can be easily parallelized as the computation of each implicit equation can be done independently. This effectively reduces its runtime by a factor equal to the number of equations required (three for space curves). For a fair comparison, we do not use parallelization in the runtimes presented here because we have not studied the potential of parallelization of the other algorithms.

The output of (CF), (GB) and (Impl) are implicit polynomials. The output of the matrix representation (Mrep) and interpolation (MI) algorithms are symbolic matrices that have the drop of rank property. It is possible to retrieve implicit equations from such matrices by computing their minors of size equal to the maximal rank of the matrix. However, this is not the aim of these algorithms, and quickly becomes too expensive. We could nonetheless compute the minors of the output of (Mrep) for space curves of parametric degree up to 7. These minors yield, up to constant factors, the same equations as the (GB) method. The minors of the matrices computed by (MI) depend on the set of monomials used as input. Computing these minors takes longer than 30 min already for space curves of parametric degree 4. Algorithm 2 outputs a different set of equations. For curves of moderate and high degrees, the Sylvester resultant computation represents more than 99% of our algorithm's computation time.

For curves parameterized by polynomials of degree up to 6 all methods are comparable in running time. For larger degrees Gröbner basis outperform our method by a factor that depends on the degree and the number of monomials in the parametric polynomials and ranges from 10 to 10 3 . In a few cases Algorithm 2 is also slower than the syzygy-based matrices. It is competitive with (Impl) in most situations.

On the upside, Algorithm 2 always outputs 3 implicit equations for space curves while their number depends on the input in the other methods and is typically much larger than 3. Also, our method works on parametric equations that contain additional formal parameters, unlike (Mrep).

The space curve example of the largest degree that we computed is of degree 18. It took about 2 min to implicitize it with (GBb), 3 min with (CF) and 4 min with (Impl). Though the runtime per equation of (Impl) is lower than that of (CF), the number of implicit equations that (Impl) outputs cannot be predicted and it seems to be only dependent on the input degree and dimensions: for space curves of the same parametric degree it always returns the same number of implicit equations. For a space curve of parametric degree larger than 8 and requested outputs of the same degree, their number is greater than 100. In contrast, Algorithm 2 can return the implicit equations one by one, and Theorem 6 guarantees that 3 of them define the curve set-theoretically. For parameterizations of degree up to 4, (Impl) computes a small number of implicit equations whose total degree is lower than the degree of those computed by our method. On the other hand, Algorithm 2 always outputs equations of the same total degree equal to the degree of the parameterization, by Lemma 5.

Examples

In the sequel we switch from the projective to the affine setting and set ξ = (ξ 1 , . . . , ξ n ) ≡ (x 1 , . . . , x n ), for emphasizing these are the implicit variables. For curves in any ambient dimension we compute the resultant directly using the Sylvester matrix, and also by interpolation, employing degree bounds relying on mixed volume. For d > 1 we use the Macaulay, or the sparse resultant matrix.

The Sylvester matrix leads to polynomials of degree twice the degree of the curve for any random points G, as expected. Interpolating the resultant leads to matrices with very large kernels: on the upside, the polynomials we obtain from the kernel vectors are of degree no greater than those of the former method. Moreover, among them we can find a number of polynomials of small degree and, often, smaller than the degree predicted by degree bounds: these polynomials define the variety set-theoretically.

Example 2. Consider the twisted cubic curve V ⊂ C 3 affinely parameterized as:

(x 1 , x 2 , x 3 ) = (t, t 2 , t 3 ), t ∈ C.

An optimal system of implicit equations for V is

x 2 1 -x 2 = x 2 2 -x 1 x 3 = x 1 x 2 -x 3 = 0. ( 9
)
Let L be the line passing through symbolic point ξ = (x 1 , x 2 , x 3 ), and generic point G ∈ V . We define two random planes H 1 (x 1 , x 2 , x 3 ), H 2 (x 1 , x 2 , x 3 ), intersecting at L by considering additional random points P 1 , P 2 / ∈ L, respectively. The Sylvester resultant of H 1 (t, t 2 , t 3 ) = H 2 (t, t 2 , t 3 ) = 0 is a polynomial of degree 6 in ξ which factors into the degree 3 polynomial

32x 2 -16x 3 + 56x 1 x 3 + 16x 1 x 2 -80x 2 2 -32x 2 1 -40x 2 3 + 42x 3 2 -2x 3 x 2 x 1 +56x 3 x 2 -5x 3 x 2 2 + 5x 3 2 x 1 -8x 3 x 1 2 -48x 2 2 x 1 + 24x 2 x 1 2
and the extraneous linear factor raised to the power 3 predicted in Lemma 5. This yields a surface containing V but not of minimal degree. Repeating the procedure 3 times, the ideal of the resulting polynomials equals the ideal defined by the polynomials in [START_REF] Dokken | Approximate implicitization[END_REF].

Alternatively, we may interpolate the Sylvester resultant above. We take as predicted support the lattice points in a 3-simplex of size 6. The 84 × 84 matrix constructed has a kernel of dimension 65. The corresponding 65 kernel polynomials are of degrees from 2 to 6. Among them, we can find the three polynomials in (9) (up to sign). The interpolation matrix built by (MI) is:

          1 1 1
         (10)
In contrast, computing the Chow form as in [6, Section 3.3] gives 16 (homogeneous) implicit equations, all of degree 3, see Example 1.

If we give up the computation of minors of the matrices given by (MI) and (Mrep), we can still check whether a point belongs to the curve with a rank computation. Computing the rank of the 2 × 3 matrix given by (Mrep) at a specified point costs about 2ms on average. For the 7 × 6 matrix given by (MI), it costs about 9ms on average. Consequently, (MI) is faster when one needs to check the ownership of a very small number of points (less than 20). If many points need to be checked or if the preimage parameters are required, then (Mrep) should be preferred over (MI).

Example 3. Consider the space curve V in Figure 1 affinely parameterized as:

(x 1 , x 2 , x 3 ) = 1 -t 2 1 + t 2 , 2t 1 + t 2 , 1 -t 2 1 + t 2 2 , t ∈ C.
It is the intersection of two cylinders:

x 2 1 -x 3 = x 2 2 + x 3 -1 = 0. (11)
Let line L be defined from the symbolic point ξ = (x 1 , x 2 , x 3 ), and "generic" point G ∈ V . Define two random planes H 1 and H 2 that intersect at L, by choosing random points P 1 , P 2 / ∈ L, respectively. Then, the Sylvester resultant of

H 1 1 -2 1 + t 2 , 2t 1 + t 2 , 1 -t 2 1 + t 2 2 , H 2 1 -t 2 1 + t 2 , 2t 1 + t 2 , 1 -t 2 1 + t 2 2
is a polynomial of degree 8 in ξ which factors into the following degree 4 polynomial:

29 + 120x 1 -56x 2 -182x 3 + 168x 1 x 2 -110x 2 1 -78x 2 2 + 90x 2 2 x 3 -12x 2 x 3 + 70x 2 1 x 3 -18x 1 x 3 -120x 1 x 2 2 + 40x 2 1 x 2 + +156x 2 3 + 56x 3 2 + 49x 4 2 + 25x 4 1 - 48x 2 3 x 2 -72x 2 3 x 1 + 12x 2 3 x 2 2 + 28x 3 x 3 2 + 12x 2 3 x 2 1 -168x 3 2 x 1 + 222x 2 2 x 2 1 - 30x 3 x 3 1 -120x 2 x 3 1 -48x 3 x 2 x 1 -102x 3 x 2 2 x 1 + 76x 3 x 2 1
x 2 and the expected extraneous linear factor raised to the power 4. This yields a surface containing V but not of minimal degree. Repeating the procedure 3 times, we obtain three surfaces that intersect on the curve.

The presence of the base points indeed reduced both the degree of the implicit equation containing the surface and the degree of the resultant. Thus we obtained more extraneous factors from the Macaulay construction.

S G V (x 0 : x 1 : x 2 : x 3 : x 4 ) = x 2 2 + x 2 4 + 3x 3 x 2 -2x 4 x 3 + 3x 2 3 -2x 0 x 2 + 2x 4 x 0 - -5x 0 x 3 + 2x 2 0 + 3x 1 x 2 -5x 4 x 1 + 9x 3 x 1 -8x 0 x 1 + 8x 2 1 .
The extraneous factor predicted by the algorithm can be computed without factoring, by using the exterior algebra formulae. For the algorithm, we chose 7 random points: G and P ij , 0 i 2, 1 j 2. The formulae is the following: So the extraneous factor is obtained by computing these four 5 × 5 determinants. Note that they can be reduced to 4 × 4 determinants since their last row consists of ones. b. Consider the following curve of degree 4 in C 3 :

(ξ ∧ G ∧ P 01 ∧ P 02 )• (ξ ∧ G ∧ P 11 ∧ P 12 )• (ξ ∧ G ∧ P 21 ∧ P 22 ) = [det(ξ,
x 1 (t) = -at 12 -t 7 -bt 6 -8t 3 -t + 9

x 2 (t) = (b -3)t 11 + 3t 10 -at 8 -8t 7 + 8t 6 -t x 3 (t) = (a -b)t 12 -2t 11 -t 9 -3t 7 -7t 6 + t 2 + b (CF) computes 3 implicit equations of conical surfaces, of degree 12 in the x i , in about 10h. None of the other algorithms tested (Impl), (GBa), (GBb) and (Mrep) was able to implicitize that curve, either because they cannot handle formal parameters ((Mrep)) or because it took too long to compute ((GB) and (Impl)).

c. Consider Viviani's curve parameterized as:

x 1 = 128 a 5 t 2 /(256 a 8 + 32 a 4 t 2 + t 4 ), x 2 = (256 a 7 t -16 a 3 t 3 )/(256 a 8 + 32 a 4 t 2 + t 4 ), x 3 = (32 a 5 -2 a t 2 ) (16 a 4 + t 2 )/(256 a 8 + 32 a 4 t 2 + t 4 ). Its implicit equations are:

Future work

Interesting mathematical questions arise in studying the complexity of our method, namely how to find tighter bounds for the degree of the computed conical implicit surfaces and how to bound their number in the general case by results similar to Theorem 6. This also requires to study in more detail the extraneous factors of the computed equations and to distinguish them into those coming from the chosen method to compute the resultant, and those inherited in our implicitization method.

As it stands, the algorithm chooses random vertices G for the hypersurfaces. A question arises of how to choose these vertices more efficiently, either by using a better random distribution or by finding particularly interesting points. The answer may depend on the tasks at hand: checking membership is faster when the coefficients are integral and of small size, computing intersections is easier when the same vertices are always used, when evaluating the distance to the variety perpendicular hypersurfaces are better, etc.

A possible improvement of the algorithm at least in 3D space, could be to always compute cylindrical hypersurfaces instead of conical hypersurfaces.

For space curves, we can reduce the computation of three conical surfaces to the computation of three cylinders by using projective linear maps, as in the proof of Theorem 6. The computation of a cylinder in P 3 is equivalent to the computation of its intersection with a perpendicular plane, i.e., a planar curve. Using cylinders would thus allow the use of methods that are specialized to planar curve implicitization.

An interesting direction for future research concerns the possibility of defining a sparse version of the Chow form. In Corollary 2, consider the polynomials in (t 0 : • • • : t d ) obtained after substituting in the equations of the hyperplanes H i , the x i 's by the parametric expressions. These d + 1 polynomial equations in d + 1 homogeneous variables may have Newton polytopes smaller than the simplex of size δ and thus be suitable for a sparse resultant computation. In the sparse resultant method, the monomial bases of the input polynomials are analyzed very closely to deduce a sparse basis of the output resultant. By analogous means, we expect that a sparse Chow form of an input variety V of codimension c in P n would be a hypersurface of a carefully picked subspace of Gr(c, n + 1). In Definition 2, the linear forms H i would have fewer degrees of freedom (by forcing a linear relation between some of the parameters u ij , for instance). It is not obvious how to correctly pick subspace D ⊂ Gr(c, n + 1) such that Z D (V ) := {L ∈ D | L ∩ V = ∅} would still satisfy Proposition 1. This potential generalization of the Chow form may lead to better algorithms for computing Chow forms.

Figure 1 :

 1 Figure 1: Left: The space curve of Example 3. Right: The two implicit surfaces defining the curve.

Algorithm 1 :

 1 Interpolation matrix of a d-dimensional variety V ⊂ C n Input : -V parameterized as in (4), or -A set of monomials S and a point cloud (x k ) k ⊂ V of at least |S| points such that S contains all monomials of a set of generators of I(V ) Output: An interpolation matrix M (x) whose rank drops on V (1) If S is not given then Let δ := d i=1 max j (deg ti f j ) and S = (m 1 , . . . , m |S| ) be the monomial basis of polynomials of degree δ

(

  CF) our main contribution, Algorithm 2 based on Chow forms, (MI) the matrix interpolation Algorithm 1 of Section 3, (GB) a Gröbner bases computation using one of Maple's native commands: (GBa) Basis from the package Groebner or (GBb) EliminationIdeal from the package PolynomialIdeals, (Mrep) a matrix representation algorithm based on syzygies [2, 3], using the Maple code of the authors, and

Figure 2 :

 2 Figure 2: Timing (sec) of various algorithms on random space curves: (CF) implements Alg. 2, (MI) matrix interpolation Alg. 1, (GBa) Gröbner bases by Maple's Basis, (GBb) Gröbner bases by Maple's EliminationIdeal, (Mrep) Maple code for matrix representation algorithm using syzygies, (Impl) Maple's implicitize implementing [4]. (Mrep-Minors) is the time of computing Mrep minors, providing implicit equations. Computing minors by (MI) takes > 20 min at degree 5 thus there is no row for (MI-Minors).

Example 9 .

 9 Algorithm 2 can also handle parametric equations containing additional formal parameters. a. Consider for instance a parameterized cubic:(x 1 , x 2 , x 3 ) = at, bt 2 , a + b 2 t 3 , t ∈ C.Algorithm (CF) computes 3 implicit equations of degree 3 in the x i and of degree 5 in {a, b}. Those equations match the ones obtained for the twisted cubic when a = b = 1.

Figure 4 :

 4 Figure 4: (Left) The curve of Example 9(c.) for a = 1. (Right) Two of its three defining equations found by Algorithm 2. A third equation is needed in order to remove the extraneous surface intersections.

x 2 1 + x 2 2 + x 2 3 = 4 a 2 , 3 -4 a 4 x 2 3 -393216 a 19 x 2 1 x 2 3 +x 2 1 -

 2323331 (x 1 -a) 2 + x 2 2 = a 2 . The 3 implicit equations computed by Algorithm 2 in 7,72 sec are of total degree 4 in x i . Since they are quite large, we show only one of them which is of degree 74 in a and omit most of its terms: 281474976710656 a 36 (-1048576 a 19 x 1 x 2 • • • + 524288 a 18 x 2 3 786432 a 18 x 3 x 3 1 ) = 0. See Fig. 4.

  03 )u 2 10 u 12 + (u 01 u 02 u 03 -u 3 02 )u 2 10 u 13 + (u 01 u 2 03 -u 2 02 u 03 )u 10 u 2 11 +(u 00 u 2 03 -u 01 u 02 u 03 )u 10 u 11 u 12 + (3 u 01 u 2 02 -2 u 2 01 u 03 -u 00 u 02 u 03 )u 10 u 11 u 13 +(u 00 u 01 u 03 -2 u 00 u 2 02 -3 u 2 01 u 02 )u 10 u 12 u 13 + (u 3 03 -u 00 u 01 u 02 )u 10 u 2 13 +(u 3 02 -u 00 u 2 03 )u 3 11 + (3 u 00 u 02 u 03 -3 u 01 u 2 02 )u 2 11 u 12 + (u 2 00 u 03 -u 3 01 )u 3 12 +(2 u 00 u 01 u 03 -2 u 00 u 2 02 )u 2 11 u 13 + (3 u 2 01 u 02 -3 u 00 u 01 u 03 )u 11 u 2 12 +(u 00 u 01 u 02 -u 2 00 u 03 )u 11 u 12 u 13 + (u 2 00 u 02 -u 00 u 2 01 )u 11 u 2 13

Table 1 :

 1 Comparison of the different features of implicitization algorithms

		(CF)	(Impl)	(GB)	(Mrep)
	Runtime
	(representation)	Medium	Medium	Fast	Medium
	Runtime
	(membership)	Immediate	Immediate	Immediate	May be slow
	# of outputs	Low and fixed	None (failure) or very high	Medium, mode-rately varying	Matrix of full rank, varying size
	Degree of outputs	Same as input	Set by user	Same as input or a bit lower	Matrix entries are linear polynomials
	Restrictions on the input	Rational* parameterization, n -d > 1	None (accepts analytic parameterizations)	Rational* parameterization, Non-parametric coefficients	Rational* parameterization, Non-parametric coefficients
	Retrieve point
	parameters
	Randomized	**			**
	Geometric
	Feature	Conical surfaces
	(*) Also accepts homogeneous parameterizations as the conversion rational ↔ homogeneous is
	straightforward
	(**) Depends on the implementation

Table 2 :

 2 Comparison of the different algorithms on the twisted cubic (Example 2)

  G, P 11 , P 12 , P 01 )(ξ ∧ G ∧ P 02 ) -det(ξ, G, P 11 , P 12 , P 02 )(ξ ∧ G ∧ P 01 )] • • (ξ ∧ G ∧ P 21 ∧ P 22 ) = [det(ξ, G, P 11 , P , P 01 ) det(ξ, G, P 21 , P 22 , P 02 )--det(ξ, G, P 11 , P 12 , P 02 ) det(ξ, G, P 21 , P 22 , P 01 )](ξ ∧ G)

http://sourceforge.net/projects/respol

A prototype version for space curves is available at: http://users.uoa.gr/~claroche/publications/ChowFormImplicitize.zip

Acknowledgements: All authors belong to team AROMATH, joint between INRIA Sophia-Antipolis (France) and NKUA. This project has received funding from the European UnionâĂŹs Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 675789 (AR-CADES). We thank L. Busé and B. Mourrain for discussions.

Independently of the resultant algorithm used, though, there is always an extraneous factor in the resultant that is similar to the one pinpointed by Lemma 5. Which leads to:

Identifying the extraneous factor in the resultant

The resultant is indeed always reducible and only one of its irreducible components is relevant for describing V . To address this issue, Lemma 5 can be generalized as follows.

Lemma 7. Let δ = deg F j (t), j = 0, . . . , n and R G be the resultant of the equations [START_REF] Dickenstein | Explicit formulas for the multivariate resultant[END_REF]. Then R G factors into a polynomial defining a hypersurface S G V which is of degree at most δ d (equality holds when there are no base points) and contains V , and an extraneous factor E p , where E is a polynomial of degree d and p δ d .

Proof. We define the generalized conical hypersurface of directrix V and vertices G = (G 1 , • • • , G n-d-1 ) as following:

) By abuse of notation, we shall also denote by S G V the square-free polynomial defining the hypersurface S G V (unique up to a constant factor). Let us first note that R G has a total degree in ξ of (d + 1)δ d . Indeed, the degree in t of every

) is δ, and the coefficients of the H i 's are linear polynomials in ξ. The resultant of these polynomials has degree in the coefficients of each H i bounded by δ d , therefore total degree (d + 1)δ d , see e.g. [START_REF] Cox | Using Algebraic Geometry[END_REF]Thm.3.1]. Now, R G vanishes if and only if the equations ( 8) have a common solution, which happens when: i. either the hyperplanes defined H i = 0, i = 0, . . . , d, intersect along a linear subspace L of dimension n -d and L ∩ V = ∅, or ii. these hyperplanes intersect along a linear subspace of dimension > n -d.

The first case is dealt with by the condition ξ ∈ S G V . It remains to prove that the second case is equivalent to ξ ∈ E with E being a hypersurface of degree d.

In what follows, we use the theory of exterior algebra ([22, Ch.10]) to compute the equations of E.

Then,

Thus,

(-1) i det(ξ, G, P 11 , . . . , P 1d , P 0i )(ξ ∧ (ΛG) ∧ P 01 ∧ . . . ∧ P 0(i-1) ∧ P 0(i+1) ∧ . . .

Alternatively, we interpolate the Sylvester resultant above using as support the lattice points in a 3-simplex of size 8. The 165 × 165 matrix constructed has a kernel of dimension 133. The degrees of the corresponding 133 kernel polynomials vary from 2 to 8. Among them, we find the two quadratic polynomials in [START_REF] Emiris | Geometric operations using sparse interpolation matrices[END_REF].

The equations above were obtained by choosing random points with integral coordinates and relatively close to the curve. When we increase the range or allow for non-integral points, in order to have better chances to avoid the nongeneric points, the size of the coefficients increase. For example, using random points with integral coordinates in a box of size 100 around the curve yields equations as below, where we omit most terms:

On this example, the algorithms (Impl), (GB) and (MI) compute the two equations in [START_REF] Emiris | Geometric operations using sparse interpolation matrices[END_REF]. This cannot be done by Algorithm 2 unless we know which special apex G for the conical surfaces should be used: for this example, points at infinity in the x 1 and x 3 directions allow to use only two equations instead of three. The algorithm (Mrep) computes a 3×5 matrix, resulting in 6 independent minors of degree 3.

Example 4. [20, Example 5.3] Consider the space curve parameterized as ( p1 q , p2 q , p3 q ), where: The denominator q has 2 real roots, -1.143 and 1.13. Consequently, the curve has 3 connected components and so do the surfaces of the implicit equations. We computed the 3 implicit equations in 0.171 sec. Their coefficients are quite large; below we show one equation where the coefficients are rounded. Example 5. We compare four implicitization algorithms on two space curves V 1 , V 2 with parameterizations: V 1 : (x 1 , x 2 , x 3 ) =(5 t 3 -2 t 2 + 4, 2 t 3 -t 2 + 3 t + 6, -9 t 3 -4 t 2 + 3 t + 4), V 2 : (x 1 , x 2 , x 3 ) =(-6 t 5 + 4 t 4 + 9 t 3 -3 t -9, -5 t 5 -4 t 4 + 9 t 2 + 7 t -6, -7 t 5 + 6 t 4 -7 t 3 -t 2 -9 t -9).

Tables 3 and4 The implicit polynomials for curve V 2 obtained as minors of the matrix by (Mrep) are all of degree 4. Since the entries of the matrices constructed by this method are always linear, the degree of the resulting implicit polynomials equals the rank of the matrix. Example 6. Consider the curve in C 4 with parameterization: The curve is defined by 5 implicit equations of degrees 1,2,2,2 and 3. Our algorithm computes 5 equations of degree 6 in 0.06 sec. These equations contain linear extraneous factors raised to the power 3. When these are divided out we obtain 5 degree 3 equations which define the curve set-theoretically.

Example 7. We tested our method on a surface in C 4 with parameterization:

It has an implicit representation defined by 7 equations all of degree 3. When computing the resultant of the hyperplane equations in step 3 of Algorithm 2, we used the Macaulay matrix. While the resultant of the equations is of degree 12, the Macaulay determinant is of degree 15 and factors into 4 polynomials:

Polynomial p 3 is irreducible of degree 4 and contains the surface; p 4 is a quadratic polynomial raised to the power 4; it is the extraneous factor predicted by the algorithm. The other two factors p 1 , p 2 , of degrees 1 and 2, respectively, are extraneous factors from the Macaulay's matrix construction. A total of 5 polynomials like p 3 define the surface set-theoretically.

Example 8. Algorithm 2 works even in the presence of base points. Consider the surface in C 4 with projective parameterization (s : t : u) → p 1 (s : t : u) q(s : t : u) , p 2 (s : t : u) q(s : t : u) , p 3 (s : t : u) q(s : t : u) , p 4 (s : t : u) q(s : t : u) ,

where: p 1 (s : t : u) = s 2 + t 2 -u 2 , p 2 (s : t : u) = u(s + t -u), p 3 (s : t : u) = su -(t -u) 2 , p 4 (s : t : u) = (s -u) 2 + (t -u) 2 -u 2 , q(s : t : u) = st.

Then (0 : 1 : 1) and (1 : 0 : 1) are two base points of the surface. The Macaulay matrix computed by the algorithm is of size 15 and of rank 13 while we expect the resultant to be of degree at most 12. Taking a non-zero minor of size 13 yields a polynomial that can be factored into the following:

• a factor S G V (x 0 : x 1 : x 2 : x 3 : x 4 ) of degree 2 that contains the surface,

• a factor degree 2 that is the extraneous factor predicted by the algorithm, raised to a power 3,

• another factor of degree 2 and two factors of degree 1, one of which is squared.