Musterlösung Aufgabe 28

Aufgabe 28: (QR-Zerlegung von 3-Band-Matrizen, technische Durchführung)

Sei A = ( \begin{matrix}β_{1} & γ_{1} & & 0 \\ α_{1} & ... & ... \\ & ... & ... & γ_{n−1} \\ 0 & & α_{n−1} & β_{n} \end{matrix} ).

a) Zeige induktiv

A^{(j)} := H^{(j−1)} ⋯ H^{(1)}A = ( \begin{matrix}u_{1} & v_{1} & w_{1} \\ & ... & ... & ... \\ & & ... & ... & ... \\ & & & u_{j−1} & v_{j−1} & w_{j−1} \\ & & & & β̃_{j} & γ̃_{j} & 0 \\ & & & & α_{j} & β_{j+1} & γ_{j+1} \\ & & & & 0 & α_{j+1} & β_{j+2} & γ_{j+2} \\ & & & & & & ... & ... & ... \end{matrix} ) ,

wobei rekursiv

H^{(j)} = ( 1I ; c_{j} s_{j} ; s_{j} −c_{j} ; 1I )\begin{matrix}← j \\ ← j + 1 \end{matrix} mit \begin{matrix}c_{j} = β̃_{j} / \sqrt{β̃_{j}^{2} + α_{j}^{2}} \\ s_{j} = α_{j} / \sqrt{β̃_{j}^{2} + α_{j}^{2}} \end{matrix}

(bzw. c_{j} = 1, s_{j} = 0 für α_{j} = 0) gewählt wird. Gestartet wird mit A^{(1)} := A, d.h.,

β̃_{1} := β_{1}, γ̃_{1} := γ_{1}. Anleitung: H^{(j)} wirkt nur auf die Zeilen j und j + 1 des Schemas.

Zeige, daß H^{(j)} orthogonal ist.

Anmerkung: es braucht nicht verifiziert zu werden, daß H^{(j)} der Householder-Konstruktion aus der Vorlesung entspricht (was –bis auf Vorzeichen– in der Tat der Fall ist).

b) Leite aus a) eine Rekursion zur Berechnung der 3 Bänder von

R = H^{(n−1)} ⋯ H^{(1)}A = ( \begin{matrix}u_{1} & v_{1} & w_{1} & & 0 \\ & ... & ... & ... \\ & & ... & ... & w_{n−2} \\ & & & ... & v_{n−1} \\ 0 & & & & u_{n} \end{matrix} )

her. Benutze dabei die Daten c_{j}, s_{j} der Matrizen H^{(j)}. Mit welchem Aufwand wird R ermittelt? Implementiere die Rekursion und gebe diese Daten für (α_{i}) = (−3, −2, …, 5), (β_{i}) = (1, 2, …, 10), (γ_{i}) = (1, 2, …, 9) an. c) Die Daten c_{j}, s_{j} der Householder-Matrizen aus a) liegen vor. Wie und mit welchem

Aufwand berechnet sich die Transformation Q^{T} b⃗ = H^{(n−1)} ⋯ H^{(1)}b⃗ eines beliebigen Vektors b⃗ ?

d) Zeige Q^{T} = H^{(n−1)} ⋯ H^{(1)} = ( \begin{matrix}∗ & ∗ & & 0 \\ ⋮ & ... & ... \\ ⋮ & ... & ... & ∗ \\ ∗ & … & … & ∗ \end{matrix} ).

Anleitung: H^{(j−1)} ⋯ H^{(1)} = ( ∗ ∗ ; ⋮ ... ... ; ∗ … ∗ ∗ 0 ; ∗ … ∗ ∗ 0 ; 0 … 0 0 1 ; I 1 )← j

Musterlösung: a) Induktion nach j. Für j = 1 ist nichts zu zeigen. Schritt j − 1 → j: sei H^{(j−1)} ⋯ H^{(1)}A von der angegebenen Form, wähle H^{(j)} wie angegeben. Da H^{(j)} nur auf die Zeilen j und j+1 des Zwischenschemas H^{(j−1)} ⋯ H^{(1)}A wirkt, bleiben alle Zeilen außer j, j + 1 in der angegebenen behaupteten Form. Die Zeilen j und j + 1 von H^{(j)} ⋯ H^{(1)}A berechnen sich zu

( \begin{matrix}c_{j} & s_{j} \\ s_{j} & −c_{j} \end{matrix} ) ( \begin{matrix}… & 0 & β̃_{j} & γ̃_{j} & 0 & 0 & … \\ … & 0 & α_{j} & β_{j+1} & γ_{j+1} & 0 & … \end{matrix} )

= ( \begin{matrix}… & 0 & c_{j}β̃_{j} + s_{j}α_{j} & c_{j}γ̃_{j} + s_{j}β_{j+1} & s_{j}γ_{j+1} & 0 & … \\ … & 0 & s_{j}β̃_{j} − c_{j}α_{j} & s_{j}γ̃_{j} − c_{j}β_{j+1} & −c_{j}γ_{j+1} & 0 & … \end{matrix} )

=: ( \begin{matrix}… & 0 & u_{j} & v_{j} & w_{j} & 0 & … \\ … & 0 & 0 & β̃_{j+1} & γ̃_{j+1} & 0 & … \end{matrix} ) .

Hierbei ist entscheidend, daß mit

c_{j} = β̃_{j} / \sqrt{β̃_{j}^{2} + α_{j}^{2}} , s_{j} = α_{j} / \sqrt{β̃_{j}^{2} + α_{j}^{2}}  (1)

das Element α_{j} eliminiert wird: s_{j}β̃_{j} − c_{j}α_{j} = \frac{α_{j}β̃_{j}}{\sqrt{β̃_{j}^{2} + α_{j}^{2}}} − \frac{β̃_{j}α_{j}}{\sqrt{β̃_{j}^{2} + α_{j}^{2}}} = 0 .

Dies beweist die in a) behauptete Form der Zwischenschemata A^{(j)} = H^{(j−1)} ⋯ H^{(1)}A. Zusätzlich ergibt sich die Rekursion

\begin{matrix}u_{j} = c_{j}β̃_{j} + s_{j}α_{j} , v_{j} = c_{j}γ̃_{j} + s_{j}β_{j+1} , w_{j} = s_{j}γ_{j+1} , \\ β̃_{j+1} = s_{j}γ̃_{j} − c_{j}β_{j+1} , γ̃_{j+1} = −c_{j}γ_{j+1} . \end{matrix}  (2)

Falls α_{j} = 0 gilt, so braucht dieses Element nicht mehr eliminiert zu werden. Man kann in diesem Fall c_{j} = 1, s_{j} = 0 setzen. In jedem Fall ist H^{(j)} orthogonal:

H^{(j)}(H^{(j)})^{T} = ( \begin{matrix}c_{j} & s_{j} \\ s_{j} & −c_{j} \end{matrix} ) ( \begin{matrix}c_{j} & s_{j} \\ s_{j} & −c_{j} \end{matrix} ) = ( \begin{matrix}c_{j}^{2} + s_{j}^{2} & 0 \\ 0 & s_{j}^{2} + c_{j}^{2} \end{matrix} ) = ( \begin{matrix}1 & 0 \\ 0 & 1 \end{matrix} ) .

b) Aus (1) und (2) ergibt sich die definierende Rekursion von u_{j}, v_{j}, w_{j}, j = 1, …, n − 2. Im letzten Schritt j = n − 1 kann man die nicht existenten w_{n−1} und γ̃_{n} vergessen. Es gilt β̃_{n} = u_{n}. Speichert man β̃_{j} als β_{j} und γ̃_{j} als γ_{j}, so ergibt sich folgende Pseudo-Implementierung:

for j = 1 .. n − 2 do begin

if α_{j} = 0 then c_{j} := 1 ; s_{j} := 0 ; else t := \sqrt{β_{j}^{2} + α_{j}^{2}} ; c_{j} := β_{j}/t; s_{j} := α_{j}/t;

end if;

u_{j} := c_{j}β_{j} + s_{j}α_{j} ; v_{j} := c_{j}γ_{j} + s_{j}β_{j+1} ; w_{j} := s_{j}γ_{j+1} ;

β_{j+1} := s_{j}γ_{j} − c_{j}β_{j+1} ; γ_{j+1} = −c_{j}γ_{j+1} ;

end for;

if α_{n−1} = 0 (∗ letzter Schritt j = n − 1: ∗)

then c_{n−1} := 1 ; s_{n−1} := 0 ; else t := \sqrt{β_{n−1}^{2} + α_{n−1}^{2}} ; c_{n−1} := β_{n−1}/t; s_{n−1} := α_{n−1}/t;

end if;

u_{n−1} := c_{n−1}β_{n−1} + s_{n−1}α_{n−1} ; v_{n−1} := c_{n−1}γ_{n−1} + s_{n−1}β_{n}; u_{n} := s_{n−1}γ_{n−1} − c_{n−1}β_{n};

Rechenaufwand: 12 (n − 2) + 10 = 12 n − 2 Multiplikationen, 5 (n − 2) + 4 = 5 n + 6 Additionen, n − 1 Wurzeln.

|\^/|     Maple V Release 4 (Universitaet GH Paderborn)
._|\|   |/|_. Copyright (c) 1981-1996 by Waterloo Maple Inc. All rights
\  MAPLE  /  reserved. Maple and Maple V are registered trademarks of
<________>  Waterloo Maple Inc.
|       Type ? for help.

> n:=10:
> for i to n do beta[i]:=1.0*i: od:
> for i to n-1 do alpha[i]:=1.0*(i-4): od:
> for i to n-1 do Gamma[i]:=1.0*i: od:
> ...
> print(A):
[ 1.0  1.0   0     0    0    0    0    0    0    0  ]
[-3.0  2.0  2.0    0    0    0    0    0    0    0  ]
[  0  -2.0  3.0   3.0   0    0    0    0    0    0  ]
[  0    0  -1.0   4.0  4.0   0    0    0    0    0  ]
[  0    0    0     0   5.0  5.0   0    0    0    0  ]
[  0    0    0     0   1.0  6.0  6.0   0    0    0  ]
[  0    0    0     0    0   2.0  7.0  7.0   0    0  ]
[  0    0    0     0    0    0   3.0  8.0  8.0   0  ]
[  0    0    0     0    0    0    0   4.0  9.0  9.0 ]
[  0    0    0     0    0    0    0    0   5.0 10.0 ]
> print(R);
[3.162277660 , -1.581138830 , -1.897366596 , 0 , 0 , 0 , 0 , 0 , 0 , 0]
[0 , 2.549509757 , -1.961161352 , -2.353393622 , 0 , 0 , 0 , 0 , 0 , 0]
[0 ,  0 ,  2.560048076 , .150237563 , -1.562470657 , 0 , 0 , 0 , 0 , 0]
[0 ,  0 ,  0 ,  -4.408964407 ,  -3.682212032 ,  0  ,  0 ,  0 ,  0 ,  0]
[0 ,  0 ,  0 , 0 , 5.099019513 , 6.079600189 , 1.176696811 , 0 , 0 , 0]
[0 ,  0 ,  0 , 0 , 0 , 5.295135649 , 8.091606465 , 2.643936044 , 0 , 0]
[0 ,  0 ,  0 , 0 , 0 , 0 , 5.209730265 , 9.905754690 , 4.606764414 , 0]
[0 ,  0 ,  0 , 0 , 0 , 0 , 0 , 4.887292319 , 11.12405771 , 7.366041900]
[0 ,  0 ,  0 , 0 , 0 , 0 , 0 ,  0 ,           5.00330508 , 9.805463395]
[0 ,  0 ,  0 , 0 , 0 , 0 , 0 ,  0 ,   0 ,                 -5.531212734]

c) Mit

H^{(j)}x⃗ = ( 1I ; c_{j} s_{j} ; s_{j} −c_{j} ; 1I ) ( ⋮ ; x_{j} ; x_{j+1} ; ⋮ ) = ( ⋮ ; c_{j}x_{j} + s_{j}x_{j+1} ; s_{j}x_{j} − c_{j}x_{j+1} ; ⋮ )

ergibt sich x⃗ = H^{(n−1)} ⋯ H^{(1)}b⃗ durch die überschreibende Rekursion x⃗ := H^{(j)}x⃗, j = 1, …, n − 1:

x⃗ := b⃗;

for j = 1 … n − 1 do begin

t := x_{j} ;

x_{j} := c_{j} t + s_{j} x_{j+1} ; x_{j+1} := s_{j} t − c_{j} x_{j+1} ;

end ; mit 4 (n − 1) Multiplikationen.

d) Induktion nach j. Start: H^{(1)} hat die behauptete Form. Schritt j − 1 → j: H^{(j)} (H^{(j−1)} ⋯ H^{(1)}) =

(  1I ; c_{j} s_{j} ; s_{j} −c_{j} ; 1I ) ( ∗ ∗ ; ⋮ ... ... ; ∗ … ∗ ∗ 0 ; ∗ … ∗ ∗ 0 ; 0 … 0 0 1 ; I 1 ) = ( ∗ ∗ ; ⋮ ... ... ; ∗ … ∗ ∗ 0 ; ∗ … ∗ ∗ ∗ ; ∗ … ∗ ∗ ∗ ; I 1 ) .