---
abstract: 'We consider the problem of modeling, estimating, and controlling the latent state of a spatiotemporally evolving continuous function using very few sensor measurements and actuator locations. Our solution to the problem consists of two parts: a predictive model of functional evolution, and feedback based estimator and controllers that can robustly recover the state of the model and drive it to a desired function. We show that layering a dynamical systems prior over temporal evolution of weights of a kernel model is a valid approach to spatiotemporal modeling that leads to systems theoretic, control-usable, predictive models. We provide sufficient conditions on the number of sensors and actuators required to guarantee observability and controllability. The approach is validated on a large real dataset, and in simulation for the control of spatiotemporally evolving function.'
author:
- 'Hassan A. Kingravi, Harshal Maske and Girish Chowdhary [^1]'
title: ' Kernel Controllers: A Systems-Theoretic Approach for Data-Driven Modeling and Control of Spatiotemporally Evolving Processes '
---

=1

Introduction {#sec:intro}
============

Modeling, control, and estimation of spatiotemporally varying systems is a challenging area in controls research. These systems are characterized by dynamic evolution in both the spatial and temporal variables. Some examples of relevant problems include active wing-shaping based control of flexible aircraft, control of heat or particulate diffusion in manufacturing processes, control of rumor spreading across a social network, and tactical asset allocation and control problems in dynamically varying battlespaces. The traditional approach to modeling and control of spatiotemporal systems have relied on Partial Differential Equations (PDEs) [@Alabau2012controlpde], solutions to which are functions that evolve in both space and time. However, PDE models can be limited in situations where exact physics based models of the functional evolution are difficult to formulate, or are inherently limited due to the physical understanding of the process or unknown spatiotemporal interactions [@cressie2011statistics]. Furthermore, the control of PDEs is fundamentally more challenging than the control of finite-dimensional state-space systems because the evolution and control spaces are infinite dimensional Hilbert spaces, as opposed to $\mathbb{R}^n$ [@Alabau2012controlpde].

Accordingly, there has been significant work in approximate modeling of spatiotemporally evolving functions using data-driven or distributed parameter based approximations of PDEs [@cressie2011statistics; @wikle2001kernel]. One way to model spatiotemporally evolving functions is to approximate the function at several sampling locations and build an autoregressive model of the evolution of the function’s output over that grid [@baker2000finite]. The fidelity of these models heavily depends on the number of sampling (equivalently Euclidean grid locations in the independent variable space) locations employed, with a large number of grid locations leading to large-scale state-space models that are difficult to manage. An alternative approach to modeling spatiotemporal functional evolution relies on modeling the correlation between any two sampling locations through a smooth covariance kernel [@cressie2011statistics]. The model of the evolution is then formed through a linear, weighted combination of the kernels, and the hyperparameters of the spatiotemporal covariance kernel and the weights are learned by solving an optimization problem. The power and flexibility of this approach lies in the fact that kernels can be defined over abstract objects, and not just Euclidean grid locations, leading to a modeling technique that is domain agnostic. For example, kernel embeddings are available for graphical models studied in decentralized control [@johansson2014global], images [@ren2012coupled], and many other domains. However, formulating control-usable kernel-based models of spatiotemporal phenomena can be challenging due to the need to take into account the spatiotemporal dependence. Many recent techniques in spatiotemporal modeling have focused on covariance kernel design and associated hyperparameter learning algorithms [@garg2012AAAI; @ma2003nonstationary; @RasmussenWilliams2005; @plagemann2008nonstationary]. The main benefit of careful design of covariance kernels over approaches that simply include time in as an additional input variable [@perez:13:gaussian; @Chowdhary13_ACC2] is that they can account for intricate spatiotemopral couplings. However, there are two key challenges with these approaches: the first challenge is in ensuring the scalability of the model to large scale phenomena. This is difficult due to the fact that the hyperparameter optimization problem is not convex in general, and because when time is used as a kernel input, it is nontrivial to restrict the number of kernels used without losing modeling fidelity [@garg2012AAAI; @ma2003nonstationary; @plagemann2008nonstationary]. The second very important challenge is concerned with the formulation of feasible control strategies utilizing predictive kernel-based models of spatiotemporal phenomena. In particular, when the spatiotemporal evolution is embedded in the design of complex covariance kernel, the resulting model of functional evolution can be highly nonlinear and difficult to utilize in control design.

In this paper, we pursue an alternative systems-theoretic approach to the modeling, control, and estimation of spatiotemporally varying functions that fuses the strengths of kernel methods with systems theory. Our main contribution is to provide a systems-theoretic formulation for approximating, with very high accuracy, spatiotemporal functional evolution by layering a linear dynamical systems prior over temporal evolution of weights of a kernel model. For a class of linearly evolving PDEs, such as the heat diffusion and the wave equation, our approach can lead to a very high-accuracy approximation. This modeling approach is also applicable to data-driven modeling of real-world phenomena, which we demonstrate on a challenging inference problem on satellite data of sea surface temperatures. One benefit of our model is that it can encode spatiotemporal evolution of complex nonlinear surfaces through an Ordinary Differential Equation (ODE) evolving in a Hilbert space induced by the specific kernel choice. Yet, the main benefit of our systems-theoretic approach is that it is highly conducive to control synthesis. To illustrate this fact, we demonstrate that feasible control strategies for a class of spatiotemporally evolving systems can be found using linear control synthesis. In particular, we derive sufficient conditions on the kernel selection to guarantee observability and controllability of the presented model. Furthermore, we demonstrate control synthesis for a diffusion PDE using simple Gaussian kernels distributed uniformly in the input domain.

The outline of this paper is as follows, Section \[sec:observers\] focuses on the development of a systems-theoretic kernel-based model of spatiotemporal evolution, Section \[sec:theory\_results\] presents the main theoretical results, Section \[sec:experimental\_results\] presents modeling results on a real-world large dataset and control synthesis results for a diffusion PDE.

Kernel Controllers {#sec:observers}
==================

This section outlines our modeling framework and presents theoretical results associated with the number of sampling locations required for monitoring functional evolution.

Problem Formulation {#sec:formulation}
-------------------

We focus on predictive inference and control over a time-varying stochastic process, whose mean $f$ is temporally evolving: $$\begin{aligned}
{\label{eq:e:bnp_model1}}
& f_{k+1}& \sim \mathbb{F}(f_k,\eta_k) $$ where $\mathbb{F}$ is a distribution varying with time $t$ and exogenous inputs $\eta$. The theory of reproducing kernel Hilbert spaces (RKHSs) provides powerful tools for generating flexible classes of functions with relative ease, and is thus a natural choice for modeling complex spatial functions [@scholkopf:bk:2002]. Therefore, our focus will be on spatiotemporally evolving kernel-based models, such as Gaussian Processes (GPs). In a kernel-based model, ${k}:{\Omega}\times{\Omega}\to{\mathbb{R}}$ is a positive definite kernel on some compact domain ${\Omega}$ that models the covariance between any two points in the input space. A Mercer kernel [@scholkopf:bk:2002] implies the existence of a smooth map ${\psi}:{\Omega}\to{\mathcal{H}}$, where ${\mathcal{H}}$ is an RKHS with the property $$\begin{aligned}
 {k}(x,y) &= {\langle {\psi}(x),{\psi}(y) \rangle_{{\mathcal{H}}}} = {\langle {\psi}({k}(x,\cdot)),{\psi}({k}(y,\cdot)) \rangle_{{\mathcal{H}}}}. \end{aligned}$$ There is a large body of literature on modeling spatiotemporal evolution in ${\mathcal{H}}$ [@wikle2002kernel; @cressie2011statistics]. A simple approach for spatiotemporal modeling is to utilize both spatial and temporal variables as inputs to the kernel [@perez:13:gaussian; @Chowdhary13_ACC2]. However, this technique leads to an ever-growing kernel dictionary, which is computationally taxing. Furthermore, constraining the dictionary size or utilizing a moving window will occlude the learning of long-term patterns. Periodic or nonstationary covariance functions and nonlinear transformations have been proposed to address this issue [@ma2003nonstationary; @RasmussenWilliams2005]. Furthermore, work in the design of nonseparable and nonstationary covariance kernels seeks to design kernels optimized to environment-specific dynamics, and optimize their hyperparameters in local regions of the input space [@garg2012AAAI; @das2014nonstationary; @plagemann2008nonstationary]. The model of spatiotemporal functional evolution proposed in this paper builds on the idea that modeling the temporal evolution of mixing weights of a kernel model is a valid approach to spatiotemporal modeling. The key idea behind our approach is that the spatiotemporal evolution of a kernel-based model can be directly modeled by tracing the evolution of the mean embedded in a RKHS using switched ordinary differential equations (ODE) when the evolution is continuous, or switched difference equations when it is discrete (Figure \[fig:hilbert\_evolution\]). The advantage of this approach is that it allows us to utilize powerful ideas from systems theory for knowing necessary conditions for functional convergence; furthermore, it offers a natural framework for designing control mechanisms as well.

![Two types of Hilbert space evolutions. Left: the model, represented by the functions $m_i$, switches discretely in the Hilbert space ${\mathcal{H}}$; Right: the evolution of the function $m_t$ is smooth, represented by a solution to an ordinary differential equation in ${\mathcal{H}}$.[]{data-label="fig:hilbert_evolution"}](model.png){width="0.8\columnwidth"}

In this paper, we restrict our attention to the class of functional evolutions $\mathbb{F}$ defined by linear Markovian transitions in an RKHS. While extension to the nonlinear case is possible (and non-trivial), it is not pursued in this paper to help ease the exposition of key ideas. Let $y\in{\mathbb{R}}^{{N}}$ be the measurements of the function available from ${N}$ sensors, ${\mathcal{A}}:{\mathcal{H}}\to{\mathcal{H}}$ be a linear transition operator in the RKHS ${\mathcal{H}}$, and ${\mathcal{K}}:{\mathcal{H}}\to{\mathbb{R}}^{{N}}$ be a linear measurement operator, the model for the infinite-dimensional functional evolution and measurement studied in this paper is: $$\begin{aligned}
{\label{eq:ideal_lin_evol}}
 f_{k+1} &= {\mathcal{A}}f_k + \eta_k\\
 y_k &= {\mathcal{K}}f_k + \zeta_k,\end{aligned}$$ where $\eta_k$ is a zero-mean stochastic process in ${\mathcal{H}}$, and $\zeta_k$ is a Wiener process in ${\mathbb{R}}^{{N}}$. For many kernels, the feature map ${\psi}$ is unknown, and therefore it is necessary to work in the dual space of ${\mathcal{H}}$. For concreteness, we work with an approximate space as follows: given points ${\mathcal{C}}= {\{c_1,\dots,c_{{M}}\}}$, $c_i\in{\Omega}$, we have a dictionary of atoms ${\mathcal{F}_{{C}}}= \begin{bmatrix}{\psi}(c_1) &\cdots & {\psi}(c_{{M}})
\end{bmatrix}$, ${\psi}(c_i)\in{\mathcal{H}}$, the span of which is a strict subspace of the RKHS generated by the kernel. Formally, we have $$\begin{aligned}
 {\mathcal{C}}\mapsto {\ensuremath{{\mathcal{H}}_{{C}}}}:= \operatorname{span}\begin{bmatrix}{\psi}(c_1) &\cdots & {\psi}(c_{{M}})
\end{bmatrix}
\subset {\mathcal{H}}.\end{aligned}$$ This regime, which trades off the flexibility of a truly nonparametric approach for computational realizability, still allows for the representation of rich phenomena. Let ${N}$ represent the number of sampling locations, and ${M}$ be the number of bases generating ${\ensuremath{{\mathcal{H}}_{{C}}}}$. Note that every function $f\in{\ensuremath{{\mathcal{H}}_{{C}}}}$ has an expansion of the form $$\begin{aligned}
{\label{eq:expansion}}
 f(x) = \sum_{i=1}^{{M}} w_i k(c_i,x).\end{aligned}$$ This expansion allows us to write the $w_i$ coordinates in the dual space as vectors $w\in{\mathbb{R}}^{{M}}$. We can show the relation of the function spaces to their Euclidean counterparts via commutative diagrams. Define ${\ensuremath{\mathcal{W}}}:{\ensuremath{{\mathcal{H}}_{{C}}}}\to{\mathbb{R}}^{{M}}$ as the operator that maps the coordinates $w_i$ in [(\[eq:expansion\])]{} to vectors ${w}\in{\mathbb{R}}^{{M}}$, and let ${\ensuremath{\mathcal{W}^{-1}}}:{\mathbb{R}}^{{M}}\to{\ensuremath{{\mathcal{H}}_{{C}}}}$. Note that for finite-dimensional spaces, this inverse map always exists. These definitions allow us to outline the relations between the dynamics operators ${\mathcal{A}}$ and $A$, and the measurement operators ${\mathcal{K}}$ and ${\ensuremath{K}}$ using the commutative diagrams in Figure \[fig:commute\_sysop\_dyn\] and Figure \[fig:commute\_sysop\_meas\] respectively.

The finite-dimensional evolution equations equivalent to [(\[eq:ideal\_lin\_evol\])]{} in the dual space can be formulated as $$\begin{aligned}
  {w}_{k+1} &= A{w}_k + \eta_k {\label{eq:k_measure}}\\
 y_{k} &= {\ensuremath{K}}_{k} w_{k} +\zeta_k, {\label{eq:k_measure1}}\end{aligned}$$ where we have matrices $A\in {\mathbb{R}}^{{M}\times{M}}, \ {\ensuremath{K}}_{k}\in {\mathbb{R}}^{{N}\times{M}}$, the vectors ${w}_k, {w}\in{\mathbb{R}}^{{M}}$, and we have slightly abused notation to let $\eta_k$ and $\zeta_k$ denote their ${\ensuremath{{\mathcal{H}}_{{C}}}}$ counterparts. Note that the measurement operator ${\mathcal{K}}$ is simply a sampling of the function $f$ at an arbitrary set of sensing locations ${\mathcal{X}}= {\{x_1, \dots, x_{{N}}\}}$, where $x_i\in{\Omega}$: we will see how this affects the structure of ${\ensuremath{K}}_k$ momentarily. The equations [(\[eq:ideal\_lin\_evol\])]{} suggest an immediate extension to functional control problems. Pick another dictionary of atoms ${\mathcal{F}_{{D}}}= \begin{bmatrix}{\psi}(d_1) &\cdots & {\psi}(d_{{\ell}})
\end{bmatrix}$, ${\psi}(d_j)\in{\mathcal{H}}$, $d_j\in{\Omega}$, the span of which, denoted by ${\ensuremath{{\mathcal{H}}_{{D}}}}$, is a strict subspace of the RKHS ${\mathcal{H}}$ generated by the kernel. The functional evolution equation is then as follows: $$\begin{aligned}
{\label{eq:ideal_lin_evol_control}}
 f_{k+1} &= {\mathcal{A}}f_k + {\mathcal{B}}{\delta}_k + \eta_k\\
 y_k &= {\mathcal{K}}_k f_k + \zeta_k,\end{aligned}$$ where the control functions ${\delta}_k$ evolve in ${\ensuremath{{\mathcal{H}}_{{D}}}}$, and ${\mathcal{B}}:{\ensuremath{{\mathcal{H}}_{{D}}}}\to{\ensuremath{{\mathcal{H}}_{{C}}}}$. To derive the finite-dimensional equivalent of ${\mathcal{B}}$, we have to work out the structure of the matrix $B$: since ${\ensuremath{{\mathcal{H}}_{{C}}}}$ is not, in general, isomorphic to ${\ensuremath{{\mathcal{H}}_{{D}}}}$, this imposes strict restrictions on $B$. We derive $B$ using least squares using the inner product of ${\mathcal{H}}$. Let ${\delta}= \sum_{j=1}^{{\ell}}{\acute{{w}}}_j{k}(d_j,x)$, and let ${\mathcal{F}_{{C}}}= 
 \begin{bmatrix}{\psi}(c_1) &\cdots & {\psi}(c_{{M}})
\end{bmatrix}$ be the basis for ${\ensuremath{{\mathcal{H}}_{{C}}}}$. Then the projection of $\delta$ onto ${\ensuremath{{\mathcal{H}}_{{C}}}}$ can be derived as $$\begin{aligned}
 \begin{bmatrix}
  \l \delta, {\psi}(c_1) \r_{{\mathcal{H}}}\\
   \vdots\\
  \l \delta, {\psi}(c_{{M}}) \r_{{\mathcal{H}}} 
 \end{bmatrix}
 &= 
 \underbrace{
 \begin{bmatrix}
  {k}(d_1,c_1) & \cdots & {k}(d_{{\ell}},c_1)\\
   \vdots  &\ddots &\vdots\\
  {k}(d_1,c_{{M}}) & \cdots & {k}(d_{{\ell}},c_{{M}})
 \end{bmatrix}}_{{{\ensuremath{K}}_{{C}{D}}}}
 \begin{bmatrix}
  {\acute{{w}}}_1\\
  \vdots\\
  {\acute{{w}}}_{{\ell}}
 \end{bmatrix},\end{aligned}$$ using the reproducing property. This derivation shows that the operator $B = {{\ensuremath{K}}_{{C}{D}}}\in{\mathbb{R}}^{{M}\times{\ell}}$, the kernel matrix between the data ${C}$ generating the atoms ${\mathcal{F}_{{C}}}$ of ${\ensuremath{{\mathcal{H}}_{{C}}}}$ and the data ${D}$ generating the atoms ${\mathcal{F}_{{D}}}$ of ${\ensuremath{{\mathcal{H}}_{{D}}}}$. Using similar arguments, it can be shown that, given sensing locations ${X}= \{x_1,x_2,\dots, x_{{N}}\}$, ${{\ensuremath{K}}_{{D}}}\in{\mathbb{R}}^{{N}\times{\ell}}$ is the kernel matrix between ${X}$ and ${D}$. Thus the finite-dimensional evolution equations equivalent to [(\[eq:ideal\_lin\_evol\_control\])]{} are $$\begin{aligned}
 {w}_{k} &= A{w}_k + {{\ensuremath{K}}_{{C}{D}}}{\acute{{w}}}_k{\label{eq:k_measure_c}}\\
 y_{k} &= {\ensuremath{K}}_{k} w_{k} {\label{eq:k_measure_c1}}.\end{aligned}$$ We pause here to point out just how flexible the kernel-based framework is. First of all, the choice of kernel completely determines the space ${\mathcal{H}}$, which may allow wildly different functional outputs for the same dynamics matrix, as shown in Figure \[fig:kernel\_variation\]. Note also that the dynamical equations [(\[eq:k\_measure\_c\])]{} and [(\[eq:k\_measure\_c1\])]{} are *independent of the choice of domain ${\Omega}$*: different domains with different kernels may result in the same sequence of matrices ${\ensuremath{K}}_k$. This allows our results to hold for any domain over which a kernel can be defined, including examples like graphs, hidden Markov models, and strings, which are not typically studied in the controls literature, at virtually no extra complexity in implementation beyond the design of the actual sensors and actuators. This remarkable fact is why we denote our method to be *domain agnostic*.

\

Since ${\ensuremath{K}}_{k+1}$ is the kernel matrix between the data points and basis vectors, its rows are of the form ${\ensuremath{K}}_{(i)} = 
 \begin{bmatrix}
  {k}(x_i, c_1) & {k}(x_i, c_2) & \cdots & {k}(x_i, c_{{M}}) 
 \end{bmatrix}$. In systems-theoretic language, each row of the kernel matrix corresponds to a *measurement* at a particular location, and the matrix itself acts as a measurement operator. We define the *generalized observability matrix* [@zhou:bk:96] as $$\begin{aligned}
{\label{eq:obs_mat}}
 {\mathcal{O}}_{{\Upsilon}} = 
 \begin{bmatrix}
  {\ensuremath{K}}_{t_1} A^{t_1}\\  
  \cdots\\
  {\ensuremath{K}}_{t_{L}} A^{t_{L}}
 \end{bmatrix},\end{aligned}$$ where ${\Upsilon}= \{t_1, t_2, \dots, t_{{L}}\}$ are the set of instances $t_i$ when we apply the measurement operators ${\ensuremath{K}}_{t_i}$. Note that ${\mathcal{O}}_{{\Upsilon}}\in{\mathbb{R}}^{{N}{L}\times{M}}$. Similarly, we can define the *generalized controllability matrix* as $$\begin{aligned}
{\label{eq:control_mat}}
 {\Psi}_{{\Upsilon}} = 
 \begin{bmatrix}
   {A^{t_1}}^T {{{\ensuremath{K}}_{{D}}}}_{t_1} & {A^{t_2}}^T {{{\ensuremath{K}}_{{D}}}}_{t_2} & \cdots {A^{t_{L}}}^T {{{\ensuremath{K}}_{{D}}}}_{t_{L}}
 \end{bmatrix},\end{aligned}$$ ${\Psi}_{{\Upsilon}}\in{\mathbb{R}}^{{M}\times{L}{\ell}}$ A linear system is said to be observable if ${\mathcal{O}}_{{\Upsilon}}$ has full column rank (i.e. $\mathrm{Rank}~{\mathcal{O}}_{{\Upsilon}}={M}$) and is controllable if ${\Psi}_{{\Upsilon}}$ has full row rank, for ${\Upsilon}= \{0, 1, \dots, {M}-1\}$ [@zhou:bk:96].

Observability guarantees that a feedback-based observer can be designed such that the estimate of $w$ denoted by $\hat {w_k}$ converges exponentially fast to the true state $w_k$. In particular, observability is the necessary condition for the existence of a unique solution to the Riccatti equation required in designing a Kalman filter. Therefore, when $\eta,\zeta$ have a zero mean Gaussian distribution, a Bayes optimal filter can be designed for estimating $w$ if and only if $\mathrm{Rank}~{\mathcal{O}}_{{\Upsilon}}={M}$. Similarly, controllability guarantees that a feedback-based controller can drive the current functional state of the system $f_{k}$ to a reference function $f_{\text{ref}}$, as long as $f_{\text{ref}}\in{\ensuremath{{\mathcal{H}}_{{C}}}}$.

We are now in a position to formally state the spatiotemporal monitoring and control problem considered: Given a spatiotemporally evolving system modeled using [(\[eq:ideal\_lin\_evol\_control\])]{}, choose a set of ${N}$ sensing locations ${\mathcal{X}}= {\{x_1, \dots, x_{{N}}\}}$ and ${\ell}$ actuating locations ${\mathcal{D}}= {\{d_1,\dots,d_{{\ell}}\}}$ such that even with ${N}\ll {M}$ and ${\ell}\ll {M}$, the functional evolution of the spatiotemporal model can be estimated robustly, and driven (controlled) to a reference function $f_{\text{ref}}$. Our approach to solve this problem relies on the design of the measurement operator ${\ensuremath{K}}$ such that the pair $(A,{\ensuremath{K}})$ is observable, and the control operator ${{\ensuremath{K}}_{{D}}}$ such that the pair $(A,{{\ensuremath{K}}_{{D}}})$ is controllable.

Theoretical Results {#sec:theory_results}
-------------------

In this section, we prove results concerning the observability of spatiotemporally varying functions modeled by the functional evolution and measurement equations [(\[eq:k\_measure\])]{} and [(\[eq:k\_measure1\])]{} formulated in Section \[sec:formulation\]. In particular, observability of the system states implies that we can recover the current state of the spatiotemporally varying function using a small number of sampling locations ${N}$, which allows us to 1) track the function, and 2) predict its evolution forward in time. It should be noted that the results are also applicable to controllability of the system in [(\[eq:k\_measure\_c1\])]{} since the structure of the control matrix $K_{CD}$ is also that of a Kernel matrix. We first show in Proposition \[prop:1\] that if $A$ has a full-rank Jordan decomposition, the kernel matrix meeting a condition called *shadedness* (to be defined below) is sufficient for the system to be observable. In Proposition \[prop:2\], we prove a lower bound on the number of sampling locations required for observability which holds for more general $A$. Finally, in Proposition \[prop:3\], we outline a method that achieves this lower bound for certain kernels. Since both ${\ensuremath{K}}$ and ${{\ensuremath{K}}_{{C}{D}}}$ are kernel matrices generated from a shared kernel, these observability results translate directly into controllability results.

To prove our results, we will leverage the spectral decomposition of $A$. Specifically, recall that any matrix $A\in{\mathbb{R}}^{{M}\times{M}}$ is similar to a unique block diagonal matrix ${\Lambda}$ (i.e. $\exists P\in{\mathbb{R}}^{{M}\times{M}}$ invertible such that $A = P{\Lambda}P^{-1}$) whose diagonal blocks are matrices of the form $$\begin{aligned}
{\Lambda}_k(\lambda_i, \lambda_i^*) := 
\begin{bmatrix}
{M}& {I_2}& \cdots &  {0}\\
\vdots & \vdots & \ddots  & {I_2}\\
{0}& {0}& \cdots & {M}{\label{eq:jor_com}}
\end{bmatrix}.\end{aligned}$$ where $(\lambda_i, \lambda_i^*)$ is a complex conjugate eigenvalue of $A$, and $
{M}= \begin{bmatrix}
\mu_1 & \mu_2\\
-\mu_2 & \mu_1
\end{bmatrix}$ and $ {I_2}= \begin{bmatrix}
1 & 0\\
0 & 1
\end{bmatrix}$. Real eigenvalues ${\lambda}_i$ correspond to the case ${M}= {\lambda}_i$ and ${I_2}= 1$. Thus the complete real Jordan form of $A$ will be the appropriate diagonal array of these blocks. If all the eigenvalues $\lambda_i$ are nonzero and real, we say the matrix has a *full-rank Jordan decomposition*.

\[def:shaded\] **(Shaded Kernel Matrix)** Let ${k}:{\Omega}\times{\Omega}\to{\mathbb{R}}$ be a positive-definite kernel on a compact domain ${\Omega}$. Let $C = [c_1,  c_2, \cdots , c_{{M}}\}$, $c_j\in{\Omega}$ be the points generating a finite-dimensional covering of the reproducing kernel Hilbert space ${\mathcal{H}}$ associated to ${k}(x,y)$, and let ${\mathcal{X}}= {\{x_1, \dots, x_{{N}}\}}$, $x_i\in{\Omega}$ Let ${\ensuremath{K}}\in{\mathbb{R}}^{{N}\times{M}}$ be the kernel matrix, where $ {\ensuremath{K}}_{ij} := {k}(x_i,c_j)$. For each row ${\ensuremath{K}}_{(i)} := [k(x_i,c_1), k(x_i,c_2),\dots, k(x_i,c_{{M}})]$, define the set ${\mathcal{I}}_{(i)} := \{\iota_1^{(i)},\iota_2^{(i)},\dots, \iota_{{M}_i}^{(i)}\}$ to be the indices in the kernel matrix row $i$ which are nonzero. Then if $$\begin{aligned}
{\label{eq:shaded_cond}}
    \bigcup_{1 \leq i \leq {N}} {\mathcal{I}}^{(i)} = \{1,2,\dots, {M}\},
    \end{aligned}$$ we denote ${\ensuremath{K}}$ as a *shaded kernel matrix* (see figure \[shaded\_matrix\]).

This condition implies that the null space of the adjoint of ${\ensuremath{K}}$ as a linear operator between Euclidean spaces, i.e. ${\ensuremath{K}}^T:{\mathbb{R}}^{{N}}\to{\mathbb{R}}^{{M}}$ is trivial. Note that, in principle, for the Gaussian kernel, a single row generates a shaded kernel matrix, although this matrix can have many entries that are extremely close to zero. With this definition in place, we can prove the following proposition, which shows that if $A$ has a full-rank Jordan decomposition, a shaded kernel matrix is sufficient to prove observability.

\[prop:1\] Let ${k}:{\Omega}\times{\Omega}\to{\mathbb{R}}$ be a positive definite kernel on a domain ${\Omega}$. Let $C = [c_1,  c_2, \cdots , c_{{M}}\}$, $c_j\in{\Omega}$ be the points generating a finite-dimensional covering of the reproducing kernel Hilbert space ${\mathcal{H}}$ associated to ${k}(x,y)$, and consider the discrete linear system on ${\mathcal{H}}$ given by the evolution and measurement equations [(\[eq:k\_measure\])]{} and [(\[eq:k\_measure1\])]{}. Let $A\in{\mathbb{R}}^{{M}\times{M}}$ be a full-rank Jordan decomposition of the form $A = {P}{\Lambda}{P}^{-1}$, where ${\Lambda}= \operatorname{diag}(\begin{bmatrix}{\Lambda}_1 & {\Lambda}_2 &\cdots & {\Lambda}_{{O}}\end{bmatrix})$, and there are no repeated eigenvalues. Given a set of time instances ${\Upsilon}= \{t_1,t_2,\dots,t_{{L}}\}$, and a set of sampling locations ${\mathcal{X}}={\{x_1, \dots, x_{{N}}\}}$, the system [(\[eq:k\_measure\])]{} is observable if the kernel matrix ${\ensuremath{K}}_{ij} := {k}(x_i,c_j)$ is shaded, ${\ensuremath{K}}^D$, the row vector generated by summing the rows of ${\ensuremath{K}}$, has all nonzero entries, ${\Upsilon}$ has distinct values, and $|{\Upsilon}| \geq {M}$.

To begin, consider a system where $A = {\Lambda}$, with Jordan blocks $\{{\Lambda}_1, {\Lambda}_2, \dots, {\Lambda}_{{O}}\}$ along the diagonal. Then $A^{t_i} = \operatorname{diag}(\begin{bmatrix}{\Lambda}_1^{t_i} & {\Lambda}_2^{t_i} & \cdots & {\Lambda}_{{O}}^{t_i}\end{bmatrix})$. We have that $${\mathcal{O}}_{{\Upsilon}} = 
 \begin{bmatrix}
  {\ensuremath{K}}A^{t_1}\\
  \cdots\\
  {\ensuremath{K}}A^{t_{L}}
 \end{bmatrix}
 =
 \underbrace{
 \begin{bmatrix}
  {\ensuremath{K}}& \cdots & {\ensuremath{K}}\end{bmatrix}}_{{\widehat{\mathbf{K}}}\in{\mathbb{R}}^{{N}\times{M}{L}}}
 \underbrace{
 \begin{bmatrix}
  {\Lambda}_1^{t_1} & \cdots & 0\\
   \vdots & \ddots & \vdots\\
  0 & \cdots & {\Lambda}_{{O}}^{t_1}\\ 
   \hline
   \vdots & \ddots & \vdots\\
   \hline
  {\Lambda}_1^{t_{{L}}} & \cdots & 0\\
   \vdots & \ddots & \vdots\\
  0 & \cdots & {\Lambda}_{{O}}^{t_{L}}   
 \end{bmatrix}}_{{\widehat{\mathbf{A}}}\in{\mathbb{R}}^{{M}{L}\times{M}}}$$ Recall that a matrix’s rank is preserved under a product with an invertible matrix. Design a matrix $U\in{\mathbb{R}}^{{N}\times{N}}$ s.t. ${\widetilde{{\ensuremath{K}}}}:= U{\ensuremath{K}}$ is a matrix with one row vector of nonzeros, and all of the remaining rows as zeros. Then $\operatorname{rank}({\widehat{\mathbf{K}}}{\widehat{\mathbf{A}}}) = \operatorname{rank}(U{\widehat{\mathbf{K}}}{\widehat{\mathbf{A}}})$. Therefore, we have that $${\widetilde{{\ensuremath{K}}}}A^{t_j} = 
 \begin{bmatrix}
  {\widetilde{{\ensuremath{K}}}}_{(1)}\\
  0\\
  \vdots\\
  0
 \end{bmatrix}
 A^{t_j}\\
 = 
 \begin{bmatrix}
  k_{11}{\lambda}_1^{t_j} & \binom{t_j}{1}{\lambda}_1^{t_j-1} + k_{12}{\lambda}_1^{t_j}  & \cdots & k_{1{M}}{\lambda}_{{O}}^{t_j}\\
  0 & 0 & \cdots & 0\\
  \vdots & \vdots & \ddots & 0\\
  0 & 0 & \cdots & 0
 \end{bmatrix}$$ Therefore, following some more elementary row operations encoded by $V\in{\mathbb{R}}^{{M}{L}\times{M}{L}}$, we get that $$V \begin{bmatrix}
    {\widetilde{{\ensuremath{K}}}}& \cdots & {\widetilde{{\ensuremath{K}}}}\end{bmatrix}
   \begin{bmatrix}
    A^{t_1}\\
    \vdots\\
    A^{t_{{L}}}
   \end{bmatrix}
   = 
   \begin{bmatrix}
    \tilde{k}_{11}{\lambda}_1^{t_1} & \cdots & \tilde{k}_{1{M}}{\lambda}_{{O}}^{t_1}\\
    \tilde{k}_{11}{\lambda}_1^{t_2} & \cdots & \tilde{k}_{1{M}}{\lambda}_{{O}}^{t_2}\\
    \vdots & \ddots & 0\\
    \tilde{k}_{11}{\lambda}_1^{t_{{L}}} & \cdots 
    & \tilde{k}_{1{M}}{\lambda}_{{O}}^{t_{{L}}}\\
    \mathbf{0} & \cdots & \mathbf{0}
   \end{bmatrix}\\
   = 
   \begin{bmatrix}
    \boldsymbol{\Phi}\\
    \widehat{\mathbf{0}}
   \end{bmatrix}.$$ If the individual entries $\tilde{k}_{1i}$ are nonzero, and the Jordan block diagonals have nonzero eigenvalues, the columns of $\boldsymbol\Phi$ become linearly independent. Therefore, if ${L}\geq {M}$, the column rank of ${\mathcal{O}}_{{\Upsilon}}$ is ${M}$, which results in an observable system.

To extend this proof to matrices $A = {P}{\Lambda}{P}^{-1}$, note that $${\mathcal{O}}_{{\Upsilon}} = 
 \begin{bmatrix}
  {\ensuremath{K}}A^{t_1}\\
  \cdots\\
  {\ensuremath{K}}A^{t_{L}}
 \end{bmatrix}
 =
  \begin{bmatrix}
  {\ensuremath{K}}{P}{\Lambda}^{t_1}{P}^{-1}\\
  \cdots\\
  {\ensuremath{K}}{P}{\Lambda}^{t_{L}}{P}^{-1}.
 \end{bmatrix}
 =
 \begin{bmatrix}
  {\ensuremath{K}}& \cdots & {\ensuremath{K}}\end{bmatrix}
 \boldsymbol{{P}}
 \boldsymbol{{\Lambda}}^t
 \boldsymbol{{P}^{-1}},$$ where $\boldsymbol{{P}}\in{\mathbb{R}}^{{M}{L}\times{M}{L}}$, $\boldsymbol{{\Lambda}}^t\in{\mathbb{R}}^{{M}{L}\times{M}{L}}$, and $\boldsymbol{{P}^{-1}}\in{\mathbb{R}}^{{M}{L}\times{M}{L}}$ are the block diagonal matrices associated with the system. Since $\boldsymbol{{P}}$ is an invertible matrix, the conclusions about the column rank drawn before still hold, and the system is observable.

When the eigenvalues of the system matrix are repeated, it is not enough for ${\ensuremath{K}}$ to be shaded. The next proposition proves a lower bound on the number of observations required.

\[prop:2\] Suppose that the conditions in Proposition \[prop:1\] hold, with the relaxation that the Jordan blocks $\begin{bmatrix}{\Lambda}_1 & {\Lambda}_2 &\cdots & {\Lambda}_{{O}}\end{bmatrix}$ may have repeated eigenvalues. Let ${r}$ be the number of unique eigenvalues of $A$, and let ${\gamma}({{\lambda}}_i)$ denote the geometric multiplicity of eigenvalue ${{\lambda}}_i$. Then there exist kernels ${k}(x,y)$ such that the lower bound ${l}$ on the number of sampling locations ${N}$ is given by the cyclic index of $A$, which can be computed as $$\begin{aligned}
{\label{eq:geom_mult}}
  {l}= \max_{1\leq i\leq{r}}{\gamma}({{\lambda}}_i).
 \end{aligned}$$

We first prove the lower bound. WLOG, let ${\mathbf{{\ensuremath{K}}}}$ have ${l}-1$ fully shaded, linearly independent rows, and write it as $$\begin{aligned}
  {\mathbf{{\ensuremath{K}}}}&= \begin{bmatrix}
                    k_{11} & k_{12} & \cdots & k_{1{M}} \\
                    \vdots & \vdots & \cdots & \vdots \\
                    k_{({l}-1)1} & k_{({l}-1)2} & \cdots & k_{({l}-1){M}} 
                   \end{bmatrix}.
 \end{aligned}$$ Since the cyclic index is ${l}$, this implies that at least one eigenvalue, say ${{\lambda}}$, has ${l}$ Jordan blocks. Define indices $j_1, j_2, \dots, j_{{l}} \in \{1,2,\dots,{M}\}$ as the columns corresponding to the leading entries of the ${l}$ Jordan blocks corresponding to ${{\lambda}}$. WLOG, let $j_1 = 1$. Using ideas similar to the last proof, we can write the observability matrix as $$\begin{aligned}
  {\mathcal{O}}_{{\Upsilon}}
   &:= 
   \begin{bmatrix}
    {k}_{11}{{\lambda}}^{t_1}  & \cdots & {k}_{1j_{{l}}}{{\lambda}}^{t_1} & \cdots\\
    \vdots & \ddots &\vdots & \ddots\\
    {k}_{11}{{\lambda}}^{t_{{L}}}  & {k}_{1j_{{l}}}{{\lambda}}^{t_{{L}}} & \cdots\\
    \vdots & \ddots & \vdots & \ddots\\
    {k}_{({l}-1)1}{{\lambda}}^{t_1}  \cdots & {k}_{({l}-1)j_{{l}}}{{\lambda}}^{t_1} & \cdots\\
    \vdots & \ddots &\vdots & \ddots\\
    {k}_{({l}-1)1}{{\lambda}}^{t_{{L}}}  & \cdots & {k}_{({l}-1)j_{{l}}}{{\lambda}}^{t_{{L}}} & \cdots
   \end{bmatrix}.
 \end{aligned}$$ Define ${\boldsymbol{{{\lambda}}}}:= \begin{bmatrix}{{\lambda}}^{t_1} & {{\lambda}}^{t_2} & \cdots {{\lambda}}^{t_{{L}}}\end{bmatrix}^T$. Then the above matrix becomes $$\begin{aligned}
  {\mathcal{O}}_{{\Upsilon}}
   &:= 
   \begin{bmatrix}
    {k}_{11}{\boldsymbol{{{\lambda}}}}& \cdots & {k}_{1j_2}{\boldsymbol{{{\lambda}}}}& \cdots & {k}_{1j_{{l}}}{\boldsymbol{{{\lambda}}}}& \cdots\\
    \vdots & \ddots & \vdots & \ddots &\vdots & \ddots\\
    {k}_{({l}-1)1}{\boldsymbol{{{\lambda}}}}& \cdots & {k}_{({l}-1)j_2}{\boldsymbol{{{\lambda}}}}& \cdots & {k}_{({l}-1)j_{{l}}}{\boldsymbol{{{\lambda}}}}& \cdots
   \end{bmatrix}.
 \end{aligned}$$ We need to show that one of the columns above can be written in terms of the others. This is equivalent to solving the linear system $$\begin{aligned}
  \begin{bmatrix}
   {k}_{1j_1}\\
   {k}_{2j_1}\\
   \vdots\\
   {k}_{({l}-1)j_1}
  \end{bmatrix}
  &=
  \begin{bmatrix}
   {k}_{1j_2} & \cdots & {k}_{1j_{{l}}}\\
   {k}_{2j_2} & \cdots & {k}_{2j_{{l}}}\\
   \vdots & \ddots & \vdots\\
   {k}_{({l}-1)j_2} & \cdots & {k}_{({l}-1)j_{{l}}}\\
  \end{bmatrix} 
  \begin{bmatrix}
   c_1\\
   c_2\\
   \vdots\\   
   c_{({l}-1)}
  \end{bmatrix}. 
 \end{aligned}$$ Suppose the kernel matrix on the RHS is generated from the Gaussian kernel. From [@micchelli1984interpolation], it’s known that every principal minor of a Gaussian kernel matrix is invertible, which implies that ${\mathcal{O}}_{{\Upsilon}}$ cannot be observable.

We now prove a sufficient condition for the observability of a system with repeated eigenvalues, but with the condition that the Jordan blocks are trivial.

\[prop:3\] Suppose that the conditions in Proposition \[prop:1\] hold, with the relaxation that the Jordan blocks $\begin{bmatrix}{\Lambda}_1 & {\Lambda}_2 &\cdots & {\Lambda}_{{O}}\end{bmatrix}$ may have repeated eigenvalues, and where ${\Lambda}_i$ are single-dimensional. Let ${l}$ be the cyclic index of $A$. We define $$\begin{aligned}
{\label{eq:empKShadFull}}
  {\mathbf{{\ensuremath{K}}}}= \begin{bmatrix}
                    {\ensuremath{K}}^{(1)}\\
                    \vdots\\
                    {\ensuremath{K}}^{({l})}
                  \end{bmatrix}
 \end{aligned}$$ as the *${l}$-shaded matrix* which consists of ${l}$ shaded matrices with the property that any subset of ${l}$ columns in the matrix are linearly independent from each other. Then system [(\[eq:k\_measure\])]{} is observable if ${\Upsilon}$ has distinct values, and $|{\Upsilon}| \geq {M}$.

A cyclic index of ${l}$ for this system implies that there exists an eigenvalue ${{\lambda}}$ that’s repeated ${l}$ times. WLOG, let ${\mathbf{{\ensuremath{K}}}}$ have ${l}$ fully shaded, linearly independent rows, and, assume that the column indices corresponding to this eigenvalue are $\{1,2,\dots,{l}\}$. Define ${\boldsymbol{{{\lambda}}}}_i:= \begin{bmatrix}{{\lambda}}_i^{t_1} & {{\lambda}}_i^{t_2} & \cdots {{\lambda}}_i^{t_{{L}}}\end{bmatrix}^T$. Then $$\begin{aligned}
  {\mathcal{O}}_{{\Upsilon}}
   &:= 
   \begin{bmatrix}
    k_{11} {\boldsymbol{{{\lambda}}}}_1 & k_{12} {\boldsymbol{{{\lambda}}}}_2 & \cdots & k_{1{M}} {\boldsymbol{{{\lambda}}}}_{{M}}\\
    \vdots & \vdots & \ddots & \vdots\\
    k_{{l}1} {\boldsymbol{{{\lambda}}}}_1 & k_{{l}2} {\boldsymbol{{{\lambda}}}}_2 & \cdots & k_{{l}{M}} {\boldsymbol{{{\lambda}}}}_{{M}}.
   \end{bmatrix}
 \end{aligned}$$ Let ${\boldsymbol{{{\lambda}}}}_1 = {\boldsymbol{{{\lambda}}}}_2 = \cdots {\boldsymbol{{{\lambda}}}}_{{l}} := {\boldsymbol{{{\lambda}}}}$. Focusing on these first ${l}$ columns of this matrix, this implies that we need to find constants $c_1,c_2,\dots, c_{{l}-1}$ s.t. $$\begin{aligned}
  \begin{bmatrix}
   k_{11}\\
   \vdots\\
   k_{{l}1}
  \end{bmatrix}
  &= 
  c_1 
  \begin{bmatrix}
   k_{12}\\
   \vdots\\
   k_{{l}2}
  \end{bmatrix}
  + \cdots + 
  c_{{l}-1} 
  \begin{bmatrix}
   k_{1{l}}\\   
   \vdots\\
   k_{{l}{l}}. 
  \end{bmatrix}
 \end{aligned}$$ However, these columns are linearly independent by assumption, and thus no such constants exist, implying that ${\mathcal{O}}_{{\Upsilon}}$ is observable.

An example of a kernel such that any subset of ${l}$ columns in ${\mathbf{{\ensuremath{K}}}}$ are linearly independent of each other is the Gaussian kernel evaluated on sampling locations ${\{x_1, \dots, x_{{N}}\}}$, where $x_i\in{\Omega}\subset{\mathbb{R}}^d$, and $x_i\neq x_j$.

We can reuse Propositions \[prop:1\], \[prop:2\], and \[prop:3\] to prove kernel controllability results, because the structure of the control matrix ${{\ensuremath{K}}_{{C}{D}}}$ in [(\[eq:k\_measure\_c\])]{} is also that of a kernel matrix.

Kernel ${k}$, basis points ${\mathcal{C}}$, final time step $T_f$. $1)$ Sample data $\{y^i_k\}_{i=1}^{{N}}$ from $f(x,k)$. $2)$ Estimate ${\widehat{w}}_k$ via standard kernel inference procedure. $3)$ Store weights ${\widehat{w}}_k$ in matrix ${\mathcal{W}}\in{\mathbb{R}}^{{M}\times T_f}$. Infer ${\widehat{A}}$ using method of choice (e.g. matrix least squares). Compute the covariance matrix ${\widehat{B}}$ of the observed weights ${\mathcal{W}}$. estimated transition matrix ${\widehat{A}}$, predictive covariance matrix ${\widehat{B}}$.

Kernel ${k}$, basis points ${\mathcal{C}}$, estimated system matrix ${\widehat{A}}$, estimated covariance matrix ${\widehat{B}}$. Compute the cyclic index ${l}$ of ${\widehat{A}}$, and compute [(\[eq:empKShadFull\])]{}, by possibly iterating over ${\mathcal{X}}= {\{x_1, \dots, x_{{N}}\}}$. Use ${\widehat{A}}$, ${\widehat{B}}$, and ${\mathbf{{\ensuremath{K}}}}$ to initialize a state-observer (e.g. Kalman filter (KF)) on ${\ensuremath{{\mathcal{H}}_{{C}}}}$. 1) Sample data $\{y^i_k\}_{i=1}^{{N}}$ from $f(x,k)$. 2) Propagate KF estimate ${\widehat{w}}_{k+1}$ forward to time $t_f$, correct using measurement feedback with $\{y^i_k\}_{i=1}^{{N}}$. 3) Output predicted function $\widehat{f}(x,k+1)$ and predictive covariance of KF.

Kernel ${k}$, basis points ${\mathcal{C}}$, estimated system matrix ${\widehat{A}}$, estimated covariance matrix ${\widehat{B}}$, and function $f_{\text{ref}}$ to drive initial function to. (see Algorithm \[alg:egp\_inf\]). Use Jordan decomposition of ${\widehat{A}}$ to obtain control locations ${\mathcal{D}}$, compute kernel matrix ${{\ensuremath{K}}_{{C}{D}}}\in{\mathbb{R}}^{{\ell}\times{M}}$ between ${\mathcal{D}}$ and ${\mathcal{C}}$, and initialize controller (e.g. LQR) utilizing $({\widehat{A}}, {\widehat{B}})$. 1) Sample data $\{y^i_k\}_{i=1}^{{N}}$ from $f(x,k)$. 2) Utilize observer to estimate ${\widehat{w}}_{k+1}$. 3) Use ${\widehat{w}}_{k+1}$ and $f_{\text{ref}}$ as input to controller to get feedback.

Experimental Results {#sec:experimental_results}
====================

We report experimental results on controlling synthetic and modeling real-world data. All experiments were performed using MATLAB on a laptop running Ubuntu 14.04 with $8$ GB of RAM, and an Intel core i7 processor.

Prediction of global ocean surface temperature
----------------------------------------------

We first analyzed the feasibility of this modeling approach on a large dataset: the $4$ km AVHRR Pathfinder project, which is a satellite monitoring global ocean surface temperature. This data was obtained from the National Oceanographic Data Center. The data consists of longitude-latitude measurements on a 2D domain ${\Omega}\subset[-180,180]\times[-90,90]$; this dataset is challenging, with measurements at over $37$ million coordinates, and several missing pieces of data. The goal was to learn the day and night temperature models $f_k(x,y)\in{\ensuremath{{\mathcal{H}}_{{C}}}}$, where ${\ensuremath{{\mathcal{H}}_{{C}}}}$ was generated using the Gaussian kernel ${k}(x,y) = e^{-(\|x-y\|^2/2{\sigma}^2)}$. We first did a search for the ideal bandwidth ${\sigma}$ for a $304$-dimensional sparse Gaussian process model with a Gaussian kernel. The set of atoms ${\mathcal{F}_{{C}}}$ was determined through a linear independence test based sparsification algorithm [@csato2002sparse]. Once the parameters were chosen, a budgeted GP was learned for each date, resulting in weight vectors ${w}_i, \ i\in\{1,2,\dots,365\}$. We used Algorithm \[alg:egp\_trans\] to infer ${\widehat{A}}$, and applied Algorithm \[alg:egp\_inf\] with ${N}\in \{280,500,1000,2000\}$ chosen randomly in the ${\Omega}$ to track the system state given a random initial condition ${w}_0$. Figures \[fig:pathfinder\_errors\_boxplots\_day\] and \[fig:pathfinder\_errors\_boxplots\] show a comparison of the deviation in percentage of the estimated values from the real data, averaged over all the days. As can be seen, the observer enables the prediction of functional evolution *without needing all the measurements (37 million)*, and performance comparable to sampling over all locations is obtained with sampling only over $2,000$ locations. Note that here, even though the system model is observable at ${N}=280$, since the dynamics are not truly linear in ${\ensuremath{{\mathcal{H}}_{{C}}}}$, we get better performance with more sampling locations. Finally, \[fig:pathfinder\_tr\_times\_boxplots\_day\] and \[fig:pathfinder\_tr\_times\_boxplots\] show that the time required to estimate the state during function tracking with kernel observer are an order of magnitude better than retraining the model every time step (“original” in the figure), with comparable performance.

\

Control of a linear PDE
-----------------------

We then employed kernel controllers for controlling an approximation to the scalar diffusion equation $u_t = bu_{xx}$ on the domain ${\Omega}=[0,1]$, with $b=0.25$. The solution to this equation is infinite-dimensional, so we chose a kernel ${k}(x,y) = e^{-(\|x-y\|^2/2{\sigma}^2)}$, and a set of atoms ${\mathcal{F}_{{C}}}={\{c_1,\dots,c_{{M}}\}}$, $c_i\in{\Omega}$, with ${M}= 25$ generating ${\ensuremath{{\mathcal{H}}_{{C}}}}$, the space approximating ${\mathcal{H}}$, and another set of atoms ${\mathcal{F}_{{D}}}=\{{\psi}(d_1),\dots,{\psi}(d_{{\ell}})\}$, $d_j\in{\Omega}$, ${\ell}=13$, generating the control space ${\ensuremath{{\mathcal{H}}_{{D}}}}$. The number of, and the location of the observations was chosen to be the same as that of the actuation locations $d_j$. First, tests (not reported here) were conducted to ensure that the solution to the diffusion equation is well approximated in ${\ensuremath{{\mathcal{H}}_{{C}}}}$. Algorithm \[alg:egp\_trans\] was then used to infer ${\widehat{A}}$. Figure \[fig:uncontrolled\_pde\] shows an example of an initial function $f_{\text{init}}$ evolving according to the PDE. A reference function $f_{\text{ref}}\in{\ensuremath{{\mathcal{H}}_{{C}}}}$ was chosen to drive $f_{\text{init}}$ to $f_{\text{ref}}$ under the action of the PDE. Finally, Algorithm \[alg:egp\_control\] was used to control the PDE. Figure \[fig:controlled\_pde\] shows $f_{\text{init}}$ being driven to $f_{\text{ref}}$, while Figure \[fig:controlled\_pde\_error\] shows the absolute value of the error between $f_k $ and $f_{\text{ref}}$ as a function of time.

Conclusions
===========

In this paper we presented a systems theoretic approach to the problem of modeling, estimating, and controlling complex spatiotemporally evolving phenomena. Our approach focused on developing a predictive model of spatiotemporal evolution by layering a dynamical systems prior over temporal evolution of weights of a kernel model. The resulting model can approximate PDE evolution, while it has the form of a finite state linear dynamical system. The lower bounds on the number of sampling and actuation locations provided in this paper are non-conservative, as such they provide direct guidance in ensuring robust real-world sensor network and actuation matrix design that must also account for fault-tolerance and reliability considerations.

[^1]: This work was supported in parts by DOE Award Number DE-FE0012173 and AFOSR Award Number FA9550-14-1-0399. Hassan Kingravi is with Pindrop Security, Harshal Maske, and Girish Chowdhary are with the Distributed Autonomous Systems (DAS) laboratory Oklahoma State University,[{hkingravi@pindropsecurity.com, maske@okstate.edu, girish.chowdhary@okstate.edu}]{}
