MATLAB p ( x ) = x n + a n − 1 x n − 1 + a n − 2 x n − 2 + · · · + a 0 = 0 monic polynomial companion matrix 0 0 · · · − a 0 1 0 0 − a 1 · · · . . ... . . A = 1 . . ... 0 − a n − 2 1 − a n − 1 balance, then . . . . . . get the zeros of p by computing the eigenvalues. David S. Watkins Happy Birthday, Volker!
MATLAB p ( x ) = x n + a n − 1 x n − 1 + a n − 2 x n − 2 + · · · + a 0 = 0 monic polynomial companion matrix 0 0 · · · − a 0 1 0 0 − a 1 · · · . . ... . . A = 1 . . ... 0 − a n − 2 1 − a n − 1 balance, then . . . . . . get the zeros of p by computing the eigenvalues. This is not always the best thing to do. David S. Watkins Happy Birthday, Volker!
Chebfun p ( x ) = T n ( x ) + b n − 1 T n − 1 ( x ) + · · · b 0 T 0 ( x ) Chebyshev polynomials David S. Watkins Happy Birthday, Volker!
Chebfun p ( x ) = T n ( x ) + b n − 1 T n − 1 ( x ) + · · · b 0 T 0 ( x ) Chebyshev polynomials colleague matrix David S. Watkins Happy Birthday, Volker!
Chebfun p ( x ) = T n ( x ) + b n − 1 T n − 1 ( x ) + · · · b 0 T 0 ( x ) Chebyshev polynomials colleague matrix This is sometimes better. David S. Watkins Happy Birthday, Volker!
What we’ve been doing David S. Watkins Happy Birthday, Volker!
What we’ve been doing companion matrix or David S. Watkins Happy Birthday, Volker!
What we’ve been doing companion matrix or companion pencil p ( x ) = a n x n + a n − 1 x n − 1 + a n − 2 x n − 2 + · · · + a 0 1 0 0 · · · 0 · · · 0 − a 0 1 0 0 · · · 1 0 0 − a 1 · · · . . ... . . ... . . . . . . − λ 1 . . ... 0 − a n − 2 1 0 1 − a n − 1 a n David S. Watkins Happy Birthday, Volker!
What we’ve been doing companion matrix or companion pencil p ( x ) = a n x n + a n − 1 x n − 1 + a n − 2 x n − 2 + · · · + a 0 1 0 0 · · · 0 · · · 0 − a 0 1 0 0 · · · 1 0 0 − a 1 · · · . . ... . . ... . . . . . . − λ 1 . . ... 0 − a n − 2 1 0 1 − a n − 1 a n . . . and variants. David S. Watkins Happy Birthday, Volker!
What we’ve been doing companion matrix or companion pencil p ( x ) = a n x n + a n − 1 x n − 1 + a n − 2 x n − 2 + · · · + a 0 1 0 0 · · · 0 · · · 0 − a 0 1 0 0 · · · 1 0 0 − a 1 · · · . . ... . . ... . . . . . . − λ 1 . . ... 0 − a n − 2 1 0 1 − a n − 1 a n . . . and variants. Today we restrict attention to the companion matrix. David S. Watkins Happy Birthday, Volker!
Cost of solving companion eigenvalue problem David S. Watkins Happy Birthday, Volker!
Cost of solving companion eigenvalue problem If structure not exploited: O ( n 2 ) storage, O ( n 3 ) flops Francis’s implicitly-shifted QR algorithm David S. Watkins Happy Birthday, Volker!
Cost of solving companion eigenvalue problem If structure not exploited: O ( n 2 ) storage, O ( n 3 ) flops Francis’s implicitly-shifted QR algorithm If structure exploited: O ( n 2 ) flops O ( n ) storage, data-sparse representation + Francis’s algorithm David S. Watkins Happy Birthday, Volker!
Cost of solving companion eigenvalue problem If structure not exploited: O ( n 2 ) storage, O ( n 3 ) flops Francis’s implicitly-shifted QR algorithm If structure exploited: O ( n 2 ) flops O ( n ) storage, data-sparse representation + Francis’s algorithm Several methods proposed David S. Watkins Happy Birthday, Volker!
Some of the Competitors Chandrasekaran, Gu, Xia, Zhu (2007) Bini, Boito, Eidelman, Gemignani, Gohberg (2010) Boito, Eidelman, Gemignani, Gohberg (2012) David S. Watkins Happy Birthday, Volker!
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 Happy Birthday, Volker!
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 Happy Birthday, Volker!
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 Happy Birthday, Volker!
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 Happy Birthday, Volker!
Our Contribution We present Yet another O ( n ) representation David S. Watkins Happy Birthday, Volker!
Our Contribution We present Yet another O ( n ) representation Francis algorithm in O ( n ) flops/iteration David S. Watkins Happy Birthday, Volker!
Our Contribution We present Yet another O ( n ) representation Francis algorithm in O ( n ) flops/iteration Fortran codes (we’re faster) David S. Watkins Happy Birthday, Volker!
Our Contribution We present Yet another O ( n ) representation Francis algorithm in O ( n ) flops/iteration Fortran codes (we’re faster) normwise backward stable (We can prove it.) David S. Watkins Happy Birthday, Volker!
Structure Companion matrix is unitary-plus-rank-one − e i θ − a 0 e i θ 0 0 0 0 · · · · · · 1 0 0 0 − a 1 + . . . . ... . . . . . . . . 1 0 0 0 − a n − 1 · · · David S. Watkins Happy Birthday, Volker!
Structure Companion matrix is unitary-plus-rank-one − e i θ − a 0 e i θ 0 0 0 0 · · · · · · 1 0 0 0 − a 1 + . . . . ... . . . . . . . . 1 0 0 0 − a n − 1 · · · preserved by unitary similarities David S. Watkins Happy Birthday, Volker!
Structure Companion matrix is unitary-plus-rank-one − e i θ − a 0 e i θ 0 0 0 0 · · · · · · 1 0 0 0 − a 1 + . . . . ... . . . . . . . . 1 0 0 0 − a n − 1 · · · preserved by unitary similarities Companion matrix is also upper Hessenberg. David S. Watkins Happy Birthday, Volker!
Structure Companion matrix is unitary-plus-rank-one − e i θ − a 0 e i θ 0 0 0 0 · · · · · · 1 0 0 0 − a 1 + . . . . ... . . . . . . . . 1 0 0 0 − a n − 1 · · · preserved by unitary similarities Companion matrix is also upper Hessenberg. preserved by Francis algorithm David S. Watkins Happy Birthday, Volker!
Structure Companion matrix is unitary-plus-rank-one − e i θ − a 0 e i θ 0 0 0 0 · · · · · · 1 0 0 0 − a 1 + . . . . ... . . . . . . . . 1 0 0 0 − a n − 1 · · · preserved by unitary similarities Companion matrix is also upper Hessenberg. preserved by Francis algorithm We exploit this structure. David S. Watkins Happy Birthday, Volker!
Structure Chandrasekaran, Gu, Xia, Zhu (2007) A = QR David S. Watkins Happy Birthday, Volker!
Structure Chandrasekaran, Gu, Xia, Zhu (2007) A = QR Q is upper Hessenberg and unitary. R is upper triangular and unitary-plus-rank-one. David S. Watkins Happy Birthday, Volker!
Structure Chandrasekaran, Gu, Xia, Zhu (2007) A = QR Q is upper Hessenberg and unitary. R is upper triangular and unitary-plus-rank-one. We do this too. David S. Watkins Happy Birthday, Volker!
The Unitary Part x x x x x x 1 1 1 x x x x x x x x = x x x 1 x x x x x x 1 1 x x David S. Watkins Happy Birthday, Volker!
The Unitary Part x x x x x x 1 1 1 x x x x x x x x = x x x 1 x x x x x x 1 1 x x � � � Q = � � � David S. Watkins Happy Birthday, Volker!
The Unitary Part x x x x x x 1 1 1 x x x x x x x x = x x x 1 x x x x x x 1 1 x x � � � Q = � � � O ( n ) storage David S. Watkins Happy Birthday, Volker!
The Upper Triangular Part R = U + xy T unitary-plus-rank-one, so David S. Watkins Happy Birthday, Volker!
The Upper Triangular Part R = U + xy T unitary-plus-rank-one, so R has quasiseparable rank 2. x · · · x x · · · x . . . ... . . . . . . x x x · · · R = x x · · · . ... . . x David S. Watkins Happy Birthday, Volker!
The Upper Triangular Part R = U + xy T unitary-plus-rank-one, so R has quasiseparable rank 2. x · · · x x · · · x . . . ... . . . . . . x x x · · · R = x x · · · . ... . . x quasiseparable generator representation ( O ( n ) storage) David S. Watkins Happy Birthday, Volker!
The Upper Triangular Part R = U + xy T unitary-plus-rank-one, so R has quasiseparable rank 2. x · · · x x · · · x . . . ... . . . . . . x x x · · · R = x x · · · . ... . . x quasiseparable generator representation ( O ( n ) storage) Chandrasekaran et. al. exploit this structure. David S. Watkins Happy Birthday, Volker!
The Upper Triangular Part R = U + xy T unitary-plus-rank-one, so R has quasiseparable rank 2. x · · · x x · · · x . . . ... . . . . . . x x x · · · R = x x · · · . ... . . x quasiseparable generator representation ( O ( n ) storage) Chandrasekaran et. al. exploit this structure. We do it differently. David S. Watkins Happy Birthday, Volker!
Our Representation Add a row/column for extra wiggle room 0 1 − a 0 1 − a 1 0 . . ... . . A = . . 1 0 − a n − 1 0 0 Extra zero root can be deflated immediately. David S. Watkins Happy Birthday, Volker!
Our Representation Add a row/column for extra wiggle room 0 1 − a 0 1 − a 1 0 . . ... . . A = . . 1 0 − a n − 1 0 0 Extra zero root can be deflated immediately. A = QR , where 0 ± 1 0 1 0 − a 1 . . ... 1 0 0 . . . . . . ... . . Q = R = . . 1 − a n − 1 0 1 0 0 ∓ 1 ± a 0 0 1 0 0 David S. Watkins Happy Birthday, Volker!
Our Representation 0 ± 1 0 1 0 0 . . ... . . Q = . . 1 0 0 0 1 Q is stored in factored form � � � Q = � � � David S. Watkins Happy Birthday, Volker!
Our Representation 0 ± 1 0 1 0 0 . . ... . . Q = . . 1 0 0 0 1 Q is stored in factored form � � � Q = � � � Q = Q 1 Q 2 · · · Q n − 1 David S. Watkins Happy Birthday, Volker!
Our Representation 1 0 − a 1 . . ... . . . . R = 1 − a n − 1 0 ∓ 1 ± a 0 0 0 R is unitary-plus-rank-one: 1 0 0 0 − a 1 0 . . . . ... ... . . . . . . . . + 1 0 0 0 0 − a n − 1 0 ∓ 1 ± a 0 0 ± 1 0 ∓ 1 0 David S. Watkins Happy Birthday, Volker!
Representation of R R = U + xy T , where David S. Watkins Happy Birthday, Volker!
Representation of R R = U + xy T , where − a 1 . . . xy T = � � 0 0 1 0 · · · − a n − 1 ± a 0 ∓ 1 David S. Watkins Happy Birthday, Volker!
Representation of R R = U + xy T , where − a 1 . . . xy T = � � 0 0 1 0 · · · − a n − 1 ± a 0 ∓ 1 Next step: Roll up x . David S. Watkins Happy Birthday, Volker!
Representation of R x x x x = x x x x David S. Watkins Happy Birthday, Volker!
Representation of R x x x x = x x � � x 0 David S. Watkins Happy Birthday, Volker!
Representation of R x x x x � = � 0 x � � x 0 David S. Watkins Happy Birthday, Volker!
Representation of R x x � � x 0 � = � 0 x � � x 0 David S. Watkins Happy Birthday, Volker!
Representation of R x x � � 0 x � = � x 0 � � 0 x C 1 · · · C n − 1 C n x = α e 1 (w.l.g. α = 1) David S. Watkins Happy Birthday, Volker!
Representation of R C 1 · · · C n − 1 C n x = e 1 David S. Watkins Happy Birthday, Volker!
Representation of R C 1 · · · C n − 1 C n x = e 1 Cx = e 1 David S. Watkins Happy Birthday, Volker!
Representation of R C 1 · · · C n − 1 C n x = e 1 Cx = e 1 C ∗ e 1 = x David S. Watkins Happy Birthday, Volker!
Representation of R C 1 · · · C n − 1 C n x = e 1 Cx = e 1 C ∗ e 1 = x R = U + xy T = U + C ∗ e 1 y T = C ∗ ( CU + e 1 y T ) David S. Watkins Happy Birthday, Volker!
Representation of R C 1 · · · C n − 1 C n x = e 1 Cx = e 1 C ∗ e 1 = x R = U + xy T = U + C ∗ e 1 y T = C ∗ ( CU + e 1 y T ) R = C ∗ ( B + e 1 y T ) David S. Watkins Happy Birthday, Volker!
Representation of R C 1 · · · C n − 1 C n x = e 1 Cx = e 1 C ∗ e 1 = x R = U + xy T = U + C ∗ e 1 y T = C ∗ ( CU + e 1 y T ) R = C ∗ ( B + e 1 y T ) B is upper Hessenberg (and unitary) David S. Watkins Happy Birthday, Volker!
Recommend
More recommend