An Expectation-Maximization Algorithm to Compute a Stochastic Factorization From Data

1 Introduction

Results from the supervised learning literature indicate that general efficient learning is only possible if one exploits some kind of regularity in the problem [Gy¨orfi et al., 2002]. When learning a model that describes a stochastic process, one can exploit, for example: a sparse or low-rank transition matrix, a smooth transition kernel, or the fact that the data lies close to a lower dimensional manifold [Farahmand, 2011]. Such regularities allow estimating parameters with a reasonable amount of data [Gy¨orfi et al., 2002].

Joelle Pineau and Doina Precup

McGill University Montreal, QC, Canada {jpineau,dprecup}@cs.mcgill.ca

One type of structural regularity that has recently received attention is stochastic factorization [Barreto and Fragoso, 2011]. In this case, a transition probability matrix of size n×n is represented as the product of two stochastic matrices, of sizes n×m and m×n, where m < n. Interestingly, swapping the factors of the multiplication yields another stochastic matrix, of size m×m, which retains some fundamental characteristics of its precursor [Barreto and Fragoso, 2011]. This makes it possible to use the new matrix instead of the original-which can lead to significant savings in terms of computational effort if m is considerably smaller than n. This strategy, dubbed the "stochastic-factorization trick," has been used to compute the stationary distribution of a Markov chain [Barreto and Fragoso, 2011], to determine the fundamental matrix of an absorbing chain [Barreto and Fragoso, 2011], and to compute decision policies in large problems through dynamic programming [Barreto et al., 2014] and reinforcement learning [Barreto et al., 2011].

So far, the stochastic factorization problem has appeared in the literature with a strict "algebraic" formulation: given a transition matrix, find two lower-order stochastic matrices whose multiplication approximates the original matrix as well as possible [Barreto and Fragoso, 2011; Barreto et al., 2014]. In this paper we are interested in computing a stochastic factorization when the transition matrix is not known, and we only have access to transitions sampled from it. One possible approach in this case is to use the samples to compute an estimate of the transition matrix, which takes us back to the scenario described above. However, this approach has the disadvantage that many entries may need to be estimated, perhaps from not much data. Instead, we would like to estimate the factors directly from data, because this entails potentially many fewer parameters. We introduce a probabilistic interpretation of a stochastic factorization and build on the resulting model to derive an expectation-maximization algorithm that does just that.

Since the order of the factorization (i.e., the value of m) is a parameter of the proposed algorithm, by changing it one can control the complexity of the resulting model, and thus the trade off between the approximation and estimation errors [Gy¨orfi et al., 2002]. This means that, if the transition matrix can be well approximated by a low-order stochastic factorization, the stochastic-factorization trick may provide benefits not only from a computational point of view but also in terms of the number of sample transitions needed in the approximation. We empirically show that, when compared to estimating the transition matrix directly via maximum likelihood (that is, by "counting" events), the proposed method is able to compute approximations of roughly the same quality using less data. In order to illustrate the practical utility of our algorithm, we use it to learn a model of the game of blackjack.

2 Background and Notation

This section introduces the notation adopted and briefly reviews some concepts that will be used throughout the paper.

Random variables are represented by capital letters; we use the notation A1:τ to refer to a sequence of variables A1,A2, ...,Aτ . Scalar variables are represented by small letters; boldface capital letters and boldface small letters are used to denote matrices and vectors, respectively. We will need the following definitions:

Definition 1 A matrix P ∈ Rn×q is called stochastic if and only if pij ≥ 0 for all i, j and ∑qj=1 pij = 1 for all i. A square stochastic matrix is called a transition matrix.

Definition 2 Given a stochastic matrix P ∈ Rn×q , the relation P = DK is called a stochastic factorization of P if D ∈ Rn×m and K ∈ Rm×q are also stochastic matrices. The integer m > 0 is the order of the factorization.

Definition 3 The stochastic rank of a stochastic matrix P ∈ Rn×q , denoted by srk(P), is the smallest possible order of the stochastic factorization P = DK.

Given a stochastic factorization of a transition matrix, P = DK, swapping the factors of the factorization yields another transition matrix ¯P = KD, potentially much smaller than the original, which retains the basic topology of P-that is, the number of recurrent classes and their respective reducibilities and periodicities (see Barreto and Fragoso's [2011] article for formal definitions). Since ¯P can be considerably smaller than P, the idea of replacing the latter with the former comes almost inevitably: this is the "stochastic-factorization trick."

3 Stochastic Factorization as a Probabilistic Model

So far, an algebraic view of the stochastic factorization has prevailed in the literature, with P, D, and K treated simply as matrices with a particular structure. In this paper we introduce a probabilistic interpretation of a stochastic factorization. In particular, we see the factorization as a Markov model that represents the joint distribution of a sequence of random variables.

3.1 Stochastic FactorizationModel

Let S ≡ {1, 2, ..., n}, H ≡ {1, 2, ...,m}, and A ≡ {1, 2, ..., na}. Consider the stochastic process (S1,A1,H1, S2,A2,H2, ...), where St ∈ S are observable states, At ∈ A are observable actions (or decisions), and Ht ∈ H are hidden states. Our probabilistic model builds on the following Markov assumptions:

Given τ > 0, it is easy to compute the probability of any sequence of length τ using (i), (ii), and (iii):

Thus, in order to define our model we only need to represent Pr(S1), Pr(St+1|Ht,At), Pr(At|St), and Pr(Ht|St,At). This leads to the following definition:

Definition 4 A stochastic factorization model (SFM) is a tuple F ≡ (Da ,Ka ,Π,µ), where Da ∈ Rn×m and Ka ∈ Rn×m are na stochastic matrices, one for each a ∈ A, Π ∈ Rn×na is a stochastic matrix, and µ ∈ Rn is a distribution. Based on Assumptions (i), (ii), and (iii), the model represents the stochastic process (S1, A1, H1, S2, A2, H2, ...) as daij = Pr(Ht = j|St = i,At = a), ka ij = Pr(St+1 = j|Ht = i,At = a), πia = Pr(At = a|St = i), and µi = Pr(S1 = i).

Note that a SFM with na = 1 is an uncontrolled system. In order to provide some intuition about our model, in Figure 1 we compare it with another uncontrolled system, the well known hidden Markov model (HMM) [Rabiner, 1989]. As shown, the main difference between the two models is the fact that in a SFM knowledge of the current observable state is sufficient to define the distribution of the next hidden state.

3.2 Computing Probabilities

Recall that in a SFM only the variables St and At are observable. Let z1:τ = (s1, a1, s2, ..., aτ-1, sτ). In order to do inference with our model, we need to be able to compute the probability of a given sequence happening, that is, Pr(Z1:τ = z1:τ |λ), where λ ≡ (Da ,Ka ,Π,µ). In this section we discuss how to compute this quantity. We base our strategy on the forward-backward procedure to compute probabilities in an HMM [Baum, 1972; Rabiner, 1989].

We start by noting that, if we can compute Pr(Z1:τ = z1:τ , Ht = i|λ), then we can compute the desired quantity by marginalizing over Ht. Given z1:τ , we define the "forward variable" αi(t) as:

In Table 1 we show that, given a SFM, we can compute αi(t) for any i and any t using the following recursive formulae:

Table 1: Computing probabilities in a SFM; the dependence on λ is omitted to improve readability.

We now define the "backward variable" βi(t) as:

As with αi(t), given a SFM, the variable βi(t) can be computed recursively, in the following way (see Table 1):

Now that we are able to compute αi(t) and βi(t) for any i and any t, we can compute the desired probability by simply multiplying these variables, that is,

The equality above is also justified in Table 1.

Scaling

Using the variables αi(t) and βi(t) we can compute Pr(Z1:τ ,Ht = i|λ), which in turn allows us to compute the probability of any sequence being generated by a SFM. However, since the variables αi(t) get very small as t →∞ (and equivalently for βi(t) as t → 0), in practice the recursive equations presented above may lead to numerical instability. We now present a simple workaround, based on Rabiner's [1989] solution for a similar problem, which consists in rescaling the variables αi(t) and βi(t) after each recursion.

Let wt = ∑i αi(t). Starting with ˆαi(1) = αi(1)/w1, we will apply the following recursion, for t ≥ 2:

where w′t = ∑i α′i (t). Using induction, it is easy to show that

and thus the variables ˆαi(t) define a distribution (and are in general representable with standard machine precision). From Table 1 we know that

which means that, as long as the variables w′j are available, we can use the scaled forward variables to compute the probability of any sequence in a SFM.

Although we do not need the backward variables to compute Pr(Z1:τ |λ), it will be convenient to also define their scaled version. Starting with ˆ βi(τ - 1) = βi(τ - 1) and ˆ βi(τ - 2) = βi(τ -2)/wτ-1, we will apply the following recursion, for t < τ - 2:

Note that, since the magnitudes of β′ i (t - 1) and α′i (t) are comparable, the computation of ˆ βi(t) will also be numerically stable. One can show that

From (6) and (9) we know that ˆαi(t)ˆ βi(t) = ∏τ-1 j=1 (w′j )-1αi(t)βi(t), which, combined with (4), implies that

Equation (10) will be useful for learning a SFM from data, as we discuss next.

4 Expectation-Maximization Algorithm

One of the motivations for introducing the SFM is to provide a framework for the application of the stochastic-factorization trick directly from data. The idea is as follows: given a sequence of observation-action pairs coming from a set of na transition matrices Pa , instead of estimating the matrices Pa directly, we estimate a SFM F ≡ (Da , Ka , Π, µ). By doing so, we automatically have approximations DaKa ≈ Pa .

The strategy above has two benefits. First, by estimating Da and Ka instead of Pa , we reduce the number of parameters we are estimating from n2 to 2nm. Second, once we have the SFM, we do not need to actually compute the multiplication DaKa : instead, we compute ¯Pa = KaDa and use this smaller matrix in place of Pa .

In this section we discuss how to compute a SFM directly from data. Since SFMs are a latent-variable model, we resort to the well know expectation-maximization (EM) algorithm [Dempster et al., 1977]. We assume the reader is familiar with the basic principles of EM; for good reviews of the subject see for example Bishop's [2006] book or Bilmes's [1998] tutorial.

4.1 Derivation

Given a sequence of observation-action pairs, z1:τ , our goal is to compute a SFM that maximizes the likelihood of the data, L(λ|z1:τ ) = Pr(z1:τ |λ). As well known, the EM algorithm does this through the function

In the "E" step we use the current parameter values λ′ to compute the expectation above; in the "M" step we maximize the expectation we computed with respect to λ. In this section we will show how to derive update equations to compute a SFM based on the above framework. In order to improve readability we will drop the subscript 1 : τ throughout the section, that is, we will use h, Z, and z to refer to h1:τ , Z1:τ , and z1:τ .

As shown in (1), in the case of a SFMwe have

Substituting (12) in (11), the Q function becomes

(recall that here H represents the space of all possible sequences h of length τ). We will focus on the third term of the equation above, tagged as "∆D," to show how to derive update equations for the elements of the matrices Da ; the update equations for the remaining parameters of a SFM are obtained in an analogous way. Our derivation closely follows that of Bilmes [1998].

Given z, suppose that st = i and at = u. Then, for each value of j ∈ {1, 2, ...,m}, the corresponding term log duij appearing in ∆D will multiply Pr(h|Z, λ′ ) if and only if the tth element in h is j-that is, duij will multiply ∑h∈H Pr(h|Z, λ′ )1{ht = j} = Pr(Ht = j|Z, λ′ ) (here 1{·} is the indicator function). Thus, we can rewrite ∆D as:

We want to maximize Q(λ, λ′ ) with respect to λ. One way to do so is to add the Lagrange multipliers κui to account for the constraints ∑j duij = 1, which leads to

Now, if we make ∂∆′D /∂duij = 0, we have a necessary condition for duij to be an extreme point of Q. Since

we know that

If we sum (13) over j, using the fact that ∑j duij = 1, we obtain κui = ∑τ-1 t=1 1{st = i, at = u}. Substituting κui back in (13), we have the following update rule for duij :

If we follow the same steps as above replacing ∆D with ∆K we arrive at the following update equation for ku ij :

As discussed in Section 3.2, we can compute Pr(Ht = i|Z, λ′ ) using (10).

Update equations (14) and (15) make intuitive sense: loosely speaking, they are the equations for computing the relative frequency of events with the occurrence of Ht = i replaced by its probability. Based on this insight, and following Rabiner's [1989] reasoning, we can extend our algorithm to the case in which we have c sequences of observations and actions, zl1:τl = {sl1 , al1 , ..., alτl-1 , slτl }, with l = 1, 2, ..., c. In this case we have the following update rules:

In principle, we could follow the steps above to derive update equations for µi and πiu. Note though that, since these parameters are based on observable quantities, the associated update equations will be simply the relative frequency of the corresponding events, that is:

Since the quantities appearing in (18) do not depend on the parameters of the SFM, the estimates of µ and Π will be fixed throughout the EM iterations. Note that in scenarios where we have control over the collection of data we generally know µ and Π, which means that these parameters do not have to be estimated at all.

To conclude, we note that for both practical and theoretical reasons the stochastic-factorization trick is sometimes applied with a single D or a single K (see Barreto et al., 2011 and Barreto et al., 2014, respectively). One can easily specialize the proposed method to these particular cases. Having a single K, for example, corresponds to replacing Assumption (ii) with Pr(Ht|At, St, Ht-1, ..., S1) = Pr(Ht|St). The resulting update equation will be expression (17) without the restriction alt = u.

4.2 Algorithm

Algorithm 1 shows a step by step description of the proposed method, which we call EMSF. The pseudo-code uses two matrices to represent the forward and backward variables: α ∈ Rτ-1×m, where αti = ˆαi(t), and β ∈ Rτ-1×m, where βti = ˆ βi(t). For each sequence zl1:τl , EMSF first computes α using (5) and then computes β using (8). Then, α and β are multiplied element-wise, giving rise to matrix C ∈ Rτ-1×m. After the rows of C have been normalized, as shown in line 39 of Algorithm 1, we have cti = Pr(Ht = i|Z1:τl = zl , λ′ ) (see (10)). These values are then used to compute the numerator of update equations (16) and (17), which are accumulated in the auxiliary matrices ˆDa ∈ Rn×m and ˆKa ∈ Rm×n . After all the sequences have been processed, the rows of ˆDa and ˆKa are normalized (lines 51 and 52 of Algorithm 1), which corresponds to completing the application of (16) and (17) to obtain the new estimates of Da

Algorithm 1 EMSF

1: Input: {z1 1:τ1 , z2 1:τ2 , ..., zc 1:τc } ▷ Data m ∈ N ▷ Order of the factorization µ ∈ Rn ▷ Initial distribution Π ∈ Rn×|A| ▷ Exploration policy ϵ > 0 ▷ Stop criterion 2: Output: DaKa ≈ Pa , for a ∈ A 3: 4: for each a ∈ A do 5: Da ← random stochastic matrix in Rn×m 6: Ka ← random stochastic matrix in Rm×n 7: 8: L← 0; L′ ←∞ ▷ L is the data log likelihood 9: while |L-L′ | > ϵ do 10: L′ ←L ; L← 0 11: 12: for each a ∈ Ado ˆDa ← 0 ∈ Rn×m; ˆKa ← 0 ∈ Rm×n 13: 14: for l ← 1 to c do 15: τ ← τl 16: s ← {sl1, sl2, ..., slτ } ▷ Extract sli:τl from zl1:τl 17: a ← {al1, al2, ..., alτ-1} ▷ Extract ali:τl-1 from zl1:τl 18: α ← 0 ∈ Rτ-1×m 19: β ← 0 ∈ Rτ-1×m 20: w′ ← 0 ∈ Rτ-1 21: ▷ Compute α 22: 23: (α)1 ← µ1πs1a1 (Da1 )s1 ▷ (α)i is α's ith row 24: w′1 ← ∑i α1i 25: (α)1 ← (α)1/w′1 26: for t ← 2 to τ - 1 do 27: (α)t ← (∑i α(t-1)ikat-1 ist ) πstat (Dat )st 28: w′t ← ∑i αti 29: (α)t ← (α)t/w′t 30: ▷ Compute β 31: 32: (β)τ-1 ← (Kaτ-1 )sτ ▷ (K)i is K's ith column 33: for t ← τ - 2 to 1 do 34: (β)t ← ∑i β(t+1)idat+1 (st+1)iπ(st+1)(at+1)(Kat )st+1 35: (β)t ← (β)t/w′t+1 36: ▷ Update ˆDa and ˆKa 37: 38: C ← α ⊗ β ▷ Element-wise multiplication 39: for i ∈ 1, 2, ..., τ - 1 do (C)i ← (C)i/ ∑j cij 40: for t ← 1 to τ - 1 do 41: ( ˆDat )st ← ( ˆDat )st + (C)t 42: ( ˆKat )st+1 ← ( ˆKat )st+1 + (C)t 43: 44: ▷ Update L based on (7) 45: 46: L ← L + ∑i log w′i + log ((α)τ-1(Kaτ-1 )sτ ) 47: 48: ▷ Update Da and Ka 49: 50: for each a ∈ A do 51: for i ∈ 1, 2, ..., n ( ˆDa )i ← ( ˆDa )i/ ∑j ˆ daij 52: for i ∈ 1, 2, ..., m ( ˆKa )i ← ( ˆKa )i/ ∑j ˆ ka ij 53: Da ← ˆDa ; Ka ← ˆKa 54: 55: returnDa and Ka , with a ∈ A

and Ka . The process continues until the decrease on the data log likelihood falls below a given threshold.

EMSF is an expectation-maximization algorithm, and thus it will only converge to a local maximum of L. The solution the algorithm converges to is determined by the initial values of Da and Ka . Besides its initialization, the approximation computed by EMSF depends only on one parameter, m. The parameter m controls the complexity of the SFM computed by EMSF, and thus it can be seen as a "regularizer" of the approximations DaKa ≈ Pa (see discussion in Section 5) [Gy¨orfi et al., 2002; Hastie et al., 2002].

As for computational complexity, each iteration of EMSF is O(τm), where τ = ∑l τl. Thus, the cost of the algorithm is not a direct function of n, the size of the matrices Pa being estimated. Of course, n will play a role in the definition of an appropriate value for τ, since EMSF estimates 2nanm parameters.

5 Experiments

In this section we present computational experiments to illustrate the properties and usefulness of our algorithm. The first experiment is a "proof of concept": we generated transition matrices P ∈ R100×100 , with srk(P) = 20, and tried to recover them with EMSF using different values for m. We then compared EMSF's results with those obtained by directly estimating P via maximum likelihood (referred to in the plots as CNT, for "counting"). The results are shown in Figure 2.

Since EMSF estimates fewer parameters than CNT, its approximation error decreases faster with the number of sample transitions. However, since the EM algorithm only converges to a local maximum of L, EMSF's approximation error stagnates at a positive value, even when m = srk(P), while CNT's approximation keeps improving as τ →∞.

These results are exactly as expected. As CNT is statistically consistent, we know that its approximation error will converge to zero, but this comes at the price of higher estimation errors at lower values of τ . EMSF reduces the estimation error by enforcing a specific structure on P-in terms of the classical bias-variance analysis, EMSF increases bias and decreases variance when m < n [Hastie et al., 2002].

This view of EMSF also helps to explain the impact of the parameter m over its performance. Since m directly controls the trade off between bias and variance, when computing a stochastic factorization from data the "optimal" value for this parameter depends not only on srk(P), but also on the amount of data available. In general, as the number of transitions increases so does the best value for m (this is, again, exactly the trend observed in Figure 2). The definition of an appropriate value for m is thus an instance of the well known model selection problem; one can address this problem using cross validation [Hastie et al., 2002].

It is clear then that EMSF can be seen as a strategy to regularize the approximation of P in order to reduce the number of transitions needed for learning. An interesting question is whether the structure imposed on P by our algorithm really arises in practice. In order to answer this question, we used EMSF to compute a model of the game of blackjack.

The game of blackjack was implemented exactly as described in Section 5.1 of Sutton and Barto's [1998] book. The experiments were carried out as follows. First, an exploration policy that selects actions uniformly at random played a number of hands of blackjack. The resulting sample transitions were then used by CNT to compute estimates of both the transition matrices Pa and the expected-reward vectors ra of the Markov decision process (MDP) describing the game. Finally, policy iteration (PI) was applied to the MDP to compute a policy, which was evaluated on 106 hands of blackjack.

In order to evaluate EMSF we used an algorithm introduced by Barreto et al. [2014] called policy iteration based on stochastic factorization (PISF). Given DaK ≈ Pa and Da¯r ≈ ra , PISF computes an approximation of the MDP's value function, ˜v ≈ v∗ , in O(n) time per iteration; Barreto et al. [2014] have shown that ‖ v∗ - ˜v ‖∞ → 0 as ‖ Pa -DaK ‖∞ → 0 and ‖ ra -Da¯r ‖∞ → 0 for all a ∈ A. EMSF computed the factorization used by PISF from the same data used by CNT. Since in blackjack the rewards are completely defined by the end state of a transition, we can compute ¯r as ¯r = Kˆr′ , where ˆ r′i is an estimate of the reward associated with transitions ending in state i, ˆ ri. This is clear when we note that ra = Paˆr ≈ DaKˆr′ .

Figure 3 compares the results obtained by CNT+PI and EMSF+PISF on blackjack. Both methods are contrasted with an agent that uses the same strategy used by the game's dealer [Sutton and Barto, 1998]. Note how the results shown in Figure 3 reproduce the trend seen in Figure 2. In particular, after playing only 3000 hands, all the EMSF+PISF agents have already outperformed the dealer's strategy, while CNT+PI needs twice as many hands to reach the same level of playing. The superior performance of EMSF+PISF is maintained up until 30000 hands played. These results show that the structural assumption underlying EMSF, a stochastic factorization, do arise in real applications-which in turn suggests that our algorithm can be useful in practice.

6 RelatedWork

Stochastic factorization is a particular case of nonnegative matrix factorization [Paatero and Tapper, 1994; Lee and Seung, 1999]. In fact, Cohen and Rothblum [1991] have shown that one can always derive a stochastic factorization from a nonnegative factorization of a stochastic matrix. In principle, thus, any algorithm designed for the latter can also be used to compute the former (see discussion in Barreto et al.'s [2014] paper). Note though that algorithms designed for nonnegative matrix factorization assume that we have access to the elements of P (or at least to a subset of them). In contrast, EMSF has no access to any element of P, only to data extracted from this matrix (this is possible precisely because of the extra structure assumed by stochastic factorization, since it is unclear how to extract data from a nonnegative matrix that is not stochastic).

As mentioned in the introduction, in the scenario considered here it is possible to start from an estimate of P and then use it to compute D and K through a nonnegative matrix factorization algorithm. Similar ideas have been explored by Finesso et al. [2010], Lakshminarayanan and Raich [2010], and Cybenko and Crespi [2011] in the context of HMM learning-except that instead of estimating P these authors estimate alternative matrix representations of the observation dynamics. A potential problem with applying these or similar approaches to learn a SFM is that the matrices representing the observation dynamics are at least n×n, which amounts to estimating O(n2 ) parameters or more-exactly what EMSF is trying to avoid.

EMSF has clear similarities with the well known BaumWelch algorithm, an EM method to compute the parameters of an HMM from data [Baum, 1972]. Given this similarity, and considering the connection between SFMs and HMMs, shown in Figure 1, it is natural to ask whether other methods to learn an HMM could also be adapted to compute a stochastic factorization (since the literature on this subject is vast, we redirect the reader to Capp´e et al.'s [2005] book to serve as a starting point).

From a broader perspective, SFMs are a particular instance of latent variable models. Recently, there has been renewed interest in the development of algorithms to learn the parameters of such models, motivated by a reduction of the problem to a singular value decomposition of some matrix representing the observation dynamics [Hsu et al., 2009; Siddiqi et al., 2010; Song et al., 2010; Bailly, 2011; Balle et al., 2012]. Unlike EM, these so-called "spectral methods" always converge to a global optimum, and are also statistically consistent (under some reasonable assumptions). As with the algorithms that learn an HMM based on a nonnegative matrix factorization, the main difficulty in adapting spectral methods to learn a SFM is the fact that they build on the decomposition of a matrix that is at least n × n, which means they require estimating at least n2 parameters. In any case, the development of alternative methods to learn a SFM that overcome the limitations of EM is a promising topic for future research.

7 Conclusion

This paper introduced SFM, a probabilistic model that encompasses a stochastic factorization, and presented an expectation-maximization algorithm, EMSF, to compute the model's parameters from data. Our algorithm has a single parameter, m, which can be seen as a practical mechanism to control the bias-variance trade off in the approximation of a transition matrix. This makes it possible to adjust the complexity of the model to the amount of data available; if the transition matrix is factorizable or nearly so, one can compute good approximations using less data than would be needed for direct maximum likelihood estimation. Our experiments suggest that factorizable transition matrices arise in practice.

One interesting direction for future research is to extend the ideas presented in this paper to continuous state spaces. Another promising line of investigation is the development of an incremental version of EMSF which makes it possible to incorporate new data into the model without the need to store the sample transitions already used [Khreich et al., 2012]. This extension would allow EMSF to scale to much larger datasets and also to be applied on-line, potentially using the model already constructed to guide the collection of new sample transitions [Kearns and Singh, 1998; Brafman and Tennenholtz, 2003]. Finally, another interesting direction for future research is the study of SFMs with factored dynamics, in which the model's conditional probabilities can be represented by a dynamic Bayesian network [Guestrin et al., 2003]. This would bridge the gap between "factorizable" and "factored" models, potentially decreasing even further the number of samples needed to learn a transition matrix [Barreto et al., 2014].

Acknowledgments

The authors would like to thank Borja Pigem for helpful discussions regarding topics related to this work. Funding for this research was provided by Conselho Nacional de Desenvolvimento Cient´ıfico e Tecnol´ogico (CNPq), grant 461739/2014-3, and by the NSERC Discovery grant program.