A global approach to the refinement of manifold data
Abstract
A refinement of manifold data is a computational process, which produces a denser set of discrete data from a given one. Such refinements are closely related to multiresolution representations of manifold data by pyramid transforms, and approximation of manifoldvalued functions by repeated refinements schemes. Most refinement methods compute each refined element separately, independently of the computations of the other elements. Here we propose a global method which computes all the refined elements simultaneously, using geodesic averages. We analyse repeated refinements schemes based on this global approach, and derive conditions guaranteeing strong convergence.
Key Words. Manifold data, geodesic average, convergence analysis.
AMS(MOS) subject classification. 65D99, 40A99, 58E10.
1 Introduction
In recent years many modern sensing devices produce data on manifolds or data that is modelled as points on a manifold. An example of such data is orientations of a rigid body as function of time, which can be regarded as data sampled from a function mapping a real interval to the Lie group of orthogonal matrices [29]. The classical methods for the approximation of a function from its samples, such as polynomial or spline interpolation, are linear, and there is no guarantee that such approximations produce always manifold values, due to the nonlinearity of manifolds. Therefore, alternative methods are required.
Contrary to the development of classical approximation methods and numerical analysis methods for realvalued functions, the development in the case of manifoldvalued functions, which is rather recent, was mainly concerned in its first stages with advanced numerical and approximation processes. Examples of such processes are geometric integration of ODE on manifolds (see e.g. [19]), subdivision schemes on manifolds (see e.g. [34, 37]) and waveletstype approximation on manifolds (see e.g. [17, 29]).
Subdivision schemes were created originally to design geometrical models [3, 23]. Later, they were recognized as methods for approximation [5, 11]. The important advantage of these schemes is their simplicity and locality. They are defined by repeatedly refining sequences of points, applying in each refinement step simple and local arithmetic averaging. This enables the extension of subdivision schemes to more abstract settings, such as matrices [32] and sets [9].
For manifold valued data, Wallner and Dyn [36] introduced the concept of adapting linear subdivision schemes to manifold data, and in particular for Lie group data. That paper initiated a new path of research on subdivision schemes for manifold data, e.g., [32, 34]. Adaptation of a linear subdivision scheme can be done in several ways, for example, by rewriting the refinement rules as repeated binary averages, and then replacing each binary average by a geodesic average, see e.g., [32, 36].
Averages play a significant role in the methods for the adaptation of linear subdivision schemes to manifold data. A natural choice of an average of two points on a geodesically complete manifold is the midpoint of the geodesic curve between the two points. In some cases, the geodesic curve is known explicitly, e.g., [14, 16, 18, 25], while in general it can be calculated numerically, e.g., [4, 15, 22, 26].
The weighted geodesic average is induced by the geodesic curve, and acts as a generalization of the weighted arithmetic average in Euclidean spaces. For a weight , it is the point on the geodesic curve, connecting the two averaged points, which divides this curve segment in the ratio . Furthermore, on several manifolds, the geodesic average can also be extended to weights outside , that is extrapolating the geodesic curve of two points beyond these points, e.g., [20]. The geodesic average is also welldefined on more general spaces known as geodesic metric spaces, e.g., [1]. Thus, in such spaces our adaptation method is also valid.
We present here a method for the adaptation of linear subdivision schemes to manifold data based on the idea of replacing weighted arithmetic averages by weighted geodesic averages in a generalized LaneRiesenfeld algorithm [23]. The refinement step in this proposed generalization consists of an elementary refinement of doubling the data, followed by several rounds of averaging. In each round of averaging the data is replaced by the same weighted average of all pairs of adjacent points in the data. Such an adaptation is discussed shortly in [8, 36]. We term such a refinement step “global refinement”.
Many results, concerning the convergence and smoothness of adapted subdivision schemes, are presented in the literature of the past few years, e.g., [34, 36, 37]. Most of these results are based on proximity conditions. A proximity condition bounds the distance between the operation of an adapted refinement step to the operation of its linear counterpart in terms of the maximal distance between adjacent data points. Such proximity conditions hold, since a manifold is locally close to a Euclidean space. Thus, the convergence results are often valid only for “dense enough data”, which is, in general, a condition that is hard to quantify and depends on properties of the manifold (such as curvature).
Recently, a progress in the convergence analysis is established in several papers which address the question of convergence from any initial data. Such a result is presented in [13] for adapted subdivision schemes to data in Hadamard spaces. Results for data on the manifold of positive definite matrices are derived in [32]. For the case of interpolatory subdivision schemes there are also results for several different metric spaces e.g., [20, 21, 35].
Here we prove convergence from all initial data, of the above adapted generalized LaneRiesenfeld algorithm, when the weighted average in each round corresponds to a weight in , and give conditions for such convergence when some averages have weights outside . In addition, we extend the above construction to a wider class of linear schemes, by introducing weighted trinary averages based on geodesic weighted averages, and give sufficient conditions for convergence from all initial data. In all these cases, and for manifolds with globally bounded curvature, the convergence guarantees that the limits are , based on the proximity analysis in [36].
Three important observations on our adaptation method:

It extends the class of linear schemes for which an adapted scheme is known to be convergent from all initial data.

It is welldefined and convergent from all data in a wide class of geodesic metric spaces.

It leads to computationally feasible subdivision schemes.
The convergence analysis introduced in this paper supplies a new tool for the analysis of linear schemes. In particular, this analysis guarantees the convergence of any linear scheme with a symbol which is a Hurwitz polynomial, up to multiplication by a monomial. The question whether this method can improve our ability to determine the convergence of linear subdivision schemes is beyond the scope of this paper and is still under investigation.
The paper is organized as follows. We start in Section 2 by providing a short survey of the required background, including a summary on the LaneRiesenfeld algorithm and a short review on geodesics and manifolds. We conclude Section 2 with a short discussion on a sufficient condition for the convergence of adapted subdivision schemes. Section 3 introduces our generalization of the LaneRiesenfeld algorithm. Then, we give conditions for the convergence of an adapted scheme based on this algorithm, from any initial manifold data, where the corresponding linear scheme has a factorizable symbol over the reals. In Section 4 we further extend the algorithm to the adaptation of general linear schemes, and conclude the paper by the convergence analysis of these schemes.
2 Preliminaries
2.1 Subdivision schemes and the LaneRiesenfeld algorithm
Linear, univariate subdivision schemes are defined on numbers (the functional setting) , and are extended to vectors by operating on each component separately. In the functional setting, these schemes are approximation operators, when the data is sampled uniformly from a continuous function . We denote the sampled data , , by . Any subdivision scheme consists of refinement rules that map to a new sequence associated with the values at , .
Let us denote by a refinement rule, defined by a finitely supported mask , as
(1) 
A (stationary) subdivision scheme with a refinement rule is a repeated application of (1) and is also denoted by .
A subdivision is termed convergent if the sequence of piecewise linear interpolants to the data converges uniformly (see e.g. [7]). By definition, the limit is a continuous function.
The LaneReisenfeld (LR) algorithm is a classical algorithm, which executes the refinement rules of a Bspline subdivision scheme [23]. This algorithm replaces each step of refinement by an elementary refinement (doubling all the data points) followed by several stages of averaging. In each stage of averaging, the data points are replaced by the midpoints of all pairs of consecutive data points. As a result, the refinement is done simultaneously to all data points. We term this refinement a global refinement, in contrary to the direct evaluation of (1), where each refined point is calculated independently of the other refined points. The refinement step of the LR algorithm is presented in Algorithm 1.
An important tool in the analysis of convergence and smoothness of subdivision schemes is the symbol, defined as the transform of the mask , that is . For example, the symbol of the Bspline subdivision scheme of degree is . A necessary condition for convergence is and implying that the subdivision scheme is invariant to a translation of the data [7, Proposition 2.1]. With the symbol the refinement rules (1) can be written algebraically as
(2) 
where the equality is in the sense of equal coefficients corresponding to the same power of . The LR algorithm is an interpretation of (2) with the symbols of the Bspline subdivision schemes. For explanation see Section 3.1 and in particular (10).
Over the years, several generalizations of the LR algorithm have been proposed. In [2] any step of the subdivision consists of a refinement step of a fixed converging subdivision scheme, followed by a fixed number of “smoothing rounds” based on another subdivision scheme (e.g., applying the insertion rule of an interpolatory scheme to each point). In [10, 31] nonlinear averages of numbers replace the arithmetic (linear) averages. A generalization based on a geodesic average goes back to [27, 28] where a corner cutting subdivision scheme based on geodesic averages is presented and analysed. In [9] the LR algorithm is adapted to compact sets based on the metric average which is a geodesic average in the metric space of compact sets with the Hausdorff metric.
In this paper we discuss the adaptation of subdivision schemes from numbers to manifold data. To distinguish between sequences of numbers (or vectors) to sequences on a manifold, we denote by and a sequence of Euclidean data and manifold, respectively.
2.2 On manifolds and geodesics
A geodesic (or a geodesic curve) is a fundamental notion in differential geometry. This notion is an extension of the shortest arc on a surface, joining two arbitrary points and on the surface. On a plane, the geodesic is simply the line segment connecting and , described by
(3) 
This line can be also characterized by its zero curvature and its endpoints. For a manifold, this property is generalized by having zero geodesic curvature (or constant velocity derived from the first fundamental form). In Riemannian manifolds, the geodesic curve is defined as the solution to the geodesic EulerLagrange equations. It turns out that any shortest path between two points is a geodesic curve.
In connected Riemannian manifolds, the HopfRinow theorem guarantees that geodesic curves connecting any two points are globally well defined and smooth, see e.g., [6]. Such manifolds are also known as geodesically complete or simply complete Riemannian manifolds. For such manifolds, one can derive the uniqueness of the geodesic curve connecting any two points, in case one point is outside the cut locus of the other. Henceforth, we will use the term geodesic curve for such shortest path curves.
The geodesic curve is of great importance in our adaptation procedures. A natural question is its availability in different manifolds. Indeed, in many cases, the geodesic curve is known explicitly. Here are several examples: on a sphere (e.g., [14]), on an ellipsoid (e.g., [16]), on the cone of positive definite matrices (e.g., [18]), in the Lie group of orthogonal matrices of the same determinant (e.g., [33, Chapter 3]), in the Heisenberg groups (e.g., [25]). Alternatively, geodesics can be calculated numerically. This can be done by directly solving the EulerLagrange equations (e.g., [15]), by fast marching methods (e.g., [22]), by exploiting heat kernels based methods (e.g., [4]), or other hypersurfaces techniques (e.g., [26]), just to name a few.
An important property of the geodesic curve is the metric property. Let be a complete Riemannian manifold with associated metric . Then, for any the geodesic curve connecting and , that is , with and , satisfies
(4) 
Since is a metric, we also have the compliment formula . In this paper, we consider data such that the geodesic curve between any two adjacent data points in is welldefined, and term such data “admissible”. Then, the geodesic curve is used as a weighted mean, that is the manifold analogue of the arithmetic mean (3). In some cases, we may need to be defined for values of outside , but close to it. Therefore, we must assume that the geodesic curve is welldefined for these “extrapolation” values. In these cases the metric property (4) is modified, replacing by .
There are some nonlinear spaces, other than Riemannian manifolds, where the geodesic curve connecting any two points is unique. These are the geodesic metric spaces, see e.g., [1]. In such spaces, the differential structure is missing and a geodesic curve is defined as the path satisfying (4). Clearly, this definition agrees with the geodesic curve on Riemannian manifolds. Note that, in general, we do not need the uniqueness of the geodesic curve, but a canonical way to choose it, see e.g., [9].
2.3 Sufficient conditions for convergence of manifoldvalued subdivision schemes
The convergence of manifoldvalued subdivision schemes can be defined intrinsically. For that, we defined for any data sequence , a piecewise geodesic interpolant , connecting any pair of consecutive points in by their geodesic curve. The manifoldvalued subdivision scheme is convergent, if the sequence , converges uniformly relative to the metric of the manifold (see [12]).
The analysis of adapted subdivision schemes in many papers is based on the method of proximity, introduced in [36]. This analysis uses conditions that indicate the proximity of the adapted refinement rule to its corresponding linear refinement rule . The simplest proximity condition is
(5) 
In [36] it is proved that if is a refinement rule of a convergent scheme that generates limits, then condition (5) implies (with additional mild assumptions on the refinement rule ) that for small enough, the adapted subdivision scheme , applied to the initial data , converges to a limit.
The weakness of the proximity method is that convergence is only guaranteed for “close enough” data points. This requirement is typically not easy to quantify and it depends on the manifold and its curvature.
For a linear subdivision schemes a contractivity factor , namely
(6) 
implies the convergence of the scheme from any initial data, see e.g. [7].
For nonlinear subdivision schemes, and in particular for schemes adapted to manifold data, contractivity is not sufficient for convergence, and an additional condition is required, see [12].
Definition 2.1 (Displacementsafe).
Let be a subdivision scheme adapted to manifold data. We say that is “displacementsafe” if
(7) 
for any sequence of manifold data , where is a constant independent of .
In [12], it is proved that
Theorem 2.2.
Let be a displacementsafe subdivision scheme for manifold data with a contractivity factor . Then, is convergent for any input manifold data.
Remark 2.3.
Two concluding remarks:

Note that interpolatory schemes satisfy (7) with by definition and thus are displacementsafe.

In [36] it is proved that any adaptation of (1) based on repeated geodesic averages satisfies (5), under mild assumptions on the manifold, such as manifolds with globally bounded curvature. This observation implies that for with , (7) is also satisfied. Thus, for such schemes, it is enough to show that the scheme has a contractivity factor to obtain convergence for any initial data and to conclude that the limit is .
3 Adaptation of generalized LR algorithms
We present an adaptation method of generalized LR algorithms, based on geodesic averages. This method is already introduced in [8, 36]. Nevertheless, the convergence result stated there is the one that follows from proximity conditions, which applies only for small enough. First, we discuss in detail our adaptation and then analyze the resulting schemes, charactering classes of schemes for which convergence from any initial data is guaranteed.
3.1 The algorithm of global refinement
Consider a linear subdivision scheme of the form (1), with a symbol . The factorization of the symbol plays an important role in the analysis of convergence and smoothness of linear subdivision schemes [7], and is also significant in our adaptation.
We start with a class of convergent linear subdivision schemes having symbols which can be factorized into real linear factors. Recall that a necessary condition for convergence is that and [7, Proposition 2.1]. Thus, we can write
(8) 
where are the nonzero roots of the symbol and is an integer. Note that cannot be a root of a symbol since . Thus, . and (8) is welldefined. We further define to be the minimizer of
(9) 
among . The reason will become clear later.
The relation between the factorization (8) and the global refinement is based on (2). For the symbol (8) we get from (2) that the linear scheme can be interpreted as
(10) 
By this interpretation, the factor indicates an initial elementary refinement step in which the data is duplicated. Then, each of the factors , implies a step of averaging, in which the current data is replaced by the weighted averages with weights on its pairs of adjacent points. A zero root of the symbol merely changes the value of . This value determines the shift of indices required to be applied, at the end of each refinement step. Note that for , , this interpretation becomes the LR algorithm. Thus, we consider the global refinement step corresponding to (8) a generalized LR algorithm.
The adaptation of the global refinement, based on geodesic averages, is summarized in Algorithm 2.
Note that for data sampled from a geodesic curve, all points generated by Algorithm 2, are on this geodesic curve.
3.2 Analysis of schemes corresponding to factorizable symbols over the reals
For our first result, we restrict the discussion to the case where the symbol (8) has a full set of real negative roots, namely , .
Theorem 3.1.
Proof.
Note that the contractivity factor of Theorem 3.1 satisfies since and , with for .
The LR algorithm satisfies the conditions of Theorem 3.1. Indeed, this theorem is a generalization of a similar result in [9, Lemma 4.1] for the adapted LR algorithm to compact sets.
Next, we show that the adapted subdivision schemes corresponding to symbols having a full set of real negative roots, are displacementsafe.
Theorem 3.2.
Proof.
The proof shows by induction that , . Denote by the linear subdivision scheme with a symbol obtained from the symbol of by retaining the first factors, , so that the adapted scheme of , , uses only iterations of the loop of Line 3 in Algorithm 2. Obviously . We use induction on . For , after the initial steps of Lines 1 and 2, Algorithm 2 inserts new points on the geodesic curves, connecting adjacent data points. Therefore, it is clear that we have , namely we get the constant for the case . The induction step assumes
for a given , with a constant , which depends on and is independent of . Then, using the triangle inequality we get
While by the metric property (4) (see Line 5 in Algorithn 2)
(11) 
Since Theorem 3.1 implies that
(12) 
we can choose and the proof follows. The shift, defined by in (8) and done in Line 9 of Algorithm 2, does not affect the above bound, since is the same for all . ∎
We conclude
Corollary 3.3.
The second case analyzed here corresponds to symbols of the form (8) with several positive roots. Positive roots mean negative weights in the averages, namely extrapolating averages in Line 5 of Algorithm 2.
Theorem 3.4.
Let be a linear convergent subdivision scheme with symbol of the form (8), such that has at least one negative root in addition to the root . Define
and renumerate the factors in (8) such that is attained at . If
(13) 
where
then the adapted scheme based on global refinement has a contractivity factor , and it converges from any admissible initial data on the manifold.
Proof.
The proof basically modifies the proofs of Theorem 3.1 and Theorem 3.2. By assumption the set is not empty, and therefore . Similarly to the proof of Theorem 3.1 the application of an averaging step in Line 5 of Algorithm 2, corresponding to , does not expand the bound on the distances between consecutive points in the data. On the other hand, an averaging step corresponding to expands the bound.
To obtain the expanding factor note that after the th step in Line 5 of Algorithm 2 we can bound the distance between consecutive points by
(14) 
Defining , , we obtain from (14)
(15) 
This together with assumption (13) shows that is a contractivity factor of the adapted scheme.
To complete the convergence proof, we observe that since , assumption (13) implies that , . Modifying the proof of Theorem 3.2, we get in its notation that (11) is replaced by
Using the same inductive argument, and the bound (15), we get
Thus, in this case . By (13) , and since implies , we finally arrive at .
Remark 3.5.
Two remarks for section 3.2:

As is proved in Theorems 3.1 and 3.2 the adaptation of Algorithm 2 leads to converging subdivision schemes when applied to linear subdivision schemes with positive mask coefficients, such that their symbols have a full set of negative roots. Theorem 3.4 extends the convergence to schemes with symbols having few positive roots in addition to at least two negative ones, which may correspond to masks with some negative coefficients.

Negative coefficients necessarily appear in the masks of smooth interpolatory schemes. However, the adaptation based on global refinement is inappropriate for interpolatory subdivision schemes, since the adapted schemes are not interpolatory any more. The commutativity of multiplication of numbers guarantees that for numbers the local refinement and the global refinement coincide.
In the next section we show that the global refinement can be interpreted as local refinements, based on a “pyramid averaging”.
3.3 interpretation of the global refinement as local refinement
Most known adaptation methods of convergent linear subdivision schemes to manifold data are based on first rewriting the average (1) in terms of repeated binary averages, and then replacing the linear averages by some manifold averages, see e.g. [34, 36, 37]. We term the so obtained refinement rules “local refinement”.
Next we show that global refinement can be interpreted as local refinement based on geodesic averages. This observation together with 2 of Remark 2.3 leads to the conclusion that the convergence of schemes adapted by global refinement guarantees limits.
We now describe how the global refinement can be interpreted as local refinement. For even, in Algorithm 2 can be calculated by a series of repeated averaging operating on . First we replace by , . We take from this sequence the first points, to form the initial level for a “pyramid averaging” of levels. In the th level of the pyramid averaging any pair of adjacent points is replaced by its geodesic average with weight , . Thus at the th level there are points. is the only value obtained at level of the pyramid averaging.
For odd, in Algorithm 2 can be calculated similarly, starting the same pyramid averaging from a different sequence. This sequence is obtained from by first replacing by , and then taking the first points. For illustrations and explanation of the pyramid averaging notion see [30].
The global refinement calculates only once each geodesic averages of adjacent points in the data, while the same average appears in the calculation of several points by local refinement. Thus, the global refinement is more efficient in terms of computational operations as compared to its local refinement interpretation. Note that it is possible to define a scheme adapted by local refinement which uses the same number of geodesic averages as the global refinement [12].
4 Adaptation based on global refinement – the general case
We extend the global refinement algorithm to converging linear schemes with general symbols. Then, instead of (8) such symbols, which are real polynomials, can be factorized into real linear factors (in addition to ) and quadratic real factors, with . Any complex root of the symbol corresponds to a real quadratic irreducible factor over the reals of the form
(16) 
where and is the real part of . The average associated with such a factor has, in the sense of the global refinement algorithm, the following weights
(17) 
Note that . Instead of (8) we have in this case the factorization
(18) 
Lemma 4.1.
For any complex ,
(19) 
Proof.
4.1 The general algorithm of global refinement
For an irreducible quadratic factor in (18) one is required to average points on the manifold at once. Motivated by the pyramid averaging of Section 3.3, we define such an average and term it a three pyramid.
Definition 4.2.
For three points with corresponding weights , the “three pyramid” is
where the following constraints must hold

.

.

.
Remark 4.3.
The global refinement of Algorithm 2 uses uniform averaging in each level. The following lemma shows that this is not possible for symbols with complex roots.
Lemma 4.4.
Proof.
4.2 Optimal choice of parameters in the three pyramid
To optimally bound the distance
(22) 
we start by setting . The reasons for this choice are presented in details in Appendix A.1. For the other parameters, we first prove the following Lemma.
Proof.
Theorem 4.6.
Proof.
Figure 1 accompanies the proof. There and correspond to and respectively, while and correspond to and respectively. and there correspond to and respectively.
We first apply the metric property (4) and the triangle inequality to get (see the schematic illustration in Figure a)
(24) 
Note that and that . Similarly we get
and since by Lemma 4.5, we conclude that is closer to than (see Figure b). Observing that these two averages lie on the geodesic curve connecting and , we conclude that
(25) 
To prove (23) we sum the following three bounds, on the lengths of the three parts of the path connecting to via and in Figure c,
The first and third bounds are obtained from Definition 4.2 by (4) and (24), the second bound is (25).
∎
Remark 4.7.
Two important conclusions, related to the parameters of the three pyramid:
Note that in the linear case, any averaging step corresponding to a complex root does not expand the distance between consecutive points as long as the weights (17) are positive, that is the real part of is positive.
4.3 Analysis of convergence
First, we consider the case of symbols of the form (18) having several complex roots and then discuss in detail the case of a single complex root.
In case of positive roots, which is analysed in Theorem 3.4, we show an initial contractivity factor induced by , associated with the negative root, followed by a series of expanding factors for , associated with the positive roots. Equipped with Theorem 4.6, the analysis of the convergence of the schemes adapted by Algorithm 3 is essentially the same.