Salta al contenuto
Note per Studenti Esercizio - residuo ed errore in un sistema ben e mal condizionato

Esercizio - residuo ed errore in un sistema ben e mal condizionato

Teoria: Norme di matrici e condizionamento dei sistemi lineariNorme di vettore $|x|1=\sum|x_i|$, $|x|2=\sqrt{\sum x_i^2}$, $|x|\infty=\max|x_i|$. Norma di matrice indotta $|A|=\max{x\ne0}\frac{|Ax|}{|x|}$: $|A|1$ = massima somma per colonne, $|A|\infty$ = massima somma per righe, $|A|2=\sqrt{\lambda{\max}(A^TA)}$ (= $\rho(A)$ se $A$ simmetrica); valgono $|Ax|\le|A||x|$, $|AB|\le|A||B|$, $\rho(A)\le|A|$. Numero di condizionamento $\kappa(A)=|A||A^{-1}|\ge1$ ($\kappa_2=\frac{\sigma_{\max}}{\sigma_{\min}}$, per SPD $\frac{\lambda_{\max}}{\lambda_{\min}}$, per ortogonali 1). Se $A\mathbf x=\mathbf b$ e $A(\mathbf x+\delta\mathbf x)=\mathbf b+\delta\mathbf b$: $\frac{|\delta\mathbf x|}{|\mathbf x|}\le\kappa(A)\frac{|\delta\mathbf b|}{|\mathbf b|}$ (D); con residuo $\mathbf r=\mathbf b-A\hat{\mathbf x}$: $\frac1\kappa\frac{|\mathbf r|}{|\mathbf b|}\le\frac{|\mathbf x-\hat{\mathbf x}|}{|\mathbf x|}\le\kappa\frac{|\mathbf r|}{|\mathbf b|}$. Si perdono circa $\log_{10}\kappa$ cifre; un residuo piccolo non garantisce un errore piccolo se $\kappa$ è grande.Norme di matrici e condizionamento dei sistemi lineari →, 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 →.

Testo. Due programmi forniscono soluzioni approssimate di due sistemi.

  • Sistema I: M=(410131012)M=\begin{pmatrix}4&1&0\\1&3&1\\0&1&2\end{pmatrix}, b=(6,10,8)T\mathbf b=(6,10,8)^T (soluzione esatta x=(1,2,3)T\mathbf x=(1,2,3)^T), soluzione calcolata x^=(1,01; 1,98; 3,00)T\hat{\mathbf x}=(1{,}01;\ 1{,}98;\ 3{,}00)^T.
  • Sistema II: A=(1224,002)A=\begin{pmatrix}1&2\\2&4{,}002\end{pmatrix}, b=(3; 6,002)T\mathbf b=(3;\ 6{,}002)^T (soluzione esatta (1,1)T(1,1)^T), soluzione calcolata x^=(3, 0)T\hat{\mathbf x}=(3,\ 0)^T.

Per ciascuno:

  1. calcolare il residuo r=b−Ax^\mathbf r=\mathbf b-A\hat{\mathbf x} e il residuo relativo ∥r∥∞∥b∥∞\frac{\|\mathbf r\|_\infty}{\|\mathbf b\|_\infty};
  2. calcolare l'errore relativo vero ∥x−x^∥∞∥x∥∞\frac{\|\mathbf x-\hat{\mathbf x}\|_\infty}{\|\mathbf x\|_\infty};
  3. verificare le disuguaglianze con κ∞\kappa_\infty (sistema I: κ∞(M)=409=4,44\kappa_\infty(M)=\frac{40}9=4{,}44; sistema II: κ∞(A)=18012\kappa_\infty(A)=18012);
  4. dire se il residuo relativo sarebbe un buon criterio d'arresto in ciascun caso.

Soluzione

Sistema I.

  • Mx^=(4⋅1,01+1,98, 1,01+3⋅1,98+3,00, 1,98+2⋅3,00)=(6,02; 9,95; 7,98)M\hat{\mathbf x}=(4\cdot1{,}01+1{,}98,\ 1{,}01+3\cdot1{,}98+3{,}00,\ 1{,}98+2\cdot3{,}00)=(6{,}02;\ 9{,}95;\ 7{,}98).
  • Residuo: r=(6−6,02, 10−9,95, 8−7,98)=(−0,02; 0,05; 0,02)\mathbf r=(6-6{,}02,\ 10-9{,}95,\ 8-7{,}98)=(-0{,}02;\ 0{,}05;\ 0{,}02), ∥r∥∞=0,05\|\mathbf r\|_\infty=0{,}05, residuo relativo 0,0510=5⋅10−3\frac{0{,}05}{10}=5\cdot10^{-3}.
  • Errore: x−x^=(−0,01; 0,02; 0)\mathbf x-\hat{\mathbf x}=(-0{,}01;\ 0{,}02;\ 0), ∥⋅∥∞=0,02\|\cdot\|_\infty=0{,}02, errore relativo 0,023=6,67⋅10−3\frac{0{,}02}3=6{,}67\cdot10^{-3}.
  • Disuguaglianze: 1κ⋅5⋅10−3=1,1⋅10−3≤6,67⋅10−3≤κ⋅5⋅10−3=2,2⋅10−2\frac1\kappa\cdot5\cdot10^{-3}=1{,}1\cdot10^{-3}\le6{,}67\cdot10^{-3}\le\kappa\cdot5\cdot10^{-3}=2{,}2\cdot10^{-2} ✓.

Residuo ed errore relativi sono dello stesso ordine (0,5%0{,}5\% e 0,7%0{,}7\%), come ci si aspetta con κ≈4\kappa\approx4.

Sistema II.

  • Ax^=(1⋅3+2⋅0, 2⋅3+4,002⋅0)=(3; 6)A\hat{\mathbf x}=(1\cdot3+2\cdot0,\ 2\cdot3+4{,}002\cdot0)=(3;\ 6).
  • Residuo: r=(3−3, 6,002−6)=(0; 0,002)\mathbf r=(3-3,\ 6{,}002-6)=(0;\ 0{,}002), ∥r∥∞=0,002\|\mathbf r\|_\infty=0{,}002, residuo relativo 0,0026,002=3,33⋅10−4\frac{0{,}002}{6{,}002}=3{,}33\cdot10^{-4}.
  • Errore: x−x^=(−2; 1)\mathbf x-\hat{\mathbf x}=(-2;\ 1), ∥⋅∥∞=2\|\cdot\|_\infty=2, errore relativo 21=2\frac21=2 (200%).
  • Disuguaglianze: 1κ⋅3,33⋅10−4=1,85⋅10−8≤2≤κ⋅3,33⋅10−4=6,0\frac1\kappa\cdot3{,}33\cdot10^{-4}=1{,}85\cdot10^{-8}\le2\le\kappa\cdot3{,}33\cdot10^{-4}=6{,}0 ✓.

Il residuo relativo è 3,3⋅10−43{,}3\cdot10^{-4} (cioè «tre cifre giuste»), ma la soluzione è completamente sbagliata. La stima superiore 6,06{,}0 è compatibile con 22; quella inferiore è inutile.

Spiegazione geometrica (sistema II). x^\hat{\mathbf x} soddisfa quasi perfettamente la prima equazione e la seconda per 0,0020{,}002; ma le due rette sono quasi parallele, e ci si può spostare lungo di esse di molto (x^\hat{\mathbf x} è a distanza 22 da x\mathbf x) con un residuo minuscolo. In termini di A−1A^{-1}: e=A−1r=(2001−1000−1000500)(00,002)=(−2; 1)\mathbf e=A^{-1}\mathbf r=\begin{pmatrix}2001&-1000\\-1000&500\end{pmatrix}\begin{pmatrix}0\\0{,}002\end{pmatrix}=(-2;\ 1) ✓, perché A−1A^{-1} moltiplica il residuo per fattori dell'ordine di 10310^3.

4. Criterio d'arresto. Per il sistema I vale: il residuo relativo ≤tol\le\text{tol} garantisce un errore relativo ≤κ⋅tol≈4,4 tol\le\kappa\cdot\text{tol}\approx4{,}4\,\text{tol}, quindi il criterio è affidabile. Per il sistema II no: per garantire un errore relativo 10−610^{-6} servirebbe un residuo relativo ≤10−618012=5,6⋅10−11\le\frac{10^{-6}}{18012}=5{,}6\cdot10^{-11}. Nei metodi iterativi (Metodi iterativi - Jacobi, Gauss-Seidel e convergenzaSi spezza $A=M-N$ con $M$ facile da invertire e si itera $M\mathbf x^{(k+1)}=N\mathbf x^{(k)}+\mathbf b$, cioè $\mathbf x^{(k+1)}=B\mathbf x^{(k)}+\mathbf c$ con $B=M^{-1}N=I-M^{-1}A$. Con $A=D+L+U$ (diagonale, parte strettamente inferiore e superiore): Jacobi $M=D$, $x_i^{(k+1)}=\frac1{a_{ii}}\big(b_i-\sum_{j\ne i}a_{ij}x_j^{(k)}\big)$; Gauss-Seidel $M=D+L$, usa subito i valori nuovi $x_j^{(k+1)}$, $j<i$. Convergenza per ogni $\mathbf x^{(0)}$ se e solo se $\rho(B)<1$ (D); basta $|B|<1$. Matrice a dominanza diagonale stretta per righe: convergono entrambi; $A$ SPD: converge Gauss-Seidel; per matrici tridiagonali $\rho(B_{GS})=\rho(B_J)^2$. Errore: $\mathbf e^{(k)}=B^k\mathbf e^{(0)}$, $|\mathbf e^{(k)}|\le|B|^k|\mathbf e^{(0)}|$ e $|\mathbf x-\mathbf x^{(k)}|\le\frac{|B|}{1-|B|}|\mathbf x^{(k)}-\mathbf x^{(k-1)}|$. Per guadagnare $d$ cifre servono $k\approx\frac d{-\log_{10}\rho(B)}$ iterazioni.Metodi iterativi - Jacobi, Gauss-Seidel e convergenza →) si usa spesso il residuo come criterio d'arresto, e questo risultato mostra che bisogna tenere conto di κ\kappa.

Risultati: sistema I residuo relativo 5⋅10−35\cdot10^{-3}, errore relativo 6,7⋅10−36{,}7\cdot10^{-3}; sistema II residuo relativo 3,3⋅10−43{,}3\cdot10^{-4}, errore relativo 22 (200%).

Errori tipici. Pensare che un residuo piccolo implichi un errore piccolo; calcolare l'errore con la norma 2 e il residuo con la norma ∞\infty nello stesso rapporto; dimenticare di dividere per ∥b∥\|\mathbf b\| e per ∥x∥\|\mathbf x\|.

Versione ripasso

(Norme di matrici e condizionamento dei sistemi lineariNorme di vettore $|x|1=\sum|x_i|$, $|x|2=\sqrt{\sum x_i^2}$, $|x|\infty=\max|x_i|$. Norma di matrice indotta $|A|=\max{x\ne0}\frac{|Ax|}{|x|}$: $|A|1$ = massima somma per colonne, $|A|\infty$ = massima somma per righe, $|A|2=\sqrt{\lambda{\max}(A^TA)}$ (= $\rho(A)$ se $A$ simmetrica); valgono $|Ax|\le|A||x|$, $|AB|\le|A||B|$, $\rho(A)\le|A|$. Numero di condizionamento $\kappa(A)=|A||A^{-1}|\ge1$ ($\kappa_2=\frac{\sigma_{\max}}{\sigma_{\min}}$, per SPD $\frac{\lambda_{\max}}{\lambda_{\min}}$, per ortogonali 1). Se $A\mathbf x=\mathbf b$ e $A(\mathbf x+\delta\mathbf x)=\mathbf b+\delta\mathbf b$: $\frac{|\delta\mathbf x|}{|\mathbf x|}\le\kappa(A)\frac{|\delta\mathbf b|}{|\mathbf b|}$ (D); con residuo $\mathbf r=\mathbf b-A\hat{\mathbf x}$: $\frac1\kappa\frac{|\mathbf r|}{|\mathbf b|}\le\frac{|\mathbf x-\hat{\mathbf x}|}{|\mathbf x|}\le\kappa\frac{|\mathbf r|}{|\mathbf b|}$. Si perdono circa $\log_{10}\kappa$ cifre; un residuo piccolo non garantisce un errore piccolo se $\kappa$ è grande.Norme di matrici e condizionamento dei sistemi lineari →)

Disuguaglianza: 1κ∥r∥∥b∥≤∥e∥∥x∥≤κ∥r∥∥b∥\frac1\kappa\frac{\|\mathbf r\|}{\|\mathbf b\|}\le\frac{\|\mathbf e\|}{\|\mathbf x\|}\le\kappa\frac{\|\mathbf r\|}{\|\mathbf b\|}, e=A−1r\mathbf e=A^{-1}\mathbf r.

  • I (MM SPD tridiagonale, κ∞=4,44\kappa_\infty=4{,}44): Mx^=(6,02;9,95;7,98)M\hat{\mathbf x}=(6{,}02;9{,}95;7{,}98), r=(−0,02;0,05;0,02)\mathbf r=(-0{,}02;0{,}05;0{,}02), residuo relativo 5⋅10−35\cdot10^{-3}; e=(−0,01;0,02;0)\mathbf e=(-0{,}01;0{,}02;0), errore relativo 6,67⋅10−36{,}67\cdot10^{-3}; limiti 1,1⋅10−31{,}1\cdot10^{-3} e 2,2⋅10−22{,}2\cdot10^{-2}.
  • II (A=(1224,002)A=\begin{pmatrix}1&2\\2&4{,}002\end{pmatrix}, κ∞=18012\kappa_\infty=18012): Ax^=(3;6)A\hat{\mathbf x}=(3;6), r=(0;0,002)\mathbf r=(0;0{,}002), residuo relativo 3,33⋅10−43{,}33\cdot10^{-4}; e=(−2;1)\mathbf e=(-2;1), errore relativo 22; limiti 1,85⋅10−81{,}85\cdot10^{-8} e 6,06{,}0. Verifica A−1r=(−2;1)A^{-1}\mathbf r=(-2;1).
  • Criterio d'arresto sul residuo: affidabile con κ\kappa piccolo; con κ=18012\kappa=18012 per errore 10−610^{-6} serve residuo 5,6⋅10−115{,}6\cdot10^{-11}.

Errori tipici: residuo piccolo ⇒\Rightarrow errore piccolo; norme diverse nello stesso rapporto; mancano ∥b∥\|\mathbf b\|, ∥x∥\|\mathbf x\|.

Teoria collegata