% 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[9pt]{beamer}
\input{hal_prefs}
\setlength{\tabcolsep}{2.5pt}
\setbeamertemplate{navigation symbols}{} %no nav symbols
\title[Algebra lineare numerica \hspace{2em}\insertframenumber/
\inserttotalframenumber]{Algebra lineare numerica}
\author{Alvise Sommariva}
\institute{Universit\`a degli Studi di Padova}
\date{\today}
\begin{document}
\frame{\titlepage}



\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

\smnsp

\small{
\begin{definizione}[Norma vettoriale (spazio vettoriale sopra $\mathbb{R}$)]
Dato uno spazio vettoriale $V$ (sopra il campo $\mathbb{R}$), la funzione $\|\cdot\|: V \rightarrow \mathbb{R}$ \mye una {\reds{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|_{\mathbb{R}} \| x\|$ per ogni $x \in V$ e scalare $\alpha \in \mathbb{R}$,
\item vale la cosidetta disuguaglianza triangolare $$\|x+y\| \leq \|x\|+\|y\|$$ per ogni $x,y \in V$.
\end{enumerate}
\end{definizione}

\smnsp

\begin{definizione}[Norma vettoriale (spazio vettoriale sopra $\mathbb{C}$)]
Dato uno spazio vettoriale $V$ (sopra il campo $\mathbb{C}$), la funzione $\|\cdot\|: V \rightarrow \mathbb{R}$ \mye una {\reds{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|_{\mathbb{C}} \| x\|$ per ogni $x \in V$ e scalare $\alpha \in \mathbb{C}$,
\item vale la cosidetta disuguaglianza triangolare $$\|x+y\| \leq \|x\|+\|y\|$$ per ogni $x,y \in V$.
\end{enumerate}
\end{definizione}

\smnsp

\begin{nota}[Norma in $\mathbb C$]
Si ricorda che se $z=a+i \cdot b$ allora $|z|_{\mathbb{C}}=\sqrt{a^2+b^2}$.
\end{nota}

}

\end{frame}




\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

\smnsp

\begin{esempio}[Alcune norme vettoriali]
Se $V=\C^n$ o $V=\R^n$ e $x=(x_1,x_2,\ldots,x_n) \in V$ alcuni esempi di norme sono
\begin{itemize}
\item per $p \in [1,+\infty)$, la {\reds{norma}} $p$ \mye definita come 
\bc
\btc[width=8cm]
$${ \|x\|_p=\left( \sum_{j=1}^n |x_j|^p \right)^{1/p} }$$ 
\etc
\ec
\item la {\reds{norma del massimo}} \mye definita per $p=\infty$ come \bc
\btc[width=8cm]$${\|x\|_{\infty}=\max_{j=1,\ldots,n} |x_j| }$$ 
\etc
\ec
\end{itemize}
\end{esempio}

\smnsp

\begin{nota}
Di seguito, non distingueremo se lo spazio vettoriale $V$, sia uguale $\C^n$ o $V=\R^n$. Qualora necessario, per semplicit\myash, si ponga $V=\R^n$.
\end{nota}

\end{frame}




\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

\smnsp

\begin{nota}[Alcuni casi notevoli della norma $p$]
Si osservi che essendo ${ \|x\|_p=( \sum_{j=1}^n |x_j|^p )^{1/p} }$
\begin{itemize}
\item per $p=1$, si definisce la importante {\reds{norma}} $1$
\bc
\btc[width=6cm]
$${\|x\|_1=\sum_{j=1}^n |x_j|}$$ 
\etc
\ec
\item per $p=2$, si definisce la importante {\reds{norma}} $2$, detta anche  {\reds{euclidea}}
\bc
\btc[width=6cm]
$${\|x\|_2=\left( \sum_{j=1}^n |x_j|^2 \right)^{1/2}}$$ 
\etc
\ec
\end{itemize}
\end{nota}

\end{frame}




\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

\smnsp
\small{
\begin{esempio}
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{esempio}

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

I due vettori sono numericamente {\reds{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{(|-0.001|,|-0.001|)}=0.001$,
\end{itemize}
risultati che sottolineano che $u$ e $v$ sono vicini rispetto alle distanze definite dalle norme.
 \end{esempio}
}
\end{frame}




\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

\smnsp

{{
\begin{definizione}
Sia $\| \cdot \|$ una norma vettoriale. Si definisce {\reds{norma vettoriale indotta}} di una matrice in $\C^{n \times n}$ la funzione che mappa ogni matrice $A \in \C^{n \times n}$ nel valore
\bc
\btc[width=6cm]
$$
\|A\|=\sup_{x \neq 0} \frac{\|Ax\|}{\|x\|}
$$
\etc
\ec
\end{definizione}

\smnsp

\begin{nota}
Tali norme si definiscono anche nel caso particolare che $A \in \R^{n \times n}$.
\end{nota}

\smnsp

\begin{definizione}[Raggio spettrale]
Si definisce {\reds{raggio spettrale}} di $A $ la quantit\mya 
\bc
\btc[width=6cm]
$${\rho(A)=\max_{k=1,\ldots,n} (|\lambda_k|)}$$
\etc
\ec
dove $\lambda_k$ \mye un {\reds{autovalore}} di $A$, ovvero esiste $x_k \neq 0$ detto {\reds{autovettore}} di $A$, tale che $Ax_k=\lambda_k x_k$.
\end{definizione}
}}
\end{frame}




\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

\smnsp

{\fns{
\begin{esempio}[Alcune norme matriciali]
Posto $A=(a_{i,j})$, alcuni esempi di norme matriciali risultano:
\begin{itemize}
\item per $p=1$, partendo dalla norma vettoriale $\|\cdot\|_1$, si dimostra che la {\reds{norma}} $1$ matriciale risulta uguale al massimo della somma dei moduli delle componenti di ogni {\purples{colonna}} ovvero

\smnsp

\bc
\btc[width=6cm]
\smnsp
$${\|A\|_1=\max_{1 \leq j \leq n} \sum_{i=1}^n |a_{i,j}|};$$
\etc
\ec
\item per $p=2$, partendo dalla norma vettoriale $\|\cdot\|_2$, si dimostra che la {\reds{norma}} $2$ matriciale risulta uguale a 

\smnsp

\bc
\btc[width=6cm]
$${\|A\|_2=\sqrt{\rho(A^* A)}};$$
\etc
\ec  

\smnsp

dove 
\begin{itemize}
\item se $A=(a_{i,j}) \in \R^{n \times n}$ allora $A^*$ \mye la trasposta $A^T$ di $A$, ovvero $a^*_{i,j}=a_{j,i}$,
\item se $A \in \C^{n \times n}$ allora $A^*$ \mye la trasposta coniugata $A^H$ di $A$ ovvero $a^*_{i,j}={\overline{a_{j,i}}}$;
\end{itemize}
\item per $p=\infty$, partendo dalla norma vettoriale $\|\cdot\|_{\infty}$, si dimostra che la {\reds{norma}} $\infty$ matriciale risulta uguale al massimo della somma dei moduli delle componenti di ogni {\purples{riga}} ovvero

\smnsp

\bc
\btc[width=6cm]
\smnsp
$${\|A\|_{\infty}=\max_{1 \leq i \leq n} \sum_{j=1}^n |a_{i,j}|}$$
\smnsp
\etc
\ec
\end{itemize}
\end{esempio}
}}
\end{frame}



\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

\smnsp

{{
\begin{esempio}[1]
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{esempio}
}}
\end{frame}




\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

{{
\begin{esempio}[2]
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 {\reds{vicine}}. Per verificarlo calcoliamo la loro distanza in termine di norme ovvero {\reds{$
{\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=\max(0.003,0.006)=0.006$,
\item $\|A-B\|_2=0.005415654779058332$,
\item $\|A-B\|_{\infty}=\max(0.002,0.007)=0.007$,
\end{itemize}
risultati che sottolineano che $A$ e $B$ sono matrici vicine rispetto alle distanze definite dalle norme.
 \end{esempio}
}}
\end{frame}





\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

\smnsp

\begin{teorema}[Disuguaglianze norme indotte]
Per le norme matriciali indotte, se $A,B \in \C^{n \times n}$, $x \in \C^{n \times 1}$,

\bc
\btc[width=6cm]
\begin{equation}\label{dis1}
\|Ax\| \leq \|A\|\|x\|,
\end{equation}
\etc
\ec

\smnsp

\bc
\btc[width=6cm]
\begin{equation}\label{dis2}
\|AB\| \leq \|A\|\|B\|.
\end{equation}
\etc
\ec


\end{teorema}

\smnsp

{\fns{
\begin{dimostrazione}[Facoltativa]
\begin{enumerate}
\item La prima disuguaglianza \'e banalmente verificata per $x=0$. Altrimenti, per $x \neq 0$,
$$
\|A\|=\sup_{y \neq 0} \frac{\|Ay\|}{\|y\|} \geq \frac{\|Ax\|}{\|x\|} \Rightarrow \|Ax\| \leq \|A\|\|x\|.
$$

\item Per quanto riguarda la seconda disuguaglianza, dal primo punto, se $y \neq 0$, allora $$
\|ABy\|=\|A(By)\| \leq \|A\| \|By\| \leq \|A\| \|B\|\|y\| 
$$ da cui $ {\|ABy\|}/{\|y\|} \leq \|A\| \|B\|$
e quindi 
$$
\|AB\|=\sup_{y \neq 0} \frac{\|ABy\|}{\|y\|}  \leq \|A\| \|B\|.
$$
\end{enumerate}
\end{dimostrazione}
}}
\end{frame}





\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare: Indice di condizionamento}

\smnsp



\begin{definizione}[Indice di condizionamento]
Sia $A \in \C^{n \times n}$ una matrice invertibile ossia tale che $\det(A) \neq 0$. Si dice {\reds{indice di condizionamento}} della matrice $A$, relativamente alla norma indotta $\|\cdot\|$,
\bc
\btc[width=6cm]
$${
k(A):=\|A\|\|A^{-1}\|.
}
$$
\etc
\ec
\end{definizione}

\begin{nota}
Si osservi che dalla definizione, se $I $ \mye la matrice identica di $\C^{n \times n}$
$$
\|I\|=\sup_{x \neq 0} \frac{\|Ix\|}{\|x\|}=\sup_{x \neq 0} \frac{\|x\|}{\|x\|}=1
$$
e quindi, se la norma $\|\cdot\|$ \mye indotta, essendo $\|ST\| \leq \|S\|\|T\| $ per ogni $S,T \in \C^{n \times n}$, ricaviamo
$$
1=\|I\|=\|A A^{-1}\|\leq \|A\| \|A^{-1}\|=k(A)
$$
ovvero {\reds{$k(A) \geq 1$}}.
\end{nota}

\end{frame}




\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

\smnsp

Valgono i seguenti risultati

\begin{teorema}[Stime dal basso dell'indice di condizionamento]
\begin{enumerate}
\item Qualsiasi sia la norma matriciale indotta, e $A$ invertibile,
\bc
\btc[width=6cm]
$$
1 \leq \frac{|\lambda_{\mbox{max}}|}{|\lambda_{\mbox{min}}|} \leq k(A)
$$
\etc
\ec
dove $\lambda_{\mbox{min}}$ e $\lambda_{\mbox{max}}$ sono rispettivamente gli autovalori di $A$ di minimo e massimo modulo.

\item Se $A$ \mye {\purples{simmetrica}}, ovvero $A=A^T$, 

\bc
\btc[width=6cm]
$$
k_2(A)=\frac{|\lambda_{\mbox{max}}|}{|\lambda_{\mbox{min}}|}.
$$
\etc
\ec
\end{enumerate}
\end{teorema}

\end{frame}






\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

\smnsp

\begin{esempio}[Matrice di Hilbert]
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{esempio}


\end{frame}




\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

\smnsp

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

\begin{teorema}[Condizionamento ed errore nel termine noto]
Se 
\begin{enumerate}
\item $A \in \C^{n \times n}$ \mye una matrice invertibile, 
\item $b \in \C^{n}$, $b \neq 0$,
\item $Ax=b$,
\item $A(x+\delta x)={\reds{b+\delta b}}$,
\item $\|\cdot\|$ \mye una norma indotta,
\end{enumerate}
allora vale la stima
\bc
\btc[width=6cm]
$$
{
\frac {\|\delta x\|}  {\|x\|} \leq k(A) \frac {\|\delta b\|}  {\|b\|}
}
$$
\etc
\ec
dove $k(A)=\|A\|\|A^{-1}\|$.
\end{teorema}

\end{frame}



\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

\smnsp
{\small{
\begin{dimostrazione}[Facoltativa]
Da $$A(x+\delta x)=b+\delta b, \quad Ax=b$$ abbiamo che
$$
A\delta x=A((x+\delta x)-x)=A(x+\delta x)-Ax=(b+\delta b)-b=\delta b
\Leftrightarrow {\reds{\delta x=A^{-1}\delta b}}.$$
Da (\ref{dis1}), $\|Su\| \leq \|S\|\|u\|$, per ogni $S \in \C^{n \times n}$, $u \in C^{n}$ e quindi per $S=A^{-1}$, $u=\delta b$ si ha $\|A^{-1}\delta b\| {{\leq \|A^{-1}\| \|\delta b\|}}$ da cui 
\begin{equation}\label{theo1.1}
0 \leq {\reds{\|\delta x\|}}=\|A^{-1}\delta b\| {\reds{\leq \|A^{-1}\| \|\delta b\|}}.
\end{equation}
D'altra parte, essendo $Ax=b$, necessariamente per  (\ref{dis1})
{{$$
\|b\|=\|Ax\| \leq \|A\|\|x\|
$$}}
e quindi, essendo $x \neq 0$ in quanto $b \neq 0$, moltiplicando i membri di $0 \leq \|b\| \leq \|A\|\|x\|$ per $\frac{1}{\|x\|\|b\|}$
\begin{equation}\label{theo1.2}
0 \leq {\reds{\frac{1}{\|x\|} \leq \frac{\|A\|}{\|b\|} }}.
\end{equation}
Da (\ref{theo1.1}), (\ref{theo1.2}), visto che se $0 \leq s \leq u$, $0 \leq t \leq v$ allora $0 \leq s t \leq uv$  e $k(A)=\|A^{-1}\|\|A\|$, ricaviamo 
$$
\frac {\|\delta x\|}{\|x\|}= {\|\delta x\|} \cdot \frac {1}{\|x\|}  \leq \|A^{-1}\|\|\delta b\|   \frac{\|A\|}{\|b\|}=\|A^{-1}\|\|A\| \frac{\|\delta b\|}{\|b\|}=k(A) \frac{\|\delta b\|}{\|b\|}.
$$
\end{dimostrazione}
}}
\end{frame}





\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

\smnsp

\begin{commento}[Relazione tra $\|\delta x\|/\|x\|$ e $\|\delta b\|/\|b\|$]
\begin{itemize}
\item Questo teorema mostra che si fa un errore relativo sul termine noto pari a $\frac{\|\delta b\|}{\|b\|}$ allora si compie un errore relativo tra soluzione $x$ e soluzione $x+\delta x$ del problema perturbato 
$$\frac {\|x-(x+\delta x) \|}{\|x\|}=\frac {\|\delta x\|}{\|x\|}$$
 tale che
$$
\frac {\|\delta x\|}{\|x\|}\leq k(A) \frac{\|\delta b\|}{\|b\|}.
$$

\item Quindi, se $k(A) \gg 1$, pu\myo accadere che la soluzione  $x+\delta x$  del problema perturbato sia molto distante da $x$ e di conseguenza la soluzione del sistema lineare $Ax=b$ molto suscettibile a piccoli errori sul termine noto.

\item In altri termini {\reds{pi\myu grande \mye $k(A)$}} pi\myu i sistemi lineari $Ax=b$ sono {\reds{difficili da trattare}} nel senso che piccoli errori sui dati $b$ possono portare a grossi errori sulla soluzione $x$.
\end{itemize}
\end{commento}

\smsp

\end{frame}



\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

\smnsp
{\small{
\begin{esempio}
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 $A(x+\delta x)=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$. 

\smsp

Questo significa che a un piccolo errore relativo sui dati abbiamo riscontrato un eccessivo errore relativo sulla soluzione.

\smsp

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



\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

\smnsp


\begin{teorema}[Condizionamento e errore sui dati della matrice]
Se 
\begin{enumerate}
\item $A \in \C^{n \times n}$ \mye una matrice invertibile, 
\item $b \in \C^{n}$, $b \neq 0$,
\item $Ax=b$,
\item ${\reds{(A+\delta A)}}(x+\delta x)=b$, con $\det(A+\delta A) \neq 0$,
\end{enumerate}
allora per una qualsiasi norma indotta $\|\cdot\|$ vale la stima
\bc
\btc[width=7cm]
$$
\frac {\|\delta x\|}  {\|x+\delta x\|} \leq k(A) \frac {\|\delta A\|}  {\|A\|}.
$$
\etc
\ec
dove $k(A)=\|A\| \|A^{-1}\|$.
\end{teorema}

\begin{commento}
Il precedente teorema dice che in presenza di un forte indice di condizionamento, piccoli errori relativi sulla matrice possono indurre grandi errori relativi nella soluzione.  Si noti che usualmente $\|x\| \approx \|x+\delta x\|$.
\end{commento}

\end{frame}




\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

\smnsp
{{
\begin{dimostrazione}[Facoltativa]
Essendo 
\begin{itemize}
\item $(A+\delta A)(x+\delta x)=b$, 
\item $Ax=b$,
\end{itemize}
 necessariamente
\begin{eqnarray}
{\reds{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={\reds{A\delta x+\delta A (x+\delta x)}} \nonumber
\end{eqnarray}
da cui {\reds{$
A \delta x=-\delta A (x + \delta x)
$}}
ovvero, moltiplicando ambo i membri per $A^{-1}$,
{{$$
\delta x=-A^{-1}\delta A (x + \delta x)
$$}}
e per (\ref{dis1}), (\ref{dis2}), da {\purples{$\|Tu\| \leq \|T\|\|u\|$}}, {\purples{$\|ST\| \leq \|S\| \|T\|$}} per ogni $S,T \in {\mathbb{C}}^{n \times n}$, $u \in {\mathbb{C}}^n$, 
\begin{equation}\label{3242}
{\reds{\|\delta x\|}}=\|-A^{-1}\delta A (x + \delta x)\| \leq\|A^{-1}\delta A\|\| x + \delta x\| {\reds{\leq \|A^{-1}\|\|\delta A\|\| x + \delta x\|}}.
\end{equation}
Essendo $\|A\| \neq 0$, dividendo ambo i membri di (\ref{3242}) per $\|x+\delta x\|$, da $ k(A)=\|A\| \|A^{-1}\|$
$$
\frac {\|\delta x\|}{\| x + \delta x\|} 
\leq \|A^{-1}\|\|\delta A\|=
\frac {\|A^{-1}\| \|A\|} {\|A\|}   \|\delta A\| = k(A) \frac {\|\delta A\|}{ \|A\|}.
$$
\end{dimostrazione}
}}
\end{frame}












\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare}

\smnsp
{\small{

\begin{esempio}
Siano
$$
A=\left[
\begin{array}{ccc}
$1.0000 $&$   0.5000 $&$   0.3333$ \\
$0.5000 $&$   0.3333  $&$  0.2500$ \\
$0.3333 $&$   0.2500  $&$  0.2000$ 
\end{array} 
\right],
A+\delta A=\left[
\begin{array}{ccc}
$1.0010 $&$   0.5030  $&$  0.3403$ \\
 $   0.5040  $&$  0.3383  $&$  0.2520$ \\
 $   0.3423  $&$  0.2520 $&$   0.2050$ 
\end{array} 
\right],
$$
Il sistema $Ax=b$, per $b=[1, \, 2, \, 3]^T$ ha soluzione approssimativamente
$$
x\approx \left[
\begin{array}{c}
$27.3232$ \\  $-193.6572$ \\ $211.5374$
\end{array} 
\right]
$$
mentre $(A+\delta A)(x+\delta x)=b$ ha soluzione
$$
x+\delta x\approx \left[
\begin{array}{c}
$5.5489 $ \\  $ -75.3599  $ \\  $ 98.0063$ \\ 
\end{array} 
\right].
$$

Di conseguenza seppure $\frac{\|\delta A\|_{2}}{\|A\|_{2}}=9.35310e-03$ abbiamo $\frac{\|\delta x\|_{2}}{\|x+\delta x\|_{2}}=1.33653e+00$. Questo significa che a un piccolo errore relativo sui dati abbiamo riscontrato un eccessivo errore relativo sulla soluzione.

\smsp

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

\end{esempio}
}}
\end{frame}













\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare: facoltativo}

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

\begin{teorema}[Facoltativo]\label{thm_3}
Se 
\begin{enumerate}
\item $A \in \C^{n \times n}$ \mye una matrice invertibile, 
\item $b \in \C^{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 $\|\cdot\|$ \mye una norma indotta,
\item $k(A) \|\delta A\| < \|A\|$,
\end{enumerate}
allora vale la stima
\bc
\btc[width=10cm]
$$
\frac {\|\delta x\|}  {\|x\|} \leq \frac{k(A)} { 1-k(A) \|\delta A\|/\|A\|}    \left( \frac {\|\delta A\|}  {\|A\|} +\frac {\|\delta b\|}  {\|b\|} \right).
$$
\etc
\ec
\end{teorema}
\end{frame}




\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare: facoltativo}

\smnsp

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

\smsp

Il valore 
$$
\mu(A):=\frac{k(A)} { 1-k(A) \|\delta A\|/\|A\|}
$$
\mye la costante di amplificazione degli errori.

\smsp

Visto che il denominatore \mye minore di $1$ (in effetti $k(A) \|\delta A\| < \|A\|$ porge $k(A) \|\delta A\|/\|A\| < 1$), deduciamo che se $k(A) \gg 1$, allora $\mu(A) \gg 1$ e quindi pu\'o essere $\frac{\| \delta x\|}{\|x\|} \gg 1 $, ovvero ${\| \delta x\|}\gg {\|x\|}$ per cui 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{commento}

\end{frame}

\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare: facoltativo}

\smnsp

{\footnotesize{
\begin{teorema}[Facoltativo]\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\|}, {\mbox{     }}
\|(I+A)^{-1}\| \leq \frac{1}{1-\|A\|}.
$$
\end{teorema}

\smnsp

\begin{dimostrazione}[Facoltativo]
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 = b+A\delta x+\delta A x +\delta A\delta x \nonumber 
\end{eqnarray}
da cui sottraendo $b$ ad ambo i membri {\reds{$
\delta b=A \delta x+\delta A x +\delta A \delta x=\delta A x + (A+\delta A) \delta x
$}}
ovvero $(A+\delta A) \delta x=\delta b-\delta A x$.

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\|=\frac{ \|A^{-1}\| \|A\|} {\|A\|}\| \delta A\|=k(A) \frac{\| \delta A\|}{\|  A\|} 
 \leq  \frac{\|  A\|}{\|  A\|} =1 \nonumber 
\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$).
\end{dimostrazione}
}}
\end{frame}



\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare: facoltativo}

\smnsp

{\footnotesize{
\begin{vuoto}
Da un teorema precedente, moltiplicando ambo i membri per $(A+\delta A)^{-1}$, visto che $(MN)^{-1}=N^{-1}M^{-1}$ 

\smnsp

\begin{eqnarray}
\delta x &=&(A+\delta A)^{-1} (\delta b-\delta A x)=(A(I+B))^{-1} (\delta b-\delta A x) = (I+B)^{-1} A^{-1}  (\delta b-\delta A x) \nonumber 
\end{eqnarray}

da cui, essendo $\| (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) \| \leq \| (I+B)^{-1}\| \| A^{-1} \| \| \delta b-\delta A x \| \nonumber \\
&\leq& \frac{1}{1-\|B\|} \| A^{-1} \| \| \delta b-\delta A x \|  =  \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 \|) \nonumber
\end{eqnarray}
e quindi riassumendo

\smnsp

\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) \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) \leq     \frac{ k(A)  }{1-k(A) \| \frac{\|\delta A\|}{\|A\|}}     \left(\frac{ \| \delta b\|}{\|b\|}+\frac{\|\delta A\|}{\|A\|} \right).
\end{eqnarray}
\end{vuoto}
}}
\end{frame}




\begin{frame}[fragile]
\frametitle{Cenni di algebra lineare: facoltativo}

\smnsp

{\footnotesize{
\begin{esempio}[Facoltativo]
Il sistema $Ax=b$, dove

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

\smnsp

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

\smnsp

\begin{itemize}
\item $(A+\delta A)u=b$ ha soluzione $u=(3,0)^T$ (perturbata matrice $A$), 
\item  $Av=b+\delta b$ ha soluzione $v=(2,0.5)^T$  (perturbato termine noto $b$),
\item $(A+\delta A)w=(b+\delta b)$ ha soluzione $w=(5,-1)^T$ (perturbati $A$ e $b$).
\end{itemize}
Ne risulta che per possibili piccole perturbazioni alla matrice e/o al termine noto abbiamo ottenuto nei diversi casi una soluzione molto {\reds{distante}} dall'originaria. Si ha
\begin{itemize}
\item $\det(A) \approx 0.003$, $\det(\delta A) \approx 0.0015$, $k_2(A) \approx 2083.666853410356$,
\item abbiamo $
 { k_2(A)  }/{(1-k_2(A) \| {\|\delta A\|_2}/{\|A\|_2})} \approx 3.05387e+04
$;
\item $\frac{ \| \delta x\|_2}{\|x\|_2} \approx 3.16228e+00$, $\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

\smnsp

\smnsp

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

\end{esempio}
}}
\end{frame}




\begin{frame}[fragile]
\frametitle{Metodo di eliminazione di Gauss e fattorizzazione LU}

\smnsp
{\small{

\begin{proposito}
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$.
\end{proposito}

Per semplificare l'esposizione cominciamo con un esempio.

\begin{esempio}
Si risolva 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.
$$
mediante eliminazione gaussiana.
\end{esempio}
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 {\reds{se moltiplichiamo la $i$-sima equazione per una costante non nulla e sottraiamo il risultato membro a membro alla $j$-sima}}, ricaviamo una equazione che sostituita alla $j$-sima, propone un nuovo problema lineare che ha le stesse soluzioni del precedente.
\end{itemize}
}}
\end{frame}




\begin{frame}[fragile]
\frametitle{Metodo di eliminazione di Gauss e fattorizzazione LU}
{\footnotesize{
\begin{svolgimento}
Diremo che due sistemi sono equivalenti, se hanno le stesse soluzioni. Da:
$$
\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.
$$

\begin{itemize}
\item moltiplichiamo la {\purples{prima equazione}} per {\blues{$m_{2,1}=6/3=2$}} e la sottraiamo alla seconda, ottenendo il sistema lineare equivalente
 $$
\left\{
\begin{array}{l}
3x_1-2x_2-x_3=0 \\
{\reds{{0}}} \cdot x_1+2x_2+4x_3=6 \\
-9x_1+7x_2+x_3=-1 \\
\end{array}
\right.
$$
\item moltiplichiamo la {\purples{prima equazione}} per {\blues{$m_{3,1}=-9/3=-3$}} e la sottraiamo alla terza, ottenendo il sistema lineare equivalente
 $$
\left\{
\begin{array}{l}
3x_1-2x_2-x_3=0 \\
{\magentas{0}} x_1+2 x_2+4x_3=6 \\
{\reds{{0}}} x_1+1 x_2-2x_3=-1 \\
\end{array}
\right.
$$
in cui le componenti nulle della prima {\it{colonna}}, precedentemente ottenute, non vengono modificate.
\item moltiplichiamo la {\purples{seconda equazione}} per {\blues{$m_{3,2}=1/2$}} e la sottraiamo alla terza, ottenendo il sistema lineare equivalente
 $$
\left\{
\begin{array}{l}
3x_1-2x_2-x_3=0 \\
{\magentas{0}} x_1+2 x_2+4x_3=6 \\
{\magentas{0}} x_1+{\reds{{0}}} x_2-4x_3=-4 \\
\end{array}
\right.
$$
in cui le componenti nulle della {\it{colonna}} precedente non vengono modificate.
\end{itemize}
\end{svolgimento}
}}
\end{frame}




\begin{frame}[fragile]
\frametitle{Metodo di eliminazione di Gauss e fattorizzazione LU}
{\small{
\begin{vuoto}
Il sistema finale
 $$
\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.
$$
 \mye facile da risolvere. 

\begin{itemize}
\item L'ultima equazione {\purples{$-4x_3=-4$}} ha una sola variabile e una sola incognita e comporta che $$x_3=1.$$
\item Inseriamo {\purples{$x_3=1$}} nella penultima equazione {\purples{$2 x_2+4x_3=6$}} e otteniamo 
$$ 
2 x_2+4\cdot 1=6
$$
da cui $2x_2=2$ e quindi $x_2=1$.
\item Inseriamo {\purples{$x_2=x_3=1$}}  nella prima equazione {\purples{$3x_1-2x_2-x_3=0 $}} e otteniamo 
$$
3x_1-2 \cdot 1-1=0
$$
da cui $$x_1=1.$$
\end{itemize}
Quindi il sistema ha soluzioni $x_1=1$, $x_2=1$, $x_3=1$.
\end{vuoto}
}}
\end{frame}




\begin{frame}[fragile]
\frametitle{Metodo di eliminazione di Gauss e fattorizzazione LU}
{\small{
\begin{nota}[Sistema matriciale]
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 $Ux=c$ dove
$$
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 {\reds{semplici}}. 

\smsp

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 {\reds{sostituzione all'indietro}}.
% % \end{esempio}
\end{nota}
}}
\end{frame}




\begin{frame}[fragile]
\frametitle{Sistemi lineari triangolari}

\smnsp
\small{
\begin{definizione}
Una matrice $A=(a_{i,j}) \in \C^{n \times n}$ si dice 
\begin{itemize}
\item {\reds{triangolare inferiore}} se $a_{i,j}=0$ qualora $i < j$;
\item {\reds{triangolare superiore}} se $a_{i,j}=0$ qualora $i > j$.
\end{itemize}
\end{definizione}

\begin{esempio}
La matrice
$$
L=\left[
\begin{array}{c c c}
3 & 0 & 0\\
4 & 2 & 0 \\
1 & 0 & -4 \\
\end{array}
\right]
$$
\mye triangolare inferiore ($a_{i,j}=0$ se l'indice di riga strettamente minore di quello di colonna).
\end{esempio}


\begin{esempio}
La matrice
$$
U=\left[
\begin{array}{c c c}
3 & -2 & -1 \\
0 & 2 & 4 \\
0 & 0 & -4 \\
\end{array}
\right]
$$
\mye triangolare superiore ($a_{i,j}=0$ se l'indice di riga strettamente maggiore di quello di colonna).
\end{esempio}
}
\end{frame}




\begin{frame}[fragile]
\frametitle{Sistemi lineari triangolari}

\smnsp
\begin{algoritmo}[Sostituzione all'indietro e in avanti]
Sia $b=(b_j) \in \C^{n}$ e supponiamo di dover risolvere il sistema $Ax=b$, con {\purples{$A_{k,k} \neq 0$}} per $k=1,\ldots,n$. 

\smsp

In tal caso, 
\begin{itemize}
\item se $A$ \mye triangolare superiore e invertibile, si usa l'algoritmo della {\reds{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}
\item se $A$ \mye triangolare inferiore e invertibile, si usa l'algoritmo della {\reds{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}
\end{itemize}
\end{algoritmo}
\end{frame}




\begin{frame}[fragile]
\frametitle{Sistemi lineari triangolari}

\smnsp

\begin{nota}[Complessit\'a risoluzione sistemi triangolari con eliminazione gaussiana]
Si osservi che la {\reds{risoluzione di sistemi lineari triangolari}} necessita approssimativamente {\reds{$n^2/2$ operazioni moltiplicative}}, dove $n$ \mye l'ordine della matrice.

\smsp

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

\end{nota}

\end{frame}




\begin{frame}[fragile]
\frametitle{Metodo di eliminazione gaussiana senza pivoting}

\smnsp


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

\smsp

Posti
\begin{itemize}
\item $A=A^{(1)}=(a^{(1)}_{i,j})$, 
\item $b=b^{(1)}=(b^{(1)}_j)$, 
\end{itemize}
possiamo immagazzinare le componenti che lo descrivono nella matrice {\reds{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].
$$

\end{frame}



\begin{frame}[fragile]
\frametitle{Metodo di eliminazione gaussiana senza pivoting}


{\footnotesize{
{\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 {\reds{moltiplicatore}}
\bc
\btc[width=8cm]
$$
m_{i,1}=\frac{a^{(1)}_{i,1}}{a^{(1)}_{1,1}}, i=2,\ldots,n
$$
\etc
\ec
e in successione moltiplichiamo la {\purples{prima riga}} per $-m_{i,1}$ e la sommiamo l'i-sima, 
$$
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)}_{1}, i=2,\ldots,n.
$$
Si vede 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].
$$


}}


\end{frame}



\begin{frame}[fragile]
\frametitle{Metodo di eliminazione gaussiana senza pivoting}

{\vspace{-0.1cm}}

{\footnotesize{

{\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'{\reds{elemento pivotale}} {\purples{$a^{(k)}_{k,k} \neq 0$}} definiamo il moltiplicatore

\bc
\btc[width=7cm]
\smnsp
$$
m_{i,k}=\frac {a^{(k)}_{i,k}} {a^{(k)}_{k,k}}, \,\,\, i=k+1,\ldots,n.
$$
\smnsp
\etc
\ec

moltiplichiamo la $k$-sima riga per $-m_{i,k}$ e la sommiamo alla $i$-sima, con $i=k+1,\ldots,n$, ovvero

\bc
\btc[width=7cm]
$$
a^{(k+1)}_{i,j}=a^{(k)}_{i,j}-m_{i,k}a^{(k)}_{k,j},\,\,\, j=k+1,\ldots,n,
$$
\smnsp
\etc
\ec


\bc
\btc[width=7cm]
$$
b^{(k+1)}_{i}=b^{(k)}_{i}-m_{i,k}b^{(k)}_{k},\,\,\, j=k+1,\ldots,n.
$$
\smnsp
\etc
\ec
Se $i \leq k$, per $j=1,\ldots,n$, poniamo $a^{(k+1)}_{i,j}=a^{(k)}_{i,j}$, $b^{(k+1)}_{i}=b^{(k)}_{i}$.
}}


\end{frame}



\begin{frame}[fragile]
\frametitle{Metodo di eliminazione gaussiana senza pivoting}

\smnsp
{\small{
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].
$$

\smsp

Con questo procedimento, dopo $n-1$ passi otteniamo un sistema $A^{(n)}x=b^{(n)}$ in cui la matrice $A^{(n)}$ \mye triangolare superiore con elementi diagonali $A^{(n)}_{k,k} \neq 0$, $k=1,\ldots,n$.

\smsp

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

}}
\end{frame}


\begin{frame}[fragile]
\frametitle{Metodo di eliminazione gaussiana senza pivoting}

\smnsp
{\fns{
\begin{definizione}[Fattorizzazione LU]
Una matrice $A \in \C^{n \times n}$ \mye fattorizzabile $LU$ se esistono
\begin{itemize}
\item $L=(l_{i,j}) \in \C^{n \times n}$ triangolare inferiore con elementi diagonali $l_{i,i}=1$, 
\item $U=(u_{i,j}) \in \C^{n \times n}$ triangolare superiore, 
\end{itemize}
tali che $A=LU$ (cf. \cite{WIKILU}).
\end{definizione}

\begin{importante}[Pivot non nulli]
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)}x=b^{(n)}$ 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)}$,
\item $L$ 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 & 0\\
m_{n,1} &m_{n,2} & \ldots & \ldots & m_{n,n-1} & 1 \\
\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.  
\end{importante}
}}
\end{frame}



\begin{frame}[fragile]
\frametitle{Metodo di eliminazione gaussiana senza pivoting}

\smnsp
{\fns{
\begin{esempio}
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 dell'eliminazione gaussiana senza pivoting abbiamo determinato $U=A^{(n)}$  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 $m_{2,1}=2$, $m_{3,1}=-3$, $m_{2,2}=0.5$,
$$
L
=
\left[
\begin{array}{c c c}
1 & 0 & 0 \\
m_{2,1} & 1 & 0 \\
m_{3,1} & m_{3,2} & 1 \\
\end{array}
\right]
=
\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. In Matlab:

\smnsp

\begin{verbatim}
>> A=[3 -2 -1; 6 -2 2; -9 7 1]; 
>> U=[3 -2 -1; 0 2 4; 0 0 -4]; L=[1 0 0; 2 1 0; -3 0.5 1];
>> L*U  % verifichiamo che A=L*U
ans =
     3    -2    -1
     6    -2     2
    -9     7     1
>> norm(A-L*U) % la norma di una matrice C nulla se e solo se C e' nulla
ans =
     0
>>
\end{verbatim}
\end{esempio}
}}
\end{frame}



\begin{frame}[fragile]
\frametitle{Metodo di eliminazione gaussiana con pivoting (parziale)}

\smnsp

\begin{commento}
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 
\bc
\btc[width=8cm]
\smnsp
$$
m_{i,k}=\frac {a^{(k)}_{i,k}} {a^{(k)}_{k,k}}, \,\,\, k=1,\ldots,n-1, \,\,\, i=k+1,\ldots,n
$$
\etc
\ec
e di seguito ridurre il problema iniziale alla risoluzione di un sistema $A^{(n)}x=b^{(n)}$, con $A^{(n)}$ triangolare superiore.

\smsp

Possono sorgere i seguenti problemi: 
\begin{itemize}
\item qualche {\reds{${a^{(k)}_{k,k}}$ risulta nullo}} e quindi il processo precedente \mye  inapplicabile;
\item qualche {\reds{${a^{(k)}_{k,k}}$ risulta 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}}).
\end{itemize}

La soluzione risulta la seguente, detta del {\reds{pivoting (parziale)}} o del {\reds{massimo pivot}}. 
\end{commento}
\end{frame}



\begin{frame}[fragile]
\frametitle{Metodo di eliminazione gaussiana con pivoting (parziale)}

\smnsp
\begin{vuoto}
Supponiamo di voler {\it{azzerare}} le componenti sotto la diagonale della colonna k-sima e determiniamo
\bc
\btc[width=6cm]
$$
c_k=\max_{k \leq i \leq n} |a^{(k)}_{i,k}|.
$$
\etc
\ec

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.

\smsp
\begin{itemize}
\item Sia {\reds{$i_k:=\phi(k) \geq k$}} il pi\myu piccolo indice di riga, per cui il massimo $c_k$ \mye ottenuto. 
\item Se $i_k > k$, {\reds{scambiamo la riga $i_k$-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. 

\smsp

Con questa tecnica, la nuova matrice aumentata avr\mya il pivot pi\myu grande possibile in modulo, con vantaggi per la propagazione degli errori.
\end{itemize}
\end{vuoto}
\end{frame}




\begin{frame}[fragile]
\frametitle{Metodo di eliminazione gaussiana con pivoting (parziale)}

\begin{definizione}[Matrice di permutazione]
Una matrice $P=(p_{i,j})$ si dice di {\reds{\it{permutazione}}} se le sue componenti sono $0$ o $1$, e ogni riga e colonna contengono una sola componente uguale a $1$.
\end{definizione}
Se calcoliamo $y=Px$ il risultato sara' un vettore con le stesse componenti di $x$ opportunamente {\it{scambiate}} (in matematica si usa dire {\it{permutate}}).

\begin{esempio}
Consideriamo la matrice di permutazione
$$
P=\left[
\begin{array}{c c c}
0 & 1 & 0 \\
1 & 0 & 0 \\
0 & 0 & 1 \\
\end{array}
\right].
$$
e il vettore $x=(x_1,x_2,x_3)^T$. Allora 
$$
\left[
\begin{array}{c c c}
0 & 1 & 0 \\
1 & 0 & 0 \\
0 & 0 & 1 \\
\end{array}
\right].
\left(
\begin{array}{c}
{\blues{x_1}}\\
{\purples{x_2 }}\\
x_3 \\
\end{array}
\right)
=
\left(
\begin{array}{c}
0\cdot x_1+1 \cdot x_2+0\cdot x_3\\
1\cdot x_1+0 \cdot x_2+0\cdot x_3\\
0\cdot x_1+0 \cdot x_2+1\cdot x_3\\
\end{array}
\right)
=
\left(
\begin{array}{c}
{\purples{x_2}}\\
{\blues{x_1}} \\
x_3 \\
\end{array}
\right).
$$
 
\end{esempio}
\end{frame}





\begin{frame}[fragile]
\frametitle{Metodo di eliminazione gaussiana con pivoting (parziale)}

\smnsp
{\footnotesize{
\begin{importante}
Dall'applicazione della procedura dell'eliminazione gaussiana con pivoting per la risoluzione di un sistema lineare, abbiamo ottenuto un sistema $A^{(n)}x=b^{(n)}$ equivalente all'originario.
Poniamo 
\begin{itemize}
\item $P$ matrice di {\purples{permutazione}} ottenuta come $P=P_{n-1}\cdot \ldots \cdot P_1$ dove $P_k$ \mye la matrice di permutazione che la scambia la riga $i_k$ con la $k$-sima, ovvero la matrice identica eccetto per le righe {\reds{$i_k$}}, {\reds{$k$}} in cui 
\bc
\btc[width=8cm]
$$
P_{i_k,i_k}=0, P_{i_k,k}=1, P_{k,i_k}=1, P_{k,k}=0;
$$
\etc
\ec

\smnsp

\item $U=A^{(n)}$ di forma {\purples{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 {\purples{triangolare inferiore}}, con componenti diagonali {\reds{$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 indipendentemente dalla risoluzione di sistemi lineari.  
\end{importante}

\smnsp

\begin{nota}[storica]
L'algoritmo di eliminazione gaussiana \mye stato scoperto in Cina circa 2000 anni fa, con pi\myu recenti contributi in lavori di Lagrange e Jacobi. La descrizione moderna  \mye stata introdotta solo nel 1930. {\cite[p.607]{TREF}}.
\end{nota}


}}
\end{frame}





\begin{frame}[fragile]
\frametitle{Metodo di eliminazione gaussiana con pivoting (parziale)}



\smnsp
Per capire perch\'e si eseguono queste operazioni matriciali, producendo una sequenza di sistemi lineari equivalenti tra loro, mostriamo le idee su un esempio. Si user\'a un linguaggio matriciale, seppure si analizzi un sistema di equazioni. Per "riga" si intender\'a "equazione", per "colonna" si intender\'a "variabile".

\begin{esempio}
Si risolva  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.
$$
mediante il metodo di eliminazione gaussiana con pivoting (parziale).
\end{esempio}
Il {\purples{termine maggiore in modulo della prima colonna}} \mye il terzo e vale $|-9|=9$. Scriviamo $\phi(1)=3$, intendendo che nella prima colonna, la componente di massimo modulo \'e la terza. Scambiamo la prima riga, ovvero equazione, 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.
$$

\end{frame}



\begin{frame}[fragile]
\frametitle{Metodo di eliminazione gaussiana con pivoting (parziale)}

Applichiamo ora i primi passi dell'eliminazione gaussiana relativamente a

$$
\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  {\purples{prima riga}}, ovvero equazione, per {\blues{$m_{2,1}=(-6/9)=-2/3$}} e la sottraiamo 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  {\purples{prima riga}}, ovvero equazione, per {\blues{$m_{3,1}=(-1/3)=-1/3$}} e la sottraiamo 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.
$$
\end{itemize}
\end{frame}


\begin{frame}[fragile]
\frametitle{Metodo di eliminazione gaussiana con pivoting (parziale)}

\smnsp
{\small{

\begin{itemize}
\item Osserviamo che il {\purples{termine maggiore in modulo della seconda colonna}}, di indice maggiore o uguale di quello di colonna,  \'e quello della seconda riga visto che $\max({8/3,1/3})=8/3$. Con le notazioni introdotte, $\phi(2)=2$ e quindi non si richiedono scambi di equazioni (si dovrebbe scambiare la seconda equazione con la seconda equazione). 

Moltiplichiamo la  {\purples{seconda riga}} per {\blues{$m_{3,2}=(1/3)/(8/3)=1/8$}} e la sottraiamo 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 mediante sostituzione all'indietro. 
\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}
}}
\end{frame}



\begin{frame}[fragile]
\frametitle{Metodo di eliminazione gaussiana con pivoting (parziale)}

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

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)}$ e dunque
$$
U=\left[
\begin{array}{c c c}
-9 & 7 & 1 \\
0 & 8/3 & 8/3 \\
0 & 0 & -1 \\
\end{array}
\right].
$$
\end{frame}



\begin{frame}[fragile]
\frametitle{Metodo di eliminazione gaussiana con pivoting (parziale)}

\smnsp

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].
$$
Nel passo relativo alla prima colonna, avendo scambiato la prima con la terza riga (era {\purples{$\phi(1)=3$}}), la matrice $P_1$, detta di permutazione, risulta la matrice identica eccetto per le righe $1$, $3$ che sono nulle ad eccezione di $(P_1)_{1,3}=(P_1)_{3,1}=1$, da cui 
$$
P_1=\left[
\begin{array}{c c c}
0 & 0 & 1 \\
0 & 1 & 0 \\
1 & 0 & 0 \\
\end{array}
\right],
$$
mentre nel passo relativo alla seconda colonna, visto che non abbiamo effettuato scambi (in effetti era {\purples{$\phi(2)=2$}}) $$P_2=\left[
\begin{array}{c c c}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1 \\
\end{array}
\right]$$ (ovvero la matrice identica). Quindi {\reds{$P=P_2 \cdot P_1=I \cdot P_1=P_1$}} e $PA=LU$.


\end{frame}


\begin{frame}[fragile]
\frametitle{Metodo di eliminazione gaussiana con pivoting (parziale)}

\smnsp

Avendo calcolato la fattorizzazione $PA=LU$, verifichiamo i risultati in Matlab.

\begin{lstlisting}[basicstyle=\small]
>> P=[0 0 1; 0 1 0; 1 0 0];
>> A=[3 -2 -1; 6 -2 2; -9 7 1];
>> PA=P*A;
>> L=[1 0 0; -2/3 1 0; -1/3 1/8 1];
>> U=[-9 7 1; 0 8/3 8/3; 0 0 -1];
>> LU=L*U;
>> norm(PA-LU,2)
ans =
   4.9651e-16
>> 
\end{lstlisting}

\smnsp

\noindent Si afferma che $\|P*A-L*U\|_2 \approx 4.9651e-16$ ovvero che anche numericamente $P*A=L*U$ (la distanza euclidea matriciale tra $PA$ e $LU$ \'e molto piccola). 

\end{frame}







\begin{frame}[fragile]
\frametitle{Alcune note sulla eliminazione gaussiana e la fattorizzazione LU}

\smnsp

{\small{

\begin{definizione}[{\white{\href{https://it.wikipedia.org/wiki/Matrice_a_diagonale_dominante}{Predominanza diagonale}}}]
La matrice $A \in \C^{n \times n}$, $A=(a_{i,j})$, \mye a {\blues{\href{https://it.wikipedia.org/wiki/Matrice\_a\_diagonale\_dominante}{predominanza diagonale per righe}}} se per ogni 
$i=1,\ldots,n$ risulta 

\smnsp

$${ { |a_{i,i}|  \geq \sum_{j=1,j \neq i}^n |a_{i,j}| }}$$
e per almeno un indice $s$ si abbia

\smnsp

$$|a_{s,s}| >\sum_{j=1,j \neq s}^n |a_{s,j}|.$$

La matrice $A \in \C^{n \times n}$ \mye a  {\blues{\href{https://it.wikipedia.org/wiki/Matrice\_a\_diagonale\_dominante}{predominanza diagonale stretta per righe}}}  se per ogni 
$i=1,\ldots,n$ risulta 

\smnsp

$${ { |a_{i,i}|  > \sum_{j=1,j \neq i}^n |a_{i,j}| }}.$$

La matrice $A \in \C^{n \times n}$ \mye a {\blues{\href{https://it.wikipedia.org/wiki/Matrice\_a\_diagonale\_dominante}{predominanza diagonale (stretta) per colonne}}}  se $A^T$  \mye a predominanza diagonale (stretta) per righe.
\end{definizione}
}}
\end{frame}



\begin{frame}[fragile]
\frametitle{Alcune note sulla eliminazione gaussiana e la fattorizzazione LU}

\smnsp
\begin{nota}[Matrici che hanno fattorizzazione $A=LU$]
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$.
\end{nota}


\begin{nota}[Fattorizzazione $PA=LU$ anche per matrici non invertibili]
Abbiamo visto che una qualsiasi matrice invertibile possiede una fattorizzazione del tipo $PA=LU$. 

\bigskip

In realt\mya una matrice arbitraria possiede una fattorizzazione di questo tipo (cf. {\cite[p.145]{BINI}}).
\end{nota}

\end{frame}


\begin{frame}[fragile]
\frametitle{Alcune note sulla eliminazione gaussiana e la fattorizzazione LU}

\smnsp

\begin{esempio}
La matrice
$$
\left(
\begin{array}{c c c}
4 & -1 & 0 \\
-1 & 4 & -1 \\
0 & -1 & 4 \\
\end{array}
\right).
$$

\noindent \mye 
\begin{itemize}
\item a predominanza diagonale stretta sia per righe che per colonne.
\item  simmetrica in quanto $A=A^T$ e definita positiva avendo per autovalori 
$$
\lambda_1 \approx  5.4142,
\lambda_2 \approx 4.0000,
\lambda_3 \approx  2.5858
$$
A tal proposito ricordiamo che una matrice \'e {\blues{\href{https://it.wikipedia.org/wiki/Matrice\_definita\_positiva}{definita positiva}}} se e solo se ha tutti gli autovalori positivi.
\end{itemize}
\end{esempio}

\end{frame}



\begin{frame}[fragile]
\frametitle{Alcune note sulla eliminazione gaussiana e la fattorizzazione LU}

\smnsp
\begin{pseudocodice}
Le seguenti righe di codice implementano la prima parte dell'algoritmo di eliminazione gaussiana senza pivoting, cui deve seguire la risoluzione di un sistema triangolare superiore mediante sostituzione all'indietro.

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

Notiamo che al termine del processo la matrice $A$ coincide con $L$ e che il codice produce errore se durante il processo qualche {\purples{\tt{a(k,k)}}} \'e nullo.

\smsp

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$. 
\end{pseudocodice}
\end{frame}



\begin{frame}[fragile]
\frametitle{Alcune note sulla eliminazione gaussiana e la fattorizzazione LU}

\smnsp

\begin{teorema}[Complessit\mya fattorizzazione PA=LU]
Il numero di operazioni moltiplicative richiesto per effettuare la fattorizzazione PA=LU \mye al pi\myu dell'ordine di {\reds{$n^3/3$}}.
\end{teorema}

\smnsp

{\scriptsize{
\begin{dimostrazione}[Facoltativa]
Dal pseudocodice, si vede che il numero di operazioni moltiplicatorie necessarie per effettuare la parte di eliminazione gaussiana, precedente alla sostituzione, \mye

\smnsp
\smnsp

\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) =
2 \sum_{p=1}^{n-1} (n-p) + \sum_{p=1}^{n-1} (n-p)^2 = 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

\smnsp
\smnsp

\begin{eqnarray}
\frac{(n-1)^3}{3} &=& \int_0^{n-1} x^2 dx \leq \sum_{k=1}^{n-1} k^2 \leq  \int_1^{n} x^2 dx = \frac{n^3-1}{3}
\end{eqnarray}

\smnsp

da cui, utilizzando la definizione di asintotico e il teorema del confronto,
{\purples{$
\sum_{k=1}^{n-1} k^2 \sim  \frac{n^3}{3}.
$}}
\end{itemize}

In definitiva, 

\smnsp

$$
2 \sum_{k=1}^{n-1} k + \sum_{k=1}^{n-1} k^2 \approx 
{n^2} + \frac{n^3}{3} \approx  {\reds{\frac{n^3}{3}}}.
$$
\end{dimostrazione}
}}
\end{frame}



\begin{frame}[fragile]
\frametitle{Alcune note sulla eliminazione gaussiana e la fattorizzazione LU}

\smnsp

{\small{
\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}
}}
\end{frame}




\begin{frame}[fragile]
\frametitle{Alcune note sulla eliminazione gaussiana e la fattorizzazione LU}

\begin{nota}[Importante]
\small{
{\reds{1}} 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$ con l'eliminazione gaussiana necessita di circa $n^3/3$ operazioni moltiplicative (si ricordi che $O(n^3/3)+2O(n^2/2)=O(n^3/3)$).	

\bigskip
 
{\reds{2}} La complessit\mya della fattorizzazione $PA=LU$ \mye essenzialmente la stessa della eliminazione gaussiana, ovvero $ \frac{n^3}{3}$.

\bigskip

{\reds{3}} In generale l'algoritmo di eliminazione gaussiana per la risoluzione di un sistema lineare $Ax=b$ con $A$ invertibile, pu\'o essere implementato come segue
\begin{enumerate}
\item {\purples{Calcola la fattorizzazione $PA=LU$}}. Osserviamo che  essendo $P$  invertibile (si dimostra che $|\det(P)|=1$)
$$Ax=b \Leftrightarrow PAx=Pb \Leftrightarrow LUx=Pb \Leftrightarrow {{Ly=\tilde{b} {\mbox{ con }} Ux=y, \tilde{b}=Pb}}.$$
\item {\purples{Calcola $\tilde{b}=Pb$}}.
\item {\purples{Risolvi il sistema triangolare inferiore $Ly=\tilde{b}$}} con la sostituzione in avanti.
\item {\purples{Risolvi il sistema triangolare superiore $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{nota}



\end{frame}



\begin{frame}[fragile]
\frametitle{Alcune note sulla eliminazione gaussiana e la fattorizzazione LU}

\smnsp






\begin{nota}[Decomposizione di Cholesky]
Se la $A$ matrice \mye {\reds{simmetrica e definita positiva}}, ovvero 
\begin{itemize}
\item $A=A^T$, 
\item tutti gli autovalori di $A$ sono positivi,
\end{itemize}
 \mye pi\myu conveniente utilizzare la {\blues{
\href
{https://it.wikipedia.org/wiki/Decomposizione\_di\_Cholesky}
{decomposizione di Cholesky}
}} \cite[p.524]{ATK}, in cui $${\reds{A=L L^T}}$$ con {\reds{$L=(l_{i,j})$ triangolare inferiore con elementi diagonali $l_{k,k} > 0$}}, $k=1,\ldots,n$. 

\smsp

La determinazione di questa fattorizzazione richiede approssimativamente $ \frac{n^3}{6}$ operazioni moltiplicative.

\smsp

Di conseguenza, per calcolare la soluzione di $Ax=b$, basta risolvere i due sistemi triangolari $Ly=b$, $L^Tx=y$. Quindi in generale la soluzione del sistema lineare si pu\myo ottenere in
$$
O \left( \frac{n^3}{6}  \right)+O \left( \frac{n^2}{2} \right)+O \left( \frac{n^2}{2} \right) = {\reds{O \left( \frac{n^3}{6} \right)}}
$$
operazioni moltiplicative.
\end{nota}
\end{frame}



\begin{frame}[fragile]
\frametitle{Alcune note sulla eliminazione gaussiana e la fattorizzazione LU}

\smnsp
{\fns{
\begin{esempio}
La matrice simmetrica
$$ 
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 definita positiva, in quanto l'autovalore minimo \mye $\lambda_{\mbox{min}} = 0.271\ldots >0$.

Di conseguenza possiede una decomposizione di Cholesky $A=L L^T$. In questo caso particolare risulta:
$$ 
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].
$$

In Matlab la routine {\tt{chol}} calcola la trasposta di tale matrice $L$. Vediamo in effetti che i risultati sono corretti.
\begin{verbatim}
>> A=gallery('minij',5);
>> LT=chol(A); L=LT';
>> L
L =
     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
>> norm(A-L*L',inf)
ans =
     0
\end{verbatim}
\end{esempio}

}}
\end{frame}



\begin{frame}[fragile]
\frametitle{Alcune applicazioni della fattorizzazione LU: calcolo del determinante di una matrice}

\smnsp

\begin{problema}
Di seguito mostriamo alcune applicazione della fattorizzazione LU, ovvero \begin{itemize}
\item il {\reds{calcolo del determinante di una matrice}},
\item la determinazione dell'inversa di una matrice quadrata non singolare,
\end{itemize}
paragonandoli con altri metodi.

\smsp
A tal proposito ricordiamo che
\bit
\item Il metodo di eliminazione di Gauss pu\myo essere usato per calcolare il determinante di una matrice quadrata.

\item Il costo computazionale della {\reds{formula ricorsiva di Laplace}} per il determinante \mye di circa {\reds{$2n!$ flops}} (cf. {\cite{WIKILAP}}).
\eit
\end{problema}
\end{frame}



\begin{frame}[fragile]
\frametitle{Alcune applicazioni della fattorizzazione LU: calcolo del determinante di una matrice}

\smnsp
\begin{vuoto}
Osserviamo che 
\begin{itemize}
\item se $PA=LU$, visto che $P^{-1}=P^T$ abbiamo che $A=P^{T}LU$ e quindi per il teorema di Binet per cui $\det(M*N)=\det(M)\cdot \det(N)$ qualora $M,N \in {\mathbb{R}}^{n \times n}$
$$
\det(A)=\det(P^{T}LU)=\det(P^{T}) \cdot \det(L) \cdot \det(U).
$$
\item se $s$  \'e il numero di scambi effettuati dal pivoting, allora {\purples{$\det(P^{T})=(-1)^s$}};
\item visto che il determinante di una matrice triangolare $T=(t_{i,j}) \in \C^{n \times n}$ \mye $$
\det(T)=\prod_{k=1}^n t_{k,k}
$$
abbiamo che 
\begin{itemize}
\item da $l_{k,k}=1$, $k=1,\ldots,n$, si ricava {\purples{$\det(L)=1$}}, 
\item {\purples{$\det(U)=\prod_{k=1}^n u_{k,k}$}}, 
\end{itemize}
e quindi
{\reds{
$$
\det(A)=(-1)^s \cdot \prod_{k=1}^n u_{k,k}.
$$
}}
\end{itemize}
\end{vuoto}
\end{frame}



\begin{frame}[fragile]
\frametitle{Alcune applicazioni della fattorizzazione LU: calcolo del determinante di una matrice}

\smnsp
\begin{vuoto}
In definitiva, visto che 
\bit 
\item il calcolo di $\prod_{k=1}^n u_{k,k}$ necessita solo di $n-1$ operazioni moltiplicative, 
\item la fattorizzazione $PA=LU$ ha complessit\mya dell'ordine di $n^3/3$, 
\eit
il calcolo del determinante pu\myo essere effettuato con {\reds{circa $n^3/3$ operazioni moltiplicative}}.

\smsp

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}

\end{vuoto}

\end{frame}



\begin{frame}[fragile]
\frametitle{Alcune applicazioni della fattorizzazione LU: calcolo dell'inversa di una matrice}

\smnsp

\begin{vuoto}
Per calcolare l'inversa di una matrice $A \in \C^{n \times n}$ con il {\reds{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$. 

\smsp

Se brutalmente per ognuna di queste $n^2+1$ matrici di dimensione $(n-1) \times (n-1)$ si usa la fattorizzazione LU, avente complessit\'a $\frac{(n-1)^3}{3}$, necessitano 
\bc
\btc[width=8cm]
$$
C_{COF} = (n^2+1) \cdot \frac{(n-1)^3}{3} \approx {{\frac{n^5}{3}}}
$$
\etc
\ec
operazioni moltiplicative.
\end{vuoto}
\end{frame}



\begin{frame}[fragile]
\frametitle{Alcune applicazioni della fattorizzazione LU: calcolo dell'inversa di una matrice}

\smnsp

{\fns{
\begin{vuoto}
Si supponga per semplicit\mya sia $A=LU$ e si risolvano i sistemi {\reds{$$
Ax^{(j)}=e^{(j)}, \,\,\, j=1,\ldots,n, \quad e^{(j)}=(0,\ldots,0,1,0\ldots,0).
$$}}


Allora {\purples{$A^{-1}$ \mye la matrice la cui $j$-sima colonna corrisponde a $x^{(j)}$}}.

\thsp

Per eseguire i calcoli necessitano
\begin{itemize}
\item una {\reds{fattorizzazione $LU$}} per cui servono circa {\reds{$n^3/3$}} operazioni,
\item risolvere $n$ sistemi lineari $Ax^{(j)}=e^{(j)}$, ognuno richiedente una sostituzione in avanti $Ly^{(j)}=e^{(j)}$ e una all'indietro $Ux^{(j)}=y^{(j)}$. 

\thsp

(a) Nel caso di {\reds{$Ly^{(j)}=e^{(j)}$}}, si vede facilmente che le prime $j-1$ componenti di $y^{(j)}$ sono nulle (si veda ad esempio (\ref{fs}) ponendo $b_1,\ldots,b_{j-1}=0$) e quindi per risolvere $Ly^{(j)}=e^{(j)}$ servono circa {\purples{$ (n-j+1)^2/2$}}
operazioni moltiplicative.


Quindi la risoluzione di tutti i sistemi $Ly^{(j)}=e^{(j)}$, $j=1,\ldots,n$ ha complessit\mya
$$
\sum_{j=1}^n (n-j+1)^2/2=\sum_{j=1}^n j^2/2=(1/2) \sum_{j=1}^n j^2 {\reds{\approx \frac{n^3}{6}}}
$$

\thnsp

(b) Nel caso di {\reds{$Ux^{(j)}=y^{(j)}$}}, il costo computazionale \mye invece di circa
$n^2/2$ operazioni moltiplicative per ogni $j=1,\ldots,n$, ovvero in totale {\reds{$n^3/2$}} operazioni.
\end{itemize}

Quindi, per calcolare l'inversa di $A$ con fattorizzazione $LU$ servono circa $$\frac{n^3}{3}+ \frac{n^3}{6}+ \frac{n^3}{2}= \frac{(2+1+3)n^3}{6}={\reds{n^3}}$$  {\reds{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$.
\end{vuoto}
}}
\end{frame}



\begin{frame}[fragile]
\frametitle{Alcune applicazioni della fattorizzazione LU: calcolo dell'inversa di una matrice}

\smnsp
\begin{esempio}
Si calcoli l'inversa della matrice
$$
A =\left[ \begin{array}{c c c }
     2  &   5   &  7 \\
     1   &  5   &  2 \\
     1    & 4   &  8 \\
     \end{array} \right].
$$
\end{esempio}

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
$$
x_2 =\left[ \begin{array}{c}
    -4.444444444444446e-01 \\
     3.333333333333334e-01\\
    -1.111111111111111e-01\\
     \end{array} \right].
$$
\end{itemize}
\end{frame}



\begin{frame}[fragile]
\frametitle{Alcune applicazioni della fattorizzazione LU: calcolo dell'inversa di una matrice}
\begin{itemize}
\item la soluzione $x_3$ di $Ax=[0,0,1]^T$ \mye
$$
x_3 =\left[ \begin{array}{c}
    -9.259259259259259e-01 \\
     1.111111111111111e-01\\
     1.851851851851852e-01\\
     \end{array} \right].
$$
\end{itemize}
e quindi essendo $A^{-1}=[x_1, \,\,\, x_2, \,\,\, x_3 ]$
$$
\scriptsize{
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{esempio}

In Matlab
\begin{Verbatim}[fontsize=\footnotesize]
>> A=[2 5 7; 1 5 2; 1 4 8]; B=inv(A); format long e
>> B
B =
     1.185185185185185e+00    -4.444444444444445e-01    -9.259259259259258e-01
    -2.222222222222222e-01     3.333333333333334e-01     1.111111111111111e-01
    -3.703703703703703e-02    -1.111111111111111e-01     1.851851851851852e-01
>> 
\end{Verbatim}


\end{frame}




\begin{frame}[fragile]
\frametitle{Metodi iterativi}

\smnsp

{{
\begin{commento}
 I metodi come l'eliminazione gaussiana calcolano in aritmetica esatta la soluzione $x^*$ di $Ax=b$, $A \in \C^{n \times n}$ in un numero finito di passi, sono detti metodi {\reds{diretti}}.

\smsp

D'altra parte
\begin{itemize}
\item la matrice {\reds{$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 {\reds{molti zeri}}, che si verifica spesso non aiutare a migliorare la performance del metodo diretto.
 \end{itemize}

\smsp

Una alternativa consiste nei {\reds{metodi iterativi}} che tipicamente non \mye detto determinino la soluzione esatta, ma possono approssimarla {\reds{arbitrariamente bene}}.
Usualmente {\reds{ogni iterazione}} calcola un prodotto matrice-vettore, e quindi ha complessit\mya di circa {\reds{$n^2$ operazioni moltiplicative}}.
 
\smsp 
 
Se i metodi iterativi forniscono l'approssimazione della soluzione richiesta in un {\reds{$k \ll n/3$}} iterazione, essendo $k n^2 \ll n^3/3$ risultano {\reds{convenienti}} rispetto alla generale fattorizzazione LU.
\end{commento}
}}
\end{frame}










\begin{frame}[fragile]
\frametitle{Metodi iterativi}

{{

Sia $A \in \C^{n \times n}$ e definiamo lo {\it{splitting}} 
\bc
\btc[width=5cm]
$$A=D-E-F$$
\etc
\ec 
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}

\begin{esempio}
Sia $A$ la matrice 
$$
\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 pu\'o essere descritta in termini di $D$, $E$, $F$ con
$$
\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{esempio}
}}
\end{frame}


\begin{frame}[fragile]
\frametitle{Metodi iterativi}

{\small{
\begin{vuoto}
Se $A \in \C^{n \times n}$ \mye una matrice quadrata, risulta importante anche lo {\it{splitting}}

\bc
\btc[width=6cm]
\begin{equation}\label{ispl}
A = P - N, {\mbox{      }}\det(P) \neq 0.
\end{equation}
\etc
\ec

Allora 
$$Ax = b \Leftrightarrow (P-N) x=b \Leftrightarrow Px-Nx=b \Leftrightarrow Px=Nx+b$$ 
da cui, per l'invertibilit\mya di $P$, moltiplicando ambo i membri per $P^{-1}$, ricaviamo
\begin{equation}\label{staz}
x=P^{-1}Px=P^{-1}(Nx+b)={\underbrace{P^{-1}N}_{B}}x+{\underbrace{P^{-1}b}_{c}}=Bx+c:=\phi(x).
\end{equation}
Quindi, posti $B:=P^{-1}N$,  $c:=P^{-1}b$, la soluzione $x^*$ di $Ax=b$ risolve pure $x=Bx+c$ (e viceversa).

\smsp

Conseguentemente, vista la struttura di punto fisso $x=\phi(x)$ con $\phi(x)=Bx+c$, si considera il {\purples{metodo iterativo (stazionario)}}  
\bc
\btc[width=6cm]
\begin{equation}\label{isucc}
x^{(k+1)}=\phi(x^{(k)})=Bx^{(k)}+c
\end{equation}
\etc
\ec
\end{vuoto}
}}
\end{frame}



\begin{frame}[fragile]
\frametitle{Metodi iterativi: alcuni esempi}

{\small{
\begin{vuoto}[Metodi di Jacobi e Gauss-Seidel]
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 {\reds{metodo di Jacobi}} corrisponde a scegliere, se $D$ \mye invertibile, 
\bc
\btc[width=5cm]
$$P=D, {\mbox{   }} N=E+F$$
\etc
\ec
\item il {\reds{metodo di Gauss-Seidel}} corrisponde a scegliere, se $D$ \mye invertibile,
\bc
\btc[width=5cm]
$$P=D-E, {\mbox{   }} N=F$$
\etc
\ec
\end{itemize}
\end{vuoto}
}}

\begin{nota}
{\fns{
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{nota}
\end{frame}






\begin{frame}[fragile]
\frametitle{Metodi iterativi: alcuni esempi}

\smnsp
{\small{
\begin{commento}[Osservazioni sul metodo di Jacobi e Gauss-Seidel]
Se $A=(a_{i,j})_{i,j=1,\ldots,n} \in \C^{n \times n}$, $det(A) \neq 0$  ed $Ax=b$, visto che $(Ax)_i=\sum_{j=1}^n a_{i,j}x_{j}$, allora 
$$\sum_{j=1}^n a_{i,j}x_{j}=b_i, {\mbox{   }} i=1,\ldots,n$$
e quindi, per $i=1,\ldots,n$, evidenziando il termine $i$-simo,
$${\purples{\sum_{j=1}^{i-1} a_{i,j}x_{j}}}+a_{i,i}x_{i}+{\purples{\sum_{j=i+1}^{n} a_{i,j}x_{j}}}=\sum_{i=1}^n a_{i,j}x_{j}=b_i, $$ 
da cui portando a secondo membro $\sum_{j=1}^{i-1} a_{i,j}x_{j}+\sum_{j=i+1}^{n} a_{i,j}x_{j}$
$$a_{i,i}x_{i}=b_i-\sum_{j=1}^{i-1} a_{i,j}x_{j}-\sum_{j=i+1}^{n} a_{i,j}x_{j} ,$$ 
e quindi nel caso $a_{i,i} \neq 0$ per $i=1,\ldots,n$
\bc
\btc[width=10cm]
$${x_{i}=\frac{b_i-\sum_{j=1}^{i-1} a_{i,j}x_{j}-\sum_{j=i+1}^{n} a_{i,j}x_{j}} {a_{i,i}}, {\mbox{   }} i=1,\ldots,n.}$$ 
\etc
\ec
\end{commento}
}}
\end{frame}








\begin{frame}[fragile]
\frametitle{Metodi iterativi: metodo Jacobi}

\begin{vuoto}
Da
{\reds{
$${x_{i}=\frac{b_i-\sum_{j=1}^{i-1} a_{i,j}x_{j}-\sum_{j=i+1}^{n} a_{i,j}x_{j} } {a_{i,i}}, {\mbox{   }} i=1,\ldots,n.}$$ 
}}
il metodo di {\purples{Jacobi}} si ricava ponendo nelle occorrenze del membro di 
\begin{itemize}
\item {\it{sinistra}} {\purples{$x^{(k+1)}_{i}$}} al posto di $x_i$,
\item  {\it{destra}} {\purples{$x^{(k)}_{j}$}} al posto di ogni occorrenza di $x_j$ (con $j=1,\ldots,n$),
\end{itemize}

\noindent ricavando

\bc
\btc
\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}	
\etc
\ec
\end{vuoto}

\begin{nota}
Il metodo ottenuto e' lo stesso descritto matricialmente come $$x^{(k+1)}=P^{-1}Nx^{(k)}+P^{-1}b$$ con $P=D$ e $N=E+F$.
\end{nota}
\end{frame}





\begin{frame}[fragile]
\frametitle{Metodo di Jacobi}

\smnsp

\begin{esempio}
{\small{
Applicando 

\smnsp

$$x^{(k+1)}_i={ \frac {1} {a_{i,i}} } \left(b_i-\left( \sum_{j<i} a_{i,j} x^{(k)}_j+\sum_{j>i} a_{i,j} x^{(k)}_j\right)\right), \,\,i=1,\ldots,n$$ al caso di $A \in {\mathbb{C}}^{3 \times 3}$, otteniamo, supponendo disponibile $x^{(k)}$:

\begin{enumerate}
\item $x^{(k+1)}_1=\frac {1} {a_{1,1}} (b_1- a_{1,2} x^{(k)}_2- a_{1,3} x^{(k)}_3)$;
\item $x^{(k+1)}_2=\frac {1} {a_{2,2}}  (b_2- a_{2,1} {\reds{x^{(k)}_1}}- a_{2,3} x^{(k)}_3)$;
\item $x^{(k+1)}_3=\frac {1} {a_{3,3}}  (b_3- a_{3,1} {\reds{x^{(k)}_1}}- a_{3,2} {\reds{x^{(k)}_2}})$.
\end{enumerate}

Si nota che 
\begin{itemize}
\item dopo aver calcolato $x^{(k+1)}_1$ tale nuova approssimazione della prima componente $x^*_1$ della soluzione $x^*$ non \'e sfruttata nel determinare $x^{(k+1)}_2$;
\item dopo aver calcolato $x^{(k+1)}_2$ tale nuova approssimazione della prima componente $x^*_2$ della soluzione $x^*$ non \'e sfruttata nel determinare $x^{(k+1)}_3$.
\end{itemize}
Ogni $x^{(k+1)}_j$, $j=1,2,3$ \'e calcolato direttamente dall'iterazione precedente $x^{(k)}=(x^{(k)}_i)_i$, senza richiedere la conoscenza di qualche $x^{(k+1)}_i$ con $i=1,\ldots,j-1$.

}}


\end{esempio}
\end{frame}








\begin{frame}[fragile]
\frametitle{Metodi iterativi: metodo di Gauss-Seidel}

\begin{vuoto}
Da
{\reds{
$${x_{i}=\frac{b_i-\sum_{j=1}^{i-1} a_{i,j}x_{j}-\sum_{j=i+1}^{n} a_{i,j}x_{j} } {a_{i,i}}, {\mbox{   }} i=1,\ldots,n.}$$ 
}}
il metodo di {\purples{Gauss-Seidel}} si ricava ponendo nelle occorrenze del membro di 
\begin{itemize}
\item {\it{sinistra}} {\purples{$x^{(k+1)}_{i}$}} al posto di $x_i$,
\item  {\it{destra}} nella prima sommatoria {\purples{$x^{(k+1)}_{i}$}} al posto di $x_i$ e nella seconda sommatoria {\purples{$x^{(k)}_{i}$}}
\end{itemize}

\noindent ricavando

\bc
\btc
\begin{equation}\label{GS}
	x_{i}^{(k+1)}=(b_{i} - \sum_{j=1}^{i-1} a_{ij} x_{j}^{(k+1)}- 
	\sum_{j=i+1}^{n} a_{ij} x_{j}^{(k)} )/a_{ii}, \,\, i=1,\ldots,n,
\end{equation}	
\etc
\ec
\end{vuoto}

\begin{nota}
{\fns{
Il metodo ottenuto e' lo stesso descritto matricialmente come $$x^{(k+1)}=P^{-1}Nx^{(k)}+P^{-1}b, \smhsp P=D-E, \smhsp N=F.$$ Una sua generalizzazione porta ai metodi SOR (cf. \cite{WIKISOR}).
}}
\end{nota}
\end{frame}











\begin{frame}[fragile]
\frametitle{Metodo di Gauss-Seidel}

\smnsp
\smnsp

\begin{esempio}
{{
Applicando $$x^{(k+1)}_i={ \frac {1} {a_{i,i}} } \left(b_i-\left( \sum_{j<i} a_{i,j} x^{(k+1)}_j+\sum_{j>i} a_{i,j} x^{(k)}_j\right)\right), \,\,i=1,\ldots,n$$ al caso di $A \in {\mathbb{C}}^{3 \times 3}$, otteniamo:

\begin{enumerate}
\item $x^{(k+1)}_1=\frac {1} {a_{1,1}} (b_1- a_{1,2} x^{(k)}_2- a_{1,3} x^{(k)}_3)$;
\item $x^{(k+1)}_2=\frac {1} {a_{2,2}}  (b_2- a_{2,1} {\reds{x^{(k+1)}_1}}- a_{2,3} x^{(k)}_3)$;
\item $x^{(k+1)}_3=\frac {1} {a_{3,3}}  (b_3- a_{3,1} {\reds{x^{(k+1)}_1}}- a_{3,2} {\reds{x^{(k+1)}_2}})$.
\end{enumerate}

Si nota che 
\begin{itemize}
\item dopo aver calcolato $x^{(k+1)}_1$ tale nuova approssimazione della prima componente $x^*_1$ della soluzione $x^*$ \'e sfruttata nel determinare $x^{(k+1)}_2$;
\item dopo aver calcolato $x^{(k+1)}_2$ tale nuova approssimazione della prima componente $x^*_2$ della soluzione $x^*$ \'e sfruttata nel determinare $x^{(k+1)}_3$.
\end{itemize}

}}

\smnsp

\end{esempio}

\smnsp


\end{frame}









\begin{frame}[fragile]
\frametitle{Metodi iterativi, convergenza}

\smnsp

{\small{
\begin{nota}[Consistenza]
Consideriamo un metodo iterativo stazionario

$$
x^{(k+1)} = B x^{(k)}+ c, \,\,\, k = 0, 1, 2, \ldots
$$
per l'approssimazione della soluzione $x^*$ del sistema lineare $Ax=b$, con $A $ invertibile, ovvero $\det(A) \neq 0$.

\begin{enumerate}
\item Vogliamo che il metodo iterativo stazionario converga a $x^*$, ovvero che {\reds{$$\lim_{k \rightarrow +\infty}\|x^{(k)} - x^*\|=0$$}} per qualche norma indotta, e che tale $x^*$ risolva $Ax=b$. 
\item Osserviamo che se $x^{(k)} \rightarrow v^*$ allora pure $x^{(k+1)} \rightarrow v^*$ da cui,
\begin{eqnarray}
{\reds{v^*}} &=& \lim_{k \rightarrow +\infty} x^{(k+1)} = \lim_{k \rightarrow +\infty} (Bx^{(k)} +c)= \lim_{k \rightarrow +\infty} Bx^{(k)} + \lim_{k \rightarrow +\infty} c \nonumber \\
&=& B \lim_{k \rightarrow +\infty} x^{(k)} +c = {\reds{Bv^*+c}}, \nonumber
\end{eqnarray}
ovvero {\reds{$v^*$ risolve il problema di punto fisso $x=\phi(x)$}} con 

$$
\phi(x)=Bx+c.
$$

Siccome desideriamo che $x^{(k)} \rightarrow x^*$, con $x^*$ soluzione di $Ax=b$, necessariamente si chiede che tale $x^*$ sia l'unica soluzione del problema di punto fisso $x=Bx+c$. Qualora $\|B\| < 1$, tale propriet\'a detta di consistenza \'e garantita tanto dal metodo di Jacobi che da quello di Gauss-Seidel.
\end{enumerate}
\end{nota}
}}
\end{frame}





\begin{frame}[fragile]
\frametitle{Metodi iterativi, convergenza}

\smnsp

\begin{teorema}[Convergenza globale dei metodi iterativi]
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$ sia l'unica soluzione di $x = Bx + c$.
\end{itemize}
Allora {\reds{qualsiasi sia}} $x_0 \in \C^n$, la successione definita da
$
x^{(k+1)} = B x^{(k)}+ c$, con $k = 0, 1, 2, \ldots$, converge a $x^*$.
\end{teorema}

{\tiny{
\begin{dimostrazione}[Facoltativa]
Da $\|B\| < 1$, abbiamo che la matrice $I-B$ \mye invertibile e quindi $x^*$ \mye l'unica soluzione di 
$(I-B)x=c$ ovvero di $x=Bx+c$. Da 
$$
x^{(k+1)} = B x^{(k)} + c, \,\,\, k = 0, 1, 2, \ldots, {\mbox{    }}
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, da $\|My\| \leq \|M\| \|y\|$, $\|M N\| \leq \|M\| \|N\|$, per $M$, $N$, $y$ arbitrari, 
\begin{eqnarray}
0 \leq \|x^{(k)}-x^*\|&=&\|B(x^{(k-1)}-x^*)\| \leq \|B\| \|x^{(k-1)}-x^*\| \leq \|B\|^2 \|x^{(k-2)}-x^*\|  \nonumber \\
&\leq&  |B\|^3 \|x^{(k-3)}-x^*\| \leq  \ldots \leq  \|B\|^{k+1} \|x^{(0)}-x^*\| \nonumber
\end{eqnarray}
e per il teorema del confronto, visto che $\|B\|^k \rightarrow 0$, 
{\reds{$
\lim_k \|x^{(k)}-x^*\|=0
$}}
ovvero $x^{(k)} \rightarrow x^*$.
\end{dimostrazione}
}}
\end{frame}



\begin{frame}[fragile]
\frametitle{Metodi iterativi, convergenza}

\smnsp

\begin{teorema}[Convergenza globale dei metodi iterativi]\label{cnorm}
Siano
\begin{itemize}
\item $B \in \C^{n \times n}$,
\item si supponga che la soluzione $x^*$ di $Ax=b$ sia l'unica soluzione di $x = Bx + c$.
\item sia  $\rho(B)=\max_{k=1,\ldots,n}|\lambda_k|$ il raggio spettrale di $B$ (dove $\lambda_1,\ldots,\lambda_n$ sono gli autovalori di $B$).
\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^*$ {\reds{se e solo se}} $\rho(B) < 1$ .
\end{teorema}

\smnsp

\begin{nota}
\begin{enumerate}
\item Si noti che in questo asserto compare, a differenza del precedente {\purples{un se e solo se}}.
\item Si noti compare il {\purples{raggio spettrale di $B$}} (e non di $A$).
\item Si noti che richiede il calcolo dell'autovalore di modulo massimo di $B$, un problema potenzialmente complicato. Di solito, noto la struttura della matrice $A$ ci sono dei teoremi che offrono risultati sul raggio spettrale di $B$.
\end{enumerate}
\end{nota}
\end{frame}







\begin{frame}[fragile]
\frametitle{Metodi iterativi, convergenza caso predominanza diagonale stretta}

A partire dai precedenti teoremi si \mye potuto dimostrare (non facile!)

\begin{teorema}[Convergenza Jacobi per matrici a predominanza diagonale stretta]
Se la matrice $A$ \mye a {\reds{predominanza diagonale (stretta) per righe}}, allora il metodo di {\purples{Jacobi}} converge.
\end{teorema}

\begin{teorema}[Convergenza Gauss-Seidel per matrici a predominanza diagonale stretta]
Se la matrice $A$ \mye a {\reds{predominanza diagonale (stretta) per righe}}, allora il metodo di {\purples{Gauss-Seidel}} converge.
\end{teorema}

\begin{teorema}[Convergenza G.-S. per matrici simmetriche e definite positive]
Se la matrice $A$ \mye {\reds{simmetrica e definita positiva}} allora il metodo di {\purples{Gauss-Seidel}} converge.
\end{teorema}

\begin{nota}[Grandezza raggio spettrale]
In generale, con stime quantitative, si vede che sono da preferire metodi con raggio spettrale $\rho(B)$ il pi\'u piccolo possibile.
\end{nota}

\end{frame}



\begin{frame}[fragile]
\frametitle{Alcuni esempi}

\smnsp
{{
\begin{esempio}[Matrice a predominanza diagonale stretta]
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}}).
Si dica se convergono i metodi di Jacobi e Gauss-Seidel.
\end{esempio}


\begin{itemize}
\item 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$ e questo si dimostra implicare che asintoticamente il metodo di Jacobi converge pi\'u rapidamente di quello di Gauss-Seidel.
\item 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{itemize}
}}
\end{frame}



\begin{frame}[fragile]
\frametitle{Alcuni esempi}

\smnsp
 \begin{table}[!ht]
\begin{center}
\begin{tabular}{| c c c |}
\hline
$n$&   $E_J$ & $E_{GS}$  \\
\hline
 	 $ 10$ & $1.6e-01$ & $2.2e-01$ \\
 	 $ 30$ & $1.5e-03$ & $6.2e-03$ \\
 	 $ 50$ & $1.4e-05$ & $1.7e-04$ \\
 	 $ 70$ & $1.3e-07$ & $4.8e-06$ \\
 	 $ 90$ & $1.2e-09$ & $1.3e-07$ \\
 	 $110$ & $1.1e-11$ & $3.8e-09$ \\
 	 $130$ & $1.0e-13$ & $1.1e-10$ \\
 	 $150$ & $1.1e-15$ & $2.9e-12$ \\
 	 $170$ & $-$ & $8.2e-14$ \\
 	 $190$ & $-$ & $2.2e-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}

\end{frame}






\begin{frame}[fragile]
\frametitle{Alcuni esempi}

\smnsp

\begin{esempio}[Matrice definita positiva]
Si studi la convergenza dei metodi di Jacobi e Gauss-Seidel, relativamente al 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$ \\ 
 	 \hline 
 	 $-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$ \\ 
 	 \hline 
 	 $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$. 
\end{esempio}

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

Le linee mostrano come $A \in \R^{9 \times 9}$ segua {\it{pattern}} particolari con {\sl{blocchi}} pari a 
$$
\footnotesize{
B = \left[
  \begin{array}{c c c}
 	 $4$ & $-1$ & $0$  \\ 
 	 $-1$ & $4$ & $-1$  \\ 
 	 $0$ & $-1$ & $4$ 
         \end{array}
\right],
\,\,\,
-I = \left[
  \begin{array}{c c c}
 	 $-1$ & $0$ & $0$  \\ 
 	 $0$ & $-1$ & $0$  \\ 
 	 $0$ & $0$ & $-1$ 
         \end{array}
\right],
\,\,\,
0_{3 \times 3} = \left[
  \begin{array}{c c c}
 	 $0$ & $0$ & $0$  \\ 
 	 $0$ & $0$ & $0$  \\ 
 	 $0$ & $0$ & $0$ 
         \end{array}
\right],
}
$$

\end{frame}



\begin{frame}[fragile]
\frametitle{Alcuni esempi}

\smnsp

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, in quanto $A$ \'e simmetrica e definita positiva.  

\smsp

Si vede che $\rho(B_J) \approx 0.7071$, mentre $\rho(B_{GS}) \approx 0.5$ e quindi ci si aspetta che asintoticamente il metodo di Gauss-Seidel converga in meno iterazioni di quello di Jacobi.

 \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$.  convergenti.}
\end{center}
\end{table}
 	 % % \end{esempio}
\end{frame}



\begin{frame}[fragile]
\frametitle{Alcuni esempi}

\smnsp
 	 
 	 
 	 
\begin{esempio}[Mancata convergenza dei metodi di Jacobi e Gauss-Seidel]
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$. Si verifichi la mancata convergenza globale dei metodi di Jacobi e Gauss-Seidel partendo da $x^{(0)}=(1-\epsilon,1-\epsilon,1-\epsilon)^T$ con $\epsilon=10^{-14}$.
\end{esempio}

\begin{itemize}
\item 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}$.

\item Entrambi i metodi citati, pur partendo da un tal $x^{(0)}$, non convergono, come si pu\myo vedere dalla relativa tabella. In effetti, si verifica che $1 < \rho(B_J) \approx 8.66$, $1 < \rho(B_{GS}) \approx 9.62$, e quindi entrambi i metodi iterativi non sono globalmente convergenti.
\end{itemize}
\end{frame}



\begin{frame}[fragile]
\frametitle{Alcuni esempi}

\smnsp
 \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{esempio}
\end{frame}




\begin{frame}[fragile]
\frametitle{Test di arresto}

\smnsp

\begin{vuoto}
Sia $\epsilon >0$ una tolleranza fissata dall'utente.
Il metodo iterativo
$$
x^{(k+1)}=Bx^{(k)}+c
$$
viene arrestato utilizzando 
\begin{itemize}
\item il {\reds{criterio dello step}}, ovvero si interrompe il processo quando
\bc
\btc[width=5cm]
$$
\|x^{(k+1)} - x^{(k)}\| \leq \epsilon,
$$
\etc
\ec
\item {\reds{criterio del residuo}}, ovvero si interrompe il processo quando
\bc
\btc[width=5cm]
$$
\|b-Ax^{(k)}\| \leq \epsilon,
$$
\etc
\ec
\item {\reds{criterio del residuo relativo}}, ovvero si interrompe il processo quando
\bc
\btc[width=5cm]
$$
{ \frac {\|b-Ax^{(k)}\|} {\|b\|} } \leq \epsilon, \smhsp b \neq 0_{{\mathbb{C}}^n}
$$
\etc
\ec
\end{itemize}
\end{vuoto}
\end{frame}
















\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati}

Si consideri il sistema lineare di 3 equazioni in due incognite
$$\left\{
\begin{array}{l} 
x_1+x_2=1 \\ 
x_1-x_2=0 \\ 
x_1+3x_2=0
\end{array}
\right.$$

Risulta chiaro che il sistema non ha soluzione. Infatti l'unica soluzione 
delle prime due equazioni \mye $x_1=x_2=1/2$, che per\myo non verifica 
$x_1+3x_2=0$.

\bigskip

In termini matriciali, il problema $Ax=b$, non ha soluzione, dove 
$$A=\left( \begin{array}{ll} 1 & 1 \\ 1  & -1 \\ 1 & 3 \end{array} 
\right), {\mbox{    }} b=(1,0,0)^T.$$

In alternativa, risulta quindi ragionevole calcolare $x^*=(x_1^*,x_2^*)$ 
tale che 
\begin{equation}
\gamma:=\min_{x \in {\mathbb{R}}^n }  \|b-Ax\|_2=\|b-Ax^*\|_2
\end{equation}
dove
$$\|u\|_2=\sqrt{\sum_i u_i^2}.$$
\end{frame}
















\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati}
{\small{
La prima questione riguarda l'esistenza e unicit\mya di tale minimo $x^*$. 
\begin{teorema}[cf. {\cite{BINI}}, p. 432]\label{th1} 
Sia $X$ l'insieme dei vettori di ${\mathbb{R}}^n$ tali che ${x}^* \in 
X$ se e solo se
\begin{equation}
\min_{x \in {\mathbb{R}}^n }  \|b-Ax\|_2=\|b-A{x}^* \|_2
\end{equation}
Supponiamo $A \in {\mathbb{R}}^{{m}\times {n}} $ con $m \geq n$ e $b \in 
{\mathbb{R}}^m$.
Valgono le seguenti propriet\mya
\begin{enumerate}
\item $x \in X$ se e solo se 
\begin{equation}\label{eq2}
A^TAx=A^Tb,
\end{equation}
cio\'e $x$ risolve il {\reds{sistema delle equazioni normali}}.
\item $X$ \mye un insieme non vuoto, chiuso e convesso (ovvero se $u,v \in X$ allora $\theta u +(1-\theta) v \in X$, $\theta \in [0,1]$ .
\item Esiste $x^* \in X$ tale che
$$\|x^*\|_2=\min_{x\in X} \|x\|_2.$$
Tale $ x^*$ \mye detto soluzione di {\reds{minima norma}}.
\item L'insieme $X$ si riduce ad {\reds{un solo elemento}} $x^*$ se e solo se la 
matrice $A$ ha {\blues{\href{https://it.wikipedia.org/wiki/Rango\_(algebra\_lineare)}{rango massimo}}}.
\end{enumerate}
\end{teorema}
Si osservi che tale teorema vale pure sostituendo $\mathbb{C}$ a ogni occorrenza di $\mathbb{R}$.
}}
\end{frame}













\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati}

\smnsp

Di conseguenza se $A \in {\reds{\R^{m \times n}}}$, con $m \geq n$ \mye una matrice di {\blues{\href{https://it.wikipedia.org/wiki/Rango\_(algebra\_lineare)}{rango massimo}}} $n$ (ovvero, le colonne di A
sono $n$ vettori linearmente indipendenti di $\R^m$ e quindi si pu\'o {\sl{estrarre}}, eliminando righe, una matrice $n \times n$ non singolare), allora la soluzione del sistema sovradeterminato $Ax=b$ 
ai minimi quadrati,
\begin{equation}\label{(20)}
\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 lineare  detto {\reds{delle equazioni normali}}
$$
A^T A x=A^Tb
$$
dove la matrice $A^T A$ \mye simmetrica e definita positiva.


\end{frame}









\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati e un problema di approssimazione polinomiale}

\smnsp
{\footnotesize{
Desideriamo calcolare il polinomio di miglior approssimazione $p^*_n \in {\mathbb{P}}_n$ tale che date $m$ coppie ($s_i, t_i$), $i=1,\ldots,m$, $m \geq n+1$, sia minima 
\begin{equation}\label{eq22}
\min_{p_n \in  {\mathbb{P}}_n} \sqrt{\sum_{i=1}^m (t_i-p_n(s_i))^2}.
\end{equation}

Visto che {\purples{$p_n(s)=\sum_{j=1}^{n+1} {\gamma}_j s^{j-1}=s_1+s_2 s+s_3 s^2+\ldots+s_{n+1} s^n$}}, da (\ref{eq22}) abbiamo che i coefficienti incogniti $\gamma=(\gamma_j)_{j=1,\ldots,n+1}$ sono tali che sono soluzione di 
\begin{equation}\label{eq23}
\min_{\gamma \in  {\mathbb{R}}^{n+1}} \sqrt{\sum_{i=1}^m (t_i-\sum_{j=1}^{n+1} \gamma_j s_i^{j-1})^2}
\end{equation}

ovvero $\gamma$ \mye soluzione ai minimi quadrati del sistema sovradeterminato $A \gamma=b$, cio\mye
\begin{equation}\label{eq24}
\min_{\gamma \in \R^{n+1}} \|b-A\gamma\|_2= \min_{\gamma \in \R^{n+1}} \sqrt{\sum_{i=1}^m (b_i-\sum_{j=1}^n a_{i,j} \gamma_j)^2}
\end{equation}

in cui 
\begin{itemize}
\item $A=(a_{i,j}) \in \R^{m \times n}$ con {\purples{$a_{i,j}=s_i^{j-1}$}}, 
\item $b=(t_i)_{i=1,\ldots,m}$, 
\item $\gamma=(\gamma_i)_{i=1,\ldots,n+1}$.
\end{itemize}
Quindi per risolvere (\ref{eq22}) si determina la soluzione $\gamma=(\gamma_i)$ del sistema sovradeterminato (\ref{eq24}) e con tali coefficienti si ottiene il polinomio di miglior approssimazione $p_n(s)=\sum_{j=1}^{n+1} {\gamma}_j s^{j-1}$.
}}
\end{frame}






\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati e fattorizzazione di Cholesky}
Se la matrice $A \in \R^{m \times n}$, $m \geq n$, ha rango massimo (ossia pari a $n$), si dimostra che la matrice $A^T A$ \mye simmetrica e definita positiva e quindi possiede una fattorizzazione di Cholesky, ossia esiste una matrice $L$ triangolare inferiore con elementi diagonali positivi tale che {\reds{$A^TA=LL^T$}}.

\smsp

Di conseguenza, la soluzione $x^*$ del problema ai minimi quadrati deve  risolvere le equazioni normali {\purples{$A^TAx=A^Tb$}} e posto $y=L^Tx$ \'e tale che
$$
A^TAx=A^Tb \Leftrightarrow LL^Tx=A^Tb  \Leftrightarrow Ly=A^Tb
$$
Quindi si risolve prima {\purples{$Ly=A^Tb$}} e, detta $y^*$ la soluzione, di seguito {\purples{$L^Tx=y^*$}}.

Il costo computazionale per la risoluzione del sistema \mye
\begin{itemize}
\item $n^2m/2$ operazioni moltiplicative per il calcolo di $A^T A$ (si usa il fatto che tale matrice \mye simmetrica!);
\item $n m$ operazioni moltiplicative per il calcolo di $A^T b$;
\item $n^3/6$ operazioni per la risoluzione del sistema lineare $A^T Ax=A^Tb$ col metodo di Cholesky.
\end{itemize}
e visto che usualmente $\frac{n^2m}{2} \gg n m$, il numero di operazioni moltiplicative risulta
$$
\frac{n^2m}{2} +n m +\frac{n^3}{6} \approx {\reds{\frac{n^2m}{2} + \frac{n^3}{6}}}.
$$

\end{frame}






\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati e fattorizzazione di Cholesky}



{\small{

\begin{esempio}
Risolvere il problema 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].
$$
mediante fattorizzazione di Cholesky.
\end{esempio}

Svolgiamo i conti in Matlab:
\begin{lstlisting}[basicstyle=\small]
>> A=[2 5; 1 1 ; 3 2];
>> b=ones(3,1);
>> R=chol(A'*A); % triangolare sup. tale che R'*R=A
>> L=R'; % triangolare inf. tale che R'*R=A
>> y=L\A'*b; % risolve Ly=A^Tb
>> x=L'\y % risolve L'x=y.
x =
     3.358778625954199e-01
     7.633587786259541e-02
>> 
\end{lstlisting}

}}
\end{frame}









\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati e fattorizzazione QR}

\smnsp
{\small{
\begin{definizione}[{\href{https://it.wikipedia.org/wiki/Decomposizione\_QR}{Fattorizzazione QR}}]
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$.
\end{definizione}

\begin{importante}
Sotto queste ipotesi, si mostra che la fattorizzazione QR esiste ed e' unica.
\end{importante}

\begin{esempio}
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{esempio}
}}
\end{frame}











\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati e fattorizzazione QR}

\smnsp
\begin{esempio}
Si determini, se esistente, la fattorizzazione $QR$ di
$$
A=\left[
\begin{array}{c c}
1 & 0 \\ 1 & 1 \\ 0 & 1
\end{array}
\right]
$$
\end{esempio}

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$).
\end{frame}







\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati e fattorizzazione QR}

\smnsp
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_3\|_{\infty} \approx 6.3642e-16$, 
\item $ \|QQ'-I_3\|_{\infty} \approx 7.3882e-16$, 
\item $\|A-QR\|_{\infty} \approx 7.3882e-16$.
\end{itemize}
dove $I_3 \in \R^{3 \times 3}$ \'e la matrice identica. 

\smsp

Il fatto che le norme infinito matriciali siano vicino a $0$ implica che il loro argomento e' circa la matrice nulla, ovvero che numericamente 
\begin{itemize}
\item $Q'Q-I_3 \approx 0$ (cio\'e $Q'Q \approx I_3$), 
\item $QQ'-I_3 \approx 0$ (cio\'e $QQ' \approx I_3$), 
\item $A-QR \approx 0$ (cio\'e $A \approx QR$).
\end{itemize}
% % \end{esempio}
\end{frame}




\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati e fattorizzazione QR}

\smnsp

Vediamo un ulteriore esempio in Matlab.

\begin{lstlisting}[basicstyle=\footnotesize]
>> A=[1 2 3; 5 3 6; 2 8 3; 1 5 9] % matrice 4 x 3 (m=4, n=3)
A =
     1     2     3
     5     3     6
     2     8     3
     1     5     9
>> [Q,R]=qr(A) % display in "format short"
Q = % matrice tale che Q'*Q=I
   -0.1796    0.1040    0.1548   -0.9659
   -0.8980   -0.4203    0.0297    0.1265
   -0.3592    0.7453   -0.5587    0.0575
   -0.1796    0.5070    0.8142    0.2185
R = % matrice rettangolare e "triangolare"
   -5.5678   -6.8250   -8.6211
         0    7.4444    4.5888
         0         0    6.2944
         0         0         0
>> norm(A-Q*R) % numericamente A=Q*R se  norm(A-Q*R) e' piccola
ans =
   4.9370e-15
>> Q'*Q % verifichiamo che Q'*Q=I
ans =
    1.0000   -0.0000   -0.0000    0.0000
   -0.0000    1.0000         0   -0.0000
   -0.0000         0    1.0000   -0.0000
    0.0000   -0.0000   -0.0000    1.0000
>> 
\end{lstlisting}


\end{frame}





\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati e fattorizzazione QR}

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

 \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{nota}
\end{frame}












\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati e fattorizzazione QR}

\smnsp

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

Supponiamo di dover calcolare la soluzione $x$ ai minimi quadrati di $Ax=b$. 

\smsp

Per il teorema precedente $x^*$ \mye pure soluzione delle equazioni normali {\purples{$A^TA x =A^Tb$}}.

Essendo 
\begin{enumerate}
\item $A=QR$, 
\item $Q^T Q =I$, 
\item $(QR)^T=R^T Q^T$, 
\end{enumerate} 
$$A^TA{\stackrel{(1)}{=}}(QR)^T QR{\stackrel{(3)}{=}}R^T Q^T Q R{\stackrel{(2)}{=}}R^T R$$
e quindi {\purples{$A^TA x =A^Tb$}} se e solo se
$$
R^T R x = A^T b.
$$
Ricordando che $R \in \R^{m \times n}$ \mye una matrice triangolare, {\reds{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} 
}}
\end{frame}






\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati e fattorizzazione QR}

{\magentas{\underline{Stadio 1}}}.  Analizziamo il {\purples{primo sistema $R^Ty=A^Tb$}} con $R \in \R^{m \times n}$, $A^Tb \in \R^{n \times 1}$. 


Se per una certa matrice $L \in \R^{n \times n}$ triangolare inferiore
$$
R^T= [{\underbrace{L}_n} {\mbox{ }} {\underbrace{\bf{0}}_{m-n}} ] \} n, {\mbox{ }} y=\left[
\begin{array}{l}
y_1\\
y_2
\end{array}
\right],
y_1 \in \R^{n}, y_2 \in \R^{m-n}
$$
si vede facilmente che essendo $$A^Tb=R^Ty= [L  \quad {\bf{0}}] \left[
\begin{array}{l}
y_1\\
y_2
\end{array}
\right]=  Ly_1+{{\bf{0}}}y_2=Ly_1$$ si ha {\reds{$$R^Ty=A^Tb \Leftrightarrow Ly_1=A^Tb$$}} 
\noindent che richiede la risoluzione di un sistema triangolare inferiore $Ly_1=A^Tb$ di dimensione $n \times n$.

\smsp

Si osservi che qualsiasi sia $y^*_2 \in \R^{m-n}$, il vettore $y=[y^*_1, y^*_2]^T$ risolve $R^Ty=A^Tb$.
\end{frame}









\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati e fattorizzazione QR}

\smnsp
{\magentas{\underline{Stadio 2}}}.  Una volta ottenuto $y^*_1 \in \R^{n \times 1}$ soluzione di $Ly_1=A^Tb$, e posto $y^*_2 \in \R^{(m-n) \times 1}$ arbitrario, visto che $R \in \R^{m \times n}$, dobbiamo risolvere il sistema sovradeterminato
$Rx=y^*$ con $y^*=[y^*_1; y^*_2] \in \R^{m \times 1}$. 


A tal proposito, visto che $x^* \in \R^{n \times 1}$ e
$$
R= \left[
\begin{array}{l}
L^T\\
{\bf{0}}
\end{array}
\right]
, {\mbox{ }} 
y^*=\left[
\begin{array}{l}
y^*_1\\
y^*_2
\end{array}
\right] \Rightarrow {\reds{\left[
\begin{array}{l}
y^*_1\\
y^*_2
\end{array}
\right]}}= y^*=Rx=\left[
\begin{array}{l}
L^T\\
{\bf{0}}
\end{array}
\right] x = {\reds{\left[
\begin{array}{l}
L^Tx\\
0_{\R^{m-n}}
\end{array}
\right]}},
$$
e quindi la soluzione $x^* \in \R^{n \times 1}$ del problema sovradeterminato $Ax=b$ pu\'o essere ricavata calcolando la soluzione del sistema quadrato triangolare superiore e nonsingolare $L^T x=y^*_1$.
\end{frame}




















\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati e fattorizzazione QR}

\smnsp
Poich\'e la soluzione ai minimi quadrati richiede la soluzione di
\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} 
e considerato per $A$ generica
\begin{itemize}
\item il costo della {\purples{fattorizzazione QR}} (si utilizza un metodo dovuto ad Householder o a Givens);
\item la {\purples{valutazione di $A^T b$}} che richiede in generale circa {\reds{$m n$}} operazioni moltiplicative, visto che $A \in \R^{m \times n}$;
\item la {\purples{risoluzione di $R^T y=A^T b$}} che richiede di determinare le prime $n$ componenti di $y$ (le altre vengono poste uguali a 0) e quindi richiede circa {\reds{$n^2/2$}} operazioni, visto che $R$ \mye triangolare;
\item la {\purples{risoluzione di $Rx=y^*$}} che richiede circa {\reds{$n^2/2$}} operazioni, visto che $R \in \R^{m \times n}$ \mye una matrice triangolare e le ultime $m-n$ componenti di $y$ sono nulle.
\end{itemize} 
Si dimostra che il numero di moltiplicazioni necessarie per la risoluzione di sistemi sovradeterminati, con il metodo QR, risulta circa
{\purples{$2mn^2-2n^3/3$}}
utilizzando il metodo di Householder e {\purples{$3mn^2-n^3$}}
utilizzando il metodo di Givens (cf. \cite[p.263]{GVL}).
\end{frame}




\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati e fattorizzazione QR}

\smnsp
\begin{esempio}
Risolvere, mediante fattorizzazione QR, il problema 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].
$$
\end{esempio}

\smnsp

\begin{lstlisting}
>> format long; A=[2 5; 1 1 ; 3 2]; b=ones(3,1); [Q,R]=qr(A);
>> Q
Q =
  -0.534522483824849   0.840626032489748   0.087370405666104
  -0.267261241912424  -0.070052169374146  -0.961074462327142
  -0.801783725737273  -0.537066631868450   0.262111216998311
>> R
R =
  -3.741657386773941  -4.543441112511214
                   0   3.058944729337695
                   0                   0
>> Rtrasposta=R';
>> L=Rtrasposta(1:2,1:2)
L =
  -3.741657386773941                   0
  -4.543441112511214   3.058944729337695
>> y=L\(A'*b)
y =
  -1.603567451474546
   0.233507231247152
>> x=L'\y
x =
   0.335877862595420
   0.076335877862595
>> 
\end{lstlisting}
\end{frame}




\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati e fattorizzazione SVD}

\begin{definizione}[Decomposizione SVD (ovvero ai valori singolari, cf. {\cite[p.444]{BINI}})]
Sia $A \in {\mathbb{R}}^{m \times n}$. Allora esistono 
\begin{itemize}
\item una matrice $U \in {\mathbb{R}}^{m \times m}$ unitaria (ovvero $U^T U=U U^T=I$);
\item una matrice $V \in {\mathbb{R}}^{n \times n}$ unitaria (ovvero $V^T V=V V^T=I$);
\item una matrice $\Sigma \in {\mathbb{R}}^{m \times n}$ che ha elementi $\sigma_{i,j}=0$ se $i \neq j$ e se $i=j$ allora $\sigma_{i,i}=\sigma_i$ con $$\sigma_1 \geq \sigma_2 \geq \ldots \geq \sigma_k \geq \sigma_{k+1}= \ldots =\sigma_p = 0, {\mbox{     }} p=\min\{m,n\}$$
tali che
$$
A=U \Sigma V^T.
$$
Gli elementi $\sigma_j$ sono detti {\purples{valori singolari}} di $A$.
\end{itemize}
\end{definizione}

\begin{itemize}
\item I valori singolari $\{\sigma_k\}_{k=1,\ldots,p}$ sono le radici quadrate degli autovalori di $A^TA$.
\item Se la matrice $A$ \'e simmetrica allora i valori singolari corrispondono con gli autovalori di $A$ in modulo.
\item L'ultimo $k$ per cui $\sigma_k > 0$ \'e il rango di $A$ (cf. \cite[p.448]{BINI}).
\end{itemize}

\end{frame}




\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati e fattorizzazione SVD}

\begin{definizione}[Decomposizione SVD (ovvero ai valori singolari, cf. {\cite[p.454]{BINI}})]
Sia $A \in {\mathbb{R}}^{m \times n}$ con $m \geq n \geq k$ di rango $k$ e {\purples{$
A=U \Sigma V^T.
$}}
la decomposizione in valori singolari di $A$. 
Allora, se
\begin{itemize}
\item $u_i$ \mye l'$i$-sima colonna di $U$;
\item $v_i$ \mye l'$i$-sima colonna di $V$;
\end{itemize}
la soluzione di {\purples{minima norma}} risulta
\begin{equation}\label{svdexp}
x^*=\sum_{i=1}^k \frac{u^T_i b}{\sigma_i} v_i.
\end{equation}
Se $k=n$ allora $x^*$ \'e l'unica soluzione nel senso dei minimi quadrati del problema $Ax=b$.
\end{definizione}

\begin{itemize}
\item Si osservi che per calcolare $x^*$ non bisogna risolvere sistemi lineari.
\item Esistono vari algoritmi per calcolare la soluzione mediante fattorizzazione SVD. Col metodo di Golub-Reinsch, la determinazione di $U$, $\Sigma$, $V$ costa {\purples{$4m^2n+8n^3$}} operazioni moltiplicative, mentre con il metodo R-SVD la complessit\'a risulta {\purples{$2m^2n+11n^3$}} (\cite[p.263]{GVL}).
\end{itemize}
\end{frame}





\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati e fattorizzazione SVD}



{\footnotesize{
\begin{esempio}
Risolvere, mediante fattorizzazione SVD, il problema 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]
$$
\end{esempio}

Con tale metodo, ricaviamo come in precedenza la soluzione $x^*$ di $Ax=b$.
\begin{lstlisting}
>> A=[2 5; 1 1 ; 3 2];
>> [U,S,V] = svd(A) % decomposizione SVD (display matrici in "format short")
U =
   -0.8288    0.5527    0.0874
   -0.2161   -0.1722   -0.9611
   -0.5161   -0.8154    0.2621
S =
    6.3866         0
         0    1.7921
         0         0
V =
   -0.5358   -0.8443
   -0.8443    0.5358
>> b=ones(3,1);
>> format long e;
>> x=sum(diag((U(:,1:size(A,2)))'*b./diag(S))*V); % soluzione calcolata con SVD.
>> x'
ans =
     3.358778625954198e-01
     7.633587786259541e-02
>> 
\end{lstlisting}

\bignsp
\smnsp

Si noti che {\purples{\tt{x=sum(diag((U(:,1:size(A,2)))'*b./diag(S))*V);}}} \'e equivalente a (\ref{svdexp}).
}}
\end{frame}







\begin{frame}[allowframebreaks] {Bibliografia}
\bibliographystyle{plain}
\bibliography{bibliography/bibliography} 
\usebibitemtemplate{\color{structure}\insertbiblabel}
\usebibliographyblocktemplate{\color{structure}}{\color{black}}{\color{structure!75}}{\color{structure!75}}


\begin{thebibliography}{99} 


\bibitem{ATK} K. Atkinson, {\it{Introduction to Numerical Analysis}}, Wiley, 1989.
\bibitem{ATKNEW} K. Atkinson, W. Han, {\it{Elementary Numerical Analysis}}, Wiley, 2003.
\bibitem{BINI} D. Bini, M. Capovani, O. Menchi, {\it{Metodi numerici per l'algebra lineare}}, Zanichelli, 1988.
%\bibitem{COM} V. Comincioli, {\reds{Analisi Numerica, metodi modelli applicazioni}}, Mc Graw-Hill, 1990.
% \bibitem{CON} S.D. Conte e C. de Boor, {\reds{Elementary Numerical Analysis, 3rd Edition}}, Mc Graw-Hill, 1980.
\bibitem{EPS} J.F. Epperson, {\it{Introduzione all'analisi numerica. Teoria, metodi, algoritmi.}}, Mc Graw-Hill, 2009.
\bibitem{GVL} G.H.Golub and C.F. Van Loan, {\em Matrix Computations\/}, 3rd edition, John Hopkins University Press, (1996).
\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, {\it{Elementi di calcolo numerico}}, Progetto Leonardo, 1999.
% \bibitem{QUA} A. Quarteroni e F. Saleri, {\reds{Introduzione al calcolo scientifico}}, Springer Verlag, 2006.
%\bibitem{QUALAR} A. Quarteroni, R. Sacco e F. Saleri, {\it{Matematica Numerica}}, Springer Verlag, 1998.
%\bibitem{SULI} A. Suli e D. Mayers, {\reds{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} {\blues{\href{https://it.wikipedia.org/wiki/Decomposizione\_di\_Cholesky}{Wikipedia, Decomposizione di Cholesky}}}.
\bibitem{WIKILU} {\blues{\href{https://it.wikipedia.org/wiki/Decomposizione\_LU}{Wikipedia, Decomposizione LU}}}.
\bibitem{WIKIQR} {\blues{\href{https://it.wikipedia.org/wiki/Decomposizione\_QR}{Wikipedia, Decomposizione QR}}}.
\bibitem{WIKIDD}  {\blues{\href{https://it.wikipedia.org/wiki/Matrice\_a\_diagonale\_dominante}{Wikipedia, Matrice a diagonale dominante}}}.
\bibitem{WIKIDP} {\blues{\href{https://it.wikipedia.org/wiki/Matrice\_definita\_positiva}{Wikipedia, Matrice definita positiva}}}.
\bibitem{WIKICOF} {\blues{\href{Matrice\_invertibile\#Metodo\_della\_matrice\_dei\_cofattori}{Wikipedia, Matrice invertibile: Metodo della matrice dei cofattori}}}.
\bibitem{WIKISOR} {\blues{\href{http://it.wikipedia.org/wiki/Successive\_Over\_Relaxation}{Wikipedia, SuccessiveOverRelaxation}}}.
\bibitem{WIKISUP}  {\blues{\href{https://it.wikipedia.org/wiki/Supercomputer}{Wikipedia, Supercomputer}}}.
\bibitem{WIKIBIN} {\blues{\href{https://it.wikipedia.org/wiki/Teorema\_di\_Binet}{Wikipedia, Teorema di Binet}}}.
\bibitem{WIKILAP}  {\blues{\href{https://it.wikipedia.org/wiki/Teorema\_di\_Laplace}{Wikipedia, Teorema di Laplace}}}.





\end{thebibliography}
\end{frame}






\end{document}











\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati}

\smnsp

{\scriptsize{
\begin{dimostrazione}[Facoltativa]
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 

\smnsp

$$\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 {\reds{$\phi(x) \leq \phi(x+v)$}}.

Da $\phi(x) \leq \phi(x+v)$ ovvero $\phi(x+v)-\phi(x) \geq 0$, sviluppando i calcoli, 
\begin{eqnarray}
\phi(x+v)&=& \|b-A(x+v)\|^2=\|(b-Ax)-Av\|^2= (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)=  2(Ax-b,Av)+(Av,Av). \nonumber
\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}
\end{dimostrazione}
}}
\end{frame}


\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati}

\smnsp

{\scriptsize{
\begin{vuoto}
Se $v=\epsilon u$ con $u$ versore, ovvero vettore tale che $\|u\|_2=1$ allora da $0\leq 2(Ax-b,Av)+(Av,Av)$ 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 $0 \leq 2(Ax-b,Av)+(Av,Av)$,

\smnsp

$$
0 \leq -2\epsilon (Ax-b,Au)+ \epsilon^2 (Au,Au)
$$
ovvero per $\epsilon > 0$ arbitrario

\smnsp

$$
0 \leq -2 (Ax-b,Au)+ \epsilon (Au,Au)
$$
cio\mye 
\begin{equation}\label{the5_4}
0 \leq -2 (Ax-b,Au).
\end{equation}

Per $u$ versore arbitrario, valgono contemporaneamente (\ref{the5_3}) e (\ref{the5_4}), per cui $0 = -2 (Ax-b,Au)$ ovvero


$$
0=(Ax-b,Au)=x^TA^TAu-b^TAu=(A^TAx)^Tu-(A^Tb)^Tu=(A^TAx-A^Tb,u)
$$
da cui essendo $u$ versore arbitrario, necessariamente
 $A^TAx-A^Tb=0$ ovvero {\reds{$
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$.
\end{vuoto}
}}
\end{frame}



\begin{frame}[fragile]
\frametitle{Sistemi sovradeterminati}

\smnsp

{\scriptsize{
\begin{vuoto}
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!).

\smsp

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{vuoto}
}}
\end{frame}


