Algebra lineare (3): calcolo differenziale con le matrici

Nelle prime due parti di questa serie abbiamo costruito gli spazi vettoriali e studiato le matrici. Restava fuori un pezzo: cosa succede quando si deriva rispetto a un vettore o a una matrice.

È il pezzo che serve davvero se si scrive codice di ottimizzazione. Ogni algoritmo di addestramento, dalla regressione lineare alla discesa del gradiente in una rete neurale, calcola la derivata di una funzione scalare — la loss — rispetto a un insieme di parametri organizzati in matrici. Le regole che seguono sono quelle che permettono di farlo senza scendere agli indici ogni volta, e quello che si scopre è che ricalcano fedelmente le regole del calcolo a una variabile: dove nel caso scalare c'è \(ax\), qui c'è \(b^T x\); dove c'è \(ax^2\), qui c'è \(x^T A x\); e le derivate si corrispondono una per una.

1. Il gradiente rispetto a una matrice

La definizione è la più naturale possibile: si deriva rispetto a ogni singolo elemento e si ripongono i risultati nella stessa posizione in cui stava l'elemento derivato.

DefinizioneGradiente. Sia \(f : \mathbb{R}^{m \times n} \to \mathbb{R}\) una funzione che prende in input una matrice di dimensione \(m \times n\) e restituisce un valore reale. Allora il gradiente di \(f\) rispetto ad \(A \in \mathbb{R}^{m \times n}\) è la matrice delle derivate parziali definita come: $$\nabla_A f(A) = \begin{bmatrix} \frac{\partial f(A)}{\partial A_{11}} & \frac{\partial f(A)}{\partial A_{12}} & \cdots & \frac{\partial f(A)}{\partial A_{1n}} \\[.4em] \frac{\partial f(A)}{\partial A_{21}} & \frac{\partial f(A)}{\partial A_{22}} & \cdots & \frac{\partial f(A)}{\partial A_{2n}} \\[.4em] \vdots & \vdots & \ddots & \vdots \\[.4em] \frac{\partial f(A)}{\partial A_{m1}} & \frac{\partial f(A)}{\partial A_{m2}} & \cdots & \frac{\partial f(A)}{\partial A_{mn}} \end{bmatrix} \ \in \mathbb{R}^{m \times n}$$ ossia una matrice \(m \times n\) con \(\big(\nabla_A f(A)\big)_{ij} = \frac{\partial f(A)}{\partial A_{ij}}\).

Il gradiente ha dunque sempre la stessa forma dell'oggetto rispetto a cui si deriva. Nel caso particolare in cui l'argomento sia un vettore \(x \in \mathbb{R}^n\), il gradiente è a sua volta un vettore colonna:

$$\nabla_x f(x) = \begin{bmatrix} \frac{\partial f(x)}{\partial x_1} \\[.4em] \frac{\partial f(x)}{\partial x_2} \\[.4em] \vdots \\[.4em] \frac{\partial f(x)}{\partial x_n} \end{bmatrix}$$

Questa corrispondenza di forma non è un dettaglio estetico: è ciò che rende possibile scrivere l'aggiornamento dei parametri come \(\theta \leftarrow \theta - \eta\, \nabla_\theta L\), sommando due oggetti della stessa dimensione.

Attenzione. Il gradiente è definito solamente se \(f\) ha valori reali. Non possiamo per esempio calcolare il gradiente di \(f(x) = Ax\) con \(x \in \mathbb{R}^n\), \(A \in \mathbb{R}^{n \times n}\), perché il risultato è un vettore e non uno scalare. Per funzioni a valori vettoriali serve la matrice jacobiana, che è un altro oggetto.

ProposizioneProprietà del gradiente. Derivano direttamente dalle proprietà delle derivate parziali:

• \(\nabla_x\big(f(x) + g(x)\big) = \nabla_x f(x) + \nabla_x g(x)\)
• \(\forall t \in \mathbb{R} : \nabla_x\big(t\, f(x)\big) = t\, \nabla_x f(x)\)

Il gradiente è cioè un operatore lineare, il che permette di spezzare il calcolo di espressioni complicate in pezzi trattabili separatamente — cosa che faremo tra poco derivando i minimi quadrati.

2. L'hessiana

Passando alle derivate seconde, le combinazioni possibili di due indici diventano \(n^2\) e il risultato è di nuovo una matrice, anche quando si parte da un vettore.

DefinizioneHessiana. Sia \(f : \mathbb{R}^n \to \mathbb{R}\). Si definisce hessiana di \(f\) rispetto a \(x\), denotata con \(\nabla_x^2 f(x)\) o semplicemente \(H\), la matrice \(n \times n\): $$\nabla_x^2 f(x) = \begin{bmatrix} \frac{\partial^2 f(x)}{\partial x_1^2} & \frac{\partial^2 f(x)}{\partial x_1 \partial x_2} & \cdots & \frac{\partial^2 f(x)}{\partial x_1 \partial x_n} \\[.4em] \frac{\partial^2 f(x)}{\partial x_2 \partial x_1} & \frac{\partial^2 f(x)}{\partial x_2^2} & \cdots & \frac{\partial^2 f(x)}{\partial x_2 \partial x_n} \\[.4em] \vdots & \vdots & \ddots & \vdots \\[.4em] \frac{\partial^2 f(x)}{\partial x_n \partial x_1} & \frac{\partial^2 f(x)}{\partial x_n \partial x_2} & \cdots & \frac{\partial^2 f(x)}{\partial x_n^2} \end{bmatrix} \ \in \mathbb{R}^{n \times n}$$ ossia \(\big(\nabla_x^2 f(x)\big)_{ij} = \frac{\partial^2 f(x)}{\partial x_i \partial x_j}\).

L'hessiana è sempre simmetrica, perché per il teorema di Schwarz l'ordine di derivazione non conta:

$$\frac{\partial^2 f(x)}{\partial x_i \partial x_j} = \frac{\partial^2 f(x)}{\partial x_j \partial x_i}$$

Ed è una simmetria fortunata, perché è esattamente l'ipotesi sotto cui valgono i teoremi più forti visti nella seconda parte: autovalori reali, base ortonormale di autovettori, caratterizzazione della definitezza tramite il segno degli autovalori. Il legame fra definitezza dell'hessiana e natura di un punto critico è quello che ci si aspetta dal caso a una variabile.

ProposizioneCondizione sufficiente per minimi e massimi relativi. Sia \(A\) un aperto di \(\mathbb{R}^n\), \(f \in \mathbf{C}^2(A)\), \(\dot x \in A\). Sia \(H_f(\dot x)\) la matrice hessiana di \(f\) in \(\dot x\). Allora:

• se \(df(\dot x) = 0\) e \(H_f(\dot x)\) è definita positiva, allora \(\dot x\) è un minimo relativo per \(f\);
• se \(df(\dot x) = 0\) e \(H_f(\dot x)\) è definita negativa, allora \(\dot x\) è un massimo relativo per \(f\).

dove \(df(\dot x)\) è il differenziale di \(f\) calcolato in \(\dot x\).

È l'analogo esatto del criterio della derivata seconda: \(f''(x) \gt 0\) diventa "hessiana definita positiva". Quando l'hessiana è indefinita — autovalori di entrambi i segni — il punto critico non è né un minimo né un massimo, ma una sella: la funzione sale lungo certe direzioni e scende lungo altre. Sono i punti che rendono difficile l'ottimizzazione non convessa, e la seconda parte ci ha già detto come riconoscerli.

Anche l'hessiana è definita solo se \(f\) ha valori reali. Nel caso \(n = 1\) si riduce a \(\frac{\partial^2 f(x)}{\partial x^2}\), che è un numero.

3. Gradiente e hessiana delle forme lineari e quadratiche

Passiamo al calcolo effettivo. Ricaviamo le derivate dei due casi che compaiono praticamente ovunque, e verifichiamo la corrispondenza con il caso scalare.

3.1 Forma lineare

Per \(x \in \mathbb{R}^n\) sia \(f(x) = b^T x\), con \(b \in \mathbb{R}^n\) vettore noto. Esplicitando,

$$f(x) = \sum_{i=1}^n b_i x_i$$

e derivando rispetto alla generica componente \(x_k\), tutti i termini della somma spariscono tranne quello che contiene effettivamente \(x_k\):

$$\frac{\partial f(x)}{\partial x_k} = \frac{\partial}{\partial x_k} \sum_{i=1}^n b_i x_i = b_k$$

Raccogliendo le \(n\) derivate parziali in un vettore si ottiene

$$\nabla_x\, b^T x = b$$

che è l'analogo diretto della derivata in una variabile \(\frac{d}{dx}(ax) = a\).

3.2 Forma quadratica: il gradiente

Consideriamo ora \(f(x) = x^T A x\) con \(A \in \mathbf{S}^n\) simmetrica — ipotesi legittima, come abbiamo visto nella seconda parte, perché la parte antisimmetrica non contribuisce alla forma quadratica. Ricordiamo che

$$f(x) = \sum_{i=1}^n \sum_{j=1}^n A_{ij} x_i x_j$$

Per calcolare la derivata parziale rispetto a \(x_k\) conviene separare i termini secondo che contengano \(x_k\) una volta, due volte o nessuna:

$$\frac{\partial f(x)}{\partial x_k} = \frac{\partial}{\partial x_k}\left[ \sum_{i \neq k}\sum_{j \neq k} A_{ij}x_i x_j + \sum_{i \neq k} A_{ik}x_i x_k + \sum_{j \neq k} A_{kj}x_k x_j + A_{kk}x_k^2 \right]$$

Il primo blocco non dipende da \(x_k\) e si annulla; il secondo e il terzo sono lineari in \(x_k\); l'ultimo è quadratico e dà un fattore 2. Restano

$$= \sum_{i \neq k} A_{ik}x_i + \sum_{j \neq k} A_{kj}x_j + 2 A_{kk}x_k = \sum_{i=1}^n A_{ik}x_i + \sum_{j=1}^n A_{kj}x_j \ \overset{(A \text{ simmetrica})}{=}\ 2\sum_{i=1}^n A_{ki}x_i$$

dove nel passaggio centrale i due termini \(A_{kk}x_k\) mancanti sono stati riassorbiti nelle sommatorie complete. L'ultimo membro è il prodotto scalare della \(k\)-esima riga di \(A\) per \(x\), moltiplicato per 2; considerando tutte le \(n\) righe si ottiene

$$\nabla_x\, x^T A x = 2Ax$$

di nuovo in perfetta analogia col caso scalare, dove \(\frac{d}{dx}(ax^2) = 2ax\).

3.3 Forma quadratica: l'hessiana

Derivando una seconda volta il risultato appena ottenuto:

$$\frac{\partial^2 f(x)}{\partial x_k \partial x_j} = \frac{\partial}{\partial x_k}\left[ 2\sum_{i=1}^n A_{ji}x_i \right] = 2A_{jk} = 2A_{kj}$$

dove l'ultima uguaglianza usa ancora la simmetria di \(A\). Dunque

$$\nabla_x^2\, x^T A x = 2A$$

analogamente a \(\frac{d^2}{dx^2}(ax^2) = 2a\). Riassumendo le tre formule:

Caso	Funzione	Derivata
Forma lineare, gradiente	\(f(x) = b^T x\)	\(\nabla_x\, b^T x = b\)
Forma quadratica, gradiente	\(f(x) = x^T A x\)	\(\nabla_x\, x^T A x = 2Ax\)  (\(A\) simmetrica)
Forma quadratica, hessiana	\(f(x) = x^T A x\)	\(\nabla_x^2\, x^T A x = 2A\)  (\(A\) simmetrica)

Vale la pena notare che l'hessiana di una forma quadratica è costante: non dipende da \(x\). Una forma quadratica con \(A\) definita positiva è quindi convessa ovunque, e ha un unico minimo globale. È la ragione per cui i problemi ai minimi quadrati si risolvono in forma chiusa, mentre l'addestramento di una rete neurale — dove la loss non è quadratica — richiede metodi iterativi.

4. Le equazioni dei minimi quadrati

Con questi strumenti possiamo derivare uno dei risultati più usati in assoluto, e chiudere il conto lasciato aperto nella seconda parte.

DefinizioneEquazioni dei minimi quadrati. Le equazioni dei minimi quadrati sono usate per trovare il migliore adattamento fra un insieme di dati osservati e un modello matematico \(f(x, \theta)\), dove \(\theta\) sono i parametri del modello. Si cercano i valori di \(\theta\) che minimizzano la somma dei quadrati delle differenze fra i dati osservati e quelli predetti: $$\min_\theta \sum_{i=1}^n \big[\, y_i - f(x_i, \theta) \,\big]^2$$

ProposizioneSoluzione in forma chiusa. Sia \(A \in \mathbb{R}^{m \times n}\) di pieno rango e \(b \in \mathbb{R}^m\) con \(b \notin \mathbf{R}(A)\). Allora il vettore \(x\) che minimizza \(\lVert Ax - b \rVert_2^2\) è $$x = (A^T A)^{-1}A^T b$$

Dimostrazione

Poiché \(b\) non appartiene al range di \(A\) — cioè, ricordando la seconda parte, \(b\) non è combinazione lineare delle colonne di \(A\) — non esiste alcun \(x\) con \(Ax = b\). Cerchiamo allora l'\(x\) che rende \(Ax\) il più vicino possibile a \(b\), misurando con la norma euclidea.

Usando \(\lVert v \rVert_2^2 = v^T v\) ed espandendo:

$$\lVert Ax - b \rVert_2^2 = (Ax-b)^T(Ax-b) = (x^T A^T - b^T)(Ax - b) = x^T A^T A x - x^T A^T b - b^T A x + b^T b$$

I due termini centrali sono scalari, e uno è la trasposta dell'altro: \(x^T A^T b = (b^T A x)^T = b^T A x\). Si sommano quindi in un unico termine:

$$= x^T A^T A x - 2\, b^T A x + b^T b$$

Ora prendiamo il gradiente rispetto a \(x\), sfruttando la linearità e riconoscendo in ciascun addendo una delle forme derivate al punto 3 — il primo è una forma quadratica nella matrice \(A^T A\) (che è simmetrica), il secondo è una forma lineare nel vettore \(A^T b\), il terzo è costante in \(x\):

$$\nabla_x\big(x^T A^T A x - 2 b^T A x + b^T b\big) = \underbrace{\nabla_x\, x^T (A^TA) x}_{= \,2A^TAx} \ -\ \underbrace{2\,\nabla_x\, (A^Tb)^T x}_{=\, 2A^Tb} \ +\ \underbrace{\nabla_x\, b^Tb}_{=\,0} = 2A^T A x - 2A^T b$$

Ponendo il gradiente uguale a zero:

$$2A^T A x - 2A^T b = 0 \quad \Rightarrow \quad A^T A x = A^T b \quad \Rightarrow \quad x = (A^T A)^{-1}A^T b$$

L'inversa esiste perché \(A\) ha pieno rango, e quindi \(A^T A\) è definita positiva per la proposizione sulla matrice di Gram vista nella seconda parte. Essendo \(A^T A\) definita positiva, l'hessiana \(2A^TA\) lo è pure e il punto critico trovato è effettivamente un minimo. \(\blacksquare\)

L'equazione intermedia \(A^T A x = A^T b\) è quella che si chiama sistema delle equazioni normali, ed è la forma in cui il problema viene effettivamente risolto in pratica — senza mai invertire esplicitamente \(A^T A\).

Coerenza con la proiezione. Nella seconda parte avevamo enunciato che la proiezione di \(b\) sul range di \(A\) vale \(A(A^TA)^{-1}A^T b\). Ora abbiamo trovato che il minimizzante è \(x = (A^TA)^{-1}A^Tb\); ma il punto del range più vicino a \(b\) è \(Ax\), cioè \(A(A^TA)^{-1}A^Tb\). I due risultati coincidono, come devono: risolvere ai minimi quadrati è proiettare ortogonalmente sul sottospazio generato dalle colonne.

5. Gradiente del determinante

Un caso meno frequente ma istruttivo, perché mostra come le proprietà algebriche viste nella seconda parte si traducano direttamente in formule di derivazione. Cominciamo con una riscrittura dell'identità che lega inversa e aggiunta.

Lemma \(\mathrm{adj}(A) = |A|\, A^{-1}\).

Dimostrazione

Dalla seconda parte sappiamo che \(A^{-1} = \frac{\mathrm{adj}(A)}{|A|}\). Moltiplicando ambo i membri per lo scalare \(|A|\) si ottiene immediatamente la tesi. \(\blacksquare\)

ProposizioneGradiente del determinante. Data \(f : \mathbb{R}^{n \times n} \to \mathbb{R}\) con \(f(A) = |A|\), si ha $$\nabla_A |A| = |A|\, A^{-T}$$

Dimostrazione

Dalla definizione di determinante per sviluppo di Laplace, scelta una colonna \(j\) qualsiasi e ricordando la definizione di complemento algebrico \(\Delta_{ij}\):

$$|A| = \sum_{i=1}^n A_{ij}\Delta_{ij}$$

La chiave è che il complemento algebrico \(\Delta_{kl}\) si ottiene cancellando la riga \(k\) e la colonna \(l\), quindi non contiene l'elemento \(A_{kl}\). Rispetto ad \(A_{kl}\), l'espressione è dunque lineare, e derivando resta il solo coefficiente:

$$\frac{\partial |A|}{\partial A_{kl}} = \Delta_{kl}$$

Il gradiente è quindi la matrice dei complementi algebrici, che è la trasposta dell'aggiunta classica (si ricordi l'inversione degli indici nella definizione di \(\mathrm{adj}\)):

$$\nabla_A |A| = \big(\mathrm{adj}(A)\big)^T \ \overset{\text{(lemma)}}{=}\ \big(|A|\, A^{-1}\big)^T = |A|\, A^{-T}$$

dove nell'ultimo passaggio si è portato fuori \(|A|\), essendo uno scalare. \(\blacksquare\)

ProposizioneDifferenziale di \(\log|A|\). Sia \(f : \mathbf{S}_{++}^n \to \mathbb{R}\) con \(f(A) = \log|A|\). Allora $$\nabla_A \log|A| = A^{-T} = A^{-1}$$

Dimostrazione

Restringiamo il dominio alle matrici simmetriche definite positive, così che il determinante sia positivo e il logaritmo restituisca valori reali. Applicando la regola di derivazione delle funzioni composte:

$$\frac{\partial \log|A|}{\partial A_{ij}} = \frac{\partial \log|A|}{\partial |A|} \cdot \frac{\partial |A|}{\partial A_{ij}} = \frac{1}{|A|} \cdot \frac{\partial |A|}{\partial A_{ij}}$$

e quindi, in forma matriciale, usando il risultato precedente:

$$\nabla_A \log|A| = \frac{1}{|A|}\, \nabla_A |A| = \frac{1}{|A|}\,|A|\, A^{-T} = A^{-T} = A^{-1}$$

dove l'ultimo passaggio usa la simmetria di \(A\). \(\blacksquare\)

Il risultato è l'analogo multidimensionale di \(\frac{d}{dx}\log x = \frac{1}{x}\): il posto del reciproco è preso dall'inversa. Non è una curiosità isolata — \(\log|A|\) compare nella log-verosimiglianza di una distribuzione gaussiana multivariata, e questa derivata è ciò che serve per stimarne la matrice di covarianza.

6. Gli autovalori come soluzione di un problema di ottimizzazione

Chiudiamo la serie con un risultato che ricollega tutto. Nella seconda parte avevamo dimostrato, con un lungo conto sulle derivate parziali, che il massimo e il minimo di \(x^T A x\) sulla sfera unitaria sono autovalori di \(A\). Lo stesso fatto si ottiene ora in poche righe, usando il metodo dei moltiplicatori di Lagrange e le formule di derivazione appena ricavate.

ProposizioneAutovalori come ottimizzatori. Consideriamo il problema di ottimizzazione vincolata $$\max_{x \in \mathbb{R}^n} x^T A x \qquad \text{soggetto a} \qquad \lVert x \rVert_2^2 = 1$$ per una matrice simmetrica \(A \in \mathbb{R}^{n \times n}\). Allora gli unici punti che possono massimizzare (o minimizzare) \(x^T A x\) sono gli autovettori di \(A\).

Dimostrazione

Il metodo standard consiste nello scrivere la lagrangiana, una funzione che incorpora il vincolo di uguaglianza:

$$\mathbf{L}(x, \lambda) = x^T A x - \lambda\,(x^T x - 1)$$

dove \(\lambda\) è il moltiplicatore di Lagrange. Si dimostra che, affinché \(x^*\) sia una soluzione, il gradiente della lagrangiana deve annullarsi in \(x^*\) — condizione necessaria ma non sufficiente. Calcoliamolo usando le formule del punto 3, osservando che anche \(x^T x = x^T I x\) è una forma quadratica, con matrice \(I\):

$$\nabla_x \mathbf{L}(x, \lambda) = \nabla_x\big(x^T A x - \lambda\, x^T x\big) = 2Ax - 2\lambda x$$

Ponendo il gradiente uguale a zero:

$$2Ax - 2\lambda x = 0 \quad \Rightarrow \quad Ax - \lambda x = 0 \quad \Rightarrow \quad Ax = \lambda x$$

che è esattamente la definizione di autovettore. \(\blacksquare\)

Il moltiplicatore di Lagrange, introdotto come artificio tecnico per gestire il vincolo, si rivela essere l'autovalore. E il valore della funzione obiettivo nel punto ottimo è \(x^T A x = x^T \lambda x = \lambda \lVert x \rVert^2 = \lambda\), dato che \(\lVert x \rVert = 1\): il massimo cercato coincide con l'autovalore più grande, il minimo con il più piccolo. Ritroviamo, per una strada completamente diversa, il teorema di min-max della seconda parte.

Questa proposizione è, in pratica, l'analisi delle componenti principali. Cercare la direzione lungo cui un insieme di dati ha varianza massima significa massimizzare \(x^T \Sigma x\) sulla sfera unitaria, con \(\Sigma\) matrice di covarianza — che è simmetrica e semidefinita positiva, come tutte le matrici di Gram. La risposta è l'autovettore associato all'autovalore più grande, e le componenti successive sono gli autovettori seguenti in ordine decrescente.

Conclusione

Con questo si chiude la serie. Ripercorrendo il cammino: nella prima parte abbiamo costruito la struttura astratta — spazi vettoriali, basi, applicazioni lineari, prodotti scalari — e visto che una matrice è un'applicazione lineare scritta in coordinate. Nella seconda abbiamo studiato le matrici per conto loro, fino a decomporle in autovalori e valori singolari per capire in quali direzioni agiscono. In questa terza abbiamo aggiunto la derivazione, che è ciò che trasforma tutto l'apparato in un metodo di calcolo effettivo.

I tre risultati da portare a casa sono probabilmente \(\nabla_x b^T x = b\), \(\nabla_x x^T A x = 2Ax\) e la loro conseguenza, le equazioni normali \(A^T A x = A^T b\). Da lì passa la regressione lineare, e con qualche complicazione in più anche tutto il resto.