Pivoting

Se, per esempio, avessimo che $a_{11}^{(1)}=0$ con il metodo precedente non potremmo proseguire, ma siamo sicuri che all'interno della stessa colonna esiste un elemento non nullo (altrimenti avremmo una colonna tutta nulla che comporterebbe la singolarità della matrice, in opposizione alle ipotesi). Cerchiamo allora un altro elemento diverso da zero nel seguente modo, limitandoci al caso reale:

\begin{displaymath}a_{k1,1}^{(1)}=max_i \vert a_{i1}^{(1)}\vert\end{displaymath}

prendiamo cioè l'elemento che ha modulo massimo sulla prima colonna. Individuato il suo indice di riga, scambiamo allora la prima riga con la k1-esima, in modo da ottenere un elemento valido sulla diagonale. Per effettuare questo scambio utiliziamo una matrice di permutazione così costruita:

\begin{displaymath}
P_1= \left(
\begin{array}{ccccccccc}
0 & \cdots & \cdots ...
...dots & 0 & \cdots & \cdots & \cdots & 1
\end{array}
\right)
\end{displaymath}

che non è altro che la matrice identità con la prima e la k1-esima riga scambiate. $P_1$ è una matrice simmetrica ed ortogonale ( $P_1=P_1^T=P_1^{-1}$) che scambia la prima riga con la k1-esima:

\begin{displaymath}
P_1 \left(
\begin{array}{c}
1\\
2\\
\vdots\\
k1-1\\...
...
k1-1\\
1\\
k1+1\\
\vdots\\
n
\end{array}
\right)
\end{displaymath}

Dopo aver costruito la matrice di permutazione elementare, possiamo eseguire $P_1A^{(1)}$ che scambia nella matrice $A$ le righe $1$ e $k1$. Ristrutturata in questo modo, la matrice $P_1A^{(1)}$ possiede tutti i requisiti richiesti dal metodo di eliminazione di Gauss, che dunque possiamo applicare; ottenuta dunque la matrice $L_1$ come visto per l'algoritmo precedente, portiamo a termine il primo passo ottenendo come risultato:

\begin{displaymath}
L_1P_1A^{(1)}= \left(
\begin{array}{cccc}
a_{k1,1}^{(1)} ...
... \cdots & a_{nn}^{(2)}
\end{array}
\right)
\equiv A^{(2)}
\end{displaymath}

Un fatto interessante da notare è che, calcolando il vettore elementare di Gauss come già visto, e cioè come

\begin{displaymath}g=\frac{1}{a_{k1,1}^{(1)}}({0,a_{21}^{(1)}, \cdots ,a_{n1}^{(1)}})^T\end{displaymath}

i suoi elementi avranno modulo massimo pari ad uno: infatti abbiamo scelto $a_{k1,1}$ come il massimo della prima colonna, e come favorevole conseguenza abbiamo che gli errori sono ben limitati in aritmetica finita; questa considerazione è valida per la matrice $L$ mentre per $U$ non possiamo dire niente.

Quello che abbiamo descritto ora è il metodo di pivoting parziale, in quanto viene scelto il pivot (l'elemento diagonale utilizzato per costruire il vettore elementare di Gauss) non in modo statico come avveniva per la fattorizzazione $A=LU$ ma viene ricercato all'interno della colonna in esame l'elemento di modulo massimo, e con quello si procede con l'algoritmo di eliminazione di Gauss.

Ottenuta la matrice $A^{(2)}$ possiamo continuare ad applicare l'algoritmo di scambio e di fattorizzazione. Ripetuto questo passo per $n-1$ volte otteniamo

\begin{displaymath}
L_{n-1}P_{n-1}\cdots L_1P_1A^{(1)}=\left(
\begin{array}{cc...
...cdots & 0 & a_{kn,n}^{(n)} \\
\end{array}
\right) \equiv U
\end{displaymath}

Applicando il pivoting, siamo dunque ancora in grado di eseguire l'eliminazione di Gauss. Se l'elemento di massimo modulo di una colonna avesse valore pari a zero, questo implicherebbe l'esistenza di un blocco diagonale con una colonna nulla e questo implicherebbe la singolarità di $A$ in opposizione alle nostre ipotesi: infatti sia

\begin{displaymath}
B= \left(
\begin{array}{c\vert c}
B_1 & C \\ \hline \mathcal{O} & B_2
\end{array}
\right)
\end{displaymath}

a blocchi quadrati, allora

\begin{displaymath}det(B)=det(B_1)det(B_2)\end{displaymath}

e

\begin{displaymath}det(B) \neq 0 \qquad \Rightarrow \qquad \left\{
\begin{array...
...& \neq & 0 \\
det(B_2) & \neq & 0 \\
\end{array}
\right.
\end{displaymath}

cioè la non singolarità di $B$ è derivata dalla non singolarità dei suoi blocchi diagonali.

Abbiamo ottenuto la matrice $U$, ma come possiamo ottenere $L$? Facciamo un esempio con $n=4$ tenendo a mente che $P_i$ è simmetrica ed ortogonale:

\begin{displaymath}
\begin{array}{l}
L_3P_3L_2P_2L_1P_1A=U \\
\\
L_3P_3L_2...
...
L_3(P_3L_2P_3)(P_3P_2L_1P_2P_3)(P_3P_2P_1)A\\
\end{array}
\end{displaymath}

da cui con opportune sostituzioni si ha

\begin{displaymath}\hat{L}_3\hat{L}_2\hat{L}_1PA=U\mbox{.}\end{displaymath}

In generale definiamo

\begin{displaymath}
\begin{array}{l}
P=P_{n-1} \cdots P_1 \\
\\
\hat{L}_i=...
...} \cdots P_{n-1} \qquad
(\hat{L}_{n-1}=L_{n-1})
\end{array}
\end{displaymath}

tali che

\begin{displaymath}L_{n-1}P_{n-1} \cdots L_2P_2L_1P_1A=\hat{L}_{n-1} \cdots
\hat{L}_1PA (=U)\end{displaymath}

Se le matrici $\hat{L}_i$ avessero la stessa struttura delle matrici elementari di Gauss avremmo trovato la nostra fattorizzazione poichè potremmo nuovamente definire $L^{-1}=\prod_{i=1}^{n-1} \hat{L}_I$ per poi ottenere la fattorizzazione $PA=LU$; cerchiamo di verificare allora che $\hat{L}_i$ ed $L_i$ hanno la stessa struttura: partendo dal vettore elementare di Gauss

\begin{displaymath}g_i={(\underbrace{0\cdots 0}_{i}, g_{i+1}^{(i)}, \cdots ,
g_n^{(i)})}^T\end{displaymath}

cerchiamo di vedere com'è fatta $\hat{L}_i$

\begin{displaymath}
\begin{array}{lcl}
\hat{L}_i & = & P_{n-1} \cdots P_{i+1}(...
...\cdots P_{i+1}g_i)(e_i^TP_{i+1} \cdots
P_{n-1})
\end{array}
\end{displaymath}

prestiamo attenzione alla forma di $e_i^TP_j$ per $j>i$. $e_i^T$ applicato ad una matrice seleziona l'i-esima riga, ma $P_{i+1}$, ad esempio, ha le prime $i$ righe dell'identità, quindi

\begin{displaymath}
\begin{array}{l}
e_i^TP_j=e_i^T \quad \forall j>i \\
\\
e_i^TP_{i+1} \cdots P_{n-1}=e_i^T
\end{array}
\end{displaymath}

Ci siamo ridotti a

\begin{displaymath}\hat{L}_i=I-(P_{n-1} \cdots P_{i+1}g_i)e_i^T\end{displaymath}

guardiamo adesso alla forma di $P_{i+1}g_i$: $P_{i+1}$ è la matrice elementare di permutazione che scambia $i+1$ con $ki(>i+1)$ che applicata a $g_i$ produce una permutazione degli elementi del vettore dal'indice $i+1$ in poi, e così per tutte le altre matrici $P_k$; chiamando $\hat{g}_i$ il nuovo vettore, questo ha forma

\begin{displaymath}\hat{g}_i={(\underbrace{0\cdots 0}_{i}, \hat{g}_{i+1}^{(i)}, \cdots ,
\hat{g}_n^{(i)})}^T\mbox{.}\end{displaymath}

L'applicazione delle matrici di permutazione al vettore $g_i$ ha prodotto un vettore con gli elementi $g_{i+1}^{(i)}, \cdots , g_n^{(i)}$ permutati, ma che mantiene ancora la struttura del vettore elementare di Gauss. Unendo i due risultati ottenuti si perviene a

\begin{displaymath}\hat{L}_i=I-\hat{g}_ie_i^T\end{displaymath}

proprio una matrice elementare di Gauss, da cui possiamo scrivere:

\begin{displaymath}L^{-1}PA=U \qquad \longrightarrow \qquad PA=LU\end{displaymath}

la fattorizzazione che stavamo cercando. Il metodo di pivoting ci consente di applicare l'eliminazione di Gauss anche in caso di elementi diagonali nulli, e quindi con qualsiasi matrice non singolare: basta applicare una matrice di permutazione che la renda conforme alle richieste della fattorizzazione $LU$.

Morpheus 2004-01-04