
\section{Predizione in RKHS e approssimazione finito-dimensionale}
\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 a dire che le risposte impulsive $h^{[i]}$ sono processi Gaussiani con funzione di covarianza $K$.\\
Il problema di stimare le risposte impulsive $h^{[i]}$ da dati misurati può essere formulato come il seguente problema di regolarizzazione di tipo Tikhonov:
\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 faccia soddisfare 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. $$
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.


\section{Stabil Spline Group LAR (SSGLAR) algoritm}
\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 alcun'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 soluzione degli ingressi come descritto prima può essere eseguito dall'algoritmo GLAR applicato al problema di regressione. Si definisce SS-Glar il seguente algotimo:

\subsection{Algoritmo SS-Glar}
\label{sec:AlgoritmoSSGlar}

\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 a \ref{SchwuleForml1};
\end{enumerate}


\subsection{Stima dei hyperparametri}
\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$ di 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 un approccio di validazione come segue: 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, i dati rimanenti l'insieme di validazione. Si esegue l'algoritmo per hyperparametri fissati usando l'insieme di identificazione e successivamente si valuta il modello ottenuto sull'insieme di validazione. Dall'insieme degli hyperparametri $(\beta\in \mathbb{R}^+ ,\gamma \in \mathbb{R}^+ )$ si sceglie un insieme con un numero finito (possibilmente piccolo) di alternative, su quale poi deve essere testato l'algoritmo. Inoltre vengono anche 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$ dell'errore di predizione a un passo:
$$ RMS_k:=\sqrt{\frac{1}{500}\sum_{i=1}^{500}(y_t-\hat{y}_{t|t-k})^2}. $$
Vengono quindi fissati gli hyperparametri e il livello di sparsità con i quali viene rifatto la stima del modello usando l'intero insieme di dati $\{y_t,u_t\}_{t=1,\ldots,N}$.


\section{Stable Spline with Exponential Hyperparametri (SSEH) Algoritm}
\label{sec:StableSplineWithExponentialHyperparametriSSEHAlgoritm}

In questo 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 i hyperparametri $\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{Hyperprobabilità a priori per i hyperparametri}
\label{sec:HyperprobabilitaAPrioriPerIHyperparametri}

La varianza di rumore $\sigma^2$ viene stimata in un passo preliminare usando un modello ARX. Quindi questo parametro viene supposto noto nella stima Bayesiana. I hyperparametri $\beta$, $\theta$ e $\lambda_i$ vengono invece modellati come variabili aleatorie indipendenti. $\beta$ ha una densità di probabilità poco informativa su $\mathbb{R}^+$, mentre $\theta$ è distribuita uniformemente. Ogni $\lambda_i$ è una variabile aleatoria esponenziale con l'inversa della media $\xi \in \mathbb{R}^+$, cioè
$$ 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 interpretato come una variabile aleatoria poco informativa su $\mathbb{R}^+$. Tutti i hyperparametri possono essere raccolti in 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$ dove $\gamma^2$ è definito in \ref{eq:Tikhonov}, possono ora essere tutti differenti, il che aumenta la flessibilità del modello.\\


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

Sia $A_{ti}\in\mathbb{R}^{N\times\infty}$ 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 i componenti non noti. Viene quindi imposto 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*}
Il $\propto$ deriva dal assunzione che il predittore di $u_t$ dato il passato di $u_t^-$ e $y_{t+1}^-$ non dipende da $\zeta$, il $\approx$ invece dall'assunzione che il passato $u_t^-$ e $y_{t}^-$ porta informazione sulla risposta impulsiva del predittore e sui hyperparametri. Questo modello statistico è descritto dalla rete Bayesiana riportata in figura \ref{fig:ReteBayesiana} (sinistra).\\
Come visto sopra, per questioni di complessità computazionale bisogna inoltre troncare la risposta impulsiva $h^{[i]}$ a una certa lunghezza $t_0$.\\

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

Il vettore degli hyperparametri $\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, \ 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^+])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 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'uso pratica dello schema matematico è l'utilizzo di un buon valore iniziale. Sotto viene descritto un schema con il quale si ottiene una soluzione sub ottimale risolvendo il problema di ottimizzazione in $\mathbb{R}^4$ in relazione al modello Bayesiano in figura \ref{fig:ReteBayesiana} (destra).

\begin{enumerate} 
  \item Si ottengono $\{\hat{\lambda_i}\}$, $\hat{\theta}$ e $\hat{\beta}$ risolvendo il problema \ref{eq:MinimizZeta} nella versione modificato seguente
  $$ 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 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 che la i-esima componente di $\bar{\zeta}$ che viene messa a zero e se $J(y^+;\bar{\zeta}) \leq J(y^+;\hat{\zeta})$ si fissa $\hat{\zeta}=\bar{\zeta}$.
\end{enumerate} 



La procedura descritta al passo 2 potrebbe non essere adatta nel caso in cui gli ingressi sono altamente correlati, perché potrebbe risultare sensitiva dal 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 2 con il passo seguente:
\begin{enumerate}
	\item 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{enumerate}
Il valore $\hat{\zeta}$ finale può essere usato come punto iniziale in problema \ref{eq:SchwuleForml2}.


\subsection{Stima delle risposte impulsive del predittore per $\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_{1if^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_{k=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}).


