Natrag: QR rastav matrice   Gore: QR RASTAV   Naprijed: Rješavanje problema najmanjih  

Numeričko računanje QR rastava

Za formiranje Householderove matrice $H$ potrebno je $O(n^2)$ računskih operacija. Za množenje nekog vektora $x$ Householderovom matricom također je potrebno $O(n^2)$ računskih operacija. No, umnožak $Hx$ može se izračunati i bez formiranja matrice $H$:

$$Hx=\left(I-2\frac{vv^T}{v^Tv}\right) x= x-v \frac{2(v^Tx)}{v^Tv}.$$

Za računanje produkta pomoću druge jednakosti potrebno je samo $O(6n)$ računskih operacija. Slično, produkt $HA$ se u praksi računa pomoću formula:

$$\beta =-\frac{2}{v^Tv},$$

$$w = \beta A^T v$$ (2.4)

$$HA = A+v w^T.$$

Koristeći prethodna poboljšanja možemo napisati potprograme za računanje QR rastava. Prvi potprogram računa vektor $v$ iz zadanog vektora $x$ prema formuli (2.3) i napomeni 2.1:

function v=House(x)
% racuna Householderov vektor v iz vektora x
sig=sign(x(1)); if sig==0, sig=1; end
v=x;
v(1)=v(1)+sig*norm(v);
end

Slijedeći potprogram računa umnožak $HA$ bez formiranja matrice $H$ prema formulama (2.4):

function B=mnozi_House(v,A)
% racuna produkt HA=(I-2*v*v'/(v'*v))*A
% bez formiranja matrice H
beta=-2/(v'*v);
w=beta*A'*v;
B=A+v*w';
end

Konačno, slijedeći program računa QR rastav matrice $A$ oprema postupku opisanom u prethodnom poglavlju:

function [Q,R]=moj_QR(A)
% QR rastav A=Q*R. A je m x n i mora biti rank(A)n.
[m,n]=size(A);
Q=eye(m);
for i=1:min(m-1,n)
   v=House(A(i:m,i));
   A(i:m,i:n)=mnozi_House(v, A(i:m,i:n));
   Q(i:m,:)=mnozi_House(v, Q(i:m,:));
end
R=A;
Q=Q';
end

Zadatak 2.2   Izračunajte QR rastav matrice $A$ iz primjera 1.2 pomoću prethodnih potprograma i usporedite rješenje s onim koje daje Matlabova naredba [Q,R]=qr(A).

Octave On-line

A = [1 1 1 4 2 1 16 4 1 25 5 1 36 6 1] % QR rastav pomocu vlastitih potprograma function v=House(x) % racuna Householderov vektor v iz vektora x sig=sign(x(1)); if sig==0, sig=1; end v=x; v(1)=v(1)+sig*norm(v); end function B=mnozi_House(v,A) % racuna produkt HA=(I-2*v*v'/(v'*v))*A % bez formiranja matrice H beta=-2/(v'*v); w=beta*A'*v; B=A+v*w'; end function [Q,R]=moj_QR(A) % QR rastav A=Q*R. A je m x n i mora biti rank(A)n. [m,n]=size(A); Q=eye(m); for i=1:min(m-1,n) v=House(A(i:m,i)); A(i:m,i:n)=mnozi_House(v, A(i:m,i:n)); Q(i:m,:)=mnozi_House(v, Q(i:m,:)); end R=A; Q=Q'; end % moj QR rastav [Q,R]=moj_QR(A) % Matlabov QR rastav [Qm,Rm]=qr(A)

     

  

[Octave On-line Home]    [Octave User's Guide]

Natrag: QR rastav matrice   Gore: QR RASTAV   Naprijed: Rješavanje problema najmanjih