Risoluzione di equazioni lineari in

Informazioni su Risoluzione di equazioni lineari

Ricordi quando hai imparato l'algebra per la prima volta, al liceo?

Di solito, ci veniva dato un sistema di equazioni e ci veniva chiesto di risolvere rispetto a x e y.

4x - 3y = -2
2x + 7y = 16

Due equazioni, due incognite: un po' di sostituzione e presto vedrai che x = 1, y = 2.

Guarda di nuovo il lato sinistro dell'esempio precedente. Sembra una matrice [4 -3; 2 7] che moltiplica un vettore incognito [x, y] per dare il vettore [-2, 16].

Nella notazione di algebra lineare: A x = b, seguendo la solita convenzione delle lettere maiuscole per le matrici e minuscole per i vettori.

A x = 0 e lo spazio nullo

Un caso speciale importante è quando tutte le nostre equazioni hanno zero nel membro destro.

4x - 3y = 0
2x + 7y = 0

Nell'esempio sopra, l'unica soluzione è quando x = y = 0, il che di solito non è interessante.

Soluzioni non banali per A x = 0 esistono solo quando A ha determinante nullo.

julia> using LinearAlgebra

julia> A = [4 -3; 2 7]
2×2 Matrix{Int64}:
 4  -3
 2   7

julia> det(A)  # not zero!
34.0

Considera invece queste equazioni:

4x - 3y = 0
8x + 6y = 0

La seconda equazione è semplicemente il doppio della prima, non aggiunge nuove informazioni. Qualsiasi valore per cui x = 0.75y sarà una soluzione.

Questi valori (x, y) giacciono lungo una linea nello spazio bidimensionale. Ci sono infinite soluzioni, e formano lo spazio nullo della matrice.

In questo caso, le righe della matrice non sono linearmente indipendenti, e il rango della matrice è inferiore al numero di righe o colonne.

julia> A = [4 -3; 8 -6]
2×2 Matrix{Int64}:
 4  -3
 8  -6

julia> det(A)
0.0

julia> size(A)  # total number of rows and columns
(2, 2)

julia> rank(A)  # number of linearly-independent rows (or columns)
1

julia> nullspace(A)
2×1 Matrix{Float64}:
 -0.6
 -0.8

julia> nullspace(A) |> norm  # Julia returns unit vectors when possible
1.0

La funzione nullspace() restituisce un singolo punto, scelto per essere normalizzato a un vettore unitario, ma qualsiasi multiplo scalare di questo vettore è anche nello spazio nullo.

Risolvere A x = b

Ora supponi di avere 1000 equazioni in 1000 incognite? Se ti sembra assurdo, ricorda che i trasduttori digitali wireless oggi sono economici e versatili (misurano deformazioni, velocità del vento, accelerazione su 3 assi, qualsiasi cosa...). Averne 1000 che monitorano una struttura moderna come un ponte sospeso è perfettamente ragionevole. Serve un sistema di controllo che interpreti il flusso di dati, pronto a segnalare un allarme se la situazione diventa allarmante.

Di nuovo, A x = b, dove abbiamo A (dal progetto ingegneristico) e b (valori misurati dai trasduttori) ma dobbiamo trovare x.

Un primo pensiero, per chi non ha familiarità, è in qualche modo «dividere per A» per spostarlo al membro destro.

In realtà, A ha un'inversa (la maggior parte, ma non tutte, le matrici quadrate ce l'hanno: il determinante deve essere diverso da zero), e il calcolo funziona (lentamente!).

julia> A = [4 -3; 2 7]
2×2 Matrix{Int64}:
 4  -3
 2   7

julia> b = [-2, 16]
2-element Vector{Int64}:
 -2
 16

# the inverse matrix
julia> inv(A)
2×2 Matrix{Float64}:
  0.205882   0.0882353
 -0.0588235  0.117647

julia> inv(A) * b
2-element Vector{Float64}:
 0.9999999999999999
 2.0

Sfortunatamente, calcolare l'inversa è lento per matrici piccole e glaciale per quelle grandi. Non è l'ideale, se il tuo ponte crolla durante il calcolo!

Fortunatamente, altri algoritmi sono molto più veloci (in questo caso l'eliminazione di Gauss).

Julia (copiando Matlab) usa semplicemente una barra rovesciata per il risolutore (tecnicamente, «divisione a sinistra»).

julia> x = A \ b
2-element Vector{Float64}:
 1.0
 2.0

Questo è un esempio banale, ma ci sono un paio di dettagli a cui fare attenzione.

• I coefficienti devono allinearsi in colonne, quindi inserisci degli zeri dove servono per i termini mancanti.
• Le equazioni (e quindi le righe della matrice) devono essere linearmente indipendenti. Se la riga 3 è semplicemente la somma delle righe 1 e 2, non aggiunge nuove informazioni e il problema è sottodeterminato.

Il comando rank() dà il numero di righe/colonne linearmente indipendenti, quindi cerca di farlo coincidere con il numero di incognite.

julia> rank(A)
2

Matrici rettangolari

Probabilmente, da preadolescente, ti è stato insegnato che serve un sistema di N equazioni per risolvere N incognite.

Come per la maggior parte delle cose in matematica, la realtà è un po' più sfumata (non solo il dettaglio dell'indipendenza lineare, menzionato nella sezione precedente).

Righe «troppo poche»

Con N-1 equazioni (righe nella matrice A), il problema è sottodeterminato. Ci arrendiamo alla disperazione?

Dipende!

Le N-1 equazioni contengono comunque molte informazioni. In termini geometrici, ci servono N per determinare la soluzione come un punto nello spazio N-dimensionale, ma N-1 ci dirà che la soluzione deve trovarsi da qualche parte su una linea.

Questo ci lascia con un numero infinito di soluzioni, ma anche una quantità infinita di spazio dei parametri con nessuna soluzione.

La tua linea di soluzioni passa attraverso una regione dello spazio dei parametri di cui un ingegnere rispettabile dovrebbe preoccuparsi? O che potrebbe giustificare la chiusura della struttura al pubblico (i ponti sospesi possono diventare vivaci con forti venti trasversali)? Forse vale la pena fare qualche calcolo in più per determinarlo!

julia> A3 = [4 -3 1; 2 2 2; 3 -4 -3]
3×3 Matrix{Int64}:
 4  -3   1
 2   2   2
 3  -4  -3

julia> b3 = [1, 12, -14]
3-element Vector{Int64}:
   1
  12
 -14

# fully-determined solution
julia> A3 \ b3
3-element Vector{Float64}:
 1.0
 2.0
 3.0

# remove row 3 from A3 and B3
julia> A2 = A3[1:2, :]
2×3 Matrix{Int64}:
 4  -3  1
 2   2  2

julia> b2 = b3[1:2]
2-element Vector{Int64}:
  1
 12

# an under-determined solution
julia> z1 = A2 \ b2
3-element Vector{Float64}:
 1.594594594594594
 2.445945945945947
 1.9594594594594592

Il risolutore con barra rovesciata ci dà una soluzione! Ufficiosamente, questa è apparentemente la soluzione con la norma più piccola (la più vicina all'origine, geometricamente).

Per ottenere le altre soluzioni, ci serve il nullspace di A2.

Aggiungere qualsiasi multiplo scalare dello spazio nullo alla soluzione originale darà un'altra soluzione valida.

# get the nummspace of A2
julia> N = nullspace(A2)
3×1 Matrix{Float64}:
 -0.46499055497527725
 -0.3487429162314578
  0.813733471206735

# add a random multiple of the nullspace to the earlier solution
julia> z = z1 + rand() * N
3×1 Matrix{Float64}:
 1.274331860650258
 2.205748895487695
 2.5199192438620472

# our random z is a valid solution, recovering b2 = [1, 12]
julia> A2 * z
2×1 Matrix{Float64}:
  0.9999999999999947
 12.0

Allo stesso modo, N-2 equazioni vincoleranno le soluzioni (infinite) a un piano, e valgono gli stessi principi.

Righe «troppe»

La situazione opposta si presenta quando ci sono più equazioni che incognite, eppure in qualche modo sono ancora linearmente indipendenti.

Nell'ingegneria del mondo reale, questo è del tutto normale, ed è visto come una cosa positiva!

I trasduttori hanno precisione limitata, il vettore b ha precisione limitata, e c'è rumore nel calcolo. Ora la soluzione è un adattamento ai minimi quadrati dei dati rumorosi.

La tecnica più semplice usa la pseudoinverse della matrice, che Julia implementa come funzione pinv().

Può essere usata in modo molto simile all'inversa di una matrice quadrata (non singolare), dando una stima ai minimi quadrati delle variabili invece di una soluzione esatta.

# create 5-row A and b for 2 variables
julia> A5 = [4 -3; 2 7; -1 2; 1 -1; 3 1]
5×2 Matrix{Int64}:
  4  -3
  2   7
 -1   2
  1  -1
  3   1

julia> rank(A5)
2

# add some random-normal noise to b5
julia> b5n = [-2, 16, 3, -1, 5] + randn(5) * 0.1
5-element Vector{Float64}:
 -1.9337349223581917
 16.024913008059457
  3.0087646177373992
 -0.8675748721176717
  4.914568047308805

# solve to get x ≈ 1, y ≈ 2, using the pseudoinverse
julia> pinv(A5) * b5n
2-element Vector{Float64}:
 1.006117942769543
 1.9962973764508272

Autovalori e autovettori

La sezione precedente descriveva le soluzioni per A x = b, equivalenti all'algebra familiare.

Questa sezione riguarda le soluzioni di A x = λ x, dove λ è uno scalare.

Questo è un concetto tipico dell'algebra lineare, ma possiamo interpretarlo geometricamente:

Per valori opportuni di x e λ, la matrice quadrata A scala il vettore non nullo x di un fattore λ in lunghezza, senza cambiarne la direzione (a parte un λ negativo che la inverte).

Sembra di nicchia, ma si scopre che è assurdamente utile!

Terminologia: I valori validi di λ sono gli autovalori di A, e i valori corrispondenti di x sono gli autovettori di A. Sfortunatamente, dobbiamo convivere con parole che cambiano lingua a metà.

Agli studenti di solito viene insegnato come calcolare autovalori/autovettori per matrici 2×2 a mano, ma usare un computer è molto più semplice (anche se resta piuttosto lento, e scala come O(n^3) per una matrice n×n nel caso generale).

julia> A = rand(-9:9, 2, 2)
2×2 Matrix{Int64}:
 9  5
 3  2

julia> F = eigen(A)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
2-element Vector{Float64}:
  0.27984674554472466
 10.720153254455274
vectors:
2×2 Matrix{Float64}:
 -0.497417  0.945605
  0.867511  0.325317

# eigenvalues
julia> F.values
2-element Vector{Float64}:
  0.27984674554472466
 10.720153254455274

# each column is an eigenvector, normalized to a unit vector
julia> F.vectors
2×2 Matrix{Float64}:
 -0.497417  0.945605
  0.867511  0.325317

In generale, una matrice n×n avrà n autovalori, anche se non sempre distinti. Pensali come le radici di un polinomio di grado n (chiamato polinomio caratteristico), che possono essere ripetute, e sono spesso complesse anche per una matrice a valori reali.

Ogni autovettore rappresenta una direzione, e qualsiasi multiplo scalare è anch'esso un autovettore valido. Per comodità nei calcoli successivi, Julia restituisce vettori unitari con norma 1.

Applicazioni

Spiegare perché gli autovettori sono importanti è idealmente un compito per un manuale di 500 pagine, piuttosto che per qualche breve paragrafo. Questo argomento permea così tanta della matematica applicata moderna.

A livello generale, gli autovettori rappresentano gli assi (direzioni) «più importanti» in un set di dati. Gli autovalori indicano l'«importanza relativa» di ciascun asse (soggetta ad alcune ipotesi sulla normalizzazione degli input).

Ciò che questo significa realmente dipende dall'applicazione.

Analisi delle componenti principali

In un tipico problema di data science, potremmo avere 100 o più «caratteristiche», memorizzate come colonne di dati. Quasi inevitabilmente, ci saranno rumore, ridondanza e correlazioni indesiderate.

Dobbiamo fare riduzione dimensionale, e la PCA è un modo per portare ordine nel caos:

• Calcola la matrice di covarianza dei dati.
• Ora abbiamo una matrice quadrata, quindi calcola gli autovalori e gli autovettori.
• Ordina gli autovettori in ordine decrescente rispetto ai loro autovalori (più precisamente, il valore assoluto |λ|).
• I primi k autovettori sono ora le principal components, dove k è significativamente più piccolo del numero originale di caratteristiche.
• Proietta l'intero set di dati su questi k assi, e inizia a cercare schemi interessanti.

In questa forma, la PCA è stata usata da gruppi diversi come scienziati medici che cercano di migliorare le proprietà molecolari nella progettazione di farmaci, e attivisti politici che cercano di capire le preferenze degli elettori dai dati dei sondaggi. Probabilmente anche da gruppi di marketing che cercano di venderti cose, ma qualsiasi tecnologia può essere usata per il bene o per il male.

La PCA non fa parte di un'installazione minima di Julia, ma (al di fuori di Exercism) il pacchetto MultivariateStats contiene ciò che serve.

Elaborazione delle immagini

Portando l'idea della PCA un passo avanti, le immagini digitali sono solo matrici di valori di pixel, e possiamo calcolare le componenti principali per esse.

Questo è stato usato per decenni nella compressione delle immagini, dove la PCA è una guida su cosa è più importante preservare quando si riduce la dimensione del file.

Sempre più spesso, la PCA è una parte vitale dei classificatori di immagini come il riconoscimento facciale. La prossima volta che attraversi un aeroporto, Grande Fratello non ti sta solo guardando: usa anche l'algebra lineare per capire ciò che vede!

Ingegneria meccanica

In modo meno controverso, gli assi principali di un componente meccanico sono gli autovettori del suo tensore dei momenti d'inerzia I (una matrice con un nome leggermente diverso).

Ogni ruota della tua auto probabilmente ha uno o più contrappesi, regolati per azzerare gli elementi fuori diagonale in questo tensore I. Il meccanico dell'officina senza dubbio evita di fare i calcoli (a differenza dei progettisti di aerei e degli ingegneri aerospaziali), ma «oscillazione» è solo una parola quotidiana per descrivere i termini incrociati nel tensore, che in questo caso renderebbero il viaggio meno confortevole e aumenterebbero l'usura meccanica sui cuscinetti.