{'fraction_non_alphanumeric': 0.051455605365274555, 'fraction_numerical': 0.01570778035685809, 'mean_word_length': 3.9099569359258775, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Total Generalized Variation for Manifold-valued Data

In this paper we introduce the notion of second-order total generalized variation (TGV) regularization for manifold-valued data in a discrete setting. We provide an axiomatic approach to formalize reasonable generalizations of TGV to the manifold setting and present two possible concrete instances that fulfill the proposed axioms. We provide well-posedness results and present algorithms for a numerical realization of these generalizations to the manifold setup. Further, we provide experimental results for synthetic and real data to further underpin the proposed generalization numerically and show its potential for applications with manifold-valued data.


Introduction
In this work we introduce and explore a generalization of second-order total generalized variation (TGV) regularization for manifold-valued data. The TGV functional has originally been introduced in [29] for the vector space setting as generalization of the total variation (TV) functional [73], which is extensively used for regularization in image processing and beyond. The advantage of TV regularization compared to, e.g., classical H 1 regularization approaches, is that jump discontinuities can be much better reconstructed. This can be seen in the function space setting since functions of bounded variation, as opposed to Sobolev functions, can have jump discontinuities. It is also reflected in numerical realizations where TV minimization allows to effectively preserve sharp interfaces. A disadvantage of TV regularization, however, is its tendency to produce piecewise constant results even in non-piecewise-constant regions, which is known as the staircasing effect. Employing a regularization with higher order derivatives, such as second-order TV regularization, overcomes this drawback, but again does not allow for jump discontinuities. As a result, a lot of recent research aims at finding suitable extensions of TV that overcome the staircasing effect, but still allow for jumps [31,75,21]. While infimal-convolutiontype approaches can be seen as the first methods to achieve this, the introduction of the TGV functional (of arbitrary order k) in 2010 finally provided a complete model for piecewise smooth data with jump discontinuities. This is achieved by an optimal balancing between first and higher order derivatives (up to the order k), which is carried out as part of the evaluation of the functional. We refer to [29] for more motivation and to [25] for a detailed analysis of TGV in the context of inverse problems. In particular, second-order TGV, which balances between first and second-order derivatives, achieves piecewise linear -as opposed to piecewise constant -image reconstructions while it still allows for jumps. This renders second-order TGV a well suited model for piecewise smooth images and can be seen as the motivation of the use of second-order TGV regularization in a plethora of applications [57,26,27,58,65]. Up to now, TGV regularization was only available for vector space data and applications are hence limited to this situation.
In various problems of applied sciences, however, the data do not take values in a linear space but in a nonlinear space such as a smooth manifold. Examples of manifold-valued data are circle and sphere-valued data as appearing in interferometric SAR imaging [62], wind directions [74], orientations of flow fields [3,79], and color image processing [33,89,56,60]. Other examples are data taking values in the special orthogonal group SO(3) expressing vehicle headings, aircraft orientations or camera positions [87], Euclidean motion group-valued data [72] as well as shapespace data [63]. Another prominent manifold is the space of positive (definite) matrices endowed with the Fisher-Rao metric [70]. This space is the data space in DTI [69]. It is a Cartan-Hadamard manifold and as such it has particularly nice differential-geometric properties. DTI allows to quantify the diffusional characteristics of a specimen non-invasively [11,54] which is helpful in the context of neurodegenerative pathologies such as schizophrenia [44] or autism [5]. Because of the natural appearance of these nonlinear data spaces quite a lot of recent work deals with them. Examples are wavelet-type multiscale transforms for manifold-valued data [87,51,91,90,92] as well as manifold-valued partial differential equations [85,35,48]. Work on statistics on Riemannian manifolds can be found in [66,23,43,24,68,42]. Optimization problems for manifold-valued data are for example the topic of [2,1], of [49] and of [52] with a view towards learning in manifolds. We also mention related work on optimization in Hadamard spaces [10,9,7] and on the Potts and Mumford-Shah models for manifold-valued data [94,81].
TV functionals for manifold-valued data have been considered from the analytical side in [46,47,45]; in particular, the existence of minimizers of certain TV-type energies has been shown. A convex relaxation based algorithm for TV regularization for S 1 -valued data was considered in [82,37]. Approaches for TV regularization for manifold-valued data are considered in [61] which proposes a reformulation as multi-label optimization problem and a convex relaxation, in [50] which proposes iteratively reweighted minimization, and in [93] which proposes cyclic and parallel proximal point algorithms. An exact solver for the TV problem for circle-valued signals has been proposed in [80]. Furthermore, [15] considers half-quadratic minimization approaches, which may be seen as an extension of [50], and [19] considers an extension of the Douglas-Rachford algorithm for manifold-valued data. Applications of TV regularization to shape spaces may be found in [13,78]. TV regularization with a data term involving an imaging type operator has been considered in [14] where a forward-backward type algorithm is proposed to solve the corresponding inverse problem involving manifold-valued data.
As with vector space data, TV regularization for manifold-valued data has a tendency to produce piecewise constant results in regions where the data is smooth. As an alternative which prevents this, second-order TV type functionals for circle-valued data have been considered in [18,20] and, for general manifolds, in [10]. However, similar to the vector space situation, regularization with second-order TV type functionals tends towards producing solutions which do not preserve the edges as desired. To address this drawback, the vector space situation guides us to considering TGV models and models based on infimal convolution to address this issue. The most difficult part in this respect is to define suitable notions in the manifold setup which have both reasonable analytic properties on the one hand and which are algorithmically realizable on the other hand. Concerning infimal-convolution type functionals such as TV-TV 2 infimal convolutions, a first effort towards a generalization to the manifold setting has been made in the recent conference proceeding [17] which has later been extended in an arXiv preprint [16]. The present manuscript [28], which was submitted to arXiv at the same time as [16], emerged independently. In this work, we focus on TGV and aim at providing a thorough study of reasonable generalizations of TGV for piecewise smooth, manifold-valued data by first investigating crucially required properties of generalizations of TGV to the manifold setting. Then we propose suitable extensions that fulfill these properties and which are, in addition, computationally feasible. In this respect, it is important to note that due to the cascadic nature of TGV (as opposed to infimal-convolution-type functionals), its generalization to manifolds requires a conceptually new approach. For this reason, both the techniques we propose for a generalization of TGV as well as our algorithmic setting are rather different from the one of [17], which uses a mid-point approach to generalize TV-TV 2 infimal convolutions and gradient descent for numerical minimization.

Contributions
The contributions of the paper are as follows: (i) we lay the foundations for -and provide concrete realizations of -a suitable variational model for second-order total generalized variation for manifold-valued data, (ii) we provide algorithms for the proposed models, (iii) we show the potential of the proposed algorithms by applying them to synthetic and real data. Concerning (i), we use an axiomatic approach. We first formalize reasonable fundamental properties of vectorvalued TGV which should be conserved in the manifold setting. Then we propose two concrete realizations which we show to fulfill the axioms. The one is based on parallel transport and the other is motivated by the Schild's approximation of parallel transport. We obtain well-posedness results for TGV-based denoising of manifold valued data for both variants. Concerning (ii) we provide the details for a numerical realization of variational regularization for general manifolds using either of the two proposed generalizations of TGV. We build on the well-established concept of cyclic proximal point algorithms (using the inexact variant). The challenging part and the main contribution consists in the computation of the proximal mappings of the TGV atoms involved. In particular, for the class of symmetric spaces, we obtain closed form expressions of the related derivatives in terms of Jacobi fields for the Schild variant; for the parallel transport variant, we obtain closed form expressions for the related derivatives for general symmetric spaces up to a derivative of a certain parallel transport term for which we can still analytically compute the derivative for the class of manifolds considered in the paper. Concerning (iii), we provide a detailed numerical investigation of the proposed generalization and its fundamental properties. Furthermore, we provide experiments with real and synthetic data and a comparison to existing methods.

Outline of the paper
The paper is organized as follows. The topic of Section 2 is the derivation of suitable TGV models for manifold-valued data. We start out with a detailed discussion of fundamental properties expected by a reasonable TGV version for manifold-valued data. To this end, we first reconsider vector-space TGV in a form suitable for our purposes in Subsection 2.1. Then we derive an axiomatic extension of TGV 2 α to the manifold setting where we -for a better understandingfirst consider the univariate case in Subsection 2.2. Next we suggest two realizations which we call the Schild variant of TGV and the parallel transport variant of TGV and show that these realizations indeed fulfill all desired properties. Then, in Subsection 2.3, we extend the axiomatic framework to the bivariate setting and provide bivariate versions of the Schild variant of TGV and the parallel transport variant of TGV, respectively, that we show to fulfill the required axioms. In Section 3 we then provide existence and lower semi-continuity results for the proposed variants of TGV from which the existence of minimizers for the TGV-regularized denoising problems of manifold-valued data follows. The topic of Section 4 is the algorithmic realization of the proposed models. We start out recalling the concept of a cyclic proximal point algorithm (CPPA). Then, in Subsection 4.1, we consider the implementation of the CPPA for manifold-valued TGV. We identify certain TGV atoms whose proximal mappings are challenging to compute. The necessary derivations needed for their computation are provided for the Schild variant of TGV in Subsection 4.2, and, for the parallel transport variant of TGV, in Subsection 4.3. The topic of Section 5 is the numerical evaluation of the proposed schemes. There, we first carry out a detailed numerical evaluation of the proposed model for synthetic data. We test extreme parameter choices and ensure consistency of the results of our numerical method with a reference implementation for vector spaces. Then we consider various application scenarios and compare to existing methods. Finally, we draw conclusions in Section 6.

Basic Notation and Concepts from Differential Geometry
Throughout this work, M ⊂ R K will always denote a complete manifold with a Riemannian structure (with its canonical metric connection, its Levi-Civita connection). Let us explain these notions more precisely. We consider a manifold M with a Riemannian metric which is a smoothly varying inner product ·, · p in the tangent space of each point p. (We note that any paracompact manifold can be equipped with a Riemannian structure.) As usual, we will frequently omit the dependence on p in the notation in the following. According to the Hopf-Rinow theorem, the complete manifold M is geodesically complete in the sense that any geodesic can be prolongated arbitrarily. Here a geodesic is a curve γ of zero acceleration, i.e., D dt d dt γ = 0, where D dt denotes the covariant derivative along the curve γ (see below for details on the covariant derivative). Geodesics are invariant under affine reparametrizations and usually identified with their image in M. (We parametrize them on [0, 1] or, more generally, on [t 0 , t 1 ] with t 0 < t 1 , depending on the context.) We always denote by d : M × M → [0, ∞) the distance on M which is induced by the Riemannian structure and note that, since M is complete, for any a, b ∈ M there is always a geodesic from a to b whose length equals d(a, b). Such geodesics are called length minimizing geodesics between a and b. By T a M we denote the tangent space of M at a ∈ M and by T M the tangent bundle of M. Further, we need the notion of parallel transport. The parallel transport of a vector v ∈ T a M with a ∈ M to a point b ∈ M along a curve γ : where the covariant derivative D dt is taken along the curve γ. The covariant derivative D dt is induced by the Levi-Civita connection ∇ X Y of the Riemannian manifold M; a connection ∇ X Y is a C ∞ -linear mapping w.r.t. the vector field X and a derivation w.r.t. the vector field Y. (Connections are needed since, in a general manifold M, there is no a priori canonical way to define the directional derivative in direction X of a vector field Y .) For the covariant derivative D dt Y along a curve γ, a direction X along the curve γ is given by d dt γ and the vector field Y is given along the curve; hence, if d dt γ = 0, D dt Y is just given by ∇ d dt γ Y. In a Riemannian manifold M, there is a (uniquely determined) canonical metric connection, its Levi-Civita connection. The Levi-Civita connection is the only connection which is symmetric (i.e.,∇ X Y = ∇ Y X for any vector fields X, Y which commute) and which is compatible with the metric (i.e., in terms of the induced covariant derivative which we need later on, d dt X, Y = D dt X, Y + X, D dt Y for any two vector fields X, Y along a curve γ, along which the covariant derivatives are taken.) For an account on Riemannian geometry we refer to the books [76,38] or to [59].
For later use, we fix some further notation. To this end, let a, b ∈ M, v ∈ T a M, and γ : [0, L] → M be a curve connecting a and b such that γ(t 0 ) = a, γ( 4. in case γ is a geodesic and length-minimizing between a and b on [t 0 , t 1 ], we denote log γ a (b) = (t 1 − t 0 ) d dt γ(t 0 ), i.e., the vector in log a (b) that is parallel to d dt γ(t 0 ) and such that exp a (log γ a (b)) = b, where z is the vector resulting from the parallel transport of v from a to b along the curve ψ that reparametrizes [t 0 , t 1 ] to [0, 1] in an affine way, i.e., ψ is a length-minimizing geodesic connecting a and b}, We extend this definition for t ∈ R \ [0, 1] by extending the corresponding geodesic.
We also note that throughout the paper we identify sets having only one element with the corresponding element.
We further need the concept of geodesic variations for the existence results in Section 3 and the derivation of the algorithms in Section 4. A variation of a curve γ : I → M defined on an interval I is a smooth mapping V : I × J → M, J an interval containing 0, such that V (s, 0) = γ(s) for all s ∈ I. A variation is called geodesic, if all curves s → V (s, t) are geodesics for any t in J.

Definition of TGV for manifold-valued data
The goal of this section is to define a discrete total generalized variation (TGV) functional of second order for manifold-valued data. To this end, we first state some fundamental properties of the TGV functional in infinite and finite dimensional vector spaces. The definition of a generalization of TGV to the manifold setting will then be driven by the goal of preserving these fundamental properties of vector-space TGV.

TGV on vector spaces
We first recall the definition of TGV on infinite-dimensional vector spaces via its minimum representation which is, according to results in [30,88,25,26] covering the second-and general order case as well as the scalar and vectorial case, equivalent to the original definition as provided in [29]. Then we present a discretization which is slightly different from the one typically used.

Definition 2.1 (Minimum representation of TGV in vector spaces
Here , · M denotes the Radon norm in the space of Radon measures M(Ω, X) K := (C 0 (Ω, X) K ) * with X ∈ {R d , S d×d }, S d×d the space of symmetric matrices and further BD(Ω, The derivatives Du ∈ D(Ω, R d ) K and Ew ∈ D(Ω, S d×d ) K are defined in the weak sense by See [25,26] for details.
Note that TGV 2 α (u) is finite if and only if u ∈ BV(Ω) K and, in this case, the minimum is actually obtained [25]. Hence the term min in the above definition is justified.
For a discretization and generalization later on, it is convenient to list some of the main properties of second-order TGV in function space (see [25,26]): α (u) = 0 if and only if u is affine. Here, TV 2 denotes a second-order TV functional which can be defined as TV 2 (u) = D 2 u M where D 2 u denotes the second-order distributional derivative of u and the above quantities are finite if and only if D 2 u ∈ M(Ω, S d×d ) K . We note that w = 0 and w = Du can trivially be obtained when u is constant and affine, respectively, but w = 0 is for instance also obtained under some symmetry conditions when u is piecewise constant or when the ratio α 0 /α 1 is large enough [29,67].
Using the minimum representation above, a discretization of TGV 2 α in vector spaces on twodimensionals grids can be given as follows.
Definition 2.2 (Discrete isotropic and anisotropic TGV in vector spaces). Set U = R K . For u = (u i,j ) i,j ∈ U N ×M with u i,j ∈ U , we define the discrete second-order TGV functional as where and with δ x+ , δ y+ and δ x− , δ y− being standard forward and backward differences, respectively. Introducing a parameter p ∈ [1, ∞), the one-norms are given as and where we extend the signals by zero to have the same size. Here, |(w 1 , w 2 )| p := |w 1 | p + |w 2 | p 1/p and Note that we incorporate a pointwise p -norm with p ∈ [1, ∞) in the definition of TGV 2 α . In the function space setting of Definition 2.1, this corresponds to choosing the appropriate dual norm on X ∈ {R d , S d×d }. There, any choice of point-wise norm (also beyond p -norms) yields equivalent functionals and function spaces, with the Euclidean norm being used in the first paper on TGV [29] and the most popular choices being p norms with p ∈ {1, 2}. The focus of this paper is the case p = 1 since this is more frequently used in the manifold context because of being more amenable to numerical realization. In the theory part we include the case of p ∈ [1, ∞) since it can be done without additional effort.
We remark that the above discretization of TGV 2 α is slightly different from the standard one as provided, e.g., in [29]. The purpose of this is to achieve consistency of the zero set of both the continuous and discrete version as follows. Also, note that, notation-wise, we do not distinguish between the continuous and the discrete version of TGV (and of TV and TV 2 ). In the following, we will always refer to discrete versions.
Proof. We provide only the basic steps: Setting w = Du we get that Now it is easy to see that if u is linear as above, EDu = 0 and hence the choice w = Du renders TGV 2 α (u) to be zero. Conversely, TGV 2 α (u) = 0 implies w = Du and hence EDu = 0. It is easy to see that δ x− δ x+ u = δ y− δ y+ u = 0 implies that u is of the form u i,j = rij + ai + bj + c with r, a, b, c ∈ U . But in this case, the last component of EDu being zero implies that r = δ y− δ x+ u = −δ x− δ y+ u = −r, hence r = 0.
The latter proposition shows that, also after discretization, the kernel of TGV 2 α consists exactly of (discrete) affine functions. In fact, this is one of the fundamental properties of TGV (corresponding to (P3)) which should also be transferred to an appropriate generalization of TGV for manifolds. Regarding appropriate counterparts of (P1) and (P2), we introduce the following definition.
In summary, using the discretized version of TV, TV 2 and TGV 2 α and the discrete notion of affine functions as in Proposition 2.3, we can conclude that the properties (P1) to (P3) of TGV in continuous vector spaces also transfer to its discretization. Figure 1: A three-point section of a signal u together with tangent vectors w represented by the endpoints y. The red vectors show the tangent vectors at the signal points, the blue dotted lines indicate the geodesics t → exp u (tw) and the y are the endpoints exp u (w). The black lines indicate geodesic interpolations between the signal points and are for visualization purposes only.

Univariate TGV 2 α on manifolds
When moving from the vector space to the manifold setting, a main difference is that vectors representing derivatives can no longer be made independent of their location, but are attached to a base point on the manifold. In other words, in the Euclidean setting, all tangent spaces at all locations can be identified, which is not possible for manifolds in general. Accordingly, when aiming to extend TGV to the manifold setting, important questions are how to introduce the vector fields (w i,j ) i,j appearing in the minimum representation, how to define Ew for such elements, and, most importantly, how to define a suitable distance-type function between such tangent vectors sitting in different points.
In the following, we describe our main ideas to resolve these questions. In order to allow the reader to understand the underlying ideas more easily, we consider the univariate setting first.
Our first goal is to suitably extend the notion of forward differences (δ x+ u) i = u i+1 − u i and introduce auxiliary variables w i which they can be compared to. To this end, the central idea is to identify tangent vectors w i in the space T ui M with point-tuples, i.e., w i [u i , y i ] with y i ∈ M. In the vector-space case, this can be done via the correspondence w i = y i − u i . For manifolds, the correspondence can be established via the exponential and the logarithmic map. That is, any w i ∈ T ui M can be assigned to a unique point-tuple [u i , y i ] such that exp ui (w i ) = y i . Conversely, any point-tuple [u i , y i ] can be assigned to (generally multiple) tangent vectors w i such that w i ∈ log ui (y i ). (Note that the ambiguities in logarithmic map result from non-uniqueness of distance-minimizing geodesics, which is a rather degenerate situation in the sense that it only occurs for a set of points with measure zero [53,34].) Figure 1 visualizes the correspondence between tangent vectors w and point-tuples [u, y].
Exploiting this correspondence, our approach is to work with a discrete "tangent bundle" that is defined as the set of point-tuples [u, y] with u, y ∈ M, rather than with the continuous one. We refer to the elements in this discrete "tangent bundle" as tangent tuples. The identification of forward differences with a tangent tuple is then naturally given via (δ The discretized vector-space version of second-order TGV (for which there is a one-to-one correspondence of tangent vectors and point tuples) can then in the univariate case be rewritten as . Note that, in this context, (w i ) i , (y i ) i always denote finite sequences of points with their length being the same as the one of δ x+ u and, whenever out of bound indices are accessed in such a summation, we assume the signals to be extended such that the evaluation of D yields zero cost. Hence, in order to extend TGV to the manifold-setting, we are left to appropriately measure the distance of two tangent tuples, i.e., generalize the expression D([x, y], [u, v] In case both tuples have the same base-point on the manifold, a simple and rather natural generalization is to measure their distance by the distance of their endpoints, that is, we set D([x, y], [x, v]) = d(y, v). In the other case, i.e., for evaluating D([u i , y i ], [u i−1 , y i−1 ]) for u i = u i−1 , there is no such simple or even unique generalization since the two tangent tuples belong to different (tangential) spaces. A quite direct approach to overcome this is to first employ the concept of parallel transport on manifolds to shift both tangent tuples to the same location, and then to measure the distance of their endpoints. One possible generalization for TGV proposed in this paper builds on this idea. This is not the only possible way of measuring the distance of two tangent tuples. Alternatively, one can for instance build on a discrete approximation of parallel transport, i.e., on a variant of the Schild approximation, to measure distance of point tuples at different locations. This second approach is simpler to realize numerically, and it can be build on simpler differential-geometric tools, while the first one is more straight-forward. Both approaches have their motivation and, given this ambiguity in possible generalizations, we take a more systematic approach to the topic. That is, we first formulate fundamental properties that we require from reasonable generalizations of TGV to the manifold setting. Then we translate these properties to requirements on admissible distance-type functions for measuring the distance of two tangent tuples. Having established the later, we propose two possible concrete realizations that fulfill the required properties. Axiomatic extension. Assuming D : M 2 × M 2 → [0, ∞) to be a distance-type function for tangent tuples, we define a generalization of second-order TGV for univariate manifold-valued data as M-TGV 2 α (u) = min In this context, we note that, while with TGV for discrete and continuous vector spaces it is known that the minimum in the above expression is attained, this is not clear a-priori for our generalization. In order not to lose focus, we shift the discussion on existence to Section 3, noting at this point that for all versions of M-TGV 2 α proposed in this work, existence can be shown. Accounting for the fact that the tangent tuples represent vectors in the tangent space, some basic identifications and requirements for the distance-type function D follow naturally. Zero elements in the discrete "tangent bundle" for instance correspond to tangent tuples of the form [u, u] with u ∈ M and the distance of [x, x] and [y, y] should be zero. Also, the distance of two identical elements should also be zero, i.e., D([x, y], [x, y]) = 0, and, in the vector space case, the distance function should reduce to the difference of the corresponding tangent vectors.
Our goal is now to further restrict possible choices of distance-type function in order to obtain a meaningful generalization of TGV. To this aim, we want to ensure that appropriate counterparts of the properties (P1) to (P3) are preserved for the resulting version of M-TGV 2 α also in the manifold setting. In order to make this more precise, we first have to generalize the involved concepts. We start with the notion of "affine".
Given that the geodesics in a manifold play the role of straight lines in vector space, a natural generalization for u = (u i ) i a finite sequence in M to be affine is to require that all points (u i ) i are on a single geodesic at equal distance. A difficulty that arises in connection with this definition is that, in general, as opposed to the vector space case, a geodesic connecting two points is not necessarily unique. As a consequence, even though all points might be at the same distance when following a joint geodesic, the distance of each single pair of points in the manifold is not necessarily equal. To account for that, we require in addition that the geodesic connecting all points also realizes the shortest connection between all involved points on the geodesic locally.
Definition 2.5. Let u = (u i ) n i=0 be a finite sequence of points on a manifold M. We say that u is geodesic if there exists a geodesic γ : [0, L] → M parametrized with unit speed such that γ(iL/n) = u i for i = 0, . . . , n and d(u i , u i+1 ) = L/n for i = 0, . . . , n − 1, i.e., γ| [iL/n,(i+1)L/n] is a geodesic of minimal length connecting u i and u i+1 for all i.
A second issue arising from non-uniqueness of geodesics is the fact that, even though every triplet of points in a signal (u i ) i might be connected by a geodesic, we cannot conclude that all points are connected by a unique geodesic. Consequently, as a reasonable generalization of TGV will typically act local, in particular will be based on three point stencils, we cannot hope to obtain more than a local assertion for signals in the zero set of TGV. To account for that, we introduce the following notion.
Definition 2.6. Let u = (u i ) n i=0 be a finite sequence of points on a manifold M. We say that u is locally geodesic if for each i ∈ {2, . . . , n − 1}, the sequence (u j ) i+1 j=i−1 is geodesic. As one might expect, in the situation that subsequent points of a signal are sufficiently close such that they admit a unique connecting geodesic we obtain equivalence of the notions of geodesic and locally geodesic signals. In this respect, we again note that, in the general non-local situation, the existence of unique minimizing geodesics is true for any set of points outside a set of measure zero [53,34]. Lemma 2.7. If a sequence of points u = (u i ) i in M is such that the length minimizing geodesic connecting each pair u i , u i with |i − i | ≤ 1 is unique, then u is locally geodesic if and only if it is geodesic.
Proof. If u is geodesic, it is obviously locally geodesic. Now assume that u is locally geodesic and that the length minimizing geodesic connecting each pair u i , u i with |i − i | ≤ 1 is unique. Then there exists a geodesic connecting u 0 , u 1 , u 2 at equal distance, i.e., the distance between two subsequent points equals the length of the geodesic segment. We proceed recursively: Now assume that u i−2 , u i−1 , u i are connected by a geodesic at equal distance. As u is locally geodesic also u i−1 , u i , u i+1 are connected by a geodesic. Now by uniqueness, the two geodesics between u i and u i−1 must coincide. Hence all points u i−2 till u i+1 are on the same geodesic at equal speed. Proceeding iteratively, the result follows.
In order to investigate the counterparts of properties (P1) and (P2), we need to define generalizations of the TV and TV 2 functionals to the manifold setting. For TV, a natural generalization is to set, for u = (u i ) i in M, see also [93]. For TV 2 , a generalization is not that immediate. In [10], second-order TV was essentially (and up to a constant) defined as TV 2 (u) = 2 i inf c∈[ui−1,ui+1] 1 2 d(c, u i ), which, in the vector space setting reduces to However, as can be deduced by considering TV 2 as special case of the different versions of TGV proposed below, this is not the only generalization which fulfills such a property. We will call a function an admissible generalization of TV 2 whenever, in the vector space setting, it reduces to TV 2 as in Equation (5) above. Using these prerequisites, we now formulate our requirements on an appropriate generalization of TGV to the manifold setting of the form (3).
(M-P1) In the vector space setting, M-TGV 2 α reduces to the univariate version of (2).
(M-P2) If the minimum in (3) is attained at (y i ) i = (u i ) i , i.e., the tangent tuples [u i , y i ] all correspond to zero vectors, then M-TGV 2 α (u) = α 1 TV(u) with TV as in (4).
(M-P4) M-TGV 2 α (u) = 0 if and only if u is locally geodesic according to Definition 2.6. As the following Proposition shows, the Properties (M-P1) to (M-P3) already follow from basic requirements on the distance function that arise from the intuition that tangent tuples represent tangent vectors. The most restrictive property is (M-P4), which requires an additional assumption on the distance function.
Proposition 2.8. Assume that the function D : Then, for M-TGV 2 α as in (3), the properties (M-P1) to (M-P3) hold. If we further assume that for any geodesic three-point signal j=i−1 is geodesic. Then, for any signal Proof. In the vector space case we get by our assumptions and since d(x, u) = |x − u|, that which coincides with the univariate version of TGV as in (2), hence (M-P1) holds. Now in case the minimum in (3) is achieved for (y i ) i = (u i ) i we get that and which is an admissible generalization of TV 2 since, by assumption, from above with the right-hand side in (7) and obtain that M-TGV 2 α (u) = 0. Conversely, in case M-TGV 2 α (u) = 0 we get that necessarily (y i ) i = (u i+1 ) i and M-TGV 2 α reduces to the right hand side in (7). By the assumption on D, this implies that u is locally geodesic and by Lemma 2.7 the result follows. Guided by the aim of fulfilling the requirements of Proposition 2.8, in the following, we propose two possible choices for D(·, ·) which result in two different concrete versions of M-TGV 2 α . Both variants follow the general idea of first transporting different tangent tuples to the same location in M and then measuring the distance there. The main difference between the two variants will be the way the transport is carried out. The first concept is motivated by the so called Schild's ladder approximation of parallel transport. Its implementation requires only basic differential geometric concepts and is therefore presented first. The second realization, which we call the parallel transport variant, requires more differential geometric concepts [76,38] and is therefore presented afterwards. However, the Schild's ladder variant can be seen as an approximation of the parallel transport variant as explained below. [u, v]) and assume for the moment that distance-minimizing geodesics are unique. Motivated by the Schild's ladder approximation of parallel transport [39,55], we consider the following construction (see Figure 2): take c to be the midpoint of the points v and , reflect u at c. Then, we claim that the distance-type function D([x, y], [u, v]) = d(y, y ) fulfills the requirements of Proposition 2.8.
The above construction may be motivated as follows. To compare the tangent vectors log x (y) and log u (v), which correspond to the tangent tuples [x, y] and [u, v], we could move them to the same point using parallel transport and then take the norm induced by the Riemannian metric in the corresponding tangent space. Using the above construction to obtain y from x and [u, v], the tangent vector log x (y ) can be seen as an approximation of the parallel transport of log u (v) to x (in the limit x → u) [39,55]. The difference of this vector to log x (y) is then approximated by d(y, y ). We notice that, as can be easily seen, this Schild's ladder approximation of parallel transport, i.e., the construction of y and log x (y ) as above, is exact in the vector space case.
Hence we propose to measure the deviation between the tangent tuples [x, y] and [u, v] by transporting [u, v] to x using the construction above to obtain the tangent tuple [x, y ] and then to compare the resulting tangent tuples sitting in the same point x by measuring the distance of their endpoints y, y . Since in general geodesics are not unique, we have to minimize over all possible constructions as above, i.e., consider all midpoints of all length-minimizing geodesics. This yields the following distance-type function It is immediate that D S is positive and zero for identical elements. Plugging Regarding the properties (M-P1) to (M-P4), we then get the following result.
Theorem 2.9. The S-TGV 2 α functional as in (9) satisfies (M-P1) and (M-P3). If a finite sequence of points u = (u i ) i is locally geodesic, then S-TGV 2 α (u) = 0. Further, if S-TGV 2 α (u) = 0 and the geodesic connecting each pair u i , u i with |i − i | ≤ 1 is unique, then u is geodesic.
Proof. We verify the assumptions of Proposition 2.8 for D = D S . It can be easily seen that Also, in the case M = R K , the approximate parallel transport as in the definition of D S above coincides with the parallel shift of vectors. 2 and, by the assumption on unique geodesics, it follows that v is geodesic.
Remark 2.10. Assuming that, for each i, u i and y i are sufficiently close such that they are connected by a unique length minimizing geodesic, we get that Consequently, in this situation, an equivalent definition of S-TGV 2 α can be given as This definition regards the mapping (u i ) i → ([u i , u i+1 ]) i as discrete gradient operator, mapping from M to the discrete tangent space, and exclusively works in the discrete tangent space with D S the canonical distance-type function.
Remark 2.11. We note that an alternative choice for D S would have been to transport [u i , y i ] to u i−1 rather than [u i−1 , y i−1 ] to u i . However, since the distance d(u i+1 , y i ) in S-TGV 2 α as in (9) can be interpreted as evaluating the difference of the forward difference (δ x+ u) i with [u i , y i ] in the center point u i , it seems natural to evaluate also the backward difference of the signal ([u i , y i ]) i at the center point.
Remark 2.12. We also note that, as opposed to the vector space setting, the distance function D S is in general not symmetric, i.e., . To obtain symmetry, alternative definitions of D S could be given asD For the sake of simplicity and in order to obtain the relation with parallel transport, however, we have defined D S as in (8).
Realization via parallel transport. The Schild's ladder construction defined above can be seen as a discrete approximation of parallel transport. Alternatively, we can use the identification of point-tuples [u, v], [x, y] with vectors in the tangent space at u and x, respectively, and use the parallel transport directly to transport the respective vectors to a common base point.
That is, assuming -for the moment -uniqueness of length-minimizing geodesics, we can identify each [u, v] with w ∈ T u M such that v = log u (w) and compare pt x (w), the vector resulting from the parallel transport of w to T x M, to log x (y) ∈ T x M. This yields a distancetype function for two point-tuples [u, v], [x, y] as where · x denotes the norm in T x M. Comparing this to D S we note that, besides using a different notion of transport, we now measure the distance with the norm in the tangent space rather than with the distance of endpoints. The reason for doing so is that, in this situation, approximating the norm in the tangent space via the distance of endpoints does not yield further simplification since the the parallel transport forces us to work in the tangent space anyway. Now in general, length minimizing geodesics are not necessarily unique. To deal with this issue, we have defined the log mapping and the parallel transport to be set-valued (see Section 1.3), i.e., for u, v ∈ M, not necessarily close, and w ∈ T v M we define pt u (w) ⊂ T u M to be the set of all vectors in T u M which can be obtained by parallel-transporting w to T u M along a length minimizing geodesic. Note that by isometry of the parallel transport, the length of all such vectors is the same but by holonomy their orientation might be different. In order to adapt the above definition of D pt to this situation, we generalize the distance function D pt for tangent Using this notation, a parallel-transport-based version of M-TGV 2 α , denoted by PT-TGV 2 α , can then be defined for the general situation as Similarly to the Schild's-version, we then get the following result. Proof. The proof is similar to the one of Theorem 2.19 below considering the bivariate setting.
Remark 2.14. Since the expression d(u i+1 , y i ) can be seen as approximation of | log ui (u i+1 ) − w i | ui , where w i ∈ log ui (y i ), an alternative definition of PT-TGV 2 α (assuming uniqueness of geodesics for simplicity) only in terms of tangent vectors is given as We believe, however, that the originally proposed version is preferable since the fact that the first term in PT-TGV 2 α only involves the standard manifold distance of two points simplifies the numerical realization. That is, with the standard distance we are able to solve certain subproblems (proximal mappings) in the algorithm explicitly whereas the other version would require additional inner loops (see Section 4 for details).

Bivariate TGV 2 α on manifolds
The goal of this section is to extend M-TGV 2 α to the bivariate case. To this aim, we first observe that the bivariate version of TGV 2 α for vector spaces as in (2) can be written as where we set all norms | · | p to be zero whenever out of bound indices are involved. The first four summands in the above definition of TGV 2 α can all be transferred to the manifold setting using the previously introduced distance-type functions for tangent tuples. The additional difficulty of the bivariate situation arises form the fifth term, which first combines both differences and sums of tangent tuples and then measures the norm. Again possible generalizations are not unique and we employ an axiomatic approach to propose reasonable choices. Axiomatic extension. Denote by D(·, ·) one of the two previously introduced distance functions for tangent tuples and assume that D sym : M 2 × M 2 × M 2 × M 2 → R generalizes the fifth term in (14), which corresponds to the mixed derivatives. A bivariate version of M-TGV 2 α can then be given as The basis for this generalization is again the representation of tangent vectors with tangent tuples, only that now for each u i,j we consider two tangent vectors w 1 i,j , w 2 i,j and corresponding points y 1 i,j , y 2 i,j in order to represent horizontal and vertical derivatives, see Figure 3. Remember that in this paper, we focus on the case p = 1, for which we derive an algorithmic realization later on. However, since it causes no additional effort in this section, we provide a formulation for all p ∈ [1, ∞) here.
We now extend the requirements (M-P1) to (M-P4) to the bivariate case. The generalization of TV to bivariate manifold-valued data is quite straightforward [93]. For u = (u i,j ) i,j we define A generalization of second-order TV is again less straightforward and we call any functional TV 2 acting on u = (u i,j ) i,j an admissible generalization of TV 2 if it reduces, in the vector space setting, to TV 2 as given in Definition 2.4. Note that our version of TV 2 differs from the definition of TV 2 as given in [10] since we use a symmetrization of the mixed derivatives.
A generalization of affine functions for the bivariate manifold setting is given as follows.
Definition 2.15. Let u = (u i,j ) i,j be a finite sequence of points on a manifold M. We say that u is (locally) geodesic if, for each (i 0 , j 0 ), the univariate signals Note that this indeed generalizes the notion of affine for vector space data. While the first two conditions ensure that u is of the form u i,j = aij + bi + cj + d, the last condition ensures that a = 0, i.e., no mixed terms occur. In the vector space case the third condition is equivalent to requiring that (u i0+k,j0+k ) k is (locally) geodesic. In the manifold setting, this is not equivalent in general and hence our definition of geodesic is somewhat anisotropic. The reason for using the anti-diagonal ((i 0 + k, j 0 − k)) and not the diagonal ((i 0 + k, j 0 + k)) direction is that the former arises naturally from the usage of forward-backward differences, as will become clear in the discussion after Proposition 2.17. An alternative would also be to require both (u i0+k,j0+k ) k and (u i0+k,j0−k ) k to be geodesic. This, however, seems rather restrictive and for general manifolds such signals might even not exist. Hence we do not impose this additional restriction.
As a direct consequence of the corresponding result in the univariate setting, we obtain equivalence of the notion of locally geodesic and geodesic if neighboring points are sufficiently close. Lemma 2.16. Let u = (u i,j ) i,j be a finite sequence in M. If the length minimizing geodesic connecting any two points u i,j and u i ,j with max{|i − i |, |j − j |} ≤ 1 is unique, then u is locally geodesic if and only if it is geodesic.
With these prerequisites, we now state our requirements for a reasonable generalization of TGV in the bivariate case.
Then, for M-TGV 2 α as in (3), for any geodesic three-by-three signal v = (v i,j ) i,j , then M-TGV 2 α (u) = 0 for any locally geodesic u. Conversely, assume that for any three-by-three k=−1 being geodesic. Then, for any u = (u i,j ) i,j such that the geodesic connecting each pair u i,j and u i ,j with max{|i − i |, |j − j |} ≤ 2 is unique, we get that M-TGV 2 α (u) = 0 implies u being geodesic.
Proof. See Section A.1 in the Appendix.
Similar to the univariate case, the most restrictive requirement for D sym is that, in case of uniqueness of length-minimizing geodesics, for (u i,j ) i,j being locally geodesic in horizontal and vertical direction, it holds that Let us discuss, again temporarily assuming uniqueness of geodesics, a general strategy to design the function D sym such that this property holds. We consider the situation around a point u •,• and denote by w 1 at which D sym is evaluated, see Figure 3. The evaluation of D sym at these points should, in the vector space case, correspond to one of the equivalent formulations Disregarding for the moment the fact that the corresponding tangent tuples live on different locations, the difficulty here is how to define algebraic operations, i.e., sum or difference operations, on two tangent tuples. For the sum, there is no direct way of doing so, the difference operation, however, quite naturally transfers to tangent tuples as Hence, in the situation that all tangent tuples live on the same location, we can define D sym by first evaluating pair-wise differences of the four point tuples and then measuring the distance of the two resulting tuples by D(·, ·). Now in the general situation, we first need to transport tangent tuples to the same location. As described in Section 2.2, we have means of doing this that build either on parallel transport or its Schild's approximation. Taking this into account, we note that the evaluation of D sym in the vector space case can be written purely in terms of distances of tangent-vector-differences in various ways, e.g., While any order of taking the differences is equivalent in the vector space case, in the manifold setting the order defines how the corresponding tangent tuples need to be transported to the same location. Looking again at Figure  The drawback of this simple solution can be found when looking at the conditions of Proposition 2.17. There, the main difficulty is to ensure that allows to conclude that the signal is geodesic in the anti-diagonal direction. The situation that is relevant for this condition is when both (u i,j0 ) i and (u i0,j ) j are locally geodesic for each (i 0 , j 0 ) and the endpoints y coincide with the respective signal points, i.e., y 1 i,j = u i+1,j and y 2 i,j = u i,j+1 (see Figure 4, left). In this situation, in order to fulfill the above condition, we need to know something about the transported tangent tuples, e.g., • . Due to holonomy however, even if the transport is carried out along a geodesic, there is, to the best of our knowledge, no way of obtaining such a result. That is, the vector corresponding to the transported tangent tuple [u •,• , x] might be arbitrarily rotated such that x is far away from u +,• .
On the other hand, a well-known fact for manifolds is that the parallel transport of a vector that is tangential to a geodesic along the geodesic itself again results in a vector that is tangential to the geodesic. As we will see, this implies an equivalent assertion for both our variants for transporting tangent tuples. That is, the transport of a tangent tuple corresponding to a tangential vector of a geodesic along the geodesic itself again results in a tangent tuple that corresponds to a tangential vector of the geodesic. In the above-described particular situation (see again Figure 4, left), this means for example that, since the points u •,+ , u •,• , u •,− are on a geodesic and w 2

•,•
[u •,• , u •,+ ] corresponds to a tangent vector that is tangential to this geodesic, we know that the transport of [u •,• , u •,+ ] to u •,− will result in the tangent tuple [u •,− , u •,• ]. More generally, for the situation that both (u i,j0 ) i and (u i0,j ) j are locally geodesic for each (i 0 , j 0 ), it means that we can always transport the tangent tuples corresponding to w 1 •,• and w 1 •,− in horizontal direction and the tangent tuples corresponding to the w 2 •,• and w 2 −,• in vertical direction and still know something about the transported tangent tuples.
In view of (M-P4'), a natural approach to define D sym in the general situation is hence to restrict ourselves to these particular transport directions. Of course, we also want to define D sym in an as-simple-as-possible way, meaning that we want to carry out as few transport operations as possible. A quick case study in Figure 3, shows that we have to transport at least three times and that, accounting for the facts that we need to generalize the expression (18) and that we can take differences but not sums of tangent tuples in the same location, there is only one possibility to achieve this. As highlighted in Figure 4 on the right for the general situation, our proposed approach is to transport the tangent tuple corresponding to w 2 and take the differences, which results in the tangent tuples [y 1 , which again requires one transport operation. This corresponds to the reformulation of In the particular situation of Figure 4, left, due to the usage of the particular transport directions, the two tangent tuples resulting from the transport of the tangent tuples corresponding to w 1 The above-described strategy is now used to extend both S-TGV 2 α and PT-TGV 2 α to the bivariate setting. As such, the only difference will be how to carry out the transport operations and which of the two functions D S , D pt is used to measure the distance of tangent tuples. Realization via Schild's-approximation. We now realize the axiomatic setting for bivariate TGV for manifold-valued data of Section 2.3 using the proposed Schild's approximation. In particular, we use the Schild's approximation of parallel transport to obtain an instance of the distance-type function D sym .
We define the following bivariate version S-TGV 2 α of bivariate TGV for manifold-valued data using the Schild's approximation by S-TGV 2 α (u) = min Here D S is given as in Equation (8) We note that, for the sake of a formally correct definition, we have introducedũ •,• , but in fact we will always chooseũ •,• = u •,• since D sym S effectively depends only on seven variables. Also note that D sym S exactly carries out the construction of D sym as described in Section 2.3, using the Schild's approximation as shown in Figure 2 for transporting tangent tuples and D S to measure their difference. In the notation of Figure 4  Realization via parallel transport. Using parallel transport in order to realize a symmetrized gradient as described in Section 2.3 yields a parallel transport version PT-TGV 2 α of TGV for bivariate manifold-valued data case as follows. We let PT-TGV 2 α (u) = min Here D pt is given as in (11) s.t. r 1 ∈ exp(pt u−,• (w 1 )) with w 1 ∈ log u•,• (y 1 •,• ) and r 2 ∈ exp(pt u•,− (w 2 )) with w 2 ∈ logũ •,• (y 2 •,• ).
Note that, as for D sym S , for the sake of a formal correctness, we have introducedũ •,• , but in fact we will always chooseũ •,• = u •,• since also D sym pt effectively depends only on seven variables. Also, we note again that D sym pt realizes exactly the construction of D sym as discussed in Section 2.3 (see also Figure 4), only that now the parallel transport is used to transport tangent tuples and the distance of two tangent-tuple-differences is again measured with D pt . Due to our careful construction, the following result follows easily. If u is such that the geodesic connecting each pair u i,j and u i ,j with max{|i−i |, |j −j |} ≤ 2 is unique and PT-TGV 2 α (u) = 0, then u is geodesic.
Proof. See Section A.1 in the Appendix.

Existence results for TGV for manifold-valued data
This section provides existence results for the minimization problem appearing the definition of M-TGV 2 α as well as for variational M-TGV 2 α -regularized denoising. For the sake of brevity, we only consider the bivariate setting and note that the univariate counterpart follows as special case. We first state the main existence results in a general setting, that builds on lower semi-continuity of the involved distance-type function. The latter is then proven in subsequent lemmata.
Also, all previously introduced distance-type functions D S , D sym S , D pt , D sym pt are lower semicontinuous and a solution to (23) exists if M-TGV 2 α is replaced by S-TGV 2 α or PT-TGV 2 α . In particular, all the models proposed in this paper have a minimizer.
Proof. The existence results in the general setting follow by standard arguments and the Hopf-Rinow theorem and we only provide a sketch of proof. Regarding existence for the evaluation of M-TGV 2 α , we note that all involved terms are positive and, for fixed u, any infimizing sequence ((y 1 ) n , (y 2 ) n ) n is bounded due to the first two distance functions involved in M-TGV 2 α as in Equation (15). Hence it admits a convergent subsequence and from lower semi-continuity of all the involved terms, it follows that the limit is a minimizer for the evaluation of M-TGV 2 α . Now take any sequence (u n ) n converging to some u for which, without loss of generality, we assume that lim inf n M-TGV 2 α (u n ) = lim n M-TGV 2 α (u n ). We can pick elements ((y 1 ) n , (y 2 ) n ) for which the minimum in M-TGV 2 α (u n ) is attained. Again by boundedness of (u n ) n and the first two distance functions in the definition of M-TGV 2 α , also ((y 1 ) n , (y 2 ) n ) n is bounded and we can extract a subsequence ((y 1 ) n k , (y 2 ) n k ) k converging to some (ŷ 1 ,ŷ 2 ). Defining E(·) such that M-TGV 2 α (u) = inf y 1 ,y 2 E(u, y 1 , y 2 ), we obtain from lower semi-continuity of all involved terms that and since lim inf k M-TGV 2 α (u n k ) = lim inf n M-TGV 2 α (u n ), M-TGV 2 α is lower semi-continuous. Given lower semi-continuity of M-TGV 2 α , existence for (23) follows by similar arguments. Regarding the particular realizations S-TGV 2 α and PT-TGV 2 α , it suffices to show lower semicontinuity of the involved distance-type functions, which is done on the sequel.
We are now left to show lower semi-continuity of the proposed distance-type functions D S , D pt , D sym S and D sym pt . To this aim, we exploit that they are all defined by minimizing a distancetype function subject to a constraint. Using this joint structure, we first provide a standard abstract lower semi-continuity result that covers this setting and reduces lower semi-continuity to a closedness condition on the constraint. Proposition 3.2. Take K, N ∈ N and functions G : M N × M K → [0, ∞) and C : M N → P(M K ) such that G is lower semi-continuous, C(x) = ∅ for all x ∈ M N and for any bounded sequence (x n ) n with elements in M N and (y n ) n such that y n ∈ C(x n ) for all n, there exist subsequences (x n k ) k and (y n k ) k converging to some x ∈ M N , y ∈ M K such that y ∈ C(x).
Then, for any x ∈ M N there exists y ∈ C(x) such that S(x) = G(x, y). Further, S is lower semi-continuous.
Proof. See Section A.2 in the Appendix.
In the application of this result to the proposed distance-type functions, the constraint y ∈ C(x) will correspond to constraining points to particular positions on distance minimizing geodesics. Consequently, in order to verify the closedness condition, the general result that the set of shortest geodesics connecting two points is closed w.r.t. perturbation of these points will be necessary and is provided as follows. Then there exist subsequences ((x n k , y n k , u n k , v n k )) k and ((c n k ,ỹ n k )) k converging to (x, y, u, v) and (c,ỹ), respectively, such that In particular, C : Proof. Since each c n is on a length-minimizing geodesic between x n and v n , we get that d(x n , c n ) ≤ d(x n , v n ) and hence (c n ) n is bounded. Also d(c n ,ỹ n ) = d(c n , u n ) and hence (ỹ n ) n is bounded. Consequently we can pick subsequences ((x n k , y n k , u n k , v n k )) k and ((c n k ,ỹ n k )) k converging to (x, y, u, v) and (c,ỹ). Now with γ n k c : [0, 1] → M being a shortest geodesic connecting x n k and v n k such that γ n k c (1/2) = c n k we get by Lemma 3.3 that, up to a subsequence, γ n k c → γ c uniformly, where γ c is a again a shortest geodesic connecting x and v. Consequently, 2 . Now pick γ n k y : [0, 1] → M to be a length minimizing geodesic between γ n k y (0) = u n k and γ n k y (1) = c n k such that γ n k y (2) =ỹ n k . By Lemma 3.3 we get that, up to subsequences, γ n k y converges uniformly to a geodesic γỹ on [0, 2] such that γỹ is length-minimizing between γỹ(0) = u and γỹ(1) = c. By uniform convergence on [0, 2] we get that y = lim k y n k = lim k γ n k y (2) = γỹ(2) and the assertion follows.
Combining this with the general assertion of Proposition 3.2, existence and lower semicontinuity results for D S and D sym S follow as direct consequences. Lemma 3.5. The minimum in the definition of D S as in Equation (8) is attained and D S is lower semi-continuous. Also, the minimum in the definition of D sym S as in Equation (20) is attained and D sym S is lower semi-continuous.
Proof. See Section A.2 in the Appendix.
Using similar techniques, existence and lower semi-continuity results for the parallel transport variants can be established as follows. Lemma 3.6. The minimum in the definition of D pt as in Equation (11) is attained and D pt is lower semi-continuous. Also, the minimum in the definition of D sym pt as in Equation (22) is attained and D sym pt is lower semi-continuous.
Proof. See Section A.2 in the Appendix.

Algorithmic approach to TGV for manifold-valued data
In order to algorithmically approach denoising problems using TGV regularization in the manifold setting, we employ the concept of cyclic proximal point algorithms (CPPAs). A reference for cyclic proximal point algorithms in vector spaces is [22]. In the context of Hadamard spaces, the concept of CPPAs was developed by [8], where it is used to compute means and medians. In the context of variational regularization methods for nonlinear, manifold-valued data, they were first used in [93]. More precisely, the reference [93] deals with TV regularization as well as classical smooth methods for manifold-valued data. The CPPA approach was later used for higher-order TV-type methods in [18] for circle-valued data and in [10] for manifold-valued data.
Principle of a CPPA. The idea of a CPPAs is to compose a functional f : M → R into basic atoms f i and then to compute the proximal mappings of the f i in a cyclic, iterative way. More precisely, assume that and consider the proximal mappings [64,41,6] prox λfi : M → M given as One cycle of a CPPA then consists of applying each proximal mapping prox λfi once in a prescribed order, e.g., prox λf1 , prox λf2 , prox λf3 , . . . , or, generally, prox λf σ(1) , prox λf σ(2) , prox λf σ(3) , . . . , where the symbol σ is employed to denote a permutation. The cyclic nature is reflected in the fact that the output of prox λf σ(i) is used as input for prox λf σ(i+1) . Since the ith update is immediately used for the (i + 1)st step, it can be seen as a Gauss-Seidel-type scheme. A CPPA now consists of iterating these cycles, i.e., in the kth cycle, we have and x (k+1) 1 is obtained by applying prox λ k f σ(n) to x (k) n . Here, n denotes the number of elements in the cycle. During the iteration, the parameter λ k of the proximal mappings is successively decreased. In this way, the penalty for deviation from the previous iterate is successively increased.
It is chosen in a way such that the sequence (λ k ) k is square-summable but not summable. Provided that this condition on the stepsize parameters holds, the cyclic proximal point algorithm can be shown to converge to the optimal solution of the underlying minimization problem, at least in the context of Hadamard manifolds and convex (f i ) i [9, Theorem 3.1].
For general manifolds, the proximal mappings (25) are not globally defined, and the minimizers are not unique, at least for general possibly far apart points; cf. [41,6]. This is a general issue in the context of manifolds that are -in a certain sense -a local concept involving objects that are often only locally well defined. In case of ambiguities, we hence consider the above objects as set-valued quantities. Furthermore, we cannot guarantee -and in fact do not expectthe second distance-type functions D in the definition of the M-TGV 2 α functional to be convex. Hence convergence of the CPPA algorithm to a globally optimal solution cannot be ensured. It thus should be seen as an approximative strategy. Nevertheless, as will be seen in the numerical experiments section, we experience a good convergence behavior of the CPPA algorithm in practice. This was also observed in previous works, where the CPPA algorithm was employed to minimize second-order TV-type functionals [18,10], which are also non-convex. In order to be precise, we further point out that we approximately compute the proximal mappings of the distance-type functions D; we hence use an inexact proximal point algorithm.
We also point out that a parallel proximal point algorithm has been proposed in [93]. Here the proximal mappings of the f i are computed in parallel and then averaged using intrinsic means. In order to reduce the costs for averaging, an approximative variant of the parallel proximal point algorithm has been proposed in [93] as well. Principally, the cyclic proximal point algorithm actually applied in this paper might be replaced by this parallel strategy; the derivations in the following provide all information necessary.

Algorithms for manifold-valued TGV
We employ a cyclic proximal point algorithm for the minimization of a TGV-regularized variational approach in the manifold context. For simplicity, we consider the univariate setting first; all aspects discussed in the univariate situation are prototypic and basic for the bivariate situation. Starting with the general version of TGV for manifolds as in (3), we aim to solve the denoising problem min (ui)i,(yi)i where (h i ) i , being a finite sequence in M, denotes the given data. We decompose (27) into the atoms Now we use this decomposition into the atoms g i , g i , g i for the decomposition (24) for the CPPA. We apply the iteration (26) to this decomposition. We remark that the data terms g i are not coupled w.r.t. the index i. The same applies to the g i . This allows the parallel computation of the proximal mappings of the g i for all i and, separately, of the g i for all i. The g i are actually coupled. However, grouping even and odd indices, we see that the g i , i even, are not coupled and that the g i , i odd, are not coupled. So we may compute their proximal mappings in parallel. Together, this results in a cycle of length four per iteration, one step for the g i , one step for the g i and two steps (even, odd) for the g i . In the following, the task is to compute the proximal mappings of the atoms g i , g i , g i of (28). To this end, we will from now on always assume that the involved points are locally sufficiently close such that there exist unique length minimizing geodesics connecting them. We note that this restriction is actually not severe: in the general (non-local) situation we have to remove the cut points -which are a set of measure zero -to end up with the corresponding setup.
We remark that, as we will see, an explicit computation of the proximal mapping is possible for g i and g i , but not for g i . We believe that it is an important feature of the proposed definition of manifold-TGV via point-tuples that the proximal mappings of the first term g i is still explicit and hence only one part of the overall problem does not allow an explicit proximal mapping. Indeed, also existing generalizations of second-order TV to the manifold setting incorporate one part with non-explicit proximal mappings [10]. Hence, from the algorithmic viewpoint, the step from second-order TV to our proposed version of TGV does not introduce essential additional computational effort.
For the data atoms g i , the proximal mappings prox λgi are given by They have been derived in [41]. We recall that the symbol [·, ·] t denotes the point reached after time t on the (non unit speed) geodesic starting at the first argument reaching the second argument at time 1. Further, the proximal mappings of the distance terms g i have a closed form representation as well (see [93]), which is given by where t = λα 1 /d(u i+1 , y i ), if λα 1 < 1 2 , and t = 1/2 else. Note that for determining the position, we view u, y as column vectors of a matrix with two columns.
The next task is to compute the proximal mappings of the g i . Unfortunately, no closed form for the proximal mappings of the g i is available. Instead, we use a subgradient descent scheme to compute these proximal mappings, i.e., to solve the problem min u,y where the u j , y j are the input of the proximal mapping. A subgradient descent scheme has already been used to compute proximal mappings in [10]. Looking at (31), the optimal u, y fulfill u j = u j , y j = y j whenever j / ∈ {i, i − 1}. Hence we may restrict ourselves to consider the four arguments u i , y i , u i−1 , y i−1 . Further the gradient of the mapping u j → 1 2 j d(u j , u j ) 2 is given as − log uj u j . (For background information on (Riemannian) differential geometry, we refer to the books [76,38].) So we have to compute the (sub)gradient of the mapping D([u i , y i ], [u i−1 , y i−1 ]) as a function of u i , y i , u i−1 , y i−1 . To this end, we have to specify the setting to the two versions of D as proposed in Section 2; the respective derivations are topics of Section 4.2 and Section 4.3.
We next consider the bivariate denoising problem. We use M-TGV 2 α regularization with the M-TGV 2 α functional given in (15) where, as pointed out, we focus on p = 1. We obtain the denoising problem min (ui,j )i,j ,(y 1 i,j )i,j ,(y 2 i,j )i,j where the (h ij ) ij denotes the bivariate data. We decompose the objective functional in (32) into the atoms i (u, y 1 , y 2 ) := α 1 d(u i,j+1 , y 1 i,j ); g (4) ij (u, y 1 , We use the decomposition w.r.t. the atoms g (k) ij as the decomposition (24) for the CPPA, and apply the iteration (26) w.r.t. the derived decomposition. The proximal mappings prox λg (1) ij are given by (29), and the proximal mappings prox λg (2) ij , prox λg (3) ij are given by (30). The proximal mappings prox λg (4) ij , prox λg (5) ij are computed as in the univariate case explained above. It remains to compute the proximal mappings of the atoms g (6) ij . As before, there is no explicit formula available and we use a subgradient descent for their computation as well. The computation of the corresponding derivatives are topics of Section 4.2 and Section 4.3.

Riemannian gradients for the Schild version
We here derive the derivatives needed to compute the proximal mappings of the mappings D and D sym in Section 4.1 when specifying to the Schild's variant. We start out to compute the derivative of the mapping D S given by Recall that, for our computations, we assume that geodesics are unique, which is, as pointed out above, the case up to a set of measure zero. We directly see that D S is symmetric w.r.t. interchanging y i−1 , u i . Hence we only have to compute the differential with respect to one of these variables y i−1 , u i . Furthermore, the gradient ∇ yi D S w.r.t. the fourth argument y i of D S for y i with ] 2 is just the well-known gradient of the Riemannian distance function which is known to be given by [38,4] where | · | is the norm in the tangent space associated with the point y i , i.e., the square root of the Riemannian scalar product. Here, denotes the result of applying the Schild's construction to the respective arguments. In order to determine the gradients w.r.t. the other variables we have to apply the adjoint of the respective differentials to the gradient of the distance function w.r.t. the first argument. More precisely, we have to calculate where T 1 is the adjoint of the differential of the mapping ] 2 , and where T 2 is the adjoint of the differential of the mapping The differential w.r.t. u i is obtained by symmetry as pointed out above. In the following we derive these adjoint mappings in terms of Jacobi fields.
To this aim, we first explain the notion of a Jacobi field and then point out the connection with geodesic variations. In a Riemannian manifold M, the Riemannian curvature tensor R is given by Z where X, Y, Z are vector fields and ∇ denotes the Levi-Civita connection, and where [X, Y ] = XY − Y X denotes the Lie bracket of the vector fields X, Y. A Jacobi field Y along a geodesic γ is the solution of the differential equation The space of Jacobi fields is a 2N -dimensional linear space where N denotes the dimension of M. The connection to geodesic variations is as follows: the "derivative" vector field Y (s) = d dt | t=0 V (s, t) of a geodesic variation V is a Jacobi field, and conversely, any Jacobi field is obtained from a geodesic variation. Further details may be found in the books [34,76,38].
We also need the notion of (Riemannian) symmetric spaces for the formulation of the next lemma. A symmetric space is a Riemannian manifold for which its curvature tensor R is invariant under parallel transport, R(pt x,y X, pt x,y Y ) pt x,y Z = pt x,y R(X, Y )Z where X, Y, Z are tangent vectors at x and pt x,y denotes the parallel transport from x to y along a curve γ (here not reflected in the notation). Then, in particular, the covariant derivative of the curvature tensor R along any curve equals zero in a symmetric space. A reference for symmetric spaces is [34]. Proof. The proof is given in Section A.3.
Regarding the fact that we consider points with D S = 0, we remark that the points u i , where T 3 is the adjoint of the derivative of the mapping m → [u i−1 , m] 2 , and where T 4 is the adjoint of the derivative of the mapping , that is, T 2 as in (38) is given as Further T 3 and T 4 can be computed using Jacobi fields.
Proof. The proof is given in Section A.3.
For the class of symmetric spaces, we make the mappings T 3 , T 4 more explicit. . We consider an orthonormal basis v n of eigenvectors of the self-adjoint Jacobi operator J → R( γ (1) |γ (1)| , J) γ (1) |γ (1)| with corresponding eigenvalues λ n , and v 1 tangent to γ. W.r.t. this basis, the mapping T 3 is given by where, using the symbol d for the distance between u i−1 and Here, the α n are the coefficients of the corresponding basis representation. Further, let ξ : [0, 1] → M be the geodesic connecting y i−1 = ξ(0) and u i = ξ(1) and w n be an orthonormal basis of eigenvectors of J → R( ξ (0) |ξ (0)| , J) ξ (0) |ξ (0)| with eigenvalues µ n , and w 1 tangent to ξ. Then, w.r.t. this basis, the mapping T 4 is given by where, using the symbol d for the distance between y i−1 and u i , Here, the β n are the coefficients of the corresponding basis representation.
Proof. The proof is given in Section A.3.
Next, we consider the Riemannian gradients for the Schild version D sym = D sym S as given in (20) which are needed for the CPPA for the bivariate TGV functional. To this end, we that note that D sym S can be expressed in terms of the function S of (36) by Hence, the mapping D sym S (which depends on seven arguments) is an expression of the Riemannian distance function and three realizations of S. We have seen how to differentiate S in Lemma 4.1 and Lemma 4.2. We derive the gradient of D sym S by iterated application of these mappings and the concepts of these lemmata. We point out the symmetry of D sym S with respect to u i−1,j , y 1 i,j , and with respect to u i,j−1 , y 2 i,j . This reduces the task to actually considering five different arguments. The derivative of D sym S is provided in the following proposition. Its proof is a direct consequence of the previous results and the decomposition of D sym S into decompositions of the function S and the distance function d as above. i,j is given by with the adjoint operators T 3 , T 4 given as in Lemma  i,j is given by with the adjoint operators T 3 , T 4 as in Lemma  i,j−1 is given by with the adjoint operator T 1 given as in Lemma 4.1 formed w.r.t. the points y 1 i,j−1 , r 2 , r 1 . The gradient with respect to the variable y 2 i−1,j is simply given by Finally, the gradient of D sym Here the adjoint operators T 3 and T 4 are given as in Lemma 4.2 formed w.r.t. the points y 1 i,j−1 , r 2 , r 1 , T 1 is given as in Lemma 4.1 formed w.r.t. the points u i,j , y 1 i,j , u i−1,j andT 1 is given as in Lemma 4.1 formed w.r.t. the points u i,j , y 2 i,j , u i,j−1 .

Riemannian gradients for the parallel transport version
We here derive the derivatives needed to compute the proximal mappings in Section 4.1 when specifying to the parallel transport variant. We use the short-hand notation For the following computations recall that, by our assumption of uniqueness of length-minimizing geodesics, the log mapping, initially defined as set-valued, always maps to single points in the tangent space.
Lemma 4.5. The function F given by (50) is symmetric with respect to interchanging (u i , y i ) with (u i−1 , y i−1 ). In particular, for points with F = 0, the gradient of F w.r.t. the third variable y i , is given by ∇ yi F (u i , u i−1 , y i , y i−1 ) = ∇ yi−1 F (u i−1 , u i , y i−1 , y i ). Further, again for points with F = 0, the gradient of the function F w.r.t. the first component variable u i , is given by Proof. The proof is given in Section A.3 by specifying F = F 0 in Lemma A.5.
Regarding the fact that we consider points with F = 0 in the above statement, we remark that the points u i , u i−1 , y i , y i−1 such that F (u i , u i−1 , y i , y i−1 ) = 0 form a set of measure zero. Only on this zero set, F is not differentiable. Further, in this case, the four-tuple consisting of the four zero-tangent vectors sitting in u i , u i−1 , y i , y i−1 belong to the subgradient of F. We note that Lemma 4.5 tells us that, in order to compute the gradient of the second order type difference F , we only need to compute the respective gradient of F w.r.t. u i−1 and y i−1 . This is done in the following. Lemma 4.6. The gradient of the function F for points with F = 0 w.r.t. the variable y i−1 is given by where T = (d yi−1 log ui−1 ) * is the adjoint of the (Fréchet) derivative (denoted by the symbol d yi−1 ) of the log mapping w.r.t. the variable y i−1 . T can be computed using Jacobi fields.
Proof. The proof is given in Section A.3.
In Riemannian symmetric spaces the above mapping T can be made more explicit.
where the α n are the coefficients of the corresponding basis representation and the function f, depending on the sign of λ n is given by Proof. The proof is given in Section A.3.
Finally, we consider the gradient of F w.r.t. the variable u i−1 .
Lemma 4.8. The gradient of the function F for points with F = 0 w.r.t. the variable u i−1 is given by where the (v n ) n form an orthonormal basis of the tangent space at u i−1 , and the coefficients α n are given by Here L n t , B n t denote the vector fields L, B (defined by L : u i−1 → log ui−1 y i−1 and B : u i−1 → pt ui,ui−1 z, where z = log ui y i ,) along the (specific) geodesic t → exp ui−1 tv n , t ∈ [0, 1], determined by v n , i.e., L n t = log exp u i−1 tvn y i−1 and B n t = pt ui,exp u i−1 tvn z, with z = log ui y i .
Proof. The proof is given in Section A.3.
The precise computation of D dt | t=0 L n t in symmetric spaces is topic of Lemma 4.9. Further, the computation of D dt | t=0 B n t is carried out for the manifolds explicitly considered in this paper: this is done for the sphere in Lemma 4.10, and for the space of positive matrices in Lemma 4.11. Lemma 4.11 and its proof may serve as a prototypic guide for deriving similar expressions for other symmetric spaces such as the rotation groups or the Grassmannians for instance. We note that the approach is by no means restricted to the two considered classes of spaces and might serve as a guide for other manifolds; we only did not derive a more explicit representation on the general level of symmetric spaces. We further note that numerical differentiation of the particular term is a second option as well. Lemma 4.9. Assume that the manifold M is a symmetric space. Let (v n ) n be an orthonormal basis of eigenvectors of the self-adjoint Jacobi operator J → R( γ (0) |γ (0)| , J) γ (0) |γ (0)| where the (constant speed) geodesic γ : [0, 1] → M connects u i−1 = γ(0) and y i−1 = γ(1), and where R denotes the Riemannian curvature tensor. For each n, we denote by λ n the eigenvalue associated with v n .
The covariate derivatives D dt | t=0 L n Here d = d(u i−1 , y i−1 ) denotes the length of the geodesic connecting u i−1 , y i−1 . (If the term √ λ n d = 0 in the denominators of the second line in (57), then u i−1 = y i−1 , and the formula is still valid since we are facing a removable singularity then.) Proof. The proof is given in Section A.3.
and v n to the left of log ui−1 u i (otherwise multiplied by −1 accounting for the change of orientation). Here the skew-symmetric matrix L ω = 0 ω −ω 0 is taken w.r.t. the basis {log ui−1 u i , where w is the vector determined by w ⊥ log ui−1 u i , |w| = 1, and w is to the left of log ui−1 u i . In other words, we have to multiply the second line of (58) with the signed length of the projection of v n to the normalized vector (log ui−1 u i ) ⊥ . If u i = u i−1 , then the differential D dt | t=0 B n t = 0 (which is consistent with letting d → 0 in the above formulae.) Proof. The proof is given in Section A.3. Lemma 4.11. Let M be the space of symmetric positive definite matrices. Then, the covariate derivative D dt | t=0 B n t (which is a tangent vector sitting in u i−1 ) is given by the following sum of matrices where   (10) and (22). Then we can use this analogue of the decomposition (33) and apply the CPPA iteration (26) to this decomposition. The proximal mappings prox λg (1) ij are given by (29), prox λg (2) ij , prox λg (3) ij are given by (30) as for the Schild variant above. The proximal mappings prox λg (4) ij , prox λg (5) ij are computed as in the univariate case considered in Section 4.3. In order to compute the proximal mappings of the atoms g (6) ij using a subgradient descent, it is necessary to differentiate the mapping D sym pt of (22) (which is the analogue of D sym S for the Schild case) w.r.t. its seven arguments. As in the Schild case, it is possible to decompose the mapping into simpler functions which we have already considered in Section 4.3. We do not carry out the rather space consuming derivation here.

Numerical results
This section provides numerical results for synthetic and real signals and images. We first describe the experimental setup.
We carry out experiments for S 1 , S 2 and Pos 3 (the manifold of symmetric positive definite 3 × 3 matrices equipped with the Fisher-Rao metric) valued data. S 1 data is visualized by the phase angle, and color-coded as hue value in the HSV color space when displaying image data.
We visualize S 2 data either by a scatter plot on the sphere as in Figure 6, or by a color coding as in Figure 9. Data on the Pos 3 manifold is visualized by the isosurfaces of the corresponding quadratic forms. More precisely, the ellipse visualizing the point f p at voxel p are the points x fulfilling (x − p) f p (x − p) = c, for some c > 0.
To quantitatively measure the quality of a reconstruction, we use the manifold variant of the signal-to-noise ratio improvement see [86,93]. Here f is the noisy data, g the ground truth, and u a regularized restoration. A higher ∆SNR value means better reconstruction quality. For adjusting the model parameters α 0 , α 1 of M-TGV 2 α , it is convenient to parametrize them by , and α 1 = r s min(s, 1 − s) , so that r ∈ (0, ∞) controls the overall regularization strength and s ∈ (0, 1) the balance of the two TGV penalties. For s → 0 we get α 0 → ∞ and α 1 = r, so that TGV minimization approximates the minimization of TV modulo a linear term which can be added at no cost. For s → 1 we have α 0 = r and α 1 → ∞ which corresponds to pure second-order TV regularization. One may think of r as the parameter mainly depending on the noise level, and of s as the parameter mainly depending on the geometry of the underlying image. Figure 5 illustrates the influence on r and s for a synthetic S 1 -valued image corrupted by von Mises noise with concentration parameter κ = 5. There and in most of the following experiments, we observed satisfactory results for the fixed value s = 0.3. Based on these observations, we suggest to use s = 0.3 as a starting point, to decrease or increase it if the image is dominated by edges or smooth parts, respectively. We have implemented the presented methods in Matlab 2016b. We use the toolbox MVIRT [10] for the basic manifold operations such as log, exp and parallel transport, and for parts of the visualization 1 . We used 100000 iterations for all experiments with univariate data and 1000 iterations for the image data. The cooling sequence (λ k ) k∈N used as stepsize in the gradient descent for computing the non-explicit proximal mappings was chosen as λ k = λ 0 k −τ with τ = 0.55. For the univariate spherical data we observed faster convergence when using a stagewise cooling, i.e., letting the sequence fixed to λ 0 for 500 iterations in the first stage, use the moderate cooling τ = 0.35 in the second stage until iteration 1000 and then the cooling τ = 0.55 afterwards.

Numerical evaluation of the proposed models
Comparison between Schild variant and parallel transport variant. First we compare the two proposed realizations of manifold TGV: the parallel transport variant and the Schild variant. Figure 6 shows the results with both variants for some typical univariate signals for the S 1 , the S 2 and the Pos 3 manifold. (The S 1 -valued signal was corrupted by von Mises noise with concentration parameter κ = 10. The S 2 -valued signal was taken from [10] and corrupted by applying the exponential mapping of Gaussian distributed tangent vectors as in [10] with σ = 0.1. The Pos 3 -valued signal was corrupted by applying the exponential mapping of Gaussian noise distributed on the tangent vectors with σ = 0.2. The r-parameters were chosen as r = 1 for the spherical signals and as r = 0.    some points, e.g. near the discontinuity of the second to fourth position of the signal and at the last position. In summary, we observe a qualitatively similar result in the sense that the geodesic parts are well reconstructed and that the edges are preserved. In the following, we focus on the Schild variant.
Comparison of manifold TGV with vector-space TGV. In order to validate both our generalizations of TGV to the manifold setting as well as numerical feasibility of our optimization algorithms, we carry out a comparison to the vector-space case. This is done for S 1 , the unit sphere in R 2 . By generating a signal with values that are strictly contained in one hemisphere, we can unroll the signal and compare to vector-space TGV-denoising in a way that the same results can be expected. We carried out this comparison for synthetic data without noise and different values of the balancing parameter s. We tested the setting of s = 0.3, which comprises a good balance between first and second order terms, as well as the rather extreme settings of s = 0.05 and s = 0.95. The regularization parameter r was fixed to r = 1. In order to approximate the ground-truth solution of second-order TGV denoising in the vector space setting, we implemented the Chambolle-Pock algorithm [32] for this situation. To ensure a close approximation of the ground truth for all parameter settings, we carried out a dedicate stepsize tuning to accelerate convergence of the algorithm, computed 2 × 10 5 iterations and ensured optimality by measuring the duality gap. The result of this evaluation can be found in Figure 7. It can be observed that the qualitative properties of the numerical solution obtained with the manifold-TGV code are similar to the ones of the approximate ground-truth of the vector space setting, confirming overall feasibility of our model and implementation.
For the case s = 0.05, one can see in particular in the right part of the bottom plot (starting with the first plateau after 180 on the x-axis), that the solution is piecewise constant up to a linear part, which is approximately the same for all four plateaus. This is what one would expect for the extreme case s → 0, as in this case, TGV minimization mimics TV minimization up to a globally linear term. On the other hand, for the case s = 0.95, one can see that the solution is piecewise linear with no jumps. This can be expected from second-order TV minimization, which coincides with second order TGV minimization for s → 1. In particular, the piecewise linear part on the left is approximated well, while the plateaus on the right are not well captured. The case s = 0.3 provides a good compromise here: The linear parts on the left are still well captured, but the solution still admits jumps on the plateau parts, as can be seen in particular outermost right jump after point.
We also note that there are still some differences of the solution obtained with the manifold code and the Chambolle-Pock result interpreted as approximate ground truth of the vector space case, in particular for the cases s = 0.05 and s = 0.95. We believe that this is mostly an issue of the algorithmic realization rather than the model itself caused by the numerical solutions obtained for the more extreme cases s = 0.05 and s = 0.95.
Basic situations. The aim of this experiment is to investigate the performance of the proposed manifold-TGV model for different basic situations and noise levels. The experiment is carried out on a univariate signal on the two-sphere S 2 , that comprises jumps in the signal as well as its derivative and is composed of piecewise constant and geodesic signals.
The results of the experiment for different noise levels (no noise, intermediate noise, strong noise), a fixed value of r = 1 and different values of s, namely s ∈ {0.05, 0.95, 0.3}, can be found in Figure 8. As can be seen there, the manifold-TGV functional regularizes the data quite well and is able to achieve good results, even in the case with relatively strong noise where it is hard to see any structure in the data. As in the previous experiment, the choice s = 0.05 promotes piecewise constant solutions, which is naturally best for regions where the signal is piecewise constant, but leads to "staircasing" in smooth regions, as can be particularly seen in the geodesic parts of the strong-noise case. The choice s = 0.95 approximates geodesics well, but does not allow for jumps and also produces some oscillations in the smooth parts of the case with strong noise. Again the choice s = 0.3 is a good compromise. Even though it does not reconstruct jumps as well as the TV-like version, it still allows for jumps and reconstructs the geodesic parts rather accurately.

Results for synthetic data
Next we investigate the denoising performance of M-TGV 2 α on synthetic images. We compare the results of the proposed algorithm with the results of manifold TV regularization [93], and  with the result of second-order manifold TV, denoted by TV 2 [10]. The model parameter of first-order TV is denoted by α and that of TV 2 by β. For all algorithms, 1000 iterations were used.
First we look at S 2 -valued images. As in [10], the noisy data f were created from the original image g by f ij = exp gij η ij where η ij is a tangent vector at g ij and both its components are Gaussian distributed with standard deviation σ = 4 45 π. For comparison with first order TV we scanned the parameter range α = 0.1, 0.2, . . . , 1 and the special parameters α = 0.22 and α = 3.5·10 −2 given in [10] and the corresponding implementation. Similarly, for TV 2 we scanned the parameter range β = 1, 5, 10, 15 . . . , 30, and the special parameters β = 8.6 and β = 29.5 suggested in [10]. As no beneficial effect of combining first and second-order TV was observed in [10], we used pure TV and TV 2 regularization. For the proposed method, we fixed s = 0.3 and report the best result among the six parameters r = 0.05, 0.1, 0.15, . . . , 0.3. The results in Figure 9 show that the second-order methods give a significantly smoother result than first order TV and that they do not suffer from staircasing effects. On the flipside, TV 2 regularization results in an undesired smoothing effect of the edges between the tiles, best seen at the bottom left tile. The proposed TGV regularization preserves these edges which results in an improved reconstruction quality.
Next we look at Pos 3 -valued images. Such images appear naturally in diffusion tensor imaging (DTI), so we briefly describe the setup. DTI is a medical imaging modality based on diffusion weighted magnetic resonance images (DWI) which for example allows to reconstruct fiber tract orientations [11]. A DWI captures the diffusivity of water molecules with respect to a direction v ∈ R 3 . The relation between a diffusion tensor f (p) and the DWIs D v (p) at the voxel p is given by the Stejskal-Tanner equation with constants b, A 0 > 0. A standard noise model for the DWIs is the Rician model which assumes a complex-valued Gaussian noise in the original frequency domain measurements [12]. This means that assuming the model (62) the actual measurement in direction v at pixel p is given by D v (p) = ((X + D v (p)) 2 + Y 2 ) 1/2 , with the Gaussian random variables X, Y ∼ N (0, σ 2 ). Building on this model, we impose the noise as follows. From the synthetic tensor image, we compute DWIs according to (62) with respect to 21 different directions. Then we impose Rician noise to all derived DWIs, and we estimate the corresponding diffusion tensor image f using a least squares approach on (62). Due to the noise, some of the fitted tensors might not be positive definite and thus have no meaningful interpretation. To handle such invalid tensor we exclude them from the data fidelity term. On the other hand, we keep the corresponding pixels in the regularizing term so that we achieve a reasonable inpainting. In Figure 10, we illustrate the effects of the regularization for a synthetic Pos 3 -valued image. For comparison with (combined) first and second-order TV, we scanned all combinations of the parameters α = 0, 0.01, 0.035, 0.1, 0.3125, 0.375, and β = 0.02, 0.05, 0.125, 0.625. As before, this comprises the parameters used in [10] and the corresponding implementation for that manifold. For the proposed method, we report the best result among the possible combinations of the parameters r = 0.06, 0.08, 0.09, 0.1, 0.12 and s = 0.3, 0.35, 0.4, 0.45. In the example we observe that combined TV/TV 2 regularization yields a similar reconstruction quality as pure first order TV, whereas S-TGV gives a significantly higher reconstruction quality.

Results for real data
We illustrate the effects of TGV regularization on real manifold-valued signals and images. First we look at smoothing a time series of wind directions. The natural data space for a signal of orientations is the unit circle S 1 . The present data shows wind directions recorded every hour in 2016 by the weather station SAUF1, St. Augustine, FL.; see Figure 11 2 . We observe that the proposed method smoothes the orientations while respecting the phase jump from −π to π and preserving linear trends in the data.
Next we look at smoothing of interferometric synthetic aperture radar (InSAR) images. Synthetic aperture radar (SAR) is a radar technique for sensing the earth's surface from measurements taken by aircrafts or satellites. InSAR images consist of the phase difference between two SAR images, recording a region of interest either from two different angles of view or at two different points in time. Important applications of InSAR are the creation of accurate digital elevation models and the detection of terrain changes; cf. [62,71]. As InSAR data consists of phase values that are defined modulo 2π, their natural data space is the unit circle. In Figure  12, we illustrate the effect of S-TGV to an InSAR image taken from [83] 3 . We observe a clear smoothing effect while sudden phase changes are preserved.
At last we consider real DTI data which was taken from the Camino project [36] 4 . In Figure 13, we see an axial slice of a human brain (slice 20 of the Camino data). We also display a zoom to the corpus callosum which connects the right and the left hemisphere. The   original tensors were computed from the diffusion weighted images by a least squares method based on the Stejskal-Tanner equation. We observe that the proposed method smoothes the tensors in the corpus callosum while it preserves the sharp edges to the adjacent tissue. Also observe the inpainting of the invalid tensors.

Conclusion
In this work, we have introduced and explored a notion of second-order total generalized variation (TGV) regularization for manifold-valued data. First, we have derived a variational model for total generalized variation for manifold-valued data. For this purpose, we have used an axiomatic approach. We have first formalized reasonable fundamental properties of vector-valued TGV that should be conserved in the manifold setting. Then we have proposed two realizations which we have shown to fulfill the required axioms. The realization based on parallel transport is rather natural -although not straight forward -from the point of view of differential geometry. The realization based on the Schild's ladder may be seen as an approximation of the parallel transport variant. It requires less differential geometric concepts and it is easier to realize numerically while yielding comparable numerical results as shown in the experimental section. Existence results for M-TGV 2 α -based denoising have been obtained for the proposed variants. Next, we have derived an algorithm for the proposed model. To this end we built on the well-established concept of a cyclic proximal algorithm. As main contribution, we have performed the challenging task to derive all quantities necessary to compute the proximal mappings of the involved atoms. Finally, we have conducted a numerical study of the proposed scheme. The experiments revealed that M-TGV regularization reliably removes noise while it preserves edges and smooth trends. A quantitative comparison on images with groundtruth indicates that TGV regularization improves upon first and second-order TV for manifold-valued data.
An interesting topic of future research are improved numerics such as earlier stopping criteria.  (15) is attained at y 1 i,j = y 2 i,j = u i,j we get, since On the other hand, in case the minimum is attained at y 1 i,j = u i+1,j and y 2 i,j = u i,j+1 we get, in the vector space situation, that which coincides with α 0 TV 2 and TV 2 as in Definition 2.4. Now suppose that (u i,j ) i,j is locally geodesic. Choosing y 1 i,j = u i+1,j and y 2 i,j = u i,j+1 for all i, j we get, as in the proof of Theorem 2.9, that M-TGV 2 α is bounded above by which is zero by assumption. Now conversely, assume that M-TGV 2 α (u) = 0 and any two points in u i,j , u i ,j with max{|i − i |, |j − j |} ≤ 2 are connected by unique geodesics. As in the proof of Proposition 2.8 this implies that both (u i,j ) i and (u i,j ) j are locally geodesic as univariate signals. Also we get that which implies by assumption, that (u i+k,j−k ) k is locally geodesic. Hence, by definition, (u i,j ) i,j is locally geodesic as bivariate signal. Finally, by Lemma 2.16, u is geodesic.
and, using Proposition 2.17, (M-P1') to (M-P3') follow. Now suppose that (u i,j ) i,j is a geodesic three-by-three signal. Hence, with the notation as in the definition of D sym S , c 1 = u i,j ∈ [u i−1,j , u i+1,j ] 1 2 and we can choose But, as shown in the proof of Theorem 2.9, the right-hand-side vanishes since (u i−k,j+k ) k is geodesic as univariate signal. Now conversely, assume that u = (u i,j ) i,j is a three-by-three signal such that both (u i, ) i and (u i,j ) j are locally geodesic and that By the assumption on uniqueness of geodesics we get, again with the notation as in the definition of D sym S , that c 1 = u i,j and hence r 1 = u i,j and also that c 2 = u i,j and hence . But the latter being zero implies that u i−1,j+1 ∈ [u i+1,j−1 , u i,j ] 2 and, consequently, u is geodesic.
Hence also the remaining assumptions of Proposition 2.17 are satisfied and the assertion follows.
The following well-known fact on the parallel transport in manifolds will be required in the proof of Theorem 2.19. For the sake of completeness, we provide a short proof.
and from Proposition 2.17 it follows that (M-P1') to (M-P3') holds. Now let u = (u i,j ) i,j be a three-by-three geodesic signal and choose y 1 i,j = u i+1,j and y 2 i,j = u i,j+1 for each (i, j). Denoting by γ 1 and γ 2 distance minimizing geodesics such that γ 1 (0) = u i−1,j , γ 1 (1/2) = u i,j , γ 1 (1) = u i+1,j and γ 2 (0) = u i,j−1 , γ 2 (1/2) = u i,j , γ 2 (1) = u i,j−1 , respectively, Lemma A.1 implies that log γ 1 ui−1,j (u i,j ) = pt ui−1,j (log γ 1 ui,j (u i+1,j )) as well as log γ 2 ui,j−1 (u i,j ) = pt ui,j−1 (log γ 2 ui,j (u i,j+1 )). Hence we can choose r 1 = r 2 = u i,j and where the last term is zero since (u i−k,j+k ) k is locally geodesic as univariate signal. Now conversely, assume that u = (u i,j ) i,j is a three-by-three signal such that the geodesics connecting each pair of points are unique and such that (u i,j ) i and (u i,j ) j are locally geodesic. By Lemma A.1 and uniqueness of geodesics we get log ui−1,j (u i,j ) = pt ui−1,j (log ui,j (u i+1,j )) as well as log ui,j−1 (u i,j ) = pt ui,j−1 (log ui,j (u i,j+1 )). Again by uniqueness of geodesics, with the notation as in the definition of D sym pt , we conclude that r 1 = r 2 = u i,j hence, From the properties of D pt it hence follows that the points u i+1,j−1 , u i,j , u i−1,j+1 are on a joint, length minimizing geodesic at equal distance. Hence u is geodesic. This ensures that all assumptions of Proposition 2.17 are fulfilled and hence the assertion follows.

A.2 Proofs for Section 3
Proof of Proposition 3.2. Take (y n ) n such that y n ∈ C(x) and S(x) = lim n G(x, y n ). By assumption (y n ) n admits a subsequence (y n k ) k converging to some y ∈ C(x). By lower semi-continuity of G it follows that S(x) = G(x, y).
For lower semi-continuity of S now take (x n ) n in M N converging to some x ∈ M N for which, without loss of generality, we assume that lim inf n S(x n ) = lim n S(x n ). Pick y n ∈ C(x n ) such that S(x n ) = G(x n , y n ). By assumption and uniqueness of limits we can obtain a subsequences (x n k ) k and (y n k ) k converging to x and y, respectively, such that y ∈ C(x). We conclude that S(x) ≤ G(x, y) ≤ lim inf k G(x n k , y n k ) = lim inf k S(x n k ) = lim inf n S(x n ) and the assertion follows.
Proof of Lemma 3.3. The proof relies on the following fact: for any bounded subset N of M, there is a constant D such that, if the length of any geodesic ψ : [r, s] → M with ψ([r, s]) ⊂ N is smaller than D, then ψ is a unique distance minimizing geodesic between ψ(r) and ψ(s), and the Jacobi fields have no zero along this geodesic ψ; see for instance [34, Chapter 1.6/1.7]. In addition, there is a constant L which is uniform for all such geodesics such that, with f : [r, s] × [0, 1] → M a geodesic variation of ψ with f (·, 0) = ψ, we have In order to obtain uniform convergence of (γ n ) n , our approach is now to use compactness arguments and subdivide the curves into small segments where these assertions hold.
Further, by construction, γ n ([0, 1]) ⊂ N for any n. With this choice of N , we now choose the constant D as above and we choose l ∈ N large enough such that l ≥ 2 d(p,q) statements hold forw n l with the same argument. So, by going to subsequencesz nr ofz n l , and w nr ofw n l , we get thatz nr → z,w nr → w, for a tangent vector z sitting in x and a tangent vector w sitting in u. We have that the limit z ∈ log x (y), and that w ∈ log u (v), since the exponential map from T M → M is differentiable in a Riemannian manifold which implies exp x z = lim r→∞ exp x nr z nr = lim r→∞ y nr = y, and exp u w = v.
We are now prepared to estimate | pt γ γ(0),γ(1) z − w| which allows us to bound D pt ([x, y], [u, v]) from above; see (64). Since the parallel transport is an isometry, we have We further have that which is a consequence of the parallel transport along nearby geodesics being differentiably dependent on the variation of the geodesics. Using (70) together with (69), we have ≤ pt γ nr γ nr (0),γ nr (1) pt x,x nrz nr − pt u,u nrw nr + |z nr − z| + |w nr − w| + ε r .
The second inequality is a consequence of the triangle inequality together with the fact that parallel transport is an isometry. We take the limit w.r.t. r on the right hand side of (71): by (68), we have thatz nr → z, and thatw nr → w and by (70)  The second inequality in (72) is a consequence of our choice of z n k , w n k , and γ n k in (66) and the definition ofz n l ,w n l as in (67), and the equality in the last line follows by our choice of subsequences in (65). Then passing to the infimum according to (64) shows (63) which shows the first assertion. To see that the infimum in (11) is attained for any x, y, u, v ∈ M, we choose the constant sequences x n := x, y n := y, u n := u, and v n := v, and apply the previous result to these sequences. This shows the second assertion and completes the proof.
Having shown Lemma A.2 we can now employ Proposition 3.2 to show existence and lower semi-continuity results for D sym pt . As a preparation we need the following lemma. Lemma A.3. Let u n → u,ũ n →ũ, and y n → y in the complete manifold M. We consider a sequence (r n ) n with r n ∈ exp(ptũn (w n )) with w n ∈ log u n y n .
Then there is a subsequence (r n k ) k which converges and a limit r = lim k→∞ r n k such that r fulfills r ∈ exp(ptũ(w)) for some w ∈ log u y.
Proof. The present proof essentially employs the techniques already used in the proof of Lemma A.2. For this reason we keep the following arguments rather short. The sequencew n := pt u (w n ) is bounded since parallel transport is an isometry. So, by going to a subsequencew n l ofw n , we get thatw n l → w for a tangent vector w sitting in u. We have that the limit w ∈ log u (y), since the exponential map T M → M is differentiable which implies exp u w = lim l→∞ exp u n l w n l = lim l→∞ y n l = y. For each l we choose a distance-minimizing geodesic γ n l connecting γ n l (0) = u n l with γ n l (1) =ũ n l . Then we use Lemma 3.3 to choose a subsequence of geodesics γ n k of γ n l (suppressing iterated subindices) which uniformly converges to a geodesic γ : [0, 1] → M connecting γ(0) = u with γ(1) =ũ. Then, since parallel transport along nearby geodesics is differentiably dependent on the variation of the geodesics and by an argument similar to the one used for the convergence in Equation (70), ptũn k ,ũ pt γ n k u n k ,ũ n k w n k → pt γ u,ũ w as k → ∞. Then, by the continuity of the exponential map, r = lim k→∞ r n k = lim k→∞ exp(pt γ n k u n k ,ũ n k (w n k )) = exp(pt γ u,ũ (w)), i.e., r n k converges and (74) holds true which was the assertion to show.
Finally, the proof of Lemma 3.6 follows as consequence.
Proof of Lemma 4.2. We observe that the concatenation of the mappings m → [u i−1 , m] 2 and y i−1 → [u i , y i−1 ] 1 2 equals the mapping y i−1 → [u i−1 , [u i , y i−1 ] 1 2 ] 2 . Hence, we conclude from the discussion leading to (38) in connection with the chain rule the validity of (40).
It remains to express T 3 , T 4 in terms of Jacobi fields. Concerning T 3 we consider the geodesic γ connecting u i−1 = γ(0) and m = γ(1). T 3 is related to the Jacobi fields along γ as follows: consider those Jacobi fields J along γ for which J(0) = 0 which means that they belong to geodesic variations leaving u i−1 fixed. Then, the adjoint of the mapping equals T 3 . Concerning T 4 , we consider the geodesic ξ connecting y i−1 = ξ(0) and u i = ξ(1). T 4 is related to the Jacobi fields J along ξ for which J(1) = 0 by the mapping Proof of Lemma 4.3. We consider the assertion of Lemma 4.3 for the mapping T 3 . By the proof of Lemma 4.2, we have to determine the adjoint of the mapping given by (77) more explicitely.
Since the covariant derivative of the curvature tensor equals zero in a symmetric space, the differential equation for the Jacobi fields J in terms of the scalar coefficient x of the vector field pt m,[m,ui−1]t v n obtained by parallel transport of the eigenvector v n along γ reads x = −λ n d 2 x, x(0) = 0. Solving this scalar ODE shows (41). Then, f (λ n ) corresponds to the value x(2) of the solution of the scalar ODE of the previous line with additional boundary condition x (1) = 1. Solving this scalar ODE yields (42). We have already derived the corresponding statement for T 4 in [10].
Our next goal is to show Lemma 4.5. To this end, we introduce the mapping F t which is a slight extension of D pt for different parameters as follows. For t ∈ [0, 1], we consider the mapping where, for the definition of F t and in the following Lemma, we use the shorthand notation pt s,t to denote the parallel transport from [u i , u i−1 ] s to [u i , u i−1 ] t and the norm on the right hand side is the one introduced by the Riemannian scalar product in the point [u i , u i−1 ] t . Note that F 0 (u i , u i−1 , y i , y i−1 ) = D pt ([u i , y i ], [u i−1 , y i−1 ]), so in order to obtain the derivative of D pt it suffices to differentiate F 0 w.r.t. its four arguments. The following lemma shows that we can also consider F t instead.
Lemma A.4. The function F t given by (79) is independent of t.
Proof. The Riemannian scalar product is invariant under parallel transport, i.e., for any tangent vectors v, v sitting in the arbitrary point x, and the parallel transport along any curve γ with any points y, z sitting on that curve, we have pt x,y v, pt x,y v y = pt x,z v, pt x,z v z . Hence, for all s, t ∈ [0, 1], pt 0,s log ui y i − pt 1,s log ui−1 y i−1 2 = pt t,s pt 0,t log ui y i − pt t,s pt 1,t log ui−1 y i−1 2 = pt 0,t log ui y i − pt 1,t log ui−1 y i−1 2 .
Since the first expression in (80) equals the square of F s , and the last expression equals the square of F t , this shows that F s = F t for all s, t ∈ [0, 1].
As pointed out below Lemma 4.5, Lemma 4.5 is obtained by specifying F = F 0 in the following lemma.
Lemma A.5. The function F t given by (79) is symmetric with respect to interchanging (u i , y i ) with (u i−1 , y i−1 ) i.e., F t (u i , u i−1 , y i , y i−1 ) = F t (u i−1 , u i , y i−1 , y i ). In particular, for points with F t = 0, the gradient of the function F t w.r.t. the third variable y i , is given by ∇ yi F t (u i , u i−1 , y i , y i−1 ) = ∇ yi−1 F t (u i−1 , u i , y i−1 , y i ), where ∇ yi−1 F t denotes the derivative of F t w.r.t. the fourth argument. Further, again for points with F t = 0, the gradient of the function F t w.r.t. the first component variable u i , is given by ∇ ui F t (u i , u i−1 , y i , y i−1 ) = ∇ ui−1 F t (u i−1 , u i , y i−1 , y i ), where ∇ ui−1 F t denotes the derivative of F t w.r.t. the second argument.
Proof. By the definition of F t in (79), we have F t (u i , u i−1 , y i , y i−1 ) = F 1−t (u i−1 , u i , y i−1 , y i ). By Lemma A.4, F t is independent of t. Together, this implies the identity F t (u i , u i−1 , y i , y i−1 ) = F t (u i−1 , u i , y i−1 , y i ) which is the first assertion of the lemma. The following two assertions on the gradients are then an immediate consequence of this symmetry property.
Proof of Lemma 4.6. We note that F agrees with the function F 0 with F t as in Equation (79) by definition and, consequently, also with F 1 by Lemma A.4. For fixed u i , u i−1 and y i , we decompose the mapping y i−1 → F 1 (u i , u i−1 , y i , y i−1 ) into the mappings G, H, i.e., F 1 = H • G, where G(y i−1 ) = log ui−1 y i−1 (81) locally maps the manifold M to the tangent space at u i−1 , and maps from the tangent space at u i−1 to the positive real numbers. The differential of H (as a map defined on the tangent bundle) is given by d ξ H(η) = ξ−z |ξ−z| , η , for ξ = z. Here ξ is a point in the tangent space at u i−1 , and η is the direction in the tangent space with respect to which the differentiation is performed. Hence, the gradient of H at ξ equals ξ−z |ξ−z| . So the gradient of H at the point ξ = log ui−1 y i−1 equals ∇H(log ui−1 y i−1 ) = log ui−1 y i−1 − z / log ui−1 y i−1 − z .
In order to get the gradient of F 1 , we have to multiply ∇H(log ui−1 y i−1 ) with the adjoint of the differential of G. The mapping G is related to the Jacobi fields along the geodesic γ : [0, 1] → M connecting the points γ(0) = u i−1 and γ(1) = y i−1 as follows. Consider the collection J of Jacobi fields J along γ for which J(0) = 0, which means that they belong to geodesic variations leaving u i−1 fixed. Then the mapping equals the differential of G.
Proof of Lemma 4.7. We use the basis {v n } to express the operator T of Lemma 4.6 (which is given as the adjoint of the derivative of G defined by (81) evaluated at y i−1 ). We use the expression (84) for the derivative of G. Since the covariant derivative of the curvature tensor equals zero in a symmetric space, the differential equation for the Jacobi fields J in terms of the spherical excess A t + B t + C t − π of the triangle ∆ t , where A t is the angle at u i , C t is the angle at u i−1 and B t is the angle at exp ui−1 tv n of the triangle ∆ t . Hence, Here, the first identity is a consequence of the chain rule combined with α 0 = 0. Since sin t t = 1 + o(1), and since C t = π/2 by the orthogonality of v n and log ui−1 u i , we get ω = lim We recall that d = d(u i , u i−1 ) and use the following identities for spherical triangles with an angle of π/2 (cf. [84]) A t = arctan tan t sin d , B t = arctan tan d sin t .
Using the Taylor expansion of the arctan function w.r.t. 0 we obtain that Now, we combine (99) with (100) and conclude, using (97), that ω = 1 sin d − 1 tan d . Together with (96), this shows (58). To see (59), we notice that the connection is linear w.r.t. the direction of differentiation. Therefore, (59) is a consequence of (58) and the expression v n , w equals the coefficient of the corresponding linear combination.
If u i = u i−1 , the mapping in (95) is constant; hence differentiation of this mapping yields zero which implies that the differential D dt | t=0 B n t = 0. This shows the last assertion and completes the proof.
Proof of Lemma 4.11. Since the space of positive matrices is a Riemannian symmetric space representable as quotient of matrix Lie groups there are explicit formulae for the objects of Riemannian geometry such as the exp mapping and the parallel transport available. The corresponding formulae may be found in, e.g., [77]. Our first task is to explicitly express the mapping B n t in the space of positive matrices which form a symmetric space. We use the notation γ t : [0, 1] → M to denote the geodesic starting in u i−1 with direction v, i.e., Here, Exp denotes the ordinary matrix exponential. Then, B n t may be expressed in the space of positive matrices by B n t = pt ui,γt z = u