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 {ℝ}^{n\times n}$ simmetrica e non singolare. $A$ è definita positiva se e solo se esiste $R\in {ℝ}^{n\times n}$ triangolare superiore tale che $A=R^{t}R$.

Dimostrazione

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

$detA=det{R}^t\ast detR=|detR{|}^2$

e quindi $detA\ne 0$ implica $R$ non singolare. Allora

$\mathrm{\forall }𝒙\ne 0,𝒚=R𝒙\ne 0,$

quindi

$|𝒚{|}_2^2\ge 0$

e siccome ho escluso il caso nullo,

$|𝒚{|}_2^2>0.$

e $A$ è definita positiva.

$1⟶2$: consideriamo $A$ simmetrica e definita positiva:

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

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

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

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

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

Il prodotto delle due matrici è:

$R^{t}R=\left(\begin{array}{cc}{\rho }^2& \rho \ast {𝒓}^t\\ \rho \ast 𝒓& 𝒓{𝒓}^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{array}{cc}\alpha ={\rho }^2& ⟶\rho =\sqrt{\alpha }\\ 𝒂=\rho 𝒓& 𝒓=\frac{1}{\sqrt{\alpha }}𝒂\\ {𝒂}^t=\rho {𝒓}^t& ⟶{𝒓}^t=\frac{{𝒂}^t}{\sqrt{\alpha }}\\ A_{\ast }=𝒓{𝒓}^t+R_{\ast }^t{R}_{\ast }& ⟶R_{\ast }^t{R}_{\ast }=A_{\ast }-𝒓{𝒓}^t=A_{\ast }-1/\alpha 𝒂{𝒂}^t\end{array}$

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 & {𝒂}^t\\ 𝒂& A_{\ast }\end{array}\right)$

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

Il risultato del prodotto è:

$MA=\left(\begin{array}{cc}\alpha & {𝒂}^t\\ 0& -1/\alpha 𝒂{𝒂}^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 & {𝒂}^t\\ 0& A_{\ast }-1/\alpha 𝒂{𝒂}^t\end{array}\right)$

$\hat{M}=\left(\begin{array}{cc}1& -{𝒂}^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 𝒂{𝒂}^t\end{array}\right)$

Osserviamo che $\hat{M}=M^t$. Ponendo $𝒚=M^t𝒙$, per ogni $𝒙\ne 0$,

${𝒙}^{t}MA\hat{M}𝒙={𝒙}^{t}MA{M}^t𝒙={𝒚}^{t}A𝒚$

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

Allora

$A_{\ast }-1/\alpha 𝒂{𝒂}^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 𝒂{𝒂}^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$.

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.

$\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{k}\mathrm{=}\mathrm{1}\mathrm{:}\mathrm{n}$

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

$\mathrm{e}\mathrm{n}\mathrm{d}$

$fork=1:n-1$

$forj=k+1:n$

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

$A(j,k)=0$

$\mathrm{e}\mathrm{n}\mathrm{d}$

$\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{=}\mathrm{k}\mathrm{+}\mathrm{1}\mathrm{:}\mathrm{n}$

$\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{j}\mathrm{=}\mathrm{k}\mathrm{+}\mathrm{1}\mathrm{:}\mathrm{n}$

$A(i,j)=A(i,j)-A(i,k)\ast A(j,k)$

$\mathrm{e}\mathrm{n}\mathrm{d}$

$\mathrm{e}\mathrm{n}\mathrm{d}$

$\mathrm{e}\mathrm{n}\mathrm{d}$

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{𝒂}{\sqrt{\alpha }}$, definizione del vettore $𝒓$)

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 }-𝒓{𝒓}^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 ${𝒓}^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\le A_{kk}$, quindi

$|R_{jk}|\le \sqrt{A_{kk}}\mathrm{\forall }k$

Chiamando $R_m$ l'elemento di modulo massimo, si ha

$R_m\le \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\ne U^t$ perché la diagonale di $L$ ha entrate uguali a 1.

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

Per ipotesi

$A=LU=A^t=(LD\overline{U}{)}^t={\overline{U}}^{t}D{L}^t$

e per l'unicità di fattorizzazione,

$LU={\overline{U}}^{t}D{L}^t⟶L={\overline{U}}^t,D{L}^t=U$

allora

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

con $D$ matrice a elementi diagonali positivi. Poniamo $𝒙=\overline{U}𝒚$, allora

${𝒙}^{t}D𝒙=(\overline{U}𝒚{)}^{t}D(\overline{U}𝒚)$

$={𝒚}^t{\overline{U}}^{t}D\overline{U}𝒚={𝒙}^{t}A𝒙>0$

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

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

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

Allora

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

cioè

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

dove $\overline{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 & {𝒂}^t\\ 𝒂& A_{\ast }\end{array}\right)$

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

$MA=\left(\begin{array}{cc}\alpha & {𝒂}^t\\ 0& A_{\ast }-1/\alpha 𝒂{𝒂}^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