Natrag: Numeričko računanje QR   Gore: QR RASTAV   Naprijed: Ekonomični QR rastav  

Rješavanje problema najmanjih kvadrata pomoću QR rastava

Koristeći svojstvo ortogonalne matrice $Q$ da je $\Vert Qx\Vert=\Vert x\Vert$ možemo lako riješiti problem najmanjih kvadrata. Neka je $A$ matrica tipa $m\times n$, $m>n$, neka je $\mathop{\mathrm{rang}}\nolimits A=n$ i neka je $A=QR$ rastav matrice $A$. Tada je

$$\Vert Ax-b\Vert^2= \Vert QRx-Q Q^T b\Vert^2= \Vert Q(Rx-Q^T b)\Vert^2 = \Vert Rx-Q^T b \Vert^2.$$

Matricu $R$ možemo zapisati kao

$$R=\begin{bmatrix}R_0 0 \end{bmatrix},$$ (2.5)

gdje je $R_0$ gornje trokutasta kvadratna matrica reda $n$, a vektor $Q^Tb$ možemo zapisati kao

$$Q^T b=\begin{bmatrix}c d \end{bmatrix},$$

gdje je $c$ dimenzije $n\times 1$ i $d$ dimenzije $(m-n)\times 1$. Sada imamo

$$\Vert Ax-b\Vert^2=\Vert Rx-Q^T b \Vert^2 = \left\Vert \begin{bmat... ... x -c -d \end{bmatrix}\right\Vert^2 = \Vert R_0 x-c\Vert^2+\Vert d\Vert^2.$$

Kako je $\mathop{\mathrm{rang}}\nolimits A=n$, iz svojstva QR3. zaključujemo da trokutasti sustav $R_0x=c$ ima jedinstveno rješenje $x$ za koje je $\Vert R_0x-c\Vert=0$. Kako $d$ ne ovisi o $x$ vrijednost $\Vert Ax-b\Vert$ se ne može više smanjiti pa je $x$ upravo (jedinstveno) rješenje problema najmanjih kvadrata. Očito je

$$\min_x \Vert Ax-b\Vert =\Vert d\Vert.$$

Slijedeći Matlab potprogram rješava problem najmanjih kvadrata pomoću QR rastava:

function x=LS(A,b)
% Rjesava problem najmanjih kvadrata || Ax-b || --> min
% za matricu punog stupcanog ranga A koristeci QR rastav.
[m,n]=size(A);
[Q,R]=moj_QR(A);
b1=Q'*b;
c=b1(1:n);
x=R(1:n,1:n)\c;
end

Zadatak 2.3   Riješite problem najmanjih kvadrata iz primjera 1.2 pomoću prethodnog potprograma i usporedite rješenje s rješenjima dobivenim korištenjem normalne jednadžbe i Matlabove naredbe x=A\b.

Octave On-line

A = [1 1 1 4 2 1 16 4 1 25 5 1 36 6 1] b = [0 1 4 8 14]' % Rjesenje problema najmanjih kvadrata pomocu QR rastava 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 function x=LS(A,b) % Rjesava problem najmanjih kvadrata || Ax-b || --> min % za matricu punog stupcanog ranga A koristeci QR rastav. [m,n]=size(A); [Q,R]=moj_QR(A); b1=Q'*b; c=b1(1:n); x=R(1:n,1:n)\c; end x=LS(A,b) % Matlabovo rjesenje xm=A\b

     

  

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

Napomena 2.2   Vidjeli smo da problem najmanjih kvadrata možemo rješavati na dva načina: pomoću normalne jednadžbe i pomoću QR rastava. Rješavanje problema najmanjih kvadarata pomoću QR rastava je otprilike dva puta sporije ali zato ima bolja numerička svojstva odnosno u određenim situacijama daje točnije rješenje.

Natrag: Numeričko računanje QR   Gore: QR RASTAV   Naprijed: Ekonomični QR rastav