Efficient Algorithms for Structured Matrices Vadim Olshevsky - - PowerPoint PPT Presentation

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Efficient Algorithms for Structured Matrices Vadim Olshevsky

www.math.uconn.edu/˜olshevsky University of Connecticut Joint works with T.Bella (U Conn), Yu.Eidelman (Tel Aviv), I.Gohberg (Tel Aviv), A.Shokrollahi (EPFL). MTNS, Kyoto, Japan, July 2006 This work was supported by the NSF grants CCR 0242518 .

– Typeset by FoilT EX – 1

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Contents

• Examples and introduction.
• Algorithms.

– Tangential Interpolation. – Displacement structure. Fast algorithms. – Superfast algorithms.

• Applications

– List Decoding.

• Quasiseparable matrices and eigenvalue problems

– Typeset by FoilT EX – 2

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Part 1. Examples of structured matrices

Toeplitz, T = ti−j

• Vandermonde, V =
• xj−1

i

• 

     t0 t−1 · · · · · · t−n+1 t1 t0 t−1 . . . . . . ... ... ... . . . . . . ... t0 t−1 tn−1 · · · · · · t1 t0             1 x1 x2

1

· · · xn−1

1

. . . . . . . . . . . . . . . . . . . . . . . . 1 xn x2

n

· · · xn−1

n

      Cauchy, C =

• 1

xi−yj

• Pick, P =

1−xix∗

j

zi+z∗

j

• 

 

1 x1−y1

· · ·

1 x1−yn

. . . . . .

1 xn−y1

· · ·

1 xn−yn

       

1−x1·x1∗ z1+z1∗

· · ·

1−x1·xn∗ z1+zn∗

. . . . . .

1−xn·x1∗ zn+z1∗

· · ·

1−xn·xn∗ zn+zn∗

    

• Exploiting the structure ⇛ Fast algorithms.

General matrices Standard algorithms Gaussian Elimination O(n3) flops Toeplitz Matrices Fast algorithms Levinson and Schur O(n2) Toeplitz Matrices Superfast algorithms Brent-Gustafson-Yan O(n log2 n)

– Typeset by FoilT EX – 3

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Growing interest

• Conferences:

Santa Barbara-1996, Boulder-1999 (AMS), Cortona-2000, South Hadley-2001 (AMS), Moscow-2003, Cortona-2004.

• Special sessions and minisymposia at SIAM, ILAS, MTNS, IWOTA (2002, 2003,

2004, 2005, 2006) and SPIE conferences.

• Recent publishing projects:

– 1999. Kailath & Sayed, Eds, SIAM Publicatiosn. – 2000. Bini, Tyrtyshnikov and Yalamov, Eds, Nova Publications. – 2001. Two-Volume set (38 papers) edited V.O. to be published by AMS. – 2002. A special issue of LAA, Dewilde, Sayed & V.O. Eds. – 2003. Fast algorithms for structured matrices, V.O., ed., AMS and SIAM. —————————————–

• 2004-2006. Semiseparable matrices. Cortona-2004, FOCM-2005, Householder-2005,

IWOTA-2005, Moscow-2005, ILAS-2006.

– Typeset by FoilT EX – 4

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Examples of efficient algorithms

• Full research cycle.

– Extraction. Applications. Mathematical models. ∗ Estimation, Signal processing, Filtering, Circuit simulations, Systems and Control, Image Processing, Error-correcting codes. – Processing. Fundamental Mathematics. ∗ Various interpolation and Approximation problems, Orthogonal polynomials, Moment problems, RKHS, Commutant-lifting theory, etc. – Injection. Algorithm development. Software design.

• A

RECENT DISCOVERY: Though fast algorithms were believed to be typically inaccurate, it turns out that there are better modified algorithms that BLEND SPEED AND ACCURACY.

– Typeset by FoilT EX – 5

Efficient Algorithms for Structured Matrices Vadim Olshevsky

• Example. Positive definite Hankel matrices (Injection)

H =       h0 h1 h2 · · · hn−1 h1 h2 ... . . . h2 ... h2n−3 . . . ... h2n−3 h2n−21 hn−1 · · · h2n−3 h2n−2 h2n−1       > 0 Conditioning

• Tyrtyshnikov (1994):

k2(H) = H2H−12 > 2n−6

– Typeset by FoilT EX – 6

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Backward Stable Fast Algorithms (Injection)

• Slow and accurate: Gaussian elimination:

– Exact arithmetic: H = LLT – Floating point arithmetic: H ≈ L LT – Cost: O(n3). – Stability: H − L LT2 < O(ǫ)H2 .

• Fast: (There are many...)

– Cost: O(n2) or O(n log2 n). – Stability: ? .

• Fast and Accurate [OS01a]:

– Cost: O(n2). – Stability: H − L LT2 < O(ǫ)H2 .

– Typeset by FoilT EX – 7

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Part 2. The confluent tangential Nevanlinna-Pick problem

The classical scalar Nevanlinna-Pick problem

• Given:

points z1, . . . , zn in the RHP: Rezk ≥ 0 values f1, . . . , fn in the unit disk: |fk| < 1

• Construct an interpolant

f(zk) = fk , such that (i) f(z) is

passive |f(z)| < 1 in the RHP

(ii) f(z) is rational The Nevanlinna algorithm (1929) Its matrix interpretation computation of the interpolant f(z) The Cholesky factorization R = LL∗, [ L is lower triangular] for the Pick matrix R = 1−fif∗

j

zi+z∗

j

• O(n2)

O(n3)

We shall show how to do faster than O(n2)

– Typeset by FoilT EX – 8

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Applications to system and circuit theory

• Many physical phenomena can be described by a

Linear Time-Invariant system

Input

✲

LTI system

✲ Output

u(z)

✲

f(z)

✲ y(z) = f(z) · u(z)

Passivity

∆

= |f(x)| ≤ 1

∆

= Conservation of energy

• Energy of input

u(z) ≥ Energy of output y(z).

• Many engineering applications have been discovered, e.g., the classical Nevanlinna

(1929) algorithm is known under the name “Darlington (1939) synthesis procedure”

– Typeset by FoilT EX – 9

Efficient Algorithms for Structured Matrices Vadim Olshevsky

MIMO systems. Tangential interpolation

• MIMO systems (multi-input-multi-output) give rise to matrix rational functions

Θ(z) = pij(z)

qij(z)

• .

x1 x2 xN

✲ ✲

. . .

✲

Θ(z) y1 y2 yM

✲ ✲

. . .

✲

• This gives rise to the tangential interpolation conditions:
• x1

· · · xN

• · Θ(zk)

=

• y1

· · · yM

• – Typeset by FoilT

EX – 10

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Confluent interpolation conditions

• Given:

Now for each interpolation point zk we are given two chains of vectors: {xk1, . . . , xk,mk} {yk1, . . . , yk,mk}

• We consider confluent interpolation conditions:

xk1 . . . xk,mk

• ×

       F(zk) F ′(zk) · · · · · ·

F (mk−1)(zk) (mk−1)!

F(zk) ... . . . ... ... . . . . . . ... ... F ′(zk) · · · · · · F(zk)        = yk,1 · · · yk,mk

• .

– Typeset by FoilT EX – 11

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Example: linear pencil and Jordan chains

• Let

F(z) = A − zI, then F ′(z) = −I, F ′′(z) = 0.

• If

uk1 uk2 uk3

• ·

  (A − zkI) −I (A − zkI) −I (A − zkI)   = .

• Hence uk1 is the (left)

eigenvector : uk1(A − zkI) = 0.

• Further, uk2 is the first

generalized eigenvector uk2(A − zkI) = uk1.

• Finally, uk3 is the second

generalized eigenvector : uk3(A − zkI) = uk2.

– Typeset by FoilT EX – 12

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Computational cost of our interpolation problem

• Given:

{zk}r

k=1

are given r points {mk}r

k=1

are their multiplicities {xk1, . . . , xk,mk}r

k=1

are the associated chains {yk1, . . . , yk,mk}r

k=1

• The cost C(n) is [OS01b] :

n log2 n Caratheodory-Fejer ≤ C(n) = O(n log2 n · [1 + r

k=1 mk n log n mk])

≤ n log3 n Nevanlinna-Pick arithmetic operations, where n = m1 + . . . + mr HOW? VIA STRUCTURED MATRICES. We present next a matrix description of the algorithm

– Typeset by FoilT EX – 13

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Part 3. Fast multiplication algorithms for structured matrices. Toeplitz matrices. Discrete Convolution.

• Matrix-vector product:

      a0 a1 a0 a2 a1 a0 . . . ... ... ... · · · a2 a1 a0             b0 b1 b2 . . .       =       a0b0 a1b0 + a0b1 a2b0 + a1b1 + a0b3 . . .      

• Polynomials:

(a0+a1z+a2z2+. . .)(b0+b1z+b2z2+. . .) = a0b0+(a1b0+a0b1)z+(a2b0+a1b1+a0b3)z2+. . .

– Typeset by FoilT EX – 14

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Vandermonde matrices. Multipoint polynomial evaluation.

• Matrix-vector product:

    1 x1 x2

1

· · · xn−1

1

1 x2 x2

2

· · · xn−1

2

. . . . . . . . . . . . 1 xn x2

n

· · · xn−1

n

        a0 a1 . . . an−1     =     f(x1) f(x2) . . . f(xn)    

• Polynomial multipoint evaluation of

f(x) = a0 + a1x + a2x2 + . . . + an−1xn−1 at x1, x2, . . . , xn.

• Matrix-vector product for the inverse of Vandermonde matrix: interpolation.
• Discrete Fourier transforms : a special case. xk = exp(2πik

n )

– Typeset by FoilT EX – 15

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Cauchy matrices. Multipoint rational evaluation.

• Matrix-vector product:

  

1 x1−y1

· · ·

1 x1−yn

. . . . . .

1 xn−y1

· · ·

1 xn−yn

     a0 . . . an−1   =   f(x1) . . . f(xn)  

• Rational multipoint evaluation of

f(x) = a0 1 x − y1 + . . . + an−1 1 x − yn . at x1, . . . , xn

• Matrix-vector product for the inverse of Cauchy matrix: rational interpolation.

– Typeset by FoilT EX – 16

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Cost for structured matrix-vector multiplication

DFT matrices (i.e., Vandermonde matrices with special nodes) Fast Fourier transform O(n log n) Toeplitz/Hankel/Bezoutian matrices convolution O(n log n) Vandermonde matrices multipoint polynomial evaluation O(n log2 n) inverse Vandermonde matrices polynomial interpolation O(n log2 n) Cauchy matrices multipoint rational evaluation O(n log2 n) inverse Cauchy matrices rational interpolation O(n log2 n)

– Typeset by FoilT EX – 17

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Displacement rank (needed to solve tangential Nevanlinna-Pick)

• Let {F, A} be fixed matrices.

α = Displacement rank of R = rank (FR - RA)

• Recall basic classed of matrices with low displacement rank:

Toeplitz-like rank (ZR − RZ)<< n Vandermonde-like rank (D−1

x R − RZT)<< n

Cauchy-like rank (DxR − RDy)<< n ZT = Jn(0), Dx = diag(x).

– Typeset by FoilT EX – 18

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Why the names? A few examples

• The displacement rank of a Cauchy matrix C =
• 1

xi−yj

• is equal to 1:

DxC − CDy = xi−yj

xi−yj

• =

    1 1 · · · 1 1 1 · · · 1 . . . . . . . . . 1 1 · · · 1    

• If the displ. rank (R) << n then R is called Cauchy-like.
• The Pick matrices P =

1−fif∗

j

zi+z∗

j

• are Cauchy-like (displ. rank = 2):

DzP + PD∗

z = [1 − fif ∗ j ] =

    1 −f1 1 −f2 . . . . . . 1 −fn    

• 1

1 · · · 1 f ∗

1

f ∗

2

· · · f ∗

n

• .
• For Toeplitz matrices the concept of displacement was introduced by Kailath et al

(1978).

Fast matrix-vector multiplication?

– Typeset by FoilT EX – 19

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Cost for structured matrix-vector multiplication

Toeplitz-like matrices O(n log n) Vandermonde-like matrices O(n log2 n) Cauchy-like matrices O(n log2 n) Toeplitz-like rank (ZR − RZ)<< n Vandermonde-like rank (D−1

x R − RZT)<< n

Cauchy-like rank (DxR − RDy)<< n ZT = Jn(0), Dx = diag(x).

– Typeset by FoilT EX – 20

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Matrices in the confluent tangential Nevanlinna-Pick problem?

• A superfast algorithm for R = L · L∗ for what we called

confluent Cauchy-like rank (AζR − RAπ)<< n where matrices Aπ, Aζ are assumed to be in the Jordan form.

• Hence confluent Cauchy-like matrices include

Toeplitz-like convolution Vandermonde-like FFT, polynomial interpolation, multipoint evaluation Cauchy-like rational interpolation rational multipoint evaluation

– Typeset by FoilT EX – 21

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Confluent Cauchy matrices.

• Let {

n

• x1, . . . , x1
• m1

, . . . , xs, . . . , xs

• ms

}{

n

• y1, . . . , y1
• k1

, . . . , yt, . . . , yt

• kt

}, and B(x, y) = b(x) − b(y) x − y , where b(x) = (x − y1)k1 · (x − y2)k2 · . . . · (x − ys)kt

• We define the block Confluent Cauchy matrix matrix

C =

• Ci,j
• 1≤i≤s,1≤j≤t ,

where each mi × kj block Cij has the form (next slide):

– Typeset by FoilT EX – 22

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Confluent Cauchy matrices. Continuation

• Cij =

        B(xi, yj) ∂1

yB(xi, yj)

· · · ∂

kj−1 y

B(xi, yj) ∂1

xB(xi, yj)

∂1

x∂1 yB(xi, yj)

· · · ∂1

x∂ kj−1 y

B(xi, yj) ∂2

xB(xi, yj)

∂2

x∂1 yB(xi, yj)

· · · ∂2

x∂ kj−1 y

B(xi, yj) . . . . . . . . . ∂mi−1

x

B(xi, yj) ∂mi−1

x

∂1

yB(xi, yj)

· · · ∂mi−1

x

∂

kj−1 y

B(xi, yj)         ,

• Fast matrix-vector product?
• Yes. Next slide.

– Typeset by FoilT EX – 23

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Cost of matrix-vector product [OS00]

• Thm. Confluent Cauchy-like matrices can be multiplied by vectors using

C(n) = O(n log n · [1 +

s

• i=1

mi n log n mi +

t

• j=1

kj n log n kj ]) flops, where – {m1, . . . , ms} are the sizes of Jordan blocks of Aπ – {k1, . . . , kt} are the sizes of Jordan blocks of Aζ

• C(n) can vary from O(n log n) to O(n log2 n) depending upon the Jordan structure.

– Typeset by FoilT EX – 24

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Part 4. Codes

Codes are used when messages are to be transmitted over noisy channels.

E S S A G M E S S A G M E A E S S A G M E E S S A G M E E S S A G M E E S S A G M E E S S A G M E S S A G M E A With coding Without coding

– Typeset by FoilT EX – 25

Efficient Algorithms for Structured Matrices Vadim Olshevsky

One idea for encoding: linear codes

MESSAGE CODEWORD A A A A − → encoding − → A A A A B B ⋔ ⋔ F4

q

F6

q

space subspace 4 is called the dimension of the code. 6 is called block-length of the code.

– Typeset by FoilT EX – 26

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Reed-Solomon Codes

• Given: x1, . . . , xn distinct elements of Fq..
• Encode

{a1, . . . , ak} →

• f(x1), . . . , f(xn)
• where

k ≤ n where f(x) = a1 + a2x + . . . + akxk−1,

– Typeset by FoilT EX – 27

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Welch-Berlekamp Method

Let (y1, . . . , yn) be close to the codeword

• f(x1), . . . , f(xn)
• .

Find g and h such that ∀ i = 1, . . . , n: h(xi)yi + g(xi) = 0. Then f(X) = −g(X) h(X).

vanishes at these points h sent word received word f Error locator polynomial

– Typeset by FoilT EX – 28

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Example

n = 6 and k = 4:            1 x1 x2

1

x3

1

x4

1

y1 y1x1 y1x2

1

1 x2 x2

2

x3

2

x4

2

y2 y2x2 y2x2

2

1 x3 x2

3

x3

3

x4

3

y3 y3x3 y3x2

3

1 x4 x2

4

x3

4

x4

4

y4 y4x4 y4x2

4

1 x5 x2

5

x3

5

x4

5

y5 y5x5 y5x2

5

1 x6 x2

6

x3

6

x4

6

y6 y6x6 y6x2

6

                           g0 g1 g2 g3 g4 h0 h1 h2                 = 0.

– Typeset by FoilT EX – 29

Efficient Algorithms for Structured Matrices Vadim Olshevsky

List-Decoding of Reed-Solomon Codes (Madhu Sudan)

What if the number of errors is larger than (n − k)/2? Compute polynomials h0, . . . , hℓ such that ∀ i = 1, . . . , n: h0(xi) + h1(xi)yi + · · · + hℓ(xi)yℓ

i = 0.

– Typeset by FoilT EX – 30

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Bivariate Interpolation

Want: Nontrivial Polynomial H(x, y) := ℓ

i=0 hi(x)yi

H := h00 + h01x + h02x2 + h10y + h11xy + h20y2. We have       1 x1 x2

1

y1 y1x1 y2

1

1 x2 x2

2

y2 y2x2 y2

2

1 x3 x2

3

y3 y3x3 y2

3

1 x4 x2

4

y4 y4x4 y2

4

1 x5 x2

5

y5 y5x5 y2

5

              h00 h01 h02 h10 h11 h20         = 0.

– Typeset by FoilT EX – 31

Efficient Algorithms for Structured Matrices Vadim Olshevsky

List decoding via displacement [OS03]

B B B B B B B B B B B B B B @ 1 x1 1 x2 1 x3 1 x4 1 x5 1 C C C C C C C C C C C C C C A B B B B B B @ 1 x1 x2 1 y1 y1x1 y2 1 1 x2 x2 2 y2 y2x2 y2 2 1 x3 x2 3 y3 y3x3 y2 3 1 x4 x2 4 y4 y4x4 y2 4 1 x5 x2 5 y5 y5x5 y2 5 1 C C C C C C A − B B B B B B @ 1 x1 x2 1 y1 y1x1 y2 1 1 x2 x2 2 y2 y2x2 y2 2 1 x3 x2 3 y3 y3x3 y2 3 1 x4 x2 4 y4 y4x4 y2 4 1 x5 x2 5 y5 y5x5 y2 5 1 C C C C C C A B B B B B B @ 1 1 1 1 1 1 C C C C C C A = B B B B B B @ 1/x1 y1/x1 − x2 1 y2 1/1 − y1x1 1/x2 y2/x2 − x2 2 y2 2/2 − y2x2 1/x3 y3/x3 − x2 3 y2 3/3 − y3x3 1/x4 y4/x4 − x2 4 y2 4/4 − y4x4 1/x5 y5/x5 − x2 5 y2 5/5 − y5x5 1 C C C C C C A @ 1 1 1 1 A .

Displacement rank = number of the codewords in the list +1.

– Typeset by FoilT EX – 32

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Tridiagonal, unitary Hessenberg, semi-separable and quasi-separable matrices

Quasi-separable

✬ ✫ ✩ ✪ ✬ ✫ ✩ ✪

Tridiagonal

✬ ✫ ✩ ✪

Unitary Hessenberg

✬ ✫ ✩ ✪

Semi- separable

• ✠

❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❘

– Typeset by FoilT EX – 33

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Two Definitions

• R is called order (rL, rU) quasiseparable if for some small rL, rU we have

rL = max rankR21, rU = max rankR12, where the maximum is taken over all symmetric partitions of the form R =

• ∗

R12 R21 ∗

• .
• R is called order (rL, rU)-semiseparable if for some small rL, rU we have

R = D + tril(RL) + triu(RU), where rankRL = rL, rankRU = rU with some RL, RU.

– Typeset by FoilT EX – 34

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Two examples

• Tridiagonal matrices:

R =           

β1 α1 γ2 α2

· · ·

1 α1 β2 α2 γ3 α3

... . . .

1 α2 β3 α3

... . . . . . .

1 α3

...

γn−1 αn−1

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

βn−1 αn−1 γn αn

· · ·

1 αn−1 βn αn

           (1)

• Unitary Hessenberg matrices:

      

−ρ1ρ∗ −ρ2µ1ρ∗ −ρ3µ2µ1ρ∗ · · · −ρn−1µn−2...µ1ρ∗ −ρnµn−1...µ1ρ∗ µ1 −ρ2ρ∗

1

−ρ3µ2ρ∗

1

· · · −ρn−1µn−2...µ2ρ∗

1

−ρnµn−1...µ2ρ∗

1

µ2 −ρ3ρ∗

2

· · · −ρn−1µn−2...µ3ρ∗

2

−ρnµn−1...µ3ρ∗

2

. . . ... µ3 . . . . . . . . . ... ... −ρn−1ρ∗

n−2

−ρnµn−1ρ∗

n−2

· · · · · · µn−1 −ρnρ∗

n−1

       (2)

– Typeset by FoilT EX – 35

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Two examples

• Tridiagonal and UH matrices are NOT SEMI-SEPARABLE.
• In fact, they are QUASI-SEPARABLE:

  · · · · · ·

γ4 α4

· · ·  

• R12 for (7)

,   −ρ4µ3µ2µ1ρ∗ · · · −ρn−1µn−2...µ1ρ∗ −ρnµn−1...µ1ρ∗ −ρ4µ3µ2ρ∗

1

· · · −ρn−1µn−2...µ2ρ∗

1

−ρnµn−1...µ2ρ∗

1

−ρ4µ3ρ∗

2

· · · −ρn−1µn−2...µ3ρ∗

2

−ρnµn−1...µ3ρ∗

2

 

• R12 for (9)
• They are order-one quasiseparable.
• More examples: companion matrices, arrow matrices, comrade matrices.

– Typeset by FoilT EX – 36

Efficient Algorithms for Structured Matrices Vadim Olshevsky

First references

• I. Gohberg, M.A. Kaashoek, L. Lerer, In Time-Variant Systems and Interpolation,

OT, 1992 (for infinite case)

• in the papers by A.-J.van der Veen starting from 1993 and in the monograph (see

also the complete list of references)

• P. Dewilde, A.-J. Van Der Veen, Time-varying Systems and Computations, Kluwer,

1998 matrices with small Hankel rank;

• in Eidelman and Gohberg, IEOT 27 (1997) as inverses to diagonal plus semiseparable

matrices. quasi-separable matrices

• Tyrtyshnikov , Mosaic ranks for weakly semiseparable matrices. (2000).

– Typeset by FoilT EX – 37

Efficient Algorithms for Structured Matrices Vadim Olshevsky

• Storage. O(N) parameters
• An order (rL, rU) quasiseparable matrix R can be represented by only

O(N) = O(rL + rU) parameters via R =

❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅

d1 ... ... dN pia×

ijqj

gib×

ijhj

(3) a×

ij = ai−1 · . . . · aj+1,

b×

ij = bi+1 · . . . · bj−1 with ai,i+1 = bi,i+1 = 1.

• The vectors or matrices {pi, qj, ak, gi, hj, bk, dk} are called generators of R.
• rL
• pi

aij× qj        rL

rU

• gi

bij× hj        rU

– Typeset by FoilT EX – 38

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Part I. Order-one quasiseparable matrices

– Typeset by FoilT EX – 39

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Divide-and-conquer eigenvalue algorithm

Quasi-separable

✬ ✫ ✩ ✪ ✬ ✫ ✩ ✪

Tridiagonal Gu,Eisenstat Bini,Pan

✬ ✫ ✩ ✪

Unitary Gragg, Reichel Ammar Sorensen Hessenberg

✬ ✫ ✩ ✪

Semi-

Chandrasecaran,

Gu

Mast., Van Camp

Van Barel separable

• ✠

❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❘

– Typeset by FoilT EX – 40

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Part 3. Divide-and-conquer algorithms

✫✪ ✬✩

N

• ✠

❅ ❅ ❅ ❅ ❘ ✫✪ ✬✩

N 2

• ✠

❅ ❅ ❅ ❘ ✫✪ ✬✩

N 2

• ✠

❅ ❅ ❅ ❅ ❘ ✫✪ ✬✩

N 4

... ...

✫✪ ✬✩

N 4

... ...

✫✪ ✬✩

N 4

... ...

✫✪ ✬✩

N 4

... ...

✫✪ ✬✩

1

✫✪ ✬✩

1

✫✪ ✬✩

1

✫✪ ✬✩

1

✫✪ ✬✩

1

✫✪ ✬✩

1

✫✪ ✬✩

1

✫✪ ✬✩

1 level 1 2 3 . . . log N Cost C 2 × C

2

4 × C

4

. . . N × C

N

Total Cost: C log N

– Typeset by FoilT EX – 41

Efficient Algorithms for Structured Matrices Vadim Olshevsky

D&C eigenvalue algorithm for order-one quasi-separable [EGKO06]

• An order-one quasiseparable R can be partitioned:

R =

• A′

m

B′

m+1

• rder-one quasiseparable

+

• Gm

Pm+1 Qm Hm+1

• rank-one

,

• The second matrix has rank = 1:

Pm+1 = col(pka×

km)N k=m+1, Qm = row(a× m,kqk)m k=1,

Gm = col(gkb×

km)m k=1, Hm+1 = row(b× mkhk)N k=m+1.

• The first matrix is order-one quasiseparable, e.g.,

A′

m(i, j) = (pi − gib× ima× m,i−1)a× ijqj,

1 ≤ j < i ≤ m. A′

m(i, i) = di − gib× ima× miqi,

i = 1, . . . , m A′

m(i, j) = gib× ij(hj − b× j−1,ma× mjqj),

1 ≤ i < j ≤ m.

– Typeset by FoilT EX – 42

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Generalized Sturm theorem [EGO05a].

• THM Let R be a Hermitian order-one-quasiseparable matrix. Let

Q0(λ) ≡ 1, Qk(λ) = det(R(1 : k, 1 : k) − λI) Let ν(λ) be the number of sign changes in the sequence QN(λ), QN−1(λ), . . . , Q1(λ), Q0(λ) (4) and let (α, β) be such that QN(α) = 0, QN(β) = 0. Then: – the difference ν(β) − ν(α) equals the number of eigenvalues of the matrix R in the interval (α, β) counting with its multiplicities.

How to evaluate (4)?

• In case {Qk} are orthogonal polynomials or the Szego polynomials the answer is

known.

– Typeset by FoilT EX – 43

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Example 1. Tridiagonal matrices and real orthogonal polynomials

• Let {Qk(x)} be real orthogonal polynomials satisfying three-term recurrence

relations: Qk(x) = (αk · x − βk) · Qk−1(x) − γk · Qk−2(x), (5)

• The relations (5) translate into the matrix form

Qk(x) = (α0 · . . . · αk) · det(xI − Rk×k) (1 ≤ k ≤ N) (6) where R =           

β1 α1 γ2 α2

· · ·

1 α1 β2 α2 γ3 α3

... . . .

1 α2 β3 α3

... . . . . . .

1 α3

...

γn−1 αn−1

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

βn−1 αn−1 γn αn

· · ·

1 αn−1 βn αn

           (7)

– Typeset by FoilT EX – 44

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Example 2. UH matrices and the Szego polynomials

• Let {Qk(x)} be the Szego polynomials satisfying two-term recurrence relations
• Gk+1(x)

Qk+1(x)

• =

1 µk+1

• 1

−ρ∗

k+1

−ρk+1 1 1 x Gk(x) Qk(x)

• .

(8)

• The relations (8) translate into the matrix form

Qk(x) = det(xI − Rk×k) µ0 · . . . · µk (1 ≤ k ≤ N) where R =       

−ρ1ρ∗ −ρ2µ1ρ∗ −ρ3µ2µ1ρ∗ · · · −ρn−1µn−2...µ1ρ∗ −ρnµn−1...µ1ρ∗ µ1 −ρ2ρ∗

1

−ρ3µ2ρ∗

1

· · · −ρn−1µn−2...µ2ρ∗

1

−ρnµn−1...µ2ρ∗

1

µ2 −ρ3ρ∗

2

· · · −ρn−1µn−2...µ3ρ∗

2

−ρnµn−1...µ3ρ∗

2

. . . ... µ3 . . . . . . . . . ... ... −ρn−1ρ∗

n−2

−ρnµn−1ρ∗

n−2

· · · · · · µn−1 −ρnρ∗

n−1

       (9)

– Typeset by FoilT EX – 45

Efficient Algorithms for Structured Matrices Vadim Olshevsky

New (quasi-separable?) polynomials

Quasi-separable New (quasi-separable?) polynomials

✬ ✫ ✩ ✪ ✬ ✫ ✩ ✪

Tridiagonal Real

• rthogonal

polynomials

✬ ✫ ✩ ✪

Unitary The Szego polynomials Hessenberg

✬ ✫ ✩ ✪

Semi- ??? ??? separable

• ✠

❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❘

• These new class of polynomials includes real orthogonal and the Szego polynomials

as special cases.

– Typeset by FoilT EX – 46

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Main results

• Three-term and two-term rr for the characteristic polynomials of submatrices of

general order-one quasi-separable.

• Eigenstructure analysis, formulas for the eigenvectors. Simple and multiple eigenvalue

cases are considered.

– Typeset by FoilT EX – 47

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Three-term recurrence relations [EGO05a]

• THM Let R be order-one quasiseparable with pihi = 0, and

Qk(λ) = det(R(1 : k, 1 : k) − λI). Then the Qk(λ) satisfy the three-term recurrence relations: Q1(λ) = d1 − λ, Q2(λ) = (d2 − λ)(d1 − λ) − p2q1g1h2; Qk(λ) = ψk(λ)Qk−1(λ) − ϕk(λ)Qk−2(λ), k = 3, . . . , N, where ψk(λ) = dk − λ + pkhk(ck−1 − ak−1bk−1λ) pk−1hk−1 , ϕk(λ) = pkhk pk−1hk−1 [(λ − dk−1)ak−1 + pk−1qk−1]·[(λ − dk−1)bk−1 + gk−1hk−1].

– Typeset by FoilT EX – 48

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Two-term recurrence relations [EGO05a]

• THM
• Q♯

1(λ)

Q1(λ)

• =
• d2λ − p2q1g1h2

λ − d1

• ;
• Q♯

k(λ)

Qk(λ)

• =
• ˜

ck(λ) − dk+1 λdk+1 + ˜ dk(λ) −1 λ Q♯

k−1(λ)

Qk−1(λ)

• ,

where ˜ ck(λ) = pk+1hk+1ck(λ) pkhk , ˜ dk(λ) = pk+1hk+1qkgk − ˜ ck(λ)dk; and QN(λ) = λQN−1(λ) − Q♯

N−1(λ).

– Typeset by FoilT EX – 49

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Part II Arbitrary Order quasiseparable matrices QR algorithm for Arbitrary Order quasiseparable matrices Application to the General eigenvalue problem

– Typeset by FoilT EX – 50

Efficient Algorithms for Structured Matrices Vadim Olshevsky

The results on QR factorization (not QR iteration) of arbitrary rank quasiseparable

• Infinite dimensional frameworks:

– P. Dewilde, A.-J. Van Der Veen, Time-varying Systems and Computations, Kluwer, 1998 QR = Inner – Outer – P. Dewilde, Matrices of small Hankel rank: A survey. In in Fast Algorithms for Structured Matrices: Theory and Applications, CONM/323, p. AMS publications, 2003.

• Y. Eidelman and I. Gohberg, A modification of the Dewilde-van der Veen method

for inversion of finite structured matrices. Linear Algebra and Application 343-344: 419-450 (2002)

What about QR iteration?

– Typeset by FoilT EX – 51

Efficient Algorithms for Structured Matrices Vadim Olshevsky

The quasiseparable structure of RQ

• The quasiseparable structure of S(1) = QR is inherited by S(2) = RQ:

S(1)

• (r,r)

= Q

• (r,r)

· R.

• (0,2r)

, S(2)

• (r,r)

= R

• (0,2r)

· Q.

• (r,r)
• Indeed,

S(1) =

• S(1)

11

S(1)

12

S(1)

21

S(1)

22

• =
• Q11

Q12 Q21 Q22

• ·
• R11

R12 R22

• S(2) =
• S(2)

11

S(2)

12

S(2)

21

S(2)

22

• =
• R11

R12 R22

• ·
• Q11

Q12 Q21 Q22

• ,

S(1)

21 = Q21R11,

S(2)

21 = R22Q21.

• The iteration was not considered in Dewilde-Van der Vein (1993) and Eidelman-

Gohberg (2002).

– Typeset by FoilT EX – 52

Efficient Algorithms for Structured Matrices Vadim Olshevsky

QR iteration

• Several papers/technical reports by Marc Van Barel with coauthors (E. Van Camp,

N.Mastronardi, R.Vanderbil, S. Van Huffel, D.Fasino....) – A QR algorithm for rank-one semiseparable. – Implicit shift technique. – Application to the general eigenvalue problem: ∗ Symmetric − → rank-one semiseparable. ∗ Symmetric − → tridiagonal − → rank-one semiseparable.

• D.Bini, L.Gemignani, V.Pan, 2004.

– A QR algorithm for rank-two quasiseparable with rank-one modification. – Application to root finding

A QR algorithm for symmetric arbitrary order quasiseparable matrices?

– Typeset by FoilT EX – 53

Efficient Algorithms for Structured Matrices Vadim Olshevsky

A QR algorithm for symmetric arbitrary order quasiseparable matrices [EGO05b]

• Starting point:

A(1)

• (r,r)

= V

• (r,∗)

· U

• (∗,r)

· R

• (0,2r)
• Computing Q:

A(1)

• (r,r)

= Q

• (r,r)

· R.

• (0,2r)
• Computing RQ:

A(2)

• (r,r)

= R.

• (0,2r)

· Q

• (r,r)

– Typeset by FoilT EX – 54

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Applications to the general eigenvalue problem

• Classical.

1. Reduce A − → T = QAQ∗ to a tridiagonal T. 2. Run tridiagonal QR iteration on T.

• Quasiseparable.

1. Reduce A − → S = QAQ∗ to a quasiseparable S. 2. Run quasiseparable QR iteration on S.

– Typeset by FoilT EX – 55

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Reduction

• First

step

• I2

f P1

• 

 

a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 a41 a42 a43 e A

  

• I2

f P1

• =

  

∗ ∗ ∗ hT ∗ ∗ ∗ hT ∗ ∗ ∗ ∗ h h ∗ ∗

   Here P1 = Householder reflection s.t. P1(

• a31

a41

• −
• a32

a42

• ) =
• ∗

· · · T .

• Second step.

  

∗ ∗ ∗ hT ∗ ∗ ∗ hT ∗ ∗ ∗ ∗ h h ∗ ∗

   − →    

∗ ∗ ∗ hTf P2

T

∗ ∗ ∗ hTf P2

T

∗ ∗ ∗ ∗ f P2h f P2h ∗ ∗

   

• Output:
• rder (2, 2) quasiseparable matrix.

Freedom: Preceed with the zero step:

• P0. The classical Householder reduction to tridiagonal is a a special case

corresponding to P0 =

• 1

f P0

• ,

such that

• P0 [ a21

a31

41 ]T = [ ∗

· · · 0 ]T .

• Conclusion: This is a family of quasiseparable algorithms INCLUDING the classical

tridiagonal as a special case. Which one is the best?

– Typeset by FoilT EX – 56

Efficient Algorithms for Structured Matrices Vadim Olshevsky

Tridiagonal, unitary Hessenberg, semi-separable and quasi-separable matrices

Quasi-separable

✬ ✫ ✩ ✪ ✬ ✫ ✩ ✪

Tridiagonal

✬ ✫ ✩ ✪

Unitary Hessenberg

✬ ✫ ✩ ✪

Semi- separable

• ✠

❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❘

– Typeset by FoilT EX – 57

Efficient Algorithms for Structured Matrices Vadim Olshevsky

References

[BEGKO07 ] T.Bella, Y.Eidelman, I.Gohberg, I.Koltracht, V.Olshevsky, A Bjorck- Pereyra-type algorithm for Szego-Vandermonde matrices based on properties of unitary Hessenberg matrices, to appear in LAA, 2007. [EGKO06a ] Yu.Eidelman, I.Gohberg, S. Kaur and V. Olshevsky A divide-and-conquer eigenvalue algorithm for quasiseparable matrices, technical report. [BEGKO06b ] T.Bella, Y.Eidelman, I.Gohberg, I.Koltracht, V.Olshevsky, A fast Bjorck-Pereyra like algorithm for solving quasiseparable-Hessenberg-Vandermonde systems, technical report. [BO06c ] T.Bella, V.Olshevsky, Fast inversion

• f

quasiseparable-Hessenberg- Vandermonde matrices, technical report. [EGO05a ] Yu.Eidelman, I.Gohberg and V.Olshevsky, Eigenstructure of Order-One- Quasiseparable Matrices. Three-term and Two-term Recurrence Relations, Linear Algebra and its Applications, Volume 405, 1 August 2005, Pages 1-40. [EGO05b ] Yu.Eidelman, I.Gohberg and V.Olshevsky, The QR iteration method for Hermitian quasiseparable matrices of an arbitrary order, Linear Algebra and its Applications, Volume 404, 15 July 2005, Pages 305-324

– Typeset by FoilT EX – 58

Efficient Algorithms for Structured Matrices Vadim Olshevsky

[OS03 ] V. Olshevsky and M. A. Shokrollahi, A displacement approach to decoding algebraic codes, in Fast Algorithms for Structured Matrices: Theory and Applications, CONM/323, p. 265 - 292, AMS publications, May 2003. [OS01a ] V.Olshevsky and Michael Stewart. Stable Factorization of Hankel and Hankel-like Matrices. Numerical Linear Algebra, Vol 8, Issue 6/7, p 401-434, 2001. [OS01b ] V.Olshevsky and Amin Shokrollahi, A Superfast Algorithm for Confluent Rational Tangential Interpolation Problem via Matrix-vector Multiplication for Confluent Cauchy-like Matrices, in “Structured Matrices in Mathematics, Computer Science, and Engineering I,” Contemporary Mathematics series, Vol 280, 32-46, AMS, 2001. [OS00 ] V.Olshevsky and A.Shokrollahi, Fast matrix-vector multiplication algorithms for confluent Cauchy-like matrices with applications, in Proc of of the Thirty Second ACM Symposium on Theory of Computing (STOC’00), p.573-581; ACM, New York, 2000.

– Typeset by FoilT EX – 59