A variational model for data fitting on manifolds by minimizing the acceleration of a Bézier curve
Ronny Bergmann Pierre-Yves Gousenbourger



Learning to Predict Crisp Boundaries
    Ruoxi Deng^1,⋆, Chunhua Shen^2, Shengjun Liu^1, Huibing Wang^3,⋆, Xinru Liu^1
    Received 22 May 2018 / Accepted 17 July 2018
=================================================================================


[1]Research Group Numerical Mathematics
  (Partial Differential Equations),
    Faculty of Mathematics,

    Technische Universität Chemnitz,
    D-09107 Chemnitz, Germany 

    E-mail: mailto:ronny.bergmann@mathematik.tu-chemnitz.de
[2]ICTEAM Institute,

  Université catholique de Louvain,
  B-1348 Louvain-la-Neuve, Belgium 

    E-mail: mailto:pierre-yves.gousenbourger@uclouvain.be





  We derive a variational model to fit a composite Bézier curve to a set of
  data points on a Riemannian manifold. The resulting curve 
  is obtained in such a way that its mean squared acceleration is minimal 
  in addition to remaining close the data points.
  We approximate the acceleration 
  by discretizing the squared second order derivative along the curve. 
  We derive a closed-form, numerically stable and efficient
  algorithm to compute the gradient of a Bézier curve on manifolds
  with respect to its control points, expressed as a concatenation of so-called
  adjoint Jacobi fields.
  Several examples illustrate the capabilites and validity of this
  approach both for interpolation and approximation. The examples also
  illustrate that the approach outperforms previous works tackling this problem.



  
AMS subject classification (2010). 65K10, 65D10, 65D25, 53C22, 49Q99


  
Key words. Riemannian manifolds, curve fitting, composite Bézier curves, Jacobi fields, variational models.




§ INTRODUCTION

  
  This papers addresses the problem of fitting a smooth curve to
  data points d_0,…,d_n lying on a Riemannian manifold  and
  associated with real-valued parameters t_0,…,t_n. The curve
  strikes a balance between a data proximity constraint and a smoothing
  regularization constraint.

  
  Several applications motivate this problem in engineering and the sciences.
  For instance, curve fitting is of high interest in
  projection-based model order reduction of one-dimensional dynamical
  systems <cit.>. In that application, the dynamical
  system depends on the Reynolds number and the model reduction is
  obtained by computing suitable projectors as points on a Grassmann
  manifold. Finding a projector is however a time- and memory-consuming
  task. Based on projectors precomputed for given parameter values,
  fitting is used to approximate the projector associated with a new
  parameter value.  In <cit.>, the same strategy is
  used to approximate wind field orientations represented as points on
  the manifold of positive semidefinite covariance matrices of size p and
  rank r. Further applications are involving rigid-body motions on
  SE(3), like Cosserat rods applications <cit.>, 
  or orientation tracking <cit.>. Sphere-valued data are also of
  interest in many applications of data analysis, for storm tracking and
  prediction, or the study of bird migration <cit.>.

  
  There exists different approaches to tackle the curve fitting problem.
  Among others, we name here the subdivision schemes
  approach <cit.> or the Lie-algebraic
  methods <cit.>.
  However, the most popular approach nowadays is probably to encapsulate
  the two above-mentioned constraints into an optimization problem
  
    min_γ∈Γ E_λ(γ) ∫_t_0^t_n‖D^2 γ(t)/ṭ^2‖^2_γ(t)ṭ
        +
        λ/2∑_i=0^n ^̣2(γ(t_i),d_i)),

  where Γ is an admissible space of curves γ
  [t_0,t_n] →, t ↦γ(t),
  D^2/ṭ^2 denotes the (Levi-Civita) second
  covariant derivative,
  ‖·‖_γ(t) denotes the Riemannian metric at γ(t)
  from which we can define a Riemannian distance (̣·,·).
  Finally, λ∈ is a parameter that strikes the balance between the
  regularizer
  ∫_t_0^t_n‖D^2
  γ(t)/ṭ^2‖^2_γ(t)ṭ and the
  fitting term ∑_i=0^n ^̣2(γ(t_i),d_i)).
  This approach leads to the remarkable property that, when the
  manifold  reduces to the Euclidean space and Γ is the
  Sobolev space H^2(t_0,t_n), the solution to (<ref>) is the natural
  cubic spline <cit.>.
  
  Optimization on manifolds has gained a lot of interest this last decade, 
  starting with the textbook <cit.> that summarizes several optimization
  methods on matrix manifolds. Recently, toolboxes have emerged, providing 
  easy access to such optimization methods, , Manopt <cit.> and
  MVIRT <cit.>.
  The former received a very positive return in many different topics of
  research, with applications for example in low-rank
  modelling in image analysis <cit.>,
  dimensionality reduction <cit.>, phase retrieval <cit.>
  or even 5G-like MIMO systems <cit.>.
  The latter stems from recent interest in manifold-valued image and data
  processing, phrased as variational models on the product manifold ^N,
  where N is the number of pixels. Starting with total variation (TV)
  regularization of phase-valued data <cit.>,
  different methods for TV
  on manifolds <cit.> have been developed as well as
  second order methods <cit.> up to infimal
  convolution <cit.> and total generalized
  variation <cit.>.
  Furthermore, different algorithms have been generalized to manifold-valued
  data, besides the previous works using gradient descent or cyclic proximal point
  methods, a Douglas–Rachford splitting <cit.>, iteratively reweighted
  least squares <cit.> and more general half-quadratic
  minimization <cit.> have been introduced.

  The curve fitting problem (<ref>), has been tackled differently the past few
  years. Samir  <cit.> considered the case where Γ is an
  infinite dimensional Sobolev space of curves, and used the Palais-metric to
  design a gradient descent algorithm for (<ref>) (see, , <cit.> 
  for an application of this approach).
  Another method consists in discretizing the curve γ in N points,
  and therefore consider Γ = ^N (see, , <cit.> for
  a result on SO(3)).
  Finally, the limit case where λ→ 0 is already well studied and
  known as the geodesic regression <cit.>.

  
  A recent topic concerns curve fitting by means of Bézier curves.
  In that approach, the search space Γ is reduced to as set of
  composite Bézier curves. Those are a very versatile tool to model
  smooth curves and surfaces for real- and vector-valued discrete data
  points (see <cit.> for a comprehensive textbook), but they
  can also be used to model smooth curves and surfaces for manifold-valued
  data <cit.>. The advantage to work with such objects,
  compared to classical approaches, are that (i) the search space is
  drastically reduced to the so-called control points of the Bézier curves
  (and this leads to better time and memory performances) and (ii) it is
  very simple to impose differentiability for the optimal curve, which is
  appreciated in several of the above-mentioned applications.
  However, while obtaining such an optimal curve reduces directly to
  solving a linear system of equations for data given on a Euclidean space,
  there is up to now no known closed form of the optimal Bézier curve for
  manifold valued data.

  
  In this work, we derive a gradient descent algorithm to compute a
  differentiable composite Bézier curve [t_0,t_n]→ that 
  satisfies (<ref>), , such that (t) has a minimal 
  mean squared acceleration, and fits the set of n+1 manifold-valued 
  data points at their associated time-parameters. 
  We consider the manifold-valued generalization of Bézier curves <cit.>
  in the same setting as in Arnould  <cit.>, 
  or more recently in <cit.>.
  We employ the (squared) second order absolute differences introduced
  in <cit.> to obtain a discrete approximation of the regularizer
  from (<ref>). The quality of the approximation depends only on the number
  of sampling points. We exploit the recursive structure of the De Casteljau 
  algorithm <cit.> to derive the gradient of the objective function 
  with respect to the control points of (t). The gradient is built
  as a recursion of Jacobi fields that, for numerical reasons, are implemented
  as a concatenation of so-called adjoint Jacobi fields. Furthermore,
  the corresponding variational model only depends on the number of control 
  points of the composite Bézier curve, and not on the number of sampling points. 
  We finally obtain an approximating model to (<ref>) that we
  solve with an algorithm only based on three tools on the manifold:
  the exponential map, the logarithmic map, and a certain Jacobi field along geodesics.
  
  
  The paper is organized as follows. We introduce the necessary
  preliminaries —Bézier curves, Riemannian manifolds and Riemannian second
  order finite differences— in Section <ref>.
  In Section <ref> we derive the 
  gradient of the discretized mean squared acceleration of the composite
  Bézier curve with respect to its control points,
  and thus of the regularizer of (<ref>).
  In Section <ref>, we present the corresponding gradient descent 
  algorithm, as well as an efficient gradient evaluation method, to
  solve (<ref>) for different values of λ. The limit case
  where λ→∞ is studied as well.
  Finally, in Section <ref>, we validate, analyze and illustrate
  the performance of the algorithm for several numerical examples on the
  sphere 𝕊^2 and on the special orthogonal group SO(3).
  We also compare our solution to existing Bézier fitting methods.
  A conclusion is given in Section <ref>.





§ PRELIMINARIES

  
  
  

 §.§ Bézier functions and composite Bézier spline

  
    Consider the Euclidean space ^m. A Bézier curve of degree
    K∈ is a function β_K[0,1]→^m parametrized by
    control points b_0,…,b_K ∈^m and taking the form
    
    β_K(t ; b_0, …, b_K)
          ∑_j=0^K b_j B_j,K(t),

    where B_j,K(t) = Kjt^j (1-t)^K-j are the Bernstein basis
    polynomials <cit.>.
    For example the linear Bézier curve β_1 is just the line segment
    (1-t)b_0 + tb_1 connecting the two control points b_0 and b_1. 
    The explicit formulae of the quadratic and
    cubic Bézier curves read
    
    β_2(t ; b_0, b_1, b_2)    = b_0(1-t)^2 + 2b_1(1-t)t + b_2 t^2,
    β_3(t ; b_0, b_1, b_2, b_3)
               = b_0(1-t)^3 + 3 b_1(1-t)^2t + 3 b_2(1-t)t^2 + b_3 t^3,

    for given control points b_0,b_1,b_2∈ℝ^m and
    an additional point b_3∈ℝ^m for the cubic case.

    A composite Bézier curve is a function
    [0,n]→^m composed of n Bézier curves and defined
    as
    
    (t) β_K_0(t; b_0^0,…,b_K_0^0)    t∈ [0,1]
    β_K_i(t-i; b_0^i,…,b_K_i^i)   
                t∈ (i,i+1],    i∈{1,…,n-1},

    where K_i∈ denotes the degree of the i^th 
    Bézier curve β_K_i of  and b_j^i, j=0,…,K_i,
    are its control points. Furthermore, (t) is continuous if the last
    and first control points of two consecutive segments
    coincide <cit.>. We introduce p_i as the point at
    the junction of two consecutive Bézier segments, , p_i
    b_K_i-1^i-1 = b_0^i. Differentiability is obtained if the control
    points of two consecutive segments are aligned.
    We introduce further b_i+1^-  b_K_i-1^i and b_i+1^+  b_1^i+1 such that the differentiability condition
    reads p_i+1 = K_i b_i+1^- + K_i+1 b^i+1_+/K_i +
    K_i+1.

    [C^1 conditions for composite cubic Bézier curves]
      We consider the case where the Bézier segments are all cubic, as
      represented in Fig. <ref>. The composite Bézier
      curve (t) is C^1 if p_i = 1/2(b_i^- +
      b_i^+), i = 1,…,4, with p_i = b_3^i-1 = b_0^i, b_i^- = b_2^i-1, and b_i^+ = b_1^i.
    

    
  
  
  

 §.§ Riemannian Manifolds

    We consider a complete m-dimensional Riemannian manifold . We refer
    to <cit.> for an introduction to Riemannian manifolds
    and to <cit.> for optimization thereon. We will use the following
    terminology and notations.

    We denote by T_a the (Euclidean) tangent space to  at a ∈;
    T∪_a∈T_a is the tangent bundle to ; ⟨·, ·⟩_a denotes the inner product in the tangent space T_a at
    a and from which we deduce the norm of v ∈ T_a denoted by ‖ v
    ‖_a = √(⟨ v, v ⟩_a). For a (not necessarily unique)
    shortest geodesic between a and b∈, we write g(·;a,b)→,
    t↦ g(t;a,b), parametrized such that g(0;a,b) = a
    and g(1;a,b) = b. This choice of parametrization also means
    that the covariant derivative D/dt of g with respect to time
    satisfies ‖Ddtg(t;a,b)‖_g(t) = d_ℳ(a,b),
    for all t∈[0,1], where d_(a,b) is the geodesic distance between a and b ∈.
    The Riemannian exponential reads a T_a→, v ↦ b = a[v]
    and we denote by r_a∈ℝ the maximal radius such
    that the exponential map is bijective on 𝒟_a { b∈ :
    d_ℳ(a,b) < r_a}. Then a𝒟_a
    → T_a, b ↦ v = a[b] is called the Riemannian logarithm
    which is (locally) the inverse of the exponential. A Riemannian manifold is
    called symmetric in x∈ℳ if the geodesic reflection s_x at
    x∈ℳ given by the mapping γ(t)↦γ(-t)
    is an isometry at least locally near x,
    for all geodesics through γ(0)=x. If ℳ is symmetric
    in every x∈ℳ, the
    manifold is called (Riemannian) symmetric space or symmetric manifold.

    In the following we assume, that both the exponential and the logarithmic
    map are available for the manifold and that they are computationally not too
    expensive to evaluate. Furthermore, we assume that the manifold is
    symmetric.
  
  
  

 §.§ Composite Bézier curves on manifolds

    One well-known way to generalize Bézier curves to a Riemannian manifold
    ℳ is via the De Casteljau algorithm <cit.>.
    This algorithm only
    requires the Riemannian exponential and logarithm and conserves the
    interpolation property of the first and last control points. Some examples
    on interpolation and fitting with Bézier curves or Bézier surfaces,
    , generalizations of tensor product Bézier curves,
    can be found in Absil  <cit.>.

    Consider β_K[0,1] →,
    t ↦β_K(t;b_0,b_1,…,b_K),
    the manifold-valued Bézier curve of order
    K driven by K+1 control points b_0,…,b_K∈.
    We introduce the points
    x_i^0 = b_i and iterate the construction of further points. 
    For i=0,…,K-k, k=1,…,K, we define
    
    x_i^kβ_k(t;b_i,…,b_i+k)
          = g(t;x_i^k-1,x_i+1^k-1)

    as the i^th point of the k^th step of the De Casteljau
    algorithm, and obtain β_K(t;b_0,…,b_K) = x_0^ K.

    The De Casteljau algorithm is illustrated on Fig. <ref> for
    a Euclidean cubic Bézier curve β_3(t;b_0,b_1,b_2,b_3). The general
    cubic Bézier curve can be explicitly expressed on a manifold  as
    
    β_3(t;x^ 0_0,x^ 0_1,x^ 0_2,x^ 0_3)
               = g
            (t;
              g(t;
                g(t;x^ 0_0,x^ 0_1)
                ,
                g(t;x^ 0_1,x^ 0_2)
              )
              , 
                       g(t;
                g(t;x^ 0_1,x^ 0_2)
                ,
                g(t;x^ 0_2,x^ 0_3)
              )
            )
            
       = g(t;
              g(t;x^ 1_0,x^ 1_1),
              g(t;x^ 1_1,x^ 1_2)
            )
            
       = g(t; x^ 2_0,x^ 2_1)
            
       =x_0^ 3.

    

    The conditions of continuity and differentiability are generalized
    to manifolds in <cit.>.
    
      
      Consider the composite Bézier
      
curve [0,2] →ℳ
      consisting of
      f  [0,1] →, t ↦
          f(t) = β_K(t;b_0^0,b_1^0,…,b_K^0) and
      g  [1,2] →, t ↦
          g(t) = β_K̅(t-1;b_0^1,b_1^1,…,b_K̅^1), i.e.
      
    (t) 
              f(t)     for  t ∈ [0,1], 
    
              g(t)     for  t ∈ (1,2].

      The composite Bézier curve (t) is continuous and
      continuously differentiable if the two following conditions hold:
      
    b_K^0 = b_0^1
                and    
            b_K^0    = g(s; b_K-1^0, b_1^1),
               s = K/K + K̅.

    

    The composite Bézier curve [0,n] →ℳ is then
    defined completely analogously to the Euclidean case from
    Eq. (<ref>). The differentiability
    conditions (<ref>) of Lemma <ref> have to be
    satisfied at each junction point p_i,
    i=1,…,n-1.

    
    [Composite cubic Bézier curves]
      The composite cubic Bézier

      curve (t)  [0,n] →
      is C^1 if p_i = g(1/2; b_i^-, b_i^+).
      See Fig. <ref> for an example on = with n=5
      segments
      and Fig. <ref> for an example on =𝕊^2
      with n=3 segments.
    
  
  
  

 §.§ Discrete approximation of the mean squared acceleration

    
    We discretize the mean squared acceleration (MSA) of a
    curve γ [0,1]→, by discretizing the corresponding integral
    
    ∫_0^1 ‖D^2γ(t)/ṭ^2‖^2_γ(t)ṭ,

    , the regularizer from (<ref>).
    We approximate the squared norm of the
    second (covariant) derivative by the second order absolute finite difference
    introduced by Bačák  <cit.>.
    Consider three points x,y,z ∈ and the set of mid-points of x and z
    
    𝒞_x,z{
            c
             | 
            c = g(12;x,z)
            },

    for all (not necessarily shorest) geodesics g connecting x and z.
    The manifold-valued second order absolute finite differences
    is defined by
    
    _̣2[x,y,z] = min_c∈𝒞_x,z2d_(c,y).

    This definition is equivalent, on the Euclidean space, to
    ‖ x - 2y + z‖ = 2‖1/2(x+z)-y‖.

    Using equispaced points t_0,…,t_N, N∈ℕ, with
    step size Δ_t = t_1-t_0 = 1/N,
    we approximate ‖D^2γ(t_i)/ṭ^2‖_γ(t_i)≈1/Δ_t^2d_2[γ(t_i-1),γ(t_i),γ(t_i+1)],
        i=1,…,N-1, and obtain by the trapezoidal rule
    
    ∫_0^1 ‖D^2γ(t)/ṭ^2‖^2_γ(t)ṭ≈∑_k=1^N-1Δ_t d_2^2[γ(t_i-1),γ(t_i),γ(t_i+1)]/Δ_t^4.

    For Bézier curves γ(t)=(t) we obtain for the regularizer
    in (<ref>) the discretized MSA A(b) that
    depends on the control points b and reads
    
    A(b) ∑_i=1^N-1^̣2_2 [
                (t_i-1),(t_i),(t_i+1)
              ]/Δ_t^3.





§ THE GRADIENT OF THE DISCRETIZED MEAN SQUARED ACCELERATION


  In order to minimize the discretized MSA A(b), we aim to employ a
  gradient descent algorithm on the product manifold ^M, where M is the number
  of elements in b. In the following, we derive a closed form of the gradient
  ∇_b A(b) of the discretized MSA (<ref>). This gradient
  is obtained by means of a recursion and the chain rule. In fact, the 
  derivative of (<ref>) is already known <cit.>,
  such that it only remains to compute the derivative of the composite Bézier curve.

  We first introduce the following notation.
  We denote by D_xf[η](x_0)∈ T_f(x_0) the
  directional derivative of fℳ→ℳ evaluated at x_0,
  with respect to its argument x and in the direction η∈ T_x.
  We use the short hand D_xf[η] = D_xf[η](x) whenever this
  directional derivative is evaluated afterwards again
  at x.
  
  We now state the two following definitions, which are crucial for the rest of this section.

  
    
    Let f→ be a real-valued function on a manifold ,
    x ∈ and η∈ T_x.

    The gradient ∇_ f(x)∈ T_xℳ of f at x
    is defined as the tangent vector that fulfills
    
    ⟨∇_f(x),η⟩_x = D_x f[η]
            for all  η∈ T_xℳ.

  
  For multivariate functions f(x,y), we denote the gradient of
  f with respect to x at (x_0,y_0) by writing
  ∇_,xf(x_0,y_0). We shorten this notation as ∇_,xf =
  ∇_ℳ,x f(x,y) when this gradient is seen as a function
  of x (and y).
   
    Let f→, h→ be two functions on a manifold  and F→,
    x↦ F(x) = (f ∘ h)(x) = f(h(x)),
    their composition.
    Letx ∈andη∈ T_x .
    The directional derivative D_xF[η] ofFwith respect toxin the directionηis given by
    
    D_x F[η] = D_h(x) f[D_x h[η]],

    whereD_x h[η] ∈ T_h(x)andD_x F[η] ∈ T_F(x).
  

  The remainder of this section is organized in four parts. We first recall the
  theory on Jacobi fields in Section <ref> and their relation to
  the differential of geodesics (with respect to start and end point). 
  In Section <ref>, we
  apply the chain rule to the composition of two geodesics, which appears
  within the De Casteljau algorithm. We use this result to build an algorithmic derivation of the differential of
  a general Bézier curve on manifolds with respect to its control points 
  (Section <ref>). We extend the result to composite Bézier
  curves in Section <ref>, including their constraints on
  junction points p_i to enforce the C^1 condition (<ref>), and finally
  gather these results
  to state the gradient ∇_ℳ A(b) of the discretized
  MSA (<ref>) with respect to the control points. 
  
  
  

 §.§ Jacobi fields as derivative of a geodesic

    In the following, we introduce a closed form of the
    differential D_x g(t;·,y) of a
    
geodesic g(t;x,y), t∈[0,1],
    with respect to its start point x∈. The differential with respect to
    the end point y∈ can be obtained by taking the
    geodesic g(t,y,x) = g(1-t; x,y).

    

    As represented in Fig. <ref>, we denote by γ_x,ξ,
    the geodesic starting in γ_x,ξ(0)=x and with
    direction D/dtγ_x,ξ(0) =ξ∈ T_x. We introduce
    ζ(s)γ_x,ξ(s)y, the tangential vector in
    T_γ_x,ξ(s) pointing towards y. Then, the
    geodesic variationΓ_g,ξ(s,t) of the geodesic
    g(·;x,y) with respect to the tangential direction ξ∈ T_x
    is given by
    
    Γ_g,ξ(s,t) γ_x,ξ(s)[tζ(s)],
                s∈(-ε,ε), t∈[0,1],

    where ε > 0. The corresponding Jacobi field J_g,ξ
    along g is then given by the vector field
    
    J_g,ξ(t)
          D/∂ sΓ_g,ξ(s,t)|_s=0

    that represents the direction of the displacement ofgifxis perturbed in a directionξ.
    
    We directly obtain J_g,ξ(0) = ξ, and J_g,ξ(1) =
    0 as well as J_g,ξ(t)∈ T_g(t;x,y). Since
    Γ_g,ξ(s,t) = g(t; γ_x,ξ(s), y) we obtain by the chain
    rule
    
    D_x g(t,·,y)[ξ]
          = d/ds  g(t; γ_x,ξ(s), y) |_s=0
          = d/dsΓ_g,ξ(s,t)|_s=0 = J_g,ξ(t).

    This relation between the derivative of geodesics and
    Jacobi fields is of high interest on symmetric spaces, where Jacobi fields
    can be computed in closed form, as summarized in the following Lemma.
  <cit.>
      Let  be a m-dimensional Riemannian manifold.
      Let g(t; x,y), t∈ [0,1], be a geodesic between x,y∈, η∈ T_x a
      tangent vector and {ξ_1,…ξ_m} be an orthonormal basis
      (ONB) of T_x that diagonalizes the curvature operator of  with eigenvalues
      κ_ℓ, ℓ=1,…,m. For details, see Ch. 4.2 and 5 (Ex. 5) of <cit.>.
      Let further denote by {Ξ_1(t),…,Ξ_m(t)} the parallel transported frame
      of {ξ_1,…,ξ_m} along g. 
      Decomposing η = ∑_ℓ=1^m η_ℓξ_ℓ∈ T_x, 
      the derivative D_x g[η] becomes
      
    D_x g(t;x,y) [η] = J_g,η(t) = ∑_ℓ=1^m η_ℓ J_g,ξ_ℓ(t),

      where the Jacobi field J_g,ξ_ℓℝ→ T_g(t;x,y) along g and in the
      direction ξ_ℓ is given by
      
    J_g,ξ_ℓ(t) =
            sinh(d_g(1-t)√(-κ_ℓ))
              /
                sinh(d_g√(-κ_ℓ))
              Ξ_ℓ(t)   κ_ℓ<0,
    sin(d_g(1-t)√(κ_ℓ)) /sin(√(κ_ℓ)d_g)Ξ_ℓ(t)   κ_ℓ>0,
    
              1-tΞ_ℓ(t)   κ_ℓ=0,

      with d_g = d(x,y) denoting the length of the
      geodesic g(t;x,y), t∈[0,1].
    

    The Jacobi field of the reversed geodesic g̅(t)  g(t;y,x) = g(1-t;x,y) is obtained using the same
    orthonormal basis and transported frame
    but evaluated at s=1-t. We thus obtain D_y
    g(t;x,y)[ξ_ℓ] = D_y g(1-t;y,x)[ξ_ℓ] = J_g̅,ξ_ℓ(1-t), where
    
    J_g̅,ξ_ℓ(1-t)
          =
          sinh(d_gt√(-κ_ℓ))
            /
              sinh(d_g√(-κ_ℓ))
            Ξ_ℓ(t)   κ_ℓ<0,
    sin(d_gt√(κ_ℓ)) /sin(d_g√(κ_ℓ))Ξ_ℓ(t)   κ_ℓ>0,
    
            tΞ_ℓ(t)   κ_ℓ=0.

  
  

 §.§ Derivative of coupled geodesics

    Letbe a symmetric Riemannian manifold.
    We use the result of Lemma <ref> to directly compute the
    derivative of coupled geodesics, , a function composed
    of g_1(t)  g(t;x,y) and g_2(t)  g(t;g_1(t),z).
    By Definition <ref>, we have
    
    D_x g_2(t)[η]
          = D_g_1(t)g(t;·,z)[D_xg_1(t)[η]]

    and by (<ref>), we obtain
    
    D_x g_2(t)[η]
          = J_g_2,D_xg_1(t)[η](t),

    where the variation direction used in the Jacobi field is now the
    derivative ofg_1(t)in directionη.
    Similarily, we compute the derivative of a reversed coupled
    geodesicg_3(t)  g(t;z,g_1(t))as
    
    D_x g_3(t)[η]
          = D_g_1(t)g(t;z,·)[D_xg_1(t)[η]]
          = J_g̅_3,D_xg_1(t)[η](1-t).

    Note that the Jacobi field is here reversed, but that its variation
    direction is the same as the one of the Jacobi field introduced forg_2(t). In a computational perspective, it means that we can use the same ONB
    for the derivatives of both g_3 and g̅_3
    Furthermore, in this case, the variation direction is also computed by
    a Jacobi field sinceD_xg_1(t)[η] = J_g_1,η(t).

    Finally the derivative ofg_2(resp.g_3) on symmetric spaces
    is obtained as follows. Let{ξ_1^[1],…,ξ_m^[1]}be an ONB ofT_xfor the
    inner Jacobi field along g_1, and{ξ_1^[2],…,ξ_m^[2]}be an ONB
    ofT_g_1(t)for the outer Jacobi
    field alongg_2(resp.g_3). Asη = ∑_ℓ=1^m η_ℓξ^[1]_ℓ∈ T_x, and statingJ_g_1,ξ^[1]_ℓ(t) = ∑_l=1^m
    μ_l ξ_l^[2]∈ T_g_1(t), the derivative ofg_2(resp.g_3)
    with respect toxin the directionη∈ T_xreads
    
    D_x g_2(t)[η]
        = ∑_l=1^m ∑_ℓ=1^m J_g_2,ξ^[2]_l(t) μ_l η_ℓ,

    and accordingly forg_3.
  
  
  

 §.§ Derivative of a Bézier curve

    Sections <ref> and <ref> introduced
    the necessary concepts to compute the derivative of a general Bézier curve
    β_K(t;b_0,…,b_K), as described in Equation (<ref>),
    with respect to its control pointsb_j.
    For readability of the recursive structure investigated in the following,
    we introduce a slightly simpler notation and the following setting.
    
    LetKbe the degree of a Bézier curveβ_K(t;b_0,…,b_K)with
    the control pointsb_0,…,b_K∈. We fix k∈{1,…,K},
    i∈{0,…,K-k} and t∈[0,1]. We introduce
    
    g_i^ k (t)  g(t; g_i^k-1(t), g_i+1^k-1(t))
            = β_i^ k (t; b_i,…,b_i+k),

    for the i^th Bézier curve of degree k in the De Casteljau
    algorithm, and g_i^0(t) = b_i. 
    
    Furthermore, givenx ∈{ b_i,…,b_i+k}, we denote by
    
    η_i^ k  D_x g_i^ k(t)[η],

    its derivative with respect to one of its control points x
    in the direction η∈ T_x. 
    
    
      Clearly any other derivative of g_i^k with respect
      to x = b_j, j<i  or j>i+k is zero. In addition we have
    η_i^0=D_xg_i^0[η] = η for x=b_i and zero otherwise.
    
      Let k∈{1,…,K} and i∈{0,…,K-k} be given.
      The derivative η_i^k = D_x g_i^ k(t)[η] of g_i^ k
      with respect to its control point x  b_j, i ≤ j ≤ i+k,
      and in the direction η∈ T_x is given by
      
    η_i^ k  D_x g_i^ k (t) [η]
            =
            
              J_g_i^k-1,η_i^k-1(t)
                    if  j = i,
    
              J_g_i^k-1,η_i^k-1(t)
                + J_g̅_i+1^k-1,η_i+1^k-1(1-t)
                    if  i<j<i+k,
    
              J_g̅_i+1^k-1,η_i+1^k-1(1-t)
                    if  j = i+k.

      Let fix t∈[0,1]
      and x=b_j, i≤ j ≤ i+k. For readability we
      set a  g_i^k-1(t), b  g_i+1^k-1(t), and f  g_i^ k(t) = g(t;a,b). Note that while
      f depends on the control points b_i,…,b_i+k and is a
      Bézier curve of degree k, both a and b are Bézier curves of
      degree k-1. The former does not depend
      on b_i+k, and the latter is independent of b_i.

      We prove the claim by induction. For k=1 the function g_i^1
      is just a geodesic. The case i<j<i+1 does not occur and the remaining
      first and third cases follow by the notation introduced for k=0 and
      Lemma <ref>.

      For k>1 we apply the chain rule (<ref>) to D_xf[η]
      and obtain
      
    D_x f[η]
            = D_a f [ D_x a [η]]
            +
            D_b f [ D_x b [η]].

      Consider the first term D_af[D_xa[η]] and j < i+k. 
      By (<ref>) and the notation
      from (<ref>), one directly has
      
    D_a f[η_i^k-1]
            = J_f, η_i^k-1(t).

      For j=i+k, clearly D_xa[η]=D_af[D_xa[η]]=0, as a does not
      depend on b_i+k.

      We proof the second term similarily. For j>i, by applying the chain
      rule and using the reversed Jacobi field formulation of
      Lemma <ref> for the derivative of a geodesic with respect to its end point, we obtain
      
    D_b f[η_i+1^k-1]
            = J_f̅, η_i+1^k-1(1-t).

      Finally, as D_xb[η]=D_bf[D_xb[η]]=0 for x=b_i, the
      assumption follows.
    
    Fig. <ref> represents one level of the
    schematic propagation tree to compute the derivative of a Bézier curve.
    [Quadratic Bézier curve]
      
      Consider the quadratic Bézier curve β_2 [0,1] →
      defined as
      
    β_2(t; b_0,b_1,b_2)
            =
            g(
            t;
            g(t; b_0,b_1),
            g(t;b_1,b_2)
            ).

      Using the notations (<ref>), we have
      
    g_0^ 1(t)     g(t; b_0, b_1),  
            g_1^ 1(t)  g(t; b_1, b_2), 
    
            g_0^ 2(t)     g(t;g_0^ 1,g_1^ 1).

      The derivative of β_2 at t with respect to b_0 in the
      direction η∈ T_b_o, is given by
      
    D_b_0β_2[η] = J_g_0^ 2, η_0^ 1(t),
             with η_0^ 1  J_g_0^ 1, η(t).

      The derivative of β_2 at t with respect to b_2 in the
      direction η∈ T_b_2 can be seen as deriving by the first
      point after inverting the Bezier curve, i.e. looking at β̅_2(t)
      = β_2(1-t), hence we have analogously to the first term
      
    D_b_2β_2[η] = J_g̅_0^ 2, η_1^ 1(1-t),
             with η_1^ 1  J_g̅_1^ 1, η(1-t).

      The case D_b_1β_2[η], η∈ T_b_1, involves a
      chain rule, where b_1 appears in both g_0^ 1 (as its starting point) and g_1^
      1 (as its end point). Using the two intermediate results (or Jacobi fields of geodesics)
      
    η_0^ 1   J_g̅_0^ 1, η(1-t)
              and  η_1^ 1   J_     g_1^ 1, η(  t),

      we obtain
      
    D_b_1β_2[η] =
            J_g_0^ 2, η_0^ 1(t)
            +
            J_g̅_0^ 2, η_1^ 1(1-t).
[Cubic Bézier curve]
      Consider the cubic Bézier curve β_3 [0,1] → defined
      as
      
    β_3(t; b_0,b_1,b_2,b_3)
            =
            g(
              t;
              β_2(t; b_0,b_1,b_2),
              β_2(t; b_1,b_2,b_3)
            ).

      As in Example <ref>, we use the
      notations (<ref>) and define
      
    g_j^ 1(t)     g(t; b_j, b_j+1),      j = 0,1,2, 
    
            g_j^ 2(t)     g(t; g_j^ 1, g_j+1^ 1),
                   j = 0,1,  and 
    
            g_0^ 3(t)     g(t;g_0^ 2,g_1^ 2).


      The derivation of β_3 with respect to b_0 or b_3 follows
      the same structure as in Example <ref>. The case of D_b_1β_3[η], however, requires two chain rules. The needed
      Jacobi fields follow the tree structure shown in Fig. <ref>:
      given η∈ T_b_1, we define at the first recursion step
      
    η_0^ 1   J_g̅_0^ 1, η(1-t),
                η_1^ 1   J_     g_1^ 1, η(  t),

      and at the second recursion step
      
    η_0^ 2   J_g_0^ 2, η_0^ 1(t)
              + J_g̅_0^ 2, η_1^ 1(1-t),
                η_1^ 2   J_g_1^ 2, η_1^ 1(t).

      Note that both η_0^ 2 and η_1^ 2 are actually the derivative of
      β_2(t;b_0,b_1,b_2) and β_2(t;b_1,b_2,b_3), respectively, with respect to b_1 and direction
      η∈ T_b_1. Finally we have
      
    D_b_1β_3[η] =
            J_g_0^ 3, η_0^ 2(t)
            +
            J_g̅_0^ 3, η_1^ 2(1-t).

      The case of D_b_2β_3[η] is obtained symmetrically, with
      arguments similar to b_1.
    
  
  

 §.§ Joining segments and deriving the gradient

    In this subsection we derive the differential of a composite Bézier
    curve (t) consisting ofnsegments and take the C^1
    conditions into account. 
    We simplify the notations from
    Section <ref> and set the degree fixed to K_i=K for
    all segments, , K=3 for a cubic composite Bézier curve. Then the
    control points are b_j^i, j=0,…,K, i=0,…,n-1. We further
    denote by p_i=b_K^i-1=b_0^i, i=1,…,n-1 the common junction
    point of the segments and p_0 and p_n the start and end points,
    respecively. For ease of notation we denote by b_i^- = b_K-1^i and
    b_i^+=b_1^i,i=1,…,n-1, the two points needed for differentiability (C^1)
    condition investigation, cf. Fig. <ref> for an illustration of
    the framework on = , withK=3.

    One possibility to enforce the C^1 condition (<ref>) is to
    include it into the composite Bézier curve by replacing b_i^+ with
    
    b_i^+ = g(2; b_i^-,p_i),   i=1,…,n-1.

    This way both the directional derivatives of(t)with respect
    to b_i^+ and p_i change due to a further (most inner) chain rule.
    
      Let  be a composite Bézier curve and p_i,b_i^+,b_i^-
      introduced as above. Replacing b_i^+= g(2;b_i^-,p_i) eliminates that
      variable from the composite Bezier curve and keeps the remaining
      differentials unchanged, despite the following which now read
      
    D_b_i^-(t)[η] =
            
          D_b_i^-β_K(t-i+1; p_i-1,b_i-1^+,…,·,p_i)[η]
               t∈ (i-1,i],
    
          D_b_i^+β_K(t-i; p_i,b_i^+,…,·,p_i+1)
                  [
                    D_b_i^-g(2; ·,p_i)[η]
                ]
                   t∈ (i,i+1],

      and
      
    D_p_i(t)[η] =
            
          D_p_iβ_K(t-i+1; p_i-1,b_i-1^+,…,b_i^-,·)[η]
               t∈ (i-1,i],
    
          D_p_iβ_K(t-i; ·,b_i^+,…,b_i+1^-,p_i+1)[η]
      +
          D_b_i^+β_K(t-i; p_i,·,…,b_i+1^-,p_i+1)
                  [
                    D_p_ig(2; b_i^-,·)[η]
                ]
               t∈ (i,i+1].

      In both cases, the first interval includes i-1=0, when i=1.
    
      Both first cases are from the derivation of a bezier curve as before,
      for both second cases replacing b_i^+ = g(2; b_i^-,p_i) yields one
      (for p_i additional) term.
    
  
    We now derive the gradient of the objective function (<ref>).
    We introduce the abbreviation _i = (t_i) ∈,
    and ^̣2_2,i = ^̣2_2 ( _i-1,_i,_i+1 ).
    
      Let  be a m-dimensional manifold, x be one of the control
      points of a composite Bézier curve , 
      and {ξ_1,…,ξ_m } be a corresponding orthonormal basis (ONB) of T_x.
      The gradient ∇_,xA(b) of the discretized mean squared
      acceleration A(b) of , w.r.t. x, discretized
      at N+1 equispaced times t_0,…,t_N is given by
      
    ∇_,xA(b) =
            ∑_ℓ=1^m
            ∑_i=1^N-1∑_j = i-1^i+1⟨∇_^̣2_2,j, D_x _j[ξ_ℓ]
              ⟩__jξ_ℓ.

      As ∇_,xA(b) ∈ T_x, we seek for the
      coefficients a_ℓ a_ℓ(x) such that
      
    ∇_ f(x) = ∑_ℓ=1^m a_ℓξ_ℓ.

      Therefore, for any tangential vector
      η∑_ℓ=1^m η_ℓξ_ℓ∈ T_x,
      we have,
      
    ⟨∇_,xA(b),η⟩_x  = ∑_ℓ=1^m a_ℓη_ℓ.

      By definition of A and Definition <ref> this yields
      
    ⟨∇_,xA(b),η⟩_x
              = ∑_i=1^N-1⟨∇_^̣2_2,i(x),η⟩_x
              = ∑_i=1^N-1 D_x ^̣2_2,i [η].

      We compute D_x^̣2_2,i [η] using the chain rule (Definition <ref>) as
      
    D_x ^̣2_2,i [η] =
              ∑_j = i-1^i+1D__j^̣2_2,j[ D_x _j[η]],

      which, by Definition <ref>, again becomes
      
    D__j^̣2_2,j[ D_x _j[η]]
              = ⟨∇_^̣2_2,j , D_x _j[η]
              ⟩__j.

      The term on the left of the inner product is given
      in <cit.> and the right term is given in Section <ref>. 
      While the former can be computed using  Jacobi fields and a logarithmic map, 
      the latter is the iteratively coupling of Jacobi fields.
      Furthermore, the differential D_x_j[η] can be written as
      
    D_x _j[η] = ∑_ℓ=1^m η_ℓ D_x _j [ξ_ℓ] ∈ T__j.

      Hence, we obtain
      
    D__j^̣2_2,j[ D_x _j[η]]
            = ∑_ℓ=1^m η_ℓ⟨∇_^̣2_2,j,
              D_x _j[ξ_ℓ]
            ⟩__j,

      and by (<ref>), (<ref>) and (<ref>), it follows
      
    ⟨∇_,xA(b),η⟩_x =
                ∑_ℓ=1^m η_ℓ∑_i=1^N-1∑_j = i-1^i+1⟨∇_^̣2_2,j,
                  D_x _j[ξ_ℓ] ⟩__j,

      which yields the assertion (<ref>).
    



§ APPLICATION TO THE FITTING PROBLEM

  The fitting problem has been tackled different
  ways this last decades. The approach with Bézier curves is more recent, and
  we refer to <cit.> for a detailed 
  overview of these methods.
  
  In this section, we present the numerical framework we use in order to
  fit a composite Bézier curveto a set of data pointsd_0,…,d_n ∈associated with time-parameterst_0,…,t_n, 
  such that we meet (<ref>). 
  For the sake of simplicity, we limit the study to the case wheret_i = i,i=0,…,n. Therefore, the fitting problem (<ref>) becomes
  
    min_b∈Γ_ E_λ(b) ∫_t_0^t_n‖D^2 (t)/ṭ^2‖^2_(t)ṭ
        +
        λ/2∑_i=0^n ^̣2(p_i,d_i),

  whereΓ_∈^Mis the set of theMcontrol points of.
  Remark that, compared to (<ref>), the optimization is now done on the
  product manifold^M. Furthermore, the fitting term now implies a distance
  betweend_iandp_i, as(t_i) = p_i.
  
  The section is divided in three parts: the product manifold^Mis
  defined in Section <ref>, where the contribution of the
  fitting term in the gradient ofEis also presented. 
  Then, we propose an efficient algorithm to compute the gradient of
  the discretized MSA, based on so-called adjoint Jacobi fields.
  We finally shortly mention the gradient descent algorithm we use as well as
  the involved Armijo rule.
  
  
  
  

 §.§ Fitting and interpolation

    Let us clarify the setΓ_∈^Mfrom (<ref>). 
    We will explicitly present the
    vector b and state its sizeM.
    The setΓ_is the set of theMremaining free control
    points to optimize,
    when the^1continuity constraints are imposed. We distinguish two cases:
    (i) the fitting case, that corresponds to formulae presented in
    Section <ref>, and (ii) the interpolation case
    (λ→∞) where the constraintd_i = p_iis imposed as well.
    
    For a given composite Bézier curve [0,n]→
    consisting of a Bézier curve of degree K on each
    segment, and given the C^1 conditions (<ref>), the segments
    are determined by the points
    
    b = (p_0,b_0^+,b_2^0,…,b_K-2^0,b_1^-,p_1,…,b_n^-,p_n)∈^M.

    We investigate the length of M. First, b_i^+ is given by p_i-1
    and b_i^- via (<ref>).
    Second, for the segments i = 2,…,n we can further omit p_i-1
    since the value also occurs in the segmentias last entry.
    Finally, the first segment contains
    the additional value of b_0^+ that is not fixed by C^1 constraints.
    The first segment is thus composed ofK+1control points, while then-1remaining
    segments are determined by K-1 points. In total we obtain M=n(K-1)+2
    control points to optimize.
    
    MinimizingA(b)alone leads to the trivial solution, for any set
    of control pointsb=(x,…,x),x∈,
    and this is why the fitting term from (<ref>) is important.

    

  
Fitting (0 < λ < ∞).
    If the segment start and end points are obstructed by noise or allowed to
    move, we employ a fitting scheme to balance the importance given to the data pointsd_0,…,d_n.
    Equation (<ref>) reads
    
    _b∈Γ_Ã(b),   Ã(b) 
            A(b) + λ/2∑_i=0^n d^2_(d_i,p_i),

    where λ∈^+ sets the priority 
    to either the data term (large λ) or the mean squared acceleration
    (small λ) within the minimization.
    The gradient of the data term is given in <cit.>,
    and the gradient of Ã is given by
    
    ∇_^M,xÃ(b)
          =
          ∇_^M,xA(b) -λlog_p_id_i
                   if  x=p_i, i=0,…,n,
    ∇_^M,xA(b)    otherwise.


  
Interpolation (λ→∞).
    For interpolation we assume that the start pointp_i-1and end pointp_iof the segmentsi=1,…,nare fixed to given datad_i-1andd_i∈ℳ, respectively.
    The optimization of the discrete mean squared accelerationA(b)reads
    
    _b∈Γ_ A(b)   s. t.  p_i=d_i, i=0,…,n.

    Since thep_iare fixed by constraint, they can be omitted from the
    vectorb. We obtain
    
    b =
          (b_0^+,b_2^0,…,b_1^-,b_2^1,…,b_n^-) ∈^M'

    Since there aren+1additional constraints, the minimization is hence performed on the product 
    manifold^M',M'=M-(n+1)=n(K-2)+1.
  
  
  

 §.§ Adjoint Jacobi fields

    In the Euclidean spaceℝ^m, the adjoint operatorT^*of a
    linear bounded operatorTℝ^m→ℝ^qis the operator
    fulfilling
    
    ⟨ T(x),y ⟩_ℝ^q = ⟨ x,T^*(y) ⟩_ℝ^m
          ,   for all x∈ℝ^m, y∈ℝ^q.

  
    The same can be defined for a linear operatorS T_x→ T_y,x,y∈ℳ, on am-dimensional Riemannian manifold . The adjoint
    operator S^* T_y→ T_x satisfies
    
    ⟨ S(η),ν⟩_y = ⟨η, S^*(ν)⟩_x,
             for all η∈ T_x, ν∈ T_y.

    
    We are interested in the case whereSis the differential
    operatorD_xof a geodesicF(x) = g(t;x,y)for some fixedt∈ℝandy∈.
    The differential D_xF T_x→ T_F(x) can be written as 
    
    D_xF[η]
          =
          J_g,η(t)
          =
          ∑_ℓ=1^m⟨η,ξ_ℓ⟩_xα_ℓΞ_ℓ(t),

    where α_ℓ are the coefficients of the Jacobi
    field (<ref>), and ξ_ℓ, Ξ_ℓ(t) are given as in
    Lemma <ref>.
    To derive the adjoint differential(D_xF)^* T_F(x)→ T_x
    we observe that for anyν∈ T_F(x)we have
    
    ⟨ D_xF[η],ν⟩_F(x)
          = ∑_ℓ=1^m ⟨η,ξ_ℓ⟩_xα_ℓ⟨Ξ_ℓ(t),ν⟩_F(x)
          = ⟨η, ∑_ℓ=1^m ⟨Ξ_ℓ(t),ν⟩_F(x)α_ℓξ_ℓ⟩_x.

    Hence the adjoint differential is given by
    
    (D_xF)^*[ν] = ∑_ℓ=1^m
            ⟨ν,Ξ_ℓ(t)⟩_F(x)α_ℓξ_ℓ,
              ν∈ T_F(x).

    We introduce
    the adjoint Jacobi field J^*_F,νℝ→
    T_xℳ, ν∈ T_F(x) as
    
    J^*_F,ν(t) = ∑_ℓ=1^m ⟨ν,Ξ_ℓ(t)⟩_F(x) J^*_F,Ξ_ℓ(t)(t)
          = ∑_ℓ=1^m ⟨ν,Ξ_ℓ(t)⟩_F(x)α_ℓξ_ℓ.

    Note that evaluating the adjoint Jacobi fieldJ^*involves
    the same transported

    frame {Ξ_1(t),…Ξ_m(t)}
    and the same coefficientsα_ℓas the Jacobi fieldJ,
    which means that the evaluation of the adjoint is in no way inferior to
    the Jacobi field itself.
    
    
    The adjointD^*of the differential is useful in particular, when computing the
    gradient ∇_(h∘ F) of the composition of Fℳ→ℳ
     with h→.
    Settingy  F(x) ∈, we obtain for any η∈ T_xℳ that
    
    ⟨∇_,x(h∘ F)(x),η⟩_x
             = D_x(h∘ F)[η] 
       = D_F(x)h[ D_x F[η] ] 
       = ⟨∇_,yh(y), D_xF[η] ⟩_y
       = ⟨ (D_xF)^*[∇_,y h(y)], η⟩_x.

    
    Especially for the evaluation of the gradient of the composite functionh∘ Fwe obtain
    
    ∇_,x(h∘ F)(x) =
          (D_xF)^*[∇_,y h(y)] = J^*_F,∇_,yh(y)(t).

    The main advantage of this technique appears in the case of composite
    functions, , of the form h∘ F_1∘ F_2 (the generalization to composition withKfunctions is straightforward). 
    The gradient
    ∇_,x (h∘ F_1 ∘ F_2)(x) now reads,
    
    ∇_,x (h∘ F_1 ∘ F_2)(x) = (D_xF_2)^*[∇_,y_2 h∘ F_1(y_2)] = J^*_F_2,∇_,y_2 h ∘ F_1(y_2))(t).

    The recursive computation ofη^ 3 =∇_,x h (x)is then given by the following algorithm
  
    η^ 1    =   ∇_,y_1h(y_1), 
    η^ 2    =    J^*_F_1,η^ 1(t), 
    η^ 3    =    J^*_F_2,η^ 2(t).

    For h=d_2^2^3→ℝ we know ∇_^3h by
    Lemma 3.1 and Lemma 3.2 from <cit.>. Let t_1,t_2,t_3 ∈ [0,1]
    be time points, and b∈^M be a given (sub)set of the control points
      of a (composite) Bézier curve .
      We define F^M→^3,b↦ F(b) =
      ((t_1),(t_2),(t_3)) as the evaluations
      of  at the three given time points.
      The composition h∘ F hence consists of (in order of evaluation) the
      geodesic evaluations of the De Casteljau algorithm, the mid point function
      for the first and third time point and a distance function.
      
      The recursive evaluation of the gradient starts with the
      gradient of the distance function.
      Then, for the first and third arguments, a mid point Jacobi field
      is applied. The result is plugged into the last geodesic evaluated within the
      De Casteljau “tree” of geodesics. At each geodesic, a tangent vector at
      a point g(t_j;a,b) is the input for two adjoint Jacobi fields, one
      mapping to T_a, the other to T_b. This information is available
      throughout the recursion steps anyways. After traversing this tree backwards,
      one obtains the required gradient of h∘ F.
    
    
    Note also that even the differentiability constraint (<ref>)
    yields only two further (most outer) adjoint Jacobi fields,
    namely J^*_g(2,b_i^-,p_i),∇_,b_i^+(t_j)(2)
    and J^*_g̃(2,b_i^-,p_i),∇_,b_i^+(t_j)(2).
    They correspond to variation of the start point b_i^- and the end
    point p_i, respectively as stated in (<ref>).
  
  
  

 §.§ A gradient descent algorithm

    To address (<ref>) or (<ref>),
    we use a gradient descent algorithm, as described in <cit.>. For completeness
    the algorithm is given in Algorithm <ref>.
    
     The step sizes are given by the Armijo line search condition 
     presented in <cit.>. Let𝒩be a Riemannian manifold, x=x^(k)∈𝒩 be an iterate of the
     gradient descent, and β,σ∈(0,1), α>0. Let m be the
     smallest positive integer such that
     
    F(x) - F(
             exp_x(-β^mα∇_𝒩F(x))
           )
           ≥σβ^mα‖∇_𝒩F(x)‖_x.

     We set the step size to s_k β^mα
     in Algorithm <ref>.

     As a stopping criterion we use a maximal number k_max of iterations or
     a minimal change per iteration d_𝒩(x^k,x^k+1)<ε. In practice, this last criterion is matched first.

     The gradient descent algorithm converges to a critical
     point if the functionFis convex <cit.>.
     The mid-points
     model (<ref>) posesses two advantages: (i) the
     complete discretized MSA (<ref>) consists of (chained)
     evaluations of geodesics and a distance function, and (ii) it reduces
     to the classical second order differences on the Euclidean space. However, this model
     is not convex on general manifolds.
     An example is given in the arXiv preprint (version 3)
     of <cit.>[see https://arxiv.org/abs/1506.02409v3
     arxiv.org/abs/1506.02409v3], Remark 4.6. Another possibility (also
     reducing to classical second order differences in Euclidean space) is the
     so-called Log model (see, , <cit.>)
     
    d_2,Log[x,y,z]  = ‖log_yx+log_yz‖_y.

     The (merely technical) disadvantage of the Log model is
     that the computation
     of the gradient involves further Jacobi fields than the one presented above,
     namely to compute the differentials of the logarithmic map both with respect
     to its argumentD_xlog_yxas well as its base pointD_ylog_yx.
     Still, these can be given in closed form for symmetric Riemannian
     manifolds <cit.><cit.>.
     To the best of our knowledge, the joint convexity of the Log model
     inx,yandzis still an open question.
     




§ EXAMPLES

  In this section, we provide several examples of our algorithm
  applied to the fitting problem (<ref>).
  
  We validate it first on the Euclidean space and verify that it retrieves
  the natural cubic smoothing spline. We then present examples on the
  sphere𝕊^2and the special orthogonal groupSO(3).
  We compare our results with the fast algorithm of Arnould  <cit.>,
  generalized to fitting, and the so-called blended
  cubic splines from <cit.>.
  The control points of the former are obtained by generalizing the optimal
  Euclidean conditions of (<ref>) (in^m, this is a 
  linear system of equations) to the manifold setting; 
  the curve is afterwards reconstructed by a classical De Casteljau algorithm.
  In the latter, the curve is obtained as a blending of solutions computed
  on carefully chosen tangent spaces, , Euclidean spaces.
  
  The following examples were implemented in MVIRT <cit.>[
  open source, available at <http://ronnybergmann.net/mvirt/>], and the
  comparison implementations from <cit.> and <cit.> use
  Manopt <cit.>[ open source, available
  at <http://www.manopt.org>]. Note that both toolboxes use a very similar
  matrix structure and are implemented in Matlab, such that the results can
  directly be compared.

  

 §.§ Validation on the Euclidean space

  As a first example, we perform a minimization on the Euclidean space ℳ=ℝ^3.
  A classical result on^mis that the curveγminimizing (<ref>) is the natural
  (^2) cubic spline whenΓis a Sobolev spaceH^2(t_0,t_n).
  We compare our approach to the tangential linear system approach derived 
  in <cit.> in order to validate our model. 
  We use the following data points
  
    d_0 = [ 0; 0; 1 ],  
        d_1 = [  0; -1;  0 ],  
        d_2 = [ -1;  0;  0 ],  
        d_3 = 1/√(82)[  0; -1; -9 ],

  as well as the control pointsp_i = d_i, and
  
    b_0^+    = exp_p_0π/8√(2)[  1; -1;  0 ],
           b_1^+    = exp_p_1 -π/2√(2)[ 1; 0; 1 ],
           b_1^-    = g(2;b_1^+,p_1),
    
        b_2^+    = exp_p_2π/2√(2)[  0;  1; -1 ],
           b_3^-    = exp_p_3π/8[ -1;  0;  0 ],
           b_2^-    = g(2;b_2^+,p_2),

  where the exponential map and geodesic on ℝ^3
  are actually the addition
  and line segments, respectively. Note that, by construction, the 
  initial curve is continuously differentiable
  but (obviously) does not minimize (<ref>). The parameter λ is set to 50.
  The MSA of this initial curve Ã(b) is
  approximately 18.8828.

  The data points are given to the algorithm of <cit.> to reconstruct
  the optimal Bézier curve(t), which  is the natural^2cubic spline.
  The result is shown in
  Fig. <ref> together with the first order
  differences along the curve in Fig. <ref>
  as a dotted curve.

   

  The set of control points (p_0,b_0^+,…,b_3^-,p_3) is
  optimized with our proposed method. We discretize the second
  order difference using N=1600 points.
  The resulting curve and the first order difference plots
  are also obtained using these sampling values and a first order forward
  difference.
  The gradient descent algorithm from Algorithm <ref>
  employs the Armijo rule (<ref>) setting β=1/2, σ=10^-4,
  and α=1.
  The stopping criteria are ∑_i=1^1600 d_ℳ(x^(k)_i,x^(k)_i) < ϵ = 10^-15 or
  ‖∇_ℳ^nÃ‖_2 < 10^-9, and the algorithm 
  stops when one of the two is met. For this example, the
  first criterion stopped the algorithm, while the norm of the gradient was of
  magnitude10^-6.

  Both methods improve the initial functional value of Ã(b) ≈
  18.8828 to a value of Ã(b_min) ≈ 4.981218.
  Both the linear system approach and the gradient descent perform equally.
  The difference of objective value is
  2.4524×10^-11 smaller for the gradient descent, and the maximal
  distance of any sampling point of the resulting curves is of size
  4.3×10^-7. Hence, in the Euclidean space, the proposed gradient descent
  yields the natural cubic spline, as  one would expect.

  
  

 §.§ Examples on the sphere 𝕊^2


  
Validation on a geodesic.
  As a second example with a known minimizer,
  we consider the manifold ℳ=𝕊^2, , the two-dimensional unit sphere embedded
  in ℝ^3, where geodesics are great arcs. We use the data points
  
    d_0 = [ 0; 0; 1 ],  
        d_1 = [ 0; 1; 0 ],  
        d_2 = [  0;  0; -1 ]

  aligned on the geodesic connecting the
  north pole p_0 and the south pole p_2, and running through a point p_1
  on the equator. We define the control points of the cubic Bézier curve as follows:
  
    x_0 = 1/√(6)[ 1; 1; 2 ],
               
        b_0^+ = p_0[3p_0[x_0]] ,
              
        b_1^- = 1/√(6)[ 1; 2; 1 ],
        
    
        x_2 = 1/√(6)[ -1;  1; -2 ], 
               
        b_1^+ = 1/√(6)[ -1;  2; -1 ],
              
        b_2^- = p_2[1/3p_0[x_2]],

  wherex_0andx_2are temporary points andp_i = d_i. We
  obtain two segments smoothly connected since log_p_1b_1^- = -log_p_1b_1^+.
  The original curve is shown in Fig. <ref>, where especially
  the tangent vectors illustrate the different speed atp_i.
  
  The control points are optimized
  with our interpolating model, , we fix the start and end
  points p_0,p_1,p_2 and minimize A(b).
  
  The curve, as well as the second and first order differences, is sampled
  with N=201 equispaced points.  
  The parameters of the Armijo rule are again set
  to β=1/2, σ=10^-4, and α=1.
  The stopping criteria are slightly relaxed to 10^-7 for the distance
  and 10^-5 for the gradient, because of the sines and cosines involved
  in the exponential map.
  

  The result is shown in Fig. <ref>. Since the minimizer is
   the geodesic running from p_0 to p_2
  through p_1, we measure the perfomance first by looking at the
  resulting first order difference, which is constant, as can be seen in Fig. <ref>.
  As a second validation, we observe that the maximal distance of the resulting
  curve to the geodesic is of 2.2357×10^-6. These evaluations again
  validate the quality of the gradient descent.


  

  
Effect of the data term.
  As a third example we investigate the effect of λ in the
  fitting model. We consider the data points
  
    d_0 = [ 0 0 1 ]^,     
        d_1 = [  0 -1  0 ]^,     
        d_2 = [ -1  0  0 ]^,     
        d_3 = [  0  0 -1 ]^,

  as well as the control pointsp_i = d_i, and
  
    b_0^+    = p_0[π/8√(2)[  1; -1;  0 ]],    
        b_1^+    = p_1[-π/4√(2)[ -1;  0;  1 ]], 
    
        b_2^+    = p_2[π/4√(2)[  0;  1; -1 ]],    
        b_3^-    = p_3[-π/8√(2)[ -1;  1;  0 ]].

  The remaining control pointsb_1^-andb_2^-are given by
  the^1conditions (<ref>).
  The corresponding cuve 𝐁(t), t∈ [0,3], is shown in
  Fig. <ref> in dashed black. When computing a minimal MSA curve that
  interpolates 𝐁(t) at the points p_i, i=0,…,3, the
  acceleration is not that much reduced, see the blue curve in Fig. <ref>.

  For fitting, we consider different values of λ
  and the same parameters as for the last example. 
  The optimized curve fits the data points closer and closer
  asλgrows, and the limit
  λ→∞ yields the interpolation case. On the other hand smaller
  values ofλyield less fitting, but also a smaller value of the
  mean squared acceleration. In the limit case, ,λ=0, the
  curve (more precisely the control points of the Bézier curve) just follows the
  gradient flow to a geodesic.
    
  The results are collected in Fig. <ref>.
  In Fig. <ref>, the original curve (dashed) is shown together
  with the solution of the interpolating model (solid blue) and the gradient
  flow result, , the solution of the fitting  model (solid black)
  with λ=0. Fig. <ref> illustrates the effect
  of λ even further with a continuous variation of λ
  from a large value ofλ=10to a small value ofλ=0.01. The image is generated by sampling this range with1000equidistant values colored in the colormap .
  Furthermore, the control points are also shown in the same color. 
  
  Several corresponding functional values, see Tab. <ref>, further
  illustrate that with smaller values of λ the discretized MSA also
  reduces more and more  to a geodesic. For λ=0 there is no coupling
  to the data points and hence the algorithm does not restrict the position of
  the geodesic. In other words, any choice for the control points that yields a
  geodesic, is a solution. Note that the gradient flow still chooses a reasonably near
  geodesic to the initial data (Fig. <ref>, solid black).
  
  

  
Comparison with the tangential solution.
  In the Euclidean space,
  the method introduced in <cit.> and the
  gradient descent proposed here yield the same curve, , the natural cubic spline.
  On Riemannian manifolds, however, all approaches approximate the optimal solution. Indeed,
  their method provides a solution by working in different tangent spaces,
  and ours minimizes a discretization of the objective functional.

  In this example we take the same data points d_i as in (<ref>) now
  interpreted as points on ℳ=𝕊^2. 
  Note that the control points constructed in (<ref>)
  still fit to the sphere since each b_i^± is built with a vector
  in T_p_iℳ and using the corresponding exponential map.
  The result is shown in Fig. <ref>.
  The norm of the first order differences is given in Fig. <ref>. 
  The initial  curve (dashed black) has an objective value of 10.9103.
  The tangent version (dotted) reduces this
  value to 7.3293. The proposed method (solid blue)
  yields a value of 2.7908, regardless of whether the starting curve is
  the initial one, or the one already computed by <cit.>.
  Note that, when p_3 = [0,0,-1]^, the tangent version is even not able to
  compute any result, since the construction is performed in the tangent space
  of p_0 and since log_p_0p_3 is not defined.

  

 §.§ An example of orientations

Finally we compare our method with the blended splines introduced in <cit.>
for orientations, , data given onSO(3).
Let

    R_xy(α) =[  cosα  sinα     0; -sinα  cosα     0;     0     0     1 ], R_xz(α) =[  cosα     0  sinα;     0     1     0; -sinα     0  cosα ], R_yz(α) =[     1     0     0;     0  cosα  sinα;     0 -sinα  cosα ]

denote the rotation matrices in thex-y,x-z, andy-zplanes,
respectively. 
We introduce the three data points

    d_0    = R_xy(4π/9)R_yz(-π/2), 
    
      d_1    = R_xz(-π/8)R_xy(π/18)R_yz(-π/2),
    
      d_2    = R_xy(5π/9)R_yz(-π/2).

These data points are shown in the first line of Fig. <ref> in cyan.

We setλ=10and discretize (<ref>) with N=401 equispaced points.
This sampling is also
used to generate the first order finite differences. The parameters of the
gradient descent algorithm are set to the same values as for the examples on the sphere.

We perform the blended spline fitting with two 
segments and cubic splines. The resulting control points are shown in the second
row of Fig. <ref>, in green. The objective value isÃ(b) = 0.6464.
We improve this solution with the minimization problem (<ref>) on
the product manifoldSO(3)^6.
We obtain the control points shown in the last line of 
Fig. <ref> and an objective value ofÃ(b̂) = 0.2909for the resulting minimizerb̂.

We further compare both curves by looking at their absolute first order differences.
In the third line of Fig. <ref>, we display17orientations
along the initial curve with its first order difference as height. The line without
samples is the absolute first order differences of the minimizerb̂,
scaled to the same magnitude. The line is straightened especially for the first
Bézier segment.
Nevertheless, the line is still bent a little bit and
hence not a geodesic. This can be seen in the fourth line which represents17samples of the minimizing curve compared to its control points in the last
line, which are drawn on a straight line.



§ CONCLUSION AND FUTURE WORK

  In this paper, we introduced a method to solve the curve fitting
  problem to data points on a manifold, using composite Bézier curves.
  We approximate the mean squared acceleration of the curve 
  by a suitable second order difference and a trapezoidal rule, 
  and derive the corresponding gradient with respect to its control points.
  The gradient is computed in closed form by exploiting the recursive structure
  of the De Casteljau algorithm. Therefore, we obtain
  a formula that reduces to a concatenation of adjoint Jacobi fields.
  
  The evaluation of Jacobi fields is the only additional requirement compared
  to previous methods, which are solely evaluating exponential and logarithmic maps.
  For these, closed forms are available on symmetric manifolds.
  
  On the Euclidean space our solution reduces to the natural smoothing spline, the unique
  acceleration minimizing polynomial curve.
  The numerical experiments further confirm that the method presented in this paper 
  outperforms the tangent space(s)-based approaches with respect to the functional value.
  
  It is still an open question whether there exists a second order absolute
  finite difference on a manifold, that is jointly convex. Then convergence would follow by standard
  arguments of the gradient descent. For such a model, another interesting
  point for future work is to find out whether only an approximate evaluation of the Jacobi fields
  suffices for convergence. This would mean that, on manifolds where the Jacobi field
  can only be evaluated approximately by solving an ODE,
  the presented approach would still converge.




 §.§ Acknowledgements

  This work was supported by the following fundings: RB gratefully acknowledges
  support by DFG grant BE 5888/2–1; PYG gratefully acknowledges support by
  the Fonds de la Recherche Scientifique – FNRS and 
  the Fonds Wetenschappelijk Onderzoek – Vlaanderen under EOS Project no 30468160, 
  and “Communauté française de Belgique - Actions de Recherche
  Concertées” (contract ARC 14/19-060).
  Fig. <ref> is based on an idea of Rozan I. Rosandi.