Calcolo numerico/Fattorizzazione di Cholesky

Esistenza della fattorizzazione $A=R^{t}R$

Questo tipo di fattorizzazione si applica solo al caso di matrici definite positive.

Teorema 2.9

Supponiamo di avere $A\in \mathbb {R} ^{n\times n}$ simmetrica e non singolare. $A$ è definita positiva se e solo se esiste $R\in \mathbb {R} ^{n\times n}$ triangolare superiore tale che $A=R^{t}R$.

Dimostrazione

$2\longrightarrow 1$: $A$ è non singolare, e supponiamo $A=R^{t}R$, allora $R$ è non singolare infatti

$\det A=\det R^{t}*\det R=|\det R|^{2}$

e quindi $\det A\neq 0$ implica $R$ non singolare. Allora

$\forall {\boldsymbol {x}}\neq 0,\quad {\boldsymbol {y}}=R{\boldsymbol {x}}\neq 0,$

quindi

$\vert {\boldsymbol {y}}\vert _{2}^{2}\geq 0$

e siccome ho escluso il caso nullo,

$\vert {\boldsymbol {y}}\vert _{2}^{2}>0.$

e $A$ è definita positiva.

$1\longrightarrow 2$: consideriamo $A$ simmetrica e definita positiva:

$\left({\begin{array}{cc}\alpha &{\boldsymbol {a}}^{t}\\{\boldsymbol {a}}&A_{\ast }\end{array}}\right)$

dove $A_{\ast }$ è una sottomatrice di dimensione $n-1\times n-1$ e a ${\boldsymbol {a}}$ è un vettore di $n-1$ componenti.

Definiamo $R,R^{t}$ nel seguente modo:

$R=\left({\begin{array}{cc}\rho &{\boldsymbol {r}}^{t}\\{\boldsymbol {0}}&R_{\ast }\end{array}}\right)$

$R^{t}=\left({\begin{array}{cc}\rho &{\boldsymbol {0}}^{t}\\{\boldsymbol {r}}&(R_{\ast })^{t}\end{array}}\right)$

Il prodotto delle due matrici è:

$R^{t}R=\left({\begin{array}{cc}\rho ^{2}&\rho *{\boldsymbol {r}}^{t}\\\rho *{\boldsymbol {r}}&{\boldsymbol {r}}{\boldsymbol {r}}^{t}+R_{\ast }^{t}R_{\ast }\end{array}}\right)$

ed eguagliamo termine a termine la decomposizione a blocchi con il risultato del prodotto ottengo:

${\begin{aligned}\alpha =\rho ^{2}&\longrightarrow \rho ={\sqrt {\alpha }}\\{\boldsymbol {a}}=\rho {\boldsymbol {r}}&{\boldsymbol {r}}={\frac {1}{\sqrt {\alpha }}}{\boldsymbol {a}}\\{\boldsymbol {a}}^{t}=\rho {\boldsymbol {r}}^{t}&\longrightarrow {\boldsymbol {r}}^{t}={\frac {{\boldsymbol {a}}^{t}}{\sqrt {\alpha }}}\\A_{\ast }={\boldsymbol {r}}{\boldsymbol {r}}^{t}+R_{\ast }^{t}R_{\ast }&\longrightarrow R_{\ast }^{t}R_{\ast }=A_{\ast }-{\boldsymbol {r}}{\boldsymbol {r}}^{t}=A_{\ast }-1/\alpha {\boldsymbol {a}}{\boldsymbol {a}}^{t}\end{aligned}}$

Definiamo la quantità

${\hat {A}}_{\ast }=R_{\ast }^{t}R_{\ast }$

Se ${\hat {A}}_{\ast }$ è ancora simmetrica e definita positiva, si può iterare la procedura, lavorando su $R_{\ast }^{t}R_{\ast }$, fino a ottenere la fattorizzazione con $R$ triangolare superiore. Quindi verifichiamo che ${\hat {A}}_{\ast }$ è simmetrica definita positiva.

Simmetria: sappiamo che

${\hat {A}}_{\ast }^{t}=(R_{\ast }^{t}R_{\ast })^{t}=R_{\ast }^{t}(R_{\ast }^{t})^{t}=R_{\ast }^{t}R_{\ast }={\hat {A}}_{\ast }$

Positività: appliciamo ad $A$ una matrice elementare $M$ per azzerare la prima sottocolonna, con

$A=\left({\begin{array}{cc}\alpha &{\boldsymbol {a}}^{t}\\{\boldsymbol {a}}&A_{\ast }\end{array}}\right)$

$M=\left({\begin{array}{cc}1&0\\-{\boldsymbol {a}}/\alpha &I\end{array}}\right)$

Il risultato del prodotto è:

$MA=\left({\begin{array}{cc}\alpha &{\boldsymbol {a}}^{t}\\0&-1/\alpha {\boldsymbol {a}}{\boldsymbol {a}}^{t}+A_{\ast }\end{array}}\right)$

(come nell'eliminazione di Gauss ordinaria, la prima riga rimane invariata e la prima colonna viene azzerata)

Per poter azzerare una sottoriga, applichiamo alla matrice $MA$ un passo di eliminazione gaussiana per righe. Eseguiamo il prodotto $MA{\hat {M}}$ con

$MA=\left({\begin{array}{cc}\alpha &{\boldsymbol {a}}^{t}\\0&A_{\ast }-1/\alpha {\boldsymbol {a}}{\boldsymbol {a}}^{t}\end{array}}\right)$

${\hat {M}}=\left({\begin{array}{cc}1&-{\boldsymbol {a}}^{t}/\alpha \\0&I\end{array}}\right)$

(ora la sottocolonna dei coefficienti è sulla prima riga)

Calcolando il prodotto

$MA{\hat {M}}=\left({\begin{array}{cc}\alpha &0\\0&A_{\ast }-1/\alpha {\boldsymbol {a}}{\boldsymbol {a}}^{t}\end{array}}\right)$

Osserviamo che ${\hat {M}}=M^{t}$. Ponendo ${\boldsymbol {y}}=M^{t}{\boldsymbol {x}}$, per ogni ${\boldsymbol {x}}\neq 0$,

${\boldsymbol {x}}^{t}MA{\hat {M}}{\boldsymbol {x}}={\boldsymbol {x}}^{t}MAM^{t}{\boldsymbol {x}}={\boldsymbol {y}}^{t}A{\boldsymbol {y}}$

e ${\boldsymbol {y}}=M^{t}{\boldsymbol {x}}\neq 0\forall {\boldsymbol {x}}\neq 0$ perché $M$ è non singolare. Siccome per ipotesi $A$ è definita positiva, ${\boldsymbol {y}}^{t}A{\boldsymbol {y}}$ è maggiore di 0, e quindi la matrice grande $MA{\hat {M}}$ è definita positiva.

Allora

$A_{\ast }-1/\alpha {\boldsymbol {a}}{\boldsymbol {a}}^{t}$

è definita positiva, perché sottomatrice principale di una definita positiva. Quindi ${\hat {A}}_{\ast }=R_{\ast }^{t}R$ è definita positiva e la procedura si può iterare sulla matrice di dimensioni $n-1\times n-1$. Lavorando su $R_{\ast }^{t}R_{\ast }$, si determinano progressivamente gli elementi di $R_{\ast }$, triangolare superiore.

cvd

Osservazione 2.4

Consideriamo

$A_{\ast }-1/\alpha {\boldsymbol {a}}{\boldsymbol {a}}^{t}.$

Alla matrice $A_{\ast }$ definita positiva, sottraiamo una matrice semidefinita positiva, e non è ovvio che il risultato sia una definita positiva. Infatti, ad esempio

$I-(0,1)(0,1)^{t}$

non è definita positiva.

Complemento di Shur

Questa dimostrazione permette di osservare la seguente proprietà: il complemento di Shur di una definita positiva è definito positivo, dove ${\hat {A}}_{\ast }$ è il complemento di Shur di $\alpha$ in $A$.

Definizione

In generale, data una matrice ${\bar {A}}$ a blocchi della forma

$\left({\begin{array}{cc}A_{11}&A_{12}\\A_{21}&A_{22}\end{array}}\right)$

dove il primo blocco è $n\times n$, il secondo $n\times m$, e così via, e dove la dimensione della matrice è $(n+m)\times (n+m)$, si chiede che $A_{11}$ sia non singolare, allora il complemento di Shur è definito come

$A_{22}-A_{21}A_{11}^{-1}A_{12}$

Algoritmo per il calcolo di R

In questa scrittura supponiamo che gli elementi di $A$ vengano progressivamente sovrascritti con quelli di $R$ trovati passo per passo.

${\rm {{for}k=1:n}}$

$A(k,k)={\sqrt {A(k,k)}}$

${\rm {end}}$

$fork=1:n-1$

$forj=k+1:n$

$A(k,j)=A(k,j)/A(k,k)$

$A(j,k)=0$

${\rm {end}}$

${\rm {{for}i=k+1:n}}$

${\rm {{for}j=k+1:n}}$

$A(i,j)=A(i,j)-A(i,k)*A(j,k)$

${\rm {end}}$

${\rm {end}}$

${\rm {end}}$

Descrizione esplicita: per $k=1:n$ calcolo

$R_{kk}=A_{kk}^{(k+1)}={\sqrt {A_{kk}}}=\rho$

e al passo $n$-esimo questa è l'unica operazione da fare.

Allora, per $k=1:n-1$, facciamo le seguenti operazioni:

1. per $j=k+1:n$,

$R_{kj}=A_{kj}^{(k+1)}={\frac {A_{kj}}{A_{kk}}}={\frac {A_{kj}}{R_{kk}}}$

(queste istruzioni corrispondono a ${\frac {\boldsymbol {a}}{\sqrt {\alpha }}}$, definizione del vettore ${\boldsymbol {r}}$)

2. per $i=k+1:n$, e per $j=k+1:n$,

${\hat {A}}_{ij,\ast }=A_{ij}^{(k+1)}=A_{ij}-A_{ki}A_{kj}=A_{ij}-R_{ki}R_{kj}$

(nell'ultimo blocco di istruzioni si definiscono i coefficienti di ${\hat {A}}_{\ast }$, cioè $A_{\ast }-{\boldsymbol {r}}{\boldsymbol {r}}^{t}$).

Ovviamente, alla fine di ogni passo $k$, le entrate $A_{kj}$ con $j<k$ vengono annullate.

Costo computazionale

È sufficiente calcolare solo metà della sottomatrice $A(k+1:n,k+1:n)$, che è simmetrica.

Le operazioni richieste sono quindi:

1. $n$ estrazioni di radici quadrate (una per ogni passo di fattorizzazione);

2. $n-k$ operazioni moltiplicative per ogni passo $k$. Infatti, al passo ${\boldsymbol {r}}^{t}$ ha indici da $k+1$ a $n$, eper ciascun elemento houna divisione;

3. per ogni elemento di ${\hat {A}}_{\ast }$, ho un'operazione moltiplicativa, e quindi per ogni passo le operazioni sono

$\sum _{i=1}^{m}i=m(m+1)/2$

con $m=n-k+1-1=n-k$, cioè $m(m+1)/2=(n-k)(n-k+1)/2$

Complessivamente il numero di operazioni è:

$n+\sum _{k=1}^{n-1}[n-k+(n-k)(n-k+1)/2]$

e raccogliendo $n-k$ otteniamo:

$n+\sum _{k=1}^{n-1}(n-k)(1+(n-k+1)/2)$

$n+1/2\sum _{k=1}^{n-1}(n-k)(2+n-k+1)$

$n+1/2\sum _{k=1}^{n-1}(n-k)(n-k+3)$

$\simeq n+1/2\sum _{k=1}^{n-1}(n-k)^{2}$

$=n^{3}/6$

Il costo computazionale è pari a quello dell'eliminazione di Gauss per le simmetriche.

Stabilità

Considerando il prodotto $A=R^{t}R$, segue che

$A_{kk}=\sum _{j=1}^{n}R_{kj}^{t}R_{jk}$

$A_{kk}=\sum _{j=1}^{n}R_{jk}^{2}$

Allora ciascun $R_{jk}$ è tale che $R_{jk}^{2}\leq A_{kk}$, quindi

$|R_{jk}|\leq {\sqrt {A_{kk}}}\,\forall k$

Chiamando $R_{m}$ l'elemento di modulo massimo, si ha

$R_{m}\leq {\sqrt {A_{m}}}$

Quindi gli elementi di $R$ hanno crescita controllata da ${\sqrt {A_{m}}}$.

L'algoritmo è stabile.

Rapporto tra fattorizzazione LU e di Cholesky

Data una matrice simmetrica definita positiva, consideriamo il rapporto tra la fattorizzazione $LU$ e quella di Cholesky. La fattorizzazione $LU$ esiste unica per una definita positiva, perché le matrici di testa di una definita positiva sono definite positive, e non singolari, quindi le ipotesi del teorema di esistenza e unicità della fattorizzazione sono soddisfatte.

Tuttavia $L\neq U^{t}$ perché la diagonale di $L$ ha entrate uguali a 1.

Poniamo $U=D*{\bar {U}}$ dove $D$ è una matrice con elementi diagonali dati dalla diagonale principale di $U$.

Per ipotesi

$A=LU=A^{t}=(LD{\bar {U}})^{t}={\bar {U}}^{t}DL^{t}$

e per l'unicità di fattorizzazione,

$LU={\bar {U}}^{t}DL^{t}\quad \longrightarrow \quad L={\bar {U}}^{t},\quad DL^{t}=U$

allora

$A={\bar {U}}^{t}D{\bar {U}}$

con $D$ matrice a elementi diagonali positivi. Poniamo ${\boldsymbol {x}}={\bar {U}}{\boldsymbol {y}}$, allora

${\boldsymbol {x}}^{t}D{\boldsymbol {x}}=({\bar {U}}{\boldsymbol {y}})^{t}D({\bar {U}}{\boldsymbol {y}})$

$={\boldsymbol {y}}^{t}{\bar {U}}^{t}D{\bar {U}}{\boldsymbol {y}}={\boldsymbol {x}}^{t}A{\boldsymbol {x}}>0$

perché $A$ è definita positiva, allora ${\bar {U}}$ è non singolare.

$A={\bar {U}}^{t}D{\bar {U}}$

Definiamo $D^{1/2}$ come la matrice che ha come elementi diagonali le radici quadrate degli elementi di $D$.

Allora

$A={\bar {U}}^{t}D^{1/2}D^{1/2}{\bar {U}}=(D^{1/2}{\bar {U}})^{t}D^{1/2}{\bar {U}}=R^{t}R$

cioè

$R=D^{1/2}{\bar {U}}$

dove ${\bar {U}}=D^{-1}U$, $R=D^{1/2}D^{-1}U=D^{-1/2}U$.

Proposizione 2.1

L'eliminazione gaussiana conserva simmetria e definita positività.

Dim.

Simmetria: sappiamo che

$A_{ij}^{(k+1)}=A_{ij}^{(k)}-M_{ik}A_{kj}^{(k)}$

con

$M_{ik}={\frac {A_{ik}^{(k)}}{A_{kk}^{(k)}}}$

cioè, sostituendo esplicitamente $M_{ik}$:

$A_{ij}^{(k+1)}=A_{ij}^{(k)}-{\frac {A_{ik}^{(k)}}{A_{kk}^{(k)}}}A_{kj}^{(k)}$

$A_{ji}^{(k)}=A_{ij}^{(k)}$ per simmetria di $A$, quindi si può scrivere:

$A_{ij}^{(k+1)}=A_{ji}^{(k)}-{\frac {A_{ki}^{(k)}}{A_{kk}^{(k)}}}A_{jk}^{(k)}$

e tenendo conto che ${\frac {A_{jk}}{A_{kk}}}=M_{jk}$ si ha

$A_{ij}^{(k+1)}=A_{ji}^{(k)}-M_{jk}A_{ki}^{(k)}=A_{ji}^{(k+1)}$

cioè la sottomatrice in basso a destra ottenuta dopo il passo $k$-esimo è ancora simmetrica.

Definita positività: consideriamo

$A=\left({\begin{array}{cc}\alpha &{\boldsymbol {a}}^{t}\\{\boldsymbol {a}}&A_{\ast }\end{array}}\right)$

$M=\left({\begin{array}{cc}1&0\\-{\boldsymbol {a}}/\alpha &I\end{array}}\right)$

$MA=\left({\begin{array}{cc}\alpha &{\boldsymbol {a}}^{t}\\0&A_{\ast }-1/\alpha {\boldsymbol {a}}{\boldsymbol {a}}^{t}\end{array}}\right)$

L'ultima entrata è il complemento di Shur di $\alpha$ in $A$. Siccome $A$ è definita positiva, anche il completamento è definito positivo. Allora l'eliminazione gaussiana è stabile sulle definite positive.

cvd