Götterdämmerung over total least squares

J. Geod. Sci. 2016; 6:43ś60 Research Article Open Access G. Malissiovas*, F. Neitzel, and S. Petrovic Götterdämmerung over total least squares DOI 10.1515/jogs-2016-0003 Received November 11, 2015; accepted April 15, 2016 1 Introduction Abstract: The traditional way of solving non-linear least squares (LS) problems in Geodesy includes a linearization of the functional model and iterative solution of a nonlinear equation system. Direct solutions for a class of nonlinear adjustment problems have been presented by the mathematical community since the 1980s, based on total least squares (TLS) algorithms and involving the use of singular value decomposition (SVD). However, direct LS solutions for this class of problems have been developed in the past also by geodesists. In this contribution we attempt to establish a systematic approach for direct solutions of non-linear LS problems from a "geodetic" point of view. Therefore, four non-linear adjustment problems are investigated: the fit of a straight line to given points in 2D and in 3D, the fit of a plane in 3D and the 2D symmetric similarity transformation of coordinates. For all these problems a direct LS solution is derived using the same methodology by transforming the problem to the solution of a quadratic or cubic algebraic equation. Furthermore, by applying TLS all these four problems can be transformed to solving the respective characteristic eigenvalue equations. It is demonstrated that the algebraic equations obtained in this way are identical with those resulting from the LS approach. As a by-product of this research two novel approaches are presented for the TLS solutions of fitting a straight line to 3D and the 2D similarity transformation of coordinates. The derived direct solutions of the four considered problems are illustrated on examples from the literature and also numerically compared to published iterative solutions. For more than two centuries mathematicians and geodesists have solved overdetermined algebraic problems that often occur in the mathematical modelling of measurement results using the method of least squares (LS), see [1]. The simplicity of the łrecipež of the LS adjustment is recognised by its wide application in geodesy, a science that traditionally deals with redundant observations whilst seeking an łoptimalž solution. In geodetic literature the method of LS is mostly applied in the form of one of the two mainly referred adjustment models, namely the Gauss-Markov Model (GMM), see [2, p.137ff.], and the Gauss-Helmert Model (GHM), see e.g. [2, p.172ff.]. The two models can be found in [3, pp.7-26] as well, under the name parametric (case) adjustment and combined (case) adjustment respectively. Both models have as common feature the linearization of the functional model and the iterative solution of a non-linear equation system. In contrast to the classic LS solutions, the so called total least squares (TLS) solution for LS problems within the errors in variables (EIV) model started to be considered in the scientific community in the last three decades. It should be mentioned that the term EIV is mainly used by the statistical community for a special case of non-linear LS problems, e.g. a definition is given in [4] or [5, p.5]. Obviously the traditional geodetic approach can easily handle problems within the EIV, as it is presented for example in [3, p.10], [6, p.102] or [7]. Moreover, an attempt to introduce the TLS approach as a generalization of (ordinary) LS was udertaken in [8]. Within literature the solutions coming from TLS are often distinguished from the classical LS by stating that TLS functions differently due to a contra- Keywords: direct solution; eigenvalue problem; nonlinear least squares; planar similarity transformation; plane fit; singular value decomposition; straight line fit; total least squares *Corresponding Author: G. Malissiovas: Institute of Geodesy and Geoinformation Science, Technische Universität Berlin, Strasse des 17. Juni 135, 10623 Berlin, Germany, E-mail: [email protected] F. Neitzel: Institute of Geodesy and Geoinformation Science, Technische Universität Berlin, Strasse des 17. Juni 135, 10623 Berlin, Germany S. Petrovic: Institute of Geodesy and Geoinformation Science, Technische Universität Berlin, Strasse des 17. Juni 135, 10623 Berlin, Germany and GFZ German Research Centre for Geosciences, Section 1.2, Telegrafenberg, 14473 Potsdam, Germany © 2016 G. Malissiovas et al., published by De Gruyter Open. This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivs 3.0 License. 44 | G. Malissiovas et al. diction adjustment model, see e.g. [9, 10], which can be formulated as l = (A − EA ) x + el , dim(A) = n × m, rank(A) = m < n, (1) where l and el are the vectors of observations and their random errors, respectively. EA are random errors in the coefficient matrix A and x is the vector containing the unknown parameters. For a full rank TLS coefficient matrix A in (1) the number of observations n is larger than the unknown parameters m. By the definition of TLS [5, p. 33] the TLS solution is based on the minimization of the objective function ^ , ^l] ||F → min, || [A, l] − [A (2) with || ||F being the Frobenius norm of a matrix. The ad^ and ^l are justed matrices A ^ , ^l] = [A, l] − [EA , el ]. [A (3) The objective function under minimization can be also expressed as in [11] by el T el + vec(EA )T vec(EA ) → min, (4) with łvecž indicating a function that stacks the columns of the error matrix EA into one vector. In this study, the terms TLS and TLS solution will refer to the solution of model (1) by minimizing the objective function (4) through SVD of the augmented matrix [A, l] (i.e. the matrix containing the coefficient matrix A and the observation vector l). In the literature there are expectations that TLS might produce a łmore realisticž solution than the classic LS as indicated for example in [11, 12]. This has been caused possibly by the work of Golub and Van Loan [9], where the solution of TLS was compared with that of LS for fitting a straight line in 2D. In the latter study, the LS solution was assumed when only the y-coordinates were regarded as observed values and the x-coordinates as error free, which led to the misleading conclusion that TLS functions differently from LS. The absence of any work from the geodetic literature is noticeable. For geodesists it has been already clear that the most important steps for the adjustment of observations is to build a correct model and minimize the errors of a correct objective function. When these requirements are fulfilled then for a linear problem the solution will be unique, regardless of the solution strategy that has been followed, e.g. TLS, GHM or some other approach. Contrary to the belief that TLS is an additional method like LS (or even a generalisation of it), Neitzel and Petrovic [13] showed that in fact TLS can be regarded as a special case of the LS method within the GHM for the example of fitting a straight line to equally weighted two-dimensional (2D) data. The iterative solution of the GHM has been proven to be equivalent to the TLS solution. Thus is any discussion unnecessary, which of the two approaches is better. It is only necessary always to model the given problem and not something completely different from it. However, in contrast to GHM, the elegant solution of TLS was derived with no need for iterations or starting values and by making use of singular value decomposition (SVD). Moreover, Felus and Schaffrin [11] attempted to derive the TLS solution for the example of a 2D similarity transformation of coordinates using SVD. However, Neitzel [14] showed that their solution needed modifications and provided the correct TLS results by evaluating the solution obtained by an iteratively linearised GHM. Additionaly, an iterative TLS solution for the same problem was presented in [15]. Concerning the two latter LS problems, it has been shown that a TLS solution can be obtained equivalently from a linearized GHM, following the traditional geodetic procedure for solving non-linear LS problems. Hence, the following questions arise out of this: ś If it is possible to solve an adjustment problem with TLS and SVD, is it also possible to obtain the same eigenvalue problem from a geodetic point of view and solve the problem directly? ś Are there additional non-linear LS problems (besides the generally well-known case of the straight line fitting to equally weighted 2D data) which can be solved directly? ś Is it possible to classify those non-linear LS problems with a direct solution and solve them by using a systematic approach? These objectives gave the motivation for deeper and further investigation of the solutions obtained by using TLS. We want to develop a clear mathematical relationship between TLS and direct solutions of non-linear LS problems for the following four individual cases: 1. Fitting a straight line in 2D, 2. Fitting a straight line in 3D, 3. Fitting a plane in 3D, Götterdämmerung over total least squares | 45 4. 2D similarity transformation of coordinates. In all four cases under investigation the coordinates in all directions are regarded as measurements. In TLS literature these problems are often distinguished as EIV model. Moreover, we always postulate a regular adjustment model and the observations are equally weighted and uncorrelated. The concept of solving non-linear LS problems applied here is directly based on [16], where the non-linear problem of fitting a straight line to a set of points in 3D space was examined. In that work, the developed LS solution was obtained by solving an eigenvalue problem, which is one of the key elements of TLS as well. In this contribution, under the title "Götterdämmerung³ over total least squares", we want to establish a systematic approach for solving the four investigated adjustment problems, based on the solution of [16]. Our mathematical approach involves a sophisticated parametrization of the problem which can be always solved by building a Lagrange function that results in a quadratic or cubic algebraic equation. The equivalence between singular values and Lagrange multipliers has been analyzed already in [18, p. 44]. A flowchart presenting both ways of solving directly nonlinear LS problems is depicted in Fig. 1. This article is arranged as follows: ś In the first section the direct LS solution for the problem of fitting a straight line in 2D is approached from geodetic point of view. A TLS solution is provided as well, which sets the foundation for the mathematical relationship between the two individual estimates (i.e. the direct LS and the TLS). ś It follows a second section which studies the fitting of a straight line in 3D. Following the same procedure as in the first section two individual solutions are obtained, one for a direct LS and a novel TLS solution for fitting a straight line in space using SVD. ś The next two following sections are dealing with fitting a plane in a 3D point cloud and the 2D similarity transformation of coordinates. Here, a novel approach is presented for the TLS solution of the 2D transformation. 3 The downfall of the gods (in Germanic mythology), see Oxford dictionary. In the context of LS adjustment, the same term has been firstly used in [17]. LS Special cases of non-linear LS problems TLS Sophisticated parametrization Augmented matrix Lagrange function Singular value decomposition Characteristic polynomial Direct solution for the unknowns Figure 1: Flowchart for two possible direct solutions of non-linear LS problems. All investigated cases are compared with already existing algorithms or models (for example with LS solutions from the GHM), numerical examples are presented at the end of each section. It should be pointed out that direct solutions for the non-linear problem of fitting a straight line in 2D as well as for the case of fitting a plane in 3D can be found already in [19] and they coincide with those presented in this contribution. However, curiously, this work of Linkwitz [19] is rarely cited. 2 Fitting a straight line in 2D One of the first attempts to solve the non-linear problem of LS for fitting a straight line to a set of points in plane (i.e. in the 2D space) non-iteratively was done by Adcock [20] who provided an elegant way of finding the direct solution to the problem. Pearson [21] investigated the same problem by minimizing the sum of the squared orthogonal distances of every point to the requested line and he extended his study to fitting a plane to a set of points in space (i.e. in the 3D space) as well. On the other hand, more recently, the work of [9] provided an analysis of the TLS problem followed by the contributions of [8, 22ś24]. These authors always comprised 46 | G. Malissiovas et al. the example of the straight line fit as the most appropriate example for illustrating the idea of TLS. At the beginning an amount of 2D data is observed, e.g. a point cloud with coordinates in x and y-direction. The question is how to fit a straight line to the measured point cloud. A straight line in 2D is represented in coordinate form [25, p. 217] by di = (5) This line passes through a point with coordinates x0 and y0 , and is parallel to a direction vector with components a and b. Assuming that the observed points lie on this straight line and that the measurements are error-free (which is never the case), x and y can be regarded as observed coordinates of a point in 2D. (y i − y0 ) b − (x i − x0 ) a . a2 + b2 2D Points LS − line vx,vy 9 8 d2 = vx2 + vy2 −> min (LS) d 7 6 5 4 vy 3 2 2.1 LS line fit in 2D (7) 10 y − direction [m] y − y0 x − x0 = . a b The normal distance of every point to the requested line can be expressed by [25, p. 218] d (orthogonal distance) vx 1 In this section we try to develop a direct LS solution for fitting a straight line in 2D when both coordinates are subject to errors. The unknown line parameters can be estimated directly by constructing and minimizing an appropriate Lagrange function and by solving a system of linear normal equations. Our goal is to show that the proposed approach leads, according to the chosen technique for solving linear equation systems, to the solution of such algebraic equations that are equivalent to TLS (the solution for fitting a straight line in 2D by TLS is presented and analysed in the next section). 0 2 4 6 8 10 x − direction [m] Figure 2: Example of fitting a straight line to 2D points with both x and y coordinates subject to measurement errors. The least squares criterion can be applied to obtain the minimum normal distances of a point cloud to the fitted line by minimizing the sum of the squared normal distances n ∑︁ d2i → min. (8) i=1 2.1.0.1 Definition of the problem From a geodetic point of view it is of great importance to clarify from the beginning which quantities are observations and hence subject to random errors. This is necessary in order to define the target (or objective) function of the problem in an appropriate way. In this investigation let us assume that both coordinates (in x and y-direction) are subject to measurement errors. Furthermore, let all measurements be uncorrelated and of the same accuracy. Therefore, the aim is to find the shortest distance of each łmeasuredž point to an adjusted straight line. As noticed already in [20] the same accuracy of all coordinate measurements corresponds to the normal distance d2i = v2x i + v2y i , 0 (6) as measure of deviations, with i = 1, ..., n (n is the number of observed points). This problem is depicted in Fig. 2. There are infinitely many choices for a condition that connects the two unknown parameters a and b for the general equation of the straight line (5). We will not restrict the problem to the usual a = 1 or b = 1, for the reason that some lines in plane are excluded with these choices (e.g if we choose a = 1, then there is no solution for lines parallel to the x-direction). From all remaining restrictions, the most appropriate for us is a2 + b2 = 1, (9) as it allows all lines in the plane to be calculated. Geometrically this restriction can be seen as a normalization of the orthogonal distances from every point to the requested line (i.e. the denominator of the orthogonal distance of Eq. (7) becomes 1). Therefore, the objective function of this Götterdämmerung over total least squares Ω(a, b, y0 , x0 ) = = n ∑︁ d2i y0 = i=1 (10) 2 n n i=1 i=1 1 ∑︁ 1 ∑︁ y i and x0 = xi , n n y − yc x − xc = . a b 2.1.0.2 Elimination of two unknowns The elimination of the unknown parameters x0 and y0 is possible by showing that they can be taken as equal to the coordinates of the centre of mass yc = 1 n y i and x c = i=1 1 n n ∑︁ xi (11) i=1 of the 2D points, which are regarded as observations. This has already been proven in [20], where it was shown that the requested line passes through the centre of mass of the 2D point set. In [16] the 3D case is covered. Here a similar proof for the case of the 2D line fit is presented in short. Multiplying all terms of the objective function of equation (10) results in Ω(a, b, y0 , x0 ) = A y20 + B y0 + C 1 n n ∑︁ i=1 +2abx0 i=1 n ∑︁ Equation (18) can be further simplified by reducing the coordinates to a coordinate system with its origin in the centre of mass of the given points. The reduced coordinates of a point can be described by y′ = y − y c and x′ = x − x c , (19) y′ b = x′ a. (20) which leads to Hence, the normal distance of every reduced point to the requested line can be rewritten as di = y′i b − x′i a . a2 + b2 (21) The reduction means that the investigated straight line passes through the centre of mass which has been translated to the centre of the coordinate system, as depicted in Fig. 3. Consequently, the objective function (10) can be rewritten (with eliminated unknowns x0 , y0 ) as A =)︃nb2 , (︃ )︃]︃ n 1 ∑︁ y i + ab x i − x0 , B = −2n b n i=1 i=1 n n n ∑︁ ∑︁ ∑︁ C = a2 (x i − x0 )2 + b2 y2i − 2ab xi yi 2 (︃ (18) 2.1.0.3 LS line fit with reduced coordinates (12) with [︃ (17) satisfy the equations (15) and (16). Thus the equation of a line in 2D can be reformulated as ((y i − y0 ) b − (x i − x0 ) a) . i=1 n ∑︁ 47 It is easy to check by substitution that problem can be expressed by n ∑︁ | (13) Ω(a, b) = i=1 yi . = b i=1 2 n ∑︁ n ∑︁ d2i = i=1 2 y′i − 2ab i=1 n ∑︁ (y′i b − x′i a)2 i=1 n ∑︁ y′i x′i + a2 i=1 n ∑︁ (22) 2 x′i . i=1 Assuming that the parameters a, b and x0 are known, functions A, B and C become constant with A being positive. Therefore Eq. (12) would represent a parabola with a minimum at B y0 = − . (14) 2A We seek for a least squares solution for the unknown parameters a and b that minimizes Eq. (22) subject to the constraint (9). Therefore, we can introduce the Lagrangian Inserting the expressions from (13) into (14) yields (︃ )︃ n n a 1 ∑︁ 1 ∑︁ yi + x i − x0 . y0 = n b n where k is the Lagrange multiplier. By differentiating function K with respect to the unknown parameters a and b and setting the partial derivatives to zero we obtain i=1 Analogously it can be shown that )︃ (︃ n n 1 ∑︁ b 1 ∑︁ x0 = xi + y i − y0 . n a n i=1 (15) i=1 i=1 K(a, b, k) = Ω(a, b) − k(a2 + b2 − 1), (︂ ∑︁ )︂ (︂ ∑︁ )︂ n n ∂K 2 =a x′i − k − b y′i x′i = 0, ∂a i=1 (16) (24) i=1 )︂ (︂ ∑︁ )︂ (︂ ∑︁ n n ∂K ′2 ′ ′ = −a y i − k = 0, yi xi + b ∂b i=1 (23) i=1 (25) 48 | G. Malissiovas et al. y − direction [m] 5 4 2D Points LS − line vx,vy 3 d 5 d2 = v 2 + v 2 −> min (LS) x 4 3 2 2 1 1 0 −5 0 0 5 −1 −1 v −2 −2 −3 −3 −4 −4 d (orthogonal distance) y −5 vx −5 −5 2.2 TLS line fit in 2D y 0 5 x − direction [m] An alternative solution for finding the line fitting to a set of points in 2D can be provided by TLS [8, 9]. According to [9] the solution of TLS can be represented geometrically by minimizing the orthogonal distances as it is depicted in Fig. 3. It is noteworthy that in the latter contribution the LS problem of fitting a straight line in 2D was regarded only when the y-coordinates are observations, in contrast to the definition and solution of the LS problem that was presented in the previous section. Thus, our target is to provide a clear insight of the TLS approach and show that it is equivalent to the developed approach of section 2.1. Rearranging the functional model of Eq. (19), which already contains the reduced coordinates of the measured points to the centre of mass, yields Figure 3: Example of fitting a straight line to a set of 2D points with coordinates reduced to the centre of mass. y′ = −β x′ (28) with a . (29) b It is to be noted that using the functional model of Eq. (28) is equivalent with employing a restriction of a = 1. Therefore, it is not possible to describe all straight lines in plane, however, these are limited cases (in this case all lines that are parallel to the y-direction). Nevertheless, it is possible to build the TLS model of Eq. (1) with the respective quantities ⎤ ⎡ ⎡ ⎤ −e x′1 −x′1 ⎢ −e ′ ⎥ ⎢ −x′ ⎥ ⎢ ⎢ x2 ⎥ 2 ⎥ ⎥ ⎢ ⎥ A = ⎢ . ⎥ , EA = ⎢ ⎢ .. ⎥ , . ⎣ . ⎦ ⎣ . ⎦ −x′ −e x′n ⎡ ′ ⎤ n ⎡ ⎤ (30) e y′1 y1 ⎢ e ′ ⎥ ⎢ y′ ⎥ )︁ (︁ ⎢ y2 ⎥ ⎢ 2 ⎥ ⎥ ⎢ ⎥,x ^ = β^ . l = ⎢ . ⎥ , el = ⎢ . ⎢ ⎥ ⎣ .. ⎦ ⎣ .. ⎦ ′ yn e y′n −β = subject to a2 + b2 = 1. If the Lagrange multiplier were known, then Eq. (24) and (25) would represent a homogeneous system of equations which is linear in the unknown line parameters a and b. For a nontrivial solution the determinant of the equation system must be equal to zero: ⃒ (︂ n )︂ ⃒ ∑︁ 2 ′ ⃒ x − k i ⃒ ⃒ i=1 ⃒ ⃒ ⃒ n ⃒ ∑︁ ⃒ − y′i x′i ⃒ ⃒ i=1 n ∑︁ − y′i x′i ⃒ ⃒ ⃒ ⃒ ⃒ i=1 ⃒ ⃒ = 0, ⃒ (︂∑︁ )︂ n ⃒ 2 ⃒ y′i − k ⃒ ⃒ (26) i=1 which leads to the quadratic equation (︂ ∑︁ n i=1 2 x′i − k )︂(︂ ∑︁ n i=1 )︂ (︂ ∑︁ )︂2 n 2 = 0, y′i − k − y′i x′i (27) i=1 with one unknown parameter k and two real and positive solutions k min and k max . It can be shown that the smaller of the two solutions for k, denoted by k min , corresponds to the minimum of the objective function (22). The solution for the adjusted unknown parameters a and b can be computed by substituting the Lagrangian factor k min into equations (24)-(25) subject to the constraint (9). An equivalent direct solution for the discussed nonlinear problem has already been presented by Linkwitz [19]. In his contribution the Hessian normal form is employed for the representation of the line in plane and leads to the same quadratic equation. The solution is obtained also by solving an eigenvalue problem. Matrix A contains the coefficients of Eq. (28) with respect to the unknown parameter β. 2.2.0.1 The minimum eigenvalue principle The TLS solution has been presented amongst others by Felus and Schaffrin [11], by employing the minimum eigenvalue principle. The first step is to construct the augmented matrix ⎡ ⎤ −x′1 y′1 ⎢ −x′ y′ ⎥ ⎢ 2 2 ⎥ [A, l] = ⎢ (31) .. ⎥ ⎢ .. ⎥. ⎣ . . ⎦ −x′n y′n Götterdämmerung over total least squares By decomposing the augmented matrix [A, l] with the help of SVD yields UΣVT = SVD([A, l]), v = [v1,m+1 , · · · , v m,m+1 , v m+1,m+1 ]T , (33) where m is the number of unknown parameters. Thus, the TLS solution for the adjusted vector of unknowns is ^= x 1 [v :v ]T . v m+1,m+1 1,m+1 m,m+1 (34) It must be mentioned that in TLS literature Eq. (34) is presented with a negative sign as in [11], which is caused by the form of the functional model. In this study the functional model (Eq. 28) has been developed in such a way that the negative sign is not necessary anymore. 2.2.0.2 Connection with the eigenvalue/eigenvector decomposition To understand deeper the operation of SVD and the derivation of the adjusted unknowns of Eq. (34) it is important to connect SVD with the eigenproblem of the symmetric non-negative definite matrices ([A, l]T [A, l]) and ([A, l][A, l]T ). According to [26, p.18] matrix V containing the right singular vectors of [A, l] can also be estimated from the eigenvalue decomposition (EVD) of the squared matrix ([A, l]T [A, l]) by VΛVT = EVD([A, l]T [A, l]), (35) where matrix Λ is a diagonal matrix carrying the eigenvalues of [A, l]. According to [27, p.427] a relationship between eigenvalues and singular values can be expressed as λi = σi 2 , By = λy ⇒ (B − λI)y = 0, (38) with I denoting an identity matrix and y an eigenvector of B. The eigenvalues of matrix B can be determined by searching for non-trivial solutions y ̸= 0, i.e. by solving the characteristic equation of the eigenvalues ⃒ n ⃒ n ⃒ ∑︁ 2 ⃒ ∑︁ ′ ′ ′ ⃒ ⃒ xi − λ − yi xi ⃒ ⃒ ⃒ i=1 ⃒ i=1 ⃒ ⃒ ⃒ ⃒ = 0, (39) ⃒ ⃒ n n ⃒ ∑︁ ⃒ ∑︁ 2 ⃒ ⃒ y′i − λ ⃒ y′i x′i ⃒ − ⃒ ⃒ i=1 i=1 i.e. (︂ ∑︁ n i=1 2 x′i − λ )︂(︂ ∑︁ n i=1 )︂2 )︂ (︂ ∑︁ n 2 = 0. y′i − λ − y′i x′i (40) i=1 This is a quadratic equation with two solutions for the unknown eigenvalues λ min and λ max . By rearranging the eigenvalues and eigenvectors appropriately, the TLS solution for the adjusted line parameter β can be found from Eq. (34). Thus, by normalizing the eigenvector that corresponds to the smallest eigenvalue yields − β= n ∑︁ y′i x′i 2 x′i − λ min i=1 n ∑︁ . (41) i=1 As expected the TLS solution for the line parameters is identical to the developed LS solution. The equivalence between the two solutions is in Eq. (29) and using the minimum Lagrange multiplier k min in Eq. (24). Furthermore, the developed characteristic equation of the eigenvalues (40) corresponds to the quadratic equation (27) from the direct LS solution. The conclusion is that the presented direct LS solution for the non-linear straight line fit in plane already provides the exact result for TLS. (36) with λ and σ being the eigenvalues and singular values, respectively. For an explicit solution of the eigenproblem for the straight line fit we build the squared matrix ⎤ ⎡ n n ∑︁ ∑︁ ′2 ′ ′ x x − y i i i ⎥ ⎢ ⎥ ⎢ i=1 i=1 ⎥ ⎢ T ⎥ = B. ⎢ (37) [A, l] [A, l] = ⎢ ⎥ n n ⎥ ⎢ ∑︁ ∑︁ 2 ⎦ ⎣ y′i − y′i x′i i=1 49 The eigenvalues and eigenvectors of matrix B can be computed from the generalised eigenvalue problem [25, p. 278] (32) where the matrices U and V are orthogonal and contain the left and right singular vectors of the augmented matrix respectively and matrix Σ is diagonal carrying the singular values. The TLS solution can be derived by normalising the right singular vector of matrix V that corresponds to the minimum singular value (this is the last column of V) [5, p.35], with | i=1 2.3 Numerical example - line fit 2D In this section a numerical example is presented for validating the results of the current investigation. The example dataset that was chosen is the one used in [13] and can be found in Table 1. The measured coordinates are uncorrelated and both in x and y-direction under the influence of random errors. 50 | G. Malissiovas et al. Table 1: Example dataset for straight line fitting in 2D point no. x-coord. [m] y-coord. [m] 1 2 3 4 0 1 2 3 0 1 4 9 For solving the non-linear problem of fitting a straight line to the four points of the presented dataset, an iterative solution of a GHM has been proposed in [13] (however, in this case only one iteration was sufficient for a LS solution). The results of the GHM are listed in Table 2. Table 2: Results from [13] GHM solution ^) −β^ (Steigung a ^ γ^ (Achsabschnitt b) 3.241804 -1.362705 Secondly, the presented direct LS solution for the unknown line parameters a and b was obtained. Building the objective and Lagrange function of equations (22) and (23) leads to a homogeneous system of linear equations with ^ that can be estimated by solvone unknown parameter k, ing ⃒ ⃒ ⃒ (5 − k) −15 ⃒⃒ ⃒ (42) ⃒ ⃒ = 0. ⃒ −15 (49 − k) ⃒ This yields a quadratic equation with the solutions for the unknown parameters k min = 0.3729460886 and k max = 53.6270539113. The direct LS solution of the line parameters can be found in Table 3. The adjusted parameters a and b were utilized to compute β. Table 3: Direct LS solution - line fit 2D Estimated parameters LS solution ^ a ^ b 0.9555698150338 0.2947648700171 Derived parameters γ^ ^ a ^ b y = −β x + γ . (43) Finally, the TLS estimate was derived for the unknown line parameters. The determinant of the generalised eigenvalue problem was built following the procedure for the eigenvalue/eigenvector solution of matrix B (Eq. 39), which results in the characteristic equation of the eigenvalues (quadratic equation) with solutions for λ min = 0.3729460886 and λ max = 53.6270539113. The TLS solution for the adjusted line parameters is presented in Table 4 and is of course equal to the results in Table 3. Both solutions coincide numerically with the iterative result of the linearized GHM. Table 4: TLS solution - line fit 2D Estimated parameters −β^ = Eq. (28) for the straight line can be rewritten as 3.2418035940925 -1.3627053911388 Parameter γ can be easily derived by bringing the coordinate system into its original position. In this sense Estimated parameters TLS solution −β^ γ^ 3.2418035940925 -1.3627053911388 3 Fitting a straight line in 3D The problem of fitting a straight line to a set of 3D points has been examined e.g. in [28ś30]. The latter contributions provided iterative algorithms for minimizing the sum of the squared orthogonal distances of each 3D point to the fitted line and thus obtaining a least squares estimate for the unknown line parameters. Non-iterative adjustment solutions for the straight line fit in space can be found in [16, 31]. Both solutions can be transformed in an eigenvalue problem which gives the motivation for investigating the relationship with the TLS solution. A straight line in 3D can be represented in coordinate form by y − y0 x − x0 z − z0 = = , (44) a b c for a line that passes through a point with coordinates x0 , y0 and z0 and is parallel to a direction vector with components a, b and c. The target is to minimize the errors in all x, y and z coordinates, which implies the non-linearity of the problem. Similarly to the previous case it is possible to reduce the number of the unknowns of the model by replacing the parameters x0 , y0 and z0 with the coordinates Götterdämmerung over total least squares In order to derive the minimum solution of the objective function under the restriction of Eq. (49) the Lagrangian of the centre of mass yc = n ∑︁ 1 n 1 n yi , xc = i=1 zc = n 1 ∑︁ n n ∑︁ | 51 xi , i=1 K(a, b, c, k) = Ω(a, b, c) − k(a2 + b2 + c2 − 1) (45) zi , i=1 of the n given points. The proof that this parameter replacement is allowed, can be found in [16]. Therefore, Eq. (44) can be rewritten as y − yc x − xc z − zc = = . (46) a b c (52) can be built. A differentiation of function K with respect to the unknown line parameters a, b and c and setting the resulting partial derivatives to zero, leads to the equations Obviously a reduction of all coordinates to the centre of mass leads to the elimination of x c , y c and z c , which results in y′ x′ z′ = = , (47) a b c with x′ , y′ and z′ being coordinates reduced to the centre of mass. (︂ ∑︁ )︂ n n n ∑︁ ∂K 2 ∑︁ ′ 2 =a y′i x′i zi − k − b x′i + ∂a i=1 i=1 i=1 n ∑︁ ′ ′ −c y i z i = 0, (53) (︂ ∑︁ )︂ n n n ∑︁ ∂K 2 ∑︁ ′ 2 = −a y′i x′i + b y′i + zi − k ∂b i=1 i=1 i=1 n ∑︁ ′ ′ −c x i z i = 0, (54) n n ∑︁ ∑︁ ∂K = −a y′i z′i − b x′i z′i ∂c i=1 i=1 (︂ ∑︁ )︂ n n ∑︁ ′2 ′2 +c yi + x i − k = 0, (55) i=1 i=1 and 3.1 LS line fit in 3D To find the best fitting line to a 3D point cloud the normal distances of all points must be minimized. The normal distance of a point to the investigated line is di = (y′i b − x′i a) + (x′i c − z′i b) + (z′i a − y′i c) . a2 + b2 + c2 (48) A LS solution can be derived by minimizing the sum of the squared distances under the restriction a2 + b2 + c2 = 1. (49) The latter restriction has been chosen as the most appropriate for this case, similarly to the investigated problem of the previous section. Thus, by utilizing Eq. (49) it is possible to derive a simplified expression of the orthogonal distances d i = (y′i b − x′i a) + (x′i c − z′i b) + (z′i a − y′i c). (50) The objective function can be described by Ω(a, b, c) = = n ∑︁ 2 x′i + i=1 +c 2 (︂ ∑︁ n n ∑︁ i=1 2 y′i i=1 −2ac + n ∑︁ i=1 z′i 2 )︂ n ∑︁ i=1 + b2 2 x′i (︂ ∑︁ n i=1 )︂ − 2ab y′i z′i − 2bc n ∑︁ i=1 p1 = n ∑︁ 2 x′i + n ∑︁ i=1 q2 = − ((y′i b − x′i a) + (x′i c − z′i b) + (z′i a − y′i c))2 (︂ ∑︁ n with the respective elements p3 = d2i y′i 2 n ∑︁ i=1 x′i z′i . n ∑︁ z′i i=1 y′i x′i 2 )︂ n ∑︁ i=1 2 y′i n ∑︁ i=1 (51) i=1 subject to a2 + b2 + c2 = 1. Similarly as in the previous investigated problem, we can interpret Eq. (53) to (55) as a homogeneous system of equations, where the solution of the unknown parameter k can be obtained by solving ⃒ ⃒ ⃒ (p − k) ⃒ q1 q2 ⃒ 1 ⃒ ⃒ ⃒ (56) q1 (p2 − k) q3 ⃒ ⃒ = 0, ⃒ ⃒ ⃒ ⃒ q2 q3 (p3 − k) i=1 i=1 i=1 = a2 n ∑︁ i=1 + z′i 2 n ∑︁ , p2 = n ∑︁ 2 y′i + i=1 2 x′i , q1 = − i=1 y′i z′i and q3 = − n ∑︁ n ∑︁ 2 z′i , i=1 y′i x′i , (57) i=1 n ∑︁ x′i z′i . i=1 Equation (56) leads to a cubic equation with the unknown parameter k. As in the previously investigated case (direct LS solution for fitting a straight line in 2D), the adjusted line parameters a, b and c can be estimated by substituting k min into Eq. (53) - (55). 52 | G. Malissiovas et al. 3.2 TLS line fit in 3D A TLS solution for fitting a straight line in 3D using SVD has not been studied yet. Our goal is to derive the TLS solution of fitting a line to a set of points in space for the first time and compare it with the developed direct LS solution from the previous section. In order to build the adjustment model of Eq. (1), it is necessary to derive an appropriate functional model. Rearranging Eq. (47) yields x′ α − y′ β = 0, 0 α − z′ β = x′ , −z′ α + 0 β = y′ , (58) where a b and β = − . (59) c c Thus, the respective matrices and vectors for building the TLS adjustment model are ⎤ ⎡ ⎡ ⎤ e x′1 −e y′1 x′1 −y′1 ⎢ 0 ⎥ ⎢ −e z′1 ⎥ −z′1 ⎥ ⎥ ⎢ ⎢ 0 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ −e z′ ⎢ −z′1 0 ⎥ 0 1 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ .. . . . .. ⎥ , EA = ⎢ .. A=⎢ . .. ⎥ , ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎢ x′ ′ ⎥ ⎢ e x′n −e y′n ⎥ ⎢ n −y n ⎥ ⎥ ⎢ ⎢ ′ ⎥ −z n ⎦ ⎣ 0 ⎣ 0 −e z′n ⎦ −z′n 0 −e z′n 0 ⎡ ⎤ ⎡ ⎤ 0 0 ⎢ e ′ ⎥ ⎢ ′ ⎥ ⎢ x1 ⎥ ⎢ x1 ⎥ ⎢ ⎢ ′ ⎥ ⎥ ⎢ e y′ ⎥ ⎢ y1 ⎥ [︃ ]︃ ⎢ 2 ⎥ ⎥ ⎢ ^ α ⎢ . ⎥ ⎢ .. ⎥ ^= l = ⎢ . ⎥ , el = ⎢ .. ⎥ and x . (60) ⎢ ⎥ ⎢ ⎥ β^ ⎢ ⎢ 0 ⎥ ⎥ ⎢ 0 ⎥ ⎥ ⎢ ⎢ ⎢ ′ ⎥ ⎥ ⎣ e x′n ⎦ ⎣ xn ⎦ y′n e y′n α=− The first column of matrix A contains the coefficients of the functional model (58) with respect to the unknown parameter α, whilst in the second column are the coefficients with respect to the unknown parameter β. The augmented matrix is ⎡ ⎤ x′1 −y′1 0 ⎢ ⎥ −z′1 x′1 ⎥ ⎢ 0 ⎢ ⎥ ⎢ −z′1 0 y′1 ⎥ ⎢ ⎥ ⎢ .. .. ⎥ . [A, l] = ⎢ ... (61) . . ⎥ ⎢ ⎥ ⎢ x′ ⎥ ′ ⎢ n −y n 0 ⎥ ⎢ ⎥ −z′n x′n ⎦ ⎣ 0 ′ ′ −z n 0 yn The right singular vectors of the augmented matrix can be estimated by the eigenvalue/eigenvector decomposition that was described by Eq. (35). Thus, the squared augmented matrix can be expressed by ⎡ ⎤ p1 q1 q2 ⎢ ⎥ [A, l]T [A, l] = ⎣ q1 p2 q3 ⎦ = B , q2 q3 p3 (62) with the respective elements p and q corresponding to those of Eq. (57). The eigenvalues and eigenvectors of matrix B can be found by employing the generalised eigenvalue problem of Eq. (38), which results in ⃒ ⃒ ⃒ ⃒ (p − λ) q1 q2 ⃒ ⃒ 1 ⃒ ⃒ (63) q1 (p2 − λ) q3 ⃒ ⃒ = 0. ⃒ ⃒ ⃒ q2 q3 (p3 − λ) ⃒ The derived determinant provides the characteristic equation of the eigenvalues. In this case this is a cubic equation with three solutions for the unknown eigen^ and β^ can be values λ. The adjusted line parameters α found by employing the minimum eigenvalue principle of Eq. (34), with the right eigenvector corresponding to the smallest eigenvalue of matrix B holding the TLS solution for the adjusted 3D line parameters. Obviously the elements of matrix B coincide with those from the direct LS solution. The determinants of Eq. (56) and (63) are equal, leading to identical characteristic equations. Therefore, the TLS solution for the nonlinear problem of the straight line fit in space is identical with the presented direct LS solution. 3.3 Numerical example - line fit 3D The dataset for the current numerical example has been presented in [31, p.64] and contains the 3D coordinates of points that are subject to fitting a straight line, see Table 5. The results from [31] for the adjusted line parameters are ^x, listed in Table 6, including the orientation components n ^ y and n ^ z , which coincide with the line parameters pren sented by Eq. (47). The results of the direct LS solution are presented in Table 7. The Lagrange multipliers derived from Eq. (56) are ⎧ ⎪ ⎨ 178.9137883045 k^ = (64) 136.5705193088 ⎪ ⎩ 73.3716754581 (k^ ) min Furthermore, the TLS solution for fitting a straight line to the same dataset can be found in Table 8. It should be mentioned that Eq. (34) would have an opposite sign in this case due to the chosen functional model. The eigen- Götterdämmerung over total least squares values derived from the determinant of Eq. (63) are ⎧ ⎪ ⎨ 178.9137883045 ^λ = 136.5705193088 ⎪ ⎩ 73.3716754581 (^λ ) min | 53 Table 8: TLS solution - line fit 3D Estimated line parameters (65) ^ Parameter α Parameter β^ 0.81246017868213 -1.32653589613424 It is directly visible that the direct LS solution for the non-linear case of fitting a straight line in 3D is identical with that of TLS. Both solutions coincide perfectly with the results from [31]. 4 Fitting a plane in 3D Table 5: 3D point cloud from [31, p.64] x-coordinates y-coordinates z-coordinates 55.7290 28.2310 18.8000 12.5270 2.1830 -5.6310 -17.2090 -30.0190 -34.7260 -55.6250 -25.834 -9.0070 -3.2250 0.6280 6.9560 11.7410 18.8330 26.6700 29.5680 42.3630 71.0470 50.3240 43.2170 38.4700 30.6860 24.7830 16.0690 6.4140 2.8520 -12.8930 The third case under investigation is the non-linear problem of fitting a plane to a 3D point cloud with all coordinates being subject to measurement errors. It is assumed also for this case that the coordinates of the points are regarded as observations, which are equally weighted and uncorrelated. Several results from various TLS algorithms were presented for this problem in [12], which resulted in a slight deviation from the LS solution. However, it was already proven that the TLS solution is identical with a direct LS solution for fitting a straight in 2D and 3D. Therefore, a mathematical relation between the TLS and LS is built following the same way as in the previous cases. The general equation of a plane in 3D reads [25, p. 214] Ax + By + Cz + D = 0, Table 6: LS solution for line fit in 3D from [31, p.64] Estimated line orientation components ^x Parameter n ^y Parameter n ^z Parameter n 0.7173305867 -0.4393417007 0.5407547498 Derived line parameters ^ n ^ = − n^ y Parameter α z Parameter β^ = − nn^^ xz 0.812460178 -1.326535896 Table 7: Direct LS solution - line fit 3D Estimated line parameters ^ Parameter a ^ Parameter b Parameter ^c 0.43934170067853 -0.71733058666835 -0.54075474984038 Derived line parameters ^ = − a^^c Parameter α ^ Parameter β^ = − b ^c 0.81246017868213 -1.32653589613424 (66) with x, y and z being the observed 3D coordinates of a point that lies in the plane. A, B , C and D are the plane parameters. A reduction of the number of the unknown plane parameters is, also in this case, possible and can be proven to be equivalent to the procedure presented in the previous sections. Therefore, by reducing the coordinates of the point cloud to the centre of mass, see Eq. (45), results in a simplified model. Parameter D can be eliminated which leads to Ax′ + By′ + Cz′ = 0. (67) 4.1 LS plane fit in 3D Fitting a plane to a set of 3D points, with all coordinates subject to measurement errors, is similar to the case of fitting a straight line in space, as the normal distance of every point to the requested plane must be minimized. The normal distance for this problem can be expressed as [25, p. 216] Ax′ + By′ + Cz′ . (68) d= √ A2 + B2 + C2 54 | G. Malissiovas et al. The LS criterion is applied and a solution is estimated by minimizing the sum of the squared normal distances n ∑︁ i=1 r1 = d2i → min, (69) s1 = A2 + B2 + C2 = 1, (70) since also in this case Eq. (67) can be scaled by an arbitrary factor, which means that only two out of the three parameters A, B or C are independent. Thus, the objective function of this LS problem is n ∑︁ d2i = i=1 n ∑︁ (Ax′i + By′i + Cz′i )2 . n ∑︁ n ∑︁ i=1 subject to the constraint Ω(A, B, C) = with 2 x′i , r2 = s2 = n ∑︁ x′i z′i i=1 2 y′i , r3 = i=1 i=1 y′i x′i , n ∑︁ n ∑︁ 2 z′i , i=1 and s3 = n ∑︁ (77) y′i z′i . i=1 Equation (76) is a cubic equation and has three solutions for the unknown k. The adjusted plane parameters A, B and C can be further estimated following the same procedure of the previously investigated cases, either by substituting parameter k min into Eq. (73) - (75) under the restriction (70) or by transforming the equation system to an eigenvalue problem. The presented direct solution for fitting a plane in 3D coincides with that from [19]. (71) i=1 We search for a LS solution for the unknown parameters A, B and C that minimizes Eq. (71) subject to the constraint (70). These two equations can be utilised to build the Lagrangian K(A, B, C, k) = Ω(A, B, C) − k(A2 + B2 + C2 − 1). (72) The differentiation of K with respect to the unknown parameters A, B and C leads, after setting the partial derivatives to zero, to the system of equations 4.2 TLS plane fit in 3D The TLS solution for fitting a plane in 3D can be estimated analogously to the investigations of [12], following however an alternative functional model from the latter contribution. Based on our approach to the TLS estimate, using the coordinates reduced to the centre of mass, the functional model of Eq. (66) can be rewritten as B A z′ = − x′ − y′ ⇒ z′i = αx′ + βy′ , C C (78) with (︂ ∑︁ )︂ (︂ ∑︁ )︂ n n ∂K ′2 ′ ′ =A xi − k + B yi xi ∂A i=1 i=1 (︂ ∑︁ )︂ n +C x′i z′i = 0, (73) )︂ (︂ ∑︁ (︂ ∑︁ )︂ n n ∂K 2 =A y′i − k y′i x′i + B ∂B i=1 )︂ i=1 (︂ ∑︁ n ′ ′ +C yi zi = 0 (74) )︂ (︂ ∑︁ )︂ (︂ ∑︁ n n ∂K =A x′i z′i + B y′i z′i ∂C (︂i=1 )︂ i=1 n ∑︁ ′2 +C z i − k = 0, (75) B A and β = − . (79) C C The following matrices and vectors can be introduced in order to build the TLS model of Eq. (1) ⎡ ′ ⎤ ⎡ ⎤ e x′1 e y′1 x1 y′1 ⎢ x′ y′ ⎥ ⎢ e ′ e ′ ⎥ ⎢ 2 ⎢ x2 y2 ⎥ 2 ⎥ ⎥ A=⎢ , E = .. ⎥ A ⎢ .. ⎥ ⎢ .. ⎢ .. ⎥, ⎣ . ⎣ . . ⎦ . ⎦ x′n y′n e x′n e y′n ⎤ ⎡ ′ ⎤ ⎡ e z′ z1 ]︃ [︃ ⎢ z′ ⎥ ⎢ e 1′ ⎥ ⎢ 2 ⎥ ⎢ z2 ⎥ ^ α ⎥ ⎢ ⎥ ⎢ ^= . (80) l = ⎢ . ⎥ , el = ⎢ . ⎥ and x β^ ⎣ .. ⎦ ⎣ .. ⎦ z′n e z′n α=− i=1 i=1 and i=1 subject to (70). The solution for the unknown parameter k can be derived by solving ⃒ ⃒ ⃒ (r − k) ⃒ s1 s2 ⃒ 1 ⃒ ⃒ ⃒ (76) s1 (r2 − k) s3 ⃒ ⃒ = 0, ⃒ ⃒ ⃒ s2 s3 (r3 − k) ⃒ The first column of the coefficient matrix A contains the coefficients of the functional model of Eq. (78) with respect to the unknown plane parameter α while in the second column are the coefficients with respect to the unknown parameter β. Further, it is possible to obtain the augmented matrix ⎡ ′ ⎤ x1 y′1 z′1 ⎢ x′ y′ z′ ⎥ ⎢ 2 2 2 ⎥ [A, l] = ⎢ (81) .. .. ⎥ ⎢ .. ⎥ ⎣ . . . ⎦ x′n y′n z′n Götterdämmerung over total least squares and the square matrix [A, l]T [A, l] being equal to ⎤ ⎡ r1 s1 s2 ⎥ ⎢ [A, l]T [A, l] = ⎣ s1 r2 s3 ⎦ = B, s2 s3 r3 (82) with the respective elements from Eq. (77). The eigenvalues and eigenvectors of matrix B can be computed through the generalised eigenvalue problem by solving ⃒ ⃒ ⃒ (r − λ) ⃒ s1 s2 ⃒ 1 ⃒ ⃒ ⃒ (83) s1 (r2 − λ) s3 ⃒ ⃒ = 0, ⃒ ⃒ ⃒ s2 s3 (r3 − λ) ⃒ This characteristic equation is a cubic equation with three solutions for the unknown parameter λ (eigenvalue). The adjusted plane parameters α and β can be estimated once again using the minimum eigenvalue principle. The presented LS solution for fitting a plane in 3D coincides perfectly with the TLS solution. This can be seen from the developed determinants in Eq. (83) and (76), that both lead to the same cubic characteristic equation. 4.3 Numerical example - plane fit 3D To show the equality of the results for the case of the plane fit in 3D, a numerical example is presented. The input data is a synthetic dataset example that was used in [12] and is listed in Table 9. The coordinates are regarded as uncorrelated observations and in all directions (x, y and z) subject to measurement errors. | 55 The solution k^ min can be used to obtain a direct LS estimation for the adjusted plane parameters A, B and C, which are presented in Table 10. Table 10: Direct LS solution for plane fit in 3D Estimated plane parameters LS solution ^ Parameter A ^ Parameter B ^ Parameter C -0.0448859450686 -0.9780188144589 -0.2036282163638 Derived plane parameters ^ ^ = − A^ Parameter α C ^ Parameter β^ = − B ^ C -0.2204308708793 -4.8029631252650 Furthermore, the unknown plane parameters α and β were estimated according to the presented TLS model. Thus, the determinant of Eq. (83) leads to a cubic characteristic equation with three solutions for the unknown eigenvalues ⎧ ⎪ ⎨ 808.2763501846 ^λ = (85) 416.6394328791 ⎪ ⎩ 142.0842169363 (^λ ) min The TLS solution for the adjusted plane parameters can be obtained by employing the minimum eigennvalue principle, using the eigenvector that corresponds to the smallest eigenvalue. The results are listed in Table 11. Table 9: Dataset for fitting a plane in 3D from [12] Table 11: TLS solution for plane fit in 3D point no. x-coord. [m] y-coord. [m] z-coord. [m] 1 2 3 4 5 6 7 8 -10 -9 -11 -10 9 10 11 10 -5.25 -1.25 1.75 6.75 -6.25 -2.25 0.75 5.75 3.75 -9.25 8.75 -5.25 4.75 -8.25 9.75 -4.25 The direct least squares solution for fitting a plane in 3D was derived by employing the presented mathematical approach. Equation (76) results in a cubic equation with three solutions for the unknown parameter ⎧ ⎪ ⎨ 808.2763501846 k^ = (84) 416.6394328791 ⎪ ⎩ 142.0842169363 (k^ ) min Estimated plane parameters TLS solution ^ Parameter α Parameter β^ -0.2204308708793 -4.8029631252650 Various TLS solutions from several algorithms for the requested plane parameters have been presented in [12]. Their eigenvalue solution coincides with our solution and can be found in Table 12. Table 12: Eigenvalue solution from [12] Estimated plane parameters TLS solution ^ (Easterly slope) Parameter α ^ Parameter β (Northerly slope) -0.22043 -4.80296 56 | G. Malissiovas et al. The results coming both from TLS and the developed direct LS solution for the non-linear case of fitting a plane in space are identical, corresponding to equal characteristic equations. It has been demonstrated that both solutions coincide numerically with that of [12]. 5 2D Similarity transformation The last case under investigation is the 2D similarity transformation of coordinates, which is one of the most frequent geodetic and photogrammetric applications. A first attempt to estimate the TLS solution of the problem using SVD was that of Felus and Schaffrin [11] by presenting a Strucured TLS (STLS) algorithm for solving the problem. Neitzel [14] has shown that this algorithm needs to be modified for estimating the correct solution. For this reason, the same problem was examined again and a modified approach has been presented in [15]. The latter solution is iterative, however, it is stated that a TLS solution using SVD could be possible. Here, a new approach is presented for a direct solution of the problem (LS, and also TLS solution via SVD). The well-known equation for the planar coordinate transformation can be given as [︃ ]︃ [︃ ]︃ [︃ ]︃ [︃ ]︃ [︃ ]︃ Xi cos ϕ − sin ϕ µ 0 xi tx = + , Yi sin ϕ cos ϕ 0 µ yi ty (86) see for example [11]. Equation (86) can be rewritten as X i = (µ cos ϕ)x i − (µ sin ϕ)y i + t x Y i = (µ sin ϕ)x i + (µ cos ϕ)y i + t y , ϕ = rotation angle µ = scale factor t x = translation in x-direction t y = translation in y-direction (88) we obtain the simplified equation system X i = ξ1 x i − ξ2 y i + t x Y i = ξ2 x i + ξ1 y i + t y . (89) For a realistic functional model the translation vector has to be present. However, the elimination of the translations (90) with x c and y c denoting the coordinates of the centre of mass of the points in the source system and X c and Y c in the target system. The latter coordinates can be computed by n n 1 ∑︁ 1 ∑︁ xi , yc = yi , n n i=1 i=1 n n 1 ∑︁ 1 ∑︁ Xi , Yc = Yi . Xc = n n xc = i=1 (91) i=1 Therefore, a reduction of all coordinates to the centre of mass leads to a simplified model for the reduced coordinates X ′i = ξ1 x′i − ξ2 y′i Y i′ = ξ2 x′i + ξ1 y′i . (92) 5.1 LS 2D similarity transformation In order to obtain a direct solution in the same manner as in the previous sections, an additional unknown parameter will be taken into consideration. The functional model of Eq. (92) can be rewritten equivalently to the overparametrized system of equations −γ X ′i = α x′i − β y′i , −γ Y i′ = β x′i + α y′i . with ξ1 = − α γ (93) β and ξ2 = − . (94) γ The enforced additional unknown parameter (γ can be seen as an additional parameter) requires to apply a restriction between the unknowns. For the purposes of this research, a meaningful constraint is chosen as α2 + β2 + γ 2 = 1. By substituting ξ1 = µ cos ϕ and ξ2 = µ sin ϕ, t x = X c − ξ1 x c + ξ2 y c , t y = Y c − ξ2 x c − ξ1 y c . (87) with i = 1, ..., n being the number of the observed homologous points between the source xy and the target XY system. The unknown transformation parameters between the two systems are: ś ś ś ś (elimination of two unknown parameters) is possible and can be proven in the same way as for the previous investigated cases by showing that (95) The coordinates of the points in both coordinate systems are subject to measurement errors. By employing the least squares criterion, the goal is to minimize the errors in all homologous points and in both directions. This is equivalent to the minimization of the euclidean distances between the points in the target system and the transformed homologous points from the source system n ∑︁ i=1 D2i → min, (96) Götterdämmerung over total least squares with the distances between two homologous points expressed as D2i = (α x′i − β y′i + γ X ′i )2 +(β x′i + α y′i + γ Y i′ )2 . (97) which leads to a cubic equation with one unknown parameter. The respective elements are v1 = n ∑︁ 2 x′i + n ∑︁ y′i 2 i=1 i=1 w1 = 0, w2 = Therefore, the objective function is , v2 = n ∑︁ = n ∑︁ D2i and w3 = (α x′i − β y′i + γ X ′i )2 i=1 ]︁ +(β x′i + α y′i + γ Y i′ )2 . (98) A LS solution for the unknown transformation parameters α, β and γ is desired, that minimizes the objective function of Eq. (98) under the constraint (95). The Lagrange function can be built as K(α, β, γ , k) = Ω(α, β, γ ) − k(α2 + β2 + γ 2 − 1). (99) A differentiation of the Lagrangian K with respect to the unknown transformation parameters α, β and γ results, when setting the partial derivatives to zero, in the following system of equations (︂ ∑︁ )︂ n n ∑︁ ∂K ′2 ′2 =α xi + yi − k ∂α i=1 i=1 )︂ (︂ ∑︁ n n ∑︁ ′ ′ ′ ′ +γ xi Xi + y i Y i = 0, (100) (︂ ∑︁ )︂ n n ∂K 2 ∑︁ ′ 2 =β x′i + yi − k ∂β i=1 i=1 )︂ (︂ ∑︁ n n ∑︁ ′ ′ +γ xi Yi − y′i X ′i = 0, (101) i=1 i=1 i=1 i=1 and )︂ (︂ ∑︁ n n ∑︁ ∂K ′ ′ ′ ′ =α yi Yi xi Xi + ∂γ i=1 i=1 )︂ (︂ ∑︁ n n ∑︁ +β x′i Y i′ − y′i X ′i +γ (︂ ∑︁ n i=1 i=1 2 x′i + n ∑︁ i=1 (102) i=1 2 y′i − k )︂ = 0, subject to α2 + β2 + γ 2 = 1. Similarly to the previous cases it is possible to estimate k by solving ⃒ ⃒ ⃒ (v − k) w1 w2 ⃒⃒ 1 ⃒ ⃒ ⃒ (103) w1 (v1 − k) w3 ⃒ ⃒ = 0, ⃒ ⃒ ⃒ w2 w3 (v2 − k) ⃒ 2 X ′i + + n ∑︁ n ∑︁ 2 Y i′ , i=1 y′i Y i′ , (104) i=1 x′i Y i′ − n ∑︁ y′i X ′i . i=1 i=1 i=1 n [︁ ∑︁ n ∑︁ n ∑︁ i=1 x′i X ′i i=1 Ω(α, β, γ ) = | 57 The solution for the adjusted transformation parameters α, β and γ can be estimated either by substituting parameter k min to Eq. (100) - (102) or by transforming the determinant of Eq. (103) into an eigenvalue problem. 5.2 TLS 2D similarity transformation In this section a new approach for the TLS solution of the 2D similarity transformation is presented. By utilizing the functional model of Eq. (92) and following the same approach as in the previous cases we can built the TLS adjustment model with the relevant matrices ⎡ x′1 y′1 .. . x′n y′n −y′1 x′1 .. . −y′n x′n ⎤ ⎡ e x′1 −e y′1 ⎢ ⎢ ⎢ ⎢ ⎥ ⎢ e y′ e x′1 ⎢ ⎥ ⎢ 1 ⎢ ⎥ ⎢ . .. ⎢ ⎥ . A=⎢ ⎥ , EA = ⎢ . ⎢ . ⎢ ⎥ ⎢ ′ ⎣ ⎦ ⎢ e x n −e y′n ⎢ ⎣ e y′n e x′n ⎤ ⎡ ⎡ ′ ⎤ X1 e X1 ⎢ e ⎥ ⎢ Y′ ⎥ [︃ ⎢ Y1 ⎥ ⎢ 1 ⎥ ⎢ . ⎥ ⎢ . ⎥ ξ^1 ⎥ ⎢ ⎥ ⎢ ^= l = ⎢ .. ⎥ , el = ⎢ .. ⎥ , x ξ^2 ⎥ ⎢ ⎢ ′ ⎥ ⎣ e Xn ⎦ ⎣ Xn ⎦ e Yn Y n′ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥, ⎥ ⎥ ⎥ ⎥ ⎦ ]︃ . (105) The augmented matrix [A, l] can be described in this case by ⎡ ′ ⎤ x1 −y′1 X1′ ′ ′ ⎥ ⎢ y′ ⎢ 1 x1 Y1 ⎥ ⎢ . .. .. ⎥ ⎥ [A, l] = ⎢ (106) . . ⎥. ⎢ .. ⎢ ′ ⎥ ′ ′ ⎣ x n −y n X n ⎦ y′n x′n Y n′ The essential right eigenvectors of the augmented matrix can then be derived by the eigenvalue/eigenvector decomposition of the squared matrix ⎡ ⎤ v1 w1 w2 ⎢ ⎥ [A, l]T [A, l] = ⎣ w1 v1 w3 ⎦ = B, (107) w2 w3 v2 58 | G. Malissiovas et al. with the respective elements from Eq. (104). The solution for the transformation parameters is determined from the generalised eigenvalue problem by solving ⃒ ⃒ ⃒ (v − λ) w1 w2 ⃒⃒ 1 ⃒ ⃒ ⃒ (108) w1 (v1 − λ) w3 ⃒ = 0. ⃒ ⃒ ⃒ ⃒ ⃒ w2 w3 (v2 − λ) As expected, the derived determinant of Eq. (108) is identical to the one of the direct solution of Eq. (103). Certainly, the resulting cubic equation for the solution of the only unknown term k is identical to the characteristic equation of the eigenvalues λ, too. 5.3 Numerical example - 2D similarity transformation The dataset of the numerical example for the case of the 2D similarity transformation contains the 2D coordinates of homologous points in two coordinate systems and can be found in Table 13. The observations respectively the coordinates of the points are assumed to be equally weighted and uncorrelated. The data originates from [32, pp.397402], while the same example has been utilised in [14]. Table 13: Example dataset for the 2D similarity transformation Point no. i Target S. X i [m] 1 2 3 4 -117.478 117.472 0.015 -0.014 Y i [m] Source S. x i [m] y i [m] 0 0 -117.41 117.451 17.856 252.637 140.089 130.40 144.794 154.448 32.326 267.027 The results of the GHM for the adjusted transformation parameters between the coordinate systems presented in Table 14 originate from the contribution of Neitzel [14]. In this case an iterative solution was employed for deriving the requested transformation parameters ξ1 and ξ2 . The developed direct LS solution for the 2D similarity transformation is presented in Table 15. The adjusted transformation parameters α, β and γ can be used to compute the parameters ξ1 and ξ2 from Eq. (94). Based on the solution of Eq. (108) and the minimum eigenvalue principle, it is possible to obtain the TLS solution for the adjusted transformation parameters of the same dataset. The TLS results are presented in Table 16. Furthermore, it is possible to derive the primal transformation parameters between the two coordinate systems. The rotational angle ϕ as well as the scale factor µ Table 14: Results from [14] Estimated parameters GHM solution Parameter ξ^1 Parameter ξ^2 ^ Scale factor µ ^ Rotational angle ϕ Translation parameter ^t x Translation parameter ^t y 0.99900748077781 -0.04109806319405 0.99985248784424 -2o 21′ 20.72′′ -141.2628 mm -143.9316 mm Table 15: Direct LS solution for the 2D similarity transformation Estimated parameters LS solution Parameter ξ^1 = − γα^^ ^ Parameter ξ^2 = − γβ^ 0.99900748077781 -0.04109806319405 Table 16: TLS solution for the 2D similarity transformation Estimated parameters TLS solution Parameter ξ^1 Parameter ξ^2 0.99900748077781 -0.04109806319405 can be computed by substituting the estimated terms ξ^1 and ξ^2 into Eq. (88). Thus, the rotational angle is (︃ )︃ ξ^2 ^ (109) ϕ = arctan ξ^1 and the scale factor ^= µ ξ^1 ^ cos ϕ . (110) Similarly, the translation parameters t x and t y can be computed by Eq. (90). The solution of the latter transformation parameters coming from the direct LS approach, as well as from TLS, is identical with the presented solution of [14] in Table 14. It has been shown also for this problem that the developed direct LS solution is identical to the TLS solution and both of them correspond numerically with the results obtained from the rigorous GHM. 6 Discussion and Conclusions Within this study a mathematical relationship between TLS and possible direct LS solutions for four individual non-linear problems is developed. The investigated cases are the fit of a straight line in 2D and 3D, the fit of plane Götterdämmerung over total least squares in 3D and the 2D similarity transformation of coordinates, with the measurements being equally weighted and uncorrelated. It has been shown that all cases investigated here can be regarded as a group of non-linear problems, which can be transformed in such a way that a direct solution is possible. The following common features were identified for all the developed direct LS solutions: 1. A non-linear and over-parametrised functional model was used to express each individual problem, see for example Eq. (7) for fitting a straight line to a 2D point set. 2. A reduction of the observed coordinates to the centre of mass was proven in any case to be admissible. Moreover, this reduction leads everytime to the elimination of some unknown parameters. 3. Choosing an appropriate restriction between the unknown parameters for the adjustment of each investigated problem, it was possible to obtain an apparently linear relationship between the observations and the unknowns. 4. The developed objective function for minimising the sum of squared residuals leads always to a homogeneous system of normal equations which is linear with respect to the unknown parameters and has a direct solution. 5. The derived direct solutions have been proven to be identical with the TLS solutions obtained by using SVD. For fitting a straight line in 2D with the coordinates in both directions subject to measurements errors, a direct solution was derived by solving an eigenvalue problem that leads to a characteristic equation (i.e. a quadratic equation), which is similar to the solution presented in [19]. An identical solution (and identical characteristic equation) was obtained also from the TLS adjustment using SVD. On a numerical example it has been shown that the results of the developed direct solution coincide exactly with the TLS solution and the iterative non-linear GHM of [13]. Also for the next two examples it is already known that a direct solution exists. In [16] the case of fitting a straight line to 3D is covered, while in [19] the problem of fitting a plane in 3D was investigated. Here the equivalence of these two direct LS solutions with those of TLS using SVD is clearly demonstrated and a novel TLS solution for fitting | 59 a straight line in a 3D point cloud is presented. Moreover, the methodology applied is the same as in the first case. The last discussed case in this contribution is the 2D similarity transformation of coordinates. The TLS solution of the problem using SVD has been presented in the first place by Felus and Schaffrin [11]. However, Neitzel [14] obtained the LS solution of the same problem by means of an iterative linearized GHM and showed that modifications were necessary in the TLS algorithm of [11]. The corresponding modified TLS solution of the same problem was presented in [15]. The latter TLS solution has been obtained iteratively and a solution via SVD is not clearly defined there, however, it is stated that such a solution is possible. In this contribution a direct solution is presented for the 2D similarity transformation of coordinates. The coordinates of the points in both the target and the source system are regarded as uncorrelated observations, subject to measurements errors. Additionally, a novel approach for a TLS solution with SVD is derived. It has been shown that both the direct and TLS solution are identical and they coincide numerically with the solution from [14]. In contrast to TLS, the developed systematic approach for solving directly non linear LS problems provides a clear description of the observed quantities and the unknown parameters of the problem. An optimal solution can be estimated by minimizing a clearly defined objective function which is based on the method of LS. The investigated class of problems leads always to a system of normal equations that can be solved with various techniques, for example SVD or EVD. Therefore, the clarity of the developed method can enable scientists to deeply understand the concept of TLS and set it under scrutiny when it comes to the solution of non-linear LS problems. Now, it is clear that a direct solution exists for a certain class of non-linear LS problems that includes the four investigated cases. The derived direct LS solution was obtained in all cases by solving a characteristic equation, which was always identical to the characteristic equation of the eigenvalues from the corresponding TLS solution. The choice between the two options for a solution is always on the user’s preference. Therefore, this contribution can serve as a systematic approach for solving directly nonlinear LS problems, even if the user has never been in contact with the "TLS-world". Acknowledgement: This research was supported by the Deutsche Forschungsgemeinschaft (DFG), grant NE 1387/2-1. 60 | G. Malissiovas et al. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] C. Gauss, Theoria motus corporum coelestium in sectionibus conicis solem anbientium, F. Perthes and I. H. Besser, Hamburg, 1809. W. Niemeier, Ausgleichungsrechnung, 2nd Edition, Walter de Gruyter, New York, 2008. E. J. Krakiwsky, A synthesis of recent advances in the method of least squares, Department of Geodesy and Geomatics Engineering, University of New Brunswick, 1975. P. J. Bickel, Y. Ritov, Eflcient estimation in the errors in variables model, Ann. Stat. 15 (2) (1987) 513ś540. S. Van Huffel, J. Vandewalle, The Total Least Squares Problem, computational aspects and analysis, SIAM, 1991. D. Wells, E. J. Krakiwsky, The method of least squares, Department of Geodesy and Geomatics Engineering, University of New Brunswick, 1971. E. M. Mikhail, F. Ackermann, Observations and Least Squares, Thomas Y. Crowell Company, Inc., 1976. P. Groen, An introduction to total least squares, Nieuw Archief voor Wiskunde 14 (1996) 237ś253. G. Golub, C. Van Loan, An analysis of the total least squares problem, SIAM 17 (1980) 883ś893. S. Van Huffel, J. Vandewalle, Algebraic connections between the least squares and total least squares problems, Numer. Math. 55 (1989) 431ś449. Y. Felus, B. Schaffrin, Performing similarity transformations using the errors-in-variables-model, Proceedings of the ASPRS Meeting, Washington DC 35 (2005) 751ś762. B. Schaffrin, I. Lee, Y. Felus, Y. Choi, Total least-squares(TLS) for geodetic straight-line and plane adjustment, Boll Geod. Sci. Aflni 65 3 (2006) 141ś168. F. Neitzel, S. Petrovic, Total Least Squares (TLS) im Kontext der Ausgleichung nach kleinsten Quadraten am Beispiel der ausgleichenden Geraden, zfv 133 (2008) 141ś148. F. Neitzel, Generalisation of total least squares on example of unweighted and weighted similarity transformation, J. Geodesy 84 (12) (2010) 751ś762. B. Schaffrin, F. Neitzel, S. Uzun, V. Mahboub, Modifying cadzow’s algorithm to generate the optimal TLS-solution for the structured EIV-model of a similarity transformation, J Geod. Sci. 2 (2012) 98ś106. D. Jovičić, M. Lapaine, S. Petrović, Fitting a straight line to a set of points in space (in croatian), Geodetski list 36 (59) (1982) 260ś266. T. Krarup, J. Juhl, K. Kubik, Götterdämmerung over Least Squares Adjustment, 14th congress of the international society of photogrammetry, Hamburg B3 (1980) 369 ś 378. K. R. Koch, Parameter estimation and hypothesis testing in linear models, second, updated and enlarged Edition, SpringerVerlag Berlin Heidelberg, 1999. K. Linkwitz, Über einige Ausgleichungsprobleme und ihre Lösung mit Hilfe Matrizen-Eigenwerten, Deutsche Geodätische Kommission, Wissenschaftliche Beiträge aus dem Kreise der Schüler von Ernst Gotthard anlässlich seiner Emeritierung, 1976. R. Adcock, A problem in least squares, The Analyst 5 (1878) 53ś 54. [21] K. Pearson, On Lines and Planes of Closest Fit to Systems of Points in Space, F.R.S, University College, London. [22] S. Van Huffel, Total least squares and errors-in-variables modeling: Bridging the gap between statistics, computational mathematics and engineering, Compstat 17 (2004) 883ś893. [23] I. Markovsky, S. Van Huffel, Overview of total least-squares methods, Signal Process. 87 (2007) 2283ś2302. [24] B. Schaffrin, Connecting the dots: The straight-line case revisited, zfv 132 (2008) 385ś394. [25] I. Bronshtein, K. Semendyayev, G. Musiol, H. Muehlig, Handbook of Mathematics, 5th Edition, Springer, Berlin Heidelberg New York, 2005. [26] C. Lawson, R. Hanson, Solving Least Squares Problems, SIAM,Philadelphia, 1974. [27] G. Golub, C. Van Loan, Matrix Computation, 2nd Edition, The Johns Hopkins University Press, Baltimore, Maryland, 1989. [28] G. Kampmann, B. Renner, Vergleich verschiedener Methoden zur Bestimmung ausgleichender Ebenen und Geraden, AVN 2 (2004) 56ś67. [29] S. Kupferer, Verschiedene Ansätze zur Schätzung einer ausgleichenden Raumgeraden, AVN 5 (2004) 162ś170. [30] H. Späth, Zur numerischen Berechnung der Trägheitsgeraden und der Trägheitsebene, AVN 7 (2004) 273ś275. [31] E. Drixler, Analyse der Form und Lage von Objekten im Raum, Vol. C, 409, Deutsche Geodätische Kommission bei der Bayerischen Akademie der Wissenschaften, 1993. [32] E. Mikhail, J. Bethel, C. McGlone, Introduction to Modern Photogrammetry, Wiley, New York Chichester, 2001.