Keywords

These keywords were added by machine and not by the authors. This process is experimental and the keywords may be updated as the learning algorithm improves.

15.1 Introductory Remarks

Algebraic techniques ofGroebner bases andMultipolynomial resultants are presented in this Chapter as efficient tools for solving explicitly systems of linear and nonlinear equations. Similar to the Gauss elimination technique applied to linear systems of equations,Groebner bases andMultipolynomial resultants are useful in eliminating several variables in multivariate systems of nonlinear equations in such a manner that the end product results intounivariate polynomial equations whose roots are readily determined using existing functions such as theroots command in MATLAB.

The capability ofGroebner bases andMultipolynomial resultants to solve explicitly nonlinear geodetic problems enables them to be used as the computational engine in theGauss–Jacobi combinatorial algorithm proposed byC.F. Gauss (published posthumously, e.g.,Appendix D-3) andC.G.I. Jacobi (1841) to solve thenonlinear Gauss–Markov model. By converting nonlinear geodetic observation equations into algebraic (polynomial) form, theGauss–Jacobi combinatorial algorithm solves an overdetermined system of equations in two steps as follows:

  1. 1.

    In the first step, all them combinations of minimal subsets of the observation equations are formed, and then unknowns of each subset rigorously solved by means of Groebner bases or Multipolynomial resultants. Them solution sets are then represented as points in ann-dimensional space \({\mathbb{R}}^{n}\).

  2. 2.

    In a second step, the solution of the overdetermined Gauß-Markov-model is constructed as a weighted arithmetic mean of thosem solution points, where the weights are obtained from the Error propagation law/variance-covariance propagation.

This chapter presents the theory of Groebner basis, multipolynomial resultants and the Gauss–Jacobi combinatorials methods, and is organized as follows. Section 15.2 provides the background to the algebraic and numerical methods used to derive exact solutions. In Sect. 15.3, the algebraic algorithms needed to solve the nonlinear system of equations are presented, and enriched with examples to help the reader understand them. The algorithms are finally applied successfully to three real-time geodetic problems in Sect. 15.4; transforming in a closed form geocentric coordinates to Gauss ellipsoidal coordinates (geodetic) through Minimum Distance Mapping (MDM); solution of GNSS pseudorange equations (using the GPS satellite vehicle as example); and to obtain the seven datum transformation parameters from two sets of coordinates. By means ofGroebner basis, thescale parameter (in the seven parameters datum transformation problem) is shown to fulfill be aunivariate algebraic equation of fourth order (quartic polynomial), while the rotation parameters are functions in scale and the coordinates differences. What we present in this chapter is a nutshell of algebraic treatment of systems of polynomial equations. For a complete coverage, we refer toJ.L. Awange and E.W. Grafarend (2005), andJ.L. Awange et al. (2010).

15.2 Background to Algebraic Solutions

In Geodesy, Photogrammetry and Computer Vision,nonlinear equations are often encountered in several applications as they often relate the observations (measurements) to the unknown parameters to be determined. In case the number of observationsn and the number of unknownsm are equal (n = m), the unknown parameters may be obtained by solving explicitly or in a closed form the system of equations relating observations to the unknown parameters. For example,D. Cox et al. (1998, pp. 28–32) illustrated that for systems of equations with exact solution, the system becomes vulnerable to small errors introduced during root findings and in case of extending the partial solution to the complete solution of the system, the errors may accumulate and thus become so large. If the partial solution is derived by iterative procedures, then the errors incurred during the root-finding may blow up during the extension of the partial solution to the complete solution (back substitution).

In some applications, symbolic rather than numerical solution are desired. In such cases, explicit procedures are usually employed. The resulting symbolic expressions often consists of univariate polynomials relating the unknown parameters (unknown variables) to the known variables (observations). By inserting known values into these univariate polynomials, numerical solutions are readily computed for the unknown variables. Advantages of explicit solutions are listed, e.g., byE. L. Merritt (1949) as;provision of satisfaction to the users (Photogrammetrist and Mathematicians) of the methods,provision of data tools for checking the iterative methods,desired by Geodesist whose task of control densification does not favor iterative procedures,provision of solace andthe requirement of explicit solutions rather thaniterativeby some applications. In Geodesy for example, theMinimum Distance Mapping problem considered byE. Grafarend andP. Lohse (1991) relates a point on the Earth’s topographical surface uniquely (one-to-one) to a point on theInternational Reference Ellipsoid. The solution of such an optimization problem requires that the equations be solved explicitly.

The draw back that was experienced with explicit solutions was that they were like rare jewels. The reason for this was partly because the methods required extensive computations for the results to be obtained, and partly because the resulting symbolic expressions were too large and required computers with large storage capacity. Until recently, the computers that were available could hardly handle large computations due to lack of faster Central Processing Unit (CPU), shortage of Random Access Memory (RAM) and limited hard disk space for storage. The other set back experienced by the explicit procedures was that some of the methods, especially those from algebraic fields, were formulated based on theoretical concepts that were hard to realize or comprehend without the help of computers. For a long time, these setbacks hampered progress of the explicit procedures. The advances made in computer technology in recent years, however, has changed the tides and led to improvements in explicit computational procedures which hitherto were difficult to achieve. Apart from the improvements in existing computational procedures, new computational techniques are continuously being added to the increasing list of computational methods with an aim of optimizing computational speed and effeiciency (see, e.g., Palancz and Awange 2012). In this category are the algebraic methods ofGroebner bases andMultipolynomial resultants.

This Chapter presents algebraic computational techniques ofGroebner bases andMultipolynomial resultants as suitable tools for solving explicitly nonlinear systems of observation equations that have been converted into algebraic (polynomial) forms in Geodesy and Geoinformatics. The algebraic techniques ofGroebner bases andSylvester resultant (resultant of two polynomials) for solving polynomial equations in Geodesy have been mentioned and examples of their applications to the two-dimensional case given in the work ofP. Lohse (1994, pp. 36–39, 71–76). We will demonstrate using examples how theGroebner bases andMultipolynomial resultants (resultant of more than two polynomials) can be used to solve problems of three-dimensional nature inherent in Geodesy.E. Grafarend (1989) already suggested the use ofGroebner bases approach in solving the perspective center problem in photogrammetry.

Other than revolutionizing computation procedures, the advances made in computer technology have also led to improvement in instrumentation for data acquisition as exemplified in the case of GPS positioning satellites. Since its inception as a positioning tool, theGlobal Positioning System (GPS) – referred to byE. Grafarend andJ. Shan (1996) as theGlobal Problem Solver – has revolutionized geodetic positioning techniques and maintained its supremacy as a positioning tool. These improvements on instrumentation for data acquisition have led to improved data collection procedures and increase in accuracy. Lengthy geodetic procedures such as triangulation that required a lot of time, intervisibility between stations, and a large manpower are rapidly being replaced by satellite positioning techniques, which require shorter observation periods, no intervisibility requirement, weather independent and less manpower leading to optimization of time and money. Such improvements can be seen by the number of Global Navigation Satellite Systems (GNSS) being launched by various countries. These includes; Chinese’s Compass/Beidou and Russian’s GLONASS (see, e.g., Wellenhof 2008; Awange 2012). Whereas improved instrumentation is applauded, it comes a long with its own difficulties. One of the difficulties is that a lot of data is collected than required to determine unknown parameters leading to redundancies. In positioning with GPS for example, due to its constellation that offer a wider coverage, more than four satellites can be observed at any point of the earth. In the minimal case, only four satellites are required to fix the receiver position and the receiver clock bias assuming that the transmitter clock bias and the transmitter-receiver bias have been neglected for thepseudo-range type of observations. More than four satellites as is nowadays the case, therefore, leads to superfluous observations. In such cases, wheren > m, the explicit solutions give way to optimization procedures such asleast squares solution which work very well forlinear models under specified conditions discussed in other chapters of this book.

In Geodesy and Geoinformatics, however, the observation equations are normallynonlinear thus requiring the use ofnonlinear Gauss–Markov model, which is normally solved either by firstlinearizing the observation equations using Taylor series expansion to the second order terms about approximate values of the unknowns then applying linear model estimation procedures, or by usingiterativeprocedures such as theGauss–Newton approach. Thelinearization approach has the disadvantage that the linearized approximation of the nonlinear models may still suffer from nonlinearity and thus resulting in the estimates of such models being far from the real estimates of thenonlinear models. This can easily be checked by re-substituting the estimates from thelinearized model into the original nonlinear model.

For theiterative procedures, the greatest undoing may be the requirement of the initial approximate values to start off the iterations, which may not be available for some applications. For simper models, the approximate initial values may be computed, for others however, the approximate values may be impossible to compute. Apart from the problem of getting the initial approximate values, there also exists the problem that poor choice of approximate values may lead to lack of convergence, or if the approximate values are far from the real solutions, then a large number of iterations may be required to get close to the real solutions thus rendering the whole procedure to be quite slow, especially where multiple roots are available. For other procedures, such as the7-parameter datum transformation that requireslinearization anditerative methods, it is not feasible to take into account the stochasticity of both systems involved.

Clearly, a procedure for solvingnonlinear Gauss–Markov model that can avoid the requirement of initial approximate starting values foriteration andlinearization approaches, and also take into consideration thestochasticity of the systems involved is the desire of modern day geodesist and geoinformatist (see, e.g.,J. L. Awange andE. W. Grafarend 2005;J. L. Awange et al. 2010). With this background, this chapter aims at answering the following questions:

  • For geodetic, geoinformatics and Earth Science related problems requiring explicit solutions, can the algebraic tools ofGroebner bases andMultipolynomial resultants that have found applications in other fields such as Robotics (for kinematic modelling of robots), Visions, Computer Aided Design (CAD), Engineering (offset surface construction in solid modelling), Computer Science (automated theorem proving) be used to solvesystems of nonlinear observation equations of algebraic (polynomial) type?

  • Is there any alternative for solving thenonlinear Gauss–Markov model without resorting tolinearization oriterative procedures that require approximate starting values?

To answer the first question, the chapter uses theGroebner bases andMultipolynomial resultants to solve explicitly three geodetic problems of;Minimum Distance Mapping,GPS four-point pseudo-ranging, and the7-parameter transformation problem in Sect. 15.4. The answer to the second problem becomes clear once the first question has been answered. Should the algebraic techniques ofGroebner bases orMultipolynomial resultants be successful in solving explicitly the selected geodetic problems, they are then used as the computational engine of thecombinatorial algorithm that was first suggested byC.F. Gauss and later on byC. G. I. Jacobi (1841) and extended byP. Werkmeister (1920). We refer to this algorithm simply as theGauss–Jacobi combinatorial algorithm. In attempting to answer the questions above, the objectives of this chapter are:

  1. (1)

    To help you, the reader, understand the algebraic computational procedures of typeGroebner bases andMultipolynomial resultants as applied to solve explicitly (in closed form)geodetic nonlinear problems. In this respect, theGroebner bases andMultipolynomial resultants techniques are used to solve explicitly (symbolically) geodetic problems ofGPS pseudo-ranging four-pointP4P,Minimum Distance Mapping and the7-parameter transformation. By converting thenonlinear observation equations of these selected geodetic problems intoalgebraic (polynomial),Groebner bases andMultipolynomial resultants techniques are used to eliminate several variables in a multivariate systems of nonlinear polynomial equations in such a manner that the end products from the initial systems of nonlinear observation equations result inunivariate polynomials. The elimination procedure is similar to the Gauss elimination approach in linear systems.

  2. (2)

    To introduce to you, the reader, an algebraic approach of Gauss–Jacobian combinatorial for solving overdetermined problems that are traditionally solved using least squares, which peg their operations on assumptions of linearity and prior information about the unknown values. From the principle of weighted arithmetic mean and using theGauss–Jacobi combinatorial lemma (C. G. I. Jacobi 1841), an adjustment procedure that neitherlinearizes thenonlinear observation equations nor usesiterative procedures to solve thenonlinear Gauss–Markov model is developed. Linearization is permitted only fornonlinear error propagation/variance-covariance propagation. Such a procedure is only realized, thanks to theunivariate polynomial generated by the algebraic computational procedures of typeGroebner bases orMultipolynomial resultants as the computing engine for itsminimal combinatorial sets.

15.3 Algebraic Methods for Solving Nonlinear Systems of Equations

In the present Section, we depart from the traditional iterative procedures for estimating the unknown fixed parameters of thenonlinear Gauss–Markov model and present acombinatorial approach that traces its roots back to the work ofC.F. Gauss.C.F. Gauss first proposed the combinatorial approach using the products of squared distances (from unknown point to known points) and the square of the perpendicular distances from the sides of the error triangle to the unknown point as the weights, see, e.g., Appendix D3. According toW.K. Nicholson (1999, pp. 272–273), the motto in Gauss seal read“pauca des matura” meaningfew but ripe. This belief ledC. F. Gauss not to publish most of his important contributions. For instance,W. K. Nicholson (1999, pp. 272–273) writes “Although not all his results were recorded in the diary (many were set down only in letters to friends), several entries would have each given fame to their author if published. Gauss new about the quartenions before Hamilton…”. The combinatorial method, like many of his works, was later to be published after his death. Several years later, the method was to be developed further byC. G. I. Jacobi (1841) who used the square of the determinants as the weights in estimating the unknown parameters from the arithmetic mean.P. Werkmeister (1920) later established the relationship between the area of the error figure formed by the combinatorial approach and the standard error of the determined point. We will refer to thiscombinatorial approach as theGauss–Jacobi combinatorial algorithm needed to solve an overdetermined system of equations algebraically.

While the solution of thelinear Gauss–Markov model by B est L inear U niformly U nbiased E stimator (BLUUE) is straight forward, the solution of thenonlinear Gauss–Markov model has no straight forward procedure owing to thenonlinearity of theinjective function (or map function) that maps \({\mathbb{R}}^{m}\) to \({\mathbb{R}}^{n}\). TheGauss–Jacobi combinatorial algorithm is presented as a possible algebraic solution to thenonlinear Gauss–Markov model. In Sect. 15.3.1, we demonstrate how the procedure can be used to solve thenonlinear Gauss–Markov model.

15.3.1 Solution of Nonlinear Gauss–Markov Model

Thenonlinear Gauss–Markov model is solved in two steps:

  • Step 1: Combinatorial minimal subsets of observations are constructed and rigorously solved by means of theMultipolynomial resultant orGroebner basis (J.L . Awange andE.W. Grafarend 2005;J. L . Awange et al. 2010).

  • Step 2: The combinatorial solution points of step 1 are reduced to their final adjusted values by means of an adjustment procedure where the Best Linear Uniformly Unbiased Estimator (BLUUE) is used to estimate the vector of fixed parameters within the linear Gauss–Markov model with the dispersion matrix of the real valued random vector of pseudo-observations fromStep 1 generated via thenonlinear error propagation law also known in this case as thenonlinear variance-covariance propagation.

15.3.1.1 Construction and Solution of the Combinatorial Subsets

Sincen > m we construct a minimal combinatorial subsets comprisingm equations solvable in closed form using eitherGroebner bases orMultipolynomial resultants which we present in Sects. 15.3.1.2 and 15.3.1.3. We begin by the following elementary definitions:

Definition 15.1 (Permutation). 

Given a setS with elements {i, j, k} ∈ S, the arrangement obtained by placing {i, j, k} ∈ S in some sequence is calledpermutation. If we choose any of the elements sayi first, then each of the remaining elementsj, k can be put in the second position, while the third position is occupied by the unused letter eitherj ork. For the setS, the following permutations can be made:

$$\left.\begin{array}{ccc} \mathit{ijk}&\mathit{ikj}&\mathit{jik}\\ \mathit{jki } &\mathit{kij } & \mathit{kji } \end{array} \right ]\!.$$
(15.1)

From (15.1), there exist three ways of filling the first position, two ways of filling the second position, and one way of filling the third position. Thus the number ofpermutations is given by 3 ×2 ×1 = 6. In general, forn different elements, the number of permutation is equal ton × ×3 ×2 ×1 = n! 

Definition 15.2 (Combination). 

If forn elements onlym elements are used for the permutation, then we have acombination of themth order. If we follow the definition above, then the first position can be filled inn ways, the second inn − 1 ways and themth in \(n - (m - 1)\) ways. In (15.1) above, the combinations are identical and contain the same elements in different sequences. If the arrangement is to be neglected, then we have forn elements, a combination ofmth order being given by

$$\left (\begin{array}{c} n\\ m\end{array} \right ) = \frac{n!} {m!(n - m)!} = \frac{n(n - 1)\ldots (n - m + 1)} {m \times \cdots\times3 \times2 \times1}.$$
(15.2)

Givenn nonlinear equations to be solved, we first form \(\left (\begin{array}{c} n\\ m\end{array} \right )\) minimal combinatorial subsets each consisting ofm elements (wherem is the number of the unknown elements). Each minimal combinatorial subset is solved using the algebraic procedures discussed in Sects. 15.3.1.2 and 15.3.1.3. In geodesy, the number of elementsn normally consist of the observations in the vector \(y\), while the number of elementsm normally consist of the unknown fixed parameters in the vector \(\xi \). In the following Sections, we present two algebraic algorithms that are used to solve the minimal combinatorial subsets (15.2) in closed form. We first present the approach based on theGroebner bases and thereafter consider theMultipolynomial resultants approach.

15.3.1.2 Groebner Basis Method

Groebner basis is the greatest common divisors of a multivariate system of equations (J.L . Awange andE.W. Grafarend 2005;J. L . Awange et al. 2010). Its direct application is the elimination of variables in nonlinear systems of equations. Let us start by the problem of finding the greatest common divisors in Example 15.1:

Given the numbers 12, 20, and 18, find their greatest common divisor. We proceed by writing the factors as

$$\left.\begin{array}{c} 12 = {2}^{2}.{3}^{1}.{5}^{0} \\ 20 = {2}^{2}.{3}^{0}.{5}^{1} \\ 18 = {2}^{1}.{3}^{2}.{5}^{0}\ \end{array} \right ] \rightarrow{2}^{1}.{3}^{0}.{5}^{0} = 2,$$
(15.3)

leading to 2 as the greatest common divisor of 12, 20 and 18. Next, let us consider the case of univariate polynomialsf 1, f 2 ∈ k[x] in (15.4).

$$\left.\begin{array}{l} {f}_{1} = 3{x}^{4} - 3{x}^{3} + 8{x}^{2} + 2x - 5 \\ {f}_{2} = 5{x}^{4} - 4{x}^{2} - 9x + 21\ \end{array} \right ] \rightarrow \mathit{Euclidean\ algorithm} = f \in k[x].$$
(15.4)

Equation (15.4) employs the Euclidean algorithm which obtains one univariate polynomial as the gcd of the two univariate polynomialsf 1 andf 2. If on the other hand expressions in (15.4) were not univariate but multivariate, e.g.,g 1, g 2 ∈ k[x, y] as in (15.5), then one applies the Buchberger algorithm which then computes the Groebner bases. For example,

$$\left.\begin{array}{c} {g}_{1} = xy + x - y - 1 \\ {g}_{2} = xy - x - y + 1\ \end{array} \right ] \rightarrow \mathit{Buchberger}\mathit{algorithm = Groebnerbasis}.$$
(15.5)

Groebner basis therefore, is the greatest common divisors of a multivariate system of polynomial equations {g 1, g 2}.

As stated earlier, it is useful for eliminating variables in nonlinear systems of equations. Gauss elimination technique on the other hand is applicable for linear cases as shown in Example 15.2.

Solve the linear system of equations

$$\left [\begin{array}{c} - x + y + 2z = 2 \\ 3x - y + z = 6 \\ - x + 3y + 4z = 4.\ \end{array} \right.$$
(15.6)

The first step is to eliminatex in the second and third expressions of (15.6). This is achieved by multiplying the first expression by 3 and adding to the second expression to give the second expression of (15.7). The third expression of (15.7) is obtained by subtracting the first expression from the third expression in (15.6).

$$\left [\begin{array}{c} - x + y + 2z = 2 \\ 2y + 7z = 12 \\ 2y + 2z = 2.\ \end{array} \right.$$
(15.7)

The second step is to eliminatey in the second and third expressions of (15.7). This is achieved by subtracting the second expression from the third expression in (15.7) to give (15.8).

$$\left [\begin{array}{c} - x + y + 2z = 2 \\ 2y + 7z = 12 \\ - 5z = -10.\ \end{array} \right.$$
(15.8)

The solution ofz = 2 in (15.8) can now be substituted back into the second equation \(2y + 7z = 12\) to give the value of \(y = -1\), which together with the value ofz = 2 are substituted into the first equation to give the value ofx = 1 to complete the Gauss elimination technique.

In many applications however, equations relating unknown variables to the measured (observed) quantities are normally nonlinear and often consist of many variables (multivariate). In such cases, the Gauss elimination technique for the univariate polynomial equations employed in Example 15.2 gives way to Groebner basis as illustrated in Examples 15.3 and 15.4. In general, the Groebner basis algorithm reduces a system of multivariate polynomial equations. This is done by employing operations “addition” and “multiplication” on a polynomial ring (see, e.g., Awange and Grafarend 2005; Awange et al. 2010) to give more simplified expressions. Given a system of polynomial equations which are to be solved explicitly for unknowns, e.g., (15.5), Groebner basis algorithm is applied to reduce the set of polynomials into another set (e.g., from a systemF(x, y, z) to another systemG(x, y, z)) of polynomials with suitable properties that allow solution. IfF(x, y, z) is a set of nonlinear system of polynomial equations, Groebner basis eliminates variables in a manner similar to Gauss eliminationtechnique for linear cases to reduce it toG(x, y, z). With Lexicographic orderingof the monomials (seeDefinition D1.2 in Appendix D), one expression inG(x, y, z) always turns out to be a univariate polynomial.Its roots are easily obtained using algebraic software of Matlab, Mathematica or Maple (see e.g., Appendix D3),and can be substituted in the other elements of the setG(x, y, z) to obtain a complete solution which also satisfy the original setF(x, y, z). Examples 15.3 and 15.4 elaborate on the application of Groebner basis.

Let us consider a simple example from Buchberger (2001). Consider a setF(x, y) = {f 1, f 2} to have as its elements

$$\left [\begin{array}{l} {f}_{1} = xy - 2y \\ {f}_{2} = 2{y}^{2} - {x}^{2},\ \end{array} \right.$$
(15.9)

where {f 1, f 2} ∈ I are the generators of the IdealI (see definition of Ideal on p. 12). We now seek a simplified set of generators of this Ideal using Buchberger algorithm. By employing operations “addition” and “multiplication”, the Groebner basis algorithm (also called Buchberger algorithm) reduces the system of nonlinear equations (15.9) into another setG ofF as

$$G :=\{ -2{x}^{2} + {x}^{3},-2y + xy,-{x}^{2} + 2{y}^{2}\}.$$
(15.10)

In Mathematica software, using the lexicographic orderx > y, i.e.,x comes beforey, the Groebner basis could simply be computed by entering the command

$$\mathit{Groebner\ Basis}[F,\{x,y\}].$$
(15.11)

The setG in (15.10) contains one univariate polynomial \(-2{x}^{2} + {x}^{3}\), which can easily be solved using roots command in Matlab for solutions \(\{x = 0,x = 0,x = 2\}\) and substituted in any of the remaining elements of the setG to solve fory. The solutions ofG, i.e., the roots \(\{x = 0,x = 0,x = 2\})\) and those ofy satisfy polynomials inF. This can be easily tested by substituting these solutions into (15.9) to give 0. Let us consider as a second example an optimization problem.

Find the minimum and maximum of \(f(x,y,z) = {x}^{3} + 2xyz - {z}^{2}\), such that \(g(x,y,z) = {x}^{2} + {y}^{2} + {z}^{2} - 1\). First, we obtain the partial derivatives of \(f - L\,g = 0\) with respect to {x, y, z, L}, whereL is the Lagrangean multiplier as

$$\frac{\partial f} {\partial \{x,y,z,L\}} := F = \left [\begin{array}{l} 3{x}^{2} + 2yz - 2xL = 0 \\ 2xz - 2yL = 0 \\ 2xy - 2z - 2zL = 0 \\ {x}^{2} + {y}^{2} + {z}^{2} - 1 = 0.\ \end{array} \right.$$
(15.12)

Groebner basis is invoked in Mathematica by

$$\mathit{Groebner\ Basis}[\{F\},\{x,y,L,z\}],$$

which leads to

$$G = \left [\begin{array}{l} L - 1.5x - 1.5yz - 43.7{z}^{6} - 62.2{z}^{4} - 17.5{z}^{2} = 0 \\ {x}^{2} + {y}^{2} + {z}^{2} - 1 = 0 \\ {y}^{2}z - 1.8{z}^{5} + 2.8{z}^{3} - z = 0 \\ {z}^{7} - 1.5{z}^{5} + 0.6{z}^{3} - 0.04z = 0.\ \end{array} \right.$$
(15.13)

The solution ofz in (15.13) can then be substituted into the third equation \({y}^{2}z - 1.8{z}^{5} + 2.8{z}^{3} - z = 0\) to give the value ofy. The obtained values ofz andy are then substituted into the second equation to give the value ofx, and thus complete the Groebner basis solution. Later in the chapter, we will introduce thereduced Groebner basis which can be used to obtain directly the last expression of (15.13), i.e., the univariate polynomial inz. The theory behind the operation of Groebner basis is however not so simple. In the remainder of this chapter, we will try to present in a simplified form the algorithm behind the computation of Groebner bases. The computation of Groebner basis is achieved by the capability to manipulate the polynomials to generateIdeals defined as

Definition 15.3 (Ideal). 

An Ideal is generated by a family of generators as consisting of the set of linear combinations of these generators with polynomial coefficients. Letf 1, , f s andc 1, , c s be polynomials inkx 1, , x n , then

$$\langle {f}_{1},\ldots,{f}_{s}\rangle ={ \sum}_{i=1}^{s}{c}_{ i}{f}_{i}.$$
(15.14)

In (15.14), ⟨f 1, , f s ⟩ is an Ideal and if a subsetI ⊂ kx 1, , x n is an Ideal, it must satisfy the following conditions (Cox 1997, p. 29);

  • 0 ∈ I

  • Iff, g ∈ I, thenf + g ∈ I (i.e.,I is an additive subgroup of the additive group of the fieldk)

  • Iff ∈ I andc ∈ kx 1, , x n , thencf  ∈ I (i.e.,I is closed under multiplication ring element).

Consider equations expressed algebraically as

$$\left [\begin{array}{l} {f}_{1} := {({x}_{1} - {x}_{0})}^{2} + {({y}_{1} - {y}_{0})}^{2} - {d}_{1}^{2} \\ {f}_{2} := {({x}_{2} - {x}_{0})}^{2} + {({y}_{2} - {y}_{0})}^{2} - {d}_{2}^{2},\ \end{array} \right.$$
(15.15)

where polynomials {f 1, f 2} belong to the polynomial ring \(\mathbb{R}[{x}_{0},{y}_{0}]\). If the polynomials

$$\left [\begin{array}{l} {c}_{1} := 4{x}_{0} + 6 \\ {c}_{2} := {x}_{0} + {y}_{0}\ \end{array} \right.$$
(15.16)

also belong to the same polynomial ring \(\mathbb{R}[{x}_{0},{y}_{0}]\), an Ideal is generated by a linear combination

$$I := \left [\begin{array}{c} \langle {f}_{1},{f}_{2}\rangle = {c}_{1}{f}_{1} + {c}_{2}{f}_{2} \\ = (4{x}_{0} + 6){f}_{1} + ({x}_{0} + {y}_{0}){f}_{2}. \end{array} \right.$$
(15.17)

In this case, {f 1, f 2} are said to be generators of theIdealI.Definition (15.3) of anIdeal can be presented in terms of polynomial equationsf 1, , f s  ∈ kx 1, , x n . This is done by expressing the system of polynomial equations as

$$\left [\begin{array}{*{20}c} {f}_{1} = 0 \\ {f}_{2} = 0\\.\\. \\ {f}_{s} = 0,\\ \end{array} \right.$$
(15.18)

and using them to derive others by multiplying each individual equation\({f}_{i}\) by another polynomialc i  ∈ kx 1, , x n and summing to get \({c}_{1}{f}_{1} + {c}_{2}{f}_{2} + \ldots+ {c}_{s}{f}_{s} = 0\) (cf., 15.14). TheIdealf 1, , f s ⟩ thus consists of a system of equations\({f}_{1} = {f}_{2} = \ldots= {f}_{s} = 0\), thus indicating that iff 1 ,…,f s ∈ kx 1 ,…,x n, then⟨f 1 ,…,f s is anIdeal generated byf 1, , f s , i.e., being thebasis of theIdeal \(I\).

In this case, a collection of thesenonlinear algebraic equations formingIdeals are referred to as the set of polynomials generating theIdeal and forms the elements of thisIdeal. Perhaps a curious reader may begin to wonder why the termIdeal is used. To quench this curiosity we refer to Nicholson (1999, p. 220) and quote from Becker (1993, p. 59) who wrote:

Consider example (15.3) with polynomials in \(\mathbb{R}\left [x,y\right ]\). TheIdeal \(I =\langle xy - 2y,2{y}^{2} - {x}^{2}\rangle\). The generators of an Ideal can be computed using thedivision algorithm defined as

Definition 15.4 (Division algorithm). 

Fix a monomial order of polynomials sayx > y for polynomialsF = (h 1, , h s ). Then everyf ∈ kx, y can be written in the form \(f = {a}_{1}{h}_{1} + {a}_{2}{h}_{2} + \ldots+ {a}_{s}{h}_{s} + r,\) wherea i , r ∈ kx, y and eitherr = 0 or a linear combination with coefficients ink of monomials, none of which is divisible by any ofLT(f 1), , LT(f s ) (see Definition D1.5 in Appendix D for leading term LT).

Divide the polynomial \(f = {x}^{3} + 2{x}^{2} + x + 5\) by \(h = {x}^{2} - 2\). We proceed as follows:

$$\left [\begin{array}{ccc} & \underline{x + 2} & \\ {x}^{2} - 2\mid &{x}^{3} + 2{x}^{2} + x + 5& \\ & \underline{{x}^{3} - 2x} & \\ & 2{x}^{2} + 3x + 5 & \\ & \underline{2{x}^{2} - 4} & \\ & 3x + 1, &\ \end{array} \right.$$
(15.19)

implying \({x}^{3} + 2{x}^{2} + x + 5 = (x + 2)({x}^{2} - 2) + (3x + 1)\), with \(a = (x + 2)\) and \(r = (3x + 1)\). Thedivision algorithm given in definition (15.4) fits well to the case ofunivariate polynomials as the remainderr can uniquely be determined. Formultivariate polynomials, the remainder may not be uniquely determined as this depends on the order of the divisors. The division of the polynomialF byf 1, f 2 wheref 1 comes beforef 2 may not necessarily give the same remainder as the division ofF byf 2, f 1 in whose case the order has been changed. This problem is overcome if we pass over to Groebner basis where the existence of everyIdeal is assured by the Hilbert Basis Theorem (Cox 1997, pp. 47–61). The Hilbert Basis Theorem assures that every IdealI ⊂ kx 1, , x n has a finite generating set, that isI = ⟨g 1, , g s ⟩ for some {g 1, , g s } ∈ I. The finite generating setG inHilbert Basis Theorem is what is known as abasis. Suppose every non-zero polynomial is written in decreasing order of its monomials:

$${\sum}_{i=1}^{n}{d}_{ i}{x}_{i},\quad {d}_{i}\neq 0,\,{x}_{i} > {x}_{i+1},$$
(15.20)

if we let the system of generators of theIdeal be in a setG, a polynomialf is reduced with respect toG if no leading monomial of an element ofG (LM (G)) divides the leading monomial off (LM(f)). The polynomialf is said to becompletely reduced with respect toG if no monomials off is divisible by the leading monomial of an element ofG (Davenport 1988, pp. 96–97). ThebasisG, which completely reduces the polynomialf and uniquely determines the remainderr is also known as the Groebner basis and is defined as follows:

Definition 15.5 (Groebner basis). 

A system of generatorsG of an IdealI is called a Groebner basis (with respect to the order < ) if every reduction off ∈ I to a reduced polynomial (with respect toG) always gives zero as a remainder. This definition is a special case of a more general definition given as: Fix a monomial order and letG = g 1, , g t  ⊂ kx 1, , x n ​. Givenf ∈ kx 1, , x n , thenf reduces to zero ModuloG, written as

$$f {\rightarrow }_{G}0,$$
(15.21)

iff can be written in the form (cf., 15.18 on p. 13)

$$f = {a}_{1}{g}_{1} + \ldots+ {a}_{t}{g}_{t},$$
(15.22)

such that whenevera i g i ≠0, we have multideg(f) ≥ multideg(a i g i ) (see Definition D1.5 in Appendix D for leading term LT, LM and Multideg).

FollowingDefinition 15.5, the reader can revisit Examples 15.3 and 15.4 which present the Groebner basisG of the original systemF of equations.

Groebner basis has become a household name in algebraic manipulations and finds application in fields such as mathematics and engineering for solving partial differential equations, e.g.,(Lidl 1998, p. 432). It has also found use as a tool for discovering and proving theorems to solving systems of polynomial equations as elaborated in publications by Buchberger and Winkler (1998). Groebner basis also give a solution to theIdeal membership problem. By reducing a given polynomialf with respect to the Groebner basisG,f is said to be a member of theIdeal if zero remainder is obtained. This implies that ifG = {g 1, , g s } is a Groebner basis of anIdealI ⊂ kx 1, , x n andf ∈ kx 1, , x n a polynomial,f ∈ I if and only if the remainder on division off byG is zero. Groebner bases can also be used to show the equivalence of polynomial equations. Two sets of polynomial equations will generate the sameIdeal if and only if their Groebner bases are equal with respect to any term ordering, e.g., the solutions of (15.10) satisfy those of (15.9). This property is important in that the solutions of the Groebner basis will satisfy the original system formed by the generating set of nonlinear equations. It implies that a system of polynomial equations \({f}_{1}({x}_{1},\ldots,{x}_{n}) = 0,\ldots,{f}_{s}({x}_{1},\ldots,{x}_{n}) = 0\) will have the same solutions with a system arising from any Groebner basis of \({f}_{1},\ldots,{f}_{s}\) with respect to any term ordering. This is the main property of Groebner basis that is used to solve systems of polynomial equations as will be explained in the next section.

15.3.1.2.1 B. Buchberger algorithm

Given polynomialsg 1, , g s  ∈ I, theB. Buchberger algorithm seeks to derive the standard generators or the Groebner basis of thisIdeal. Systems of equations \({g}_{1} = 0,\ldots,{g}_{s} = 0\) to be solved in practise are normally formed by these same polynomials which here generating theIdeal. TheB. Buchberger algorithm computes theGroebner basis by making use of pairs of polynomials from the original polynomialsg 1, , g s  ∈ I and computes the subtraction polynomial known as theS-polynomial explained inD. Cox et al. (1997; p. 81) as follows:

Definition 15.6 (S-polynomial). 

Letf, g ∈ kx 1, …x n be two non-zero polynomials. If multideg (f) = α and multideg (g) = β, then let γ = γ1, , γ n , where γ i  = maxα i , β i for eachi.x γ is called theLeast Common Multiple (LCM) of LM(f) and LM(g) expressed asx γ = LCMLM(f), LM(g). TheS-polynomial off andg is given as

$$S(f,g) = \frac{{x}^{\gamma }} {LT(f)}f - \frac{{x}^{\gamma }} {LT(g)}g.$$
(15.23)

The expression above givesS as a linear combination of the monomials \(\frac{{x}^{\gamma }} {LT(f)}\), \(\frac{{x}^{\gamma }} {LT(g)}\) with polynomial coefficientsf andg and thus belongs to theIdeal generated byf andg.

Definition 15.7 (Groebner basis in terms ofS-polynomial). 

A basisG isGroebner basis if and only if for every pair of polynomialsf andg ofG,S(f, g) reduces to zero with respect toG. More generally a basisG = g 1, , g s for anIdealI is aGroebner basis if and only ifS(f, g) →  G 0, ij.

The implication ofDefinition 15.7 is the following: Given two polynomialsf, g ∈ G such thatLCM(LM(f), LM(g)) = LM(f). LM(g), the leading monomials off andg are relatively prime leading toS(f, g) →  G 0. The concept of prime integer is clearly documented inK. Ireland andM. Rosen (1990 pp. 1–17).

Consider the two polynomials

$$\left [\begin{array}{l} {g}_{1} = {x}_{1}^{2} + 2{a}_{12}{x}_{1}{x}_{2} + {x}_{2}^{2} + {a}_{oo} \\ {g}_{2} = {x}_{2}^{2} + 2{b}_{23}{x}_{2}{x}_{3} + {x}_{3}^{2} + {b}_{oo}.\\ \end{array} \right.$$
(15.24)

TheS-polynomial can then be computed as follows: First we choose alexicographic orderingx 1 > x 2 > x 3 then

$$\left [\begin{array}{c} LM({g}_{1}) = {x}_{1}^{2},\,LM({g}_{2}) = {x}_{2}^{2},LT({g}_{1}) = {x}_{1}^{2},\,LT({g}_{2}) = {x}_{2}^{2} \\ LCM(LM({g}_{1}),LM({g}_{2})) = {x}_{1}^{2}{x}_{2}^{2} \\ S({g}_{1},{g}_{2}) = \frac{{x}_{1}^{2}{x}_{ 2}^{2}} {{x}_{1}^{2}} ({x}_{1}^{2} + 2{a}_{12}{x}_{1}{x}_{2} + {x}_{2}^{2} + {a}_{oo}) -\frac{{x}_{1}^{2}{x}_{ 2}^{2}} {{x}_{2}^{2}} ({x}_{2}^{2} + 2{b}_{23}{x}_{2}{x}_{3} + {x}_{3}^{2} + {b}_{oo}) \\ = {x}_{2}^{2}{x}_{1}^{2} + 2{a}_{12}{x}_{1}{x}_{2}^{3} + {x}_{2}^{4} + {a}_{oo}{x}_{2}^{2} - {x}_{1}^{2}{x}_{2}^{2} - 2{b}_{23}{x}_{1}^{2}{x}_{2}{x}_{3} - {x}_{1}^{2}{x}_{3}^{2} - {b}_{oo}{x}_{1}^{2}) \\ = -{b}_{oo}{x}_{1}^{2} - 2{b}_{23}{x}_{1}^{2}{x}_{2}{x}_{3} - {x}_{1}^{2}{x}_{3}^{2} + 2{a}_{12}{x}_{1}{x}_{2}^{3} + {x}_{2}^{4} + {a}_{oo}{x}_{2}^{2}\\ \end{array} \right.$$
(15.25)

Consider two polynomial equations given as

$$\left [\begin{array}{l} {g}_{3} = {x}_{1}^{2} - 2{a}_{0}{x}_{1} + {x}_{2}^{2} - 2{b}_{0}{x}_{2} + {x}_{3}^{2} - 2{c}_{0}{x}_{3} - {x}_{4}^{2} + 2{d}_{0}{x}_{4} + {a}_{0}^{2} + {b}_{0}^{2} + {c}_{0}^{2} + {d}_{0}^{2} \\ {g}_{4} = {x}_{1}^{2} - 2{a}_{1}{x}_{1} + {x}_{2}^{2} - 2{b}_{1}{x}_{2} + {x}_{3}^{2} - 2{c}_{1}{x}_{3} - {x}_{4}^{2} + 2{d}_{1}{x}_{4} + {a}_{1}^{2} + {b}_{1}^{2} + {c}_{1}^{2} + {d}_{1}^{2}.\\ \end{array}\right.$$
(15.26)

By choosing thelexicographic orderingx 1 > x 2 > x 3 > x 4, theS-polynomial is computed as follows

$$\left [\begin{array}{l} LM({g}_{3}) = {x}_{1}^{2},LM({g}_{4}) = {x}_{1}^{2},LT({g}_{3}) = {x}_{1}^{2},LT({g}_{4}) = {x}_{1}^{2} \\ LCM(LM({g}_{3}),LM({g}_{4})) = {x}_{1}^{2} \\ S({g}_{3},{g}_{4}) = \frac{{x}_{1}^{2}} {{x}_{1}^{2}} ({g}_{3}) -\frac{{x}_{1}^{2}} {{x}_{1}^{2}} ({g}_{4}) = {g}_{3} - {g}_{4} \\ S({g}_{3},{g}_{4}) = 2({a}_{1} - {a}_{o}){x}_{1} + 2({b}_{1} - {b}_{o}){x}_{2} + 2({c}_{1} - {c}_{o}){x}_{3}+ \\ + 2({d}_{o} - {d}_{1}){x}_{4} + {a}_{o} - {a}_{1} + {b}_{o} - {b}_{1} + {c}_{o} - {c}_{1} + {d}_{0} - {d}_{1}\\ \end{array} \right.$$
(15.27)

As an additional example on the computation ofS-polynomial, let us consider the minimum distance mapping problem. The last two polynomials equations are given as

$$\left [\begin{array}{l} {g}_{5} = {x}_{3} + {b}^{2}{x}_{3}{x}_{4} - Z \\ {g}_{6} = {b}^{2}{x}_{1}^{2} + {b}^{2}{x}_{2}^{2} + {a}^{2}{x}_{3}^{2} - {a}^{2}{b}^{2}.\\ \end{array} \right.$$
(15.28)

By choosing thelexicographic orderingx 1 > x 2 > x 3 > x 4, theS-polynomial is computed as follows

$$\left [\begin{array}{l} LM({g}_{5}) = {x}_{3},LM({g}_{6}) = {x}_{1}^{2},LT({g}_{5}) = {x}_{3},LT({g}_{6}) = {b}^{2}{x}_{1}^{2} \\ LCM(LM({g}_{5}),LM({g}_{6})) = {x}_{1}^{2}{x}_{3} \\ S({g}_{5},{g}_{6}) = \frac{{x}_{1}^{2}{x}_{ 3}} {{x}_{3}} ({x}_{3} + {b}^{2}{x}_{3}{x}_{4} - Z) -\frac{{x}_{1}^{2}{x}_{ 3}} {{b}^{2}{x}_{1}^{2}} ({b}^{2}{x}_{1}^{2} + {b}^{2}{x}_{2}^{2} + {a}^{2}{x}_{3}^{2} - {a}^{2}{b}^{2}) \\ S({g}_{5},{g}_{6}) = -Z{x}_{{}_{1}}^{2} + {b}^{2}{x}_{1}^{2}{x}_{3}{x}_{4} - {x}_{2}^{2}{x}_{3} -\frac{{a}^{2}} {{b}^{2}} {x}_{3}^{3} + {a}^{2}{x}_{3}\\ \end{array} \right.$$
(15.29)

By means of an Example given byJ. H. Davenport et al. (1988, pp. 101–102), we illustrate how theB. Buchberger algorithm works. Let us consider theIdeal generated by the polynomial equations

$$\left [\begin{array}{l} {g}_{1} = {x}^{3}yz - x{z}^{2} \\ {g}_{2} = x{y}^{2}z - xyz \\ {g}_{3} = {x}^{2}{y}^{2} - z\\ \end{array} \right.$$
(15.30)

with thelexicographic orderingx > y > z adopted. TheS-polynomials to be considered areS(g 1, g 2),S(g 2, g 3) andS(g 1, g 3). We consider firstS(g 2, g 3) and show that the result is used to suppressg 1 upon which any pairS(g 1, g i ) (e.g.S(g 1, g 2) andS(g 1, g 3)) containingg 1 will not be considered. \(LT({g}_{2}) = x{y}^{2}z,\,\,LT({g}_{3}) = {x}^{2}{y}^{2},\,\,\) thenLCM(g 2, g 3) = x 2 y 2 z respectively

$$\left [\begin{array}{c} S({g}_{2},{g}_{3}) = \frac{{x}^{2}{y}^{2}z} {x{y}^{2}z} {g}_{2} -\frac{{x}^{2}{y}^{2}z} {{x}^{2}{y}^{2}} {g}_{3} \\ = ({x}^{2}{y}^{2}z - {x}^{2}yz) - ({x}^{2}{y}^{2}z - {z}^{2}) \\ = -{x}^{2}yz + {z}^{2}\\ \end{array} \right.$$
(15.31)

We immediately note that the leading term of the resulting polynomial LT(S(g 2, g 3)) is not divisible by any of the leading terms of the elements ofG, and thus the remainder upon the division ofS(g 2, g 3) by the polynomials inG is not zero (i.e. when reduced with respect toG) thusG isnot aGroebnerbasis. This resulting polynomial after the formation of theS-polynomial is denotedg 4 and its negative (to make calculations more reliable) added to the initial set ofG leading to

$$\left [\begin{array}{l} {g}_{1} = {x}^{3}yz - x{z}^{2} \\ {g}_{2} = x{y}^{2}z - xyz \\ {g}_{3} = {x}^{2}{y}^{2} - z \\ {g}_{4} = {x}^{2}yz - {z}^{2}.\\ \end{array} \right.$$
(15.32)

TheS-polynomials to be considered are nowS(g 1, g 2), S(g 1, g 3), S(g 1, g 4), S(g 2, g 4) andS(g 3, g 4). In the set ofG, one can writeg 1 = xg 4 leading without any change inG to the suppression ofg 1 leaving onlyS(g 2, g 4) andS(g 3, g 4) to be considered. Then

$$\left [\begin{array}{c} S({g}_{2},{g}_{4}) = x{g}_{2} - y{g}_{4} \\ = -{x}^{2}yz + y{z}^{2}\\ \end{array} \right.$$
(15.33)

which can be reduced by addingg 4 to give \({g}_{5} = y{z}^{2} - {z}^{2}\), a non zero value thus the setG isnot aGroebnerbasis. This value is added to the set ofG to give

$$\left [\begin{array}{l} {g}_{2} = x{y}^{2}z - xyz, \\ {g}_{3} = {x}^{2}{y}^{2} - z, \\ {g}_{4} = {x}^{2}yz - {z}^{2}, \\ {g}_{5} = y{z}^{2} - {z}^{2},\\ \end{array} \right.$$
(15.34)

theS-polynomials to be considered are nowS(g 3, g 4), S(g 2, g 5), S(g 3, g 5) andS(g 4, g 5). We then compute

$$\left [\begin{array}{c} S({g}_{3},{g}_{4}) = z{g}_{3} - y{g}_{4} \\ = y{z}^{2} - {z}^{2}\\ \end{array} \right.$$
(15.35)
$$\mathrm{whichuponsubtractionfrom}{g}_{5}\mathrm{reducestozero.Further,}$$
$$\left [\begin{array}{c} S({g}_{2},{g}_{5}) = z{g}_{2} - xy{g}_{5} \\ = -xy{z}^{2} + xy{z}^{2} \\ = 0\\ \end{array} \right.$$
(15.36)
$$\mbox{ and}$$
$$\left [\begin{array}{c} S({g}_{4},{g}_{5}) = z{g}_{4} - {x}^{2}y{g}_{5} \\ = {x}^{2}{z}^{2} - {z}^{3},\\ \end{array} \right.$$
(15.37)
$$\mathrm{whichisaddedto}G\mathrm{as}{g}_{6}\mbox{ giving}$$
$$\left [\begin{array}{l} {g}_{2} = x{y}^{2}z - xyz, \\ {g}_{3} = {x}^{2}{y}^{2} - z, \\ {g}_{4} = {x}^{2}yz - {z}^{2}, \\ {g}_{5} = y{z}^{2} - {z}^{2}, \\ {g}_{6} = {x}^{2}{y}^{2} - {z}^{3},\\ \end{array} \right.$$
(15.38)

TheS polynomials to be considered are nowS(g 3, g 5), S(g 2, g 6), S(g 3, g 6), S(g 4, g 6) andS(g 5, g 6). We now illustrate that all thisS-polynomials reduces to zero as follows

$$\left [\begin{array}{l} S({g}_{3},{g}_{5}) = {z}^{2}{g}_{3} - {x}^{2}y{g}_{5} = {x}^{2}y{z}^{2} - {z}^{3} - z{g}_{4} = 0 \\ S({g}_{2},{g}_{6}) = xz{g}_{2} - {y}^{2}{g}_{6} = -{x}^{2}{y}^{2}{z}^{2} + {y}^{2}{z}^{3} + {y}^{2}{g}_{4} = 0 \\ S({g}_{3},{g}_{6}) = {z}^{2}{g}_{3} - {y}^{2}{g}_{6} = {y}^{2}{z}^{3} - {z}^{3} - (yz - z){g}_{5} = 0 \\ S({g}_{4},{g}_{6}) = z{g}_{4} - y{g}_{6} = y{z}^{3} - {z}^{3} - z{g}_{5} = 0 \\ S({g}_{5},{g}_{6}) = {x}^{2}{g}_{5} - y{g}_{6} = -{x}^{2}{z}^{2} + y{z}^{3} + {g}_{6} - z{g}_{5} = 0.\\ \end{array} \right.$$
(15.39)

Thus (15.39) comprise theGroebner basis of the original set in (15.30).

The importance of theS-polynomials is that they lead to the cancelation of the leading terms of the polynomial pairs involved. In so doing the polynomial variables are systematically eliminated according to the polynomial ordering chosen. For example if thelexicographic orderingx > y > z is chosen,x will be eliminated first, followed byy and the final expression may consist only of the variablez.D. Cox et al. (1998, p.15) has indicated the advantage oflexicographic ordering as being the ability to produceGroebner basis with systematic elimination of variables.Graded lexicographic ordering on the other hand has the advantage of minimizing the amount of computational space needed to produce theGroebner basis. The procedure is thus ageneralisation of theGauss elimination procedure for linear systems of equations. If we now put our system of polynomial equations to be solved in a setG,S-pair combinations can be formed from the set ofG as explained in thedefinitions above. Thetheorem, known as theBuchberger’sS-pairpolynomial criterion, gives the criterion for deciding whether a given basis is aGroebner basis or not. It suffices to compute all theS-polynomials and check whether they reduce to zero. Should one of the polynomials not reduce to zero, then the basis fails to be aGroebner basis. Since the reduction is a linear combination of the elements ofG, it can be added to the setG without changing theIdeal generated.B. Buchberger (1979) gives anoptimisation criterion that reduces the number of theS-polynomials already considered in the algorithm. The criterion states that if there is an elementp ofG such that the leading monomial ofp (LM(p )) divides the LCM(f, g ∈ G), and ifS(f, p) , S(p, g) have already been considered, then there is no need of consideringS(f, g) as this reduces to zero.

The essential observation in using theGroebner bases to solve a system of polynomial equations is that the variety (simultaneous solution of system of polynomial equations) does not depend on the original system of the polynomialsF : = f 1, , f s but instead on theIdealI generated byF. This therefore means that for the varietyV = V (I), one makes use of the special generating set (Groebner basis) instead of the actual systemF. Since theIdeal is generated byF, the solution obtained by solving for the affine variety of thisIdeal satisfies the original systemF of equations.B. Buchberger (1970) proved thatV (I) is void, and thus giving a test as to whether the system of polynomialF can be solved, if and only if the computedGroebner basis of polynomial IdealI has 1 as its element.B. Buchberger (1970) further gives the criterion for deciding ifV (I) is finite. If the system has been proved to be solvable and finite thenF. Winkler (1996,theorem 8.4.4, p.192) gives a theorem for deciding whether the system has finitely or infinitely many solutions. Thetheorem states that ifG is aGroebner basis, then a solvable system of polynomial equations has finitely many solutions if and only if for everyx i ,  1 ≤ i ≤ n, there is a polynomialg i  ∈ G such thatLM(g i ) is a pure power ofx i . The process of addition of the remainder after the reduction by theS-polynomials and thus expanding the generating set is shown byB. Buchberger (1970),D. Cox et al. (1997 p.88) andJ. H. Davenport et al. (1988 p.101) to terminate.

TheB. Buchberger algorithm (see Appendix D2), more or less a generalization of theGauss elimination procedure, makes use of the subtraction polynomials known as theS-polynomials inDefinition 15.7 to eliminate the leading terms of a pair of polynomials. In so doing and iflexicographic ordering is chosen, the process end up with one of the computedS-polynomials being a univariate polynomial which can be solved and substituted back in the otherS-polynomials using theExtension Theorem (D. Cox et al. 1998, pp. 25–26) to obtain the other variables. TheGroebner bases approach adds to the treasures of methods that are used to solvenonlinear algebraic systems of equations in Geodesy, Photogrammetry, Machine Vision, Robotics and Surveying.

Having defined the termGroebner basis and illustrated how theB. Buchberger algorithm computes theGroebner bases, we briefly mention here how theGroebner bases can be computed using algebraic softwares of Mathematica and Maple. In Mathematica Versions 8, theGroebner basis command is executed by writingIn[1]:=GroebnerBasis[{polynomoials},{variables},{elims}] (where In[1]:= is the mathematica prompt) which computes theGroebner basis for theideal generated by thepolynomials with respect to themonomial order specified bymonomial order options with thevariables specified as in the executable command giving the reducedGroebner basis. Without specifying the elims part that indicates the elimination order, one gets too many elements of theGroebner basis which may not be relevant. In Maple Version 5 the command is accessed by typing> with (grobner); (where > is the Maple prompt and the semicolon ends the Maple command). Once theGroebner basis package has been loaded, the execution command then becomes> g basis (polynomials, variables, termorder) which computes theGroebner basis for theideal generated by thepolynomials with respect to themonomial ordering specified bytermorder andvariables in the executable command. Following suggestions fromB. Buchberger (1999), Mathematica software is applied to the examples presented in Sect. 15.4.

15.3.1.3 Multipolynomial Resultants Method

Whereas the resultant of two polynomials is well known and algorithms for computing it are well incorporated into computer algebra packages such asMaple, theMultipolynomial resultant, i.e. the resultant of more than two polynomials still remain an active area of research. In Geodesy, the use of the two polynomial resultant also known as theSylvester resultant is exemplified in the work ofP. Lohse (1994, pp. 72–76). This section therefore extends on the use ofSylvester resultants to resultants of more than two polynomials of multiple variables(Multipolynomial resultant).J.L. Awange andE.W. Grafarend (2005) andJ. L. Awange et al. (2010) illustrate how the tool can be exploited in Geodesy and Geoinformatics to solvenonlinear system of equations.

The necessity ofMultipolynomial resultant method in Geodesy is due to the fact that many geodetic problems involve the solution of more than two polynomials of multiple variables. This is true since we are living in a three-dimensional world. We shall therefore understand the termmultipolynomial resultants to mean resultants of more than two polynomials. We treat it as a tool besides theGroebner bases and perhaps more powerful to eliminate variables in solution of polynomial systems. Publications on the subject can be found in the works ofG. Salmon (1876),F. Macaulay (1902, 1916),A. L. Dixon (1908), C. Bajaj et al. (1988),J. F. Canny (1988),J. F. Canny et al. (1989),I. M. Gelfand et al. (1990, 1994), D. Manocha (1992, 1993, 1994a,b,c),D. Manocha andJ.F. Canny (1991, 1992, 1993),G. Lyubeznik (1995), S. Krishna and D. Manocha (1995),J. Guckenheimer et al.(1997), B. Sturmfels (1994, 1998) andE. Cattani et al. (1998). Text books on the subject have been written byI. Gelfand et al. (1994), byD. Cox et al. (1998, pp.71–122) and, more recently, byJ.L. Awange andE.W. Grafarend (2005), andJ. L. Awange et al. (2010) who provides interesting material.

In order to understand theMultipolynomial resultants technique, we first present the simple case; the resultant of two polynomial also known as theSylvester resultant.

15.3.1.3.1 Resultant of two polynomials

Definition 15.8 (Homogeneous polynomial). 

If monomials of a polynomialp with non zero coefficients have the sametotal degree, the polynomialp is said to behomogeneous.

A homogeneous polynomial equation of total degree 2 is \(p = {x}^{2} + {y}^{2} + {z}^{2} + xy + xz + yz\) since the monomials {x, y, z, xy, xz, yz} all have the sum of their powers (total degree) being 2.

To set the ball rolling, we examine next the resultant of two univariate polynomialsp, q ∈ k[x] ofpositive degree as

$$\left.\begin{array}{c} p = {k}_{0}{x}^{i} + \ldots. + {k}_{i},\,{k}_{0}\neq 0,\,i > 0 \\ q = {l}_{0}{x}^{j} + \ldots. + {l}_{j},\,{l}_{0}\neq 0,\,j > 0 \end{array} \right \}$$
(15.40)

the resultant ofp andq, denoted Res(p, q), is the \((i + j) \times(i + j)\) determinant

$$\mathit{Res}\,(p,\,q) = \mathit{det}\left [\begin{array}{c@{\;\;}c@{\;\;}c@{\;\;}c@{\;\;}c@{\;\;}c@{\;\;}c@{\;\;}c@{\;\;}c@{\;\;}c@{\;\;}c@{\;\;}c} {k}_{0}\;\;&{k}_{1}\;\;&{k}_{2}\;\;&. \;\;&. \;\;&. \;\;&{k}_{i}\;\;& 0 \;\;& 0 \;\;& 0 \;\;& 0 \;\;& 0 \\ 0 \;\;&{k}_{0}\;\;&{k}_{1}\;\;&{k}_{2}\;\;&. \;\;&. \;\;&. \;\;&{k}_{i}\;\;& 0 \;\;& 0 \;\;& 0 \;\;& 0 \\ 0 \;\;& 0 \;\;&{k}_{0}\;\;&{k}_{1}\;\;&{k}_{2}\;\;&. \;\;&. \;\;&. \;\;&{k}_{i}\;\;& 0 \;\;& 0 \;\;& 0 \\ 0 \;\;& 0 \;\;& 0 \;\;&{k}_{0}\;\;&{k}_{1}\;\;&{k}_{2}\;\;&. \;\;&. \;\;&. \;\;&{k}_{i}\;\;& 0 \;\;& 0 \\ 0 \;\;& 0 \;\;& 0 \;\;& 0 \;\;&{k}_{0}\;\;&{k}_{1}\;\;&{k}_{2}\;\;&. \;\;&. \;\;&. \;\;&{k}_{i}\;\;& 0 \\ 0 \;\;& 0 \;\;& 0 \;\;& 0 \;\;& 0 \;\;&{k}_{0}\;\;&{k}_{1}\;\;&{k}_{2}\;\;&. \;\;&. \;\;&. \;\;&{k}_{i} \\ {l}_{0}\;\;&{l}_{1}\;\;&{l}_{2}\;\;&. \;\;&. \;\;&. \;\;&{l}_{j}\;\;& 0 \;\;& 0 \;\;& 0 \;\;& 0 \;\;& 0 \\ 0 \;\;&{l}_{0}\;\;&{l}_{1}\;\;&{l}_{2}\;\;&. \;\;&. \;\;&. \;\;&{l}_{j}\;\;& 0 \;\;& 0 \;\;& 0 \;\;& 0 \\ 0 \;\;& 0 \;\;&{l}_{0}\;\;&{l}_{1}\;\;&{l}_{2}\;\;&. \;\;&. \;\;&. \;\;&{l}_{j}\;\;& 0 \;\;& 0 \;\;& 0 \\ 0 \;\;& 0 \;\;& 0 \;\;&{l}_{0}\;\;&{l}_{1}\;\;&{l}_{2}\;\;&. \;\;&. \;\;&. \;\;&{l}_{j}\;\;& 0 \;\;& 0 \\ 0 \;\;& 0 \;\;& 0 \;\;& 0 \;\;&{l}_{0}\;\;&{l}_{1}\;\;&{l}_{2}\;\;&. \;\;&. \;\;&. \;\;&{l}_{j}\;\;& 0 \\ 0 \;\;& 0 \;\;& 0 \;\;& 0 \;\;& 0 \;\;&{l}_{0}\;\;&{l}_{1}\;\;&{l}_{2}\;\;&. \;\;&. \;\;&. \;\;&{l}_{j} \end{array} \right ]$$
(15.41)

where the coefficients of the first polynomialp of (15.40) occupiesj rows while those of the second polynomialq occupiesi rows. The empty spaces are occupied by zeros as shown above such that a square matrix is obtained. This resultant is also known as theSylvester resultant and has the following properties (B. Sturmfels 1998,D. Cox et al. 1998, §3.5)

  1. 1.

    Res(p, q) is a polynomial in \({k}_{0},\ldots,{k}_{i},\,{l}_{0},\ldots,{l}_{j}\) with integer coefficients

  2. 2.

    Res(p, q) = 0 if and only ifp(x) andq(x) have a common factor inQ[x]

  3. 3.

    There exist a polynomialr, s ∈ Q[x] such that \(rp + sq = \mathit{Res}(p,\,q)\)

Sylvester resultants can be used to solve two polynomial equations as shown in the example below

Consider the two equations

$$\left.\begin{array}{c} p := xy - 1 = 0 \\ q := {x}^{2} + {y}^{2} - 4 = 0. \end{array} \right \}$$
(15.42)

In order to eliminate one variable, we use thehide variable technique i.e. we consider one variable sayy as a constant (of degree zero). We then have theSylvester resultant from (15.41) as

$$\mathit{Res}\,(p,\,q,x) = \mathit{det}\left [\begin{array}{ccc} y& - 1& 0 \\ 0& y & - 1 \\ 1& 0 &{y}^{2} - 4 \end{array} \right ] = {y}^{4}-4{y}^{2}+1$$
(15.43)

which can be readily solved for the variabley and substituted back in any of the equations in (15.42) to get the values of the other variablex. For two polynomials, the construction of resultant is relatively simpler and algorithms for the execution are incorporated in computer algebra algorithms. For the resultant of more than 2 polynomials of multiple variables, we turn to theMultipolynomial resultants.

15.3.1.3.1.1 Multipolynomial resultants

In defining the termmultipolynomial resultant,D. Manocha (1994c) writes:

“Elimination theory, a branch of classical algebraic geometry, deals with conditions for common solutions of a system of polynomial equations. Its main result is the construction of a single resultant polynomial ofn homogeneous polynomial equations inn unknowns, such that the vanishing of the resultant is anecessary andsufficient condition for the given system to have a non-trivial solution. We refer to this resultant as themultipolynomial resultant and use it in the algorithm presented in the paper”.

We present here two approaches for the formation of the design matrix whose determinant we need; first the approach based onF. Macaulay (1902) formulation and then a more modern approach based onB. Sturmfels (1998) formulation.

Approach based on F. Macaulay (1902) formulation: Withn polynomials, the construction of the matrix whose entries are the coefficients of the polynomials \({f}_{1},\ldots,{f}_{n}\) can be done in five steps as illustrated in the following approach ofF. Macaulay (1902):

Step 1::

The given polynomials \({f}_{1} = 0,\ldots,{f}_{n} = 0\) are considered to be homogeneous equations in the variables \({x}_{1},\ldots,{x}_{n}\) and if not, they are homogenized. Let the degree of the polynomialf i bed i . The first step involves the determination of thecritical degree given byC. Bajaj et al. (1988) as

$$d = 1 + \sum\nolimits ({d}_{i} - 1).$$
(15.44)
Step 2::

Once thecritical degree has been established, the given monomials of the polynomial equations are multiplied with each other to generate a setX whose elements consists of monomials whose total degree equals thecritical degree. Thus if we are given the polynomial equations \({f}_{1} = 0,\ldots,{f}_{n} = 0\), then each monomial off 1is multiplied by those of \({f}_{2},\ldots,{f}_{n}\), those off 2 are multiplied by those of \({f}_{3},\ldots,{f}_{n}\) until those off n − 1 are multiplied by those off n . The setX of monomials generated in this form is

$${X}^{d} =\{ {x}^{d}\mid d = {\alpha }_{ 1} + {\alpha }_{2} + \cdots+ {\alpha }_{n}\}$$
(15.45)

and the variable \({x}^{d} = {x}_{1}^{{\alpha }_{1}}\ldots {x}_{n}^{{\alpha }_{n}}\).

Step 3::

The setX containing the monomials each of total degreed is now partitioned according to the following criteria (J. F. Canny 1988, p. 54)

$$\left \{\begin{array}{lrc} {X}_{1}^{d} & =& \{{x}^{\alpha } \in {X}^{d}\mid {\alpha }_{1} \geq{d}_{1}\} \\ {X}_{2}^{d} & =& \{{x}^{\alpha } \in {X}^{d}\mid {\alpha }_{2} \geq{d}_{2}\ \mathit{and}\ {\alpha }_{1} < {d}_{1}\}\\.&.&.\\.&.&.\\.&.&. \\ {X}_{n}^{d}& =&\{{x}^{\alpha } \in {X}^{d}\mid {\alpha }_{n} \geq{d}_{n}\ \mathit{and}\ {\alpha }_{i} < {d}_{i},\,\mathit{for}\,i = 1,\ldots,n - 1\}.\end{array} \right.$$
(15.46)

The resulting sets ofX i d are disjoint and every element ofX d is contained in exactly one of them.

Step 4::

From the resulting subsetsX i d ⊂ X d, a set of polynomialsF i which are homogeneous inn variables are defined as follows

$${F}_{i} = \frac{{X}_{i}^{d}} {{x}_{i}^{{d}_{i}}}{f}_{i}$$
(15.47)

from which asquare matrix \(A\) is now formed with the row elements being thecoefficients of the monomials of the polynomialsF i i = 1, , n and the columns correspond to theN monomials of the setX d. The formed square matrix is of the order \(\left (\begin{array}{c} d + n - 1\\ d\end{array} \right )\times \left (\begin{array}{c} d + n - 1\\ d\end{array} \right )\)and is such that for a given polynomialF i in (15.47), the row is made up of the symbolic coefficients of each polynomial. The square matrix \(A\) has a special property that the non trivial solution of the homogeneous equationsF i which also form the solution of the original equationsf i are in its null space. This implies that the matrix must be singular or its determinant, \(det(A)\), must be zero. For the determinant to vanish therefore, the original equationsf i and their homogenized counterpartsF i must have the same non trivial solutions.

Step 5::

After computing the determinant of the square matrix \(A\) above,F. Macaulay (1902) suggests the computation ofextraneous factor in order to obtain the resultant.D. Cox et al. (1998,Proposition 4.6, p.99) explains theextraneous factors to be integer polynomials in the coefficients of \(\bar{{F}}_{0},\ldots,\bar{{F}}_{n-1}\), where \(\bar{{F}}_{i} = {F}_{i}({x}_{0},\ldots,{x}_{n-1},0)\) and is related to the determinant via

$$\mathit{determinant} = \mathit{Res}({F}_{1},\ldots,{F}_{n}).\mathit{Ext}$$
(15.48)

with the determinant computed as in step 4,Res(F 1, , F n ) being theMultipolynomial resultant andExt the extraneous factor. This expression was established as early as 1902 byF. Macaulay (1902) and this procedure ofresultant formulation named after him.F. Macaulay (1902) determines the extraneous factor from the sub-matrix of theN ×N square matrix \(A\) and calls it a factor of minor obtained by deleting rows and columns of theN ×N matrix \(A\). A monomialx α of total degreed is said to be reduced if \({x}_{i}^{{d}_{i}}\) dividesx α for exactly onei. Theextraneous factor is obtained by computing the determinant of the sub-matrix of the coefficient matrix \(A\) obtained by deleting rows and columns corresponding to reduced monomialsx α. 

From the relationship of step 5, it suffices for our purpose to solve for the unknown variable hidden in the coefficients of the polynomialsf i by obtaining the determinant of theN ×N square matrix \(A\) and equating it to zero neglecting the extraneous factor. This is because the extraneous factor is an integer polynomial and as such not related to the variable in the determinant of \(A\). The existence of the non-trivial solution provides thenecessary andsufficient conditions for the vanishing of the determinant. It should be pointed out that there exists several procedures for computing the resultants as exemplified in the works ofD. Manocha (1992, 1993, 1994a,b,c) who solves theMultipolynomial resultants using the eigenvalue-eigenvector approach,J. F. Canny (1988) who solves the resultant using the characteristic polynomial approach andB. Sturmfels (1994, 1998) who proposes a more compact approach for solving the resultants of a ternary quadric using the Jacobian determinant approach which we present below.

Approach based on B. Sturmfels (1998, p.26) formulation: Given three homogeneous equations of degree two as follows

$$\begin{array}{ccc} {F}_{1} := {a}_{11}{x}^{2} + {a}_{12}{y}^{2} + {a}_{13}{z}^{2} + {a}_{14}xy + {a}_{15}xz + {a}_{16}yz& =&0 \\ {F}_{2} := {a}_{21}{x}^{2} + {a}_{22}{y}^{2} + {a}_{23}{z}^{2} + {a}_{24}xy + {a}_{25}xz + {a}_{26}yz& =&0 \\ {F}_{3} := {a}_{31}{x}^{2} + {a}_{32}{y}^{2} + {a}_{33}{z}^{2} + {a}_{34}xy + {a}_{35}xz + {a}_{36}yz& =&0\end{array}$$
(15.49)

we compute the Jacobian determinants of (15.49) by

$$J = det\left [\begin{array}{c@{\;\;}c@{\;\;}c} \frac{\partial {F}_{1}} {\partial x} \;\;&\frac{\partial {F}_{1}} {\partial y} \;\;&\frac{\partial {F}_{1}} {\partial z} \\ \frac{\partial {F}_{2}} {\partial x} \;\;&\frac{\partial {F}_{2}} {\partial y} \;\;&\frac{\partial {F}_{2}} {\partial z} \\ \frac{\partial {F}_{3}} {\partial x} \;\;&\frac{\partial {F}_{3}} {\partial y} \;\;&\frac{\partial {F}_{3}} {\partial z} \end{array} \right ]$$
(15.50)

which is a cubic polynomial in the coefficients {x, y, z}. Since the determinant polynomialJ in (15.50) is a cubic polynomial, its partial derivatives will be quadratic polynomials in variables {x, y, z} and can be written in the form

$$\begin{array}{rcl} \frac{\partial J} {\partial x} := {b}_{11}{x}^{2} + {b}_{ 12}{y}^{2} + {b}_{ 13}{z}^{2} + {b}_{ 14}xy + {b}_{15}xz + {b}_{16}yz& =& 0 \\ \frac{\partial J} {\partial y} := {b}_{21}{x}^{2} + {b}_{ 22}{y}^{2} + {b}_{ 23}{z}^{2} + {b}_{ 24}xy + {b}_{25}xz + {b}_{26}yz& =& 0 \\ \frac{\partial J} {\partial z} := {b}_{31}{x}^{2} + {b}_{ 32}{y}^{2} + {b}_{ 33}{z}^{2} + {b}_{ 34}xy + {b}_{35}xz + {b}_{36}yz& =& 0.\end{array}$$
(15.51)

The coefficientsb ij in (15.51) are quadratic polynomials ina ij of equation (15.49). The final step in computing the resultant of the initial system (15.49) involves the computation of the determinant of a 6 ×6 matrix given by

$${ \mathit{Res}}_{222}({F}_{1},{F}_{2},{F}_{3}) = \mathit{det}\left [\begin{array}{c@{\;\;}c@{\;\;}c@{\;\;}c@{\;\;}c@{\;\;}c} {a}_{11}\;\;&{a}_{12}\;\;&{a}_{13}\;\;&{a}_{14}\;\;&{a}_{15}\;\;&{a}_{16} \\ {a}_{21}\;\;&{a}_{22}\;\;&{a}_{23}\;\;&{a}_{24}\;\;&{a}_{25}\;\;&{a}_{26} \\ {a}_{31}\;\;&{a}_{32}\;\;&{a}_{33}\;\;&{a}_{34}\;\;&{a}_{35}\;\;&{a}_{36} \\ {b}_{11}\;\;&{b}_{12}\;\;&{b}_{13}\;\;&{b}_{14}\;\;&{b}_{15}\;\;&{b}_{12} \\ {b}_{21}\;\;&{b}_{22}\;\;&{b}_{23}\;\;&{b}_{24}\;\;&{b}_{25}\;\;&{b}_{26} \\ {b}_{31}\;\;&{b}_{32}\;\;&{b}_{33}\;\;&{b}_{34}\;\;&{b}_{35}\;\;&{a}_{36} \end{array} \right ]\!.$$
(15.52)

The resultant (15.52) vanishes if and only if (15.49) have a common solution {x, y, z}, where {x, y, z} are complex numbers or real numbers not all equal zero. In Sect. 15.4, we use the procedure to solve the GPS pseudo-range equations.

15.3.2 Adjustment of the combinatorial subsets

Once thecombinatorial minimal subsets have been solved using either theGroebner bases or theMultipolynomial resultant approach, the resulting sets of solution are considered as pseudo-observations. For each combinatorial, the obtained minimal subset solutions considered as pseudo-observations are used as the approximate values to generate the dispersion matrix via the nonlinear error propagation law/variance-covariance propagation (e.g. E. Grafarend andB. Schaffrin, 1993, pp. 469–471) as follows:

From the nonlinear geodetic observation equations that have been converted into its algebraic (polynomial) form, thecombinatorial minimal subsets will consist of polynomials \({f}_{1},\ldots,{f}_{m} \in k[{x}_{1},\ldots,{x}_{m}]\) with \(\{{x}_{1},\ldots,{x}_{m}\}\) being the unknown variables (fixed parameters) to be determined and \(\{{y}_{1},\ldots,{y}_{n}\}\) the known variables comprising the observations or pseudo-observations. We write the polynomials as

$$\left.\begin{array}{c} {f}_{1} := g({x}_{1},\ldots,{x}_{m},{y}_{1},\ldots,{y}_{n}) = 0 \\ {f}_{2} := g({x}_{1},\ldots,{x}_{m},{y}_{1},\ldots,{y}_{n}) = 0\\.\\.\\. \\ {f}_{m} := g({x}_{1},\ldots,{x}_{m},{y}_{1},\ldots,{y}_{n}) = 0 \end{array} \right \}$$
(15.53)

which is expressed in matrix form as

$$f := g(x,y) = 0$$
(15.54)

where the unknown variables \(\{{x}_{1},\ldots,{x}_{m}\}\) are placed in a vector \(x\) and the known variables \(\{{y}_{1},\ldots,{y}_{n}\}\) are placed in the vector \(y\). Next, we implement the error propagation from the observations (pseudo-observations) \(\{{y}_{1},\ldots,{y}_{n}\}\) to the parameters \(\{{x}_{1},\ldots,{x}_{m}\}\) that are to be explicitly determined namely characterized by thefirst moments, the expectation \(E\{x\} ={ \mu }_{x}\) and \(E\{y\} ={ \mu }_{y}\), as well as thesecond moments, the variance-covariance matrices/dispersion matrices \(D\{x\} ={ \Sigma }_{x}\) and \(D\{y\} ={ \Sigma }_{y}\). FromE. Grafarend andB. Schaffrin (1993, pp. 470–471), we have up to nonlinear terms

$$D\{x\} ={ J}_{x}^{-1}{J}_{ y}{\Sigma }_{y}{J}_{y}^{{\prime}}{({J}_{ x}^{-1})}^{{\prime}}$$
(15.55)

with \({J}_{x},{J}_{y}\) being the partial derivatives of (15.54) with respect to \(x\),\(y\) respectively at the Taylor points \(({\mu }_{x},{\mu }_{y})\). The approximate values of unknown parameters \(\{{x}_{1},\ldots,{x}_{m}\} \in x\) appearing in the Jacobi matrices \({J}_{x},{J}_{y}\) are obtained fromGroebner bases orMultipolynomial resultants solution of thenonlinear system of equations (15.53).

Given \({J}_{i} ={ J}_{{x}_{i}}^{-1}{J}_{{y}_{i}}\) from theith combination and \({J}_{j} ={ J}_{{x}_{j}}^{-1}{J}_{{y}_{j}}\) from thejth combinatorials, the correlation between theith andjth combination is given by

$${ \Sigma }_{ij} ={ J}_{j}{\Sigma }_{{y}_{j}{y}_{i}}{J}_{i}^{{\prime}}$$
(15.56)

The sub-matrices variance-covariance matrix for the individual combinatorials \({\Sigma }_{1},{\Sigma }_{2},{\Sigma }_{3},\ldots,{\Sigma }_{k}\) (wherek is the number of combinations ) obtained via (15.55) and the correlations between combinatorials obtained from (15.56) form the variance-covariance/dispersion matrix

$$\Sigma= \left [\begin{array}{c@{\;\;}c@{\;\;}c@{\;\;}c@{\;\;}c@{\;\;}c} {\Sigma }_{1} \;\;&{\Sigma }_{12}\;\;&. \;\;&.\;\;&.\;\;&{\Sigma }_{1k} \\ {\Sigma }_{21}\;\;& {\Sigma }_{2} \;\;&. \;\;&.\;\;&.\;\;&{\Sigma }_{2k} \\. \;\;& \;\;&{\Sigma }_{3}\;\;& \;\;& \;\;& \\. \;\; & \;\; & \;\; &.\;\; & \;\; &\\. \;\; & \;\; & \;\; & \;\; &.\;\; & \\ {\Sigma }_{k1}\;\;&. \;\;&. \;\;&.\;\;&\;\;& {\Sigma }_{k} \end{array} \right ]$$
(15.57)

for the entirek combinations. The obtained dispersion matrix \(\Sigma \) is then used in thelinear Gauss–Markov model to obtain the estimates \(\hat{\xi }\) of the unknown parameters \(\xi \) with the combinatorial solution points in a polyhedron considered as pseudo-observations in the vector \(y\) of observations while the design matrix \(A\) comprises of integer values 1 being the coefficients of the unknowns as in (15.60). In order to understand the adjustment process, we consider the following example.

Consider that one has four observations in three unknowns in a nonlinear problem. Let the observations be given by {y 1, y 2, y 3, y 4} leading to four combinations giving the solutions \({z}_{I}({y}_{1},{y}_{2},{y}_{3})\), \({z}_{II}({y}_{2},{y}_{3},{y}_{4})\), \({z}_{III}({y}_{1},{y}_{3},{y}_{4})\) and \({z}_{IV }({y}_{1},{y}_{2},{y}_{4})\). If the solutions are placed in a vector \({z}_{J} = {[\begin{array}{c@{\;\,}c@{\;\,}c@{\;\,}c} {z}_{I}\;\,&{z}_{II}\;\,&{z}_{III}\;\,&{z}_{IV } \end{array} ]}^{{\prime}},\) the adjustment model is then defined as

$$E\{{z}_{J}\} ={ I}_{12\times 3}{\xi }_{3\times 1},D\{{z}_{J}\}\,\mathit{from\ Variance/Covariance\ propagation.}$$
(15.58)

Let

$${ \xi }^{n} = L{z}_{ J}\ \mathit{subject\ to}\ {z}_{J} := \left [\begin{array}{c} {z}_{I} \\ {z}_{\mathit{II}} \\ {z}_{\mathit{III}} \\ {z}_{\mathit{IV}} \end{array} \right ] \in {\mathbb{R}}^{12\times 1}$$
(15.59)

such that the postulations \(\mathit{trD}\{{\xi }^{n}\} = \mathit{min}\) i.e. “best” and \(E\{{\xi }^{n}\} = \xi \,\mathit{for\,all}\,{\xi }^{n} \in {\mathbb{R}}^{m}\) i.e. “uniformlyunbiased” holds. We then have from (15.57), (15.58) and (15.59) the result

$$\begin{array}{rcl} \hat{\xi }& =& ({I}_{3\times 12}^{{\prime}}{\Sigma }_{{ z}_{J}}{I}_{12\times 3}){I}_{3\times 12}^{{\prime}}{\Sigma }_{{ z}_{J}}^{-1}{z}_{ J} \\ \hat{L}& =& \mathit{arg}\{trD\{{\xi }^{n}\} = \mathit{tr}\,L{\Sigma }_{ y}{L}^{{\prime}} = \mathit{min}\mid \mathit{UUE}\}\end{array}$$
(15.60)

The dispersion matrix \(D\{\hat{\xi }\}\) of the estimates \(\hat{\xi }\) is obtained. The shift from arithmetic weighted mean to the use oflinear Gauss Markov model is necessitated as we do not readily have the weights of the minimal combinatorial subsets but their dispersion which we obtain viaerror propagation/Variance-Covariance propagation. If we employ the equivalencetheorem ofE. Grafarend andB. Schaffrin (1993, pp. 339–341), an adjustment using linear Gauss markov model instead of weighted arithmetic mean is permissible.

Let us consider a simple case of the planar ranging problem. From an unknown point \(P(X,Y ) \in {\mathbb{E}}^{2}\), let distancesS 1 andS 2 be measured to two known points \({P}_{1}({X}_{1},{Y }_{1}) \in {\mathbb{E}}^{2}\) and \({P}_{2}({X}_{2},{Y }_{2}) \in {\mathbb{E}}^{2}\) respectively. We have the distance equations expressed as

$$\left [\begin{array}{c} {S}_{1}^{2} = {({X}_{1} - X)}^{2} + {({Y }_{1} - Y )}^{2} \\ {S}_{2}^{2} = {({X}_{2} - X)}^{2} + {({Y }_{2} - Y )}^{2} \end{array} \right.$$
(15.61)
$$\mathrm{whichweexpressinalgebraicform(15.53)as}$$
$$\left [\begin{array}{c} {f}_{1} := {({X}_{1} - X)}^{2} + {({Y }_{1} - Y )}^{2} - {S}_{1}^{2} = 0 \\ {f}_{2} := {({X}_{2} - X)}^{2} + {({Y }_{2} - Y )}^{2} - {S}_{2}^{2} = 0 \end{array} \right..$$
(15.62)
$$\mathrm{Ontakingtotaldifferentialof(15.62)wehave}$$
$$\left [\begin{array}{c} d{f}_{1} := 2({X}_{1} - X)d{X}_{1} - 2({X}_{1} - X)dX + 2({Y }_{1} - Y )d{Y }_{1}- \\ - 2({Y }_{1} - Y )dY - 2{S}_{1}d{S}_{1} = 0 \\ d{f}_{2} := 2({X}_{2} - X)d{X}_{2} - 2({X}_{2} - X)dX + 2({Y }_{2} - Y )d{Y }_{2}- \\ - 2({Y }_{2} - Y )dY - 2{S}_{2}d{S}_{2} = 0 \end{array} \right.$$
(15.63)

which on arranging the differential vector of the unknown terms \(\{X,Y \} =\{ {x}_{1},{x}_{2}\} \in x\) on the left hand side and that of the known terms \(\{{X}_{1},{Y }_{1},{X}_{2},{Y }_{2},{S}_{1},{S}_{2}\} =\{ {y}_{1},{y}_{2},{y}_{3},{y}_{4},{y}_{5},{y}_{6}\} \in y\) on the right hand side leads to

$${ J}_{x}\left [\begin{array}{c} dX\\ dY \end{array} \right ] ={ J}_{y}\left [\begin{array}{c} d{S}_{1} \\ d{X}_{1} \\ d{Y }_{1} \\ d{S}_{2} \\ d{X}_{2} \\ d{Y }_{2} \end{array} \right ]$$
(15.64)
$$\mbox{ with}$$
$${ J}_{x} = \left [\begin{array}{c@{\;\;\,}c} \frac{\partial {f}_{1}} {\partial X} \;\;\,&\frac{\partial {f}_{1}} {\partial Y }\\ \;\;\, & \\ \frac{\partial {f}_{2}} {\partial X} \;\;\,&\frac{\partial {f}_{2}} {\partial Y } \end{array} \right ] = \left [\begin{array}{cc} - 2({X}_{1} - X)& - 2({Y }_{1} - Y ) \\ - 2({X}_{2} - X)& - 2({Y }_{2} - Y ) \end{array} \right ]$$
(15.65)
$$\mbox{ and}$$
$$\left.\begin{array}{c} {J}_{y} = \left [\begin{array}{c@{\;\;\,}c@{\;\;\,}c@{\;\;\,}c@{\;\;\,}c@{\;\;\,}c} \frac{\partial {f}_{1}} {\partial {S}_{1}} \;\;\,& \frac{\partial {f}_{1}} {\partial {X}_{1}} \;\;\,& \frac{\partial {f}_{1}} {\partial {Y }_{1}} \;\;\,& 0 \;\;\,& 0 \;\;\,& 0\\ \;\;\, & \;\;\, & \;\;\, & \;\;\, & \;\;\, & \\ 0 \;\;\,& 0 \;\;\,& \;\;\,&\frac{\partial {f}_{2}} {\partial {S}_{2}} \;\;\,& \frac{\partial {f}_{2}} {\partial {X}_{2}} \;\;\,& \frac{\partial {f}_{2}} {\partial {Y }_{2}} \end{array} \right ] = \\ \\ = \left [\begin{array}{c@{\;\;}c@{\;\;}c@{\;\;}c@{\;\;}c@{\;\;}c} 2{S}_{1}\;\;& - 2({X}_{1} - X)\;\;& - 2({Y }_{1} - Y )\;\;& 0 \;\;& 0 \;\;& 0 \\ 0 \;\;& 0 \;\;& 0 \;\;&2{S}_{2}\;\;& - 2({X}_{2} - X)\;\;& - 2({Y }_{2} - Y ) \end{array} \right ] \end{array} \right ]$$
(15.66)
$$\mathrm{Ifweconsiderthat}$$
$$D\{x\} ={ \Sigma }_{x} = \left [\begin{array}{c@{\;\;}c} {\sigma }_{X}^{2}\;\;&{\sigma }_{XY } \\ {\sigma }_{Y X}\;\;&{\sigma }_{Y }^{2} \end{array} \right ]$$
(15.67)
$$\mbox{ and}$$
$$D\{y\} ={ \Sigma }_{y} = \left [\begin{array}{c@{\;\;}c@{\;\;}c@{\;\;}c@{\;\;}c@{\;\;}c} {\sigma }_{{S}_{1}}^{2} \;\;&{\sigma }_{{S}_{1}{X}_{1}}\;\;&{\sigma }_{{S}_{1}{Y }_{1}}\;\;&{\sigma }_{{S}_{1}{X}_{2}}\;\;& {\sigma }_{{S}_{1}{S}_{2}} \;\;&{\sigma }_{{S}_{1}{Y }_{2}} \\ {\sigma }_{{X}_{1}{S}_{1}}\;\;& {\sigma }_{{X}_{1}}^{2} \;\;&{\sigma }_{{X}_{1}{Y }_{1}}\;\;&{\sigma }_{{X}_{1}{S}_{2}}\;\;&{\sigma }_{{X}_{1}{X}_{2}}\;\;&{\sigma }_{{X}_{1}{Y }_{2}} \\ {\sigma }_{{Y }_{1}{S}_{1}}\;\;&{\sigma }_{{Y }_{1}{X}_{1}}\;\;& {\sigma }_{{Y }_{1}}^{2} \;\;&{\sigma }_{{Y }_{1}{S}_{2}}\;\;&{\sigma }_{{Y }_{1}{X}_{2}}\;\;&{\sigma }_{{Y }_{1}{Y }_{2}} \\ {\sigma }_{{S}_{2}{S}_{1}}\;\;&{\sigma }_{{S}_{2}{X}_{1}}\;\;&{\sigma }_{{S}_{2}{Y }_{1}}\;\;& {\sigma }_{{S}_{2}}^{2} \;\;&{\sigma }_{{S}_{2}{X}_{2}}\;\;&{\sigma }_{{S}_{2}{Y }_{2}} \\ {\sigma }_{{X}_{2}{S}_{1}}\;\;&{\sigma }_{{X}_{2}{X}_{1}}\;\;&{\sigma }_{{X}_{2}{Y }_{1}}\;\;&{\sigma }_{{X}_{2}{S}_{2}}\;\;& {\sigma }_{{X}_{2}}^{2} \;\;&{\sigma }_{{X}_{2}{Y }_{2}} \\ {\sigma }_{{Y }_{2}{S}_{1}}\;\;&{\sigma }_{{Y }_{2}{X}_{1}}\;\;&{\sigma }_{{Y }_{2}{Y }_{1}}\;\;&{\sigma }_{{Y }_{2}{S}_{2}}\;\;&{\sigma }_{{Y }_{2}{X}_{2}}\;\;& {\sigma }_{{Y }_{2}}^{2} \end{array} \right ]$$
(15.68)

we obtain with (15.64), (15.65) and (15.66) the dispersion (15.55) of the unknown variables \(\{X,Y \} =\{ {x}_{1},{x}_{2}\} \in x\). The unknown values of \(\{X,Y \} =\{ {x}_{1},{x}_{2}\} \in x\) appearing in the Jacobi matrices (15.65) and (15.66) are obtained fromGroebner bases orMultipolynomial resultants solution of thenonlinear system of equations (15.61).

15.4 Examples

In this section, we present three examples fromJ.L. Awange andE.W. Grafarend (2005) andJ. L. Awange et al. (2010) based on geodetic problems ofMinimum Distance Mapping,GPS pseudo-ranging four-pointP4P, and the7-parameter transformation.

In order to relate a pointP on the Earth’s topographic surface to a point on the international reference ellipsoid \({\mathbb{E}}_{a,a,b}^{2}\), a bundle of half straight lines so calledprojection lines which depart from P and intersect \({\mathbb{E}}_{a,a,b}^{2}\) either not at all or in two points are used. There isone projection line which is at minimum distance relating P to p. Figure 15.1 is an illustration of such a minimum distance mapping. Let us formulate such an optimization problem by means of the Lagrangean \(\mathit{\pounds }({x}_{1},{x}_{2},{x}_{3},{x}_{4})\) in Solution 15.1.

In the first case, theEuclidean distance between points P and p in terms ofCartesian coordinates ofP(X, Y, Z) and ofp(x 1, x 2, x 3) is represented. TheCartesian coordinates (x 1, x 2, x 3) of the projection pointP are unknown. The constraint that the point p is an element of the ellipsoid of revolution

$${\mathbb{E}}_{a,a,b}^{2} :=\{ \mathbf{x} \in {\mathbb{R}}^{3}\vert{b}^{2}({x}_{ 1}^{2} + {x}_{ 2}^{2}) + {a}^{2}{x}_{ 2}^{3} -{a}^{2}{b}^{2} = 0, {\mathbb{R}}^{+} \ni a > b \in{\mathbb{R}}^{+}\}$$

is substituted into the Lagrangean by means of the Lagrange multiplierx 4, which is unknown too. \(\{({x}_{1}^{\wedge },{x}_{2}^{\wedge },{x}_{3}^{\wedge },{x}_{4}^{\wedge }) = arg\{\mathit{\pounds }({x}_{1},{x}_{2},{x}_{3},{x}_{4}) = min\}\) is the argument of the minimum of theconstrained Lagrangean \(\mathit{\pounds }({x}_{1},{x}_{2},{x}_{3},{x}_{4})\). The result of the minimization procedure is presented by Lemma 15.10. Equation (15.71) provides thenecessary conditions to constitute an extremum: The normal equations are of bilinear type.Products of the unknowns for instancex 1 x 4, x 2 x 4, x 3 x 4 andsquares of the unknowns, for instancex 1 2, x 2 2, x 3 2 appear. Finally the matrix of second derivatives \({\vec{H}}_{3}\) in (15.73) which ispositive definite constitutes thesufficient condition to obtain a minimum. Fortunately the matrix of second derivatives \({\vec{H}}_{3}\) isdiagonal. Using (15.72i–15.72iv), together with (15.75) leads to (15.76), which are the eigenvalues of the Hesse matrix \({\vec{H}}_{3}\). These values are \({\Lambda }_{1} = {\Lambda }_{2} = X\setminus {x}_{1}^{\wedge },{\Lambda }_{3} = Z\setminus {x}_{3}^{\wedge }\) and must be positive.

Lemma 15.10 (Constrained minimum distance mapping). 

The functional \(\mathit{\pounds }({x}_{1},{x}_{2},{x}_{3},{x}_{4})\) is minimal, if the conditions (15.71) and (15.73) hold.

On taking partial derivatives with respect tox i , we have

$$\left [\begin{array}{l} (i)\,\,\,\,\,\, \frac{\partial \mathit{\pounds }} {\partial ({x}_{1}^{\wedge })} = -(X - {x}_{1}^{\wedge }) + {b}^{2}{x}_{ 1}^{\wedge }{x}_{ 4}^{\wedge } = 0 \\ \\ (\mathit{ii})\,\,\,\,\,\, \frac{\partial \mathit{\pounds }} {\partial ({x}_{2}^{\wedge })} = -(Y - {x}_{2}^{\wedge }) + {b}^{2}{x}_{ 2}^{\wedge }{x}_{ 4}^{\wedge } = 0 \\ \\ (\mathit{iii})\,\,\,\,\,\, \frac{\partial \mathit{\pounds }} {\partial ({x}_{3}^{\wedge })} = -(Z - {x}_{3}^{\wedge }) + {a}^{2}{x}_{ 3}^{\wedge }{x}_{ 4}^{\wedge } = 0 \\ \\ (\mathit{iv})\,\,\,\,\,\, \frac{\partial \mathit{\pounds }} {\partial ({x}_{4}^{\wedge })} = \frac{1} {2}[{b}^{2}({x}_{ 1}^{\wedge 2} + {x}_{ 2}^{\wedge 2})] + {a}^{2}{x}_{ 3}^{\wedge 2} - {a}^{2}{b}^{2} = 0 \end{array} \right.$$
(15.72)
$$\begin{array}{rcl}{ \mathbf{H}}_{3}& :=& \left [ \frac{{\partial }^{2}\mathit{\pounds }} {\partial {x}_{i}\partial {x}_{j}}({\mathbf{x}}^{\wedge })\right ] \\ & =& \left [\begin{array}{ccc} 1 + {b}^{2}{x}_{4}^{\wedge }& 0 \\ 0 &1 + {b}^{2}{x}_{4}^{\wedge }& 0 \\ 0 & 0 &1 + {a}^{2}{x}_{4}^{\wedge } \end{array} \right ]\quad \in \quad {\mathbb{R}}^{3\times 3}\end{array}$$
(15.74)
$$\textquotedblleft \mathit{e}igenvalues \prime \prime$$
$$\vert {\mathbf{H}}_{3} - \Lambda {\mathbf{I}}_{3}\vert= 0\quad \Longleftrightarrow$$
(15.75)
$$\left [\begin{array}{c} {\Lambda }_{1} = {\Lambda }_{2} := 1 + {b}^{2}{x}_{4}^{\wedge } = \frac{X} {{x}_{1}^{\wedge }} = \frac{Y } {{x}_{2}^{\wedge }} \\ {\Lambda }_{3} := 1 + {a}^{2}{x}_{4}^{\wedge } = \frac{Z} {{x}_{3}^{\wedge }}\ \end{array} \right.$$
(15.76)

Without the various forward and backward reduction steps, we could automatically generate an equivalent algorithm for solving the normal equations (15.72i)–(15.72iv) in a closed form by means of Groebner basis approach.

Fig. 15.1
figure 1

Minimum distance mapping of a pointP on the Earth’s topographic surface to a pointp on theinternational reference ellipsoid \({\mathbb{E}}_{a,a,b}^{2}\)

Groebner basis computation of the normal equations (15.72i)–(15.72iv) leads to 14 elements presented in Solution 15.1 interpreted as follows: Thefirst expression is a univariate polynomial of order four (quartic) in the Lagrange multiplier, i.e.,

With the admissible valuesx 4 substituted in the linear equations (4),(8),(12) of the computed Groebner basis, i.e.,

$$\left [\begin{array}{l} (1 + {a}^{2}{x}_{4}){x}_{3} - Z \\ (1 + {b}^{2}{x}_{4}){x}_{2} - Y \\ (1 + {b}^{2}{x}_{4}){x}_{1} - X,\\ \end{array} \right.$$
(15.78)

the values (x 1, x 2, x 3) = (x, y, z) are finally produced.

Solution 15.1 (Groebner basis MDM solution; see Awange et al. 2010 for further details and examples). 

(1)

$$\left [\begin{array}{l} {a}^{2}{b}^{2}{x}_{4}^{4} + (2{a}^{6}{b}^{4} + 2{a}^{4}{b}^{6}){x}_{4}^{3} + ({a}^{6}{b}^{2} + 4{a}^{4}{b}^{4} + {a}^{2}{b}^{6} - {a}^{4}{b}^{2}{X}^{2} - {a}^{4}{b}^{2}{Y }^{2}- \\ {a}^{2}{b}^{4}{Z}^{2}){x}_{4}^{2} + +(2{a}^{4}{b}^{2} + 2{a}^{2}{b}^{4} - 2{a}^{2}{b}^{2}{X}^{2} - 2{a}^{2}{b}^{2}{Y }^{2} - 2{a}^{2}{b}^{2}{Z}^{2}){x}_{4} \\ + ({a}^{2}{b}^{2} - {b}^{2}{X}^{2} - {b}^{2}{Y }^{2} - {a}^{2}{Z}^{2}). \end{array} \right.$$

(2)

$$\left [\begin{array}{l} ({a}^{4}Z - 2{a}^{2}{b}^{2}Z + {b}^{4}Z){x}_{3} - {a}^{6}{b}^{6}{x}_{4}^{3} - (2{a}^{6}{b}^{4} + {a}^{4}{b}^{6}){x}_{4}^{2} \\ - ({a}^{6}{b}^{2} + 2{a}^{4}{b}^{4} - {a}^{4}{b}^{2}{X}^{2} - {a}^{4}{b}^{2}{Y }^{2} - {a}^{2}{b}^{4}{Z}^{2}){x}_{4} \\ - {a}^{2}{b}^{4} + {a}^{2}{b}^{2}{X}^{2} + {a}^{2}{b}^{2}{Y }^{2} + 2{a}^{2}{b}^{2}{Z}^{2} - {b}^{4}{Z}^{2}. \end{array} \right.$$

(3)

$$\left [\begin{array}{l} (2{b}^{2}Z + {b}^{4}{x}_{4}Z - {a}^{2}Z){x}_{3} + {a}^{4}{b}^{6}{x}_{4}^{3} + (2{a}^{4}{b}^{4} + {a}^{2}{b}^{6}){x}_{4}^{2} \\ + ({a}^{4}{b}^{2} + 2{a}^{2}{b}^{4} - {a}^{2}{b}^{2}{X}^{2} - {a}^{2}{b}^{2}{Y }^{2} - {b}^{4}{Z}^{2}){x}_{4} \\ + {a}^{2}{b}^{2} - {b}^{2}{X}^{2} - {b}^{2}{Y }^{2} - 2{b}^{2}{Z}^{2}. \end{array} \right.$$

(4)

$$(1 + {a}^{2}{x}_{ 4}){x}_{3} - Z$$

(5)

$$\left [\begin{array}{l} ({a}^{4} - 2{a}^{2}{b}^{2} + {b}^{4}){x}_{ 3}^{2} + (2{a}^{2}{b}^{2}Z - 2{b}^{4}Z){x}_{ 3} \\ - {a}^{4}{b}^{6}{x}_{4}^{2} - 2{a}^{4}{b}^{4}{x}_{4} - {a}^{4}{b}^{2} + {a}^{2}{b}^{2}{X}^{2} + {a}^{2}{b}^{2}{Y }^{2} + {b}^{4}{Z}^{2}). \end{array} \right.$$

(6)

$$\left [\begin{array}{l} (2{b}^{2} - {a}^{2} + {b}^{4}{x}_{4}){x}_{3}^{2} - {a}^{2}Z{x}_{3} + {a}^{4}{b}^{6}{x}_{4}^{3} + (2{a}^{4}{b}^{4} + 2{a}^{2}{b}^{6}){x}_{4}^{2} \\ + +({a}^{4}{b}^{2} + 4{a}^{2}{b}^{4} - {a}^{2}{b}^{2}{X}^{2} - {a}^{2}{b}^{2}{Y }^{2} - {b}^{4}{Z}^{2}){x}_{4} \\ + 2{a}^{2}{b}^{2} - 2{b}^{2}X - 2b{Y }^{2} - 2{b}^{2}{Z}^{2}. \end{array} \right.$$

(7)

$$\left [\begin{array}{l} ({X}^{2} + {Y }^{2}){x}_{ 2} + {a}^{2}{b}^{4}Y {x}_{ 4}^{2} + Y ({a}^{2}{b}^{2} - {b}^{2}{x}_{ 3}^{2} - {b}^{2}Z{x}_{ 3}){x}_{4} \\ + Y {x}_{3}^{2} - {Y }^{3} - Y Z{x}_{3} - Y {X}^{2}. \end{array} \right.$$

(8)

$$(1 + {b}^{2}{x}_{ 4}){x}_{2} - Y$$

(9)

$${a}^{2}{x}_{ 3} - {b}^{2}{x}_{ 3} + {b}^{2}Z){x}_{ 2} - {a}^{2}{x}_{ 3}Y$$

(10)

$$Y {x}_{1} - X{x}_{2}$$

(11) \(X{x}_{1} + {a}^{2}{b}^{4}{x}_{4}^{2} + ({a}^{2}{b}^{2} + {b}^{2}{x}_{3}^{2} -{b}^{2}Z{x}_{3}){x}_{4} + {x}_{3}^{2} -Z{x}_{3} + Y {x}_{2} -{X}^{2} -{Y }^{2}.\) (12)

$$(1 + {b}^{2}{x}_{ 4}){x}_{1} - X$$

(13)

$$({a}^{2}{x}_{ 3} - {b}^{2}{x}_{ 3} + {b}^{2}Z){x}_{ 1} - {a}^{2}X{x}_{ 3}$$

(14)\({x}_{1}^{2} +{a}^{2}{b}^{4}{x}_{4}^{2} +(2{a}^{2}{b}^{2} +{b}^{2}{x}_{3}^{2} -{b}^{2}Z{x}_{3}){x}_{4} +2{x}_{3}^{2} -2Z{x}_{3} +{x}_{2}^{2} -{X}^{2} -{Y }^{2}.\)

E. Grafarend andJ. Shan (1996) have defined theGPS pseudo-ranging four-point problem (“pseudo 4P”) as the problem of determining the four unknowns comprising thethree components of the receiver position and the stationary receiverrange bias from four observed pseudo-ranges to four satellite transmitter of given geocentric position. Geometrically, the four unknowns are obtained from the intersection of four spherical cones given by the pseudo-ranging equations. Several procedures have been put forward for obtaining closed form solution of the problem. Amongst the procedures include the vectorial approach evidenced in the works ofS. Bancroft (1985), P. Singer et al. (1993),H. Lichtenegger (1995) andA. Kleusberg (1994,1999).E. Grafarend andJ. Shan (1996) propose two approaches.

One approach is based on the inversion of a \(3 \times3\) coefficient matrix of a linear system formed by differencing of thenonlinear pseudo-ranging equations in geocentric coordinates, while the other approach uses the coefficient matrix from the linear system to solve the same equations in barycentric coordinates. In this example we present both the approaches ofGroebner bases andMultipolynomial resultants (B. Sturmfel 1998 approach) to solve the same problem. We demonstrate our algorithms by solving theGPS Pseudo-ranging four-point problem already solved byA. Kleusberg (1994) andE. Grafarend andJ. Shan (1996). Both theGroebner basis and theMultipolynomial resultant solve the same linear equations as those ofE. Grafarend andJ. Shan (1996) and lead to identical results (see alsoJ. L . Awange andE. Grafarend 2002). We start with the pseudo-ranging equations written algebraically as

$$\left [\begin{array}{c} {({x}_{1} - {a}_{0})}^{2} + {({x}_{2} - {b}_{0})}^{2} + {({x}_{3} - {c}_{0})}^{2} - {({x}_{4} - {d}_{0})}^{2} = 0 \\ {({x}_{1} - {a}_{1})}^{2} + {({x}_{2} - {b}_{1})}^{2} + {({x}_{3} - {c}_{1})}^{2} - {({x}_{4} - {d}_{1})}^{2} = 0 \\ {({x}_{1} - {a}_{2})}^{2} + {({x}_{2} - {b}_{2})}^{2} + {({x}_{3} - {c}_{2})}^{2} - {({x}_{4} - {d}_{2})}^{2} = 0 \\ {({x}_{1} - {a}_{3})}^{2} + {({x}_{2} - {b}_{3})}^{2} + {({x}_{3} - {c}_{3})}^{2} - {({x}_{4} - {d}_{3})}^{2} = 0 \\ \mathit{where}\,\,\,\,{x}_{1},{x}_{2},{x}_{3},{x}_{4} \in \\ ({a}_{0},{b}_{0},{c}_{0}) = ({x}^{0},{y}^{0},{z}^{0}) \sim{P}^{0} \\ ({a}_{1},{b}_{1},{c}_{1}) = ({x}^{1},{y}^{1},{z}^{1}) \sim{P}^{1} \\ ({a}_{2},{b}_{2},{c}_{2}) = ({x}^{2},{y}^{2},{z}^{2}) \sim{P}^{2} \\ ({a}_{3},{b}_{3},{c}_{3}) = ({x}^{3},{y}^{3},{z}^{3}) \sim{P}^{3}\\ \end{array} \right.$$
(15.79)

withP 0, P 1, P 2, P 3 being the position of the four GPS satellites, their ranges to the stationary receiver atP given byd 0, d 1, d 2, d 3. The parameters

$$(\left \{{a}_{0},{b}_{0},{c}_{0}\right \},\left \{{a}_{1},{b}_{1},{c}_{1}\right \},\left \{{a}_{2},{b}_{2},{c}_{2}\right \},\left \{{a}_{3},{b}_{3},{c}_{3}\right \},\left \{{d}_{0},{d}_{1},{d}_{2},{d}_{3}\right \})$$

are elements of the spherical cone that intersect atP to give the coordinatesx 1, x 2, x 3 of the receiver and the stationary receiver range biasx 4. The equations above is expanded and arranged in thelexicographic orderx 1 > x 2 > x 3 > x 4 and re-written with the linear terms on one side and the nonlinear terms on the other as (Awange and Grafarend 2005; Awange et al. 2010)

$$\left [\begin{array}{l} {x}_{1}^{2} + {x}_{2}^{2} + {x}_{3}^{2} - {x}_{4}^{2} = 2{a}_{0}{x}_{1} + 2{b}_{0}{x}_{2} + 2{c}_{0}{x}_{3} - 2{d}_{0}{x}_{4} + {d}_{0}^{2} - {a}_{0}^{2} - {b}_{0}^{2} - {c}_{0}^{2} \\ {x}_{1}^{2} + {x}_{2}^{2} + {x}_{3}^{2} - {x}_{4}^{2} = 2{a}_{1}{x}_{1} + 2{b}_{1}{x}_{2} + 2{c}_{1}{x}_{3} - 2{d}_{1}{x}_{4} + {d}_{1}^{2} - {a}_{1}^{2} - {b}_{1}^{2} - {c}_{1}^{2} \\ {x}_{1}^{2} + {x}_{2}^{2} + {x}_{3}^{2} - {x}_{4}^{2} = 2{a}_{2}{x}_{1} + 2{b}_{2}{x}_{2} + 2{c}_{2}{x}_{3} - 2{d}_{2}{x}_{4} + {d}_{2}^{2} - {a}_{2}^{2} - {b}_{2}^{2} - {c}_{2}^{2} \\ {x}_{1}^{2} + {x}_{2}^{2} + {x}_{3}^{2} - {x}_{4}^{2} = 2{a}_{3}{x}_{1} + 2{b}_{3}{x}_{2} + 2{c}_{3}{x}_{3} - 2{d}_{3}{x}_{4} + {d}_{3}^{2} - {a}_{3}^{2} - {b}_{3}^{2} - {c}_{3}^{2}.\\ \end{array}\right.$$
(15.80)

On subtracting (15.80 iv) from (15.80 i), (15.80 ii), and (15.80 iii), to three equations which are linear with four unknowns are derived (see, e.g., Awange et al. 2010, p. 177). Treating the unknown variable \({x}_{4}\) as a constant, bothGroebner bases and theMultipolynomial resultant techniques are applied to solve the linear system of equation for \({x}_{1} = g({x}_{4}),{x}_{2} = g({x}_{4}),{x}_{3} = g({x}_{4})\), where \(g({x}_{4})\) is a univariate polynomial (see Awange and Grafarend 2005 or Awange et al 2010, pp. 177–180 for further details).

Consider a case where coordinates have been given in two systems, A and B. For clarity purposes, let us assume the two coordinate systems to be e.g., photo image coordinates in system A and ground coordinates in system B. The ground coordinates {X i , Y i , Z i  | i, , n} of the objects are obtained from, say, GPS measurements. Given the photo coordinates \(\{{a}_{i} = {x}_{i},{b}_{i} = {y}_{i},{c}_{i} = -f\vert i,\ldots,n\}\) and their equivalent ground coordinates {X i , Y i , Z i  | i, , n}, the 7-parameter datum transformation problem concerns itself with determining;

  1. (1)

    The scale parameter \({x}_{1} \in \mathbb{R}\)

  2. (2)

    Three translation parameters \({\vec{x}}_{2} \in {\mathbb{R}}^{3}\)

  3. (3)

    The rotation matrix \({\vec{X}}_{3} \in {\mathbb{R}}^{3\times 3}\) comprising three rotation elements.

Once the listed unknowns have been determined, coordinates can subsequently be transformed from one system onto another. The nonlinear equations relating these unknowns and coordinates from both systems are given by

$$\left [\begin{array}{c} {a}_{i} \\ {b}_{i} \\ {c}_{i} \end{array} \right ] = {x}_{1}{\vec{X}}_{3}\left [\begin{array}{c} {X}_{i} \\ {Y }_{i} \\ {Z}_{i} \end{array} \right ]+{\vec{x}}_{2}\mid i = 1,2,3,\ldots,n,$$
(15.81)
$$\mathrm{subjectto}$$

In (15.81), {a i , b i , c i } and {X i , Y i , Z i } are the coordinates of the same points, e.g., in both photo and ground coordinate systems respectively. The determination of the unknowns \({x}_{1} \in \mathbb{R}\), \({\vec{x}}_{2} \in {\mathbb{R}}^{3}\), \({\vec{X}}_{3} \in {\mathbb{R}}^{3\times 3}\) require a minimum ofthree points in both systems whose coordinates are known. Owing to the nonlinearity of (15.81), the solutions have always been obtained using a least squares approach iteratively. With this approach, (15.81) is first linearized and some initial approximate starting values of the unknown parameters used. The procedure then iterates, each time improving on the solutions of the preceding iteration step. This is done until a convergence criteria is achieved.

Where the rotation angles are small e.g., in photogrammetry, the starting values of zeros are normally used. In other fields such as geodesy, the rotation angles are unfortunately not small enough to be initialized by zeros, thereby making the solution process somewhat difficult and cumbersome. Bad choices of initial starting values often lead to many iterations being required before the convergence criteria is achieved. In some cases, where the initial starting values are far from those of the unknown parameters, iteration processes may fail to converge. With these uncertainties in the initial starting values, the cumbersomeness of the linearization and iterations, procedures that would offer an exact solution to the 7-parameter datum transformation problem would be desirable. To answer this challenge, we propose algebraic approaches whose advantages over the approximate numerical methods have already been mentioned.

Apart from the computational difficulties associated with numerical procedures, the 7-parameter datum transformation problem poses another challenge to existing algorithms. This is, the incorporation of the variance-covariance (weight) matrices of the two systems involved. In practice, users have been forced to rely on iterative procedures and linearized least squares solutions which are incapable of incorporating the variance-covariance matrices of both systems in play. This challenge is addressed algebraically inJ.L. Awange and E.W. Grafarend (2005) orJ.L. Awange et al. (2010, Chap. 17).

15.5 Notes

The current known techniques for solvingnonlinear polynomial equations can be classified intosymbolic,numeric andgeometric methods (D. Manocha 1994c). Symbolic methods, discussed in this chapter for solving closed form geodetic problems, apply theGroebner bases and theMultipolynomial resultants techniques to eliminate several variables in a multivariate system of equations in such a manner that the end product often consist ofunivariate polynomials whose roots can be determined by existing programs, however, such as the roots command in MATLAB. The current available programs however are efficient only for sets of low degree polynomial systems consisting of upto three to four polynomials due to the fact that computing the roots of theunivariate polynomials can be ill conditioned for polynomials of degree greater than 14 or 15 (D. Manocha 1994c).

Elaborate literature onGroebner bases can be found in the works ofB. Buchberger (1965, 1970),J. H. Davenport et al. (1988, pp. 95–103),F. Winkler (1996),D. Cox et al. (1997, pp.47-99),H. M. Mller (1998),W. V. Vasconcelos (1998),T. Becker andV. Weispfenning (1993,1998),B. Sturmfels (1996),G. Pistone andH.P. Wynn (1996),D. A. Cox (1998) andD. Cox et al.(1998, pp. 1–50), while literature onMultipolynomial resultants procedure include the works ofG. Salmon (1876),F. Macaulay (1902, 1921),A. L. Dixon (1908), B. L. van Waerden (1950),C. Bajaj et al. (1988),J. F. Canny (1988),J. F. Canny et al. (1989),I. M. Gelfand et al. (1990), J. Weiss (1993),D. Manocha (1992, 1993, 1994a,b,c, 1998),D. Manocha andJ. F. Canny (1991, 1992, 1993),I. M. Gelfand et al. (1994), G. Lyubeznik (1995), S. Krishna and D. Manocha (1995),J. Guckenheimer et al.(1997), B. Sturmfels (1994, 1998),E. Cattani et al. (1998) andD. Cox et al. (1998, pp.71–122). Besides theGroebner bases andresultant techniques, there exists another approach for variable elimination developed by WU Wen Tsün (W. T. Wu 1984) using the ideas proposed byJ. F. Ritt (1950). This approach is based on Ritts characteristic set construction and was successfully applied to automated geometric theorem by Wu. This algorithm is referred byX. S. Gao andS. C. Chou (1990) as theRitt-Wu’s algorithm (D. Manocha andF. Canny 1993).C. L. Cheng andJ. W. Van Ness (1999) have presented polynomial measurement error models.

Numeric methods for solving polynomials can be grouped intoiterative andhomotopy methods. Forhomotopy we refer toA. P. Morgan (1992). Also in this category are geometric methods which have found application in curve and surface intersection whose convergence are however said to be slow (D. Manocha 1994c). In general, for low degree curve intersection, the algebraic methods have been found to be the fastest in practice. In Sect. 15.3.1.2 and 15.3.1.3 we present in a nut shell the theories ofGroebner bases andMultipolynomial resultants.

The problem ofnonlinear adjustment in Geodesy as in other fields continues to attract more attention from the modern day researchers as evidenced in the works ofR. Mautz (2001) andL. Guolin (2000) who presents a procedure that tests using the F-distribution whether anonlinear model can be linearized or not. The solution to the minimization problem of thenonlinear Gauss–Markov model unlike its linear counter part does not have a direct method for solving it and as such, always relies on the iterative procedures such as the Steepest-descent method, Newton’s method and the Gauss-Newton’s method discussed byP. Teunissen (1990). In particular,P. Teunissen (1990) recommends the Gauss-Newton’s method as it exploits the structure of the objective function (sum of squares) that is to be minimized.P. Teunissen andE. H. Knickmeyer (1988) considers in a statistical way how the nonlinearity of a function manifests itself during the various stages of adjustment.E. Grafarend andB. Schaffrin (1989, 1991) while extending the work ofT. Krarup (1982) onnonlinear adjustment with respect to geometric interpretation have presented thenecessary andsufficient conditions for least squares adjustment ofnonlinear Gauss–Markov model and provided the geometrical interpretation of these conditions.

Other geometrical approaches include the works ofG. Blaha andR. P. Besette (1989) andK. Borre andS. Lauritzen (1983) while non geometrically treatment of nonlinear problems have been presented byK. K. Kubik (1967),T. Saito (1973),H. J. Schek andP. Maier (1976),A. Pope (1982) andH. G. Bähr (1988). A comprehensive review to the iterative procedures for solving thenonlinear equations is presented in the work ofP. Lohse (1994).M. Gullikson andI. Sderkvist (1995) have developed algorithms for fitting surfaces which have been explicitly or implicitly defined to some measured points with negative weights being acceptable by the algorithm.

Our approach in the present chapter goes back to the work ofC. F. Gauss (Published posthumously e.g.Awange et al. 2010) andC. G. I. Jacobi (1841). Within the framework of arithmetic mean,C. F. Gauss andC. G. I. Jacobi (1841) suggest that givenn linear(ized) observation equations withm unknowns (n > m), σ combinations, each consisting ofm equations be solved for the unknown elements and the weighted arithmetic mean be applied to get the final solution. WhereasC.F. Gauss proposes weighting by using the products of the square of the measured distances from the unknown point to known points and the distances of the side of the error triangles,C. G. I. Jacobi (1841) proposed the use of the square of the determinants as weights. In tracing the method of least squares to the arithmetic mean,A. T. Hornoch (1950) shows that the weighted arithmetic mean proposed byC. G. I. Jacobi (1841) leads to the least squares solution only if the weights used are the actual weights and the pseudo-observations formed by the combinatorial pairs are uncorrelated. Using a numerical example,S. Wellisch (1910, pp. 41–49) has shown the results of least squares solution to be identical to those of theGauss–Jacobi combinatorialalgorithm once proper weighting is applied.

P. Werkmeister (1920) illustrated that for planar resection case with three directional observations from the unknown point to three known stations, the area of the triangle (error figure) formed by the resulting three combinatorial coordinates of the new point is proportional to the determinant of the dispersion matrix of the coordinates of the new station. In this chapter, these concepts of the combinatoriallinear adjustment were extended tononlinear adjustment. Awange and Grafarend (2005) illustrated using a leveling network and planar ranging problem that the results ofGauss–Jacobi combinatorialalgorithm are identical to those oflinear Gauss–Markov model if the actual variance-covariance matrices are used.

To test the algebraic computational tools ofGroebner bases andMultipolynomial resultants, geodetic problems ofMinimum Distance Mapping andGNSS pseudo-ranging four-point P4P, and7-parameter transformation were solved. The procedures are, however, not limited only to these three examples but have been shown, e.g., inJ.L . Awange andE.W. Grafarend (2005) andJ. L . Awange et al. (2010) to be applicable to many problems in geodesy and geoinformatics. An example is the classical problem of resection. In general, the search towards the solution of the three-dimensional resection problem traces its origin to the work of a German mathematicianJ. A. Grunert (1841) whose publication appeared in the year 1841.J. A. Grunert (1841) solved the three-dimensional resection problem – what was then known as the“Pothenot’s” problem – in a closed form by solving an algebraic equation of degree four. The problem had hitherto been solved by iterative means mainly in Photogrammetry and Computer Vision. Procedures developed later for solving the three-dimensional resection problem revolved around the improvements of the approach ofJ. A. Grunert (1841) with the aim of searching for the optimal means of distances determination. WhereasJ. A. Grunert (1841) solves the problem by substitution approach inthree steps, the more recent desire has been to solve the distance equations in lesser steps as exemplified in the works ofS. Finsterwalder andW. Scheufele (1937),E. L . Merritt (1949),M. A. Fischler andR. C. Bolles (1981),S. Linnainmaa et al. (1988) andE. Grafarend, P. Lohse andB. Schaffrin (1989). Other research done on the subject of resection include the works ofF. J. Mller (1925),E. Grafarend andJ. Kunz (1965), R. Horaud et al. (1989),P. Lohse (1990), andE. Grafarend andJ. Shan (1997a, 1997b). An extensive review of some of the procedures above are presented byF. J. Müller (1925) andR. M. Haralick et al. (1991, 1994).

R. M. Haralick et al. (1994) reviewed the performance of six direct solutions (J. A. Grunert 1841,S. Finsterwalder andW. Scheufele 1937,E. L. Merritt 1949,M. A. Fischler andR. C. Bolles 1981,Linnainmaa et al. 1988, andE. Grafarend, P. Lohse andB. Schaffrin 1989) with the aim of evaluating the numerical stability of the solutions and the effect of permutation within the various solutions. All the six approaches follow the outline adopted byJ. A. Grunert (1841) with the major differences being the change in variables, the reduction of these variables, and the combination of different pairs of equations. The study revealed that the higher the order of the polynomials, the more complex the computations became and thus the less accurate the solutions were due numerical instabilities. Consequently, S. Finsterwalder’s (SF) procedure which solves a third order polynomial is ranked first, J. A Grunert (JG), Fischler and Bolles (FB), and Grafarend et al. (GLS) solutions are ranked second, Linnainmaa et al. solution which generates an eighth order polynomial is ranked third. Though it does not solve the eight order polynomial, the complexity of the polynomial is still found to be higher than those of the other procedures. An amazing result is that of Merritt’s procedure which is ranked last despite the fact that it is a fourth order polynomial and is similar to Grunert’s approach except for the pairs of equations used.R. M. Haralick et al. (1994) attributes the poor performance of Merritt’s procedure to the conversion procedure adopted byE. L Merritt (1949) in reducing the equations from forth to third order. For planar resection problem, solutions have been proposed e.g. byD. Werner (1913),G. Brandsttter (1974) andJ. van Mierlo (1988).

Overdetermined planar resection problem has been treated graphically byE. Hammer (1896),C. Runge (1900),P. Werkmeister (1916) andP. Werkmeister (1920).E. Gotthardt (1940,1974) dealt with the overdetermined two-dimensional resection where more than four points were considered with the aim of studying the critical configuration that would yield a solution. The work was later to be extended byK. Killian (1990). A special case of an overdetermined two-dimensional resection has also been considered byH. G. Bähr (1991) who uses six known stations and proposes the measuring of three horizontal angles which are related to the two unknown coordinates by nonlinear equations. By adopting approximate coordinates of the unknown point, an iterative adjustment procedure is performed to get the improved two-dimensional coordinates of the unknown point. It should be noted that the procedure is based on the coordinate system of the six known stations.K. Rinner (1962) has also contributed to the problem of overdetermined two-dimensional resection.

In order to relate a point on the Earth’s topographical surface uniquely (one-to-one) to a point on theInternational Reference Ellipsoid,E. Grafarend andP. Lohse (1991) have proposed the use of theMinimum Distance Mapping (i.e., Example 15.16). Other procedures that have been proposed are either iterative, approximate “closed”, closed or higher order equations. Iterative procedures include the works ofN. Bartelme andP. Meissl. (1975),W. Benning (1987),K. M. Borkowski (1987, 1989),N. Croceto (1993),A. Fitzgibbon et al. (1999),T. Fukushima (1999), W. Gander et al. (1994),B. Heck (1987),W. Heiskannen andH. Moritz (1976),R. Hirvonen andH. Moritz (1963),P. Loskowski (1991),K. C. Lin andJ. Wang (1995),M. K. Paul (1973), L. E. Sjoeberg (1999),T. Soler andL. D. Hothem (1989),W. Torge (1976),T. Vincenty (1978) andR. J. You (2000).

Approximate “closed” procedures includeB. R. Bowring (1976, 1985),A. Fotiou (1998),B. Hofman-Wellenhof et al. (1992),M. Pick (1985) andT. Vincenty (1980), while closed procedures includeW. Benning (1987),H. Frhlich andH. H. Hansen (1976),E. Grafarend et al. (1995) andM. Heikkinen (1982). Procedures that required the solution of higher order equations includeM. Lapaine (1990),M. J. Ozone (1985),P. Penev (1978), H. Snkel (1976), andP. Vaniceck and E. Krakiwski (1982). In this chapter, theMinimum Distance Mapping problem is solved using theGroebner bases approach. Aunivariate polynomial of fourth order is obtained together with 13 other elements of theGroebner basis. The obtainedunivariate polynomial and the linear terms are compared to those ofE. Grafarend andP. Lohse (1991). Other reference includeE. Grafarend andW. Keller (1995) andMukherjee, K. (1996)

TheGNSS four-point pseudo-ranging problem presented in Example 15.17 is concerned with the determination of the coordinates of a stationary receiver together with its range bias. Several closed form procedures have been put forward for obtaining closed form solution of the problem. Amongst the procedures include the vectorial approach as evidenced in the works ofL. O. Krause (1987),J. S. Abel andJ. W. Chaffee (1991),P. Singer et al. (1993),J. W. Chaffee andJ. S. Abel (1994), H. Lichtenegger (1995) andA. Kleusberg (1994, 1999).E. Grafarend andJ. Shan (1996) propose two approaches; one approach is based on a closed form solution of the nonlinear pseudo-ranging equations in geocentric coordinates while the other approach solves the same equations in barycentric coordinates.

S. C. Han et al. (2001) have developed an algorithm for very accurate absolute positioning through Global Positioning System (GPS) satellite clock estimation whileS. Bancroft (1985) provides an algebraic closed form solution of the overdetermined GPS pseudo-ranging observations. In this chapter we solved usingGroebner basis andMultipolynomial resultants GPS four-point pseudo-ranging problem. For literature on three-dimensional positioning models, we refer toE. W. Grafarend (1975),V. Ashkenazi andS, Grist (1982),G. W. Hein (1982a, 1982b),J. Zaiser (1986),F. W. O. Aduol (1989),U. Klein (1997) andS. M. Musyoka 1999).

InJ.L. Awange andE.W. Grafarend (2005) andJ. L. Awange et al. (2010), two case studies; the conversion of the geocentric GPS points for the Baltic Sea level Project into the Gauss–Jacobi ellipsoidal coordinates and the determination of the seven datum transformation parameters are presented. Datum transformation models have been dealt with e.g. inE. Grafarend, F. Okeke (1998),F. Krumm andF. Okeke (1995), E. Grafarend andF. Okeke (1998) andE. Grafarend andR. Syffus (1998).

In summary, this chapter has demonstrated the power of the algebraic computational tools ofGroebner bases and theMultipolynomial resultants in solving selected geodetic problems. In the case of theminimal combinatorial set, we have demonstrated how the problems ofminimum distance mapping (Example 15.16) and thepseudo-ranging four-point problem (Example 15.17) could be solved in a closed form using eitherGroebner bases approach or theMultipolynomial resultants approach. We have succeeded in demonstrating that by converting thenonlinear observation equations of the selected geodetic problems above into algebraic (polynomials), themultivariate system of polynomial equations relating the unknown variables (indeterminate) to the known variables can be reduced into polynomial equations consisting ofa univariate polynomial oncelexicographic monomial ordering is specified. We have therefore managed to providesymbolicsolutions to the problems of 7-parameter transformation,minimum distance mapping and theGPS pseudo-ranging four-point P4P by obtaining in each casea univariate polynomial that can readily be solved numerically once the observations are available. Although the algebraic techniques ofGroebner bases andMultipolynomial resultants were applied to these selected geodetic problems, the tools can be used to solve explicitly any closed form problem in Geodesy. The only limitation may be the storage and computational speed of the computers when compounded with problems involving many variables. The ability of the algebraic tools (Groebner bases and theMultipolynomial resultants) to solve closed form solutions gave theGauss–Jacobi combinatorial algorithm the required operational engine as evidenced in the solution of theoverdetermined GPS pseudo-ranging problem inJ.L . Awange andE.W. Grafarend (2005) andJ. L . Awange et al. (2010) where the solutions are obtained without reverting to iterative or linearization procedures. The results ofJ.L. Awange andE.W. Grafarend (2005) andJ. L . Awange et al. (2010) compared well with the solutions obtained using thelinearized Gauss–Markov model giving legitimacy to theGauss–Jacobi combinatorial procedure.

For the case studies (see e.g.,J.L. Awange andE.W. Grafarend (2005) andJ. L . Awange et al. (2010)), theGroebner basesalgorithm successfully determines explicitly the 7-parameters of the datum transformation problem and shows the scale parameter to be represented by aunivariate polynomial offourth order while the rotation parameters are represented by linear functions comprising the coordinates of the two systems and the scale parameter. The admissible value of the scale is readily obtained by solving theunivariate polynomial and restricting the scale to a positive real value \({x}_{1} \in \mathbb{R}{.}^{+}\) This eliminates the negative components of the roots and the complex values. The admissible value \({x}_{1} \in {\mathbb{R}}^{+}\) of the scale parameter is chosen and substituted in the linear functions that characterize the three elements of theskew-symmetric matrix \(S\) leading to the solution of the elements of the rotation matrix \({X}_{3}\). The translation elements are then deduced from the transformation equation. The advantage of using theGroebner basesalgorithm is the fact that there exists no requirement for prior knowledge of the approximate values of the 7-transformation parameters as is usually the case.

TheGroebner basesalgorithm managed to solve theMinimum Distance Mapping problem and in so doing, enabling the mapping of points from thetopographical surface to theInternational Reference Ellipsoid of Revolution. Theunivariate polynomial obtained was identical to that obtained by Grafarend and Lohse (1991). This implies that the algebraic tools ofGroebner bases and theMultipolynomial resultants can also be used to check the validity of existing closed form procedures in Geodesy.

TheGauss–Jacobi combinatorial algorithm highlighted one important fact while solving the overdetermined 7 parameter transformation problem; that the stochasticity of both systems involved can be taken into account. This has been the bottleneck of the problem of7-datum transformation parameters.