% Copyright (C) 2007 Alvise Sommariva. This is a free document; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% at your option) any later version. This document is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA}

\documentclass[10pt]{elsarticle}

\input{myprefs.tex}

\begin{document}
\begin{frontmatter}

\title{\bf Algebra lineare numerica \thanks{Ultima revisione: \today}}

\author{A. Sommariva\thanks{Dipartimento di Matematica,
Universit\mya degli Studi di  Padova, stanza 419, via Trieste 63, 35121 Padova, Italia ({\tt 
alvise@math.unipd.it}).} %Telefono: +39-049-8271350.
}

\begin{keyword}
Norme di vettori e matrici, condizionamento di matrici e sistemi; metodi 
diretti: metodo di eliminazione di Gauss con pivoting e fattorizzazione 
LU, costo computazionale, calcolo del determinante, calcolo della matrice 
inversa, effetto dell'arrotondamento sull'accuratezza della fattorizzazione 
e sulla soluzione di sistemi (matrice di Hilbert); soluzione ai minimi 
quadrati di sistemi sovradeterminati e fattorizzazione QR; cenni ai metodi 
iterativi di tipo splitting, il metodo di Jacobi.
\end{keyword}

\end{frontmatter}

\noindent {\it{Revisione}}: \today

\noindent\rule{12.15cm}{0.4pt}





\section{Condizionamento di matrici e sistemi}

\begin{definition}
Dato uno spazio vettoriale $V$, la funzione $\|\cdot\|: V \rightarrow \R$ \mye una {\sl{norma vettoriale}} se
\begin{enumerate}
\item $\|x\| \geq 0$ per ogni $x \in V$ e $\|x\| = 0$ se e solo se $x=0_V$,
\item $\|\alpha x\| = |\alpha| \| x\|$ per ogni $x \in V$ e scalare $\alpha \in \R$,
\item vale la cosidetta disuguaglianza triangolare $$\|x+y\| \geq \|x\|+\|y\|$$ per ogni $x,y \in V$.
\end{enumerate}
\end{definition}

\begin{example}
Se $V=\C^n$ e $x=(x_1,x_2,\ldots,x_n) \in \C^n$ alcuni esempi di norme sono
\begin{itemize}
\item per $p \in [1,+\infty)$, la {\sl{norma}} $p$ \mye definita come $$\|x\|_p=\left( \sum_{j=1}^n |x_j|^p \right)^{1/p};$$ 
\item la {\sl{norma del massimo}} \mye definita come $$\|x\|_p=\max_{j=1,\ldots,n} |x_j|;$$ 
\end{itemize}
\end{example}

\begin{remark}
Lo stesso tipo di norme si definisce nel caso particolare i vettori appartengano a $\R^{n}$.
\end{remark}

\begin{remark}
Si osservi che
\begin{itemize}
\item per $p=1$, si definisce la importante {\sl{norma}} $1$
$$\|x\|_1=\sum_{j=1}^n |x_j| ;$$ 
\item per $p=2$, si definisce la importante {\sl{norma}} $2$, detta anche  {\sl{euclidea}}
$$\|x\|_2=\left( \sum_{j=1}^n |x_j|^2 \right)^{1/2}.$$ 
\end{itemize}
\end{remark}

\begin{example}
Si consideri il vettore $v=(4,-3) \in \R^{2}$.
Allora
\begin{itemize}
\item $\|v\|_1=|4|+|-3|=7$,
\item $\|v\|_2=\sqrt{|4|^2+|-3|^2}=5$,
\item $\|v\|_{\infty}=\max{|4|,|-3|}=4$.
\end{itemize}
\end{example}

\begin{example}
Si considerino i vettori $u=(4,-3), v=(4.001,-2.999) \in \R^{2}$.

I due vettori sono numericamente {\sl{vicini}}. Per verificarlo calcoliamo la loro distanza in termine di norme ovvero
$$
{\mbox{dist}}(u,v):=\|u-v\|_p, \,\,\, p=[1,\infty].
$$

Allora, visto che $u-v=[-0.001   -0.001]$
\begin{itemize}
\item $\|u-v\|_1=0.002$,
\item $\|u-v\|_2=0.0014$,
\item $\|u-v\|_{\infty}=\max{|4|,|-3|}=0.001$,
\end{itemize}
risultati che sottolineano che $u$ e $v$ sono vicini rispetto alle distanze definite dalle norme.
\end{example}

\begin{definition}
Sia $\| \cdot \|$ una norma vettoriale. Si definisce {\sl{norma vettoriale indotta}} di una matrice in $\C^{n \times n}$ la funzione che mappa ogni matrice $A \in \R^{n \times n}$ nel valore
$$
\|A\|=\sup_{x \neq 0} \frac{\|Ax\|}{\|x\|}.
$$
\end{definition}

\begin{remark}
Lo stesso tipo di norme si definisce nel caso particolare le matrici siano in $\R^{n \times n}$.
\end{remark}

\begin{definition}
Si definisce {\sl{raggio spettrale}} di $A \in \C^{n \times n}$ la quantit\mya $$\rho(A)=\max_{k=1,\ldots,n} (|\lambda_k|)$$ dove $\lambda_k$ \mye un {\sl{autovalore}} di $A$, ovvero esiste $x_k \neq 0$ detto {\sl{autovettore}} di $A$, tale che $Ax_k=\lambda_k x_k$.
\end{definition}

Alcuni esempi di norme matriciali risultano:
\begin{itemize}
\item per $p=1$, partendo dalla norma vettoriale $\|\cdot\|_1$, si dimostra che la {\sl{norma}} $1$ matriciale risulta uguale a 
$$\|A\|_1=\max_{1 \leq j \leq n} \sum_{i=1}^n |a_{i,j}|,$$ dove $A=(a_{i,j})$,
\item per $p=2$, partendo dalla norma vettoriale $\|\cdot\|_2$, si dimostra che la {\sl{norma}} $2$ matriciale risulta uguale a 
$$\|x\|_2=\sqrt{\rho(A^* A)}$$  dove $A=(a_{i,j})$ e se $A \in \C^{n \times n}$ allora $A^*$ \mye la trasposta coniugata di $A$ mentre se $A \in \R^{n \times n}$ allora $A^*$ \mye la trasposta di $A$;
\item per $p=\infty$, partendo dalla norma vettoriale $\|\cdot\|_{\infty}$, si dimostra che la {\sl{norma}} $\infty$ matriciale risulta uguale a 
$$\|A\|_{\infty}=\max_{1 \leq i \leq n} \sum_{j=1}^n |a_{i,j}|,$$ dove $A=(a_{i,j})$.
\end{itemize}

\begin{example}
Si consideri la matrice $A \in \R^{2 \times 2}$
$$
A=\left[
\begin{array}{cc}
$1$ & $-2$ \\
$-3$ & $4$ \\
\end{array}
\right]
$$
Allora 
\begin{itemize}
\item $\|A\|_1=\max(|1|+|-3|,|-2|+|4|)=6$;
\item osservato che $A^*=A^T$ abbiamo
$$
A^T A=\left[
\begin{array}{cc}
$1$ & $-3$ \\
$-2$ & $4$ \\
\end{array} 
\right]
\left[
\begin{array}{cc}
$1$ & $-2$ \\
$-3$ & $4$ \\
\end{array}
\right]
=
\left[
\begin{array}{cc}
$10$ & $-14$ \\
$-14$ & $20$ \\
\end{array}
\right]
$$
i cui autovalori valgono rispettivamente 
$$
\lambda_1 \approx 0,1339312526814949, \,\,\,\, \lambda_2 \approx 29.86606874731851
$$
e di conseguenza 
$$
\|A\|_2=\sqrt{\max(|\lambda_1|,|\lambda_2|)} \approx 5.464985704219044;
$$
\item $\|A\|_{\infty}=\max(|1|+|-2|,|-3|+|4|)=7$.
\end{itemize}
\end{example}

\begin{example}
Si considerino le matrici 
$$
A=\left[
\begin{array}{cc}
$1.001$ & $-2.999$ \\
$-2.002$ & $4.005$ \\
\end{array} 
\right]
{\hspace{0.5cm}}
B=\left[
\begin{array}{cc}
$1$ & $-3$ \\
$-2$ & $4$ \\
\end{array} 
\right]
$$

Le due matrici sono numericamente {\sl{vicine}}. Per verificarlo calcoliamo la loro distanza in termine di norme ovvero
$$
{\mbox{dist}}(A,B):=\|A-B\|_p, \,\,\, p=[1,\infty].
$$

Allora, visto che 
$$
A-B=\left[
\begin{array}{cc}
$0.0010$&$0.0010$ \\
$-0.0020$&$0.0050$\\
\end{array} 
\right]
$$
ricaviamo
\begin{itemize}
\item $\|A-B\|_1=0.006$,
\item $\|A-B\|_2=0.005415654779058332$,
\item $\|A-B\|_{\infty}=\max{(|4|,|-3|)}=0.007$,
\end{itemize}
risultati che sottolineano che $A$ e $B$ sono matrici vicine rispetto alle distanze definite dalle norme.
\end{example}

\begin{theorem}
Per le norme matriciali indotte, valgono le seguenti disuguaglianze
\begin{equation}\label{dis1}
\|Ax\| \leq \|A\|\|x\|,
\end{equation}
e
\begin{equation}\label{dis2}
\|AB\| \leq \|A\|\|B\|.
\end{equation}
\end{theorem}

\begin{proof}
Mostriamo la prima disuguaglianza. 
E' banalmente verificata per $x=0$.
Altrimenti, essendo
$$
\|A\|=\sup_{y \neq 0} \frac{\|Ay\|}{\|y\|} \geq \frac{\|Ax\|}{\|x\|}
$$
abbiamo che per $x \neq 0$,
$$
\|Ax\| \leq \|A\|\|x\|.
$$

\smsp

Per quanto riguarda la seguente disuguaglianza, essendo per il primo punto, se $y \neq 0$
$$
\|ABy\|=\|A(By)\| \leq \|A\| \|By\| \leq \|A\| \|B\|\|y\| 
$$
da cui 
$$
\frac{\|ABy\|}{\|y\|} \leq \|A\| \|B\|
$$
e quindi 
$$
\|AB\|=\sup_{y \neq 0} \frac{\|ABy\|}{\|y\|}  \leq \|A\| \|B\|.
$$
\end{proof}




\begin{definition}
Sia $A \in \R^{n \times n}$ una matrice invertibile ossia tale che $\det(A) \neq 0$. Si dice {\sl{indice di condizionamento}} della matrice $A$, relativamente alla norma indotta $\|\cdot\|$,
$$
k(A):=\|A\|\|A^{-1}\|.
$$
\end{definition}

\begin{remark}
Si osservi che dalla definizione, se $I $ \mye la matrice identica di $\R^{n \times n}$
$$
\|I\|=\sup_{x \neq 0} \frac{\|Ix\|}{\|x\|}=\sup_{x \neq 0} \frac{\|x\|}{\|x\|}=1
$$
e quindi
$$
1=\|I\|=\|A A^{-1}\|\leq \|A\| \|A^{-1}\|=k(A)
$$
ovvero $k(A) \geq 1$.
\end{remark}

Vale il seguente risultato

\begin{theorem}
Qualsiasi sia la norma matriciale indotta, e $A$ invertibile,
$$
1 \leq \frac{|\lambda_{\mbox{max}}|}{|\lambda_{\mbox{min}}|} \leq k(A)
$$
dove $\lambda_{\mbox{min}}$ e $\lambda_{\mbox{max}}$ sono rispettivamente gli autovalori di $A$ di minimo e massimo modulo.

Se $A$ \mye simmetrica 
$$
k_2(A)=\frac{|\lambda_{\mbox{max}}|}{|\lambda_{\mbox{min}}|}.
$$
\end{theorem}

\begin{example}
La matrice $H=(h_{i,j})_{i,j=1,\ldots,n}$ con $h_{i,j}=\frac{1}{i+j-1}$ \mye nota come matrice di Hilbert di ordine $n$. Per queste matrici gli indici di condizionamento crescono molto rapidamente, come si pu\myo vedere dalla tabella.

 \begin{table}[!ht]
\begin{center}
\begin{tabular}{| c c  | c c | c c |}
\hline
$n$&   $c_n$ & $n$ &   $c_n$ & $n$ &   $c_n$   \\
\hline
 	 $ 1$ & $1.00e+00$ & $ 5$ & $4.77e+05$ & $ 9$ & $4.93e+11$ \\
 	 $ 2$ & $1.93e+01$ & $ 6$ & $1.50e+07$ & $10$ & $1.60e+13$ \\
 	 $ 3$ & $5.24e+02$ & $ 7$ & $4.75e+08$ & $11$ & $5.22e+14$ \\
 	 $ 4$ & $1.55e+04$ & $ 8$ & $1.53e+10$ & $12$ & $1.62e+16$ \\
  \hline
\end{tabular}
\caption{Indice di condizionamento $c_n={\mbox{cond}}_2(H_n)$ della matrici di Hilbert $H_n$ di ordine $n$}
\end{center}
\end{table}


\end{example}


Vale la seguente risposta della soluzione di un sistema lineare agli errori del vettore termine noto.

\begin{theorem}
Se 
\begin{enumerate}
\item $A \in \R^{n \times n}$ \mye una matrice invertibile, 
\item $b \in \R^{n}$, $b \neq 0$,
\item $Ax=b$,
\item $A(x+\delta x)=b+\delta b$,
\end{enumerate}
allora vale la stima
$$
\frac {\|\delta x\|}  {\|x\|} \leq k(A) \frac {\|\delta b\|}  {\|b\|}.
$$
\end{theorem}

\begin{proof}
Da $A(x+\delta x)=b+\delta b$ e $Ax=b$ abbiamo che
$$
A\delta x=A(x+\delta x)-Ax=(b+\delta b)-b
$$
ovvero $A\delta x=\delta b$ e quindi $\delta x=A^{-1}\delta b$.
Ne consegue che
\begin{equation}\label{theo1.1}
\|\delta x\|=\|A^{-1}\delta b\| \leq \|A^{-1}\| \|\delta b\|.
\end{equation}
D'altra parte, essendo $Ax=b$
$$
\|b\|=\|Ax\| \leq \|A\|\|x\|
$$
e quindi, essendo $x \neq 0$ in quanto $b \neq 0$, ricaviamo
\begin{equation}\label{theo1.2}
\frac{1}{\|x\|} \leq \frac{\|A\|}{\|b\|}.
\end{equation}
Da (\ref{theo1.1}), (\ref{theo1.2}), ricaviamo
$$
\frac {\|\delta x\|}{\|x\|}\leq \|A^{-1}\| \|A\|   \frac{\|\delta b\|}{\|b\|}=k(A) \frac{\|\delta b\|}{\|b\|}.
$$
\end{proof}

\begin{comment}
Questo teorema mostra che si fa un errore relativo sul termine noto pari a $\frac{\|\delta b\|}{\|b\|}$ allora si compie un errore relativo sulla soluzione $\frac {\|\delta x\|}{\|x\|}$ tale che
$$
\frac {\|\delta x\|}{\|x\|}\leq k(A) \frac{\|\delta b\|}{\|b\|}.
$$
Quindi, se $k(A) \gg 1$, pu\myo accadere che la soluzione perturbata $x+\delta x$ sia molto distante da $x$ e di conseguenza la soluzione del sistema lineare $Ax=b$ molto suscettibile a piccoli errori sul termine noto. In altri termini pi\myu grande \mye $k(A)$ pi\myu i sistemi lineari $Ax=b$ sono difficili da trattare nel senso che piccoli errori sui dati $b$ possono portare a grossi errori sulla soluzione $x$.
\end{comment}

\begin{example}
Siano
$$
A=\left[
\begin{array}{cc}
$7$&$10$ \\
$5$&$7$\\
\end{array} 
\right],
\hspace{0.25cm}
b=\left[
\begin{array}{c}
$1$ \\
$0.7$\\
\end{array} 
\right],
\hspace{0.25cm}
b+\delta b=\left[
\begin{array}{c}
$1.01$ \\
$0.69$\\
\end{array} 
\right].
$$
Il sistema $Ax=b$ ha soluzione
$$
x=\left[
\begin{array}{c}
$0$ \\
$0.1$\\
\end{array} 
\right].
$$
mentre $Ax=b+\delta b$ ha soluzione
$$
x+\delta x=\left[
\begin{array}{c}
$-0.17$ \\
$0.22$\\
\end{array} 
\right].
$$
Di conseguenza seppure $\frac{\|\delta b\|_{\infty}}{\|b\|_{\infty}}=0.01$ abbiamo $\frac{\|\delta x\|_{\infty}}{\|x\|_{\infty}}=1.7$. Questo significa che a un piccolo errore relativo sui dati abbiamo riscontrato un eccessivo errore relativo sulla soluzione.

Osserviamo che in questo caso $k(A)=\|A\|_{\infty}\|A^{-1}\|_{\infty}=289$ e quindi non troppo vicino a $1$.
\end{example}

Mostriamo ora la risposta della soluzione di un sistema lineare agli errori sulla matrice.

\begin{theorem}
Se 
\begin{enumerate}
\item $A \in \R^{n \times n}$ \mye una matrice invertibile, 
\item $b \in \R^{n}$, $b \neq 0$,
\item $Ax=b$,
\item $(A+\delta A)(x+\delta x)=b$, con $\det(A+\delta A) \neq 0$,
\end{enumerate}
allora vale la stima
$$
\frac {\|\delta x\|}  {\|x+\delta x\|} \leq k(A) \frac {\|\delta A\|}  {\|A\|}.
$$
\end{theorem}

\begin{proof}
Essendo $(A+\delta A)(x+\delta x)=b$ e $Ax=b$, necessariamente
\begin{eqnarray}
0= (A+\delta A)(x+\delta x)-b&=& Ax+A\delta x+\delta A x +\delta A \delta x -b \nonumber \\
&=&A\delta x+\delta A x +\delta A \delta x=A\delta x+\delta A (x+\delta x) \nonumber
\end{eqnarray}
da cui
$$
A \delta x=-\delta A (x + \delta x)
$$
ovvero
$$
\delta x=-A^{-1}\delta A (x + \delta x)
$$
e di conseguenza 
$$
\|\delta x\|=\|-A^{-1}\delta A (x + \delta x)\| \leq \|A^{-1}\|\|\delta A\|\| x + \delta x\|.
$$
Quindi, essendo $\|A\| \neq 0$,
$$
\frac {\|\delta x\|}{\| x + \delta x\|} 
\leq \|A^{-1}\|\|\delta A\|=
\frac {\|A^{-1}\| \|A\|} {\|A\|}   \|\delta A\| 
$$
da cui
$$
\frac {\|\delta x\|} {\| x + \delta x\|}  \leq k(A) \frac {\|\delta A\|}{ \|A\|}.
$$
\end{proof}

Vogliamo ora studiare il caso in cui tanto la matrice quanto il termine noto siano soggetti ad errori.

Per prima cosa introduciamo il seguente importante teorema.

\begin{theorem}\label{thm3lemma}
Se
\begin{enumerate}
\item $\|\cdot \|$ \mye una norma matriciale indotta da una vettoriale,
\item $\|A\| < 1$,
\end{enumerate}
allora "$I-A$" e "$I+A$" sono invertibili e si ha
$$
\|(I-A)^{-1}\| \leq \frac{1}{1-\|A\|},
$$
$$
\|(I+A)^{-1}\| \leq \frac{1}{1-\|A\|}.
$$
\end{theorem}

Finalmente,

\begin{theorem}\label{thm_3}
Se 
\begin{enumerate}
\item $A \in \R^{n \times n}$ \mye una matrice invertibile, 
\item $b \in \R^{n}$, $b \neq 0$,
\item $Ax=b$,
\item $(A+\delta A)(x+\delta x)=b+\delta b$, con $\det(A+\delta A) \neq 0$,
\item $k(A) \|\delta A\| < \|A\|$,
\end{enumerate}
allora vale la stima
$$
\frac {\|\delta x\|}  {\|x\|} \leq \frac{k(A)} { 1-k(A) \|\delta A\|/\|A\|}    \left( \frac {\|\delta A\|}  {\|A\|} +\frac {\|\delta b\|}  {\|b\|} \right).
$$
\end{theorem}

\begin{comment}
Il teorema stima l'errore relativo $\frac {\|\delta x\|}  {\|x\|}$ relativamente alle soluzioni, in funzione degli errori relativi sulla perturbazione della matrice $\frac {\|\delta A\|}  {\|A\|}$ e sul termine noto $\frac {\|\delta b\|}  {\|b\|}$.
Il valore 
$$
\frac{k(A)} { 1-k(A) \|\delta A\|/\|A\|}
$$
\mye la costante di amplificazione degli errori.

Visto che il denominatore \mye maggiore di $1$ visto che $k(A) \|\delta A\|/\|A\| < 1$, deduciamo che se $k(A) \gg 1$, la soluzione del sistema perturbato $(A+\delta A)(x+\delta x)=b+\delta b$ pu\myo essere molto distante da quella di $Ax=b$.
\end{comment}

\begin{proof}
Essendo $(A+\delta A)(x+\delta x)=b+\delta x$ e $Ax=b$ deduciamo 
\begin{eqnarray}
b+\delta b&=& (A+\delta A)(x+\delta x) =Ax+A\delta x+\delta A x +\delta A\delta x \nonumber \\
&=& b+A\delta x+\delta A x +\delta A\delta x \nonumber 
\end{eqnarray}
da cui sottraendo $b$ ad ambo i membri
$$
\delta b=A \delta x+\delta A x +\delta A \delta x=\delta A x + (A+\delta A) \delta x
$$ 
ovvero
\begin{equation}\label{thm3.1}
(A+\delta A) \delta x=\delta b-\delta A x.
\end{equation}
Posto $B=A^{-1} \delta A$, ricaviamo 
$$
A+\delta A=A+AA^{-1}\delta A=A(I+A^{-1}\delta A)=A(I+B)
$$
e da $k(A) \|\delta A\| < \|A\|$
\begin{eqnarray}
\|B\|&=&\|A^{-1} \delta A\| \leq \|A^{-1}\|\| \delta A\| \nonumber \\
&=&\frac{ \|A^{-1}\| \|A\|} {\|A\|}\| \delta A\|=k(A) \frac{\| \delta A\|}{\|  A\|} \nonumber \\
& \leq & \frac{\|  A\|}{\|  A\|} =1
\end{eqnarray}
ovvero $
\|B\| < 1.
$

Osserviamo che la matrice $A+\delta A$ \mye invertibile in quanto lo \mye $A$, lo \mye $I+B$ per il teorema \ref{thm3lemma} e da
$$
A+\delta A=A(I+A^{-1}\delta A)=A(I+B)
$$
lo \mye pure $A+\delta A$ perch\mye il prodotto di matrici invertibili \mye invertibile (ricordare che $\det(MN)=\det(M) \det(N)$, e quindi se $\det(M), \det(N) \neq 0$ allora $\det(MN) \neq 0$).

Dal teorema \ref{thm3.1}, moltiplicando ambo i membri per $(A+\delta A)^{-1}$, visto che $(MN)^{-1}=N^{-1}M^{-1}$ ricaviamo
\begin{eqnarray}
\delta x &=&(A+\delta A)^{-1} (\delta b-\delta A x)=(A(I+B))^{-1} (\delta b-\delta A x) \nonumber \\
&=& (I+B)^{-1} A^{-1}  (\delta b-\delta A x) \nonumber \\
\end{eqnarray}
da cui, essendo nuovamente per il teorema \ref{thm3.1} $\| (I+B)^{-1}\| \leq 1/(1-\|B\|)$, $B=A^{-1} \delta A$
\begin{eqnarray}
\|\delta x\| &=& \| (I+B)^{-1} A^{-1}  (\delta b-\delta A x) \| \nonumber \\
&\leq& \| (I+B)^{-1}\| \| A^{-1} \| \| \delta b-\delta A x \| \nonumber \\
&\leq& \frac{1}{1-\|B\|} \| A^{-1} \| \| \delta b-\delta A x \|  \nonumber \\
&=&  \frac{\| A^{-1} \|}{1-\|A^{-1} \delta A\|}  \| \delta b-\delta A x \|  \nonumber \\
&\leq &  \frac{\| A^{-1} \|}{1-\|A^{-1}\| \| \delta A\|} ( \| \delta b\|+\|\delta A\|\| x \|)
\end{eqnarray}
e quindi riassumendo
\begin{equation}\label{thm3_4}
\|\delta x\| \leq  \frac{\| A^{-1} \|}{1-\|A^{-1}\| \| \delta A\|} ( \| \delta b\|+\|\delta A\|\| x \|).
\end{equation}
Essendo $b=Ax$ abbiamo $\|b\|=\|Ax\| \leq \|A\| \|x\|$ da cui
\begin{equation}\label{thm3_5}
\frac{1}{\|A\|\|x\|} \leq \frac{1}{\|b\|}
\end{equation}
e moltiplicando membro a membro i risultati di (\ref{thm3_4}), (\ref{thm3_5}), 
\begin{eqnarray}\label{thm3_6}
\frac{\|\delta x\|}{\|x\|} 
&\leq&     \frac{\| A^{-1} \|}{1-\|A^{-1}\| \| \delta A\|} \left(\frac{ \| \delta b\|}{\|  x \|}+\|\delta A\|\| x \| \frac{ 1}{\|  x \|} \right) \nonumber \\
&\leq&     \frac{ \| A^{-1} \|  \| A \|   }{1-\|A^{-1}\| \| \delta A\|}  \frac{ \| 1 \|}{  \| A \|   }   \left(\frac{ \| \delta b\|}{\|  x \|}+\|\delta A\| \right) \nonumber \\
&\leq&     \frac{ k(A)  }{1-\|A^{-1}\| \|A\| \| \frac{\|\delta A\|}{\|A\|}}     \left(\frac{ \| \delta b\|}{\|  A \|\|  x \|}+\frac{\|\delta A\|}{\|A\|} \right) \nonumber \\
&\leq&     \frac{ k(A)  }{1-k(A) \| \frac{\|\delta A\|}{\|A\|}}     \left(\frac{ \| \delta b\|}{\|b\|}+\frac{\|\delta A\|}{\|A\|} \right).
\end{eqnarray}

\end{proof}

\begin{example}
Il sistema $Ax=b$, dove 
$$
A=\left[
\begin{array}{cc}
$1$&$2$ \\
$0.499$&$1.001$\\
\end{array} 
\right],
b=\left[
\begin{array}{c}
$3$ \\
$1.5$\\
\end{array} 
\right],
$$
ha soluzione $x=(1,1)^T$.

Se consideriamo  
$$
A+\delta A=\left[
\begin{array}{cc}
$1$&$2$ \\
$0.5$&$1.0015$\\
\end{array} 
\right],
\hspace{0.25cm}
b+\delta b=\left[
\begin{array}{c}
$3$ \\
$1.4985$\\
\end{array} 
\right].
$$
il sistema 
\begin{itemize}
\item $(A+\delta A)u=b$ ha soluzione $u=(3,0)^T$ (perturbata matrice), 
\item  $Av=b+\delta b$ ha soluzione $v=(2,0.5)^T$  (perturbato termine noto),
\item $(A+\delta A)w=(b+\delta b)$ ha soluzione $w=(5,-1)^T$ (perturbati matrice e termine noto).
\end{itemize}
Ne risulta che per possibili piccole perturbazioni alla matrice e/o al termine noto abbiamo ottenuto nei diversi casi una soluzione molto {\sl{distante}} dall'originaria.

Con riferimento al teorema \ref{thm_3}
\begin{itemize}
\item $\det(A) \approx 0.003$, $\det(\delta A) \approx 0.0015$,
\item $k_2(A) \approx 2083.666853410356$,
\item abbiamo
$$
 \frac{ k_2(A)  }{1-k_2(A) \| \frac{\|\delta A\|_2}{\|A\|_2}} \approx 3.05387e+04,
$$
\item $\frac{ \| \delta x\|_2}{\|x\|_2} \approx 3.16228e+00$,
\item $\frac{ \| \delta b\|_2}{\|b\|_2}  \approx 4.47214e-04$ e $\frac{\|\delta A\|_2}{\|A\|_2} \approx 4.47178e-04$.
\end{itemize}

La stima asserisce in effetti che
\begin{eqnarray}
3.16228e+00 &\leq& 3.05387e+04 \cdot  (4.47214e-04+4.47178e-04) \nonumber \\ &\approx  & 2.73136e+01.
\end{eqnarray}

\end{example}

\section{Metodo di eliminazione di Gauss e fattorizzazione LU}

Il proposito di questa sezione \mye di introdurre il metodo dell'eliminazione gaussiana per risolvere sistemi lineari $Ax=b$ dove supponiamo $A$ invertibile, cio\mye tale che $\det(A) \neq 0$.

Per semplificare l'esposizione cominciamo con un esempio.

\begin{example}
Si consideri il sistema lineare
$$
\left\{
\begin{array}{l}
3x_1-2x_2-x_3=0 \\
6x_1-2x_2+2x_3=6 \\
-9x_1+7x_2+x_3=-1 \\
\end{array}
\right.
$$
Sappiamo che 
\begin{itemize}
\item se moltiplichiamo ambo i membri di una stessa equazione per una costante non nulla, il nuovo problema ha le stesse soluzioni del precedente;
\item se  scambiamo due equazioni,  il nuovo problema ha le stesse soluzioni del precedente;
\item in generale se moltiplichiamo una equazione per una costante e sommiamo il risultato membro a membro a un'altra equazione, otteniamo un nuovo problema che ha le stesse soluzioni del precedente.
\end{itemize}

Diremo che due sistemi sono equivalenti, se hanno le stesse soluzioni.

\begin{itemize}
\item moltiplichiamo la prima riga per $-2$ e la sommiamo alla seconda, ottenendo
 $$
\left\{
\begin{array}{l}
3x_1-2x_2-x_3=0 \\
0 x_1+2x_2+4x_3=6 \\
-9x_1+7x_2+x_3=-1 \\
\end{array}
\right.
$$
\item moltiplichiamo la prima riga per $3$ e la sommiamo alla terza, ottenendo
 $$
\left\{
\begin{array}{l}
3x_1-2x_2-x_3=0 \\
0 x_1+2 x_2+4x_3=6 \\
0 x_1+1 x_2-2x_3=-1 \\
\end{array}
\right.
$$
\item moltiplichiamo la seconda riga per $-1/2$ e la sommiamo alla terza, ottenendo
 $$
\left\{
\begin{array}{l}
3x_1-2x_2-x_3=0 \\
0 x_1+2 x_2+4x_3=6 \\
0 x_1+0 x_2-4x_3=-4 \\
\end{array}
\right.
$$
\end{itemize}
Il sistema finale, \mye facile da risolvere. 
\begin{itemize}
\item L'ultima equazione ha una sola variabile e una sola incognita e comporta che $$x_3=1.$$
\item Inseriamo questo risultato nella penultima equazione e otteniamo 
$$ 
2 x_2+4\cdot 1=6
$$
da cui $2x_2=2$ e quindi $x_2=1$.
\item Inseriamo questi risultati nella prima equazione e otteniamo 
$$
3x_1-2 \cdot 1-1=0
$$
da cui $$x_1=1.$$
\end{itemize}
Il sistema pu\myo essere scritto matricialmente come $Ax=b$ dove
$$
A=\left[
\begin{array}{c c c}
3 & -2 & -1 \\
6 & -2 & 2 \\
-9 & 7 & 1 \\
\end{array}
\right],
\hspace{0.5cm}
b=\left[
\begin{array}{c}
0 \\
6 \\
-1 \\
\end{array}
\right],
$$
e abbiamo mostrato essere equivalente a
$$
U=\left[
\begin{array}{c c c}
3 & -2 & -1 \\
0 & 2 & 4 \\
0 & 0 & -4 \\
\end{array}
\right],
\hspace{0.5cm}
c=\left[
\begin{array}{c}
0 \\
6 \\
-4\\
\end{array}
\right].
$$
Dal punto di vista matriciale, abbiamo generato una sequenza di sistemi lineari equivalenti, ma via via pi\myu {\sl{semplici}}. L'ultimo \mye un sistema $Ux=c$ con $U$ triangolare superiore, ovvero tale che $U_{i,j}=0$ qualora $i > j$, che abbiamo risolto facilmente mediante {\sl{sostituzione all'indietro}}.
\end{example}

\subsection{Sistemi lineari triangolari}

\begin{definition}
Una matrice $A=(a_{i,j}) \in \R^{n \times n}$ si dice 
\begin{itemize}
\item {\sl{triangolare inferiore}} se $A_{i,j}=0$ qualora $i < j$;
\item {\sl{triangolare superiore}} se $A_{i,j}=0$ qualora $i > j$.
\end{itemize}
\end{definition}

Sia $b=(b_j) \in \R^{n}$ e supponiamo di dover risolvere il sistema $Ax=b$. 

In tal caso, 
\begin{itemize}
\item se $A$ \mye triangolare superiore e invertibile, si usa l'algoritmo della {\sl{sostituzione all'indietro}} ovvero
\begin{equation}\label{bs}
\left\{
\begin{array}{l}
x_n=\frac{b_n}{a_{n,n}}\\
x_k=\frac {b_k-(a_{k,k+1}x_{k+1} + \ldots +a_{k,n}x_{n})}  {a_{k,k}}, \,\,\,k=n-1,\ldots,1.
\end{array}
\right.
\end{equation}
Si noti che essendo $A$ invertibile, necessariamente $a_{k,k} \neq 0$ per $k=1,\ldots,n$.
\item se $A$ \mye triangolare inferiore e invertibile, si usa l'algoritmo della {\sl{sostituzione in avanti}} ovvero
\begin{equation}\label{fs}
\left\{
\begin{array}{l}
x_1=b_1/a_{1,1}, \\
x_k=\frac {b_k-(a_{k,1}x_{1} + \ldots +a_{k,k-1}x_{k-1})}  {a_{k,k}}, \,\,\,k=2,\ldots,n.
\end{array}
\right.
\end{equation}
Si noti che essendo $A$ invertibile, necessariamente $a_{k,k} \neq 0$ per $k=1,\ldots,n$.
\end{itemize}

\begin{remark}
Si osservi che la risoluzione di sistemi lineari triangolari necessita approssimativamente $n^2/2$ operazioni moltiplicative, dove $n$ \mye l'ordine della matrice.

Infatti il numero di operazioni moltiplicative, nel caso della sostituzione all'indietro risulta
$$
1+2+\ldots+n = n(n+1)/2 \approx n^2/2.
$$

\end{remark}

\subsection{Metodo di eliminazione gaussiana senza pivoting}

In questa sezione trattiamo la risoluzione di sistemi lineari $Ax=b$ dove in generale si suppone 
\begin{itemize}
\item $A \in \R^{n \times n}$, invertibile e non triangolare, 
\item $b \in \R^{n}$.
\end{itemize}

Posto $A=A^{(1)}=(a^{(1)}_{i,j})$, $b=b^{(1)}=(b^{(1)}_j)$, possiamo riscrivere il sistema nella formula nota come {\sl{aumentata}}

$$
[A^{(1)} | b^{(1)}]=\left[
\begin{array}{c c c c | c}
a^{(1)}_{1,1} & a^{(1)}_{1,2} & \ldots & a^{(1)}_{1,n} & b^{(1)}_{1}   \\
a^{(1)}_{2,1} & a^{(1)}_{2,2} & \ldots & a^{(1)}_{2,n} & b^{(1)}_{2}   \\
\ldots & \ldots & \ldots & \ldots & \ldots   \\
a^{(1)}_{n,1} & a^{(1)}_{n,2} & \ldots & a^{(1)}_{n,n} & b^{(1)}_{n}   \\
\end{array}
\right].
$$

\begin{itemize}
\item {\bf{Primo passo}}. Supponiamo sia $a^{(1)}_{1,1} \neq 0$. Intendiamo definire un sistema $[A^{(2)}|b^{(2)}]$ tale che $A^{(2)}x=b^{(2)}$ sia equivalente a $A^{(1)}x=b^{(1)}$, ovvero con la stessa soluzione, ma in cui  $A^{(2)}=(a^{(2)}_{i,j})$ sia tale che 
$$a^{(2)}_{i,1}=0, \,\, i=2,\ldots,n.$$
A tal proposito definiamo il {\sl{moltiplicatore}}
$$
m_{i,1}=\frac{a^{(1)}_{i,1}}{a^{(1)}_{1,1}}, i=2,\ldots,n
$$
e in successione moltiplichiamo la prima riga per $m_{i,1}$ e la sommiamo all'i-sima, ottenendo
$$
a^{(2)}_{i,j}=a^{(1)}_{i,j}-m_{i,1} a^{(1)}_{1,j}, j=2,\ldots,n,
$$
$$
b^{(2)}_{i}=b^{(1)}_{i}-m_{i,1} b^{(1)}_{i}, i=2,\ldots,n.
$$
Si vede immediatamente che si raggiunge lo scopo prefisso visto che per definizione dei moltiplicatori
$$
a^{(2)}_{i,1}=a^{(1)}_{i,1}-m_{i,1}a^{(1)}_{1,1}=0, \,\,\, i=2,\ldots,n.
$$
Di conseguenza la matrice aumentata $[A^{(2)}|b^{(2)}]$ avr\mya la forma
$$
[A^{(2)} | b^{(2)}]=\left[
\begin{array}{c c c c | c}
a^{(2)}_{1,1} & a^{(2)}_{1,2} & \ldots & a^{(2)}_{1,n} & b^{(2)}_{1}   \\
0 & a^{(2)}_{2,2} & \ldots & a^{(2)}_{2,n} & b^{(2)}_{2}   \\
\ldots & \ldots & \ldots & \ldots & \ldots   \\
0 & a^{(2)}_{n,2} & \ldots & a^{(2)}_{n,n} & b^{(2)}_{n}   \\
\end{array}
\right].
$$
\item {\bf{Passo k-simo}}. Supponiamo $k=2,\ldots,n-1$. Assumiamo che sia 
$$
[A^{(k)} | b^{(k)}]=\left[
\begin{array}{c c c c c c c | c}
a^{(k)}_{1,1} & a^{(k)}_{1,2} & \ldots & \ldots & \ldots & \ldots & a^{(k)}_{1,n} & b^{(k)}_{1}   \\
0 & a^{(k)}_{2,2} & \ldots & \ldots & \ldots & \ldots & a^{(k)}_{2,n} & b^{(k)}_{2}   \\
\ldots & \ddots & \ddots & \ldots & \ldots & \ldots & \ldots & \ldots   \\
0 & \ldots &0 & a^{(k)}_{k,k} & \ldots & \ldots & a^{(k)}_{k,n} & b^{(k)}_{k}   \\
0 & \ldots &0 & a^{(k)}_{k+1,k} & \ldots & \ldots & a^{(k)}_{k+1,n} & b^{(k)}_{k+1}   \\
0 & \ldots & 0 & \ldots & \ldots & \ldots & \ldots & \ldots   \\
0 & \ldots &0 & a^{(k)}_{n,k} & \ldots & \ldots & a^{(k)}_{n,n} & b^{(k)}_{n}   \\
\end{array}
\right].
$$
Se l'{\sl{elemento pivotale}} $a^{(k)}_{k,k} \neq 0$ definiamo
$$
m_{i,k}=\frac {a^{(k)}_{i,k}} {a^{(k)}_{k,k}}, \,\,\, i=k+1,\ldots,n.
$$
e similmente al passo 1, moltiplichiamo la $k$-sima riga per $m_{i,k}$ e la sommiamo alla $i$-sima, $i=k+1,\ldots,n$ ovvero
$$
a^{(k+1)}_{i,j}=a^{(k)}_{i,j}-m_{i,k}a^{(k)}_{k,j},\,\,\, j=k+1,\ldots,n,
$$
$$
b^{(k+1)}_{i}=b^{(k)}_{i}-m_{i,k}b^{(k)}_{k},\,\,\, j=k+1,\ldots,n.
$$
Nuovamente, come nel primo passo, dalla definizione dei moltiplicatori, otteniamo 
$$
a^{(k+1)}_{i,k}=0, \,\,\, i=k+1,\ldots,n,
$$
da cui
$$
[A^{(k+1)} | b^{(k+1)}]=\left[
\begin{array}{c c c c c c c | c}
a^{(k+1)}_{1,1} & a^{(k+1)}_{1,2} & \ldots & \ldots & \ldots & \ldots & a^{(k+1)}_{1,n} & b^{(k+1)}_{1}   \\
0 & a^{(k+1)}_{2,2} & \ldots & \ldots & \ldots & \ldots & a^{(k+1)}_{2,n} & b^{(k+1)}_{2}   \\
\ldots & \ddots & \ddots & \ldots & \ldots & \ldots & \ldots & \ldots   \\
0 & \ldots &0 & a^{(k+1)}_{k,k} & \ldots & \ldots & a^{(k+1)}_{k,n} & b^{(k)}_{k}   \\
0 & \ldots &0 & 0 & a^{(k+1)}_{k+1,k+1} & \ldots & a^{(k)}_{k+1,n} & b^{(k)}_{k+1}   \\
0 & \ldots & 0 & 0 & a^{(k+1)}_{n-1,k+1} & \ldots & \ldots & \ldots   \\
0 & \ldots &0 & 0 &a^{(k+1)}_{n,k+1}  & \ldots & a^{(k)}_{n,n} & b^{(k)}_{n}   \\
\end{array}
\right].
$$
\end{itemize}
Con questo procedimento, dopo $n-1$ passi otteniamo un sistema $A^{(n-1)}x=b^{(n-1)}$ in cui la matrice $A^{(n-1)}$ \mye triangolare superiore con elementi diagonali $A^{(n-1)}_{k,k} \neq 0$, $k=1,\ldots,n$.

Di conseguenza il sistema $A^{(n-1)}x=b^{(n-1)}$ pu\myo essere facilmente risolto con l'algoritmo di sostituzione all'indietro.

\begin{definition}
Una matrice $A \in \R^{n \times n}$ \mye fattorizzabile $LU$ se esistono
\begin{itemize}
\item $L=(l_{i,j}) \in \R^{n \times n}$ triangolare inferiore con elementi diagonali $l_{i,i}=1$, 
\item $u=(l_{i,j}) \in \R^{n \times n}$ triangolare superiore, 
\end{itemize}
tali che $A=LU$ (cf. \cite{WIKILU}).
\end{definition}

\begin{important}
Dall'applicazione della procedura precedente per la risoluzione di un sistema lineare, qualora tutti i pivot $a^{(k)}_{k,k}$ siano non nulli, otteniamo un sistema $A^{(n-1)}x=b^{(n-1)}$ equivalente all'originario e durante il calcolo introduciamo i moltiplicatori $m_{i,k}$, $k=1,\ldots,n-1$, $i=k+1,\ldots,n$. 
Siano 
\begin{itemize}
\item $U=A^{(n-1)}$,
\item $L$ la matrice
$$
L=\left[
\begin{array}{c c c c c}
1 & 0 & 0 & \ldots & 0  \\
m_{2,1} & 1 & 0 & \ldots & 0 \\
\vdots & \vdots & \vdots & \ldots & \vdots \\
m_{n,1} &m_{n,2} & \ldots & \ldots & \ \\
\end{array}
\right].
$$
\end{itemize}
Allora $A=LU$ (cf. {\cite[p.511]{ATK}, \cite[p.143, p.158]{BINI}}).

\smsp

Questo risultato \mye rilevante, perch\mye tale fattorizzazione viene molto utilizzata, come vedremo in seguito con qualche esempio, indipendentemente dalla risoluzione di sistemi lineari.  

\begin{example}
Sia, come nell'esempio iniziale,
$$
A=\left[
\begin{array}{c c c}
3 & -2 & -1 \\
6 & -2 & 2 \\
-9 & 7 & 1 \\
\end{array}
\right],
\hspace{0.5cm}
$$
Alla fine del processo di eliminazione gaussiana senza pivoting abbiamo determinato $U=A^{(n-1)}$  uguale a
$$
U=\left[
\begin{array}{c c c}
3 & -2 & -1 \\
0 & 2 & 4 \\
0 & 0 & -4 \\
\end{array}
\right].
$$
Tenendo conto dei moltiplicatori
$$
L=\left[
\begin{array}{c c c}
1 & 0 & 0 \\
2 & 1 & 0 \\
-3 & 0.5 & 1 \\
\end{array}
\right],
$$
e si verifica facilmente che $A=LU$, con $L,U$ triangolari e della forma richiesta.
\end{example}
\end{important}

\begin{remarkh}
L'algoritmo di eliminazione gaussiana \mye stato scoperto in Cina circa 2000 anni fa, sebbene pi\myu recenti contributi includano i lavori di Lagrange e Jacobi. La metodologia moderna per descrivere questi algoritmi \mye stata introdotta solo nel 1930. {\cite[p.607]{TREF}}.
\end{remarkh}

\subsection{Metodo di eliminazione gaussiana con pivoting (parziale)}

Nella precedente sezione, uno dei punti cruciali \mye risultato che i pivot $a^{(k)}_{k,k}$  fossero non nulli. Questo ci ha permesso di calcolare i moltiplicatori 
$$
m_{i,k}=\frac {a^{(k)}_{i,k}} {a^{(k)}_{k,k}}, \,\,\, k=1,\ldots,n-1, \,\,\, i=k+1,\ldots,n
$$
e di seguito ridurre il problema iniziale alla risoluzione di un sistema $A^{(n-1)}x=b^{(n-1)}$, con $A^{(n-1)}$ triangolare superiore.

Il problema \mye che qualche ${a^{(k)}_{k,k}}$ potrebbe risultare nullo e quindi il processo precedente risulterebbe inapplicabile. 

Ma anche se cos\myi non fosse, potrebbe accadere che ${a^{(k)}_{k,k}}$ sia molto piccolo in modulo, questione che rende il metodo soggetto a una cattiva propagazione degli errori (cf. {\cite[p.516]{ATK}, \cite[p.174]{BINI}}).

La soluzione risulta la seguente, detta del {\sl{pivoting (parziale)}} o del {\sl{massimo pivot}}. 

Al passo $k$-simo dell'eliminazione gaussiana, determiniamo
$$
c_k=\max_{k \leq i \leq n} |a^{(k)}_{i,k}|.
$$
Certamente $c_k \neq 0$, altrimenti, dopo qualche conto, si vedrebbe che la matrice originaria $A$ di cui studiamo $Ax=b$ sarebbe non invertibile, diversamente da quanto richiesto.

Sia $i:=\phi(k)$ il pi\myu piccolo indice di riga, con $i \geq k$, per cui il massimo $c_k$ \mye ottenuto. Se $i > k$, scambiamo la riga $i$-sima con la $k$-sima della matrice aumentata
$[A^{(k)} | b^{(k)}]$ ottenendo una nuova matrice aumentata, diciamo $[{{\hat{A}}}^{(k)} | {{\hat{b}}}^{(k)}]$, su cui possiamo applicare la tecnica introdotta al $k$-simo passo dell'eliminazione gaussiana senza pivoting.



\begin{important}
Dall'applicazione della procedura dell'eliminazione gaussiana con pivoting per la risoluzione di un sistema lineare, abbiamo ottenuto un sistema $A^{(n-1)}x=b^{(n-1)}$ equivalente all'originario.
Poniamo 
\begin{itemize}
\item $P=P_{n-1}\cdot \ldots \cdot P_1$ dove $P_k=((p_k)_{i,j})$ \mye la matrice identica, ad eccezione delle righe $k$ e $\phi(k)$ in cui sono uguali a $1$ esclusivamente le componenti $(p_k)_{k,\phi(k)}$ e $(p_k)_{\phi(k),k}$,
\item $U=A^{(n-1)}$ di forma triangolare superiore,
\item $L=(l_{i,j})$ la matrice
$$
L=\left[
\begin{array}{c c c c c c}
1 & 0 & 0 & \ldots & 0 & 0  \\
m_{2,1} & 1 & 0 & \ldots & 0 & 0 \\
\vdots & \vdots & \vdots & \ldots & \vdots & \vdots \\
m_{n,1} &m_{n,2} & \ldots & \ldots &m_{n,n-1} & 1 \ \\
\end{array}
\right].
$$
 di forma triangolare inferiore, con componenti diagonali $l_{k,k}=1$, $k=1,\ldots,n$.
\end{itemize}
Allora $PA=LU$ (cf. {\cite[p.511]{ATK}, \cite[p.143, p.158]{BINI}}).

\smsp

Questo risultato \mye rilevante, perch\mye tale fattorizzazione viene molto utilizzata, come vedremo in seguito con qualche esempio, indipendentemente dalla risoluzione di sistemi lineari.  
\end{important}

\begin{example}
Si consideri il sistema lineare
$$
\left\{
\begin{array}{l}
3x_1-2x_2-x_3=0 \\
6x_1-2x_2+2x_3=6 \\
-9x_1+7x_2+x_3=-1 \\
\end{array}
\right.
$$

Il termine maggiore in modulo della prima colonna \mye il terzo e vale $9$. Con le notazioni introdotte, $\phi(1)=3$. Scambiamo la prima riga con la terza e ricaviamo il sistema equivalente

$$
\left\{
\begin{array}{l}
-9x_1+7x_2+x_3=-1 \\
6x_1-2x_2+2x_3=6 \\
3x_1-2x_2-x_3=0 \\
\end{array}
\right.
$$

\begin{itemize}
\item Moltiplichiamo la prima riga per $-m_{2,1}=-(-6/9)=2/3$ e la sommiamo alla seconda, ottenendo
 $$
\left\{
\begin{array}{l}
-9x_1+7x_2+x_3=-1 \\
0 x_1+(8/3)x_2+(8/3)x_3=16/3 \\
3x_1-2x_2-x_3=0  \\
\end{array}
\right.
$$
\item Moltiplichiamo la prima riga per $-m_{3,1}=-(-1/3)=1/3$ e la sommiamo alla terza, ottenendo
 $$
\left\{
\begin{array}{l}
-9x_1+7x_2+x_3=-1 \\
0 x_1+(8/3)x_2+(8/3)x_3=16/3 \\
0 x_1+(1/3) x_2-(2/3)x_3=-1 /3\\
\end{array}
\right.
$$
\item Osserviamo che $\max{8/3,1/3}=8/3$ Con le notazioni introdotte, $\phi(2)=2$ e quindi non si richiedono scambi di equazioni. 

Moltiplichiamo la seconda riga per $-m_{3,2}=-(1/3)/(8/3)=-1/8$ e la sommiamo alla terza, ottenendo
 $$
\left\{
\begin{array}{l}
-9x_1+7x_2+x_3=-1 \\
0 x_1+(8/3)x_2+(8/3)x_3=16/3 \\
0 x_1+0 x_2- x_3=-1 \\
\end{array}
\right.
$$
\end{itemize}
Il sistema finale, \mye facile da risolvere. 
\begin{itemize}
\item L'ultima equazione ha una sola variabile e una sola incognita e comporta che $$x_3=1.$$
\item Inseriamo questo risultato nella penultima equazione e otteniamo 
$$ 
(8/3)x_2+(8/3) \cdot 1=16/3
$$
da cui $x_2=1$.
\item Inseriamo questi risultati nella prima equazione e otteniamo 
$$
-9x_1+7 \cdot 1+1=-1
$$
da cui $x_1=1.$
\end{itemize}

\begin{remark}
L'esempio precedente \mye stato scritto in forma di sistema di equazioni, ma pu\myo essere facilmente riscritto direttamente in forma matriciale.
\end{remark}

Passiamo alla fattorizzazione LU della matrice 
 $$
A=\left[
\begin{array}{c c c}
3 & -2 & -1 \\
6 & -2 & 2 \\
-9 & 7 & 1 \\
\end{array}
\right]
$$
utilizzata nell'esempio precedente per risolvere $Ax=b$, con $b=(0,6,-1)^T$.
Per quanto visto, analizzando il sistema finale del processo, $U=A{(n-1)}$ e dunque
$$
U=\left[
\begin{array}{c c c}
-9 & 7 & 1 \\
0 & 8/3 & 8/3 \\
0 & 0 & -1 \\
\end{array}
\right].
$$
mentre ricordando i moltiplicatori $m_{2,1}=-2/3$, $m_{3,1}=1/3$, $m_{3,2}=1/8$, abbiamo
$$
L=\left[
\begin{array}{c c c}
1 & 0 & 0 \\
-2/3 & 1 & 0 \\
-1/3 & 1/8 & 1 \\
\end{array}
\right].
$$
Avendo scambiato la prima con la terza riga, la matrice $P$, detta di permutazione, risulta
$$
P_1=\left[
\begin{array}{c c c}
0 & 0 & 1 \\
0 & 1 & 0 \\
1 & 0 & 0 \\
\end{array}
\right],
$$
mentre al secondo passo, visto che non abbiamo effettuato scambi $P_2=I_n$, con $I_n \in \R^{n \times n}$ matrice identica. Quindi $P=P_2 \cdot P_1=P_1$.
\end{example}

\section{Alcune note sulla eliminazione gaussiana e la fattorizzazione LU}


\begin{remark}

\begin{definition}
La matrice $A \in \C^{n \times n}$ \mye a predominanza diagonale per righe se per ogni 
$i=1,\ldots,n$ risulta 
$${ { |a_{i,i}|  \geq \sum_{j=1,j \neq i}^n |a_{i,j}| }}$$
e per almeno un indice $s$ si abbia
$$|a_{s,s}| >\sum_{j=1,j \neq s}^n |a_{s,j}|.$$

La matrice $A \in \C^{n \times n}$ \mye a predominanza diagonale stretta per righe se per ogni 
$i=1,\ldots,n$ risulta 
$${ { |a_{i,i}|  > \sum_{j=1,j \neq i}^n |a_{i,j}| }}.$$

La matrice $A \in \C^{n \times n}$ \mye a predominanza diagonale (stretta) per colonne  se $A^T$  \mye a predominanza diagonale (stretta) per righe.
\end{definition}

\smsp

Esistono classi di matrici come quelle 
\begin{enumerate}
\item a predominanza diagonale stretta per righe (cf. {\cite{WIKIDD}}), o
\item a predominanza diagonale stretta per colonne (cf. {\cite{WIKIDD}}), o
\item simmetriche e definite positive,
\end{enumerate}
per cui \mye sempre possibile la fattorizzazione $A=LU$ (cf. \cite[p.76]{QUALAR}).
\end{remark}



\begin{remark}
Abbiamo visto che una qualsiasi matrice invertibile possiede una fattorizzazione del tipo $PA=LU$. In realt\mya una matrice arbitraria possiede una fattorizzazione di questo tipo (cf. {\cite[p.145]{BINI}}).
\end{remark}






\begin{remark}
 Il seguente pseudocodice implementa la prima parte dell'algoritmo di eliminazione gaussiana senza pivoting, cui segue la risoluzione di un sistema triangolare superiore.

\begin{verbatim}
for k=1:n-1
  for i=k+1:n
       m(i,k)=a(i,k)/a(k,k);
       for j=k+1:n
         a(i,j)=a(i,j)-m(i,k)*a(k,j);
       end
       b(i)=b(i)-m(i,k)*b(k);
  end 
end
\end{verbatim}

Il costo della eliminazione gaussiana con pivoting \mye essenzialmente lo stesso di quello senza pivoting, visto che a parte il calcolo di un massimo, si tratta esclusivamente di scambiare al pi\myu due righe per ogni indice $k=1,\ldots,n-1$. 

Dal pseudocodice, si vede che il numero di operazioni moltiplicatorie necessarie per effettuare la parte di eliminazione gaussiana, precedente alla sostituzione, \mye
\begin{eqnarray}\label{223}
\sum_{p=1}^{n-1} \sum_{k=p+1}^n (n-p+2) &=&
\sum_{p=1}^{n-1} \sum_{k=p+1}^n 2 + \sum_{p=1}^{n-1} \sum_{k=p+1}^n (n-p) \nonumber \\ &=&
2 \sum_{p=1}^{n-1} (n-p) + \sum_{p=1}^{n-1} (n-p)^2 \nonumber \\
&=& 2 \sum_{k=1}^{n-1} k + \sum_{k=1}^{n-1} k^2 \nonumber \\
\end{eqnarray}

Nell'ultima riga, abbiamo osservato che posto $k=n-p$, 
\begin{itemize}
\item [a)] $\sum_{p=1}^{n-1} (n-p)=\sum_{k=1}^{n-1} k$,
\item [b)] $\sum_{p=1}^{n-1} (n-p)^2=\sum_{k=1}^{n-1} k^2$.
\end{itemize}

Inoltre
\begin{itemize}
\item si ha $2 \sum_{k=1}^{n-1} k=2 \frac{(n-1)n}{2} \approx {n^2}$
\item per quanto riguarda $ \sum_{k=1}^{n-1} k^2$ abbiamo che
\begin{eqnarray}
\frac{(n-1)^3}{3} &=& \int_0^{n-1} x^2 dx \leq \sum_{k=1}^{n-1} k^2 \nonumber \\
&\leq &  \int_1^{n} x^2 dx = \frac{n^3-1}{3}
\end{eqnarray}
da cui, utilizzando la definizione di asintotico e il teorema del confronto,
$$
\sum_{k=1}^{n-1} k^2 \sim  \frac{n^3}{3}.
$$
\end{itemize}

In definitiva, 
$$
2 \sum_{k=1}^{n-1} k + \sum_{k=1}^{n-1} k^2 \approx 
{n^2} + \frac{n^3}{3} \approx  \frac{n^3}{3}.
$$

\begin{itemize}
\item Visto che risolvere ognuno dei sistemi triangolari  costa circa $n^2/2$ mentre la parte precedente necessita di circa $n^3/3$ operazioni moltiplicative, deduciamo che risolvere il sistema $Ax=b$ necessita di circa $n^3/3$ operazioni moltiplicative;

\begin{table}[!ht]
\begin{center}
\begin{tabular}{| c c  | c c |}
\hline
$n$&   cputime  \\
\hline
 	 $    1$ & $2.28081e-04$ \\
 	 $    2$ & $1.14614e-03$\\
 	 $    4$ & $1.14304e-04$\\
 	 $    8$ & $4.85808e-05$\\
 	 $   16$ & $4.79512e-04$\\
 	 $   32$ & $1.15395e-04$\\
 	 $   64$ & $2.23282e-03$\\
 	 $  128$ & $2.85307e-03$\\
 	 $  256$ & $2.95903e-03$\\
 	 $  512$ & $1.71150e-02$\\
 	 $ 1024$ & $4.20015e-02$\\
 	 $ 2048$ & $2.45661e-01$\\
 	 $ 4096$ & $1.55579e+00$\\
 	 $ 8192$ & $9.08754e+00$\\
 	 $16384$ & $6.18414e+01$\\
  \hline
\end{tabular}
\caption{Secondi necessari per risolvere un sistema $Ax=b$, con una generica $A \in \R^{n \times n}$, mediante eliminazione gaussiana su un Mac Book Pro, con processore 2,7 GHz Intel Core i5 e 16 GB di memoria.}
\end{center}
\end{table}
 	 
\item La complessit\mya della fattorizzazione $PA=LU$ \mye essenzialmente lo stesso della eliminazione gaussiana, ovvero $ \frac{n^3}{3}$.
\end{itemize}
\end{remark}












\begin{remark}
In generale l'algoritmo di eliminazione gaussiana per la risoluzione di un sistema lineare $Ax=b$ con $A$ invertibile, \mye implementato come segue
\begin{enumerate}
\item Calcola la fattorizzazione $PA=LU$. A questo punto siccome $P$ \mye invertibile $Ax=b$ se e solo se $LUx=PAx=Pb$. 
\item Calcola $\tilde{b}=Pb$ .
\item Risolvi il sistema triangolare inferiore $Ly=\tilde{b}$ con la sostituzione in avanti.
\item Risolvi il sistema triangolare inferiore $Ux=y$ con la sostituzione all'indietro.
\end{enumerate}
La complessit\mya dell'algoritmo cos\myi implementato \mye di $ \frac{n^3}{3}$ operazioni moltiplicative, ma offre il vantaggio di avere a disposizione la fattorizzazione $PA=LU$.
\end{remark}










\begin{remark}
Se la $A$ matrice \mye simmetrica e definita positiva {\cite{WIKIDP}}, ovvero $A=A^T$ e tutti gli autovalori di $A$ sono positivi, \mye pi\myu conveniente utilizzare la cosidetta decomposizione di Cholesky \cite[p.524]{ATK}, \cite{WIKICHOL}, in cui $A=L L^T$ con $L=(l_{i,j})$ triangolare inferiore con elementi diagonali $l_{k,k} > 0$, $k=1,\ldots,n$. La determinazione di questa fattorizzazione richiede approssimativamente $ \frac{n^3}{6}$ operazioni moltiplicative.

Di conseguenza, per calcolare la soluzione di $Ax=b$, basta risolvere i due sistemi triangolari $Ly=b$, $L^Tx=y$. Quindi la soluzione del sistema lineare si pu\myo ottenere in circa
$$
\frac{n^3}{6}+\frac{n^2}{2}+\frac{n^2}{2} \approx \frac{n^3}{6}
$$
operazioni moltiplicative.

\begin{example}
La matrice
$$ 
A=\left[
\begin{array}{c c c c c}
     1   &  1 &    1  &   1  &   1 \\
     1    & 2  &   2  &   2  &   2\\
     1   &  2    & 3  &   3   &  3\\
     1   &  2   &  3  &   4  &   4\\
     1   &  2   &  3   &  4  &   5\\
\end{array}
\right]
$$
\mye simmetrica e definita positiva, in quanto l'autovalore minimo \mye $\lambda_{\mbox{max}} = 0.271\ldots >0$.

Inoltre $A=L L^T$ con
$$ 
L=\left[
\begin{array}{c c c c c}
     1   &  0 &    0  &   0  &   0 \\
     1    & 1  &   0  &   0  &   0\\
     1   &  1    & 1  &   0   &  0\\
     1   &  1   &  1  &   1  &   0\\
     1   &  1   &  1   &  1  &   1\\
\end{array}
\right].
$$
\end{example}

\end{remark}







\subsection{Alcune applicazioni della fattorizzazione LU}
Di seguito mostriamo alcune applicazione della fattorizzazione LU, ovvero il calcolo del determinante di una matrice e quello della sua inversa, paragonandoli con altri metodi.

\subsubsection{Calcolo del determinante di una matrice}
Il metodo di eliminazione di Gauss pu\myo essere usato per calcolare il determinante di una matrice quadrata.

Il costo computazionale della formula ricorsiva di Laplace per il determinante \mye di circa $2n!$ flops (cf. {\cite{WIKILAP}}.

Se $PA=LU$, visto che $P^{-1}=P^T$ abbiamo che $A=P^{T}LU$ e quindi per il teorema di Binet 
$$
\det(A)=\det(P^{T}LU)=\det(P^{T}) \cdot \det(L) \cdot \det(U).
$$
Osserviamo che 
\begin{itemize}
\item a seconda degli scambi effettuati $\det(P^{T})=(-1)^s$ con $s$ numero di scambi effettuati dal pivoting,
\item visto che il determinante di una matrice triangolare $T=(t_{i,j}) \in \R^{n \times n}$ \mye $$
\det(T)=\prod_{k=1}^n t_{k,k}
$$
abbiamo che da $l_{k,k}=1$, $k=1,\ldots,n$, si ricava $\det(L)=1$, e quindi
$$
\det(A)=(-1)^s \cdot \prod_{k=1}^n u_{k,k}.
$$
con $\det(P^{T})=(-1)^s$ che \mye calcolabile senza operazioni moltiplicative, ma solo contando gli scambi effettuati.
\end{itemize}
In definitiva, visto che $\prod_{k=1}^n u_{k,k}$ necessita solo di $n-1$ operazioni moltiplicative, mentre la fattorizzazione $PA=LU$ ha complessit\mya dell'ordine di $n^3/3$, il calcolo del determinante pu\myo essere effettuato con circa $n^3/3$ operazioni moltiplicative.

Nella tabella paragoniamo il tempo di calcolo necessario a un supercomputer {\cite{WIKISUP}} come Summit che eroga 200 petaflops, ovvero $2 \cdot 10^{17}$ operazioni al secondo, per calcolare il determinante di una matrice generica di ordine $n$. Si tenga conto che si presume che l'et\mya dell'universo sia di circa "$13.82e+9$" anni.


 \begin{table}[!ht]
\begin{center}
\begin{tabular}{| c c c  | c c c |}
\hline
$n$&   $CPU_L$ & $CPU_{LU}$&   $n$&   $CPU_L$ & $CPU_{LU}$  \\
\hline
 	 $  5$ & $6.0e-16$ sec & $2.1e-16$  sec & $ 75$ & $1.2e+92$ anni & $7.0e-13$  sec \\
 	 $ 10$ & $1.8e-11$  sec & $1.7e-15$  sec & $100$ & $4.7e+140$ anni & $1.7e-12$  sec \\
 	 $ 25$ & $9.0e+02$ anni & $2.6e-14$  sec & $125$ & $9.4e+191$ anni & $3.3e-12$  sec \\
 	 $ 50$ & $1.8e+42$ anni & $2.1e-13$  sec & $150$ & $2.9e+245$ anni & $5.6e-12$  sec \\
 	   \hline
\end{tabular}
\caption{In questa tabella paragoniamo il tempo di calcolo $CPU_L$, $CPU_{LU}$ necessario a un supercomputer come Summit che eroga 200 petaflops, ovvero $2 \cdot 10^{17}$ operazioni al secondo, per calcolare il determinante di una matrice generica di ordine $n$, rispettivamente con la regola di Laplace e via fattorizzazione $PA=LU$.}
\end{center}
\end{table}

\subsubsection{Calcolo dell'inversa di una matrice}

Per calcolare l'inversa di una matrice $A \in \R^{n \times n}$ con il {\sl{metodo dei cofattori}} {\cite{WIKICOF}} serve calcolare $n^2+1$ determinanti di matrici estratte da $A$ eliminando la $i$-sima riga e $j$-sima colonna, con $i,j=1,\ldots,n$, pi\myu il determinante di $A$. Se brutalmente si usa la fattorizzazione LU per ognuna di esse, necessitano 
$$
C_{COF} = (n^2+1) \cdot \frac{n^3}{3} \approx \frac{n^5}{3}
$$
operazioni moltiplicative.

Si supponga sia $PA=LU$ e si risolvano i sistemi $$
Ax^{(j)}=e^{(j)}, \,\,\, j=1,\ldots,n
$$
dove $e^{(j)}=(0,\ldots,0,1,0\ldots,0)$ \mye il $j$-simo vettore della base canonica. Allora
$A^{-1}$ \mye la matrice la cui $j$-sima colonna corrisponde a $x^{(j)}$.
Per eseguire i calcoli necessitano
\begin{itemize}
\item una fattorizzazione $LU$ per cui servono circa $n^3/3$ operazioni,
\item risolvere $n$ sistemi lineari con una sostituzione in avanti e una all'indietro. 

Nel caso di $Ly^{(k)}=e^{(k)}$, si vede facilmente che le prime $k$ cifre di $y^{(k)}$ sono nulle (si veda ad esempio (\ref{fs}) ponendo $b_1,\ldots,b_{k-1}=0$) e quindi per risolvere $Ly^{(k)}=e^{(k)}$ servono
$$
1+\ldots+(n-k+1)=(n-k+1)(n-k+2)/2 \approx (n-k+1)^2/2
$$
operazioni moltiplicative.
Quindi la risoluzione di tutti i sistemi $Ly^{(k)}=e^{(k)}$, $k=1,\ldots,n$ ha complessit\mya
$$
\sum_{k=1}^n (n-k+1)^2/2=\sum_{j=1}^n j^2/2=(1/2) \sum_{j=1}^n j^2 \approx \frac{n^3}{6}
$$
Nel caso di $Lz^{(k)}=y^{(k)}$, il costo computazionale \mye invece di 
$n^2/2$ operazioni moltiplicative per ogni $k$, ovvero in totale $n^3/2$ operazioni.
\end{itemize}
Quindi, per calcolare l'inversa di $A$ con fattorizzazione $LU$ servono $$\frac{n^3}{3}+ \frac{n^3}{6}+ \frac{n^3}{2}= \frac{(2+1+3)n^3}{6}=n^3$$  operazioni moltiplicative, e quindi con complessit\mya inferiore a quella richiesta da una implementazione brutale del metodo dei cofattori che \mye risultata essere $n^5/3$.

\begin{example}
Sia
$$
A =\left[ \begin{array}{c c c }
     2  &   5   &  7 \\
     1   &  5   &  2 \\
     1    & 4   &  8 \\
     \end{array} \right].
$$
Essendo 
\begin{itemize}
\item la soluzione $x_1$ di $Ax=[1,0,0]^T$ \mye
$$
x_1 =\left[ \begin{array}{c}
     1.185185185185185e+00 \\
    -2.222222222222222e-01 \\
    -3.703703703703703e-02 \\
     \end{array} \right].
$$
\item la soluzione $x_2$ di $Ax=[0,1,0]^T$ \mye
$$
A =\left[ \begin{array}{c}
    -4.444444444444446e-01 \\
     3.333333333333334e-01\\
    -1.111111111111111e-01\\
     \end{array} \right].
$$
\item la soluzione $x_3$ di $Ax=[0,0,1]^T$ \mye
$$
A =\left[ \begin{array}{c}
    -9.259259259259259e-01 \\
     1.111111111111111e-01\\
     1.851851851851852e-01\\
     \end{array} \right].
$$
\end{itemize}
e quindi
$$
\footnotesize{
A^{-1} =\left[ \begin{array}{c c c }
     1.185185185185185e+00  &  -4.444444444444446e-01   & -9.259259259259259e-01 \\
    -2.222222222222222e-01  &   3.333333333333334e-01  &   1.111111111111111e-01 \\
    -3.703703703703703e-02  &  -1.111111111111111e-01   &  1.851851851851852e-01 \\
     \end{array} \right].
     }
$$
\end{example}

\smsp

\section{Sistemi sovradeterminati, minimi quadrati e fattorizzazione QR}

\begin{theorem}
Sia $A \in \R^{m \times n}$, con $m>n$, una matrice di rango massimo $n$ (ovvero, le colonne di A
sono $n$ vettori linearmente indipendenti di $\R^m$).

Allora la soluzione del sistema sovradeterminato $Ax=b$ 
ai minimi quadrati,
\begin{equation}\label{QR1}
\min_{x \in \R^n} \|b-Ax\|_2= \min_{x \in \R^n} \sqrt{\sum_{i=1}^m (b_i-\sum_{j=1}^n a_{i,j} x_j)^2}
\end{equation}
\mye equivalente a risolvere il sistema con matrice simmetrica, non singolare, definita positiva
$$
A^T A x=A^Tb
$$
detto {\sl{sistema delle equazioni normali}}.
\end{theorem}

\begin{remark}
Essendo il sistema sovradeterminato, se $m > n$ allora ha pi\myu equazioni che incognite e quindi pu\myo non esistere $x \in \R^n$ tale che $Ax=b$ ma comunque (\ref{QR1}) pu\myo avere soluzione (che \mye appunto quella nel senso dei minimi quadrati).
\end{remark}

\begin{proof}
Sia $(u,v)=u^Tv \in \R$ il prodotto scalare in $\R^n$. Ricordiamo che 
\begin{itemize}
\item $(\alpha u,v)=\alpha u^Tv=\alpha (u,v)$ e similmente $( u,\beta v)=\beta u^Tv=\beta (u,v)$;
\item $(u,v)=u^Tv=(v^Tu)=(v,u)$ in quanto il trasposto di un numero \mye il numero stesso;
\item $(u+v,u+v)=(u,u)+(u,v)+(v,u)+(v,v)=(u,u)+(u,v)+(u,v)+(v,v)=(u,u)+2(u,v)+(v,v)$.
\end{itemize}

Si osservi che $x$ minimizza il polinomio quadratico 
$$\phi(x) = \|b-Ax\|^2=(b-Ax,b-Ax)=\sum_{i=1}^m (b_i-\sum_{j=1}^n a_{i,j} x_j)^2$$ 
se e solo se per ogni $v\in \R^n$ si ha 
$$\phi(x) \leq \phi(x+v)$$
ovvero $\phi(x+v)-\phi(x) \geq 0$ e sviluppando i calcoli, 
\begin{eqnarray}
\phi(x+v)&=& \|b-A(x+v)\|^2=\|(b-Ax)-Av\|^2 \nonumber \\
&=& (b-Ax,b-Ax)-2(b-Ax,Av)+(Av,Av) \nonumber \\
&=& \phi(x)-2(b-Ax,Av)+(Av,Av)
\end{eqnarray}
deduciamo che 
\begin{eqnarray}\label{the5_1}
0 &\leq& \phi(x+v)-\phi(x)=-2(b-Ax,Av)+(Av,Av)\nonumber \\ &=& 2(Ax-b,Av)+(Av,Av).
\end{eqnarray}


Notiamo che essendo $v$ arbitrario, (\ref{the5_1}) vale pure per $v^*=-v$ e quindi
\begin{equation}\label{the5_2}
0 \leq \phi(x+v^*)-\phi(x)=2(b-Ax,Av^*)+(Av^*,Av^*)=-2(Ax-b,Av)+(Av,Av).
\end{equation}

Se $v=\epsilon u$ con $u$ versore, ovvero vettore tale che $\|u\|_2=1$ allora da (\ref{the5_1}) abbiamo con facili conti
$$
0 \leq 2\epsilon (Ax-b,Au)+ \epsilon^2 (Au,Au)
$$
ovvero per $\epsilon > 0$ arbitrario
$$
0 \leq 2 (Ax-b,Au)+ \epsilon (Au,Au)
$$
cio\mye 
\begin{equation}\label{the5_3}
0 \leq 2 (Ax-b,Au).
\end{equation}
Similmente da (\ref{the5_1})
$$
0 \leq -2\epsilon (Ax-b,Au)+ \epsilon^2 (Au,Au)
$$
ovvero per $\epsilon > 0$ arbitrario
$$
0 \leq -2 (Ax-b,Au)+ \epsilon (Au,Au)
$$
cio\mye 
\begin{equation}\label{the5_4}
0 \leq -2 (Ax-b,Au).
\end{equation}
Di conseguenza, per $u$ versore arbitrario, dovendo valere contemporaneamente (\ref{the5_3}) e (\ref{the5_4}) , ricaviamo $0 = -2 (Ax-b,Au)$ ovvero
\begin{equation}\label{the5_5}
0=(Ax-b,Au)=x^TA^TAu-b^TAu=(A^TAx)^Tu-(A^Tb)^Tu=(A^TAx-A^Tb,u)
\end{equation}
da cui essendo $u$ versore arbitrario, necessariamente
 $A^TAx-A^Tb=0$ ovvero
$$
A^TAx=A^Tb.
$$

In effetti, se fosse $A^TAx-A^Tb \neq 0$, posto $u=(A^TAx-A^Tb)/\|A^TAx-A^Tb\|$ avremmo
\begin{eqnarray}
0&=&(A^TAx-A^Tb,(A^TAx-A^Tb)/\|A^TAx-A^Tb\|) \nonumber \\
&=&\|A^TAx-A^Tb\|^2/\|A^TAx-A^Tb\|=\|A^TAx-A^Tb\| \nonumber 
\end{eqnarray}
ovvero essendo il vettore nullo l'unico di norma nulla
$$
A^TAx-A^Tb=0
$$
il che sarebbe assurdo perch\mye avevamo supposto $A^TAx-A^Tb \neq 0$.

Per dimostrare che $A^TA$ \mye non singolare, si osservi
che se $A^TAv = 0$ allora 
$$0=(0,v)=(A^TAv, v) = v^TA^TAv= (Av,Av) =\|Av\|^2_2$$ 
da cui siccome $\|x\|=0$ se e solo $x=0$, deduciamo che $Av = 0$. Poich\mye il rango di $A$ \mye $n$, le colonne di $A$ sono linearmente indipendenti e quindi $Av = 0$ implica $v=0$.
Ne comporta che $A^TAv=0$ se e solo $v=0$, ovvero $A^TA$ \mye non singolare (il suo nucleo ha solo il vettore nullo!).

Osserviamo che $A^TA \in \R^{n \times n}$ \mye simmetrica essendo
$$
(A^TA)^T=A^T (A^T)^T=A^TA.
$$

Infine visto che $A^TA$ \mye definita positiva se e solo se 
$$
x^TA^TAx > 0, {\mbox{ per }} x \neq 0
$$
essendo $Ax=0$ se e solo se $x=0$ e 
$$
x^TA^TAx=(Ax)^TAx=(Ax,Ax)=\|Ax\|^2 \geq 0,
$$
deduciamo che $A^TA$ \mye definita positiva.
\end{proof}

\begin{remark}
L'approssimazione polinomiale di grado $k$ ai minimi quadrati pu\myo essere reinterpretata
come soluzione ai minimi quadrati del sistema sovradeterminato $$ Va = y$$ dove
\begin{itemize}
\item $V=(x^{j}_i)_{i=1,\ldots,m,j=0,\ldots,k} \in \R^{m \times (k+1)}$, $m \geq k+1$, \mye una matrice di Vandermonde
\item $y=(y_j)_{j=1,\ldots,m}$ \mye tale che $y_j=f(x_j)$ per $j=1,\ldots,m$.
\end{itemize}

Il teorema precedente asserisce che $a=(a_j)_{j=0,\ldots,k}$ \mye la soluzione del sistema delle equazioni normali
$$
V^T V a = V^T y.
$$
\end{remark}

\begin{definition}
Data una matrice $A \in \R^{m \times n}$, $m \geq n$, avente rango massimo $n$, si dice che $A$ \mye fattorizzabile QR se e solo se esistono
\begin{itemize}
\item la matrice quadrata $Q \in \R^{m \times m}$ unitaria ovvero $QQ^T=Q^TQ=I_m$,
\item la matrice rettangolare $R \in \R^{m \times n}$ triangolare superiore, ovvero $r_{i,j} >0$ se $i > j$.
\end{itemize}
tali che $A=QR$ (cf. {\cite{WIKIQR}}).
\end{definition}

\begin{remark}
Il calcolo della fattorizzazione QR \mye in sostanza equivalente all'ortogonalizzazione di Gram-Schmidt delle colonne di $A$ (cf. {\cite[p.9]{BINI}}, {\cite[p.86]{QUALAR}}).

 \begin{table}[!ht]
\begin{center}
\begin{tabular}{| c c  | c c |}
\hline
$n$&   cputime & $n$&   cputime \\
\hline
  	 $    2$ & $1e-05$ & $  128$ & $3e-03$  \\
 	 $    4$ & $2e-05$ & $  256$ & $5e-03$ \\
 	 $    8$ & $5e-05$ & $  512$ & $2e-02$ \\
 	 $   16$ & $7e-05$ & $ 1024$ & $1e-01$  \\
 	 $   32$ & $1e-04$ & $ 2048$ & $6e-01$ \\
 	 $   64$ & $2e-03$ & $ 4096$ & $4e+00$ \\
  \hline
\end{tabular}
\caption{Secondi necessari per calcolare la fattorizzazione QR  di una generica $A \in \R^{n \times n}$, mediante eliminazione gaussiana su un Mac Book Pro, con processore 2,7 GHz Intel Core i5 e 16 GB di memoria.}
\end{center}
\end{table}


\end{remark}

\begin{example}
Sia 
$$
A=\left[
\begin{array}{c c c}
72 & -144 & -144 \\ -144 & -36 & -360 \\ -144 & -360 & 450
\end{array}
\right]
$$
Il determinante della matrice \mye $\det(A)=-34012224$ e quindi la matrice $A$ ha rango massimo, ovvero $rk(A)=3$.

Si vede che (cf. {\cite[p. 184]{BINI}})
$$
Q=\frac{1}{6} \cdot \left[
\begin{array}{c c c}
-2 & 4 & 4 \\  4 & -2 & 4 \\ 4 & 4 & -2
\end{array}
\right],
\,\,\,\
R=\left[
\begin{array}{c c c}
-216 & -216 & 108 \\  0 & -324 & 324 \\ 0 & 0 & -486
\end{array}
\right].
$$
\end{example}

\begin{example}
Sia 
$$
A=\left[
\begin{array}{c c}
1 & 0 \\ 1 & 1 \\ 0 & 1
\end{array}
\right]
$$
La matrice $A \in \R^{3 \times 2}$ ha rango $2$ perch\mye la sottomatrice $$
A_{1,2}=\left[
\begin{array}{c c}
1 & 0 \\ 1 & 1 
\end{array}
\right]
$$ composta dalle prime due righe \mye non singolare (si ha $\det(A_{1,2})=1$).

Se fattorizziamo $A=QR$ in Matlab otteniamo
$$
\footnotesize{
Q = \left[
  \begin{array}{c c c}
    -7.071067811865472e{\hspace{-0.05cm}}-01  &   4.082482904638630e{\hspace{-0.05cm}}-01  &   5.773502691896258e{\hspace{-0.05cm}}-01\\
    -7.071067811865475e{\hspace{-0.05cm}}-01  &  -4.082482904638630e{\hspace{-0.05cm}}-01  &  -5.773502691896258e{\hspace{-0.05cm}}-01 \\
                         0  &  -8.164965809277261e{\hspace{-0.05cm}}-01  &   5.773502691896256e{\hspace{-0.05cm}}-01 \\
         \end{array}
\right]
}
$$

$$
R = \left[
  \begin{array}{c c }
    -1.414213562373095e+00  &  -7.071067811865475e-01 \\
                         0  &  -1.224744871391589e+00 \\
                         0         &                0 \\
         \end{array}
\right]
$$
In particolare, si ha che 
\begin{itemize}
\item $\|Q'Q-I_m)\|_{\infty} \approx 6.3642e-16$, 
\item $ \|QQ'-I_m)\|_{\infty} \approx 7.3882e-16$, 
\item $\|A-QR)\|_{\infty} \approx 7.3882e-16$.
\end{itemize}
\end{example}

\subsubsection{Applicazione della fattorizzazione QR alla soluzione di sistemi sovradeterminati}

Sia $A \in \R^{m \times n}$, $m \geq n$, una matrice avente rango massimo, e sia $A=QR$.

Se $b \in \R^n$, e supponiamo di dover calcolare la soluzione $x$ ai minimi quadrati di $Ax=b$. Per il teorema precedente $x^*$ \mye pure soluzione delle equazioni normali $$A^TA x =A^Tb.$$

Essendo $A=QR$, e $Q^T Q =I$ ricaviamo 
$$A^TA=(QR)^T QR=R^T Q^T Q R=R^T R$$
e quindi
$$
R^T R x = A^T b.
$$
Ricordando che $R \in \R^{m \times n}$ \mye una matrice triangolare, si risolve
\begin{itemize}
\item $R^T y=A^T b$ ($n$ equazioni e $m$ incognite),
\item $R x=y$ ($m$ equazioni e $n$ incognite, ma solo $n$ sono significative).
\end{itemize} 

\begin{example}
Supponiamo di voler risolvere ai minimi quadrati $Ax=b$ dove 
$$
A = \left[
  \begin{array}{c c }
     2   &  5 \\
     1   &  1 \\
     3   &  2 \\
         \end{array}
\right],
\,\,\,\
b = \left[
  \begin{array}{c }
     1 \\
     1 \\
     1 \\
         \end{array}
\right].
$$

Visto che $A=QR$ con

$$
\footnotesize{
Q = \left[
  \begin{array}{c c c}
    -5.345224838248486e-01  &   8.406260324897477e-01   &  8.737040566610366e-02 \\
    -2.672612419124244e-01  &  -7.005216937414555e-02  &  -9.610744623271416e-01 \\
    -8.017837257372731e-01  &  -5.370666318684501e-01  &   2.621112169983115e-01 \\
         \end{array}
\right]
}
$$

$$
\footnotesize{
R = \left[
  \begin{array}{c c}
    -3.741657386773941e+00   & -4.543441112511214e+00 \\
                         0  &   3.058944729337695e+00 \\
                         0    &                     0 \\
         \end{array}
\right]
}
$$

e 

$$
\footnotesize{
A^T b = \left[
  \begin{array}{c }
6 \\
8\\
         \end{array}
\right]
}
$$

\smsp

\noindent per prima cosa risolviamo il sistema $R^T y=A^Tb$ ovvero

$$
\left\{
\begin{array}{l}
   -3.741657386773941 \cdot y_1+0\cdot y_2+0\cdot y_3 =6  \\
 -4.543441112511214 \cdot y_1 +     3.058944729337695 y_2+0\cdot y_3=8   \\
\end{array}
\right.
$$

\noindent che ha soluzione

$$
{
y = \left[
\begin{array}{c }
    -1.603567451474547e+00 \\
     2.335072312471520e-01 \\
                         0 \\
                         \end{array}
\right].
}
$$

\noindent e quindi $Rx=y$ ovvero
$$
\left\{
\begin{array}{l}
  -3.741657386773941 \cdot x_1 -4.543441112511214 \cdot x_2 =-1.603567451474547  \\
 0 \cdot x_1 +     3.058944729337695 \cdot x_2 =0.2335072312471520   \\
  0 \cdot x_1 +     0 \cdot x_2 =0.   \\
\end{array}
\right.
$$

Alla fine otteniamo

$$
{
x = \left[
\begin{array}{c }
     3.358778625954200e-01 \\
     7.633587786259531e-02 \\
                         \end{array}
\right].
}
$$

\end{example}


\begin{remark} [Facoltativo]
Dal punto di vista del condizionamento \mye meglio risolvere $R^TRx = Q^Tb$ invece di $A^TAx = b$,
perch\mye l'indice di condizionamento di $R$ in norma $2$ \mye molto pi\myu piccolo dell'indice di
condizionamento di $A^TA$. Si pu\myo intuire osservando che se A \mye quadrata e simmetrica, 
\begin{itemize}
\item $k_2(A^TA) = k_2(A^2) = (k_2(A))^2$,
\item $k_2(R) \leq k_2(A)$ visto che $R = Q^TA$ e $\|Q\|_2 = \|Q^T\|_2= \|Q^{-1}\| = 1$ da cui $k_2(Q^T)=1$. In realt\mya si pu\myo mostrare che $k_2(R) = k_2(A)$.
\end{itemize}
\end{remark} 

\section{Metodi iterativi}

\smsp
I metodi che in aritmetica esatta calcolano la soluzione $x^*$ di $Ax=b$, $A \in \R^{n \times n}$ in un numero finito di passi, sono detti metodi {\sl{diretti}}.

Un esempio \mye l'eliminazione gaussiana, che per matrici generali determina la soluzione in circa $n^3/3$ operazioni moltiplicative.

Nei problemi determinati dalle scienze applicate, pu\myo accadere che 
\begin{itemize}
\item la matrice $n$ \mye di grandi dimensioni e quindi ci sono sia problemi di immagazzinamento dati che di complessit\mya computazionale (ovvero $n^3/3$ pu\myo essere proibitivo in termini di costo);
\item la matrice $A$ pu\myo possedere molti zeri, che si verifica spesso non aiutare a migliorare la performance del metodo diretto.
 \end{itemize}
 
 I metodi iterativi tipicamente non \mye detto determinino la soluzione esatta, ma possono approssimarla {\sl{arbitrariamente bene}}. Tipicamente, ogni iterazione \mye basata sull'applicazione di un prodotto matrice-vettore, e quindi ha complessit\mya di circa $n^2$ operazioni moltiplicative.
 
 Di conseguenza, se raggiungono l'approssimazione della soluzione richiesta in un numero di iterazioni $k \ll n/3$, essendo $k n^2 \ll n^3/3$ risultano convenienti rispetto la fattorizzazione LU.


\begin{theorem}
Siano
\begin{itemize}
\item $\| \cdot \|$ una norma matriciale indotta da una norma vettoriale,
\item $B \in \C^{n \times n}$, con $\|B\| < 1$,
\item si supponga che la soluzione $x^*$ di $Ax=b$ verifichi $x^* = Bx^* + c$.
\end{itemize}
Allora qualsiasi sia $x_0 \in \C^n$, la successione $\{x^{(k)}\}$ definita da
$$
x^{(k+1)} = B x^{(k)}+ c, \,\,\, k = 0, 1, 2, \ldots
$$
converge a $x^*$.
\end{theorem}

\begin{proof}
Essendo $\|B\| < 1$, per quanto visto nel teorema \ref{thm3lemma}, abbiamo che la matrice $I-B$ \mye unica e quindi $x^*$ \mye l'unica soluzione di 
$(I-B)x=c$ ovvero di $x=Bx+c$.

Essendo 
$$
x^{(k+1)} = B x^{(k)} + c, \,\,\, k = 0, 1, 2, \ldots
$$
e
$$
x^*= B x^* + c, \,\,\, k = 0, 1, 2, \ldots
$$
necessariamente, sottraendo membro a membro, per la linearit\mya di $B$,
$$
x^{(k+1)}-x^*=(B x^{(k)} + c)-(B x^* + c)=B(x^{(k)}-x^*).
$$
Quindi
\begin{eqnarray}
0 \leq \|x^{(k)}-x^*\|&=&\|B(x^{(k-1)}-x^*)\| \leq \|B\| \|x^{(k-1)}-x^*\| \nonumber \\
&\leq& \ldots \leq  \|B\|^{k+1} \|x^{(0)}-x^*\| \nonumber
\end{eqnarray}
e per il teorema del confronto
$$
\lim_k \|x^{(k)}-x^*\|=0
$$
ovvero $x^{(k)} \rightarrow x^*$.
\end{proof}

\begin{remark}
I metodi della forma
$$
x^{(k+1)} = B x^{(k)}+ c, \,\,\, k = 0, 1, 2, \ldots
$$
per la soluzione del sistema lineare $Ax=b$, sono detti usualmente {\sl{iterativi stazionari}}, perch\mye $B$ e $c$ non dipendono dall'indice di iterazione $k$.
\end{remark}

\begin{theorem}\label{cnorm}
Siano
\begin{itemize}
\item $\| \cdot \|$ una norma matriciale indotta da una norma vettoriale,
\item $B \in \C^{n \times n}$,
\item si supponga che la soluzione $x^*$ di $Ax=b$ verifichi $x^* = Bx^* + c$.
\end{itemize}
Allora qualsiasi sia $x^{(0)} \in \C^n$, la successione $\{x^{(k)}\}$ definita da
$$
x^{(k+1)} = B x^{(k)} + c, \,\,\, k = 0, 1, 2, \ldots
$$
converge a $x^*$ se e solo se $\rho(B) < 1$.
\end{theorem}

\begin{proof}
Dimostriamo l'asserto nel caso $B=(b_{i,j}) \in \C^{n \times n}$ sia diagonalizzabile, ovvero in cui i suoi autovettori $\{v_k\}$ formino una base di $\C^n$.

\smsp

($\Leftarrow$) Poniamo di seguito $e^{(k)}:=x^{(k)} - x^*$ e supponiamo $\rho(B) < 1$.

Visto che $x^{(0)}, x^* \in \C^n$, abbiamo $$e^{(0)}:=x^{(0)}-x^=\sum_{s=1}^n c_s v_s$$ con $Bv_s=\lambda_s v_s$, $s=1,\ldots,n$.

Ma allora, per la linearit\mya di $B$, 
\begin{eqnarray}
e^{(1)}&:=&x^{(1)} - x^*= (B x^{(0)} + c) -(B x^{*} + c) \nonumber \\
&=& B (x^{(0)}-x^{*}) =B e^{(0)} = B \sum_{s=1}^n c_s v_s \nonumber \\
&=& \sum_{s=1}^n c_s B v_s = \sum_{s=1}^n c_s \lambda_s v_s,
\end{eqnarray}
come pure
\begin{eqnarray}
e^{(2)}&:=&x^{(2)} - x^*= (B x^{(1)} + c) -(B x^{*} + c) \nonumber \\
&=& B e^{(1)} = B \sum_{s=1}^n c_s \lambda_s v_s \nonumber \\
&=& \sum_{s=1}^n c_s \lambda_s  B v_s = \sum_{s=1}^n c_s \lambda^2_s v_s,
\end{eqnarray}
e pi\myu in generale, per $k=1,2,3,\ldots$
\begin{eqnarray}
e^{(k)}&:=& \sum_{s=1}^n c_s \lambda^k_s v_s,
\end{eqnarray}
Ma allora 
\begin{eqnarray}
0 \leq \|e^{(k)}\|&:=& \| \sum_{s=1}^n c_s \lambda^k_s v_s \| \leq  \sum_{s=1}^n |c_s| |\lambda^k_s| \|v_s\| \nonumber \\
& \leq & \max_{s=1,\ldots,n}  |\lambda^k_s| \sum_{s=1}^n |c_s| \|v_s\|  \nonumber \\
&=& (\rho(B))^k \sum_{s=1}^n |c_s| \|v_s\|
\end{eqnarray}
ed essendo $\rho(B) < 1$ concludiamo per il teorema del confronto che
$$
\lim_k \|e^{(k)}\|=0,
$$
e quindi $x^{(k)} \rightarrow x^*$.

\smsp

($\Rightarrow$) Mostriamo il fatto che se il metodo converge per qualsiasi $x^{(0)} \in \C^n$ allora $\rho(B) < 1$. Supponiamo per assurdo sia $\rho(B) \geq 1$. Quindi esiste un autovalore $\lambda$ con autovettore $v$ tale che $|\lambda| \geq 1$.

Sia $x^{(0)} \in \C^n$ scelto cosicch\mye sia $e^{(0)}=x^*-x^{(0)}=v$, ovvero $x^{(0)}=x^{*}-v$. Allora
\begin{eqnarray}
e^{(1)}&:=&x^{(1)} - x^*= (B x^{(0)} + c) -(B x^{*} + c) \nonumber \\
&=& B (x^{(0)}-x^{*}) =B e^{(0)} = B v=\lambda v,
\end{eqnarray}
come pure
\begin{eqnarray}
e^{(2)}&:=&x^{(2)} - x^*= (B x^{(1)} + c) -(B x^{*} + c) \nonumber \\
&=& B (x^{(1)}-x^{*}) =B e^{(1)} = B (\lambda v)=\lambda Bv =\lambda^2 v 
\end{eqnarray}
e pi\myu in generale
$$
e^{(k)}=\lambda^k v.
$$
Ma allora
$$
\|e^{(k)}\|=|\lambda^k| \|v\|=|\lambda|^k \|v\|.
$$
Se 
\begin{itemize}
\item $|\lambda| > 1$ ne consegue che $\|e^{(k)}\| \rightarrow +\infty$,
\item se $|\lambda|=1$ abbiamo che
$
\|e^{(k)}\|= \|v\| \neq 0
$
e quindi non \mye infinitesima, perch\mye costantemente uguale a una costante strettamente positiva.
\end{itemize}
\end{proof}


Sia 
\begin{itemize}
\item $A \in \R^{n \times n}$ una matrice quadrata, 
\item $A = P - N$ un cosidetto {\sl{splitting della matrice}} $A$, con $\det(P) \neq 0$
\end{itemize}
Allora il sistema $(P-N)x=Ax = b$ e quindi $Px=Nx+b$ da cui, per l'invertibilit\mya di $P$ ricaviamo
\begin{equation}\label{staz}
x=P^{-1}Nx+P^{-1}b.
\end{equation}
Quindi, posto $B:=P^{-1}N$, $c:=P^{-1}b$, la soluzione $x^*$ di $Ax=b$ risolve pure $x=Bx+c$.

Questa formulazione suggerisce di considerare metodi le cui iterazioni siano fornite dalle iterazioni successive 
\begin{equation}\label{isucc}
x^{(k+1)})=P^{-1}Nx^{(k)})+P^{-1}b.
\end{equation}

Sia $A=D-E-F$ con
\begin{enumerate}
\item $D$ matrice diagonale,
\item $E$ triangolare inferiore,
\item $F$ triangolare superiore.
\end{enumerate}
Allora
\begin{itemize}
\item il metodo di Jacobi corrisponde a scegliere, se $D$ \mye invertibile, $P=D$, $N=E+F$,
\item il metodo di Gauss-Seidel corrisponde a scegliere, se $D$ \mye invertibile, $P=D-E$, $N=F$.
\end{itemize}

\begin{remarkh}
Il metodo di Jacobi fu pubblicato dall'autore nel 1845 in un giornale di astronomia, mentre quello che conosciamo come metodo di Gauss-Seidel fu il risultato di una collaborazione di Seidel con Jacobi nel 1874. Gauss aveva gi\mya descritto un metodo analogo nel 1845 \cite[p.466]{EPS}. 
\end{remarkh}

\begin{remark}
Esiste la seguente generalizzazione del metodo di Gauss-Seidel:
\begin{definition}[Successive over relaxation, (SOR), 1884-1950 (Frankel-Young)]
Sia $A=D-E-F \in \C^{n \times n}$ dove
\begin{itemize}
\item $D$ \mye una matrice diagonale, {{non singolare}}
\item $E$ \mye triangolare inferiore con elementi diagonali nulli,
\item $F$ \mye triangolare superiore con elementi diagonali nulli.
\end{itemize}
Sia fissato $x^{(0)} \in \C^n$ e $\omega \in \C \backslash 0$.
Il {\sl{metodo SOR}} genera $\{x^{(k)}\}_{k}$ con
$$
x^{(k+1)}=\phi(x^{(k)})=M^{-1}Nx^{(k)}+M^{-1}b
$$
dove 
$${{M=    { \frac {D} {\omega} }   -E  , \,\,\, N=\left({ \frac 
{1} 
{\omega}  }-1 \right)D+F}}$$
\end{definition}
Si nota che per il metodo SOR con $w=1$ corrisponde a quello di Gauss-Seidel.
\end{remark}

\begin{remark}
Per matrici simmetriche e definite positive esistono metodi pi\myu adatti che non verranno trattati, come quelli della famiglia del gradiente. Per dettagli su questi importanti metodi e su altri che rappresentano lo stato dell'arte si veda ad esempio \cite{KEL}.
\end{remark}

\begin{example}
Mostriamo un esempio di come una matrice $A$ pu\myo essere riscritta in termini di $D$, $E$, $F$ come descritti sopra.

Se
$$
\footnotesize{
A = \left[
  \begin{array}{c c c c}
     1   &  2   &  3   &  4 \\ 
     5   &  6  &   7  &   2 \\
     8    & 9  &   1  &   2 \\
     3    & 4  &   5 &    1 \\
         \end{array}
\right],}
$$
allora
$$
\footnotesize{
D = \left[
  \begin{array}{c c c c}
     1   &  0   &  0   &  0 \\ 
     0   &  6  &   0  &   0 \\
     0    & 0  &   1  &   0 \\
     0    & 0  &   0 &    1 \\
         \end{array}
\right],
\,\,
E = \left[
  \begin{array}{c c c c}
     0   &  0   &  0   &  0 \\ 
     -5   &  0  &   0  &   0 \\
     -8    & -9  &   0  &   0 \\
     -3    & -4  &   -5 &    0 \\
         \end{array}
\right],
\,\,\
F = \left[
  \begin{array}{c c c c}
     0   &  -2   &  -3   &  -4 \\ 
     0   &  0  &   -7  &   -2 \\
     0    & 0  &   0  &   -2 \\
     0    & 0  &   0 &    0 \\
         \end{array}
\right]
}
$$
\end{example}

\begin{remark}
I metodi di Jacobi e Gauss-Seidel sono riscrivibili componente per componente rispettivamente come
\begin{equation}\label{Jac}
	x_{i}^{(k+1)}=(b_{i} - \sum_{j=1}^{i-1} a_{ij} x_{j}^{(k)}- 
	\sum_{j=i+1}^{n} a_{ij} x_{j}^{(k)} )/a_{ii}, \,\, i=1,\ldots,n,
\end{equation}	
\begin{equation}\label{GS}
	x_{i}^{(k+1)}=\left( b_{i} - \sum_{j=1}^{i-1} a_{ij} x_{j}^{(k+1)}- 
	\sum_{j=i+1}^{n} a_{ij} x_{j}^{(k)} \right)/a_{ii}, \,\, i=1,\ldots,n.
\end{equation}
\end{remark}

\begin{theorem}
Se la matrice $A$ \mye a predominanza diagonale (stretta) per righe, allora i metodi di Jacobi e Gauss-Seidel convergono.
\end{theorem}

\begin{proof}
Dimostriamo esclusivamente la parte relativa al metodo di Jacobi.

Essendo la matrice a predominanza diagonale stretta per righe, si ha 
$$
|a_{i,j}| > \sum_{i=1}^n |a_{i,j}|, \,\,\, j \neq i,
$$
ovvero
$$
\sum_{i=1}^n \frac{|a_{i,j}|}{|a_{i,i}|} < 1, \,\,\, j \neq i.
$$
Poich\mye $P=D$, $N=E+F$, allora le iterazioni sono del tipo $x^{(k+1)}=B_J x^{(k)}+c$ con
$$
B_J=P^{-1}N=D^{-1}(E+F).
$$
In particolare, visto che $(D^{-1}(E+F)|)_{i,j}=-a_{i,j}/a_{i,i}$ se $i \neq j$, e $0$ altrimenti,
\begin{eqnarray}
\|B_J\|_{\infty} &=&  \max_{1 \leq i \leq n} \sum_{j=1, j\neq i}^n |-a_{i,j}/a_{i,i}|  \nonumber \\
&=&   \max_{1 \leq i \leq n} \sum_{j=1, j\neq i}^n |a_{i,j}|/|a_{i,i}| < 1
\end{eqnarray}
e quindi per (\ref{cnorm}), il metodo converge per ogni $x^{(0)} \in \C^n$.
\end{proof}
Inoltre,
\begin{theorem}
Se la matrice $A$ \mye simmetrica e definita positiva, allora il metodo di Gauss-Seidel converge.
\end{theorem}

Per una dimostrazione si veda \cite[p.132]{QUALAR}.

\subsubsection{Alcuni esempi}

\begin{example}
Consideriamo il sistema lineare $Ax=b$ dove
$$
\footnotesize{
A = \left[
  \begin{array}{c c c }
     11   & -5   &  5 \\ 
     5   &  12  &   6 \\
     6    & -4  &   11 \\
         \end{array}
\right],
\,\,\,
b = \left[
  \begin{array}{c}
     1 \\ 
    23 \\
     13 \\
         \end{array}
\right],
}
$$
che ha soluzione $x^*=[1,1,1]^T$ (cf. {\cite[p.252]{BINI}}).

 \begin{table}[!ht]
\begin{center}
\begin{tabular}{| c c c |}
\hline
$n$&   $E_J$ & $E_{GS}$  \\
\hline
  	  	 $ 10$ & $1.1e-01$ & $1.8e-01$ \\
 	 $ 30$ & $1.1e-03$ & $5.1e-03$ \\
 	 $ 50$ & $1.1e-05$ & $1.4e-04$ \\
 	 $ 70$ & $1.0e-07$ & $4.0e-06$ \\
 	 $ 90$ & $9.6e-10$ & $1.1e-07$ \\
 	 $110$ & $9.4e-12$ & $3.2e-09$ \\
 	 $130$ & $8.7e-14$ & $8.8e-11$ \\
 	 $150$ & $8.9e-16$ & $2.5e-12$ \\
 	 $170$ & $-$ & $6.9e-14$ \\
 	 $190$ & $-$ & $1.9e-15$ \\
  \hline
\end{tabular}
\caption{Errori assoluti $\|x^{(k)}-x^*\|_{\infty}$ compiuti dalle $n$-sime iterazioni di Jacobi e Gauss-Seidel, nel risolvere un certo problema $Ax=b$ con $A$ a predominanza diagonale stretta, con vettore iniziale $x_0=0 \in \R^3$.}
\end{center}
\end{table}


La matrice $A$ \mye a predominanza diagonale in senso stretto e quindi tanto il metodo di Jacobi, quanto quello di Gauss-Seidel, risultano convergenti. In effetti $\rho(B_J) \approx 0.79$, $\rho(B_{GS}) \approx 0.83$.

Gli errori assoluti $\|x^{(k)}-x^*\|_{\infty}$ compiuti rispettivamente dal metodo di Jacobi e Gauss-Seidel sono esposti nella tabella che segue, da cui si evince la convergenza dei metodi.

\end{example}

\begin{example}
Consideriamo il sistema lineare $Ax=b$ dove
$$
\footnotesize{
A = \left[
  \begin{array}{c c c c c c c c c}
 	 $4$ & $-1$ & $0$ & $-1$ & $0$ & $0$ & $0$ & $0$ & $0$ \\ 
 	 $-1$ & $4$ & $-1$ & $0$ & $-1$ & $0$ & $0$ & $0$ & $0$ \\ 
 	 $0$ & $-1$ & $4$ & $0$ & $0$ & $-1$ & $0$ & $0$ & $0$ \\ 
 	 $-1$ & $0$ & $0$ & $4$ & $-1$ & $0$ & $-1$ & $0$ & $0$ \\ 
 	 $0$ & $-1$ & $0$ & $-1$ & $4$ & $-1$ & $0$ & $-1$ & $0$ \\ 
 	 $0$ & $0$ & $-1$ & $0$ & $-1$ & $4$ & $0$ & $0$ & $-1$ \\ 
 	 $0$ & $0$ & $0$ & $-1$ & $0$ & $0$ & $4$ & $-1$ & $0$ \\ 
 	 $0$ & $0$ & $0$ & $0$ & $-1$ & $0$ & $-1$ & $4$ & $-1$ \\ 
 	 $0$ & $0$ & $0$ & $0$ & $0$ & $-1$ & $0$ & $-1$ & $4$ \\
         \end{array}
\right],
\,\,\,
b = \left[
  \begin{array}{c}
          2 \\
     1\\
     2\\
     1\\
     0\\
     1\\
     2\\
     1\\
     2\\
         \end{array}
\right],
}
$$
che ha soluzione $x^*=[1,\ldots,1]^T \in \R^9$. La matrice $A$ \mye un esempio di matrice di Poisson ed \mye definita positiva, in quanto il suo autovalore pi\myu piccolo vale $$\lambda_{\mbox{min}} \approx 1.171572875253810e+00 >0.$$


Gli errori assoluti $\|x^{(k)}-x^*\|_{\infty}$ compiuti rispettivamente dal metodo di Jacobi e Gauss-Seidel partendo da $x_0=[0,\ldots,0]^T \in \R^9$. Entrambi i metodi convergono, e per quanto concerne il metodo di Gauss-Seidel, ci\myo era garantito dal teorema di convergenza.  

 \begin{table}[!ht]
\begin{center}
\begin{tabular}{| c c c |}
\hline
$n$&   $E_J$ & $E_{GS}$  \\
\hline
 	 $ 10$ & $4.7e-02$ & $1.8e-03$ \\
 	 $ 30$ & $4.6e-05$ & $1.7e-09$ \\
 	 $ 50$ & $4.5e-08$ & $1.8e-15$ \\
 	 $ 70$ & $4.4e-11$ & $-$ \\
 	 $ 90$ & $4.3e-14$ & $-$ \\
 	   \hline
\end{tabular}
\caption{Errori assoluti $\|x^{(k)}-x^*\|_{\infty}$ compiuti dalle $n$-sime iterazioni di Jacobi e Gauss-Seidel, nel risolvere un certo problema $Ax=b$ con $A$ matrice di Poisson, con vettore iniziale $x_0=0 \in \R^9$.}
\end{center}
\end{table}
 	 \end{example}
 	 
 	 
 	 
 	 
 	\begin{example}
 	Consideriamo il sistema $Ax=b$ con
 	 $$
A=\left[
\begin{array}{c c c}
3 & -2 & -1 \\
6 & -2 & 2 \\
-9 & 7 & 1 \\
\end{array}
\right],
\,\,
b=\left[
\begin{array}{c}
0 \\
6 \\
-1 \\
\end{array}
\right],
$$
avente soluzione $x^*=(1,1,1)^T$.

Utilizziamo i metodi di Jacobi e Gauss-Seidel con dato iniziale $$x^{(0)}=(1-10^{-14},1-10^{-14},1-10^{-14})^T$$
e quindi molto prossimo al vettore soluzione $x^*$ in quanto $\|x^*-x^{(0)}\|_{\infty}=10^{-14}$.

Entrambi i metodi citati, pur partendo da un tal $x^{(0)}$, non convergono, come si pu\myo vedere dalla relativa tabella.

 \begin{table}[!ht]
\begin{center}
\begin{tabular}{| c c c |}
\hline
$n$&   $E_J$ & $E_{GS}$  \\
\hline
 	 $ 10$ & $2.6e-12$ & $1.4e-05$ \\
 	 $ 30$ & $6.7e-09$ & $6.3e+14$ \\
 	 $ 50$ & $1.8e-05$ & $2.9e+34$ \\
 	 $ 70$ & $4.5e-02$ & $1.4e+54$ \\
 	 $ 90$ & $1.2e+02$ & $6.3e+73$ \\
 	 $110$ & $2.9e+05$ & $2.9e+93$ \\
 	 $130$ & $7.3e+08$ & $1.4e+113$ \\
 	 $150$ & $1.8e+12$ & $6.3e+132$ \\
 	 $170$ & $4.4e+15$ & $2.9e+152$ \\
 	   \hline
\end{tabular}
\caption{Errori assoluti $\|x^{(k)}-x^*\|_{\infty}$ compiuti dalle $n$-sime iterazioni di Jacobi e Gauss-Seidel, nel risolvere un certo problema $Ax=b$ con $A$ matrice utilizzata come esempio per l'eliminazione gaussiana, con vettore iniziale $x^{(0)}=(1-10^{-14},1-10^{-14},1-10^{-14})^T  \in \R^3$.}
\end{center}
\end{table}

 	\end{example}


\subsubsection{Metodi di Richardson stazionari}

Da (\ref{staz}), essendo $A=P-N$ ovvero $N=P-A$ con $P$ invertibile, deduciamo 
$$
x=P^{-1}Nx+P^{-1}b=P^{-1}(P-A)x+P^{-1}b= (I - P^{-1 }A)x+P^{-1}b
$$
e quindi il metodo delle approssimazioni successive (\ref{isucc}) si pu\myo riscrivere come
\begin{eqnarray}\label{itrichst}
x^{(k+1 )}&=& (I - P^{-1 }A)x^{(k)}+P^{-1}b=x^{(k)}+P^{-1 }(b-Ax^{(k)}) \nonumber \\
&=&x^{(k)}+P^{-1 }r(x^{(k)}),
\end{eqnarray}
dove $r(x^{(k )} ) = b - Ax^{(k )} $ \mye il vettore detto {\sl{residuo}} al passo $k$-esimo.

Un metodo le cui iterazioni siano del tipo (\ref{itrichst}), si dice di {\sl{Richardson stazionario}}.

\smsp

Osserviamo che essendo $x^*=Bx^*+c$, con 
\begin{itemize}
\item $B=P^{-1}N$, dove $A=P-N$ ovvero $N=P-A$,
\item $c=P^{-1}b$,
\end{itemize}
allora 
\begin{eqnarray}
P^{-1}Ax^*&=&(I-I+P^{-1}A))x^* =(I-P^{-1}(P-A))x^* \nonumber \\
&=&(I-P^{-1}N)x^*=(I-B)x^*=c=P^{-1}b. \nonumber
\end{eqnarray}

Quindi
\begin{equation}\label{prpr}
P^{-1}Ax^* =P^{-1}b,
\end{equation}
ovvero invece di $Ax=b$ si studia il problema equivalente (\ref{prpr}) in cui ci si {\sl{augura}} che $k(P^{-1}A) \ll k(A)$. Una tal matrice si chiama di {\sl{precondizionamento}}. Di conseguenza la matrice $P$ dello splitting $A=P-N$ pu\myo essere vista come matrice di precondizionamento.

Si mostra che l'azione di $P^{-1 }$\mye efficace quando $P^{-1 } \approx A^{-1 }$, nel senso che gli autovalori di $P^{-1 } A$ si accumulano intorno ad $1$ e quindi la convergenza diventa pi\myu rapida. 

\begin{remark}
Si noti che per quanto concerne i metodi di Jacobi e di Gauss-Seidel, dato un vettore $v$, il calcolo di $z = P^{−1}v$, ovvero la soluzione del sistema $Pz = v$, ha un basso costo computazionale, essendo $P$ diagonale o triangolare superiore. 
\end{remark}

\begin{remark}
Vari metodi introducono un parametro $\alpha$ di rilassamento, 
$$x^{(k+1 )}=x^{(k )}+\alpha P^{-1} r(x^{(k)})$$
che aumenti l'efficacia del precondizionatore, cercando di diminuire o addirittura minimizzare il raggio spettrale di $B_{\alpha}=I-\alpha P^{-1}A$, che ha l'effetto di aumentare la velocit\mya di convergenza.
\end{remark}

\subsubsection{Test di arresto}

Sia $\epsilon >0$ una tolleranza fissata dall'utente.
Il metodo iterativo
$$
x^{(k+1)}=Bx^{(k)}+c
$$
viene arrestato utilizzando 
\begin{itemize}
\item il {\sl{criterio dello step}} (cf. {\cite[p.159]{QUALAR}}), ovvero si interrompe il processo quando
$$
\|x^{(k+1)} - x^{(k)}\| \leq \epsilon,
$$
\item {\sl{criterio del residuo}} , ovvero si interrompe il processo quando
$$
\|b-Ax^{(k)}\| \leq \epsilon,
$$
\item {\sl{criterio del residuo relativo}} , ovvero si interrompe il processo quando
$$
{ \frac {\|b-Ax^{(k)}\|} {\|b\|} } \leq \epsilon,
$$
\end{itemize}
Si dimostra che vale la seguente stima dell'errore relativo 
$$
\frac {\|x^*-x^{(k)}\|}  {x^*} \leq  {k(A)}  { \frac {\|b-Ax^{(k)}\|} {\|b\|} } .
$$
Quindi, se $k(A) \gg 1$, nonostante ${ \frac {\|b-Ax^{(k)}\|} {\|b\|} } < \epsilon$, ovvero possa essere piccolo, possiamo avere che l'errore relativo $\frac {\|x^*-x^{(k)}\|}  {x^*}$ sia molto grande (cf. {\cite[p.161]{QUALAR}}), mentre la stima \mye buona qualora $k(A) \approx 1$.



\smsp

{\bf{Bibliografia}}

\begin{thebibliography}{99}
\bibitem{ATK} K. Atkinson, {\sl{Introduction to Numerical Analysis}}, Wiley, 1989.
\bibitem{ATKNEW} K. Atkinson, W. Han, {\sl{Elementary Numerical Analysis}}, Wiley, 2003.
\bibitem{BINI} D. Bini, M. Capovani, O. Menchi, {\sl{Metodi numerici per l'algebra lineare}}, Zanichelli, 1988.
%\bibitem{COM} V. Comincioli, {\sl{Analisi Numerica, metodi modelli applicazioni}}, Mc Graw-Hill, 1990.
% \bibitem{CON} S.D. Conte e C. de Boor, {\sl{Elementary Numerical Analysis, 3rd Edition}}, Mc Graw-Hill, 1980.
\bibitem{EPS} J.F. Epperson, {\sl{Introduzione all'analisi numerica. Teoria, metodi, algoritmi.}}, Mc Graw-Hill, 2009.
\bibitem{KEL} C.T. Kelley, {\em Iterative Methods for Linear and Nonlinear Equations\/}, SIAM Frontiers in Applied Mathematics. SIAM, Philadelphia, (1995).
\bibitem{QUASM} A. Quarteroni e F. Saleri, {\sl{Elementi di calcolo numerico}}, Progetto Leonardo, 1999.
% \bibitem{QUA} A. Quarteroni e F. Saleri, {\sl{Introduzione al calcolo scientifico}}, Springer Verlag, 2006.
\bibitem{QUALAR} A. Quarteroni, R. Sacco e F. Saleri, {\sl{Matematica Numerica}}, Springer Verlag, 1998.
%\bibitem{SULI} A. Suli e D. Mayers, {\sl{An Introduction to Numerical Analysis}}, Cambridge University Press, 2003.
\bibitem{TREF} L.N. Trefethen, {\it{Numerical Analysis}}, Princeton Companion to Mathematics, 2008.
%\bibitem{WIKINC} Wikipedia, Formule di Newton-Cotes,\\   {\href{https://it.wikipedia.org/wiki/Formule\_di\_Newton-Cotes}{https://it.wikipedia.org/wiki/Formule\_di\_Newton-Cotes}}.

\bibitem{WIKICHOL} Wikipedia, Decomposizione di Cholesky,\\   {\href{https://it.wikipedia.org/wiki/Decomposizione\_di\_Cholesky}{https://it.wikipedia.org/wiki/Decomposizione\_di\_Cholesky}}
\bibitem{WIKILU} Wikipedia, Decomposizione LU,\\   {\href{https://it.wikipedia.org/wiki/Decomposizione\_LU}{https://it.wikipedia.org/wiki/Decomposizione\_LU}}
\bibitem{WIKIQR} Wikipedia, Decomposizione QR,\\   {\href{https://it.wikipedia.org/wiki/Decomposizione\_QR}{https://it.wikipedia.org/wiki/Decomposizione\_QR}}
\bibitem{WIKIDD} Wikipedia, Matrice a diagonale dominante,\\   {\href{https://it.wikipedia.org/wiki/Matrice\_a\_diagonale\_dominante}{https://it.wikipedia.org/wiki/Matrice\_a\_diagonale\_dominante}}
\bibitem{WIKIDP} Wikipedia, Matrice definita positiva,\\   {\href{https://it.wikipedia.org/wiki/Matrice\_definita\_positiva}{https://it.wikipedia.org/wiki/Matrice\_definita\_positiva}}
\bibitem{WIKICOF} Wikipedia, Metodo della matrice dei cofattori,\\   {\href{Matrice\_invertibile\#Metodo\_della\_matrice\_dei\_cofattori}{Matrice\_invertibile\#Metodo\_della\_matrice\_dei\_cofattori}}.
\bibitem{WIKISOR} Wikipedia, Successive Over Relaxation\\   {\href{http://it.wikipedia.org/wiki/Successive\_Over\_Relaxation}{http://it.wikipedia.org/wiki/Successive\_Over\_Relaxation}}
\bibitem{WIKISUP} Wikipedia, Supercomputer,\\   {\href{https://it.wikipedia.org/wiki/Supercomputer}{https://it.wikipedia.org/wiki/Supercomputer}}.
\bibitem{WIKIBIN} Wikipedia, Teorema di Binet,\\   {\href{https://it.wikipedia.org/wiki/Teorema\_di\_Binet}{https://it.wikipedia.org/wiki/Teorema\_di\_Binet}}.
\bibitem{WIKILAP} Wikipedia, Teorema di Laplace,\\   {\href{https://it.wikipedia.org/wiki/Teorema\_di\_Laplace}{https://it.wikipedia.org/wiki/Teorema\_di\_Laplace}}.


\end{thebibliography}


\end{document}

 \begin{figure}[ht!]\label{fig1aaaahjtfty}
 \centering
   {\includegraphics[scale=0.2725,clip]{....eps}}
 \caption{...}
 \end{figure}


 \begin{table}[!ht]
\begin{center}
\begin{tabular}{| c c  | c c |}
\hline
$n$&   $E_n$ & $n$&   $E_n$  \\
\hline
  	 $       2$ & $9.62e-01$ & $     512$ & $1.58e-07$ \\
  \hline
\end{tabular}
\caption{....}
\end{center}
\end{table}






\begin{lstlisting}%[frame=single]

function [x,L,U,P]=lupivot(A,b)
%------------------------------------
% INPUT:
% A matrice.
% b termine noto.
% OUTPUT:
% x soluzione del sistema Ax=b.
% A matrice elaborata dalla funzione.
% r vettore degli scambi.
%-------------------------------------

[n,n]=size(A); x=zeros(n,1); y=zeros(n,1); c=zeros(1,n); 
r=1:n;

for p=1:n-1

    % Calcola il pivot relativo alla "p"-sima colonna.
    [max1,j]=max(abs(A(p:n,p)));

    % Scambia righe "p" e "j".
    c=A(p,:);A(p,:)=A(j+p-1,:); A(j+p-1,:)=c;
    d=r(p); r(p)=r(j+p-1); r(j+p-1)=d;

    if A(p,p) == 0, fprintf('\n \t [WARNING]: A NON INVERTIBILE'); break; end

    % Eliminazione. 
    for k=p+1:n
        mult=A(k,p)/A(p,p); A(k,p)=mult; 
        A(k,p+1:n)=A(k,p+1:n)-mult*A(p,p+1:n);
    end
    
end

% Calcola "y".
y(1)=b(r(1));
for k=2:n, y(k)=b(r(k))-A(k,1:k-1)*y(1:k-1); end

% Calcola "x".
x(n)=y(n)/A(n,n);
for k=n-1:-1:1, x(k)=(y(k)-A(k,k+1:n)*x(k+1:n))/A(k,k); end

% Calcoliamo P, L, U.
U=(tril(A'))'; L=A-U+eye(size(A)); 

P=zeros(size(A));
for k=1:n, P(k,r(k))=1; end
\end{lstlisting}
%\end{verbatim}
