Il determinante: quanto una trasformazione gonfia lo spazio

Il determinante è più vecchio della matrice. Nel 1683 il matematico giapponese Seki Takakazu calcola quelli che oggi chiamiamo determinanti, fino al quinto ordine, per risolvere equazioni; dieci anni dopo, in una lettera a l’Hôpital, Leibniz scrive la condizione perché tre equazioni lineari in due incognite abbiano una soluzione comune, cioè l’annullarsi di un determinante \(3\times 3\). Nessuno dei due ha una parola per la tabella di numeri da cui quel valore si ricava: James Joseph Sylvester conierà il termine «matrice» nel 1850, e la teoria arriverà con Arthur Cayley nel 1858.

Per un secolo e mezzo, quindi, si è calcolato il determinante di oggetti che non avevano nome. Quel numero nasce per rispondere alla domanda che la sezione sui sistemi lineari ha lasciato aperta: una soluzione esiste, ed è una sola? La risposta geometrica è ancora più netta, e più vecchia della parola «matrice»: già nel 1773 Lagrange trova che il volume del tetraedro che ha per vertici l’origine e tre punti è un sesto di un determinante \(3\times 3\) delle loro coordinate. In generale quel numero misura di quanto la trasformazione moltiplica i volumi (le aree, nel piano), e se vale zero li annulla: lo spazio è stato schiacciato.

Un numero che misura un’area

Fig. 3.8 Il quadretto di partenza, lo stesso quadretto dopo la trasformazione, e il caso in cui la trasformazione lo appiattisce. Il valore assoluto del determinante è il rapporto fra l’area finale e quella iniziale, e il segno dice se la figura è stata anche specchiata; il determinante vale zero esattamente quando la figura collassa.

Elementare

Prendi un foglio di gomma con sopra disegnato un quadretto di un centimetro di lato, e tiralo. Il quadretto diventa un parallelogramma, e la domanda naturale è di quante volte la sua area sia cambiata. Quel numero, per lo stiramento che hai fatto, è il determinante, e Fig. 3.8 lo disegna.

Quel numero non dipende dal quadretto che hai scelto. Qualunque figura tu avessi disegnato sul foglio (un cerchio, una lettera, la sagoma di un gatto), la sua area sarebbe cambiata dello stesso fattore, perché lo stiramento è lo stesso dappertutto. Un solo numero riassume l’effetto della trasformazione su tutte le aree del piano.

Il determinante ha quattro proprietà, tutte conseguenze del vederlo come un fattore di area.

Può essere negativo. Se oltre a tirare il foglio lo rovesci, l’area c’è ancora ma il disegno è specchiato: ciò che prima girava in senso orario ora gira in senso antiorario. Il segno meno registra il ribaltamento, e il valore assoluto resta il fattore di area.

Se vale zero, hai schiacciato tutto. Uno stiramento che appiattisce il foglio su una retta manda ogni figura in un segmento, e l’area diventa zero. Poiché da un segmento non si risale al quadretto di partenza, la trasformazione non si può disfare: è il caso che la sezione sui sistemi lineari chiamava avere un nucleo più grande del solo zero.

Due stiramenti in sequenza moltiplicano i fattori. Se il primo, ad esempio, raddoppia le aree e il secondo le triplica, insieme le moltiplicano per sei. Ne segue che tornare indietro divide: se una trasformazione moltiplica le aree per cinque, la sua inversa le divide per cinque.

In tre dimensioni il fattore riguarda il volume. Il cubetto di un centimetro di lato diventa un solido obliquo, e il determinante dice di quante volte il volume è cambiato. Oltre le tre dimensioni non c’è più nulla da disegnare, ma il calcolo vale e ha lo stesso significato.

Superiore

Per \(\mathbf{A}\in\mathbb{R}^{2\times 2}\) il determinante è

\[\begin{split} \det\mathbf{A} = \det\begin{pmatrix} a & b\\ c & d\end{pmatrix} = ad - bc , \end{split}\]

ed è l’area con segno del parallelogramma generato dalle due colonne. Più in generale, in \(\mathbb{R}^n\) è il volume con segno del parallelepipedo generato dalle \(n\) colonne, e il fattore vale per ogni figura: se \(S\subset\mathbb{R}^n\) ha volume, \(\operatorname{vol}(\mathbf{A}S)=|\det\mathbf{A}|\,\operatorname{vol}(S)\) (si approssima \(S\) con cubetti, e ciascuno va in un parallelepipedo di volume \(|\det\mathbf{A}|\) volte il suo).

La definizione che rende tutto immediato è quella assiomatica: il determinante è l’unica funzione \(\mathbb{R}^{n\times n}\to\mathbb{R}\) che sia multilineare nelle colonne, alternante (cambia segno scambiando due colonne) e normalizzata (\(\det\mathbf{I}=1\)).

Da questi tre assiomi discende tutto il resto. Per \(n=2\), ad esempio, scrivendo le colonne come \(a\mathbf{e}_1+c\mathbf{e}_2\) e \(b\mathbf{e}_1+d\mathbf{e}_2\), la multilinearità dà quattro termini; due contengono una colonna ripetuta e si annullano per l’alternanza, e restano \(ad\det(\mathbf{e}_1,\mathbf{e}_2)+cb\det(\mathbf{e}_2,\mathbf{e}_1)=ad-bc\). Valgono inoltre le proprietà che si usano davvero:

\[ \det(\mathbf{A}\mathbf{B}) = \det\mathbf{A}\,\det\mathbf{B}, \qquad \det(\mathbf{A}^{-1}) = \frac{1}{\det\mathbf{A}}, \qquad \det(\mathbf{A}^\top) = \det\mathbf{A}, \]

e per uno scalare \(\alpha\), in dimensione \(n\), \(\det(\alpha\mathbf{A}) = \alpha^n\det\mathbf{A}\), che è una delle cose che si sbagliano più spesso scrivendo a memoria (il fattore va all’\(n\)-esima potenza perché scala tutte e \(n\) le direzioni).

La caratterizzazione che lega il determinante al nucleo è

\[ \det\mathbf{A} \neq 0 \;\Longleftrightarrow\; \mathbf{A} \text{ invertibile} \;\Longleftrightarrow\; \ker(\mathbf{A}) = \{\mathbf{0}\} \;\Longleftrightarrow\; \operatorname{rank}(\mathbf{A}) = n , \]

quattro modi di dire la stessa cosa.

Vale infine il legame con lo spettro: \(\det\mathbf{A} = \prod_{i=1}^{n}\lambda_i\), il prodotto degli autovalori contati con la loro molteplicità, e \(|\det\mathbf{A}| = \prod_i \sigma_i\), il prodotto dei valori singolari. La seconda scrittura rende evidente il senso geometrico. La decomposizione ai valori singolari, infatti, dice che ogni trasformazione è una rotazione (se serve con una riflessione), poi una dilatazione lungo assi ortogonali di fattori \(\sigma_i\), poi un’altra rotazione (anche qui, se serve, con una riflessione); queste non toccano i volumi, perché il loro determinante vale \(\pm 1\), quindi il fattore di volume è il prodotto delle dilatazioni.

Come diagnosi, \(|\det|\) non dice nulla sul condizionamento. Una matrice \(\operatorname{diag}(10^{-3},10^{3})\) ha determinante \(1\) e conserva le aree, ma è pessimamente condizionata, perché schiaccia mille volte in una direzione e stira mille volte nell’altra. Il determinante è un prodotto, e in un prodotto un fattore quasi nullo e uno molto grande si compensano.

Poiché si annulla esattamente quando la trasformazione perde delle direzioni, il determinante è la risposta più compatta alla domanda dei sistemi lineari. Vale però soltanto per le matrici quadrate: un sistema con più equazioni che incognite non ha un determinante da guardare, e lì la domanda si sposta sul rango.

import numpy as np

A = np.array([[3.0, 1.0],
              [1.0, 2.0]])
print(np.linalg.det(A))              # 3*2 - 1*1 -> 5.000000000000001

specchio = np.array([[0.0, 1.0],     # scambia le due coordinate
                     [1.0, 0.0]])
print(np.linalg.det(specchio))       # ribalta senza deformare -> -1.0

# due trasformazioni di fila: i fattori si moltiplicano
print(round(np.linalg.det(A @ specchio), 6),
      round(np.linalg.det(A) * np.linalg.det(specchio), 6))   # -> -5.0 -5.0

# seconda colonna = doppio della prima: il foglio viene schiacciato
S = np.array([[2.0, 4.0],
              [1.0, 2.0]])
print(np.linalg.det(S), np.linalg.matrix_rank(S))             # -> 0.0 1

Come si calcola, e come non si calcola

Lo sviluppo di Laplace, che per una matrice \(2\times 2\) o \(3\times 3\) si fa a mano, vale in ogni dimensione, ma costa un numero di operazioni che cresce come \(n!\), e per questo è la strada sbagliata.

Elementare

Per una tabella due per due il calcolo si fa a mano. Con le righe \((3, 1)\) e \((1, 2)\) si moltiplicano i due numeri sulla diagonale che scende, \(3\cdot 2 = 6\), e si toglie il prodotto dei due sull’altra diagonale, \(1\cdot 1 = 1\). Il determinante è quindi \(5\), cioè l’area del parallelogramma di Fig. 3.8.

Per tabelle più grandi c’è un procedimento che si insegna a mano, lo sviluppo, che spezza la tabella in tante tabelle più piccole e ripete il calcolo su ciascuna. Funziona benissimo per due o tre righe, ma diventa assurdo appena dopo: per una tabella di venti righe per venti chiederebbe più di quattro miliardi di miliardi di moltiplicazioni, e una macchina che ne fa un miliardo al secondo ci metterebbe più di un secolo.

La strada buona è la stessa dei sistemi lineari, l’eliminazione: si porta la tabella a scaletta, e a quel punto il determinante è il prodotto dei numeri sulla diagonale, con un cambio di segno per ogni scambio di righe fatto per strada. In questo modo servono circa ottomila operazioni (venti al cubo) invece di quattro miliardi di miliardi, e il calcolo finisce prima che tu abbia alzato gli occhi.

Resta un problema pratico: la grandezza del risultato. Moltiplicare fra loro qualche centinaio di numeri dà un valore assurdamente grande o assurdamente piccolo. Su una tabella di quattrocento righe presa a caso, ad esempio, il determinante esce dalla scala dei numeri che la macchina sa scrivere, e al posto della risposta compare la parola «infinito».

La cura è antica e si usa dappertutto: invece del prodotto si sommano i logaritmi. Il logaritmo conta quante volte si deve moltiplicare un numero fisso per ottenerne un altro (quanti raddoppi portano da uno a otto? tre), quindi trasforma i prodotti in somme; il logaritmo di un numero enorme è un numero comodo, e si somma senza problemi. Le librerie offrono la funzione apposita, che restituisce due cose separate, il segno e il logaritmo del valore assoluto.

Superiore

Lo sviluppo di Laplace lungo una riga o una colonna è

\[ \det\mathbf{A} = \sum_{j=1}^{n} (-1)^{i+j} A_{ij}\,\det \mathbf{A}_{(ij)}, \]

dove \(i\) è una riga qualsiasi, tenuta fissa, e \(\mathbf{A}_{(ij)}\) è la sottomatrice ottenuta cancellando la riga \(i\) e la colonna \(j\) (il suo determinante è il minore di posto \((i,j)\)). Ha costo \(\Theta(n!)\) e serve solo a dimostrare teoremi.

La via praticabile è la fattorizzazione \(\mathbf{P}\mathbf{A}=\mathbf{L}\mathbf{U}\) prodotta dall’eliminazione con pivoting: poiché \(\det\mathbf{L}=1\) e il determinante di una matrice triangolare è il prodotto della diagonale,

\[ \det\mathbf{A} = (-1)^{s}\prod_{i=1}^{n} U_{ii}, \]

con \(s\) il numero di scambi di riga. Il costo scende a \(\Theta(n^3)\), lo stesso dell’eliminazione, e il determinante è un sottoprodotto della risoluzione di un sistema.

Il prodotto di \(n\) pivot ha però una dinamica esponenziale: per \(\mathbf{M}\in\mathbb{R}^{400\times 400}\) con entrate normali standard il valore vero è dell’ordine di \(e^{998}\), cioè fuori dai numeri in virgola mobile a doppia precisione, che si fermano attorno a \(1{,}8\cdot 10^{308}\). Per questo esiste numpy.linalg.slogdet, che restituisce la coppia \((\operatorname{sign}, \log|\det|)\) sommando i logaritmi dei pivot invece di moltiplicarli.

Si noti che nelle applicazioni la quantità che serve è quasi sempre \(\log|\det|\) e non \(\det\), quindi non si perde niente: è il termine che compare nelle log-verosimiglianze, e ottenerlo direttamente evita di calcolare un numero enorme per poi prenderne il logaritmo.

import numpy as np
rng = np.random.default_rng(0)
M = rng.normal(size=(400, 400))

print(np.linalg.det(M))          # -> -inf, con un avviso di overflow

segno, log_det = np.linalg.slogdet(M)
print(segno, round(log_det, 4))  # -> -1.0 998.4819

# su una triangolare il conto è il prodotto della diagonale, e si vede
T = np.tril(rng.normal(size=(5, 5)))
np.fill_diagonal(T, [1.2, 0.7, 2.0, 0.5, 1.5])
print(round(np.linalg.det(T), 8), round(np.prod(np.diag(T)), 8))   # -> 1.26 1.26

Il determinante che va in inf spiega perché nel codice si incontra quasi sempre sotto forma di logaritmo. Il calcolo sulla triangolare, invece, è il principio su cui si fonda un’intera famiglia di modelli generativi.

Il cambio di variabile: come cambia una densità quando lo spazio si deforma

Fin qui il determinante ha misurato aree e volumi. La ragione per cui torna di continuo nel machine learning è un’altra, e discende da questa.

Elementare

Hai un litro di vernice steso su una parete di dieci metri quadri, quindi un decimo di litro per metro quadro. Se la parete si allarga fino a venti metri quadri e la vernice resta quella, lo strato si dimezza: la quantità totale non cambia, quindi quanta ce n’è per metro quadro deve scendere nella stessa misura in cui la parete si è allargata.

In generale, ogni volta che una trasformazione deforma lo spazio, qualunque cosa sia distribuita su quello spazio si concentra dove lo spazio si è ristretto e si dirada dove si è allargato, e il fattore è il determinante.

Dove la trasformazione non è la stessa dappertutto (un foglio tirato più da una parte che dall’altra) si guarda punto per punto: attorno a ogni punto, per uno spostamento piccolo, ogni deformazione somiglia a uno stiramento uniforme, e il determinante di quello stiramento locale dice quanto la vernice si assottiglia lì.

Questo calcolo permette a certi generatori di immagini di dire quanto è probabile ciò che hanno prodotto, invece di limitarsi a produrlo. Partono da una nuvola di punti semplicissima e la deformano, con una trasformazione che si può disfare, finché non somiglia ai dati veri. Per sapere quanto è denso il risultato in un punto bisogna sapere quanto la deformazione ha allargato lo spazio in quel punto, e questo richiede il determinante in ogni punto e a ogni passo.

Il calcolo, però, ha un costo: richiede il cubo del numero di coordinate, e su un’immagine di centomila pixel sarebbe fuori discussione. Allora si costruisce la trasformazione apposta per renderlo semplice. In particolare, si impone che la prima coordinata nuova dipenda solo dalla prima vecchia, la seconda dalle prime due, la terza dalle prime tre.

Con questo vincolo la tabella viene a scaletta, e il determinante torna a essere il prodotto della sua diagonale: un fattore per coordinata, quindi un costo proporzionale al numero di coordinate invece che al suo cubo. Una proprietà matematica decide così la forma di un’architettura: la trasformazione è scelta in modo che il calcolo sia semplice.

Superiore

Sia \(\mathbf{f}:\mathbb{R}^n\to\mathbb{R}^n\) un diffeomorfismo (cioè invertibile, differenziabile con continuità e con inversa differenziabile) e sia \(\mathbf{y}=\mathbf{f}(\mathbf{x})\). In modo equivalente, per una \(\mathbf{f}\) invertibile di classe \(C^1\), si chiede che la matrice jacobiana \(\mathbf{J}_{\mathbf{f}}\) (definita qui sotto) abbia \(\det\mathbf{J}_{\mathbf{f}}\neq 0\) in ogni punto. L’ipotesi sull’inversa serve. La funzione \(f(x)=x^3\), ad esempio, è invertibile e derivabile, ma in \(0\) ha derivata nulla, e la densità trasformata vi diverge.

La formula del cambio di variabile per una densità è

\[ p_Y(\mathbf{y}) = p_X(\mathbf{x})\, \left|\det \mathbf{J}_{\mathbf{f}}(\mathbf{x})\right|^{-1}, \qquad \big(\mathbf{J}_{\mathbf{f}}\big)_{ij} = \frac{\partial f_i}{\partial x_j}, \]

dove \(\mathbf{J}_{\mathbf{f}}\) è la matrice jacobiana, cioè la migliore approssimazione lineare di \(\mathbf{f}\) attorno a \(\mathbf{x}\), e le derivate parziali che la compongono arrivano nella sezione su analisi e ottimizzazione.

La formula esprime la conservazione della massa: la probabilità contenuta in un volumetto attorno a \(\mathbf{x}\) è la stessa che sta nel volumetto immagine attorno a \(\mathbf{y}\), che è \(\lvert\det\mathbf{J}_{\mathbf{f}}(\mathbf{x})\rvert\) volte più grande, quindi la densità si divide per quel fattore.

In forma logaritmica, che è quella che si usa,

\[ \log p_Y(\mathbf{y}) = \log p_X(\mathbf{x}) - \log\left|\det \mathbf{J}_{\mathbf{f}}(\mathbf{x})\right| , \]

e componendo \(K\) trasformazioni i termini si sommano. Per la regola della catena, infatti, la jacobiana della composizione è il prodotto delle jacobiane calcolate nei punti intermedi, e il determinante di un prodotto è il prodotto dei determinanti. È l’ossatura dei flussi normalizzanti, i modelli generativi fatti di trasformazioni invertibili messe in fila, che la sezione sul flusso che si può invertire costruisce proprio su questa formula.

Il vincolo di progetto discende dal costo. Un determinante generico costa \(\Theta(n^3)\) per passo, insostenibile per \(n\) dell’ordine di \(10^5\). Le architetture si costruiscono quindi perché la jacobiana sia triangolare: imponendo che \(f_i\) dipenda solo da \(x_1,\dots,x_i\), tutti gli elementi sopra la diagonale sono nulli e

\[ \log\left|\det \mathbf{J}\right| = \sum_{i=1}^{n} \log\left|\frac{\partial f_i}{\partial x_i}\right| , \]

cioè \(\Theta(n)\).

Il determinante compare in altri due posti. Nella densità della gaussiana multivariata, il fattore di normalizzazione contiene \((\det\boldsymbol{\Sigma})^{-1/2}\) e misura il volume dell’ellissoide di covarianza; nella verosimiglianza dei processi gaussiani, \(\log\det\boldsymbol{\Sigma}\) e \(\mathbf{y}^\top\boldsymbol{\Sigma}^{-1}\mathbf{y}\) si ricavano dalla stessa fattorizzazione di Cholesky della matrice di covarianza, che costa \(\Theta(N^3)\) nel numero \(N\) di osservazioni [RW06].

Una rete neurale ordinaria, invece, in genere non è invertibile: la funzione di attivazione, cioè quella che piega i numeri fra uno strato e l’altro, nella forma più diffusa (la ReLU) manda a zero tutti i valori negativi, e da uno zero non si risale al numero di partenza. Il determinante che misura quanto la rete deforma lo spazio attorno a ogni punto, dilatandolo o comprimendolo, vale zero su intere regioni, e quindi la formula del cambio di variabile non si applica. Quando poi la rete manda i dati in uno spazio di dimensione diversa, quel determinante non esiste affatto.

I flussi calcolano la probabilità esatta con la formula del cambio di variabile, e per riuscirci accettano una trasformazione di forma obbligata. Altre famiglie arrivano alla probabilità per altre vie (i modelli autoregressivi la scompongono in un prodotto di probabilità condizionate, come racconta la sezione sul generare un pixel alla volta) o ne danno soltanto un limite inferiore, e altre ancora, come le GAN, non la calcolano affatto.

In pratica, con NumPy

import numpy as np

A = np.array([[3.0, 1.0], [1.0, 2.0]])

np.linalg.det(A)          # il fattore di area, con segno
np.linalg.slogdet(A)      # (segno, log|det|): la forma da usare quasi sempre
np.linalg.matrix_rank(A)  # la domanda giusta se la matrice non è quadrata

np.linalg.det(2 * A)      # -> 4 volte det(A): il fattore va alla potenza n

Raddoppiare una matrice \(2\times 2\) ne quadruplica il determinante, e su una \(n\times n\) lo moltiplica per \(2^n\), perché si raddoppiano \(n\) lati del parallelepipedo.

Elementare

Da ricordare

• Il determinante è di quante volte una trasformazione cambia le aree (o i volumi), ed è lo stesso numero per qualunque figura, perché lo stiramento è lo stesso dappertutto.

• Il segno dice se la trasformazione ha anche ribaltato il foglio. Il valore zero dice che ha schiacciato tutto su una retta: da lì non si torna indietro, ed è lo stesso caso in cui un sistema non ha una sola soluzione.

• Due trasformazioni di fila moltiplicano i loro fattori, e disfarne una divide per il suo.

• Lo sviluppo che si fa a mano, su venti righe, chiederebbe più di un secolo di conti: si calcola con l’eliminazione, portando la tabella a scaletta e moltiplicando la diagonale. Poiché quel prodotto diventa subito enorme, nel codice si usa il logaritmo del determinante invece del determinante.

• Ciò che è distribuito nello spazio si assottiglia dove lo spazio si allarga, nella misura data dal determinante. Questo calcolo permette a certi generatori di dire quanto è probabile ciò che producono, e li obbliga a una forma che lo rende poco oneroso.

Superiore

Da ricordare

• \(\det\mathbf{A}\) è il volume con segno del parallelepipedo generato dalle colonne; è l’unica funzione multilineare, alternante e con \(\det\mathbf{I}=1\). Valgono \(\det(\mathbf{A}\mathbf{B})=\det\mathbf{A}\det\mathbf{B}\), \(\det(\mathbf{A}^{-1})=1/\det\mathbf{A}\) e \(\det(\alpha\mathbf{A})=\alpha^n\det\mathbf{A}\).

• \(\det\mathbf{A}\neq 0\) equivale a invertibilità, a nucleo banale e a rango pieno. Inoltre \(\det\mathbf{A}=\prod_i\lambda_i\) e \(|\det\mathbf{A}|=\prod_i\sigma_i\). Il determinante non dice niente sul condizionamento: \(\operatorname{diag}(10^{-3},10^{3})\) ha determinante \(1\).

• Laplace costa \(\Theta(n!)\) e serve per le dimostrazioni; si calcola con \(\mathbf{P}\mathbf{A}=\mathbf{L}\mathbf{U}\) in \(\Theta(n^3)\), e nel codice si usa slogdet perché il prodotto dei pivot va in overflow già a poche centinaia di righe.

• Cambio di variabile: \(\log p_Y(\mathbf{y}) = \log p_X(\mathbf{x}) - \log|\det\mathbf{J}_{\mathbf{f}}(\mathbf{x})|\), e componendo trasformazioni i termini si sommano. È l’ossatura dei flussi normalizzanti, i modelli generativi fatti di trasformazioni invertibili in fila, e la ragione per cui le loro jacobiane si costruiscono triangolari, portando il costo da \(\Theta(n^3)\) a \(\Theta(n)\).

• Lo stesso \(\log\det\) compare nella normalizzazione della gaussiana multivariata e nella verosimiglianza dei processi gaussiani, dove la fattorizzazione da cui si ricava rende il metodo cubico nel numero di osservazioni.

Con il determinante i vettori, le matrici e le trasformazioni che le matrici rappresentano sono al completo: si sa mettere i dati in fila, trasformarli, chiedersi che cosa una trasformazione può produrre e che cosa perde, trovare la risposta migliore quando quella esatta non esiste, e misurare come cambiano i volumi. Serve ora un altro linguaggio, che non guarda una trasformazione fissa ma il modo in cui una quantità cambia quando se ne muove un’altra: le derivate, da cui si ricava come migliorare un modello.