Awange J.,Grafarend E., Palancz B., Zaletnyik P. Algebraic Geodesy and Geoinformatics

Подождите немного. Документ загружается.

5-3 Multipolynomial resultants 55

5-33 The Dixon resultant

5-331 Basic concepts

Consider a system of polynomials

F = {f

, f

, ..., f

} (5.16)

where

α∈A

i,α

(5.17)

and

= x

...x

(5.18)

for each i = 0,...,d, and with the polynomial coeﬃcient c

i,α

, which are also some-

times referred to as the parameters of the system. In general, we consider a col-

lection of polynomials in x

,...,x

with coeﬃcients in the ﬁeld of an arbitrary ring

⊂ N

[113]. Elimination theory tells us how to construct a single polynomial of

d homogeneous polynomial equations in d unknowns, such that its vanishing is a

necessary and suﬃcient condition for a given system to have a non-trivial solution.

This single polynomial is called a resultant polynomial [288].

One way to compute the resultant of a given polynomial system is to construct

a matrix with the property that whenever the polynomial system has a solution,

such a matrix has a deﬁcient rank, thereby implying that the determinant of any

maximal minor is a multiple of the resultant. A simple way to construct a resultant

matrix is to use the dialytic method [76], i.e. multiply each polynomial with a ﬁnite

set of monomials, and rewrite the resulting s ystem in matrix form. We call such

matrix the dialytic matrix.

This alone, however, does not guarantee that a matrix so constructed is a re-

sultant matrix. Note that such matrices are usually quite sparse: matrix entries are

either zero or coeﬃcients of the polynomials in the original system. Goo d exam-

ples of re sultant dialytic matrices are Sylvester [377] for the univariate case, and

Macaulay [281], as well as the Newton sparse matrices of Canny and Emiris [102]

for the multivariate case, where they all diﬀer only in the selection of multiplier

monomial sets.

In contrast to dialytic matrices, the Dixon matrix is dense since its entries are

determinants of the coeﬃcients of the polynomials in the original system. It has

the advantage of being an order of magnitude smaller in comparison to a dialytic

matrix, which is important as the computation of the symbolic determinant of

a matrix is sensitive to its size. The Dixon matrix is constructed through the

computation of the Dixon polynomial, which is expressed in matrix form.

Below, the generalized multivariate Dixon formulation for simultaneously elim-

inating many variables from a polynomial system and computing its resultant is

brieﬂy reviewed [114].

Let π

) = σ

...σ

i+1

...x

where i ∈ {0, 1, ..., d} and the σ

s are new

variables; π

) = x

. Now π

is extended to polynomials in a natural way as

56 5 Polynomial resultants

, ..., x

)) = f

(σ

, ...σ

, x

i+1

, ..., x

) (5.19)

Given a polynomial system F = {f

, f

,...f

}, the Dixon polynomial is deﬁned

θ (f

, f

, ...f

) = δ (x

, . . . , x

, σ

, . . . , σ

) =

i=1

− x

Det







) π

) . . . π

)

) π

) . . . π

)

. . . . . .

) π

) . . . π

)







(5.20)

The order in which the original variables in x are replaced by new variables in

σ is signiﬁcant in the sense that the Dixon polynomial computed using two distinct

variable orderings may diﬀer.

A Dixon polynomial δ( x

x, σ

σ) can be written in bilinear form as [114],

δ (x

, ..., x

, σ

, ...σ

) = ΞΘX

ΞΘX

, (5.21)

where Ξ

Ξ = (σ

,...,σ

) and X

X = (x

,...,x

) are row vectors. The k × s matrix

Θ is called the Dixon matrix.

5-332 Formulation of the Dixon resultant

Cayley’s Formulation of B´ezout’s Method

Let us recall here the Cayley’s [104] formulation of Bezout’s method [309] for

solving two polynomial equations in the univariate case. However, this method is

actually due to Euler (e.g., Salmon [352]).

Let us consider two univariate polynomials f (x) and g(x ), and let

deg = max(degree(f), degree(g)) (5.22)

and σ be an auxiliary variable. The quantity

δ(x, σ) =

x − σ

det



f(x) g(x)

f(σ) g(σ)



f(x)g(σ) − f(σ)g(x)

x − σ

(5.23)

is a symmetric polynomial in x and σ of deg-1 which is called the Dixon polynomial

of f and g. Every common zero of f and g is a zero of δ(x,σ) for all values of σ.

Example 5.3. Let us consider two univariate polynomials with parameter π,

f(x) = (x − 1)(x + 3)(x −4) = 12 − 11x − 2x

+ x

g(x) = (x − π)(x + 4) = −4π + 4x − πx + x

(5.24)

The Dixon polynomial is then expressed as

5-3 Multipolynomial resultants 57

δ =

x − σ

det



f(x) g(x)

f(σ) g(σ)



, (5.25)

which on substituting Eq.(5.24) gives

δ =

x − σ

det



12 − 11x −2x

+ x

−4π + 4x − πx + x

12 − 11σ − 2σ

+ σ

−4π + 4σ − πσ + σ



, (5.26)

and evaluating the determinant and simplifying leads to −48 + 56π −12x + 8πx −

4πx

− 12σ + 8πσ + 3xσ − 2πxσ + 4x

σ − πx

σ − 4πσ

+ 4xσ

− πxσ

+ x

Hence, at a common zero, each coeﬃcient of σ

in δ vanishes, leading to

: −4π + 4x − πx + x

= 0

: −12 + 8π + 3x − 2πx + 4x

− πx

= 0

: −48 + 56π − 12x + 8πx − 4πx

= 0.

(5.27)

It is a homogeneous system in variables x

, x

and x





−4π 4 − π 1

−12 + 8π 3 − 2π 4 − π

−48 + 56π −12 + 8π −4π





















(5.28)

This system has non-trivial solutions if and only if its determinant D is zero.

D is called the Dixon resultant of f and g. The matrix of the system M is the

Dixon matrix,

M =





−4π 4 − π 1

−12 + 8π 3 − 2π 4 − π

−48 + 56π −12 + 8π −4π





(5.29)

and its determinant is

D = −480 + 440π + 80π

− 40π

(5.30)

The polynomials have common zeros if D vanishes. Indeed, solving the equation

D = 0 we get

= −3, π

= 1and π

= 4

5-333 Dixon’s generalization of the Cayley-B´ezout method

Dixon [128] generalized Cayley’s approach to Bezout’s method to systems of three

polynomials equations in three unknowns. Let

f(x, y, z) = 0

g(x, y, z) = 0

h(x, y, z) = 0.

(5.31)

Now, the Dixon polynomial is deﬁned by

δ(x, y, z, σ, ξ) =

(x − σ)(y − ξ)

det





f(x, y, z) g(x, y, z) h(x, y, z)

f(σ, y, z) g(σ, y, z) h(σ, y, z)

f(σ, ξ, z) g(σ, ξ, z) h(σ, ξ, z)





(5.32)

58 5 Polynomial resultants

Example 5.4. Let us consider the three-variate polynomials

f = x

+ y

− 1

g = x

+ z

− 1

h = y

+ z

− 1.

(5.33)

The Dixon polynomial is

δ(x, y, z, σ, ξ) =

(x − σ)(y − ξ)

det





+ y

− 1 x

+ z

− 1 y

+ z

− 1

+ y

− 1 σ

+ z

− 1 y

+ z

− 1

+ ξ

− 1 σ

+ z

− 1 ξ

+ z

− 1





(5.34)

Equation (5.34) leads to yz − 2x

yz + yξ − 2x

yξ + zσ − 2x

zσ + ξσ − 2x

ξσ

Now we eliminate the variables y and z. Considering that

δ(x, y, z, σ, ξ) = σξ



1 − 2x



+ξy



1 − 2x



+σz



1 − 2x



+yz



1 − 2x



, (5.35)

and the system of equations to be

. y

: 0 0 0 1 − 2x

: 0 0 1 − 2x

: 0 1 − 2x

0 0

: 1 − 2x

0 0 0

, (5.36)

The Dixon matrix is then given as

M =







0 0 0 1 − 2x

0 0 1 − 2x

0 1 − 2x

0 0

1 − 2x

0 0 0







, (5.37)

and its determinant, the Dixon resultant, as D = 1 − 8x

+ 24x

− 32x

+ 16x

Dixon (1908) proved that for the three polynomials of degree two, the vanishing

of D is a necessary condition for the existence of a common zero. Furthermore, D

is not identically zero, namely, it is not zero for e very value of x.

5-334 Improved Dixon resultant - Kapur-Saxena-Yang method

Dixon’s method and proofs easily generalize to a system of n+1 generic n degree

polynomials in n unknowns (for more details, see, e.g., [239]). Recall that a poly-

nomial is generic if all its coeﬃcients are independent parameters, unrelated to

each other. A polynomial in n variables is n-degree if all powers to the maximum

of each variable appear in it. Dixon’s method only applies to generic n degree

polynomials. If this condition fails, then one can face the following problems:

5-3 Multipolynomial resultants 59

(a) The Dixon matrix, M, may be singular;

(b) After removal of the rows and columns containing zeros, the vanishing of the

determinant of the Dixon matrix, D, may not give a necessary condition for

the existence of a common zero;

not even be square. Hence, its determinant cannot be deﬁned.

Kapur et al. [239] addressed all three problems successfully, provided a certain

precondition holds. Namely, assuming that the column that corresponds to the

monomial 1 = x

. . ., the Dixon matrix is not a linear combination of the re-

maining ones. If this precondition is true, then D = 0 is a necessary condition of

common zeros. This theorem yields a simple algorithm for obtaining the necessary

condition D = 0, which is called the Kapur-Saxena-Yang-Dixon resultant [308].

The algorithm involves the following:

(1) Compute the Dixon matrix M. If the precondition holds, continue;

(2) Reduce the M rows without scaling to row echelon form, M’.

(3) Compute the product D of the pivots of M’.

This algorithm was implemented into Mathematica by Nakos and Williams [309,

308]. We e mploy this to illustrate how the last example could be solved. In Math-

ematica 5.2 [386], the Dixon resultant prompt is called by typing

<< ResultantDixon.

In solving Eq.(5.34), we ﬁrst compute the Dixon polynomial using

DixonP olynomial[{f, g, h}, {y, z}, {σ, ξ}],

followed by the Dixon matrix using

M = DixonMatrix [{f, g, h}, {y, z}, {σ, ξ}]

which leads to







0 0 0 1 − 2x

0 0 1 − 2x

0 1 − 2x

0 0

1 − 2x

0 0 0







. (5.38)

The Dixon resultant is then computed via DixonResultant[{f, g, h}, {y, z}, {σ, ξ}]

which gives 1 − 8x

+ 24x

− 32x

+ 16x

. Since all these functions are in-built,

one can proceed to compute the Dixon resultant directly without going through

Dixon polynomial and matrix.

5-335 Heuristic methods to accelerate the Dixon resultant

The basic idea of the Dixon method is to construct a square matrix M whose

determinant D is a multiple of the resultant. Usually M is not unique, it is obtained

60 5 Polynomial resultants

as a maximal minor, in a larger matrix we shall call M

, and there are usually

many maximal minors - any one of which will do. The entries in M are polynomials

in parameters. The factors of D that are not the res ultant are called the spurious

factors, and their product is sometimes referred to as the spurious factor.

The naive way to proceed is to compute D, factor it, and separate the spurious

factor from the actual resultant. But there are problems. On one the hand, the

determinant may be so large as for it to be impractical or even impossible to

compute, even though the resultant is relatively small, the spurious factor is huge.

On the other hand, the determinant may be so large that factoring it is impractical.

Lewis developed three heuristic methods to overcome these problems [267, 268].

The ﬁrst may be used on any polynomial system. It uses known factors of D to

compute other factors. The second also may be used on any p olynomial system and

it discovers factors of D so that the complete determinant is never produced. The

third applies only when the resultant appears as a factor of D in certain exponential

patterns. These methods were discovered by experimentation and may apply to

other resultant formulations, such as the Macaulay.

5-336 Early discovery of factors: the EDF method

This method is to exploit the observed fact that D has many factors. In other

words, we try to turn the existence of spurious factors to our advantage. By ele-

mentary row and column manipulations (Gaussian elimination) we discover prob-

able factors of D and extract them from M

≡ M . This produces a smaller matrix

, still with polynomial entries, and a list of discovered numerators and denom-

inators. Here is a very simple example.

Example 5.5. Given initially



9 2

4 4



numerators : denominators :

We factor a 2 out of the second column, then a 2 from the second row. Thus



9 1

2 1



numerators : 2, 2 denominators :

Note that 9 × 4 − 2 ×4 = 2 × 2 × (9 ×1 −2 × 1).

We change the second row by subtracting 2/9 of the ﬁrst



9 1

0 7/9



numerators : 2, 2 denominators :

We pull out the denominator 9 from the second row, and factor out 9 from the

ﬁrst column:



1 1

0 7



numerators : 2, 2, 9 denominators : 9

5-4 Concluding remarks 61

Note that 9 × 7/9 − 1 ×0 =

2×2×9

× (1 × 7 −1 ×0).

We “clean up” or consolidate by dividing out the common factor of 9 from the

numerator and denominator lists; any one that occurs may be erased and the list

compacted since the ﬁrst column is canonically simple. We have hence ﬁnished one

step of the algorithm, and have produced a smaller M

= (7) numerators : 2, 2 denominators :

The algorithm terminates by pulling out the 7:

numerators : 2, 2, 7 denominators :

Realize, however, that the det(M

) = det



9 2

4 4



= 2 ∗ 2 ∗ 7 = 28

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.

The accelerated Dixon resultant by the Early Discovery Factors (Dixon - EDF)

algorithm, was suggested and implemented in the computer algebra system Fermat

by Lewis [267, 269].

The Dixon resultant is a very attractive tool for solving system of multivariate

polynomial geodetic equations (see [319]). Com paring it to other multipolynomial

resultants like Strumfels’s method, it has the advantages of (i) its small size and

high density of the Dixon matrix, (ii) faster computational speed and (iii) being

robust.

5-4 Concluding remarks

With modern computers, the polynomial resultant approaches discussed can eas-

ily be used to develop algorithms for solving systems of nonlinear equations.

Compared to Groebner basis, these have the advantage of not computing ex-

tra parameters, thus requiring less computer memory. Its shortcoming, how-

ever, lies in the formation of the design matrix which become more compli-

cated and cumbersome as the number of polynomials and variables increases.

Unless Groebner basis fails, we recommend it for solving geodetic and geoinfor-

matics nonlinear systems of equations. On the other hand, the polynomial re-

sultants approach comes in handy when the computer’s space is limited. With

modern computer storage capacity though, most problems requiring algebraic

solutions in the ﬁelds mentioned above can easily be handled by Groebner

basis without fear of a computer breakdown. Publications on the subject in-

clude: [22, 28, 34, 49, 100, 101, 103, 118, 128, 148, 149, 195, 254, 279, 280, 281,

282, 284, 285, 286, 287, 288, 352, 289, 290, 291, 292, 305, 373, 375, 401, 404].

Besides Groebner bases and polynomial resultants techniques, there exists an-

other approach for eliminating variables developed by W. WU [418] using the

62 5 Polynomial resultants

ideas proposed by [347]. This approach is based on Ritt’s characteristic set con-

struction and was successfully applied to automated geometric theorem by Wu.

This algorithm is referred as the Ritt-Wu’s algorithm [292].

6 Linear homotpy

6-1 Introductory remarks

A fundamental task in geodesy is the solving of systems of equations. Many geode-

tic problems are represented as systems of multivariate polynomials. A common

problem in solving such systems is improper initial starting values for iterative

methods, leading to the convergence to solutions with no physical meaning, or con-

vergence that requires global method. Although symbolic methods such as Groeb-

ner bases or resultants have been shown to be very eﬃcient, i.e., providing solutions

for determined systems such as 3-point problem of 3D aﬃne transformation, the

symbolic algebra can be very time consuming, even with special Computer Algebra

Systems (CAS). This Chapter proposes the Linear Homotopy method that can be

implemented easily in high level computer languages like C++ and Fortran, which

are faster than CAS by at least two orders of magnitude. Using Mathematica, the

power of Homotopy is demonstrated by solving three nonlinear geodetic problems:

resection, GPS positioning and aﬃne transformation. The method enlarging the

domain of convergence is found to be eﬃcient, less sensitive to rounding errors, and

has a lower complexity compared to other local methods like Newton-Raphson.

The Chapter is organized as follows: In Section 6-2, a background to linear

homotopy is presented. Section 6-3 presents the deﬁnition and basic concepts of

homotopy required to understand the s olution of nonlinear equations. For simplic-

ity, a quadratic and a third degree polynomial equations are used to illustrate the

approach. Although the suggested technique can be applied to square systems,

overdetermined problems can also be solved after transforming the N-point prob-

lem into a square one in a least squares sense. In part II of the book, it will be

demonstrated how the technique can be used to solve three nonlinear geodetic

problems.

6-2 Background to homotopy

Solving nonlinear s ystems , especially algebraic polynomial systems, is a fundamen-

tal problem that occurs frequently in various ﬁelds of science and engineering, like

robotics, computational chemistry, computer vision, computational geometry, and

signal processing. Using numerical algorithms to solve polynomial systems with

tools originating from algebraic geometry is the main activity in so-called Numer-

ical Algebraic Geometry, see e.g., [366, 367]. This is a new and developing ﬁeld

on the crossroads of algebraic geometry, numerical analysis, computer science and

engineering.

The homotopy continuation method is a global numerical method for solving

nonlinear systems in general, and also polynomial systems, see [193]. This method

is used to locate all geometrically isolated roots as well as decompose positive di-

mensional solution sets into irreducible comp onents [394]. Homotopy continuation

J.L. Awange et al., Algebraic Geodesy and Geoinformatics, 2nd ed.,

DOI 10.1007/978-3-642-12124-1 6,

 Springer-Verlag Berlin Heidelberg 2010

64 6 Linear homotpy

has been established as a reliable and eﬃcient method for solving polynomial sys-

tems over the last two decades, originating from the works of [129, 146]. Some com-

puter codes have been developed for homotopy methods, for example Bertini [59],

HOM4PS [263], HOMPACK [402], PHCpack [393], and PHoM [200]. Other codes

are written in Maple (see e.g., [270]) and Matlab [127]. Recently, codes for ﬁxed

point and Newton homotopy were developed in Mathematica see [80, 81].

In geodesy, several algebraic procedures have been put forward for solving

nonlinear systems of equations (see, e.g., [41, 320]). The procedures suggested

in these studies, such as Groebner bases and resultant approaches, are, however,

normally restricted by the size of the nonlinear systems of equations involved.

In most cases, the symbolic computations are time consuming. These symbolic

computations were necessitated by the failure to obtain suitable starting values

for numerical iterative procedures. In situations where large systems of equations

are to be solved (e.g., aﬃne transformation), and where symbolic methods are

insuﬃcient, there exists the need to investigate the suitability of other alternatives.

One such alternative is linear homotopy.

In this study linear homotopy continuation methods are introduced (see [249]).

The deﬁnitions and basic ideas are considered and illustrated. The general algo-

rithms of the linear homotopy m ethod are considered: iterated solution of homo-

topy equations via Newton-Raphson method; and as an initial value problem of

a system of ordinary diﬀerential equations. The eﬃc iency of these methods are

illustrated here by solving polynomial equation systems in geodesy, namely the

solution of 3D resection, GPS navigation and 3D aﬃne transformation problems.

Our computations were carried out with Mathematica and Fermat computer alge-

bra systems on a HP workstation xw 4100 with XP operation system, 3 GHz P4

Intel processor and 1 GB RAM. The details can be found in MathSource (Wolfram

Research Inc.) in [320], while the mathematical background of the algorithms is

discussed in this chapter, and in the literature, e.g. [367]. The application of the

implemented Mathematica functions is illustrated in the Appendix A-3.

6-3 Deﬁnition and basic concepts

The continuous deformation of an object to another object is known as homotopy.

Let us consider two univariate polynomials p

and q

of the same degree,

(x) = x

− 3 (6.1)

(x) = −x

− x + 1. (6.2)

The linear convex function, H, expressed by the variables x and λ, give s the

homotopy function deﬁned as

H(x, λ) = (1 − λ)p

(x) + λq

(x) , (6.3)

where λ ∈ [0, 1].