Geometric subdivision and multiscale transforms

TL;DR

This paper explores how to perform subdivision and multiscale transforms on data with geometric structures, ensuring operations respect data symmetries and intrinsic properties.

Contribution

It provides a framework for geometric subdivision and multiscale transforms on various structured spaces, emphasizing intrinsic operations and convergence properties.

Findings

Analysis of geometric structures like metric spaces and manifolds.

Discussion of intrinsic averaging and refinement operations.

Current understanding of convergence and smoothness in multiscale transforms.

Abstract

Any procedure applied to data, and any quantity derived from data, is required to respect the nature and symmetries of the data. This axiom applies to refinement procedures and multiresolution transforms as well as to more basic operations like averages. This chapter discusses different kinds of geometric structures like metric spaces, Riemannian manifolds, and groups, and in what way we can make elementary operations geometrically meaningful. A nice example of this is the Riemannian metric naturally associated with the space of positive definite matrices and the intrinsic operations on positive definite matrices derived from it. We disucss averages first and then proceed to refinement operations (subdivision) and multiscale transforms. In particular, we report on the current knowledge as regards convergence and smoothness.

Equations

$$x={\mathrm{avg}}_{j\in \mathbb{Z}}(a_j,x_j):=\sum a_j{x}_j,\text{where}\sum a_j=1.$$

$$x=\mathrm{avg}(a_j,x_j)$$

$$⟺x=y+\sum a_j(x_j-y)\text{for any}y$$

$$⟺x\text{minimizes}\sum a_j\Vert x-x_j{\Vert }^2.$$

$${\mathrm{avg}}_F(a_j,x_j):=\arg \underset{x\in \mathcal{M}}{\min }\sum a_j{\text{d}}_{\mathcal{M}}(x,x_j{)}^2.$$

$${\mathrm{avg}}_F(\mu )=\arg \underset{x\in \mathcal{M}}{\min }\int {\text{d}}_{\mathcal{M}}(x,\cdot {)}^{2}d\mu .$$

$$p\oplus v$$

$$\nabla {\text{d}}_{\mathcal{M}}(\cdot ,a)(x)$$

$$v\in \mathfrak{g}⟹e\oplus v=\exp (v)={\sum }_{k\ge 0}\frac{1}{k!}{v}^k.$$

$$p\oplus v$$

$$\displaystyle d(a,b)=\|\log(a^{-1/2}ba^{-1/2})\|_{2}=\Big{(}\sum\nolimits_{\lambda_{1},\ldots,\lambda_{n}\in\sigma(a^{-1}b)}\log^{2}\lambda_{j}\Big{)}^{1/2}.$$

$$x={\mathrm{avg}}_E(a_j,x_j):⟺\sum a_j(x_j\ominus x)=0.$$

$${\mathrm{avg}}_y(a_j;x_j):=y\oplus \sum a_j(x_j\ominus y).$$

$$y={\mathrm{avg}}_F(a_j;x_j)⟹{\mathrm{avg}}_y(a_j;x_j)=y\oplus \sum a_j(x_j\ominus y)=y\oplus 0=y,$$

$$L^{N}S=SL.$$

$$p,\text{Sp},S^{2}p,S^{3}p,\dots ,$$

$${\text{Sp}}_i={\sum }_j{a}_{i-Nj}{p}_j.$$

$${\text{Sp}}_i=\mathrm{avg}(a_{i-Nj};p_j).$$

$${\text{Tp}}_i={\mathrm{avg}}_{m_i}(a_{i-Nj};p_j),$$

$${\text{Tp}}_i=\pi ({\text{Sp}}_i).$$

$$\pi :{\text{GL}}_n\to {\text{O}}_n,\pi (g)=(g{g}^T{)}^{-1/2}g.$$

$${\text{Sp}}_{2i}$$

$${\text{Sp}}_{2i+1}$$

$$T$$

$$\delta (p)=\underset{|i-j|\le 1}{\sup }{\text{d}}_{\mathcal{M}}(p_i,p_j),$$

$$\delta (\text{Sp})\le \gamma \delta (p),\text{for some}\gamma <1,\text{resp.}{\sup }_{i\in {\mathbb{Z}}^s}{\text{d}}_{\mathcal{M}}({\text{Sp}}_{Ni},p_i)\le \lambda \delta (p).$$

$${\text{d}}_{\mathcal{M}}(f_k(\xi ),f_{k+1}(\xi ))\le \lambda \delta (S^{k}p)\le \lambda {\gamma }^k\delta (p).$$

$${\text{d}}_{\mathcal{M}}(f(\xi ),f(\eta ))$$

$$\le \sum {\text{d}}_{\mathcal{M}}(f_j(a_j),f_{j+1}(a_j))+{\text{d}}_{\mathcal{M}}(f_{j+1})(a_j),f_{j+1}(a_{j+1}))+(\ast ).$$

$${\text{d}}_{\mathcal{M}}(f(\xi ),f(\eta ))$$

$$\le 2(\lambda +s\gamma )\delta (p){\sum }_{j\ge k}{\gamma }^j\le C\delta (p)\frac{{\gamma }^k}{1-\gamma }.$$

$${\text{d}}_{\mathcal{M}}(f(\xi ),f(\eta ))\le C^{\prime }{\gamma }^{-\log |\xi -\eta |/\log N}=C^{\prime }|\xi -\eta {|}^{-\log \gamma /\log N}.$$

$$a^{\ast }(z)\Delta p(z^N)$$

Full text

1. Computing averages in nonlinear geometries

The line of research presented in this chapter was first suggested by a 2001 presentation by D. Donoho on multiscale representations of discrete data [11]. A subsequent Ph.D. thesis and accompanying publication appeared a few years later [48]. Multiscale representations are intimately connected with refinement procedures (prediction operators). These are in themselves an interesting topic with applications, e.g. in computer graphics. Iterative refinement a.k.a. subdivision in turn is based on the notion of average. Consequently this chapter is structured into the following parts: Firstly a discussion of averages, in particular averages in metric spaces and in manifolds. Secondly, subdivision rules and the limits generated by them. Thirdly, multiresolution representations.

We start with the affine average w.r.t. weights $a_{j}$ of data points $x_{j}$ contained in a vector space. It is defined by

[TABLE]

In this chapter we stick to finite averages, but we allow negative coefficients. For data whose geometry is not that of a vector space, but that of a surface contained in some Euclidean space, or that of a group, or that of a Riemannian manifold, this affine average often does not make sense. In any case it is not natural. Examples of such data are, for instance, unit vectors, positions of a rigid body in space, or the 3 by 3 symmetric positive definite matrices which occur in diffusion-tensor MRI. In the following paragraphs we show how to extend the notation of affine average to nonlinear situations in a systematic way. We start by pointing out equivalent characterizations of the affine average:

[TABLE]

The Fréchet mean

Each of (2)–(4) has been used to generalize the notion of weighted average to nonlinear geometries. Some of these generalizations are conceptually straightforward. For example, Equation (4) has an analogue in any metric space $({\mathord{\mathcal{M}}},\operatorname{\text{\it d}_{\!{\mathord{\mathcal{M}}}}})$, namely the weighted Fréchet mean defined by

[TABLE]

It is a classical result that in case of nonnegative weights, the Fréchet mean exists and is unique, if ${\mathord{\mathcal{M}}}$ is a Hadamard metric space. This property means ${\mathord{\mathcal{M}}}$ is complete, midpoints exist uniquely, and triangles are slim, cf. [1].111More precisely, for all $a,b\in{\mathord{\mathcal{M}}}$ there is a unique midpoint $x=m(a,b)$ defined by $\operatorname{\text{\it d}_{\!{\mathord{\mathcal{M}}}}}(x,a)=\operatorname{\text{\it d}_{\!{\mathord{\mathcal{M}}}}}(x,b)=\operatorname{\text{\it d}_{\!{\mathord{\mathcal{M}}}}}(a,b)/2$, and for any $a,b,c\in{\mathord{\mathcal{M}}}$ and points $a^{\prime},b^{\prime},c^{\prime}\in{\mathbb{R}}^{2}$ which have the same pairwise distances as $a,b,c$, the inequality $\operatorname{\text{\it d}_{\!{\mathord{\mathcal{M}}}}}(c,m(a,b))\leq d_{{\mathbb{R}}^{2}}(c^{\prime},m(a^{\prime},b^{\prime}))$ holds.

The Fréchet mean in Riemannian manifolds. In a surface resp. Riemannian manifold ${\mathord{\mathcal{M}}}$, the Fréchet mean locally exists uniquely. A main reference here is the paper [39] by H. Karcher. He considered the more general situation that $\mu$ is a probability measure on ${\mathord{\mathcal{M}}}$, where the mean is defined by

[TABLE]

In this chapter we stick to the elementary case of finite averages with possibly negative weights. The Fréchet mean exists uniquely if the manifold is Hadamard – this property is usually called “Cartan-Hadamard” and is characterized by completeness, simple connectedness, and nonpositive sectional curvature. For unique existence of $\operatorname{avg}_{F}$, we do not even have to require that weights are nonnegative [37, Th. 6].

The Fréchet mean in the non-unique case. If the Cartan-Hadamard property is not fulfilled, the Fréchet mean does not have to exist at all, e.g. if the manifold is not complete (cutting a hole in ${\mathord{\mathcal{M}}}$ exactly where the mean should be makes it nonexistent). If the manifold is complete, the mean exists, but possibly is not unique.

If ${\mathord{\mathcal{M}}}$ is complete with nonpositive sectional curvature, but is not simply connected, there are situations where a unique Fréchet mean of given data points can still be defined, e.g. if the data are connected by a path $c\colon[a,b]\to{\mathord{\mathcal{M}}}$ with $c(t_{j})=x_{j}$. This will be the case e.g. if data represent a time series. Existence or maybe even canonical existence of such a path depends on the particular application. We then consider the simply connected covering $\widetilde{\mathord{\mathcal{M}}}$, find a lifting $\widetilde{c}\colon I\to\widetilde{\mathord{\mathcal{M}}}$ of $c$, compute the Fréchet mean $\operatorname{avg}_{F}(a_{j},\widetilde{c}(t_{j}))$, and project it back to ${\mathord{\mathcal{M}}}$. This average does not only depend on the data points and the weights, but also on the homotopy class of $c$. In fact instead of a path, any mapping $c\colon I\to{\mathord{\mathcal{M}}}$ can be used for such purposes as long as its domain $I$ is simply connected [37].

Finally, if ${\mathord{\mathcal{M}}}$ is complete but has positive sectional curvatures, a unique Fréchet mean is only defined locally. The size of neighbourhoods where uniqueness happens has been discussed by [15, 16, 34]. This work plays a role in investigating convergence of subdivision rules in Riemannian manifolds, see Section 2.4.

The exponential mapping

From the different expressions for the affine average, (2) and (3) seem to be specific to linear spaces, because they involve the $+$ and $-$ operations. However, it turns out that there is a big class of nonlinear geometries where natural analogues $\oplus$ and $\ominus$ of these operations exist, namely the exponential mapping and its inverse. We discuss this construction in surfaces resp. Riemannian manifolds, in groups, and in symmetric spaces.

The exponential mapping in Riemannian geometry. In a Riemannian manifold ${\mathord{\mathcal{M}}}$, for any $p\in{\mathord{\mathcal{M}}}$ and tangent vector $v\in\mathord{T\hskip-1.72218pt_{p}\hskip-0.43057pt{\mathord{\mathcal{M}}}}$, the point $\exp_{p}(v)$ is the endpoint of the geodesic curve $c(t)$ which starts in $p$, has initial tangent vector $v$, and whose length equals $\|v\|$. We let

[TABLE]

One property of the exponential mapping is the fact that curves of the form $t\mapsto p\oplus tv$ are shortest paths with initial tangent vector $v$. The mapping $v\mapsto p\oplus v$ is a diffeomorphism locally around $v=0$. Its differential equals the identity.

Properties of the Riemannian exponential mapping. For complete Riemannian manifolds, $p\oplus v$ is always well defined. Also $q\ominus p$ exists by the Hopf-Rinow theorem, but it does not have to be unique. Uniqueness happens if $\operatorname{\text{\it d}_{\!{\mathord{\mathcal{M}}}}}(p,q)$ does not exceed the injectivity radius $\rho_{\text{inj}}(p)$ of $p$. In Cartan-Hadamard manifolds, injectivity radii are infinite and the exponential mapping does not decrease distances, i.e., $\operatorname{\text{\it d}_{\!{\mathord{\mathcal{M}}}}}(p\oplus v,p\oplus w)\geq\|v-w\|_{\mathord{T\hskip-1.20552pt_{p}\hskip-0.3014pt{\mathord{\mathcal{M}}}}}$. The injectivity radius can be small for topological reasons (e.g.  a cylinder of small radius which is intrinsically flat, can have arbitrarily small injectivity radius), but even in the simply connected case, one cannot expect $\rho_{\text{inj}}$ to exceed $\pi K^{-1/2}$, if $K$ is a positive upper bound for sectional curvatures.

Further, the $\ominus$ operation and the Riemannian distance are related by

[TABLE]

if $v=a\ominus x$ refers to the smallest solution $v$ of $x\oplus a=v$. For more properties of the exponential mapping we refer to [39] and to differential geometry textbooks like [9].

The exponential mapping in groups. In Lie groups, which we describe only in the case of a matrix group $G$, a canonical exponential mapping is defined: With the notation ${\mathord{\mathfrak{g}}}=T_{e}G$ for the tangent space in the identity element, we let

[TABLE]

The curve $t\mapsto e\oplus tv$ is the unique one-parameter subgroup of $G$ whose tangent vector at $t=0$ is the vector $v\in{\mathord{\mathfrak{g}}}$. Again, $v\mapsto e\oplus v$ is locally a diffeomorphism whose differential is the identity mapping.

An inverse log of exp is defined locally around $e$. Transferring the definition of $\oplus$ to the entire group by left translation, the defining relation $g\oplus gv:=g(e\oplus v)$ yields

[TABLE]

Addition is always globally well defined, but the difference $q\ominus p$ might not exist. For example, in $\mathord{{\text{\rm GL}}}_{n}$, the mapping $v\mapsto e\oplus v$ is not onto. The difference exists always, but not uniquely, in compact groups. See e.g. [2].

The exponential mapping in symmetric spaces. Symmetric spaces have the form $G/H$, where $H$ is a Lie subgroup of $G$. There are several definitions which are not entirely equivalent. We use the one that the tangent spaces ${\mathord{\mathfrak{g}}}=T_{e}G$, ${\mathord{\mathfrak{h}}}=T_{e}H$ obey the condition that ${\mathord{\mathfrak{h}}}$ is the $+1$ eigenspace of an involutive Lie algebra automorphism $\sigma$ of ${\mathord{\mathfrak{g}}}$.222i.e., $\sigma$ obeys the law $\sigma([v,w])=[\sigma(v),\sigma(w)]$, where in the matrix group case, the Lie bracket operation is given by $[v,w]=vw-wv$. The tangent space ${\mathord{\mathfrak{g}}}/{\mathord{\mathfrak{h}}}$ of $G/H$ in the point $eH\in G/H$ is naturally identified with the $-1$ eigenspace ${\mathord{\mathfrak{s}}}$ of the involution, and is transported to all points of $G/H$ by left translation. The exponential mapping in $G$ is projected onto $G/H$ in the canonical way and yields the exponential mapping in the symmetric space.

We do not go into more details but refer to the comprehensive classic [35] instead. Many examples of well-known manifolds fall into this category, e.g. the sphere $S^{n}$, hyperbolic space $H^{n}$, and the Grassmannians. We give an important example:

• Example 1.

The Riemannian symmetric space of positive-definite matrices. The space $\mathord{{\text{\rm Pos}}}_{n}$ of positive definite $n\times n$ matrices is made a metric space by letting

[TABLE]

Here $\|\cdot\|_{2}$ means the Frobenius norm, and $\sigma(m)$ means the eigenvalues of a matrix.

The metric (7) is actually that of a Riemannian manifold. $\mathord{{\text{\rm Pos}}}_{n}$, as an open subset of the set $\mathord{{\text{\rm Sym}}}_{n}$ of symmetric matrices, in each point has a tangent space $T_{a}\mathord{{\text{\rm Pos}}}_{n}$ canonically isomorphic to $\mathord{{\text{\rm Sym}}}_{n}$ as a linear space. The Riemannian metric in this space is defined by $\|v\|=\|a^{-1/2}va^{-1/2}\|_{2}$.

$\mathord{{\text{\rm Pos}}}_{n}$ is also a symmetric space: We know that any $g\in\mathord{{\text{\rm GL}}}_{n}$ can be uniquely written as a product $g=au$, with $a=\sqrt{gg^{T}}\in\mathord{{\text{\rm Pos}}}_{n}$ and $u\in\mathord{{\text{\rm O}}}_{n}$. Thus $\mathord{{\text{\rm Pos}}}_{n}=G/H$, with $G=\mathord{{\text{\rm GL}}}_{n}$, $H=\mathord{{\text{\rm O}}}_{n}$, and the canonical projection $\pi(x)=\sqrt{xx^{T}}$.

The respective tangent spaces ${\mathord{\mathfrak{g}}},{\mathord{\mathfrak{h}}}$ of $G,H$ are given by ${\mathord{\mathfrak{g}}}={\mathbb{R}}^{n\times n}$ and ${\mathord{\mathfrak{h}}}={\mathfrak{so}}_{n}$, which is the set of skew-symmetric $n\times n$ matrices. The involution $\sigma(x)=-x^{T}$ in ${\mathord{\mathfrak{g}}}$ obeys $[\sigma(v),\sigma(w)]=\sigma([v,w])$, and ${\mathord{\mathfrak{h}}}$ is its $+1$ eigenspace. We have thus recognized $\mathord{{\text{\rm Pos}}}_{n}$ as a symmetric space. It turns out that $a\oplus v=a\exp(a^{-1}v)$, where exp is the matrix exponential function.

The previous paragraphs define two different structures on $\mathord{{\text{\rm Pos}}}_{n}$, namely that of a Riemannian manifold, and that of a symmetric space. They are compatible in the sense that the $\oplus$, $\ominus$ operations derived from either structure coincide. For more information we refer to [40, 24, 58]. Subdivision in particular is treated by [38]. ∎

Averages defined in terms of the exponential mapping

If $\oplus$ and $\ominus$ are defined as discussed in the previous paragraphs, it is possible to define a weighted affine average implicitly by requiring that

[TABLE]

Any Fréchet mean in a Riemannian manifold is also an average in this sense, which follows directly from (5) together with (6). Locally, $\operatorname{avg}_{E}$ is well defined and unique. As to the size of neighbourhoods where this happens, in the Riemannian case the proof given by [15, 16] for certain neighbourhoods enjoying unique existence of $\operatorname{avg}_{F}$ shows that the very same neighbourhoods also enjoy unique existence of $\operatorname{avg}_{E}$.

Affine averages with respect to a base point. From the different expressions originally given for the affine average, $x=y+\sum a_{j}(x_{j}-y)$ is one we have not yet defined a manifold analogue for. With $\ominus$ and $\oplus$ at our disposal, this can be done by

[TABLE]

We call this the log/exp average with respect to the base point $y$. It has the disadvantage of a dependence on the base point, but for the applications we have in mind, there frequently is a natural choice of base point. Its advantages lie in the easier analysis compared to the Fréchet mean. One should also appreciate that the Fréchet mean is a log/exp mean w.r.t. to a basepoint, if that basepoint is the Fréchet mean itself:

[TABLE]

because of (8). This may be a trivial point, but it has been essential in proving smoothness of limit curves for manifold-based subdivision processes (see Th. 3.2 and [29]).

The possibility to define averages w.r.t. basepoints rests on the possibility of defining $\ominus$, which has been discussed above.

2. Subdivision

2.1. Defining stationary subdivision

Subdivision is a refinement process acting on input data lying in some set ${\mathord{\mathcal{M}}}$, which in the simplest case are indexed over the integers and are interpreted as samples of a function $f\colon{\mathbb{R}}\to{\mathord{\mathcal{M}}}$. A subdivision rule refines the input data, producing a sequence $\mathord{{\text{\it Sp}}}$ which is thought of denser samples of either $f$ itself, or of a function approximating $f$.

One mostly considers binary rules, whose application “doubles” the number of data points. The dilation factor of the rule, generally denoted by the letter $N$, then equals $2$. We require that the subdivision rule is invariant w.r.t. index shift, which by means of the left shift operator $L$ can be formalized as

[TABLE]

We require that each point $\mathord{{\text{\it Sp}}}_{i}$ depends only on finitely many data points $p_{j}$. Together with shift invariance this means that there is $s>0$ such that $p_{i}$ influences only $\mathord{{\text{\it Sp}}}_{Ni-s},\ldots,\mathord{{\text{\it Sp}}}_{Ni+s}$.

Subdivision rules are to be iterated: We create finer and finer data

[TABLE]

which we hope approach a continuous limit (the proper definition of which is given below).

Subdivision was invented by G. de Rham [7], who considered the process of iteratively cutting corners from a convex polygon contained in ${\mathord{\mathcal{M}}}={\mathbb{R}}^{2}$, and asked for the limit shape. If cutting corners is done by replacing each edge $p_{i}p_{i+1}$ by the shorter edge with vertices $\mathord{{\text{\it Sp}}}_{2i}=(1-t)p_{i}+tp_{i+1}$, $\mathord{{\text{\it Sp}}}_{2i+1}=tp_{i}+(1-t)p_{i+1}$, this amounts to a subdivision rule. In de Rham’s example, only two data points $p_{i}$ contribute to any individual ${\mathord{{\text{\it Sp}}}_{\hskip-0.3014ptj}}$.

Primal and dual subdivision rules. The corner-cutting rules mentioned above are invariant w.r.t. reordering indices according to $\ldots,0\mapsto 1,1\mapsto 0,2\mapsto-1,\ldots$. With inversion $U$ defined by $(Up)_{i}=p_{-i}$ we can write this invariance as $(LU)S=S(LU)$. An even simpler kind of symmetry is enjoyed by subdivision rules with obey $US=SU$. The latter are called primal rules, the former dual ones. The reason why we emphasize these properties is that they give guidance for finding manifold analogues of linear subdivision rules.

Subdivision of multivariate data. It is not difficult to generalize the concept of subdivision to multivariate data $p\colon{\mathbb{Z}}^{s}\to{\mathord{\mathcal{M}}}$ indexed over the grid ${\mathbb{Z}}^{s}$. A subdivision rule $S$ must fulfill $L_{v}^{N}S=SL_{v}$, for all shifts $L_{v}$ w.r.t. a vector $v\in{\mathbb{Z}}^{s}$.

Data with combinatorial singularities have to be treated separately, cf. Sec. 3.4. Here basically only the bivariate case is studied, but this has been done extensively, mostly because of applications in Computer Graphics [43].

Linear subdivision rules and their nonlinear analogues

A linear subdivision rule acting on data $p\colon{\mathbb{Z}}^{2}\to{\mathbb{R}}^{d}$ has the form

[TABLE]

If the sum $\sum_{j}a_{i-Nj}$ of coefficients contributing to $\mathord{{\text{\it Sp}}}_{i}$ equals $1$, the application of the rule amounts to computing a weighted average:

[TABLE]

Subdivision rules not expressible in this way might occur as auxiliary tools in proofs, but are not meant to be applied to data which are points of an affine space. This is because if $\sum a_{i-Nj}\neq 1$, then the linear combination $\sum a_{i-Nj}p_{j}$ is not translation-invariant, and the rule depends on the choice of origin of the coordinate system.

Besides, the iterated application of rules not expressible as weighted averages either leads to divergent data $S^{k\!}p$, or alternatively, to data approaching zero. For this reason, one exclusively considers linear rules of the form (11). A common definition of convergent subdivision rule discounts the case of zero limits and recognizes translation invariance as a necessary condition for convergence, cf. [17].

For a thorough treatment of linear subdivision rules, conveniently done via $S$ acting as a linear operator in $\ell^{\infty}({\mathbb{Z}}^{s},{\mathbb{R}})$ and using the appropriate tools of approximation theory, see, e.g. [4].

In the following we discuss some nonlinear, geometric, versions of subdivision rules. We use the various nonlinear versions of averages introduced above, starting with the Fréchet mean in metric spaces.

$\bullet$*   Subdivision using the Féchet mean*. A natural analogue of (11) is found by replacing the affine average by the Fréchet mean. This procedure is particularly suited for Hadamard metric spaces and also in complete Riemannian manifolds.

$\bullet$*   Log/exp subdivision*. In a manifold equipped with an exponential mapping, an analogue of (11) is defined by

[TABLE]

where $m_{i}$ is a base point computed in a meaningful manner from the input data, e.g. $m_{i}=p_{\lfloor i/N\rfloor}$. In case of combinatorial symmetries of the subdivision rule, it makes sense to make the choice of $m_{i}$ conform to these symmetries.

$\bullet$*   Subdivision using projections*. If ${\mathord{\mathcal{M}}}$ is a surface embedded in a vector space and $\pi$ is a projection onto ${\mathord{\mathcal{M}}}$, we might use the subdivision rule

[TABLE]

If the intrinsic symmetries of ${\mathord{\mathcal{M}}}$ extend to symmetries of ambient space, then this projection analogue of a linear subdivision rule is even intrinsic – see Example 2.1.

• Example 2.

Subdivision in the motion group. The groups $\mathord{{\text{\rm O}}}_{n}$ and $\mathord{{\text{\rm SO}}}_{n}$ are ${1\over 2}n(n-1)$-dimensional surfaces in the linear space ${\mathbb{R}}^{n\times n}$. A projection onto $\mathord{{\text{\rm O}}}_{n}$ is furnished by singular value decomposition, or in an alternate way of expressing it, by the polar decomposition of Example 1:

[TABLE]

This projection is $\mathord{{\text{\rm O}}}_{n}$-equivariant in the sense that for $u\in\mathord{{\text{\rm O}}}_{n}$, we have both $\pi(ug)=u\pi(g)$ and $\pi(gu)=\pi(g)u$. The same invariance applies to application of a linear subdivision rule acting in ${\mathbb{R}}^{n\times n}$. So for any given data in $\mathord{{\text{\rm O}}}_{n}$, and a linear subdivision rule $S$, the subdivision rule $\pi\circ S$ produces data in $\mathord{{\text{\rm O}}}_{n}$ in a geometrically meaningful way, as long as we do not exceed the bounds of $\mathord{{\text{\rm GL}}}_{n}$. Since $\mathord{{\text{\rm GL}}}_{n}$ is a rather big neighbourhood of $\mathord{{\text{\rm O}}}_{n}$, this is in practice no restriction. Figure 1 shows an example. ∎

2.2. Convergence of subdivision processes

Definition of convergence. When discrete data $p$ are interpreted as samples of a function, then refined data $\mathord{{\text{\it Sp}}}$, $S^{2}p$ etc. are interpreted as the result of sampling which is $N$ times, $N^{2}$ times etc. as dense as the original. We therefore define a convergent refinement rule as follows.

• Definition 3

Discrete data $S^{k\!}p\colon{\mathbb{Z}}^{s}\to{\mathord{\mathcal{M}}}$ at the $k$-th iteration of refinement determine a function $f_{k}\colon N^{-k}{\mathbb{Z}}^{s}\to{\mathord{\mathcal{M}}}$, whose values are the given data points: For any $N$-adic point $\xi$, we have $(S^{k\!}p)_{N^{k}\xi}=f_{k}(\xi)$, provided $N^{k}\xi$ is an integer. For all such $\xi$, the sequence $(f_{k}(\xi))_{k\geq 0}$ is eventually defined and we let $f(\xi)=\lim_{k\to\infty}f_{k}(\xi)$. We say $S$ is convergent for input data $p$, if the limit function $f$ exists for all $\xi$ and is continuous. It can be uniquely extended to a continuous function $S^{\infty}p\colon{\mathbb{R}}^{s}\to{\mathord{\mathcal{M}}}$.

Another way of defining the limit is possible if data $p_{i},\mathord{{\text{\it Sp}}}_{i},\ldots$ lie in a vector space. We linearly interpolate them by functions $g_{0},g_{1},\ldots$ with $g_{k}(N^{-k}i)=S^{k}p_{i}$. Then the limit of functions $g_{k}$ agrees with the limit of Def. 2.2 (which is pointwise, but in fact convergence is usually uniform on compact sets.)

The following lemma is the basis for investigating convergence of subdivision rules in metric spaces. The terminology is that of [21, 20].

• Lemma 4

Let ${\mathord{\mathcal{M}}}$ be a complete metric space, and let the subdivision rule $S$ operate with dilation $N$ on data $p\colon{\mathbb{Z}}^{s}\to{\mathord{\mathcal{M}}}$. We measure the density of the data by

[TABLE]

where we use the 1-norm on the indices. $S$ is contractive, resp. displacement-safe, if

[TABLE]

If these two conditions are met, any input data with bounded density have a limit $S^{\infty}p$, which is Hölder continuous with exponent $-{\log\gamma\over\log N}$.

Proof.

Contractivity implies $\delta(S^{k\!}p)\leq\gamma^{k}\delta(p)$. For any $N$-adic rational point $\xi\in N^{-r}{\mathbb{Z}}^{s}$, the sequence $f_{k}(\xi)=(S^{k\!}p)_{N^{k}\xi}$ is defined for all $k\geq r$. It is Cauchy, since

[TABLE]

Thus the limit function $S^{\infty}p\equiv f$ is defined for all $N$-adic points.

Consider now two $N$-adic points $\xi,\eta$. Choose $k$ such that $N^{-(k+1)}\leq|\xi-\eta|\leq N^{-k}$. For all $j\geq k$, approximate $\xi$ resp. $\eta$ by $N$-adic points $a_{j},b_{j}\in N^{-j}{\mathbb{Z}}^{s}$, such that none of $|a_{j}-a|$, $|b_{j}-b|$, $|a_{j}-a_{j+1}|$ $|b_{j}-b_{j+1}|$ exceeds $sN^{-j}$. One can choose $a_{k}=b_{k}$. The sequence $a_{j}$ is eventually constant with limit $\xi$, and similarly the sequence $b_{j}$ is eventually constant with limit $\eta$. Using the symbol $(*)$ for “similar terms involving $b_{j}$ instead of $a_{j}$”, we estimate

[TABLE]

Using the contractivity and displacement-safe condition, we further get

[TABLE]

The index $k$ was chosen such that $k\leq-\log|\xi-\eta|/\log N$, so in particular $\gamma^{k}\leq\gamma^{-\log|\xi-\eta|/\log N}$. We conclude that

[TABLE]

Thus $f$ is continuous with Hölder exponent $-{\log\gamma\over\log N}$ on the $N$-adic rationals, and so is the extension of $f$ to all of ${\mathbb{R}}^{s}$. ∎

The scope of this lemma can be much expanded by some obvious modifications.

$\bullet$*   Input data with unbounded density $d(p)$*. Since points ${\mathord{{\text{\it Sp}}}_{\hskip-0.3014ptj}}$ only depend on finitely many $p_{i}$’s, there is $m>0$ such that $p_{i}$ only influences $\mathord{{\text{\it Sp}}}_{Ni+j}$ with $|j|<m$. By iteration, $p_{i}$ influences $S^{2}p_{N^{2}i+j}$ with $|j|<Nm+m$, and so on. It follows that $p_{i}$ influences the value $S^{\infty}p(i+\xi)$ of the limit function only for $|\xi|<{m\over N}+{m\over N^{2}}+\cdots={m\over N-1}$. We can therefore easily analyze the restriction of the limit function to some box by re-defining all data points away from that box in a manner which makes $d(p)$ finite.

$\bullet$*   Partially defined input data*. If data are defined not in all of ${\mathbb{Z}}^{s}$ but only in a subset, the limit function is defined for a certain subset of ${\mathbb{R}}^{s}$. Finding this subset goes along the same lines as the previous paragraph – we omit the details.

$\bullet$*   Convergence for special input data*. In order to check convergence for particular input data $p$, it is sufficient that the contractivity and displacement-safe conditions of Lemma 2.2 hold for all data $S^{k\!}p$ constructed by iterative refinement from $p$. A typical instance of this case is that contractivity can be shown only if $\delta(p)$ does not exceed a certain threshold $\delta_{0}$. It follows that neither does $\delta(S^{k\!}p)$, and Lemma 2.2 applies to all $p$ with $\delta(p)\leq\delta_{0}$.

$\bullet$*   Powers of subdivision rules*. A subdivision rule $S$ might enjoy convergence like a contractive rule without being contractive itself. This phenomenon is analogous to a linear operator $A$ having norm $\|A\|\geq 1$ but spectral radius $\rho(A)<1$, in which case some $\|A^{m}\|<1$. In that case we consider some power $S^{m}$ as a new subdivision rule with dilation factor $N^{m}$. If $S^{m}$ is contractive with factor $\gamma^{m}<1$, Lemma 2.2 still applies, and limits enjoy Hölder smoothness with exponent $-{\log\gamma^{m}\over\log N^{m}}=-{\log\gamma\over\log N}$.

• Example 5.

Convergence of linear subdivision rules. Consider a univariate subdivision rule $S$ defined by finitely many nonzero coefficients $a_{j}$ via (11). $S$ acts as a linear operator on sequences $p\colon{\mathbb{Z}}\to{\mathbb{R}}^{d}$. The norm $\|p\|=\sup_{i}\|p_{i}\|_{{\mathbb{R}}^{d}}$ induces an operator norm $\|S\|$ which obeys $\|\mathord{{\text{\it Sp}}}\|\leq\|S\|\|p\|$. It is an exercise to check $\|S\|=\max_{i}\sum_{j}|a_{i-Nj}|$. Equality is attained for suitable input data with values in $\{-1,0,1\}$.

With $(\Delta p)_{i}=p_{i+1}-p_{i}$ we express the density of the data as $\delta(p)=\sup\|\Delta p_{i}\|$. Contractivity means that $\sup\|\Delta\mathord{{\text{\it Sp}}}_{i}\|\leq\gamma\sup\|\Delta p_{i}\|$ for some $\gamma<1$.

Analysis of this contractivity condition uses a trick based on the generating functions $p(z)=\sum p_{j}z^{j}$ and $a(z)=\sum a_{j}z^{j}$. Equation (11) translates to the relation $(\mathord{{\text{\it Sp}}})(z)=a(z)p(z^{N})$ between generating functions, and we also have $\Delta p(z)=(z^{-1}-1)p(z)$. The trick consists in introducing the derived subdivision rule $S^{*}$ with coefficients $a_{j}^{*}$ which obeys $S^{*}\Delta=N\Delta S$. The corresponding relation between generating functions reads

[TABLE]

This division is possible in the ring of Laurent polynomials, because for all $i$, $\sum_{j}a_{i-Nj}=1$. The contractivity condition now reads $\sup\|\Delta\mathord{{\text{\it Sp}}}_{i}\|={1\over N}\sup\|S^{*}\Delta p_{i}\|\leq{1\over N}\|S^{*}\|\sup\|\Delta p_{i}\|$, i.e., the contractivity factor of the subdivision rule $S$ is bounded from above by ${1\over N}\|S^{*}\|$. The “displacement-safe” condition of Lemma 2.2 is fulfilled also, which we leave as an exercise (averages of points $p_{i}$ are not far from the $p_{i}$’s).

The above computation leads to a systematic procedure for checking convergence: we compute potential contractivity factors ${1\over N}\|S^{*}\|$, ${1\over N^{2}}\|S^{2*}\|$, and so on, until one of them is $<1$. The multivariate case is analogous but more complicated [19, 17, 4]. ∎

• Example 6.

Convergence of geodesic corner-cutting rules. Two points $a,b$ of a complete Riemannian manifold ${\mathord{\mathcal{M}}}$ are joined by a shortest geodesic path $t\mapsto a\oplus tv$, $v=b\ominus a$, $t\in[0,1]$. The difference vector $v$ and thus the path are generically unique, but do not have to be, if the distance between $a$ and $b$ exceeds both injectivity radii $\rho_{\text{inj}}(a),\rho_{\text{inj}}(b)$. The point $x=a\oplus tv$ has $\operatorname{\text{\it d}_{\!{\mathord{\mathcal{M}}}}}(a,x)=t\operatorname{\text{\it d}_{\!{\mathord{\mathcal{M}}}}}(a,b)$, $\operatorname{\text{\it d}_{\!{\mathord{\mathcal{M}}}}}(b,x)=(1-t)\operatorname{\text{\it d}_{\!{\mathord{\mathcal{M}}}}}(a,b)$. It is a Fréchet mean of points $a,b$ w.r.t. weights $(1-t),t$.

With these preparations, we consider two elementary operations on sequences, namely averaging $A_{t}$ and corner cutting $S_{t,s}$:

[TABLE]

The distance of $A_{t}p_{i}$ from $A_{t}p_{i+1}$ is bounded by the length of the broken geodesic path which connects the first point with $p_{i+1}$ and continues on to the second; its length is bounded by $\delta(p)$. Similarly, the distance of successive points of the sequence $S_{ts}p$, for $0\leq t<s\leq 1$ is estimated by $\max(1-(s-t),s-t)\delta(p)$. It follows immediately that a concatenation of operations of this kind is a subdivision rule where Lemma 2.2 applies, if at least one $S_{t,s}$ with $0<s-t<1$ is involved. Any such concatenation therefore is a convergent subdivision rule in any complete Riemannian manifold. A classical example are the rules $S^{(k)}=(A_{1/2})^{k}\circ S_{0,1/2}$, which insert midpoints,

[TABLE]

and then compute $k$ rounds of averages. E.g.,

[TABLE]

The rule $S^{(2)}$ (Chaikin’s rule, see [5]) is one of de Rham’s corner cutting rules. In the linear case, $S^{(k)}$ has coefficients $a_{j}={1\over 2^{k}}{k\choose j}$, apart from an index shift. Its limit curves are the B-spline curves whose control points are the initial data $p_{j}$ [45].

The corner-cutting rules discussed above are well defined and convergent in any Hadamard metric space – those spaces have geodesics in much the same way as Riemannian manifolds. Subdivision rules based on geodesic averaging (not necessarily restricted to values $t,s\in[0,1]$) have been treated by [51, 54, 21, 20]. We should also mention that adding a round $A_{1/2}$ to a subdivision increases smoothness of limit curves, which was recently confirmed in the manifold case [14]. ∎

• Example 7.

Convergence of interpolatory rules. A subdivision rule $S$ with dilation factor $N$ is called interpolatory if $\mathord{{\text{\it Sp}}}_{Ni}=p_{i}$, i.e., the old data points are kept and new data points are inserted in between. In the linear case, a very well studied subdivision rule of this kind is the four-point rule proposed by Dyn, Gregory and Levin [18]. We let $\mathord{{\text{\it Sp}}}_{2i}=p_{i}$ and

[TABLE]

In the special case $\omega={1\over 16}$, the point $\mathord{{\text{\it Sp}}}_{2i+1}$ is found by evaluating the cubic Lagrange polynomial interpolating $p_{i-1},\ldots,p_{i+2}$, which accounts for the high approximation order of $S$. There is in fact a whole series of interpolatory rules based on the idea of evaluating Lagrange interpolation polynomials (the Dubuc-Deslauriers subdivision schemes, see [8]).

$S$ is a binary “dual” subdivision rule with combinatorial symmetry about edges. Thus it makes sense to define a Riemannian version of $S$ by means of averages w.r.t. geodesic midpoints of $p_{i},p_{i+1}$ as base points, cf. Equ. (9). Using $m_{p_{i},p_{i+1}}=p_{i}\oplus{1\over 2}(p_{i+1}\ominus p_{i})$, we let

[TABLE]

The distance of successive points $\mathord{{\text{\it Tp}}}_{2i}$ and $\mathord{{\text{\it Tp}}}_{2i+1}$ is bounded by half the geodesic distance of $p_{i},p_{i+1}$ plus the length of the vector added to the midpoint in the previous formula. This yields the inequality $\delta(Tp)\leq{1\over 2}\delta(p)+2|\omega|{3\over 2}\delta(p)=({1\over 2}+3|\omega|)\delta(p)$. Lemma 2.2 thus shows convergence, if $|\omega|<1/6$.

We cannot easily extend this “manifold” four-point rule to more general metric spaces. The reason is that we used the linear structure of the tangent space. A general discussion of univariate interpolatory rules is found in [50]. ∎

2.3. Probabilistic interpretation of subdivision

in metric spaces

O. Ebner in [22, 23] gave a probabilistic interpretation of subdivision. This goes as follows. Consider a linear subdivision rule as in (11), namely

[TABLE]

acting on data $p\colon{\mathbb{Z}}^{s}\to{\mathbb{R}}^{d}$. Consider a stochastic process $J_{0},J_{1},\ldots$ defined as the random walk on ${\mathbb{Z}}^{s}$ with transition probabilities

[TABLE]

Then the expected value of $p_{J_{n+1}}$, conditioned on $J_{n}=j$ is given by

[TABLE]

by definition of the expected value. Now the expectation $\mathop{{\mathbb{E}\,}}(X)$ of an ${\mathbb{R}}^{d}$-valued random variable $X$ has a characterization via distances: $\mathop{{\mathbb{E}\,}}(X)$ is that constant $c\in{\mathbb{R}}^{d}$ which is closest to $X$ in the sense of $\mathop{{\mathbb{E}\,}}(d(X,c)^{2})\to\min$. A similar characterization works for the conditional expectation $\mathop{{\mathbb{E}\,}}(X|Y)$ which is the random variable $f(Y)$ closest to $X$ in the $L^{2}$ sense. These facts inspired a theory of random variables with values in Hadamard metric spaces developed by K.-T. Sturm [46, 47]. The minimizers mentioned above can be shown to still exist if ${\mathbb{R}}^{d}$ is replaced by ${\mathord{\mathcal{M}}}$.

Since the way we compute subdivision by Fréchet means is compatible with the distance-based formula for expected values, Equation (13) holds true also in the case that both the expectation and the subdivision rule are interpreted in the Hadamard space sense. On that basis, O. Ebner could show a remarkable statement on convergence of subdivision rules:

• Theorem 8

[23, Th. 1]* Consider a binary subdivision rule $\mathord{{\text{\it Tp}}}_{i}=\operatorname{avg}_{F}(a_{i-2j};p_{j})$ with nonnegative coefficients $a_{i}$. It produces continuous limits for any data $p_{j}$ in any Hadamard space ${\mathord{\mathcal{M}}}$ if and only if it produces a continuous limit function when acting on real-valued data.*

Sketch of proof.

With the random walk $(J_{i})_{i=0,1,\ldots}$ defined above, (13) directly implies

[TABLE]

Unlike for ${\mathbb{R}}^{d}$-valued random variables, there is no tower property for iterated conditioning, so in general $(T^{n}p)_{J_{0}}\neq\mathop{{\mathbb{E}\,}}(p_{J_{n}}|J_{0})$. That expression has a different interpretation: $T$ is analogous to the linear rule $S$ of (12), which is nothing but the restriction of the general rule $T$ to data in Euclidean spaces. Its $n$-th power $S^{n}$ is a linear rule of the form $(S^{n}q)_{i}=\sum a^{[n]}_{i-2^{n}j}q_{j}$, and we have

[TABLE]

This follows either directly (computing the coefficients of the $n$-th iterate $S^{n}$ corresponds to computing transition probabilites for the $n$-iterate of the random walk), or by an appeal to the tower property in (14).

Sturm [46] showed a Jensen’s inequality for continuous convex functions $\Psi$,

[TABLE]

We choose $\Psi=\operatorname{\text{\it d}_{\!{\mathord{\mathcal{M}}}}}(\cdot,x)$ and observe that $q_{J_{n}}=\operatorname{\text{\it d}_{\!{\mathord{\mathcal{M}}}}}(p_{J_{n}},x)$ is a real-valued random variable. Combining Jensen’s inequality with (14) and (15) yields

[TABLE]

To continue, we need some information on the coefficients $a^{[n]}_{i}$. For that, we use the limit function $\phi\colon{\mathbb{R}}^{s}\to[0,1]$ generated by applying $T$ (or rather, $S$), to the delta sequence. By construction (see Lemma 2.2), $|a^{[n]}_{j}-\phi(2^{-n}j)|\to 0$ as $n\to\infty$. These ingredients allow us to show existence of $n$ with $T^{n}$ contractive. ∎

As a corollary we get, for instance, that subdivision with nonnegative coefficients works in $\mathord{{\text{\rm Pos}}}_{n}$ in the same way as in linear spaces, as far as convergence is concerned. Since $\mathord{{\text{\rm Pos}}}_{n}$ is not only a Hadamard metric space, but even a smooth Riemannian manifold, also the next section will yield a corollary regarding $\mathord{{\text{\rm Pos}}}_{n}$.

2.4. The convergence problem in manifolds

The problem of convergence of subdivision rules in manifolds (Riemannian manifolds, groups, and symmetric spaces) was at first treated by means of so-called proximity inequalities which compare linear rules with their analogous counterparts in manifolds. This approach was successful in studying smoothness of limits (see Section 3 below), but less so for convergence. Unless subdivision rules are of a special kind (interpolatory, corner-cutting$\,,\ldots$) convergence can typically be shown only for “dense enough” input data, with very small bounds on the maximum allowed density. On the other hand numerical experiments demonstrate that a manifold rule analogous to a convergent linear rule usually converges. This discrepancy between theory and practice is of course unsatisfactory from the viewpoint of theory, but is not so problematic from the viewpoint of practice. The reason is the stationary nature of subdivision — if $\delta(p)$ is too big to infer existence of a continuous limit $S^{\infty}p$, we can check if $\delta(S^{k\!}p)$ is small enough instead. As long as $S$ converges, this leads to an a-posteriori proof of convergence.

More recently, convergence of subdivision rules of the form $\mathord{{\text{\it Sp}}}_{i}=\operatorname{avg}_{F}(a_{i-Nj};p_{j})$ in Riemannian manifolds has been investigated along the lines of Lemma 2.2. This work is mainly based on the methods of H. Karcher’s seminal paper [39]. So far, only the univariate case of data $p\colon{\mathbb{Z}}\to{\mathord{\mathcal{M}}}$ has been treated successfully, cf. [53, 37, 36].

There are two main cases to consider. In Cartan-Hadamard manifolds (curvature $\leq 0$) the Fréchet mean is well defined and unique also if weights are allowed to be negative [37, Th. 6]. Subdivision rules are therefore globally and uniquely defined. We have the following result:

• Prop. 9

[37, Th. 11]* Consider a univariate subdivision rule $\mathord{{\text{\it Sp}}}_{i}=\operatorname{avg}_{F}(a_{i-Nj};p_{j})$ acting on sequences in a Cartan-Hadamard manifold ${\mathord{\mathcal{M}}}$. Consider also the norm $\|S^{*}\|$ of its linear derived subdivision rule according to Example 2.2. If*

[TABLE]

then $S$ meets the conditions of Lemma 2.2 (with contractivity factor $\gamma$) and produces continuous limits.

This result is satisfying because it allows us to infer convergence from a condition which is well known in the linear case, cf. [17]. If ${1\over N}\|S^{*}\|\geq 1$, we can instead check if one of ${1\over N^{n}}\|S^{*n}\|$, $n=2,3,\ldots$ is smaller than $1$. If this is the case, then the manifold subdivision rule analogous to the linear rule $S^{n}$ converges.

Subdivision in Riemannian manifolds with positive curvature. Recent work [36] deals with spaces of positive curvature, and initial results have been achieved on the unit sphere, for subdivision rules of the form $\mathord{{\text{\it Sp}}}_{i}=\operatorname{avg}_{F}(a_{i-2j};p_{j})$. Figure 3 shows two examples. One aims at finding a bound $\delta_{0}$ such that for all data $p$ with $\delta(p)<\delta_{0}$, $S$ acts in a contractive way so that Lemma 2.2 shows convergence.

Rules defined in a different way are sometimes much easier to analyze. E.g. the Lane-Riesenfeld subdivision rules defined by midpoint insertion, followed by $k$ rounds of averaging, can be transferred to any complete Riemannian manifold as a corner-cutting rule and will enjoy continuous limits, see Example 2.2. Similarly, the interpolatory four-point rule can be generalized to the manifold case in the manner described by Example 2.2, and will enjoy continuous limits. The generalization via the Fréchet mean (Fig. 3) on the other hand, is not so easy to analyze. The approach by [36] is to control $\delta(\mathord{{\text{\it Sp}}})$ by introducing a family $S^{(t)}$, $0\leq t\leq 1$, of rules where $S^{(0)}$ is easy to analyze, and $S^{(1)}=S$. If one manages to show $\delta(S^{(0)}p)<\gamma_{1}\delta(p)$ and $\|{d\over dt}S^{(t)}p_{i}\|\leq C\delta(p)$, then the length of each curve $t\mapsto S^{(t)}p_{i}$ is bounded by $C\delta(p)$, and

[TABLE]

Contractivity is established if $\gamma_{1}+2C<1$, in which case Lemma 2.2 shows convergence. The bounds mentioned in Fig. 3 have been found in this way. Estimating the norm of the derivative mentioned above involves estimating the eigenvalues of the Hessian of the right hand side of (5).

The state of the art regarding convergence of refinement schemes. Summing up, convergence of geometric subdivision rules is treated in a satisfactory manner for special rules (interpolatory, corner-cutting), for rules in special spaces (Hadamard spaces and Cartan-Hadamard manifolds), and in the very special case of the unit sphere and univariate rules. General manifolds with positive curvature have not been treated. Multivariate data are treated only in Hadamard metric spaces and for subdivision rules with nonnegative coefficients. In other situations, we know that convergence happens only for “dense enough” input data, where the required theoretical upper bounds on $\delta(p)$ are very small compared to those inferred from numerical evidence.

3. Smoothness analysis of subdivision rules

For linear subdivision rules, the question of smoothness of limits can be considered as largely solved, the derived rule $S^{*}$ introduced in Example 2.2 being the key to the question if limits are smooth. Manifold subdivision rules do not always enjoy the same smoothness as the linear rules they are derived from. The constructions mentioned in Section 2 basically yield manifold rules whose limits enjoy $C^{1}$ resp. $C^{2}$ smoothness if the original linear rule has this property, but this general statement is no longer true if $C^{3}$ or higher smoothness is involved. Manifold rules generated via Fréchet means or via projection [26, 59] retain the smoothness of their linear counterparts. Others, e.g. constructed by means of averages w.r.t. basepoints in general do not. This is to be expected, since the choice of basepoint introduces an element of arbitrariness into manifold subdivision rules. The following paragraphs discuss the method of proximity inequalities which was successfully employed in treating the smoothness of limits.

3.1. Derivatives of limits

A subdivision rule $S$ acting on a sequence $p$ in ${\mathbb{R}}^{d}$ converges to the limit function $S^{\infty}p$, if the refined data $S^{k\!}p$, interpreted as samples of functions $f_{k}$ at the finer grid $N^{-k}{\mathbb{Z}}$, approach that limit function (see Definition 2.2):

[TABLE]

whenever $N^{k}\xi$ is an integer. A similar statement holds for derivatives, which are approximated by finite differences. With $h=N^{-k}$, we get

[TABLE]

Here $S^{*}$ is the derived rule defined by the relation $S^{*}\Delta=N\Delta S$, see Ex. 2.2. For the $r$-th derivative of the limit function we get

[TABLE]

These relations, except for references to derived rules, are valid even if $S$ does not act linearly. $S$ could be a manifold rule expressed in a coordinate chart, or it could be acting on a surface contained in ${\mathbb{R}}^{d}$.

If $S$ does act linearly, one proves that $S$ has $C^{1}$ smooth limits, if $S^{*}$ has continuous ones, and in that case $(S^{\infty}p)^{\prime}=S^{*\infty}\Delta p$. To treat higher order derivatives, this statement can be iterated. For multivariate data $p_{i}$, $i\in{\mathbb{Z}}^{s}$, the situation is analogous but more complicated to write down. For the exact statements, see [17, 4].

3.2. Proximity inequalities

Smoothness from proximity. Manifold subdivision rules were first systematically analyzed with regard to derivatives by [51]. The setup is a linear rule $S$ and a nonlinear rule $T$ both acting on data contained in the same space ${\mathbb{R}}^{d}$. $T$ could be a manifold version of $S$, with ${\mathbb{R}}^{d}$ being a coordinate chart of the manifold; or $T$ could act on points of a surface contained in ${\mathbb{R}}^{d}$. Then $S$, $T$ are in proximity, if

[TABLE]

This formula is motiviated by a comparison of the shortest path between two points within in a surface (which is a geodesic segment), with the shortest path in Euclidean space (which is a straight line). These two paths differ by exactly the amount stated in (16). Two statements were shown in [51]:

• (1)

Certain manifold subdivision rules $T$ derived from a convergent linear rule $S$ obey the proximity inequality (16) whenever data are dense enough (i.e., $\delta(p)$ is small enough).

• (2)

in that case, if limit curves of $S$ enjoy $C^{1}$ smoothness, then $T$ produces continuous limit curves for data with $d(p)$ small enough; and all continuous limit curves enjoy $C^{1}$ smoothness.

To demonstrate how proximity inequalities work, we prove a convergence statement like the ones given by [54, Th. 1] or [51, Th. 2+3] (with slightly different proofs).

• Prop. 10

Assume the setting of Equ. (16), with a subdivision rule $T$ being in proximity with a linear subdivision rule $S$. We also assume ${1\over N}\|S^{*}\|<1$.333implying convergence of the linear rule $S$. $N$ is the dilation factor, $S^{*}$ is the derived rule, cf. Ex. 2.2. Then $T$ produces continuous limit curves from data $p$ with $\delta(p)$ small enough.

Proof.

Generally $\sup_{i}\|p_{i}-q_{i}\|\leq K\implies\delta(p)\leq\delta(q)+2K$. Thus (16) implies

[TABLE]

Choose $\epsilon>0$ with $\lambda:={1\over N}\|S^{*}\|+2C\epsilon<1$. If $\delta(p)<\epsilon$, then $T$ is contractive:

[TABLE]

By recursion, $\delta(T^{k+1}p)\leq\lambda\delta(T^{k}p)$. As to the displacement-safe condition of Lemma 2.2, recall from Example 2.2 that $S$ has it. For $T$, observe that

[TABLE]

Now Lemma 2.2 shows convergence. ∎

The convergence of vectors $N^{k}\Delta T^{k}p$ to derivatives of the limit function $T^{\infty}p$ is proved in a way which is analogous in principle. The method was extended to treat $C^{2}$ smoothness by [49], using the proximity condition

[TABLE]

A series of publications treated $C^{2}$ smoothness of Lie group subdivision rules based on log/exp averages [52, 31], the same in symmetric spaces [53], $C^{1}$ smoothness in the multivariate case [25], higher order smoothness of interpolatory rules in groups [27, 62], and higher order smootheness of projection-based rules [26, 59]. The proximity conditions involving higher order smoothness become rather complex, especially in the multivariate case.

Smoothness equivalence. If a manifold subdivision rule $T$ is created on basis of a linear rule $S$, it it interesting to know if the limit functions of $T$ enjoy the same smoothness as the limits of $S$. For $C^{1}$ and $C^{2}$ smoothness, $T$ can basically be constructed by any of the methods described above, and it will enjoy the same smothness properties as $S$ (always assuming that convergence happens, and that the manifold under consideration is itself as smooth as the intended smoothness of limits). This smoothness equivalence breaks down for $C^{k}$ with $k\geq 3$.

A manifold subdivision rule based on the log/exp construction, using averages w.r.t. basepoints,

[TABLE]

does not enjoy $C^{k}$ smoothness equivalence for $k\geq 3$ unless the base points $m_{i}$ obey a technical condition which can be satisfied e.g. if they themselves are produced by certain kinds of subdivision [60, 29]. Necessary and sufficient conditions for smoothness equivalence are discussed by [13]. We pick one result whose proof is based on this method (using (10) for a “base point” interpretation of Fréchet means):

• Theorem 11

[29, Th. 4.3]* Let $S$ be a stable444“Stable” means existence of constants $C_{1},C_{2}$ with $C_{1}\|p\|\leq\|S^{\infty}p\|_{\infty}\leq C_{2}\|p\|$ for all input data where $\|p\|:=\sup_{i}\|p_{i}\|$ is bounded. Stable rules with $C^{n}$ limits generate polynomials of degree $\leq n$, which is a property used in the proof. subdivision rule $\mathord{{\text{\it Sp}}}_{i}=\operatorname{avg}(a_{i-Nj}$; $p_{j})$ acting on data $p\colon{\mathbb{Z}}^{s}\to{\mathbb{R}}^{d}$, which is convergent with $C^{n}$ limits. Then all continuous limits of its Riemannian version $\mathord{{\text{\it Tp}}}_{i}$ $=$ $\operatorname{avg}_{F}(a_{i-Nj}$; $p_{j})$ likewise are $C^{n}$.*

We should also mention that proximity conditions relevant to the smoothness analysis of manifold subdivision rules can take various forms, cf. the “differential” proximity condition of [30, 13, 12].

Finally we point out a property which manifold rules share with linear ones: For any univariate linear rule $S$ which has $C^{k}$ limits, the rule $A_{1/2}^{k}\circ S$ has limits of smoothess $C^{n+k}$, where $A_{1/2}$ is midpoint-averaging as described by Ex. 2.2. It has been shown in [14] that an analogous statement holds true also in the manifold case, for a general class of averaging operators.

3.3. Subdivision of Hermite data

Hermite subdivision is a refinement process acting not on points, but on tangent vectors, converging to a limit and its derivative simultaneously. In the linear case, data $(p,v)\colon{\mathbb{Z}}\to{\mathbb{R}}^{d}\times{\mathbb{R}}^{d}$ undergo subdivision by a rule $S$ which obeys basic shift invariance $SL=L^{N}S$. The interpretation of $p_{i}$ as points and $v_{i}$ as vectors leads to

[TABLE]

$S$ is invariant w.r.t. translations, which act via $p\mapsto p+x$ on points, but act identically on vectors. Iterated refinement creates data $S^{k}\mathord{\textstyle{p\choose v}}$ converging to a limit $f\colon{\mathbb{R}}\to{\mathbb{R}}^{d}$,

[TABLE]

We say that $S$ converges, if the limit $(f,f^{\prime})$ exists and $f$ enjoys $C^{1}$ smoothness, with $f^{\prime}$ then being continuous. A manifold version of $S$, operating on data

[TABLE]

faces the difficulty that each $v_{i}$ is contained in a different vector space. One possibility to overcome this problem is to employ parallel transport $\mathop{{\text{\rm pt}\hskip 0.43057pt}}\nolimits_{p}^{q}\colon T_{p}{\mathord{\mathcal{M}}}\to T_{q}{\mathord{\mathcal{M}}}$ between tangent spaces. In Riemannian manifolds, a natural choice for $\mathop{{\text{\rm pt}\hskip 0.43057pt}}\nolimits_{p}^{q}$ is parallel transport w.r.t. the canonical Levi-Civita connection along the shortest geodesic connecting $p$ and $q$, cf. [9]. In groups, we can simply choose $\mathop{{\text{\rm pt}\hskip 0.43057pt}}\nolimits_{p}^{q}$ as left translation by $qp^{-1}$ resp. the differential of this left translation. Then the definition

[TABLE]

is meaningful (provided the base point $m_{i}$ is chosen close to $p_{\lfloor i/N\rfloor}$). In a linear space, this expression reduces to (19). C. Moosmüller could show $C^{1}$ smoothness of limits of such subdivision rules, by methods in the spirit of Section 3.2, see [41, 42].

3.4. Subdivision with irregular combinatorics

A major application of subdivision is in Computer Graphics, where it is ubiquitously used as a tool to create surfaces from a finite number of handle points whose arrangement is that of the vertices of a 2D discrete surface. That surface usually does not have the combinatorics of a regular grid.

Two well known subdivision rules acting on such data are the Catmull-Clark rule and the Doo-Sabin rule, see [3, 6]. Such subdivision rules create denser and denser discrete surfaces which are mostly regular grids but retain a constant number of combinatorial singularities. This implies that the limit surface is locally obtained via Definition 2.2, but with a nontrivial overlapping union of several such limits as one approaches a combinatorial singularity. A systematic way of analyzing convergence and smoothness was found by U. Reif [44], see also the monograph [43]. There is a wealth of contributions to this topic, mostly because of its relevance for Graphics.

A. Weinmann in [55, 57, 56] studied intrinsic manifold versions of such subdivision rules. They are not difficult to define, since the linear subdivision rules which serve as a model are defined in terms of averages. We do not attempt to describe the methods used for establishing convergence and smoothness of limits other than to say that a proximity condition which holds between a linear rule $S$ and a nonlinear rule $T$ eventually guarantees that in the limit, smoothness achieved by $S$ carries over to $T$ — the perturbation incurred by switching from a linear space to a manifold is not sufficient to destroy smoothness. Figure 5 illustrates a result obtained by [55].

4. Multiscale transforms

4.1. Definition of intrinsic multiscale transforms

A natural multiscale representation of data, which does not suffer from distortions caused by the choice of more or less arbitrary coordinate charts, is required to be based on operations which are themselves adapted to the geometry of the data. This topic is intimately connected to subdivision, since upscaling operations may be interpreted as subdivision.

A high-level introduction of certain kinds of multiscale decompositions is given by [33]. We start with an elementary example.

• Example 12.

A geometric Haar decomposition and reconstruction procedure. Consider data $p\colon{\mathbb{Z}}\to{\mathord{\mathcal{M}}}$, and the upscaling rule $S$ and downscaling rule $D$,

[TABLE]

The use of $\oplus$ and $\ominus$ refers to the exponential mapping, as a means of computing differences of points, and adding vectors to points. $D$ is a left inverse of $S$ but not vice versa: $SDp\neq p$ in general. However, if we store the difference between $p$ and $SDp$ as detail vectors $q$:

[TABLE]

then the reconstruction procedure

[TABLE]

recovers the information destroyed by downsampling. ∎

More systematically, we have employed two upscaling rules $S,R$ and two downscaling rules $D,Q$ which obey

[TABLE]

($L$ is left shift). We have data $p^{(j)}$ at level $j$, $j=0,\dots,M$, where we interpret the data at the highest (finest) level as given, and the data at lower (coarser) level computed by downscaling. We also store details $q^{(j)}$ at each level:

[TABLE]

We require that upscaled level $j-1$ data and level $j$ details can restore level $j$ data:

[TABLE]

Generally, $S$, $D$ compute points from points, so they are formulated via averages:

[TABLE]

In Example 4.1, averages are computed w.r.t. base points $p_{\lfloor i/2\rfloor}$ for $S$ resp. $p_{i}$ for $D$, and coefficients $a_{j}$ and $b_{j}$ vanish except $a_{0}=a_{1}=1$, $b_{0}=b_{1}={1\over 2}$.

The downscaling operator $Q$ acts on tangent vectors $v_{i}=p_{i}\ominus(SDp)_{i}\in T_{p_{i}}{\mathord{\mathcal{M}}}$, so it has to deal with vectors potentially contained in different vector spaces. In our special case, $Q$ simply forgets one half of the data:

[TABLE]

Finally, the upscaling operator $R$ takes the vectors stored in $q^{(j)}$ and converts them into vectors which can be added to upscaled points $Sp^{(j-1)}$. Thus $R$ potentially has to deal with vectors contained in different tangent spaces. In our special case, the points $(Sp^{(j-1)})_{2i}$ $(Sp^{(j-1)})_{2i+1}$ both coincide with $p^{(j-1)}_{i}$, and that is also the point where the detail coeffient $q^{(j)}_{i}$ is attached to. We therefore might be tempted to write $(Rq)_{2i}=q_{i}$, $(Rq)_{2i+1}=-q_{i}$. This simple rule however does not take into account that along reconstruction, data and details could have been modified, and no longer fit together. We therefore use parallel transport to move the vector to the right tangent space, just in case:

[TABLE]

The symbol $\mathop{{\text{\rm pt}\hskip 0.43057pt}}\nolimits^{(p^{(j-1)})_{i}}(q_{i})$ refers to transporting $q_{i}$ to a tangent vector attached to $({p^{(j-1)}})_{i}$, see Section 3.3.

The operations $S,R,D,Q$ must be compatible, in the sense that reconstruction is a left inverse of downscaling plus computing details. While in the linear case, where $S,D,R,Q$ are linear operators on $\ell^{\infty}({\mathbb{R}}^{d})$, one usually requires $QR={\mathord{\text{\rm id}}}$ and $QS=0$ as well as $SD+RQ={\mathord{\text{\rm id}}}$, in the geometric case we must be careful not to mix operations on points with operations on tangent vectors. We therefore require

[TABLE]

• Example 13.

Interpolatory wavelets. Consider an interpolatory subdivision rule $S$ with dilation factor $2$, i.e., $\mathord{{\text{\it Sp}}}_{2i}=p_{i}$, and the forgetful downscaling operator $p^{(j-1)}_{i}=(Dp^{(j)})_{i}=p^{(j)}_{2i}$. If we store as details the difference vectors between $SDp$ and $p$ for odd indices, the data points $p_{2i+1}$ can be easily reconstructed:

[TABLE]

This procedure fits into the general scheme described above if we we let $Q=DL$ ($L$ is left shift) and define the upscaling of details by $(Rq)_{2i}=0$, $(Rq)_{2i+1}=q_{i}$. To admit the possibility that before reconstruction, data and details have been changed, we define

[TABLE]

in order to account for the possibility that $q_{i}^{(j)}$ is not yet contained in the “correct” tangent space. The decimated data $p^{(j-1)}$ together with details $q^{(j)}$ ($j\leq M$) may be called a geometric interpolatory-wavelet decomposition of the data at the finest level $p^{(M)}$. That data itself comes e.g. from sampling a function, cf. [10]. ∎

Definability of multiscale transforms without redundancies. The previous examples use upscaling and downscaling operations which are rather simple, except that in Example 4.1 one may use any interpolatory subdivision rule. It is also possible to extend Example 4.1 to the more general case of a midpoint-interpolating subdivision rule $S$, which is a right inverse of the decimation operator $D$. In [33] it is argued that it is highly unlikely that in the setup described above, which avoids redundancies, more general upscaling and downscaling rules will manage to meet the compatibility condition (22) needed for perfect reconstruction. In the linear case, where all details $q_{i}^{(j)}$ are contained in the same vector space, (22) is merely an algebraic condition on the coefficents involved in the definition of $S,D,Q,R$ which can be solved. In the geometric case, the usage of parallel transport makes a fundamental difference.

4.2. Properties of multiscale transforms

Charcterizing smoothness by coefficient decay. One purpose of a multiscale decomposition of data is to read off properties of the original data. Classically, the faster the magnitude of detail coefficients $q^{(j)}_{i}$ decays as $j\to\infty$, the smoother the original data. A corresponding result for the interpolatory wavelets of Example 4.1 in the linear case is given by [10, Th. 2.7]. To state a result in the multivariate geometric case, let us first introduce new notation for interpolatory wavelets, superseding Example 4.1.

We consider an interpolatory subdivision rule $S$ acting with dilation factor $N$ on data $p\colon{\mathbb{Z}}^{s}\to{\mathord{\mathcal{M}}}$. We define data $p^{(j)}$ at level $j$ as samples of a function $f\colon{\mathbb{R}}^{s}\to{\mathord{\mathcal{M}}}$, and construct detail vectors similar to Example 4.1:

[TABLE]

This choice is consistent with the decimation operator $Dp_{i}:=p_{Ni}$. The difference to Example 4.1 is firstly that here we allow multivariate data, and secondly that we do not “forget” redundant information such as $q^{(j)}_{Ni}=0$.

The result below uses the notation ${\operatorname{Lip\,}}\gamma$ for functions which are $C^{k}$ with $k=\lfloor\gamma\rfloor$ and whose $k$-th derivatives are Hölder continuous of exponent $\gamma-k$. The critical Hölder regularity of a function $f$ is the supremum of $\gamma$ such that $f\in{\operatorname{Lip\,}}\gamma$.

• Theorem 14

[32, Th. 8]* Assume that the interpolatory upscaling rule $S$, when acting linearly on data $p\colon{\mathbb{Z}}^{s}\to{\mathbb{R}}$, reproduces polynomials of degree $\leq d$ and has limits of critial Hölder regularity $r$.*

Consider a continous function $f\colon{\mathbb{R}}^{s}\to{\mathord{\mathcal{M}}}$, and construct detail vectors $q^{(j)}$ at level $j$ for the function $x\mapsto f(\sigma\cdot s)$ for some $\sigma>0$ (whose local existence is guaranteed for some $\sigma>0$).

Then $f\in{\operatorname{Lip\,}}\alpha$, $\alpha<d$ implies that detail vectors decay with $\sup_{i}\|q^{(j)}_{i}\|\leq C\cdot N^{-\alpha j}$ as $j\to\infty$. Conversely, that decay rate together with $\alpha<r$ implies $f\in{\operatorname{Lip\,}}\alpha$. The constant is understood to be uniform in a compact set.

The manifold ${\mathord{\mathcal{M}}}$ can be any of the cases we defined $\oplus$ and $\ominus$ operations for. Of course, smoothness of $f\colon{\mathbb{R}}^{s}\to{\mathord{\mathcal{M}}}$ is only defined up to the intrinsic smoothness of ${\mathord{\mathcal{M}}}$ as a differentiable manifold. An example of an upscaling rule $S$ is the four-point scheme with parameter $1/16$ mentioned in Example 2.2, which reproduces cubic polynomials and has critical Hölder regularity 2, cf. [18].

The proof is conducted in a coordinate chart (it does not matter which), and uses a linear vision of the theorem as an auxiliary tool. It further deals with the extensive technicalities which surround proximity inequalities in the multivariate case.

It is worth noting that A. Weinmann in [56] succeeded in transferring these ideas to the combinatorially irregular setting. The results are essentially the same, with the difference that one can find upscaling rules only up to smoothness $2-\epsilon$.

Stability. Compression of data is a main application of multiscale decompositions, and it is achieved e.g. by thresholding or quantizing detail vectors. It is therefore important to know what effect these changes have when reconstruction is performed. What we bascially want to know is whether reconstruction is Lipschitz continuous. In the linear case the problem does not arise separately, since the answer is implicitly contained in norms of linear operators. For the geometric multiscale transforms defined by upscaling operations $S,R$ and downscaling operations $D,Q$ according to (22), this problem is discussed by [33]. Consider data $p^{(j)}$ at level $j$ with $p^{(j-1)}=Dp^{(j)}$ such that $\delta(p^{(j)})\leq C\mu^{j}$, for some $\mu<1$. Consider recursive reconstruction of data $p^{(j)}$ from $p^{(0)}$ and details $q^{(1)},\ldots,q^{(j)}$ according to Equation (21). Then there are constants $C_{k}$ such that for modified details $\tilde{q}^{(j)}$, leading to modified data $\tilde{p}^{(j)}$, we have the local Lipschitz-style estimate

[TABLE]

It refers to a coordinate chart of the manifold ${\mathord{\mathcal{M}}}$ (it does not matter which).

Approximation Order. For an interpolatory upscaling operator $S$, and data $p_{i}\in{\mathord{\mathcal{M}}}$ defined by sampling, $p_{i}=f(h\cdot i)$, we wish to know to what extent the original function differs from the limit created by upscaling the sample. We say that $S$ has approximation order $r$, if there are $C>0$, $h_{0}>$ such that for all $h<h_{0}$

[TABLE]

It was shown by [61] that a manifold subdivision rule has in general the same approximation order as the linear rule we get by restricting $S$ to linear data.

This question is directly related to stability as discussed above: Both $f$ and $S^{\infty}p$ can be reconstructed from samples $p^{(j)}$, if $h=N^{-j}$: Detail vectors $q^{(k)}$, $k>j$, according to (23) reconstruct $f$, whereas details $\tilde{q}^{(k)}=0$ reconstruct $S^{\infty}p$. Stability of reconstruction and knowledge of the asymptotic magnitude of details $q^{(k)}_{i}$, $k>j$ directly corresponds to approximation order. On basis of this relationship one can again show an approximation order equivalence result, cf. [28].

Conclusion

The preceding pages give an account of averages, subdivision, and multiscale transforms defined via geometric operations which are intrinsic for various geometries (metric spaces, Riemannian manifolds, Lie groups, and symmetric spaces). We reported on complete solutions in special cases (e.g. convergence of subdivision rules in Hadmard metric spaces) and on other results with much more general scope as regards the spaces and subdivision rules involved, but with more restrictions on the data they apply to.

Acknowledgments

The author gratefully acknowledges the support of the Austrian Science fund through grant No. W1230, and he also wants to thank all colleagues whose work he was able to present in this survey.