Salta al contenuto
Note per Studenti Esercizio - pivot piccolo in aritmetica a quattro cifre e costo della fattorizzazione

Esercizio - pivot piccolo in aritmetica a quattro cifre e costo della fattorizzazione

Teoria: Fattorizzazione LU, pivoting e CholeskyRisolvere $A\mathbf x=\mathbf b$ con Gauss equivale a fattorizzare $A=LU$ ($L$ triangolare inferiore con $1$ sulla diagonale e i moltiplicatori $l_{ik}=\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}$ sotto la diagonale, $U$ la matrice a scala finale) e poi risolvere $L\mathbf y=\mathbf b$ (in avanti) e $U\mathbf x=\mathbf y$ (all'indietro). Esiste (ed è unica) se i minori principali di testa $\det A_k$, $k<n$, sono $\ne0$. Costo $\frac23n^3$ flop per la fattorizzazione, $2n^2$ per ogni termine noto; $\det A=\prod u_{kk}$. Con un pivot piccolo l'algoritmo è instabile: il pivoting parziale (scambio di righe che porta in posizione di pivot l'elemento di modulo massimo della colonna) dà $PA=LU$ con $|l_{ik}|\le1$ ed esiste per ogni $A$ invertibile. Se $A$ è simmetrica definita positiva esiste $A=LL^T$ (Cholesky), con $l_{jj}=\sqrt{a_{jj}-\sum_{k<j}l_{jk}^2}$ e $l_{ij}=\frac{a_{ij}-\sum_{k<j}l_{ik}l_{jk}}{l_{jj}}$; costo $\frac13n^3$, senza pivoting.Fattorizzazione LU, pivoting e Cholesky →, Aritmetica di macchina - floating point, precisione di macchina e cancellazioneUn calcolatore rappresenta solo un insieme finito di numeri F(β,t,L,U); ogni numero reale viene arrotondato con errore relativo al più u = ½β^(1−t) (u ≈ 1,1·10⁻¹⁶ in doppia precisione) e ogni operazione è esatta più un errore (1+δ), |δ| ≤ u. Conseguenze: overflow e underflow, somma non associativa, assorbimento, e soprattutto la cancellazione numerica nella differenza di numeri quasi uguali, che si evita riformulando l'espressione.Aritmetica di macchina - floating point, precisione di macchina e cancellazione →.

Testo. Si risolve il sistema (10−5111)(x1x2)=(12)\begin{pmatrix}10^{-5}&1\\1&1\end{pmatrix}\begin{pmatrix}x_1\\x_2\end{pmatrix}=\begin{pmatrix}1\\2\end{pmatrix} con un calcolatore che lavora con 4 cifre significative decimali e arrotonda ogni operazione.

  1. Calcolare la soluzione esatta (a 6 cifre).
  2. Eseguire la LU senza scambi in aritmetica a 4 cifre e dare la soluzione ottenuta.
  3. Ripetere con il pivoting parziale.
  4. Un calcolatore esegue 10910^9 flop/s. Quanto tempo serve per risolvere un sistema 2000×20002000\times2000 con LU, e quanto per 100 termini noti diversi con la stessa matrice, riusando la fattorizzazione o rifacendola ogni volta?

Soluzione

1. Soluzione esatta. Sottraendo la prima equazione dalla seconda: x1(1−10−5)=2−1=1x_1(1-10^{-5})=2-1=1, quindi x1=11−10−5=1,00001000,x2=2−x1=0,99998999≈0,99999.x_1=\frac1{1-10^{-5}}=1{,}00001000,\qquad x_2=2-x_1=0{,}99998999\approx0{,}99999 . La matrice è ben condizionata: ∥A∥∞=2\|A\|_\infty=2, A−1=110−5−1(1−1−110−5)A^{-1}=\frac1{10^{-5}-1}\begin{pmatrix}1&-1\\-1&10^{-5}\end{pmatrix}, ∥A−1∥∞=21−10−5≈2\|A^{-1}\|_\infty=\frac{2}{1-10^{-5}}\approx2, quindi κ∞≈4\kappa_\infty\approx4. Piccole perturbazioni dei dati non fanno danni: un eventuale errore grande è colpa dell'algoritmo.

2. Senza scambi (4 cifre).

  • Moltiplicatore: l21=110−5=105=1,000⋅105l_{21}=\frac{1}{10^{-5}}=10^5=1{,}000\cdot10^5.
  • u22=1−l21⋅1=1−100000=−99999u_{22}=1-l_{21}\cdot1=1-100000=-99999, che a 4 cifre diventa −1,000⋅105-1{,}000\cdot10^5 (il 11 è sparito).
  • y2=2−l21⋅1=2−100000=−99998→−1,000⋅105y_2=2-l_{21}\cdot1=2-100000=-99998\to-1{,}000\cdot10^5.
  • x2=y2u22=−1,000⋅105−1,000⋅105=1,000x_2=\frac{y_2}{u_{22}}=\frac{-1{,}000\cdot10^5}{-1{,}000\cdot10^5}=1{,}000 (accettabile: esatto 0,999990{,}99999).
  • x1=1−1⋅x210−5=1−1,00010−5=0x_1=\frac{1-1\cdot x_2}{10^{-5}}=\frac{1-1{,}000}{10^{-5}}=\mathbf 0: sbagliato del 100% (esatto 1,000011{,}00001). Il piccolo errore su x2x_2 (di ordine 10−510^{-5}, invisibile) viene diviso per il pivot 10−510^{-5} e amplificato di 10510^5: errore assoluto ≈1\approx1 su x1x_1.

3. Con il pivoting. Si scambiano le righe: (1110−51)\begin{pmatrix}1&1\\10^{-5}&1\end{pmatrix} con termini noti (2,1)(2,1).

  • l21=10−51=10−5l_{21}=\frac{10^{-5}}1=10^{-5}; u22=1−10−5⋅1=0,99999→1,000u_{22}=1-10^{-5}\cdot1=0{,}99999\to1{,}000.
  • y2=1−10−5⋅2=0,99998→1,000y_2=1-10^{-5}\cdot2=0{,}99998\to1{,}000.
  • x2=1,0001,000=1,000x_2=\frac{1{,}000}{1{,}000}=1{,}000; x1=2−1⋅1,0001=1,000x_1=\frac{2-1\cdot1{,}000}{1}=1{,}000. Soluzione (1,000; 1,000)(1{,}000;\ 1{,}000), con errore ≤10−5\le10^{-5} come ci si aspetta da un calcolo a 4 cifre (l'errore di arrotondamento di ciascuna operazione è ≤0,5⋅10−3\le0{,}5\cdot10^{-3} relativo).

4. Costo. Per n=2000n=2000: LU 23n3=23⋅8⋅109=5,3⋅109\frac23n^3=\frac23\cdot8\cdot10^9=5{,}3\cdot10^9 flop ⇒5,3\Rightarrow5{,}3 s. Ogni coppia di sostituzioni: 2n2=8⋅1062n^2=8\cdot10^6 flop ⇒8\Rightarrow8 ms. Per 100 termini noti:

  • riusando la fattorizzazione: 5,3+100⋅0,008=6,15{,}3+100\cdot0{,}008=6{,}1 s;
  • rifacendo tutto ogni volta: 100⋅(5,3+0,008)≈531100\cdot(5{,}3+0{,}008)\approx531 s. Il risparmio è di un fattore ≈87\approx87. Se la matrice è anche SPD, Cholesky dimezza il costo di fattorizzazione (≈2,7\approx2{,}7 s).

Risultati: esatto (1,00001; 0,99999)(1{,}00001;\ 0{,}99999); senza pivoting (0; 1)(0;\ 1) (errore 100%100\%); con pivoting (1,000; 1,000)(1{,}000;\ 1{,}000); tempo 5,35{,}3 s + 88 ms per termine noto.

Errori tipici. Attribuire l'errore al condizionamento del problema (κ≈4\kappa\approx4 è piccolo: il problema è buono, l'algoritmo senza pivoting no); dimenticare di arrotondare a 4 cifre ogni risultato intermedio; calcolare il costo di LU come n3n^3 senza il fattore 23\frac23.

Versione ripasso

(Fattorizzazione LU, pivoting e CholeskyRisolvere $A\mathbf x=\mathbf b$ con Gauss equivale a fattorizzare $A=LU$ ($L$ triangolare inferiore con $1$ sulla diagonale e i moltiplicatori $l_{ik}=\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}$ sotto la diagonale, $U$ la matrice a scala finale) e poi risolvere $L\mathbf y=\mathbf b$ (in avanti) e $U\mathbf x=\mathbf y$ (all'indietro). Esiste (ed è unica) se i minori principali di testa $\det A_k$, $k<n$, sono $\ne0$. Costo $\frac23n^3$ flop per la fattorizzazione, $2n^2$ per ogni termine noto; $\det A=\prod u_{kk}$. Con un pivot piccolo l'algoritmo è instabile: il pivoting parziale (scambio di righe che porta in posizione di pivot l'elemento di modulo massimo della colonna) dà $PA=LU$ con $|l_{ik}|\le1$ ed esiste per ogni $A$ invertibile. Se $A$ è simmetrica definita positiva esiste $A=LL^T$ (Cholesky), con $l_{jj}=\sqrt{a_{jj}-\sum_{k<j}l_{jk}^2}$ e $l_{ij}=\frac{a_{ij}-\sum_{k<j}l_{ik}l_{jk}}{l_{jj}}$; costo $\frac13n^3$, senza pivoting.Fattorizzazione LU, pivoting e Cholesky →)

(10−5111)x=(12)\begin{pmatrix}10^{-5}&1\\1&1\end{pmatrix}\mathbf x=\begin{pmatrix}1\\2\end{pmatrix}, 4 cifre significative.

  • Esatto: x1=11−10−5=1,00001x_1=\frac1{1-10^{-5}}=1{,}00001, x2=0,99999x_2=0{,}99999; κ∞≈4\kappa_\infty\approx4, problema ben condizionato.
  • Senza scambi: l21=105l_{21}=10^5, u22=1−105→−1,000⋅105u_{22}=1-10^5\to-1{,}000\cdot10^5, y2=2−105→−1,000⋅105y_2=2-10^5\to-1{,}000\cdot10^5, x2=1,000x_2=1{,}000, x1=1−1,00010−5=0x_1=\frac{1-1{,}000}{10^{-5}}=0: errore 100%100\% (colpa dell'algoritmo).
  • Con pivoting: l21=10−5l_{21}=10^{-5}, u22=1,000u_{22}=1{,}000, y2=1,000y_2=1{,}000, x2=1,000x_2=1{,}000, x1=2−1,000=1,000x_1=2-1{,}000=1{,}000.
  • n=2000n=2000, 10910^9 flop/s: LU 23n3=5,3⋅109→5,3\frac23n^3=5{,}3\cdot10^9\to5{,}3 s; sostituzioni 2n2=8⋅106→82n^2=8\cdot10^6\to8 ms. 100 termini noti: 6,16{,}1 s riusando LU contro 531531 s rifacendo (fattore 8787).

Errori tipici: errore attribuito al condizionamento; arrotondamenti intermedi saltati; n3n^3 senza 23\frac23.

Teoria collegata