Basi di algebra lineare in

Informazioni su Basi di algebra lineare

Che cos'è l'«algebra lineare»?

Esistono molte definizioni tecniche, su Wikipedia, nei libri di testo e su siti eccellenti come 3Blue1Brown, ma per i nostri scopi possiamo essere più informali.

I maintainer di Julia (e di Python) adorano questo genere di cose. La dirigenza di Exercism, non altrettanto.

Questo concetto sarà limitato agli aspetti più semplici di un argomento enorme, ma attenzione: si tratta, inevitabilmente, di un concetto piuttosto matematico.

Che cos'è un vettore?

Dipende da chi lo chiedi e da come vuoi visualizzare qualcosa di piuttosto astratto.

• In Julia, Vector{T} è solo un alias per Array{T, 1}, e l'eltype T può essere qualsiasi cosa.
• In (buona parte della) fisica, un vettore è una freccia con una length e una direction in uno spazio a N dimensioni, ma senza una posizione fissa.
• In algebra lineare, la freccia dei fisici è fissata nello spazio, con la coda nell'origine. Questo rende superflua l'asta della freccia, quindi possiamo rappresentare il vettore come un point in uno spazio a N dimensioni, dove gli N elementi rappresentano la distanza dall'origine lungo ciascuno degli N assi.

Ai fini di questa trattazione, ignora i vettori di stringhe o di caratteri. In questo concetto lavoreremo con tipi numerici: Int, Float o Complex.

Che cos'è una matrice?

Anche qui troverai risposte diverse.

• Un array bidimensionale rettangolare di numeri (senza bordi irregolari). Una matrice quadrata è un caso speciale comune.
• Una linear combination di vettori colonna, affiancati uno accanto all'altro.
• Una transformation che può essere applicata a un vettore, proprio come una function viene applicata ad altri tipi di input.

Alcuni tipi di matrice quadrata sono abbastanza comuni da avere dei nomi speciali.

Matrice diagonale: Tutti gli elementi diversi da zero si trovano sulla main diagonal (dall'angolo in alto a sinistra a quello in basso a destra).

julia> [1 0 0; 0 2 0; 0 0 3]
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

# a convenient shortcut
julia> diagm(1:3)
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

Matrice identità: Una matrice diagonale con soltanto uni sulla diagonale (per motivi che diventeranno più chiari in seguito). Spesso abbreviata con I.

julia> Matrix{Float64}(I, 3, 3)
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

Triangolare superiore: Valori diversi da zero sulla diagonale e sopra di essa, zeri sotto. Triangolare inferiore puoi immaginarlo.

Trasporre una matrice

Se vuoi davvero scambiare le righe con le colonne, la funzione permutedims() lo fa e restituisce una nuova matrice.

Questo comporta una copia, che è lenta e avida di memoria quando si lavora con matrici di grandi dimensioni.

Ai fini dell'algebra lineare, la funzione transpose() è più utile, perché crea rapidamente un wrapper lazy attorno alla matrice originale.

Ancora più utile è la funzione adjoint(), che cambia anche il segno della parte immaginaria dei numeri complessi (se ti chiedi perché sia così importante, una risposta rapida è la meccanica quantistica). Questa operazione è abbastanza comune da poter semplicemente aggiungere un apostrofo ' al nome della variabile per creare l'aggiunto.

julia> m
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

# new, full copy
julia> permutedims(m)
3×2 Matrix{Int64}:
 1  4
 2  5
 3  6

# lazy version
julia> transpose(m)
3×2 transpose(::Matrix{Int64}) with eltype Int64:
 1  4
 2  5
 3  6

julia> mc = [1+2im 2+3im; 3+2im 1+2im]
2×2 Matrix{Complex{Int64}}:
 1+2im  2+3im
 3+2im  1+2im

# lazy conjugate transpose
julia> adjoint(mc)
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

# syntactic sugar for adjoint
julia> mc'
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

Moltiplicazione

Il resto di questo documento include funzionalità del modulo LinearAlgebra. Una riga using LinearAlgebra lo porterà nel namespace, ma non continueremo a ripeterlo negli esempi (troppo disordine visivo).

Vettori, prodotto elemento per elemento

Ne abbiamo parlato nel concetto Operazioni sui vettori. L'operatore è .*, che agisce sugli input a coppie di elementi, restituendo un output della stessa dimensione e dello stesso tipo degli input.

julia> [1, 2] .* [3, 4]
2-element Vector{Int64}:
 3
 8

Vettori, prodotto scalare

Questa operazione estremamente comune viene scritta nei libri di testo come u ⋅ v ed è chiamata «prodotto scalare».

Un prodotto scalare equivale alla somma del prodotto elemento per elemento.

Due vettori si possono moltiplicare con il solito operatore *, ma solo se il vettore di sinistra è un adjoint: comodamente scritto u' * v. Questo dovrebbe diventare più chiaro nella sezione successiva sulla moltiplicazione tra matrici.

Il punto sopraelevato (centrato) è disponibile in Julia (si digita \cdot e poi tab) come zucchero sintattico per la funzione dot(). Non c'è bisogno di specificare l'aggiunto, perché questo dettaglio viene gestito automaticamente.

julia> using LinearAlgebra

# with element-wise syntax
julia> sum([1, 2] .* [3, 4])
11

# with u' * v syntax
julia> [1, 2]' * [3, 4]
11

# with dot()
julia> dot([1, 2], [3, 4])
11

# with \cdot syntax
julia> [1, 2] ⋅ [3, 4]
11

Per i vettori a valori complessi, il vettore di sinistra deve essere il coniugato (segno invertito sulla parte immaginaria). La funzione dot lo fa automaticamente, la sintassi u' lo fa esplicitamente, ma sum(u .* v) fallirà se u e v sono complessi.

Vettori, prodotto vettoriale

Un po' di nostalgia per chiunque abbia seguito in passato un corso di elettricità e magnetismo! E anche per gli ingegneri che hanno familiarità con il calcolo del momento torcente o dei vettori del momento angolare.

Mentre il prodotto scalare trasforma due vettori in uno scalar, il prodotto vettoriale trasforma due vettori a 3 componenti in un terzo vettore a 3 componenti.

Nella rappresentazione geometrica dello spazio vettoriale, il nuovo vettore è perpendicolare al piano che contiene i due vettori di input. Se i due input sono paralleli (a meno del segno), non definiscono un piano, quindi l'output sarà [0, 0, 0].

L'ordine conta: u × v == -(v × u). Questa è la famosa regola della mano destra, che ha lasciato molti di noi a fissare il pollice e due dita mentre li ruotavamo nello spazio (spesso con un'espressione perplessa).

# with \times syntax
julia> [1, 2, 3] × [3, 4, 5]
3-element Vector{Int64}:
 -2
  4
 -2
# with cross()
julia> cross([1, 2, 3], [3, 4, 5])
3-element Vector{Int64}:
 -2
  4
 -2

Nota che questa operazione è limitata ai vettori di lunghezza 3 (equivalenti allo spazio euclideo con assi x, y, z ortogonali). I matematici sono felici di definire l'operazione per altre dimensioni, ottenendo come risultato un'incomprensibile insalata di parole (tensori di ordine superiore, prodotto wedge, prodotto esterno, multivettori...). A quel punto, la maggior parte di noi scappa o cambia rapidamente argomento!

Norme

Quanto è «grande» un vettore?

La norm è un tentativo di catturare questo, riducendo il vettore a uno scalare appropriato.

Esiste tutta una famiglia di norme, ma di gran lunga la più comune è la norma 2, che è √(v ⋅ v).

Questa operazione di radice della media dei quadrati è la distanza pitagorica dall'origine (in uno spazio a N dimensioni). Se visualizziamo il vettore come una freccia, con la coda nell'origine, la norma 2 è la lunghezza della freccia.

Qualsiasi norma p si può calcolare fornendo p come secondo argomento. La norma 1 a volte è utile: è semplicemente la somma dei valori assoluti degli elementi (quindi molto rapida e facile da calcolare).

# defaults to the 2-norm
julia> norm([1, 2, 3])
3.7416573867739413
# the 1-norm
julia> norm([1, -2, 3], 1)
6.0

Moltiplicazione tra matrici

A volte sembra che ogni matematico applicato del mondo abbia passato gran parte degli ultimi 80 anni a trasformare ogni calcolo in una serie di moltiplicazioni tra matrici.

È un'operazione in cui i computer sono molto bravi:

• È altamente ripetitiva.
• Può essere parallelizzata in modo efficiente.
• È stato sviluppato molto hardware specializzato per renderla più veloce, dal Cray-1 da 80 milioni di dollari negli anni '70 alla GPU che probabilmente è integrata nel tuo portatile.

I dettagli sono piuttosto semplici, anche se non molto intuitivi a prima vista.

Considera una matrice A che moltiplica un vettore v per dare un output w (per convenzione dell'algebra lineare, useremo lettere maiuscole per le matrici e minuscole per i vettori).

La prima riga di A viene moltiplicata scalarmente per v per ottenere il primo elemento di w, la seconda riga dà il secondo elemento, e così via.

julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4
julia> v = [5, 6]
2-element Vector{Int64}:
 5
 6
julia> A * v
2-element Vector{Int64}:
 17  # equals [1, 2] ⋅ [5, 6]
 39  # equals [3, 4] ⋅ [5, 6]

Una moltiplicazione matrice * matrice estende questo a tutte le colonne della matrice di destra.

Per C = A * B, possiamo considerare che A moltiplica ogni colonna di B per dare la colonna corrispondente di C: una serie di moltiplicazioni matrice-vettore.

In modo equivalente, potremmo dire che la prima riga di A viene moltiplicata scalarmente per ogni colonna di B per dare la prima riga di C, la seconda riga dà la seconda, e così via.

Esistono diverse rappresentazioni mentali di questo tipo, e gli studenti interessati possono guardare un'intera lezione del MIT che ne parla.

julia> A
2×2 Matrix{Int64}:
 1  2
 3  4

julia> B = [5 7; 6 8]
2×2 Matrix{Int64}:
 5  7
 6  8

# left column is the same as A*v previously
julia> A * B
2×2 Matrix{Int64}:
 17  23
 39  53

julia> B * A
2×2 Matrix{Int64}:
 26  38
 30  44

Come mostrato nell'esempio sopra, la moltiplicazione tra matrici non è commutativa: non c'è una relazione semplice tra A*B e B*A.

La moltiplicazione tra matrici è probabilmente difficile da visualizzare per chi la incontra per la prima volta, solo leggendo le parole. YouTube ha molti video che la mostrano graficamente, quindi cerca «matrix multiplication» e scegli quello con lo stile, il livello di dettaglio e la lingua che preferisci.

Dimensioni

Il prodotto scalare di due vettori si basa sul fatto che abbiano la stessa lunghezza.

Di conseguenza, per la moltiplicazione tra matrici il numero di colonne della matrice di sinistra deve corrispondere al numero di righe della matrice di destra.

Esprimendo le dimensioni come tuple (nrows, ncols), come restituite da size(A), abbiamo (a, b) * (b, c) -> (a, c). Le dimensioni «interne», qui b e b, sono compatibili con l'esecuzione dei prodotti scalari. Le dimensioni «esterne», qui a e c, determinano le dimensioni dell'output.

Un esempio con matrici rettangolari:

julia> D = reshape(1:6, 2, 3)
2×3 reshape(::UnitRange{Int64}, 2, 3) with eltype Int64:
 1  3  5
 2  4  6

julia> E = reshape(1:12, 3, 4)
3×4 reshape(::UnitRange{Int64}, 3, 4) with eltype Int64:
 1  4  7  10
 2  5  8  11
 3  6  9  12

julia> D * E
2×4 Matrix{Int64}:
 22  49   76  103
 28  64  100  136

# (3, 4) * (2, 3) not possible
julia> E * D
ERROR: DimensionMismatch: matrix A has axes (Base.OneTo(3),Base.OneTo(4)), matrix B has axes (Base.OneTo(2),Base.OneTo(3))

Moltiplicare coppie di vettori in stile matriciale ha due possibilità.

Per convenzione, useremmo u' * v come equivalente del prodotto scalare. Le dimensioni sono (1, 3) * (3, 1), e Julia semplifica l'output (1, 1) in uno scalare (a differenza, per esempio, di R).

In alternativa, potremmo usare u * v', con dimensioni (3, 1) * (1, 3) -> (3, 3), per moltiplicare tutte le possibili coppie di elementi ed espandere i vettori in una matrice.

julia> u = [1, 2]
2-element Vector{Int64}:
 1
 2
julia> v = [3, 4]
2-element Vector{Int64}:
 3
 4
julia> u' * v
11
julia> u * v'
2×2 Matrix{Int64}:
 3  4
 6  8

u' * v è talvolta chiamato prodotto interno.

Il (meno comune) u * v' è corrispondentemente il prodotto esterno. Questo è legato al prodotto tensoriale, anche se va ben oltre i nostri scopi.

Rotazioni

Un tipo particolarmente comune di moltiplicazione tra matrici coinvolge le matrici di rotazione.

In 2D, esiste una matrice relativamente semplice per ruotare un vettore in senso antiorario di θ radianti.

julia> rot2d(θ, vec) = [cos(θ) -sin(θ); sin(θ) cos(θ)] * vec
rot2d (generic function with 1 method)

# unit vector in the x direction
julia> i_hat = [1, 0]
2-element Vector{Int64}:
 1
 0

# rotate 45 degrees
julia> rot2d(π/4, i_hat)
2-element Vector{Float64}:
 0.7071067811865476
 0.7071067811865475

# rotate 90 degrees -> unit vector in the y direction, j_hat
julia> rot2d(π/2, i_hat)
2-element Vector{Float64}:
 6.123233995736766e-17  # zero, within numerical error
 1.0

Le rotazioni in 3D richiedono una matrice più complessa, con due angoli. È difficile da mostrare in un documento Markdown senza incorporare molta LaTeX, quindi controlla Wikipedia per la formula.

Per la grafica 3D, come OpenGL e i suoi successori, la convenzione è lavorare con le coordinate omogenee.

Ogni vertice è rappresentato da un vettore a 4 componenti (o equivalentemente da una colonna di matrice): [x, y, z, 1.0] per un punto in (x, y, z).

Questo consente matrici di trasformazione più elaborate. La rotazione è ancora in A[1:3, 1:3], le traslazioni sono [Δx, Δy, Δz] in A[1:3, 4], il ridimensionamento è sulla diagonale, e ci sono altre possibilità per inclinazione, prospettiva, ecc.

Prestazioni

Nei dipartimenti di informatica di tutto il mondo ci sono molte tesi di dottorato con uno o più capitoli su piccoli miglioramenti incrementali negli algoritmi per la moltiplicazione tra matrici. Questo è davvero importante, e varie organizzazioni sono disposte a finanziare la ricerca per il proprio tornaconto.

Le prestazioni sono un argomento vasto che non possiamo approfondire, ma per illustrare possiamo provare a moltiplicare matrici casuali di varie dimensioni.

Il codice qui sotto è una stima rapida e approssimativa (usa BenchmarkTools.jl per un approccio migliore).

Il sistema usato era un piccolo PC da 430 $ (US): processore Ryzen 9, 32 GB di RAM, Linux Mint 22.1, Julia 1.11.6.

julia> mmul(A) = A * A
mmul (generic function with 1 method)
julia> A = rand(Float64, 1_000, 1_000);
julia> @time mmul(A);
  0.011055 seconds (3 allocations: 7.629 MiB)
julia> A = rand(Float64, 10_000, 10_000);
julia> @time mmul(A);
  7.236284 seconds (1.61 k allocations: 763.027 MiB, 17.92% gc time, 0.15% compilation time)

All'incirca, una coppia di matrici da un milione di elementi ha richiesto qualche millisecondo, matrici da 100 milioni di elementi hanno richiesto qualche secondo. Aggiungere diversi altri ordini di grandezza richiederà hardware migliore...