Matrici di Householder

L' algoritmo parte con

$$A= \begin{pmatrix} a_{1 \: 1 } ^0 & \ldots & \ldots & a_{1 \: n} ... ...vdots \\ a_{n \: 1 } ^0 & \ldots & \ldots & a_{n \: n} ^0 \end{pmatrix}= A^0$$

e al primo passo procede premoltiplicando a sinistra la matrice A per una matrice $P^1$ in modo che la prima colonna di $A$ sia un multiplo di $e_1$, $\det(A) \neq 0 \; \rightarrow \; Ae_1 \neq \underline{0}$ ed esiste $v_1 \in R^n$ tale che

$$P^1 A e_1 = ( I_n -2 \frac{v_1 v_1 ^T}{v_1 ^T v_1}) \begin{pmatri... ...end{pmatrix} = \begin{pmatrix} a_{1\:1} ^1 \\ 0 \\ \vdots \\ 0 \end{pmatrix}$$

e questo $v_1$ è appunto il vettore di Householder o di riflessione associato alla prima colonna della matrice A. Dopo questa moltiplicazione la situazione diventa la seguente:

$$P^1 A^0 = A^1 = \begin{pmatrix} a_{1\:1}^1 & a_{1\:2} ^1 & \ldots... ...& \vdots \\ 0 & a_{n\:2} ^1 & \ldots & \ldots & a_{n \:n} ^1 \\ \end{pmatrix}$$

dove si nota che, a differenza della fattorizzazione LU, al primo passo anche la prima riga viene alterata.
Il secondo passo inizia considerando la matrice $A^1$ come triangolare a blocchi

$$A^1 = \left( \begin{array}{c\vert ccc} a_{1\:1}^1 & a_{1\:2... ...e 0 & & & \\ \vdots & & A_1 & \\ 0 & & & \end{array}\right)$$

con $A_1$ non singolare, ed esiste una matrice di riflessione $P_2$ che premoltiplicata ad $A_1$ trasforma la sua prima colonna in un multiplo di $e_1$ (con $n-1$ componenti); questa matrice non può essere moltiplicata per $A^1$ per incompatibilità di dimensione, quindi è necessario trovare una matrice ortogonale che sia moltiplicabile per $A^1$, che trasformi la prima colonna di $A_1$ nel modo desiderato e che non alteri il resto della matrice. Questa matrice esiste ed è definita come

$$P^2= \left( \begin{array}{c\vert c c c c } 1 & & \bigcirc ^T... ...circ& & & & \\ & &P_2 A_1 & & \\ & & & & \end{array}\right)$$

Gli altri passi sono identici, via via si lavora su sottomatrici più piccole tutte non singolari, si trova una matrice di riflessione opportuna e la si ``inserisce'' moltiplicata per la matrice corrente dentro una matrice più grande con tutti 1 sulla diagonale. In questo modo si ottiene $P^{n-1} \cdots P^1 A = R$ dove tutte le $P$ sono ortogonali ed il loro pordotto vale $Q^T$, quindi $Q^T A = R\; \rightarrow \; A=QR$.

[x,A]=QRhouse(A,b) - Data la matrice $A$ ed il vettore $\underline{b}$ la funzione fattorizza $A$ come il prodotto di due matrici $QR$ con $Q$ ortogonale ed $R$ triangolare superiore, secondo il metodo di Householder e risolve il sistema $R\underline{x}=Q^{T}\underline{b}$. La matrice $R$ viene memorizzata nella parte triangolare superiore di $A$, ed al posto della $Q$ vengono memorizzate nella parte inferiore le parti significative dei vettori di $Householder$ associati a ciascuna colonna della matrice originaria.

%QRHOUSE
%[x,A]=QRHOUSE(A,b)
%
%Pre: A è una matrice quadrata non singolare n x n; 
% b vettore di lunghezza n
%
% La procedura fornisce la fattorizzazione QR di una matrice 
% A e risolve
%  il sistema lineare Ax=b;  Q è una matrice 
% ortogonale (Q'Q=I) e R è
% triangolare superiore.
% In x viene restituita la soluzione del sistema lineare mentre
% in A viene 
%restituita una matrice che contiene nella sua 
% parte triangolare superiore
% la matrice R e nella parte 
% strettamente triangolare inferiore le componenti
%  essenziali
% dei vettori di householder che generano la matrice Q
% See also HOUSE
function [x,A]=QRhouse(A,b)
n=length(b);
% beta = 2/(v' * v) con v vettore di houseolder
beta=zeros(1,n-1);
gamma=zeros(1,n);

for i=1:n-1
   % adesso il vettore v si trova in
   % A(i:n,i)
   A(i:n,i)=house(A(i:n,i))';
   beta(i)=2/([1;A(i+1:n,i)]'*[1;A(i+1:n,i)]);
   % calcola v'*y con y vettore colonna corrente
   for j=i+1:n 
      gamma(j)=([1;A(i+1:n,i)]'*A(i:n,j));
   end
   % aggiorna la riga i-esima
   for t=i+1:n
      A(i,t)=A(i,t)-gamma(t)*beta(i);
   end
   % viene aggiornata la sottomatrice in
   % esame per colonne, moltiplicando per 
   % una matrice di Householder
   for j=i+1:n
      % prodotto vettoriale con fattore moltiplicativo
      % aggiuntivo -gamma(k)*beta(i)
      for k=i+1:n
         A(k,j)=A(k,j)-gamma(j)*A(k,i)*beta(i);
      end
   end
end
% aggiornamento termine noto
x=b(:);
for i=1:n-1
   temp=beta(i)*[1;A(i+1:n,i)]'*x(i:n);
   x(i:n)=x(i:n)-temp*[1;A(i+1:n,i)];
end
x=solveU(A,x);
return


   
   
   
   

>> A

A =

3.0000000000 1.0000000000 0.5000000000 0.8000000000
1.5000000000 9.0000000000 0.5000000000 0.6000000000
3.2000000000 4.0000000000 12.7000000000 3.2000000000
0.7000000000 3.2000000000 4.1000000000 11.3000000000

>> b

b =

   0.50000000000000
   1.20000000000000
   3.00000000000000
   6.20000000000000

>>  inv(A)*b

ans =

  -0.01092762955907
   0.09772671951632
   0.08447073456824
   0.49102600234596

>> QRHouse(A,b)

ans =

  -0.01092762955907
   0.09772671951632
   0.08447073456824
   0.49102600234596