8.5 Beispiel Prototyping:
Update der QR-Zerlegung

In der Statistik kommt es oft vor, dass eine Serie von Ausgleichsproblemen gelöst werden muss, wobei die Matrix $A \in {\mathbb{R}}^{m\times m}$ jeweils nur wenig ändert. Wir betrachten dazu die folgende Problemstellung:

Gegeben sei die QR-Zerlegung $A=Q{R\choose 0}$ mit $Q \in {\mathbb{R}}^{m\times m}$. Gegeben seien ferner die Vektoren $u \in {\mathbb{R}}^m$ und $v \in {\mathbb{R}}^n$.

Kann man die QR-Zerlegung der Matrix

$$A' = A + uv^T = Q' R'$$

einfacher aus der QR-Zerlegung von $A$ berechen?

Es ist

$$Q^T A' = Q^TA + Q^Tuv^T = {R\choose 0} + w v^T$$   mit$$\quad w = Q^Tu.$$

Sei $J_k = G(k,k+1,\phi_k)$ eine (orthogonale) Givensrotationsmatrix, welche bei der Multiplikation $J_k w$ die Elemente $w_k$ und $w_{k+1}$ verändert.

Der Algorithmus besteht aus drei Schritten:

1. Wir wählen Rotationen $J_k$ und Winkel $\phi_k$ so, dass die Elemente $n+2:m$ von $w$ zu Null rotiert werden:

   $$\underbrace{J_{n+1}^T \ldots J_{m-2}^T J_{m-1}^T Q^T}_{Q_1^T} A' = {R\choose 0} + J_{n+1}^T \ldots J_{m-2}^T J_{m-1}^T w v^T$$

   $$= {R\choose 0} + {w'\choose 0} v^T$$

   
   
   Dabei wird nur die Matrix $Q$ zu $Q_1$ und der Vektor $w$ geändert.

2. Nun werden weitere Givensrotationen durchgeführt, welche die Komponenten von $w'$ bis auf die erste $\alpha$ annullieren. Dabei ändert $Q_1$ zu $Q_2$ und auch die Matrix $R$ zu einer Hessenbergmatrix $H_1$.
   $$\underbrace{J_1^T \ldots J_{n}^T Q_1^T}_{Q_2^T} A' = H_1 + {\alpha \choose 0} v^T = H$$
   Der Beitrag der Rang-1 Matrix ändert nun in der letzten Gleichung nur die erste Zeile von $H_1$; man erhält eine neue Hessenbergmatrix $H$.

3. Nun wird im letzten Schritt durch weitere Givensrotationen die untere Nebendiagonale von $H$ annulliert und man erhält die gesuchte neue QR-Zerlegung.
   $$\underbrace{J_{n}^T \ldots J_{1}^T Q_2^T}_{Q'^T} A' = J_{n}^T \ldots J_{1}^T H = R'$$

Um diesen Algorithmus zu implementieren, konstruieren wir zunächst eine Demonstrationsversion. Dazu benötigt man die Givensrotationsmatrizen $J_k$:

function [G ] = rot(m,i,k,x,y);
%  [G ] = rot(m,i,k,x,y);
% konstruiert eine orthogonale Matrix G mxm,  welche y zu null rotiert
G = eye(m,m);
if y~=0,
   cot = x/y; si = 1/sqrt(1+cot^2); co = si*cot;
   G(i,i) = co; G(k ,k) = co; G(i,k) = si; G(k,i) = -si;
end;

Damit kann das folgende Programm geschrieben werden. Nach jeder Rotation wurde ein Pausenbefehl eingefügt, man kann damit den Ablauf des Algorithmus leicht verfolgen.

A= rand(7,4);
u=rand(7,1);
v=[1 2 3 4]';
[m n] = size(A);
[Q R] = qr(A);
QS =Q;
RS = R;
%
w = Q'*u
disp(' 1. Phase:  Anullieren  w(n+2:m) und aendern nur Q')
for i=m:-1:n+2,
  G= rot(m,i-1,i,w(i-1), w(i))
  w = G*w
  QS = QS*G';
  pause
end;
disp(' 2. Phase  Annullieren w(2:n+1) und aendern R und Q')
for i= n+1:-1:2,
  G= rot(m,i-1,i,w(i-1), w(i))
  w = G*w
  RS = G*RS
  QS = QS*G';
  pause
end
disp(' Addieren jetzt die Rang 1 Mod. zur ersten Zeile von R')
RS = RS + w*v'
pause
disp(' 3.Phase:  RS ist nun Hessenberg und wird reduziert')
for i = 1:n,
  G= rot(m,i,i+1,RS(i,i), RS(i+1,i))
 RS = G*RS
  QS = QS*G';
  pause
end
disp(' Kontrolle')
AS = A + u*v'; [qs rs] = qr(AS); 
qs
QS
rs
RS

Nun ist es einfach die Demonstrationsversion zu einem effektiven Programm umzuschreiben. Wir benötigen dazu eine Funktion giv, um die Rotationswinkel, bzw. den $\sin$ und $\cos$ davon zu berechnen:

function [s,c] = giv(x,y);
% bestimmt eine Givensrotation,
% welche y zu null rotiert
if y~=0,
   cot = x/y; s = 1/sqrt(1+cot^2); c = s*cot;
else c=1; s=0;
end

Im nachfolgenden Programm werden die Transformationen nicht mehr durch Multiplikation mit vollen Matrizen ausgeführt. Es werden nur noch jene Elemente neu berechnet, die von den Givensrotationen betroffen sind.

Man beachte, dass in MATLAB Vektoroperationen sehr gut dargestellt werden können, z. B. wird durch

     rs(i,i-1:n) = -si*rs(i-1,i-1:n) + co*rs(i,i-1:n);

die $i$-te Zeile der Matrix $R$ neu berechnet. Das Programm kann dadurch mit wenig Aufwand auf einen Vektorrechner übertragen werden.

A= round(10*rand(7,4))
 u=rand(7,1)
 v=[1 2 3 4]'

[m n] = size(A);  
[Q R] = qr(A); QS =Q; RS = R;
%
w = Q'*u
%  anullieren w(n+2:m) und aendern nur Q
for i=m:-1:n+2,
  [si , co]= giv(w(i-1), w(i));
  w(i-1) = co*w(i-1) + si*w(i);
%  QS = QS*G';
  h = QS(:,i-1);
  QS(:,i-1) = co*h + si*QS(:,i);
  QS(:,i) = -si*h + co*QS(:,i)
end
%  annullieren w(2:n+1) und aendern R und Q
for i= n+1:-1:2,
  [si , co]= giv(w(i-1), w(i));
%  w = G*w
  w(i-1) = co*w(i-1) + si*w(i);
%  RS = G*RS
  h = co*RS(i-1,i-1:n) + si*RS(i,i-1:n);
  RS(i,i-1:n) = -si*RS(i-1,i-1:n) + co*RS(i,i-1:n);
  RS(i-1,i-1:n) = h
%  QS = QS*G';
  h = QS(:,i-1);
  QS(:,i-1) = co*h + si*QS(:,i);
  QS(:,i) = -si*h + co*QS(:,i)
end
% Addieren jetzt die Rang 1 Mod. zur eRSten Zeile von R
RS(1,:) = RS(1,:) + w(1)*v'
% RS ist nun Hessenberg und wird reduziert
for i = 1:n,
  [si , co]= giv(RS(i,i), RS(i+1,i))
  % RS = G*RS
  h = co*RS(i,i:n) + si*RS(i+1,i:n);
  RS(i+1,i:n) = -si*RS(i,i:n) + co*RS(i+1,i:n);
  RS(i,i:n) = h
  %  QS = QS*G';
  h = QS(:,i);
  QS(:,i) = co*h + si*QS(:,i+1);
  QS(:,i+1) = -si*h + co*QS(:,i+1)
end
% Kontrolle
AS = A + u*v'; [qs rs] = qr(AS);
qs
QS
rs
RS