Contour integration methods for self adjoint operators OTKR - - PowerPoint PPT Presentation

▶
contour integration methods for self adjoint operators
SMART_READER_LITE
LIVE PREVIEW

Contour integration methods for self adjoint operators OTKR - - PowerPoint PPT Presentation

Sep 03, 2023 136 likes •923 views

Prirodoslovno-matemati cki fakultet Matemati cki odsjek Sveu cili ste u Zagrebu Contour integration methods for self adjoint operators OTKR December 1922 2019 TU Vienna Luka Grubi si c Department of Mathematics,

slide-1
SLIDE 1

Prirodoslovno-matematiˇ cki fakultet Matematiˇ cki odsjek Sveuˇ ciliˇ ste u Zagrebu

Contour integration methods for self adjoint operators

OTKR December 19–22 2019 – TU Vienna Luka Grubiˇ si´ c Department of Mathematics, University of Zagreb luka.grubisic@math.hr

Joint work with Jay Goplakrishnan and Jeffrey Ovall

Luka Grubiˇ si´ c 1 / 36

slide-2
SLIDE 2

Outline

1 About the problem 2 About the method 3 Discrete resolvent estimates 4 Numerical experiments

Luka Grubiˇ si´ c 2 / 36

slide-3
SLIDE 3

About this talk

FEAST idea (Polizzi et al.)

  • Substitute numerical integration formula for Cauchy resolvent integral
  • Then subspace iterate the so obtained operator
  • Use an appropriate subspace extraction (eg. RR) to determine

eigenvalue/vector approximation

  • Other contour integration approaches: Beyn, Sakurai & Sugiura, RIM
  • f Sun et al.

In this talk:

Luka Grubiˇ si´ c 3 / 36

slide-4
SLIDE 4

About this talk

FEAST idea (Polizzi et al.)

  • Substitute numerical integration formula for Cauchy resolvent integral
  • Then subspace iterate the so obtained operator
  • Use an appropriate subspace extraction (eg. RR) to determine

eigenvalue/vector approximation

  • Other contour integration approaches: Beyn, Sakurai & Sugiura, RIM
  • f Sun et al.

In this talk:

  • Looping contours is naturally dictated by perturbation theory

Luka Grubiˇ si´ c 3 / 36

slide-5
SLIDE 5

About this talk

FEAST idea (Polizzi et al.)

  • Substitute numerical integration formula for Cauchy resolvent integral
  • Then subspace iterate the so obtained operator
  • Use an appropriate subspace extraction (eg. RR) to determine

eigenvalue/vector approximation

  • Other contour integration approaches: Beyn, Sakurai & Sugiura, RIM
  • f Sun et al.

In this talk:

  • Looping contours is naturally dictated by perturbation theory
  • It will work – for self adjoint operators – with coarse integration

formulae

Luka Grubiˇ si´ c 3 / 36

slide-6
SLIDE 6

About this talk

FEAST idea (Polizzi et al.)

  • Substitute numerical integration formula for Cauchy resolvent integral
  • Then subspace iterate the so obtained operator
  • Use an appropriate subspace extraction (eg. RR) to determine

eigenvalue/vector approximation

  • Other contour integration approaches: Beyn, Sakurai & Sugiura, RIM
  • f Sun et al.

In this talk:

  • Looping contours is naturally dictated by perturbation theory
  • It will work – for self adjoint operators – with coarse integration

formulae

  • With very few iterations per refinement level (asymptotically

speaking)

Luka Grubiˇ si´ c 3 / 36

slide-7
SLIDE 7

About this talk

FEAST idea (Polizzi et al.)

  • Substitute numerical integration formula for Cauchy resolvent integral
  • Then subspace iterate the so obtained operator
  • Use an appropriate subspace extraction (eg. RR) to determine

eigenvalue/vector approximation

  • Other contour integration approaches: Beyn, Sakurai & Sugiura, RIM
  • f Sun et al.

In this talk:

  • Looping contours is naturally dictated by perturbation theory
  • It will work – for self adjoint operators – with coarse integration

formulae

  • With very few iterations per refinement level (asymptotically

speaking)

  • Strong resolvent convergence with“contraction property”is enough

Luka Grubiˇ si´ c 3 / 36

slide-8
SLIDE 8

Model operator

Model operator Au = −△u − Vu, u ∈ H1

∗(Ω).

  • Ω ⊂ Rd open and boundˇ

zed

  • A self-adjoint and unbounded
  • (ξ − A)−1 operator valued function
  • Spec(A) is countable without finite accumulation points.
  • Notation: Aψi = λiψi, and −V ≤ λ1 ≤ λ2 ≤ · · ·
  • We are counting eigenvalues with multiplicity.
  • Variational (energy) scalar product (u, v)V = a[u, v], u, v ∈ H1

∗(Ω).

About the problem Luka Grubiˇ si´ c 4 / 36

slide-9
SLIDE 9

Model operator

Model operator Au = −△u − Vu, u ∈ H1

∗(Ω).

  • Ω= Rd
  • A self-adjoint and unbounded
  • (ξ − A)−1 operator valued function
  • It holds for the discrete spectrum Specdisc(A) ⊂ [inf V , 0
  • Notation: Aψi = λiψi, and −V ≤ λ1 ≤ λ2 ≤ · · · < 0 are only

isolated eigenvalues call them Λ

  • We are counting eigenvalues with multiplicity.
  • Variational (energy) scalar product (u, v)V = a[u, v], u, v ∈ H1

∗(Ω).

  • CT estimates: χBx(ξ − A)−1χBy ≤ C exp(−τx − y)

About the problem Luka Grubiˇ si´ c 5 / 36

slide-10
SLIDE 10

Bounded and extended states

Quantum states from chebfun V (t)= −50

  • exp(−t2)+0.7 exp(−(t − 3)2)
  • About the problem

Luka Grubiˇ si´ c 6 / 36

slide-11
SLIDE 11

Bounded and extended states

Quantum states from chebfun V (t)= −50

  • exp(−t2)+0.7 exp(−(t − 3)2)
  • We consider AR on finite ΩR.

About the problem Luka Grubiˇ si´ c 6 / 36

slide-12
SLIDE 12

Bounded and extended states

Quantum states from chebfun V (t)= −50

  • exp(−t2)+0.7 exp(−(t − 3)2)
  • We consider AR on finite ΩR.
  • Perturbation argument: put a contour on [−V , 0, unwanted

clustered eigenvalues in [0, ∞.

About the problem Luka Grubiˇ si´ c 6 / 36

slide-13
SLIDE 13

Also in 2D (thanks to NG Solve)

About the problem Luka Grubiˇ si´ c 7 / 36

slide-14
SLIDE 14

However!

  • We would like to spice things up
  • essential spectrum

About the problem Luka Grubiˇ si´ c 8 / 36

slide-15
SLIDE 15

However!

  • We would like to spice things up
  • essential spectrum

extreme clustering

  • Enter

About the problem Luka Grubiˇ si´ c 8 / 36

slide-16
SLIDE 16

However!

  • We would like to spice things up
  • essential spectrum

extreme clustering

  • Enter

1 Periodic potentials

About the problem Luka Grubiˇ si´ c 8 / 36

slide-17
SLIDE 17

However!

  • We would like to spice things up
  • essential spectrum

extreme clustering

  • Enter

1 Periodic potentials 2 First order block operators, e.g. Dirac operators

About the problem Luka Grubiˇ si´ c 8 / 36

slide-18
SLIDE 18

However!

  • We would like to spice things up
  • essential spectrum

extreme clustering

  • Enter

1 Periodic potentials 2 First order block operators, e.g. Dirac operators 3 Strong resolvent convergence.

About the problem Luka Grubiˇ si´ c 8 / 36

slide-19
SLIDE 19

Model operator – periodic potential e.g. V = cos(·)

Model operator Au = −△u − Vu, u ∈ H1

∗(Ω).

  • Ω= Rd
  • A self-adjoint and unbounded
  • (ξ − A)−1 operator valued function
  • Only essential spectrum with band gaps.
  • Notation: Aψi = λiψi, for isolated eigenvalues

call them Λ

  • We are counting eigenvalues with multiplicity.
  • Variational (energy) scalar product (u, v)V = a[u, v], u, v ∈ H1

∗(Ω).

About the problem Luka Grubiˇ si´ c 9 / 36

slide-20
SLIDE 20

Quantum pendulum

About the problem Luka Grubiˇ si´ c 10 / 36

slide-21
SLIDE 21

Quantum pendulum – has band gaps

About the problem Luka Grubiˇ si´ c 11 / 36

slide-22
SLIDE 22

Model operator – V = cos(x) − exp(−x2)

Model operator Au = −△u − Vu, u ∈ H1

∗(Ω).

  • Ω= Rd
  • A self-adjoint and unbounded
  • (ξ − A)−1 operator valued function
  • The same essential spectrum with band gaps and isolated eigenvalues

in the middle.

  • Notation: Aψi = λiψi, for isolated eigenvalues

call them Λ

  • We are counting eigenvalues with multiplicity.
  • Variational (energy) scalar product (u, v)V = a[u, v], u, v ∈ H1

∗(Ω).

About the problem Luka Grubiˇ si´ c 12 / 36

slide-23
SLIDE 23

Quantum pendulum – with a bump

About the problem Luka Grubiˇ si´ c 13 / 36

slide-24
SLIDE 24

Quantum pendulum – with a bump

About the problem Luka Grubiˇ si´ c 13 / 36

slide-25
SLIDE 25

Cascading spectrum

The bump is a relatively compact perturbation

About the problem Luka Grubiˇ si´ c 14 / 36

slide-26
SLIDE 26

Cascading spectrum

The bump is a relatively compact perturbation

About the problem Luka Grubiˇ si´ c 14 / 36

slide-27
SLIDE 27

Cascading spectrum

The bump is a relatively compact perturbation

About the problem Luka Grubiˇ si´ c 14 / 36

slide-28
SLIDE 28

Cascading spectrum

The bump is a relatively compact perturbation

About the problem Luka Grubiˇ si´ c 14 / 36

slide-29
SLIDE 29

Cascading spectrum

The bump is a relatively compact perturbation

About the problem Luka Grubiˇ si´ c 14 / 36

slide-30
SLIDE 30

Cascading spectrum

The bump is a relatively compact perturbation

About the problem Luka Grubiˇ si´ c 14 / 36

slide-31
SLIDE 31

Cascading spectrum

The bump is a relatively compact perturbation

About the problem Luka Grubiˇ si´ c 14 / 36

slide-32
SLIDE 32

Cascading spectrum

The bump is a relatively compact perturbation

About the problem Luka Grubiˇ si´ c 14 / 36

slide-33
SLIDE 33

First order operators are good candidates

Dirac Op with radial symmetry and the Culomb potential

  • Defined by the block operator matrix

A = mc2 − Z

r

−∂r + κ

r

∂r + κ

r

−mc2 − Z

r

  • with Dom(A) = H1(R) ⊕ H1(R).
  • Z = 1, 2, ..., 137 is the electric charge number
  • κ is the spin orbit coupling parameter
  • See L. Boulton’s web page

About the problem Luka Grubiˇ si´ c 15 / 36

slide-34
SLIDE 34

Associated norms – not used in the standard analysis so far

  • L2 norm ·

About the problem Luka Grubiˇ si´ c 16 / 36

slide-35
SLIDE 35

Associated norms – not used in the standard analysis so far

  • L2 norm ·
  • Energy norm A1/2 ·

About the problem Luka Grubiˇ si´ c 16 / 36

slide-36
SLIDE 36

Associated norms – not used in the standard analysis so far

  • L2 norm ·
  • Energy norm A1/2 ·
  • Graph norm
  • · 2 + A · 2.

About the problem Luka Grubiˇ si´ c 16 / 36

slide-37
SLIDE 37

The task

1 Start from sa self adjoint operator A and its spectrum Σ(A)

About the method Luka Grubiˇ si´ c 17 / 36

slide-38
SLIDE 38

The task

1 Start from sa self adjoint operator A and its spectrum Σ(A) 2 Target Λ ⊂ Σ(A) consisting of m isolated points, where we count the

eigenvalues according to multiplicity.

About the method Luka Grubiˇ si´ c 17 / 36

slide-39
SLIDE 39

The task

1 Start from sa self adjoint operator A and its spectrum Σ(A) 2 Target Λ ⊂ Σ(A) consisting of m isolated points, where we count the

eigenvalues according to multiplicity.

3 Also, let E be the (m-dimensional) span of the corresponding

eigenvectors,

About the method Luka Grubiˇ si´ c 17 / 36

slide-40
SLIDE 40

The task

1 Start from sa self adjoint operator A and its spectrum Σ(A) 2 Target Λ ⊂ Σ(A) consisting of m isolated points, where we count the

eigenvalues according to multiplicity.

3 Also, let E be the (m-dimensional) span of the corresponding

eigenvectors,

4 And so the heroes are:

Λ = {λ1, . . . , λm}, E = span{e1, . . . , em}, Aei = λiei.

About the method Luka Grubiˇ si´ c 17 / 36

slide-41
SLIDE 41

The tools

  • Let Γ be a simple closed contour enclosing Λ
  • Define

r(z) = 1 2π 1i

  • Γ

(ξ − z)−1 =

  • 1

, z ∈ G , z ∈ C \ (Λ ∪ Γ)

  • Spectral projection of A for Γ is

S = r(A) = 1 2π 1i

  • Γ

(ξ − A)−1 dξ .

  • If Au = λu then r(A)u = r(λ)u.

About the method Luka Grubiˇ si´ c 18 / 36

slide-42
SLIDE 42

The tools

  • Let Γ be a simple closed contour enclosing Λ
  • Define

r(z) = 1 2π 1i

  • Γ

(ξ − z)−1 =

  • 1

, z ∈ G , z ∈ C \ (Λ ∪ Γ)

  • Spectral projection of A for Γ is

S = r(A) = 1 2π 1i

  • Γ

(ξ − A)−1 dξ .

  • If Au = λu then r(A)u = r(λ)u.
  • rN(·) is“ripping”Λ from Σ(A)!

About the method Luka Grubiˇ si´ c 18 / 36

slide-43
SLIDE 43

Numerical analyst’s take

  • Rational approximation of r:

r(ξ) = 1 2π 1i

  • Γ

(ξ − z)−1 dξ ≈ rN(ξ) =

N−1

  • i=1

ωi(zi − ξ)−1 .

About the method Luka Grubiˇ si´ c 19 / 36

slide-44
SLIDE 44

Numerical analyst’s take

  • Rational approximation of r:

r(ξ) = 1 2π 1i

  • Γ

(ξ − z)−1 dξ ≈ rN(ξ) =

N−1

  • i=1

ωi(zi − ξ)−1 .

  • If Au = λu then rN(A)u = rN(λ)u.

About the method Luka Grubiˇ si´ c 19 / 36

slide-45
SLIDE 45

Numerical analyst’s take

  • Rational approximation of r:

r(ξ) = 1 2π 1i

  • Γ

(ξ − z)−1 dξ ≈ rN(ξ) =

N−1

  • i=1

ωi(zi − ξ)−1 .

  • If Au = λu then rN(A)u = rN(λ)u.
  • We use (SN :=)

rN(A) = N−1

i=1 ωi(zi − A)−1

About the method Luka Grubiˇ si´ c 19 / 36

slide-46
SLIDE 46

Numerical analyst’s take

  • Rational approximation of r:

r(ξ) = 1 2π 1i

  • Γ

(ξ − z)−1 dξ ≈ rN(ξ) =

N−1

  • i=1

ωi(zi − ξ)−1 .

  • If Au = λu then rN(A)u = rN(λ)u.
  • We use (SN :=)

rN(A) = N−1

i=1 ωi(zi − A)−1

  • rN(·) is filtering Λ from Σ(A)!

About the method Luka Grubiˇ si´ c 19 / 36

slide-47
SLIDE 47

Note!

The spectrum of rN(A) looks like Ideal for subspace iteration towards the dominant spectral component

About the method Luka Grubiˇ si´ c 20 / 36

slide-48
SLIDE 48

Note!

The spectrum of rN(A) looks like Ideal for subspace iteration towards the dominant spectral component We know which singular values of rN(A) arebig and which are small

About the method Luka Grubiˇ si´ c 20 / 36

slide-49
SLIDE 49

Spectral separation SVD picture — a semibounded example

About the method Luka Grubiˇ si´ c 21 / 36

slide-50
SLIDE 50

Spectral separation SVD picture — a semibounded example

About the method Luka Grubiˇ si´ c 21 / 36

slide-51
SLIDE 51

Randomized rank estimator based on contour integration

  • Starting from Halko, Martinsson and Tropp Finding structure with

randomness, SIREW.

  • the twist is the use of
  • Π as a projection operator onto a“finite element”space based on the

regularity theory

  • approximate range SVD of ΠP will yield the rank estimator
  • Trapeze rule converges exponentially in the number of nodes N for

periodic functions and yields rN(A) = PN.

About the method Luka Grubiˇ si´ c 22 / 36

slide-52
SLIDE 52

Randomized rank estimator — the statements

  • We prove – let k be the rank of P and p an oversampling parameter

then if we generate a random subspace S ⊂ (Ran(Π)) of dimension k + p then P − SP = (I − S)P = O(interpolation) + O(

  • k + p)σk+1(ΠP)

with probability (with respect to Gaussian measure in Ran(Π))

  • f at least 1 − 6p−p.
  • svd of SΠP yields almost optimal rank k approximation of P.
  • even a small oversampling with say p = 5 already yields highly

satisfactory estimates, e.g. statements with probability of 1 − 6 · 5−5 ≈ 99.8%

  • Compressive sampling.

About the method Luka Grubiˇ si´ c 23 / 36

slide-53
SLIDE 53

FEAST = filtered subspace iteration + RR extraction

  • Let W be a (complex) Hilbert space and let B : W → W be a

bounded linear operator.

  • Let Υ be a finite set of eigenvalues of finite multiplicity of B that are

isolated from the rest of Σ(B).

  • Consider ˜

Bℓ = B + Dℓ where Dℓ : W → W is a bounded linear

  • perator representing perturbations at step l of the iteration.
  • The iterations are started using a given initial finite-dimensional

subspace Q0 ⊂ W.

  • At step ℓ, the inexact subspace iteration computes the subspace

Qℓ = ˜ BℓQℓ−1, ℓ = 1, 2, . . . . (1)

About the method Luka Grubiˇ si´ c 24 / 36

slide-54
SLIDE 54

Assumptions

1 Consider the subspace iteration with perturbed filter ˜

B(ℓ) set to ˜ S(ℓ)

N = SN + Dℓ

assuming Dℓ → 0 strongly.

2 The perturbation size is measured by

˜ S(ℓ)

N = SN + Dℓ,

Dlq ≤ ηlq, q ∈ Ql (2)

About the method Luka Grubiˇ si´ c 25 / 36

slide-55
SLIDE 55

Assumptions

1 Consider the subspace iteration with perturbed filter ˜

B(ℓ) set to ˜ S(ℓ)

N = SN + Dℓ

assuming Dℓ → 0 strongly.

2 The perturbation size is measured by

˜ S(ℓ)

N = SN + Dℓ,

Dlq ≤ ηlq, q ∈ Ql (2)

3 If encircled eigenvalues have H1+α(Ω) regularity, then

ηℓ = O(hα) where h is the discretization parameter.

About the method Luka Grubiˇ si´ c 25 / 36

slide-56
SLIDE 56

Assumptions

1 Consider the subspace iteration with perturbed filter ˜

B(ℓ) set to ˜ S(ℓ)

N = SN + Dℓ

assuming Dℓ → 0 strongly.

2 The perturbation size is measured by

˜ S(ℓ)

N = SN + Dℓ,

Dlq ≤ ηlq, q ∈ Ql (2)

3 If encircled eigenvalues have H1+α(Ω) regularity, then

ηℓ = O(hα) where h is the discretization parameter.

4 Possible to include further computational errors in ηℓ as well!

About the method Luka Grubiˇ si´ c 25 / 36

slide-57
SLIDE 57

Assumptions

1 Consider the subspace iteration with perturbed filter ˜

B(ℓ) set to ˜ S(ℓ)

N = SN + Dℓ

assuming Dℓ → 0 strongly.

2 The perturbation size is measured by

˜ S(ℓ)

N = SN + Dℓ,

Dlq ≤ ηlq, q ∈ Ql (2)

3 If encircled eigenvalues have H1+α(Ω) regularity, then

ηℓ = O(hα) where h is the discretization parameter.

4 Possible to include further computational errors in ηℓ as well! 5 The ability to shrink ηl adaptively is key to convergence

About the method Luka Grubiˇ si´ c 25 / 36

slide-58
SLIDE 58

Assumptions

1 Consider the subspace iteration with perturbed filter ˜

B(ℓ) set to ˜ S(ℓ)

N = SN + Dℓ

assuming Dℓ → 0 strongly.

2 The perturbation size is measured by

˜ S(ℓ)

N = SN + Dℓ,

Dlq ≤ ηlq, q ∈ Ql (2)

3 If encircled eigenvalues have H1+α(Ω) regularity, then

ηℓ = O(hα) where h is the discretization parameter.

4 Possible to include further computational errors in ηℓ as well! 5 The ability to shrink ηl adaptively is key to convergence 6 Finite number of contour integration points allows to get away with

strong resolvent convergence and a sufficiently rich initial space

About the method Luka Grubiˇ si´ c 25 / 36

slide-59
SLIDE 59

Measuring the subspace error

Following H. Elman (2011–)and Y. Saad (2016).

  • For subspaces E = SW and Q of the same dimension set

t(E, Q) = max

x∈Q

(I − S)x Sx .

  • The measure t(E, Q) is obviously related to a tangent of the angle

between subspaces, but this interpretation is not necessary for the result that follows.

  • By a direct computation

gap(E, Q) = max

q∈Q

(I − S)q q ≤ max

q∈Q

(I − S)q Sq = t(E, Q) .

About the method Luka Grubiˇ si´ c 26 / 36

slide-60
SLIDE 60

Inexact subspace iteration – use a priori bounds on rN!

Theorem

Assume B, P and Q0 satisfy the assumptions. Set ˜ Bl = B + Dl, Ql = ˜ BlQ0 and suppose Dl, l ∈ N are bounded operators such that Dlq ≤ ηlq, q ∈ Ql. Then with the notation tl = t(E, Ql) we have tl ≤ µ2tl−1 + ηl

  • 1 + t2

l−1

µ1 − ηl

  • 1 + t2

l−1

.

About the method Luka Grubiˇ si´ c 27 / 36

slide-61
SLIDE 61

Inexact subspace iteration – use a priori bounds on rN!

Theorem

.... tl ≤ µ2tl−1 + ηl

  • 1 + t2

l−1

µ1 − ηl

  • 1 + t2

l−1

.

Proof:

max

q∈Q0

(I − P)˜ Bq P ˜ Bq ≤ µ2 µ1 |1 − δ| max

q∈Q0

(I − P)q Pq + 1 µ1 |1 − δ| max

q∈Q0

(I − P)Dq Pq .

About the method Luka Grubiˇ si´ c 27 / 36

slide-62
SLIDE 62

Finally

Theorem

Assuming perturbations are diminishing“quickly enough”and the spectral separation parameters are µ2 = rad(B(I − P)) and µ1 = inf{BPx : x ∈ Ran(P)} such that ˜ κ = µ2/(µ1 + γ). Then t(E, Ql) ≤ l˜ κl .

Theorem

Let Dℓ be converging strongly to zero, then there is a subsequence Dℓj such that so defined subspaces Qj satisfy t(E, Qj) ≤ j˜ κj .

About the method Luka Grubiˇ si´ c 28 / 36

slide-63
SLIDE 63

Rational approx. + restriction = computable estimate

  • Rational approximation of r:

r(ξ) = 1 2π 1i

  • Γ

(ξ − z)−1 dξ ≈ rN(ξ) =

N−1

  • i=1

ωi(zi − ξ)−1 .

  • Take Rh(z) : L2(Ω) → Vh, and define

SN,h =

N−1

  • i=1

ωiRh(zi) .

Discrete resolvent estimates Luka Grubiˇ si´ c 29 / 36

slide-64
SLIDE 64

Rational approx. + restriction = computable estimate

  • Rational approximation of r:

r(ξ) = 1 2π 1i

  • Γ

(ξ − z)−1 dξ ≈ rN(ξ) =

N−1

  • i=1

ωi(zi − ξ)−1 .

  • Take Rh(z) : L2(Ω) → Vh, and define

SN,h =

N−1

  • i=1

ωiRh(zi) .

  • Rh(.) is a fin. dim. realization of the resolvent.

Discrete resolvent estimates Luka Grubiˇ si´ c 29 / 36

slide-65
SLIDE 65

Discrete resolvent calculus

The ratio κ = sup

µ∈Σ(A)\Λ

|rN(µ)|/ min

λ∈Λ |rN(λ)|

(3) is a measure of filter quality (for separating spectral components). Let y, γ, be as before, then Note W = N−1

k=0 |wk|

Theorem

Let zk ∈ Γ, k = 1, ..., N for Γ ⊂ C \ Σ(A), then SN − Sh

NV ≤ W

max

k=1,...,N Rh(zk) − R(zk)V

gapV(E, Eh) ≤ CNW max

k=1,...,N (Rh(zk) − R(zk))

  • E
  • V .

Constant CN depends on the filter function.

Discrete resolvent estimates Luka Grubiˇ si´ c 30 / 36

slide-66
SLIDE 66

All together

the estimates

Discrete resolvent estimates Luka Grubiˇ si´ c 31 / 36

slide-67
SLIDE 67

All together

the estimates

  • for eigenvectors set V = W(= L2), then

gapW(E, Eh) ≤ CNW max

k=1,...,N

  • R(zk) − Rh(zk)
  • E
  • W ,

(4)

Discrete resolvent estimates Luka Grubiˇ si´ c 31 / 36

slide-68
SLIDE 68

All together

the estimates

  • for eigenvectors set V = W(= L2), then

gapW(E, Eh) ≤ CNW max

k=1,...,N

  • R(zk) − Rh(zk)
  • E
  • W ,

(4)

  • and eigenvalues (using |A|1/2 · as energy)

dist(Λ, Λh) ≤ Ca gapV(E, Eh)2 + CbgapH(E, Eh)2. follow

Discrete resolvent estimates Luka Grubiˇ si´ c 31 / 36

slide-69
SLIDE 69

The algorithm

Pick X on an a priori generated mesh for a coarse projection

1 Compute Y = SNX 2 compute Ritz values eig(Y ∗AY ) and the residual

R = AY − Y (Y ∗AY )

3 accept Ran(Y ) basis vectors with a small residual and prune the rest

and refine the mesh for the accepted R-values

4 form new Y from the refined projection and go to ONE until

convergence

Numerical experiments Luka Grubiˇ si´ c 32 / 36

slide-70
SLIDE 70

Once is enough

Study the convergence history of a refinement sequence where on each level we apply the filtered operator just once!

1000 2000 3000 4000 5000

  • 8
  • 6
  • 4
  • 2

#dof log10(error) Benchmark h

2

Practical λ1 Eigs λ1 Practical λ2 Eigs λ2

Numerical experiments Luka Grubiˇ si´ c 33 / 36

slide-71
SLIDE 71

Graph norm convergence of eigenspaces

a Matlab implementation of FEAST + Chebyshev polynomial discretization – the growth of the algebraic error

Numerical experiments Luka Grubiˇ si´ c 34 / 36

slide-72
SLIDE 72

When to use numerical rational functions calculuss

  • eigenvalues of infinite multiplicity

Numerical experiments Luka Grubiˇ si´ c 35 / 36

slide-73
SLIDE 73

When to use numerical rational functions calculuss

  • eigenvalues of infinite multiplicity
  • indefinite operators (Maxwell,Dirac,...)

Numerical experiments Luka Grubiˇ si´ c 35 / 36

slide-74
SLIDE 74

When to use numerical rational functions calculuss

  • eigenvalues of infinite multiplicity
  • indefinite operators (Maxwell,Dirac,...)
  • nonlocal operators

Numerical experiments Luka Grubiˇ si´ c 35 / 36

slide-75
SLIDE 75

When to use numerical rational functions calculuss

  • eigenvalues of infinite multiplicity
  • indefinite operators (Maxwell,Dirac,...)
  • nonlocal operators

◮ fractional differential operators ◮ semi-groups ◮ etc Numerical experiments Luka Grubiˇ si´ c 35 / 36

slide-76
SLIDE 76

When to use numerical rational functions calculuss

  • eigenvalues of infinite multiplicity
  • indefinite operators (Maxwell,Dirac,...)
  • nonlocal operators

◮ fractional differential operators ◮ semi-groups ◮ etc

  • Work in progress: Maxwell

M =

  • 1i 1

εcurl

− 1i 1

µcurl

  • with Dom(M) = H(curl) ⊕ H(curl).

Numerical experiments Luka Grubiˇ si´ c 35 / 36

slide-77
SLIDE 77

When to use numerical rational functions calculuss

  • eigenvalues of infinite multiplicity
  • indefinite operators (Maxwell,Dirac,...)
  • nonlocal operators

◮ fractional differential operators ◮ semi-groups ◮ etc

  • Work in progress: Maxwell

M =

  • 1i 1

εcurl

− 1i 1

µcurl

  • with Dom(M) = H(curl) ⊕ H(curl).
  • Robust end efficient resolvent approximations, such as those by

discontinuous Petrov – Galerkin methods.

Numerical experiments Luka Grubiˇ si´ c 35 / 36

slide-78
SLIDE 78

Acknowledgment

Research supported in part by the Croatian Science Foundation grant nr. HRZZ 9345

Acknowledgment

University of Zagreb, Faculty of Sci- ence research grant ” Numerical and analytical methods for parameter de- pendent eigenvalue problems” and a bilateral NSF-MSO grant. Thank you for your attention!

Numerical experiments Luka Grubiˇ si´ c 36 / 36