Year: 2020

Abstract:
Model-based reinforcement learning has been empirically demonstrated as a successful strategy to improve sample efficiency. In particular, Dyna is an elegant model-based architecture integrating learning and planning that provides huge flexibility of using a model. One of the most important components in Dyna is called search-control, which refers to the process of generating state or state-action pairs from which we query the model to acquire simulated experiences. Searchcontrol is critical in improving learning efficiency. In this work, we propose a simple and novel search-control strategy by searching high frequency regions of the value function. Our main intuition is built on Shannon sampling theorem from signal processing, which indicates that a high frequency signal requires more samples to reconstruct. We empirically show that a high frequency function is more difficult to approximate. This suggests a search-control strategy: we should use states from high frequency regions of the value function to query the model to acquire more samples. We develop a simple strategy to locally measure the frequency of a function by gradient and hessian norms, and provide theoretical justification for this approach. We then apply our strategy to search-control in Dyna, and conduct experiments to show its property and effectiveness on benchmark domains.

heading: 1 INTRODUCTION
text: Model-based reinforcement learning (MBRL) (Lin, 1992; Sutton, 1991b; Daw, 2012; Sutton & Barto, 2018) methods have successfully been applied to many benchmark domains (Gu et al., 2016; Ha, David and Schmidhuber, Jürgen, 2018; Kaiser et al., 2020). The Dyna architecture, introduced by Sutton (1991a), is one of the classical MBRL architectures, which integrates model-free and model-based policy updates in an online RL setting (Algorithm 2 in Appendix A.3). At each time step, a Dyna agent uses the real experience to learn a model and to perform model-free policy update, and during the planning stage, simulated experiences are acquired from the model to further improve the policy. A closely related method in model-free learning setting is experience replay (ER) (Lin, 1992; Adam et al., 2012), which utilizes a buffer to store experiences. An agent using the ER buffer randomly samples the recorded experiences at each time step to update the policy. Though ER can be thought of as a simplified form of MBRL (van Seijen & Sutton, 2015), a model provides more flexibility in acquiring simulated experiences.
A crucial aspect of the Dyna architecture is the search-control mechanism. It is the mechanism for selecting states or state-action pairs to query the model in order to generate simulated experiences (cf. Section 8.2 of Sutton & Barto 2018). We call the corresponding data structure for storing those states or state-action pairs the search-control queue. Search-control is of vital importance in Dyna, as it can significantly affect the model-based agent’s sample efficiency. A simple approach to searchcontrol is to sample visited states or state-action pairs, i.e., use the initial state-action pairs stored in the ER buffer as the search-control queue. This approach, however, does not lead to an agent that outperforms a model-free agent that uses ER. To see this, consider a deterministic environment, and assume that we have the exact model. If we simply sample visited state-action pairs for searchcontrol, the next-state and reward would be the same as those in the ER buffer. In practice, we have
∗Equal contribution.
model errors too, which causes some performance deterioration (Talvitie, 2014; 2017). Without an elegant search-control mechanism, we are not likely to benefit from the flexibility given by a model.
Several search-control mechanisms have already been explored. Prioritized sweeping (Moore & Atkeson, 1993) is one such method that is designed to speed up the value iteration process: the simulated transitions are updated based on the absolute temporal difference error. It has been adopted to continuous domains with function approximation too (Sutton et al., 2008; Pan et al., 2018; Corneil et al., 2018). Gu et al. (2016) utilizes local linear models to generate optimal trajectories through iLQR (Li & Todorov, 2004). Pan et al. (2019) suggest a method to generate states for the searchcontrol queue by hill climbing on the value function estimate.
This paper proposes an alternative perspective to design search-control strategy: we can sample more frequently from the state space where the value function is more difficult to estimate. We first review some basic background in MBRL (Section 2). Afterwards, we review some concepts in signal processing and conduct experiments in the supervised learning setting to show that a high frequency function is more difficult to approximate (Section 3). In order to quantify the difficulty of estimation, we borrow a crucial idea from the signal processing literature: a signal with higher frequency terms requires more samples for accurate reconstruction. We then propose a method to locally measure the frequency of a point in a function’s domain and provide a theoretical justification for our method (Theorem 1 in Section 3.2). We use the hill climbing approach of Pan et al. (2019) to adapt our method to design a search-control mechanism for the Dyna architecture (Section 4). We conduct experiments on benchmark and challenging domains to illustrate the properties and utilities of our method (Section 5).

heading: 2 BACKGROUND
text: Reinforcement learning (RL) problems are typically formulated as Markov Decision Processes (MDPs) (Sutton & Barto, 2018; Szepesvári, 2010). An MDP (S,A,P, R, γ) is determined by state space S, action space A, transition function P, reward function R : S × A × S → R, and discount factor γ ∈ [0, 1]. At each step t, an agent observes a state st ∈ S , and takes an action at ∈ A. The environment receives at, and transits to the next state st+1 ∼ P(·|st, at). The agent receives a reward scalar rt+1 = R(st, at, st+1). The agent maintains a policy π : S × A → [0, 1] that determines the probability of choosing an action at a given state. For a given state-action pair (s, a), the action-value function of policy π is defined as Qπ(s, a) = E[Gt|St = s,At = a;At+1:∞ ∼ π] where Gt def = ∑∞ t=0 γ
tR(st, at, st+1) is the return of s0, a0, s1, a1, ... following the policy π and transition P. Value-based RL methods learn the action-value function (Watkins & Dayan, 1992), and act greedily w.r.t. the action-value function. Policy-based RL methods perform gradient update of parameters to learn policies with high expected rewards (Sutton et al., 1999). Both value and policy-based RL methods can be easily adopted in the Dyna framework.
Model-based RL. A model is a mapping that takes a state-action pair as its input and outputs some projection of the future state. A model can be local (Tassa et al., 2012; Gu et al., 2016) or global (Ha, David and Schmidhuber, Jürgen, 2018; Pan et al., 2018), deterministic (Sutton et al., 2008) or stochastic (Deisenroth & Rasmussen, 2011; Ha, David and Schmidhuber, Jürgen, 2018), feature-to-feature (Corneil et al., 2018; Ha, David and Schmidhuber, Jürgen, 2018) or observationto-observation (Gu et al., 2016; Pan et al., 2018; Kaiser et al., 2020), one-step (Gu et al., 2016; Pan et al., 2018), or multi-step (Sorg & Singh, 2010; Oh et al., 2017), or decision-aware (Joseph et al., 2013; Farahmand et al., 2017; Silver et al., 2017). Modelling the environment dynamics through a reproducing kernel Hilbert space (RKHS) embedding has been also studied (Grunewalder et al., 2012), where the Bellman operator is approximated in an RKHS. The model we consider in this work is a one-step environment dynamics model, which takes a state-action pair as its input and returns the next-state and reward. Our proposed search-control approach, however, can be naturally used for different types of models.
The most relevant work to ours is hill climbing Dyna (Pan et al., 2019). Pan et al. (2019) proposes a search-control mechanism based on hill climbing on the value estimates (see Algorithm 3 in Appendix A.3). We briefly review the key steps of their algorithm, which is called (Hill Climbing)HCDyna, as it helps to understand ours. HC-Dyna maintains an ER buffer. At each step, a state is randomly sampled from the ER buffer and is used as the initial state to perform hill climbing (i.e.
gradient ascent) on the learned value function. The states along the trajectory are stored in the search-control queue.1 During the planning stage, states are sampled from the search-control queue and are paired with their corresponding on-policy actions (i.e., actions selected by the current Q network at the sampled states). Afterwards, the model is queried for each of the state-action pairs to get the next-state and reward. These simulated transitions are then mixed with samples from the ER buffer, which are observed by the agent during its interaction with the real environment, to train the value function estimator, e.g., a deep neural network.
The heuristic idea behind the search-control mechanism of HC-Dyna is that the magnitude of the value function provides useful information for guiding where to query the model. This heuristic can intuitively be understood by noticing that an RL agent tends to move towards high-value regions of the state space; by performing gradient ascent on the (estimated) value function, we provide the agent with more samples from regions where it may move towards in the future. Even if the estimated value function is incorrect and the samples are indeed from the low-value regions of the state space, these extra samples lead to the fast correction of the estimated value in those regions. Nevertheless, the magnitude of the value function is only one source of extra information from which we can design a search-control mechanism. This work suggests a different perspective: we should sample more from the regions of the state space where learning the value function is more difficult.

heading: 3 UNDERSTANDING THE DIFFICULTY OF FUNCTION APPROXIMATION
text: In a regular regression setting, we illustrate that high frequency regions of a function is difficult to approximate. We show that by assigning more training data to those regions, the learning performance considerably improves. To make this insight practically useful, we employ the sum of gradient and hessian norms of a function as a measure of the local frequency of a function. We establish a theoretical connection between our proposed criterion and the local frequency of a function. This would be the foundation of our frequency-based search-control method in Section 4.

heading: 3.1 WHAT TYPE OF FUNCTION IS DIFFICULT TO APPROXIMATE?
text: Consider the standard regression problem with the mean square loss. Given a training set D = {(xi, yi)}i=1:n, our goal is to learn an unknown target function f∗(x) = E[Y |X = x] by empirical risk minimization. Formally, we aim to solve
f = arg min f∈H
1
n n∑ i=1 (f(xi)− yi)2,
where H is some hypothesis space. Suppose that we can choose the distributions of samples {xi}. How should we select them in order to improve the quality of the learned function? One intuitive heuristic is that if we know the regions in the domain of f∗ that are more difficult to approximate, we can assign more training data there in order to help the learning process. The important question is how to quantify the difficulty of approximating a function. We borrow an idea from the field of signal processing to suggest a method.
The Nyquist-Shannon sampling theorem in signal processing states that given a band-limited function (or signal) f : R → R with the highest frequency (in the Fourier domain) of ωbandwidth, we can perfectly reconstruct it based on regular samples (in the time domain) obtained at the sampling rate of 2ωbandwidth (Zayed, 1993).2 Therefore, if the Fourier transform of a function has high frequency terms, more samples are required to reconstruct it accurately. We note that the sampling theory has been applied in the sample complexity analysis of machine learning algorithms (Smale & Zhou, 2004; 2005; Jiang, 2019). Although the problem setting in machine learning is somewhat different from this result in signal processing, it still provides a high-level intuition for us: regions with more high frequency signals require more learning data.
1According to the original paper, natural gradient is used to guarantee a certain level of coordinate invariance property, so it can handle state variables with vastly different numerical scales.
2Sampling rate refers to number of samples per second used to reconstruct continuous signals.
To make this high-level intuition concrete, we consider the following function:
fsin(x) = { sin(8πx) x ∈ [−2, 0), sin(πx) x ∈ [0, 2]. (1)
It is easy to check that the regions [−2, 0) and [0, 2] contain signals with frequency ratio 8 : 1. Based on the intuition from the sampling theorem, the [−2, 0) interval requires more training data than the [0, 2] interval. Given the same amount of training data, and the same learning algorithm, we would expect that assigning more fraction of the training data on [−2, 0) to perform better than distributing them uniformly or assigning more samples to the [0, 2] interval.
An illustrative experiment. To empirically verify the intuition, we conduct a simple regression task, with fsin as the target function. The training set D = {(xi, yi)}i=1:n is generated by sampling x ∈ [−2, 2], and adding Gaussian noise N(0, σ2) on Eq. (1), where the standard deviation is set to be σ = 0.1. We present the `2 regression learning curves of training datasets with different biased sampling ratios pb ∈ {60%, 70%, 80%}, as shown in Fig. 1 (a)-(c). We observe that biased training data sampling ratios towards high frequency region clearly speeds up learning. This is consistent with the intuitive insight and suggests that our heuristic to assign more data to high frequency regions leads to faster learning.

heading: 3.2 IDENTIFYING HIGH FREQUENCY REGIONS OF A FUNCTION
text: Identifying the high frequency region of fsin in the previous toy problem was easy, as each region contained a signal with a constant known frequency. In practice, we face two main difficulties to identify the high frequency regions of a function. The first is that we do not have access to the underlying target function, but only to data or possibly an approximate function that is estimated using data, e.g., a trained neural network. The second is that frequency is a global property rather than a local one. The value of the function at each (non-zero measure) region of the domain has impact on its global frequency representation. To make the high frequency heuristic practically useful, we need a simple criterion that (a) uses function approximation, (b) characterizes local frequency information, and (c) can be efficiently calculated. Inspired by the function fsin in Eq. (1), a natural idea is to calculate the first order f ′(x) def= df(x)dx or second order derivative f ′′(x) def = d 2f(x) dx2 because they both satisfy (a) and (c). As a “sanity check” for property (b), consider the following examples.
Example 1. For fsin defined in Eq. (1), calculate the integrals of squared first order derivative f ′sin on high frequency region [−2, 0) and low frequency region [0, 2], respectively:∫ 0
−2 |f ′sin(x)| 2 dx = 64π2, ∫ 2 0 |f ′sin(x)| 2 dx = π2.
Example 2. Let f : [−π, π]→ R be a band-limited real valued function defined as
f(x) = a0 2 + N∑ n=1 an cos (nx) + bn sin (nx),
where a0, an, bn ∈ R, n = 1, 2, . . . , N are Fourier coefficients of frequency n2π . Then,∫ π −π |f ′(x)|2 dx = π · N∑ n=1 n2 ( a2n + b 2 n ) , ∫ π −π |f ′′(x)|2 dx = π · N∑ n=1 n4 ( a2n + b 2 n ) .
Example 1 shows that the integral of squared first order derivative ratio is 64 : 1 (the frequency ratio is 8 : 1), and the region with large derivative magnitude is indeed the high frequency region. Moreover, Example 2 indicates that for one dimensional real-valued functions over a bounded domain, the integral of a derivative magnitude is closely related to the frequency information. For the squared derivative, the integral is the same as weighting the frequency terms an and bn proportional to n2, and for the squared second-order derivative, the integral is the same as weighting the frequency terms proportional to n4. The weighting schemes n2 or n4 emphasize the higher frequency terms.
Empirical demonstration. Our calculation in the above examples implies that regions with large gradient and Hessian norm correspond to high frequency regions. Based on the same spirit of the l2 regression task in Section 3.2, we empirically verify this insight. Our expectation is that biasing training dataset towards high gradient norm and Hessian norm would achieve better learning results. In Fig. 2(a), Biased-GradientNorm corresponds to uniformly sampling x ∈ [−2, 2] for 60% of training data and sampling proportional to gradient norm (i.e., p(x) ∝ |f ′sin(x)|) for the remaining 40%; while Biased-HessianNorm corresponds to sampling proportional to Hessian norm (i.e., p(x) ∝ |f ′′sin(x)|) for the remaining 40% of training data. In Fig. 2(b)(c), we visualize the two types of biased training points. Sampling according to the gradient norm or the Hessian norm leads to denser point distribution in the high frequency region [−2, 0): there are 65.35%, 68.97% of training points fall in [−2, 0] in Fig. 2(b), (c) respectively. An important difference between Fig. 2(b) and (c) is that, sampling according to Hessian norm leads to denser points around spikes: there are 18.17% points fall in the yellow area in (b) and 27.45% such points in (c). Those areas around spikes should be more difficult to approximate as the underlying function changes sharply, which explains the superior performance on the data set biased by Hessian norm. Fig. 2(a) shows that such biased training datasets provide fast learning, similar to the high frequency biased training datasets in Fig. 1.
Given a function f : X 7→ Y and a point x ∈ X , we propose to measure frequency of f around a small neighborhood of x (we call this local frequency) using the following function:
g(x) def = ‖∇xf(x)‖2 + ‖Hf (x)‖2F , (2)
where ‖∇xf(x)‖ is the `2-norm of the gradient at x, and Hf (x)F is the Frobenius norm of the Hessian matrix of f at x. We claim that local frequency of f around x is proportional to g(x).We theoretically justify this claim. For real-valued functions in Euclidean spaces, our theory connects local gradient and Hessian norms, local function energy 3, to local frequency distribution. The proof of our theorem and its connection to the well-known uncertainty principle are in Appendix A.2. Theorem 1. Given any function f : Rn → R, for any frequency vector k ∈ Rn, define its local Fourier transform around x ∈ Rn,
f̂(k) def = ∫ y∈B(x,1) f(y) exp { −2πi · y>k } dy,
for local function around x, i.e., y ∈ B(x, 1) def= {y : ‖y − x‖ < 1}. Assume the local function “energy” is finite, ∫
y∈B(x,1) [f(y)]
2 dy = ∫ Rn ‖f̂(k)‖2dk <∞, ∀x ∈ Rn. (3)
3We consider the notion of energy in signal processing terminology: the energy of a continuous time signal x(t) is defined as ∫ x(t)2dt. In our theory, function f is the signal.
Define “local frequency distribution” of f(x) as:
πf̂ (k) def = ‖f̂(k)‖2∫ Rn ‖f̂(k̃)‖2dk̃ , ∀k ∈ Rn. (4)
Then, for any x ∈ Rn, we have: 1) The first order connection:∫
y∈B(x,1) ‖∇f(y)‖2 dy = 4π2 · [∫ y∈B(x,1) [f(y)] 2 dy ] · [∫ Rn πf̂ (k) · ‖k‖ 2 dk ] , (5)
2) The second order connection:∫ y∈B(x,1) ||Hf (y)||2F dy = 16π4 [∫ y∈B(x,1) [f(y)] 2 dy ] · [∫ Rn πf̂ (k) · ‖k‖ 4 dk ] (6)
Remark 1. Note that πf̂ defined in Eq. (4) is a probability distribution over R n as:∫
k∈Rn πf̂ (k)dk = 1, and πf̂ (k) ≥ 0, ∀k ∈ R n.
We use such a distribution to characterize local frequency behaviour for reasons. First, comparing frequencies of regions is more naturally captured by a distribution than one single scalar, since signals usually are within a range of frequencies. Second, to eliminate the impact of the function energy Eq. (3), we normalize the Fourier coefficient f̂ to get πf̂ .
Remark 2. For a frequency vector k ∈ Rn, the larger its norm ‖k‖ is, the higher its frequency is. Given any x and its local function (i.e., f(·) around x), πf̂ (k) is the proportion/percentage that frequency k occupies. Therefore, the integral of πf̂ (k) · ‖k‖
2 reflects the contribution of high frequency terms in the local frequency distribution of a function. Remark 3. Consider f as a value function in reinforcement learning setting. Theorem 1 indicates that regions with large gradient norm can either have large absolute value function, or high local frequency, or both. To prevent finding regions that only have large negative value function, our theory implies that it is reasonable to take both gradient norm and value function into account, as our proposed method does in the next section.

heading: 4 FREQUENCY-BASED SEARCH-CONTROL IN DYNA
text: We present the Dyna architecture with the frequency-based search-control (Algorithm 1) in this section. It combines the idea that samples from high-frequency regions of the state space is important,
as discussed in the previous section, and the hill climbing process to effectively draw samples from those regions, as introduced by Pan et al. (2019). We omit implementation details such as preconditioning, noisy gradient for the hill climbing process, and refer readers to Appendix A.6 and A.7.
Our goal is to query the model more often from the regions of the state space where the local frequency of the value function is higher. The intuition behind this search-control mechanism, as discussed in the previous section in the context of supervised learning, is that those regions correspond to where learning the (value) function is more difficult, hence more samples from the model might be helpful. To populate the search-control queue with states from those regions, we can do hill climbing on g(s) = ‖∇sV (s)‖2 + ‖Hv(s)‖2F . Theorem 1, however, suggests that states with large gradient norm can either have large absolute value, or high local frequency, or both. We want to avoid many samples from regions with large negative value states, as those states may be rarely visited under the optimal policy anyway. A sensible strategy to get around this problem is to combine the proposed hill climbing method with the previous hill climbing on the value function (Pan et al., 2019), as the latter tends to generate samples from high value states.
We propose the following method for combining those approaches. At each time step, with certain probability, we perform hill climbing by either
s← s+ α { ∇sg(s) with probability of p (7a) ∇sV (s) with probability of 1− p (7b)
and store states along the gradient trajectory in the search-control queue. When hill climbing on the value function (7b), we sample the initial state from the ER buffer as suggested by the previous work (Pan et al., 2019). This populates the search-control queue with states from the high value regions of the state space. When hill climbing on g(s) (7a), however, we sample the initial state from the search-control queue itself (instead of the ER buffer). This way ensures that the initial state for searching high frequency region has relatively high value. Hill climbing on g(s) from an initial state with a high value populates the search-control queue with high frequency samples around high value regions of the state space. We discuss some other intuitive mechanisms that we have tested in Appendix A.4.
Similar to the previous work (Pan et al., 2019), we obtain the state-value function in both (7a) and (7b) by taking the maximum of the estimated action-value, i.e. V (s) = maxaQ(s, a) ≈ maxaQθ(s, a) where θ is the parameter of the Q-network. Similar to the Dyna architecture (Algorithm 2), during planning stage, we sample multiple mixed mini-batches to update the parameters (i.e. we call multiple planning steps/updates). The mixed mini-batch was also used in the work by Gu et al. (2016) and can alleviate off-policy sampling issue as studied by Pan et al. (2019).

heading: 5 EXPERIMENTS
text: In the experiments, we carefully study the properties of our algorithm on the MountainCar benchmark domain. Then we illustrate the utility of our algorithm on a challenging self-designed MazeGridWorld domain, by which we illustrate the practical implication of having samples from the high frequency regions. Though we mainly focuses on search-control instead of how to learn a model, we include the result of using an online learned model for our algorithm. We refer readers to Appendix A.5 for additional experiments and Appendix A.6 for the reproducibility detail.

heading: 5.1 UTILITY OF FREQUENCY-BASED SEARCH-CONTROL
text: The MountainCar (Brockman et al., 2016) domain is well-studied, and it is known that the value function under the optimal value function has sharp changing regions (Sutton & Barto, 2018), which is the setting where our algorithm should be more effective. The agent needs to learn to reach the goal state within as few steps as possible since the reward is −1 per time step. The purposes of experimenting on this domain are: 1) verify that our search-control can outperform several natural competitors under different number of planning updates; 2) show that our search-control is robust to environment noise.
We use the following competitors. Dyna-Frequency is the Dyna architecture using the proposed search-control strategy (Algorithm 1); Dyna-Value is Algorithm 3 from the previous work (Pan et al., 2019); PrioritizedER is DQN with prioritized experience replay (Schaul et al., 2016); ER
Algorithm 1 Dyna architecture with Frequency-based search-control B: the ER buffer, Bs: search-control queue M : S ×A → S × R, the model outputs the next-state and reward m: number of states we want to fetch by hill climbing, d: number of planning steps a: the threshold for accepting a state Q,Q′: current and target Q networks, respectively b: the mini-batch size, β ∈ (0, 1): the proportion of simulated samples in a mini-batch τ : update target network Q′ every τ updates to Q t← 0 is the time step, nτ is the number of parameter updates while true do
Observe st, take action at (i.e. -greedy w.r.t. Q) Observe st+1, rt+1, add (st, at, st+1, rt+1) to B // Gradient ascent hill climbing With probability p, 1− p, choose hill climbing rule (7a) or (7b) respectively; sample s from Bs if choose rule (7a), or from B otherwise; set c← 0, s̃← s while c < m do
update s by executing the chosen hill climbing rule if s is out of state space then: // resample the initial state and hill climbing rule
With probability p, 1− p, choose hill climbing rule (7a) or (7b) respectively; sample s from Bs if choose (7a), or from B otherwise; set c← 0, s̃← s continue
if ||s− s̃||2/ √ n > a then: // n is the state dimension, i.e. S ⊂ Rn
add s to Bs, s̃← s, c← c+ 1 for d times do // d planning updates: sample d mini-batches
draw βb sample states from the search-control queue Bs, pair them with their corresponding on-policy action, and queryM to get the corresponding next-states and rewards draw (1 − β)b sample transitions from the ER buffer B and add them to the simulated transitions
use the mixed mini-batch for parameter update of the estimator, e.g., DQN nτ ← nτ + 1 if mod (nτ , τ) == 0 then:
Q′ ← Q t← t+ 1
is simply DQN with experience replay (ER) (Mnih et al., 2015). Figure 3 shows the learning curves of all those algorithms using 10 planning updates (a)(b) and 30 planning updates (c)(d) under different stochasticity. In Figure 3(b)(d), the reward is sampled from the Gaussian distribution N(−1, σ2), σ ∈ {0.0, 0.1}. We make several important observations: 1) With increased number of planning updates, these algorithms do not necessarily perform better, as shown in Figure 3(c). The proposed algorithm, however, appears to gain more through more number of updates since the difference between DynaFrequency and Dyna-Value seems to be clearer in Figure 3(c) than in Figure 3(a). 2) Since both Dyna-Value and our algorithm fetch the same number of states (i.e. m = 20) by hill climbing, the superior performance of our algorithm indicates the advantage of using samples from the high frequency regions. 3) PrioritizedER clearly performs worse than our algorithm and Dyna-Value, which probably implies the utility of the generalization power of the value function to acquire additional samples. 4) Our algorithm maintains superior performance in the presence of noise. One reason is that, noisy perturbation leads to more “energy” in all frequencies. When taking derivative, those high frequency terms are amplified. Hence, even with perturbation, high frequency region remains while the value estimate itself may get affected in an unpredictable manner.

heading: 5.2 A CASE STUDY: MAZEGRIDWORLD
text: We now illustrate the utility of our method on a challenging MazeGridWorld domain as shown in Figure 4(a). The domain has continuous state space S = [0, 1]2 and four discrete actions {up,down, left, right}. There are three walls in the middle, each of which has a hole for the agent
to go through. Each episode starts from the bottom left and ends at top right and the agent receives a reward of −1 at each time step, hence the agent should learn to use as few steps as possible to reach the goal area. On this domain, we mainly study our algorithm and the Dyna-Value algorithm.
Figure 4(b) shows the evaluation curves of the two algorithms. An important difference between our algorithm and the previous work is in the variance of the evaluation curve, which implies a robust policy learned by our method. In Figure 5, we further investigate the state distribution in searchcontrol queues of the two algorithms by uniformly sampling 1000 states from the two queues. Notice that a very important difference between the two distributions is that our search-control queue has a clearly high density around the bottleneck area, i.e., the hole areas where the agent can go across the walls. Learning a stable policy around such areas is extremely important: the agent simply fails to reach the goal state if they cannot pass any one of the holes. This distinguishes our algorithm with the previous work, which appears to acquire states near the goal area.

heading: 6 DISCUSSION
text: We motivated and studied a new category of methods for search-control by considering the approximation difficulty of a function. We provided a method for identifying the high frequency regions of a function, and justified it theoretically. We conducted experiments to illustrate our theory. We incorporated the proposed method into the Dyna architecture and empirically investigated its benefits. The method achieved competitive learning performances on a difficult domain.
There are several promising future research directions. First, it is worth exploring the combination of different search-control strategies. Second, we can explore the use of active learning methods (Settles, 2010; Hanneke, 2014) in the design of search-control mechanisms, since active learning con-
cerns about learning a function with as few samples as possible. This directly corresponds to our main purpose of using a smart search-control method in Dyna: to improve policy by using as few planning steps as possible.

heading: ACKNOWLEDGMENTS
text: We would like to thank the anonymous reviewers for their helpful feedback. Amir-massoud Farahmand acknowledges the funding from the the Canada CIFAR AI Chairs program. Yangchen Pan acknowledges the funding from Amii.

heading: A APPENDIX
text: We provide calculations for Example 1 and Example 2 in Section A.1; theoretical proof of Theorem 1 in Section A.2. Background of Dyna is reviewed in Section A.3 and additional discussions regarding search-control are provided in Section A.4. Additional experiments regarding separate criterion of hill climbing, and on continuous control problems are shown in Section A.5. Experimental details for reproducing our empirical results are in Section A.6.

heading: A.1 CALCULATIONS FOR EXAMPLE 1 AND EXAMPLE 2
text: Example 1. For fsin defined in Eq. (1), calculate the integrals of squared first order derivative f ′sin on high frequency region [−2, 0) and low frequency region [0, 2], respectively:∫ 0
−2 [f ′sin(x)] 2 dx = 64π2, ∫ 2 0 [f ′sin(x)] 2 dx = π2.
Proof. Taking derivative and integral,∫ 0 −2 [f ′(x)] 2 dx = 64π2 ∫ 0 −2
[cos (8πx)] 2 dx = 64π2,∫ 2
0
[f ′(x)] 2 dx = π2 ∫ 2 0 [cos (πx)] 2 dx = π2.
Example 2. Let f : [−π, π]→ R be a band-limited real valued function defined as
f(x) = a0 2 + N∑ n=1 an cos (nx) + bn sin (nx),
where a0, an, bn ∈ R, n = 1, 2, . . . , N are Fourier coefficients of frequency n2π . Then,∫ π −π |f ′(x)|2 dx = π · N∑ n=1 n2 ( a2n + b 2 n ) , ∫ π −π |f ′′(x)|2 dx = π · N∑ n=1 n4 ( a2n + b 2 n ) .
Proof. Taking derivative of f ,
f ′(x) = N∑ n=1 [−nan sin (nx)] + N∑ n=1 [nbn cos (nx)].
Taking square of f ′,
[f ′(x)] 2 = N∑ n=1 N∑ m=1 [nmanam sin (nx) sin (mx)]− N∑ n=1 N∑ m=1 [nmanbm sin (nx) cos (mx)]
− N∑ n=1 N∑ m=1 [mnambn sin (mx) cos (nx)] + N∑ n=1 N∑ m=1 [nmbnbm cos (nx) cos (mx)].
Taking integral,∫ π −π [f ′(x)] 2 dx = ∫ π −π N∑ n=1 N∑ m=1 [nmanam sin (nx) sin (mx)]dx− ∫ π −π N∑ n=1 N∑ m=1 [nmanbm sin (nx) cos (mx)]dx
− ∫ π −π N∑ n=1 N∑ m=1 [mnambn sin (mx) cos (nx)]dx+ ∫ π −π N∑ n=1 N∑ m=1 [nmbnbm cos (nx) cos (mx)]dx
= N∑ n=1 N∑ m=1 [nmanamπδn,m − 0− 0 + nmbnbmπδn,m]
= π · N∑ n=1 n2 ( a2n + b 2 n ) ,
where
δn,m def = { 1, if n = m, 0, otherwise.
Using similar arguments, taking derivative of f ′(x),
f ′′(x) = N∑ n=1 [ −n2an cos (nx) ] + N∑ n=1 [ −n2bn sin (nx) ] .
Taking integral, ∫ π −π [f ′′(x)] 2 dx = π · N∑ n=1 n4 ( a2n + b 2 n ) .

heading: A.2 PROOF FOR THEOREM 1
text: Notations. For any vector norm || · ||, we mean l2 norm and we ignore the subscript unless clarification is needed. We use Frobenius norm || · ||F for matrix. We use subscript yl to denote the lth element in vector y. Let Hf (y) be the Hessian matrix of f(y). We write H for short unless clarification is needed. Let Hl,: be the lth row of the Hessian matrix.
Proof description. We establish the connection between local gradient norm, Hessian norm and local frequency. To build such connection, we introduce a definition of πf̂ as shown below and we call it “local frequency distribution” of f(x). πf̂ is a probability distribution over R
n, i.e.,∫ k∈Rn πf̂ (k)dk = 1, and πf̂ (k) ≥ 0, ∀k ∈ R
n. Within an open subset of domain (an unit ball), this distribution characterizes the proportion of a particular frequency component occupies. The proof can be described by three key steps:1) We use a local Fourier transform to express a function locally (i.e. within an unit ball). 2) we calculate the gradient/Hessian norm based on this local Fourier transform; 3) we take integration over the unit ball of the gradient/Hessian norm to build the connection with the local frequency distribution πf̂ and function energy.
Theorem 1. Given any function f : Rn → R, for any frequency vector k ∈ Rn, define its local Fourier transform of x ∈ Rn,
f̂(k) def = ∫ y∈B(x,1) f(y) exp { −2πi · y>k } dy,
for local function around x, i.e., y ∈ B(x, 1) def= {y : ‖y − x‖ < 1}. Assume the local function “energy” is finite, ∫
y∈B(x,1) [f(y)]
2 dy = ∫ Rn ‖f̂(k)‖2dk <∞, ∀x ∈ Rn.
Define “local frequency distribution” of f(x) as:
πf̂ (k) def = ‖f̂(k)‖2∫ Rn ‖f̂(k̃)‖2dk̃ , ∀k ∈ Rn.
Then, ∀x ∈ Rn, we have: 1) the first order connection:∫
y∈B(x,1) ‖∇f(y)‖2 dy = 4π2 · [∫ y∈B(x,1) [f(y)] 2 dy ] · [∫ Rn πf̂ (k) · ‖k‖ 2 dk ] ,
2) the second order connection:∫ y∈B(x,1) ||H(y)||2F dy = 16π4 [∫ y∈B(x,1) [f(y)] 2 dy ] · [∫ Rn πf̂ (k) · ‖k‖ 4 dk ]
Proof. 1) We first prove the first order connection. Consider the following function defined locally around x,
fx(y) def = { f(y), if y ∈ B(x, 1), 0, otherwise.
By definition, the Fourier transform of fx is
f̂(k) = ∫ y∈B(x,1) f(y) exp { −2πi · y>k } dy
= ∫ Rn fx(y) exp { −2πi · y>k } dy.
And the inverse Fourier transform of fx(y), ∀y ∈ B(x, 1) is,
fx(y) = ∫ Rn f̂(k) exp { 2πi · y>k } dk,
and then the gradient ∀y ∈ B(x, 1) is ∇f(y) = ∇fx(y) = ∫ Rn f̂(k) exp { 2πi · y>k } (2πi · k) dk. (8)
To calculate gradient norm, we use complex conjugate, ∇f∗(y) = ∫ Rn f̂∗(k′) exp { −2πi · y>k′ } (−2πi · k′) dk′,
where
f̂∗(k′) = ∫ Rn fx(y ′) exp { 2πi · y′>k′ } dy′
is the complex conjugate of f̂(k). Therefore,
‖∇f(y)‖2 = 〈∇f(y),∇f∗(y)〉
= ∫ Rn ∫ Rn f̂(k)f̂∗(k′) exp { 2πi · y> (k − k′) } ( 4π2k>k′ ) dkdk′.
(9)
Taking integral of ‖∇f(y)‖2 within the unit ball centered at x,∫ y∈B(x,1) ‖∇f(y)‖2 dy = ∫ Rn ‖∇fx(y)‖2 dy, by function definition (10a)
= ∫ Rn ∫ Rn f̂(k)f̂∗(k′) [∫ Rn exp { 2πi · y> (k − k′) } dy ] ( 4π2k>k′ ) dkdk′ (10b)
= ∫ Rn ∫ Rn f̂(k)f̂∗(k′)δk−k′,0 ( 4π2k>k′ ) dkdk′ (10c)
= 4π2 ∫ Rn ‖f̂(k)‖2 · ‖k‖2 dk. (10d)
Recall the definition of local function “energy” around x,∫ Rn ‖f̂(k)‖2dk = ∫ Rn 〈 f̂(k), f̂∗(k) 〉 dk (11a)
= ∫ y∈Rn ∫ y∈Rn fx(y)fx(y ′) [∫ Rn exp { 2πi · k> (y′ − y) } dk ] dydy′ (11b)
= ∫ y∈Rn ∫ y∈Rn fx(y)fx(y ′)δy′−y,0dydy ′ (11c)
= ∫ y∈Rn f2x(y)dy (11d)
= ∫ y∈B(x,1) f2(y)dy <∞. (by definition of fx(y) and finite energy assumption)
(11e)
For y ∈ B(x, 1), the local gradient information is related to local energy and frequency distribution,
∫ y∈B(x,1) ‖∇f(y)‖2 dy = 4π2 ∫ Rn ‖f̂(k)‖2 · ‖k‖2 ∫ Rn ‖f̂(k̃)‖ 2dk̃∫ Rn ‖f̂(k̃)‖2dk̃ dk (12a)
= 4π2 ∫ Rn πf̂ (k) ‖k‖ 2 ∫ Rn ‖f̂(k̃)‖2dk̃dk (12b)
= 4π2 · [∫ y∈B(x,1) f2(y)dy ] · [∫ Rn πf̂ (k) · ‖k‖ 2 dk ] , (12c)
where the last equality follows by ∫ Rn ‖f̂(k̃)‖ 2dk̃ = ∫ y∈B(x,1) f
2(y)dy which is established in the derivation (11).
2) Now we prove the second order connection. To show the second order connection, we start from Eq. (8). Then the lth row of the Hessian matrix Hl,: can be written as:
Hl,: = ∂∇f(y) ∂yl
>
where we use the notation ∂∇f(y)∂yl to denote the vector formed by taking partial derivative of each element in the gradient vector∇f(y) w.r.t. yl. Then,
∂∇f(y) ∂yl = ∫ Rn f̂(k) exp { 2πi · y>k } ( 4π2i2(e>l k)k ) dk,
where el is standard basis vector where the lth element is one. To calculate the norm of the vector Hl,: =
∂∇f(y) ∂yl
> , we use complex conjugate again and follow the similar derivation as done
in Eq. (9):
||Hl,:||22 = 〈Hl,:, Hl,:〉
= ∫ Rn ∫ Rn f̂(k)f̂∗(k′) exp { 2πi · y>(k − k′) } ( 16π4i4(e>l k)(e > l k ′)k>k′ ) dkdk′
Note that the square of Frobenius norm of the Hessian matrix can be written as ||H||2F =∑ i,j H 2 i,j = ∑n l=1 ||Hl,:||22. Then,
||H||2F = n∑ l=1 ||Hl,:||22
= ∫ Rn ∫ Rn f̂(k)f̂∗(k′) exp { 2πi · y>(k − k′) }( 16π4i4 n∑ l=1 (e>l k)(e > l k ′)k>k′ ) dkdk′
= 16π4 ∫ Rn ∫ Rn f̂(k)f̂∗(k′) exp { 2πi · y>(k − k′) } (k>k′)2dkdk′
Taking the integration of ||H||2F over y variable within a ball with center x and unit radius, we acquire: ∫
y∈B(x,1) ||H(y)||2F dy = 16π4 ∫ Rn ||f̂(k)||2||k||4dk
= 16π4 [∫ y∈B(x,1) [f(y)] 2 dy ] · [∫ Rn πf̂ (k) · ‖k‖ 4 dk ] where the derivation process for the first equation is a simple modification from the derivation (10) and the second equation follows the same derivation (12).
Algorithm 2 Generic Dyna Architecture: Tabular Setting Initialize Q(s, a) and modelM(s, a), ∀(s, a) ∈ S ×A while true do
observe s, take action a by -greedy w.r.t Q(s, ·) execute a, observe reward R and next state s′ Q-learning update for Q(s, a) update modelM(s, a) (i.e. by counting) store (s, a) into search-control queue for i=1:d do
sample (s̃, ã) from search-control queue (s̃′, R̃)←M(s̃, ã) // simulated transition Q-learning update for Q(s̃, ã) // planning update
Discussion with Uncertainty Principle. We now provide an intuitive interpretation of our theorem from the well-known uncertainty principle. The Uncertainty Principle says that a function cannot be simultaneously concentrated in both spatial and frequency space. That is, the more concentrated a function is, the more spread out its Fourier transform must be, indicating that more high frequency terms are needed to express the function. For example, by uncertainty principle (Stein & Shakarchi, 2003), ∀x ∈ R we have[∫
‖y−x‖≤1 (y − x)2 · [fx(y)]2 dy
] · [∫
Rn ‖f̂(k)‖2 · ‖k‖2 dk
] ≥ 1
16π2 . (13)
Note that the first term [∫ ‖y−x‖≤1 (y − x) 2 · [fx(y)]2 dy ]
measures the dispersion of the function around the point x. Hence the smaller this term is, the more concentration of the function is on the specified domain. Plug into our Eq. (5) and replace the function energy term by Fourier transform and cancel out the normalization term, we get:[∫
‖y−x‖≤1 (y − x)2 · [fx(y)]2 dy
] · [∫ ‖y−x‖≤1 ‖∇f(y)‖2 dy ] ≥ 1 4 ,
which means the more concentrated f is locally around x, the larger the lower bound of the local gradient norm must be.

heading: A.3 BACKGROUND IN DYNA
text: In this section, we provide the vanilla (tabular) Dyna (Sutton, 1991a; Sutton & Barto, 2018) in Algorithm 2, and Hill Climbing Dyna by Pan et al. (2019) in Algorithm 3. Dyna is a classic modelbased reinforcement learning architecture. As described in Algorithm 2, at each time step, the real experience is used to directly improve policy/value estimates, and is also used to learn the environment dynamics model. During planning stage, simulated experiences are acquired from the learned model and are used to further improve the policy. The critical component in Dyna is the search-control mechanism, which decides what simulated experiences to use during planning. This area is relatively unexplored, though abundant literature is available regarding how to learn a model.

heading: A.4 A DISCUSSION ON SEARCH-CONTROL DESIGN BASED ON HILL CLIMBING
text: There are different ways to combine hill climbing strategies. Here are some unsuccessful trials. For example, climbing on direct combinations of V (s) (value function) and g(s) (frequency criterion), such as V (s) + g(s), or V (s)g(s), did not work well. The reasons are as following. First, such combination can lead to unpredictable gradient behaviour. It can alter the trajectory solely based on either g(s) or V (s), and the effect is unclear. It may lead to states with neither high value or high frequency. Last, and probably the most important, hill climbing on V (s) and on g(s) have fundamentally different insights. The former is based on the intuition that the value information should be propagated from the high value region to low value region; as a result, it requires to store states along the whole trajectory, including those in low value region. This is empirically verified by Pan et al. (2019). However, the latter is based on the insight that the function value
Algorithm 3 HC-Dyna architecture Bs: search-control queue, B: the experience replay buffer m: number of states to fetch by search-control b: the mini-batch size while true do
Observe st, take action at (i.e. -greedy w.r.t. action value function) Observe st+1, rt+1, add (st, at, st+1, rt+1) to B sample s from visited states, i.e. ER buffer B // Hill climbing by gradient ascent while get less than m states do
s← s+∇sV (s), V (s) = maxsQθ(s, a) store s into search control queue Bs
// Planning stage for d times do
// sample states from Bs and pair them with on-policy actions, query the model to get next states and rewards // mix simulated and real experiences into a mini-batch and use it to update parameters (i.e. DQN update)
t← t+ 1
in high frequency region is more difficult to approximate and needs more samples, while there is no obvious reason to propagate those information back to low frequency region. As a result, this approach does not emphasize on recording states throughout the whole hill climbing trajectory.

heading: A.5 ADDITIONAL EXPERIMENTS
text: In this section, we briefly study the effect of doing hill climbing on only gradient norm or Hessian norm. Then we demonstrate that our search-control strategy can be also used for continuous control algorithms.
Hill climbing on only gradient norm or Hessian norm. Throughout our paper, we use the form of g(s) = ‖∇sV (s)‖2 + ‖Hv(s)‖2F to search states from high (local) frequency region of the value function. Besides the theoretical reason, there is a practical demand of such design. On value function surface, regions which have low (or even zero) gradient magnitude may have high Hessian magnitude, and vice versa. Hence, it can help move along the gradient trajectory in case that one of the term vanished at some point. Such cases can be a result of function approximation (smoothness/differentiability), or of the nature of the task, or both. In Fig. 6, we show the results of using only either gradient norm or Hessian norm. The reason we choose MountainCar and GridWorld (the same domain as described by Pan et al. (2019)) is that, the former has a value function surface with lots of variations; while the latter’s value function increases smoothly from the initial state to the goal state, which indicates a small magnitude second-order derivative. Indeed, we empirically observe that the term ∇s‖Hv(s)‖2F frequently gives a zero vector on GridWorld. This explains the bad performance of Dyna-HessNorm in Fig. 6(b). In contrast, Fig. 6(a) shows slightly better performance of Dyna-HessNorm and Dyna-GradNorm. Notice that, an intuitive and more general form of g(x) can be g(s) = η1‖∇sV (s)‖2 + η2‖Hv(s)‖2F , at the cost that additional meta-parameters are introduced.
Continuous Control. In this section, we show a simple demonstration where our method is adapted to two continuous control tasks: Hopper-v2 and Walker2d-v2 from Mujoco (Todorov et al., 2012) by using a continuous Q learning algorithm called NAF (Normalized Advantage Function) (Gu et al., 2016). The algorithm parameterizes the action value function as Q(s, a) = V (s) − (a − µ(s))TP (a − µ(s)) where P is a positive semi-definite matrix and hence the action with maximum value can be easily found: arg maxaQ(s, a) = µ(s). Our search-control strategy naturally applies here by utilizing the value function V (s). From Fig. 7, one can see that our algorithm (DynaNAF-Frequency) finds a better policy comparing with the model-free NAF.

heading: A.6 REPRODUCIBILITY DETAIL
text: All of our implementations are based on tensorflow with version 1.13.0 (Abadi et al., 2015). For DQN update, we use Adam optimizer (Kingma & Ba, 2014). We use mini-batch size b = 32 except on the supervised learning experiment where we use 128. For reinforcement learning experiment, we use buffer size 100k. All activation functions are tanh except the output layer of the Q-value is linear. Except the output layer parameters which were initialized from a uniform distribution [−0.003, 0.003], all other parameters are initialized using Xavier initialization (Glorot & Bengio, 2010). For model learning, we use a 64 × 64 relu units neural network to predict s′ − s given a state-action pair with mini-batch size 128 and learning rate 0.0001.
For the supervised learning experiment shown in Section 3, we use 16 × 16 tanh units neural network, with learning rate 0.001 for all algorithms. The learning curve is plotted by computing the testing error every 20 iterations. When generating Fig. 2, in order to sample points according to p(x) ∝ |f ′(x)| or p(x) ∝ |f ′′(x)|, we use 10, 000 even spaced points on the domain [−2, 2] and the probabilities are computed by normalization across the 10k points.
The experiment on MountainCar is based on the implementation from OpenAI (Brockman et al., 2016), we use 32 × 32 tanh layer, with target network moving rate 1000 and learning rate 0.001. Exploration noise is 0.1 without decaying. For all algorithms, we use warm up steps = 5000 (i.e. random action is taken in the first 5k time steps). Prioritized experience replay (PrioritizedER) is implemented as the proportional version with sum tree data structure. We use prioritized ER without importance ratio but half of mini-batch samples are uniformly sampled from the buffer as a strategy for bias correction. For Dyna-Value and Dyna-Frequency, we use: gradient ascent step size (in
search-control) 0.01, mixing rate β = 0.5 and m = 20, i.e., at each environment time step we fetch 20 states by hill climbing. We fix p = 0.5 across all experiments, hence the hill climbing rules (7a) and (7b) are selected with equal probability. We use natural projected gradient ascent for hill climbing as introduced by Pan et al. (2019).
For the experiment on MazeGridWorld, each wall’s width is 0.1 and each hole has height 0.1. The left-top point of the hole in the first wall (counting from left to right) has coordinate (0.2, 0.5); the hole in the second wall has coordinate (0.4, 1.0) and the third one is 0.7, 0.2. Each action leads to 0.05 unit move perturbed by a Gaussian noise from N(0, 0.01). On this domain, for both DynaValue and Dyna-Frequency, all parameters are set the same with that used on MountainCar except that we use 64 × 64 tanh units for Q network, and number of search-control samples is set as m = 50, number of planning updates is 30. As a supplement to the Section 5.2, we also provide the state distribution from ER buffer in Figure 8. One can see that ER buffer has very different state distribution with search-control queue.

heading: A.7 ALGORITHMIC DETAILS
text: We provide the pseudo-code in Algorithm 4 with sufficient details to recreate our experimental results. The hill climbing rules we used is the same as introduced by Pan et al. (2019). Define
vs def = ∇s max
a Q(s, a),
then
gs def = ∇sg(s) = ∇s(||∇s max a Q(s, a)||22 + ||Hv(s)||2F ) = ∇s(||vs||22 + ||∇svs||2F ).
Note that we use a squared norm to ensure numerical stability when taking gradient. Then for value-based search-control, we use
s← s+ α ||Σ̂svs|| Σ̂svs +Xi, Xi∼ N(0, ηΣ̂s) (14)
and for frequency-based search-control, we use
s← s+ α ||Σ̂sgs|| Σ̂sgs +Xi, Xi∼ N(0, ηΣ̂s) (15)
where Σ̂s is empirical covariance matrix estimated from visited states, and we set η = 0.01, α = 0.01 across all experiments. Notice that comparing with the previous work, we omitted the projection step as we found it is unnecessary in our experiments.
Algorithm 4 Dyna architecture with Frequency-based search-control with additional details Bs: search-control queue, B: the experience replay buffer M : S ×A → S × R, the environment model m: number of search-control samples to fetch at each step p: probability of choosing value-based hill climbing rule (we set p = 0.5 for all experiments) β ∈ [0, 1]: mixing factor in a mini-batch, i.e. βb samples in a mini-batch are simulated from model n: number of state variables, i.e. S ⊂ Rn a: empirically learned threshold as sample average of ||st+1 − st||2/ √ n
d: number of planning steps Q,Q′: current and target Q networks, respectively b: the mini-batch size τ : update target network Q′ every τ updates to Q t← 0 is the time step nτ ← 0 is the number of parameter updates // Gradient ascent hill climbing With probability p, 1− p, choose hill climbing Eq. (14) o Eq. (15) respectively; sample s from Bs if choose rule Eq. (14), or from B otherwise; set c← 0, s̃← s while c < m do
update s by executing the chosen hill climbing rule if s is out of state space then: // resample the initial state and hill climbing rule
With probability p, 1− p, choose hill climbing rule Eq. (14) or Eq. (15) respectively; sample s from Bs if choose Eq. (7), or from B otherwise; set c← 0, s̃← s continue
if ||s− s̃||2/ √ n > a then:
add s to Bs, s̃← s, c← c+ 1 // d planning updates: sample d mini-batches for d times do // d planning updates
sample βb states from Bs and pair them with on-policy actions, and query M to get next states and rewards
sample b(1− β) transitions from B an stack these with the simulated transitions use the mixed mini-batch for parameter (i.e. DQN) update nτ ← nτ + 1 if mod(nτ , τ) == 0 then:
Q′ ← Q t← t+ 1

Review summaries:
This paper basically built upon [1]. The authors propose to do sampling in the high-frequency domain to increase the sample efficiency.  They first argue that the high-frequency part of the function is hard to approximate (i.e., needs more sample points) in section 3.1. They argue that the gradient and Hessian can be used to identify the high-frequency region. And then they propose to use g(x)=||gradient||+||Hessian || as the sampling metric as illustrated in Algorithm 1. To be noticed that, they actually hybrid the proposed metric (6) and the value-based metric (7, proposed in [1]) in their algorithm.

This paper proposes a new way to select states from which do do transitions in dyna algorithm (which trains policy from model experience as if it was a real experience). It proposes to look for states where frequency of value function as a function of a real valued state is large, because these are the states where the function is harder to approximate. The paper also shows that such frequency is large where the gradient of the function is large in magnitude which allows for finding such states in practice. In more detail, similar to previous algorithms, this algorithm keeps both an experience replay buffer as well as another buffer of states (search-control queue) and uses a hill climbing strategy to find states with both higher frequency and higher value. The paper tests the algorithm on toy domains - the mountain car and a maze with doors.