Fast computation of eigenvalues of companion, comrade, and related - - PowerPoint PPT Presentation

▶
fast computation of eigenvalues of companion comrade and
SMART_READER_LITE
LIVE PREVIEW

Fast computation of eigenvalues of companion, comrade, and related - - PowerPoint PPT Presentation

Apr 08, 2024 32 likes •1.4k views

Fast computation of eigenvalues of companion, comrade, and related matrices David S. Watkins Department of Mathematics Washington State University Providence, June 2013 David S. Watkins companion, comrade, and related matrices Collaborators

slide-1
SLIDE 1

Fast computation of eigenvalues of companion, comrade, and related matrices

David S. Watkins

Department of Mathematics Washington State University

Providence, June 2013

David S. Watkins companion, comrade, and related matrices

slide-2
SLIDE 2

Collaborators

This is joint work with Raf Vandebril (KU Leuven) Jared Aurentz (WSU)

David S. Watkins companion, comrade, and related matrices

slide-3
SLIDE 3

Zeros of a polynomial (monomial basis)

p(z) = zn − c1zn−1 − c2zn−2 − · · · − cn = 0 companion matrix A =        c1 c2 · · · cn 1 1 ... 1        Get the zeros of p by computing the eigenvalues of A. Example: MATLAB’s roots command

David S. Watkins companion, comrade, and related matrices

slide-4
SLIDE 4

Zeros of a polynomial (monomial basis)

p(z) = zn − c1zn−1 − c2zn−2 − · · · − cn = 0 companion matrix A =        c1 c2 · · · cn 1 1 ... 1        Get the zeros of p by computing the eigenvalues of A. Example: MATLAB’s roots command

David S. Watkins companion, comrade, and related matrices

slide-5
SLIDE 5

Zeros of a polynomial (monomial basis)

p(z) = zn − c1zn−1 − c2zn−2 − · · · − cn = 0 companion matrix A =        c1 c2 · · · cn 1 1 ... 1        Get the zeros of p by computing the eigenvalues of A. Example: MATLAB’s roots command

David S. Watkins companion, comrade, and related matrices

slide-6
SLIDE 6

Zeros of a polynomial (monomial basis)

p(z) = zn − c1zn−1 − c2zn−2 − · · · − cn = 0 companion matrix A =        c1 c2 · · · cn 1 1 ... 1        Get the zeros of p by computing the eigenvalues of A. Example: MATLAB’s roots command

David S. Watkins companion, comrade, and related matrices

slide-7
SLIDE 7

Zeros of a polynomial (monomial basis)

p(z) = zn − c1zn−1 − c2zn−2 − · · · − cn = 0 companion matrix A =        c1 c2 · · · cn 1 1 ... 1        Get the zeros of p by computing the eigenvalues of A. Example: MATLAB’s roots command

David S. Watkins companion, comrade, and related matrices

slide-8
SLIDE 8

Zeros of a polynomial (“arbitrary” basis)

p(z) = pn(z) − c1pn−1(z) − c2pn−2(z) − · · · − cnp0(z) = 0 confederate matrix (Barnett 1983) A =         × × × × × × × × × × × × × × × × × × × × × × × × × ×         pk(z)bk+1,k = zpk−1(z) −

k

  • j=1

pj−1(z)bjk Short recurrence implies sparse matrix.

David S. Watkins companion, comrade, and related matrices

slide-9
SLIDE 9

Zeros of a polynomial (“arbitrary” basis)

p(z) = pn(z) − c1pn−1(z) − c2pn−2(z) − · · · − cnp0(z) = 0 confederate matrix (Barnett 1983) A =         × × × × × × × × × × × × × × × × × × × × × × × × × ×         pk(z)bk+1,k = zpk−1(z) −

k

  • j=1

pj−1(z)bjk Short recurrence implies sparse matrix.

David S. Watkins companion, comrade, and related matrices

slide-10
SLIDE 10

Zeros of a polynomial (“arbitrary” basis)

p(z) = pn(z) − c1pn−1(z) − c2pn−2(z) − · · · − cnp0(z) = 0 confederate matrix (Barnett 1983) A =         × × × × × × × × × × × × × × × × × × × × × × × × × ×         pk(z)bk+1,k = zpk−1(z) −

k

  • j=1

pj−1(z)bjk Short recurrence implies sparse matrix.

David S. Watkins companion, comrade, and related matrices

slide-11
SLIDE 11

Zeros of a polynomial (“arbitrary” basis)

p(z) = pn(z) − c1pn−1(z) − c2pn−2(z) − · · · − cnp0(z) = 0 confederate matrix (Barnett 1983) A =         × × × × × × × × × × × × × × × × × × × × × × × × × ×         pk(z)bk+1,k = zpk−1(z) −

k

  • j=1

pj−1(z)bjk Short recurrence implies sparse matrix.

David S. Watkins companion, comrade, and related matrices

slide-12
SLIDE 12

Zeros of a polynomial (“arbitrary” basis)

p(z) = pn(z) − c1pn−1(z) − c2pn−2(z) − · · · − cnp0(z) = 0 confederate matrix (Barnett 1983) A =         × × × × × × × × × × × × × × × × × × × × × × × × × ×         pk(z)bk+1,k = zpk−1(z) −

k

  • j=1

pj−1(z)bjk Short recurrence implies sparse matrix.

David S. Watkins companion, comrade, and related matrices

slide-13
SLIDE 13

Zeros of a polynomial (3-term recurrence)

pk(z)αk = (z − βk)pk−1(z) − γkpk−2(z) p(z) = pn(z) − c1pn−1(z) − c2pn−2(z) − · · · − cnp0(z) = 0 comrade matrix A =         × × × × × × × × × × × × × × × × × × × ×         tridiagonal plus spike Examples: Chebyshev, Legendre, Hermite, . . . Example: Chebfun’s roots command

David S. Watkins companion, comrade, and related matrices

slide-14
SLIDE 14

Zeros of a polynomial (3-term recurrence)

pk(z)αk = (z − βk)pk−1(z) − γkpk−2(z) p(z) = pn(z) − c1pn−1(z) − c2pn−2(z) − · · · − cnp0(z) = 0 comrade matrix A =         × × × × × × × × × × × × × × × × × × × ×         tridiagonal plus spike Examples: Chebyshev, Legendre, Hermite, . . . Example: Chebfun’s roots command

David S. Watkins companion, comrade, and related matrices

slide-15
SLIDE 15

Zeros of a polynomial (3-term recurrence)

pk(z)αk = (z − βk)pk−1(z) − γkpk−2(z) p(z) = pn(z) − c1pn−1(z) − c2pn−2(z) − · · · − cnp0(z) = 0 comrade matrix A =         × × × × × × × × × × × × × × × × × × × ×         tridiagonal plus spike Examples: Chebyshev, Legendre, Hermite, . . . Example: Chebfun’s roots command

David S. Watkins companion, comrade, and related matrices

slide-16
SLIDE 16

Zeros of a polynomial (3-term recurrence)

pk(z)αk = (z − βk)pk−1(z) − γkpk−2(z) p(z) = pn(z) − c1pn−1(z) − c2pn−2(z) − · · · − cnp0(z) = 0 comrade matrix A =         × × × × × × × × × × × × × × × × × × × ×         tridiagonal plus spike Examples: Chebyshev, Legendre, Hermite, . . . Example: Chebfun’s roots command

David S. Watkins companion, comrade, and related matrices

slide-17
SLIDE 17

Zeros of a polynomial (2-term recurrence)

Newton basis pk(z)αk = (z − βk)pk−1(z) p(z) = pn(z) − c1pn−1(z) − c2pn−2(z) − · · · − cnp0(z) = 0 A =         × × × × × × × × × × × × × × × ×         bidiagonal plus spike

David S. Watkins companion, comrade, and related matrices

slide-18
SLIDE 18

Zeros of a polynomial (2-term recurrence)

Newton basis pk(z)αk = (z − βk)pk−1(z) p(z) = pn(z) − c1pn−1(z) − c2pn−2(z) − · · · − cnp0(z) = 0 A =         × × × × × × × × × × × × × × × ×         bidiagonal plus spike

David S. Watkins companion, comrade, and related matrices

slide-19
SLIDE 19

Zeros of a polynomial (2-term recurrence)

Newton basis pk(z)αk = (z − βk)pk−1(z) p(z) = pn(z) − c1pn−1(z) − c2pn−2(z) − · · · − cnp0(z) = 0 A =         × × × × × × × × × × × × × × × ×         bidiagonal plus spike

David S. Watkins companion, comrade, and related matrices

slide-20
SLIDE 20

Zeros of a polynomial (1-term recurrence)

monomial basis pk(z) = zpk−1(z) p(z) = pn(z) − c1pn−1(z) − c2pn−2(z) − · · · − cnp0(z) = 0 A =         × × × × × × × × × × ×        

  • ne diagonal plus spike

companion matrix

David S. Watkins companion, comrade, and related matrices

slide-21
SLIDE 21

Zeros of a polynomial (1-term recurrence)

monomial basis pk(z) = zpk−1(z) p(z) = pn(z) − c1pn−1(z) − c2pn−2(z) − · · · − cnp0(z) = 0 A =         × × × × × × × × × × ×        

  • ne diagonal plus spike

companion matrix

David S. Watkins companion, comrade, and related matrices

slide-22
SLIDE 22

How do we compute the eigenvalues?

narrow band plus spike A =         × × × × × × × × × × × × × × × × × × × ×         Francis’s bulge-chasing algorithm (“implicitly-shifted QR”) O(n3) time, O(n2) space Is there a more economical approach?

David S. Watkins companion, comrade, and related matrices

slide-23
SLIDE 23

How do we compute the eigenvalues?

narrow band plus spike A =         × × × × × × × × × × × × × × × × × × × ×         Francis’s bulge-chasing algorithm (“implicitly-shifted QR”) O(n3) time, O(n2) space Is there a more economical approach?

David S. Watkins companion, comrade, and related matrices

slide-24
SLIDE 24

How do we compute the eigenvalues?

narrow band plus spike A =         × × × × × × × × × × × × × × × × × × × ×         Francis’s bulge-chasing algorithm (“implicitly-shifted QR”) O(n3) time, O(n2) space Is there a more economical approach?

David S. Watkins companion, comrade, and related matrices

slide-25
SLIDE 25

How do we compute the eigenvalues?

narrow band plus spike A =         × × × × × × × × × × × × × × × × × × × ×         Francis’s bulge-chasing algorithm (“implicitly-shifted QR”) O(n3) time, O(n2) space Is there a more economical approach?

David S. Watkins companion, comrade, and related matrices

slide-26
SLIDE 26

How do we compute the eigenvalues?

banded plus spike A =         × × × × × × × × × × × × × × × × × × × ×         Exploit rank structure (generators) We pursue a different approach.

David S. Watkins companion, comrade, and related matrices

slide-27
SLIDE 27

How do we compute the eigenvalues?

banded plus spike A =         × × × × × × × × × × × × × × × × × × × ×         Exploit rank structure (generators) We pursue a different approach.

David S. Watkins companion, comrade, and related matrices

slide-28
SLIDE 28

How do we compute the eigenvalues?

banded plus spike A =         × × × × × × × × × × × × × × × × × × × ×         Exploit rank structure (generators) We pursue a different approach.

David S. Watkins companion, comrade, and related matrices

slide-29
SLIDE 29

Fast Methods We have developed great methods that are . . .

lightning fast unstable

David S. Watkins companion, comrade, and related matrices

slide-30
SLIDE 30

Fast Methods We have developed great methods that are . . .

lightning fast unstable

David S. Watkins companion, comrade, and related matrices

slide-31
SLIDE 31

Fast Methods We have developed great methods that are . . .

lightning fast unstable

David S. Watkins companion, comrade, and related matrices

slide-32
SLIDE 32

A useful factorization (by Gaussian elimination)

      × × × × × × × × × × × × × × × ×      

David S. Watkins companion, comrade, and related matrices

slide-33
SLIDE 33

A useful factorization (by Gaussian elimination)

     × × × × × × × × × × × × × × × ×       =       × × × × × × × × × × × × × × × ×      

David S. Watkins companion, comrade, and related matrices

slide-34
SLIDE 34

A useful factorization (by Gaussian elimination)

     × × × × × × × × × × × × × × × ×       =       × × × × × × × × × × × × × × ×      

David S. Watkins companion, comrade, and related matrices

slide-35
SLIDE 35

A useful factorization (by Gaussian elimination)

     × × × × × × × × × × × × × × × ×       =       × × × × × × × × × × × × × × ×      

David S. Watkins companion, comrade, and related matrices

slide-36
SLIDE 36

A useful factorization (by Gaussian elimination)

     × × × × × × × × × × × × × × × ×       =       × × × × × × × × × × × × × ×      

David S. Watkins companion, comrade, and related matrices

slide-37
SLIDE 37

A useful factorization (by Gaussian elimination)

     × × × × × × × × × × × × × × × ×       =       × × × × × × × × × × × × × ×      

David S. Watkins companion, comrade, and related matrices

slide-38
SLIDE 38

A useful factorization (by Gaussian elimination)

     × × × × × × × × × × × × × × × ×       =       × × × × × × × × × × × × ×      

David S. Watkins companion, comrade, and related matrices

slide-39
SLIDE 39

A useful factorization (by Gaussian elimination)

     × × × × × × × × × × × × × × × ×       =       × × × × × × × × × × × × ×       and so on . . . to (banded) triangular form.

David S. Watkins companion, comrade, and related matrices

slide-40
SLIDE 40

A useful factorization (by Gaussian elimination)

We get the factorization A =

     × × × × × × × × × × ×       core transformations times banded upper triangular band width = length of recurrence

David S. Watkins companion, comrade, and related matrices

slide-41
SLIDE 41

What core transformations look like

  • =

   × × × × 1 1       1 × × × × 1       1 1 × × × ×   

David S. Watkins companion, comrade, and related matrices

slide-42
SLIDE 42

Bulge Chasing Algorithm

     × × × × × × × × × × ×       We will develop a nonunitary bulge-chasing algorithm that preserves this factorization. Note: Storage is O(n) (about 6n for this case). We have single-shift and double-shift variants.

David S. Watkins companion, comrade, and related matrices

slide-43
SLIDE 43

Bulge Chasing Algorithm

     × × × × × × × × × × ×       We will develop a nonunitary bulge-chasing algorithm that preserves this factorization. Note: Storage is O(n) (about 6n for this case). We have single-shift and double-shift variants.

David S. Watkins companion, comrade, and related matrices

slide-44
SLIDE 44

Bulge Chasing Algorithm

     × × × × × × × × × × ×       We will develop a nonunitary bulge-chasing algorithm that preserves this factorization. Note: Storage is O(n) (about 6n for this case). We have single-shift and double-shift variants.

David S. Watkins companion, comrade, and related matrices

slide-45
SLIDE 45

Bulge Chasing Algorithm

     × × × × × × × × × × ×       We will develop a nonunitary bulge-chasing algorithm that preserves this factorization. Note: Storage is O(n) (about 6n for this case). We have single-shift and double-shift variants.

David S. Watkins companion, comrade, and related matrices

slide-46
SLIDE 46

Bulge Chasing Algorithm (single-shift variant)

     × × × × × × × × × × ×      

David S. Watkins companion, comrade, and related matrices

slide-47
SLIDE 47

Bulge Chasing Algorithm (single-shift variant)

     × × × × × × × × × × ×      

  • David S. Watkins

companion, comrade, and related matrices

slide-48
SLIDE 48

Bulge Chasing Algorithm (single-shift variant)

     × × × × × × × × × × ×      

  • David S. Watkins

companion, comrade, and related matrices

slide-49
SLIDE 49

Passing Gauss transform through triangular matrix

Blue core transformations are Gauss transforms. × × × × ×

  • =

× × × × × ×

  • =
  • ×

× × × ×

  • We do this at every opportunity.

From this point on, we suppress the triangular matrix.

David S. Watkins companion, comrade, and related matrices

slide-50
SLIDE 50

Passing Gauss transform through triangular matrix

Blue core transformations are Gauss transforms. × × × × ×

  • =

× × × × × ×

  • =
  • ×

× × × ×

  • We do this at every opportunity.

From this point on, we suppress the triangular matrix.

David S. Watkins companion, comrade, and related matrices

slide-51
SLIDE 51

Passing Gauss transform through triangular matrix

Blue core transformations are Gauss transforms. × × × × ×

  • =

× × × × × ×

  • =
  • ×

× × × ×

  • We do this at every opportunity.

From this point on, we suppress the triangular matrix.

David S. Watkins companion, comrade, and related matrices

slide-52
SLIDE 52

Bulge Chasing Algorithm (single-shift variant)

     × × × × × × × × × × ×      

  • David S. Watkins

companion, comrade, and related matrices

slide-53
SLIDE 53

Bulge Chasing Algorithm (single-shift variant)

     × × × × × × × × × × ×      

David S. Watkins companion, comrade, and related matrices

slide-54
SLIDE 54

Bulge Chasing Algorithm (single-shift variant)

  • David S. Watkins

companion, comrade, and related matrices

slide-55
SLIDE 55

Bulge Chasing Algorithm (single-shift variant)

  • David S. Watkins

companion, comrade, and related matrices

slide-56
SLIDE 56

Bulge Chasing Algorithm (single-shift variant)

  • David S. Watkins

companion, comrade, and related matrices

slide-57
SLIDE 57

Bulge Chasing Algorithm (single-shift variant)

  • David S. Watkins

companion, comrade, and related matrices

slide-58
SLIDE 58

Turnover operation

Can we do this?

  • ⇔
  • crucial question

Pessimistic answer: No! (not always) Optimistic answer: Yes! (almost always)

David S. Watkins companion, comrade, and related matrices

slide-59
SLIDE 59

Turnover operation

Can we do this?

  • ⇔
  • crucial question

Pessimistic answer: No! (not always) Optimistic answer: Yes! (almost always)

David S. Watkins companion, comrade, and related matrices

slide-60
SLIDE 60

Turnover operation

Can we do this?

  • ⇔
  • crucial question

Pessimistic answer: No! (not always) Optimistic answer: Yes! (almost always)

David S. Watkins companion, comrade, and related matrices

slide-61
SLIDE 61

Turnover operation

Can we do this?

  • ⇔
  • crucial question

Pessimistic answer: No! (not always) Optimistic answer: Yes! (almost always)

David S. Watkins companion, comrade, and related matrices

slide-62
SLIDE 62

Turnover operation

  • =

× × × × × × × ×

  • =

× × × × × × × × ×

  • Red submatrix still has rank one.

David S. Watkins companion, comrade, and related matrices

slide-63
SLIDE 63

Turnover operation

  • =

× × × × × × × ×

  • =

× × × × × × × × ×

  • Red submatrix still has rank one.

David S. Watkins companion, comrade, and related matrices

slide-64
SLIDE 64

Turnover operation

So . . . × × × × × × × × ×

  • =
  • ×

× × × × × × ×

  • =
  • Gauss transform does not disturb 2 × 2 submatrix.

David S. Watkins companion, comrade, and related matrices

slide-65
SLIDE 65

Turnover operation

So . . . × × × × × × × × ×

  • =
  • ×

× × × × × × ×

  • =
  • Gauss transform does not disturb 2 × 2 submatrix.

David S. Watkins companion, comrade, and related matrices

slide-66
SLIDE 66

Turnover operation

So . . . × × × × × × × × ×

  • =
  • ×

× × × × × × ×

  • =
  • Gauss transform does not disturb 2 × 2 submatrix.

David S. Watkins companion, comrade, and related matrices

slide-67
SLIDE 67

Turnover operation

Complete turnover looks like this:

  • =
  • David S. Watkins

companion, comrade, and related matrices

slide-68
SLIDE 68

Bulge Chasing Algorithm (single-shift variant)

  • David S. Watkins

companion, comrade, and related matrices

slide-69
SLIDE 69

Bulge Chasing Algorithm (single-shift variant)

  • David S. Watkins

companion, comrade, and related matrices

slide-70
SLIDE 70

Bulge Chasing Algorithm (single-shift variant)

  • David S. Watkins

companion, comrade, and related matrices

slide-71
SLIDE 71

Bulge Chasing Algorithm (single-shift variant)

  • David S. Watkins

companion, comrade, and related matrices

slide-72
SLIDE 72

Bulge Chasing Algorithm (single-shift variant)

  • David S. Watkins

companion, comrade, and related matrices

slide-73
SLIDE 73

Bulge Chasing Algorithm (single-shift variant)

  • David S. Watkins

companion, comrade, and related matrices

slide-74
SLIDE 74

Bulge Chasing Algorithm (single-shift variant)

  • David S. Watkins

companion, comrade, and related matrices

slide-75
SLIDE 75

Bulge Chasing Algorithm (single-shift variant)

  • David S. Watkins

companion, comrade, and related matrices

slide-76
SLIDE 76

Bulge Chasing Algorithm (single-shift variant)

  • David S. Watkins

companion, comrade, and related matrices

slide-77
SLIDE 77

Bulge Chasing Algorithm (single-shift variant)

  • David S. Watkins

companion, comrade, and related matrices

slide-78
SLIDE 78

Bulge Chasing Algorithm (single-shift variant)

  • David S. Watkins

companion, comrade, and related matrices

slide-79
SLIDE 79

Bulge Chasing Algorithm (single-shift variant)

  • David S. Watkins

companion, comrade, and related matrices

slide-80
SLIDE 80

Bulge Chasing Algorithm (single-shift variant)

  • Done!

Now repeat again, again, and again

David S. Watkins companion, comrade, and related matrices

slide-81
SLIDE 81

Bulge Chasing Algorithm (single-shift variant)

  • Done!

Now repeat again, again, and again

David S. Watkins companion, comrade, and related matrices

slide-82
SLIDE 82

Bulge Chasing Algorithm (single-shift variant)

  • Done!

Now repeat again, again, and again

David S. Watkins companion, comrade, and related matrices

slide-83
SLIDE 83

Bulge Chasing Algorithm (single-shift variant)

  • Done!

Now repeat again, again, and again

David S. Watkins companion, comrade, and related matrices

slide-84
SLIDE 84

Bulge Chasing Algorithm (single-shift variant)

  • Done!

Now repeat again, again, and again

David S. Watkins companion, comrade, and related matrices

slide-85
SLIDE 85

Computational complexity

O(n) storage (everything in cache) O(n) flops per iteration O(n) iterations O(n2) flops total

David S. Watkins companion, comrade, and related matrices

slide-86
SLIDE 86

Computational complexity

O(n) storage (everything in cache) O(n) flops per iteration O(n) iterations O(n2) flops total

David S. Watkins companion, comrade, and related matrices

slide-87
SLIDE 87

Computational complexity

O(n) storage (everything in cache) O(n) flops per iteration O(n) iterations O(n2) flops total

David S. Watkins companion, comrade, and related matrices

slide-88
SLIDE 88

Computational complexity

O(n) storage (everything in cache) O(n) flops per iteration O(n) iterations O(n2) flops total

David S. Watkins companion, comrade, and related matrices

slide-89
SLIDE 89

Computational complexity

O(n) storage (everything in cache) O(n) flops per iteration O(n) iterations O(n2) flops total

David S. Watkins companion, comrade, and related matrices

slide-90
SLIDE 90

Stability issue

This is a structured, implicitly-shifted LR algorithm. Unstable! Can we trust the results? No, but we can do a posteriori tests. (λI − A)v = |p(λ)| | p(λ)

p′(λ) | is estimate of error.

Newton correction Cost: O(n2) flops for n roots.

David S. Watkins companion, comrade, and related matrices

slide-91
SLIDE 91

Stability issue

This is a structured, implicitly-shifted LR algorithm. Unstable! Can we trust the results? No, but we can do a posteriori tests. (λI − A)v = |p(λ)| | p(λ)

p′(λ) | is estimate of error.

Newton correction Cost: O(n2) flops for n roots.

David S. Watkins companion, comrade, and related matrices

slide-92
SLIDE 92

Stability issue

This is a structured, implicitly-shifted LR algorithm. Unstable! Can we trust the results? No, but we can do a posteriori tests. (λI − A)v = |p(λ)| | p(λ)

p′(λ) | is estimate of error.

Newton correction Cost: O(n2) flops for n roots.

David S. Watkins companion, comrade, and related matrices

slide-93
SLIDE 93

Stability issue

This is a structured, implicitly-shifted LR algorithm. Unstable! Can we trust the results? No, but we can do a posteriori tests. (λI − A)v = |p(λ)| | p(λ)

p′(λ) | is estimate of error.

Newton correction Cost: O(n2) flops for n roots.

David S. Watkins companion, comrade, and related matrices

slide-94
SLIDE 94

Stability issue

This is a structured, implicitly-shifted LR algorithm. Unstable! Can we trust the results? No, but we can do a posteriori tests. (λI − A)v = |p(λ)| | p(λ)

p′(λ) | is estimate of error.

Newton correction Cost: O(n2) flops for n roots.

David S. Watkins companion, comrade, and related matrices

slide-95
SLIDE 95

Stability issue

This is a structured, implicitly-shifted LR algorithm. Unstable! Can we trust the results? No, but we can do a posteriori tests. (λI − A)v = |p(λ)| | p(λ)

p′(λ) | is estimate of error.

Newton correction Cost: O(n2) flops for n roots.

David S. Watkins companion, comrade, and related matrices

slide-96
SLIDE 96

Stability issue

This is a structured, implicitly-shifted LR algorithm. Unstable! Can we trust the results? No, but we can do a posteriori tests. (λI − A)v = |p(λ)| | p(λ)

p′(λ) | is estimate of error.

Newton correction Cost: O(n2) flops for n roots.

David S. Watkins companion, comrade, and related matrices

slide-97
SLIDE 97

Stability issue

This is a structured, implicitly-shifted LR algorithm. Unstable! Can we trust the results? No, but we can do a posteriori tests. (λI − A)v = |p(λ)| | p(λ)

p′(λ) | is estimate of error.

Newton correction Cost: O(n2) flops for n roots.

David S. Watkins companion, comrade, and related matrices

slide-98
SLIDE 98

Stability issue

This is a structured, implicitly-shifted LR algorithm. Unstable! Can we trust the results? No, but we can do a posteriori tests. (λI − A)v = |p(λ)| | p(λ)

p′(λ) | is estimate of error.

Newton correction Cost: O(n2) flops for n roots.

David S. Watkins companion, comrade, and related matrices

slide-99
SLIDE 99

Stability issue

This is a structured, implicitly-shifted LR algorithm. Unstable! Can we trust the results? No, but we can do a posteriori tests. (λI − A)v = |p(λ)| | p(λ)

p′(λ) | is estimate of error.

Newton correction Cost: O(n2) flops for n roots.

David S. Watkins companion, comrade, and related matrices

slide-100
SLIDE 100

Performance

First example: companion matrices, complex, single-shift code

Contestants

LAPACK code ZHSEQR (O(n3)) BBEGG (quasiseparable, unitary) AVW

  • ur code

David S. Watkins companion, comrade, and related matrices

slide-101
SLIDE 101

Speed Comparison, companion case

10

1

10

2

10

3

10

−4

10

−3

10

−2

10

−1

10 10

1

10

2

Degree Elapsed time LAPACK BBEGG AVW David S. Watkins companion, comrade, and related matrices

slide-102
SLIDE 102

Speed Comparison, companion case

AVW times include a posteriori test times. No crossover point Slope = 1.84 At n = 1210: method time LAPACK 28.2 BBEGG 6.5 AVW 0.5

David S. Watkins companion, comrade, and related matrices

slide-103
SLIDE 103

Speed Comparison, companion case

AVW times include a posteriori test times. No crossover point Slope = 1.84 At n = 1210: method time LAPACK 28.2 BBEGG 6.5 AVW 0.5

David S. Watkins companion, comrade, and related matrices

slide-104
SLIDE 104

Speed Comparison, companion case

AVW times include a posteriori test times. No crossover point Slope = 1.84 At n = 1210: method time LAPACK 28.2 BBEGG 6.5 AVW 0.5

David S. Watkins companion, comrade, and related matrices

slide-105
SLIDE 105

Accuracy Comparison, companion case

Maximum of residual (λI − A)v /Av degree 17 72 295 1210 LAPACK 1.1 × 10−14 3.1 × 10−14 6.8 × 10−14 2.3 × 10−13 BBEGG 1.5 × 10−14 1.1 × 10−13 4.1 × 10−12 2.4 × 10−11 AVW 7.4 × 10−12 5.0 × 10−11 3.7 × 10−10 1.4 × 10−09

David S. Watkins companion, comrade, and related matrices

slide-106
SLIDE 106

Accuracy Comparison, companion case

Error estimate |p(λ)/p′(λ)| degree 17 72 295 1210 LAPACK 4.3 × 10−15 8.5 × 10−15 1.4 × 10−14 2.4 × 10−14 BBEGG 2.4 × 10−14 1.8 × 10−13 1.3 × 10−12 1.2 × 10−11 AVW 4.7 × 10−12 7.5 × 10−11 5.4 × 10−11 1.7 × 10−10

David S. Watkins companion, comrade, and related matrices

slide-107
SLIDE 107

Accuracy Comparison, companion case

After Newton correction Maximum of residual (λI − A)v /Av degree 17 72 295 1210 LAPACK 8.1 × 10−16 1.8 × 10−15 2.6 × 10−15 5.2 × 10−15 BBEGG 7.5 × 10−16 1.5 × 10−15 3.5 × 10−15 5.5 × 10−15 AVW 8.6 × 10−16 1.4 × 10−15 2.6 × 10−15 5.3 × 10−15

David S. Watkins companion, comrade, and related matrices

slide-108
SLIDE 108

Accuracy Comparison, companion case

After Newton correction Error estimate |p(λ)/p′(λ)| degree 17 72 295 1210 LAPACK 3.2 × 10−16 3.4 × 10−16 3.2 × 10−16 3.1 × 10−16 BBEGG 5.2 × 10−16 3.1 × 10−16 2.9 × 10−16 3.8 × 10−16 AVW 3.3 × 10−16 3.1 × 10−16 3.1 × 10−16 3.1 × 10−16

David S. Watkins companion, comrade, and related matrices

slide-109
SLIDE 109

Speed Comparison, colleague case

Second example: Chebyshev basis, complex, single-shift code

3-term recurrence comrade matrix, tridiagonal plus spike colleague matrix (Good 1961)

David S. Watkins companion, comrade, and related matrices

slide-110
SLIDE 110

Speed Comparison, colleague case

Second example: Chebyshev basis, complex, single-shift code

3-term recurrence comrade matrix, tridiagonal plus spike colleague matrix (Good 1961)

David S. Watkins companion, comrade, and related matrices

slide-111
SLIDE 111

Speed Comparison, colleague case

Second example: Chebyshev basis, complex, single-shift code

3-term recurrence comrade matrix, tridiagonal plus spike colleague matrix (Good 1961)

David S. Watkins companion, comrade, and related matrices

slide-112
SLIDE 112

Speed Comparison, colleague case

Second example: Chebyshev basis, complex, single-shift code

3-term recurrence comrade matrix, tridiagonal plus spike colleague matrix (Good 1961)

David S. Watkins companion, comrade, and related matrices

slide-113
SLIDE 113

Speed Comparison, colleague case

10

1

10

2

10

3

10

−4

10

−3

10

−2

10

−1

10 10

1

10

2

Degree Elapsed Time LAPACK AVW David S. Watkins companion, comrade, and related matrices

slide-114
SLIDE 114

Speed Comparison, colleague case

No crossover point Slope = 1.76 At n = 1210: method time LAPACK 31.4 AVW 1.1

David S. Watkins companion, comrade, and related matrices

slide-115
SLIDE 115

Speed Comparison, colleague case

No crossover point Slope = 1.76 At n = 1210: method time LAPACK 31.4 AVW 1.1

David S. Watkins companion, comrade, and related matrices

slide-116
SLIDE 116

Speed Comparison, colleague case

No crossover point Slope = 1.76 At n = 1210: method time LAPACK 31.4 AVW 1.1

David S. Watkins companion, comrade, and related matrices

slide-117
SLIDE 117

Accuracy Comparison, colleague case

Maximum of residual (λI − A)v /Av degree 17 72 295 1210 LAPACK 1.7 × 10−13 1.1 × 10−12 6.1 × 10−12 2.8 × 10−11 AVW 7.7 × 10−12 2.8 × 10−10 1.3 × 10−08 3.3 × 10−07

David S. Watkins companion, comrade, and related matrices

slide-118
SLIDE 118

Accuracy Comparison, colleague case

Error estimate |p(λ)/p′(λ)| degree 17 72 295 1210 LAPACK 1.9 × 10−14 2.7 × 10−14 3.1 × 10−14 5.1 × 10−14 AVW 3.8 × 10−13 7.9 × 10−12 3.9 × 10−11 1.2 × 10−10

David S. Watkins companion, comrade, and related matrices

slide-119
SLIDE 119

Accuracy Comparison, colleague case

After Newton correction Maximum of residual (λI − A)v /Av degree 17 72 295 1210 LAPACK 7.7 × 10−15 2.7 × 10−14 1.7 × 10−13 9.0 × 10−13 AVW 7.4 × 10−15 2.7 × 10−14 1.4 × 10−13 9.9 × 10−13

David S. Watkins companion, comrade, and related matrices

slide-120
SLIDE 120

Accuracy Comparison, colleague case

After Newton correction Error estimate |p(λ)/p′(λ)| degree 17 72 295 1210 LAPACK 9.8 × 10−16 6.6 × 10−16 4.9 × 10−16 2.1 × 10−16 AVW 1.7 × 10−15 3.2 × 10−16 2.1 × 10−16 3.4 × 10−16

David S. Watkins companion, comrade, and related matrices

slide-121
SLIDE 121

Speed Comparison, real colleague case

Third example: Chebyshev basis, real, double-shift code

Contestants

LAPACK code DHSEQR (O(n3)) EGG (symmetric plus rank-one, backward stable) AVW

  • ur code (double-shift)

David S. Watkins companion, comrade, and related matrices

slide-122
SLIDE 122

Speed Comparison, real colleague case

10

1

10

2

10

3

10

−4

10

−3

10

−2

10

−1

10 10

1

10

2

Degree Elapsed Time LAPACK AVW EGG David S. Watkins companion, comrade, and related matrices

slide-123
SLIDE 123

Speed Comparison, real colleague case

We are just a bit faster than EGG At n = 1210: method time LAPACK 9.00 EGG 0.64 AVW 0.49

David S. Watkins companion, comrade, and related matrices

slide-124
SLIDE 124

Speed Comparison, real colleague case

We are just a bit faster than EGG At n = 1210: method time LAPACK 9.00 EGG 0.64 AVW 0.49

David S. Watkins companion, comrade, and related matrices

slide-125
SLIDE 125

Speed Comparison, real colleague case

We are just a bit faster than EGG At n = 1210: method time LAPACK 9.00 EGG 0.64 AVW 0.49

David S. Watkins companion, comrade, and related matrices

slide-126
SLIDE 126

Accuracy Comparison, real colleague case

Maximum of residual (λI − A)v /Av degree 17 72 295 1210 LAPACK 1.7 × 10−13 1.7 × 10−12 6.3 × 10−12 5.4 × 10−11 EGG 1.9 × 10−13 2.1 × 10−12 7.4 × 10−12 1.1 × 10−10 AVW 4.3 × 10−09 1.2 × 10−06 4.6 × 10−07 2.4 × 10−03

David S. Watkins companion, comrade, and related matrices

slide-127
SLIDE 127

Accuracy Comparison, real colleague case

Maximum of residual (λI − A)v /Av degree 17 72 295 1210 LAPACK 1.7 × 10−13 1.7 × 10−12 6.3 × 10−12 5.4 × 10−11 EGG 1.9 × 10−13 2.1 × 10−12 7.4 × 10−12 1.1 × 10−10 AVW 4.3 × 10−09 1.2 × 10−06 4.6 × 10−07 2.4 × 10−03

Ouch!

David S. Watkins companion, comrade, and related matrices

slide-128
SLIDE 128

Accuracy Comparison, real colleague case

After Newton correction Maximum of residual (λI − A)v /Av degree 17 72 295 1210 LAPACK 8.2 × 10−15 3.5 × 10−14 1.0 × 10−13 1.6 × 10−12 EGG 6.5 × 10−15 3.3 × 10−14 1.3 × 10−13 2.1 × 10−12 AVW 5.8 × 10−15 5.8 × 10−13 5.0 × 10−12 4.3 × 10−05 A second Newton correction does the trick.

David S. Watkins companion, comrade, and related matrices

slide-129
SLIDE 129

Accuracy Comparison, real colleague case

After Newton correction Maximum of residual (λI − A)v /Av degree 17 72 295 1210 LAPACK 8.2 × 10−15 3.5 × 10−14 1.0 × 10−13 1.6 × 10−12 EGG 6.5 × 10−15 3.3 × 10−14 1.3 × 10−13 2.1 × 10−12 AVW 5.8 × 10−15 5.8 × 10−13 5.0 × 10−12 4.3 × 10−05 A second Newton correction does the trick.

David S. Watkins companion, comrade, and related matrices

slide-130
SLIDE 130

Final Remarks

Our code is lightning fast, but it is unstable. Newton correction is helpful. Double-shift code disappointing. These are benign examples. Hybrid methods?

David S. Watkins companion, comrade, and related matrices

slide-131
SLIDE 131

Final Remarks

Our code is lightning fast, but it is unstable. Newton correction is helpful. Double-shift code disappointing. These are benign examples. Hybrid methods?

David S. Watkins companion, comrade, and related matrices

slide-132
SLIDE 132

Final Remarks

Our code is lightning fast, but it is unstable. Newton correction is helpful. Double-shift code disappointing. These are benign examples. Hybrid methods?

David S. Watkins companion, comrade, and related matrices

slide-133
SLIDE 133

Final Remarks

Our code is lightning fast, but it is unstable. Newton correction is helpful. Double-shift code disappointing. These are benign examples. Hybrid methods?

David S. Watkins companion, comrade, and related matrices

slide-134
SLIDE 134

Final Remarks

Our code is lightning fast, but it is unstable. Newton correction is helpful. Double-shift code disappointing. These are benign examples. Hybrid methods?

David S. Watkins companion, comrade, and related matrices

slide-135
SLIDE 135

Final Remarks

Our code is lightning fast, but it is unstable. Newton correction is helpful. Double-shift code disappointing. These are benign examples. Hybrid methods?

David S. Watkins companion, comrade, and related matrices

slide-136
SLIDE 136

Final Remarks

Our code is lightning fast, but it is unstable. Newton correction is helpful. Double-shift code disappointing. These are benign examples. Hybrid methods?

David S. Watkins companion, comrade, and related matrices

slide-137
SLIDE 137

Final Remarks

Our code is lightning fast, but it is unstable. Newton correction is helpful. Double-shift code disappointing. These are benign examples. Hybrid methods?

Thank you for your attention.

David S. Watkins companion, comrade, and related matrices