Computational Study of  the 3D Affine Transformation 
Part II. N-point Problem

Bela (Palancz^1), Piroska Zaletnyi(k^2) ,  Robert H. (Lewis^3) and Joseph Awang(e^4)

1Department of Photogrammetry and Geoinformatics
Budapest University of Technology and Economy,H-1521,Hungary
e-mail:palancz@epito.bme.hu

 2Department of Geodesy and Surveying and
Research Group of Physical Geodesy and Geodynamics of the Hungarian Academy of Sciences
Budapest University of Technology and Economy,H-1521,Hungary
e-mail: zaletnyikp@hotmail.com

(^3)Department of Mathematics
Fordham University, Bronx, NY 10458, USA
e-mail: rlewis @fordham.edu

4Western Australian Centre for Geodesy
Department of Spatial Sciences,Division of Resource and Environmental
Curtin University of Technology,GPO Box U1987, Perth,WA 6845 Australia
e-mail: J.awange@curtin.edu.au

Abstract

In case of  considerable nonlinearity. i.e. in geodesy, photogrammetry, robotics, it is difficult to find proper initial values to solve the parameter estimation problem of 3D affine transformation with 9 parameters  via linearization and/or iteration. In this paper we developed a symbolic - numeric method to achieve the solution without initial guess. Our method employs explicit analytical expressions developed by computer algebra technique  via Dixon resultant as well as reduced Groebner basis for solving 3 points problem, see in Part I. Then this solution can be used as initial value for a Newton -Raphson- Krylov numerical  method to solve the N points problem in a form of determined system of  9 polynomials developed by symbolic computation. 
The suggested method is fast, robust and it has low complexity comparing with other global as well as local numerical methods,  like direct minimization via genetic algorithm, linear homotopy continuation method, Newton-Raphson method with deflation employing SVD algorithm. Numerical illustration is presented with real world geodesical data representing 1138 points of  Hungarian Datum, OGPSH network.

Keywords:   9-parameter 3D affine transformation, solution of polynomial system,  local numerical methods, Newton-Raphson method,  deflation,  global minimization, genetic algorithm.

1. Introduction

In case of  considerable nonlinearity. i.e. in geodesy, photogrammetry, robotics, it is difficult to find proper initial values to solve the parameter estimation problem of 3D affine transformation with 9 parameters  via linearization and/or iteration.  The purpuse of this work, is to suggest and demonstrate a  symbolic-numeric  method to solve the 9 parameter 3D affine transformation problem  without initial guess for the case when N stations are known in both systems,i.e. Hooijberg (2008).

In geodesy, to solve the 3D Helmert transformation problem with 7 parameters,  there are different types of methods. The traditional numerical methods use linearization and linear regression subsequently, which needs inital values and iteration. One of the recent methods solves the 3 point problem in symbolic form  with reduced Groebner basis and the employing Gauss -Jacobi combinatoric algorithm with complicated weighting technique,i.e. Awange and Grafarend (2005). This method does not need initial guess, but it has two important disanvantages in case of real world problems: 1) increasing number of data points, N causes computational run-away, i.e. in case of N =100 we need 161 700 evalutions of the 3 point models; 2) real data are not well-behaved, different 3 points give very different parameter values and their proper weighting is quite hopeless, see in Part I. (Palancz et al 2008). Another, somewhat new and effective method the Procrustean transformation, which also alleviates the need for linearization and starting values, see Awange and Grafarend (2005). Recently, this method has been improved and extended to 9-parameter affine transformation, Awange, Bae and Claessens (2008), but its validity is restricted. 

 N points problem is basically an overdetermined problem consisting of  3N equations. and because of  the size of real world problems,  it is  not proper for symbolic solution, although some new research results have been published recently, for example Szanto(2007),  Therefore the solution of the overdetermined model was investigated via two numerical methods, with  global minimization using genetic algorithm and with a local method,  the Newton -Raphson method with deflation employing singular value decomposition to compute pseudoinverse of the Jacobian. This latest can be considered as the modern version of the linearization followed by linear regression method labeled as traditional method above, see Quoc-Nam Tran (1994) and Yhao (2007). 

In order to reduce the size and make the system square, the overdetermined modell can be transformed into a determined one by employing symbolic evalution of the least square objective and its symbolic derivation providing the necessary condition of its minimum. However, in this case one will get a higher order and more complex polynomial system than the original overdetermined one, but in that case we always have a system of  9 equations independently on the number of the data points.
The symbolic solution of this determined model proved to be hopeless because of its complexity, even the solution with coefficient in form of rational numbers seems to be very difficult. 

In case of the determined model of the N points problem,  the gobal numerical methods like  the genetic algorithm and  the Jenkins-Traub method failed. Linear homotopy could be used only with fixed point homotopy, while no natural starting system can be extracted from the original target system. Consequently the local numerical methods can be considered as candidates. That is why Newton-Raphson method with Krylov iteration algorithm has been chosen in order to recondition the Jacobian in case of ill-conditioned situation, Stoer and Bulirsch (2005).
 
Our suggested method employs explicit, analytical expressions developed by computer algebra technique  via Dixon resultant and reduced Groebner basis for solving 3 point problemin Part I. Then the result of the 3 point problem can be used as a good initial guess for a Newton -Raphson- Krylov numerical  method to solve the N point problem in a form of determined system of 9 polynomials developed by symbolic computation. 

This  method is very fast, robust and has low complexity comparing with other global as well as local numerical methods,  like direct minimization via genetic algorithm, linear homotopy continuation method and Newton-Raphson method with deflation employing pseudoinverse computed via singular value decomposition. Numerical illustration is presented with real world geodetic coordinates representing 1138 points of  Hungarian Datum, OGPSH network.
                                                                    
Numerical computations were carried out with hp workstation xw 4100 with XP operation system, 3 GHz P4 Intel processor and 1 GB RAM.

2. Definition of the N points problem

The (Subscript[C, 9])(3,3) 3D affine transformation is one possible generalization of the (Subscript[C, 7])(3,3) Helmert transformation, using three different scale ((Subscript[s, 1]), (Subscript[s, 2]), (Subscript[s, 3])) parameters instead of a single one. In this case the scale factors can be modeled by a diagonal matrix (W).

(({{(Subscript[x, i])}, {(Subscript[y, i])}, {(Subscript[z, i])}}))=W R (({{(Subscript[X, i])}, {(Subscript[Y, i])}, {(Subscript[Z, i])}}))+(({{(Subscript[X, 0])}, {(Subscript[Y, 0])}, {(Subscript[Z, 0])}}))

where W=(({{(Subscript[s, 1]), 0, 0}, {0, (Subscript[s, 2]), 0}, {0, 0, (Subscript[s, 3])}}))- scale matrix, (Subscript[X, 0]), (Subscript[Y, 0]), (Subscript[Z, 0]) - 3 translation parameters and R the rotation matrix

The rotation matrix can be expressed with the skew-symmetric matrix (S) (see Awange and Grafarend (2003))

S=(({{0, -c, b}, {c, 0, -a}, {-b, a, 0}}));

The rotation matrix is

R=((((Subscript[I, 3])-S))^(-1))((Subscript[I, 3])+S),

where (Subscript[I, 3]) is a 3×3 identity matrix.

(Subscript[I, 3])=IdentityMatrix[3];

R=Inverse[((Subscript[I, 3])-S)] .((Subscript[I, 3])+S)//Simplify; MatrixForm[R]

for which the next restriction is true:

R.Transpose[R]//Simplify//MatrixForm

let us call (Subscript[σ, i]) =(1/(Subscript[s, i]))and Ω = (W^(-1))

Ω=(({{(Subscript[σ, 1]), 0, 0}, {0, (Subscript[σ, 2]), 0}, {0, 0, (Subscript[σ, 3])}}));

For the determination of the 9 parameters of the transformation (a, b, c, (Subscript[X, 0]), (Subscript[Y, 0]), (Subscript[Z, 0]), (Subscript[s, 1]), (Subscript[s, 2]), (Subscript[s, 3])) we need at least 3 non-collinear points with known coordinates in both coordinate systems. In further,  instead of the scale parameters ((Subscript[s, 1]), (Subscript[s, 2]), (Subscript[s, 3])), we will use ((Subscript[σ, 1]), (Subscript[σ, 2]), (Subscript[σ, 3])) to get more simple equations. However in case of  N > 3 points in both systems, expressing the rotation matrix with the skew-symmetric matrix and using the inverse of the scale matrix (Ω) the nonlinear  system to be solved for ((Subscript[f, i]) =0 ) determining the 9 parameters leads to a system of  3 N polynomial equations. Now, in our numerical example, we have N = 1138 points.

𝒩=1138;

EQm=Table[(({{(Subscript[f, 1+3(j-1)])}, {(Subscript[f, 2+3(j-1)])}, {(Subscript[f, 3+3(j-1)])}}))=(((Subscript[I, 3])-S).Ω.(({{(Subscript[x, i])}, {(Subscript[y, i])}, {(Subscript[z, i])}}))- ((Subscript[I, 3])+S).(({{(Subscript[X, i])}, {(Subscript[Y, i])}, {(Subscript[Z, i])}}))-((Subscript[I, 3])-S).Ω.(({{(Subscript[X, 0])}, {(Subscript[Y, 0])}, {(Subscript[Z, 0])}})))/.i->j//Expand,{j,1,𝒩}];

sys=Table[(Subscript[f, i]),{i,3 𝒩}];Short[sys,10]

The number of the equations is

Length[sys]

This model is an overdetermined polynomial system with 9 unknowns and  3414 equations.

2.1 Numerical values of 1138 Hungarian points

Let us consider the numerical values of  1138  Hungarian points in the system of ETRS89 ((Subscript[x, 1]), (Subscript[y, 1]), (Subscript[z, 1]), ... , (Subscript[z, 3414]) ) and in the local Hungarian system OGPSH (Hungarian Datum ) ((Subscript[X, !]), (Subscript[Y, 1]), (Subscript[Z, 1]), ..., (Subscript[Z, 3414])).

The numerical values of the local and global coordinates are the following :

dataxyz=ReadList["G:\dataezer\IUGG67_ellipszoidi.txt",{Number,Number,Number}];

Short[dataxyz,10]

dataXYZ=ReadList["G:\dataezer\WGS84_ellipszoidi.txt",{Number,Number,Number}];

Short[dataXYZ,10]

dataxyzS=Drop[dataxyz,0];

dataXYZS=Drop[dataXYZ,0];

nxyz=Length[dataxyzS]

Now, we prepare the proper form of the data pairs

dataxyzT=Transpose[dataxyzS];

numericalValuesxyz=Table[{(Subscript[x, i])->dataxyzT[[1,i]],(Subscript[y, i])->dataxyzT[[2,i]],(Subscript[z, i])->dataxyzT[[3,i]]},{i,nxyz}]//Flatten;

dataXYZT=Transpose[dataXYZS];

numericalValuesXYZ=Table[{(Subscript[X, i])->dataXYZT[[1,i]],(Subscript[Y, i])->dataXYZT[[2,i]],(Subscript[Z, i])->dataXYZT[[3,i]]},{i,nxyz}]//Flatten;

numericalValues=Join[numericalValuesxyz,numericalValuesXYZ];

Short[numericalValues,10]

3. Solutions of the overdetermined model

In this section we solve this  system of  3 N = 3414  equations for the unknown parameters a, b, c, (Subscript[σ, 1]), (Subscript[σ, 2]),(Subscript[σ, 3]) , (Subscript[X, 0]),(Subscript[Y, 0])and (Subscript[Z, 0]). Two methods will be considered here.  Global minimization of the sum of the  square errors of the equations, which does not required initial guess values. Newton method with deflation using pseudoinverse of the Jacobian computed with singular value decomposition. It needs start values for the iteration, i.e. it is a local method.

3.1 Direct Numerical Solution via Global Minimization

The objective function can be constructed by considering the least square residiums of the equations, namely

obj=Apply[Plus,Map[((#)^2)&,sys]];Short[obj,10]

which should be minimized.

To carry out global minimization, the built -in function in Mathematica  NMinimize  based on genetic algorithm, can be employed.

nobj=obj/.numericalValues;

solN=NMinimize[nobj,{a,b,c,(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3]),(Subscript[X, 0]),(Subscript[Y, 0]),(Subscript[Z, 0])},Method->"DifferentialEvolution"]//Timing

3.2 Newton-Raphson with Deflation

This method solves the numerical Jacobi matrix in very steps with singular value decomposition method, therefore it can solve overdetermined system in least squares sense. However, this method requires initial guess values. These can be computed with the symbolic formulas developed for the 3 points problem. We use now the numerical values computed  on the basis of the first three points of the data set of 1138 points, with the analytic expressions developed  by computer algebra method  in Part I, see Table 1.

Table 1.
Result of the 3-point problem in case of N = 1138

({{a, -6.99099212136712*(10^(-6))}, {b, 8.020666005031663*(10^(-6))}, {c, 1.595856535443034*(10^(-6))}, {(Subscript[σ, 1]), 0.9999741620583953}, {(Subscript[σ, 2]), 0.9999995727567995}, {(Subscript[σ, 3]), 1.0000271383148969}, {(Subscript[X, 0]), -209.45964972433358752}, {(Subscript[Y, 0]), -1.39086256794054952}, {(Subscript[Z, 0]), 248.531507793257077310}})

ic={a->-6.99099212136712`*^-6,b->8.020666005031663`*^-6,c->1.595856535443034`*^-6,(Subscript[σ, 1])->0.9999741620583953`,(Subscript[σ, 2])->0.9999995727567995`,(Subscript[σ, 3])->1.0000271383148969`,(Subscript[X, 0])->-209.45964972433358752823551185429,(Subscript[Y, 0])->-1.39086256794054952123929069785,(Subscript[Z, 0])->248.53150779325707731004513334483}

X00={a,b,c,(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3]),(Subscript[X, 0]),(Subscript[Y, 0]),(Subscript[Z, 0])}/.ic;

The list of variables is

XX={a,b,c,(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3]),(Subscript[X, 0]),(Subscript[Y, 0]),(Subscript[Z, 0])};

Although built- in function FindMinimum works with Newton Raphson, but we should use here pseudoinverse, the algorithm of the Newton-Raphson with Deflation is implemented in Mathematica as , NewtonExtended , see Paláncz (2008).

nsys=sys/.numericalValues;

The implementation is a modified adoptation of  the function written by Ruskeepää (2004).

NewtonExtended[f_List,x_List,x0_List,eps_:10^-12,n_:100]:=With[{jac=N[Outer[D,f,x]]},FixedPointList[(#+PseudoInverse[jac/.Thread[x->N[#]]].(-f/.Thread[x->N[#]]))&,N[x0],n,SameTest->(Sqrt[Abs[(#1-#2).(#1-#2)]]<eps&)]]

NewtonExtended[nsys,XX,X00,10^-8,100]//Last//Flatten//Timing

Remark : This algorithm is faster by about two orders of magnitude than global minimization and it is quite robust, see the considerable deviations between the result and the initial guess values!

4. The determined model

Determined model can be achived by employing symbolic evalution of the least square objective and its symbolic derivation providing the necessary condition of its minimum. However, in this case one will get a higher order and more complex polynomial system than the original overdetermined one. In symbolic form most of these equations have many thousands of terms, therefore it is useful to collect terms corresponding to the same multivariate expression.

First let us write the 3 definable equations for the i - th point

eq1i=-(Subscript[X, i])+c (Subscript[Y, i])-b (Subscript[Z, i])+(Subscript[x, i]) (Subscript[σ, 1])-(Subscript[X, 0]) (Subscript[σ, 1])+c (Subscript[y, i]) (Subscript[σ, 2])-c (Subscript[Y, 0]) (Subscript[σ, 2])-b (Subscript[z, i]) (Subscript[σ, 3])+b (Subscript[Z, 0]) (Subscript[σ, 3]);

eq2i=-c (Subscript[X, i])-(Subscript[Y, i])+a (Subscript[Z, i])-c (Subscript[x, i]) (Subscript[σ, 1])+c (Subscript[X, 0]) (Subscript[σ, 1])+(Subscript[y, i]) (Subscript[σ, 2])-(Subscript[Y, 0]) (Subscript[σ, 2])+a (Subscript[z, i]) (Subscript[σ, 3])-a (Subscript[Z, 0]) (Subscript[σ, 3]);

eq3i=b (Subscript[X, i])-a (Subscript[Y, i])-(Subscript[Z, i])+b (Subscript[x, i]) (Subscript[σ, 1])-b (Subscript[X, 0]) (Subscript[σ, 1])-a (Subscript[y, i]) (Subscript[σ, 2])+a (Subscript[Y, 0]) (Subscript[σ, 2])+(Subscript[z, i]) (Subscript[σ, 3])-(Subscript[Z, 0]) (Subscript[σ, 3]);

The objective function for the i-th point is

obji=(eq1i^2)+(eq2i^2)+(eq3i^2)//Expand

The list of variables is

XX={a,b,c,(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3]),(Subscript[X, 0]),(Subscript[Y, 0]),(Subscript[Z, 0])};

The 9 equations of the determined model can be obtained by partial derivating the objective function with respect to the variables

{Eqai,Eqbi,Eqci,EqX0i,EqY0i,EqZ0i,Eqσ1i,Eqσ2i,Eqσ3i}=Map[D[obji,#]&,XX];

The first from the nine equations for the i-th point is obtained by partial differentiating the objective function with respect to the variable a ,

Eqai

To determine the nine equations to be solved for N points, the coefficients of the above determined nine equations need to be summarized. To do this we need to find the coefficients of the terms corresponding to the 9 variables (and their products) and to summarize them.

To find the coefficients we can use the next functions,

vari={a,b,c,(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3]),(Subscript[X, 0]),(Subscript[Y, 0]),(Subscript[Z, 0])};

Expon[equ_]:=Union[Table[Exponent[equ[[j]],Reverse[vari]],{j,Length[equ]}]]

Koeff[poly_,{i_,j_,k_,l_,m_,n_,o_,p_,q_}]:=CoefficientList[poly,Reverse[vari]][[i+1,j+1,k+1,l+1,m+1,n+1,o+1,p+1,q+1]]

KoeffList[poly_]:=Table[Koeff[poly,Expon[poly][[j]]],{j,Length[Expon[poly]]}]

ExpVari[{a1_,a2_,a3_,a4_,a5_,a6_,a7_,a8_,a9_},{i_,j_,k_,l_,m_,n_,o_,p_,q_}]:=a1^i a2^j a3^k a4^l a5^m a6^n a7^o a8^p a9^q

ExpVariList[poly_]:=Table[ExpVari[Reverse[vari],Expon[poly][[j]]],{j,Length[Expon[poly]]}]

For the first equation the summation and simplification can be done in the following steps:

First the equation can be simplified

gr1=Eqai/2//Simplify//Expand

The corresponding terms are

d1=ExpVariList[gr1]

and their coefficients

K1=KoeffList[gr1]

It is easy to check our result

gr1-K1.d1//Simplify

Now the summation of the coefficient can be done

C100=Map[(∑_(i=1)^n)#&,K1]

And the first equation of the determined model after summation for N point is

Eqa=C100.d1

All the equations of the determined model can be defined similarly.

Eqb=b n (Subscript[X, 0]^2) (Subscript[σ, 1]^2)-a n (Subscript[X, 0]) (Subscript[Y, 0]) (Subscript[σ, 1]) (Subscript[σ, 2])-c n (Subscript[Y, 0]) (Subscript[Z, 0]) (Subscript[σ, 2]) (Subscript[σ, 3])+b n (Subscript[Z, 0]^2) (Subscript[σ, 3]^2)+b (Subscript[X, 0]) (Subscript[σ, 1]^2) (∑_(i=1)^n)-2 (Subscript[x, i])+a (Subscript[Y, 0]) (Subscript[σ, 1]) (Subscript[σ, 2]) (∑_(i=1)^n)(Subscript[x, i])+b (Subscript[σ, 1]^2) (∑_(i=1)^n)(Subscript[x, i]^2)+b (Subscript[X, 0]) (Subscript[σ, 1]) (∑_(i=1)^n)-2 (Subscript[X, i])+(Subscript[Z, 0]) (Subscript[σ, 3]) (∑_(i=1)^n)-2 (Subscript[X, i])+a (Subscript[Y, 0]) (Subscript[σ, 2]) (∑_(i=1)^n)(Subscript[X, i])+b (Subscript[σ, 1]) (∑_(i=1)^n)2 (Subscript[x, i]) (Subscript[X, i])+a (Subscript[X, 0]) (Subscript[σ, 1]) (Subscript[σ, 2]) (∑_(i=1)^n)(Subscript[y, i])+c (Subscript[Z, 0]) (Subscript[σ, 2]) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[y, i])+a (Subscript[σ, 1]) (Subscript[σ, 2]) (∑_(i=1)^n)-(Subscript[x, i]) (Subscript[y, i])+a (Subscript[σ, 2]) (∑_(i=1)^n)-(Subscript[X, i]) (Subscript[y, i])+a (Subscript[X, 0]) (Subscript[σ, 1]) (∑_(i=1)^n)(Subscript[Y, i])+c (Subscript[Z, 0]) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[Y, i])+a (Subscript[σ, 1]) (∑_(i=1)^n)-(Subscript[x, i]) (Subscript[Y, i])+a (∑_(i=1)^n)-(Subscript[X, i]) (Subscript[Y, i])+b (Subscript[Z, 0]) (Subscript[σ, 3]^2) (∑_(i=1)^n)-2 (Subscript[z, i])+c (Subscript[Y, 0]) (Subscript[σ, 2]) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[z, i])+(Subscript[σ, 3]) (∑_(i=1)^n)2 (Subscript[X, i]) (Subscript[z, i])+c (Subscript[σ, 2]) (Subscript[σ, 3]) (∑_(i=1)^n)-(Subscript[y, i]) (Subscript[z, i])+c (Subscript[σ, 3]) (∑_(i=1)^n)-(Subscript[Y, i]) (Subscript[z, i])+b (Subscript[σ, 3]^2) (∑_(i=1)^n)(Subscript[z, i]^2)+b (Subscript[Z, 0]) (Subscript[σ, 3]) (∑_(i=1)^n)-2 (Subscript[Z, i])+c (Subscript[Y, 0]) (Subscript[σ, 2]) (∑_(i=1)^n)(Subscript[Z, i])+(Subscript[X, 0]) (Subscript[σ, 1]) (∑_(i=1)^n)2 (Subscript[Z, i])+(Subscript[σ, 1]) (∑_(i=1)^n)-2 (Subscript[x, i]) (Subscript[Z, i])+c (Subscript[σ, 2]) (∑_(i=1)^n)-(Subscript[y, i]) (Subscript[Z, i])+c (∑_(i=1)^n)-(Subscript[Y, i]) (Subscript[Z, i])+b (Subscript[σ, 3]) (∑_(i=1)^n)2 (Subscript[z, i]) (Subscript[Z, i])+b (∑_(i=1)^n)((Subscript[X, i]^2)+(Subscript[Z, i]^2));

Eqc=c n (Subscript[X, 0]^2) (Subscript[σ, 1]^2)+c n (Subscript[Y, 0]^2) (Subscript[σ, 2]^2)-a n (Subscript[X, 0]) (Subscript[Z, 0]) (Subscript[σ, 1]) (Subscript[σ, 3])-b n (Subscript[Y, 0]) (Subscript[Z, 0]) (Subscript[σ, 2]) (Subscript[σ, 3])+c (Subscript[X, 0]) (Subscript[σ, 1]^2) (∑_(i=1)^n)-2 (Subscript[x, i])+a (Subscript[Z, 0]) (Subscript[σ, 1]) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[x, i])+c (Subscript[σ, 1]^2) (∑_(i=1)^n)(Subscript[x, i]^2)+c (Subscript[X, 0]) (Subscript[σ, 1]) (∑_(i=1)^n)-2 (Subscript[X, i])+a (Subscript[Z, 0]) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[X, i])+(Subscript[Y, 0]) (Subscript[σ, 2]) (∑_(i=1)^n)2 (Subscript[X, i])+c (Subscript[σ, 1]) (∑_(i=1)^n)2 (Subscript[x, i]) (Subscript[X, i])+c (Subscript[Y, 0]) (Subscript[σ, 2]^2) (∑_(i=1)^n)-2 (Subscript[y, i])+b (Subscript[Z, 0]) (Subscript[σ, 2]) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[y, i])+(Subscript[σ, 2]) (∑_(i=1)^n)-2 (Subscript[X, i]) (Subscript[y, i])+c (Subscript[σ, 2]^2) (∑_(i=1)^n)(Subscript[y, i]^2)+(Subscript[X, 0]) (Subscript[σ, 1]) (∑_(i=1)^n)-2 (Subscript[Y, i])+c (Subscript[Y, 0]) (Subscript[σ, 2]) (∑_(i=1)^n)-2 (Subscript[Y, i])+b (Subscript[Z, 0]) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[Y, i])+(Subscript[σ, 1]) (∑_(i=1)^n)2 (Subscript[x, i]) (Subscript[Y, i])+c (Subscript[σ, 2]) (∑_(i=1)^n)2 (Subscript[y, i]) (Subscript[Y, i])+c (∑_(i=1)^n)((Subscript[X, i]^2)+(Subscript[Y, i]^2))+a (Subscript[X, 0]) (Subscript[σ, 1]) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[z, i])+b (Subscript[Y, 0]) (Subscript[σ, 2]) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[z, i])+a (Subscript[σ, 1]) (Subscript[σ, 3]) (∑_(i=1)^n)-(Subscript[x, i]) (Subscript[z, i])+a (Subscript[σ, 3]) (∑_(i=1)^n)-(Subscript[X, i]) (Subscript[z, i])+b (Subscript[σ, 2]) (Subscript[σ, 3]) (∑_(i=1)^n)-(Subscript[y, i]) (Subscript[z, i])+b (Subscript[σ, 3]) (∑_(i=1)^n)-(Subscript[Y, i]) (Subscript[z, i])+a (Subscript[X, 0]) (Subscript[σ, 1]) (∑_(i=1)^n)(Subscript[Z, i])+b (Subscript[Y, 0]) (Subscript[σ, 2]) (∑_(i=1)^n)(Subscript[Z, i])+a (Subscript[σ, 1]) (∑_(i=1)^n)-(Subscript[x, i]) (Subscript[Z, i])+a (∑_(i=1)^n)-(Subscript[X, i]) (Subscript[Z, i])+b (Subscript[σ, 2]) (∑_(i=1)^n)-(Subscript[y, i]) (Subscript[Z, i])+b (∑_(i=1)^n)-(Subscript[Y, i]) (Subscript[Z, i]);

EqX0=n (Subscript[X, 0]) (Subscript[σ, 1]^2)+(b^2) n (Subscript[X, 0]) (Subscript[σ, 1]^2)+(c^2) n (Subscript[X, 0]) (Subscript[σ, 1]^2)-a b n (Subscript[Y, 0]) (Subscript[σ, 1]) (Subscript[σ, 2])-a c n (Subscript[Z, 0]) (Subscript[σ, 1]) (Subscript[σ, 3])+(Subscript[σ, 1]^2) (∑_(i=1)^n)-(Subscript[x, i])+(b^2) (Subscript[σ, 1]^2) (∑_(i=1)^n)-(Subscript[x, i])+(c^2) (Subscript[σ, 1]^2) (∑_(i=1)^n)-(Subscript[x, i])+(b^2) (Subscript[σ, 1]) (∑_(i=1)^n)-(Subscript[X, i])+(c^2) (Subscript[σ, 1]) (∑_(i=1)^n)-(Subscript[X, i])+(Subscript[σ, 1]) (∑_(i=1)^n)(Subscript[X, i])+a b (Subscript[σ, 1]) (Subscript[σ, 2]) (∑_(i=1)^n)(Subscript[y, i])+c (Subscript[σ, 1]) (∑_(i=1)^n)-2 (Subscript[Y, i])+a b (Subscript[σ, 1]) (∑_(i=1)^n)(Subscript[Y, i])+a c (Subscript[σ, 1]) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[z, i])+a c (Subscript[σ, 1]) (∑_(i=1)^n)(Subscript[Z, i])+b (Subscript[σ, 1]) (∑_(i=1)^n)2 (Subscript[Z, i]);

EqY0=-a b n (Subscript[X, 0]) (Subscript[σ, 1]) (Subscript[σ, 2])+n (Subscript[Y, 0]) (Subscript[σ, 2]^2)+(a^2) n (Subscript[Y, 0]) (Subscript[σ, 2]^2)+(c^2) n (Subscript[Y, 0]) (Subscript[σ, 2]^2)-b c n (Subscript[Z, 0]) (Subscript[σ, 2]) (Subscript[σ, 3])+a b (Subscript[σ, 1]) (Subscript[σ, 2]) (∑_(i=1)^n)(Subscript[x, i])+a b (Subscript[σ, 2]) (∑_(i=1)^n)(Subscript[X, i])+c (Subscript[σ, 2]) (∑_(i=1)^n)2 (Subscript[X, i])+(Subscript[σ, 2]^2) (∑_(i=1)^n)-(Subscript[y, i])+(a^2) (Subscript[σ, 2]^2) (∑_(i=1)^n)-(Subscript[y, i])+(c^2) (Subscript[σ, 2]^2) (∑_(i=1)^n)-(Subscript[y, i])+(a^2) (Subscript[σ, 2]) (∑_(i=1)^n)-(Subscript[Y, i])+(c^2) (Subscript[σ, 2]) (∑_(i=1)^n)-(Subscript[Y, i])+(Subscript[σ, 2]) (∑_(i=1)^n)(Subscript[Y, i])+b c (Subscript[σ, 2]) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[z, i])+a (Subscript[σ, 2]) (∑_(i=1)^n)-2 (Subscript[Z, i])+b c (Subscript[σ, 2]) (∑_(i=1)^n)(Subscript[Z, i]);

EqZ0=-a c n (Subscript[X, 0]) (Subscript[σ, 1]) (Subscript[σ, 3])-b c n (Subscript[Y, 0]) (Subscript[σ, 2]) (Subscript[σ, 3])+n (Subscript[Z, 0]) (Subscript[σ, 3]^2)+(a^2) n (Subscript[Z, 0]) (Subscript[σ, 3]^2)+(b^2) n (Subscript[Z, 0]) (Subscript[σ, 3]^2)+a c (Subscript[σ, 1]) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[x, i])+b (Subscript[σ, 3]) (∑_(i=1)^n)-2 (Subscript[X, i])+a c (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[X, i])+b c (Subscript[σ, 2]) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[y, i])+b c (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[Y, i])+a (Subscript[σ, 3]) (∑_(i=1)^n)2 (Subscript[Y, i])+(Subscript[σ, 3]^2) (∑_(i=1)^n)-(Subscript[z, i])+(a^2) (Subscript[σ, 3]^2) (∑_(i=1)^n)-(Subscript[z, i])+(b^2) (Subscript[σ, 3]^2) (∑_(i=1)^n)-(Subscript[z, i])+(a^2) (Subscript[σ, 3]) (∑_(i=1)^n)-(Subscript[Z, i])+(b^2) (Subscript[σ, 3]) (∑_(i=1)^n)-(Subscript[Z, i])+(Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[Z, i]);

Eqσ1=n (Subscript[X, 0]^2) (Subscript[σ, 1])+(b^2) n (Subscript[X, 0]^2) (Subscript[σ, 1])+(c^2) n (Subscript[X, 0]^2) (Subscript[σ, 1])-a b n (Subscript[X, 0]) (Subscript[Y, 0]) (Subscript[σ, 2])-a c n (Subscript[X, 0]) (Subscript[Z, 0]) (Subscript[σ, 3])+(Subscript[X, 0]) (Subscript[σ, 1]) (∑_(i=1)^n)-2 (Subscript[x, i])+(b^2) (Subscript[X, 0]) (Subscript[σ, 1]) (∑_(i=1)^n)-2 (Subscript[x, i])+(c^2) (Subscript[X, 0]) (Subscript[σ, 1]) (∑_(i=1)^n)-2 (Subscript[x, i])+a b (Subscript[Y, 0]) (Subscript[σ, 2]) (∑_(i=1)^n)(Subscript[x, i])+a c (Subscript[Z, 0]) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[x, i])+(Subscript[σ, 1]) (∑_(i=1)^n)(Subscript[x, i]^2)+(b^2) (Subscript[σ, 1]) (∑_(i=1)^n)(Subscript[x, i]^2)+(c^2) (Subscript[σ, 1]) (∑_(i=1)^n)(Subscript[x, i]^2)+(b^2) (Subscript[X, 0]) (∑_(i=1)^n)-(Subscript[X, i])+(c^2) (Subscript[X, 0]) (∑_(i=1)^n)-(Subscript[X, i])+(Subscript[X, 0]) (∑_(i=1)^n)(Subscript[X, i])+(∑_(i=1)^n)-(Subscript[x, i]) (Subscript[X, i])+(b^2) (∑_(i=1)^n)(Subscript[x, i]) (Subscript[X, i])+(c^2) (∑_(i=1)^n)(Subscript[x, i]) (Subscript[X, i])+a b (Subscript[X, 0]) (Subscript[σ, 2]) (∑_(i=1)^n)(Subscript[y, i])+a b (Subscript[σ, 2]) (∑_(i=1)^n)-(Subscript[x, i]) (Subscript[y, i])+c (Subscript[X, 0]) (∑_(i=1)^n)-2 (Subscript[Y, i])+a b (Subscript[X, 0]) (∑_(i=1)^n)(Subscript[Y, i])+a b (∑_(i=1)^n)-(Subscript[x, i]) (Subscript[Y, i])+c (∑_(i=1)^n)2 (Subscript[x, i]) (Subscript[Y, i])+a c (Subscript[X, 0]) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[z, i])+a c (Subscript[σ, 3]) (∑_(i=1)^n)-(Subscript[x, i]) (Subscript[z, i])+a c (Subscript[X, 0]) (∑_(i=1)^n)(Subscript[Z, i])+b (Subscript[X, 0]) (∑_(i=1)^n)2 (Subscript[Z, i])+b (∑_(i=1)^n)-2 (Subscript[x, i]) (Subscript[Z, i])+a c (∑_(i=1)^n)-(Subscript[x, i]) (Subscript[Z, i]);

Eqσ2=-a b n (Subscript[X, 0]) (Subscript[Y, 0]) (Subscript[σ, 1])+n (Subscript[Y, 0]^2) (Subscript[σ, 2])+(a^2) n (Subscript[Y, 0]^2) (Subscript[σ, 2])+(c^2) n (Subscript[Y, 0]^2) (Subscript[σ, 2])-b c n (Subscript[Y, 0]) (Subscript[Z, 0]) (Subscript[σ, 3])+a b (Subscript[Y, 0]) (Subscript[σ, 1]) (∑_(i=1)^n)(Subscript[x, i])+a b (Subscript[Y, 0]) (∑_(i=1)^n)(Subscript[X, i])+c (Subscript[Y, 0]) (∑_(i=1)^n)2 (Subscript[X, i])+(Subscript[Y, 0]) (Subscript[σ, 2]) (∑_(i=1)^n)-2 (Subscript[y, i])+(a^2) (Subscript[Y, 0]) (Subscript[σ, 2]) (∑_(i=1)^n)-2 (Subscript[y, i])+(c^2) (Subscript[Y, 0]) (Subscript[σ, 2]) (∑_(i=1)^n)-2 (Subscript[y, i])+a b (Subscript[X, 0]) (Subscript[σ, 1]) (∑_(i=1)^n)(Subscript[y, i])+b c (Subscript[Z, 0]) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[y, i])+a b (Subscript[σ, 1]) (∑_(i=1)^n)-(Subscript[x, i]) (Subscript[y, i])+c (∑_(i=1)^n)-2 (Subscript[X, i]) (Subscript[y, i])+a b (∑_(i=1)^n)-(Subscript[X, i]) (Subscript[y, i])+(Subscript[σ, 2]) (∑_(i=1)^n)(Subscript[y, i]^2)+(a^2) (Subscript[σ, 2]) (∑_(i=1)^n)(Subscript[y, i]^2)+(c^2) (Subscript[σ, 2]) (∑_(i=1)^n)(Subscript[y, i]^2)+(a^2) (Subscript[Y, 0]) (∑_(i=1)^n)-(Subscript[Y, i])+(c^2) (Subscript[Y, 0]) (∑_(i=1)^n)-(Subscript[Y, i])+(Subscript[Y, 0]) (∑_(i=1)^n)(Subscript[Y, i])+(∑_(i=1)^n)-(Subscript[y, i]) (Subscript[Y, i])+(a^2) (∑_(i=1)^n)(Subscript[y, i]) (Subscript[Y, i])+(c^2) (∑_(i=1)^n)(Subscript[y, i]) (Subscript[Y, i])+b c (Subscript[Y, 0]) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[z, i])+b c (Subscript[σ, 3]) (∑_(i=1)^n)-(Subscript[y, i]) (Subscript[z, i])+a (Subscript[Y, 0]) (∑_(i=1)^n)-2 (Subscript[Z, i])+b c (Subscript[Y, 0]) (∑_(i=1)^n)(Subscript[Z, i])+b c (∑_(i=1)^n)-(Subscript[y, i]) (Subscript[Z, i])+a (∑_(i=1)^n)2 (Subscript[y, i]) (Subscript[Z, i]);

Eqσ3=-a c n (Subscript[X, 0]) (Subscript[Z, 0]) (Subscript[σ, 1])-b c n (Subscript[Y, 0]) (Subscript[Z, 0]) (Subscript[σ, 2])+n (Subscript[Z, 0]^2) (Subscript[σ, 3])+(a^2) n (Subscript[Z, 0]^2) (Subscript[σ, 3])+(b^2) n (Subscript[Z, 0]^2) (Subscript[σ, 3])+a c (Subscript[Z, 0]) (Subscript[σ, 1]) (∑_(i=1)^n)(Subscript[x, i])+b (Subscript[Z, 0]) (∑_(i=1)^n)-2 (Subscript[X, i])+a c (Subscript[Z, 0]) (∑_(i=1)^n)(Subscript[X, i])+b c (Subscript[Z, 0]) (Subscript[σ, 2]) (∑_(i=1)^n)(Subscript[y, i])+b c (Subscript[Z, 0]) (∑_(i=1)^n)(Subscript[Y, i])+a (Subscript[Z, 0]) (∑_(i=1)^n)2 (Subscript[Y, i])+(Subscript[Z, 0]) (Subscript[σ, 3]) (∑_(i=1)^n)-2 (Subscript[z, i])+(a^2) (Subscript[Z, 0]) (Subscript[σ, 3]) (∑_(i=1)^n)-2 (Subscript[z, i])+(b^2) (Subscript[Z, 0]) (Subscript[σ, 3]) (∑_(i=1)^n)-2 (Subscript[z, i])+a c (Subscript[X, 0]) (Subscript[σ, 1]) (∑_(i=1)^n)(Subscript[z, i])+b c (Subscript[Y, 0]) (Subscript[σ, 2]) (∑_(i=1)^n)(Subscript[z, i])+a c (Subscript[σ, 1]) (∑_(i=1)^n)-(Subscript[x, i]) (Subscript[z, i])+a c (∑_(i=1)^n)-(Subscript[X, i]) (Subscript[z, i])+b (∑_(i=1)^n)2 (Subscript[X, i]) (Subscript[z, i])+b c (Subscript[σ, 2]) (∑_(i=1)^n)-(Subscript[y, i]) (Subscript[z, i])+a (∑_(i=1)^n)-2 (Subscript[Y, i]) (Subscript[z, i])+b c (∑_(i=1)^n)-(Subscript[Y, i]) (Subscript[z, i])+(Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[z, i]^2)+(a^2) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[z, i]^2)+(b^2) (Subscript[σ, 3]) (∑_(i=1)^n)(Subscript[z, i]^2)+(a^2) (Subscript[Z, 0]) (∑_(i=1)^n)-(Subscript[Z, i])+(b^2) (Subscript[Z, 0]) (∑_(i=1)^n)-(Subscript[Z, i])+(Subscript[Z, 0]) (∑_(i=1)^n)(Subscript[Z, i])+(∑_(i=1)^n)-(Subscript[z, i]) (Subscript[Z, i])+(a^2) (∑_(i=1)^n)(Subscript[z, i]) (Subscript[Z, i])+(b^2) (∑_(i=1)^n)(Subscript[z, i]) (Subscript[Z, i]);

The list of the equations of the determined system

sys9={Eqa,Eqb,Eqc,EqX0,EqY0,EqZ0,Eqσ1,Eqσ2,Eqσ3};

The numerical form of the determined model substituting the 1138 Hungarian coordinates

numSys9=sys9/.n->nxyz/.numericalValues;

For example the first equation is,

numSys9[[1]]

4.1 Direct Numerical Solution via Global Minimization

The objective function can be constructed by considering the least square residiums of the equations, namely

nobjD=Apply[Plus,Map[((#)^2)&,numSys9]];

which should be minimized.

To carry out global minimization, the built -in function in Mathematica  NMinimize (#) based on genetic algorithm, can be employed.

solND=NMinimize[nobjD,XX,Method->"DifferentialEvolution"]//Timing

solND=NMinimize[nobjD,XX,Method->"DifferentialEvolution",MaxIterations->200]//Timing

Remark: The global minimization method without initial guess values failed to converge.

4.2 General polynomial solver

TimeConstrained[NSolve[numSys9,XX],300,"Not converging after 300 seconds"]

Remark: It seems to be very difficult to solve it with global method! We need to use local method for determined model, too.

4.3 Newton-Raphson Method

Here we use our own function newtonSolveSystem employing Krylov method for solving the linear system of the Jacobian, in order to be able to compare the running time with that of the overdetermined model solved by Newton deflation method. This function is a slightly modified version of the function written by Ruskeepää (2004).

newtonSolveSystem[f_List,x_List,x0_List,eps_:10^-8,n_:20]:=With[{jac=N[Outer[D,f,x]]},FixedPointList[(#+LinearSolve[N[jac/.Thread[x->N[#]],30],N[-f/.Thread[x->N[#]],30],Method->"Krylov"])&,N[x0],n,SameTest->(Sqrt[(#1-#2).(#1-#2)]<eps&)]]

Off[LinearSolve::"krymit"];

Employing initial guess computed for  3 points section, see Table 1.

newtonSolveSystem[numSys9,XX,X00,10^-8,100]//Last//Flatten//N//Timing

Let us summarize our results in case of the N point problem

Table 2.
Results in case of N point Problem

({{Model, Method, Initial guess value, Running Time
(s)}, {overdetermined, Global minimization
 with genetic algorithm, no, 130.91}, {overdetermined, Newton method with Deflation
using SVD, (SuperStar[yes]), 3.03}, {determined, Global minimization
 with genetic algorithm, no, failed}, {determined, General Polynomial Solver
, no, failed}, {determined, Newton Method
with Krylov , (SuperStar[yes]), 0.75}})

(SuperStar[])initial starting values are based on the result of the 3-point problem, taken the first three points of the 1138 points data set.

The traditional Newton - Raphson method without Krylov iteration is faster, but much more sensitive for the initial values. In our case it fails,

ini=SetPrecision[MapThread[{#1,#2}&,{XX,X00}],20]

FindRoot[numSys9,ini,WorkingPrecision->20]//Timing

5. Complexity study of the algorithms

In order to get more information about the algorithms can be considered as numerical methods for solving N point problem, the computations have been carried out for different N values, too. We consider here only the two Newton methods, namely Newton -Raphson method with deflation for overdetermined system and the Newton-Raphson with Krylov iteration for the determined model.

This implies that for N points  problem probably the best choice is employing the determined model with Newton -Raphson method with Krylov iteration. 

The initial values for this method can be computed from the 3 point model, but one need to be caushes and not to do it blindly. In the last section of Part I. (see Palancz et al 2008) we discussed how one can properly select 3 points from the N ones, in order to compute good initial values ensuring  fast convergence of  the Newton method employed for the solution of the N points problem.

Table 3.
Results  of  the complexity study

({{Number of points, Time
in case A
(s), Time 
in case B
(s)}, {200, 0.47, 0.78}, {400, 1.67, 0.81}, {600, 1.89, 0.69}, {800, 2.5, 0.73}, {1000, 2.78, 0.73}, {1138, 3.03, 0.75}})

where 
case A  -  stands for the Newton -Raphson method with deflation for overdetermined system
case B -  stands for the Newton -Raphson method with Krylov iteration for  the determined system

From Table 3. it is clear that Newton -Raphson method with Krylov iteration for  the determined system has zero complexity as it was expected. This means that computation time remains practically unchanged with increasing problem size.

Conclusions

According to our computational study of  different algorithms,  global minimization with genetic algorithm, Newton Raphson with deflation, Newton -Raphson method with Krylov iteration  as well as study of  different models, overdetermined and determined one, the determined model solved by Newton -Raphson method with Krylov iteration proved to be the most efficient technique, which being fast, robuts and almost zero complexity.

Acknowledgement

The authors are indebted to professor Daniel Lichtblau (Wolfram Research) for his valueable remarks and suggestions, especially calling up their attention to the proper implementation of NSolve as well as to the application of the reduced Groebner basis.

References

Awange J.L. and Grafarend E.W.: (2005) Solving Algebraic Computational Problems in Geodesy and Geoinformatics, Springer, Berlin, 

Awange, Bae and Claessens (2008) Procrustean solution of the 9-parameter transformation problem

Hooijberg M. (2008) Geometrical Geodesy, Springer, Berlin

Palancz B. (2008) Extended Newton-Raphson's Method, MathSource, http://library.wolfram.com/infocenter/MathSource/7118/

Palancz B., Lewis R.H., Zaletnyik P., Awange J.L. (2008): Computational Study of the 3D Affine Transformation - Part I. 3-point Problem, MathSource, http://library.wolfram.com/infocenter/MathSource/7090/

Ruskeepää H. (2004). Mathematical Navigator, Second Edition, Elsevier Academic Press, Amsterdam

Quoc-Nam Tran (1994) Extended Newton's method  for  finding the roots of an arbitrary system of nonlinear equations. In Hamza, M.H. ed. Proc. 12th IASTED Int. Conf.on Applied Informatics, IASTED, Anaheim, CA. 

Stoer and Bulirsch (2005) Numerische Mathematik 2. 5.Auflage Springer, Berlin

Szanto A. (2008) Solving over-determined systems by the subresultant method (with an appendix by Marc Chardin), J.of Symbolic Computation 43. pp.46-74  

Yhao (2007) Newton's Method with Deflation for Isolated Singularities of Polynomial Systems, Ph.D. Thesis. University of Illinois at Chicago, USA