Queste note sono basate sugli appunti fatti con Gianmarco Manzini negli anni 1995-2005

Fattorizzazione di Cholesky

Appunti di calcolo numerico

Autore/Autrice

Affiliazione

Enrico Bertolazzi

University of Trento, Department of Industrial Engineering

Sia \boldsymbol{A}\in \Bbb{R}^{n \times n} una matrice non singolare e \boldsymbol{b}\in \Bbb{R}^n. Vogliamo determinare \boldsymbol{x}\in \Bbb{R}^n tale che

\boldsymbol{A}\boldsymbol{x}= \boldsymbol{b}.

Sappiamo che è possibile risolvere il problema tramite la fattorizzazione di Gauss \boldsymbol{A}= \boldsymbol{L}\boldsymbol{U}, dove, per semplicità di esposizione, assumiamo che non sia necessario il pivoting. Questo approccio richiede di risolvere due sistemi lineari con matrici triangolari, uno inferiore (\boldsymbol{L}) e uno superiore (\boldsymbol{U}):

\boldsymbol{L}\boldsymbol{y}= \boldsymbol{b}, \qquad \boldsymbol{U}\boldsymbol{x}= \boldsymbol{y}.

Il costo computazionale di questa procedura non risiede soltanto nella soluzione dei due sistemi, ma anche nel processo di calcolo dei fattori \boldsymbol{L} e \boldsymbol{U}. Inoltre, il calcolatore deve disporre di sufficiente memoria per memorizzare entrambi i fattori. Questo non rappresenta un problema per matrici piene, poiché la memoria già occupata dalla matrice originale può essere riutilizzata. Tuttavia, per matrici sparse, il processo di fill-in (creazione di nuovi elementi non nulli durante la fattorizzazione) può incrementare notevolmente i requisiti di memoria, rendendo difficile la memorizzazione dei due fattori.

Per queste ragioni, è utile considerare situazioni in cui è possibile una fattorizzazione alternativa del tipo:

\boldsymbol{A}= \boldsymbol{R}^T \boldsymbol{R},

dove \boldsymbol{R} è una matrice triangolare superiore, detta fattorizzazione di Cholesky.

Poiché questa fattorizzazione richiede il calcolo di un solo fattore anziché due, possiamo aspettarci una riduzione sia del costo computazionale sia dell’occupazione di memoria. Sebbene il guadagno non sia necessariamente dimezzato, si ottiene comunque un notevole risparmio¹.

Il seguente teorema descrive le matrici per le quali la fattorizzazione di Cholesky è possibile.

Teorema 1 (Fattorizzazione di Cholesky e matrici SPD) Le matrici quadrate non-singolari \boldsymbol{A} per cui è possibile una fattorizzazione della forma \boldsymbol{A}=\boldsymbol{R}^T\boldsymbol{R} con \boldsymbol{R} matrice triangolare superiore sono tutte e sole le matrici simmetriche e definite positive, (SPD).

Dimostrazione. Mostriamo innanzitutto la necessità, cioè che la fattorizzazione di Cholesky è possibile solo per matrici simmetriche e definite positive. Mostreremo in seguito la sufficienza, cioè che data una matrice che ha tali proprietà è sempre possibile calcolare il suo fattore di Cholesky \boldsymbol{R}.

• Necessità
  Sia \boldsymbol{A},\,\boldsymbol{R}\in\Bbb{R}^{n\times n} con \boldsymbol{A} non singolare ed \boldsymbol{R} triangolare superiore tale che valga la decomposizione \boldsymbol{A}=\boldsymbol{R}^T\boldsymbol{R}.

  • Simmetria
    \boldsymbol{A}^T = (\boldsymbol{R}^T\boldsymbol{R})^T = \boldsymbol{R}^T\boldsymbol{R}= \boldsymbol{A}.

  • Positività
    Se \boldsymbol{A} è non-singolare, allora anche \boldsymbol{R} è non-singolare perché dal teorema di Binet segue
    \left|\begin{matrix}\boldsymbol{A}\end{matrix}\right|=\left|\begin{matrix}\boldsymbol{R}^T\boldsymbol{R}\end{matrix}\right|=\left|\begin{matrix}\boldsymbol{R}^T\end{matrix}\right|\left|\begin{matrix}\boldsymbol{R}\end{matrix}\right|=\left|\begin{matrix}\boldsymbol{R}\end{matrix}\right|^2.

    Inoltre, se \boldsymbol{R} è non-singolare, si ha \boldsymbol{R}\boldsymbol{x}\neq\boldsymbol{0} per ogni vettore \boldsymbol{x}\neq\boldsymbol{0}. Quindi,
    \boldsymbol{x}^T\boldsymbol{A}\boldsymbol{x}= \boldsymbol{x}^T\boldsymbol{R}^T\boldsymbol{R}\boldsymbol{x}= \left|\left|\boldsymbol{R}\boldsymbol{x}\right|\right|^2_{2} > 0

• Sufficienza
  La dimostrazione procede per induzione sull’ordine della matrice, ma è anche in parte costruttiva. Illustreremo un algorimo di riduzione, che permette di ricondurre il problema della fattorizzazione di una matrice SPD di ordine n a quello di una matrice SPD di ordine n-1. Per induzione, si troverà che tutte le matrici SPD sono fattorizzabili.

  Per facilitare la costruzione dell’algoritmo di riduzione, riscriviamo la matrice simmetrica \boldsymbol{A}\in\Bbb{R}^{n\times n} in una forma equivalente a blocchi

  \boldsymbol{A}= \left(\begin{array}{c|ccc} A_{1,1} & A_{1,2} & \cdots & A_{1,n} \\ \hline A_{2,1} & A_{2,2} & \cdots & A_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ A_{n,1} & A_{n,2} & \cdots & A_{n,n} \end{array}\right) \begin{pmatrix} \alpha & \boldsymbol{a}^T \\ \boldsymbol{a}& \tilde{\boldsymbol{A}} \end{pmatrix}

  introducendo la sottomatrice

  \tilde{\boldsymbol{A}}\in\Bbb{R}^{(n-1)\times (n-1)}, il vettore \boldsymbol{a}\in\Bbb{R}^{n-1} e lo scalare \alpha\in\Bbb{R}^n \tilde{\boldsymbol{A}} = \begin{pmatrix} A_{2,2} & \cdots & A_{2,n} \\ \vdots & \ddots & \vdots \\ A_{n,2} & \cdots & A_{n,n} \end{pmatrix},\qquad \boldsymbol{a}= \begin{pmatrix} A_{2,1} \\ \vdots \\ A_{n,1} \end{pmatrix}, \qquad \alpha = A_{1,1}. Scriviamo la matrice \boldsymbol{R} a blocchi in modo simile ad \boldsymbol{A} \boldsymbol{R}= \left(\begin{array}{c|c} \rho & \boldsymbol{r}^T \\ \boldsymbol{0}& \boldsymbol{R}^\prime \end{array}\right) e sviluppiamo il prodotto \boldsymbol{R}^T\boldsymbol{R} \boldsymbol{R}^T\boldsymbol{R}= \begin{pmatrix} \rho & \boldsymbol{0}^{T} \\ \boldsymbol{r}& (\boldsymbol{R}^\prime)^T \end{pmatrix} \begin{pmatrix} \rho & \boldsymbol{r}^T \\ \boldsymbol{0}& \boldsymbol{R}^\prime \end{pmatrix} = \begin{pmatrix} \rho^2 & \rho\boldsymbol{r}^T \\ \rho\boldsymbol{r}& \boldsymbol{r}\boldsymbol{r}^T + (\boldsymbol{R}^\prime)^T\boldsymbol{R}^\prime \end{pmatrix}. La fattorizzazione \boldsymbol{A}=\boldsymbol{R}^T\boldsymbol{R} si può esprimere uguagliando i blocchi nelle posizioni corrispondenti di \boldsymbol{R}^{T}\boldsymbol{R} ed \boldsymbol{A} \begin{pmatrix} \alpha & \boldsymbol{a}^T \\ \boldsymbol{a}& \tilde{\boldsymbol{A}} \end{pmatrix} = \begin{pmatrix} \rho^2 & \rho\boldsymbol{r}^T \\ \rho\boldsymbol{r}& \boldsymbol{r}\boldsymbol{r}^T + (\boldsymbol{R}^\prime)^T\boldsymbol{R}^\prime \end{pmatrix}, per cui si ottengono immediatamente le relazioni seguenti \alpha = \rho^2,\qquad \boldsymbol{a}= \rho\boldsymbol{r}, o equivalentemente \boldsymbol{a}^T= \rho\boldsymbol{r}^T,\qquad \tilde{\boldsymbol{A}} = (\boldsymbol{R}^\prime)^T\boldsymbol{R}^\prime + \boldsymbol{r}\boldsymbol{r}^T. \tag{1} Le prime due relazioni forniscono una espressione diretta per determinare \rho ed \boldsymbol{r}, cioè \rho = \sqrt{\alpha},\qquad\textrm{e}\qquad \boldsymbol{r}= \dfrac{1}{\sqrt{\alpha}}\boldsymbol{a}, dove supponiamo per il momento² che \alpha sia un numero reale strettamente positivo.

  L’ultima relazione in Equazione 1, riscritta come (\boldsymbol{R}^\prime)^T\boldsymbol{R}^\prime = \tilde{\boldsymbol{A}} - \dfrac{1}{\alpha}\boldsymbol{a}\boldsymbol{a}^T, suggerisce che la matrice \tilde{\boldsymbol{R}} sia a sua volta il fattore triangolare superiore – ammesso ovviamente che esista – della fattorizzazione di Cholesky della matrice ridotta \boldsymbol{A}^\prime\in\Bbb{R}^{(n-1)\times (n-1)} \boldsymbol{A}^\prime= \tilde{\boldsymbol{A}} - \dfrac{1}{\alpha}\boldsymbol{a}\boldsymbol{a}^T. Procedendo in questo modo potremmo scrivere il fattore triangolare superiore \boldsymbol{R} di Cholesky come \boldsymbol{R}= \begin{pmatrix} \sqrt{\alpha} & \frac{1}{\sqrt{\alpha}}\boldsymbol{a}^T \\ \boldsymbol{0}& \boldsymbol{R}^\prime \end{pmatrix}, in cui la prima riga sarebbe completamente determinata e resterebbe da calcolare il fattore triangolare superiore \boldsymbol{R}^\prime della matrice \boldsymbol{A}^\prime di ordine n-1. Di fatto, avremmo ridotto di uno la dimensione del problema ed è evidente che se si potesse ripetere questo procedimento n volte, arriveremmo alla fine ad una matrice di ordine 1 la cui fattorizzazione sarebbe banalmente \begin{pmatrix}A_{1,1}\end{pmatrix} = \begin{pmatrix}\sqrt{A_{1,1}}\end{pmatrix}\begin{pmatrix}\sqrt{A_{1,1}}\end{pmatrix}.

  Mostreremo che l’ipotesi che \boldsymbol{A} sia SPD garantisce che tutto il procedimento ipotizzato è in realtà ben definito. Innanzitutto, osserviamo che valgono le seguenti proprietà.

  • (a) Se \boldsymbol{A} è SPD, allora \alpha>0
    Il risultato si ottiene applicando la definizione di matrice definita positiva ad \boldsymbol{e}_1, primo vettore della base canonica:
    \alpha = A_{1,1} = \boldsymbol{e}_1^T\boldsymbol{A}\boldsymbol{e}_1 > 0.

  • (b) Se \boldsymbol{A} è SPD, allora \boldsymbol{A}^\prime è SPD
    La simmetria di \boldsymbol{A}^\prime è evidente dalla sua definizione. Vogliamo mostrare inoltre che per ogni vettore \boldsymbol{y}\in\Bbb{R}^{n-1} non nullo si ha \boldsymbol{y}^T\boldsymbol{A}^\prime\boldsymbol{y}> 0 \qquad\textrm{cioè}\qquad \boldsymbol{y}^T\left(\tilde{\boldsymbol{A}} - \alpha^{-1}\boldsymbol{a}^T\boldsymbol{a}\right)\boldsymbol{y}> 0, nell’ipotesi che \boldsymbol{A} sia definita positiva. Dato un vettore generico \boldsymbol{y}\in\Bbb{R}^{n-1}, introduciamo un vettore ausiliario \boldsymbol{z}= \begin{pmatrix}\eta & \boldsymbol{y}^T\end{pmatrix}^{T}\in\Bbb{R}^{n}, dove \eta è a sua volta un qualsiasi numero reale non nullo. La positività di \boldsymbol{A} implica che \boldsymbol{z}^T\boldsymbol{A}\boldsymbol{z}>0, cioè \begin{pmatrix} \eta & \boldsymbol{y}^T\end{pmatrix}\begin{pmatrix} \alpha & \boldsymbol{a}^T \\ \boldsymbol{a}& \boldsymbol{A}^\prime \end{pmatrix} \begin{pmatrix} \eta \\ \boldsymbol{y}\end{pmatrix} = \alpha\,\eta^2 + 2\boldsymbol{a}^T\boldsymbol{y}\,\eta + \boldsymbol{y}^T\tilde{\boldsymbol{A}}\boldsymbol{y}> 0. La positività di questo polinomio di secondo grado in \eta implica che il discriminante sia negativo (ovvero non ci sono radici reali), per cui sicuramente vale \Delta = \Big(\boldsymbol{a}^T\boldsymbol{y}\Big)^2 - \alpha\boldsymbol{y}^T\tilde{\boldsymbol{A}}\boldsymbol{y}<0, e questa è proprio la condizione che dobbiamo verificare!

Sfruttando le due proprietà (a) e (b) possiamo concludere la dimostrazione procedendo per induzione sull’ordine della matrice. Assumiamo che \boldsymbol{A}\in\Bbb{R}^{n\times n} sia SPD, e che si abbia

• per n=1 una matrice SPD si fattorizza come già riportato prima, e cioè \begin{pmatrix}A_{1,1}\end{pmatrix} = \begin{pmatrix}\sqrt{A_{1,1}}\end{pmatrix}\begin{pmatrix}\sqrt{A_{1,1}}\end{pmatrix}.

• per ogni matrice \boldsymbol{A}^\prime SPD di ordine n-1 esiste una matrice triangolare superiore \boldsymbol{R} che la fattorizza nella forma di Cholesky \boldsymbol{A}^\prime = \boldsymbol{R}^\prime(\boldsymbol{R}^\prime)^T.

Allora, le proprietà (a) e (b) implicano che se \boldsymbol{A} è SPD è possibile operare la riduzione, e dalle ipotesi induttive segue subito che la matrice è fattorizabile in forma di Cholesky.

Algoritmo di calcolo

In questo paragrafo, descriviamo il processo logico della fattorizzazione mediante il metodo di Cholesky. Analizziamo in dettaglio il primo e il secondo step, lo step generico k e

Step 1

Il teorema appena dimostrato ci permette di scrivere la matrice \boldsymbol{A} con una decomposizione del tipo:

\begin{aligned} \boldsymbol{A}&= \begin{pmatrix} \alpha & \boldsymbol{a}^T \\ \boldsymbol{a}& \tilde{\boldsymbol{A}} \end{pmatrix} \\ &= \begin{pmatrix} \sqrt{\alpha} & 0 \\ \frac{1}{\sqrt{\alpha}} \boldsymbol{a}& \boldsymbol{I}_{n-1} \end{pmatrix} \begin{pmatrix} 1 & 0^T \\ 0 & \boldsymbol{A}^\prime \end{pmatrix} \begin{pmatrix} \sqrt{\alpha} & \frac{1}{\sqrt{\alpha}} \boldsymbol{a}^T \\ 0 & \boldsymbol{I}_{n-1} \end{pmatrix} \\ & = \boldsymbol{R}_1^T \boldsymbol{A}_1 \boldsymbol{R}_1 \end{aligned}

dove le matrici \boldsymbol{A}_1 e \boldsymbol{R}_1 sono definite in modo naturale. Inoltre, abbiamo introdotto la matrice ridotta di ordine (n-1):

\boldsymbol{A}^\prime = \tilde{\boldsymbol{A}} - \dfrac{1}{\alpha} \boldsymbol{a}\boldsymbol{a}^T,

che risulta essere SPD (definita positiva), poiché la matrice originale \boldsymbol{A} è anch’essa SPD. La matrice \boldsymbol{A}^\prime è chiamata complemento di Schur della matrice \tilde{\boldsymbol{A}} rispetto all’elemento pivotale \alpha.

Step 2

Poiché la matrice ridotta \boldsymbol{A}^\prime è ancora SPD (definita positiva), possiamo ripetere su di essa lo stesso processo del primo step. Iniziamo riscrivendo \boldsymbol{A}^\prime in forma a blocchi:

\boldsymbol{A}^\prime = \begin{pmatrix} \alpha^\prime & (\boldsymbol{a}^\prime)^T \\ \boldsymbol{a}^\prime & \tilde{\boldsymbol{A}}^\prime \end{pmatrix}

Applichiamo quindi la riduzione, decomponendo \boldsymbol{A}^\prime come segue:

\begin{pmatrix} \alpha^\prime & (\boldsymbol{a}^\prime)^T \\ \boldsymbol{a}^\prime & \tilde{\boldsymbol{A}}^\prime \end{pmatrix} = \begin{pmatrix} \sqrt{\alpha^\prime} & 0 \\ \frac{1}{\sqrt{\alpha^\prime}} \boldsymbol{a}^\prime & \boldsymbol{I}_{n-2} \end{pmatrix} \begin{pmatrix} 1 & 0^T \\ 0 & \boldsymbol{A}^{\prime\prime} \end{pmatrix} \begin{pmatrix} \sqrt{\alpha^\prime} & \frac{1}{\sqrt{\alpha^\prime}} (\boldsymbol{a}^\prime)^T \\ 0 & \boldsymbol{I}_{n-2} \end{pmatrix}.

La matrice ridotta \boldsymbol{A}^{\prime\prime} \in \Bbb{R}^{(n-2)\times(n-2)} è data da:

\boldsymbol{A}^{\prime\prime} = \tilde{\boldsymbol{A}}^\prime - \dfrac{1}{\alpha^\prime} \boldsymbol{a}^\prime (\boldsymbol{a}^\prime)^T,

che rappresenta il complemento di Schur di \tilde{\boldsymbol{A}}^\prime rispetto al pivot \alpha^\prime. Dato che \boldsymbol{A}^\prime è SPD, anche \boldsymbol{A}^{\prime\prime} sarà SPD.

A questo punto abbiamo decomposto la matrice \boldsymbol{A}_1, ottenuta nel primo step, come segue:

\begin{aligned} \boldsymbol{A}_1 &= \begin{pmatrix} 1 & 0^T \\ 0 & \begin{pmatrix} \alpha^\prime & (\boldsymbol{a}^\prime)^T \\ \boldsymbol{a}^\prime & \tilde{\boldsymbol{A}}^\prime \end{pmatrix} \end{pmatrix} \\ &= \begin{pmatrix} 1 & 0^T \\ 0 & \begin{pmatrix} \sqrt{\alpha^\prime} & 0 \\ \frac{1}{\sqrt{\alpha^\prime}} \boldsymbol{a}^\prime & \boldsymbol{I}_{n-2} \end{pmatrix} \end{pmatrix} \begin{pmatrix} 1 & 0^T \\ 0 & \begin{pmatrix} 1 & 0^T \\ 0 & \boldsymbol{A}^{\prime\prime} \end{pmatrix} \end{pmatrix} \begin{pmatrix} 1 & 0^T \\ 0 & \begin{pmatrix} \sqrt{\alpha^\prime} & \frac{1}{\sqrt{\alpha^\prime}} (\boldsymbol{a}^\prime)^T \\ 0 & \boldsymbol{I}_{n-2} \end{pmatrix} \end{pmatrix} \\ &= \boldsymbol{R}_2^T \boldsymbol{A}_2 \boldsymbol{R}_2, \end{aligned}

dove \boldsymbol{R}_2 e \boldsymbol{A}_2 sono definite in modo naturale. Dopo questi due step, abbiamo ottenuto la seguente decomposizione per la matrice \boldsymbol{A}:

\boldsymbol{A}= \boldsymbol{R}_1^T \boldsymbol{A}_1 \boldsymbol{R}_1 = \boldsymbol{R}_1^T \boldsymbol{R}_2^T \boldsymbol{A}_2 \boldsymbol{R}_2 \boldsymbol{R}_1.

Step k generico

Dopo k-1 step, abbiamo decomposto la matrice \boldsymbol{A} nel prodotto dei fattori:

\boldsymbol{A}= \boldsymbol{R}_1^T \boldsymbol{R}_2^T \ldots \boldsymbol{R}_{k-1}^T \boldsymbol{A}_{k-1} \boldsymbol{R}_{k-1} \ldots \boldsymbol{R}_2 \boldsymbol{R}_1,

dove

\boldsymbol{A}_{k-1} = \begin{pmatrix} \boldsymbol{I}_{k-1} & \boldsymbol{0}^T \\ \boldsymbol{0}& \boldsymbol{A}^{(k-1)} \end{pmatrix}.

Dal teorema generale, si deduce che il blocco di matrice \boldsymbol{A}^{(k-1)} \in \Bbb{R}^{(n-k+1) \times (n-k+1)} è dato da:

\boldsymbol{A}^{(k-1)} = \begin{pmatrix} \alpha^{(k-1)} & (\boldsymbol{a}^{(k-1)})^T \\ \boldsymbol{a}^{(k-1)} & \tilde{\boldsymbol{A}}^{(k-1)} \end{pmatrix},

ed è una matrice SPD, il che consente di continuare con la decomposizione.

Procediamo ora a decomporre a blocchi la matrice \boldsymbol{A}_{k-1}:

\begin{aligned} \boldsymbol{A}_{k-1} &= \begin{pmatrix} \boldsymbol{I}_{k-1} & \boldsymbol{0}^T \\ \boldsymbol{0}& \begin{pmatrix} \alpha^{(k-1)} & (\boldsymbol{a}^{(k-1)})^T \\ \boldsymbol{a}^{(k-1)} & \tilde{\boldsymbol{A}}^{(k-1)} \end{pmatrix} \end{pmatrix} \\ &= \begin{pmatrix} \boldsymbol{I}_{k-1} & \boldsymbol{0}^T \\ \boldsymbol{0}& \begin{pmatrix} \sqrt{\alpha^{(k-1)}} & \boldsymbol{0}\\ \frac{1}{\sqrt{\alpha^{(k-1)}}}\boldsymbol{a}^{(k-1)} & \boldsymbol{I}_{n-k} \end{pmatrix} \end{pmatrix} \\ &\qquad \begin{pmatrix} \boldsymbol{I}_{k-1} & \boldsymbol{0}^T \\ \boldsymbol{0}& \begin{pmatrix} 1 & \boldsymbol{0}^T \\ \boldsymbol{0}& \boldsymbol{A}^{(k)} \end{pmatrix} \end{pmatrix} \\ &\qquad\qquad \begin{pmatrix} \boldsymbol{I}_{k-1} & \boldsymbol{0}^T \\ \boldsymbol{0}& \begin{pmatrix} \sqrt{\alpha^{(k-1)}} & \frac{1}{\sqrt{\alpha^{(k-1)}}}(\boldsymbol{a}^{(k-1)})^T \\ \boldsymbol{0}& \boldsymbol{I}_{n-k} \end{pmatrix} \end{pmatrix} \end{aligned}

Dopo k step, abbiamo quindi decomposto la matrice \boldsymbol{A} nella forma seguente:

\boldsymbol{A}= \boldsymbol{R}_1^T \boldsymbol{R}_2^T \ldots \boldsymbol{R}_k^T \boldsymbol{A}_k \boldsymbol{R}_k \ldots \boldsymbol{R}_2 \boldsymbol{R}_1.

Step n

Procedendo in questo modo, è chiaro che dopo n passi otterremo la seguente decomposizione della matrice \boldsymbol{A} originale:

\boldsymbol{A}= \boldsymbol{R}_1^T \boldsymbol{R}_2^T \ldots \boldsymbol{R}_n^T \boldsymbol{I}_n \boldsymbol{R}_n \boldsymbol{R}_{n-1} \ldots \boldsymbol{R}_1.

Le matrici \boldsymbol{R}_i^T per i = 1, 2, \ldots, n non sono solo triangolari inferiori, ma anche matrici di Frobenius. In particolare, ogni matrice \boldsymbol{R}_i^T è una matrice di Frobenius di indice i. Poiché il prodotto di queste matrici (e delle loro trasposte) è disposto nell’ordine corretto, esso risulta in una matrice triangolare inferiore che, nella colonna i-esima, contiene l’i-esima colonna di \boldsymbol{R}_i^T.

Possiamo definire la matrice \boldsymbol{R}^T come:

\boldsymbol{R}^T = \boldsymbol{R}_1^T \boldsymbol{R}_2^T \ldots \boldsymbol{R}_n^T,

e concludere che, ad ogni passo dell’algoritmo, stiamo effettivamente calcolando una colonna (o una riga) del fattore di Cholesky (o del suo trasposto).

Osservazione

Si osservi che la simmetria della matrice \boldsymbol{A} consente di ottimizzare il calcolo degli elementi dei vari complementi di Schur, poiché anch’essi risulteranno simmetrici in base al teorema generale. Nel metodo di fattorizzazione di Cholesky, quindi, si lavora esclusivamente sulla diagonale e sulla metà superiore (o inferiore) della matrice, mentre gli altri elementi sono noti per trasposizione.

Connessioni con la fattorizzazione di Gauss

Osservazione

La fattorizzazione di Cholesky rappresenta una variante intelligente della fattorizzazione di Gauss per matrici SPD. Procediamo decomponendo la matrice \boldsymbol{A} a blocchi nel seguente modo:

\begin{pmatrix} \alpha\beta & \boldsymbol{r}^T \\ \boldsymbol{c}& \boldsymbol{A}^\prime \end{pmatrix} = \begin{pmatrix} \alpha & \boldsymbol{0}\\ \dfrac{1}{\beta}\boldsymbol{c}& \boldsymbol{L}^\prime \end{pmatrix} \begin{pmatrix} \beta & \dfrac{1}{\alpha}\boldsymbol{r}^T \\ \boldsymbol{c}& \boldsymbol{U}^\prime \end{pmatrix},

dove \alpha\beta rappresenta l’elemento pivotale, mentre \boldsymbol{r} e \boldsymbol{c} indicano la prima riga e colonna di \boldsymbol{A}, escluse le prime voci. Supponiamo inoltre che esistano due matrici triangolari, \boldsymbol{L}^\prime e \boldsymbol{U}^\prime, tali che:

\boldsymbol{A}^\prime = \boldsymbol{L}^\prime\boldsymbol{U}^\prime + \dfrac{1}{\alpha\beta}\boldsymbol{c}\boldsymbol{r}^T.

• Se \alpha=1 e \textrm{diag}(\boldsymbol{L}^\prime)=1, otteniamo la fattorizzazione di Gauss (senza pivot) \boldsymbol{A}=\boldsymbol{L}\boldsymbol{U}.

• Se \alpha=\beta, \boldsymbol{L}^\prime=\boldsymbol{U}^\prime, e \boldsymbol{c}=\boldsymbol{r}, allora si ha la fattorizzazione di Cholesky \boldsymbol{A}=\boldsymbol{R}^T\boldsymbol{R}, con \boldsymbol{R}=\boldsymbol{U}^\prime.

In quest’ultimo caso, la matrice deve necessariamente essere simmetrica.

In entrambi i casi, la matrice \boldsymbol{A}^{\prime\prime} ottenuta durante il processo è il complemento di Schur del blocco di elementi corrispondenti alle righe e colonne i,j=2,\ldots,n di \boldsymbol{A}, calcolato rispetto all’elemento pivotale.

Osservazione

Dato che il teorema generale garantisce che l’elemento pivotale è sempre strettamente positivo, si potrebbe pensare che non è necessario operare alcun pivoting. Valgono tuttavia le stesse considerazioni fatte a proposito del pivot numerico e del pivot simbolico per l’algoritmo di Gauss. Nel caso di matrici sparse è sempre consigliabile operare un pivot simbolico per controllare la formazione dei non-zeri.

Note

1. In sintesi, il risparmio è considerevole!

2. Se \alpha fosse strettamente negativo, il ragionamento che segue avrebbe ancora senso ammettendo di calcolare la radice quadrata in campo complesso. Questo tuttavia non produrrebbe una algoritmo numerico efficiente. Se \alpha fosse nullo, il ragionamento che segue non sarebbe possibile. Il problema, in realtà, non si pone, perché dimostreremo tra un attimo che l’ipotesi SPD assicura la positività stretta di questo scalare.