Fattorizzazione $A=QR$

Quello che facciamo adesso è cercare una fattorizzazione di $A$ come un fattore ortogonale ed un fattore triangolare superiore.

Data $A \in {\mathbb{R}}^{m \times n} \mbox{ } m \geq n$ (in generale ha più righe che colonne) e richiediamo che abbia rango massimo, cioè $rank(A)=n$ (e quindi se la matrice è quadrata, sarà anche non singolare). Allora ciò che cerchiamo è

$$A=QR \quad Q \in {\mathbb{R}}^{m \times m} \mbox{, } R \in {\mathbb{R}}^{m \times n}$$

$$R= \left( \begin{array}{c} R_1 \\ \hline \mathcal{O} \en... ... \in {\mathbb{R}}^{n \times n} \mbox{ triangolare superiore}$$

Come osservazione possiamo notare che se $m=n$, $R$ è triangolare superiore ed $R=R_1$.

Mostreremo l'esistenza di questa fattorizzazione mediante una dimostrazione costruttiva che ci guiderà nella costruzione dell'algoritmo.

Dalle proprietà del rango abbiamo che

$$n=rank(A)=rank(QR)=rank(R)$$

visto che $Q$ ha rango massimo; $R$ ha rango massimo, e per costruzione

$$rank(R)=rank(R_1)$$

ed essendo $R_1 \in {\mathbb{R}}^{n \times n}$ questo implica che $R_1$ è non singolare.

Dunque, come già visto durante l'algoritmo di eliminazione di Gauss, dato un vettore $x$ cerchiamo una matrice ortogonale $P \in {\mathbb{R}}^{n \times n}$ tale che:

$$Px=\left( \begin{array}{c} \alpha \\ 0\\ \vdots \\ 0 \end{array} \right) = \alpha e_1\mbox{,}$$

che azzeri tutti gli elementi di $x$ da un certo punto in poi. Inoltre, prendendo il quadrato della norma euclidea di $Px$ si ottiene:

$$\begin{array}{l} x^TP^TPx=x^Tx={\Vert x\Vert}_2^2 \\ \mb... ...e_1= {\alpha}^2 {\Vert e_1\Vert}_2^2= {\alpha}^2 \end{array}$$

da cui otteniamo $\alpha= \pm {\Vert x\Vert}_2$.

Cerchiamo dunque una matrice $P$ nella forma

$$P=I-2\frac{vv^T}{v^Tv} \qquad v \in {\mathbb{R}}^n \mbox{ assegnato}$$

con $v$ vettore da determinare; notiamo inoltre che $\frac{2}{v^Tv}$ è uno scalare. $P$ è una matrice simmetrica, perché sia $I$ che $vv^T$ lo sono ( $vv^T={(vv^T)^T}={(v^T)}^Tv^T=vv^T$) ed è detta matrice di Householder; sapendo che $P=P^T$ verifichiamo l'ortogonalità di P:

$$\begin{eqnarray*} P^2 & = & P^TP=I \\ & = & \left( I-2\frac{vv^T}{v^Tv} \rig... ...T}{{(v^Tv)}^2} = I - 4\frac{vv^T}{v^Tv}+ 4\frac{vv^T}{v^Tv}=I \end{eqnarray*}$$

Come per la matrice elementare di Gauss, anche $P$ è ottenuta come la matrice identità a cui viene sommata una matrice di correzione di rango $1$, e sommare matrici di rango basso significa eseguire pochi conti.

Eravamo rimasti che il vettore $v$ doveva essere determinato, verifichiamo allora che

$$v=x-\alpha e_1$$

è proprio ciò di cui abbiamo bisogno, cioè che:

$$\begin{eqnarray*} \left( I-2\frac{vv^T}{v^Tv} \right) x & = & \alpha e_1 \\ \... ...ac{v^Tx}{v^Tv}\right)+ \alpha 2\frac{v^Tx}{v^Tv} e_1 \mbox{.} \end{eqnarray*}$$

Se si riuscisse a dimostrare che il coefficiente del vettore $x$ è pari a zero, e quello del vettore $e_1$ è pari ad $\alpha$ avremmo terminato; per fare questo bisogna dimostrare che

$$2\frac{v^Tx}{v^Tv}=1 \qquad \quad 2v^Tx=v^Tv$$

$$\begin{eqnarray*} 2v^Tx & = & \left( \begin{array}{l} v=x - \alpha e_1 \\ \... ...overbrace{e_1^Te_1}^{1} = \\ & = & 2{\alpha}^2 - 2 \alpha x_1 \end{eqnarray*}$$

Le due quantità sono uguali e quindi abbiamo verificato che una matrice $P$ così costruita soddisfa le nostre richieste; il vettore $v$ visto è detto vettore di Householder.

La scelta del segno di $\alpha$, cioè se scegliere $\alpha = + {\Vert x\Vert}_2$ oppure $\alpha = - {\Vert x\Vert}_2$, è sottesa alla riduzione degli errori nell'aritmetica finita del calcolatore (la matrice ortogonale non pone grossi problemi di calcolo poichè l'imposizione che la somma dei quadrati sia pari ad uno la rende tale che gli elementi siano ben limitati in modulo e quindi che anche gli errori sono limitati). Come detto

$$v=x- \alpha e_1 = \left( \begin{array}{c} x_1- \alpha \\ x_2 \\ \vdots \\ x_n \end{array} \right)$$

è quindi opportuno che $x_1$ ed $\alpha$ siano di segno concorde così da eliminare il problema della cancellazione (soprattutto nel caso $x_1$ fosse l'elemento di modulo massimo), dunque

$$\begin{array}{rcl} \mbox{se }x_1 \geq 0 & , & \alpha = - {... ...se }x_1 < 0 & , & \alpha = +{\Vert x\Vert}_2 \\ \end{array}$$

Siamo ora in grado di dimostrare il teorema di esistenza della fattorizzazione $QR$

Teorema 1   Se $A$ è una matrice ${\mathbb{R}}^{m \times n}$, di rango massimo, allora $\exists Q \in {\mathbb{R}}^{m \times m}$ ortogonale ed $\exists R= \left( \begin{array}{c} R_1 \\ \hline \mathcal{O} \end{array} \right) \in {\mathbb{R}}^{m \times n}, R_1 \in {\mathbb{R}}^{n \times n}$ triangolare superiore e non singolare tali che $A=QR$

La dimostrazione di questo teorema definirà un algoritmo che ci consentirà di generare la fattorizzazione. Ragioniamo come abbiamo fatto con l'eliminazione di Gauss: quindi poniamo

$$A= \left( \begin{array}{ccc} a_{11}^{(0)} & \cdots & a_{1... ...cdots & a_{mn}^{(0)} \\ \end{array} \right) \equiv A^{(0)}$$

considerando adesso

$$x=\left( \begin{array}{c} a_{11}^{(0)} \\ \vdots \\ a_{m1}^{(0)} \end{array} \right)$$

possiamo costruire la matrice $P_1$ tale che

$$P_1\left( \begin{array}{c} a_{11}^{(0)} \\ \vdots \\ a_... ... a_{11}^{(1)} \\ 0 \\ \vdots \\ 0 \end{array} \right)$$

con $a_{11}^{(0)} \neq a_{11}^{(1)}$; in questo caso non abbiamo bisogno di pivoting come con l'algoritmo di Gauss: per come è costruito il vettore di Householder questo sarà sempre diverso dal vettore nullo ( $x- \alpha e_1 = \underline{0} \Longleftrightarrow x= \underline{0}$) e dunque la matrice $P$ sarà sempre ben definita.

Ottenuta la matrice $P_1$ possiamo operare come già esaminato con Gauss

$$P_1 A^{(0)} = \left( \begin{array}{cccc} a_{11}^{(1)} & \cd... ... & \cdots & a_{mn}^{(1)} \end{array} \right) \equiv A^{(1)}$$

come si può notare gli elementi sottodiagonali della prima colonna vengono annullati ma gli altri vengono comunque modificati, anche quelli della prima riga.

Consideriamo la seconda colonna, dall'elemento diagonale in poi e definiamo la matrice $P^{(2)}$ in modo tale che

$$P^{(2)}\left( \begin{array}{c} a_{22}^{(1)} \\ \vdots \\ ... ...ay} \right) \qquad P^{(2)} \in {\mathbb{R}}^{m-1 \times m-1}$$

possiamo adesso definire la matrice $P_2$ come

$$P_2 = \left( \begin{array}{c\vert c} 1 & \\ \hline & ... ...} \end{array} \right) \qquad P_2 \in {\mathbb{R}}^{m \times m}$$

ed è ancora simmetrica ed ortogonale. Otteniamo infine

$$P_2P_1 A = \left( \begin{array}{ccccc} a_{11}^{(1)} & \cdot... ...s & a_{mn}^{(2)} \end{array} \right) \equiv A^{(2)}\mbox{.}$$

La prima riga non viene modificata poichè la prima riga di $P_2$ è la prima riga dell'identità, vengono modificate solo le $m-1$ righe sotto la prima.

In generale, al passo i-esimo

$$P^{(i)}\left( \begin{array}{c} a_{ii}^{(i-1)} \\ \vdots \... ... \right) \qquad P^{(i)} \in {\mathbb{R}}^{m-i+1 \times m-i+1}$$

$$P_i = \left( \begin{array}{c\vert c} I_{i-1} & \\ \hlin... ...} \end{array} \right) \qquad P_i \in {\mathbb{R}}^{m \times m}$$

ortogonale e simmetrica. Alla fine, dopo n passi, otteniamo:

$$\underbrace{P_n\cdots P_2P_1}_{Q^T}A= \left( \begin{array}... ...1 \\ \hline \mathcal{O} \end{array} \right) \equiv R \mbox{.}$$

Il prodotto di matrici ortogonali è ancora una matrice ortogonale, da cui

$$\begin{array}{l} Q^TQ=(P_n \cdots P_1)(P_1 \cdots P_n)=I \\ \\ QQ^TA=QR \\ \\ A=QR \end{array}$$

Nel caso di matrici quadrate, se la matrice in esame è non singolare, allora è anche fattorizzabile $QR$, mentre per la fattorizzazione $LU$, come ci ricordiamo, ci sono ben altri vincoli, molto più stringenti che la semplice non singolarità.

Soffermiamoci sullo spazio di memoria di cui abbiamo bisogno. Ci riferiremo al primo passo, dal momento che gli altri sono simili. Per conoscere $P_1$ l'unica cosa di cui abbiamo bisogno è il vettore $v$, che ha dimensioni $m$: infatti noto $v$, possiamo ricavare $P_1$.

Dato che $rank(A)=n$ è massimo

$$\begin{array}{l} x\neq \underline{0} \qquad {\alpha}^2={\V... ... x_2 \\ \vdots\\ x_m \end{array} \right) \end{array}$$

e sicuramente $v_1 \neq 0$ poichè $x_1$ ed $\alpha$ sono di segno concorde. Possiamo dunque definire

$$\tilde{v}=\frac{v}{v_1} \qquad {\tilde{v}}_1 \neq 0$$

e quindi possiamo evitare di memorizzare la prima componente del vettore $\tilde{v}$ in quanto sappiamo che è pari ad uno. Vediamo cosa comporta questa scelta di vettore per la matrice di Householder:

$$\tilde{P}=I-2\frac{\tilde{v}{\tilde{v}}^T}{{\tilde{v}}^T\tild... ...^T(\frac{v}{v_1})}= I-2\frac{v_1^2}{v_1^2} \frac{vv^T}{v^Tv}=P$$

la matrice di Householder è perciò invariante per moltiplicazioni per scalari. In conclusione, invece di dover utilizzare $m-i$ componenti, ne abbiamo bisogno solo di $m-i-1$, proprio quelle che andiamo ad azzerare; possiamo riscrivere la matrice $A$ con tutte e sole le informazioni necessarie alla fattorizzazione: nella parte triangolare superiore metteremo $R_1$ mentre nelle parti che andiamo ad azzerare metteremo i vettori $\tilde{v}$ a partire dalla seconda componente in poi.