Risoluzione di sistemi lineari

Un sistema di equazioni lineari può essere scritto in forma matriciale come Ax = b, dove A è la matrice dei coefficienti, x è il vettore delle incognite e b è il vettore dei termini noti. Ad esempio, due equazioni con due incognite possono essere espresse come un'equazione matriciale 2×2. Risolvere il sistema significa trovare x tale che A moltiplicata per x sia uguale a b. np.linalg.solve(A, b) di NumPy gestisce questa operazione in modo efficiente e accurato.

Dal punto di vista matematico, la soluzione di Ax = b è x = A⁻¹b. Tuttavia, calcolare np.linalg.inv(A) @ b è più lento e meno preciso che chiamare direttamente np.linalg.solve(A, b). solve() usa internamente la decomposizione LU, che richiede meno operazioni e accumula meno errori di arrotondamento. Regola generale: se sta risolvendo x in Ax=b, usi sempre solve, mai inv.

np.linalg.solve(A, b) richiede che A sia quadrata e non singolare. Se A è singolare (il determinante è zero), la funzione solleva np.linalg.LinAlgError: Singular matrix. A è singolare quando le righe sono linearmente dipendenti, ad esempio quando un'equazione è semplicemente un multiplo di un'altra. Verifichi sempre che il sistema sia ben posto prima di chiamare solve. Se A non è quadrata, usi invece np.linalg.lstsq().

np.linalg.solve applica internamente la decomposizione LU: fattorizza A in una matrice triangolare inferiore L e una matrice triangolare superiore U, tali che A = LU. La risoluzione diventa quindi la risoluzione di due sistemi triangolari più semplici (sostituzione in avanti e all'indietro), molto più efficienti dell'eliminazione di Gauss ripetuta da zero. È lo stesso algoritmo usato da MATLAB, Fortran LAPACK e praticamente tutte le librerie di calcolo scientifico.

np.linalg.solve(A, B) può gestire una matrice B con più colonne, risolvendo simultaneamente il sistema per ogni colonna. Se B ha forma (n, k), il risultato X ha forma (n, k), dove ogni colonna di X è la soluzione corrispondente alla colonna di B. Questo è utile nelle simulazioni fisiche e nel machine learning, quando è necessario risolvere la stessa matrice dei coefficienti con molti termini noti diversi.

Quando il sistema è sovradeterminato (con più equazioni che incognite), non esiste una soluzione esatta. np.linalg.lstsq(A, b) trova il valore di x che minimizza la somma dei quadrati dei residui ||Ax - b||²: è esattamente ciò che calcola la regressione lineare. Restituisce la soluzione x, i residui, il rango e i valori singolari. L'argomento rcond=None evita un avviso di deprecazione e imposta la soglia per trattare i valori singolari come zero.

Le equazioni normali dei minimi quadrati ordinari sono A.T @ A @ x = A.T @ b, un sistema quadrato risolvibile con np.linalg.solve. Il risultato coincide con quello di lstsq quando A.T @ A è ben condizionata. Nella pratica, lstsq è preferibile perché usa la SVD ed è stabile anche quando le colonne sono quasi collineari; tuttavia, le equazioni normali aiutano a comprendere la matematica alla base della regressione.

Prima di risolvere il sistema, è buona pratica controllare il numero di condizionamento con np.linalg.cond(A). Un numero di condizionamento vicino a 1 indica che il sistema è ben condizionato e che la soluzione è affidabile. Un numero di condizionamento molto grande (ad esempio, 1e12) indica che la matrice è quasi singolare: piccole variazioni in b producono grandi variazioni in x e la soluzione potrebbe essere numericamente inaffidabile. In tal caso, valuti la regolarizzazione (regressione ridge) o l'uso di un risolutore più robusto.

Per sistemi molto grandi in cui A contiene prevalentemente valori nulli (una matrice sparsa), memorizzare l'intera matrice spreca memoria e usare np.linalg.solve denso è lento. SciPy fornisce scipy.sparse e scipy.sparse.linalg.spsolve, che sfruttano la sparsità per gestire sistemi enormi, come quelli dei modelli a elementi finiti o dei problemi di rete. Capire quando una matrice è sparsa è un aspetto importante delle prestazioni nel calcolo scientifico.

Una volta ottenuto il vettore soluzione x da np.linalg.solve(A, b), lo verifichi sempre calcolando A @ x e confrontandolo con b tramite np.allclose(A @ x, b). Nell'aritmetica in virgola mobile, raramente si ottiene un'uguaglianza esatta: allclose verifica che i valori rientrino in una piccola tolleranza. Il residuo ||Ax - b|| dovrebbe essere prossimo all'epsilon di macchina moltiplicato per la scala del problema, se il sistema è ben condizionato.

I risolutori lineari sono presenti in tutto il machine learning: la regressione lineare risolve le equazioni normali, i processi gaussiani risolvono sistemi con matrici kernel e i filtri di Kalman risolvono le equazioni di aggiornamento dello stato. Comprendendo np.linalg.solve e lstsq, può implementare questi algoritmi da zero e analizzarli quando i risultati non sono quelli attesi. Il modulo linalg di NumPy offre gli stessi strumenti numerici utilizzati nelle librerie scientifiche di produzione.

Verifichi la sua comprensione dei concetti di analisi dei dati trattati in questa lezione.

In questa lezione ha imparato che: np.linalg.solve(A, b) risolve in modo efficiente i sistemi lineari quadrati utilizzando la decomposizione LU, np.linalg.lstsq() trova la soluzione ai minimi quadrati per i sistemi sovradeterminati e il numero di condizionamento misura la stabilità numerica: valori elevati segnalano matrici quasi singolari. Ora esamineremo gli autovalori e la SVD, alla base della PCA e di molti algoritmi di ML.