Lecture II Univariate subdivision schemes and their analysis by the matrix approach

(1)

Lecture II

Univariate subdivision schemes and their analysis by the matrix approach

Ulrich Reif

Technische Universit¨at Darmstadt

Bertinoro, May 18, 2010

(2)

Overview

Univariate subdivision as a binary tree of products of square matrices.

Contractivity of matrices ↔ continuity

Joint spectral radius matrices ↔ H¨older continuity Strategies for computing the JSR.

Examples

(3)

Approach 1: Laurent Series

Sequence of data points at level k ∈N0, f^k ={f_i^k}_i∈_Z.

Identify with Laurent series f^k(z) :=X

i∈Z

f_i^kzⁱ.

Subdivision represented by symbol a, f^k+1(z) =a(z)f^k(z²).

Study properties of product function

a(z)a(z²)a(z⁴)· · ·a(z²^k) → exponential growth in k.

(4)

Approach 2: Matrices (global)

Sequence of data points at level k ∈N0, f^k ={f_i^k}_i∈_Z

Subdivision represented byinfinte matrixA∈R^Z^×^Z, f^k =Af^k−1 =A²f^k−2=· · ·=A^kf⁰

Problem: How to study properties ofA^k?

(5)

Approach 2: Matrices (local)

Idea: Reduce infinite sequencef⁰ of initial data points to the vector F⁰ defining the limit functionf on [0,1],

(6)

Approach 2: Matrices (local)

k = 0

k = 1

k = 2

k = 3

k = 4

(7)

Approach 2: Matrices (local)

k = 0 k = 1 k = 2 k = 3 k = 4

(8)

Approach 2: Matrices (local)

0 1

k = 0 k = 1 k = 2 k = 3 k = 4

(9)

Approach 2: Matrices (local)

F⁰ ={f_i⁰}_i=1:N₀, N0 =N F¹ ={f_i¹}_i=1:N₁, N₁ =N+ 1 F² ={f_i²}_i=1:N₂, N₂ =N+ 3 F^k ={f_i^k}_i=1:N_k, N_k =N+ 2^k −1 Local subdivision represented by finitematrices A^k,

F^k =A^kF^k⁻¹=A^kA^k−1F^k−2 =· · ·=A^kA^k−1· · ·A¹F⁰ Problem: How to study properties of the product matrix?

(10)

Approach 2: Matrices (local, square)

F⁰ ={f_i⁰}_i=1:N₀, N₀ =N F¹ ={f_i¹}_i=1:N₁, N1 =N+ 1

Idea: Partition F¹ into two sub-vectors of length N, F¹ = [f₁¹,f₂¹, . . . ,f_N¹,f_N+1¹ ]

F_`¹ = [f₁¹,f₂¹, . . . ,f_N¹] F_r¹ = [f₂¹, . . . ,f_N¹,f_N+1¹ ]

Local subdivision represented by a pair (S_`,Sr) of square matrices, F_`¹ =S_`F⁰ definingf on [0,1/2]

F_r¹ =S_rF⁰ definingf on [1/2,1]

(11)

Approach 2: Matrices (local, square)

F⁰ ={f_i⁰}_i=1:N₀, N0 =N F¹ ={f_i¹}_i=1:N₁, N1 =N+ 1

Idea: Partition F¹ into two sub-vectors of length N, F¹ = [f₁¹,f₂¹, . . . ,f_N¹,f_N+1¹ ]

F_`¹ = [f₁¹,f₂¹, . . . ,f_N¹] F_r¹ = [f₂¹, . . . ,f_N¹,f_N+1¹ ]

Local subdivision represented by a pair (S0,S1) of square matrices, F₀¹ =S0F⁰ definingf on [0,1/2]

F₁¹ =S₁F⁰ definingf on [1/2,1]

(12)

Approach 2: Matrices (local, square)

F=F⁰

S₀ F S₁ F

S₀ S₀ F S₁ S₀ F S₀ S₁ F S₁ S₁ F

S0S 0S

0F S 1S

0S 0F S

0S 1S

0F S 1S

1S 0F S

0S 0S

1F S 1S

0S 1F S

0S 1S

1F S 1S

1S 1F

(13)

Approach 2: Matrices (local, square)

At levelk, there are 2^k sub-intervals, indexed by I = [i1, . . . ,ik]∈ {0,1}^k.

The binary number 0.i₁· · ·i_k is the left end-point of the sub-interval corresponding toI.

The vectorF_I^k of data defining the limit functionf on the sub-interval with index I is given by

F_I^k =SIF⁰, where SI :=Si_k· · ·Si1. Analyze binary tree of products of (S₀,S₁).

(14)

The difference scheme

Let

∆ :=







−1 1 0 · · · 0 . .. . ..

0 · · · 0 −1 1







denote thedifference matrix.

The matrices D = (D₀,D₁) representing the difference scheme

∆F_i¹ =D_i∆F⁰ must satisfy

∆S_i =D_i∆.

A solution exists and is unique iff the rows of S₀,S₁ sum up to 1.

(15)

The difference scheme

Let

∆ :=







−1 1 0 · · · . .. ...

· · · 0 −1 1





, ∆⁻¹ :=







0 1 · · · 1 . .. ...

0 · · · 0 1







T

denote thedifference matrix and thesummation matrix, resp.

The matrices D = (D0,D1) representing the difference scheme

∆F_i¹ =S_i∆F⁰ must satisfy

∆S_i =D_i∆.

A solution exists and is unique iff the rows of S₀,S₁ sum up to 1, Di = ∆Si∆⁻¹.

(16)

The difference scheme

The differences at level k are given by

∆F_I^k =D_I∆F⁰, D_I :=D_i_k· · ·D_i₁, whereI = [i₁, . . . ,i_k]∈ {0,1}^k.

The subdivision schemeS = (S0,S1) generates aC⁰-limit function iff the difference scheme D= (D₀,D₁) is contractive, i.e., iff

DI = 0 for any infinite sequenceI ∈ {0,1}^N.

The difference scheme is not contractiveif

%(D_I)≥1 for some index vector I.

The difference scheme is contractiveif there existsN ∈Nsuch that kD_Ik<1 for all index vectorsI of length N.

(17)

The difference scheme

The difference scheme is contractiveif there existsN ∈Nsuch that kD_Ik<1 for all index vectorsI of length N.

(18)

The difference scheme

The difference scheme is contractiveif there existsN ∈Nsuch that kD k

(19)

The four-point scheme I

0 1 0 0 −w w’ w’ −w

Subdivision scheme S:

S₀ =







0 1 0 0 0 0

−ω ω⁰ ω⁰ −ω 0 0

0 0 1 0 0 0

0 −ω ω⁰ ω⁰ −ω 0

0 0 0 1 0 0

0 0 −ω ω⁰ ω⁰ −ω





 , S₁=







−ω ω⁰ ω⁰ −ω 0 0

0 0 1 0 0 0

0 −ω ω⁰ ω⁰ −ω 0

0 0 0 1 0 0

0 0 −ω ω⁰ ω⁰ −ω

0 0 0 0 0 1







(20)

The four-point scheme I

0 1 0 0 −w w’ w’ −w

Difference scheme D:

D0 =







−ω ¹₂ ω 0 0

ω ¹₂ −ω 0 0

0 −ω ¹₂ ω 0 0 ω ¹₂ −ω 0 0 0 −ω ¹₂ ω







, D1 =







ω ¹₂ −ω 0 0

0 −ω ¹₂ ω 0 0 ω ¹₂ −ω 0 0 0 −ω ¹₂ ω 0 0 ω ¹₂ −ω







(21)

The four-point scheme I

0 1 0 0 −w w’ w’ −w

Divided difference scheme ¯D:= 2D:

D¯0 = 2







−ω ¹₂ ω 0 0

ω ¹₂ −ω 0 0

0 −ω ¹₂ ω 0 0 ω ¹₂ −ω 0 0 0 −ω ¹₂ ω







, D¯1 = 2







ω ¹₂ −ω 0 0

0 −ω ¹₂ ω 0 0 ω ¹₂ −ω 0 0 0 −ω ¹₂ ω 0 0 ω ¹₂ −ω







(22)

The four-point scheme I

0 1 0 0 −w w’ w’ −w

Difference scheme D² of divided difference scheme:

D₀² = 2







2ω 2ω 0 0

−ω ω⁰−ω −ω 0

0 2ω 2ω 0

0 −ω ω⁰−ω −ω







, D₁² = 2







−ω ω⁰−ω −ω 0

0 2ω 2ω 0

0 −ω ω⁰−ω −ω

0 0 2ω 2ω







(23)

The four-point scheme I

The FPS generatesC¹-limit functions iff the difference schemeD² of the divided difference scheme is contractive.

Determining the maximal set (0, ω_sup) providingC¹ is a challenge:

I ’87, based on leveln= 2,ωsup≥ ¹₈=.125

I ’92, based on leveln= 3,ωsup≥

√5−1 8 ≈.155

I ’96, based on leveln= 22,ω_sup≥.188

I ’06, based on refined analytic approach,

ωsup:= 27 + 3√ 105^2/3

−6 12 27 + 3√

105^1/3

≈0.192729.

(24)

The four-point scheme I

The FPS generatesC¹-limit functions iff the difference schemeD² of the divided difference scheme is contractive.

Determining the maximal set (0, ω_sup) providingC¹ is a challenge:

I ’87, based on leveln= 2,ωsup≥ ¹₈=.125

I ’92, based on leveln= 3,ωsup≥

√5−1 8 ≈.155

I ’96, based on leveln= 22,ω_sup≥.188

I ’06, based on refined analytic approach,

ω_sup:= 27 + 3√ 105^2/3

−6 12 27 + 3√

105^1/3

≈0.192729.

(25)

Breadth-first vs. depth-first search

0 0.2 0.4 0.6 0.8 1

0 1 2 3 4 5 6

D is contractive if there exists a level N with contractive nodes, kD_Ik<1 for allI of lengthN,

i.e., if abreadth first search terminates.

(26)

Breadth-first vs. depth-first search

0 0.2 0.4 0.6 0.8 1

0 1 2 3 4 5

6 0 0.2 0.4 0.6 0.8 1

0 1 2 3 4 5 6

D is contractive if there exists a proper subtreeT with contractive nodes, i.e., if adepth first search terminates.

For ω= 0.188, reduction from 4,000,000to 159matrices.

(27)

Bad news I

The problem of checking a pair of matrices (D₀,D₁) for contractivity is undecidable.

(28)

H¨ older continuity and the joint spectral radius

The rate of convergence of kD_Ik →0 as #I → ∞determines the H¨older continuity of the limit functionf.

Thejoint spectral radius of (D0,D1) is defined by jsr(D0,D1) := sup

n∈N

sup

I∈{0,1}ⁿ

pn

%(DI).

The limit function f isC⁰ if jsr(D₀,D₁)<1.

The limit function f isC^0,α if jsr(D₀,D₁)<2^−α.

(29)

Basic observations

For anyI ∈ {0,1}ⁿ, pn

%(D_I)≤jsr(D₀,D₁).

For any norm,

jsr(D0,D1)≤ max

I∈{0,1}ⁿ

pn

kD_Ik.

There exists a normk · k∗ onR^d with jsr(D₀,D₁) = max

kD₀k_∗,kD₁k_∗ .

For symmetric subdivision schemes and (2×2)-matricesD₀,D₁, jsr(D0,D1) = max

%(D0),p

%(D0D1) .

(30)

Corner cutting

1−w w w 1−w

S0 =





1−ω ω 0

ω 1−ω 0

0 1−ω ω



, S1 =





ω 1−ω 0

0 1−ω ω

0 ω 1−ω





(31)

Corner cutting

1−w w w 1−w

D₀ =

1−2ω 0

ω ω

, D₁ =

ω ω 0 1−2ω

The joint spectral radius is given by jsr(D₀,D₁) = maxn

1−2ω, ω/2+

q

ω−7ω²/4o

0 0.1 0.2 0.3 0.4 0.5

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

(32)

The four-point scheme II

For the four-point scheme with parameter

ω∈[0,1/20], it is

jsr(D₀²,D₁²) =%(D₀²).

ω∈[1/10,2/10], it is jsr(D₀²,D₁²) =

q

%(D₀²D₁²).

0 0.05 0.1 0.15 0.2

0.4 0.5 0.6 0.7 0.8 0.9 1 1.1

The upper bound ωsup is obtained from solving %(D₀²D₁²) = 1.

(33)

Bad news II

Thefiniteness conjecture jsr(D₀,D₁) =pⁿ

%(D_I) for someI ∈ {0,1}ⁿ was disproven.

The jsr-problem is np-complete with respect to accuracy and dimension.

In general, the numerical computation of jsr(D0,D1) with accuracyε requires

O (dimD0)^1/ε operations.

(34)

Bad news II

Thefiniteness conjecture jsr(D₀,D₁) =pⁿ

%(D_I) for someI ∈ {0,1}ⁿ was disproven.

The jsr-problem is np-complete with respect to accuracy and dimension.

In general, the numerical computation of jsr(D0,D1) with accuracyε requires

O (dimD0)^1/ε operations.

(35)

Good news

So far, no counterexamples to the finiteness conjecture have been encountered in practice.

Robust algorithm for verifying jsr(D₀,D₁) =pⁿ

%(D_I)

for given I is available (implementation in progress).

(36)

Exact evaluation of the JSR

Determine a candidateI for the finiteness conjecture, jsr(D0,D1) =pⁿ

%(DI).

Let

D˜i = Di

jsr(D₀,D₁).

Verifying the conjecture is equivalent to showing jsr( ˜D₀,D˜₁) = 1.

(37)

Exact evaluation of the JSR - Method 1

Start with unit cube M⁰. If the recursion

M^k⁺¹= conv( ˜D0M^k,D˜1M^k) with stopping criterion

M^k⁺¹⊆M^k

terminates, then the conjecture is verified.

The setM^k defines the unit ball wrt. the optimal norm k · k_∗.

(38)

Exact evaluation of the JSR - Method 2

Conisder tree of matrix products with

I edges of type

D0,D1,{D_Iⁿ:n∈N0}

I set-valued nodes of type D={DJD_IⁿDK :n∈N0}

If a depth-first search with stopping criterion maxD∈DkDk<1

terminates, then the conjecture is verified.

(39)

Conclusion

The matrix approach provides an alternative to the Laurent series formalism.

From a theoretical point of view, both methods are equivalent.

For special purposes, one approach may be more efficient than the other.

In general, sharp results are beyondreach, and even good estimates may bevery hard to determine.

In practise, sharp results arewithinreach.