Calcolo numerico/Fattorizzazione QR

Motivazioni

$Q$ è una matrice unitaria, cioè $Q{Q}^h=Q^{h}Q=1$, allora $Q^{-1}=Q^h$, e $R$ è una triangolare superiore.

Ci chiediamo per quale motivo serve un'alternativa alla fattorizzazione $PA=LU$.

Il vantaggio della fattorizzazione $QR$ è una maggior stabilità rispetto alla $PA=LU$. Infatti la $QR$ è una fattorizzazione con matrici elementari unitarie, che conservano la norma 2.

Lo svantaggio è il raddoppio del costo computazionale.

Non occorre la tecnica del pivot.

Ci sono situazioni più sofisticate in cui si usa una tecnica di pivot per colonne, per garantire che gli elementi diagonali della $R$ abbiano una certa proprietà di ordinamento.

Questa fattorizzazione si usa solo in casi limite, ad esempio nel calcolo degli autovalori (spettro completo) e del rango, e della s.v.d.

Data una matrice $A$ hermitiana, allora si può diagonalizzare con trasformazioni unitarie, scrivendo $Q\lambda Q^h$.

La decomposizione in fattori singolari è l'estensione di quest'idea anche ai casi in cui la matrice non è simmetrica e non è quadrata.

La fattorizzazione $QR$ esiste sempre ma non è unica.

Risoluzione del sistema

Con la fattorizzazione $QR$ $A$ viene fattorizzata nel prodotto di una matrice unitaria e di una triangolare superiore.

$Q$ è unitaria quindi $Q^{-1}=Q^h$.

Nel sistema $A𝒙=𝒃$, cioè $QR𝒙=𝒃$, poniamo $R𝒙=𝒚$, risolvo prima

$Q𝒚=𝒃⟶𝒚=Q^{-1}𝒃=Q^h𝒃\text{prodotto matrice per vettore}$

quindi si risolve

$R𝒙=𝒚$

con $R$ matrice triangolare.

Il costo di risoluzione è $2/3n^3$.

Matrici elementari di Hausolder

La matrice della diade ha rango uguale a 1.

Nella fattorizzazione $QR$ consideriamo le matrici elementari di Hausolder (o matrici di riflessione). Queste matrici sono della forma:

$P=I-\beta 𝒗{𝒗}^h,\beta \in ℝ,𝒗\in {ℂ}^n$

Proposizione 2.2

Le matrici di riflessione sono hermitiane, e sono unitarie se scegliamo $\beta =\frac{2}{|𝒗{|}_2}$.

Dim.

Hermitiane:

$P^h=I-\beta (𝒗{𝒗}^h{)}^h=I-\beta 𝒗{𝒗}^h=P$

Unitarie: sappiamo che $I=P{P}^h$, e per hermitianità

$P{P}^h=PP=(I-\beta 𝒗{𝒗}^h)(I-\beta 𝒗{𝒗}^h)$

$=I-\beta 𝒗{𝒗}^h-\beta 𝒗{𝒗}^h+{\beta }^2𝒗{𝒗}^h𝒗{𝒗}^h$

$=I-2\beta 𝒗{𝒗}^h+{\beta }^2𝒗(|𝒗|{)}^2{𝒗}^h$

(infatti ${𝒗}^h𝒗=(|𝒗{|}_2{)}^2$)

quindi

$I=I-2\beta 𝒗{𝒗}^h+\beta (|𝒗{|}_2{)}^2𝒗{𝒗}^h$

implica

$\beta 𝒗{𝒗}^h\ast [-2+\beta |𝒗{|}_2^2]=0$

cioè, se $\beta \ne 0$, deve valere $-2+\beta |𝒗{|}_2^2=0$, cioè

$\beta =\frac{2}{(|𝒗{|}_2{)}^2}$

cvd

Altre proprietà sulle matrici di Hausolder

Sia $P=P(𝒗)$ (la matrice di Hausolder è univocamente determinata da $𝒗$), con $P=I-\beta 𝒗{𝒗}^h$, e $\beta =\frac{2}{(|𝒗{|}_2{)}^2}$. Allora per ogni $𝒙\in {𝒗}_{\perp }$, la matrice di Hausolder si comporta come l'identità. Infatti:

$𝒙(I-\beta 𝒗{𝒗}^h)=𝒙-\beta 𝒗{𝒗}^h𝒙=𝒙$

infatti $𝒙𝒗=0$ per $𝒙\in {𝒗}_{\perp }$.

$P$ applicata nel sottospazio $\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\{𝒗\}$ è:

$(I-\beta 𝒗{𝒗}^h)(\gamma 𝒗)=\gamma 𝒗-\beta \gamma 𝒗{𝒗}^h𝒗$

$=\gamma 𝒗-\beta \gamma |𝒗{|}^2𝒗$

e per definizione di $\beta$:

$=\gamma 𝒗-\frac{2}{(|𝒗{|}_2{)}^2}\gamma |𝒗{|}^2𝒗$

$=\gamma 𝒗-2\gamma 𝒗=-\gamma 𝒗$

cioè si ha una riflessione.

Esistenza della matrice elementare P

Teorema 2.10

Per ogni $𝒙\in {ℂ}^n\setminus 0$ esiste $P\in {ℂ}_{n\times n}$ matrice elementare di Hausolder tale che $P𝒙=\alpha {𝒆}_1$.

Dimostrazione

Dobbiamo determinare il parametro $\alpha \in ℂ$ e la matrice elementare di Hausolder $P$.

Siccome $\alpha \in ℂ$, $\alpha =|\alpha |\ast e^{i\theta }$ (relazione $\circ$). Deve essere $P𝒙=\alpha {𝒆}_1$, allora

$|P𝒙|=|\alpha ||{𝒆}_1{|}_2=|\alpha |$

Inoltre $P$ è un'isometria quindi conserva la norma 2, cioè:

$(|P𝒙{|}_2{)}^2=P^2{𝒙}^h𝒙={𝒙}^h𝒙=(|𝒙{|}_2{)}^2$

allora, eguagliando le due espressioni trovate per $(|P𝒙{|}_2{)}^2$, ottengo:

$|\alpha |=|𝒙{|}_2$

Per determinare l'angolo $\theta$, dev'essere

${𝒙}^{h}P𝒙={𝒙}^h\alpha {𝒆}_1=\alpha {𝒙}^h{𝒆}_1=\alpha {\overline{x}}_1\text{ è uno scalare}$

Per hermitianità di $P$:

$\overline{(}{𝒙}^{h}P𝒙)=({𝒙}^{h}P𝒙{)}^h={𝒙}^h{P}^h𝒙={𝒙}^{h}P𝒙$

allora siccome lo scalare è uguale al suo coniugato, è reale, quindi $\alpha {\overline{x}}_1\in ℝ$.

${\overline{x}}_1\in ℂ$, allora possiamo scriverlo come $|x_1|e^{i\varphi }$ (relazione $\circ \circ$).

$\alpha {\overline{x}}_1=|\alpha |e^{i\theta }|x_1|e^{i\varphi }$

e sostituendo l'espressione di $\alpha$:

$=|𝒙{|}_2|x_1|e^{i(\theta -\varphi )}$

$\alpha {\overline{x}}_1\in ℝ$, allora

$\theta -\varphi =k\pi$

ovvero $\theta =\varphi +\delta \pi$ con $\delta =0,1$.

In conclusione:

$\alpha =|\alpha |e^{i\theta }=|𝒙{|}_2\ast e^{i(\varphi +\delta \pi )}$

$e^{i\delta \pi }=\pm 1⟶\alpha =\pm |𝒙{|}_2{e}^{i\varphi }$

e per la $\circ \circ$,

$=\pm |𝒙{|}_2\ast \frac{x_1}{|x_1|}$

Se $x_1=0$, allora poniamo $\alpha =-|𝒙{|}_2$.

Vale che

$P𝒙=(I-\beta 𝒗{𝒗}^h)(𝒙)=𝒙-\beta 𝒗{𝒗}^h𝒙,$

e per definizione di $\beta$ otteniamo:

$P𝒙=𝒙-\frac{2}{(|𝒗{|}_2{)}^2}{𝒗}^h𝒙𝒗$

Chiamiamo $\hat{\beta }$ il coefficiente che moltiplica $𝒗$, e mettiamo in evidenza $\hat{\beta }𝒗$:

$\hat{\beta }𝒗=𝒙-P𝒙=𝒙-\alpha {𝒆}_1$

Osservazione 2.5

A patto di escludere che $\hat{\beta }=0$, è lecito porre $\hat{\beta }=1$, in quanto

$P(\gamma 𝒗)=I-\frac{2}{|\gamma 𝒗{|}^2}\ast (\gamma 𝒗)\ast (\gamma {𝒗}^h)$

$=I-\frac{2}{|\gamma {|}^2|𝒗{|}_2}|\gamma {|}^2\ast 𝒗{𝒗}^h$

$=I-\frac{2}{|𝒗{|}_2}\ast 𝒗{𝒗}_h=P(𝒗)$

allora per il fatto che $P(\gamma 𝒗)=P(𝒗)$, possiamo porre $\hat{\beta }=1$.

Tornando alla formula:

$\hat{\beta }𝒗=𝒙-\alpha {𝒆}_1$

e ponendo $\hat{\beta }=1$, otteniamo:

$𝒗=𝒙-\alpha {𝒆}_1$

Nota bene: possiamo escludere il caso $\hat{\beta }=0$, infatti, per definizione, $\hat{\beta }=0$ se e solo se $𝒙\in {𝒗}_{\perp }$, e possiamo scegliere il vettore $𝒙$ al di fuori del complemento ortogonale.

Calcoliamo:

${𝒗}^h𝒗={𝒗}^h(𝒙-\alpha {𝒆}_1)=(|𝒙|{)}^2+|x_1||𝒙{|}_2$

e la quantità è positiva se $𝒙\ne 0$ come nell'ipotesi.

cvd

Riassumendo, poniamo

$\alpha =\{\begin{array}{cc}\pm |𝒙{|}_2\frac{x_1}{|x_1|}& \text{se}{x}_1\ne 0\\ -|𝒙{|}_2& \text{se}{x}_1=0\end{array}$

$𝒗=𝒙-\alpha {𝒆}_1$

${𝒆}_1$ modifica solo la prima componente del vettore $𝒗$. Quindi:

$v_1=x_1-\alpha ,=x_1-(-|𝒙{|}_2\frac{x_1}{|x_1|})=$

$=x_1+|𝒙{|}_2\frac{x_1}{|x_1|})=$

$=x_1\ast (\frac{|x_1|+|𝒙{|}_2}{|x_1|})$

e non c'è possibilità di errori di cancellazione.

Se $𝒙\in {ℝ}^n$, $\alpha \in ℝ,𝒗\in {ℝ}^n$ per costruzione, allora $P$ è una simmetrica ortogonale.

Il calcolo di $𝒗$ costa $n$ operazioni per il calcolo della norma di $𝒙$, 2 operazioni moltiplicative, un'estrazione di radice, cioè un numero dell'ordine di $n$ operazioni.

Algoritmo della fattorizzazione QR

Utilizziamo il teorema dimostrato per determinare la fattorizzazione $QR$.

Supponiamo di essere al passo $k$-esimo della fattorizzazione $QR$. Allora chiamiamo ${𝒂}_k=(A_{kk},\dots ,A_{nk})$ la $k$-esima sottocolonna della matrice $A^{(k)}$, che dobbiamo azzerare (corrisponde al vettore $𝒙$ del teorema). Scegliamo il vettore ${𝒗}_k=𝒗$ che ha 0 nelle prime $k-1$ componenti, e le componenti successive sono:

${𝒗}_k={\gamma }_k(|{𝒂}_k{|}_2+|A_{kk}^{(k)}|)=\mathrm{s}\mathrm{g}\mathrm{n}{a}_{kk}\ast (|{\sigma }_k|)$

${𝒗}_{k+1}=a_{k+1}^{(k)}$

${𝒗}_{k+2}=a_{k+2}^{(k)}$

$\dots$

${𝒗}_n=a_n^{(k)}$

dove definiamo

${\gamma }_k=\{\begin{array}{cc}\frac{A_{kk}}{|A_{kk}|}& ⟺A_{kk}\ne 0\\ 1& \text{se}{A}_{kk}=0\end{array}$

(per la definizione di ${𝒗}_k$ sfruttiamo l'espressione ricavata nel teorema, cioè ${𝒗}_k=𝒙-\alpha {𝒆}_1$, e sfruttiamo anche l'espressione ricavata per $v_1$)

${\beta }_k=\frac{2}{(|{𝒗}_k{|}_2{)}^2}$

Allora calcoliamo:

$|{𝒗}_k{|}^2=|{\gamma }_k{|}^2(|{𝒂}_k{|}_2+|A_{kk}|{)}^2+|{𝒂}_k{|}_2^2-|A_{kk}{|}^2$

$=(|{𝒂}_k{|}_2{)}^2+2|A_{kk}||{𝒂}_k{|}_2+(|A_{kk}|{)}^2+|{𝒂}_k{|}_2^2-|A_{kk}{|}^2$

$=2(|{𝒂}_k{|}_2{)}^2+2|A_{kk}||{𝒂}_k{|}_2$

e raccogliendo $2|{𝒂}_k{|}_2$ otteniamo:

$=2|{𝒂}_k{|}_2\ast [|{𝒂}_k{|}_2+|A_{kk}|]$

cioè

${\beta }_k=\frac{2}{2(|{𝒗}_k{|}_2{)}^2}$

$=\frac{2}{2|{𝒂}_k{|}_2\ast [|{𝒂}_k|+A_{kk}]}$

${\beta }_k=\frac{1}{|{𝒂}_k{|}_2\ast [|{𝒂}_k|+A_{kk}]}$

Allora abbiamo determinato univocamente la matrice

$P_k=I-\beta 𝒗{𝒗}^h$

Se il vettore ${𝒂}_k$ è tutto nullo, consideriamo ${\gamma }_k=1$ e ${𝒗}_k=0$, in modo che a quel passo la matrice di Hausolder sia l'identità.

Applicando la matrice $P$ trovata al vettore $𝒙$ cioè alla colonna ${𝒂}_k$, tutte le entrate di ${𝒂}_k$ da $k+1$ a $n$ vengono azzerate, mentre si ha $v_k=A_{kk}=-|𝒙{|}_2$ (infatti $P𝒙=\alpha {𝒆}_1$). Applichiamo poi la matrice di Hausolder alle colonne rimanenti di $A$.

Riepilogando, data la matrice $A$ di partenza, con il primo passo $P_{1}A$ azzero la prima colonna e modifico $A_{11}$, e poi ripeto il procedimento, per $n-1$ passi, e otteniamo la matrice triangolare superiore $R$.

Il prodotto di matrici unitarie è unitario, e sappiamo che

$R=P_{n-1}\dots P_2{P}_{1}A$

e per hermitianità:

$A=P_1{P}_2\dots P_{n-1}R$

quindi

$Q=P_1{P}_2\dots P_{n-1}$

Non unicità della fattorizzazione

Proposizione 2.3

Per ogni $A$ esiste la fattorizzazione $A=QR$, ma non è unica.

Dim. Consideriamo una matrice di fase $S$, cioè una matrice diagonale con elementi ${\theta }_1,\dots ,{\theta }_n$ con $|{\theta }_i|=1\mathrm{\forall }i$, ed è una matrice unitaria. Allora $S{S}_h=I$, cioè $A=QR=QS{S}_{h}R$, e definiamo $\overline{Q},\overline{R}$ dove $\overline{Q}=QS$ e $\overline{R}=S_{h}R$, dove $\overline{Q}$ è ancora unitaria e $\overline{R}$ è triangolare superiore. Allora abbiamo comunque una fattorizzazione $QR$ della matrice.

cvd

Questa proprietà può essere sfruttata per scegliere la $\overline{R}$ in modo che gli elementi sulla diagonale siano ordinati in una certa maniera.

Fattorizzazione e risoluzione dei sistemi

Per risolvere il sistema lineare, non è necessario trovare esplicitamente la matrice $Q$, e nemmeno le $P_h$.

Passo 1: poniamo $A_1=A,{𝒃}_1=𝒃$, e chiamiamo $T_1$ la matrice che accosta ${𝒃}_1$.

Passo 2: chiamiamo $T_2=P_1{T}_1$, con

$P_1=I-{\beta }_1{𝒗}_1{𝒗}_1^h$

Chiamiamo poi ${𝒚}_1^h$ il vettore riga ${𝒗}_1^h{T}_1$, in modo da scrivere:

$T_2=P_1{T}_1=T_1-{\beta }_1{𝒗}_1{𝒚}_1^h$

Passo $k+1$:

$T_{k+1}=P_k{T}_k=T_k-{\beta }_k{𝒗}_k{𝒚}_k^h$

con

${𝒚}_k^h={𝒗}_k^h{T}_k$

Infine

Passo $n$:

$T_n=R⟶{𝒚}_k^h={𝒗}_{k}R$

il vettore $𝒚$ è stato calcolato senza calcolare esplicitamente le matrici $Q$, ed è tale che $q𝒚=𝒃$.

Costo computazionale

Per analizzare il costo del $k$-esimo passo, tengo conto che ${𝒗}_k$ ha $k-1$ componenti nulle e $n-k+1$ componenti diverse da 0. Per calcolare ${𝒗}_k$, le operazioni sono quindi $n-k+1+2=n-k+3$.

$\dots \dots \dots$

Stabilità della fattorizzazione QR

Questa non è una norma matriciale indotta, infatti

$|I{|}_f=\sqrt{n}$

ed è diversa da 1, mentre con le norme matriciali indotte l'identità a norma 1.

Teorema 2.11 (analisi dell'errore all'indietro)

Sia $\overline{𝒙}$ la soluzione effettivamente calcolata con l'applicazione delle ${\overline{P}}_k$ effettive alla matrice $A$ per calcolare la $\overline{R}$. Sia $\overline{𝒚}$ il vettore che compare nella risoluzione backward del sistema $R𝒙=𝒚$. Allora $\overline{𝒙}$ è soluzione esatta di un sistema

$(A+{\delta }_A)\overline{𝒙}=𝒃+{\delta }_b$

dove

$|{\delta }_A{|}_f<u\gamma n^2|A{|}_f+n|\overline{R}{|}_f+o(u^2)$

$|{\delta }_b{|}_2<u\gamma n^2|𝒃{|}_2+O(u^2)$

L'entità della perturbazione dipende dai dati iniziali, dalle dimensioni della matrice, e dalla matrice di fattorizzazione. Non dipende dalla norma di Frobenius di $Q$ che è 1.

Consideriamo le matrici $A^{(k)}$ al $k$-esimo passo di fattorizzazione, e calcolo il prodotto

$A_{k+1}^h{A}_{k+1}=A_k^h{P}_k^h{P}_k{A}_k$

siccome $P_k^h{P}_k=I$ allora $A_k^h{A}_k$ non varia con i passi di fattorizzazione.

Allora

$\mathrm{t}\mathrm{r}(A_{k+1}^h{A}_{k+1})=\mathrm{t}\mathrm{r}(A_k^h{A}_k)$

$⟶|A_{k+1}{|}_f^2=|A_k{|}_f^2$

allora

$|A_k{|}_f=|A{|}_f$

inoltre

$|R{|}_f=|A_n{|}_f=|A{|}_f$

allora non c'è possibilità di esplosione dei coefficienti.

La $QR$ è stabile e può essere fatta senza pivot.