QR-algorithms for eigenvalue computation of structured matrices. - - PowerPoint PPT Presentation

qr algorithms for eigenvalue computation of structured
SMART_READER_LITE
LIVE PREVIEW

QR-algorithms for eigenvalue computation of structured matrices. - - PowerPoint PPT Presentation

QR-algorithms for eigenvalue computation of structured matrices. The case of low-rank corrections of unitary matrices L. Gemignani University of Pisa gemignan@dm.unipi.it http://www.dm.unipi.it/gemignan Cortona, 2004 p.1/19


slide-1
SLIDE 1

QR-algorithms for eigenvalue computation of structured matrices. The case of low-rank corrections of unitary matrices

  • L. Gemignani

University of Pisa gemignan@dm.unipi.it http://www.dm.unipi.it/˜gemignan

Cortona, 2004 – p.1/19

slide-2
SLIDE 2

The Companion Matrix

  • ✁✄✂
☎ ✆ ✝ ✞✠✟ ✡
✂ ✞ ☛ ☞ ✆ ✌ ✍✎✍✎✍✎✍ ✏ ✑ ✒ ✒ ✒ ✑ ✓

... . . .

✓ ✗✖ ✔

...

. . .

✕ ✓
  • ✝✙✘
✖ ✔
✚ ✛✎✛✎✛✎✛ ✜

Eigenvalues of

== Zeros of

✂ ☎ ☞

is upper Hessenberg MATLAB uses the shifted

✢ ✣

eigenvalue algorithm for the computation of polynomial zeros MATLAB only exploits the Hessenberg structure

Cortona, 2004 – p.2/19

slide-3
SLIDE 3

The (Shifted) QR Algorithm

☞ ✡ ✆ ☞ ☞✥✤ ✓ ✦ ✤ ✧ ✝ ✆ ✢ ✤ ✣ ✤ ☞✥✤ ★ ✖ ✩ ✆ ✣ ✤ ✢ ✤ ✪ ✦ ✤ ✧ ✝ ☛

Classical Theory: The Hessenberg structure is invariant under

✢ ✣

(

☞ ✡

Hessenberg

✫ ☞ ✖

Hessenberg)

✬ ✁✄✭ ✮ ☎

flops for each step of

✢ ✣

applied to a Hessenberg matrix The MATLAB implementation generally requires

✬ ✁✄✭ ✯ ☎

flops and

✬ ✁✄✭ ✮ ☎

storage for computing all the zeros

Cortona, 2004 – p.3/19

slide-4
SLIDE 4

The Research Problem

To design a

✢ ✣

eigenvalue algorithm for companion matrices which requires

✬ ✁✄✭ ☎

flops and

✬ ✁✄✭ ☎

storage per iteration Recent Contributions:

  • 1. [Bini, Daddi, G.]. (To appear in ETNA) FAST but

UNSTABLE

  • 2. [Chandrasekaran, Gu]. (Talk at the Workshop in Banff,

November 2003)

www.pims.math.ca/birs/workshops/2003/03w5008 FAST but

(probably) UNSTABLE The main goal: To achieve both efficiency and stability

Cortona, 2004 – p.4/19

slide-5
SLIDE 5

The Novel Approach

☞ ✆ ✌ ✍✎✍✎✍✎✍ ✏ ✑ ✒ ✒ ✒ ✑ ✕ ✕

... . . .

...

. . .

✑ ✕ ✑ ✚ ✛✎✛✎✛✎✛ ✜ ✪✱✰ ✲ ✳ ✝ ✆ ✴ ✪✱✰ ✲ ✳ ✝ ✴

is a very special unitary Hessenberg matrix Exploit the representation of

as a unitary Hessenberg plus a rank-one matrix Unitary Hessenberg matrices belong to the class of rank-structured matrices

Cortona, 2004 – p.5/19

slide-6
SLIDE 6

Unitary Hessenberg Matrices

  • 1. Let

be unitary Hessenberg and

✵ ✆ ✢ ✣

its

✢ ✣

factorization, where

has nonnegative diagonal entries

  • 2. We have
✣ ✆ ✧ ✝

and, therefore,

✵ ✆ ✢ ✶✸✷ ✹ ✺✻✺✼✺✽✺✻✺✼✺✽✺✻✺✼✺✽✾ ✿ ❀❂❁❄❃ ❁❄❅ ✿ ❀❂❁❄❃ ❆ ❅ ❁❄❇ ❈ ❈ ❈ ❈ ❈ ❈ ✿ ❀ ❁❄❃ ❆ ❅ ❉ ❉ ❉ ❆❋❊✗● ❅ ❁ ❊ ❀✱❆ ❅ ✿ ❀ ❁❍❅ ❁❍❇ ❈ ❈ ❈ ❈ ❈ ❈ ✿ ❀ ❁❍❅ ❆ ❇ ❉ ❉ ❉ ❆■❊✗● ❅ ❁ ❊ ❀✱❆ ❇

... . . . ... ... . . .

❀ ❆ ❊✗● ❅ ✿ ❀ ❁ ❊✗● ❅ ❁ ❊ ❏ ❑✻❑✼❑✽❑✻❑✼❑✽❑✻❑✼❑✽▲

[Gill, Golub, Murray, Saunders , 1974]

Cortona, 2004 – p.6/19

slide-7
SLIDE 7

The Rank Structure

▼ ✁ ✵ ☎ ✆ ◆ ❖ P ✖ ◗ ❘ ◗ ✝✙✘ ✖

rank

✁ ✵ ❙ ❚ ✪ ✕ ✩ ✭ ☛ ✕ ✩ ❚ ❯ ☎ ❱ ✁ ✵ ☎ ✆ ◆ ❖ P ✖ ◗ ❘ ◗ ✝✙✘ ✖

rank

✁ ✵ ❙ ✕ ✩ ❚ ☛ ❚ ✪ ✕ ✩ ✭ ❯ ☎ ▼ ✁ ✵ ☎❳❲ ✕

and

❱ ✁ ✵ ☎❳❲ ✕

for

unitary Hessenberg Depending on some small differences in the definition and in the representation, we may have:

  • 1. Quasiseparable matrices [Eidelman, Gohberg]
  • 2. Matrices with low Hankel rank [Dewilde, Van der Veen]
  • 3. Weakly semiseparable matrices [Tyrtyshnikov]
  • 4. Sequentially/Hierarchically semiseparable matrices

[Chandrasekaran, Gu]

Cortona, 2004 – p.7/19

slide-8
SLIDE 8

Companion Matrices Under QR

Theorem 1 Let

❨ ☞ ✤ ❩

be the sequence of matrices generated by the shifted

✢ ✣

algorithm applied to the companion matrix

☞ ✆ ☞ ✡

. Then

▼ ✁ ☞❬✤ ☎ ❲ ✕

and

❱ ✁ ☞❬✤ ☎ ❲ ❭

for all

Proof: From

☞ ✤ ✆ ✢ ✤ ☞ ✤ ✘ ✖ ✢ ❫ ✤

and

☞ ✡ ✆ ✴ ✡ ✪✱❴ ❵ ❛ ❜❞❝ ❵ ❛ ❜ ❫

, one gets

☞ ✤ ✆ ✴ ✤ ✪ ❴ ❵ ❡ ❜❞❝ ❵ ❡ ❜ ❫ ☛ ✴ ✤ ❢❣ ❤ ✐ ❖ ❥ ❦ ❧♥♠♦ ♦ ♠ ❣ ♣ ♠ ❥ q

Since

☞✥✤

is Hessenberg, it also holds

r s ✤ t ✞✈✉ ❘ ✆ ✓ ✇ s ✤ t ✞ ①③② s ✤ t ❘ ☛ ④⑥⑤ ❥ ⑦ ✓ ❚ ⑧ ⑨

To conclude: Characterize the rank structure of such a matrix

✴ ✤

Cortona, 2004 – p.8/19

slide-9
SLIDE 9

Unitary Rank-Structured Matrices

Theorem 2 Let

be a unitary matrix such that tril

✁ ✴ ☛ ✓ ⑨ ☎ ✆

tril

✁ ❴ ❝ ❫ ☛ ✓ ⑨ ☎

. Then,

❱ ✁ ✴ ☎❳❲ ⑨

. In particular, there exist

✭ ✓ ✕

vectors

⑩ ❘ ❶ ❷ ✮

,

✕ ❲ ❚ ❲ ✭ ✓ ✕

,

✭ ✓ ✕

vectors

❸ ❘ ❶ ❷ ✮

,

⑨ ❲ ❚ ❲ ✭

and

✭ ✓ ⑨

lower triangular matrices

❹ ❘ ❶ ❷ ✮❂❺ ✮

,

⑨ ❲ ❚ ❲ ✭ ✓ ✕

, such that

r ✞✈✉ ❘ ✆ ⑩ ✳ ✞ ❹ ✞ ★ ✖ ❻ ❻ ❻ ❹ ❘ ✘ ✖ ❸ ❘ ☛ ④ ⑤ ❥ ⑦ ✪ ✕ ❲ ❚ ❲ ✭

Proof: Compute the QR factorization of

. The unitary factor

is the product of

✬ ✁✄✭ ☎
  • rotations. The rank structure of

coincides with the rank structure of

Cortona, 2004 – p.9/19

slide-10
SLIDE 10

A Representation for

Theorem 3 For each matrix

☞❽✤

generated by shifted QR algorithm applied to the companion matrix

☞ ✆ ☞ ✡

, there exist

✭ ✓ ✕

vectors

⑩ ❵ ❡ ❜ ❾ ❶ ❷ ✮

,

✭ ✓ ✕

vectors

❸ ❵ ❡ ❜ ❾ ❶ ❷ ✮

and

✭ ✓ ⑨

lower triangular matrices

❹ s ✤ t ❘ ❶ ❷ ✮❂❺ ✮

such that

  • s
✤ t ✞ ✉ ❘ ✆ ⑩ ❵ ❡ ❜ ❿ ✳ ❹ s ✤ t ✞ ★ ✖ ❻ ❻ ❻ ❹ s ✤ t ❘ ✘ ✖ ❸ ❵ ❡ ❜ ❾ ✪ ✇ s ✤ t ✞ ①③② s ✤ t ❘ ☛ ④ ⑤ ❥ ⑦ ✪ ✕ ❲ ❚ ❲ ✭ ➀
  • s
✤ t ✞ ✉ ✞ ✆
  • s
✤ t ✞ ☛ ④⑥⑤ ❥ ✕ ❲ ⑦ ❲ ✭ ➀
  • s
✤ t ✞ ★ ✖ ✉ ✞ ✆ ➁ s ✤ t ✞ ☛ ④ ⑤ ❥ ✕ ❲ ⑦ ❲ ✭ ✓ ✕ ➀

Each iterate

☞✥✤

can be represented by means of

➂ ✕ ✕ ✭

parameters.

Cortona, 2004 – p.10/19

slide-11
SLIDE 11

The Structured QR Iteration

Let

➃ ✁ ❨ ⑩ ❵ ❡ ❜ ❿ ❩ ☛ ❨ ❸ ❵ ❡ ❜ ❿ ❩ ☛ ❨ ❹ s ✤ t ✞ ❩ ☛ ✰ ❵ ❡ ❜ ☛ ➄ ❵ ❡ ☛ ❴ ❵ ❡ ❜ ☛ ❝ ❵ ❡ ❜ ☎ ✆ ☞ ✤

The structured QR iteration takes in input

✁ ❨ ⑩ ❵ ❡ ❜ ❿ ❩ ☛ ❨ ❸ ❵ ❡ ❜ ❿ ❩ ☛ ❨ ❹ s ✤ t ✞ ❩ ☛ ✰ ❵ ❡ ❜ ☛ ➄ ❵ ❡ ❜ ☛ ❴ ❵ ❡ ❜ ☛ ❝ ❵ ❡ ❜ ☛ ➅ ➆ ➇ ➈ ➉ ➊ ➋➌➎➍ ✦ ✤ ☎

and returns in output

✁ ❨ ⑩ ❵ ❡ ➏ ➐ ❜ ❿ ❩ ☛ ❨ ❸ ❵ ❡ ➏ ➐ ❜ ❿ ❩ ☛ ❨ ❹ s ✤ ★ ✖ t ✞ ❩ ☛ ✰ ❵ ❡ ➏ ➐ ❜ ☛ ➄ ❵ ❡ ➏ ➐ ❜ ☛ ❴ ❵ ❡ ➏ ➐ ❜ ☛ ❝ ❵ ❡ ➏ ➐ ❜ ☎

such that

☞✥✤ ★ ✖

is equal to

➃ ✁ ❨ ⑩ ❵ ❡ ➏ ➐ ❜ ❿ ❩ ☛ ❨ ❸ ❵ ❡ ➏ ➐ ❜ ❿ ❩ ☛ ❨ ❹ s ✤ ★ ✖ t ✞ ❩ ☛ ✰ ❵ ❡ ➏ ➐ ❜ ☛ ➄ ❵ ❡ ➏ ➐ ❜ ☛ ❴ ❵ ❡ ➏ ➐ ❜ ☛ ❝ ❵ ❡ ➏ ➐ ❜ ☎

Cortona, 2004 – p.11/19

slide-12
SLIDE 12

A Fast and Stable Structured Iteration

INPUT:

✁ ❨ ⑩ ❵ ❡ ❜ ❿ ❩ ☛ ❨ ❸ ❵ ❡ ❜ ❿ ❩ ☛ ❨ ❹ s ✤ t ✞ ❩ ☛ ✰ ❵ ❡ ❜ ☛ ➄ ❵ ❡ ☛ ❴ ❵ ❡ ❜ ☛ ❝ ❵ ❡ ❜ ☛ ➅ ➆ ➇ ➈ ➉ ➊ ➋➌➎➍ ✦ ✤ ☎ ➃ ✁ ❨ ⑩ ❵ ❡ ❜ ❿ ❩ ☛ ❨ ❸ ❵ ❡ ❜ ❿ ❩ ☛ ❨ ❹ s ✤ t ✞ ❩ ☛ ✰ ❵ ❡ ❜ ☛ ➄ ❵ ❡ ❜ ☛ ❴ ❵ ❡ ❜ ☛ ❝ ❵ ❡ ❜ ☎ ✆ ☞ ✤ ➑ ✁ ❨ ⑩ ❵ ❡ ❜ ❿ ❩ ☛ ❨ ❸ ❵ ❡ ❜ ❿ ❩ ☛ ❨ ❹ s ✤ t ✞ ❩ ☛ ✰ ❵ ❡ ❜ ☛ ➄ s ✤ t ☛ ❴ ❵ ❡ ❜ ☛ ❝ ❵ ❡ ❜ ☎ ✆ ✵ ✤
  • 1. Compute
✢ ✤

such that

☞ ✤ ✓ ✦ ✤ ✧ ✝ ✆ ✢ ✤ ✣ ✤
  • 2. (Expansion Phase) Compute
✵ ✤ ★ ✖

: (a)

✵ ✤ ★ ✖ ✆ ✢ ❫ ✤ ✵ ✤ ✢ ✤

(b)

☞ ✤ ★ ✖ ✩ ✆ ✣ ✤ ✢ ✤ ✪ ✦ ✤ ✧ ✝

,

✵ ✤ ★ ✖ ✆ ☞✥✤ ★ ✖ ✓ ❴ ❵ ❡ ➏ ➐ ❜❞➒ ❵ ❡ ➏ ➐ ❜ ❫
  • 3. (Compression Phase) Compute a QR factorization of
✵ ✤ ★ ✖

OUTPUT:

✁ ❨ ⑩ ❵ ❡ ➏ ➐ ❜ ❿ ❩ ☛ ❨ ❸ ❵ ❡ ➏ ➐ ❜ ❿ ❩ ☛ ❨ ❹ s ✤ ★ ✖ t ✞ ❩ ☛ ✰ ❵ ❡ ➏ ➐ ❜ ☛ ➄ ❵ ❡ ➏ ➐ ❜ ☛ ❴ ❵ ❡ ➏ ➐ ❜ ☛ ❝ ❵ ❡ ➏ ➐ ❜ ☎

Cortona, 2004 – p.12/19

slide-13
SLIDE 13

A Representation for

❼ ➓

Theorem 4 At the end of the expansion phase the matrix

✵ ✤ ★ ✖

admits the following representation. There exist vectors

➔ ❵ ❡ ➏ ➐ ❜ ❿ ❶ ❷ ✯

, vectors

→ ❵ ❡ ➏ ➐ ❜ ❿ ❶ ❷ ✯

and matrices

➣ s ✤ ★ ✖ t ✞ ❶ ❷ ✯ ❺ ✯

such that

↔ s ✤ ★ ✖ t ✞ ✉ ❘ ✆ ➔ ❵ ❡ ➏ ➐ ❜ ❿ ✳ ➣ s ✤ ★ ✖ t ✞ ★ ✖ ❻ ❻ ❻ ➣ s ✤ ★ ✖ t ❘ ✘ ✖ → ❵ ❡ ➏ ➐ ❜ ❾ ☛ ④ ⑤ ❥ ❚ ✓ ⑦ ⑧ ✕ ↔ s ✤ ★ ✖ t ✞ ✉ ❘ ✆ ✓ ✇ s ✤ ★ ✖ t ✞ ① ② s ✤ ★ ✖ t ❘ ④⑥⑤ ❥ ⑦ ✓ ❚ ⑧ ⑨ ✒ ❱ ✁ ✵ ✤ ★ ✖ ☎ ❲ ❭ ✫

redundant representation QR factorization computed in linear time [Eidelman, Gohberg]

Cortona, 2004 – p.13/19

slide-14
SLIDE 14

Computational Analysis

Computational cost:

➂ ⑨ ✑ ✑ ✭
  • flops. Maybe can be

slightly reduced Linear storage: all the matrices are stored in linear data structures No divisions

absolute errors are small as possible. Perturbation results for eigenvalues depend on the magnitude of the absolute perturbation Unitary transformations

no significant growth of intermediate results. NUMERICALLY STABLE BEHAVIOR

Cortona, 2004 – p.14/19

slide-15
SLIDE 15

A Numerical Experiment

Companion matrix of size

✭ ✆ ↕ ➙

with random entries

20 40 60 80 100 120 140 10

−16

10

−15

10

−14

10

−13

Plot of || eig QR −eig QRS||

Cortona, 2004 – p.15/19

slide-16
SLIDE 16

Topics for Future Researches

  • 1. Incorporate scaling and balancing techniques
  • 2. Compare theoretically and experimentally with the fast

QR algorithms for unitary Hessenberg matrices designed by [Gragg & coworkers]

  • 3. Investigate the problem of reconstructing a

rank-structured matrix given a partial description and some additional requirements [Dym, Gohberg for invertible linear

  • perators]
  • 4. Backward stability analysis for the structured QR

algorithm [Tisseur for the standard QR algorithm and Stewart for the fast unitary

Hessenberg QR algorithms]

  • 5. Block companion matrices (Matrix Polynomials)

Cortona, 2004 – p.16/19

slide-17
SLIDE 17

Scaling and Balancing

Backward stability for companion matrices gives poor bounds for polynomials

  • 1. Balancing can not be used
  • 2. Consider the QZ algorithm applied to the companion

matrix pencil

☞ ✓ ➛ ❹

suitably normalized. From [Van Dooren & Dewilde]

backward stability: The computed roots are exact for

➝➟➞ ✁ ✂ ☎

with

➠ ➡♥➢ ✓ ➢ ➠ ✮ ❲ ❪ ✁✄✭ ☎ ➠➤➢ ➠ ✮ ➥

Experimentally: The rank structure is invariant under the QZ algorithm

Cortona, 2004 – p.17/19

slide-18
SLIDE 18

Software: The State of Art

Only MATLAB implementations are available Need Fortran and C++ routines to compare timing and accuracy with the LAPACK programs Software development is time consuming but I think This should be the time to create a software package concerning fast methods for rank-structured eigenvalue problems Collaborations are welcome!!!

Cortona, 2004 – p.18/19

slide-19
SLIDE 19

References

  • 1. D. A. Bini, L. Gemignani, V. Y. Pan.
➦ ➧
  • like algorithms for generalized

semiseparable matrices. Submitted to Numer. Math, 2003.

  • 2. D. Fasino, L. Gemignani, N. Mastronardi, M. Van Barel. Orthogonal

rational functions and structured matrices. To appear in SIAM J. Matrix Anal., 2003.

  • 3. D. A. Bini, F

. Daddi, L. Gemignani. The shifted QR algorithm for companion matrices. To appear in ETNA, 2004.

  • 4. D. A. Bini, L. Gemignani, V. Y. Pan. Computing polynomial roots in
➨ ➩ ➫ ➭

memory and

➨ ➩ ➫ ❇ ➭

time by a numerically reliable algorithm based on

➦ ➧
  • iteration. To appear in ETNA, 2004.
  • 5. L. Gemignani. A
➦ ➧
  • based method for solving the unitary

Hessenberg eigenvalue problem via semiseparable matrices. Submitted to Comput. Appl. Math.

Cortona, 2004 – p.19/19