An Efficient, Expressive and Local Minima-free Method for Learning Controlled Dynamical Systems

TL;DR

This paper introduces RFF-PSRs, a novel method for modeling controlled dynamical systems that combines kernel embeddings, moment-matching, and local optimization for efficient and expressive learning of complex dynamics.

Contribution

It presents a new Predictive State Representation framework using Random Fourier Features, enabling non-parametric, expressive, and theoretically sound modeling of controlled systems.

Findings

Effective modeling of continuous non-linear dynamics.

Efficient learning algorithm with theoretical guarantees.

Experimental results demonstrate superior performance.

Abstract

We propose a framework for modeling and estimating the state of controlled dynamical systems, where an agent can affect the system through actions and receives partial observations. Based on this framework, we propose the Predictive State Representation with Random Fourier Features (RFFPSR). A key property in RFF-PSRs is that the state estimate is represented by a conditional distribution of future observations given future actions. RFF-PSRs combine this representation with moment-matching, kernel embedding and local optimization to achieve a method that enjoys several favorable qualities: It can represent controlled environments which can be affected by actions; it has an efficient and theoretically justified learning algorithm; it uses a non-parametric representation that has expressive power to represent continuous non-linear dynamics. We provide a detailed formulation, a theoretical analysis and an experimental evaluation that demonstrates the effectiveness of our method.

Tables

Table 1. Table 1: Comparison between proposed RFF-PSR and existing system identification methods in terms of the type of systems they can model as well as their computational efficiency and statistical consistency. The table should be interpreted as follows: for each method there exists an instantiation that simultaneously satisfies all properties marked with ✓ ✓ \checkmark but there is no instantiation that is guaranteed to satisfy the properties marked with × \times . A method is scalable if computational and memory costs scale at most linearly with the number of training examples. For RFF-based methods, consistency is up to an approximation error that is controllable by the number of features  (Rahimi and Recht, 2008 ) .

Method	Actions	Continuous	Non-linear	Partially observable	Scalable	Consistent
Non-linear ARX	✓	✓	✓	$\times$	✓	✓
N4SID for Kalman Filter	✓	✓	$\times$	✓	✓	✓
Non-convex optimization (e.g. EM)	✓	✓	✓	✓	✓	$\times$
Gram-Matrix (e.g. HSE-PSR)	✓	✓	✓	✓	$\times$	✓
Spectral PSR/POMDP	✓	$\times$	✓	✓	✓	✓
Reduction to Supervised Learning	$\times$	✓	✓	✓	✓	✓
RFF-PSR	✓	✓	✓	✓	✓	✓

Equations

$$\mathbb{E}[A|B=b,C=c]=V_{A\mid B;c}b\forall b,c$$

$$V_{A,B\mid C,D}{\times }_{C}c{\times }_{D}d$$

$$P_t=W_{\mathrm{sys}}(Q_t)$$

$$Q_{t+1}=f_{\mathrm{filter}}(W_{\mathrm{sys}}(Q_t),o_t,a_t)$$

$${\bar{P}}_t=W_{\mathrm{sys}}({\bar{Q}}_t)$$

$$\arg \underset{f\in \mathcal{F}}{\min }\sum _{t=1}^{T}l(f({\psi }_t^h),{\psi }_t^{oa})+R(f)$$

$$\arg \underset{f\in \mathcal{F}}{\min }\sum _{t=1}^{T}l(f({\psi }_t^h){\psi }_t^a,{\psi }_t^o)+R(f)$$

$${\bar{P}}_{n,t}\approx W_{\mathrm{sys}}({\bar{Q}}_{n,t})\forall n,t$$

$$\Vert ({\hat{W}}_{\mathrm{sys}}-W_{\mathrm{sys}})(Q)\Vert \le$$

$$O\left({\eta }_{\delta ,N}\left((1/{\lambda }_2)+(1/{\lambda }_2^{\frac{3}{2}})\sqrt{1+\sqrt{\frac{\log (1/\delta )}{N}}}\right)\right)$$

$$+O\left(\frac{\log (1/\delta )}{\sqrt{N}}\left(\frac{1}{{\lambda }_2}+\frac{1}{{\lambda }_2^{\frac{3}{2}}}\right)\right)+O\left(\sqrt{{\lambda }_2}\right),$$

$${\eta }_{\delta ,N}=O_p\left(\frac{1/\sqrt{N}+{\lambda }_1}{c+{\lambda }_1}\right),$$

$$Q_t$$

$$P_t$$

$$V_{{\psi }_{t+1}^o\mid {\psi }_{t+1}^a,{\varphi }_t^o,{\varphi }_t^a;h_t^{\infty }}$$

$$=V_{{\psi }_{t+1}^o\otimes {\varphi }_t^o\mid {\psi }_{t+1}^a,{\varphi }_t^a;h_t^{\infty }}{\times }_{{\varphi }_t^o}(C_{o_t{o}_t\mid h_t^{\infty },a_t}+\lambda I{)}^{-1}$$

$$Q_{t+1}$$

$$=V_{{\psi }_{t+1}^o\mid {\psi }_{t+1}^a,{\varphi }_t^o,{\varphi }_t^a;h_t^{\infty }}{\times }_{{\varphi }_t^o}{\varphi }_t^o{\times }_{{\varphi }_t^a}{\varphi }_t^a$$

$${\bar{Q}}_t=W_{{\psi }_t^o|{\psi }_t^a;{\psi }_t^h}=C_{{\psi }_t^o{\psi }_t^a|{\psi }_t^h}(C_{{\psi }_t^a{\psi }_t^a|{\psi }_t^h}+\lambda I{)}^{-1}$$

$$T^{\ast }=\arg \underset{T\in \mathcal{F}}{\min }\Vert (T{\times }_{{\psi }_t^h}{\psi }_t^h{\times }_{{\psi }_a^t}{\psi }_a^t)-{\psi }_o^t{\Vert }^2+\lambda \Vert T{\Vert }_{HS}^2,$$

$${\bar{Q}}_t=T^{\ast }{\times }_{{\psi }_t^h}{\psi }_t^h$$

$$q_{t+1}$$

$$\mathbb{E}[o_t|q_t]$$

$${\dot{x}}_1(t)$$

$$-0.5\cos (x_1(t))a(t)$$

$${\dot{x}}_2(t)$$

$$o(t)$$

$$O^k{\times }_s{e}_i{\times }_a(e_j\otimes e_l)=\mathrm{vec}(O_{ij}\otimes (O^{k-1}{\times }_a{e}_l{\times }_s{T}_{ij}))$$

$$Q_t=O^k{\times }_s{s}_t$$

$$P_t=O^{k+1}{\times }_s{s}_t$$

$$\mathrm{vec}(Q_t)={\tilde{O}}^k{s}_t$$

$$P_t=O^{k+1}{\times }_s{s}_t=O^{k+1}{\times }_s(({\tilde{O}}^k{)}^{+}\mathrm{vec}(Q_t)),$$

$$x_t$$

$$o_t$$

$$\mathbb{E}[o_{t:t+k-1}\mid s_t,a_{t:t+k-1}]={\Gamma }_k{s}_t+U_k{a}_{t:t+k-1},$$

$${\Gamma }_k$$

Full text

1 Introduction

Controlled dynamical systems, where an agent can influence an environment through actions and receive partial observations, emerge in numerous applications in robotics and automatic control. Modeling and learning these systems from data is of great importance in these fields.

The general problem of learning dynamical systems from data (also known as system identification) has been extensively studied and several methods were proposed to tackle it. However, having an expressive, efficient and consistent method for non-linear controlled systems remains an open problem.

Many system identification methods rely on likelihood-based optimization or sampling using EM, MCMC or gradient descent. which makes them prone to poor local optima. There is another class of methods that alleviates the local optima problem and offers a tractable and statistically consistent approach to system identification. These methods, usually referred to as spectral algorithms, have two key properties in common: predictive representation and method of moments. Instead of the state being a latent variable, they represent the estimated state by the expectation of sufficient statistics (or features) of future observations; and they use method of moments to learn model parameters from data.111There is a class of spectral algorithms that maintains the latent variable view. This is exemplified by tensor decomposition methods (Anandkumar et al., 2014).

Initially introduced for linear-Gaussian systems (van Overschee and de Moor, 1996), these algorithms have been extended to discrete systems (Hsu, Kakade, and Zhang, 2009; Siddiqi, Boots, and Gordon, 2010; Boots, Siddiqi, and Gordon, 2011) and then to general smooth continuous systems (Boots, Gretton, and Gordon, 2013). More recently, it has been shown that a wide class of spectral learning algorithms for uncontrolled systems are instances of a two-stage regression framework (Hefny, Downey, and Gordon, 2015), where system identification is reduced to solving a set of regression problems. This framework allows for seamless integration of compressing non-linearities, sparsity (Xia, 2016) and online learning (Venkatraman et al., 2016) into system identification, and for establishing theoretical guarantees by leveraging the rich literature on supervised regression.

Unfortunately, the formulation in (Hefny, Downey, and Gordon, 2015) is limited to uncontrolled systems. On the contrary, we are interested in controlled systems, where the user can affect the system through actions. This gives rise to a key issue: the policy that determines the actions can change at test time. For this reason, the representation of the predictive state must be independent of the training policy and therefore must encode a conditional distribution of future observations given future actions. To adopt such a representation into a practical method that retains the benefits of the two-stage regression formulation, there are a number of challenges that need to be tackled.

First, we need a suitable state representation and dynamics model that can be used to represent a wide class of controlled dynamical systems while ensuring the learning problem remains tractable. Second, we would like to benefit from the two-stage regression view of (Hefny, Downey, and Gordon, 2015) to facilitate model formulation. However, a key assumption in that work is that future observations provide an unbiased estimate of the predictive state, which is not true when the state is a conditional distribution. Third, having a different state representation and having action policy playing a key role on determining the training data require a different theoretical analysis than the one in (Hefny, Downey, and Gordon, 2015). Fourth, because they are based on method of moments, two stage regression models are statistically inefficient. Having the ability to refine the model using local optimization can lead to significant gains in predictive performance.

In this work we address these challenges by combining ideas from two-stage regression, kernel embedding and approximation and gradient descent with backpropagation through time to develop RFF-PSRs. Overall, RFF-PSRs enjoy a number of advantages that, to our knowledge, are not attained by existing system identification methods. We summarize these advantages in Table 1.

In summary, the contributions of this work are as follows: (1) We develop a two-stage regression framework for controlled dynamical systems that admits tractable learning (Sections 3-4). (2) Through the two-stage regression view, we provide theoretical guarantees on learning the parameters of a controlled system (Section 4.4). (3) We use the extended formulation to construct RFF-PSRs, an efficient approximation of kernel-based predictive state representations (HSE-PSRs) (Section 5). (4) We provide a means to refine the parameters of a controlled dynamical system and apply it to our proposed RFF-PSR model (Section 5.5). (5) We demonstrate the advantages of our proposed method through synthetic and robot simulation experiments (Section 6).

2 Related Work

Developing tractable and consistent algorithms for latent state dynamical systems dates back to spectral subspace identification algorithms for Kalman filters(van Overschee and de Moor, 1996). At their heart, these algorithms represent the state as a prediction of the future observations conditioned on history and future actions, and use matrix factorization to obtain a basis for the state.

This notion of the state as a prediction is the basis of predictive state representations (PSRs) (Singh, James, and Rudary, 2004), where the state is represented by the success probabilities of a number of tests. A test succeeds if a specified sequence of test observations is observed when administering a specified sequence of test actions.

Noting that the state and parameters of a PSR are defined up to a similarity transformation has led to a family of tractable and consistent spectral algorithms for learning PSRs (Rosencrantz and Gordon, 2004). More recently, Boots, Gretton, and Gordon (2013) proposed a generalization of PSRs in a reproducing kernel Hilbert space (RKHS). This Hilbert space embedding of PSRs (HSE-PSRs) is able to represent systems with continuous observations and actions while still offering a tractable and consistent learning algorithm. HSE-PSRs, however, use a Gram matrix formulation, whose computational and storage requirements can grow rapidly with the size of training data. A finite dimensional approximation for non-linear PSRs was proposed by Boots and Gordon (2011). However, it can be thought of as an approximation of HSE-HMMs (Song et al., 2010) with actions, a method that has poor theoretical guarantees (Boots, Gretton, and Gordon, 2013). In addition, Boots and Gordon (2011) did not provide examples of how to apply the proposed model to controlled processes with continuous actions. In contrast, the model we propose is an approximation of HSE-PSRs, which is a more principled generalization of PSRs as it performs true Bayesian inference in the RKHS. In addition, our proposed learning algorithm incorporates a local optimization procedure that we demonstrate to be very effective.

We use a reduction of system identification to supervised regression. Similar reductions has been proposed in the literature (Langford, Salakhutdinov, and Zhang, 2009; Hefny, Downey, and Gordon, 2015; Boots and Gordon, 2011; Venkatraman et al., 2016; Sun et al., 2016). These reductions, however, assume uncontrolled systems, where future observation statistics constitute an unbiased representation of the predictive state.222implicit reductions do exist in the system identification literature (van Overschee and de Moor, 1996) but they assume linear systems. Modeling controlled systems is more subtle since the the state of the system is a conditional distribution of observations given actions.

Another related work is the spectral learning algorithm for POMDPs proposed by Azizzadenesheli, Lazaric, and Anandkumar (2016). This method uses tensor factorization to recover POMDP parameters from examples collected by a non-blind memoryless policy. However, this method is limited to discrete POMDPs. Also, PSRs have more representational capacity than POMDPs and can compactly represent more sophisticated systems (Singh, James, and Rudary, 2004). There are other classes of dynamical system learning algorithms that are based on local optimization or sampling approaches (Fox et al., 2009; Frigola et al., 2013) but they do not offer consistency guarantees.

3 Formulation

We define a class of models that extends predictive state models of Hefny, Downey, and Gordon (2015) to controlled systems. We first introduce some notation: We denote by $\mathrm{Pr}[x\mid{\mathbf{do}}(Y=y)]$ the probability of $x$ given that we intervene by setting $Y$ to $y$. This is different from $\mathrm{Pr}[x\mid Y=y]$ which denotes conditioning on observing $Y=y$; in the former case, we ignore all effects on $Y$ by other variables. We denote by $V_{A\mid B;c}$ the linear operator that satisfies

[TABLE]

In other words for each $c$, $V_{A\mid B;c}$ is a conditional expectation operator from $B$ to $A$. In the discrete case, $V_{A\mid B;c}$ is just a conditional probability table.

When dealing with multiple variables, we will use tensor notation; e.g., $V_{A,B\mid C,D}$ is a 4-mode tensor. We will use

[TABLE]

to denote multiplying $V_{A,B\mid C,D}$ by $c$ along the mode corresponding to $C$ and by $d$ along the mode corresponding to $D$. If $c$ is a matrix then the multiplication is performed along the first dimension of $c$.

We will also use $\|\cdot\|_{F}$ to denote Frobenius norm, $a\otimes b$ to denote Kronecker product of two vectors and $A\star B$ to denote the Khatri-Rao product of two matrices (columnwise Kronecker product).

3.1 Model Definition

We will consider $k$-observable systems, where the posterior belief state given all previous observations and actions is uniquely identified by the conditional distribution $\mathrm{Pr}[o_{t:t+k-1}\mid{\mathbf{do}}(a_{t:t+k-1})]$.

Following Hefny, Downey, and Gordon (2015), we denote by ${\psi}^{o}_{t}$, ${\psi}^{a}_{t}$, ${\xi}^{o}_{t}$ and ${\xi}^{a}_{t}$ sufficient features of future observations $o_{t:t+k-1}$, future actions $a_{t:t+k-1}$, extended future observations $o_{t:t+k}$ and extended future actions $a_{t:t+k}$ at time $t$ respectively.

We also use $h^{\infty}_{t}\equiv o_{1:t-1},a_{1:t-1}$ to denote the entire history of observations and actions at time $t$ and use ${\psi}^{h}_{t}\equiv{\psi}^{h}(o_{1:t-1},a_{1:t-1})$ to denote finite features of previous observations and actions before time $t$.333Often but not always, ${\psi}^{h}_{t}$ is a computed from fixed-size window of previous observations and actions ending at $t-1$.

We are now ready to define the class of systems we are interested in.

Definition 1**.**

A dynamical system is said to conform to a predictive state controlled model (PSCM) if it satisfies the following properties:

• •

For each time $t$, there exists a linear operator $Q_{t}=V_{{\psi}^{o}_{t}\mid{\mathbf{do}}({\psi}^{a}_{t});h^{\infty}_{t}}$ (referred to as predictive state) such that $\mathbb{E}[{\psi}^{o}_{t}\mid{\mathbf{do}}(a_{t:t+k-1}),h^{\infty}_{t}]=Q_{t}{\psi}^{a}_{t}$

• •

For each time $t$, there exists a linear operator $P_{t}=V_{{\xi}^{o}_{t}\mid{\mathbf{do}}({\xi}^{a}_{t});h^{\infty}_{t}}$ (referred to as extended state) such that $\mathbb{E}[{\xi}^{o}_{t}\mid{\mathbf{do}}(a_{t:t+k}),h^{\infty}_{t}]=P_{t}{\xi}^{a}_{t}$

• •

There exists a linear map $W_{\mathrm{sys}}$ (referred to as system parameter map), such that, for each time $t$,

[TABLE]

• •

There exists a filtering function $f_{\mathrm{filter}}$ such that, for each time $t$, $Q_{t+1}=f_{\mathrm{filter}}(P_{t},o_{t},a_{t})$. $f_{\mathrm{filter}}$ is typically non-linear but known in advance.

It follows that a PSCM is specified by the tuple $(Q_{0},W_{\mathrm{sys}},f_{\mathrm{filter}})$, where $Q_{0}$ denotes the initial belief state.

There are a number of aspects of PSCMs that warrant discussion. First, unlike latent state models, the state $Q_{t}$ is represented by a conditional distribution of observed quantities. Second, $Q_{t}$ is a deterministic function of the history $h^{\infty}_{t}$. It represents the belief state that one should maintain after observing the history to make optimal predictions. Third, a PSCM specifies a recursive filter where given an action $a_{t}$ and an observation $o_{t}$, the state update equation is given by

[TABLE]

This construction allows us to have a linear map $W_{\mathrm{sys}}$ and still use it to build models with non-linear state updates, including IO-HMMs (Bengio and Frasconi, 1995), Kalman filters with inputs (van Overschee and de Moor, 1996) and HSE-PSRs (Boots, Gretton, and Gordon, 2013). As we see in Section 4, avoiding latent variables and having a linear $W_{\mathrm{sys}}$ enable the formulation of a consistent learning algorithm.

4 Learning A Predictive State Controlled Model

We assume that the extended features ${\xi}^{o}_{t}$ and ${\xi}^{a}_{t}$ are chosen such that $f_{\mathrm{filter}}$ is known. The parameters to learn are thus $W_{\mathrm{sys}}$ and $Q_{0}$. We also assume that a fixed blind (open-loop) policy is used to collect training data, and so we can treat causal conditioning on action ${\mathbf{do}}(a_{t})$ as ordinary conditioning on $a_{t}$.444One way to deal with non-blind training policies is to assign importance weights to training examples to correct the bias resulting from non-blindness (Bowling et al., 2006; Boots, Siddiqi, and Gordon, 2011). This, however, requires knowledge of the data collection policy and can result in a high variance of the estimated parameters. We defer the case of unknown non-blind policy to future work. It is possible, however, that a different (possibly non-blind) policy is used at test time.

To learn model parameters, we will adapt the two-stage regression method of Hefny, Downey, and Gordon (2015). Let $\bar{Q}_{t}\equiv\mathbb{E}[Q_{t}\mid{\psi}^{h}_{t}]$ (resp. $\bar{P}_{t}\equiv\mathbb{E}[P_{t}\mid{\psi}^{h}_{t}]$) be the expected state (resp. expected extended state) conditioned on finite history features ${\psi}^{h}_{t}$. For brevity, we might refer to $\bar{Q}_{t}$ simply as the (predictive) state when the distinction from $Q_{t}$ is clear. It follows from linearity of expectation that $\mathbb{E}[{\psi}^{o}_{t}\mid{\psi}^{a}_{t},{\psi}^{h}_{t}]=\bar{Q}_{t}{\psi}^{a}_{t}$ and $\mathbb{E}[{\xi}^{o}_{t}\mid{\xi}^{a}_{t},{\psi}^{h}_{t}]=\bar{P}_{t}{\xi}^{a}_{t}$; and it follows from the linearity of $W_{\mathrm{sys}}$ that

[TABLE]

So, we train regression models (referred to S1 regression models) to estimate $\bar{Q}_{t}$ and $\bar{P}_{t}$ from ${\psi}^{h}_{t}$. Then, we train another (S2) regression model to estimate $W_{\mathrm{sys}}$ from $\bar{Q}_{t}$ and $\bar{P}_{t}$. Being conditional distributions, estimating $\bar{Q}_{t}$ and $\bar{P}_{t}$ from ${\psi}^{h}_{t}$ is more subtle compared to uncontrolled systems, since we cannot use observation features as estimates of the state. We describe two methods to construct an S1 regression model to estimate $\bar{Q}_{t}$. The same methods apply to $\bar{P}_{t}$. As we show below, instances of both methods exist in the literature of system identification.

4.1 Joint S1 Approach

Let ${\psi}^{oa}_{t}$ denote a sufficient statistic of the joint observation/action distribution $\mathrm{Pr}({\psi}^{o}_{t},{\psi}^{a}_{t}\mid{\psi}^{h}_{t})$. This distribution is fixed for each value of ${\psi}^{h}_{t}$ since we assume a fixed model and policy. We use an S1 regression model to learn the map $f:{\psi}^{h}_{t}\mapsto\mathbb{E}[{\psi}^{ao}_{t}\mid{\psi}^{h}]$ by solving the optimization problem

[TABLE]

for some suitable Bregman divergence loss $l$ (e.g., square loss) and regularization $R$.

Once we learn $f$, we can estimate $\bar{Q}_{t}$ by first estimating the joint distribution $\mathrm{Pr}({\psi}^{o}_{t},{\psi}^{a}_{t}\mid{\psi}^{h}_{t})$ and then deriving the conditional operator $\bar{Q}_{t}$. By the continuous mapping theorem, a consistent estimator of $f$ results in a consistent estimator of $\bar{Q}_{t}$. An example of applying this method is using kernel Bayes rule (Fukumizu, Song, and Gretton, 2013) to estimate states in HSE-PSR (Boots, Gretton, and Gordon, 2013).

4.2 Conditional S1 Approach

In this method, instead of estimating the joint distribution represented by $\mathbb{E}[{\psi}^{oa}_{t}\mid{\psi}^{h}_{t}]$, we directly estimate the conditional distribution $\bar{Q}_{t}$. We exploit the fact that each training example ${\psi}^{o}_{t}$ is an unbiased estimate of $\bar{Q}_{t}{\psi}^{a}_{t}=\mathbb{E}[{\psi}^{o}_{t}\mid{\psi}^{a}_{t},{\psi}^{h}_{t}]$. We can formulate the S1 regression problem as learning a function $f:{\psi}^{h}_{t}\mapsto\bar{Q}_{t}$ that best matches the training examples, i.e., we solve the problem

[TABLE]

for some suitable Bregman divergence loss $l$ (e.g., square loss) and regularization $R$. An example of applying this method is the oblique projection method used in spectral system identification (van Overschee and de Moor, 1996). It is worth emphasizing that both the joint and conditional S1 approaches assume the state to be a conditional distribution. They only differ in the way to estimate that distribution.

4.3 S2 Regression and Learning Algorithm

Given S1 regression models to estimate $\bar{Q}_{t}$ and $\bar{P}_{t}$, learning a controlled dynamical system proceeds as shown in Algorithm 1.

4.4 Theoretical Guarantees

It is worth noting that Algorithm 1 is still an instance of the two stage regression framework described in (Hefny, Downey, and Gordon, 2015) and hence retains its theoretical guarantees: mainly that we can bound the error in estimating the dynamics matrix $W_{\mathrm{sys}}$ in terms of S1 regression error bounds, assuming that we collect examples from the stationary distribution of a blind policy with sufficient exploration.

A blind policy provides sufficient exploration if it has a stationary distribution that (1) visits a sufficient history set such that the set of equations $\mathbb{E}[P_{t}|{\psi}^{h}_{t}]=W_{sys}(\mathbb{E}[Q_{t}|{\psi}^{h}_{t}])$ are sufficient for estimating $W_{sys}$ and (2) provides training data to estimate $\mathbb{E}[Q_{t}|{\psi}^{h}_{t}]$ and $\mathbb{E}[P_{t}|{\psi}^{h}_{t}]$ with increasing accuracy.

Theorem 2**.**

Let $\pi$ be a blind data collection policy with a stationary distribution. If history, action and observation features have bounded norms, $\pi$ provides sufficient exploration, and ridge regression is used with $\lambda_{1}$ and $\lambda_{2}$ regularization parameter for S1 and S2 regression respectively, then for all valid states $Q$ the following is satisfied with probability at least $1-\delta$.

[TABLE]

where

[TABLE]

where $c>0$ is a problem-dependent constant.

We provide proofs and discussion of sufficient exploration condition in the supplementary material.

5 Predictive State Controlled Models With Random Fourier Features

Having a general framework for learning controlled dynamical systems, we now focus on HSE-PSR (Boots, Gretton, and Gordon, 2013) as a non-parametric instance of that framework using Hilbert space embedding of distributions. We first describe HSE-PSR learning as a two-stage regression method. Then we demonstrate how to obtain a finite dimensional approximation using random Fourier features (RFF) (Rahimi and Recht, 2008). Before describing HSE-PSR we give some necessary background on Hilbert space embedding and random Fourier features.

5.1 Hilbert Space Embedding of Distributions

We will briefly describe the concept of Hilbert space embedding of distributions. We refer the reader to (Smola et al., 2007) for more details on this topic. Hilbert space embedding of distributions provide a non-parametric generalizations of marginal, joint and conditional probability tables of discrete variables to continuous domains: namely, mean maps, covariance operators and conditional operators.

Let $k$ be a kernel associated with a feature map $\phi(x)$ such that $k(x_{1},x_{2})=\langle\phi(x_{1}),\phi(x_{2})\rangle$. A special case for discrete variables is the delta kernel where $\phi(x)$ maps $x$ to an indicator vector. For a random variable $X$, the mean map $\mu_{X}$ is defined as $\mathbb{E}[\phi_{\cal X}(X)]$. Note that $\mu_{X}$ is an element of the reproducing kernel Hilbert space (RKHS) associated with $k$.

The uncentered covariance operator of two variables $X$ and $Y$ is $C_{XY}=\mathbb{E}[\phi_{\cal X}(X)\otimes\phi_{\cal Y}(Y)]$. For universal kernels $k_{\cal X}$ and $k_{\cal Y}$, $C_{XY}$ is a sufficient representation of the joint distribution $\mathrm{Pr}(X,Y)$. In this paper, we will use $C_{XY\mid z}$ to denote the covariance of $X$ and $Y$ given that $Z=z$.

Under smoothness assumptions, (Song et al., 2009) show that $V_{\phi_{\cal X}(X)\mid\phi_{\cal Y}(Y)}=C_{XY}C_{XX}^{-1},$ where the conditional operator $V$ is as defined in Section 3. More generally, $V_{\phi_{\cal X}(X)\mid\phi_{\cal Y}(Y);z}=C_{XY\mid z}C_{XX\mid z}^{-1}$.

5.2 HSE-PSR as a predictive state controlled model

HSE-PSR is a generalization of IO-HMM that has proven to be successful in practice (Boots, Gretton, and Gordon, 2013; Boots and Fox, 2013). It is suitable for high dimensional and continuous observations and/or actions. HSE-PSR uses kernel feature maps as sufficient statistics of observations and actions. We define four kernels $k_{O},k_{A},k_{o},k_{a}$ over future observation features, future action features, individual observations and individual actions respectively.

We can then define ${\psi}^{o}_{t}=\phi_{O}(o_{t:t+k-1})$ and similarly ${\psi}^{a}_{t}=\phi_{A}(a_{t:t+k-1})$. We will also use $\phi^{o}_{t}$ and $\phi^{a}_{t}$ as shorthands for $\phi_{o}(o_{t})$ and $\phi_{a}(a_{t})$. The extended future is then defined as ${\xi}^{o}_{t}={\psi}^{o}_{t}\otimes\phi^{o}_{t}$ and ${\xi}^{a}_{t}={\psi}^{a}_{t}\otimes\phi^{a}_{t}$

Under the assumption of a blind learning policy, the operators $Q_{t}$ and $P_{t}$ are defined to be

[TABLE]

Therefore, $Q_{t}$ specifies the state of the system as a conditional distribution of future observations given future actions while $P_{t}$ is a tuple of two operators that allow us to condition on the pair $(a_{t},o_{t})$ to obtain $Q_{t+1}$. In more detail, filtering in an HSE-PSR is carried out as follows

• •

From $o_{t}$ and $a_{t}$, obtain $\phi^{o}_{t}$ and $\phi^{a}_{t}$.

• •

Compute $C_{o_{t}o_{t}\mid h^{\infty}_{t},a_{t}}=V_{\phi^{o}_{t}\otimes\phi^{o}_{t}\mid\phi^{a}_{t};h^{\infty}_{t}}\phi^{a}_{t}$

• •

Multiply by inverse observation covariance to change “predicting $\phi^{o}_{t}$” into “conditioning on $\phi^{o}_{t}$”:

[TABLE]

• •

Condition on $\phi^{o}_{t}$ and $\phi^{a}_{t}$ to obtain shifted state

[TABLE]

Thus, in HSE-PSR, the parameter $W_{\mathrm{sys}}$ is composed of two linear maps; $f_{o}$ and $f_{\xi}$ such that $P_{t}^{\xi}=f_{\xi}(Q_{t})$ and $P_{t}^{o}=f_{o}(Q_{t})$. In the following section we show how to estimate $\bar{Q}_{t}$ and $\bar{P}_{t}$ from data. Estimation of $f_{\xi}$, $f_{o}$ can then be carried out using kernel regression.

Learning and filtering in an HSE-PSR can be implicitly carried out in the RKHS using a Gram matrix formulation. We will describe learning in terms of the RKHS elements and refer the reader to (Boots, Gretton, and Gordon, 2013) for details on the Gram matrix formulation. As we mention in Section 5, random Fourier features, provides a scalable approximation to operating in the RKHS.

5.3 S1 Regression for HSE-PSR

As discussed in section 4 we can use a joint or conditional approach for S1 regression. We now demonstrate how these two approaches apply to HSE-PSR.

Joint S1 Regression for HSE-PSR

This is the method used in (Boots, Gretton, and Gordon, 2013). In this approach we exploit the fact that

[TABLE]

So, we learn two linear maps $T_{oa}$ and $T_{a}$ such that $T_{oa}({\psi}^{h}_{t})\approx C_{{\psi}^{o}_{t}{\psi}^{a}_{t}\mid{\psi}^{h}_{t}}$ and $T_{a}({\psi}^{h}_{t})\approx C_{{\psi}^{a}_{t}{\psi}^{a}_{t}\mid{\psi}^{h}_{t}}$. The training examples for $T_{oa}$ and $T_{a}$ consist of pairs $({\psi}^{h}_{t},{\psi}^{o}_{t}\otimes{\psi}^{a}_{t})$ and $({\psi}^{h}_{t},{\psi}^{a}_{t}\otimes{\psi}^{a}_{t})$ respectively.

Once we learn this map, we can estimate $C_{{\psi}^{o}_{t}{\psi}^{a}_{t}|{\psi}^{h}_{t}}$ and $C_{{\psi}^{a}_{t}{\psi}^{a}_{t}|{\psi}^{h}_{t}}$ and consequently estimate $\bar{Q}_{t}$.

Conditional S1 Regression for HSE-PSR

It is also possible to apply the conditional S1 regression formulation in Section 4.2. Specifically, let $\cal{F}$ be the set of 3-mode tensors, with modes corresponding to ${\psi}^{o}_{t}$, ${\psi}^{o}_{t}$ and ${\psi}^{h}_{t}$. We estimate a tensor $T^{*}$ by optimizing

[TABLE]

where $\|.\|_{HS}^{2}$ is the Hilbert-Schmidt norm, which translates to Frobenius norm in finite-dimensional Euclidan spaces. We can then use

[TABLE]

For both regression approaches, the same procedure can be used to estimate the extended state $\bar{P}_{t}$ by replacing features ${\psi}^{o}_{t}$ and ${\psi}^{a}_{t}$ with their extended counterparts ${\xi}^{o}_{t}$ and ${\xi}^{a}_{t}$.

5.4 Approximating HSE-PSR with Random Fourier Features

A Gram matrix formulation of the HSE-PSR has computational and memory requirements that grow rapidly with the number of training examples. To alleviate this problem, we resort to kernel approximation—that is, we replace RKHS vectors such as $\psi_{t}^{o}$ and $\psi_{t}^{a}$ with finite dimensional vectors that approximately preserve inner products. We use random Fourier features (RFF) (Rahimi and Recht, 2008) as an approximation but it is possible to use other approximation methods. Unfortunately RFF approximation can typically require $D$ to be prohibitively large. Therefore, we apply principal component analysis (PCA) to the feature maps to reduce their dimension to $p\ll D$. We apply PCA again to quantities that require $p^{2}$ space such as extended features ${\xi}^{o}_{t}$, ${\xi}^{a}_{t}$ and states $\bar{Q}_{t}$, reducing them to $p$ dimensions. We map them back to $p^{2}$ dimensions when needed (e.g., for filtering). We also employ randomized SVD (Halko, Martinsson, and Tropp, 2011) for fast computation of PCA, resulting in an algorithm that scales linearly with $N$ and $D$. 44footnotetext: We provide pseudo-code in the supplementary material. MATLAB source code is available at: https://github.com/ahefnycmu/rffpsr

5.5 Model refinement by local optimization

A common practice is to use the output of a moment-based algorithm to initialize a non-convex optimization algorithm such as EM (Belanger and Kakade, 2015) or gradient descent (Jiang, Kulesza, and Singh, 2016). Since EM is not directly applicable to RFF-PSRs, we propose a gradient descent approach. We can observe that filtering in an RFF-PSR defines a recurrent structure given by.

[TABLE]

where $W_{\mathrm{pred}}$ is a linear operator that predicts the next observation.555The linearity of $W_{\mathrm{pred}}$ is a valid assumption for a universal kernel. If $f_{\mathrm{filter}}$ is differentiable, we can improve our estimates of $W_{\mathrm{sys}}$ and $W_{\mathrm{pred}}$ using backpropagation through time (BPTT) (Werbos, 1990). In particular, we optimize the error in predicting (features of) a window of observations. In our experiments, we learn to predict $o_{t:t+k-1}$ given $a_{t:t+k-1}$.

6 Experiments

6.1 Synthetic Data

We use the benchmark synthetic non-linear system used by (Boots, Gretton, and Gordon, 2013) :

[TABLE]

The input $a$ is generated as zero-order hold white noise, uniformly distributed between $-0.5$ and 0.5. We collected 20 trajectories of 100 observations and actions at 20Hz and we split them into 10 training, 5 validation and 5 test trajectories. The prediction target for this experiment is $o(t)$.

6.2 Predicting windshield view

In this experiment we used the TORCS car simulation server, which outputs 64x64 images (see Figure 2). The observations are produced by converting the images to greyscale and projecting them to 200 dimensions via PCA. The car is controlled by a built-in controller that controls acceleration while the external actions control steering. We collected 50 trajectories by applying a sine wave with random starting phase to the steering control and letting the simulator run until the car goes off the track. We used 40 trajectories for training, 5 for validation and 5 for testing. The prediction target is the projected image.

6.3 Predicting the nose position of a simulated swimmer robot

We consider the 3-link simulated swimmer robot from the open-source package RLPy (Geramifard et al., 2013). The 2-d action consists of torques applied on the two joints of the links. The observation model returns the angles of the joints and the position of the nose (in body coordinates). The measurements are contaminated with Gaussian noise whose standard deviation is 5% of the true signal standard deviation. To collect the data, we use an open-loop policy that selects actions uniformly at random. We collected 25 trajectories of length 100 each and use 24 for training and 1 for validation. We generate test trajectories using a mixed policy: with probability $p_{\mathrm{blind}}$, we sample a uniformly random action, while with probability $1-p_{\mathrm{blind}}$, we sample an action from a pre-specified deterministic policy that seeks a goal point. We generate two sets of 10 test trajectories each, one with $p_{\mathrm{blind}}=0.8$ and another with $p_{\mathrm{blind}}=0.2$. The prediction target is the position of the nose.

6.4 Tested Methods and Evaluation Procedure

We tested three different initializations of RFF-PSR (with Gaussian RBF kernel): random initialization, two-stage regression with joint S1, and two-stage regression with conditional S1 (Section 5.3). For each initialization, we tested the model before and after refinement. For refinement we used BPTT with a decreasing step size: the step size is reduced by half if validation error increases. Early stopping occurs if the step size becomes too small ($10^{-5}$) or the relative change in validation is insignificant ($10^{-3}$). We also test the following baselines.

HSE-PSR: We implemented the Gram matrix HSE-PSR as described in (Boots, Gretton, and Gordon, 2013).

N4SID: We used MATLAB’s implementation of subspace identification of linear dynamical systems.

Non-linear Auto Regression (RFF-ARX): We implemented a version of auto regression where the predictor variable is the RFF representation of future actions together with a finite history of previous observations and actions, and the target variable is future observations.

Models were trained with future length of 10 and history length of 20. For RFF-PSR and RFF-ARX we used 10000 random features and applied PCA to project features onto 20 dimensions. Kernel bandwidths were set to the median of the distance between training points (median trick). For evaluation, we perform filtering on the data and estimate the prediction target of the experiment at test time $t$ given the history $o_{1:t-H},a_{1:t}$, where $H$ is the prediction horizon. We report the mean square error across all times $t$ for each value of $H\in\{1,2,\dots,10\}$.

6.5 Results and Discussion

The results are shown in Figure 1. There are a number of important observations.

• •

In general, joint S1 training closely matches or outperforms conditional S1 training, with and without refinement.

• •

Local refinement significantly improves predictive performance for all initialization methods.

• •

Local refinement, on its own, is not sufficient to produce a good model. The two stage regression provides a good initialization of the refinement procedure.

• •

Even without refinement, RFF-PSR outperforms HSE-PSR. This could be attributed to the dimensionality reduction step, which adds appropriate inductive bias.

• •

Compared to other methods, RFF-PSR has better performance with non-blind test policies.

7 Conclusion

We proposed a framework to learn controlled dynamical systems using two-stage regression. We then applied this framework to develop a scalable method for controlled non-linear system identification: using RFF approximation of HSE-PSR together with a refinement procedure to enhance the model after a two-stage regression initialization. We have demonstrated promising results for the proposed method in terms of predictive performance. As future work, we would like to use this framework for further tasks such as imitation learning and reinforcement learning.

Acknowledgements

The authors gratefully acknowledge support from ONR (grant number N000141512365), DARPA (grant number FA87501720152) and NSF EAGER (grant number IIS1450543). The authors would like to thank Wen Sun and Yuxiang Wang for the helpful discussions.

Appendix A RFF-PSR Learning Algorithm

For ease of exposition, we assume that RFF features are computed prior to PCA. In our implementation, we compute the RFF features on the fly while performing PCA to reduce the required memory footprint. Here we use $A\star B$ to denote the Khatri-Rao product of two matrices (columnwise Kronecker product).

Appendix B Examples of Predictive State Controlled Models

Here we discuss IO-HMM and Kalman filter with inputs, showing that they are instances of PSCMs. We do this for each model by defining the predictive state, showing that it satisfies the condition $P_{t}=WQ_{t}$ and describing an S1 regression method.

B.1 IO-HMM

Let $T$ be the transition tensor such that $T\times_{s}s_{t}\times_{a}a_{t}=\mathbb{E}[s_{t+1}|a_{t},s_{t}]$ and $O$ be the observation tensor such that $O\times_{s}s_{t}\times_{a}a_{t}=\mathbb{E}[o_{t}|a_{t},s_{t}]$.

Define $O^{k}$ to be the extended observation tensor where $O_{k}\times_{s}s_{t}\times_{a}a_{t:t+k-1}=\mathbb{E}[o_{t:t+k-1}|a_{t:t+k-1},s_{t}]$

As a shortcut, we will denote by $T_{ij}$ the product $T\times_{s}e_{i}\times_{a}e_{j}$.

For $k=1$, we have $O^{1}=O$.

For $k>1$ we can think of $a_{t:t+k-1}$ as the outer product $a_{t}\otimes a_{t+1:t+k}$. So we can define $O^{k}$ such that

[TABLE]

In words, starting from state $e_{i}$ and applying an action $e_{j}$ followed by a sequence of $k-1$ actions denoted by indicator $e_{l}$. The expected indicator of the next $k$ observations is the outer product of expected observation $o_{t}$ (given by $O_{ij}$) with the expected indicator of observations $o_{t+1:t+k-1}$ as predicted by $O^{k-1}$. Note that the two expectations being multiplied are conditionally independent given the state $e_{i}$ and the action sequence.

Given the tensor $O^{k}$ the predictive states $Q_{t}$ and $P_{t}$ are defined to be

[TABLE]

Now to show that (1) holds, let $\tilde{O}^{k}$ be a reshaping of $O^{k}$ into a matrix such that

[TABLE]

It follows that

[TABLE]

which is linear in $Q_{t}$.

S1 Regression

Let $s_{t}=s(h^{\infty}_{t})$ be the belief state at time $t$. Note that $s_{t}$ is a deterministic function of the entire history.

Under a fixed policy assumption, an indicator vector of the joint observation and action assignment is an unbiased estimate of the joint probability table $\mathbb{P}[{\psi}^{a}_{t},{\xi}^{a}_{t}\mid h^{\infty}_{t}]$. An S1 regression model can be used to learn the mapping ${\psi}^{h}_{t}\mapsto\mathbb{P}[{\psi}^{a}_{t},{\xi}^{a}_{t}\mid{\psi}^{h}_{t}]$. It is then easy to estimate the conditional probability table $\bar{Q}_{t}$ from the joint probability table $\mathbb{P}[{\psi}^{a}_{t},{\xi}^{a}_{t}\mid{\psi}^{h}_{t}]$.

We can also use the conditional S1 approach. By exploiting the fact that ${\psi}^{o}_{t}$ is an unbiased estimate of a single column of $Q_{t}$ corresponding to ${\psi}^{a}_{t}$. We can use (3) to learn a function $f:h_{t}\mapsto\bar{Q}_{t}$ that best matches the training examples.

B.2 Kalman Filter with inputs

The Kalman filter is given by

[TABLE]

Given a belief state $s_{t}\equiv\mathbb{E}[x_{t-1}|h^{\infty}_{t}]$ we can write the predictive state as

[TABLE]

where

[TABLE]

The extended predictive state have similar form with $\Gamma_{k}$ and $U_{k}$ replaced with $\Gamma_{k+1}$ and $U_{k+1}$. Since $U$ is fixed, keeping track of the state amounts to keeping track of $Q_{t}\equiv\Gamma_{k}s_{t}$. It follows that

[TABLE]

If $h_{t}$ is a linear projection of $h^{\infty}_{t}$ (e.g. stacking of a finite window of observations and actions), it can also be shown van Overschee and de Moor (1996) that

[TABLE]

for some matrix $\tilde{\Gamma}_{k}$.

S1 Regression

Let ${\cal F}$ be the set of functions that take the form

[TABLE]

The oblique projection method van Overschee and de Moor (1996) uses linear regression to estimate $\Gamma$ and $B$ (essentially solving (3)). Having a fixed $B$, the conditional operator is determined by $\Gamma h_{t}$ through an affine transformation. Therefore we can use $\bar{Q}_{t}=\Gamma h_{t}$.

Appendix C Theoretical Analysis

Let ${\cal H}=\{h_{i}\}_{i=1}^{N}$ be a set of histories generated from an i.i.d distribution. 666The i.i.d property is achieved if we can restart the system or if the data collection policy induces an ergodic process with a stationary distribution. In the latter case, we assume the examples are sufficiently spaced in time to that allow the process to mix. However, in practice, we use all examples as this makes the error only smaller.

We use $Q({\psi}^{h})$ to denote $\mathbb{E}[Q|{\psi}^{h}]$.

The main theorem in Hefny, Downey, and Gordon (2015) bounds parameter estimation error in terms of S1 regression error. This implies that we need to analyze the properties of S1 regression to prove Theorem 2. We will look at multiple scenarios where in each scenario we develop sufficient exploration conditions and provide an S1 error bound for these conditions.

Definition C.1** (Sufficient history set).**

Consider a PSCM that satisfies

[TABLE]

A set of histories ${\cal H}=\{h_{i}\}_{i=1}^{M}$ is called a sufficient history set if it is sufficient to estimate $W_{\mathrm{sys}}$ using $\mathbb{E}[{Q}_{t}|{\psi}^{h}_{t}=h]$ and $\mathbb{E}[{P}_{t}|{\psi}^{h}_{t}=h]$ for each $h\in{\cal H}$.

Note that $W_{\mathrm{sys}}$ may not be unique, we care about estimating $W_{\mathrm{sys}}Q$ for any valid $Q$. From the above definition, it follows that a data collection policy provides sufficient exploration if it allows for estimating $\mathbb{E}[{Q}|{\psi}^{h}_{t}=h]$ and $\mathbb{E}[{P}|{\psi}^{h}_{t}=h]$ for a sufficient history set with increasing accuracy.

C.1 Case 1: Discrete Observations and Actions

Consider a discrete system where ${\cal H}$, ${\cal A}$, ${\cal A^{+}}$, ${\cal O}$, ${\cal O^{+}}$ are the set of all possible histories, future action sequences, extended future action sequences, future observation sequences and extended future observation sequences respectively.

Theorem C.2**.**

Assume a discrete system where the data collection policy induces an i.i.d distribution over histories. If the policy generates each possible extended future action sequence starting from each possible history $M$ times, then it generates an S2 training dataset of size $N=M|{\cal H}||{\cal A^{+}}|$ with S1 error bound $\eta_{\delta,N}=\sqrt{\frac{|{\cal H}||{\cal A^{+}}||{\cal O^{+}}|}{2M}\log\left(\frac{2|{\cal H}||{\cal A^{+}}||{\cal O^{+}}|}{\delta}\right)}$

Proof.

The proof follows immediately from Heoffding’s inequality which bounds the error in estimating the probability of an event by averaging.

Note that we need to estimate $|{\cal H}||{\cal A}||{\cal O}|$ probabilities to estimate $Q$ and $|{\cal H}||{\cal A^{+}}||{\cal O^{+}}|$ probabilities to estimate $P$. Therefore we divide $\delta$ by $2|{\cal H}||{\cal A^{+}}||{\cal O^{+}}|$ to correct for multiple probability estimates. ∎

Remark C.3**.**

Assume the system to be 1-observable, where the history and future are of length 1. Then a consistent estimate of $Q$ and $P$ can be obtained by a consistent estimate of the joint probability table $P(o_{t-1:t+1},a_{t-1:t+1})$.

C.2 Case 2: Continuous System

Definition C.4** (Range and span of a policy).**

Let $\pi$ be a data collection policy with a stationary distribution. For a random vector $X_{t}=f(h^{\infty}_{t},o_{t:\infty},a_{t:\infty})$, the range of $\pi$ on $X$ is the support of the stationary distribution of $X_{t}$ induced by the policy $\pi$ (i.e. the set of all possible values of $X_{t}$ that can be generated by the stationary distribution).

The span of $\pi$ on $X$ is the subspace spanned by the range of $\pi$ on $X$.

When referring to the policy range or span, we may omit the variable name when it is clear in the context.

Condition C.5** (Action span for joint S1).**

Let $\pi$ be data collection policy and let $\cal H$ be the range of $\pi$ on histories. The action span condition for joint S1 is defined as the requirement to satisfy the following:

${\cal H}$* is a sufficient history set.* 2. 2.

For any ${\psi}^{h}\in{\cal H}$, the conditional covariance $\Sigma_{{\psi}^{a}|{\psi}^{h}}$ is full rank.

Condition C.6** (Action span for conditional S1).**

Let $\pi$ be data collection policy and let $\cal H$ be the range of $\pi$ on histories. The action span condition for conditional S1 is defined as the requirement to satisfy the following:

${\cal H}$* is a sufficient history set.* 2. 2.

For any $\psi^{h}\in{\cal H}$ and any future action feature vector $\psi^{a}$, the quantity $(\psi^{h}\otimes\psi^{a})$ is in the policy span.

Remark C.7**.**

Condition C.5 implies Condition C.6.

Assumption C.8** (Bounded features).**

We assume that $\|\psi^{h}\|<c_{h}$ for all $h\in{\cal H}$. Also, we assume that $\|{\psi}^{o}\|\leq c_{O}$ and $\|{\psi}^{a}\|\leq c_{A}$ for any valid future observation sequence and action sequence respectively.

Theorem C.9**.**

Let $\pi$ be a data collection policy and let $\cal H$ be the range of $\pi$ on histories. If Assumption C.8 and Condition C.6 are satisfied and conditional S1 regression is used with a liner model as the correct model, then $\pi$ provides sufficient exploration and, for all $h\in{\cal H}$ and any $\delta\in(0,1)$ such that $N>\frac{c^{2}\log(2d_{h}d_{A}/\delta)}{\lambda_{\min}(\Sigma_{{\psi}^{h}\otimes{\psi}^{a}})}$, the following holds with probability at least $1-\delta$

[TABLE]

where

[TABLE]

In the following section we provide a proof sketch for the asymptotic form in Theorem 2 for joint S1.

Remark C.10** (Conditioning).**

It is known that linear regression converges faster if the problem is well-conditioned. In the two stage regression we need the good conditioning of both stages– that is,

• •

The set of training histories result in a problem $\bar{P}_{t}=W\bar{Q}_{t}$ that is well conditioned (S2 conditioning).

• •

The S1 regression problem is well conditioned.

The second requirement ensures that we converge fast to good estimates of $\bar{Q}_{t}$ and $\bar{P}_{t}$. Designing exploration policies that result in well conditioned two stage regression problems is an interesting direction for future work.

Appendix D Proofs of theorems

In this section we provide proofs for Theorem C.9. The asymptotic statement in Theorem 2 follows directly from the main theorem in (Hefny, Downey, and Gordon, 2015). We also provide a proof sketch for the joint S1 case.

The proof strategy is as follows: First, we use matrix concentration bounds to analyze the effect of using estimated covariance matrices. Then, we analyze the effect of error in covariance matrix on regression weights. By combining the results of both analyses, we prove the desired theorems.

Lemma D.1** (Matrix Chernoff Inequality (Tropp, 2015)).**

Consider a finite sequence $\{S_{k}\}$ of independent, random, Hermitian matrices with common dimension $d$. Assume that

[TABLE]

Introduce the random matrix

[TABLE]

Define

[TABLE]

Then, for any $\epsilon\in[0,1)$

[TABLE]

Corollary D.2** (Minimum eigenvalue of empirical covariance).**

Let $X$ be a random variable of dimensionality $d$ such that $\|X\|<c$. Let $\{x_{k}\}_{k=1}^{N}$ be $N$ i.i.d samples of the distribution of $X$.

Define

[TABLE]

For any $\delta\in(0,1)$ such that $N>\frac{c^{2}\log(2d/\delta)}{\lambda_{\min}(\Sigma_{X})}$ the following holds with probability at least $1-\delta$

[TABLE]

Proof.

Define $S_{k}=\frac{1}{N}x_{k}x_{k}^{\top}$. It follows that $\lambda_{\max}(S_{k})\leq L=c^{2}/N$ and $\mu_{\min}=\lambda_{\min}(\Sigma_{X})$. Define

[TABLE]

which implies that

[TABLE]

It follows from Lemma D.1 that $\mathrm{Pr}(\lambda_{\min}(\hat{\Sigma}_{X})\leq(1-\epsilon)\mu_{\min})\leq\delta$

∎

Lemma D.3** (Matrix Bernstein Inequality (Tropp, 2015)).**

Consider a finite sequence $\{S_{k}\}$ of independent, random matrices with common dimensions $d_{1}\times d_{2}$. Assume that

[TABLE]

Introduce the random matrix

[TABLE]

Let $v(Z)$ be the matrix variance statistic of the sum:

[TABLE]

Then

[TABLE]

Corollary D.4** (Error in empirical cross-covariance).**

With probability at least $1-\delta$

[TABLE]

where

[TABLE]

Proof.

Define $S_{k}=y_{k}x_{k}^{\top}-\Sigma_{YX}$, it follows that

[TABLE]

Applying Lemma D.3 we get

[TABLE]

and hence

[TABLE]

This quadratic inequality implies

[TABLE]

Using the fact that $\sqrt{a^{2}+b^{2}}\leq|a|+|b|$ we get

[TABLE]

∎

Corollary D.5** (Normalized error in empirical covariance).**

With probability at least $1-\delta$

[TABLE]

where

[TABLE]

Proof.

Define $S_{k}=\Sigma_{X}^{-1/2}x_{k}x_{k}^{\top}-\Sigma_{X}^{1/2}$, it follows that

[TABLE]

Applying Lemma D.3 we get

[TABLE]

and similar to the proof of Corollary D.4, we can show that

[TABLE]

∎

Lemma D.6**.**

Let $\hat{\Sigma}_{YX}=\Sigma_{YX}+\Delta_{YX}$ and $\hat{\Sigma}_{X}=\Sigma_{X}+\Delta_{X}$ where $\mathbb{E}[\Delta_{YX}]$ and $\mathbb{E}[\Delta_{YX}]$ are not necessarily zero and $\hat{\Sigma}_{X}$ is symmetric positive semidefinite. Define $W=\Sigma_{YX}\Sigma_{X}^{-1}$ and $\hat{W}=\hat{\Sigma}_{YX}(\hat{\Sigma}_{X}+\lambda)^{-1}$. It follows that

[TABLE]

Proof.

[TABLE]

It follows that

[TABLE]

As for $T_{1}$, using the matrix inverse Lemma $B^{-1}-A^{-1}=B^{-1}(A-B)A^{-1}$ and the fact that $\Sigma_{YX}=\Sigma_{Y}^{1/2}V\Sigma_{X}^{1/2}$, where $V$ is a correlation matrix satisfying $\|V\|\leq 1$ we get

[TABLE]

and hence

[TABLE]

∎

Corollary D.7**.**

Let ${x_{k}}_{k=1}^{N}$ and ${y_{k}}_{k=1}^{N}$ be i.i.d samples from two random variables $X$ and $Y$ with dimensions $d_{X}$ and $d_{Y}$ and (uncentered) covariances $\Sigma_{X}$ and $\Sigma_{Y}$ respectively. Assume $\|X\|\leq c_{x}$ and $\|Y\|\leq c_{y}$. Let $\hat{\Sigma}_{YX}=\frac{1}{N}\sum_{k=1}^{N}y_{k}x_{k}^{\top}$ and $\hat{\Sigma}_{X}=\frac{1}{N}\sum_{k=1}^{N}x_{k}x_{k}^{\top}$. Define $W=\Sigma_{YX}\Sigma_{X}^{-1}$ and $\hat{W}=\hat{\Sigma}_{YX}(\hat{\Sigma}_{X}+\lambda)^{-1}$.

For any $\delta\in(0,1)$ such that $N>\frac{c_{x}^{2}\log(2d_{X}/\delta)}{\lambda_{\min}(\Sigma_{X})}$ the following holds with probability at least $1-3\delta$

[TABLE]

where

[TABLE]

Proof.

This corollary follows simply from applying Corollaries D.2, D.4 and D.5 to Lemma D.6. The $1-3\delta$ bound follows from union bound; since we have three probabilitic bounds each of which holds with probability $1-\delta$. ∎

Lemma D.8**.**

Let $\hat{\Sigma}_{YX}=\Sigma_{YX}+\Delta_{YX}$ and $\hat{\Sigma}_{X}=\Sigma_{X}+\Delta_{X}$ where $\mathbb{E}[\Delta_{YX}]$ and $\mathbb{E}[\Delta_{YX}]$ is not necessarily zero and $\hat{\Sigma}_{X}$ is symmetric but not necessarily positive semidefinite. Define $W=\Sigma_{YX}\Sigma_{X}^{-1}$ and $\hat{W}=\hat{\Sigma}_{YX}\hat{\Sigma}_{X}(\hat{\Sigma}_{X}^{2}+\lambda)^{-1}$. It follows that

[TABLE]

Proof.

[TABLE]

Using the matrix inverse Lemma $B^{-1}-A^{-1}=B^{-1}(A-B)A^{-1}$ and the fact that $\Sigma_{YX}=\Sigma_{Y}^{1/2}V\Sigma_{X}^{1/2}$, where $V$ is a correlation matrix satisfying $\|V\|\leq 1$ we get

[TABLE]

∎

D.1 Proof of Theorem C.9

Proof.

In the linear case, we estimate a tensor $T$ with modes corresponding to ${\psi}^{h}$, ${\psi}^{a}$ and ${\psi}^{o}$ by solving the minimization problem in Section 5.3. Equivalently, we estimate a matrix $T_{r}$ of size $d_{O}\times d_{h}d_{A}$ where an input ${\psi}^{h}\otimes\psi^{a}$ is mapped to an output $\mathbb{E}[\psi^{o}\mid h,\psi^{a}]$. Note that

[TABLE]

For any history $h\in{\cal H}$ and future action feature vector ${\psi}^{a}$ we have

[TABLE]

Note that Condition C.6 implies that ${\psi}^{h}\otimes{\psi}^{a}$ will eventually be in the span of training examples. This rules out the case where the inequality is satisfied only because $({\psi}^{h}\otimes{\psi}^{a})$ is incorrectly in the null space of $\hat{T}_{r}$ and $T_{r}$.

The theorem is proven by applying Corollary D.7 to bound $\|\hat{T}_{r}-T_{r}\|$. ∎

D.2 Sketch Proof for Joint S1

Let $T_{A}$ be a tensor such that $\Sigma_{{\psi}^{a}\mid{\psi}^{h}}=T_{A}\times_{h}{\psi}^{h}$ In order to prove Theorem 2 for joint S1, note that

[TABLE]

From Lemma D.6, we obtain a high probability bound on $\|\hat{T}_{A}-T_{A}\|$ and $\|\hat{T}_{OA}-T_{OA}\|$. Then we apply these bounds to Lemma D.8 to obtain an error in $Q({\psi}^{h})$.