\section{Nuclei per l'identificazione}


\subsection{Definizione del problema}
\label{sec:DefinizioneDelProblema}

Sia $\{z_t\}_{t\in \mathbb{Z}},\, z_t\in \mathbb{R}^m $ un processo stocastico, stazionario che descrive l'evoluzione temporale di qualche variabile di interesse. Con un piccolo abuso di notazione in questo contesto $z_t$ indica sia la variabile aleatoria che la sua realizzazione. Come mostrato in figure \ref{fig:rete} ogni componente del processo $\{z_t\}_{t\in \mathbb{Z}}$ può essere pensato come un nodo di una grafo. L'obiettivo di questi algoritmi di identificazione è di trovare modelli lineari dinamici che descrivono tutti i componenti di $\{z_t\}$ in funzione degli altri componenti. A questo scopo si definisce $y_t:=z_t^{[i]}$ (la componente i-esima di $z_t$) come l'uscita e tutte le altre componenti $u_t:=z_t^{[-i]}\in \mathbb{R}^{m-1}$ come gli ingressi. Ovviamente il procedimento può essere ripetuto per ogni $i=1,\ldots,m$ ottenendo in questo modo una descrizione completa di tutte le variabili di $z_t$ in funzione delle rimanenti. In questa sezione però per questioni di chiarezza si sceglie un valore specifico per $i$, per esempio $i=1$. In questo caso $z_t$ può essere rappresentato come
\begin{equation*}
z_t = \begin{bmatrix} y_t \\ u_t \end{bmatrix}
\end{equation*}

Si definiscono gli insiemi delle variabili passate all'istante $t$
$$ Y^t=[y_{t-1},y_{t-2},\ldots], \qquad U^t=[u_{t-1},u_{t-2},\ldots] $$

Dalla stazionarietà del processo $\{z_t\}_{t\in \mathbb{Z}}$ segue che anche i processi $\{y_t\}_{t\in \mathbb{Z}}$ e $\{u_t\}_{t\in \mathbb{Z}}$ sono processi stocastici, congiuntamente stazionari che possono quindi essere pensati rispettivamente come ingressi e uscita di un sistema dinamico tempo-invariante. Possono essere quindi descritti dalla seguente relazione:
\begin{equation}
y_t = \sum_{k=1}^{\infty}f_k u_{t-k}+\sum_{k=0}^{\infty}g_k e_{t-k}
\label{eq:SisDin}
\end{equation}
dove $f_k\in\mathbb{R}^{1\times m}$ (avendo assunto assenza di retroazione diretta, cioè $f_0=0$) e $g_k\in\mathbb{R}$ sono coefficienti incogniti della risposta impulsiva e $e_t$ rappresenta un processo di innovazione. Esso è equivalente all'errore di predizione di un passo
\begin{align}
	e_t &:= y_t - \hat{y}_{t|t-1} = y_t - E[y_t|Y^t,U^t] \\
	E[y_t|Y^t,U^t] &:= \sum_{j=1}^{m-1}\left[ \sum_{k=1}^{\infty}h_k^{j}u_{t-k}^{[j]} \right] +\sum_{k=1}^{\infty}h_k^{[m]}y_{t-k}.
\label{predittore_gay}
\end{align}
Le sequenze $h_k:=[h_k^{[1]},h_k^{[2]},\ldots,h_k^{[m]}]\in\mathbb{R}^{1\times m}, k\in \mathbb{Z}^+$ rappresentano i coefficienti della risposta impulsiva del predittore e servono a descrivere sistemi (BIBO) stabili e si ha quindi $h^{[m]}\in l_1(\mathbb{Z}^+)$.\\
Nell'ambito della PEM la predizione di sistemi dinamici può essere vista come predizione dei coefficienti della risposta impulsiva $h_k$ del predittore. Si vuole prestare particolare attenzione a situazioni in cui non tutte le variabili sono necessarie per predire $y_t$, il che in termini matematici equivale ad avere $h_k^{[i]}=0$ per ogni $k\in\mathbb{Z}^+$. In una rappresentazione grafica c'è una arco diretto $u_k^{[i]}$ al nodo $y_k$ se e solo se esiste $k\in\mathbb{Z}^+$ tale che $h_k^{[i]}\neq0$ per $i=1,\ldots,m-1$; inoltre non c'è un capo se e solo se esiste $k\in\mathbb{Z}^+$ tale che $h_k^{[m]}\neq0$. Per esempio il grafo rappresentato in Figura \ref{fig:rete} ha $h_k^{[5]}=h_k^{[1]}=0, \, \forall k\in\mathbb{Z}+$, mentre $h_k^{[2]}=h_k^{[3]}=h_k^{[4]}$ e $h_k^{[6]}$ non sono identicamente nulle. Questo significa che per predire $y_t$ serve il passato di $u^{[2]},\ u^{[3]},\ u^{[4]}$ e di $y$ stesso.\\
\begin{figure}[t]
    \begin{center}
      \includegraphics[width=0.7\textwidth]{figure/rete}
    \caption[rete]
    {\label{fig:rete}
    Grafo che rappresenta l'interconnessione di $m=6$ variabili. Gli archi solidi rappresentano le connessioni nel modello dinamico con il nodo $y_t:=z_t^{[1]}$ da tutti gli altri nodi. L'assenza degli archi da $u_t^{[i]},\ i=1,5$ a $y_t$ indica che $h_k^{[i]}=0,\ i=1,5,\ \forall k \in \mathbb{Z}^+$. Il nodo $y_t$ ha l'arco entrante indicando la dipendenza dall'errore di predizione a un passo $e_t$. Gli archi punteggiati indicano dipendenze tra gli ingressi $u_t$.}
    \end{center}
\end{figure}
In pratica non si sa a priori quali variabili saranno significative per predire l'uscita $y_t$. I metodi PEM e N4SID non effettuano alcuna selezione di ingressi e stimano un modello ``completo'' che tiene conto di tutte le misure. Questo è un motivo per cui questi due algoritmi a volte producono stime scarse.\\
Qui di seguito si presentano degli algoritmi che favoriscono la sparsità e che sono quindi capaci di individuare una struttura del sistema dinamico simile a quella riportata in Figura \ref{fig:rete} ed allo stesso momento riescono a stimare le risposte impulsive $h_k^{[i]}$ nulle.





\subsection{La regolarizzazione}
Un approccio ampiamente utilizzato per ricostruire una funzione da un insieme di misure $\{y_t\}$ consiste nel minimizzare un funzionale di regolarizzazione appartenente ad un RKHS (Reproducing Kernel Hilbert Space) $\mathcal{H}$ associato ad un nucleo $K$ simmetrico e definito positivo.\\
Date $N$ misure, la regolarizzazione ai minimi quadrati in $\mathcal{H}$ stima la funzione sconosciuta come
\begin{equation}
\hat{h}=arg \min_{h}\sum_{t=1}^{N}(y_t-\Gamma_t[h])^2+\eta \|h\|^2_{\mathcal{H}} 
\label{kernel based regularization}
\end{equation}
dove $\{\Gamma_t\}$ sono dei funzionali lineari e limitati in $\mathcal{H}$ legati al modello di misura mentre il parametro $\eta$ è uno scalare positivo che permette di settare il trade-off fra l'errore di predizione e la regolarità della soluzione. La funzione che minimizza (\ref{kernel based regularization}) può essere interpretata da un punto di vista Bayesiano come la stima a minima varianza di $h$ quando $h$ è visto come un processo Gaussiano a media nulla con autocovarianza $K$ e $\{y_t-\Gamma_t[h] \}$ un rumore bianco indipendente da $h$. Spesso la conoscenza a priori è limitata al fatto che i segnali in gioco siano continui e a energia limitata. In questo caso la $h$ viene modellata con un integrale $p$-dimensionale di un rumore bianco di intensità unitaria. L'autocorrelazione $W_p$ di $h$ si può così scrivere come:
\begin{equation}
W_p(s,t)=\int_0^1G_p(s,u)G_p(t,u)du\,
\label{stable spline kernel}
\end{equation}
\begin{equation}
G_p(r,u)=\frac{(r-u)_+^{p-1}}{(p-1)!},\quad (u)_+=
\left\{ \begin{array}{l}
u\quad \text{if} \, u\geq 0 \\
0\quad \text{if} \, u < 0
\end{array}\right.
\label{id spline kernel}
\end{equation}
dove $G_p$ è l'autocovarianza associata all'interpretazione Bayesiana di una spline di regolarizzazione di ordine $p$. Se $p=2$ si ottiene un nucleo a spline cubica.\\
Il problema principale nell'utilizzo di un nucleo come (\ref{id spline kernel}) è che non tiene conto della stabilità della risposta impulsiva. Infatti se l'autocovarianza di $h$ è proporzionale a $W_p$ (come nel caso cubico) la varianza di $h(t)$ è nulla nell'istante iniziale e poi diverge all'infinito. Per questo motivo è stato definito un nucleo apposta per l'identificazione di sistemi lineari in cui $h$ rappresenti la risposta impulsiva di un sistema stabile (con covarianza non nulla nell'instante iniziale che converge esponenzialmente a zero). Questo viene garantito prendendo la covarianza proporzionale alla classe di nuclei
\begin{equation}
K_p(s,t) = W_p(e^{-\beta s},e^{-\beta t})\, \quad s,t\in \mathbb{R}^+
\label{classe di nuclei stabili}
\end{equation}
dove $\beta$ è uno scalare positivo che regola il decadimento della varianza; essendo in pratica un parametro sconosciuto sarà trattato come un iperparametro da stimare a partire dai dati. Tenendo presente le formule precedenti si vede che con $p=2$ l'autocovarianza diventa ciò che in letteratura viene chiamato $Stable\,Spline\,kernel$:
\begin{equation}
K_2(t,\tau)=\frac{e^{-\beta(t+\tau)}e^{-\beta \max(t,\tau)} }{2}-\frac{e^{-3\beta \max(t,\tau)}}{6}
\label{stable spline kernel}
\end{equation}
E' stato dimostrato (vedi[??]) che se $h$ è un processo Gaussiano a media nulla con autocovarianza $K_2$ allora con probabilità uno le realizzazioni di $h$ sono risposte impulsive di sistemi BIBO stabili.\\
E' chiaro che lavorando in pratica con segnali a tempo discreto si parlerà di versioni campionate $h_k$ di $h(t)$ così che $\{h_k\}_{k\in \mathbb{Z}^+}$ verranno viste come realizzazioni campionate del nucleo $K(i,j)$.\\

Si considera ora il problema di stimare il parametro $\theta \in \mathbb{R}^m$ dal modello lineare
\begin{equation}
Y=X\theta +W
\label{arx modello}
\end{equation}
dove si ha $Y\in \mathbb{R}^N$ è il vettore dei dati di uscita, $X\in \mathbb{R}^{N\times m}$ è il vettore di regressione e $W\in \mathbb{R}^N$ è un rumore bianco vettoriale a media nulla con $\mathbb{E}[WW^T]=\sigma^2I$. Ciò che verrà fatto è considerare ogni componente di $\theta$ appartenente ad uno spazio infinito dimensionale RKHS con nucleo $K(s,t)$ come in (\ref{classe di nuclei stabili}). In altre parole $\theta_i$ viene visto come un processo gaussiano con covarianza $K(s,t)$. La matrice di regressione $X$ sarà un operatore lineare contente la storia passata di $u$ e $y$.\\

Il problema di selezionare le variabili più significative può essere affrontato come un problema di ricerca sparsità nella soluzione con alcune importanti differenze. In particolare si tratta ora di un problema a gruppi, poichè per escludere dalla procedura di identificazione delle variabili è necessario che vengano poste a zero delle funzioni di trasferimento. Ogni gruppo di parametri descrive dunque una risposta impulsiva $h^{[i]}$. Per evitare il problema relativo alla selezione del modello e della parametrizzazione non lineare si preferisce non utilizzare algoritmi a gruppi parametrici scegliendo di lavorare in un contesto non parametrico. I parametri vivono dunque in un RKHS che come già detto possono essere trattati mediante probabilità a priori e regolarizzazione.\\
I due approcci che verranno utilizzati sono:
\begin{enumerate}
\item SS-GLAR: Versione a gruppi del problema posto dal Lasso
\footnote{Nel Lasso la scelta del regressore avviene risolvendo il problema nella forma:
\begin{equation}
\hat{\theta}:= arg \min_{\theta}\|Y-X\theta\|_2^2+\gamma_1\|\theta\|_1
\end{equation}
con $\gamma_1$ parametro di regolarizzazione compreso fra 0 e 1.
}
portato in un contesto non-parametrico dove i gruppi $h^{[i]}$ appartengono ad uno spazio infinito dimensionale. Per includere il vincolo del problema è necessario risolvere un problema di regolarizzazione misto $l_1-l_2$.
\item SSEH: un modello gerarchico dove $h^{[i]}$ è un processo Gaussiano con covarianza $\lambda_iK(s,t)$ e il fattore di scala $\lambda_i$ ha una distribuzione esponenziale per favorire la sparsità sullo spazio dei fattori di scala.

\end{enumerate}

A questo punto si rende necessario introdurre della notazione aggiuntiva per poter spiegare gli algoritmi che verranno utilizzati. Si definisce dunque
\begin{equation}
y_t^+:=
\begin{bmatrix}
y_t \\ 
\vdots \\
y_{t+N-1}
\end{bmatrix}
\quad
e_t^+:=
\begin{bmatrix}
e_t\\
\vdots  \\
e_{t+N-1}
\end{bmatrix}
\quad
h:=
\begin{bmatrix}
h^{[1]} \\
\vdots \\
h^{[m]}
\end{bmatrix}
,
\label{def yeh}
\end{equation}
dove $h^{[i]}\in l_1(\mathbb{Z}^+)$, $i=1,\ldots,m$. Il predittore (\ref{predittore_gay}) può essere così riscritto
\footnote{Il prodotto tra matrici semi-infinite va inteso come il limite della sequenza finita, che sotto l'ipotesi $h^{[i]}\in l_1{\mathbb{Z}^+}$, è ben definito.}
come:
\begin{equation}
y_t^+=\underbrace{[A_{t1}\,\ldots \, A_{tm}]}_{:=A_t}h+e_t^+
\label{predittore riscritto}
\end{equation}
con $A_{ti}\in \mathbb{R}^{N \times \infty}$, $i=1,\ldots,m$ definiti da:
\begin{equation}
\begin{array}{rcl}
A_{ti}^{[jk]} & := & u_{t-j-k}^{[i]}\,, \quad i=1,\ldots,m-1\ ,\quad j,k \in \mathbb{Z}^+ \\
A_{tm}^{[jk]} &:=& y_{t-j-k},\quad j\,,k \in \mathbb{Z}^+
\end{array}
\label{gli A}
\end{equation}
e il problema di identificazione può essere dunque visto come la stima di $h$ in (\ref{def yeh}),(\ref{predittore riscritto}) soggetto al vincolo che $h^{[i]}\in l_1(\mathbb{Z}^+)$, $i=1,\ldots,m$. Come già detto lo stimatore dovrà automaticamente selezionare tra $\{u^{[1]},\ldots,u^{[m]},y\}$ le variabili necessarie alla predizione di $y$ e quali no e questo farà in modo che alcune risposte impulsive $\hat{h}^{[i]}$ siano poste esattamente a zero.





\newpage
\section{Predizione in RKHS}
\label{sec:PredizioneInRKHSEApprossimazioneFinitoDimensionale}


Si assume che le risposte impulsive $h^{[i]}$ siano (le versioni campionate) di funzioni in $\mathcal{H}_K$, lo spazio di Hilbert riproduttivo associato al nucleo campionato $K$. Questo corrisponde al fatto che le risposte impulsive $h^{[i]}$ sono processi Gaussiani con funzione di covarianza $K$.\\
Il problema di stimare le risposte impulsive $h^{[i]}$ dai dati misurati può essere formulato come il seguente problema di tipo Tikhonov di regolarizzazione:
\begin{equation}
\{\hat{h}^{[i]}\}=\arg \min_{h^{[i]}\in\mathcal{H}_K}\sum_{t=t_0}^N(y_t-\hat{y}_{t|t-1})^2+ \gamma_2^2\left(\sum_{i=1}^{m}\|h^{[i]}\|^2_{\mathcal{H}_K}\right)
\label{eq:Tikhonov}
\end{equation}
soggetto a
\begin{equation}
\hat{y}_{t|t-1}=\sum_{i=1}^{m-1}\left[\sum_{k=1}^{\infty}h^{[i]}_k u^{[i]}_{t-k}\right] +\sum_{k=1}^{\infty}h^{[m]}_k y_{t-k}.
\label{eq:PredictorTikhonov}
\end{equation}
Il parametro $\gamma_2$ è il cosiddetto parametro di regolarizzazione che deve regolare il trade-off tra fit e regolarità di $h^{[i]}$.\\

Per portare questo problema in dimensione finita si modifica il nucleo $K$ in modo che la realizzazione del nucleo modificato soddisfi la risposta impulsiva alla condizione $h^{[i]}_k=0, \ \forall k>J$. È facile vedere che il nucleo (la covarianza, indicata con $K_J$) che soddisfa la condizione è $K_J(h,k)=0, \  \forall(h,k): \ k>J \ \text{oppure}  \ h>J$.\\
Raccogliendo i coefficienti della risposta impulsiva in un vettore colonna si ha
$$ \underline{h}^{[i]}:=[h^{[i]}_1,h^{[i]}_2,\ldots,h^{[i]}_{t_0}]^T, $$
dove $t_0=J$. Il numero $J$ non ha influenza sul trade-off tra bias e varianza, ma serve solo per ragioni di complessità computazionale. Con un piccolo abuso di notazione si indica con $K_J \in \mathbb{R}^{J\times J}$ la matrice simmetrica, definita positiva formato dai nuclei campionati fino al ritardo $J$. Sotto queste condizione la norma $\|h^{[i]}\|^2_{\mathcal{H}_{K_J}}$ può essere scritta come
$$ \|h^{[i]}\|^2_{\mathcal{H}_{K_J}}=(\underline{h}^{[i]})^T K_J^{-1}\underline{h}^{[i]}. $$
Il predittore a un passo $\hat{y}_{t|t-1}$ diventa
\begin{equation}
\hat{y}_{t|t-1}=\sum_{i=1}^{m-1}\left[\sum_{k=1}^{J}h_k^{[i]}u_{t-k}^{[i]}\right] +\sum_{k=1}^{J}h_k^{[m]}y_{t-k}.
\label{eq:PredictorTikhonovApprox}
\end{equation}
Risulta che sotto la restrizione $h^{[i]}\in\mathcal{H}_{K_J}$ la soluzione del problema \ref{eq:Tikhonov} è
\begin{equation}
\arg \min_{\underline{h}^{[i]}\in\mathbb{R}^{t_0}}\sum_{t=t_0}^N(y_t-\hat{y}_{t|t-1})^2+ \gamma_2^2\sum_{i=1}^{m}\left[(\underline{h}^{[i]})^T K_J^{-1}\underline{h}^{[i]}\right]
\label{eq:TikhonovApprox}
\end{equation}
soggetto a \ref{eq:PredictorTikhonovApprox}. Si vede quindi che il problema può essere riformulato come quello di stimare un modello ARX con un peso di regolarizzazione incluso nel problema della stima.


\newpage
\section{L'algoritmo SSGLAR}
\label{sec:StabilSplineGroupLARSSGLARAlgoritm}

Viene ora discusso come si può risolvere il problema \ref{eq:TikhonovApprox} invocando sparsità nei gruppi $h_{t_0}^{[i]}$ usando l'algoritmo GLAR. Per fare questo si assume che i parametri $\theta$ siano stati fissati; senza alcuna informazione a priori questo parametro viene posto uguale a zero.\\
Il problema \ref{eq:TikhonovApprox} è un problema di regressione lineare con una funzione peso in $l_2$ che può essere riscritto nella forma
\begin{equation}
\overline{y}^+_{t_0}=\sum_{i=1}^{m}\overline{A}_{Ji}h_{t_0}^{[i]}+W 
\label{SchwuleForml1}
\end{equation}
avendo definito
\begin{align*}
\overline{y}^+_{t_0} &:= \begin{bmatrix} y^+_{t_0} \\ 0_{1\times (t_0\cdot m)} \end{bmatrix} \\
\overline{A}_{Ji} &:= \left[A_{Ji}\quad v_i\otimes \Lambda\right] \quad \Lambda:=\gamma_2 K_J^{-1/2} \\
v_i &:= [\underbrace{0\ldots 0}_{i-1} 1 \underbrace{0\ldots 0}_{m-i}]^T	
\end{align*}

La selezione degli ingressi come descritto prima può essere eseguito dall'algoritmo GLAR applicato al problema di regressione. Nella prossima sezione si definisce nel dettaglio questo algoritmo che viene chiamato SS-GLAR.

\paragraph{I passi:}
\label{sec:IPassi}
I passi che compie l'algoritmo SS-GLAR sono:
\begin{enumerate}
	\item stabilire il parametro $\beta$ del nucleo;
	\item stabilire il parametro $\gamma_2$ di \ref{eq:TikhonovApprox} e formare $\overline{A}_{Ji}$ in \ref{SchwuleForml1};
	\item stima di $h^{[i]}$ applicando l'algoritmo GLAR al problema \ref{SchwuleForml1};
\end{enumerate}


\subsection{Stima degli iperparametri}
\label{sec:StimaDeiHyperParametri}
Prima di poter eseguire l'algoritmo appena visto bisogna fissare i seguenti parametri:
\begin{enumerate}
	\item la funzione peso $\gamma_2$ in $l_2$ nell'equazione \ref{SchwuleForml1} (regolarità di $[h^{[i]}_1,h^{[i]}_2,\ldots,h^{[i]}_{t_0}]^T$ nello spazio $\mathcal{H}_{K_{t_0}}$);
	\item il parametro $\beta$ del nucleo (velocità di decrescita del nucleo);
	\item Il numero di blocchi diversi da zero stimati dall'algoritmo GLAR;
\end{enumerate}
Questi parametri possono essere stimati usando il seguente approccio di validazione: disponendo delle misure $\{y_t,u_t\}_{t=1,\ldots,N}$ si suddivide l'insieme in due parti. Il primo $\{y_t,u_t\}_{t=1,\ldots,2N/3}$ viene chiamato l'insieme di identificazione, il secondo con i dati rimanenti viene chiamato l'insieme di validazione. Si esegue l'algoritmo per iperparametri fissati usando l'insieme di identificazione e successivamente si valuta il modello ottenuto sull'insieme di validazione. Dall'insieme degli iperparametri $(\beta\in \mathbb{R}^+ ,\gamma \in \mathbb{R}^+ )$ si sceglie un insieme con un numero finito (possibilmente piccolo) di alternative, sul quale poi si deve testare l'algoritmo. Vengono inoltre testati diversi livelli di sparsità. Alla fine viene scelto il modello identificato che ha dato le miglior prestazioni sul set di validazione in termini di $RMS_1$ (Root Mean Square) dell'errore di predizione a un passo, dove il generico $RMS_k$ ($k$ numero di passi) è definito da:
$$ RMS_k:=\sqrt{\frac{1}{N}\sum_{t=1}^{N}(y_t-\hat{y}_{t|t-k})^2}. $$
Vengono quindi fissati gli iperparametri e il livello di sparsità con i quali viene rifatta la stima del modello usando l'intero insieme di dati $\{y_t,u_t\}_{t=1,\ldots,N}$.

\newpage
\section{L'algoritmo SSEH}
\label{sec:StableSplineWithExponentialHyperparametriSSEHAlgoritm}

In questa sezione il predittore di equazione \ref{eq:PredictorTikhonov} viene formulato in un contesto Bayesiano assumendo che le risposte impulsive $h^{[i]}$ siano processi Gaussiani con covarianza $\lambda_i K(t,s)$. Serve specificare gli iperparametri $\beta$, $\theta$ e $\lambda_i$ e la varianza di rumore $\sigma^2$ che descrivono interamente le densità di probabilità delle misure di $y$ e le risposte impulsive dei predittori $h^{[i]}$.\\


\subsection{Iperprobabilità a priori per gli iperparametri}
\label{sec:HyperprobabilitaAPrioriPerIHyperparametri}

La varianza di rumore $\sigma^2$ viene stimata in un passo preliminare usando un modello ARX per poi supporla nota nella stima Bayesiana. Gli iperparametri $\beta$, $\theta$ e $\lambda_i$ vengono invece modellati come variabili aleatorie indipendenti. $\beta$ ha una densità di probabilità non informativa (cioè inizialmente non si hanno informazioni a priori sul parametro) su $\mathbb{R}^+$, mentre $\theta$ è distribuita uniformemente. Ogni $\lambda_i$ è una variabile aleatoria esponenziale con media $1/\xi \in \mathbb{R}^+$, quindi vale
$$ p(\lambda_i)=\xi exp (-\xi\lambda_i)\chi(\lambda_i\geq0), \qquad i=1,\ldots,m $$
con $\chi$ la funzione indicatrice. $\xi$ può essere interpretata come una variabile aleatoria non informativa su $\mathbb{R}^+$. Tutti gli iperparametri possono essere raccolti in un vettore aleatorio indicato con $\zeta:=[\lambda_1,\ldots,\lambda_m,\theta_1,\theta_2,\beta,\xi]$.\\
Si noti che i parametri $\lambda_i$, che hanno lo stesso ruolo di $\sigma^2 / \gamma^2$ con $\gamma^2$ definito in \ref{eq:Tikhonov}, possono ora essere tutti differenti aumentando così la flessibilità del modello.\\


\subsection{Il modello Bayesiano intero}
\label{sec:IlModelloBayesianoIntero}

Sia $A_{ti}\in\mathbb{R}^{N\times\infty}$ come definito in \ref{gli A} in modo che per le uscite future $y^+:=y_1^+$ risulta
$$ y^+=\left[\sum_{i=1}^{m}A_{1i}h^{[i]}\right]+e^+, $$
dove $e^+ := e^+_1$.\\
In pratica $y^-:=[y_0,y_{-1},y_{-2},\ldots]$ non è mai interamente noto e una soluzione consiste nel porre a zero le componenti non note. Viene quindi imposta la seguente approssimazione:
\begin{align*}
	p(y^+,\{h^{[i]}\},y^-,u|\zeta) &\propto \left[\prod_{t=1}^{N}p(y_t|\{h^{[i]}\},y_t^-,\zeta,u_t^-)\right]p(y^-,\{h^{[i]}\},u^-|\zeta)\\
	&\approx \left[\prod_{t=1}^{N}p(y_t|\{h^{[i]}\},y_t^-,\zeta,u_t^-)\right]p(\{h^{[i]}\}|\zeta)p(y^-,u^-)
\end{align*}
La proporzionalità deriva dall'assunzione che il predittore di $u_t$ dato il passato $u_t^-$ e $y_{t+1}^-$ non dipenda da $\zeta$ mentre l'approssimazione deriva dall'assunzione che il passato $u_t^-$ e $y_{t}^-$ porti informazione sulla risposta impulsiva del predittore e sugli iperparametri. Questo modello statistico è descritto dalla rete Bayesiana riportata in Figura \ref{fig:ReteBayesiana} (sinistra).\\
Come visto sopra, per questioni di complessità computazionale bisogna troncare la risposta impulsiva $h^{[i]}$ a una certa lunghezza $t_0$.\\

\subsection{Stima degli iperparametri}
\label{sec:StimaDegliHyperparametri}

Il vettore degli iperparametri $\zeta$ viene stimato minimizzando la probabilità marginale, cioè la densità di $y^+$, $\zeta$ e $\{h^{[i]}\}$. Può essere dimostrato (si guardi \cite{Chiuso10}) che la stima si ottiene da
\begin{equation}
 \hat{\zeta}=arg \min_{\zeta}J(y^+;\zeta)
 \label{eq:MinimizZeta}
\end{equation}
con $\xi,\beta>0$, $\lambda_i>0, \text{ per } i=1,\ldots,m$ e $J$ dato da
\begin{align}
	J(y^+;\zeta) &=\frac{1}{2} log (det[2\pi V[y^+]])+\frac{1}{2}(y^+)^T(V[y^+])^{-1}y^++\xi\sum_{i=1}^{m}\lambda_i-log(\xi)\\
 V[y^+] &=\sigma^2I_N+\sum_{k=1}^{m}\lambda_i A_{1i} KA_{1i}^T\,.
 \label{eq:SchwuleForml2}
\end{align}

\begin{figure}[t]
    \begin{center}
      \includegraphics[width=0.7\textwidth]{figure/ReteBayesiana}
    \caption[ReteBayesiana]
    {\label{fig:ReteBayesiana}
    Rete Bayesiana che descrive il nuovo modello non parametrico per l'identificazione di sistemi dinamici sparsi, dove $y^-:=[y_{l-1},y_{l-2},\ldots]$ e nel modello ridotto $\lambda:=\lambda_1,\ldots,\lambda_m$.}
    \end{center}
\end{figure}

Una questione molto importante per l'utilizzo pratico di questo procedimento matematico è l'utilizzo di un buon valore iniziale. Nel seguito viene descritto un schema con il quale è possibile ottenere una soluzione sub-ottimale risolvendo il problema di ottimizzazione in $\mathbb{R}^4$ in relazione al modello Bayesiano di Figura \ref{fig:ReteBayesiana} (destra).

\begin{description} 
  \item[(i)] Si ottengono $\{\hat{\lambda_i}\}$, $\hat{\theta}$ e $\hat{\beta}$ risolvendo il problema \ref{eq:MinimizZeta} modificato nella seguente versione:
  $$ \arg \min_{\zeta}[J(y^+;\zeta)-\xi\sum_{i=1}^{m}\lambda_i+log(\xi)] $$
  con $\beta > 0, \ \lambda_1=\ldots=\lambda_m\geq0$
  \item[(ii)] Si fissano $\hat{\xi}=1/\hat{\lambda_1}$ e $\hat{\zeta}:=[\hat{\lambda}_1,\ldots,\hat{\lambda}_m,\hat{\theta}_1,\hat{\theta}_2,\hat{\beta},\hat{\xi}]$. Poi, per $i=1,\ldots,m$ si fissa $\bar{\zeta}=\hat{\zeta}$ tranne la i-esima componente di $\bar{\zeta}$ che viene messa a zero. Se invece $J(y^+;\bar{\zeta}) \leq J(y^+;\hat{\zeta})$ si fissa $\hat{\zeta}=\bar{\zeta}$.
\end{description} 



La procedura descritta al passo (ii) potrebbe non essere adatta nel caso in cui gli ingressi sono altamente correlati, perché potrebbe risultare sensibile all'ordine con cui le componenti di $\hat{\zeta}$ convergono a zero. Per aggirare questo problema si propone di usare l'algoritmo SS-GLAR per effettuare questa selezione. In particolare si può sostituire il passo (ii) con il passo seguente:
\begin{description}
	\item[(ii')] Si fissa $\gamma_2^2=\sigma^2/\hat{\lambda}_1$ e si aggiusta il nucleo visto in sezione \ref{sec:StabilSplineGroupLARSSGLARAlgoritm} con $K(:,:,\hat{\theta},\hat{\beta})$. Si esegue poi l'algoritmo SS-GLAR selezionando il numero di componenti non nulli usando l'approccio di validazione. A questo punto i componenti $\{j_1,\ldots,j_k\}\subseteq\{1,\ldots,m\}$ di $\hat{\zeta}$ vengono annullati, il che equivale a annullare $\hat{h}^{[j]}=0,\ j\in\{j_1,\ldots,j_k\}$.
\end{description}
Il valore $\hat{\zeta}$ finale può essere usato come punto iniziale in problema \ref{eq:SchwuleForml2}.


\subsection{Stima delle risposte impulsive con $\zeta$ noto}
\label{sec:StimaDiRisposteImpulsivePerZetaNoto}

Sotto certi condizioni (vedi \cite{Chiuso10}) si può dimostrare:
$$	\{\hat{h^{[i]}}\}^{m}_{i=1}=arg \min_{\{f^i\in\mathcal{H}\}^{m}_{i=1}} \|y^+ -\sum_{i=1}^{m}A_{1i}f^i\|^2 + \sigma^2\sum_{i=1}^{m}\frac{\|f^i\|^2_{\mathcal{H}_K}}{\lambda^2_i} $$
Inoltre per $k=1,\ldots,m+1$ si ha
$$ \hat{h^{[i]}} = \lambda^2_i K A_{1i}^T c, \quad c=\left(\sigma ^2I_n+\sum_{i=1}^{m+1}\lambda_iA_{1i}KA_{1i}^T\right)^{-1}y^+ $$
Dopo aver ottenute le stime di $\{h^{[i]}\}$ semplice formule possono essere usate per ottenere le risposte impulsive $f$ e $g$ in
\ref{eq:SisDin} e quindi anche il predittore a $k$ passi (per dettagli vedere \cite{Ljung99}).





