% -----------------------------*- LaTeX -*------------------------------
\documentclass[12pt]{report}
\usepackage{../PSC_class}
\usepackage[latin1]{inputenc}
\usepackage[italian]{babel}

\DeclareMathOperator{\tr}{tr}

\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{Ottobre 13, 2010}       % 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 simmetrica $\Sigma = \Sigma^T \in\mathbb{R}^{n\times n}$ si dice semidefinita
positiva, $\Sigma \geq 0$, se $x^T\Sigma x\geq 0\; \forall
x\in\mathbb{R}^n$.
\end{definizione}

\begin{definizione}
Una matrice simmetrica $\Sigma= \Sigma^T \in\mathbb{R}^{n\times n}$ si dice definita
positiva, $\Sigma > 0$, se $x^T\Sigma x\geq 0\; \forall
x\in\mathbb{R}^n$ e $x^T\Sigma x= 0 \Leftrightarrow x= 0$.
\end{definizione}

\begin{definizione}
Una matrice $U \in\mathbb{R}^{n\times n}$ si dice unitaria se $UU^T=U^TU=I$.
\end{definizione}

\begin{proposizione}
Data una matrice $\Sigma \in\mathbb{R}^{n\times n}$ definita posiva $\Sigma>0$, esiste una matrice unitaria $U\in\mathbb{R}^{n\times n}$ che la diagonanalizza, cioe' $U^T\Sigma U=\Lambda$ dove $\Lambda=\mathrm{diag}\{\lambda_1,\lambda_2,\ldots, \lambda_n\}$ con $\lambda_i>0,\forall i$. La colonna $u_i$ della matrice $U=[u_1\ u_2 \ \cdots \ u_n]$ corrisponde all' autovettore relativo all'autovalore $\lambda_i$ di $\Sigma$. Tali colonne rappresentano una base ortonormale di $\mathbb{R}^n$, cioe'  $u_i^Tu_j=0, i\neq j$ ($u_i \perp u_j$) e $||u_i||^2_2=1$.
\end{proposizione}

\begin{proposizione}
Data una matrice $\Sigma \in\mathbb{R}^{n\times n}$ definita posiva $\Sigma>0$, il luogo dei punti $\mathcal{E}=\{x\in\mathbb{R}^n\,|\, x^T\Sigma^{-1} x=1\}$ corrisponde ad una ellisse di dimensione $n$ le cui direzioni degli assi sono dati dagli autovettori di $\Sigma$ e la dimensione lungo tali assi e' data da $\sqrt{\lambda_i}$ dell'autovalore corrispondente.
\end{proposizione}

\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 $\tr(A)=\sum_{i=1}^n A_{ii}$ dove $A_{ii}$ è
elemento sulla diagonale della matrice $A$.\\ La traccia gode delle seguenti
proprietà:
\begin{enumerate}
    \item $\tr(ABC)=\tr(BCA)=\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=\tr(x^TQx)=\tr(Qxx^T)$;
    \item Se $Q\geq 0$ e $P_1\geq P_2$, allora
    $\tr(QP_1)\geq\tr(QP_2)$;
    \item $P_1\geq P_2\Longrightarrow \tr(P_1)\geq
    \tr(P_2)$, ma non viceversa.
\end{enumerate}
\end{definizione}




\begin{definizione}
Una variabile aleatoria si dice gaussiana ed è 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 densità di
probabilità è data da
$$p(x)=\frac{1}{\sqrt{(2\pi)^n\det{\Sigma}}}e^{-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu)}.$$
È 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 (cioè le variabili di incertezza nulla hanno una d.d.p. assimilabile ad un impulso di Dirac nell'origine).
\end{definizione}

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

Elenchiamo ora una serie di semplici proprietà 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 densità di probabilità $p(x,y)$ e 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 indichiamo con $\mathbb{E}_x$
per indicare chiaramente rispetto a quale variabile aleatoria è
calcolata. In genere tale pedice non viene indicato ed è
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\in\mathbb{R}^{n}$ e $y\in\mathbb{R}^{m}$ due variabili congiuntamente gaussiane, cioè:
\[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 aleatorie gaussiane è ancora una
variabile aleatoria gaussiana. Infatti se $z=Ax+By$ dove $A,B$
sono matrici di dimensioni opportune, allora $p(z)\sim
\mathcal{N}(\mu_z,\Sigma_z)$, dove
$$\begin{aligned}
\mu_z &=\mathbb{E}[Ax+By] = A\mu_x +B\mu_y \\
\Sigma_z &= \mathbb{E}\left[\left(Ax+By -A\mu_x -B\mu_y\right)\left(Ax+By -A\mu_x -B\mu_y\right)^T\right] \\ &= \mathbb{E}\left[\left([A \ \ B]\left[
\begin{array}{c}
  x-\mu_x \\
  y-\mu_y
\end{array}\right]\right)\left([A \ \ B]\left[
\begin{array}{c}
  x-\mu_x \\
  y-\mu_y
\end{array}\right]\right)^T\right] \\
&= [A \ \ B]\mathbb{E}\left[\left[
\begin{array}{c}
  x-\mu_x \\
  y-\mu_y
\end{array}\right]\left[
\begin{array}{c}
  x-\mu_x \\
  y-\mu_y
\end{array}\right]^T\right][A \ \ B]^T \\
&= [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]
\end{aligned}$$

\item Il condizionamento di variabili aleatorie gaussiane è ancora una variabile aleatoria gaussiana di d.d.p. $p(x|y)=\mathcal{N}(\mu_{x|y},\Sigma_{x|y})$.\\
Se $\Sigma_{yy}>0$ allora
$$\mu_{x|y} = \mu_x +\Sigma_{xy}\Sigma_{yy}^{-1}(y-\mu_y) \qquad \Sigma_{x|y} =
  \Sigma_{xx} -\Sigma_{xy}\Sigma_{yy}^{-1}\Sigma_{yx} $$
Se $\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) \qquad \Sigma_{x|y} =
  \Sigma_{xx} -\Sigma_{xy}\Sigma_{yy}^\dag\Sigma_{yx}.$$
\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 più 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 sia $g:\mathbb{R}^m
\rightarrow \mathbb{R}^n$ una qualsiasi funzione misurabile.
Definiamo $\tilde{x}_g=g(y)$ stimatore di $x$ ed
$e_g=x-g(y)=x-\tilde{x}_g$ il corrispondente errore di stima, a loro
volta variabili aleatorie. Lo stimatore
$\hat{x}=\hat{g}(y)=\mathbb{E}[x|y]$ è detto \textbf{stimatore ottimo
a minima varianza} o stimatore ottimo in media quadratica, poiché
possiede la seguente proprietà:
\begin{equation}
 \mathbb{E}[
(x-\hat{x})(x-\hat{x})^T ]\leq \mathbb{E}[
(x-\tilde{x}_g)(x-\tilde{x}_g)^T ] \quad \forall g(\cdot).
\end{equation}
La precedente proprietà implica che
\begin{equation}\mathbb{E}[||x-\hat{x}||_Q^2]\leq
\mathbb{E}[||x-\tilde{x}_g||_Q^2] \quad \forall g(\cdot) \quad \forall Q\geq 0.
\end{equation}
Se indichiamo con $e=x-\hat{x}$ l'errore dello stimatore ottimo,
si ha che
\begin{equation}
\mathbb{E}[e\hat{g}(y)^T]=0
\end{equation}
cioè l'errore dello stimatore ottimo è scorrelato dalla sua stima.\\
Lo stimatore ottimo è inoltre corretto (unbiased) in quanto
\begin{equation}
\mathbb{E}[\hat{x}]=\mathbb{E}_y[\hat{g}(y)]=\mathbb{E}[x].
\end{equation}
\end{teorema}
\begin{dimostrazione}
Questo teorema è simile al Teorema~2.2 in \cite{Picci:06}.
Definiamo per semplicità $\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(\cdot)$. Questo dimostra (1.1). Per dimostrare (1.2) si utilizzano le proprietà della traccia, definita in~\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(\tr\big( \mathbb{E}[ Q(x-\hat{x})(x-\hat{x})^T
]\big)\leq \big(\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}[e\hat{g}(y)^T]&=\mathbb{E}_{xy}[(x-\hat{x})\hat{g}(y)^T]=\mathbb{E}_{y}[\mathbb{E}_{x|y}[(x-\hat{x})\hat{g}(y)^T|y]]=\\
&=\mathbb{E}_{y}[(\mathbb{E}_{x|y}[x|y]-\hat{x})\hat{g}(y)^T]=\mathbb{E}_{y}[(\hat{x}-\hat{x})\hat{g}(y)^T]=0\end{array}\]
Infine dimostriamo che lo stimatore ottimo è 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 $x$, rispetto ai dati misurati
$y$. In genere lo stimatore ottimo è una funzione non-lineare dei
dati e quindi di poca utilità pratica. Tuttavia, nel caso di
processi gaussiani, lo stimatore ottimo si può calcolare in modo
esplicito e risulta essere una funzione lineare\footnote{Formalmente la funzione è \textit{affine} anche se è 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 è altro che
un modo efficiente per calcolare l'aspettazione, condizionata rispetto ad una
serie di misure, di
una variabile aleatoria gaussiana.


\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 incorrelate\footnote{Per semplicità non verrà trattato il caso di variabili aleatorie correlate e/o tempo varianti. E' possibile tuttavia generalizzare i risultati che seguono per tali variabili.} 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 quantità:
\begin{eqnarray}
\hat{x}_{k|h} &=& \mathbb{E}[x_k \,|\, y_0,\ldots, y_k] \\
P_{k|h} &=& \mathbb{E}[(x_k-\hat{x}_{k|h})(x_k-\hat{x}_{k|h})^T\,
|\, y_0,\ldots, y_k]
\end{eqnarray}
Lo stimatore $\hat{x}_{k|h}$ viene interpretato come un filtro $\hat{x}_{k|k}$ se $k=h+N$, un predittore a $N$ passi $\hat{x}_{k+N|k}$ se $k=h+N$ e un interpolatore (\emph{smoother}) a $N$ passi $\hat{x}_{k|k+N}$ se $k=h-N$.
La matrice semidefinita positiva $P_{k|h}\geq 0$ corrisponde invece alla varianza dell'errore di stima.

Prima di derivare le usuali equazioni del filtro di Kalman è interessante notare come in realtà i risultati descritti nella sezione precedente siano sufficienti per ottenere lo stimatore ottimo date 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-1} \\
  0 & C & CA & \ldots & CA^{k-2} \\
  0 & 0 & C & \ldots & CA^{k-3} \\
  \vdots & & & & \vdots \\
  0 & 0 & 0 & \ldots & C \\
  0 & 0 & 0 & \ldots & 0
\end{array}%
\right) \textbf{w}_k+ \left(%
\begin{array}{c}
  CA^{k} \\
  CA^{k-1} \\
  \vdots \\
  CA^2 \\
  CA \\
  C
\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 è 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$ è 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 (\emph{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 necessita l'inversione di matrici
estremamente elevate, che crescono di dimensione all'aumentare di $k$.
Tuttavia, come si può notare, le precedenti equazioni hanno una
struttura particolare: quella più evidente è che alcune matrici
sono triangolari inferiori. Le equazioni del filtro di Kalman sono un
algoritmo ricorsivo numericamente efficiente che sfrutta appunto
la particolare struttura del problema, riducendo i tempi e le difficoltà di computazione.

\bibliographystyle{plain}
\bibliography{../PSC_bibliografia}

\end{document}
