- Email: [email protected]
Comparing acceleration techniques for the Dixon and Macaulay resultants
Available online at www.sciencedirect.com
Mathematics and Computers in Simulation 80 (2010) 1146–1152
Comparing acceleration techniques for the Dixon and Macaulay resultants Robert H. Lewis ∗ Fordham University, New York, NY 10458, USA Received 28 November 2007; received in revised form 24 March 2008; accepted 16 April 2008 Available online 16 May 2008
Abstract The Bezout–Dixon resultant method for solving systems of polynomial equations lends itself to various heuristic acceleration techniques, previously reported by the present author, which can be extraordinarily effective. In this paper we will discuss how well these techniques apply to the Macaulay resultant. In brief, we find that they do work there with some difficulties, but the Dixon method is greatly superior. That they work at all is surprising and begs theoretical explanation. © 2008 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Resultant; Polynomial system; System of equations; Bezout; Dixon
1. Introduction Using resultants is a well-known method of solving systems of polynomial equations, such as f1 (x, y, z, . . . , a, b, . . .) = 0 f2 (x, y, z, . . . , a, b, . . .) = 0 f3 (x, y, z, . . . , a, b, . . .) = 0 .. . We do not assume homogeneity. Thus, we have n equations in n − 1 variables x, y, z, . . . and a number of parameters a, b, c, . . .. The ground ring is Z, Q, or Z/p. All computations are exact. Solutions are understood to be in an algebraic closure. A resultant is a single polynomial derived from a system of polynomial equations that encapsulates the solution (common zero) to the system. We want to eliminate the variables and have a polynomial in the parameters. A common variation is to have n equations in n variables xi , i = 1, n. Then one of them, say x1 , is considered a “parameter” to bring this into the previous form. The resultant is therefore an equation in that one variable; the others have been eliminated. ∗
Tel.: +1 7188173226. E-mail address: [email protected].
0378-4754/$36.00 © 2008 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2008.04.020
R.H. Lewis / Mathematics and Computers in Simulation 80 (2010) 1146–1152
1147
The Sylvester matrix is the best-known technique to compute resultants [3]. However, it is usually not realistic for more than one variable. In that case the technique due to Bezout, Cayley, and Dixon, modified by Kapur, Saxena, and Yang, is much better [7,1,4]. As refined and extended to avoid the spurious factor problem, this method has succeeded in solving many polynomial systems that others have found intractable [13]. We shall call this simply the Dixon Resultant. Also well known is the method of Macaulay [3], which is a generalization of Sylvester’s. All three of these methods involve the construction of a matrix containing (usually) multivariate polynomials, and then the determinant of that matrix. In previous talks [11,12] and several papers [9,14,13], we have discussed two heuristic techniques to greatly speed up the Dixon Resultant. Here we apply these ideas to Macaulay and compare the times of Dixon and Macaulay on some test problems. We will discuss: • • • • • • •
Dixon and the spurious factor problem. The EDF method (Early Detection of Factors). Problems I, II: geometric. Systems with symmetry: the inner and outer exponent methods. Problem III: a system due to Kotseiras. Problem IV: a system due to Emiris–Canny. Solutions to the above using Dixon and Macaulay.
We thank Manfred Minimair for providing Maple procedures to compute the Macaulay matrices and insight into this method [15]. 2. The Dixon and Macaulay resultant methods We want to decide if there is a common root of n polynomial equations in n − 1 variables x, y, z, . . . and k parameters a, b, . . . fi (x, y, z, . . . , a, b, . . .) = 0, i = 1, . . . , n Here is a brief description of the Bezout–Dixon method. More details are in [7]. • By a two-step process, a certain matrix is created. Its dimensions are usually larger than n. Its entries are polynomials in the parameters. • Ideally, let det = the determinant of the matrix. If the system has a common solution, then det = 0. det involves only the parameters. • Problem: the matrix need not be square, or might have det = 0 identically. Then the method appears to fail. • However, we may continue [7,2]: Find any maximal rank submatrix, say M. Let det = determinant [M]. Existence of a solution implies det = 0. • One must be aware of spurious factors, i.e., the true resultant is usually only a factor of det. So: the problem is to find the right factor in the determinant of a matrix of polynomials. Here is a brief description of the Macaulay method. More details are in [3]. • A sparse square matrix is constructed. Its dimensions are much larger than n. Its entries are polynomials in the parameters. As with Dixon, its determinant is desired. • Sometimes the determinant of a second smaller submatrix must be constructed and divided into the first. Thus, one may think of this second determinant as a “spurious factor” of the first. In this paper we are not concerned with the second matrix. So: the problem is again to find (a factor of) the determinant of a matrix of polynomials. All computations below were done with the author’s computer algebra system Fermat [10], which excels at polynomial and matrix computations of this sort. Recent independent tests [16] show that Fermat is extremely good for gcd
1148
R.H. Lewis / Mathematics and Computers in Simulation 80 (2010) 1146–1152
of multivariate polynomials. Fermat is available for Windows, Mac, Linux, and Unix. We used a desktop Macintosh new in 2004 with 2 GB of RAM. 3. The acceleration techniques 3.1. EDF: Early Detection of Factors We wish to exploit the observed fact that det = determinant [M] has many factors. By elementary row and column manipulations (Gaussian elimination), we discover probable factors of det and pluck them out of M0 ≡ M. Any denominators that form in the matrix are plucked out. This produces a smaller matrix M1 still with polynomial entries, and a list of discovered numerators and denominators. Iterate. Here is a very simple example with integers:
We factor a 2 out of the second column, then a 2 from the second row. Thus:
We change the second row by subtracting 2/9 of the first:
We pull out the denominator 9 from the second row, and factor out 9 from the first column:
We “clean up” or consolidate by dividing out the common factor of 9 from the numerator and denominator lists; any 1 that occurs may be erased and the list compacted. Since the first column is canonically simple, we are finished with one step of the algorithm, and have produced a one-smaller M1 for the next step. The algorithm terminates by pulling out the 7: As expected (since the original matrix contained all integers) the denominator list is empty. The product of all the entries in the numerator list is the determinant, but we never needed to deal with any number larger than 9. This method is described in detail in [9]. • Note 1: There is no guarantee that det will be completely factored this way. • Note 2: This is a poor way to compute the determinant of a random matrix. But the Dixon and Macaulay matrices are far from random. Let us apply this technique to some test problems with both Dixon and Macaulay. Problem I (A parameter elimination problem). Let the point (x, y) be defined in terms of parameters (u, v, w): u2 − 2xu + v2 − 2yv + x2 + y2 − 1 = 0, v2 − u3 = 0, −3u2 v + 3yu2 − 2vu + 2xv = 0, 6w2 u2 v − 3u2 w − 2wv + 1 = 0
R.H. Lewis / Mathematics and Computers in Simulation 80 (2010) 1146–1152
1149
The parameter elimination problem is to eliminate (u, v, w) and have one equation in only x and y. The Dixon matrix is 16 × 16. 55% of the entries are 0. EDF takes 0.28 s and produces 27 factors having number of terms: 3, 1, 1, 1, 1, 3, 1, 3, 1, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 21, 21, 3 The answer is the piece of 21 terms: 729x8 + 216x7 + 729y2 x6 − 2900x6 − 1458y2 x5 − 2376x5 − 2619y2 x4 + 3870x4 − 1458y4 x3 − 4892y2 x3 + 4072x3 + 729y4 x2 − 297y2 x2 − 1188x2 − 4158y4 x + 5814y2 x − 1656x + 729y6 − 1685y4 + 427y2 + 529 But straight Gaussian elimination on the matrix gives the det in 0.1 s, with 211 terms. Factoring it takes 0.06 s, so EDF loses (but it does take less space). The Macaulay matrix is 286 × 286. It is extremely sparse. However, its determinant is 0! This happens under certain situations. Instead of giving up, let us try the Dixon idea and find a maximal minor. This is not typically done with Macaulay [15]. The rank is 268. EDF takes 16 s, and works. The numerator list is 3, 3, 3, 3, 3, 2, 3, 6, 6, 3, 6, 3, 6, 3, 1, 1, 1, 21, 21, 9, 3, 34, 14, 45 Note that the answer again appears twice (the polynomial of 21 terms). Gaussian elimination plus factoring takes 10.4 s but more space than EDF. Problem II (Heron’s formula in three and four dimensions). Heron’s formula is the classical relation between the area A and the sides a, b, c of a triangle, A2 = s(s − a)(s − b)(s − c), where s is the semi-perimeter. In dimension three (or higher dimension n) one can similarly ask for a formula relating the volume and sides of a tetrahedron (or n-simplex). There is a well-known solution in terms of the Cayley–Menger determinant [6], but here we will use the problem to generate systems of polynomial equations. In three dimensions, put one vertex at the origin, another on the x-axis at (a, 0, 0), a third at (x, y, 0), and the fourth at (x1, y1, z1). If the volume is v and the lengths of the sides are a, b, c, d, e, f , one easily derives six equations. x2 + y 2 − c 2 = 0 (x − a)2 + y2 − b2 = 0 x12 + y12 + z12 − f 2 = 0 (x1 − a)2 + y12 + z12 − d 2 = 0 (x1 − x)2 + (y1 − y)2 + z12 − e2 = 0 6v − a y z1 = 0 We seek to eliminate the coordinates x, y, x1, y1, z1 and get an equation for v in terms of the sides. The Dixon matrix for the above system is 13 × 13. Gaussian elimination plus factoring takes 0.04 s. EDF takes 0.06 s. Both use 9 MB of RAM. The solution polynomial has 23 terms. The corresponding Macaulay matrix is 792 × 792. The rank is 662. We extract a maximal minor. EDF takes 3.5 min, 10 MB of RAM, and produces the correct answer. But Gaussian elimination takes about an hour and uses 1.1 GB of RAM. The determinant has 7,000,000 terms and is too large to factor. In four dimensions, the idea generalizes immediately to a set of ten equations. The Dixon matrix is 63 × 63. EDF takes 3.5 s and 15 MB of RAM. The numerator list has 150 entries ending in polynomials with number of terms: 131, 1, 1, 6, 22, 1, 6, 615, 6, 1, 22, 615, 1, 1471. The answer polynomial, of 131 terms, is a factor of three of these entries. Gaussian elimination was killed after 5 min and 500 MB of RAM. • The corresponding Macaulay matrix is too large to compute. • The five dimensional problem is also solvable by Dixon. 3.2. Symmetry and exponential patterns in the determinant In some problems the Dixon resultant f (x) appears in det as either f (xn ) or f (x)n . These exponential patterns often occur in systems of equations that exhibit cyclic symmetry. We do not have a theoretical explanation for why
1150
R.H. Lewis / Mathematics and Computers in Simulation 80 (2010) 1146–1152
the symmetry produces this effect. In [9] we describe a technique that uses these patterns to efficiently compute f (x). Here we summarize the method and explore two motivating problems with both Dixon and Macaulay. Each problem is really a family of problems following a simple pattern, indexed by an integer n > 1. It is easy to solve either problem for small n, say less than 5. We use the experience gained in the small cases to solve the case of interest. The ground ring is Z. The basic technique is to use the Chinese Remainder Theorem, working mod p with carefully selected primes p. We call these the exponent methods. Problem III (A problem from dynamical systems theory). At ACA 2003, from dynamical systems theory, I. Kotsireas discussed a system of n + 1 equations in variables x1 , . . . , xn , mu [8]. For example, when n = 4 the equations are (1 − x2 )/2 − mu(1 − x12 )/4 = 0 (1 − x3 )/2 − mu(1 − x22 )/4 = 0 (1 − x4 )/2 − mu(1 − x32 )/4 = 0 (1 − x1 )/2 − mu(1 − x42 )/4 = 0 mu4 x1 x2 x3 x4 + 1 = 0 The system can be greatly simplified with a few substitutions to yield an equivalent system: −2x2 + x12 − mu = 0 −2x3 + x22 − mu = 0 −2x4 + x32 − mu = 0 −2x1 + x42 − mu = 0 x1 x2 x3 x4 + 1 = 0 The system (but not each individual equation) is invariant under cyclic permutation (x1 , x2 , x3 , x4 ) of the variables to be eliminated. The ground ring is Z. The case of interest was n = 8. We will work with the simplified system. We want to eliminate the xi and have a resultant in mu. Experimentation with small n quickly reveals a pattern: if the Dixon matrix M is r × r, its determinant det has degree r − 1 in mu, is monic, and factors as f (mu)n s(mu). s(mu) is spurious and the desired answer is f (mu). Let D(mu) be the determinant, f (mu) the answer, s(mu) the spurious factor. Then D(mu) = f (mu)n s(mu), or D(mu) 1/n f (mu) = (*) s(mu) By solving the problem first modulo some p, we know the spurious factor s(mu) over Z, since its coefficients are small. The main idea: With the Chinese Remainder Theorem and interpolation, we will compute f (mu), not D(mu). Values v in Z/pZ will be plugged in for mu in (∗). It is easy to compute D(v), and since we have s(mu), s(v) is trivial to compute. The crux of the matter is the root extraction. The point is to use primes p for which we can compute the n th root. Using elementary number theory, we easily find such primes [9] and root extraction is very easy. Features have been added to Fermat to aid the process. Relative to a straightforward call of Det(M) over Z, execution time is divided by approximately n2 since the interpolation degree is divided by n and the integer coefficient length is divided by n. Further details are in [9]. This solves the n = 8 case much more efficiently than the method of Kotsireas. We call this the outer exponent method. Does the Macaulay matrix also lend itself to this technique? For example, when n = 4, the Dixon matrix is 33 × 33. The determinant is f4 (mu)4 times a spurious factor of eight terms, where f4 (mu) is mu6 − 12mu5 + 47mu4 − 188mu3 + 527mu2 + 4913. When n = 5 the matrix is 81 × 81. The determinant is f5 (mu)5 times a spurious factor of four terms. f5 has degree 15. Amazingly, a similar pattern does happen with Macaulay. The matrix for n = 4 is 495 × 495 and has exactly the same determinant as Dixon. However, when n = 5, the matrix is 3003 × 3003 but its rank is 2787. It is necessary to find a maximal minor, and the determinant varies by minor. The same f5 (mu)5 from Dixon occurs and so does the same
R.H. Lewis / Mathematics and Computers in Simulation 80 (2010) 1146–1152
1151
spurious factor of four terms. However, the existence of other spurious factors and a numerical content would make it difficult to actually run the exponent method with Macaulay. The larger size of the matrices impacts the determinant computation time, but only by about a factor of two. The Macaulay matrix here is so much larger than the Dixon that even computing it is problematical for n > 5. Problem IV (Cyclic-n). This is the well-known system due to Emiris–Canny, of n equations in variables x1 , . . . , xn [5,17]. For example, when n = 4 the equations are x4 + x3 + x2 + x1 = 0, x3 x4 + x1 x4 + x2 x3 + x1 x2 = 0, x2 x3 x4 + x1 x3 x4 + x1 x2 x4 + x1 x2 x3 = 0, x1 x2 x3 x4 − 1 = 0 We will eliminate all variables but x1 . Each equation is invariant under cyclic permutation of all the variables (not just the ones being eliminated). The ground ring is Z. We find a situation similar to the previous problem, but now D(x1 ) = s(x1 ) ∗ f (x1 n ). Let t = x1 n . Then f (t) =
D(t 1/n ) s(t 1/n )
Similar to Problem III, we choose only primes p where x → xn is a bijection mod p. When interpolating, for each value v plug in D(v1/n ). Execution time now is divided by approximately n since the degree is divided by n and the numerical coefficients do not change much (they shrink a bit). We call this the inner exponent method. Dixon and Macaulay again show the same pattern. The Dixon matrix is 4 × 4. The Macaulay matrix is 84 × 84 2 and has rank 80. Both determinants compute in less than .01 s. The Dixon answer is x1 4 (x1 4 − 1) and the Macaulay 2 is x1 6 (x1 4 − 1) . The latter has a spurious factor of x1 2 , but we see the pattern f (x1 n ). When n = 5, the Dixon matrix is 22 × 22. The determinant can be computed in about 0.1 s and is 10
2
x1 20 (x1 5 − 1) (x1 10 + 123x1 5 + 1)
The Macaulay matrix is 1365 × 1365 but has rank 1354. The determinant of a maximal minor can be computed in about 0.7 s and is exactly the same as the Dixon, except x1 is raised to the 35. 4. Conclusions • The EDF method is effective on the Macaulay matrix, as that determinant also exhibits many factors. • Apparently, it is not well known that one can extract a maximal minor from a Macaulay matrix and work with it in the same way as with Dixon. • Exponential patterns show up in the Macaulay matrix as they do with Dixon. This is surprising, and begs a theoretical explanation. • We find that on problems of interest, Macaulay is always slower than Dixon and needs much more storage space. Sometimes the Macaulay matrix is so large it cannot be computed. There would seem to be no practical reason to use Macaulay instead of Dixon. References [1] P. Bikker, On Bezout’s method for computing the resultant. RISC-Linz Report Series, Johannes Kepler University, Linz, Austria, 1995. [2] L. Buse, M. Elkadi, B. Mourrain, Generalized resultants over unirational algebraic varieties, Journal of Symbolic Computations 29 (2000) 515–526. [3] D. Cox, J. Little, D. O’Shea, Using Algebraic Geometry. Graduate Texts in Mathematics, 185, Springer-Verlag, New York, 1998. [4] A.L. Dixon, The eliminant of three quantics in two independent variables, Proceedings of the London Mathematics Society 6 (1908) 468–478.
1152
R.H. Lewis / Mathematics and Computers in Simulation 80 (2010) 1146–1152
[5] I.Z. Emiris, J.F. Canny, Efficient incremental algorithms for the sparse resultant and the mixed volume, Journal of Symbolic Computations 20 (2) (1995) 117–149. [6] T.F. Havel, Some examples of the use of distances as coordinates for Euclidean geometry, Journal of Symbolic Computations 11 (1991) 579–593. [7] D. Kapur, T. Saxena, L. Yang, Proc. of the International Symposium on Symbolic and Algebraic Computation Algebraic and geometric reasoning using Dixon resultants, A.C.M. Press, 1994. [8] I.S. Kotsireas, K. Karamanos, Exact computation of the bifurcation point B4 of the logistic map and the Bailey–Broadhurst conjectures, International Journal of Bifurcation and Chaos 14 (7) (2004) 2417–2423. [9] R.H. Lewis, Heuristics to accelerate the Dixon resultant, Mathematics and Computers in Simulation 77 (4) (2008) 400–407. [10] R.H. Lewis, Computer Algebra System Fermat. http://home.bway.net/lewis/ (2008). [11] R.H. Lewis, Using the Dixon resultant on big problems, in: CBMS, Conference, Texas, A.M., University, May 2002. http://www.math.tamu.edu/conferences/cbms/abs.html. [12] R.H. Lewis, Exploiting Symmetry in a Polynomial System with the Dixon Resultant, International Conference on Applications of Computer Algebra (ACA), Lamar University, Texas, 2004, July. [13] R.H. Lewis, S. Bridgett, Conic tangency equations arising from Apollonius problems in biochemistry, Mathematics and Computers in Simulation 61 (2) (2003) 101–114. [14] R.H. Lewis, E.A. Coutsias, Algorithmic search for flexibility using resultants of polynomial systems, in: F. Botana, T. Recio (Eds.), Automated Deduction in Geometry: ADG 2006, Pontevedra, Spain, vol. 4869 of Lecture Notes in Computer Science, Springer, Berlin/Heidelberg, 2007. [15] M. Minimair, Macaulay Resultant Package. http://minimair.org/software/ (2007). [16] D. Robertz, V. Gerdt, comparison of software systems, http://home.bway.net/lewis/fermat/gcdcomp (2005). [17] A. Sommese, J. Verschelde, C. Wampler, Numerical decomposition of the solution sets of polynomial systems into irreducible components, SINUM 38 (6) (2001) 2022–2046.
Mathematics and Computers in Simulation 80 (2010) 1146–1152
Comparing acceleration techniques for the Dixon and Macaulay resultants Robert H. Lewis ∗ Fordham University, New York, NY 10458, USA Received 28 November 2007; received in revised form 24 March 2008; accepted 16 April 2008 Available online 16 May 2008
Abstract The Bezout–Dixon resultant method for solving systems of polynomial equations lends itself to various heuristic acceleration techniques, previously reported by the present author, which can be extraordinarily effective. In this paper we will discuss how well these techniques apply to the Macaulay resultant. In brief, we find that they do work there with some difficulties, but the Dixon method is greatly superior. That they work at all is surprising and begs theoretical explanation. © 2008 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Resultant; Polynomial system; System of equations; Bezout; Dixon
1. Introduction Using resultants is a well-known method of solving systems of polynomial equations, such as f1 (x, y, z, . . . , a, b, . . .) = 0 f2 (x, y, z, . . . , a, b, . . .) = 0 f3 (x, y, z, . . . , a, b, . . .) = 0 .. . We do not assume homogeneity. Thus, we have n equations in n − 1 variables x, y, z, . . . and a number of parameters a, b, c, . . .. The ground ring is Z, Q, or Z/p. All computations are exact. Solutions are understood to be in an algebraic closure. A resultant is a single polynomial derived from a system of polynomial equations that encapsulates the solution (common zero) to the system. We want to eliminate the variables and have a polynomial in the parameters. A common variation is to have n equations in n variables xi , i = 1, n. Then one of them, say x1 , is considered a “parameter” to bring this into the previous form. The resultant is therefore an equation in that one variable; the others have been eliminated. ∗
Tel.: +1 7188173226. E-mail address: [email protected].
0378-4754/$36.00 © 2008 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2008.04.020
R.H. Lewis / Mathematics and Computers in Simulation 80 (2010) 1146–1152
1147
The Sylvester matrix is the best-known technique to compute resultants [3]. However, it is usually not realistic for more than one variable. In that case the technique due to Bezout, Cayley, and Dixon, modified by Kapur, Saxena, and Yang, is much better [7,1,4]. As refined and extended to avoid the spurious factor problem, this method has succeeded in solving many polynomial systems that others have found intractable [13]. We shall call this simply the Dixon Resultant. Also well known is the method of Macaulay [3], which is a generalization of Sylvester’s. All three of these methods involve the construction of a matrix containing (usually) multivariate polynomials, and then the determinant of that matrix. In previous talks [11,12] and several papers [9,14,13], we have discussed two heuristic techniques to greatly speed up the Dixon Resultant. Here we apply these ideas to Macaulay and compare the times of Dixon and Macaulay on some test problems. We will discuss: • • • • • • •
Dixon and the spurious factor problem. The EDF method (Early Detection of Factors). Problems I, II: geometric. Systems with symmetry: the inner and outer exponent methods. Problem III: a system due to Kotseiras. Problem IV: a system due to Emiris–Canny. Solutions to the above using Dixon and Macaulay.
We thank Manfred Minimair for providing Maple procedures to compute the Macaulay matrices and insight into this method [15]. 2. The Dixon and Macaulay resultant methods We want to decide if there is a common root of n polynomial equations in n − 1 variables x, y, z, . . . and k parameters a, b, . . . fi (x, y, z, . . . , a, b, . . .) = 0, i = 1, . . . , n Here is a brief description of the Bezout–Dixon method. More details are in [7]. • By a two-step process, a certain matrix is created. Its dimensions are usually larger than n. Its entries are polynomials in the parameters. • Ideally, let det = the determinant of the matrix. If the system has a common solution, then det = 0. det involves only the parameters. • Problem: the matrix need not be square, or might have det = 0 identically. Then the method appears to fail. • However, we may continue [7,2]: Find any maximal rank submatrix, say M. Let det = determinant [M]. Existence of a solution implies det = 0. • One must be aware of spurious factors, i.e., the true resultant is usually only a factor of det. So: the problem is to find the right factor in the determinant of a matrix of polynomials. Here is a brief description of the Macaulay method. More details are in [3]. • A sparse square matrix is constructed. Its dimensions are much larger than n. Its entries are polynomials in the parameters. As with Dixon, its determinant is desired. • Sometimes the determinant of a second smaller submatrix must be constructed and divided into the first. Thus, one may think of this second determinant as a “spurious factor” of the first. In this paper we are not concerned with the second matrix. So: the problem is again to find (a factor of) the determinant of a matrix of polynomials. All computations below were done with the author’s computer algebra system Fermat [10], which excels at polynomial and matrix computations of this sort. Recent independent tests [16] show that Fermat is extremely good for gcd
1148
R.H. Lewis / Mathematics and Computers in Simulation 80 (2010) 1146–1152
of multivariate polynomials. Fermat is available for Windows, Mac, Linux, and Unix. We used a desktop Macintosh new in 2004 with 2 GB of RAM. 3. The acceleration techniques 3.1. EDF: Early Detection of Factors We wish to exploit the observed fact that det = determinant [M] has many factors. By elementary row and column manipulations (Gaussian elimination), we discover probable factors of det and pluck them out of M0 ≡ M. Any denominators that form in the matrix are plucked out. This produces a smaller matrix M1 still with polynomial entries, and a list of discovered numerators and denominators. Iterate. Here is a very simple example with integers:
We factor a 2 out of the second column, then a 2 from the second row. Thus:
We change the second row by subtracting 2/9 of the first:
We pull out the denominator 9 from the second row, and factor out 9 from the first column:
We “clean up” or consolidate by dividing out the common factor of 9 from the numerator and denominator lists; any 1 that occurs may be erased and the list compacted. Since the first column is canonically simple, we are finished with one step of the algorithm, and have produced a one-smaller M1 for the next step. The algorithm terminates by pulling out the 7: As expected (since the original matrix contained all integers) the denominator list is empty. The product of all the entries in the numerator list is the determinant, but we never needed to deal with any number larger than 9. This method is described in detail in [9]. • Note 1: There is no guarantee that det will be completely factored this way. • Note 2: This is a poor way to compute the determinant of a random matrix. But the Dixon and Macaulay matrices are far from random. Let us apply this technique to some test problems with both Dixon and Macaulay. Problem I (A parameter elimination problem). Let the point (x, y) be defined in terms of parameters (u, v, w): u2 − 2xu + v2 − 2yv + x2 + y2 − 1 = 0, v2 − u3 = 0, −3u2 v + 3yu2 − 2vu + 2xv = 0, 6w2 u2 v − 3u2 w − 2wv + 1 = 0
R.H. Lewis / Mathematics and Computers in Simulation 80 (2010) 1146–1152
1149
The parameter elimination problem is to eliminate (u, v, w) and have one equation in only x and y. The Dixon matrix is 16 × 16. 55% of the entries are 0. EDF takes 0.28 s and produces 27 factors having number of terms: 3, 1, 1, 1, 1, 3, 1, 3, 1, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 21, 21, 3 The answer is the piece of 21 terms: 729x8 + 216x7 + 729y2 x6 − 2900x6 − 1458y2 x5 − 2376x5 − 2619y2 x4 + 3870x4 − 1458y4 x3 − 4892y2 x3 + 4072x3 + 729y4 x2 − 297y2 x2 − 1188x2 − 4158y4 x + 5814y2 x − 1656x + 729y6 − 1685y4 + 427y2 + 529 But straight Gaussian elimination on the matrix gives the det in 0.1 s, with 211 terms. Factoring it takes 0.06 s, so EDF loses (but it does take less space). The Macaulay matrix is 286 × 286. It is extremely sparse. However, its determinant is 0! This happens under certain situations. Instead of giving up, let us try the Dixon idea and find a maximal minor. This is not typically done with Macaulay [15]. The rank is 268. EDF takes 16 s, and works. The numerator list is 3, 3, 3, 3, 3, 2, 3, 6, 6, 3, 6, 3, 6, 3, 1, 1, 1, 21, 21, 9, 3, 34, 14, 45 Note that the answer again appears twice (the polynomial of 21 terms). Gaussian elimination plus factoring takes 10.4 s but more space than EDF. Problem II (Heron’s formula in three and four dimensions). Heron’s formula is the classical relation between the area A and the sides a, b, c of a triangle, A2 = s(s − a)(s − b)(s − c), where s is the semi-perimeter. In dimension three (or higher dimension n) one can similarly ask for a formula relating the volume and sides of a tetrahedron (or n-simplex). There is a well-known solution in terms of the Cayley–Menger determinant [6], but here we will use the problem to generate systems of polynomial equations. In three dimensions, put one vertex at the origin, another on the x-axis at (a, 0, 0), a third at (x, y, 0), and the fourth at (x1, y1, z1). If the volume is v and the lengths of the sides are a, b, c, d, e, f , one easily derives six equations. x2 + y 2 − c 2 = 0 (x − a)2 + y2 − b2 = 0 x12 + y12 + z12 − f 2 = 0 (x1 − a)2 + y12 + z12 − d 2 = 0 (x1 − x)2 + (y1 − y)2 + z12 − e2 = 0 6v − a y z1 = 0 We seek to eliminate the coordinates x, y, x1, y1, z1 and get an equation for v in terms of the sides. The Dixon matrix for the above system is 13 × 13. Gaussian elimination plus factoring takes 0.04 s. EDF takes 0.06 s. Both use 9 MB of RAM. The solution polynomial has 23 terms. The corresponding Macaulay matrix is 792 × 792. The rank is 662. We extract a maximal minor. EDF takes 3.5 min, 10 MB of RAM, and produces the correct answer. But Gaussian elimination takes about an hour and uses 1.1 GB of RAM. The determinant has 7,000,000 terms and is too large to factor. In four dimensions, the idea generalizes immediately to a set of ten equations. The Dixon matrix is 63 × 63. EDF takes 3.5 s and 15 MB of RAM. The numerator list has 150 entries ending in polynomials with number of terms: 131, 1, 1, 6, 22, 1, 6, 615, 6, 1, 22, 615, 1, 1471. The answer polynomial, of 131 terms, is a factor of three of these entries. Gaussian elimination was killed after 5 min and 500 MB of RAM. • The corresponding Macaulay matrix is too large to compute. • The five dimensional problem is also solvable by Dixon. 3.2. Symmetry and exponential patterns in the determinant In some problems the Dixon resultant f (x) appears in det as either f (xn ) or f (x)n . These exponential patterns often occur in systems of equations that exhibit cyclic symmetry. We do not have a theoretical explanation for why
1150
R.H. Lewis / Mathematics and Computers in Simulation 80 (2010) 1146–1152
the symmetry produces this effect. In [9] we describe a technique that uses these patterns to efficiently compute f (x). Here we summarize the method and explore two motivating problems with both Dixon and Macaulay. Each problem is really a family of problems following a simple pattern, indexed by an integer n > 1. It is easy to solve either problem for small n, say less than 5. We use the experience gained in the small cases to solve the case of interest. The ground ring is Z. The basic technique is to use the Chinese Remainder Theorem, working mod p with carefully selected primes p. We call these the exponent methods. Problem III (A problem from dynamical systems theory). At ACA 2003, from dynamical systems theory, I. Kotsireas discussed a system of n + 1 equations in variables x1 , . . . , xn , mu [8]. For example, when n = 4 the equations are (1 − x2 )/2 − mu(1 − x12 )/4 = 0 (1 − x3 )/2 − mu(1 − x22 )/4 = 0 (1 − x4 )/2 − mu(1 − x32 )/4 = 0 (1 − x1 )/2 − mu(1 − x42 )/4 = 0 mu4 x1 x2 x3 x4 + 1 = 0 The system can be greatly simplified with a few substitutions to yield an equivalent system: −2x2 + x12 − mu = 0 −2x3 + x22 − mu = 0 −2x4 + x32 − mu = 0 −2x1 + x42 − mu = 0 x1 x2 x3 x4 + 1 = 0 The system (but not each individual equation) is invariant under cyclic permutation (x1 , x2 , x3 , x4 ) of the variables to be eliminated. The ground ring is Z. The case of interest was n = 8. We will work with the simplified system. We want to eliminate the xi and have a resultant in mu. Experimentation with small n quickly reveals a pattern: if the Dixon matrix M is r × r, its determinant det has degree r − 1 in mu, is monic, and factors as f (mu)n s(mu). s(mu) is spurious and the desired answer is f (mu). Let D(mu) be the determinant, f (mu) the answer, s(mu) the spurious factor. Then D(mu) = f (mu)n s(mu), or D(mu) 1/n f (mu) = (*) s(mu) By solving the problem first modulo some p, we know the spurious factor s(mu) over Z, since its coefficients are small. The main idea: With the Chinese Remainder Theorem and interpolation, we will compute f (mu), not D(mu). Values v in Z/pZ will be plugged in for mu in (∗). It is easy to compute D(v), and since we have s(mu), s(v) is trivial to compute. The crux of the matter is the root extraction. The point is to use primes p for which we can compute the n th root. Using elementary number theory, we easily find such primes [9] and root extraction is very easy. Features have been added to Fermat to aid the process. Relative to a straightforward call of Det(M) over Z, execution time is divided by approximately n2 since the interpolation degree is divided by n and the integer coefficient length is divided by n. Further details are in [9]. This solves the n = 8 case much more efficiently than the method of Kotsireas. We call this the outer exponent method. Does the Macaulay matrix also lend itself to this technique? For example, when n = 4, the Dixon matrix is 33 × 33. The determinant is f4 (mu)4 times a spurious factor of eight terms, where f4 (mu) is mu6 − 12mu5 + 47mu4 − 188mu3 + 527mu2 + 4913. When n = 5 the matrix is 81 × 81. The determinant is f5 (mu)5 times a spurious factor of four terms. f5 has degree 15. Amazingly, a similar pattern does happen with Macaulay. The matrix for n = 4 is 495 × 495 and has exactly the same determinant as Dixon. However, when n = 5, the matrix is 3003 × 3003 but its rank is 2787. It is necessary to find a maximal minor, and the determinant varies by minor. The same f5 (mu)5 from Dixon occurs and so does the same
R.H. Lewis / Mathematics and Computers in Simulation 80 (2010) 1146–1152
1151
spurious factor of four terms. However, the existence of other spurious factors and a numerical content would make it difficult to actually run the exponent method with Macaulay. The larger size of the matrices impacts the determinant computation time, but only by about a factor of two. The Macaulay matrix here is so much larger than the Dixon that even computing it is problematical for n > 5. Problem IV (Cyclic-n). This is the well-known system due to Emiris–Canny, of n equations in variables x1 , . . . , xn [5,17]. For example, when n = 4 the equations are x4 + x3 + x2 + x1 = 0, x3 x4 + x1 x4 + x2 x3 + x1 x2 = 0, x2 x3 x4 + x1 x3 x4 + x1 x2 x4 + x1 x2 x3 = 0, x1 x2 x3 x4 − 1 = 0 We will eliminate all variables but x1 . Each equation is invariant under cyclic permutation of all the variables (not just the ones being eliminated). The ground ring is Z. We find a situation similar to the previous problem, but now D(x1 ) = s(x1 ) ∗ f (x1 n ). Let t = x1 n . Then f (t) =
D(t 1/n ) s(t 1/n )
Similar to Problem III, we choose only primes p where x → xn is a bijection mod p. When interpolating, for each value v plug in D(v1/n ). Execution time now is divided by approximately n since the degree is divided by n and the numerical coefficients do not change much (they shrink a bit). We call this the inner exponent method. Dixon and Macaulay again show the same pattern. The Dixon matrix is 4 × 4. The Macaulay matrix is 84 × 84 2 and has rank 80. Both determinants compute in less than .01 s. The Dixon answer is x1 4 (x1 4 − 1) and the Macaulay 2 is x1 6 (x1 4 − 1) . The latter has a spurious factor of x1 2 , but we see the pattern f (x1 n ). When n = 5, the Dixon matrix is 22 × 22. The determinant can be computed in about 0.1 s and is 10
2
x1 20 (x1 5 − 1) (x1 10 + 123x1 5 + 1)
The Macaulay matrix is 1365 × 1365 but has rank 1354. The determinant of a maximal minor can be computed in about 0.7 s and is exactly the same as the Dixon, except x1 is raised to the 35. 4. Conclusions • The EDF method is effective on the Macaulay matrix, as that determinant also exhibits many factors. • Apparently, it is not well known that one can extract a maximal minor from a Macaulay matrix and work with it in the same way as with Dixon. • Exponential patterns show up in the Macaulay matrix as they do with Dixon. This is surprising, and begs a theoretical explanation. • We find that on problems of interest, Macaulay is always slower than Dixon and needs much more storage space. Sometimes the Macaulay matrix is so large it cannot be computed. There would seem to be no practical reason to use Macaulay instead of Dixon. References [1] P. Bikker, On Bezout’s method for computing the resultant. RISC-Linz Report Series, Johannes Kepler University, Linz, Austria, 1995. [2] L. Buse, M. Elkadi, B. Mourrain, Generalized resultants over unirational algebraic varieties, Journal of Symbolic Computations 29 (2000) 515–526. [3] D. Cox, J. Little, D. O’Shea, Using Algebraic Geometry. Graduate Texts in Mathematics, 185, Springer-Verlag, New York, 1998. [4] A.L. Dixon, The eliminant of three quantics in two independent variables, Proceedings of the London Mathematics Society 6 (1908) 468–478.
1152
R.H. Lewis / Mathematics and Computers in Simulation 80 (2010) 1146–1152
[5] I.Z. Emiris, J.F. Canny, Efficient incremental algorithms for the sparse resultant and the mixed volume, Journal of Symbolic Computations 20 (2) (1995) 117–149. [6] T.F. Havel, Some examples of the use of distances as coordinates for Euclidean geometry, Journal of Symbolic Computations 11 (1991) 579–593. [7] D. Kapur, T. Saxena, L. Yang, Proc. of the International Symposium on Symbolic and Algebraic Computation Algebraic and geometric reasoning using Dixon resultants, A.C.M. Press, 1994. [8] I.S. Kotsireas, K. Karamanos, Exact computation of the bifurcation point B4 of the logistic map and the Bailey–Broadhurst conjectures, International Journal of Bifurcation and Chaos 14 (7) (2004) 2417–2423. [9] R.H. Lewis, Heuristics to accelerate the Dixon resultant, Mathematics and Computers in Simulation 77 (4) (2008) 400–407. [10] R.H. Lewis, Computer Algebra System Fermat. http://home.bway.net/lewis/ (2008). [11] R.H. Lewis, Using the Dixon resultant on big problems, in: CBMS, Conference, Texas, A.M., University, May 2002. http://www.math.tamu.edu/conferences/cbms/abs.html. [12] R.H. Lewis, Exploiting Symmetry in a Polynomial System with the Dixon Resultant, International Conference on Applications of Computer Algebra (ACA), Lamar University, Texas, 2004, July. [13] R.H. Lewis, S. Bridgett, Conic tangency equations arising from Apollonius problems in biochemistry, Mathematics and Computers in Simulation 61 (2) (2003) 101–114. [14] R.H. Lewis, E.A. Coutsias, Algorithmic search for flexibility using resultants of polynomial systems, in: F. Botana, T. Recio (Eds.), Automated Deduction in Geometry: ADG 2006, Pontevedra, Spain, vol. 4869 of Lecture Notes in Computer Science, Springer, Berlin/Heidelberg, 2007. [15] M. Minimair, Macaulay Resultant Package. http://minimair.org/software/ (2007). [16] D. Robertz, V. Gerdt, comparison of software systems, http://home.bway.net/lewis/fermat/gcdcomp (2005). [17] A. Sommese, J. Verschelde, C. Wampler, Numerical decomposition of the solution sets of polynomial systems into irreducible components, SINUM 38 (6) (2001) 2022–2046.