Core-Chasing Algorithms for the Eigenvalue Problem David S. Watkins - - PowerPoint PPT Presentation

▶
core chasing algorithms for the eigenvalue problem
SMART_READER_LITE
LIVE PREVIEW

Core-Chasing Algorithms for the Eigenvalue Problem David S. Watkins - - PowerPoint PPT Presentation

Apr 07, 2023 864 likes •2.55k views

Core-Chasing Algorithms for the Eigenvalue Problem David S. Watkins Department of Mathematics Washington State University July, 2016 David S. Watkins Core-Chasing Algorithms Our International Research Group This is joint work with (Oxford

slide-1
SLIDE 1

Core-Chasing Algorithms for the Eigenvalue Problem

David S. Watkins

Department of Mathematics Washington State University

July, 2016

David S. Watkins Core-Chasing Algorithms

slide-2
SLIDE 2

Our International Research Group

This is joint work with Jared Aurentz (Oxford → Madrid) Thomas Mach (KU Leuven → Astana) Raf Vandebril (KU Leuven)

David S. Watkins Core-Chasing Algorithms

slide-3
SLIDE 3

Today’s Topic

David S. Watkins Core-Chasing Algorithms

slide-4
SLIDE 4

Today’s Topic

The matrix eigenvalue problem

David S. Watkins Core-Chasing Algorithms

slide-5
SLIDE 5

Today’s Topic

The matrix eigenvalue problem A ∈ Cn×n

David S. Watkins Core-Chasing Algorithms

slide-6
SLIDE 6

Today’s Topic

The matrix eigenvalue problem A ∈ Cn×n Find the eigenvalues (. . . vectors, invariant subspaces)

David S. Watkins Core-Chasing Algorithms

slide-7
SLIDE 7

Today’s Topic

The matrix eigenvalue problem A ∈ Cn×n Find the eigenvalues (. . . vectors, invariant subspaces) Possible structures

unitary unitary-plus-rank-one (companion) symmetric symmetric-plus-rank-one (colleague)

David S. Watkins Core-Chasing Algorithms

slide-8
SLIDE 8

Today’s Topic

The matrix eigenvalue problem A ∈ Cn×n Find the eigenvalues (. . . vectors, invariant subspaces) Possible structures

unitary unitary-plus-rank-one (companion) symmetric symmetric-plus-rank-one (colleague) . . . or no special structure

David S. Watkins Core-Chasing Algorithms

slide-9
SLIDE 9

John Francis

photo: Frank Uhlig, 2009

David S. Watkins Core-Chasing Algorithms

slide-10
SLIDE 10

John Francis

photo: Frank Uhlig, 2009 invented the winning algorithm in 1959.

David S. Watkins Core-Chasing Algorithms

slide-11
SLIDE 11

John Francis

photo: Frank Uhlig, 2009 invented the winning algorithm in 1959. implicitly shifted QR algorithm

David S. Watkins Core-Chasing Algorithms

slide-12
SLIDE 12

John Francis

photo: Frank Uhlig, 2009 invented the winning algorithm in 1959. implicitly shifted QR algorithm Our algorithms are all variants of this.

David S. Watkins Core-Chasing Algorithms

slide-13
SLIDE 13

Francis’s algorithm . . .

David S. Watkins Core-Chasing Algorithms

slide-14
SLIDE 14

Francis’s algorithm . . . . . . is a bulge chasing algorithm.

David S. Watkins Core-Chasing Algorithms

slide-15
SLIDE 15

Francis’s algorithm . . . . . . is a bulge chasing algorithm. We turn it into a core chasing algorithm.

David S. Watkins Core-Chasing Algorithms

slide-16
SLIDE 16

Francis’s algorithm . . . . . . is a bulge chasing algorithm. We turn it into a core chasing algorithm. Instead of chasing bulges, we chase core transformations.

David S. Watkins Core-Chasing Algorithms

slide-17
SLIDE 17

Core Transformations

What is a core transformation?

David S. Watkins Core-Chasing Algorithms

slide-18
SLIDE 18

Core Transformations

What is a core transformation? It’s a unitary matrix, and

David S. Watkins Core-Chasing Algorithms

slide-19
SLIDE 19

Core Transformations

What is a core transformation? It’s a unitary matrix, and it’s essentially 2 × 2 C2 =       1 × × × × 1 1      

David S. Watkins Core-Chasing Algorithms

slide-20
SLIDE 20

Core Transformations

What is a core transformation? It’s a unitary matrix, and it’s essentially 2 × 2 C2 =       1 × × × × 1 1       Ex: Givens rotator, reflector, . . .

David S. Watkins Core-Chasing Algorithms

slide-21
SLIDE 21

Core Transformations

What is a core transformation? It’s a unitary matrix, and it’s essentially 2 × 2 C2 =       1 × × × × 1 1       Ex: Givens rotator, reflector, . . . We just wanted a generic term.

David S. Watkins Core-Chasing Algorithms

slide-22
SLIDE 22

Core Transformations

  1 × × × ×     × × × × × × ×   =   × × × × × ×  

David S. Watkins Core-Chasing Algorithms

slide-23
SLIDE 23

Core Transformations

  1 × × × ×     × × × × × × ×   =   × × × × × ×   Abbreviated notation

David S. Watkins Core-Chasing Algorithms

slide-24
SLIDE 24

Core Transformations

  1 × × × ×     × × × × × × ×   =   × × × × × ×   Abbreviated notation

  • ×

=

David S. Watkins Core-Chasing Algorithms

slide-25
SLIDE 25

Hessenberg QR decomposition

× × × × × = × × × × ×

David S. Watkins Core-Chasing Algorithms

slide-26
SLIDE 26

Hessenberg QR decomposition

  • ×

× × × × = × × × ×

David S. Watkins Core-Chasing Algorithms

slide-27
SLIDE 27

Hessenberg QR decomposition

  • ×

× × × × = × × ×

David S. Watkins Core-Chasing Algorithms

slide-28
SLIDE 28

Hessenberg QR decomposition

  • ×

× × × × = × ×

David S. Watkins Core-Chasing Algorithms

slide-29
SLIDE 29

Hessenberg QR decomposition

  • ×

× × × × = ×

David S. Watkins Core-Chasing Algorithms

slide-30
SLIDE 30

Hessenberg QR decomposition

  • ×

× × × × =

David S. Watkins Core-Chasing Algorithms

slide-31
SLIDE 31

Hessenberg QR decomposition

  • ×

× × × × = Now invert the core transformations.

David S. Watkins Core-Chasing Algorithms

slide-32
SLIDE 32

Hessenberg QR decomposition

× × × × × =

  • Q

R

David S. Watkins Core-Chasing Algorithms

slide-33
SLIDE 33

Our algorithms operate on the matrix in QR decomposed form. A = QR =

  • This is not inefficient.

We apply Francis’s algorithm to this factored form.

David S. Watkins Core-Chasing Algorithms

slide-34
SLIDE 34

Operating on Core Transformations

Fusion

  • ⇒
  • David S. Watkins

Core-Chasing Algorithms

slide-35
SLIDE 35

Operating on Core Transformations

Turnover

  • ⇔

∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗

  • ⇔
  • David S. Watkins

Core-Chasing Algorithms

slide-36
SLIDE 36

Operating on Core Transformations

Turnover as a shift-through operation

  • David S. Watkins

Core-Chasing Algorithms

slide-37
SLIDE 37

Operating on Core Transformations

Turnover as a shift-through operation

  • David S. Watkins

Core-Chasing Algorithms

slide-38
SLIDE 38

Operating on Core Transformations

Turnover as a shift-through operation

  • David S. Watkins

Core-Chasing Algorithms

slide-39
SLIDE 39

Operating on Core Transformations

Turnover as a shift-through operation

  • David S. Watkins

Core-Chasing Algorithms

slide-40
SLIDE 40

Operating on Core Transformations

Turnover as a shift-through operation

  • David S. Watkins

Core-Chasing Algorithms

slide-41
SLIDE 41

Operating on Core Transformations

Turnover as a shift-through operation

  • David S. Watkins

Core-Chasing Algorithms

slide-42
SLIDE 42

Operating on Core Transformations

Turnover as a shift-through operation

  • David S. Watkins

Core-Chasing Algorithms

slide-43
SLIDE 43

Operating on Core Transformations

Turnover as a shift-through operation : Abbreviated notation

  • David S. Watkins

Core-Chasing Algorithms

slide-44
SLIDE 44

Operating on Core Transformations

Turnover as a shift-through operation : Abbreviated notation

  • David S. Watkins

Core-Chasing Algorithms

slide-45
SLIDE 45

Operating on Core Transformations

Passing a core transformation through a triangular matrix

   ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗    ⇔    ∗ ∗ ∗ ∗ ∗ ∗ ∗ + ∗ ∗ ∗    ⇔

  ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗   

Cost is O(n) flops.

David S. Watkins Core-Chasing Algorithms

slide-46
SLIDE 46

Operating on Core Transformations

Passing a core transformation through a triangular matrix Abbreviated Notation:

  • David S. Watkins

Core-Chasing Algorithms

slide-47
SLIDE 47

Core Chasing

  • David S. Watkins

Core-Chasing Algorithms

slide-48
SLIDE 48

Core Chasing

  • David S. Watkins

Core-Chasing Algorithms

slide-49
SLIDE 49

Core Chasing

  • David S. Watkins

Core-Chasing Algorithms

slide-50
SLIDE 50

Core Chasing

  • David S. Watkins

Core-Chasing Algorithms

slide-51
SLIDE 51

Core Chasing

  • David S. Watkins

Core-Chasing Algorithms

slide-52
SLIDE 52

Core Chasing

  • David S. Watkins

Core-Chasing Algorithms

slide-53
SLIDE 53

Core Chasing

  • David S. Watkins

Core-Chasing Algorithms

slide-54
SLIDE 54

Core Chasing

  • David S. Watkins

Core-Chasing Algorithms

slide-55
SLIDE 55

Core Chasing

  • David S. Watkins

Core-Chasing Algorithms

slide-56
SLIDE 56

Core Chasing

  • David S. Watkins

Core-Chasing Algorithms

slide-57
SLIDE 57

Core Chasing

  • David S. Watkins

Core-Chasing Algorithms

slide-58
SLIDE 58

Core Chasing

  • David S. Watkins

Core-Chasing Algorithms

slide-59
SLIDE 59

Core Chasing

  • David S. Watkins

Core-Chasing Algorithms

slide-60
SLIDE 60

Core Chasing

  • David S. Watkins

Core-Chasing Algorithms

slide-61
SLIDE 61

Core Chasing

  • David S. Watkins

Core-Chasing Algorithms

slide-62
SLIDE 62

Cost

David S. Watkins Core-Chasing Algorithms

slide-63
SLIDE 63

Cost

Most arithmetic in passing-through operation

David S. Watkins Core-Chasing Algorithms

slide-64
SLIDE 64

Cost

Most arithmetic in passing-through operation O(n2) flops per iteration . . . O(n3) total flops . . . about the same as for standard Francis iteration.

David S. Watkins Core-Chasing Algorithms

slide-65
SLIDE 65

Are there any advantages?

David S. Watkins Core-Chasing Algorithms

slide-66
SLIDE 66

Are there any advantages?

unitary case

David S. Watkins Core-Chasing Algorithms

slide-67
SLIDE 67

Are there any advantages?

unitary case companion case (unitary-plus-rank-one)

David S. Watkins Core-Chasing Algorithms

slide-68
SLIDE 68

Unitary Case

David S. Watkins Core-Chasing Algorithms

slide-69
SLIDE 69

Unitary Case

A = QR =

  • David S. Watkins

Core-Chasing Algorithms

slide-70
SLIDE 70

Unitary Case

A = QR =

  • David S. Watkins

Core-Chasing Algorithms

slide-71
SLIDE 71

Unitary Case

A = QR =

  • David S. Watkins

Core-Chasing Algorithms

slide-72
SLIDE 72

Unitary Case

A = QR =

  • Cost is O(n) flops per iteration,

David S. Watkins Core-Chasing Algorithms

slide-73
SLIDE 73

Unitary Case

A = QR =

  • Cost is O(n) flops per iteration, O(n2) flops total.

David S. Watkins Core-Chasing Algorithms

slide-74
SLIDE 74

Unitary Case

A = QR =

  • Cost is O(n) flops per iteration, O(n2) flops total.

Storage requirement is O(n).

David S. Watkins Core-Chasing Algorithms

slide-75
SLIDE 75

Unitary Case

A = QR =

  • Cost is O(n) flops per iteration, O(n2) flops total.

Storage requirement is O(n). Gragg (1986)

David S. Watkins Core-Chasing Algorithms

slide-76
SLIDE 76

Unitary Case

A = QR =

  • Cost is O(n) flops per iteration, O(n2) flops total.

Storage requirement is O(n). Gragg (1986) Ammar, Reichel, M. Stewart, Bunse-Gerstner, Elsner, He, W, . . .

David S. Watkins Core-Chasing Algorithms

slide-77
SLIDE 77

Companion Case

David S. Watkins Core-Chasing Algorithms

slide-78
SLIDE 78

Companion Case

p(x) = xn + an−1xn−1 + an−2xn−2 + · · · + a0 = 0 monic polynomial

David S. Watkins Core-Chasing Algorithms

slide-79
SLIDE 79

Companion Case

p(x) = xn + an−1xn−1 + an−2xn−2 + · · · + a0 = 0 monic polynomial companion matrix A =         · · · −a0 1 · · · −a1 1 ... . . . . . . ... −an−2 1 −an−1         . . . get the zeros of p by computing the eigenvalues.

David S. Watkins Core-Chasing Algorithms

slide-80
SLIDE 80

Companion Case

p(x) = xn + an−1xn−1 + an−2xn−2 + · · · + a0 = 0 monic polynomial companion matrix A =         · · · −a0 1 · · · −a1 1 ... . . . . . . ... −an−2 1 −an−1         . . . get the zeros of p by computing the eigenvalues. MATLAB’s roots command

David S. Watkins Core-Chasing Algorithms

slide-81
SLIDE 81

Cost of solving companion eigenvalue problem

David S. Watkins Core-Chasing Algorithms

slide-82
SLIDE 82

Cost of solving companion eigenvalue problem

If structure not exploited:

O(n2) storage, O(n3) flops Francis’s algorithm

David S. Watkins Core-Chasing Algorithms

slide-83
SLIDE 83

Cost of solving companion eigenvalue problem

If structure not exploited:

O(n2) storage, O(n3) flops Francis’s algorithm

If structure exploited:

O(n) storage, O(n2) flops several methods proposed data-sparse representation + Francis’s algorithm

David S. Watkins Core-Chasing Algorithms

slide-84
SLIDE 84

Representation of R

We store the QR decomposed form.

A = QR =

  • David S. Watkins

Core-Chasing Algorithms

slide-85
SLIDE 85

Representation of R

We store the QR decomposed form.

A = QR =

  • where

R =      1 · · · −a1 1 −a2 ... . . . −a0      .

David S. Watkins Core-Chasing Algorithms

slide-86
SLIDE 86

Representation of R

We store the QR decomposed form.

A = QR =

  • where

R =      1 · · · −a1 1 −a2 ... . . . −a0      . This is unitary-plus-rank-one.

David S. Watkins Core-Chasing Algorithms

slide-87
SLIDE 87

Representation of R

We store the QR decomposed form.

A = QR =

  • where

R =      1 · · · −a1 1 −a2 ... . . . −a0      . This is unitary-plus-rank-one. How do we store it?

David S. Watkins Core-Chasing Algorithms

slide-88
SLIDE 88

Representation of R

Jared covered this yesterday.

David S. Watkins Core-Chasing Algorithms

slide-89
SLIDE 89

Representation of R

Jared covered this yesterday. Add a row and column to R.

David S. Watkins Core-Chasing Algorithms

slide-90
SLIDE 90

Representation of R

Jared covered this yesterday. Add a row and column to R. ˜ R =        1 · · · −a1 1 −a2 ... . . . . . . −a0 1 · · ·        .

David S. Watkins Core-Chasing Algorithms

slide-91
SLIDE 91

Representation of R

Jared covered this yesterday. Add a row and column to R. ˜ R =        1 · · · −a1 1 −a2 ... . . . . . . −a0 1 · · ·        . This is still unitary-plus-rank-one.

David S. Watkins Core-Chasing Algorithms

slide-92
SLIDE 92

Representation of R

˜ R =

David S. Watkins Core-Chasing Algorithms

slide-93
SLIDE 93

Representation of R

˜ R = C ∗

n · · · C ∗ 1 (B1 · · · Bn + e1yT)

David S. Watkins Core-Chasing Algorithms

slide-94
SLIDE 94

Representation of R

˜ R = C ∗

n · · · C ∗ 1 (B1 · · · Bn + e1yT)

  • +

· · ·

David S. Watkins Core-Chasing Algorithms

slide-95
SLIDE 95

Representation of R

˜ R = C ∗

n · · · C ∗ 1 (B1 · · · Bn + e1yT)

  • +

· · ·

. . . and we don’t have to store the rank-one part!

David S. Watkins Core-Chasing Algorithms

slide-96
SLIDE 96

Representation of R

˜ R = C ∗

n · · · C ∗ 1 (B1 · · · Bn + e1yT)

  • +

· · ·

. . . and we don’t have to store the rank-one part! Storage is O(n).

David S. Watkins Core-Chasing Algorithms

slide-97
SLIDE 97

Passing a core transformation through R

  • David S. Watkins

Core-Chasing Algorithms

slide-98
SLIDE 98

Passing a core transformation through R

  • C ∗

n · · · C ∗ 1

  • B1 · · · Bn
  • David S. Watkins

Core-Chasing Algorithms

slide-99
SLIDE 99

Passing a core transformation through R

  • C ∗

n · · · C ∗ 1

  • B1 · · · Bn
  • Cost: O(n) flops per iteration, O(n2) flops total.

David S. Watkins Core-Chasing Algorithms

slide-100
SLIDE 100

The rank-one part

Recovering yT from the core transformations

David S. Watkins Core-Chasing Algorithms

slide-101
SLIDE 101

The rank-one part

Recovering yT from the core transformations = eT

n+1 ˜

R

David S. Watkins Core-Chasing Algorithms

slide-102
SLIDE 102

The rank-one part

Recovering yT from the core transformations = eT

n+1 ˜

R = eT

n+1C ∗ n · · · C ∗ 1 (B1 · · · Bn + e1yT)

David S. Watkins Core-Chasing Algorithms

slide-103
SLIDE 103

The rank-one part

Recovering yT from the core transformations = eT

n+1 ˜

R = eT

n+1C ∗ n · · · C ∗ 1 (B1 · · · Bn + e1yT)

= eT

n+1C ∗(B + e1yT)

David S. Watkins Core-Chasing Algorithms

slide-104
SLIDE 104

The rank-one part

Recovering yT from the core transformations = eT

n+1 ˜

R = eT

n+1C ∗ n · · · C ∗ 1 (B1 · · · Bn + e1yT)

= eT

n+1C ∗(B + e1yT)

= eT

n+1C ∗B + (eT n+1C ∗e1)yT

David S. Watkins Core-Chasing Algorithms

slide-105
SLIDE 105

The rank-one part

Recovering yT from the core transformations = eT

n+1 ˜

R = eT

n+1C ∗ n · · · C ∗ 1 (B1 · · · Bn + e1yT)

= eT

n+1C ∗(B + e1yT)

= eT

n+1C ∗B + (eT n+1C ∗e1)yT

so yT = −ρ−1eT

n+1C ∗B,

where ρ = eT

n+1C ∗e1.

David S. Watkins Core-Chasing Algorithms

slide-106
SLIDE 106

The rank-one part

Recovering yT from the core transformations = eT

n+1 ˜

R = eT

n+1C ∗ n · · · C ∗ 1 (B1 · · · Bn + e1yT)

= eT

n+1C ∗(B + e1yT)

= eT

n+1C ∗B + (eT n+1C ∗e1)yT

so yT = −ρ−1eT

n+1C ∗B,

where ρ = eT

n+1C ∗e1.

Note: |ρ| = y −1

2

= 0.

David S. Watkins Core-Chasing Algorithms

slide-107
SLIDE 107

The rank-one part

Recovering yT from the core transformations = eT

n+1 ˜

R = eT

n+1C ∗ n · · · C ∗ 1 (B1 · · · Bn + e1yT)

= eT

n+1C ∗(B + e1yT)

= eT

n+1C ∗B + (eT n+1C ∗e1)yT

so yT = −ρ−1eT

n+1C ∗B,

where ρ = eT

n+1C ∗e1.

Note: |ρ| = y −1

2

= 0. This gives us a way to compute y at any time.

David S. Watkins Core-Chasing Algorithms

slide-108
SLIDE 108

The rank-one part

Recovering yT from the core transformations = eT

n+1 ˜

R = eT

n+1C ∗ n · · · C ∗ 1 (B1 · · · Bn + e1yT)

= eT

n+1C ∗(B + e1yT)

= eT

n+1C ∗B + (eT n+1C ∗e1)yT

so yT = −ρ−1eT

n+1C ∗B,

where ρ = eT

n+1C ∗e1.

Note: |ρ| = y −1

2

= 0. This gives us a way to compute y at any time. In fact we never need to use it, . . .

David S. Watkins Core-Chasing Algorithms

slide-109
SLIDE 109

The rank-one part

Recovering yT from the core transformations = eT

n+1 ˜

R = eT

n+1C ∗ n · · · C ∗ 1 (B1 · · · Bn + e1yT)

= eT

n+1C ∗(B + e1yT)

= eT

n+1C ∗B + (eT n+1C ∗e1)yT

so yT = −ρ−1eT

n+1C ∗B,

where ρ = eT

n+1C ∗e1.

Note: |ρ| = y −1

2

= 0. This gives us a way to compute y at any time. In fact we never need to use it, . . . . . . but it comes in handy for the backward error analysis.

David S. Watkins Core-Chasing Algorithms

slide-110
SLIDE 110

Other things we can do

We can also handle

David S. Watkins Core-Chasing Algorithms

slide-111
SLIDE 111

Other things we can do

We can also handle generalized eigenvalue problem companion pencil

David S. Watkins Core-Chasing Algorithms

slide-112
SLIDE 112

Other things we can do

We can also handle generalized eigenvalue problem companion pencil matrix polynomial eigenvalue problems (with L. Robol)

David S. Watkins Core-Chasing Algorithms

slide-113
SLIDE 113

Other things we can do

We can also handle generalized eigenvalue problem companion pencil matrix polynomial eigenvalue problems (with L. Robol) symmetric-plus-rank-one (via Cayley transform)

David S. Watkins Core-Chasing Algorithms

slide-114
SLIDE 114

Other things we can do

We can also handle generalized eigenvalue problem companion pencil matrix polynomial eigenvalue problems (with L. Robol) symmetric-plus-rank-one (via Cayley transform) completely unstructured problems

David S. Watkins Core-Chasing Algorithms

slide-115
SLIDE 115

Other things we can do

We can also handle generalized eigenvalue problem companion pencil matrix polynomial eigenvalue problems (with L. Robol) symmetric-plus-rank-one (via Cayley transform) completely unstructured problems generalizations of Hessenberg form.

David S. Watkins Core-Chasing Algorithms

slide-116
SLIDE 116

Other things we can do

We can also handle generalized eigenvalue problem companion pencil matrix polynomial eigenvalue problems (with L. Robol) symmetric-plus-rank-one (via Cayley transform) completely unstructured problems generalizations of Hessenberg form.

Monograph in progress (100+ pp?)

David S. Watkins Core-Chasing Algorithms

slide-117
SLIDE 117

Beyond Hessenberg form

Form of Q can be upper Hessenberg, so-called CMV, . . .

  • David S. Watkins

Core-Chasing Algorithms

slide-118
SLIDE 118

Beyond Hessenberg form

Form of Q can be upper Hessenberg, so-called CMV, . . .

  • r “random” twisted . . .
  • David S. Watkins

Core-Chasing Algorithms

slide-119
SLIDE 119

Backward Stability

David S. Watkins Core-Chasing Algorithms

slide-120
SLIDE 120

Backward Stability

All of these algorithms are normwise backward stable.

David S. Watkins Core-Chasing Algorithms

slide-121
SLIDE 121

Backward Stability

All of these algorithms are normwise backward stable. This is “obvious”,

David S. Watkins Core-Chasing Algorithms

slide-122
SLIDE 122

Backward Stability

All of these algorithms are normwise backward stable. This is “obvious”, but our storage schemes are unorthodox.

David S. Watkins Core-Chasing Algorithms

slide-123
SLIDE 123

Backward Stability

All of these algorithms are normwise backward stable. This is “obvious”, but our storage schemes are unorthodox. So let’s take a closer look.

David S. Watkins Core-Chasing Algorithms

slide-124
SLIDE 124

Backward Stability

All of these algorithms are normwise backward stable. This is “obvious”, but our storage schemes are unorthodox. So let’s take a closer look. In exact arithmetic . . . ˆ A = U∗AU

David S. Watkins Core-Chasing Algorithms

slide-125
SLIDE 125

Backward Stability

All of these algorithms are normwise backward stable. This is “obvious”, but our storage schemes are unorthodox. So let’s take a closer look. In exact arithmetic . . . ˆ A = U∗AU In floating-point arithmetic . . . ˆ A = U∗(A + E)U Is it true that E ≈ uA? (u is unit roundoff.)

David S. Watkins Core-Chasing Algorithms

slide-126
SLIDE 126

Backward Stability

All of these algorithms are normwise backward stable. This is “obvious”, but our storage schemes are unorthodox. So let’s take a closer look. In exact arithmetic . . . ˆ A = U∗AU In floating-point arithmetic . . . ˆ A = U∗(A + E)U Is it true that E ≈ uA? (u is unit roundoff.) If so, the algorithm is normwise backward stable.

David S. Watkins Core-Chasing Algorithms

slide-127
SLIDE 127

Backward Stability

A = QR =

  • David S. Watkins

Core-Chasing Algorithms

slide-128
SLIDE 128

Backward Stability

A = QR =

  • ˆ

A = ˆ Q ˆ R = U∗QRU = (U∗QV )(V ∗RU)

David S. Watkins Core-Chasing Algorithms

slide-129
SLIDE 129

Backward Stability

A = QR =

  • ˆ

A = ˆ Q ˆ R = U∗QRU = (U∗QV )(V ∗RU) ˆ Q = U∗QV ˆ R = V ∗RU

David S. Watkins Core-Chasing Algorithms

slide-130
SLIDE 130

Backward Stability

A = QR =

  • ˆ

A = ˆ Q ˆ R = U∗QRU = (U∗QV )(V ∗RU) ˆ Q = U∗QV ˆ R = V ∗RU ˆ Q = U∗(Q + Eq)V ˆ R = V ∗(R + Er)U

David S. Watkins Core-Chasing Algorithms

slide-131
SLIDE 131

Backward Stability

A = QR =

  • ˆ

A = ˆ Q ˆ R = U∗QRU = (U∗QV )(V ∗RU) ˆ Q = U∗QV ˆ R = V ∗RU ˆ Q = U∗(Q + Eq)V ˆ R = V ∗(R + Er)U Eq ≈ u Er ≈ u R = u A . . . assuming R is stored in the conventional way.

David S. Watkins Core-Chasing Algorithms

slide-132
SLIDE 132

Backward Stability

ˆ A = ˆ Q ˆ R = U∗QVV ∗RU = U∗QRU = U∗AU

David S. Watkins Core-Chasing Algorithms

slide-133
SLIDE 133

Backward Stability

ˆ A = ˆ Q ˆ R = U∗QVV ∗RU = U∗QRU = U∗AU ˆ A = ˆ Q ˆ R = U∗(Q + Eq)(R + Er)U = U∗(QR + EqR + QEr + EqEr)U = U∗(A + E)U, where E = EqR + QEr + O(u2), so E ≈ u R .

David S. Watkins Core-Chasing Algorithms

slide-134
SLIDE 134

Backward Stability

Now consider the unitary-plus-rank-one case: R = C ∗(B + e1yT).

David S. Watkins Core-Chasing Algorithms

slide-135
SLIDE 135

Backward Stability

Now consider the unitary-plus-rank-one case: R = C ∗(B + e1yT). Is it still true that ˆ R = V ∗(R + Er)U with Er ≈ uR ?

David S. Watkins Core-Chasing Algorithms

slide-136
SLIDE 136

Backward Stability

Now consider the unitary-plus-rank-one case: R = C ∗(B + e1yT). Is it still true that ˆ R = V ∗(R + Er)U with Er ≈ uR ? R = C ∗(B + e1yT) = C ∗B + C ∗e1yT

David S. Watkins Core-Chasing Algorithms

slide-137
SLIDE 137

Backward Stability

Now consider the unitary-plus-rank-one case: R = C ∗(B + e1yT). Is it still true that ˆ R = V ∗(R + Er)U with Er ≈ uR ? R = C ∗(B + e1yT) = C ∗B + C ∗e1yT Ru = C ∗B Ro = C ∗e1yT

David S. Watkins Core-Chasing Algorithms

slide-138
SLIDE 138

Backward Stability

Now consider the unitary-plus-rank-one case: R = C ∗(B + e1yT). Is it still true that ˆ R = V ∗(R + Er)U with Er ≈ uR ? R = C ∗(B + e1yT) = C ∗B + C ∗e1yT Ru = C ∗B Ro = C ∗e1yT ˆ Ru = V ∗RuU ˆ Ro = V ∗RoU

David S. Watkins Core-Chasing Algorithms

slide-139
SLIDE 139

Backward Stability

Let’s consider the unitary part first.

David S. Watkins Core-Chasing Algorithms

slide-140
SLIDE 140

Backward Stability

Let’s consider the unitary part first. ˆ Ru = V ∗RuU = V ∗C ∗BU = (V ∗C ∗W )(W ∗BU) = ˆ C ∗ ˆ B

David S. Watkins Core-Chasing Algorithms

slide-141
SLIDE 141

Backward Stability

Let’s consider the unitary part first. ˆ Ru = V ∗RuU = V ∗C ∗BU = (V ∗C ∗W )(W ∗BU) = ˆ C ∗ ˆ B ˆ C = W ∗CV ˆ B = W ∗BU

David S. Watkins Core-Chasing Algorithms

slide-142
SLIDE 142

Backward Stability

Let’s consider the unitary part first. ˆ Ru = V ∗RuU = V ∗C ∗BU = (V ∗C ∗W )(W ∗BU) = ˆ C ∗ ˆ B ˆ C = W ∗CV ˆ B = W ∗BU ˆ C = W ∗(C + Ec)V ˆ B = W ∗(B + Eb)U Ec ≈ u Eb ≈ u

David S. Watkins Core-Chasing Algorithms

slide-143
SLIDE 143

Backward Stability

Let’s consider the unitary part first. ˆ Ru = V ∗RuU = V ∗C ∗BU = (V ∗C ∗W )(W ∗BU) = ˆ C ∗ ˆ B ˆ C = W ∗CV ˆ B = W ∗BU ˆ C = W ∗(C + Ec)V ˆ B = W ∗(B + Eb)U Ec ≈ u Eb ≈ u Put them together!

David S. Watkins Core-Chasing Algorithms

slide-144
SLIDE 144

Backward Stability

Now consider the rank-one part.

David S. Watkins Core-Chasing Algorithms

slide-145
SLIDE 145

Backward Stability

Now consider the rank-one part. Ro = C ∗e1yT,

David S. Watkins Core-Chasing Algorithms

slide-146
SLIDE 146

Backward Stability

Now consider the rank-one part. Ro = C ∗e1yT, where yT = −ρ−1eT

n+1C ∗B

David S. Watkins Core-Chasing Algorithms

slide-147
SLIDE 147

Backward Stability

Now consider the rank-one part. Ro = C ∗e1yT, where yT = −ρ−1eT

n+1C ∗B

so Ro = −ρ−1C ∗e1eT

n+1C ∗B

David S. Watkins Core-Chasing Algorithms

slide-148
SLIDE 148

Backward Stability

Now consider the rank-one part. Ro = C ∗e1yT, where yT = −ρ−1eT

n+1C ∗B

so Ro = −ρ−1C ∗e1eT

n+1C ∗B

ˆ Ro = V ∗RoU = −ρ−1V ∗C ∗e1eT

n+1C ∗BU

David S. Watkins Core-Chasing Algorithms

slide-149
SLIDE 149

Backward Stability

Now consider the rank-one part. Ro = C ∗e1yT, where yT = −ρ−1eT

n+1C ∗B

so Ro = −ρ−1C ∗e1eT

n+1C ∗B

ˆ Ro = V ∗RoU = −ρ−1V ∗C ∗e1eT

n+1C ∗BU

ˆ Ro = V ∗RoU = −ρ−1V ∗(C + Ec)∗e1eT

n+1(C + Ec)∗(B + Eb)U

David S. Watkins Core-Chasing Algorithms

slide-150
SLIDE 150

Backward Stability

Now consider the rank-one part. Ro = C ∗e1yT, where yT = −ρ−1eT

n+1C ∗B

so Ro = −ρ−1C ∗e1eT

n+1C ∗B

ˆ Ro = V ∗RoU = −ρ−1V ∗C ∗e1eT

n+1C ∗BU

ˆ Ro = V ∗RoU = −ρ−1V ∗(C + Ec)∗e1eT

n+1(C + Ec)∗(B + Eb)U

Backward stability follows.

David S. Watkins Core-Chasing Algorithms

slide-151
SLIDE 151

Advantages

Advantages of our approach

David S. Watkins Core-Chasing Algorithms

slide-152
SLIDE 152

Advantages

Advantages of our approach

Fast algorithms for certain special structures

David S. Watkins Core-Chasing Algorithms

slide-153
SLIDE 153

Advantages

Advantages of our approach

Fast algorithms for certain special structures

unitary unitary-plus-rank-one symmetric symmetric-plus-rank-one

David S. Watkins Core-Chasing Algorithms

slide-154
SLIDE 154

Advantages

Advantages of our approach

Fast algorithms for certain special structures

unitary unitary-plus-rank-one symmetric symmetric-plus-rank-one

. . . but also recommended for problems with no special structure.

David S. Watkins Core-Chasing Algorithms

slide-155
SLIDE 155

Advantages

Advantages of our approach

Fast algorithms for certain special structures

unitary unitary-plus-rank-one symmetric symmetric-plus-rank-one

. . . but also recommended for problems with no special structure.

chase multiple bulges for parallelism, BLAS3 speed (in theory)

David S. Watkins Core-Chasing Algorithms

slide-156
SLIDE 156

Advantages

Advantages of our approach

Fast algorithms for certain special structures

unitary unitary-plus-rank-one symmetric symmetric-plus-rank-one

. . . but also recommended for problems with no special structure.

chase multiple bulges for parallelism, BLAS3 speed (in theory) better deflation strategy

David S. Watkins Core-Chasing Algorithms

slide-157
SLIDE 157

Deflation

David S. Watkins Core-Chasing Algorithms

slide-158
SLIDE 158

Deflation

Conventional deflation: Decouple when |aj+1,j | < u (|ajj | + |aj+1,j+1 |).

David S. Watkins Core-Chasing Algorithms

slide-159
SLIDE 159

Deflation

Conventional deflation: Decouple when |aj+1,j | < u (|ajj | + |aj+1,j+1 |). Our deflation:

  • Q1Q2 · · · Qn−1

Qj = cj −sj sj cj

  • Decouple when |sj | < u.

David S. Watkins Core-Chasing Algorithms

slide-160
SLIDE 160

Deflation

Conventional deflation: Decouple when |aj+1,j | < u (|ajj | + |aj+1,j+1 |). Our deflation:

  • Q1Q2 · · · Qn−1

Qj = cj −sj sj cj

  • Decouple when |sj | < u.

Stronger backward stability (Mach and Vandebril 2014)

David S. Watkins Core-Chasing Algorithms

slide-161
SLIDE 161

Quick Summary

David S. Watkins Core-Chasing Algorithms

slide-162
SLIDE 162

Quick Summary

We took a new look at Francis’s algorithm.

David S. Watkins Core-Chasing Algorithms

slide-163
SLIDE 163

Quick Summary

We took a new look at Francis’s algorithm. considered QR decomposed form.

David S. Watkins Core-Chasing Algorithms

slide-164
SLIDE 164

Quick Summary

We took a new look at Francis’s algorithm. considered QR decomposed form. chased cores instead of bulges.

David S. Watkins Core-Chasing Algorithms

slide-165
SLIDE 165

Quick Summary

We took a new look at Francis’s algorithm. considered QR decomposed form. chased cores instead of bulges. We demonstrated some advantages.

David S. Watkins Core-Chasing Algorithms

slide-166
SLIDE 166

Quick Summary

We took a new look at Francis’s algorithm. considered QR decomposed form. chased cores instead of bulges. We demonstrated some advantages.

Thank you for your attention.

David S. Watkins Core-Chasing Algorithms