Spherical designs and QMC design sequences Valuable tools for - - PowerPoint PPT Presentation

Spherical designs and QMC design sequences – Valuable tools for numerical analysis

Ian H. Sloan

University of New South Wales, Sydney, Australia Joint work with J Brauchart, J Dick, EB Saff, YG Wang and R Womersley

QMC?

QMC = Quasi-Monte Carlo Quasi-Monte Carlo rules are integration rules in which all weights are equal. An important special case: spherical designs

Spherical designs

Definition: A spherical t-design on Sd ⊂ Rd+1 is a set XN := {x1, . . . , xN} ⊂ Sd such that 1 N

N

• j=1

p(xj) =

• Sd

p(x) d σd(x) ∀p ∈ Pt.

Here d σd(x) is normalised measure on Sd.

So XN is a spherical t-design if the QMC (i.e. equal weight) cubature rule with these points integrates exactly all polynomials of degree ≤ t.

A spherical 50-design

This is a Womersley spherical 50-design with 1302 ≈ 1

2512 points.

Spherical designs are good for integration

Spherical designs are tools for numerical integration. The following theorem (Hesse & IHS, 2005, 2006) shows a good rate of convergence for sufficiently smooth functions f: Theorem. Given s > d/2, there exists C(s, d) > 0 such that for every spherical t-design XN on Sd there holds

• 1

N

• x∈XN

f(x) −

• Sd f(x) d σd(x)
• ≤ C(s, d)

ts fHs. The result holds more generally, for all positive-weight rules that are exact for polynomials of degree up to t.

How many points for a spherical t-design?

It is known (Seymour & Zaslavsky, 1984) that for every t ≥ 1 (and for every

dimension of the sphere) there always exists a spherical design.

But how many points does a spherical t-design need?

There is no possible upper bound because ...

Delsarte, Goethals, Seidel (1977) established lower bounds of exact

• rder td:

N ≥       

• d+t/2

d

• +
• d+t/2 − 1

d

• ,

if t is even

• d+⌊t/2⌋

d

• ,

if t is odd Yudin (1997) established larger lower bounds, still of exact order td.

A constant ×td is enough points

It has long been conjectured that cdtd points is enough, for some cd > 0, but until very recently there was no proof. Recently Bondarenko, Radchenko and Viazovska (Annals of Mathematics, 2013) proved this important existence result. But what is the constant? The Bondarenko et al. constant is huge. For S2 Chen, Frommer and Lang (2011) proved that (t + 1)2 points is enough for all t up to 100. For S2, we believe that (t + 1)2 points is enough for all t. Even N ≈ 1

2(t + 1)2 seems to be enough (R. Womersley, private

communication).

QMC designs

• Definition. A sequence of point sets (XN) ⊂ Sd with N → ∞ is a

sequence of QMC designs for the Sobolev space Hs(Sd), for some s > d/2, if there exists c(s, d) > 0, such that for all f ∈ Hs(Sd)

• 1

N

• x∈XN

f(x) −

• Sd f(x) d σd(x)
• ≤ c(s, d)

N s/d fHs. This is the optimal rate of convergence in Hs(Sd)

K Hesse and IHS 2005 for d = 2, Hesse 2006 for general d.

The idea grew from properties of spherical designs.

Efficient spher. designs are QMC designs

Clearly, every sequence of spherical t-designs is a sequence of QMC designs for Hs(Sd), for all s > d/2, iff N ≍ td as t → ∞, since Theorem. Given s > d/2, there exists C(s, d) > 0 such that for every spherical t-design XN on Sd there holds

• 1

N

• x∈XN

f(x) −

• Sd f(x) d σd(x)
• ≤ C(s, d)

ts fHs. If N ≍ td this gives

• 1

N

• x∈XN

f(x) −

• Sd f(x) d σd(x)
• ≤ c(s, d)

N s/d fHs.

Are there other QMC designs?

There are many. Here’s one:

• Theorem. (J Brauchart, EB Saff, IH Sloan, R Womersley, Math Comp, 2014)

A sequence of N-point sets X∗

N that maximize the sum of pairwise

Euclidean distances is a sequence of QMC designs for H(d+1)/2(Sd). Thus for S2 the points that maximize the sum of Euclidean distances form a sequence of QMC designs for H3/2. To prove this and other things we need some machinery. But first:

The nested property of QMC designs

• Theorem. (Brauchart, Saff, IHS, Womersley, op. cit.)

Given s > d/2, let (XN) be a sequence of QMC designs for Hs(Sd). Then (XN) is a sequence of QMC designs for all coarser Hs′(Sd), i.e. for all s′ satisfying d/2 < s′ ≤ s.

This result isn’t trivial – for the smaller set Hs(Sd) we demand faster convergence.

So there is some upper bound on the admissible values of s: s∗ := sup{s : (XN) is a sequence of QMC designs for Hs}. We call s∗ the strength of the sequence of QMC designs (XN).

Generic QMC designs

If s∗ = +∞, we say the sequence of QMC designs is “generic”. Every sequence of spherical t-designs with N ≍ td as t → ∞ is a generic sequence of QMC designs.

We don’t know if there are other interesting examples.

The Sobolev space Hs(Sd)

With λℓ := ℓ (ℓ + d − 1) (λℓ is the ℓth eigenvalue of −∆∗

d),

Hs(Sd) =

• f ∈ L2(Sd) :

∞

• ℓ=0

Z(d,ℓ)

• k=1

(1 + λℓ)s

• fℓ,k
• 2

< ∞

• .

Thus H0(Sd) = L2(Sd). Here Laplace-Fourier coefficients

• fℓ,k = (f, Yℓ,k)L2(Sd) =
• Sd f(x)Yℓ,k(x) d σd(x).

Yℓ,k for k = 1, . . . , Z(d, ℓ) is an orthonormal set of spherical harmonics of degree ℓ: ∆∗

d Yℓ,k = −λℓYℓ,k.

Norms for Hs(Sd)

It is useful to allow also other equivalent norms for Hs(Sd): Let (a(s)

ℓ

)ℓ≥0 satisfy a(s)

ℓ

≍ (1 + λℓ)−s ≍ (1 + ℓ)−2s. Inner product and norm for f, g ∈ Hs(Sd) (f, g)Hs:=

∞

• ℓ=0

Z(d,ℓ)

• k=1

1 a(s)

ℓ

• fℓ,k

gℓ,k, fHs:=

• (f, f)Hs.

It is easily seen that Hs(Sd) is embedded in C(Sd) iff s > d/2.

Worst case error

Let Q[XN](f) := 1 N

N

• j=1

f(xj) ≈

• Sd f(x) d σd(x).

Then the worst case error in Hs(Sd) is defined by: wce(Q[XN]; Hs(Sd)) := sup

• Q[XN](f) − I(f)
• : f ∈ Hs(Sd), fHs ≤ 1
• ,

QMC designs defined in terms of WCE: A sequence (XN) of N point configurations on Sd is a sequence of QMC designs for Hs iff there exists c(s, d) > 0 such that wce(Q[XN]; Hs(Sd)) ≤ c(s, d) N s/d .

Reproducing kernel for Hs(Sd), s > d/2

K(s)(x, y)=

∞

• ℓ=0

Z(d,ℓ)

• k=1

a(s)

ℓ

Yℓ,k(x)Yℓ,k(y) =

∞

• ℓ=0

a(s)

ℓ

Z(d, ℓ)P (d)

ℓ

(x · y), where P (d)

ℓ

Legendre polynomial associated with Sd. It is a zonal kernel: i.e. K(s)(x, y) = K(s)(x · y). The reproducing kernel properties are easily verified: K(s)(·, x)∈ Hs(Sd), x ∈ Sd, (f, K(s)(·, x))Hs= f(x), x ∈ Sd, f ∈ Hs(Sd).

WCE in terms of the reproducing kernel

Recall: For a(s)

ℓ

≍ (1 + ℓ)−2s,

K(s)(x · y) :=

∞

• ℓ=0

a(s)

ℓ Z(d,ℓ)

• k=1

Yℓ,k(x)Yℓ,k(y) =

∞

• ℓ=0

a(s)

ℓ

Z(d, ℓ)P (d)

ℓ

(x · y)

is a reproducing kernel for Hs(Sd). Theorem

• wce(Q[XN]; Hs(Sd))

2 = 1 N 2

N

• j=1

N

• i=1

K(s)(xj · xi) . Here K(s)(x · y) :=

∞

• ℓ=1

a(s)

ℓ

Z(d, ℓ)P (d)

ℓ

(x · y) = K(s)(x · y) − a(s)

0 .

Optimal QMC designs

Recall: K(s)(x · y) := K(s)(x · y) − a(s)

0 .

Let X∗

N = {x∗ 1, · · · , x∗ N}, for N = 2, 3, 4, . . . be a sequence of

minimizers of

N

• j=1

N

• i=1

K(s)(xj · xi),

• Theorem. (op cit) For all s.d/2 there exists c(s, d) > 0 such that for all

N ≥ 2 wce(Q[X∗

N]; Hs(Sd)) ≤ c(s, d)

N s/d . Consequently, (X∗

N) is a sequence of QMC designs for Hs(Sd).

Proof: Spherical designs with N ≍ td exist, and satisfy the bounds in the theorem for all s > d/2. The minimizer for a particular s > d/2 must be at least as good.

QMC designs from distance kernels

Let V (Sd) :=

• Sd
• Sd |x − y| d σd(x) d σd(y). It can be shown that

K(d+1)/2(x, y) := 2V (Sd) − |x − y| is a reproducing kernel for H(d+1)/2(Sd), and that correspondingly K(d+1)/2(x, y) = V (Sd) − |x − y| . Thus wce(Q[XN]; H(d+1)/2(Sd)) =  V (Sd) − 1 N 2

N

• j=1

N

• i=1

|xj − xi|  

1/2

is the corresponding WCE in H(d+1)/2(Sd).

• Corollary. A sequence of N-point sets X∗

N that maximize the sum of

Euclidean distances is a sequence of QMC designs for H(d+1)/2(Sd).

Generalized distance kernels

In the same way, a kernel for Sobolev spaces with s ∈ (d/2, d/2 + 1) is given by (Hubbert & Baxter 2012, Brauchart & Womersley 201?) K(s)

gd (x, y) := 2V2s−d(Sd) − |x − y|2s−d

where V2s−d := V2s−d(Sd) :=

• Sd
• Sd |x − y|2s−d d σd(x) d σd(y).

A reproducing kernel for s > d/2 + 1

For d/2 + L < s < d/2 + L + 1, with L = 1, 2, · · · , a version of the reproducing kernel for Hs(Sd) is (Brauchart & Womersley 201?) K(s)

gd (x, y) :=

• 1 − (−1)L+1

Vd−2s(Sd) +QL(x · y) + (−1)L+1 |x − y|2s−d, where QL(x · y) :=

L

• ℓ=1

βℓ P (d)

ℓ

(x · y), βℓ = · · · . This is important because it gives us a practical tool for computing worst-case errors for (almost) all values of s.

QMC designs are better than average

We might wonder: perhaps ALL sequences of point sets are QMC designs? No, because: Theorem Given s > d/2,

• E
• wce(W (

Q(S[), PN)XN]; Hs(Sd))

2

• = b(s, d)

N 1/2 for some constant b(s, d) > 0, where the points x1, . . . , xN are independently and uniformly distributed on Sd. Hence for s > d/2 optimal order quadrature error O(N −s/d) is not attained by random points. On the other hand ...

Randomized equal area points

. Suppose we have a sequence (DN = {Dj,N, . . . , DN,N}) of equal area partitions of Sd with small diameter: ∪N

j=1Dj,N= Sd,

Dj,N ∩ Dk,N = ∅, j = k, σd(Dj,N) = 1/N ∀j, diamDj,N ≤ c/N 1/d ∀j. This time we distribute our N random points in a stratified way, with

• ne point to each Dj,N. The result is very different:

Theorem (op cit.) Let XN = {x1,N, . . . , xN,N}, where xj,N is chosen randomly from Dj,N with respect to uniform measure on Dj,N. Then, for d/2 < s < d/2 + 1, β′ N s/d ≤

• E
• wce(W (

Q(S[), PN)XN]; Hs(Sd))

2 ≤ β N s/d . N ≍ td. Thus on average, randomized equal area points will form a QMC design sequence for d/2 < s < d/2 + 1.

(But not for s > d/2 + 1.)

Candidates for sequences of QMC designs

For the sphere S2, some plausible candidates are:

Random points, uniformly distributed on the sphere. NO Equal area points (cf. Rakhmanov-S-Zhou). Fekete points which maximize the determinant (cf.

Sloan-Womersley).

• Log. energy points and Coulomb energy points, which minimize

N

• j=1

N

• i=1

log 1 |xj − xi|,

N

• j=1

N

• i=1

1 |xj − xi|.

(continued)

Generalized spiral points (cf. Rakhmanov-S-Zhou; Bauer), with

spherical coordinates (θj, φj) given by θj = cos−1(1 − 2j − 1 N ), φj = 1.8 √ Nθj mod 2π, j = 1, . . . ,

Distance points, which maximize

N

• j=1

N

• i=1

|xj − xi|.

Spherical t-designs with N = ⌈(t + 1)2/2⌉ + 1 points.

WCE for Hs(S2) and s = 1.5

10

1

10

2

10

3

10

4

10

−4

10

−3

10

−2

10

−1

10 Number of points N Random 1.18 N−0.52 Fekete 0.90 N−0.75 Equal area 0.93 N−0.75 Coulomb energy 0.90 N−0.75 Log energy 0.90 N−0.75 Generalized spiral 0.91 N−0.75 Distance 0.90 N−0.75 Spherical design 0.91 N−0.75

There are many overlapping curves

Random 1.18 N −0.52 Fekete 0.90 N −0.75 Equal area 0.93 N −0.75 Coulomb energy 0.90 N −0.75 Log energy 0.90 N −0.75 Generalized spiral 0.91 N −0.75 Distance 0.90 N −0.75 Spherical design 0.91 N −0.75

WCE for Hs(S2) and s = 2.5

10

1

10

2

10

3

10

4

10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10 Number of points N Random 2.25 N−0.53 Fekete 0.42 N−1.04 Equal area 0.82 N−0.98 Coulomb energy 1.03 N−1.23 Log energy 1.18 N−1.26 Generalized spiral 1.24 N−1.26 Distance 1.17 N−1.26 Spherical design 1.22 N−1.26

WCE for Hs(S2) and s = 3.5

10

1

10

2

10

3

10

4

10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10 Number of points N Random 4.92 N−0.52 Fekete 0.13 N−0.82 Equal area 1.80 N−0.95 Coulomb energy 0.17 N−1.06 Log energy 1.50 N−1.58 Generalized spiral 2.51 N−1.55 Distance 3.50 N−1.77 Spherical design 3.59 N−1.77

WCE for Hs(S2) and s = 4.5

10

1

10

2

10

3

10

4

10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10 10

1

Number of points N Random 10.22 N −0.52 Fekete 0.35 N −0.79 Equal area 4.09 N −0.95 Coulomb energy 0.29 N −1.01 Log energy 1.29 N −1.48 Generalized spiral 4.82 N −1.53 Distance 7.00 N −2.00 Spherical design 18.49 N −2.31

Integrating a smooth function

Franke function in C∞(S2) f

• x, y, z
• := 0.75 exp(−(9x − 2)2/4 − (9y − 2)2/4 − (9z − 2)2/4)

+0.75 exp(−(9x + 1)2/49 − (9y + 1)/10 − (9z + 1)/10) +0.5 exp(−(9x − 7)2/4 − (9y − 3)2/4 − (9z − 5)2/4) −0.2 exp(−(9x − 4)2 − (9y − 7)2 − (9z − 5)2)

• S2 f(x)dσ2(x) = 0.532865250084389...

And note: |Q[XN](f) − I(f)| ≤ wce(Q[XN]; Hs(Sd)) fHs. Thus we should see in the error a rate of decay of order N −s∗/2.

The Franke function

Integration errors for the Franke function

10

1

10

2

10

3

10

4

10

−13

10

−12

10

−11

10

−10

10

−9

10

−8

10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10 Number of points N Random 0.28 N−0.51 Fekete 0.03 N−0.90 Equal area 0.24 N−0.96 Coulomb energy 0.04 N−1.17 Log energy 0.47 N−1.68 Generalized spiral 0.68 N−1.71 Distance 3.23 N−2.14 Spherical design

Estimates of strength s∗ for d = 2

s∗ := sup{s : (XN) is a sequence of QMC designs for Hs(S2)}.

Table 1: Estimates of s∗ for d = 2

Point set s∗ Fekete 1.5 Equal area 2 Coulomb energy 2 Log energy 3 Generalized spiral 3 Distance 4 Spherical designs ∞

Thus for S2:

All the well known point set sequences we have looked at, except for random points, appear to be approx. spherical designs for Hs, for a significant range of values of s. But so far this has been proved only for point sets that maximize the sum of distances, and only for s up to 3/2. (And we have proved similar results for generalized sums of distances for 1 < s < 2 ). For the sum of distances points we have therefore proved that s∗ ≥ 3/2. But the experimental s∗ in this case is s∗ = 4 – there’s a big gap in our knowledge! QMC designs can be valuable tools for numerical integration, especially if s∗ is large.

Covering

Now for a new part of the story: Sequences of QMC designs have a good covering property. J Brauchart, J Dick, E Saff, I Sloan, Y Wang and R Womersley, Covering of spheres by spherical caps and worst-case error for equal weight cubature in Sobolev spaces, arXiv:1407.831v

Covering

the covering radius or mesh norm, or fill distance ρ(XN) := max

x∈Sd

min

1≤k≤N arccos(x · xk).

Covering by spherical caps

Now move just one point

Covering by spherical caps

Trivial lower bound

ρ(XN) ≥ cdN −1/d

because a spherical cap of this radius has surface measure of order

1 N .

We will say that the sequence XN) has the optimal covering property if ρ(XN) ≍ N −1/d

Covering and spherical designs

Reimer 2000, extending Yudin 1995, showed that a spherical t-design yields a covering radius ρ(XN) ≤ cdt−1. (The result holds also for any positive weight cubature rule.) Thus a sequence (XN) of spherical t-designs with N ≍ td yields

• ptimal covering:

ρ(XN) ≤ cdN −1/d.

Function spaces – generalise to p = 2

Recall: for f ∈ L2(Sd), f(x) =

∞

• ℓ=0

Z(d,ℓ)

• k=1
• fℓ,kYℓ,k(x),

where The Laplace-Fourier coefficients are defined by

• fℓ,k = (f, Yℓ,k)L2(Sd) =
• Sd f(x)Yℓ,k(x) d σd(x).

Yℓ,k for k = 1, . . . , Z(d, ℓ) is an L2 orthonormal set of spherical harmonics of degree ℓ: ∆∗

d Yℓ,k = −λℓYℓ,k.

with λℓ := ℓ (ℓ + d − 1) the ℓth eigenvalue of −∆∗

d.

The Sobolev space W s

p(Sd) Then f ∈ W s

p(Sd) if (1 − ∆∗ d)s/2f ∈ Lp(Sd), that is, if

fW s

p (Sd) :=

• ∞
• ℓ=0

Z(d,ℓ)

• k=0

(1 + λℓ)s/2 fℓ,kYℓ,k

• Lp(Sd)

< ∞. The case s = 0 gives W s

p(Sd) = Lp(Sd).

The main result

“The covering radius of a point set XN on Sd is upper bounded by a power of the worst-case error in a Sobolev space”: Theorem Let d ≥ 1, p ∈ [1, ∞] and s > d/p. For a positive integer N, let XN be an N-point set on Sd. Then ρ(XN) ≤ cs,d

• wce(Q[XN]; W s

p(Sd))

1/(s+d/q) , where 1/q + 1/p = 1.

Idea of the proof

The result is equivalent to wce(Q[XN]; W s

p(Sd)) ≥ const × ρ(XN)s+d/q,

We construct a “fooling function” f ∈ W s

p(Sd), whose support is a

spherical cap of radius ρ(XN) with no points of XN in its interior. Then obviously Q[XN]f = 0, implying wce(Q[XN]; W s

p(Sd)) ≥

|If| fW s

p (Sd)

. Finally, we show |If| ≥ cρd, and fW s

p (Sd) ≤ cρ(XN)−s+d/p.

Extending the definition of QMC designs

ρ(XN) ≤ cs,d

• wce(Q[XN]; W s

p(Sd))

1/(s+d/q) , So how fast can the WCE in W s

p(Sd) decay?

• Definition. For p ∈ [1, ∞] and s > d/p, a sequence of point sets

(XN) ⊂ Sd with N → ∞ is a sequence of QMC designs for the Sobolev space W s

p(Sd) if there exists c(s, d) > 0, such that

• 1

N

• x∈XN

f(x) −

• Sd f(x) d σd(x)
• ≤ c(s, d)

N s/d fW s

p (Sd)

∀f ∈ W s

p (Sd).

That is, it is a QMC design for W s

p(Sd) iff

wce(Q[XN]; W s

p(Sd)) ≤ c(s, d)

N s/d

QMC designs have optimal convergence

(XN) is a sequence of QMC designs iff wce(Q[XN]; W s

p(Sd)) ≤ c(s, d)

N s/d This is the optimal rate of convergence in W s

p(Sd)

– for p = 2 proved by Hesse & IHS 2005 for d = 2, then Hesse 2006 for general d; – for general p ∈ (0, ∞) proved by Brandolini et al 2013.

QMC designs lead to good covering

Corollary of main theorem Let d ≥ 1 and p ∈ [1, ∞] . Assume that (XN) with N → ∞ is a QMC design sequence for W s

p(Sd) for some s > d/p. Then there

exist a constant c > 0 depending on d, s, p and the sequence (XN) but not on N such that ρ(XN) ≤ c N −βs/d for all N, where βs := 1/(1 + d/qs) with 1/q + 1/p = 1. In particular, if p = 1 (and thus β = 1), then the sequence of covering radii ρ(XN) has the optimal covering property, ρ(XN) ≍ N −1/d.

Thus a sequence of QMC designs for W s

1 (Sd) and some s > d

necessarily has the optimal covering property.

Compare with spherical designs

Reimer 2000, extending Yudin 1995, showed that a spherical t-design yields a covering radius ρ(XN) ≤ cdt−1. (The result holds also for any positive weight cubature rule.) Thus a sequence (XN) of spherical t-designs with N ≍ td yields

• ptimal covering:

ρ(XN) ≤ cdN −1/d. The difference in the present work is that we have replaced the assumption of polynomial exactness with an assumption of optimal rate of convergence of wce(Q[XN]; W1s(Sd)), for some s > d.

References

J Brauchart, E Saff, I Sloan, R Womersley, Optimal order quasi-Monte Carlo integration schemes on the sphere. Math Comp 2014. J Brauchart, J Dick, E Saff, I Sloan, Y Wang and R Womersley, Covering of spheres by spherical caps and worst-case error for equal weight cubature in Sobolev spaces, arXiv:1407.831v