Franciss Algorithm as a Core-Chasing Algorithm David S. Watkins - - PowerPoint PPT Presentation

▶
francis s algorithm as a core chasing algorithm
SMART_READER_LITE
LIVE PREVIEW

Franciss Algorithm as a Core-Chasing Algorithm David S. Watkins - - PowerPoint PPT Presentation

Aug 30, 2023 629 likes •2.53k views

Franciss Algorithm as a Core-Chasing Algorithm David S. Watkins Department of Mathematics Washington State University PNWNAS, November 12, 2016 David S. Watkins Core-Chasing Algorithm Todays Topic David S. Watkins Core-Chasing

slide-1
SLIDE 1

Francis’s Algorithm as a Core-Chasing Algorithm

David S. Watkins

Department of Mathematics Washington State University

PNWNAS, November 12, 2016

David S. Watkins Core-Chasing Algorithm

slide-2
SLIDE 2

Today’s Topic

David S. Watkins Core-Chasing Algorithm

slide-3
SLIDE 3

Today’s Topic

The matrix eigenvalue problem

David S. Watkins Core-Chasing Algorithm

slide-4
SLIDE 4

Today’s Topic

The matrix eigenvalue problem A ∈ Cn×n

David S. Watkins Core-Chasing Algorithm

slide-5
SLIDE 5

Today’s Topic

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

David S. Watkins Core-Chasing Algorithm

slide-6
SLIDE 6

Today’s Topic

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

David S. Watkins Core-Chasing Algorithm

slide-7
SLIDE 7

Today’s Topic

The matrix eigenvalue problem A ∈ Cn×n Find the eigenvalues (. . . vectors, invariant subspaces) Many applications Interest dates back to the very beginning of the electronic computing era.

David S. Watkins Core-Chasing Algorithm

slide-8
SLIDE 8

Today’s Topic

The matrix eigenvalue problem A ∈ Cn×n Find the eigenvalues (. . . vectors, invariant subspaces) Many applications Interest dates back to the very beginning of the electronic computing era. Nobody knew how to do it.

David S. Watkins Core-Chasing Algorithm

slide-9
SLIDE 9

John Francis

David S. Watkins Core-Chasing Algorithm

slide-10
SLIDE 10

John Francis

invented the winning algorithm in 1959.

David S. Watkins Core-Chasing Algorithm

slide-11
SLIDE 11

John Francis

invented the winning algorithm in 1959. commonly called: QR algorithm

David S. Watkins Core-Chasing Algorithm

slide-12
SLIDE 12

John Francis

invented the winning algorithm in 1959. commonly called: QR algorithm more precisely: implicitly shifted QR algorithm

David S. Watkins Core-Chasing Algorithm

slide-13
SLIDE 13

John Francis

invented the winning algorithm in 1959. commonly called: QR algorithm more precisely: implicitly shifted QR algorithm better yet: Francis’s algorithm,

David S. Watkins Core-Chasing Algorithm

slide-14
SLIDE 14

John Francis

invented the winning algorithm in 1959. commonly called: QR algorithm more precisely: implicitly shifted QR algorithm better yet: Francis’s algorithm, bulge-chasing algorithm.

David S. Watkins Core-Chasing Algorithm

slide-15
SLIDE 15

Francis’s algorithm (superficial description)

David S. Watkins Core-Chasing Algorithm

slide-16
SLIDE 16

Francis’s algorithm (superficial description)

upper Hessenberg form A =       ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗       unitary similarity transformation direct method (O(n3) flops)

David S. Watkins Core-Chasing Algorithm

slide-17
SLIDE 17

Francis’s algorithm (superficial description)

upper Hessenberg form A =       ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗       unitary similarity transformation direct method (O(n3) flops) Francis: Iterate

David S. Watkins Core-Chasing Algorithm

slide-18
SLIDE 18

Francis’s algorithm (superficial description)

upper Hessenberg form A =       ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗       unitary similarity transformation direct method (O(n3) flops) Francis: Iterate Drive toward triangular form.

David S. Watkins Core-Chasing Algorithm

slide-19
SLIDE 19

Francis’s algorithm (superficial description)

upper Hessenberg form A =       ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗       unitary similarity transformation direct method (O(n3) flops) Francis: Iterate Drive toward triangular form. (Galois theory)

David S. Watkins Core-Chasing Algorithm

slide-20
SLIDE 20

Francis’s algorithm (superficial description)

Chasing the bulge       ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      

David S. Watkins Core-Chasing Algorithm

slide-21
SLIDE 21

Francis’s algorithm (superficial description)

Chasing the bulge       ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ + ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      

David S. Watkins Core-Chasing Algorithm

slide-22
SLIDE 22

Francis’s algorithm (superficial description)

Chasing the bulge       ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ + ∗ ∗ ∗ ∗ ∗      

David S. Watkins Core-Chasing Algorithm

slide-23
SLIDE 23

Francis’s algorithm (superficial description)

Chasing the bulge       ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ + ∗ ∗      

David S. Watkins Core-Chasing Algorithm

slide-24
SLIDE 24

Francis’s algorithm (superficial description)

Chasing the bulge       ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      

David S. Watkins Core-Chasing Algorithm

slide-25
SLIDE 25

Francis’s algorithm (superficial description)

Chasing the bulge       ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗       iteration complete!

David S. Watkins Core-Chasing Algorithm

slide-26
SLIDE 26

Francis’s algorithm (superficial description)

Chasing the bulge       ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗       iteration complete! repeated iterations ⇒ triangular form

David S. Watkins Core-Chasing Algorithm

slide-27
SLIDE 27

Francis’s algorithm (superficial description)

Chasing the bulge       ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗       iteration complete! repeated iterations ⇒ triangular form This is the single-shift algorithm.

David S. Watkins Core-Chasing Algorithm

slide-28
SLIDE 28

Francis’s algorithm (superficial description)

Chasing the bulge       ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗       iteration complete! repeated iterations ⇒ triangular form This is the single-shift algorithm. Double-shift algorithm chases a bigger bulge.

David S. Watkins Core-Chasing Algorithm

slide-29
SLIDE 29

Computational Cost

Computational Cost

O(n2) flops per iteration O(n) total iterations O(n3) total flops

David S. Watkins Core-Chasing Algorithm

slide-30
SLIDE 30

For details see . . .

David S. Watkins Core-Chasing Algorithm

slide-31
SLIDE 31

For details see . . .

Golub and Van Loan, Matrix Computations, 4th Ed.

David S. Watkins Core-Chasing Algorithm

slide-32
SLIDE 32

For details see . . .

Golub and Van Loan, Matrix Computations, 4th Ed. Watkins, Fundamentals of Matrix Computations, 3rd Ed.

David S. Watkins Core-Chasing Algorithm

slide-33
SLIDE 33

For details see . . .

Golub and Van Loan, Matrix Computations, 4th Ed. Watkins, Fundamentals of Matrix Computations, 3rd Ed.

David S. Watkins Core-Chasing Algorithm

slide-34
SLIDE 34

My History with this Topic

David S. Watkins Core-Chasing Algorithm

slide-35
SLIDE 35

My History with this Topic

Understanding the QR algorithm, SIAM Rev., 1982

David S. Watkins Core-Chasing Algorithm

slide-36
SLIDE 36

My History with this Topic

Understanding the QR algorithm, SIAM Rev., 1982 Fundamentals of Matrix Computations, Wiley, 1991

David S. Watkins Core-Chasing Algorithm

slide-37
SLIDE 37

My History with this Topic

Understanding the QR algorithm, SIAM Rev., 1982 Fundamentals of Matrix Computations, Wiley, 1991 Some perspectives on the eigenvalue problem, 1993 QR-like algorithms—an overview of convergence theory and practice, AMS proceedings, 1996 QR-like algorithms for eigenvalue problems, JCAM, 2000 The Matrix Eigenvalue Problem: GR and Krylov Subspace Methods, SIAM, 2007 The QR algorithm revisited, SIAM Rev., 2008 Fundamentals of Matrix Computations, 3rd Ed., 2010 Francis’s Algorithm, Amer. Math. Monthly, 2011

David S. Watkins Core-Chasing Algorithm

slide-38
SLIDE 38

My History with this Topic

Understanding the QR algorithm, SIAM Rev., 1982 Fundamentals of Matrix Computations, Wiley, 1991 Some perspectives on the eigenvalue problem, 1993 QR-like algorithms—an overview of convergence theory and practice, AMS proceedings, 1996 QR-like algorithms for eigenvalue problems, JCAM, 2000 The Matrix Eigenvalue Problem: GR and Krylov Subspace Methods, SIAM, 2007 The QR algorithm revisited, SIAM Rev., 2008 Fundamentals of Matrix Computations, 3rd Ed., 2010 Francis’s Algorithm, Amer. Math. Monthly, 2011 . . . but we’re still not done!

David S. Watkins Core-Chasing Algorithm

slide-39
SLIDE 39

Our International Research Group

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

David S. Watkins Core-Chasing Algorithm

slide-40
SLIDE 40

A new look at an old algorithm

David S. Watkins Core-Chasing Algorithm

slide-41
SLIDE 41

A new look at an old algorithm

Store in QR decomposed form A = QR Q is unitary, R is upper triangular

David S. Watkins Core-Chasing Algorithm

slide-42
SLIDE 42

A new look at an old algorithm

Store in QR decomposed form A = QR Q is unitary, R is upper triangular looks inefficient!

David S. Watkins Core-Chasing Algorithm

slide-43
SLIDE 43

A new look at an old algorithm

Store in QR decomposed form A = QR Q is unitary, R is upper triangular looks inefficient! but it’s not!

David S. Watkins Core-Chasing Algorithm

slide-44
SLIDE 44

A new look at an old algorithm

     ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      =      ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗     

David S. Watkins Core-Chasing Algorithm

slide-45
SLIDE 45

A new look at an old algorithm

    ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      =      ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗     

David S. Watkins Core-Chasing Algorithm

slide-46
SLIDE 46

A new look at an old algorithm

    ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      =      ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗     

David S. Watkins Core-Chasing Algorithm

slide-47
SLIDE 47

A new look at an old algorithm

    ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      =      ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗     

David S. Watkins Core-Chasing Algorithm

slide-48
SLIDE 48

A new look at an old algorithm

    ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      =      ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗     

David S. Watkins Core-Chasing Algorithm

slide-49
SLIDE 49

A new look at an old algorithm

    ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      =      ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗     

David S. Watkins Core-Chasing Algorithm

slide-50
SLIDE 50

A new look at an old algorithm

    ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      =      ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      Def: Core Transformation

David S. Watkins Core-Chasing Algorithm

slide-51
SLIDE 51

A new look at an old algorithm

    ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      =      ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      Def: Core Transformation Now invert the core transformations to move them to the other side.

David S. Watkins Core-Chasing Algorithm

slide-52
SLIDE 52

A new look at an old algorithm

     ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      =

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

David S. Watkins Core-Chasing Algorithm

slide-53
SLIDE 53

A new look at an old algorithm

     ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      =

    ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      A = QR

David S. Watkins Core-Chasing Algorithm

slide-54
SLIDE 54

A new look at an old algorithm

     ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      =

    ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗      A = QR Q =

  • Q requires only O(n) storage space.

David S. Watkins Core-Chasing Algorithm

slide-55
SLIDE 55

A new look at an old algorithm

Manipulating core transformations

David S. Watkins Core-Chasing Algorithm

slide-56
SLIDE 56

A new look at an old algorithm

Manipulating core transformations

Fusion

  • ⇒
  • David S. Watkins

Core-Chasing Algorithm

slide-57
SLIDE 57

A new look at an old algorithm

Manipulating core transformations

Fusion

  • ⇒
  • Turnover

(aka shift through, Givens swap, . . . )

  • ⇔
  • David S. Watkins

Core-Chasing Algorithm

slide-58
SLIDE 58

A new look at an old algorithm

Manipulating core transformations

Fusion

  • ⇒
  • Turnover

(aka shift through, Givens swap, . . . )

  • ⇔
  • Passing a core transformation through a triangular matrix

(cost O(n))    ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗    ⇔    ∗ ∗ ∗ ∗ ∗ ∗ ∗ + ∗ ∗ ∗    ⇔

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

David S. Watkins Core-Chasing Algorithm

slide-59
SLIDE 59

A new look at an old algorithm

David S. Watkins Core-Chasing Algorithm

slide-60
SLIDE 60

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

David S. Watkins Core-Chasing Algorithm

slide-61
SLIDE 61

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

  • David S. Watkins

Core-Chasing Algorithm

slide-62
SLIDE 62

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

  • David S. Watkins

Core-Chasing Algorithm

slide-63
SLIDE 63

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

  • David S. Watkins

Core-Chasing Algorithm

slide-64
SLIDE 64

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

David S. Watkins Core-Chasing Algorithm

slide-65
SLIDE 65

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

David S. Watkins Core-Chasing Algorithm

slide-66
SLIDE 66

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

David S. Watkins Core-Chasing Algorithm

slide-67
SLIDE 67

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

David S. Watkins Core-Chasing Algorithm

slide-68
SLIDE 68

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

  • David S. Watkins

Core-Chasing Algorithm

slide-69
SLIDE 69

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

David S. Watkins Core-Chasing Algorithm

slide-70
SLIDE 70

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

David S. Watkins Core-Chasing Algorithm

slide-71
SLIDE 71

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

David S. Watkins Core-Chasing Algorithm

slide-72
SLIDE 72

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

David S. Watkins Core-Chasing Algorithm

slide-73
SLIDE 73

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

  • David S. Watkins

Core-Chasing Algorithm

slide-74
SLIDE 74

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

David S. Watkins Core-Chasing Algorithm

slide-75
SLIDE 75

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

David S. Watkins Core-Chasing Algorithm

slide-76
SLIDE 76

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

David S. Watkins Core-Chasing Algorithm

slide-77
SLIDE 77

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

David S. Watkins Core-Chasing Algorithm

slide-78
SLIDE 78

A new look at an old algorithm

Francis’s algorithm on the QR decomposed form (a core chasing algorithm)

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

Done!

David S. Watkins Core-Chasing Algorithm

slide-79
SLIDE 79

A new look at an old algorithm Cost

David S. Watkins Core-Chasing Algorithm

slide-80
SLIDE 80

A new look at an old algorithm Cost

Most arithmetic in passing-through operation

David S. Watkins Core-Chasing Algorithm

slide-81
SLIDE 81

A new look at an old algorithm 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 Algorithm

slide-82
SLIDE 82

A new look at an old algorithm Are there any advantages?

David S. Watkins Core-Chasing Algorithm

slide-83
SLIDE 83

A new look at an old algorithm Are there any advantages?

unitary case

David S. Watkins Core-Chasing Algorithm

slide-84
SLIDE 84

A new look at an old algorithm Are there any advantages?

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

David S. Watkins Core-Chasing Algorithm

slide-85
SLIDE 85

A new look at an old algorithm Are there any advantages?

unitary case companion case (unitary-plus-rank-one) general case: efficient cache use

David S. Watkins Core-Chasing Algorithm

slide-86
SLIDE 86

Unitary Case

David S. Watkins Core-Chasing Algorithm

slide-87
SLIDE 87

Unitary Case

A = QR =

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

David S. Watkins Core-Chasing Algorithm

slide-88
SLIDE 88

Unitary Case

A = QR =

  • David S. Watkins

Core-Chasing Algorithm

slide-89
SLIDE 89

Unitary Case

A = QR =

  • David S. Watkins

Core-Chasing Algorithm

slide-90
SLIDE 90

Unitary Case

A = QR =

  • Cost is O(n) flops per iteration,

David S. Watkins Core-Chasing Algorithm

slide-91
SLIDE 91

Unitary Case

A = QR =

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

David S. Watkins Core-Chasing Algorithm

slide-92
SLIDE 92

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 Algorithm

slide-93
SLIDE 93

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 Algorithm

slide-94
SLIDE 94

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 Algorithm

slide-95
SLIDE 95

Companion Case

David S. Watkins Core-Chasing Algorithm

slide-96
SLIDE 96

Companion Case

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

David S. Watkins Core-Chasing Algorithm

slide-97
SLIDE 97

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 Algorithm

slide-98
SLIDE 98

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 Algorithm

slide-99
SLIDE 99

Cost of solving companion eigenvalue problem

David S. Watkins Core-Chasing Algorithm

slide-100
SLIDE 100

Cost of solving companion eigenvalue problem

If structure not exploited:

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

David S. Watkins Core-Chasing Algorithm

slide-101
SLIDE 101

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 data-sparse representation + Francis’s algorithm

David S. Watkins Core-Chasing Algorithm

slide-102
SLIDE 102

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 data-sparse representation + Francis’s algorithm several methods proposed

David S. Watkins Core-Chasing Algorithm

slide-103
SLIDE 103

Some of the Competitors

Chandrasekaran, Gu, Xia, Zhu (2007) Bini, Boito, Eidelman, Gemignani, Gohberg (2010) Boito, Eidelman, Gemignani, Gohberg (2012)

David S. Watkins Core-Chasing Algorithm

slide-104
SLIDE 104

Some of the Competitors

Chandrasekaran, Gu, Xia, Zhu (2007) Bini, Boito, Eidelman, Gemignani, Gohberg (2010) Boito, Eidelman, Gemignani, Gohberg (2012) Fortran codes available

David S. Watkins Core-Chasing Algorithm

slide-105
SLIDE 105

Some of the Competitors

Chandrasekaran, Gu, Xia, Zhu (2007) Bini, Boito, Eidelman, Gemignani, Gohberg (2010) Boito, Eidelman, Gemignani, Gohberg (2012) Fortran codes available evidence of backward stability

David S. Watkins Core-Chasing Algorithm

slide-106
SLIDE 106

Some of the Competitors

Chandrasekaran, Gu, Xia, Zhu (2007) Bini, Boito, Eidelman, Gemignani, Gohberg (2010) Boito, Eidelman, Gemignani, Gohberg (2012) Fortran codes available evidence of backward stability quasiseparable generator representation

David S. Watkins Core-Chasing Algorithm

slide-107
SLIDE 107

Some of the Competitors

Chandrasekaran, Gu, Xia, Zhu (2007) Bini, Boito, Eidelman, Gemignani, Gohberg (2010) Boito, Eidelman, Gemignani, Gohberg (2012) Fortran codes available evidence of backward stability quasiseparable generator representation We will do something else.

David S. Watkins Core-Chasing Algorithm

slide-108
SLIDE 108

Some of the Competitors

Chandrasekaran, Gu, Xia, Zhu (2007) Bini, Boito, Eidelman, Gemignani, Gohberg (2010) Boito, Eidelman, Gemignani, Gohberg (2012) Fortran codes available evidence of backward stability quasiseparable generator representation We will do something else. Our method is faster,

David S. Watkins Core-Chasing Algorithm

slide-109
SLIDE 109

Some of the Competitors

Chandrasekaran, Gu, Xia, Zhu (2007) Bini, Boito, Eidelman, Gemignani, Gohberg (2010) Boito, Eidelman, Gemignani, Gohberg (2012) Fortran codes available evidence of backward stability quasiseparable generator representation We will do something else. Our method is faster, and we can prove backward stability.

David S. Watkins Core-Chasing Algorithm

slide-110
SLIDE 110

Structure

David S. Watkins Core-Chasing Algorithm

slide-111
SLIDE 111

Structure

Companion matrix is unitary-plus-rank-one      · · · 1 1 ... . . . 1      +      · · · −a0 − 1 −a1 . . . . . . . . . · · · −an−1      We exploit this structure.

David S. Watkins Core-Chasing Algorithm

slide-112
SLIDE 112

Structure

. . . but we store the QR decomposed form

A = QR

David S. Watkins Core-Chasing Algorithm

slide-113
SLIDE 113

Structure

. . . but we store the QR decomposed form

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

David S. Watkins Core-Chasing Algorithm

slide-114
SLIDE 114

Structure

. . . but we store the QR decomposed form

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

  • ...

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

David S. Watkins Core-Chasing Algorithm

slide-115
SLIDE 115

Structure

How do we store R compactly?

David S. Watkins Core-Chasing Algorithm

slide-116
SLIDE 116

Structure

How do we store R compactly? R is unitary-plus-rank one.

David S. Watkins Core-Chasing Algorithm

slide-117
SLIDE 117

Structure

How do we store R compactly? R is unitary-plus-rank one. Adjoin a row and column for wiggle room. (not obvious)

David S. Watkins Core-Chasing Algorithm

slide-118
SLIDE 118

Structure

How do we store R compactly? R is unitary-plus-rank one. Adjoin a row and column for wiggle room. (not obvious) ˆ R =        1 −a1 ... . . . . . . 1 −an−1 −a0 1       

David S. Watkins Core-Chasing Algorithm

slide-119
SLIDE 119

Structure

How do we store R compactly? R is unitary-plus-rank one. Adjoin a row and column for wiggle room. (not obvious) ˆ R =        1 −a1 ... . . . . . . 1 −an−1 −a0 1        =        1 ... . . . . . . 1 1 1        +        −a1 ... . . . . . . −an−1 −a0 −1       

David S. Watkins Core-Chasing Algorithm

slide-120
SLIDE 120

Representation of R

R = PT ˆ RP

David S. Watkins Core-Chasing Algorithm

slide-121
SLIDE 121

Representation of R

R = PT ˆ RP ˆ R = U + xyT, where

David S. Watkins Core-Chasing Algorithm

slide-122
SLIDE 122

Representation of R

R = PT ˆ RP ˆ R = U + xyT, where xyT =        −a1 . . . − an−1 −a0 −1       

  • · · ·

1

  • David S. Watkins

Core-Chasing Algorithm

slide-123
SLIDE 123

Representation of R

R = PT ˆ RP ˆ R = U + xyT, where xyT =        −a1 . . . − an−1 −a0 −1       

  • · · ·

1

  • Next step: Roll up x.

David S. Watkins Core-Chasing Algorithm

slide-124
SLIDE 124

Representation of R

   x x x x    =    x x x x   

David S. Watkins Core-Chasing Algorithm

slide-125
SLIDE 125

Representation of R

  x x x x    =    x x x   

David S. Watkins Core-Chasing Algorithm

slide-126
SLIDE 126

Representation of R

  x x x x    =    x x   

David S. Watkins Core-Chasing Algorithm

slide-127
SLIDE 127

Representation of R

  x x x x    =    x   

David S. Watkins Core-Chasing Algorithm

slide-128
SLIDE 128

Representation of R

  x x x x    =    x    C1 · · · Cn−1Cnx = αe1 (w.l.g. α = 1)

David S. Watkins Core-Chasing Algorithm

slide-129
SLIDE 129

Representation of R

C1 · · · Cn−1Cnx = e1

David S. Watkins Core-Chasing Algorithm

slide-130
SLIDE 130

Representation of R

C1 · · · Cn−1Cnx = e1 Cx = e1

David S. Watkins Core-Chasing Algorithm

slide-131
SLIDE 131

Representation of R

C1 · · · Cn−1Cnx = e1 Cx = e1 C ∗e1 = x

David S. Watkins Core-Chasing Algorithm

slide-132
SLIDE 132

Representation of R

C1 · · · Cn−1Cnx = e1 Cx = e1 C ∗e1 = x ˆ R = U + xyT = U + C ∗e1yT = C ∗(CU + e1yT)

David S. Watkins Core-Chasing Algorithm

slide-133
SLIDE 133

Representation of R

C1 · · · Cn−1Cnx = e1 Cx = e1 C ∗e1 = x ˆ R = U + xyT = U + C ∗e1yT = C ∗(CU + e1yT) ˆ R = C ∗(B + e1yT)

David S. Watkins Core-Chasing Algorithm

slide-134
SLIDE 134

Representation of R

C1 · · · Cn−1Cnx = e1 Cx = e1 C ∗e1 = x ˆ R = U + xyT = U + C ∗e1yT = C ∗(CU + e1yT) ˆ R = C ∗(B + e1yT) B is upper Hessenberg (and unitary)

David S. Watkins Core-Chasing Algorithm

slide-135
SLIDE 135

Representation of R

C1 · · · Cn−1Cnx = e1 Cx = e1 C ∗e1 = x ˆ R = U + xyT = U + C ∗e1yT = C ∗(CU + e1yT) ˆ R = C ∗(B + e1yT) B is upper Hessenberg (and unitary) so B = B1 · · · Bn.

David S. Watkins Core-Chasing Algorithm

slide-136
SLIDE 136

Representation of R

C1 · · · Cn−1Cnx = e1 Cx = e1 C ∗e1 = x ˆ R = U + xyT = U + C ∗e1yT = C ∗(CU + e1yT) ˆ R = C ∗(B + e1yT) B is upper Hessenberg (and unitary) so B = B1 · · · Bn. R = PTC ∗(B + e1yT)P = PTC ∗

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

David S. Watkins Core-Chasing Algorithm

slide-137
SLIDE 137

Representation of R

C1 · · · Cn−1Cnx = e1 Cx = e1 C ∗e1 = x ˆ R = U + xyT = U + C ∗e1yT = C ∗(CU + e1yT) ˆ R = C ∗(B + e1yT) B is upper Hessenberg (and unitary) so B = B1 · · · Bn. R = PTC ∗(B + e1yT)P = PTC ∗

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

O(n) storage Bonus: Redundancy! No need to keep track of y.

David S. Watkins Core-Chasing Algorithm

slide-138
SLIDE 138

Representation of R

R = PTC ∗

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

David S. Watkins Core-Chasing Algorithm

slide-139
SLIDE 139

Representation of R

R = PTC ∗

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

=

   

  • +

· · ·     

David S. Watkins Core-Chasing Algorithm

slide-140
SLIDE 140

Representation of A Altogether we have

A = QR =

   

  • +

· · ·     

David S. Watkins Core-Chasing Algorithm

slide-141
SLIDE 141

Francis Iterations

We have complex single-shift code . . . real double-shift code.

David S. Watkins Core-Chasing Algorithm

slide-142
SLIDE 142

Francis Iterations

We have complex single-shift code . . . real double-shift code. We describe single-shift case for simplicity.

David S. Watkins Core-Chasing Algorithm

slide-143
SLIDE 143

Francis Iterations

We have complex single-shift code . . . real double-shift code. We describe single-shift case for simplicity. ignoring rank-one part . . . A =

  • David S. Watkins

Core-Chasing Algorithm

slide-144
SLIDE 144

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-145
SLIDE 145

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-146
SLIDE 146

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-147
SLIDE 147

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-148
SLIDE 148

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-149
SLIDE 149

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-150
SLIDE 150

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-151
SLIDE 151

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-152
SLIDE 152

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-153
SLIDE 153

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-154
SLIDE 154

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-155
SLIDE 155

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-156
SLIDE 156

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-157
SLIDE 157

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-158
SLIDE 158

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-159
SLIDE 159

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-160
SLIDE 160

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-161
SLIDE 161

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-162
SLIDE 162

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-163
SLIDE 163

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-164
SLIDE 164

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-165
SLIDE 165

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-166
SLIDE 166

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-167
SLIDE 167

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-168
SLIDE 168

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-169
SLIDE 169

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-170
SLIDE 170

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-171
SLIDE 171

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-172
SLIDE 172

The Core Chase

  • David S. Watkins

Core-Chasing Algorithm

slide-173
SLIDE 173

Done!

iteration complete!

David S. Watkins Core-Chasing Algorithm

slide-174
SLIDE 174

Done!

iteration complete! Cost: 3n turnovers/iteration, so O(n) flops/iteration

David S. Watkins Core-Chasing Algorithm

slide-175
SLIDE 175

Done!

iteration complete! Cost: 3n turnovers/iteration, so O(n) flops/iteration Double-shift iteration is similar. (Chase two core transformations instead of one.)

David S. Watkins Core-Chasing Algorithm

slide-176
SLIDE 176

Performance

10 10

1

10

2

10

3

10

4

10

5

10

−5

10

−4

10

−3

10

−2

10

−1

10 10

1

10

2

10

3

degree time (seconds) LAPACK BEGG AMVW David S. Watkins Core-Chasing Algorithm

slide-177
SLIDE 177

Performance

At degree 1000

method time LAPACK 7.2 BEGG 1.2 AMVW 0.2

David S. Watkins Core-Chasing Algorithm

slide-178
SLIDE 178

See our paper for . . .

Paper in SIAM J. Matrix Anal. Appl. has

David S. Watkins Core-Chasing Algorithm

slide-179
SLIDE 179

See our paper for . . .

Paper in SIAM J. Matrix Anal. Appl. has . . . more timings,

David S. Watkins Core-Chasing Algorithm

slide-180
SLIDE 180

See our paper for . . .

Paper in SIAM J. Matrix Anal. Appl. has . . . more timings, . . . accuracy comparisons,

David S. Watkins Core-Chasing Algorithm

slide-181
SLIDE 181

See our paper for . . .

Paper in SIAM J. Matrix Anal. Appl. has . . . more timings, . . . accuracy comparisons, . . . proof of backward stability.

David S. Watkins Core-Chasing Algorithm

slide-182
SLIDE 182

Summary

David S. Watkins Core-Chasing Algorithm

slide-183
SLIDE 183

Summary

We took a new look at Francis’s algorithm

David S. Watkins Core-Chasing Algorithm

slide-184
SLIDE 184

Summary

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

David S. Watkins Core-Chasing Algorithm

slide-185
SLIDE 185

Summary

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

David S. Watkins Core-Chasing Algorithm

slide-186
SLIDE 186

Summary

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

unitary case

David S. Watkins Core-Chasing Algorithm

slide-187
SLIDE 187

Summary

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

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

David S. Watkins Core-Chasing Algorithm

slide-188
SLIDE 188

Summary

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

unitary case unitary-plus-rank-one case (companion) efficient cache use (not demonstrated today)

David S. Watkins Core-Chasing Algorithm

slide-189
SLIDE 189

Summary

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

unitary case unitary-plus-rank-one case (companion) efficient cache use (not demonstrated today)

Thank you for your attention.

David S. Watkins Core-Chasing Algorithm