23 - Algebra lineare 1: vettori, matrici, sistemi

Corso Machine Learning Engineer · Modulo 02 · Matematica e statistica · Lezione 8

Obiettivo. Vedere l'algebra lineare per quello che è: geometria dello spazio, non tabelle di numeri. Alla fine saprai leggere un vettore come una freccia e come una riga di dati, calcolare e interpretare il prodotto scalare e la similarità del coseno, capire una matrice come una trasformazione dello spazio (con il determinante come fattore d'area), risolvere un sistema lineare con l'eliminazione di Gauss scritta da zero, riconoscere un problema mal condizionato, e derivare le equazioni normali dei minimi quadrati con molte variabili.

Falso mito da sfatare. «Per risolvere Ax = b si calcola l'inversa: x = A⁻¹b.» È il metodo più lento e meno accurato. Si risolve il sistema direttamente (lo faremo noi, a mano, e poi con np.linalg.solve).

Introduzione

Un dataset è una matrice: righe per esempi, colonne per caratteristiche. Un modello lineare è un prodotto X @ w. Una rete neurale è una successione di moltiplicazioni per matrici e di non linearità. Per capire che cosa fanno queste operazioni (e quando falliscono) serve la geometria che le sta sotto. Questa lezione la costruisce; la prossima aggiunge gli autovalori, la SVD e la PCA.

Il percorso parte dai vettori e dal prodotto scalare, passa alle matrici come trasformazioni dello spazio e al determinante, risolve i sistemi lineari con l'eliminazione di Gauss e affronta rango, inversa e condizionamento, per arrivare alle equazioni normali dei minimi quadrati. Il prodotto di matrici riprende la composizione di funzioni degli insiemi e la regola della catena di Analisi 1; la vettorizzazione in NumPy è il suo ambiente di lavoro.

Vettori

L'analogia: lo spostamento

Un vettore è uno spostamento: «tre passi a est e due a nord» è il vettore (3, 2), indipendentemente da dove parti. Si disegna come una freccia dall'origine al punto (3, 2). Con d componenti è uno spostamento in uno spazio a d dimensioni.

Lo stesso oggetto è anche una riga di dati: un cliente con età 40, reddito 35 000, 3 acquisti è il vettore (40, 35000, 3) nello spazio delle caratteristiche (lezione 1: il prodotto cartesiano ℝ³). Qui si fonde la geometria con il machine learning: un esempio è un punto, un modello è qualcosa che agisce su punti.

Due operazioni:

• Somma: si sommano le componenti, (3, 2) + (1, 4) = (4, 6): due spostamenti in sequenza.
• Moltiplicazione per uno scalare c: allunga (o accorcia, o inverte se c < 0) la freccia.

Una combinazione lineare c₁v₁ + c₂v₂ + … è la ricetta per costruire nuovi vettori con queste due operazioni. L'insieme di tutte le combinazioni lineari dei vettori v₁, …, v_k è il loro span (il sottospazio che generano): per un vettore è una retta, per due vettori non paralleli un piano.

import numpy as np

v = np.array([3.0, 2.0])
w = np.array([1.0, 4.0])
print("somma           :", v + w)
print("2.5 * v         :", 2.5 * v)
print("combinazione 2v - w:", 2 * v - w)

# un cliente come vettore di caratteristiche
cliente = np.array([40, 35_000, 3])
print("forma:", cliente.shape, "| spazio R^", cliente.size)

Output ottenuto:

somma           : [4. 6.]
2.5 * v         : [7.5 5. ]
combinazione 2v - w: [5. 0.]
forma: (3,) | spazio R^ 3

Cosa fa il codice. Esegue le due operazioni elementari con gli array NumPy (che le fanno componente per componente, come nella lezione 5) e mostra che un cliente è un vettore in ℝ³. L'espressione 2 * v - w è una combinazione lineare con coefficienti 2 e −1.

Il prodotto scalare, la norma e la similarità del coseno

Definizione e significato

Il prodotto scalare (dot product) di a e b è la somma dei prodotti delle componenti:

a · b = Σ aᵢ·bᵢ.

Ha un'interpretazione geometrica: a · b = |a|·|b|·cos θ, con |a| = √(a · a) la lunghezza (norma) di a e θ l'angolo tra i due vettori. Proviamo la corrispondenza nel piano: scrivi a = |a|(cos α, sin α) e b = |b|(cos β, sin β). Allora a · b = |a||b|(cos α cos β + sin α sin β) = |a||b| cos(α − β), per la formula di sottrazione del coseno, e α − β è l'angolo tra i vettori. ∎

Conseguenze:

• Perpendicolari ⇔ a · b = 0 (il coseno di 90° è 0).
• Proiezione: la «ombra» di b sulla direzione di a ha lunghezza (a · b)/|a|.
• Disuguaglianza di Cauchy-Schwarz: poiché |cos θ| ≤ 1, vale |a · b| ≤ |a|·|b|.
• Similarità del coseno: cos θ = (a · b)/(|a||b|) misura quanto due vettori puntano nella stessa direzione, ignorando la loro lunghezza: valore 1 stessa direzione, 0 perpendicolari, −1 opposti.

Le norme più usate: ‖x‖₂ = √Σ xᵢ² (euclidea), ‖x‖₁ = Σ |xᵢ| (somma dei valori assoluti, «distanza dei taxi»), ‖x‖∞ = max |xᵢ|. La distanza tra due punti è la norma della differenza.

a = np.array([1.0, 2.0, 3.0])
b = np.array([4.0, 5.0, 6.0])
print("a·b a mano   :", sum(x * y for x, y in zip(a, b)), "| con @:", a @ b, "| np.dot:", np.dot(a, b))
print("norme di a   : L2 =", round(float(np.linalg.norm(a)), 4), "| L1 =", np.linalg.norm(a, 1), "| Linf =", np.linalg.norm(a, np.inf))

cos = (a @ b) / (np.linalg.norm(a) * np.linalg.norm(b))
print("cos θ =", round(float(cos), 4), "| angolo in gradi =", round(float(np.degrees(np.arccos(cos))), 2))
print("perpendicolari (1,2)·(-2,1) =", np.array([1, 2]) @ np.array([-2, 1]))

# Cauchy-Schwarz su 1000 coppie casuali
rng = np.random.default_rng(0)
A, B = rng.normal(size=(1000, 5)), rng.normal(size=(1000, 5))
lhs = np.abs((A * B).sum(axis=1))
rhs = np.linalg.norm(A, axis=1) * np.linalg.norm(B, axis=1)
print("Cauchy-Schwarz vale per tutte le 1000 coppie:", bool(np.all(lhs <= rhs)))

Output ottenuto:

a·b a mano   : 32.0 | con @: 32.0 | np.dot: 32.0
norme di a   : L2 = 3.7417 | L1 = 6.0 | Linf = 3.0
cos θ = 0.9746 | angolo in gradi = 12.93
perpendicolari (1,2)·(-2,1) = 0
Cauchy-Schwarz vale per tutte le 1000 coppie: True

Cosa fa il codice. Calcola il prodotto scalare in tre modi equivalenti (somma esplicita, operatore @, np.dot), le tre norme del vettore a, e l'angolo tra a e b invertendo il coseno. Poi verifica che (1, 2) e (−2, 1) sono perpendicolari (prodotto 0) e controlla la disuguaglianza di Cauchy-Schwarz su 1000 coppie casuali: (A * B).sum(axis=1) è il prodotto scalare riga per riga, e np.linalg.norm(..., axis=1) calcola le norme di tutte le righe in un colpo.

Un'applicazione: documenti simili

Rappresentiamo tre piccoli testi come vettori che contano le parole, e confrontiamoli con la similarità del coseno. Il vantaggio rispetto alla distanza: un testo lungo il doppio, che ripete le stesse parole, deve risultare identico in direzione.

vocabolario = ["gatto", "cane", "pesce", "mangia", "dorme"]
doc1 = np.array([3, 0, 0, 2, 1], dtype=float)         # parla di un gatto
doc2 = doc1 * 2                                        # lo stesso testo, ripetuto due volte
doc3 = np.array([0, 3, 1, 1, 2], dtype=float)         # parla di un cane e di un pesce

def coseno(x, y):
    return float(x @ y / (np.linalg.norm(x) * np.linalg.norm(y)))

print("coseno(doc1, doc2) =", round(coseno(doc1, doc2), 4), "| distanza euclidea =", round(float(np.linalg.norm(doc1 - doc2)), 3))
print("coseno(doc1, doc3) =", round(coseno(doc1, doc3), 4), "| distanza euclidea =", round(float(np.linalg.norm(doc1 - doc3)), 3))

Output ottenuto:

coseno(doc1, doc2) = 1.0 | distanza euclidea = 3.742
coseno(doc1, doc3) = 0.276 | distanza euclidea = 4.583

Cosa fa il codice. doc2 è doc1 con tutte le frequenze raddoppiate: ha lunghezza diversa ma stessa direzione. La similarità del coseno lo riconosce (1.0), la distanza euclidea li considera lontani (ha lunghezza |doc1|). Questo è il motivo per cui nel recupero di documenti e negli embedding (vettori che rappresentano parole o frasi) si usa il coseno.

Trade-off. Il coseno ignora l'ampiezza: due clienti, uno che spende 10 e uno che spende 1000 negli stessi proporzioni, risultano identici. Se l'ampiezza conta, serve la distanza (o una normalizzazione diversa).

Matrici

L'analogia: la macchina che trasforma lo spazio

Una matrice A di dimensioni m × n (m righe, n colonne) è una macchina che prende un vettore x con n componenti e ne restituisce uno con m componenti: y = A x. In particolare, se A è quadrata (n × n) trasforma lo spazio in sé: ruota, allunga, schiaccia, inclina.

La regola per calcolarla: ogni componente di y è il prodotto scalare di una riga di A con x:

yᵢ = Σⱼ Aᵢⱼ·xⱼ.

Ma c'è un secondo modo di vederla, molto più illuminante, il punto di vista delle colonne:

A x = x₁·(colonna 1 di A) + x₂·(colonna 2 di A) + …

cioè A x è una combinazione lineare delle colonne di A, con i coefficienti dati da x. In particolare le colonne di A sono le immagini dei vettori di base: se e₁ = (1, 0) e e₂ = (0, 1), allora A e₁ è la prima colonna e A e₂ la seconda. Conoscendo dove vanno i vettori di base, conosci l'intera trasformazione.

Il determinante

Il quadrato unitario (lati e₁ ed e₂, area 1) diventa il parallelogramma di lati A e₁ e A e₂. La sua area è il determinante (in valore assoluto). Per A = [[a, b], [c, d]] le colonne sono (a, c) e (b, d), e l'area del parallelogramma che generano è |ad − bc|: quindi det A = ad − bc. (Il segno dice se la trasformazione ribalta l'orientamento.) Nella figura: 2·1.5 − 1·0.5 = 2.5.

Se det A = 0 il parallelogramma si schiaccia su una retta (area zero): la trasformazione perde un'intera dimensione, e non è invertibile (non puoi ricostruire l'ingresso).

A = np.array([[2.0, 1.0],
              [0.5, 1.5]])
x = np.array([3.0, -1.0])

print("A @ x                        :", A @ x)
print("combinazione delle colonne   :", 3 * A[:, 0] + (-1) * A[:, 1])
print("det A (formula ad - bc)      :", A[0, 0] * A[1, 1] - A[0, 1] * A[1, 0], "| np.linalg.det:", round(float(np.linalg.det(A)), 6))

quadrato = np.array([[0, 1, 1, 0], [0, 0, 1, 1]], dtype=float)       # i 4 vertici, uno per colonna
print("vertici del quadrato dopo A:\n", A @ quadrato)

B = np.array([[1.0, 2.0], [2.0, 4.0]])                                # seconda colonna = 2 × prima
print("det di una matrice con colonne dipendenti:", np.linalg.det(B))

Output ottenuto:

A @ x                        : [5. 0.]
combinazione delle colonne   : [5. 0.]
det A (formula ad - bc)      : 2.5 | np.linalg.det: 2.5
vertici del quadrato dopo A:
 [[0.  2.  3.  1. ]
 [0.  0.5 2.  1.5]]
det di una matrice con colonne dipendenti: 0.0

Cosa fa il codice. A @ x e la combinazione delle colonne (3·colonna1 − 1·colonna2) danno lo stesso vettore: sono due letture della stessa operazione. Il determinante calcolato con la formula coincide con quello di NumPy. Moltiplicando A per la matrice dei quattro vertici del quadrato (uno per colonna) si trasformano tutti insieme: si ottengono i vertici del parallelogramma. L'ultima matrice ha la seconda colonna pari al doppio della prima: le colonne sono dipendenti, lo span è una sola retta e il determinante è zero.

Il prodotto di matrici: comporre trasformazioni

Applicare prima B e poi A a un vettore: A(Bx). È ancora una trasformazione lineare, e la sua matrice è il prodotto AB. Per questo (AB)x = A(Bx): il prodotto di matrici è la composizione di funzioni (lezione 1, dove abbiamo introdotto g ∘ f; e la regola della catena della lezione 6). Gli elementi sono (AB)ᵢⱼ = Σₖ Aᵢₖ Bₖⱼ, e le dimensioni devono combaciare: (m × n)(n × p) = (m × p).

Il prodotto non è commutativo: ruotare e poi allungare non è come allungare e poi ruotare.

theta = np.radians(90)
R = np.array([[np.cos(theta), -np.sin(theta)],
              [np.sin(theta),  np.cos(theta)]])          # rotazione di 90° antioraria
S = np.array([[2.0, 0.0],
              [0.0, 1.0]])                                # allunga di 2 la direzione x
v = np.array([1.0, 1.0])

print("prima ruoto, poi allungo  S(Rv):", np.round(S @ (R @ v), 4))
print("prima allungo, poi ruoto  R(Sv):", np.round(R @ (S @ v), 4))
print("S@R == R@S ?", np.allclose(S @ R, R @ S))
print("(S R) v = S (R v):", np.allclose((S @ R) @ v, S @ (R @ v)))

M, N = np.arange(6).reshape(2, 3), np.arange(12).reshape(3, 4)
print("forme: (2,3)@(3,4) ->", (M @ N).shape)
print("(MN)^T = N^T M^T:", np.array_equal((M @ N).T, N.T @ M.T))

Output ottenuto:

prima ruoto, poi allungo  S(Rv): [-2.  1.]
prima allungo, poi ruoto  R(Sv): [-1.  2.]
S@R == R@S ? False
(S R) v = S (R v): True
forme: (2,3)@(3,4) -> (2, 4)
(MN)^T = N^T M^T: True

Cosa fa il codice. R è la matrice di rotazione di 90° (i vettori di base finiscono su (0, 1) e (−1, 0)), S allunga di 2 lungo x. Applicate in ordine diverso danno vettori diversi: (−2, 1) contro (−1, 2) (a meno di arrotondamenti), e np.allclose(S @ R, R @ S) è falso. La composizione è associativa: (SR)v = S(Rv). Le ultime righe mostrano la regola delle dimensioni e la proprietà della trasposta: (MN)ᵀ = NᵀMᵀ (l'ordine si inverte).

I sistemi lineari

Che cosa sono

Un sistema lineare è un insieme di equazioni in cui le incognite compaiono solo moltiplicate per numeri e sommate. Esempio con tre incognite:

 2x +  y −  z =   8
−3x −  y + 2z = −11
−2x +  y + 2z =  −3

In forma compatta: A x = b, con A la matrice dei coefficienti, x il vettore delle incognite, b i termini noti. Letto con le colonne: «trova la combinazione delle colonne di A che produce b». Geometricamente ogni equazione è un piano nello spazio, e la soluzione è il punto in cui i tre piani si incontrano.

L'eliminazione di Gauss

L'idea (che usava già Gauss e prima ancora i matematici cinesi) è di combinare le equazioni per far sparire le incognite una alla volta, fino ad avere una forma triangolare che si risolve «dal basso»:

1. Si usa la prima equazione per eliminare x dalle altre: si sottrae da ciascuna un multiplo adatto della prima.
2. Si usa la seconda equazione (già senza x) per eliminare y dalle successive.
3. Si arriva a un sistema triangolare superiore (U x = c): l'ultima equazione dà direttamente l'ultima incognita; poi si sostituisce all'indietro (back substitution).

Il numero che si usa per eliminare, nella riga di lavoro, si chiama pivot. Se il pivot è zero o molto piccolo si crea un problema: si scambiano le righe in modo da usare come pivot il valore più grande in modulo della colonna (pivoting parziale), per evitare divisioni per numeri minuscoli che amplificano gli errori di arrotondamento.

def gauss(A, b):
    """Risolve A x = b con eliminazione di Gauss e pivoting parziale. Restituisce x."""
    A = np.array(A, dtype=float)            # copie: non modifichiamo gli input
    b = np.array(b, dtype=float)
    n = len(b)
    for k in range(n):
        # pivoting: porto in posizione k la riga con il valore più grande in modulo nella colonna k
        p = k + int(np.argmax(np.abs(A[k:, k])))
        if abs(A[p, k]) < 1e-12:
            raise ValueError("matrice singolare (o quasi): nessuna soluzione unica")
        if p != k:
            A[[k, p]], b[[k, p]] = A[[p, k]], b[[p, k]]
        # eliminazione: azzero la colonna k sotto il pivot
        for i in range(k + 1, n):
            m = A[i, k] / A[k, k]           # moltiplicatore
            A[i, k:] -= m * A[k, k:]
            b[i] -= m * b[k]
    # sostituzione all'indietro
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (b[i] - A[i, i + 1:] @ x[i + 1:]) / A[i, i]
    return x

A = [[2, 1, -1], [-3, -1, 2], [-2, 1, 2]]
b = [8, -11, -3]
x = gauss(A, b)
print("soluzione con Gauss:", np.round(x, 6))
print("np.linalg.solve    :", np.round(np.linalg.solve(A, b), 6))
print("verifica A x - b   :", np.round(np.array(A) @ x - np.array(b), 12))

try:
    gauss([[1, 2], [2, 4]], [3, 6])
except ValueError as e:
    print("sistema con righe dipendenti:", e)

Output ottenuto:

soluzione con Gauss: [ 2.  3. -1.]
np.linalg.solve    : [ 2.  3. -1.]
verifica A x - b   : [0. 0. 0.]
sistema con righe dipendenti: matrice singolare (o quasi): nessuna soluzione unica

Cosa fa il codice, passo per passo. Il ciclo esterno sceglie, per ogni colonna k, la riga con il pivot più grande (np.argmax sul valore assoluto della parte sotto la diagonale) e la scambia in posizione k se serve. Poi, per ogni riga sotto, calcola il moltiplicatore m = A[i,k]/A[k,k] e sottrae m volte la riga del pivot (da entrambe le matrici A e b): la colonna k sotto il pivot diventa zero. Alla fine la sostituzione all'indietro calcola x dall'ultima incognita alla prima. Il risultato (2, 3, −1) coincide con np.linalg.solve e il residuo A x − b è zero. L'ultimo esempio ha la seconda riga uguale al doppio della prima: nessun pivot utilizzabile, e la funzione lo segnala invece di restituire spazzatura.

Costo. L'eliminazione richiede circa n³/3 moltiplicazioni: raddoppiare il numero di incognite moltiplica il tempo per 8. Per sistemi con migliaia di incognite è ancora fattibile; per milioni serve sfruttare la struttura (matrici sparse, metodi iterativi come la discesa del gradiente della lezione scorsa).

Rango, inversa, e quando un sistema non ha una soluzione unica

Il rango di una matrice è il numero di colonne linearmente indipendenti (colonne che non si possono ottenere come combinazione delle altre); coincide con il numero di righe indipendenti ed è la «dimensione» dello spazio che la matrice può raggiungere. Per una matrice quadrata n × n:

• rango n: invertibile, det ≠ 0, Ax = b ha una sola soluzione per ogni b.
• rango < n: singolare, det = 0, Ax = b ha o nessuna soluzione o infinite.

L'inversa A⁻¹ è la matrice che «disfa» A: A⁻¹A = I, con I la matrice identità (la trasformazione che non fa nulla). Esiste solo se A è invertibile. Ma calcolarla per risolvere un sistema è uno spreco: costa circa tre volte più dell'eliminazione ed è meno accurato. In pratica np.linalg.solve(A, b), mai np.linalg.inv(A) @ b.

C = np.array([[1.0, 2.0, 3.0],
              [4.0, 5.0, 6.0],
              [7.0, 8.0, 9.0]])
print("rango di C:", np.linalg.matrix_rank(C), "| det C =", round(float(np.linalg.det(C)), 10))
print("colonna 3 = 2·colonna 2 − colonna 1:", np.allclose(C[:, 2], 2 * C[:, 1] - C[:, 0]))

rng = np.random.default_rng(1)
G = rng.normal(size=(300, 300))
g = rng.normal(size=300)
x1 = np.linalg.solve(G, g)
x2 = np.linalg.inv(G) @ g
print("residuo con solve:", f"{np.linalg.norm(G @ x1 - g):.2e}", "| residuo con inv:", f"{np.linalg.norm(G @ x2 - g):.2e}")

Output ottenuto:

rango di C: 2 | det C = 0.0
colonna 3 = 2·colonna 2 − colonna 1: True
residuo con solve: 1.40e-12 | residuo con inv: 2.06e-12

Cosa fa il codice. La matrice C (i numeri da 1 a 9) sembra «generica» ma ha rango 2: la terza colonna è una combinazione delle prime due (2·col2 − col1). Il determinante è zero (a meno di arrotondamento) e quindi non è invertibile. Nella seconda parte risolve un sistema casuale 300 × 300 con solve e con l'inversa: i due residui sono confrontabili su questa matrice ben condizionata (qui solve è un po' migliore); la differenza si sente sul tempo di calcolo (l'inversa costa circa il triplo in numero di operazioni, secondo i conteggi standard: tempo non misurato in questo listato) e, soprattutto, su matrici mal condizionate, dove l'inversa esplicita amplifica gli errori.

Il condizionamento

Un sistema può essere teoricamente risolvibile (det ≠ 0) e praticamente irrisolvibile: è il caso di una matrice mal condizionata, cioè quasi singolare. Piccole perturbazioni dei dati (anche solo l'arrotondamento di un float) producono enormi variazioni nella soluzione.

L'analogia: la leva. Spostare un masso con una leva lunga richiede una forza piccola, ma un piccolo errore di mano alla estremità si traduce in un grande spostamento del masso. Una matrice mal condizionata è una leva lunghissima: amplifica gli errori.

Il numero di condizionamento κ(A) è il rapporto tra il massimo e il minimo «stiramento» che A applica ai vettori (lo calcoleremo con la SVD nella prossima lezione). Regola empirica: si perdono circa log₁₀ κ cifre decimali di accuratezza nella soluzione.

# una matrice quasi singolare: due righe quasi uguali
A_buona = np.array([[1.0, 1.0], [1.0, -1.0]])
A_cattiva = np.array([[1.0, 1.0], [1.0, 1.0001]])
b = np.array([2.0, 2.0])

for nome, A in (("ben condizionata", A_buona), ("mal condizionata", A_cattiva)):
    x = np.linalg.solve(A, b)
    x_pert = np.linalg.solve(A, b + np.array([0.0, 0.0001]))        # perturbo b di 1e-4
    print(f"{nome}: κ = {np.linalg.cond(A):.3e} | x = {np.round(x, 4)} | con b perturbato: {np.round(x_pert, 4)}")

# il polinomio di grado 14 della lezione su Matplotlib: la matrice di Vandermonde
for grado in (3, 5, 9, 14):
    punti = np.linspace(0, 1, 15)
    V = np.vander(punti, grado + 1)                                    # colonne: x^grado, ..., x, 1
    print(f"grado {grado:>2}: numero di condizionamento della matrice dei dati = {np.linalg.cond(V):.2e}")

Output ottenuto:

ben condizionata: κ = 1.000e+00 | x = [2. 0.] | con b perturbato: [ 2.e+00 -1.e-04]
mal condizionata: κ = 4.000e+04 | x = [2. 0.] | con b perturbato: [1. 1.]
grado  3: numero di condizionamento della matrice dei dati = 1.05e+02
grado  5: numero di condizionamento della matrice dei dati = 3.16e+03
grado  9: numero di condizionamento della matrice dei dati = 4.63e+06
grado 14: numero di condizionamento della matrice dei dati = 4.03e+11

Cosa fa il codice. Risolve due sistemi quasi identici. Per la matrice «buona» perturbare b di 10⁻⁴ cambia la soluzione di quell'ordine; per quella «cattiva» (due righe quasi uguali) la stessa perturbazione sposta la soluzione di ordini di grandezza di più. La seconda parte ritrova il misterioso errore del polinomio di grado 14 incontrato in Visualizzare i dati con Matplotlib: costruisce la matrice di Vandermonde (le colonne sono le potenze di x) usata per adattare un polinomio ai dati; il suo numero di condizionamento cresce rapidamente col grado, fino a valori enormi per il grado 14. Adattare quel polinomio significa risolvere un sistema con quella matrice, e gli errori di arrotondamento vengono amplificati del fattore κ: ecco perché il polinomio di grado 14 «esplodeva». È il collegamento tra algebra lineare e aritmetica in virgola mobile.

Rimedi. Standardizzare le caratteristiche; rimuovere variabili quasi collineari; usare la regolarizzazione (aggiungere un termine che penalizza i pesi grandi, che migliora il condizionamento); usare metodi numericamente stabili (QR, SVD) invece delle equazioni normali.

I minimi quadrati con molte variabili

Il problema

Con n esempi e d caratteristiche, X è una matrice n × d, w un vettore di d pesi, y il vettore degli n valori da prevedere. Vogliamo w tale che X w ≈ y. Di solito n > d (più equazioni che incognite): il sistema X w = y non ha soluzione esatta, e si cerca la migliore in senso quadratico, minimizzando

L(w) = ‖y − X w‖² = (y − Xw)ᵀ(y − Xw).

Il punto di vista geometrico

Xw è una combinazione lineare delle colonne di X: sta nel piano (sottospazio) generato dalle colonne. y in generale sta fuori da quel piano. Il punto del piano più vicino a y è la sua proiezione ortogonale ŷ, e quello che resta (il residuo r = y − ŷ) è perpendicolare al piano, cioè perpendicolare a ciascuna colonna di X:

Xᵀ (y − X w) = 0 ⇒ XᵀX · w = Xᵀy.

Queste sono le equazioni normali (normali nel senso di «perpendicolari»).

La derivazione con il gradiente

Verifichiamo con gli strumenti della lezione 7. Sviluppando: L(w) = yᵀy − 2wᵀXᵀy + wᵀXᵀXw. Il gradiente rispetto a w è ∇L = −2Xᵀy + 2XᵀXw (si ricava come per le derivate parziali, termine per termine: la derivata di wᵀa è a, quella di wᵀMw per M simmetrica è 2Mw). Imporre ∇L = 0 dà di nuovo XᵀXw = Xᵀy: due strade diverse (geometria e analisi) portano alla stessa equazione. Se le colonne di X sono indipendenti, XᵀX è invertibile e w = (XᵀX)⁻¹Xᵀy.

rng = np.random.default_rng(5)
n, d = 500, 3
X = rng.normal(size=(n, d))
w_vero = np.array([2.0, -1.0, 0.5])
y = X @ w_vero + rng.normal(0, 0.3, size=n)

# 1) equazioni normali
w_normali = np.linalg.solve(X.T @ X, X.T @ y)
# 2) risolutore stabile (usa la SVD/QR)
w_lstsq = np.linalg.lstsq(X, y, rcond=None)[0]
# 3) scikit-learn
from sklearn.linear_model import LinearRegression
w_sk = LinearRegression(fit_intercept=False).fit(X, y).coef_

print("pesi veri         :", w_vero)
print("equazioni normali :", np.round(w_normali, 4))
print("np.linalg.lstsq   :", np.round(w_lstsq, 4))
print("scikit-learn      :", np.round(w_sk, 4))

residuo = y - X @ w_lstsq
print("Xᵀ r (deve essere ≈ 0, residuo ortogonale alle colonne):", np.round(X.T @ residuo, 8))
print("R² =", round(1 - (residuo ** 2).sum() / ((y - y.mean()) ** 2).sum(), 4))

Output ottenuto:

pesi veri         : [ 2.  -1.   0.5]
equazioni normali : [ 1.9888 -1.0144  0.5059]
np.linalg.lstsq   : [ 1.9888 -1.0144  0.5059]
scikit-learn      : [ 1.9888 -1.0144  0.5059]
Xᵀ r (deve essere ≈ 0, residuo ortogonale alle colonne): [-0. -0. -0.]
R² = 0.981

Cosa fa il codice. Genera 500 esempi con 3 caratteristiche e pesi veri (2, −1, 0.5), più rumore. Calcola i pesi in tre modi (equazioni normali, lstsq, scikit-learn) e ottiene lo stesso risultato, vicino ai pesi veri. L'ultima parte verifica la condizione di perpendicolarità: Xᵀr è quasi zero (a meno di arrotondamenti), cioè il residuo è ortogonale a tutte le colonne, esattamente come nel disegno. fit_intercept=False perché qui il modello non ha l'intercetta (per includerla basta aggiungere a X una colonna di uni).

Un avvertimento: le equazioni normali raddoppiano il condizionamento

Il condizionamento di XᵀX è il quadrato di quello di X: κ(XᵀX) = κ(X)². Se κ(X) = 10⁶, κ(XᵀX) = 10¹² e si perdono quasi tutte le cifre di un float64. Per questo, quando le colonne sono quasi collineari, si preferisce lstsq (basato su fattorizzazioni stabili).

Xc = np.column_stack([X[:, 0], X[:, 1], X[:, 0] + 1e-6 * rng.normal(size=n)])     # terza colonna ≈ prima
print("κ(X)   =", f"{np.linalg.cond(Xc):.2e}", "| κ(XᵀX) =", f"{np.linalg.cond(Xc.T @ Xc):.2e}", "(circa il quadrato)")
yc = Xc @ np.array([1.0, 1.0, 1.0]) + rng.normal(0, 0.3, size=n)
print("coefficienti (lstsq):", np.round(np.linalg.lstsq(Xc, yc, rcond=None)[0], 3))
print("somma dei due coefficienti collineari:", round(float(np.linalg.lstsq(Xc, yc, rcond=None)[0][[0, 2]].sum()), 3), "(la somma è ben determinata, i singoli no)")

Output ottenuto:

κ(X)   = 1.82e+06 | κ(XᵀX) = 3.33e+12 (circa il quadrato)
coefficienti (lstsq): [ 1.700241e+03  1.004000e+00 -1.698237e+03]
somma dei due coefficienti collineari: 2.004 (la somma è ben determinata, i singoli no)

Cosa fa il codice. Costruisce una matrice in cui la terza colonna è la prima più un rumore minuscolo (collinearità): le due colonne portano quasi la stessa informazione. Il numero di condizionamento esplode, e quello delle equazioni normali è circa il quadrato. I coefficienti di lstsq sulle due colonne collineari risultano instabili (cambiano molto con piccole variazioni dei dati), ma la loro somma è ben determinata: i dati possono dire quanto pesa «la direzione comune» ma non come ripartirla. Chi interpreta i coefficienti come «importanza» di una variabile deve tenerne conto.

Errori tipici (sintomo, causa, rimedio)

• * al posto di @. Moltiplica elemento per elemento invece del prodotto matriciale (lezione 5). Rimedio: @, e controllare le forme.
• Vettori (n,) e (n, 1) confusi. Il broadcasting produce matrici (n, n) senza errore. Rimedio: assert sulle forme.
• Invertire per risolvere. inv(A) @ b è lento e impreciso. Rimedio: solve.
• Testare det == 0 per la singolarità. Un determinante 1e-17 è «zero» ma non esattamente zero. Rimedio: rango (matrix_rank) o condizionamento.
• Ignorare la collinearità. Coefficienti enormi e instabili. Rimedio: rimuovere variabili ridondanti, regolarizzare.
• Ordine delle moltiplicazioni. AB ≠ BA; (AB)ᵀ = BᵀAᵀ. Rimedio: ragionare sulle dimensioni.
• Equazioni normali con matrici mal condizionate. Rimedio: lstsq, QR, SVD.

Tre errori su sette riguardano le forme degli array (* contro @, (n,) contro (n, 1), ordine dei fattori) e si prevengono scrivendo la shape accanto a ogni variabile. Gli altri quattro (inversa, determinante, collinearità, equazioni normali) hanno la stessa radice: trattare una matrice quasi singolare come se fosse ben condizionata, e per scoprirlo basta calcolare il numero di condizionamento prima di fidarsi del risultato.

Come lavorare nella pratica

Pensa in forme: scrivi accanto a ogni array la sua shape. Usa i risolutori (solve, lstsq) e non le inverse, controlla il condizionamento con np.linalg.cond prima di fidarti di una soluzione e standardizza le caratteristiche per migliorarlo. Verifica con il residuo (‖Ax − b‖ piccolo e Xᵀr ≈ 0 per i minimi quadrati) e disegna in 2D prima di passare a dimensioni alte.

Il lessico essenziale. Un vettore è insieme uno spostamento e una lista ordinata di numeri; la combinazione lineare è una somma di vettori scalati e lo span l'insieme di tutte quelle possibili. Il prodotto scalare è Σ aᵢbᵢ = |a||b|cos θ e la norma la lunghezza. La matrice rappresenta una trasformazione lineare, il determinante il fattore con cui moltiplica le aree, il rango il numero di colonne indipendenti, l'inversa la matrice che disfa A. Il pivot è il valore usato per eliminare un'incognita, il numero di condizionamento misura quanto la soluzione amplifica gli errori sui dati, le equazioni normali sono XᵀXw = Xᵀy e la collinearità è la quasi proporzionalità fra caratteristiche.

Mini-progetto

1. Implementa cosine_similarity_matrix(X) che restituisca la matrice n × n delle similarità del coseno tra tutte le righe, senza cicli (suggerimento: normalizza le righe e usa un solo prodotto matriciale X_norm @ X_norm.T). Verificala su 5 documenti inventati.
2. Estendi gauss perché restituisca anche il determinante (il prodotto dei pivot, con il segno cambiato a ogni scambio di righe) e confrontalo con np.linalg.det.
3. Con i dati del tumore al seno di scikit-learn, risolvi la regressione lineare del primo indice di caratteristiche sulle altre con le equazioni normali e con lstsq: confronta il numero di condizionamento di X e di XᵀX, con e senza standardizzazione.

Il primo punto è un'applicazione immediata del prodotto scalare: normalizzare le righe trasforma tutte le similarità in un solo prodotto matriciale, ed è il trucco che si usa per confrontare migliaia di embedding. Il secondo estende il codice di Gauss con il determinante, e il terzo ti fa misurare con i tuoi occhi che il condizionamento delle equazioni normali è il quadrato di quello di X.

Challenge

• Easy. Dimostra che due vettori non nulli sono perpendicolari se e solo se |a + b|² = |a|² + |b|² (teorema di Pitagora generalizzato).
• Medium. Dimostra che det(AB) = det(A)·det(B) per matrici 2 × 2 e interpretalo geometricamente (aree che si moltiplicano).
• Hard. Dimostra che κ(XᵀX) = κ(X)² usando il fatto che le «lunghezze di stiramento» di XᵀX sono i quadrati di quelle di X (le giustificheremo con la SVD).
• Expert. Implementa la proiezione ortogonale di y sullo span delle colonne di X con il metodo di Gram-Schmidt (ortogonalizzazione) e confrontala con X @ lstsq(X, y).
• Enterprise. Prendi un modello lineare della tua organizzazione con almeno 10 caratteristiche. Calcola κ(X) prima e dopo la standardizzazione, individua le coppie di caratteristiche con correlazione > 0.9, e riporta quanto cambiano i coefficienti rimuovendone una.

Le prime due sfide sono dimostrazioni brevi che si controllano con un caso numerico; la terza si chiude solo con la SVD della lezione successiva, e la quarta ti obbliga a scrivere la proiezione ortogonale invece di chiamare lstsq. L'ultima è la più utile al lavoro: confronta i coefficienti con e senza una variabile quasi duplicata, e vedrai quanto poco vale l'«importanza» letta dai pesi.

Che cosa sai fare adesso

Sai leggere un vettore come freccia e come riga di dati, usare il coseno per confrontare direzioni, interpretare una matrice dalle sue colonne e il determinante come fattore d'area, e risolvere un sistema con l'eliminazione di Gauss con pivoting invece che con l'inversa. Sai riconoscere un sistema mal condizionato con np.linalg.cond e derivare le equazioni normali sia con la geometria sia con il gradiente. La lezione seguente, autovalori, SVD e PCA, dà il significato preciso del numero di condizionamento; i minimi quadrati con molte variabili tornano in Regressione multivariata.

Nella stessa serie trovi anche: 22 - Analisi matematica 2: gradiente, integrali, discesa del gradiente (la lezione precedente), 25 - Calcolo matriciale e differenziazione automatica: dal Jacobiano a un autograd scritto da zero (che richiama questa lezione).