Perplektički QR rastav centrosimetričnih matrica (seminar, 2012)

Perplektički QR rastav centrosimetričnih matrica Konrad Burnik PMF - Matematički odjel Seminar za numeričku matematiku i računarstvo

28. lipnja 2012.

Uvod

 I

1 ..

 Definiramo reverz matrice identiteta: Rn := 

.

  

1 I

Za svaki n ∈ N vrijedi Rn2 = In .

I

Definiramo (indefinitni) skalarni produkt obzirom na Rn sa hx, y iRn := x1 yn + x2 yn−1 + · · · + xn y1 = x T Rn y kojeg zovemo perplektički skalarni produkt.

Uvod

I

Skalarnom produktu hx, y iRn pridružen je funkcional („norma“): q(x) := hx, xiRn

I

Za vektor x ∈ Rn kažemo da je izotropan ako i samo ako je q(x) = 0.

I

Primjetimo da iz q(x) = 0 ne slijedi x = 0, dakle postoje izotropni vektori koji nisu jednaki nul vektoru!

Perplektičke matrice

Posebno će nas zanimati matrice koje čuvaju perplektički skalarni produkt tj. matrice G ∈ Rn×n takve da vrijedi: hGx, Gy iRn = hx, y iRn , ∀x, y ∈ Rn Takve matrice zovemo perplektičke matrice.

Theorem Matrica G ∈ Rn je perplektička ako i samo ako G T Rn G = Rn . Napomena: Perplektičke matrice čuvaju izotropnost vektora!

Perplektički QR rastav

I

Zadanu matricu A ∈ Rn×n želimo rastaviti na faktore A = QR pri čemu je matrica Q perplektička. Osnovno pitanje je kakav je „oblik“ od R?

I

Promatramo sve najprije u prostoru parne dimenzije tj. za n = 2k.

Perplektički QR rastav I

Želimo da R ima što više nula! Možemo pokušati uzeti npr. gornjetrokutastu matricu za R kao u klasičnom QR tj.   ××××××  × × × × ×    × × × × ?   R=  × × ×    × × ×

I

Problem: Prva polovica stupaca u R su izotropni vektori, pa stoga da bi postojala faktorizacija sa ovakvim R prva polovica stupaca u početnoj matrici A također moraju svi biti izotropni!

I

Slično zaključujemo da niti blok-gornjetrokutasta forma nije dobra.

Perplektički QR rastav

I

Promotrimo A koji se sastoji od samo jednog stupca.   × ×   ×  A= ×   × × Koje elemente poništiti ?

Perplektički QR rastav I

Ako poništimo sve elemente osim jednog, tada uvijek dobivamo izotropan vektor. Npr.     × ×   ×       ×  →    ×       × ×

I

Općenito, vektori kanonske baze ej u prostoru parne dimenzije su svi izotropni vektori pa ne možemo preslikavati neizotropne u njih!

I

Moramo dopustiti i preslikavanje neizotropnih vektora!

Perplektički QR rastav

I

Poništavamo dakle sve osim dva elementa.

I

Problem: ako ti elementi nisu na simetrični oko sredine vektora onda ponovo imamo izotropan vektor!     × × ×         ×   6→   ×       × × ×

Perplektički QR rastav

I

Poništavamo tako da elementi koji ostaju budu na položajima simetričnim oko sredine!

I

Imamo dvije prirodne mogućnosti za ovakvo poništavanje, od kojih obje mogu biti slika izotropnih i neizotropnih vektora, te koje su definirane za svaku duljinu n = 2k.         × × × ×   ×           ×         →   ili  ×  →  ×  ×   × ×         ×   ×   × × ×

Perplektički QR rastav

I

Pravokutna matrica A ∈ Rn×m je po definiciji centrosimetrična ako vrijedi ARm = Rn A Mi ćemo se baviti centrosimetričnim kvadratnim matricama tj. matricama A ∈ Rn×n takvima da vrijedi ARn = Rn A.

I

Napomena: Perplektičke matrice čuvaju svojstvo centrosimetrije!

Dijamantna forma (♦-forma) matrice Definition Matrica A ∈ R2k×2k je u dijamantnoj formi (skraćeno ♦-formi) ako se može prikazati u obliku   A11 A12 A= A21 A22 Gdje je I

A11 ∈ Rk×k donjetrokutasta obzirom na sporednu dijagonalu tj. aij = 0 za i > k − j

I

A12 ∈ Rk×k donjetrokutasta tj. aij = 0 za i > j

I

A21 ∈ Rk×k gornjetrokutasta obzirom na sporednu dijagonalu tj. aij = 0 za i < j

I

A22 ∈ Rk×k gornjetrokutasta tj. aij = 0 za i < k − j

Dijamantna forma (♦-forma) matrice

Example Za n = 6 sljedeća matrica A je u ♦-formi. Simbol “×” označava ne-nul elemente.   ××  ××××    × × × × × ×   A=  × × × × × ×    ××××  ××

Perplektički QR rastav

Definition Perplektički QR rastav centrosimetrične matrice A ∈ Rn×n je prikaz matrice A kao produkt oblika A = QR pri čemu je Q perplektička, a R je u ♦-formi.

Osnovni korak redukcije I

Želimo odjednom poništavati po dva stupca matrice (najjednostavnije za centrosimetrične matrice)

I

Neka je A ∈ R2k×2k centrosimetrična matrica. Tada je A oblika   A = a1 a2 . . . ak akR . . . a2R a1R

I

Za zadani stupac x ∈ {a1 , . . . , ak } tražimo perplektičko preslikavanje oblika:     x1 xn  x2 xn−1   α1 α2      .. ..  →   α2 α1   . .  xn x1 pri čemu su α1 , α2 parametri koje još treba odrediti.

I

Dodatno zahtjevamo da preslikavanje čuva normu stupaca (do na predznak)!

Osnovni korak redukcije

I

Iz zadanih uvjeta na traženo preslikavanje, parametre α1 ,α2 možemo odrediti iz sljedećeg sustava: ( α12 + α22 = kxk22 2α1 α2 = q(x)

Osnovni korak redukcije (određivanje parametara α1 , α2 )

I

U slučaju q(x) = 0 (izotropan x) imamo α1 = 0 ili α2 = 0 (bez smanjena općenitosti uzmimo α2 = 0) pa preslikavanje izgleda:     x1 xn  x2 xn−1    α1     .. ..  →   α1   . .  xn x1 pri čemu je α1 = ± kxk2 .

Osnovni korak redukcije (određivanje parametara α1 , α2 )

I

U slučaju q(x) 6= 0 (neizotropan x) za rješenje sustava uzimamo: r q kxk22 + kxk22 − q(x)2 √ α1 = 2 2 2(kxk2 α1 − α13 ) α2 = q(x)

I

Napomena: Sada znamo što se treba preslikati u što, sada još treba odrediti matricu traženog preslikavanja!

Problem preslikavanja (općeniti)

Problem preslikavanja: Za zadane vektore x, y ∈ Rn pronaći matricu A ∈ Rn×n takvu da Ax = y .

Strukturirani problem preslikavanja

I

Želimo da matrica A koja rješava problem preslikavanja ima dodatnu strukturu tj. da bude perplektička!

I

Poništavanje elemenata raznim vrstama perplektičkih matrica (dvostruke Givensove rotacije, dvostruki Householder, G-reflektori) opisano je u [2].

I

Mi ćemo koristiti (blok) G-reflektore iz [4]!

G-reflektori

I

Problem preslikavanja je djelomično riješen pomoću G-reflektora [3, 1]. Neka su x, y ∈ Rn i u := y − x. Tada za G-reflektor 2uu T Rn G = In − q(u) vrijedi G ∈ G i Gx = y . (dokaz se može pronaći u [3])

I

Problem 1: Kada je u izotropan tj. q(u) = 0 tada G nije definiran.

I

Problem 2: Nakon što smo pronašli G1 takav da G1 x = y , tada moramo pronaći G2 takav da G2 G1 x = G1 x i G2 x R = y R što je dodatan korak.

Blok G-reflektori I

Spomenute probleme rješavamo blok G-reflektorima.

I

Za zadane x, y ∈ Rn neka su X , Y ∈ Rn×2 oblika     X = x xR i Y = y yR te stavimo V := Y − X .

I

Blok G-reflektor [1][4] je matrica oblika: G := In − 2V [V T Rn V ]† V T Rn pri čemu je sa [V T Rn V ]† označen Moore-Penroseov generalizirani inverz od V T Rn V .

I

Vrijedi G ∈ G i GX = Y (dokaz se može pronaći u [4]) čime smo riješili problem preslikavanja!

Blok G-reflektori

I

Mi ćemo uzeti da su matrice X , Y , V ∈ Rn×2 oblika:   X := x x R   Y := α1 em + α2 em+1 α1 em+1 + α2 em V :=Y − X

Blok G-reflektori

I

Problem: Kako računati [V T Rn V ]† ?

I

Primjetimo V T Rn V =



q(v ) kv k22 kv k22 q(v )



dakle V T Rn V je 2 × 2 (centrosimetrična) matrica!

Blok G-reflektori

I

U slučaju detV T Rn V 6= 0 imamo [V T Rn V ]† = [V T Rn V ]−1 što znamo izračunati.

I

Vrijedi detV T Rn V = 0 ako i samo ako q(v ) = ± kv k22 , pa razlikujemo tri slučaja.

Blok G-reflektori I

Ako q(v ) = kv k22 > 0 onda 

11 [V Rn V ] = q(v ) 11 †

T

I

  1 11 = 4q(v ) 1 1

†

  1 1 −1 = 4q(v ) −1 1

Ako q(v ) = kv k22 < 0 onda 

1 −1 [V Rn V ] = q(v ) −1 1 T

I

†

†

Ako q(v ) = kv k22 = 0 onda 

00 [V Rn V ] = 00 T

†

†



00 = 00



Blok G-reflektori

Vektori v ∈ Rn za koje vrijedi q(v ) = ± kv k22 su točno oni koji su centrosimetrični tj. vrijedi sljedeća činjenica.

Fact Neka je v ∈ Rn . Tada je v = ±v R ako i samo ako q(v ) = ± kv k22 .

Osnovni korak redukcije

× ×  ×  ×  × × 

× × × × × ×

× × × × × ×

× × × × × ×

× × × × × ×

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

× × × × × ×

× × × × × ×

× × × × × ×

 × ×   × ×  × ×  ×  ×

Perplektički QR rastav

I

Kako se nakon primjene jednog reflektora osigurati da ne pokvarimo redukciju u prethodom koraku ?

I

Treba nam pojam križnog ulaganja matrice.

Križno ulaganje matrice (cross-embedding)

Definition 2k×2k , zovemo križnim ulaganjem matrice Matricu  X ∈ R A11 A12 A= ∈ R2k×2k u matricu In za n ≥ k ako i samo ako je A21 A22 X oblika   A11 A12  In−2k X = A21 A22

pri čemu su Aij kvadratne matrice reda k.

Perplektički QR rastav (nastavak)

I

Radimo blok G-reflektora  križno-ulaganje  G11 G12 G= ∈ R2k×2k . G21 G22 

G11

˜ = G

In−2k G21

I

G12

 

G22

Križno ulaganje perplektičke (ortogonalne) matrice je ponovo perplektička (ortogonalna)!

Redukcija pomoću križnog ulaganja

×  ×  × ×  × ×   × × 

× × × × × ×

× × × × × ×

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

× × × × × ×

 ×   × ×  × ×  × 

Perplektički QR rastav (nastavak)

I

U svakom idućem koraku redukcije djelujemo sa križno uloženim osnovnim blok reflektorom za k = 1, 2, . . . , n koji djeluje na stupce {(2, n − 1), (3, n − 2), . . . , (k − 1, k + 1)} redom.

Perplektički QR rastav (nastavak)

Na matricu A primjenjujemo niz (križno uloženih) blok G-reflektora i dolazimo do ♦-forme Gs . . . G2 G1 A = ♦ Stavimo Q := G1T G2T . . . GsT i R := ♦.Vrijedi A = QR Pri čemu je Q perplektička ortogonalna, a R je u ♦-formi.

Algoritam za perplektički QR rastav (slučaj n = 2k) 1: 2: 3: 4: 5:

BlockPerplecticQR(A) ˜ ← A;Q ← In ;d ← n A for k = 1 . . . n2 − 1 do x ←˜ a1 ;m ← d/2 if hx, xiRn 6= 0 then

procedure

r

6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20:

kxk2 2+

q

2 kxk2 2 −q(x)

. if x is non-isotropic 2(kxk2 α −α3 )

2 1 1 √ ; α2 ← ; α1 ← q(x)   2   R X ← x x ; Y ← α1 em + α2 em+1 α2 em + α1 em+1 ; else . if x is isotropic α1 ←kxk2    X ← x x R ; Y ← α1 em αem+1 ; end if V ← Y − X; Find a block-G-reflector G ← In − 2V [V T Rn V ]† V T Rn ;   G11 G12 ← G. G21 G22   G11 G12  Q In−2k Q← G21 G22  ˜←A ˜\ ˜ A a1 , ˜ ad , ˜ ad /2 , ˜ ad /2+1 d ←d −2 end for R ← QA return {Q, R} . Q is perplectic orthogontal, R is in ♦-form and A = Q T R

Perplektički QR rastav

Example Neka je A sljedeća centrosimetrična matrica (Toeplitzova matrica): 

1 2  3 A= 4  5 6

2 1 2 3 4 5

3 2 1 2 3 4

4 3 2 1 2 3

5 4 3 2 1 2

 6 5  4  3  2 1

Perplektički QR rastav

 0.711252 0.140987  −0.19544 0.542221 0.196768 −0.511667     −0.120788 0.373191 −0.0704935 0.827895   Q =  0.827895 0.20416 −0.337755 −0.120788     −0.511667 0.0351289 −0.605016 −0.19544  −0.133902 0.67551   −0.390879 −1.02333   1.35102 1.49201 0.957484 0.281974     9.02022 7.02036 5.09077 3.56949 2.79459 3.10414   R =  3.10414 2.79459 3.56949 5.09077 7.02036 9.02022     0.281974 0.957484 1.49201 1.35102 −1.02333 −0.390879 

0.67551 −0.605016 −0.337755 −0.0704935 0.196768 0.140987

−0.133902 0.0351289 0.20416 0.373191 0.542221 0.711252

Primjena Primjena (rješavanje centrosimetričnih sustava)

I

Perplektički QR rastav možemo primjeniti na rješavanje centrosimetričnih sustava.

Primjena Primjena (rješavanje centrosimetričnih sustava)

  (2)  r2,m−1   ..  ... .   (m)  rm,1 ×   r (m) ×  m+1,1  ..  ..  . .  (2)  rn−1,m−1 

(1)

(1)



r1,m r1,m+1 × .. . ×

× .. . ×

(2) r2,m+1

× .. . ×

× .. . ×

× .. .

(1)

(1)

rn,m rn,m+1

.. . ×

(2) rn−1,m+1

..

.

(m) rm,m (m) rm+1,m

..

.

(m)

x1   (m−1)  x  2  ..  .   (1)   xm    x (1)   m+1  ..   .   (m−1)  x n−1  (m) xm

 (1) b   1(2)   b2     ..   .     b (m)   m  =  (m)   bm+1     .   ..     (2)   bn−1 

(1)

                

bn

Figure: The scheme for solving a diamond form system Rx = b. The elements that are used in forming a 2 × 2 subsystem are marked by the number of the stage they are considered.

Primjena (rješavanje centrosimetričnih sustava) procedure SolveDiamondSystem(R,b) 2: m := n/2; 3: for k = 1 . . . m do

1:

. R is in ♦-form

4:

5: 6: 7: 8: 9: 10:

 sh := (xm−k−1 , . . . , xm+k+1 ), (rh,m−k−1 , . . . , rh,m+k+1 )   rk,m−k rk,m+k+1 Rk := ; r rn−k+1,m+k+1  n−k+1,m−k    xm−k bk − sk Xk := ; Bk := ; xm−k+1 bn−k+1 − sn−k+1 Xk ← Rk−1 Bk ; . Solve a 2 × 2 system end for return x end procedure

Slučaj n = 2k + 1 U slučaju neparne dimenzije prostora n = 2k + 1 najprije moramo definirati ♦-formu.

Definition Matrica A ∈ R2k+1×2k+1 je u dijamantnoj formi (skraćeno ♦-formi) ako se može prikazati u obliku   A11 × A12 A= × × ×  A21 × A22 Gdje je I

I I

I

A11 ∈ Rk×k donjetrokutasta obzirom na sporednu dijagonalu tj. aij = 0 za i > k − j A12 ∈ Rk×k donjetrokutasta tj. aij = 0 za i > j A21 ∈ Rk×k gornjetrokutasta obzirom na sporednu dijagonalu tj. aij = 0 za i < j A22 ∈ Rk×k gornjetrokutasta tj. aij = 0 za i < k − j

Slučaj n = 2k + 1

Example Za n = 7 sljedeća matrica A je u ♦-formi. Simbol “×” označava ne-nul elemente.   ×××  ×××××    × × × × × × ×    A= × × × × × × × × × × × × × ×    ×××××  ×××

Perplektički QR rastav za n = 2k + 1 I

˜ ∈ R2k+1×2k+1 se može particionirati Matrica A   A11 × A12 ˜= × × ×  A A21 × A22 pri čemu su blokovi Aij kvadratne matrice reda k.

I

I

I

Primjenimo algoritam za slučaj n = 2k na matrici   A11 A12 A= i dobijemo faktore Q i R reda 2k. A21 A22     Q11 Q12 R11 × R12 ˜ = ˜ =  × × ×  (pritom R 1 Sada stavimo Q Q21 Q22 R21 × R22 ˜ su sa„ד označeni prepisani elemeni iz A) ˜=Q ˜R ˜ pri čemu je Q ˜ perplektička ortogonalna, a R ˜ je Vrijedi A u ♦-formi.

Zaključak

I

Pokazali smo da vrijedi sljedeći teorem.

Theorem Neka je A ∈ Rn×n centrosimetrična matrica. Tada postoji perplektička faktorizacija A = QR pri čemu je Q perplektička ortogonalna, a R je u ♦-formi.

References

A. Bunse-Gerstner, An analysis of the HR algorithm for computing the eigenvalues of a matrix, Linear Algebra Appl. 35 (1981) 155–173. D. S. Mackey, N. Mackey, F. Tisseur, Structured tools for structured matrices, Electron. J. Linear Algebra 10 (2003) 106–145. D. S. Mackey, N. Mackey, F. Tisseur, G–reflectors: Analogues of Householder transformations in scalar product spaces, Linear Algebra Appl. 385 (2004) 187–213. S. Singer, S. Singer, Orthosymmetric block reflectors, Linear Algebra Appl. 429 (5–6) (2008) 1354–1385.