Fattorizzazione LU, pivoting e Cholesky

I sistemi lineari $A\mathbf x=\mathbf b$ sono il cuore del calcolo scientifico: l'analisi di una struttura agli elementi finiti ha migliaia di incognite (gli spostamenti dei nodi), un circuito o una rete idraulica hanno una equazione per nodo, una discretizzazione di un'equazione differenziale ne ha una per punto della griglia. L'Eliminazione di GaussCon tre operazioni elementari sulle righe (scambio, moltiplicazione per uno scalare non nullo, somma di un multiplo di un'altra riga) ogni matrice si riduce a scala senza cambiare il rango; serve a calcolare ranghi, risolvere sistemi, trovare relazioni di dipendenza e matrici che riducono a scala.Eliminazione di Gauss → risolve il sistema, ma nella forma vista in algebra «lavora» sulla matrice e sul termine noto insieme: se il termine noto cambia (stessa struttura, carichi diversi) bisogna rifare tutto. La fattorizzazione LU separa il lavoro costoso (sulla matrice) da quello economico (sul termine noto).

Prerequisiti: Eliminazione di GaussCon tre operazioni elementari sulle righe (scambio, moltiplicazione per uno scalare non nullo, somma di un multiplo di un'altra riga) ogni matrice si riduce a scala senza cambiare il rango; serve a calcolare ranghi, risolvere sistemi, trovare relazioni di dipendenza e matrici che riducono a scala.Eliminazione di Gauss →, Operazioni tra matriciLe matrici m×n formano uno spazio vettoriale (somma e prodotto per scalare elemento per elemento); il prodotto righe per colonne corrisponde alla composizione di funzioni lineari, è associativo ma non commutativo; la trasposta scambia righe e colonne e (AB)^T = B^T A^T.Operazioni tra matrici →, Matrice inversaL'inversa di una matrice quadrata A è la matrice A⁻¹ con A A⁻¹ = A⁻¹ A = I; esiste se e solo se rango(A) = n e si calcola con Gauss-Jordan riducendo (A | I) fino a (I | A⁻¹).Matrice inversa → e Proprietà del determinante, Binet e LaplaceIl determinante è lineare in ogni riga, cambia segno scambiando due righe, non cambia sommando a una riga un multiplo di un'altra: così si calcola con Gauss riducendo a triangolare. Binet: det(AB) = det A · det B. Laplace: sviluppo lungo una riga o colonna con i complementi algebrici.Proprietà del determinante, Binet e Laplace →; per la perdita di cifre Aritmetica di macchina - floating point, precisione di macchina e cancellazioneUn calcolatore rappresenta solo un insieme finito di numeri F(β,t,L,U); ogni numero reale viene arrotondato con errore relativo al più u = ½β^(1−t) (u ≈ 1,1·10⁻¹⁶ in doppia precisione) e ogni operazione è esatta più un errore (1+δ), |δ| ≤ u. Conseguenze: overflow e underflow, somma non associativa, assorbimento, e soprattutto la cancellazione numerica nella differenza di numeri quasi uguali, che si evita riformulando l'espressione.Aritmetica di macchina - floating point, precisione di macchina e cancellazione →.

1. Gauss come fattorizzazione LU

Si trasforma $A$ in una matrice triangolare superiore $U$ con $n-1$ passi. Al passo $k$ si sottrae dalla riga $i>k$ la riga $k$ moltiplicata per il moltiplicatore $$l_{ik}=\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}\qquad(i=k+1,\dots,n),$$ dove $a^{(k)}_{ij}$ sono gli elementi dopo $k-1$ passi e $a_{kk}^{(k)}$ è il pivot. Ogni passo equivale a moltiplicare a sinistra per una matrice elementare $M_k$ (identità con $-l_{ik}$ nella colonna $k$ sotto la diagonale), quindi $M_{n-1}\cdots M_1A=U$. Le inverse $M_k^{-1}$ si moltiplicano senza interferenze e danno $$A=LU,\qquad L=\begin{pmatrix}1&&&\\l_{21}&1&&\\l_{31}&l_{32}&1&\\\vdots&&&\ddots\end{pmatrix},\qquad U=\begin{pmatrix}u_{11}&u_{12}&\cdots\\&u_{22}&\cdots\\&&\ddots\end{pmatrix}.$$

Procedura (risolvere $A\mathbf x=\mathbf b$ con LU).

1. Si fattorizza $A=LU$ (una volta sola, $\approx\frac23n^3$ flop).
2. Sostituzione in avanti $L\mathbf y=\mathbf b$: $y_i=b_i-\sum_{j<i}l_{ij}y_j$ (la diagonale di $L$ è 1, nessuna divisione).
3. Sostituzione all'indietro $U\mathbf x=\mathbf y$: $x_i=\dfrac{y_i-\sum_{j>i}u_{ij}x_j}{u_{ii}}$, $i=n,n-1,\dots,1$.

I passi 2 e 3 costano $n^2$ operazioni ciascuno: con la stessa matrice e $r$ termini noti diversi si fa una fattorizzazione sola e $r$ coppie di sostituzioni.

Esempio (a mano). $A=\begin{pmatrix}2&1&1\\4&3&3\\8&7&9\end{pmatrix}$, $\mathbf b=(4,10,24)^T$ (cioè $\mathbf x=(1,1,1)^T$).

• Passo 1, pivot $2$: $l_{21}=\frac42=2$, $l_{31}=\frac82=4$. Riga 2 $-2\cdot$ riga 1: $(0,\,3-2,\,3-2)=(0,1,1)$. Riga 3 $-4\cdot$ riga 1: $(0,\,7-4,\,9-4)=(0,3,5)$.
• Passo 2, pivot $1$: $l_{32}=\frac31=3$. Riga 3 $-3\cdot$ riga 2: $(0,0,5-3)=(0,0,2)$.

$$L=\begin{pmatrix}1&0&0\\2&1&0\\4&3&1\end{pmatrix},\qquad U=\begin{pmatrix}2&1&1\\0&1&1\\0&0&2\end{pmatrix}.$$ Verifica del prodotto: riga 3 di $LU$ = $4\cdot(2,1,1)+3\cdot(0,1,1)+1\cdot(0,0,2)=(8,\,4+3,\,4+3+2)=(8,7,9)$ ✓. Sostituzione in avanti: $y_1=4$; $y_2=10-2\cdot4=2$; $y_3=24-4\cdot4-3\cdot2=2$. All'indietro: $x_3=\frac22=1$; $x_2=\frac{2-1\cdot1}1=1$; $x_1=\frac{4-1\cdot1-1\cdot1}2=1$ ✓. Il determinante è $\det A=\det L\cdot\det U=1\cdot(2\cdot1\cdot2)=4$.

Quando esiste. Ogni minore principale di testa $A_k$ (le prime $k$ righe e colonne) si fattorizza come $A_k=L_kU_k$, quindi $\det A_k=u_{11}\cdots u_{kk}$.

Teorema. Se $\det A_k\ne0$ per $k=1,\dots,n-1$, esiste un'unica fattorizzazione $A=LU$ con $L$ a diagonale unitaria. I pivot sono $u_{kk}=\dfrac{\det A_k}{\det A_{k-1}}$.

Se un pivot si annulla l'algoritmo si ferma anche per una matrice invertibile: $\begin{pmatrix}0&1\\1&0\end{pmatrix}$ non ha LU ($a_{11}=0$). Serve scambiare le righe.

Costo. Al passo $k$ si calcolano $n-k$ moltiplicatori e si aggiornano $(n-k)^2$ elementi: circa $(n-k)^2$ moltiplicazioni e altrettante addizioni. Sommando, $\sum_{k=1}^{n-1}(n-k)^2\approx\frac{n^3}3$ moltiplicazioni e $\frac{n^3}3$ addizioni: $\frac23n^3$ flop. Per $n=1000$: $6{,}7\cdot10^8$ flop per la fattorizzazione e $2\cdot10^6$ per ogni risoluzione successiva. Lo stesso sistema con la formula di Cramer costerebbe $\approx(n+1)!$ operazioni: impraticabile già per $n=20$.

Non calcolare mai l'inversa per risolvere un sistema: $A^{-1}$ costa $\approx2n^3$ flop (si risolvono $n$ sistemi) ed è meno accurata di $LU$ con sostituzioni.

2. Instabilità e pivoting

Il moltiplicatore $l_{ik}=\frac{a_{ik}}{a_{kk}}$ è grande quando il pivot è piccolo. In aritmetica di macchina (Aritmetica di macchina - floating point, precisione di macchina e cancellazioneUn calcolatore rappresenta solo un insieme finito di numeri F(β,t,L,U); ogni numero reale viene arrotondato con errore relativo al più u = ½β^(1−t) (u ≈ 1,1·10⁻¹⁶ in doppia precisione) e ogni operazione è esatta più un errore (1+δ), |δ| ≤ u. Conseguenze: overflow e underflow, somma non associativa, assorbimento, e soprattutto la cancellazione numerica nella differenza di numeri quasi uguali, che si evita riformulando l'espressione.Aritmetica di macchina - floating point, precisione di macchina e cancellazione →) questo distrugge l'informazione dei numeri «ordinari».

Esempio. $\varepsilon=10^{-20}$ in doppia precisione ($\approx16$ cifre decimali), sistema $$\begin{pmatrix}\varepsilon&1\\1&1\end{pmatrix}\begin{pmatrix}x_1\\x_2\end{pmatrix}=\begin{pmatrix}1\\2\end{pmatrix},\qquad\text{soluzione esatta }x_1=\tfrac1{1-\varepsilon}\approx1,\ x_2=\tfrac{1-2\varepsilon}{1-\varepsilon}\approx1 .$$ Senza scambi: $l_{21}=\frac1\varepsilon=10^{20}$, $u_{22}=1-10^{20}\cdot1=-10^{20}$ (il $1$ viene perso: $10^{20}+1$ in doppia precisione vale $10^{20}$), $y_2=2-10^{20}\cdot1=-10^{20}$, $x_2=\frac{y_2}{u_{22}}=1$, $x_1=\frac{1-1\cdot1}{\varepsilon}=\mathbf 0$. Il risultato $x_1=0$ è completamente sbagliato (esatto $\approx1$), anche se il sistema è ben condizionato. Scambiando le righe, il pivot è $1$, $l_{21}=\varepsilon$, $u_{22}=1-\varepsilon\approx1$, $y_2=1-2\varepsilon\approx1$ e si ottiene $x=(1,1)$ corretto.

Pivoting parziale. Al passo $k$ si cerca fra le righe $i\ge k$ quella con $|a_{ik}^{(k)}|$ massimo e la si scambia con la riga $k$. Il risultato è $PA=LU$, con $P$ matrice di permutazione e $|l_{ik}|\le1$.

Per ogni matrice invertibile esiste almeno un elemento non nullo in colonna sotto (o sulla) diagonale, quindi il procedimento non si blocca mai. Si risolve $LU\mathbf x=P\mathbf b$: si permuta il termine noto e poi si procede come prima. Il determinante è $\det A=(-1)^s\prod u_{kk}$, con $s$ numero di scambi.

Fattore di crescita. $\rho=\dfrac{\max_{ijk}|a_{ij}^{(k)}|}{\max_{ij}|a_{ij}|}$. Con pivoting parziale $\rho\le2^{n-1}$ (limite teorico raggiunto solo da matrici costruite apposta); in pratica $\rho$ è piccolo e l'algoritmo è stabile all'indietro: la soluzione calcolata risolve esattamente $(A+\delta A)\hat{\mathbf x}=\mathbf b$ con $\frac{\|\delta A\|}{\|A\|}\lesssim n\rho\,u$, $u$ precisione di macchina. L'accuratezza finale dipende poi dal condizionamento di $A$ (Norme di matrici e condizionamento dei sistemi lineariNorme di vettore $|x|1=\sum|x_i|$, $|x|2=\sqrt{\sum x_i^2}$, $|x|\infty=\max|x_i|$. Norma di matrice indotta $|A|=\max{x\ne0}\frac{|Ax|}{|x|}$: $|A|1$ = massima somma per colonne, $|A|\infty$ = massima somma per righe, $|A|2=\sqrt{\lambda{\max}(A^TA)}$ (= $\rho(A)$ se $A$ simmetrica); valgono $|Ax|\le|A||x|$, $|AB|\le|A||B|$, $\rho(A)\le|A|$. Numero di condizionamento $\kappa(A)=|A||A^{-1}|\ge1$ ($\kappa_2=\frac{\sigma_{\max}}{\sigma_{\min}}$, per SPD $\frac{\lambda_{\max}}{\lambda_{\min}}$, per ortogonali 1). Se $A\mathbf x=\mathbf b$ e $A(\mathbf x+\delta\mathbf x)=\mathbf b+\delta\mathbf b$: $\frac{|\delta\mathbf x|}{|\mathbf x|}\le\kappa(A)\frac{|\delta\mathbf b|}{|\mathbf b|}$ (D); con residuo $\mathbf r=\mathbf b-A\hat{\mathbf x}$: $\frac1\kappa\frac{|\mathbf r|}{|\mathbf b|}\le\frac{|\mathbf x-\hat{\mathbf x}|}{|\mathbf x|}\le\kappa\frac{|\mathbf r|}{|\mathbf b|}$. Si perdono circa $\log_{10}\kappa$ cifre; un residuo piccolo non garantisce un errore piccolo se $\kappa$ è grande.Norme di matrici e condizionamento dei sistemi lineari →).

Esempio (a mano, con pivoting). $A=\begin{pmatrix}1&3&2\\4&2&1\\2&2&3\end{pmatrix}$, $\mathbf b=(5,7,3)^T$ ($\mathbf x=(1,2,-1)^T$).

• Passo 1: il massimo in valore assoluto della colonna 1 è $4$ (riga 2): si scambiano le righe 1 e 2 ($P$ scambia 1 e 2): $\begin{pmatrix}4&2&1\\1&3&2\\2&2&3\end{pmatrix}$. $l_{21}=\frac14$, $l_{31}=\frac24=\frac12$. Riga 2 $-\frac14$ riga 1: $(0,\ 3-\frac12,\ 2-\frac14)=(0,2{,}5,1{,}75)$. Riga 3 $-\frac12$ riga 1: $(0,\ 2-1,\ 3-\frac12)=(0,1,2{,}5)$.
• Passo 2: fra $2{,}5$ e $1$ vince $2{,}5$: nessuno scambio. $l_{32}=\frac1{2{,}5}=0{,}4$. Riga 3 $-0{,}4$ riga 2: $(0,0,\ 2{,}5-0{,}4\cdot1{,}75)=(0,0,1{,}8)$.

$$P=\begin{pmatrix}0&1&0\\1&0&0\\0&0&1\end{pmatrix},\quad L=\begin{pmatrix}1&0&0\\0{,}25&1&0\\0{,}5&0{,}4&1\end{pmatrix},\quad U=\begin{pmatrix}4&2&1\\0&2{,}5&1{,}75\\0&0&1{,}8\end{pmatrix}.$$ Termine noto permutato: $P\mathbf b=(7,5,3)^T$. Avanti: $y_1=7$; $y_2=5-0{,}25\cdot7=3{,}25$; $y_3=3-0{,}5\cdot7-0{,}4\cdot3{,}25=-1{,}8$. Indietro: $x_3=\frac{-1{,}8}{1{,}8}=-1$; $x_2=\frac{3{,}25-1{,}75\cdot(-1)}{2{,}5}=2$; $x_1=\frac{7-2\cdot2-1\cdot(-1)}4=1$ ✓. Determinante: $\det A=(-1)^1\cdot4\cdot2{,}5\cdot1{,}8=-18$ (un solo scambio).

3. Matrici simmetriche definite positive: Cholesky

Molte matrici dell'ingegneria sono simmetriche definite positive (SPD): la matrice di rigidezza di una struttura, la matrice delle equazioni normali dei minimi quadrati, la matrice di una conduzione termica discretizzata. Una matrice $A=A^T$ è definita positiva se $\mathbf v^TA\mathbf v>0$ per ogni $\mathbf v\ne\mathbf 0$ (si veda Forme quadratiche e criteri di definitezzaUna forma quadratica q(h) = h^T M h (M simmetrica) è definita positiva se q>0 per h ≠ 0, semidefinita positiva se q ≥ 0, definita/semidefinita negativa con segni opposti, indefinita se cambia segno. In due variabili, con M = [[a,b],[b,c]]: definita positiva ⇔ det M > 0 e a > 0, definita negativa ⇔ det M > 0 e a < 0, indefinita ⇔ det M < 0 (D, completando il quadrato). In generale: segno degli autovalori di M (D, teorema spettrale): tutti > 0 definita positiva, tutti < 0 definita negativa, di segno diverso indefinita; inoltre q(h) ≥ λ_min |h|². Per n = 3 si usano i minori principali di Sylvester. È lo strumento per classificare i punti critici.Forme quadratiche e criteri di definitezza →): equivalentemente tutti gli autovalori sono positivi, oppure tutti i minori principali di testa $\det A_k$ sono positivi.

Teorema (Cholesky). Se $A$ è SPD esiste un'unica matrice triangolare inferiore $L$ con elementi diagonali positivi tale che $A=LL^T$.

Uguagliando gli elementi di $A=LL^T$ (con $j\le i$) $a_{ij}=\sum_{k=1}^jl_{ik}l_{jk}$ si ricavano per colonne $$\boxed{\,l_{jj}=\sqrt{a_{jj}-\sum_{k=1}^{j-1}l_{jk}^2},\qquad l_{ij}=\frac1{l_{jj}}\Big(a_{ij}-\sum_{k=1}^{j-1}l_{ik}l_{jk}\Big)\ (i>j)\,}$$ Il radicando $a_{jj}-\sum l_{jk}^2$ è positivo proprio perché $A$ è definita positiva (è il rapporto $\frac{\det A_j}{\det A_{j-1}}$): se compare un radicando $\le0$, la matrice non è SPD, ed è anche il test più economico per controllarlo.

Esempio. $A=\begin{pmatrix}4&2&2\\2&5&3\\2&3&6\end{pmatrix}$. Minori di testa: $4$, $4\cdot5-4=16$, $\det A=64$: tutti positivi, A è SPD.

• Colonna 1: $l_{11}=\sqrt4=2$, $l_{21}=\frac22=1$, $l_{31}=\frac22=1$.
• Colonna 2: $l_{22}=\sqrt{5-1^2}=2$, $l_{32}=\frac{3-l_{31}l_{21}}{l_{22}}=\frac{3-1\cdot1}2=1$.
• Colonna 3: $l_{33}=\sqrt{6-1^2-1^2}=\sqrt4=2$.

$$L=\begin{pmatrix}2&0&0\\1&2&0\\1&1&2\end{pmatrix},\qquad\det A=(l_{11}l_{22}l_{33})^2=8^2=64 .$$ Per risolvere $A\mathbf x=\mathbf b$ si usano $L\mathbf y=\mathbf b$ e $L^T\mathbf x=\mathbf y$.

Vantaggi rispetto a LU. Costo $\frac13n^3$ flop (metà, si memorizza solo $L$); nessun pivoting necessario e la stabilità è garantita ($|l_{jk}|^2\le a_{jj}$: gli elementi non crescono); le matrici a banda mantengono la banda. Si può evitare la radice con la forma $A=LDL^T$ ($L$ a diagonale unitaria, $D$ diagonale). Nel caso di una matrice tridiagonale (catena di molle, differenze finite 1D) tutto costa $O(n)$.

4. Implementazione

import numpy as np

def lu_fattorizza(A, pivoting=True):
    """PA = LU con pivoting parziale (o A = LU senza). Restituisce P, L, U."""
    A = np.array(A, dtype=float)
    n = len(A)
    P = np.eye(n)
    L = np.eye(n)
    for k in range(n - 1):
        if pivoting:
            r = k + np.argmax(abs(A[k:, k]))          # riga del pivot massimo in colonna
            if r != k:
                A[[k, r]] = A[[r, k]]
                P[[k, r]] = P[[r, k]]
                L[[k, r], :k] = L[[r, k], :k]
        for i in range(k + 1, n):
            L[i, k] = A[i, k] / A[k, k]               # moltiplicatore
            A[i, k:] -= L[i, k] * A[k, k:]            # riga i - l_ik * riga k
    return P, L, np.triu(A)

def risolvi(P, L, U, b):
    n = len(b)
    y = np.zeros(n)
    c = P @ b
    for i in range(n):                                # sostituzione in avanti: Ly = Pb
        y[i] = c[i] - L[i, :i] @ y[:i]
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):                    # sostituzione all'indietro: Ux = y
        x[i] = (y[i] - U[i, i + 1:] @ x[i + 1:]) / U[i, i]
    return x

A = np.array([[1, 3, 2], [4, 2, 1], [2, 2, 3]])
b = A @ np.array([1, 2, -1])
P, L, U = lu_fattorizza(A)
print("L =\n", L, "\nU =\n", U, "\nP =\n", P)
print("x =", risolvi(P, L, U, b), " det(A) =", round(np.linalg.det(A), 6), " det(P)*prod(U_ii) =", round(np.linalg.det(P) * np.prod(np.diag(U)), 6))

# pivot piccolo: eps = 1e-20, soluzione esatta quasi (1, 1)
eps = 1e-20
B = np.array([[eps, 1.0], [1.0, 1.0]])
c = np.array([1.0, 2.0])
for piv in (False, True):
    print("pivoting", piv, "->", risolvi(*lu_fattorizza(B, piv), c))

# Cholesky
M = np.array([[4, 2, 2], [2, 5, 3], [2, 3, 6]], dtype=float)
print("chol =\n", np.linalg.cholesky(M))

Output:

L =
 [[1.   0.   0.  ]
 [0.25 1.   0.  ]
 [0.5  0.4  1.  ]] 
U =
 [[4.   2.   1.  ]
 [0.   2.5  1.75]
 [0.   0.   1.8 ]] 
P =
 [[0. 1. 0.]
 [1. 0. 0.]
 [0. 0. 1.]]
x = [ 1.  2. -1.]  det(A) = -18.0  det(P)*prod(U_ii) = -18.0
pivoting False -> [0. 1.]
pivoting True -> [1. 1.]
chol =
 [[2. 0. 0.]
 [1. 2. 0.]
 [1. 1. 2.]]

Equivalenti MATLAB: [L,U,P] = lu(A) (pivoting parziale, $PA=LU$); x = A\b (sceglie da solo LU, Cholesky o altro in base alla matrice: non calcola inv(A)*b); R = chol(A) restituisce il fattore superiore $R=L^T$ con $A=R^TR$; det(A). In numpy np.linalg.solve(A, b) usa LAPACK con pivoting parziale.

Errori da evitare

• Dimenticare di scambiare anche le righe già calcolate di $L$ (e il termine noto) quando si cambia pivot a metà procedimento: si scambia la riga intera, compresi i moltiplicatori già memorizzati.
• Segno dei moltiplicatori: in $L$ compaiono i moltiplicatori con segno $+$ (sono quelli sottratti).
• Usare $\det A=\prod u_{kk}$ con pivoting: serve il fattore $(-1)^s$.
• Applicare Cholesky senza controllare che $A$ sia SPD: la simmetria da sola non basta (esempio: $\begin{pmatrix}1&2\\2&1\end{pmatrix}$ ha $l_{22}=\sqrt{1-4}$, autovalori $3,-1$).
• Calcolare $A^{-1}\mathbf b$: più lento e meno accurato di LU con sostituzioni.
• Confondere $R$ (MATLAB chol) con $L$: $R=L^T$.

Esercizi svolti: Esercizio - fattorizzazione LU di una matrice 3x3, sistema e determinante, Esercizio - pivoting parziale su una matrice 3x3 con scambio di righe, Esercizio - Cholesky della matrice di rigidezza di una catena di tre masse, Esercizio - pivot piccolo in aritmetica a quattro cifre e costo della fattorizzazione.

Domande d'esame

1. Quando una matrice ammette la fattorizzazione LU senza scambi? Dare un esempio di matrice invertibile che non la ammette. Traccia. $\det A_k\ne0$ per $k=1,\dots,n-1$ (pivot $u_{kk}=\det A_k/\det A_{k-1}$); $\begin{pmatrix}0&1\\1&0\end{pmatrix}$.
2. Spiegare perché si fa il pivoting e che cosa garantisce. Traccia. Pivot piccolo $\Rightarrow$ moltiplicatori grandi $\Rightarrow$ perdita di cifre ($\varepsilon=10^{-20}$ dà $x_1=0$); pivoting parziale: $|l_{ik}|\le1$, $PA=LU$, esiste per ogni $A$ invertibile.
3. Costo di LU e delle sostituzioni; come si risolvono 100 sistemi con la stessa matrice. Traccia. $\frac23n^3$ una volta, poi $2n^2$ per sistema.
4. Enunciare il teorema di Cholesky e ricavare $l_{jj}$ e $l_{ij}$. Traccia. Uguagliare $a_{ij}=\sum_{k\le j}l_{ik}l_{jk}$; $a_{jj}=\sum_{k<j}l_{jk}^2+l_{jj}^2$; radicando positivo $\iff$ SPD.