In this paper we address the problem of finding a geodesic curve that best fits a given set of time-labeled points  on a sphere.
Since the corresponding normal equations are highly non-linear, we formulate the problem as a constrained nonlinear optimization problem and solve it using the routine fmincon from MATLAB with the SQP (Sequential Quadratic Programming) algorithm.Foundation for Science and Technology in Projects scope: FCOMP-01-0124-FEDER-022674 and PTDC/EEACRO/
122812/201

Machado, L.

Monteiro, M. Teresa T.

Universidade do Minho: RepositoriUM

Proceedings of the 13th International Conferenceon Computational and Mathematical Methodsin Science and Engineering, CMMSE 201324{27 June, 2013.Geodesic regression on spheres:a numerical optimization approachL. Machado1 and T. Monteiro21 Department of Mathematics & ISR, University of Tras-os-Montes and Alto Douro &University of Coimbra2 R&D Centre Algoritmi, Department of Production and Systems, University of Minho,Portugalemails: lmiguel@utad.pt, tm@dps.uminho.ptAbstractIn this paper we address the problem of nding a geodesic curve that best ts agiven set of time-labeled points on a sphere. Since the corresponding normal equationsare highly non-linear, we formulate the problem as a constrained nonlinear optimiza-tion problem and solve it using the routine fmincon from MATLAB with the SQP(Sequential Quadratic Programming) algorithm.Key words: Manifolds, geodesics, geodesic distance, normal equations, constrainednonlinear optimization, SQP.1 IntroductionThe astounding development of mechanical and robotics industry in the past few yearshas required the generalization of classical methods to more general curved spaces. Thisis mainly due to the fact that the conguration systems of the most part of mechanicalsystems are particular manifolds, like Lie groups or symmetric spaces. Although obtainingsuch generalizations is not as straightforward as we might expect.In this paper we show how to generalize the linear regression problem on Euclideanspaces to the n dimensional unit sphere. This technique of approximating data is a commonprocedure in several applications from a wide range of elds including statistics, computervision, signal processing, fuzzi control, air trac and aeronautics control.In the classical problem, [6], we are given a collection of points and the same numberof instants of time and the objective is to nd a parameterized straight line that best tscCMMSE ISBN: 978-84-616-2723-3Geodesic regression on spheresthe data. However, nding the best approximant curve is not a trivial task especially whenthe points belong to some curved space. We refer to [5], where it has been developed avariational approach to generate tting curves on the sphere and to [8] and [9], where thisapproach has been generalized to the more general context of Riemannian manifolds. Themajor diculty that we face is the fact that, in general, explicit formulas to the analoguesof straight lines, the so-called geodesics, are not available. Such is not the case of the unitn sphere, when equipped with the metric induced by the metric in the embedding space,where geodesics are the great arc circles, [4].To get some insight, we start, in Section 2, by recalling the classical linear regressionproblem on Euclidean spaces. In Section 3, we formulate the corresponding problem on then dimensional unit sphere and then present the rst order necessary optimality conditions.Unlike the Euclidean case, this system of equations is extremely non-linear and numericaloptimization methods are put to use in Section 4. In this section, the problem is formu-lated as a constrained nonlinear optimization problem and numerical experiments using theoptimization toolbox from MATLAB will be provided. Some conclusions and future workideas are carried out in Section 5.2 Linear regression problemAlthough in the most part of the literature, [6, 11, 3], the linear least squares problem isaddressed for data in R, the approach to more general Euclidean spaces is straightforward,[7].Let us consider the Euclidean space Rn endowed with the usual inner product and letus denote by d the metric induced by the l2-norm kak = ha; ai 12 .In the linear regression problem, we are given a nite set of points in Rn, p0; : : : ; pN ,and a set of instants of time, t0; : : : ; tN , and wish to nd a parameterized straight linet 7 ! (t) = a0 + a1t 2 Rn, that best ts the given data, in the sense that the functionalE, dened byE() =NXi=0d2(pi; (ti));should be as small as possible.Theorem 2.1. For each N  1, the parameterized straight line t 7! (t) = a0 + a1t thatbest ts the given data is unique and is the solution of the following system of equations8>>>><>>>>:NXi=0(ti) =NXi=0piNXi=0ti(ti) =NXi=0tipi: (1)cCMMSE ISBN: 978-84-616-2723-3L. Machado, T. MonteiroThe linear system of equations (1), known in the literature as the normal equations,can be solved explicitly and its unique solution is given by(t) =NXi=0t2iNXi=0pi  NXi=0tiNXi=0tipi(N+1)NXi=0t2i    NXi=0ti2 +(N+1)NXi=0tipi  NXi=0tiNXi=0pi(N+1)NXi=0t2i    NXi=0ti2 t:This classical problem can be naturally generalized to more general curved spaces aslong as explicit formulas for geodesics are available, [8]. Nevertheless, even in those casesobtaining exact solutions is not an easy task mainly because the counterpart of the normalequations give rise to nonlinear systems of equations.In the next section we will present the generalization of this classical problem to then dimensional unit sphere Sn.3 Regression problem on spheresLet us consider the unit n sphere Sn as an embedded submanifold of the Euclidean spaceRn+1. Since the tangent space of Sn at a point p 2 Sn isTpSn =v 2 Rn+1 : hv; pi = 0	;let us dene an inner product in TpSn byhu; vi = u>v; u; v 2 TpSn:With this inner product, Sn can be considered a Riemannian manifold whose metric isthe one induced by the Euclidean inner product in Rn+1.Geodesics with respect to this metric are the solutions of the second order dierentialequation   h; i = 0:The unique geodesic t 7! (t) with initial conditions (0) = p 2 Sn and _(0) = v 2 TpSnis given by(t) = p cos(tkvk) + vkvk sin(tkvk): (2)If we settled on the 2 dimensional sphere, S2, geodesics are simply the great arc circles.In order to nd out the geodesic distance with respect to the above metric, one justneeds to compute the velocity vector of the geodesic that joins two points on the sphere.Let us assume that we are given two points p and q on Sn such that the angle betweencCMMSE ISBN: 978-84-616-2723-3Geodesic regression on spheresthem lies on the interval ]0; [. Then, the geodesic that joins p, at t = 0, to q, at t = 1, canbe parameterized explicitly as(t) = p cos(t) +q   p cossinsin(t);where t 2 [0; 1] and  = arccoshp; qi, [2]. Therefore, the geodesic distance between p and qis given by the length of the velocity vector of  at t = 0, that is,d(p; q) = arccoshp; qi: (3)3.1 Problem's formulationSince we have already dened the distance function, we are now in conditions to formulatethe analogous to the linear regression problem on Sn. Let us consider a collection of N +1points, p0; : : : ; pN , on Sn and a monotone increasing sequence of instants of time t0 <t1    < tN , that we assume, for simplicity, that they form a partition of the unit timeinterval [0; 1]. Our main goal is to nd a geodesic on Sn parameterized explicitly by(t) = p cos(tkvk) + vkvk sin(tkvk); (4)where p 2 Sn and v 2 TpSn, that best ts the given data in the sense that it yields theminimum value for the functionalE() =NXi=0d2(pi; (ti));where d is the geodesic distance on Sn dened by (3).Notice that nding  is equivalent to nd p 2 Sn and v 2 TpSn that minimize thefunctionF (p; v) =NXi=0arccos2pi; (ti): (5)Theorem 3.1. ([8]) A necessary condition for t 7! (t) = p cos (kvkt) + vkvk sin (kvkt) to bethe geodesic that best ts the given data (points and instants of time) is that the pair (p; v)satises the following system of equations:8>>>><>>>>:NXi=0isinicos kvkti pi   hpi; pip = 0NXi=0i sin(kvkti)sini pi   hpi; pi(tiv + p)  hpi;vikvk2 v=NXi=0 i cos(kvkti)sinihpi;vikvk tiv; (6)where i = arccospi; p cos kvkti+ vkvk sin kvkti, for i = 0; : : : ; N .cCMMSE ISBN: 978-84-616-2723-3L. Machado, T. MonteiroAnalogously to what happens in the Euclidean case, we call to the system of equations(6) the normal equations for the geodesic regression problem on Sn.Since the above system of equations is highly nonlinear, a numerical optimization ap-proach to nd out the approximate solutions of the proposed optimization problem will beconsidered in the next section.4 Numerical tests using an optimization approachLet us consider the constrained optimization problem:minp;vF (p; v) =NXi=0arccos2hpi; (ti)is.t. kpk = 1hp; vi = 0(7)whose objective function is dened in (5). The equality constraints mean that the point pmust belong to the unit n-sphere Sn and v must be orthogonal to p.4.1 Computational experimentsThe computational experiments were made on a 2.0 GHz Intel Core i7 with 8GB of RAM,Windows 7 64-bit operating system. The MATLAB version used was 7.13.0.564 (R2011b).The problem was codied in MATLAB and solved with the fmincon routine from theoptimization toolbox. This routine, fmincon, attempts to nd a constrained minimum of ascalar function of several variables starting at an initial estimate. This is generally referredto as constrained nonlinear optimization or nonlinear programming. fmincon has fouralgorithm options: interior-point, SQP, active-set and trust-region-reective (default). Inthese experiments the SQP algorithm is used, [10]. SQP is an iterative method for nonlinearoptimization used on problems for which the objective function and the constraints are twicecontinuously dierentiable. SQP methods solve a sequence of optimization subproblems,each of which optimizes a quadratic model of the objective objective function subject to alinearization of the constraints.Several numerical tests were performed on the two dimensional unit sphere S2. Table1 reports information from two of them. In the rst column it is indicated the number ofpoints (N + 1) and in the second column the set of instants of time (t). The third columnreports the N + 1 points randomly generated on S2. p, v and F  are the optimal valuesof p, v and F , respectively, and iter denotes the number of iterations carried out by thefmincon routine.cCMMSE ISBN: 978-84-616-2723-3Geodesic regression on spheresFigure 1 shows the geodesic that best ts two (1a) and ve given data (1b), respectively,corresponding to the tests mentioned on Table 1. Notice that, as expected, the geodesicthat best ts two points is exactly the one that joins them (interpolation case).Table 1: Numerical testsN + 1 t Data p v F iter2 f0,1g0.1585 -0.5605 0.1585 -0.9351 10 14 14-0.8624 0.4575 -0.8624 0.56960.4807 0.6903 0.4807 1.33025f0, 0.25, -0.89 0.8243 -0.9164 -0.9383 0.0192 -0.8006 4.63663.4203 470.5, 0.75, 1g -0.3941 -0.2085 -0.158 0.3022 -0.7635 -0.2176 -32,542-0.2292 0.5264 -0.3678 0.1678 0.6455 0.5583 5.3801(a) Two data points (b) Five data pointsFigure 1: Geodesics on S25 Conclusions and future workIn this paper we formulated the analogous linear regression problem to the unit n sphere.In contrast to what happens in the Euclidean case this optimization problem cannot besolved analytically due to nonlinearity of the counterpart of the normal equations (6).To overcome this diculty we have successfully used the MATLAB optimization toolboxroutine fmincon.Since in the Euclidean case higher order degree polynomials can be used to t a givendata set of points, as future work we aim to nd more general tting curves on spheres.cCMMSE ISBN: 978-84-616-2723-3L. Machado, T. MonteiroAcknowledgementsThis work is funded by FEDER funds through Operational Programme for Competitive-ness Factors - COMPETE and National Funds through FCT - Foundation for Scienceand Technology in Projects scope: FCOMP-01-0124-FEDER-022674 and PTDC/EEA-CRO/122812/2010.References[1] S. R. Buss and J. P. Fillmore, Spherical Averages and Applications to SphericalSplines and Interpolation, ACM Transactions on Graphics, 2, no. 20, (2001), 95{126.[2] P. Crouch, G. Kun and F. Silva Leite, The De Casteljau Algorithm on Lie Groupsand Spheres, J. of Dyn. Control Sys., 5, no. 3, (1999) 397{429.[3] G. Farin, Curves and Surfaces for Computer Aided Geometric Design, AcademicPress, 1990.[4] J. Jost Riemannian geometry and geometric analysis (6 Ed.), Universitext, Springer,2011.[5] P. E. Jupp and J. T. Kent, Fitting Smooth Paths to Spherical Data, Appl. Statist.,36, no. 1, (1987), 34{46.[6] P. Lancaster and K. Salkauskas, Curve and Surface Fitting, Academic Press,1990.[7] L. Machado, Least squares Problem on Riemannian manifolds, PhD Thesis, Univer-sity of Coimbra, 2006.[8] L. Machado and F. Silva Leite, Fitting Smooth Paths on Riemannian Manifolds,Int. J. Appl. Math. Stat., 4, no. J06, (2006), 25{53.[9] L. Machado, F. Silva Leite and K. Krakowski, Higher-order smoothing splinesversus least squares problems on Riemannian manifolds, J. of Dyn. Control Sys., 16,no. 1, (2010), 121{148.[10] J. Nocedal and S. J. Wright, Numerical Optimization, Second Edition, SpringerSeries in Operations Research, Springer Verlag, 2006.[11] M. J. D. Powell, Approximation Theory and Methods, Cambridge University Press,1981.cCMMSE ISBN: 978-84-616-2723-3

Geodesic regression on spheres : a numerical optimization approach

Abstract

Similar works

Full text

Available Versions

Universidade do Minho: RepositoriUM

Universidade do Minho: RepositoriUM