% -----------------------------*- LaTeX -*------------------------------
\documentclass[12pt,A4]{report}
\usepackage{LC2_class}
\begin{document}

% \course{CS 497}           % optional
% \coursetitle{Geometric Data Structures} % optional
% \semester{Fall 1998}          % optional
% \lecturer{Jeff Erickson}      % optional
\scribe{Luca Schenato}        % required
\lecturenumber{1}           % required, must be a number
\lecturedate{Aprile 20}       % required, omit year

\maketitle

% ----------------------------------------------------------------------


%\begin{danger}
%This is the danger environment.
%\end{danger}

\section{Variabili aleatorie e nozioni preliminari}

Le definizioni di questa sezione si trovano in qualsiasi testo di
processi stocastici lineari ed in particolare nel Capitolo 2 in
\cite{Picci:06}.

\begin{definizione}
Una matrice $\Sigma\in\mathbb{R}^{n\times n}$ si dice semidefinita
positiva se $\Sigma=\Sigma^T$ e $x^T\Sigma x\geq 0, \forall
x\in\mathbb{R}^n$, e viene indicata con $\Sigma\geq 0$. Una matrice
semidefinita positiva si dice positiva se $\Sigma> 0$.
\end{definizione}

\begin{definizione}\label{def:traccia}
La traccia $\textrm{tr}:\mathbb{R}^{n\times n}\rightarrow
\mathbb{R}$ di una matrice $A\in\mathbb{R}^{n\times n}$ e definita
come $\textrm{tr}(A)=\sum_{i=1}^n A_{ii}$ dove $A_{ij}$ e'
elemento $i-j$ della matrice $A$. La traccia gode delle seguenti
proprieta':
\begin{enumerate}
    \item $\textrm{tr}(ABC)=\textrm{tr}(BCA)=\textrm{tr}(CAB)$,
    dove le matrici $A,B,C$ hanno dimensione opportuna
    \item Sia $x\in\mathbb{R^n}$ e $Q\geq 0$, allora
    $||x||_Q^2=x^TQx=\textrm{tr}(x^TQx)=\textrm{tr}(Qxx^T)$
    \item Se $Q\geq 0$ e $P_1\geq P_2$, allora
    $\textrm{tr}(QP_1)\geq\textrm{tr}(QP_2)$
    \item $P_1\geq P_2\Longrightarrow \textrm{tr}(P_1)\geq
    \textrm{tr}(P_2)$, ma non viceversa.
\end{enumerate}
\end{definizione}




\begin{definizione}
Una variabile aleatoria si dice gaussiana ed e' indicata con
$x\sim\mathcal{N}(\mu,\Sigma)$, dove $\mu\in \mathbb{R}^n$ e
$\Sigma\in\mathbb{R}^{n\times n},\Sigma>0$, se la sua densita' di
probabilita' e' data da
$p(x)=\frac{1}{\sqrt{2^n\pi^n\det{\Sigma}}}e^{(x-\mu)^T\Sigma^{-1}(x-\mu)}$.
E' facile verificare che $\mathbb{E}[x]=\mu$ e
$\textrm{Var}(x)=\mathbb{E}[(x-\mathbb{E}[x])(x-\mathbb{E}[x])^T]=\Sigma$.
La condizione $\Sigma>0$ puo essere rilassata a $\Sigma\geq 0$
tramite l'introduzione delle distribuzioni di Dirac.
\end{definizione}

\begin{definizione}
Sia $A\in\mathbb{R}^{n\times m}$. Di definisce pseudoinversa della
matrice $A$, l'unica matrice $A^\dag \in\mathbb{R}^{m\times n}$
che soddisfa le seguenti proprieta:
\[ A^\dag A A^\dag = A^\dag, \ \ A A^\dag A = A, \ \
(AA^\dag)^T=AA^\dag, \ \ (A^\dag A)^T = A^\dag A \] Dalla
definizione segue che
\begin{enumerate}
    \item La pseudoinversa della matrice di soli zeri $0_{nm}\in 0^{n \times m}$,
    e' data da $0_{nm}^\dag= 0_{mn}$.
    \item Se la matrice $A\in\mathbb{R}^{n\times n}$ e'
    invertibile allora $A^\dag = A^{-1}$
    \item Data la matrice $A=\left[
\begin{array}{cc}
  A_{11} &0_{nm} \\
  0_{lr} & 0_{lm}
\end{array}%
\right]$ la sua pseudoinversa e' data da $A^\dag=\left[
\begin{array}{cc}
  A_{11}^\dag &0_{rl} \\
  0_{mn} & 0_{ml}
\end{array}%
\right]$.
\end{enumerate}
\end{definizione}

Elenchiamo ora una serie di semplici proprieta' che risulteranno
molto utili per derivare in maniera chiara e semplice il filtro di
Kalman:

\begin{proposizione}
Siano $x,y$ due variabili aleatorie (non necessariamente gaussiane)
con densita' di probabilita' $p(x,y)$. Sia $g(x,y)$ una generica
funzione misurabile. Allora \[ \mathbb{E}_{xy}[g(x,y)]=
\mathbb{E}_y[\,\mathbb{E}_{x|y}[g(x,y)|y]\,]\]
\end{proposizione}
\begin{dimostrazione}
Con un piccolo abuso di notazione usuamo $\mathbb{E}_x$ per indicare
chiaramente rispetto a quale variabile aleatoria e' calcolata
l'aspettazione. In genere tale pedice non viene indicato ed e'
sottointeso. Infatti
\[ \begin{array}{ll}\mathbb{E}[g(x,y)]&=\mathbb{E}_{xy}[g(x,y)]=\int_x\int_y
g(x,y)p(x,y)dx\,dy =\int_x\int_y g(x,y)p(x|y)p(y)dx\,dy \\
&= \int_y\Big(\int_x  g(x,y)p(x|y)dx\Big)p(y)dy =
\int_y\mathbb{E}_{x|y}[g(x,y)|y]p(y)dy =
\mathbb{E}_y[\,\mathbb{E}_{x|y}[g(x,y)|y]\,] \end{array}\]
\end{dimostrazione}


\begin{proposizione}\label{prop:gaussian}
Siano $x,y$ due variabili congiuntamente gaussiane, cioe':
\[p(x,y)\sim\mathcal{N}\left(\left[
\begin{array}{c}
  \mu_x \\
  \mu_y
\end{array}%
\right],\left[
\begin{array}{cc}
  \Sigma_{xx} &\Sigma_{xy} \\
  \Sigma_{yx} &\Sigma_{yy}
\end{array}%
\right] \right)\]

\begin{enumerate}
\item La combinazione lineare di variabili gaussiane e' ancora una
variabile aleatoria gaussiana. Infatti se $z=Ax+By$ dove $A,B$
sono matrici di dimensioni opportune, allora $z\sim
\mathcal{N}(\mu_z,\Sigma_z)$, dove
\[\mu_z = A\mu_z +B\mu_y, \ \ \Sigma_z = [A \ \ B]\left[
\begin{array}{cc}
  \Sigma_{xx} &\Sigma_{xy} \\
  \Sigma_{yx} &\Sigma_{yy}
\end{array}%
\right] \left[
\begin{array}{c}
  A^T \\
  B^T
\end{array}%
\right] \] \item Il condizionamento di variabili aleatorie
gaussiane e' ancora una variabile aleatoria gaussiana. Infatti
$p(x|y)=\mathcal{N}(\mu_{x|y},\Sigma_{x|y})$. Se $\Sigma_{yy}>0$
\[\mu_{x|y} = \mu_x +\Sigma_{xy}\Sigma_{yy}^{-1}(y-\mu_y), \ \ \Sigma_{x|y} =
  \Sigma_{xx} -\Sigma_{xy}\Sigma_{yy}^{-1}\Sigma_{yx}  \]
Nel caso sia semplicemente semidefinita positiva $\Sigma_{yy}\geq
0$, la media e la varianza sono date da:
\[\mu_{x|y} = \mu_x +\Sigma_{xy}\Sigma_{yy}^\dag(y-\mu_y), \ \ \Sigma_{x|y} =
  \Sigma_{xx} -\Sigma_{xy}\Sigma_{yy}^\dag\Sigma_{yx}  \] dove
  $\Sigma_{yy}^\dag$ e' la pseudoinversa di $\Sigma_{yy}$.
\end{enumerate}
\end{proposizione}
\begin{dimostrazione} Vedere Teorema 2.3 e 2.6 in \cite{Picci:06}.
\end{dimostrazione}

Passiamo ora ad enunciare uno dei piu' importanti teoremi che
riguardano la stima di variabili aleatorie.

\begin{teorema}
Siano $x\in\mathbb{R}^n,y\in\mathbb{R}^m$ due variabili aleatorie
(non necessariamente gaussiane), e la funzione $g:\mathbb{R}^m
\rightarrow \mathbb{R}^n$ una qualsiasi funzione misurabile.
Definiamo $\tilde{x}_g=g(y)$ un generico stimatore di $x$, e sia
$e_g=x-\tilde{x}_g$ il corrispondente errore di stima, che a loro
volta sono variabili aleatorie. Lo stimatore
$\hat{x}=\hat{g}(y)=\mathbb{E}[x|y]$ e' definito \textbf{stimatore a
minima varianza} o stimatore ottimo in media quadratica, poiche'
possiede la seguente proprieta': \[ \mathbb{E}[
(x-\hat{x})(x-\hat{x})^T ]\leq \mathbb{E}[
(x-\tilde{x}_g)(x-\tilde{x}_g)^T ], \forall g()
\]
La precedente proprieta' implica che
\[\mathbb{E}[||x-\hat{x}||_Q^2]\leq
\mathbb{E}[||x-\tilde{x}_g||_Q^2], \forall g(),\forall Q\geq 0.\]
Se indichiamo con $e=x-\hat{x}$ l'errore dello stimatore ottimo,
si ha che \[\mathbb{E}[eg(y)^T]=0,\forall g()\] cioe' la stimatore
e' incorrelato con il le "misure" $y$. Infine lo stimatore ottimo
e' corretto (unbiased) in quanto
\[\mathbb{E}[\hat{x}]=\mathbb{E}[x]\]
\end{teorema}
\begin{dimostrazione}
Questo teorema e' simile al Teorema~2.2 in \cite{Picci:06}.
Definiamo per semplicita' $\Delta g(y)=\hat{g}(y)-g(y)$.
\[ \begin{array}{l}\mathbb{E}[ (x-\tilde{x}_g)(x-\tilde{x}_g)^T
]=\mathbb{E}_{xy}[ (x-\tilde{x}_g)(x-\tilde{x}_g)^T
]=\mathbb{E}_{xy}[ \big(x-g(y)\big)\big(x-g(y)\big)^T ]=
\\=\mathbb{E}_{xy}[
\big(x-g(y)+\hat{g}(y)-\hat{g}(y)\big)\big(x-g(y)+\hat{g}(y)-\hat{g}(y)\big)^T
]\\ =
\mathbb{E}_{y}\big[\mathbb{E}_{x|y}[\big(x-g(y)+\hat{g}(y)-\hat{g}(y)\big)\big(x-g(y)+\hat{g}(y)-\hat{g}(y)\big)^T
|y]\big] \\ =
\mathbb{E}_{y}\big[\mathbb{E}_{x|y}[\big(x-\hat{g}(y)\big)\big(x-\hat{g}(y)\big)^T+\big(x-\hat{g}(y)\big)\Delta
g(y)^T+ \Delta g(y)\big(x-\hat{g}(y)\big)^T+ \Delta g(y)\Delta
g(y)^T|y]\big] \\
=
\mathbb{E}_{y}\big[\mathbb{E}_{x|y}[\big(x-\hat{g}(y)\big)\big(x-\hat{g}(y)\big)^T
|y]+\big(\mathbb{E}_{x|y}[x|y]-\hat{g}(y)\big)\Delta g(y)^T+\Delta
g(y)\big(\mathbb{E}_{x|y}[x|y]-\hat{g}(y)\big)^T + \\
\hspace{5mm}+ \Delta g(y)\Delta g(y)^T \big] \\
= \mathbb{E}[(x-\hat{x})(x-\hat{x})^T]+ \mathbb{E}[\Delta
g(y)\Delta g(y)^T \big] = \mathbb{E}[(x-\hat{x})(x-\hat{x})^T] + M
\geq \mathbb{E}[(x-\hat{x})(x-\hat{x})^T]
\end{array}\]
dove abbiamo usato il fatto che $\mathbb{E}_{x|y}[x|y]=\hat{g}(y)$ e
$M\geq 0, \forall g()$. Questo dimostra la prima disuguaglianza del
teorema. La seconda disuguaglianza segue utilizzando le proprieta'
della traccia indicata nella Definizione~\ref{def:traccia} come
segue
\[ \begin{array}{l}  \mathbb{E}[
(x-\hat{x})(x-\hat{x})^T ]\leq \mathbb{E}[
(x-\tilde{x}_g)(x-\tilde{x}_g)^T ] \Longrightarrow \mathbb{E}[
Q(x-\hat{x})(x-\hat{x})^T ]\leq \mathbb{E}[
Q(x-\tilde{x}_g)(x-\tilde{x}_g)^T ]\\ \Longrightarrow
\big(\textrm{tr}\big( \mathbb{E}[ Q(x-\hat{x})(x-\hat{x})^T
]\big)\leq \big(\textrm{tr}\big(  \mathbb{E}[
Q(x-\tilde{x}_g)(x-\tilde{x}_g)^T ]\big) \Longrightarrow
\mathbb{E}[ ||x-\hat{x}||_Q^2]\leq \mathbb{E}[
||x-\tilde{x}||_Q^2]
\end{array}\]
Passiamo ora a dimostrare l'incorrelazione tra lo stimatore e una
qualsiasi funzione dei dati:
\[ \begin{array}{ll}\mathbb{E}[eg(y)^T]&=\mathbb{E}_{xy}[(x-\hat{x})g(y)^T]=\mathbb{E}_{y}[\mathbb{E}_{x|y}[(x-\hat{x})g(y)^T|y]]=\\
&=\mathbb{E}_{y}[(\mathbb{E}_{x|y}[x|y]-\hat{x})g(y)^T]=\mathbb{E}_{y}[(\hat{x}-\hat{x})g(y)^T]=0\end{array}\]
Infine dimostriamo che lo stimatore ottimo e' corretto (unbiased):
\[\mathbb{E}[\hat{x}]=\mathbb{E}_y[\hat{x}]=\mathbb{E}_y[\mathbb{E}_{x|y}[ x|y]]=\mathbb{E}_{xy}[x]=\mathbb{E}_{x}[x]=\mathbb{E}[x]\]
\end{dimostrazione}

Questo teorema dimostra che lo stimatore ottimo in media
quadratica coincide con l'aspettazione condizionata della
variabile incognita di interesse $x$ rispetto ai dati misurati
$y$. In genere lo stimatore ottimo e' una funzione non-lineare dei
dati e quindi di poca utilita' pratica. Tuttavia, nel caso di
processi gaussiani, lo stimatore ottimo si puo' calcolare in modo
esplicito e risulta essere una funzione lineare\footnote{Da un
punto di vista formale la funzione e' \textit{affine} anche se e'
uso comune usare impropriamente il termine \textit{lineare}.} dei
dati come indicato nella Proposizione~\ref{prop:gaussian}-2 in
quanto $\hat{x}=\mu_{x|y}$. Il filtro di Kalman non e' altro che
un modo efficiente per calcolare l'aspettazione condizionata di
una variabile aleatoria gaussiana condizionata rispetto ad una
serie di misure.


\section{Processi lineari stocastici e stimatori ottimi}

Consideriamo ora un sistema lineare stocastico a tempo discreto
del tipo:
\begin{eqnarray}
x_{k+1} &=& Ax_k +w_k \\
  y_k &=& Cx_k + v_k
\end{eqnarray}
dove $x_k\in\mathbb{R}^{n},y_k\in\mathbb{R}^{m}$,
$A\in\mathbb{R}^{n \times n}$, $C\in \mathbb{R}^{m \times n}$, e
$w_k\sim\mathcal{N}(0,Q)$, $v_k\sim\mathcal{N}(0,R)$ e
\mbox{$x_0\sim\mathcal{N}(\bar{x}_0,P_0)$} sono variabili
aleatorie gaussiane inccorrelate\footnote{E' possibile
generalizzare i risultati che seguono per variabili aleatorie
correlate, ma per semplicita' non verranno considerate.} ed
identicamente distribuite (i.i.d), dove $Q\geq 0,R\geq 0, P_0\geq
0$. Per semplificare la notazione abbiamo indicato $x_k=x(k)$ dove
$k=0,1,2,\ldots$ sono gli istanti temporali. Definiamo anche le
seguenti quantita':
\begin{eqnarray}
\hat{x}_{k|h} &=& \mathbb{E}[x_k \,|\, y_0,\ldots, y_h] \\
P_{k|h} &=& \mathbb{E}[(x_k-\hat{x}_{k|h})(x_k-\hat{x}_{k|h})^T\,
|\, y_0,\ldots, y_h]
\end{eqnarray}
Lo stimatore $\hat{x}_{k|k}$ viene definito come fitro di kalman.
Se $k=h+N$ allora $\hat{x}_{k+N|k}$ viene indicato come predittore
a $N$ passi, mentre se $k=h-N$ allora $\hat{x}_{k+N|k}$ viene
indicato come interpolatore (smoother) a $N$ passi. La matrice
semidefinita positiva $P_{k|h}\geq 0$ corrisponde alla varianza
dell'errore di stima.

Prima di derivare le usuali equazioni del filtro di Kalman e'
interessante notare come il realta' i risultati descritti nella
sezione precedente siano sufficienti per ottenere lo stimatore
ottimo dato le misure fino all'istante $k$. A questo scopo
definiamo le nuove variabili aleatorie
$\textbf{x}_k=(x_{k},x_{k-1},\ldots,x_1)\in\mathbb{R}^{kn}$,$\textbf{y}_k=(y_k,y_{k-1},\ldots,y_1,y_0)\in\mathbb{R}^{(k+1)m}$,
$\textbf{v}_k=(v_k,v_{k-1},\ldots,v_1,v_0)\in\mathbb{R}^{(k+1)m}$
e $\textbf{w}_k=(w_{k-1},w_{k-2},\ldots,w_0)\in\mathbb{R}^{kn}$
per poter riscrivere la dinamica del sistema stocastico in maniera
compatta come segue:\begin{small}
\begin{eqnarray}
\textbf{x}_k &=& \left(%
\begin{array}{ccccc}
  I & A & A^2 &\ldots & A^{k-1} \\
  0 & I & A & \ldots & A^{k-2} \\
  0 & 0 & I & \ldots & A^{k-3} \\
  \vdots & & & & \vdots \\
  0 & 0 & 0 & \ldots & I \\
\end{array}%
\right)\textbf{w}_k+ \left(%
\begin{array}{c}
  A^{k} \\
  A^{k-1} \\
  \vdots \\
  A^2 \\
  A \\
\end{array}%
\right)x_0 = F \textbf{w}_k + G x_0 \\
  \textbf{y}_k &=& \left(%
\begin{array}{ccccc}
  C & CA & CA^2 &\ldots & CA^k \\
  0 & C & CA & \ldots & CA^{k-1} \\
  0 & 0 & C & \ldots & CA^{k-2} \\
  \vdots & & & & \vdots \\
  0 & 0 & 0 & \ldots & C \\
\end{array}%
\right) \textbf{w}_k+ \left(%
\begin{array}{c}
  CA^{k} \\
  CA^{k-1} \\
  \vdots \\
  CA^2 \\
  CA \\
\end{array}%
\right)x_0 + \textbf{v}_k = H \textbf{w}_k + L x_0 +  \textbf{v}_k
\end{eqnarray}
\end{small}
Si vede chiaramente che $\textbf{x}_k,\textbf{y}_k$ sono variabili
aleatorie gaussiane essendo combinazioni lineari di variabili
aleatorie gaussiane. Da un punto di vista teorico e' quindi facile
ottenere lo stimatore ottimo $\hat{x}_k$ calcolando medie e varianze
congiunte delle variabili $\textbf{x}_k,\textbf{y}_k$ come indicato
nella Proposizione \ref{prop:gaussian}-1 e poi calcolando la stima
condizionata $\mathbb{E}[\textbf{x}_k|\textbf{y}_k]$ come indicato
nella Proposizione \ref{prop:gaussian}-2. Lo stimatore ottimo
$\hat{x}_k$ e' dato dai primi $n$ elementi del vettore
$\mathbb{E}[\textbf{x}_k|\textbf{y}_k]$. Si noti come implicitamente
si sia calcolato anche l'interpolazione (smoothing) degli stati
precedenti $\hat{x}_{h|k}$ per $h=1,\ldots,k-1$. Sebbene molto
intuitivo e matematicamente compatto, il calcolo dello stimatore
ottimo necessitano l'inversione di matrici estremamente elevate,
anzi che crescono al crescere di $k$. Tuttavia come si puo' notare
che le precedenti equazioni hanno una struttura particolare. Quella
piu' evidente e' che alcune matrici sono triangolari e quindi
possiedono una struttura che puo' essere utilizzata per un calcolo
numerico piu' efficiente della semplice inversione. Inotre nel caso
tempo-invariante le matrici sono Toepliz e questo permette di
ottenere molti risultati che riguardano la prestazione a regime,
cioe' per $k\rightarrow \infty$. Le equazioni del filtro di Kalman
sono infatti un algoritmo ricorsivo numericamente efficiente che
appunto sfrutta la particolare struttura del problema.

%\subsection{Il filtro di Kalman}

\begin{thebibliography}{30}

\bibitem{Picci:06}
Giorgio Picci.
\newblock {\em Fitraggio Statistico (Wiener, Levinson, Kalman) e Applicazioni}.
\newblock Libreria Progetto, 2006.

\end{thebibliography}

\end{document}
