Tecniche del pivot

Tecnica del pivot parziale

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

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

Esempio in dimensione n=4

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

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

$$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$$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,$$

$${\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}$${\bar {M}}_{1}=P_{2}M_{1}P_{2}$, al passo successivo ottengo:

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

$$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}$${\bar {M}}_{2}=P_{3}M_{2}P_{3}$, e ${\hat {M}}_{1}=P_{3}{\bar {M}}_{1}P_{3}$${\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$$

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

e pongo $P_{3}P_{2}P_{1}A=PA$$P_{3}P_{2}P_{1}A=PA$.
Voglio verificare che anche le matrici segnate sono ancora matrici elementari di Gauss, e quindi invertibili.

Considero, ad esempio, ${\bar {M}}_{2}=P_{3}M_{2}P_{3}$${\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}}$$

$$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}}$$

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

Allora calcolo il prodotto $P_{3}M_{2}P_{3}$$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}}$$

$$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}}$$

(scambio la terza riga con la quarta.) Invece, moltiplicando a destra per $P_{3}$$P_{3}$, scambio la terza e la quarta colonna, e ottengo 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}}$$

$${\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}$$M_{2}$.
Consideriamo invece la matrice $P_{2}$$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}}$$

$$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}}$$

$$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}$:

$$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}}$$

$$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}$$P_{2}$ non è più un'operazione ammissibile nell'eliminazione di Gauss al passo $k=2$$k=2$, perché si può operare solo nella sottomatrice $2\times 2$$2\times 2$ in basso a destra, quindi non ci sono problemi.

Fattorizzazione

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

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

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

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

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

Allora devo considerare il termine noto permutato ${\bar {\boldsymbol {b}}}$${\bar {\boldsymbol {b}}}$, e risolvo prima il sistma lineare

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

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

e poi

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

$$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$$L$ è una matrice triangolare inferiore con 1 sulla diagonale principale e $M_{ik}$$M_{ik}$ nella parte inferiore.
2. I moltiplicatori hanno espressione
   $$M_{ik}={\frac {|A_{ik}|}{\max _{j=k}^{n}|A_{jk}|}}$$

   $$M_{ik}={\frac {|A_{ik}|}{\max _{j=k}^{n}|A_{jk}|}}$$
   quindi, per ogni $i=k+1,\dots ,n$$i=k+1,\dots ,n$ i coefficienti sono $\leq 1$$\leq 1$.
3. La crescita degli elementi di $U$$U$ è limitata, infatti i coefficienti della "nuova" matrice sono:
   $$A_{ij}^{(k+1)}=A_{ij}^{(k)}-M_{ik}A_{kj}^{(k)}$$

   $$A_{ij}^{(k+1)} = A_{ij}^{(k)}-M_{ik} A_{kj}^{(k)}$$
   mentre per $i,j=1,\dots k$$i,j=1,\dots k$ gli elementi della matrice rimangono invariati, e quindi
   $$A_{ij}^{(k+1)}=A_{ij}^{(k)}$$

   $$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}|$$

   $$|A_{ij}|^{(k+1)}\leq |A_{ij}|^{(k)}+|M_{ik}||A_{kj}|^{(k)}\leq |A_{ij}|+|A_{kj}|$$
   Chiamo $A_{m}^{(k)}$$A_{m}^{(k)}$ il massimo modulo degli elementi della matrice $A^{(k)}$$A^{(k)}$. Allora, passando ai massimi nella disuguaglianza sopra,
   $$|A_{m}^{(k+1)}|\leq |A_{m}^{(k)}|+|A_{m}^{(k)}|=2|A_{m}^{(k)}|$$

   $$|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)}|$$

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

Errore algoritmico

Teorema 2.7

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

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

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

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

 

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

Teorema 2.8

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

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

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

con

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

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

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

 

Esempio 2.7

Considero la matrice

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

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

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

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

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

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

$$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}$$

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

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

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

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

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

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

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

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

$${\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$$U$.

 

Tecnica del pivot totale

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

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

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

assumo di poter fare permutazioni di colonne ${\hat {P}}_{k}$${\hat {P}}_{k}$ oltre che permutazioni di righe $P_{k}$$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$$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}$$

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

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

 

Come prima

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

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

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

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

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

devo inserire la ${\hat {P}}$${\hat {P}}$, allora pongo ${\hat {P}}{\boldsymbol {z}}={\boldsymbol {x}}$${\hat {P}}{\boldsymbol {z}}={\boldsymbol {x}}$

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

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

Calcolo ${\bar {\boldsymbol {b}}}$${\bar {\boldsymbol {b}}}$ permutazione del termine noto, e risolvo quindi:

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

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

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

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

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

$${\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}}}$$

$$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$$

$$\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}}}$$

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

Invece la fattorizzazione $LU$$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]$$

$$\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$$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).