Algoritmo Gram-Schmidt

Fondamentale nell'algebra lineare, che permette di trasformare un una base qualsiasi in una base ortonormale, è l'algoritmo di Gram-Schmidt.
Nonostante porti i nomi di Jørgen Pedersen Gram e Erhard Schmidt, l'algoritmo ha radici più antiche. Versioni primordiali del processo apparvero già nei lavori di Laplace e Cauchy. Inizialmente, serviva a risolvere problemi di minimi quadrati in astronomia e geodesia. Jørgen Gram (1883), matematico attuariale, pubblicò il metodo in un saggio sulla scomposizione delle serie numeriche. Per lui era uno strumento per approssimare dati statistici. Erhard Schmidt (1907), allievo di Hilbert, ebbe il merito di generalizzare l'algoritmo agli spazi di funzioni a dimensione infinita, gli Spazi di Hilbert. Fu Schmidt a capire che questo processo era la chiave per rendere rigorosa la teoria di Fourier.

L'idea centrale è la sottrazione delle proiezioni. Se voglio rendere un vettore \(v_2\) perpendicolare a \(v_1\), devo togliere da \(v_2\) la sua "ombra" (proiezione) lungo \(v_1\).
Sia \(\{v_1, v_2, \dots, v_n\}\) una base di partenza. Vogliamo ottenere una base ortogonale \(\{u_1, u_2, \dots, u_n\}\).

Il primo vettore non ha nulla a cui essere perpendicolare, quindi lo prendiamo così com'è: \[u_1 = v_1\]
per rendere il secondo ortogonale gli sottraiamo la sua proiezione su \(u_1\): \[u_2 = v_2 - \text{proj}_{u_1}(v_2) = v_2 - \frac{\langle v_2, u_1 \rangle}{\langle u_1, u_1 \rangle} u_1\]
Per il terzo vettore dobbiamo sottrarre le proiezioni su tutti i vettori ortogonali già trovati \(u_1\) e \(u_2\): \[u_3 = v_3 - \frac{\langle v_3, u_1 \rangle}{\langle u_1, u_1 \rangle} u_1 - \frac{\langle v_3, u_2 \rangle}{\langle u_2, u_2 \rangle} u_2\]
Per ogni altro vettore \(k\): \[u_k = v_k - \sum_{j=1}^{k-1} \text{proj}_{u_j}(v_k)\]
Dopo aver ottenuto la base ortogonale \(\{u_1, \dots, u_n\}\), la rendiamo ortonormale dividendo ogni vettore per la sua norma (lunghezza):\[e_k = \frac{u_k}{\|u_k\|}.\]

def gram_schmidt(A):
    """
    A: sympy Matrix (n x n), colonne = vettori indipendenti
    Ritorna: sympy Matrix (n x n) con colonne ortonormali
    """
    n = A.cols
    ortho = []

    for k in range(n):
        v = A.col(k)

        # sottrai proiezioni sui vettori precedenti
        for u in ortho:
            v = v - (v.dot(u) / u.dot(u)) * u

        ortho.append(v)

    # normalizzazione
    ortho_norm = [u / sqrt(u.dot(u)) for u in ortho]

    # ricostruzione matrice
    Q = Matrix.hstack(*ortho_norm)
    return simplify(Q)
    
A = Matrix([
    [1, 1, 0],
    [1, 0, 1],
    [0, 1, 1]
])

gram_schmidt(A)

\(\left[\begin{matrix}\frac{\sqrt{2}}{2} & \frac{\sqrt{6}}{6} & - \frac{\sqrt{3}}{3}\\\frac{\sqrt{2}}{2} & - \frac{\sqrt{6}}{6} & \frac{\sqrt{3}}{3}\\0 & \frac{\sqrt{6}}{3} & \frac{\sqrt{3}}{3}\end{matrix}\right]\)

L’algoritmo è di seguito applicato alla matrice
\(A =\)
ottenendo

Sebbene elegante sulla carta, il processo di Gram-Schmidt "classico" ha un punto debole: gli errori di arrotondamento nei computer.
Sottraendo proiezioni molto simili tra loro, si possono accumulare piccoli errori che portano i vettori finali a non essere perfettamente perpendicolari (perdita di ortogonalità).
Per ovviare a questo, in informatica e ingegneria si usa il Gram-Schmidt Modificato, che ricalcola le proiezioni in modo più stabile, o altri metodi come le Riflessioni di Householder.

def gram_schmidt_modified(A):
    """
    A: sympy Matrix (n x n), colonne = vettori indipendenti
    Ritorna: sympy Matrix (n x n) con colonne ortonormali
    Implementazione: Modified Gram-Schmidt (MGS)
    """
    n = A.cols
    
    # Creiamo una lista di vettori modificabili partendo dalle colonne di A
    # Usiamo una lista perché aggiorneremo questi vettori iterativamente
    V = [A.col(i) for i in range(n)]
    
    # Lista per i vettori finali ortonormali (Q)
    ortho_norm = []

    for i in range(n):
        # 1. Prendi il vettore corrente
        v = V[i]

        # 2. Normalizzazione IMMEDIATA
        # Calcoliamo la norma e dividiamo. Ora abbiamo il vettore q_i unitario.
        norm_val = sqrt(v.dot(v))
        q = v / norm_val
        ortho_norm.append(q)

        # 3. Aggiornamento dei vettori FUTURI
        # Invece di guardare indietro, proiettiamo il nuovo q su tutti i 
        # vettori rimanenti e sottraiamo la componente subito.
        for j in range(i + 1, n):
            # Proiezione: (V[j] . q) * q
            # Nota: non serve dividere per q.dot(q) perché q è già normalizzato (è 1)
            V[j] = V[j] - V[j].dot(q) * q

    # Ricostruzione matrice finale
    Q = Matrix.hstack(*ortho_norm)
    return simplify(Q)

Bibliografia

Gram-Schmidt Orthogonalization and Projection, on Python for Linear Algebra
Matrices and Linear Algebra, Maxima 5.17.1 Manual

❮❮ ❯❯