Nonstationary Subdivision Schemes

• Nonstationary subdivision schemes are level-dependent refinement processes that vary their masks across iterations to exactly reproduce exponential-polynomial spaces.
• They enable adaptive geometric design, multiresolution analysis, and wavelet construction by using tension parameters and variable mask support for enhanced shape control.
• Analytical tools like approximate sum rules, spectral radius criteria, and symbol calculus ensure convergence, smoothness, and stability in both uniform and non-uniform settings.

Searching arXiv for recent and foundational papers on nonstationary subdivision schemes to ground the article. Searching arXiv for regularity, reproduction, and geometric applications of nonstationary subdivision. Nonstationary subdivision schemes are level-dependent refinement processes in which the refinement rule changes with the iteration index. In the standard scalar setting, one starts from discrete data on $\mathbb Z^s$ and refines by masks $a^{[k]}$ or $a(r)$ that vary with the level, rather than by a single stationary mask. This simple change enlarges the class of exactly reproducible functions from polynomials to exponential polynomials, allows tension or shape parameters to vary with scale, supports non-uniform and adaptive constructions, and creates a direct link with geometric design, wavelets, multiresolution analysis, and surface refinement on irregular topologies (Charina et al., 2015, Charina et al., 2013).

1. Formal setting and basic concepts

A scalar multivariate non-stationary subdivision scheme with dilation matrix $M\in\mathbb Z^{s\times s}$ is defined by the recursion

$$f^{[k+1]}_\alpha=\sum_{\beta\in \mathbb Z^s} a^{[k]}_{\alpha-M\beta}\, f^{[k]}_\beta, \qquad \alpha\in\mathbb Z^s,\ k\ge 0,$$

or, equivalently,

$c^{(r+1)} = S_{a(r)}c^{(r)}.$

Each mask is finitely supported, and its symbol is the Laurent polynomial

$$a^{[k]}(z)=\sum_{\alpha\in\mathbb Z^s} a^{[k]}_\alpha z^\alpha.$$

A natural parametrization is

$t_\alpha^{[k]} = M^{-k}(\alpha + \tau),$

where $\tau\in\mathbb R^s$ is a shift parameter that later enters reproduction theorems (Charina et al., 2013).

The stationary case is recovered when $a^{[k]}$ is independent of $a^{[k]}$0. In classical non-stationary theory, the dilation matrix is fixed and only the mask varies. A broader generalization, called multiple subdivision, allows both the mask and the dilation matrix to vary with the level through a finite family

$a^{[k]}$1

so that a refinement path is a sequence $a^{[k]}$2. This setting contains stationary schemes as the special case $a^{[k]}$3, and certain non-stationary schemes as the special case of a fixed dilation matrix with varying masks (Charina et al., 2018).

Convergence is usually formulated by requiring that refined data sample a continuous or $a^{[k]}$4 limit function. In the regular tensor-product setting, one asks for

$a^{[k]}$5

For schemes on arbitrary-topology meshes, the definition near extraordinary vertices or faces is localized to rings around the irregular element and asks that the corresponding surface patches converge to a point $a^{[k]}$6 (Conti et al., 2017).

2. Generation, reproduction, and symbol calculus

A central distinction in nonstationary subdivision is between generation and reproduction of function spaces. For exponential-polynomial spaces

$a^{[k]}$7

generation means that sampled data from $a^{[k]}$8 are mapped to a limit in the same space, whereas reproduction means that the original function is recovered exactly in the limit. For non-singular schemes, step-wise reproduction and reproduction are equivalent (Charina et al., 2013).

The multivariate characterization of reproduction is algebraic. For a non-singular scheme, $a^{[k]}$9-reproduction holds if and only if there exists a shift $a(r)$0 such that

$a(r)$1

for all $a(r)$2 and all $a(r)$3. Here $a(r)$4, and $a(r)$5 is the multivariate falling-factorial polynomial

$a(r)$6

At the nontrivial points $a(r)$7, the derivatives must vanish; at $a(r)$8, they must satisfy a normalization determined by the parametrization shift (Charina et al., 2013).

In the univariate binary case, analogous conditions are formulated at the level-dependent nodes $a(r)$9. A non-singular scheme generates an exponential-polynomial space if

$M\in\mathbb Z^{s\times s}$0

and reproduces it if, in addition,

$M\in\mathbb Z^{s\times s}$1

with $M\in\mathbb Z^{s\times s}$2 determined by the primal or dual parametrization (Conti et al., 2014).

A major constructive theme is the conversion of approximating non-stationary schemes into interpolatory ones. Given a levelwise approximating symbol $M\in\mathbb Z^{s\times s}$3, one seeks a polynomial $M\in\mathbb Z^{s\times s}$4 solving the generalized Bezout equation

$M\in\mathbb Z^{s\times s}$5

and then defines

$M\in\mathbb Z^{s\times s}$6

The resulting interpolatory scheme reproduces the same exponential-polynomial space $M\in\mathbb Z^{s\times s}$7 as the original approximating scheme whenever it converges (Conti et al., 2010).

This symbol calculus also underlies the family of exponential pseudo-splines. There one writes

$M\in\mathbb Z^{s\times s}$8

where $M\in\mathbb Z^{s\times s}$9 is a normalized exponential B-spline symbol and $$f^{[k+1]}_\alpha=\sum_{\beta\in \mathbb Z^s} a^{[k]}_{\alpha-M\beta}\, f^{[k]}_\beta, \qquad \alpha\in\mathbb Z^s,\ k\ge 0,$$0 is a minimal-support correction determined by Hermite interpolation conditions. In the limit $$f^{[k+1]}_\alpha=\sum_{\beta\in \mathbb Z^s} a^{[k]}_{\alpha-M\beta}\, f^{[k]}_\beta, \qquad \alpha\in\mathbb Z^s,\ k\ge 0,$$1, these symbols converge to the classical polynomial pseudo-spline symbols, so exponential pseudo-splines are genuinely a non-stationary counterpart of stationary pseudo-splines (Conti et al., 2014).

3. Convergence, approximate sum rules, and matrix criteria

In stationary subdivision, polynomial reproduction is tied to exact sum rules. In the non-stationary case, the corresponding notion is that of approximate sum rules. For a univariate scheme, approximate sum rules of order $$f^{[k+1]}_\alpha=\sum_{\beta\in \mathbb Z^s} a^{[k]}_{\alpha-M\beta}\, f^{[k]}_\beta, \qquad \alpha\in\mathbb Z^s,\ k\ge 0,$$2 require

$$f^{[k+1]}_\alpha=\sum_{\beta\in \mathbb Z^s} a^{[k]}_{\alpha-M\beta}\, f^{[k]}_\beta, \qquad \alpha\in\mathbb Z^s,\ k\ge 0,$$3

where

$$f^{[k+1]}_\alpha=\sum_{\beta\in \mathbb Z^s} a^{[k]}_{\alpha-M\beta}\, f^{[k]}_\beta, \qquad \alpha\in\mathbb Z^s,\ k\ge 0,$$4

These conditions play, in the non-stationary theory, the role played by exact sum rules in the stationary theory (Conti et al., 2014).

Approximate sum rules interact with asymptotic comparison principles. Two schemes are asymptotically equivalent if the levelwise differences are summable, while they are asymptotically similar if

$$f^{[k+1]}_\alpha=\sum_{\beta\in \mathbb Z^s} a^{[k]}_{\alpha-M\beta}\, f^{[k]}_\beta, \qquad \alpha\in\mathbb Z^s,\ k\ge 0,$$5

For bounded local univariate schemes reproducing constants, asymptotical similarity to a convergent stationary scheme is sufficient for convergence; asymptotical equivalence is not required. The proof proceeds through the difference scheme

$$f^{[k+1]}_\alpha=\sum_{\beta\in \mathbb Z^s} a^{[k]}_{\alpha-M\beta}\, f^{[k]}_\beta, \qquad \alpha\in\mathbb Z^s,\ k\ge 0,$$6

and a contraction condition called Condition A (Conti et al., 2014).

The matrix approach packages regularity into finite-dimensional spectral data. For level-dependent masks satisfying sum rules of order $$f^{[k+1]}_\alpha=\sum_{\beta\in \mathbb Z^s} a^{[k]}_{\alpha-M\beta}\, f^{[k]}_\beta, \qquad \alpha\in\mathbb Z^s,\ k\ge 0,$$7, one forms finite matrices

$$f^{[k+1]}_\alpha=\sum_{\beta\in \mathbb Z^s} a^{[k]}_{\alpha-M\beta}\, f^{[k]}_\beta, \qquad \alpha\in\mathbb Z^s,\ k\ge 0,$$8

and restricts them to a common invariant subspace $$f^{[k+1]}_\alpha=\sum_{\beta\in \mathbb Z^s} a^{[k]}_{\alpha-M\beta}\, f^{[k]}_\beta, \qquad \alpha\in\mathbb Z^s,\ k\ge 0,$$9 orthogonal to polynomial left eigenvectors. The key spectral quantity is the joint spectral radius

$c^{(r+1)} = S_{a(r)}c^{(r)}.$0

which is norm-independent. For level- and parameter-dependent schemes with linear dependence on a parameter $c^{(r+1)} = S_{a(r)}c^{(r)}.$1, every restricted matrix is a convex combination of endpoint matrices, and

$c^{(r+1)} = S_{a(r)}c^{(r)}.$2

Under the sum-rule assumptions, the scheme is $c^{(r+1)} = S_{a(r)}c^{(r)}.$3-convergent if

$c^{(r+1)} = S_{a(r)}c^{(r)}.$4

and its Hölder exponent satisfies

$c^{(r+1)} = S_{a(r)}c^{(r)}.$5

This reduces an infinite non-stationary analysis to a finite endpoint problem (Charina et al., 2015).

A parallel multivariate theory with $c^{(r+1)} = S_{a(r)}c^{(r)}.$6 combines approximate sum rules, asymptotic similarity, and joint spectral radius. If the symbols satisfy approximate sum rules of order $c^{(r+1)} = S_{a(r)}c^{(r)}.$7 and

$c^{(r+1)} = S_{a(r)}c^{(r)}.$8

then the scheme is $c^{(r+1)} = S_{a(r)}c^{(r)}.$9-convergent, with Hölder exponent bounded below by

$$a^{[k]}(z)=\sum_{\alpha\in\mathbb Z^s} a^{[k]}_\alpha z^\alpha.$$0

Under stability of the stationary limit function and sufficiently fast decay of the defects, the bound is exact; this yields, for example, the result that generalized Daubechies-type wavelets have the same Hölder regularity as the corresponding classical Daubechies wavelets (Charina et al., 2014).

4. Nonlinear, non-uniform, and topology-sensitive extensions

Nonstationarity is not confined to linear uniform refinement. On arbitrary-topology meshes, a general convergence and normal-continuity theory near extraordinary vertices and faces compares a rotationally symmetric non-stationary scheme with a stationary comparison scheme. In regular regions one assumes asymptotic equivalence or asymptotic equivalence of order $$a^{[k]}(z)=\sum_{\alpha\in\mathbb Z^s} a^{[k]}_\alpha z^\alpha.$$1; in irregular regions one requires exponential decay

$$a^{[k]}(z)=\sum_{\alpha\in\mathbb Z^s} a^{[k]}_\alpha z^\alpha.$$2

and for normal continuity the stronger condition

$$a^{[k]}(z)=\sum_{\alpha\in\mathbb Z^s} a^{[k]}_\alpha z^\alpha.$$3

Under these hypotheses, convergence and normal continuity at extraordinary elements follow (Conti et al., 2017).

A different direction is non-uniform corner-cutting. The proposed non-uniform corner-cutting scheme defines level- and location-dependent masks through exponential-polynomial reproduction, with shape parameters $$a^{[k]}(z)=\sum_{\alpha\in\mathbb Z^s} a^{[k]}_\alpha z^\alpha.$$4 chosen from the data. The masks converge to Chaikin’s quadratic B-spline mask

$$a^{[k]}(z)=\sum_{\alpha\in\mathbb Z^s} a^{[k]}_\alpha z^\alpha.$$5

so the method is asymptotically equivalent to the classical scheme. It is proved stable, convergent, and $$a^{[k]}(z)=\sum_{\alpha\in\mathbb Z^s} a^{[k]}_\alpha z^\alpha.$$6, and it attains approximation order $$a^{[k]}(z)=\sum_{\alpha\in\mathbb Z^s} a^{[k]}_\alpha z^\alpha.$$7, whereas the classical methods attain second order accuracy (Jeong et al., 2020).

Nonlinearity can also be combined with nonstationarity to reproduce polynomial data on non-uniform grids without explicit grid information. A uniform non-linear non-stationary 4-point scheme infers local grid ratios through annihilation operators, approaches the linear uniform 4-point Lagrange scheme as $$a^{[k]}(z)=\sum_{\alpha\in\mathbb Z^s} a^{[k]}_\alpha z^\alpha.$$8, reproduces $$a^{[k]}(z)=\sum_{\alpha\in\mathbb Z^s} a^{[k]}_\alpha z^\alpha.$$9 on suitable initial grids, and is $t_\alpha^{[k]} = M^{-k}(\alpha + \tau),$0 for every $t_\alpha^{[k]} = M^{-k}(\alpha + \tau),$1 (López-Ureña, 2024).

At the opposite extreme, non-stationary schemes with growing mask support can be used to generate compactly supported $t_\alpha^{[k]} = M^{-k}(\alpha + \tau),$2 refinable functions. Using box-spline masks of increasing smoothing content, one obtains multivariate Up-like functions as basic limit functions of a non-stationary scheme. In one dimension the support can be reduced to

$t_\alpha^{[k]} = M^{-k}(\alpha + \tau),$3

which approaches length $t_\alpha^{[k]} = M^{-k}(\alpha + \tau),$4 as $t_\alpha^{[k]} = M^{-k}(\alpha + \tau),$5 increases, while preserving $t_\alpha^{[k]} = M^{-k}(\alpha + \tau),$6 regularity (Charina et al., 2022).

A useful comparison point is provided by stationary nonlinear alternatives. One such interpolatory stationary nonlinear scheme reproduces

$t_\alpha^{[k]} = M^{-k}(\alpha + \tau),$7

for suitable $t_\alpha^{[k]} = M^{-k}(\alpha + \tau),$8 without prior knowledge of the conic parameter, whereas linear non-stationary schemes typically require masks tailored to $t_\alpha^{[k]} = M^{-k}(\alpha + \tau),$9. This does not replace non-stationary theory, but it clarifies one of its design motivations: exact reproduction of conics and exponential-polynomial data with geometry-dependent refinement (Donat et al., 2018).

5. Multiresolution, attractors, and scale-dependent geometry

Nonstationary subdivision also supplies the prediction machinery for multiresolution representations. In a Harten-type framework with spaces $\tau\in\mathbb R^s$0, prediction operators $\tau\in\mathbb R^s$1, decimation operators $\tau\in\mathbb R^s$2, and reconstructions $\tau\in\mathbb R^s$3, subdivision is naturally interpreted as prediction: $\tau\in\mathbb R^s$4 A central result is that one may reverse the classical order of construction: choose the prediction first and then define a compatible decimation by minimization,

$\tau\in\mathbb R^s$5

This is especially useful for nonstationary predictors such as exponential B-spline schemes, which do not naturally fit the point-value or cell-average settings (López-Ureña, 2019).

The connection with fractal geometry is obtained through sequences of function systems. For a non-stationary scheme with masks $\tau\in\mathbb R^s$6, one constructs a sequence of function systems $\tau\in\mathbb R^s$7 and studies forward and backward trajectories. The decisive result is that backward trajectories recover non-stationary subdivision limits: for a non-stationary $\tau\in\mathbb R^s$8-convergent scheme, the backward trajectories of the associated sequence of function systems converge to a unique attractor whose first components constitute the subdivision limit curve. Unlike stationary IFS attractors, these attractors may have different structures at different scales (Dyn et al., 2016).

A recent multiscale formulation is the pyramid transform based on nonstationary subdivision. At level $\tau\in\mathbb R^s$9, one uses a refinement mask $a^{[k]}$0 and a matching decimation mask $a^{[k]}$1: $a^{[k]}$2 The detail coefficients satisfy a scale decay estimate of order $a^{[k]}$3 for sampled $a^{[k]}$4 data, and both decomposition and reconstruction are stable under levelwise operator bounds. In geometric examples, schemes reproducing $a^{[k]}$5 yield near-zero details on exact circles and localized large details on perturbed arcs, making the transform effective for circle detection and anomaly localization in planar data (Landau et al., 13 Jul 2025).

6. Limitations, parameter effects, and structural constraints

Nonstationarity enlarges the design space, but it does not remove structural constraints. A necessary criterion for a compactly supported function $a^{[k]}$6 to be generated by a level-dependent scheme is that the zero set of its Fourier transform satisfies

$a^{[k]}$7

together with the scaling-periodicity relation

$a^{[k]}$8

This excludes certain compactly supported functions; for example,

$a^{[k]}$9

cannot be generated by any nonstationary subdivision scheme because its Fourier transform is the Bessel function $a^{[k]}$00, whose zeros are incompatible with the structured zero sets above (Charina et al., 2015).

A second important limitation is analytical rather than existential. Many convergence and smoothness theorems require either approximate sum rules, asymptotic similarity to a stationary comparison scheme, or explicit spectral conditions on restricted transition matrices. The non-stationary theory therefore remains strongly tied to stationary reference schemes, even when the masks do not converge in a summable way (Charina et al., 2014).

At the same time, nonstationarity can improve classical behavior. In parameter-dependent families, the actual parameter range may materially change the joint spectral radius and hence the guaranteed Hölder exponent. In geometric design, non-stationary schemes reproduce circles, conics, exponential splines, and tension-controlled shapes that stationary polynomial schemes cannot match as directly. In surface refinement, sufficiently rapid decay of level-dependent irregular stencils preserves convergence and normal continuity near extraordinary elements. In wavelet construction, growing-support schemes generate compactly supported $a^{[k]}$01 refinable functions with unusually small supports (Charina et al., 2015, Conti et al., 2017, Charina et al., 2022).

Taken together, these developments define nonstationary subdivision schemes less as a single technique than as a broad level-dependent paradigm: exact reproduction is encoded by symbol conditions, convergence and regularity are controlled by approximate sum rules and spectral radii, and geometric flexibility arises from the ability to vary masks, parameters, support sizes, and even dilation matrices across levels.