Calcolo numerico/Tecniche del pivot

Tecnica del pivot parziale

Supponiamo di essere al $k$-esimo passo di fattorizzazione, in cui bisogna azzerare i coefficienti della $k$-esima colonna a partire dall'entrata $k+1$. Cerchiamo l'indice ${\bar {i}}$ tale che $A_{{\bar {i}},k}=\max _{i=k}^{n}|a_{ik}|$, e scambiamo quindi la ${\bar {i}}$-esima riga con la $k$-esima, applicando alla matrice $A^{(k)}$ la matrice di permutazione $p_{k}$, che è l'identità in cui scambio la ${\bar {i}}$-esima riga con la $k$-esima.

Bisogna verificare che anche moltiplicando per la matrice di permutazione, la fattorizzazione rimane del tipo $LU$.

Esempio in dimensione n=4

In dimensione $n=4$, all'ultimo passo dell'eliminazione si ottiene

$U=M_{3}P_{3}M_{2}P_{2}M_{1}P_{1}A,$

(le matrici di permutazione sono alternate a quelle elementari)

Le matrici di permutazione hanno la proprietà che $P_{k}P_{k}=I$, allora, inserisco il doppio prodotto nell'uguaglianza nelle posizioni opportune, per far comparire vicine tra loro le matrici di permutazione.

Ad esempio:

${\hbox{passo 1}},\quad U=M_{3}P_{3}M_{2}P_{2}M_{1}P_{2}P_{2}P_{1}A,$

Ponendo ${\bar {M}}_{1}=P_{2}M_{1}P_{2}$, al passo successivo otteniamo:

$U=M_{3}P_{3}M_{2}P_{3}P_{3}{\bar {M}}_{1}P_{3}P_{3}P_{2}P_{1}A$

e ponendo ${\bar {M}}_{2}=P_{3}M_{2}P_{3}$, e ${\hat {M}}_{1}=P_{3}{\bar {M}}_{1}P_{3}$, ottengo:

$U=M_{3}{\bar {M}}_{2}{\hat {M}}_{1}P_{3}P_{2}P_{1}A$

e si pone $P_{3}P_{2}P_{1}A=PA$.

Vogliamo verificare che anche le matrici segnate sono ancora matrici elementari di Gauss, e quindi invertibili.

Consideriamo, ad esempio, ${\bar {M}}_{2}=P_{3}M_{2}P_{3}$

$M_{2}={\begin{array}{cccc}1&0&0&0\\0&1&0&0\\0&M_{32}&1&0\\0&M_{42}&0&1\end{array}}$

Invece

$P_{3}={\begin{array}{cccc}1&0&0&0\\0&1&0&0\\0&0&0&1\\0&0&1&0\end{array}}$

Allora calcoliamo il prodotto $P_{3}M_{2}P_{3}$:

$P_{3}M_{2}={\begin{array}{cccc}1&0&0&0\\0&1&0&0\\0&M_{42}&0&1\\0&M_{32}&1&0\end{array}}$

(scambiiamo la terza riga con la quarta). Invece, moltiplicando a destra per $P_{3}$, scambiamo la terza e la quarta colonna, e otteniamo la matrice

${\bar {M}}_{2}={\begin{array}{cccc}1&0&0&0\\0&1&0&0\\0&M_{42}&1&0\\0&M_{32}&0&1\end{array}}$

che ha la stessa forma di $M_{2}$.

Consideriamo invece la matrice $P_{2}$ che scambia la seconda riga con la terza:

$p_{2}={\begin{array}{cccc}1&0&0&0\\0&0&1&0\\0&1&0&0\\0&0&0&1\end{array}}$

Allora

$P_{2}M_{2}={\begin{array}{cccc}1&0&0&0\\0&M_{32}&1&0\\0&1&0&0\\0&M_{42}&0&1\end{array}}$

e moltiplicando per $P_{2}$:

$P_{2}M_{2}P_{2}={\begin{array}{cccc}1&0&0&0\\0&1&M_{32}&0\\0&0&1&0\\0&0&M_{42}&1\end{array}}$

che non è più una matrice elementare di Gauss.

In ogni caso, la moltiplicazione per $P_{2}$ non è più un'operazione ammissibile nell'eliminazione di Gauss al passo $k=2$, perché si può operare solo nella sottomatrice $2\times 2$ in basso a destra, quindi non ci sono problemi.

Fattorizzazione

Dato il sistema $A{\boldsymbol {x}}={\boldsymbol {b}}$, moltiplichiamo per $P$ ambo i membri per far comparire la matrice di permutazione, e otteniamo:

$PA{\boldsymbol {x}}=P{\boldsymbol {b}}$

$PA=LU\quad \longrightarrow \quad LU{\boldsymbol {x}}=P{\boldsymbol {b}}$

Allora bisogna considerare il termine noto permutato ${\bar {\boldsymbol {b}}}$, e risolviamo prima il sistema lineare

$L{\boldsymbol {y}}={\bar {\boldsymbol {b}}}$

e poi

$U{\boldsymbol {x}}={\boldsymbol {y}}$

Effetti del pivot parziale

Il pivot parziale ha un effetto stabilizzante, ovvero riduce l'errore algoritmico. In particolare:

1. Come prima, $L$ è una matrice triangolare inferiore con 1 sulla diagonale principale e $M_{ik}$ nella parte inferiore.

2. I moltiplicatori hanno espressione

$M_{ik}={\frac {|A_{ik}|}{\max _{j=k}^{n}|A_{jk}|}}$

quindi, per ogni $i=k+1,\dots ,n$ i coefficienti sono $\leq 1$.

3. La crescita degli elementi di $U$ è limitata, infatti i coefficienti della "nuova" matrice sono:

$A_{ij}^{(k+1)}=A_{ij}^{(k)}-M_{ik}A_{kj}^{(k)}$

mentre per $i,j=1,\dots k$ gli elementi della matrice rimangono invariati, e quindi

$A_{ij}^{(k+1)}=A_{ij}^{(k)}$

in definitiva

$|A_{ij}|^{(k+1)}\leq |A_{ij}|^{(k)}+|M_{ik}||A_{kj}|^{(k)}\leq |A_{ij}|+|A_{kj}|$

Chiamiamo $A_{m}^{(k)}$ il massimo modulo degli elementi della matrice $A^{(k)}$. Allora, passando ai massimi nella disuguaglianza sopra,

$|A_{m}^{(k+1)}|\leq |A_{m}^{(k)}|+|A_{m}^{(k)}|=2|A_{m}^{(k)}|$

e applicando a ritroso la relazione

$|A_{m}|^{(n)}\leq 2^{n-1}|A_{m}^{(1)}|$

Errore algoritmico

Teorema 2.7

Sia $A\in \mathbb {C} _{n,n}$ tale che i suoi elementi siano numeri macchina (focus sull'errore algoritmico). Chiamiamo ${\bar {L}},{\bar {U}}$ le matrici effettivamente calcolate, allora ${\bar {L}}*{\bar {U}}=A+\delta _{A}$ dove vale la relazione

$|\delta _{A}|\leq 2nu(|A|+|{\bar {L}}|*|{\bar {U}}|)+O(u^{2})$

(a livello notazionale, chiamiamo $|B|$ la matrice che ha come componenti i moduli degli elementi di $B$)

In questa relazione, oltre alla dipendenza dalla dimensione, da $A$ e dalla precisione di macchina, il fatto significativo è che compaiono ${\bar {L}},{\bar {U}}$.

Teorema 2.8

Siano $A,B\in \mathbb {C} _{n,n}$ e, considerando ${\bar {L}},{\bar {U}}$ effettivamente calcolate, chiamiamo ${\bar {\boldsymbol {y}}}$ la soluzione del sistema ${\bar {L}}{\boldsymbol {y}}={\boldsymbol {b}}$ e ${\bar {\boldsymbol {x}}}$ soluzione di ${\bar {U}}{\boldsymbol {x}}={\boldsymbol {y}}$. Allora ${\bar {\boldsymbol {x}}}$ è la soluzione di

$(A+\delta _{A}){\bar {\boldsymbol {x}}}={\boldsymbol {b}}$

con

$\delta _{A}\leq 4nu(|A|+|{\bar {L}}||{\bar {U}}|+o(u^{2}))$

(interpretiamo l'errore algoritmico come errore inerente del problema perturbato)

Esempio 2.7

Consideriamo la matrice

$A=\left({\begin{array}{cc}\varepsilon &1\\1&0\end{array}}\right)$

${\boldsymbol {b}}=(1+\varepsilon ,1)$

dove $\varepsilon$ è scelto in modo che la soluzione del sistema lineare $A{\boldsymbol {x}}={\boldsymbol {b}}$ sia $(1,1)$. Il problema è ben condizionato, infatti

$A^{-1}={\begin{array}{cc}0&1\\1&-\varepsilon \end{array}}$

$K(A)=\vert A\vert _{\infty }*\vert A^{-1}\vert _{\infty }=(1+\varepsilon )^{2}$

e quindi se $\varepsilon$ è piccolo, $K(A)\simeq 1$.

Per quanto riguarda la stabilità, senza tecnica del pivot:

$L={\begin{array}{cc}1&0\\1/\varepsilon &1\end{array}}$

$U={\begin{array}{cc}\varepsilon &1\\0&-1/\varepsilon \end{array}}$

${\bar {\boldsymbol {x}}}=(0,1)$

mentre la soluzione esatta ha entrambe le componenti uguali a 1, e il problema è dato dalla crescita dei coefficienti di $U$.

Tecnica del pivot totale

Per ogni passo $k$, cerchiamo ${\bar {i}},{\bar {j}}$ tale che

$|A_{{\bar {i}},{\bar {j}}}|=\max _{i,j=k}^{n}|A_{i,j}|$

assumiamo di poter fare permutazioni di colonne ${\hat {P}}_{k}$ oltre che permutazioni di righe $P_{k}$.

Il problema che sorge è che, facendo permutazioni di colonne, si permuta l'ordine delle incognite.

Esempio 2.8

Nella catena di moltiplicazioni di matrici elementari, nel caso $n=4$,

$U=M_{3}P_{3}M_{2}P_{2}M_{1}P_{1}A{\hat {P}}_{1}{\hat {P}}_{2}{\hat {P}}_{3}$

ma poniamo ${\hat {P}}_{1}{\hat {P}}_{2}{\hat {P}}_{3}={\hat {P}}$. L'altro pezzo di $U$ viene trattato come prima.

Come prima

$PA{\hat {P}}=LU$

e per risolvere il sistema $A{\boldsymbol {x}}={\boldsymbol {b}}$ permutiamo a sinistra,

$PA{\boldsymbol {x}}=P{\boldsymbol {b}}$

bisogna inserire la ${\hat {P}}$, allora poniamo ${\hat {P}}{\boldsymbol {z}}={\boldsymbol {x}}$

$PA{\hat {P}}{\boldsymbol {z}}=P{\boldsymbol {b}}$

Calcoliamo ${\bar {\boldsymbol {b}}}$ permutazione del termine noto, e risolviamo quindi:

$L{\boldsymbol {y}}={\bar {\boldsymbol {b}}}$

$U{\boldsymbol {z}}={\boldsymbol {y}}$

${\hat {P}}{\boldsymbol {z}}={\boldsymbol {x}},\,\quad {\hbox{è una moltiplicazione matrice per vettore}}$

Questo procedimento ha un effetto più stabilizzante rispetto al pivot parziale, infatti si ha:

$A_{k}\leq {\sqrt {k\prod _{j=2}^{n}j^{j-1}}}$

e non sono state trovate matrici per cui vale l'uguaglianza.

Costo computazionale

La risoluzione forward o backward, considerando il numero di operazioni moltiplicative, ha come costo computazionale

$\sum _{i=1}^{n}i$

infatti, per ogni incognita, si ha l'espressione

$x_{i}={\frac {b_{i}-\sum _{j=1}^{i-1}a_{ij}x_{j}}{a_{ii}}}$

Invece la fattorizzazione $LU$ ha costo

$\sum _{k=1}^{n-1}[(\sum _{i=k+1}^{n}1)+\sum _{i=k+1}^{n}\sum _{j=k+1}^{n}1]$

che è dell'ordine di $n^{3}/3$ (il primo addendo serve per il calcolo dei moltiplicatori, e il secondo per il calcolo del blocco quadrato in basso a destra).