Queste note sono basate sugli appunti fatti con Gianmarco Manzini negli anni 1995-2005

Metodo di eliminazione di Gauss

Appunti di calcolo numerico

Autore/Autrice

Affiliazione

Enrico Bertolazzi

University of Trento, Department of Industrial Engineering

Introduzione

Il metodo di eliminazione di Gauss¹ è una procedura sistematica per risolvere sistemi di equazioni lineari. Si basa sulla trasformazione del sistema originale in uno equivalente che abbia una struttura più semplice, detta forma triangolare, in cui tutte le equazioni successive contengono sempre meno incognite.

L’obiettivo principale del metodo è annullare gradualmente i coefficienti al di sotto della diagonale principale della matrice associata al sistema, tramite operazioni elementari sulle righe (somma, sottrazione, moltiplicazione per uno scalare). Una volta ottenuta la forma triangolare, le equazioni risultanti possono essere risolte facilmente mediante il processo di sostituzioni all’indietro, partendo dall’ultima equazione, che contiene una sola incognita, e risalendo fino alla prima.

Il metodo di Gauss è particolarmente efficiente e costituisce la base per molte altre tecniche numeriche utilizzate per la risoluzione di sistemi lineari.

Il metodo di eliminazione di Gauss, applicato ad un sistema lineare della forma:

\bgroup \left\{\begin{aligned} A_{1,1}x_1 + A_{1,2}x_2 + \cdots + A_{1,n}x_n &= b_1 \\ A_{2,1}x_1 + A_{2,2}x_2 + \cdots + A_{2,n}x_n &= b_2 \\ &\vdots\\ A_{n,1}x_1 + A_{n,2}x_2 + \cdots + A_{n,n}x_n &= b_n \end{aligned}\right.\egroup

lo trasforma in un sistema equivalente (cioè con la stessa soluzione) della forma:

\bgroup \left\{\begin{aligned} C_{1,1}x_1 + C_{1,2}x_2 + \cdots + C_{1,n}x_n &= d_1 \\ C_{2,2}x_2 + \cdots + C_{2,n}x_n &= d_2 \\ &\vdots\\ C_{n,n}x_n &= d_n \end{aligned}\right.\egroup \tag{1}

detta forma triangolare.

Un sistema in forma triangolare si può risolvere immediatamente tramite sostituzioni all’indietro, cioè partendo dall’ultima incognita, il cui valore, si osservi, è già noto. Il processo è il seguente:

\begin{aligned} x_n &= \dfrac{d_n}{C_{n,n}} \\ x_{n-1} &= \dfrac{d_{n-1} - C_{n-1,n}x_n}{C_{n-1,n-1}} \\ x_{n-2} &= \dfrac{d_{n-2} - C_{n-2,n}x_n - C_{n-2,n-1}x_{n-1}}{C_{n-2,n-2}} \\ &\vdots \\ x_1 &= \dfrac{d_1 - C_{1,n}x_n - C_{1,n-1}x_{n-1} - \cdots - C_{1,2}x_2}{C_{1,1}} \\ \end{aligned} \tag{2}

Esempio

La soluzione del sistema in forma triangolare:

\bgroup \left\{\begin{aligned} x+y+z &= 1 \\ 4y+2z &= 0 \\ z &= 2 \end{aligned}\right.\egroup

(nel quale abbiamo posto x=x_1, y=x_2 e z=x_3), si ottiene immediatamente sostituendo all’indietro le variabili come mostrato in Equazione 2,

\begin{aligned} z &= \dfrac{2}{1} = 2, \\ y &= \dfrac{0-2z}{4} = \dfrac{0-2\cdot 2}{4}=-1,\\ x &= \dfrac{1-1z-1y}{1} = \dfrac{1-1(-1)-1\cdot 2}{1} = 0. \end{aligned}

L’algoritmo di eliminazione di Gauss procede formalmente esprimendo la variabile che si vuole eliminare come combinazione lineare delle variabili non ancora eliminate. A tal scopo, si utilizza una equazione del sistema, che successivamente non sarà più presa in considerazione. L’eliminazione ha luogo per sostituzione della variabile da eliminare nelle rimanenti equazioni. L’algoritmo può essere espresso efficacemente, come vedremo, in termini matriciali poiché, nella pratica, si procede costruendo opportune combinazioni lineari delle equazioni del sistema, quindi di righe della matrice dei coefficienti.

Per fissare le idee, illustriamo il procedimento con un esempio dettagliato.

Esempio

Sia dato il seguente sistema lineare

\bgroup \left\{\begin{aligned} x+y+z+w &= 1 \qquad\qquad & [1] \\ -x-y+4z &= 0 & [2]\\ x-z-w &= 2 & [3]\\ x+w &= 0 & [4] \end{aligned}\right.\egroup

dove abbiamo numerato le equazioni a destra in parentesi quadre per facilità di riferimento. Trasformeremo questo sistema in forma triangolare eliminando, passo dopo passo, le incognite x, y, z², utilizzando combinazioni lineari delle equazioni.

Per eliminare l’incognita x dalle equazioni [2], [3] e [4], eseguiremo le seguenti operazioni:

• Eliminare x da [2]: Sommiamo l’equazione [1] all’equazione [2] per ottenere la nuova equazione [2']: [2'] \gets [2] + [1]

• Eliminare x da [3]: Sottraiamo l’equazione [1] dall’equazione [3] per ottenere la nuova equazione [3']: [3'] \gets [3] - [1]

• Eliminare x da [4]: Sottraiamo l’equazione [1] dall’equazione [4] per ottenere la nuova equazione [4']: [4'] \gets [4] - [1]

Si ottiene:

\bgroup \left\{\begin{aligned} x+y+z+w &= 1 \qquad\qquad & [1] \\ 5z+w &= 1 & [2'] \\ -y-2z-2w &= 1 & [3'] \\ -y-z &= -1 & [4'] \end{aligned}\right.\egroup

Osserviamo che l’equazione [2'] non contiene la variabile y, mentre l’equazione [3'] la contiene. Per semplificare il sistema, scambiamo l’equazione [2'] con la [3'] prima di procedere all’eliminazione successiva:

\bgroup \left\{\begin{aligned} x+y+z+w &= 1 \qquad\qquad & [1] \\ -y-2z-2w &= 1 & [3']\\ 5z+w &= 1 & [2']\\ -y-z &= -1 & [4'] \end{aligned}\right.\egroup

Per eliminare l’incognita y dall’equazione [4'], eseguiremo la seguente operazione:

• Eliminare y da [4']: Sommiamo l’equazione [3'] all’equazione [4'] per ottenere la nuova equazione [4'']: [4''] \gets [4'] - [3']

Si ottiene:

\bgroup \left\{\begin{aligned} x+y+z+w &= 1 \qquad\qquad & [1] \\ -y-2z-2w &= 1 & [3']\\ 5z+w &= 1 & [2']\\ z+2w &= -2 & [4''] \end{aligned}\right.\egroup

Per eliminare l’incognita z dall’equazione [4''], eseguiremo la seguente operazione:

• Eliminare z da [4'']: Sottraiamo \frac{1}{5} volte l’equazione [2'] dall’equazione [4''] per ottenere la nuova equazione [4''']: [4'''] \gets [4''] - \frac{1}{5}[2']

Si ottiene:

\bgroup \left\{\begin{aligned} x+y+z+w &=1 \qquad\qquad & [1] \\ -y-2z-2w &=1 & [3']\\ 5z+w &=1 & [2']\\ \frac{9}{5}w &=-\frac{11}{5} & [4'''] \end{aligned}\right.\egroup

Il sistema è ora in forma triangolare e, tramite il procedimento di sostituzioni all’indietro Equazione 2, possiamo calcolare la soluzione:

\begin{aligned} w &= \frac{-\frac{11}{5}}{\frac{9}{5}} = -\frac{11}{9} \\ z &= \frac{1 - \left(-\frac{11}{9}\right)}{5} = \frac{4}{9} \\ y &= \frac{1 + 2\left(-\frac{11}{9}\right) + 2 \frac{4}{9}}{-1} = \frac{5}{9} \\ x &= \frac{1 - \left(-\frac{11}{9}\right) - \frac{4}{9} - \frac{5}{9}}{1} = \frac{11}{9} \end{aligned}

Osservazione

Per trasformare un sistema lineare da una forma generica a una forma triangolare, sono sufficienti solo due tipi di operazioni:

1. Somma di un’equazione con un’altra moltiplicata per un opportuno scalare: Questa operazione consente di eliminare progressivamente le incognite da una o più equazioni, portando il sistema verso una forma triangolare.

2. Scambio di due equazioni: Questa operazione è utile per posizionare le equazioni in modo che il sistema possa essere facilmente risolto tramite sostituzioni all’indietro.

Si noti che la soluzione di un sistema non cambia se:

• Si scambiano due equazioni tra loro.
• Si sostituisce un’equazione con la stessa equazione sommata ad un’altra equazione moltiplicata per un opportuno scalare.

Queste operazioni mantengono l’equivalenza del sistema, preservando le soluzioni originali e facilitando il processo di riduzione a forma triangolare.

L’algoritmo di Gauss si può quindi descrivere come segue:

Metodo di Eliminazione di Gauss

Consideriamo un sistema lineare:

\boldsymbol{A}\boldsymbol{x}= \boldsymbol{b}

dove

\boldsymbol{A}= \begin{pmatrix} A_{1,1} & A_{1,2} & \cdots & A_{1,n} \\ A_{2,1} & A_{2,2} & \ddots & \vdots \\ \vdots & \ddots & \ddots & A_{n-1,n} \\ A_{n,1} & \cdots & A_{n,n-1} & A_{n,n} \end{pmatrix},

\boldsymbol{x}= \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix}, \qquad \boldsymbol{b}= \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{pmatrix}.

Algoritmo

1. Eliminazione Progressiva:

   • For i from 1 to n-1:
     • Trova il primo indice k \geq i tale che A_{k,i} \neq 0.
     • If k \neq i then
       • Scambia la i-esima equazione con la k-esima.
     • End if
     • For k from i+1 to n:
       • Calcola \beta \gets \frac{A_{k,i}}{A_{i,i}}.
       • Aggiorna b_k \gets b_k - \beta b_i.
       • For j from i+1 to n:
         • Aggiorna A_{k,j} \gets A_{k,j} - \beta A_{i,j}.
       • End for
     • End for
   • End for
   • Alla fine di queste operazioni, la matrice \boldsymbol{A} sarà in forma triangolare superiore:

   \boldsymbol{A}= \begin{pmatrix} A_{1,1} & A_{1,2} & \cdots & A_{1,n} \\ 0 & A_{2,2} & \cdots & A_{2,n} \\ 0 & 0 & \ddots & \vdots \\ \vdots & \vdots & & A_{n-1,n} \\ 0 & 0 & \cdots & A_{n,n} \end{pmatrix}, \qquad \boldsymbol{b}= \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_{n-1} \\ b_n \end{pmatrix}.

2. Sostituzione all’Indietro:

   • Calcola x_n \gets \frac{b_n}{A_{n,n}}.
   • For k from n-1 down to 1:
     • Calcola \beta\gets b_k.
     • For j from k+1 to n:
       • Aggiorna \beta \gets \beta - A_{k,j} x_j.
     • End for
     • Calcola x_k = \frac{\beta}{A_{k,k}}.
   • End for

3. Restituisci la Soluzione:

   • Restituisci il vettore \boldsymbol{x} contenente le soluzioni x_1, x_2, \ldots, x_n.

Osservazione

Nel metodo di eliminazione di Gauss, l’operazione di scambio di equazioni è fondamentale per garantire che la matrice venga trasformata correttamente in una forma triangolare superiore. L’algoritmo standard prevede di scambiare le righe per trovare un elemento diverso da zero nella colonna corrente, ma esiste un criterio alternativo chiamato pivoting che può migliorare la stabilità numerica del metodo.

Pivoting

Il criterio di pivoting consiste nel selezionare, tra tutte le righe dalla riga corrente in poi, quella con il valore assoluto massimo nella colonna corrente. Precisamente, si procede come segue:

• Sia k l’indice tale che \left\lvert A^{(i)}_{k,i}\right\rvert \geq \left\lvert A^{(i)}_{j,i}\right\rvert, \quad \forall j \geq i. In altre parole, k è l’indice della riga con il valore assoluto più grande nella colonna i, a partire dalla riga i.

• If k \neq i then

  • Scambia la i-esima riga con la k-esima riga nella matrice \boldsymbol{A}^{(i)}.
  • Scambia anche la componente i-esima con la componente k-esima del vettore \boldsymbol{b}^{(i)}.

Motivazione

Il pivoting ha lo scopo di migliorare l’accuratezza della soluzione calcolata. Questo è particolarmente importante quando si lavora con numeri molto piccoli. Infatti, se si divide per un numero molto vicino a zero, si rischia di ottenere risultati numerici imprecisi a causa dell’amplificazione degli errori di arrotondamento. Scegliendo l’elemento massimo nella colonna come pivot, si riduce il rischio di divisioni per numeri molto piccoli e si migliora la stabilità numerica dell’algoritmo.

In sintesi, il pivoting è una tecnica che aiuta a prevenire problemi numerici e a ottenere soluzioni più accurate per i sistemi di equazioni lineari.

Forma Matriciale del Metodo di Gauss

Il metodo di Gauss per la risoluzione di sistemi lineari può essere interpretato usando il formalismo delle matrici. Consideriamo un sistema lineare nella forma matriciale:

\boldsymbol{A}\boldsymbol{x}= \boldsymbol{b}

dove \boldsymbol{A} è una matrice quadrata n \times n, \boldsymbol{x} è un vettore colonna di incognite, e \boldsymbol{b} è un vettore colonna di termini noti.

Trasformazione in Forma Triangolare Superiore

L’obiettivo del metodo di Gauss è trasformare il sistema \boldsymbol{A}\boldsymbol{x}= \boldsymbol{b} in un sistema equivalente della forma:

\boldsymbol{U}\boldsymbol{x}= \boldsymbol{c}

dove \boldsymbol{U} è una matrice triangolare superiore e \boldsymbol{c} è un vettore modificato di termini noti. La matrice \boldsymbol{U} è ottenuta applicando una serie di operazioni elementari al sistema originale.

Passaggi per la Trasformazione

La trasformazione dal sistema originale \boldsymbol{A}\boldsymbol{x}= \boldsymbol{b} al sistema triangolare superiore \boldsymbol{U}\boldsymbol{x}= \boldsymbol{c} avviene attraverso i seguenti passaggi:

1. Operazioni Elementari: Utilizziamo matrici elementari, chiamate matrici di Frobenius, per applicare operazioni elementari alle righe della matrice \boldsymbol{A}. Queste operazioni nel caso del metodo di Gauss si riducono a:

   • Addizione di un multiplo di una riga ad un’altra riga.

2. Matrici di Scambio: Se necessario, scambiamo due righe della matrice \boldsymbol{A} utilizzando matrici di permutazione, che sono matrici identità con due righe scambiate. Queste matrici assicurano che l’algoritmo non fallisca a causa di zeri nei pivots (elementi diagonali).

Matrici di Frobenius³

Le matrici di Frobenius sono matrici elementari utilizzate per eseguire operazioni sui sistemi lineari. Ogni matrice di Frobenius rappresenta una trasformazione elementare specifica:

1. Addizione di Multipli di Righe ad Altre Righe: Una matrice di Frobenius per aggiungere multipli \beta_{k+1}, \beta_{k+2}, \ldots, \beta_{n} della riga k-esima alle righe k+1, k+2 fino a n è una matrice identità con gli \beta_j posizionati nella k-esima colonna. Per scriverla esplicitamente considerimao il vettore \boldsymbol{v}=\begin{pmatrix}0,\ldots,0,\beta_{k+1},\ldots,\beta_n\end{pmatrix}^{T} cioè tale che v_i = 0, \qquad i=1,2,\ldots,k in questo modo la matrice di Frobenius \boldsymbol{F} diventa \boldsymbol{F}= \boldsymbol{I}+ \boldsymbol{v}\boldsymbol{e}_k^T, o più esplicamente: \boldsymbol{F}= \left( \begin{array}{ccc|c|ccc} 1 & & & 0 & 0 & \cdots & 0 \\ & \ddots & & \vdots & \vdots & & \vdots \\ & & 1 & 0 & 0 & \cdots & 0 \\ \hline 0 & \cdots & 0 & 1 & 0 & \cdots & 0 \\ \hline 0 & \cdots & 0 & \beta_{k+1} & 1 & & \\ \vdots & & \vdots & \vdots & & \ddots & \\ 0 & \cdots & 0 & \beta_n & & & 1 \end{array} \right)

2. Scambio di Due Righe: Una matrice di permutazione per scambiare la i-esima e la j-esima riga (denotata con \boldsymbol{S}_{i,j}) è una matrice identità con le righe i e j scambiate. Esplicitamente prende la forma: \boldsymbol{S}_{i,j} = \left(\begin{array}{ccc|c|ccc|c|ccc} 1 & & & & & & & & & & \\ & \ddots & & & & & & & & & \\ & & 1 & & & & & & & & \\ \hline & & & 0 & & & & 1 & & & \\ \hline & & & & 1 & & & & & & \\ & & & & & \ddots & & & & & \\ & & & & & & 1 & & & & \\ & & & 1 & & & & 0 & & & \\ \hline & & & & & & & & 1 & & \\ \hline & & & & & & & & & \ddots & \\ & & & & & & & & & & 1 \end{array}\right) \tag{3}

   La stessa matrice ha una rappresentazione piu compatta come segue: \boldsymbol{S}_{i,j} = \boldsymbol{I}- (\boldsymbol{e}_i-\boldsymbol{e}_j)(\boldsymbol{e}_i-\boldsymbol{e}_j)^{T}

Proprietà delle matrici di Frobenius

Consideriamo la matrice di Frobenius \boldsymbol{F} definita in Equazione 3 e la matrice \boldsymbol{A},

\boldsymbol{A}= \begin{pmatrix} A_{1,1} & A_{1,2} & \cdots & A_{1,n} \\ \vdots & \vdots & & \vdots \\ A_{k,1} & A_{k,2} & \cdots & A_{k,n} \\ A_{k+1,1} & A_{k+1,2} & \cdots & A_{k+1,n} \\ \vdots & \vdots & & \vdots \\ A_{n,1} & A_{n,2} & \cdots & A_{n,n} \end{pmatrix},

il prodotto \boldsymbol{F}\boldsymbol{A} è

\boldsymbol{F}\boldsymbol{A}= \begin{pmatrix} A_{1,1} & A_{1,2} & \cdots & A_{1,n} \\ \vdots & \vdots & & \vdots \\ A_{k,1} & A_{k,2} & \cdots & A_{k,n} \\ A_{k+1,1}+\beta_{k+1}A_{k,1} & A_{k+1,2}+\beta_{k+1}A_{k,2} & \cdots & A_{k+1,n}+\beta_{k+1}A_{k,n} \\ \vdots & \vdots & & \vdots \\ A_{n,1}+\beta_n A_{k,1} & A_{n,2}+\beta_n A_{k,2} & \cdots & A_{n,n}+\beta_n A_{k,n} \end{pmatrix}.

L’effetto della moltiplicazione di \boldsymbol{F} sulla matrice \boldsymbol{A} può essere descritto come segue:

• La riga i-esima della matrice \boldsymbol{A} viene sostituita dalla somma della riga i-esima originale con un multiplo della riga k-esima, dove il multiplo è dato dallo scalare v_i.

In altre parole, stiamo effettuando le seguenti modifiche alle righe della matrice \boldsymbol{A}:

• Sostituiamo la riga (k+1)-esima con la riga

  A_{k+1,\bullet} \gets A_{k+1,\bullet} + \beta_{k+1} \cdot A_{k,\bullet}

• Sostituiamo la riga (k+2)-esima con la riga

  A_{k+2,\bullet} \gets A_{k+2,\bullet} + \beta_{k+2} \cdot A_{k,\bullet}

• E così via fino a:

  A_{n,\bullet} \gets A_{n,\bullet} + \beta_{n} \cdot A_{k,\bullet}

Lemma 1 (Prodotto di Matrici di Frobenius) Siano \boldsymbol{F}_i = \boldsymbol{I}+ \boldsymbol{v}_i \boldsymbol{e}_i^T le matrici di Frobenius, dove \boldsymbol{v}_i è un vettore con zeri nelle prime i componenti, ossia della forma:

\boldsymbol{v}_i = \begin{pmatrix} 0 \\ \vdots \\ 0 \\ \beta_{i+1} \\ \beta_{i+2} \\ \vdots \\ \beta_n \end{pmatrix}

Allora, per ogni 2 \leq k \leq n-1, si ha che:

\boldsymbol{F}_1 \boldsymbol{F}_2 \cdots \boldsymbol{F}_k = \boldsymbol{I}+ \boldsymbol{v}_1 \boldsymbol{e}_1^T + \boldsymbol{v}_2 \boldsymbol{e}_2^T + \cdots + \boldsymbol{v}_k \boldsymbol{e}_k^T

Dimostrazione. La dimostrazione procede per induzione sull’indice k.

Base dell’Induzione (k=2):

Calcoliamo il prodotto delle due matrici:

\begin{aligned} \boldsymbol{F}_1 \boldsymbol{F}_2 &= (\boldsymbol{I}+ \boldsymbol{v}_1 \boldsymbol{e}_1^T)(\boldsymbol{I}+ \boldsymbol{v}_2 \boldsymbol{e}_2^T) \\ &= \boldsymbol{I}+ \boldsymbol{v}_1 \boldsymbol{e}_1^T + \boldsymbol{v}_2 \boldsymbol{e}_2^T + \boldsymbol{v}_1 (\boldsymbol{e}_1^T \boldsymbol{v}_2) \boldsymbol{e}_2^T \\ &= \boldsymbol{I}+ \boldsymbol{v}_1 \boldsymbol{e}_1^T + \boldsymbol{v}_2 \boldsymbol{e}_2^T, \end{aligned}

poiché \boldsymbol{e}_1^T \boldsymbol{v}_2 = (\boldsymbol{v}_2)_1 = 0.

Passo Induttivo:

Supponiamo che il lemma sia vero per un indice generico k-1, cioè:

\boldsymbol{F}_1 \boldsymbol{F}_2 \cdots \boldsymbol{F}_{k-1} = \boldsymbol{I}+ \boldsymbol{v}_1 \boldsymbol{e}_1^T + \boldsymbol{v}_2 \boldsymbol{e}_2^T + \cdots + \boldsymbol{v}_{k-1} \boldsymbol{e}_{k-1}^T.

Ora, calcoliamo:

\begin{aligned} \left(\boldsymbol{F}_1 \boldsymbol{F}_2 \cdots \boldsymbol{F}_{k-1}\right) \boldsymbol{F}_k &= \left(\boldsymbol{I}+ \boldsymbol{v}_1 \boldsymbol{e}_1^T + \boldsymbol{v}_2 \boldsymbol{e}_2^T + \cdots + \boldsymbol{v}_{k-1} \boldsymbol{e}_{k-1}^T\right) \left(\boldsymbol{I}+ \boldsymbol{v}_k \boldsymbol{e}_k^T\right) \\ &= \boldsymbol{I}+ \boldsymbol{v}_1 \boldsymbol{e}_1^T + \boldsymbol{v}_2 \boldsymbol{e}_2^T + \cdots + \boldsymbol{v}_{k-1} \boldsymbol{e}_{k-1}^T + \boldsymbol{v}_k \boldsymbol{e}_k^T \\ &\quad + \boldsymbol{v}_1 (\boldsymbol{e}_1^T \boldsymbol{v}_k) \boldsymbol{e}_k^T + \cdots + \boldsymbol{v}_{k-1} (\boldsymbol{e}_{k-1}^T \boldsymbol{v}_k) \boldsymbol{e}_k^T \\ &= \boldsymbol{I}+ \boldsymbol{v}_1 \boldsymbol{e}_1^T + \boldsymbol{v}_2 \boldsymbol{e}_2^T + \cdots + \boldsymbol{v}_{k-1} \boldsymbol{e}_{k-1}^T + \boldsymbol{v}_k \boldsymbol{e}_k^T, \end{aligned}

poiché tutti i termini misti sono nulli, essendo

\boldsymbol{e}_i^T \boldsymbol{v}_k = (\boldsymbol{v}_k)_i = 0, \qquad i = 1, 2, \ldots, k-1.

Lemma 2 (Inversa della Matrice di Frobenius) Consideriamo la matrice di Frobenius \boldsymbol{F} definita come:

\boldsymbol{F}= \boldsymbol{I}+ \boldsymbol{v}\boldsymbol{e}_k^T

dove \boldsymbol{I} è la matrice identità, \boldsymbol{v} è un vettore colonna e \boldsymbol{e}_k^T è il vettore riga unitario con 1 nella k-esima posizione e zeri altrove.

Quando v_k = 0, l’inversa di \boldsymbol{F} è data dalla matrice:

\boldsymbol{F}^{-1} = \boldsymbol{I}- \boldsymbol{v}\boldsymbol{e}_k^T.

Dimostrazione. Per verificare che \boldsymbol{F}^{-1} = \boldsymbol{I}- \boldsymbol{v}\boldsymbol{e}_k^T sia effettivamente l’inversa di \boldsymbol{F}, dobbiamo mostrare che:

\boldsymbol{F}\cdot \boldsymbol{F}^{-1} = \boldsymbol{I}

Calcoliamo il prodotto:

\begin{aligned} \left(\boldsymbol{I}+ \boldsymbol{v}\boldsymbol{e}_k^T\right) \left(\boldsymbol{I}- \boldsymbol{v}\boldsymbol{e}_k^T\right) &= \boldsymbol{I}\cdot \boldsymbol{I}- \boldsymbol{I}\cdot \boldsymbol{v}\boldsymbol{e}_k^T + \boldsymbol{v}\boldsymbol{e}_k^T \cdot \boldsymbol{I}- \boldsymbol{v}\boldsymbol{e}_k^T \cdot \boldsymbol{v}\boldsymbol{e}_k^T \\ &= \boldsymbol{I}- \boldsymbol{v}\boldsymbol{e}_k^T + \boldsymbol{v}\boldsymbol{e}_k^T - \boldsymbol{v}(\boldsymbol{e}_k^T \boldsymbol{v}) \boldsymbol{e}_k^T \\ &= \boldsymbol{I}- \boldsymbol{v}(\boldsymbol{e}_k^T \boldsymbol{v}) \boldsymbol{e}_k^T \end{aligned}

Notiamo che i termini \boldsymbol{v}\boldsymbol{e}_k^T si cancellano:

-\boldsymbol{v}\boldsymbol{e}_k^T + \boldsymbol{v}\boldsymbol{e}_k^T = 0

Pertanto, rimane:

\left(\boldsymbol{I}+ \boldsymbol{v}\boldsymbol{e}_k^T\right) \left(\boldsymbol{I}- \boldsymbol{v}\boldsymbol{e}_k^T\right) = \boldsymbol{I}- \boldsymbol{v}(\boldsymbol{e}_k^T \boldsymbol{v}) \boldsymbol{e}_k^T

Per dimostrare che \boldsymbol{F}^{-1} = \boldsymbol{I}- \boldsymbol{v}\boldsymbol{e}_k^T è effettivamente l’inversa di \boldsymbol{F}, dobbiamo verificare che il termine \boldsymbol{v}(\boldsymbol{e}_k^T \boldsymbol{v}) \boldsymbol{e}_k^T sia nullo. Poiché \boldsymbol{e}_k^T \boldsymbol{v} è il k-esimo elemento del vettore \boldsymbol{v}, denotato come v_k, abbiamo:

\boldsymbol{e}_k^T \boldsymbol{v}= v_k

Se v_k = 0 (cioè l’k-esimo elemento di \boldsymbol{v} è zero), allora:

\left(\boldsymbol{I}+ \boldsymbol{v}\boldsymbol{e}_k^T\right) \left(\boldsymbol{I}- \boldsymbol{v}\boldsymbol{e}_k^T\right) = \boldsymbol{I}

Poiché v_k = 0 allora \boldsymbol{F}^{-1} = \boldsymbol{I}- \boldsymbol{v}\boldsymbol{e}_k^T è l’inversa della matrice di Frobenius.

Osservazione

• Data la matrice \boldsymbol{A} dove \boldsymbol{A}_{i,\bullet} è la riga i-esima della matrice \boldsymbol{A}, l’effetto di \boldsymbol{S}_{i,j}\boldsymbol{A} è il seguente:

  \boldsymbol{A}= \begin{pmatrix} \vdots \\ \boldsymbol{A}_{i,\bullet} \\ \vdots \\ \boldsymbol{A}_{j,\bullet} \\ \vdots \end{pmatrix}, \qquad \boldsymbol{S}\boldsymbol{A}= \begin{pmatrix} \vdots \\ \boldsymbol{A}_{j,\bullet} \\ \vdots \\ \boldsymbol{A}_{i,\bullet} \\ \vdots \end{pmatrix},

  cioè scambia la riga i-esima con la riga j-esima.

• Scambiando due volte le stesse righe di una matrice riotteniamo la matrice iniziale. Quindi, se \boldsymbol{S}_{i,j} è una matrice di scambio allora \boldsymbol{S}_{i,j}\boldsymbol{S}_{i,j}=\boldsymbol{I} cioè \boldsymbol{S}_{i,j} coincide con la sua inversa, cioè è una matrice ortoginale.

• Il prodotto di due matrici di scambio è una matrice di permutazione, che modifica l’ordinamento delle righe di una matrice.

• Le matrici di permutazione si possono scrivere come il prodotto di matrici di scambio. In generale il prodotto di matrici ortogonali è una matrice ortogonale⁴, e quindi l’inversa della matrice di permutazione coincide con la sua trasposta.

Primo passo del metodo di Gauss

Dato il sistema \boldsymbol{A}\boldsymbol{x}=\boldsymbol{b} poniamo \boldsymbol{A}^{(1)}=\boldsymbol{A} e \boldsymbol{b}^{(1)}=\boldsymbol{b},

\boldsymbol{A}^{(1)}=\begin{pmatrix} A^{(1)}_{1,1} & A^{(1)}_{1,2} & \cdots & A^{(1)}_{1,n} \\ A^{(1)}_{2,1} & A^{(1)}_{2,2} & \ddots & \vdots \\ \vdots & \ddots & \ddots & A^{(1)}_{n-1,n} \\ A^{(1)}_{n,1} & \cdots & A^{(1)}_{n,n-1} & A^{(1)}_{n,n} \end{pmatrix}.

Supponendo A^{(1)}_{1,1}\neq 0 definiamo

\begin{aligned} v^{(1)}_i &= \begin{cases} 0 & i=1 \\ \dfrac{A^{(1)}_{i,1}} {A^{(1)}_{1,1}} & i=2,3,\ldots,n \end{cases}, \\ %\vv^{(1)} &= \MAT{ % 0 \\ v^{(1)}_2 \\ v^{(1)}_3 \\ \vdots \\ v^{(1)}_n %}, %\\ \boldsymbol{L}_1 &= \boldsymbol{I}-\boldsymbol{v}^{(1)}\boldsymbol{e}_1^T = \begin{pmatrix} 1 & 0 & \cdots & & 0 \\ -v^{(1)}_2 & 1 & & & \vdots \\ \vdots & & \ddots & & 0 \\ -v^{(1)}_{n-1} & & & 1 & 0 \\ -v^{(1)}_n & 0 & \cdots & 0 & 1 \end{pmatrix}. \end{aligned}

Moltiplicando il sistema \boldsymbol{A}^{(1)}\boldsymbol{x}=\boldsymbol{b}^{(1)} a sinistra per la matrice di Frobenius \boldsymbol{L}_1 otteniamo il sistema trasformato \boldsymbol{A}^{(2)}\boldsymbol{x}= \boldsymbol{b}^{(2)} dove,

\boldsymbol{A}^{(2)} = \boldsymbol{L}_1\boldsymbol{A}^{(1)} = \left(\begin{array}{c|cccc} A^{(1)}_{1,1} & A^{(1)}_{1,2} & A^{(1)}_{1,3} & \cdots & A^{(1)}_{1,n} \\ \hline 0 & A^{(2)}_{2,2} & A^{(2)}_{2,3} & \cdots & A^{(2)}_{2,n} \\ 0 & A^{(2)}_{3,2} & A^{(2)}_{3,3} & \cdots & A^{(2)}_{3,n} \\ \vdots & \vdots & & & \vdots \\ 0 & A^{(2)}_{n,2} & A^{(2)}_{n,3} & \cdots & A^{(2)}_{n,n} \end{array}\right),

\boldsymbol{b}^{(2)} = \boldsymbol{L}_1\boldsymbol{b}^{(1)} = \left(\begin{array}{c} b^{(1)}_1 \\ b^{(2)}_2 \\ b^{(2)}_3 \\ \vdots \\ b^{(2)}_n \end{array}\right).

ed inoltre si ha

\begin{aligned} A^{(2)}_{i,j} &= A^{(1)}_{i,j}-v^{(1)}_i A^{(1)}_{1,j}, \qquad & i,j &= 2,3,\ldots,n \\ b^{(2)}_i &= b^{(1)}_i-v^{(1)}_i b^{(1)}_1 & i &= 2,3,\ldots,n. \end{aligned}

k-esimo passo del metodo di Gauss

Supponiamo di aver effettuato k-1 passi dell’algoritmo di Gauss (senza permutazione di righe). Siamo quindi nella situazione

\boldsymbol{A}^{(k)}=\boldsymbol{L}_{k-1}\cdots\boldsymbol{L}_2\boldsymbol{L}_1\boldsymbol{A}^{(1)}

e

\boldsymbol{b}^{(k)}=\boldsymbol{L}_{k-1}\cdots\boldsymbol{L}_2\boldsymbol{L}_1\boldsymbol{b}^{(1)}

dove

\boldsymbol{A}^{(k)}= \left(\begin{array}{ccc|ccc} A^{(1)}_{1,1} & \cdots & A^{(1)}_{1,k-1} & ^{(1)}_{1,k} & \cdots & A^{(1)}_{1,n} \\ & \ddots & \vdots & \vdots & & \vdots \\ 0 & & A^{(k-1)}_{k-1,k-1} & A^{(k-1)}_{k-1,k} & \cdots & A^{(k-1)}_{k-1,n} \\ \hline 0 & \cdots & 0 & A^{(k)}_{k,k} & \cdots & A^{(k)}_{k,n} \\ \vdots & & \vdots & \vdots & & \vdots \\ 0 & \cdots & 0 & A^{(k)}_{n,k} & \cdots & A^{(k)}_{n,n} \end{array}\right)

\boldsymbol{b}^{(k)} = \left(\begin{array}{c} b^{(1)}_1 \\ \vdots \\ b^{(k-1)}_{k-1} \\ b^{(k)}_k \\ \vdots \\ b^{(k)}_n \end{array}\right).

Supponendo A^{(k)}_{k,k}\neq 0 e ponendo

\begin{aligned} v^{(k)}_i &= \begin{cases} 0 & i=1,2,\ldots,k \\ \dfrac{A^{(k)}_{i,k}} {A^{(k)}_{k,k}} & i=k+1,k+2,\ldots,n \end{cases}, \\ % \vv^{(k)} = % \MAT{ 0 \\ \vdots \\ 0 \\ v^{(k)}_{k+1} \\ \vdots %\\ v^{(k)}_n }, \\ \boldsymbol{L}_k &= \boldsymbol{I}-\boldsymbol{v}_k\boldsymbol{e}_k^T = \left(\begin{array}{ccc|c|ccc} 1 & & & 0 & 0 & \cdots & 0 \\ & \ddots & & \vdots & \vdots & & \vdots \\ & & 1 & 0 & 0 & \cdots & 0 \\ \hline 0 & \cdots & 0 & 1 & 0 & \cdots & 0 \\ \hline 0 & \cdots & 0 & -v^{(k)}_{k+1} & 1 & & \\ \vdots & & \vdots & \vdots & & \ddots & \\ 0 & \cdots & 0 & -v^{(k)}_n & & & 1 \end{array}\right), \end{aligned}

otteniamo

\boldsymbol{L}_k\boldsymbol{A}^{(k)}= \left(\begin{array}{ccc|cccc} A^{(1)}_{1,1} & \cdots & A^{(1)}_{1,k-1} & A^{(1)}_{1,k} & A^{(1)}_{1,k+1} & \cdots & A^{(1)}_{1,n} \\ & \ddots & \vdots & \vdots & \vdots & & \vdots \\ 0 & & A^{(k-1)}_{k-1,k-1} & A^{(k-1)}_{k-1,k} & A^{(k-1)}_{k-1,k+1} & \cdots & A^{(k-1)}_{k-1,n} \\ \hline 0 & \cdots & 0 & A^{(k)}_{k,k} & A^{(k)}_{k,k+1} & \cdots & A^{(k)}_{k,n} \\ 0 & \cdots & 0 & 0 & A^{(k+1)}_{k+1,k+1} & \cdots & A^{(k+1)}_{k+1,n} \\ \vdots & & \vdots & \vdots & \vdots & & \vdots \\ 0 & \cdots & 0 & 0 & A^{(k+1)}_{n,k+1} & \cdots & A^{(k+1)}_{n,n} \end{array}\right),

\boldsymbol{L}_k\boldsymbol{b}^{(k)}=\left(\begin{array}{c} b^{(1)}_1 \\ \vdots \\ b^{(k-1)}_{k-1} \\ b^{(k)}_{k} \\ b^{(k+1)}_{k+1} \\ \vdots \\ b^{(k+1)}_n \end{array}\right),

dove

\begin{aligned} A^{(k+1)}_{i,j} &= A^{(k)}_{i,j}-v^{(k)}_i A^{(k)}_{k,j},\qquad &i,j&=k+1,k+2,\ldots,n \\ b^{(k+1)}_i &= b^{(k)}_i-v^{(k)}_i b^{(k)}_i & i&=k+1,k+2,\ldots,n. \end{aligned}

L’algoritmo di Gauss applica lo stesso processo iterativo su matrici sempre più piccole fino a ottenere la triangolarizzazione della matrice originale. Ad ogni passo, le operazioni elementari di riga riducono progressivamente la dimensione della matrice fino a raggiungere una forma triangolare superiore.

Osservazione

Consideriamo una matrice quadrata \boldsymbol{A} di dimensione n \times n e supponiamo che ogni passo dell’algoritmo di Gauss sia applicabile, cioè A^{(k)}_{k,k} \neq 0 per ogni k. Dopo n-1 passi dell’algoritmo, otteniamo:

\boldsymbol{L}_{n-1} \boldsymbol{L}_{n-2} \cdots \boldsymbol{L}_1 \boldsymbol{A}= \boldsymbol{U},

ponendo

\boldsymbol{L}= \boldsymbol{L}_1^{-1} \boldsymbol{L}_2^{-1} \cdots \boldsymbol{L}_{n-1}^{-1}.

possiamo scrivere

\boldsymbol{A}= \boldsymbol{L}\boldsymbol{U}.

Inoltre:

1. Le matrici \boldsymbol{L}_k = \boldsymbol{I}- \boldsymbol{v}_k \boldsymbol{e}_k^T sono matrici di Frobenius.

2. Le matrici \boldsymbol{L}_k^{-1} = \boldsymbol{I}+ \boldsymbol{v}_k \boldsymbol{e}_k^T sono ancora matrici di Frobenius.

3. Applicando il lemma Lemma 1, abbiamo:

\begin{aligned} \boldsymbol{L}&= \boldsymbol{L}_1^{-1} \boldsymbol{L}_2^{-1} \cdots \boldsymbol{L}_{n-1}^{-1}, \\ &= (\boldsymbol{I}+ \boldsymbol{v}_1 \boldsymbol{e}_1^T) (\boldsymbol{I}+ \boldsymbol{v}_2 \boldsymbol{e}_2^T) \cdots (\boldsymbol{I}+ \boldsymbol{v}_{n-1} \boldsymbol{e}_{n-1}^T), \\ &= \boldsymbol{I}+ \boldsymbol{v}_1 \boldsymbol{e}_1^T + \boldsymbol{v}_2 \boldsymbol{e}_2^T + \cdots + \boldsymbol{v}_{n-1} \boldsymbol{e}_{n-1}^T. \end{aligned}

Quindi, \boldsymbol{L} è una matrice triangolare inferiore.

La matrice \boldsymbol{A} è stata quindi decomposta come il prodotto di una matrice triangolare inferiore \boldsymbol{L} e una matrice triangolare superiore \boldsymbol{U}. Questa decomposizione è nota come decomposizione LU o fattorizzazione LU.

Osservazione

Dalla decomposizione LU della matrice \boldsymbol{A}, abbiamo:

\begin{aligned} \boldsymbol{L}&= \begin{pmatrix} 1 & 0 & 0 & \cdots & 0 \\ v^{(1)}_2 & 1 & 0 & \cdots & 0 \\ v^{(1)}_3 & v^{(2)}_3 & \ddots & \ddots & \vdots \\ \vdots & \vdots & & 1 & 0 \\ v^{(1)}_n & v^{(2)}_n & \cdots & v^{(n-1)}_n & 1 \end{pmatrix}, \quad \text{dove} \quad v^{(j)}_i = \frac{A^{(j)}_{i,j}}{A^{(j)}_{j,j}} \\ \boldsymbol{U}&= \begin{pmatrix} A^{(1)}_{1,1} & A^{(1)}_{1,2} & A^{(1)}_{1,3} & \cdots & A^{(1)}_{1,n} \\ 0 & A^{(2)}_{2,2} & A^{(1)}_{2,3} & \cdots & A^{(2)}_{2,n} \\ 0 & 0 & \ddots & \ddots & \vdots \\ \vdots & \vdots & & A^{(n-1)}_{n-1,n-1} & A^{(n-1)}_{n-1,n} \\ 0 & 0 & \cdots & 0 & A^{(n)}_{n,n} \end{pmatrix} \end{aligned}

Da quanto sopra, è evidente che le matrici \boldsymbol{L} e \boldsymbol{U} possono essere memorizzate all’interno della matrice originale \boldsymbol{A}, senza memorizzare esplicitamente la diagonale di \boldsymbol{L}. Questo approccio consente di implementare un algoritmo che calcola la decomposizione LU e salva direttamente il risultato nella stessa matrice \boldsymbol{A}.

Algoritmo Decomposizione LU

L’algoritmo modifica la matrice \boldsymbol{A} in-place per calcolare la decomposizione LU, dove \boldsymbol{L} è una matrice triangolare inferiore con 1 sulla diagonale e \boldsymbol{U} è una matrice triangolare superiore.

• INPUT: La matrice \boldsymbol{A} di dimensione n \times n.

• OUTPUT: La matrice \boldsymbol{A} modificata contenente i fattori della decomposizione LU.

• For i \gets 1 to n-1:

  • For k \gets i+1 to n:
    • Compute \beta \gets \frac{A_{k,i}}{A_{i,i}}.
    • For j \gets i+1 to n:
      • Update A_{k,j} \gets A_{k,j} - \beta A_{i,j}.
    • End For
    • Set A_{k,i} \gets \beta.
  • End For

• End For

Durante l’esecuzione, le righe della matrice \boldsymbol{A} vengono modificate per ottenere una matrice triangolare superiore \boldsymbol{U}, mentre i coefficienti di eliminazione vengono memorizzati nella parte inferiore della matrice \boldsymbol{A}, formando così la matrice triangolare inferiore \boldsymbol{L}.

Algoritmo di Gauss in presenza di pivoting

Nel processo di eliminazione di Gauss di base, l’algoritmo esegue operazioni sulle righe per convertire una matrice in una forma triangolare superiore. Tuttavia, questo metodo può incontrare instabilità numerica, specialmente quando la matrice presenta elementi di pivot molto piccoli o è mal condizionata. Per migliorare la stabilità numerica e l’accuratezza, si utilizza il pivoting parziale.

Il pivoting parziale consiste nel selezionare l’elemento di valore assoluto massimo nella colonna corrente come pivot e scambiare le righe per portare questo elemento massimo nella posizione del pivot. Questa strategia aiuta a minimizzare gli errori di arrotondamento e migliora la robustezza dell’algoritmo.

Per analizzare l’algoritmo di eliminazione di Gauss con pivoting parziale, è necessario gestire anche le matrici di scambio che vengono introdotte. Queste matrici sono utilizzate per scambiare righe e devono essere propagate attraverso le operazioni di fattorizzazione. Introduciamo quindi il seguente lemma per capire questo processo:

Lemma 3 (Frobenius e Scambi) Sia \boldsymbol{F}= \boldsymbol{I}- \boldsymbol{v}\boldsymbol{e}_k^T una matrice di Frobenius e \boldsymbol{S}_{i,j} una matrice di scambio che scambia le righe i e j. Se i e j sono entrambi maggiori di k, allora la matrice risultante \boldsymbol{F}\boldsymbol{S}_{i,j} può essere scritta come \boldsymbol{S}_{i,j} \boldsymbol{F}', dove \boldsymbol{F}' è anch’essa una matrice di Frobenius. In particolare:

\boldsymbol{F}' = \boldsymbol{I}- \boldsymbol{w}\boldsymbol{e}_k^T, \quad \text{dove} \quad \boldsymbol{w}= \boldsymbol{S}_{i,j} \boldsymbol{v}.

Dimostrazione. Calcoliamo direttamente il prodotto \boldsymbol{F}\boldsymbol{S}_{i,j} e confrontiamolo con la forma \boldsymbol{S}_{i,j} \boldsymbol{F}'.

1. Osservazione sui vettori unitari: Nota che \boldsymbol{e}_k^T \boldsymbol{S}_{i,j} \boldsymbol{x}= x_k=\boldsymbol{e}_k^T\boldsymbol{x} se i \neq k e j \neq k. Questo perché, quando i e j sono diversi da k, la matrice di scambio \boldsymbol{S}_{i,j} non sposta la k-esima componente. Quindi possiamo scrivere \boldsymbol{e}_k^T \boldsymbol{S}_{i,j} = \boldsymbol{e}_k^T.

2. Calcolo di \boldsymbol{F}\boldsymbol{S}_{i,j}: Consideriamo il prodotto: \begin{aligned} \boldsymbol{F}\boldsymbol{S}_{i,j} &= (\boldsymbol{I}- \boldsymbol{v}\boldsymbol{e}_k^T) \boldsymbol{S}_{i,j} \\ &= \boldsymbol{S}_{i,j} - \boldsymbol{v}\boldsymbol{e}_k^T \boldsymbol{S}_{i,j} \\ &= \boldsymbol{S}_{i,j} - \boldsymbol{v}\boldsymbol{e}_k^T\qquad[\text{per il punto 1}] \end{aligned}

3. Calcolo di \boldsymbol{S}_{i,j} \boldsymbol{F}': Sostituendo \boldsymbol{F}' = \boldsymbol{I}- \boldsymbol{w}\boldsymbol{e}_k^T, otteniamo:

   \begin{aligned} \boldsymbol{S}_{i,j} \boldsymbol{F}' &= \boldsymbol{S}_{i,j} (\boldsymbol{I}- \boldsymbol{w}\boldsymbol{e}_k^T) \\ &= \boldsymbol{S}_{i,j} - \boldsymbol{S}_{i,j} \boldsymbol{w}\boldsymbol{e}_k^T, \\ &= \boldsymbol{S}_{i,j} - \boldsymbol{S}_{i,j}\boldsymbol{S}_{i,j} \boldsymbol{v}\boldsymbol{e}_k^T .\\ &= \boldsymbol{S}_{i,j} - \boldsymbol{v}\boldsymbol{e}_k^T. \qquad[\text{per la proprietà $\boldsymbol{S}_{i,j}\boldsymbol{S}_{i,j}=\boldsymbol{I}$}] \end{aligned}

e quindi \boldsymbol{F}\boldsymbol{S}_{i,j}=\boldsymbol{S}_{i,j} \boldsymbol{F}'.

L’algoritmo di Gauss con pivoting può essere formulato come segue:

\boldsymbol{L}_{n-1}\boldsymbol{P}_{n-1}\boldsymbol{L}_{n-2}\boldsymbol{P}_{n-2}\cdots\boldsymbol{L}_2\boldsymbol{P}_2\boldsymbol{L}_1\boldsymbol{P}_1\boldsymbol{A}= \boldsymbol{U},

dove \boldsymbol{P}_k è una matrice di scambio, utilizzata se è necessario scambiare righe al k-esimo passo dell’algoritmo di Gauss. Quando non è necessario scambiare righe, \boldsymbol{P}_k è la matrice identità. Si osserva che \boldsymbol{P}_k = \boldsymbol{S}_{i,j} per i,j \geq k, e quindi, applicando il lemma Lemma 3, otteniamo che:

\boldsymbol{P}_k \boldsymbol{L}_{k-1} = \boldsymbol{L}_{k-1}' \boldsymbol{P}_k,

dove \boldsymbol{L}_{k-1}' = \boldsymbol{I}- \boldsymbol{w}_k \boldsymbol{e}_k^T e \boldsymbol{w}_k = \boldsymbol{P}_k \cdots \boldsymbol{P}_2 \boldsymbol{P}_1 \boldsymbol{v}_k. Spostando le matrici di scambio a destra nell’equazione precedente, otteniamo:

\left(\boldsymbol{L}_{n-1}' \boldsymbol{L}_{n-2}' \cdots \boldsymbol{L}_2' \boldsymbol{L}_1'\right) \left(\boldsymbol{P}_{n-1} \boldsymbol{P}_{n-2} \cdots \boldsymbol{P}_2 \boldsymbol{P}_1\right) \boldsymbol{A}= \boldsymbol{U},

dove la matrice \boldsymbol{U} è triangolare superiore e le matrici \boldsymbol{L}_k' sono matrici di Frobenius date da:

\boldsymbol{L}_k' = \boldsymbol{I}- \boldsymbol{w}_k \boldsymbol{e}_k^T, \qquad \boldsymbol{w}_k = \boldsymbol{P}_k \cdots \boldsymbol{P}_2 \boldsymbol{P}_1 \boldsymbol{v}_k.

Definiamo ora:

\boldsymbol{L}= (\boldsymbol{L}_1')^{-1} (\boldsymbol{L}_2')^{-1} \cdots (\boldsymbol{L}_{n-1}')^{-1}, \qquad \boldsymbol{P}= \boldsymbol{P}_{n-1} \boldsymbol{P}_{n-2} \cdots \boldsymbol{P}_2 \boldsymbol{P}_1,

e otteniamo:

\boldsymbol{P}\boldsymbol{A}= \boldsymbol{L}\boldsymbol{U}.

Quindi, l’algoritmo di eliminazione consente di decomporre la matrice \boldsymbol{P}\boldsymbol{A} in un prodotto \boldsymbol{L}\boldsymbol{U}, dove \boldsymbol{P} rappresenta la sequenza di scambi di righe applicati alla matrice \boldsymbol{A}, e \boldsymbol{L} e \boldsymbol{U} sono rispettivamente la matrice triangolare inferiore e superiore ottenute tramite la decomposizione LU.

Osservazione

Si è illustrato il funzionamento dell’algoritmo di fattorizzazione con pivoting, noto come pivot parziale, che prevede solo lo scambio delle righe. Questo tipo di pivoting equivale a riordinare le equazioni del sistema lineare.

Esiste però un pivoting più generale, chiamato pivot totale, che consente sia scambi di righe che di colonne. In questo caso, lo scambio di colonne corrisponde a una diversa numerazione delle variabili nel sistema lineare.

Il metodo di fattorizzazione LU con pivot parziale si estende facilmente al pivot totale, e si ottiene la seguente decomposizione:

\boldsymbol{P}\boldsymbol{A}\boldsymbol{Q}= \boldsymbol{L}\boldsymbol{U},

dove \boldsymbol{P} è la matrice di permutazione delle righe, ottenuta come prodotto delle matrici elementari di scambio delle righe, e \boldsymbol{Q} è la matrice di permutazione delle colonne, ottenuta come prodotto delle matrici elementari di scambio delle colonne.

Algoritmo: Decomposizione LU con Pivoting

• INPUT: La matrice \boldsymbol{A}

• OUTPUT: La matrice \boldsymbol{A} contenente la decomposizione LU e il vettore \boldsymbol{p} contenente la permutazione

1. Inizializza il vettore che memorizza le permutazioni
   • For i \gets 1 to n
     • Set p_i \gets i
   • End for
2. Algoritmo di fattorizzazione
   • For i \gets 1 to n-1
     • Set i_p \gets p_i
     • Pivoting: Trova k tale che \left\lvert A_{k_p,i}\right\rvert \geq \left\lvert A_{j,i}\right\rvert per j \geq i
     • Set k_p \gets i_p
     • For j \gets i+1 to n
       • Set j_p \gets p_j
       • If \left\lvert A_{j_p,i}\right\rvert > \left\lvert A_{k_p,i}\right\rvert then
         • Set k \gets j
         • Set k_p \gets j_p
       • End if
     • End for
     • Scambia la riga i con la riga k
       • Swap p_i con p_k
     • Elimination
     • For k \gets i+1 to n
       • Set k_p \gets p_k
         • Compute \beta \gets A_{k_p,i} / A_{i_p,i}
       • For j \gets i+1 to n
       • Update A_{k_p,j} \gets A_{k_p,j} - \beta A_{i_p,j}
       • End for
       • Set A_{k_p,i} \gets \beta
     • End for
   • End For
3. Return \boldsymbol{A}, \boldsymbol{p}

Osservazione: Pivot Numerico

Idealmente, se l’elemento pivotale è non nullo, non sarebbe necessario eseguire uno scambio di righe. Tuttavia, se l’elemento pivotale è molto piccolo in valore assoluto, l’algoritmo di fattorizzazione può produrre complementi di Schur molto grandi, il che potrebbe compromettere la stabilità numerica dell’algoritmo.

Per affrontare questo problema, nella pratica si adotta una strategia di controllo: si sceglie l’elemento pivotale con il valore assoluto maggiore nella colonna — noto come pivot parziale — o nella sottomatrice — noto come pivot totale — composta dagli elementi che devono ancora essere fattorizzati, quando l’elemento pivotale corrente è al di sotto di una soglia prestabilita.

Questa tecnica, chiamata pivot numerico, è progettata per garantire la stabilità numerica dell’algoritmo di fattorizzazione, prevenendo problemi che potrebbero derivare da piccoli valori degli elementi pivotali.

Osservazione: Pivot Simbolico

Oltre alla stabilità numerica, c’è un’altra considerazione importante nella scelta del pivot: la preservazione della sparsità della matrice.

Nella pratica, spesso si lavora con matrici sparse, cioè matrici in cui solo un numero limitato di elementi è diverso da zero, tipicamente dell’ordine della dimensione della matrice piuttosto che del suo quadrato. È cruciale mantenere questa sparsità durante il processo di fattorizzazione.

Il meccanismo del complemento di Schur tende infatti a generare molti elementi non nulli, il che può portare a una crescita esponenziale nel numero di elementi non nulli e, conseguentemente, a un eccessivo consumo di memoria per memorizzare le matrici \boldsymbol{L} e \boldsymbol{U} risultanti.

Per controllare la crescita di elementi non nulli e preservare la sparsità, è utile eseguire scambi di righe, che corrispondono a scambi delle equazioni nel sistema. Questo approccio può aiutare a mantenere la matrice più sparsa possibile.

Mostriamo questo effetto con il seguente esempio:

Esempio

Consideriamo la fattorizzazione LU della matrice \boldsymbol{A}\in\Bbb{R}^{10\times 10}

\boldsymbol{A}= \begin{pmatrix} \boldsymbol{1}& \boldsymbol{1}& \boldsymbol{1}& \boldsymbol{1}& \boldsymbol{1}& \boldsymbol{1}& \boldsymbol{1}& \boldsymbol{1}& \boldsymbol{1}& \boldsymbol{1}\\ \boldsymbol{1}& \boldsymbol{1}& 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \boldsymbol{1}& 0 & \boldsymbol{1}& 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \boldsymbol{1}& 0 & 0 & \boldsymbol{1}& 0 & 0 & 0 & 0 & 0 & 0 \\ \boldsymbol{1}& 0 & 0 & 0 & \boldsymbol{1}& 0 & 0 & 0 & 0 & 0 \\ \boldsymbol{1}& 0 & 0 & 0 & 0 & \boldsymbol{1}& 0 & 0 & 0 & 0 \\ \boldsymbol{1}& 0 & 0 & 0 & 0 & 0 & \boldsymbol{1}& 0 & 0 & 0 \\ \boldsymbol{1}& 0 & 0 & 0 & 0 & 0 & 0 & \boldsymbol{1}& 0 & 0 \\ \boldsymbol{1}& 0 & 0 & 0 & 0 & 0 & 0 & 0 & \boldsymbol{1}& 0 \\ \boldsymbol{1}& 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \boldsymbol{1} \end{pmatrix}

e della matrice che si ottiene prendendo righe e colonne nell’ordine inverso, cioè

\boldsymbol{P}\boldsymbol{A}\boldsymbol{P}= \begin{pmatrix} \boldsymbol{1}& 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \boldsymbol{1}\\ 0 & \boldsymbol{1}& 0 & 0 & 0 & 0 & 0 & 0 & 0 & \boldsymbol{1}\\ 0 & 0 & \boldsymbol{1}& 0 & 0 & 0 & 0 & 0 & 0 & \boldsymbol{1}\\ 0 & 0 & 0 & \boldsymbol{1}& 0 & 0 & 0 & 0 & 0 & \boldsymbol{1}\\ 0 & 0 & 0 & 0 & \boldsymbol{1}& 0 & 0 & 0 & 0 & \boldsymbol{1}\\ 0 & 0 & 0 & 0 & 0 & \boldsymbol{1}& 0 & 0 & 0 & \boldsymbol{1}\\ 0 & 0 & 0 & 0 & 0 & 0 & \boldsymbol{1}& 0 & 0 & \boldsymbol{1}\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \boldsymbol{1}& 0 & \boldsymbol{1}\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \boldsymbol{1}& \boldsymbol{1}\\ \boldsymbol{1}& \boldsymbol{1}& \boldsymbol{1}& \boldsymbol{1}& \boldsymbol{1}& \boldsymbol{1}& \boldsymbol{1}& \boldsymbol{1}& \boldsymbol{1}& \boldsymbol{1} \end{pmatrix} Abbiamo considerato gli elementi non nulli di \boldsymbol{A} come aventi valore 1, ma l’argomento che desideriamo illustrare riguarda esclusivamente il processo di formazione di nuovi elementi non nulli durante la fattorizzazione, noto con il termine inglese fill-in. Questo è distinto dal valore che tali elementi assumono nei fattori di Gauss \boldsymbol{L} e \boldsymbol{U}.

In seguito, analizzeremo la struttura degli elementi non nulli nella matrice \boldsymbol{A}, nella matrice permutata simmetricamente \boldsymbol{P}\boldsymbol{A}\boldsymbol{P}, e nei rispettivi fattori triangolari superiori e inferiori ottenuti dall’algoritmo di Gauss, assumendo che non sia necessario eseguire pivoting numerico.

È importante notare che nella forma originale della matrice \boldsymbol{A}, i fattori risultanti dall’algoritmo di Gauss sono pieni, il che implica un riempimento totale.

Invece, nella forma permutata, non si generano nuovi elementi non nulli. Inoltre, poiché la diagonale di \boldsymbol{L} è composta esclusivamente da 1, non è necessario memorizzarla. Di conseguenza, non viene richiesta ulteriore memoria rispetto a quella necessaria per memorizzare gli elementi non nulli della matrice.

Figura 1: In alto mostriamo la struttura dei non zeri di \boldsymbol{A} e di \boldsymbol{P}\boldsymbol{A}\boldsymbol{P}, in mezzo quella dei fattori di \boldsymbol{A} ed in basso quella dei fattori di \boldsymbol{P}\boldsymbol{A}\boldsymbol{P}.

Note

1. Johann Carl Friedrich Gauss (1777-1855)

2. Perché non si elimina l’incognita w?

3. Ferdinand Georg Frobenius (1849-1917)

4. Se \boldsymbol{A} e \boldsymbol{B} sono matrici ortognali allora (\boldsymbol{A}\boldsymbol{B})^T (\boldsymbol{A}\boldsymbol{B}) = \boldsymbol{B}^T (\boldsymbol{A}^T \boldsymbol{A}) \boldsymbol{B}= \boldsymbol{B}^T \boldsymbol{I}\boldsymbol{B}= \boldsymbol{B}^T \boldsymbol{B}= \boldsymbol{I} e quindi il prodotto è una matrice ortogonale.