Interpretable Function Approximation with Gaussian Processes in Value-Based Model-Free Reinforcement Learning

Venue: NLDL_org_2025_Conference
Published: 2024-08-23 15:48:45

Abstract: Estimating value functions in Reinforcement Learning (RL) for continuous spaces is challenging. While traditional function approximators, such as linear models, offer interpretability, they are limited in their complexity. In contrast, deep neural networks can model more complex functions but are less interpretable. Gaussian Process (GP) models bridge this gap by offering interpretable uncertainty estimates while modeling complex nonlinear functions. This work introduces a Bayesian nonparametric framework using GPs, including Sparse Variational (SVGP) and Deep GPs (DGP), for off-policy and on-policy learning. Results on popular classic control environments show that SVGPs/DGPs outperform linear models but converge slower than their neural network counterparts. Nevertheless, they do provide valuable insights when it comes to uncertainty estimation and interpretability.

Main paper: 
 model architectures: 
 The computational backend for all the models is PyTorch [26] . For the GP models, we used GPytorch [21] , a high-performance GPU-accelerated library for GP modelling, in conjunction with Botorch [27] , a Bayesian optimization library that extendsGPytorch. An NVIDIA RTX 4090 with 24 GB of VRAM was used for GPU acceleration.We were also interested in measuring VRAM usage on the GPU, energy consumption, and execution time. For VRAM usage, we used pynvml, a Python interface for the NVIDIA Management Library. For energy consumption and execution time, we used the Zeus library [28] .For the choice of GP model, we primarily considered DGPs, which reduce to an SVGP when using a single unit.In Section 2.2, the derivations assumed an observation noise variance parameter σ 2 y . In the RL setting, it is hard to estimate the exact value of σ 2 y .To address this, instead of setting σ 2 y a priori, we inferred it alongside the hyperparameters by including σ 2 y in the optimization process when minimizing the negative MLL in (18) or the negation of the ELBO in (23) .Regarding the choice of the kernel for the states, existing literature on GPs and RL have used the RBF (16) or Matern (15) kernel [9, 14] . This implies that those authors had a prior belief that the actionvalue function is reasonably smooth. We assumed that this assumption was reasonable. Among the two, the RBF kernel was used.Regarding [8] 's suggestion to use a separate kernel for the states and action is sensible, but this causes issues with sparse variational methods when learning the inducing points Z ∈ R M d using gradient methods [7] . The actions are discrete, but the problem is that the optimizer that optimizes the inducing points z i ∈ S × A does not take this information into account. Using a kernel for categorical features introduces issues since the actions may not be exact integers anymore. For this reason, we used the approach by [9, 14] , and simply used one kernel for the states and actions.Using a constant mean function in DGPs makes each GP mapping highly non-injective, leading to issues with the DGP prior [29] . Following [15] , we Note the difference in how the q-function is mod-1014 elled compared to a GP. With a neural network, 1015 we have as input the state and an output neuron 1016 for each action, containing q(s, a). This is differ-1017 ent for a GP, where the input is a state-action pair 1018 z ∈ S × A, and the output is a Gaussian with a 1019 mean and variance for the action-value. For GP-Q/GP-SARSA, due to the experimental 1030 nature of the algorithm and the computational com-1031 plexity of certain hyperparameter settings, informal 1032 testing was performed on both CartPole and Lunar 1033 Lander to identify effective hyperparameter settings. 1034 This involved examining the reward curve on runs 1035 with different configurations. We tuned GP fitting 1036 settings, such as the initial learning rate of the Adam 1037 optimizer and how optimization over the dataset was 1038 performed at each GP update.The hyperparameters used in the experimental 1040 setup are summarized in Tables D. 3 and D.2. A 1041 moderate dataset size of 10,000 or 20,000, with each 1042 GP update involving a random subset of approxi-1043 mately 3,500 samples, works well when optimized 1044 in mini-batches. Using a small learning rate (0.001 1045 or 0.005) for the Adam optimizer ensures that hy-1046 perparameters are not too tightly fitted on early 1047 inaccurate q-value estimates. However, the learning 1048 rate should not be too low to allow the agent to 1049 learn a meaningful policy. The discount factor γ 1050 was set to 0.99, the same as with the DQN algo-1051 rithm. The selected behavioural policy chosen was 1052 Thompson sampling, but a small comparison was 1053 made to ϵ-greedy and UCB in Lunar Lander.What was tested more rigorously was the perfor-1055 mance difference going from an SVGP to a DGP. 1056As GP-Q and GP-SARSA are the same, except for the computation of the TD target, such measurements were omitted from GP-SARSA for brevity. 
gaussian processes: 
 A GP is a collection of RVs {f GP (x) | x ∈ X }, any finite number of which has a joint Gaussian distribution [6] . In the context of regression, the index set X is related to the input of some function f : X → R that we want to approximate, e.g., the action-value function. It can be interpreted as: At each point x ∈ X , the output of the GP regression model is an RV denoted f GP (x).A GP, denoted f GP (x) ∼ GP(m(x), k(x, x ′ )), is fully specified by a mean function m : X → R and covariance or kernel function k : X × X → R which are defined as:m(x) = E[f GP (x)], k(x, x ′ ) = E[(f GP (x) -m(x))(f GP (x ′ ) -m(x ′ ))].(5) A GP offers a Bayesian approach to nonparametric regression. Without any data, the kernel function represents our prior belief about the function we are trying to model, by abuse of notation denoted p(f ), as it encodes similarity between data points, with closer points having higher covariance. A GP can then be conditioned on a dataset, D, to get a posterior GP f GP (x) ∼ GP(m post (x), k post (x, x ′ )), which is our posterior belief about the function, p(f |D).The Appendix B provides a comprehensive discussion of GPs, and their extensions: Sparse Variational (SVGP) and Deep GPs (DGP).Using GPs for value-based RL is done by capturing a distribution over possible q-functions, q π ∼ GP(m(z), k(z, z ′ )), where the input domain is z ∈ S × A [9] . The core idea behind the GP-Q algorithm (see Algorithm 1) is to perform a type of Bayesian optimization on the TD error by setting the target values as described in Algorithm B.1 to the TD(0) target. Using different acquisition functions, we now have more options in designing the behavioural policy other than just using ϵ-greedy. [9] for example used a variant of GP-Q using UCB.The GP-Q algorithm is off-policy and usesEQUATIONas the TD target, where q is the mean output of the GP [9] . An on-policy variant, GP-SARSA, can be created by setting the target toEQUATIONcertainty [7] . This means that the predictive vari- of that particular part of the state space [16] .This has two key applications. Choose A t from S t using behavioral policy (e.g., ϵ-greedy)Take action A t , observe R t+1 , S t+1 6: let Z t = (S t , A t ) ⊤ and y t = R t+1 + γ max a ′ q(S t+1 , a ′ ), where q is the posterior mean functionAdd state-action pair to dataset, D t+1 = D t ∪ {(Z t , y t )} 8: if |D t+1 | > Budget 
training environments: 
 Since we are comparing value-based algorithms that involve taking an (arg-)max over the action space, we are constrained to environments with discrete action spaces, however, the state spaces are continuous as we are considering the realm of approximate RL. For our experiments, we rely on two popular classic control environments vastly used within the literature: CartPole and Lunar Lander from the Gymnasium library [18] . Note that when it comes to the state representation of the environments we normalize the state vector, by relying on standardizationEQUATIONwhere x and σ 2x are the mean and variance of x and η is a small constant. One exception is that for the Lunar Lander environment, we only normalized the first 6 dimensions as some state components are booleans.The GP-Q [9] and GP-SARSA [8] algorithms were adopted for the off-policy and on-policy cases respectively with the following modifications: The updates were done in batches, instead of per timestep. Meaning, that for a batch size B, we wait B timesteps before updating the GP. This was done to reduce computational complexity and improve sampling efficiency. Furthermore, this aimed to prevent the inaccurate estimate of the q-function from changing too much per timestep. All the modifications made to the base GP-Q/GP-SARSA algorithm are summarized in algorithm 2.Initially, D 0 = ∅. Updating the GP, i.e., computing p(q|D n+1 ), should be interpreted as performing type-II Maximum Likelihood Estimation (MLE) of the hyperparameters and observation noise variance θ ′ = [θ, σ 2 y ] ⊤ , and then computing the predictive posterior distribution where the test inputs are the state-action pairs of the current batch. In the case of SVGP/DPGs, we also have the Evidence Lower Bound (ELBO) with variational parameters Z, m, S as hyperparameters. We utilized the Adam optimizer [19] to estimate θ ′ opt by minimizing the negative Marginal Log Likelihood (MLL) (18) or the negative ELBO (23) . It was implicitly assumed that any hyperparameter values after optimization carry over to the next GP update. As explained in Appendix B.4, SVGPs allow for mini-batching in ELBO computation, which helps manage memory usage and enables handling significantly larger datasets. Two viable approaches to mini-batching in terms of time complexity are: using mini-batches over the entire dataset while assuming a relatively small dataset and hoping the variational parameters retain information from previously discarded Algorithm 2 Adjusted GP-Q/GP-SARSA 1: Initialize:2: Collect initial dataset D 0 = (z i , y i ) i=1,...,n0 with z i ∈ S × A, y i ∈ R 3:Initialize GP with p(q|D 0 ), noise variance σ 2 y , and hyperparameters θ 4: Choose A t from S t using a behavioral policy (e.g., ϵ-greedy, UCB, Thompson sampling) if Off-Policy then 10:EQUATIONy t = R t+1 + γ max a ′ q(S t+1 , a ′ ) 11:else {On-Policy}A t+1 = π(S t+1 ) 13:y t = R t+1 + γ q(S t+1 , A t+1 ) 14: end if 15:Update dataset:16: Form state-action pair Z t = [S t , A t ] ⊤ 17: D t+1 = D t ∪ {(Z t , y t )} 18: if |D t+1 | > Budget then 19: Delete oldest (z i , y i ) from D t+1 20: end if 21: if b == batch size then 22:Update GP model, p(q|D n+1 ):Minimize negative MLL in (18) or negative ELBO in ( 23 action from the sample q-function for a fixed state s ∈ S. This was done by considering the points {(s, a) | a ∈ A} sampling from the GP latent RVs, {f GP (s, a) | a ∈ A}, and selecting the action a for for which the sampled q-value is highest.EQUATIONWe compared GP-Q with DQN using both a Multi-Layer Perceptron (MLP) and a linear model, as well as a random policy during both training and evaluation. In Lunar Lander, we also included GP-SARSA. Since GP-Q and GP-SARSA were expected to perform similarly, we focused mainly on GP-Q. All algorithms were trained for 1000 episodes on CartPole and 3000 episodes on Lunar Lander with evaluation over 30 episodes. The episode numbers were chosen such that all algorithms would converge. This comparison involved examining the reward curve/return graph for policy convergence and performance, including the maximum return obtained, and measuring computational complexity and energy usage during training and evaluation1 . Time complexity was measured by considering execution time during training and evaluation. Space complexity was measured by examining VRAM usage during training. Energy usage was measured in kilojoules (kJ) during training and evaluation. Additionally, we measured the average negative ELBO for each environment over the number of GP updates.At evaluation time, we also compared the means of returns. A two-sample one-tailed Wilcoxon ranksum test was performed to compare the average return between a specific RL agent and a random policy after training. The null hypothesis wasH 0 : E πagent [G] = E π random [G]A : E πagent [G] > E π random [G],where G is the undiscounted sum of rewards over an episode.We also measured the stability of the GPQ with a DGP by examining how many times out of C r runs of at least C e episodes the algorithms performed better than random. For CartPole C r = 10, C e = 200 and Lunar Lander C r = 7, C e = 1000. If an algorithm is stable, we expect the same hyperparameter configuration and environment to produce consistent results w.r.t. random in most runs.Finally, we recorded the agent environment interaction after training and described the general behaviour in cases where it was insightful. For details about model architecture, experimental design, and hyperparameter tuning, refer to Appendix D. 
results: 
 The following showcases the results from the experiment described in Section 3. It takes around 1000 more episodes for GP-Q/GP-SARSA to reach a close-to-optimal policy in Lunar Lander compared to CartPole, which may be due to the larger state space. This raises questions about how well GP-Q and GP-SARSA converge as the state and action spaces grow larger. However, benchmarks from [15] show that DGPs can handle high-dimensional data, as demonstrated by their application to a 784-dimensional image dataset.Lunar Lander 393In the case of Lunar Lander, we visualized the GP 394 posterior distribution by fixing a state s ∈ S and 395 plotting a Gaussian with mean q(s, a) and variance 396 σ 2 (s, a) for each a ∈ A. The work aimed to compare a GP-based RL al-407 gorithm with linear and neural network function 408 approximators, specifically against DQN using lin-409 ear and MLP models. The goal was to assess the 410 strengths and limitations of GPs in RL.Findings from simulations align with the hypothe-412 sis: GP-based algorithms (GP-Q/GP-SARSA), par-413 ticularly using SVGPs and DGPs, outperform linear 414 function approximation in CartPole and Lunar Lan-415 der. However, they do not match the stability of 416 DQN with an MLP in CartPole or the overall per-417 formance in Lunar Lander. Additionally, GP-Q and 418 GP-SARSA are more computationally expensive, 419 even during inference. Given a train set of noisy observations D = (x i , y i ) i=1,...,N , where x i ∈ X and y i = f (x i ) + ϵ i with ϵ i ∼ N (0, σ 2 y ). Suppose we are interested in get-ting predictions f * = [f GP (x * 1 ), . . . , f GP (x * N * )] ⊤2 for test inputs X * = (x * 1 , . . . , x * N * ). By the definitions of a GP, it follows that the prior joint distribution p(y, f * |X, X * ) has the following form:EQUATIONwhereµ X = [(m(x 1 ), ..., m(x N )] ⊤ , µ * = [m(x * 1 ), . . . , m(x * N * )] ⊤ , K X,X = k(X, X) ∈ R N ×N , K σ = K X,X + σ 2 y I, K X, * = k(X, , X * ) ∈ R N ×N * ,and K * , * = k(X * , X * ) ∈ R N * ×N * are matrices of all the covariances between relevant datapoints.One can then condition on the observations to get the Bayesian predictive distribution for the test 2 f * is implicitly assumed to be a realization of the random vector.points. By Gaussian identities, we can get this in 696 closed form:EQUATION) 700Assuming a zero mean function, this reduces for a 701 single test point to:702EQUATION) 703 where k * = [k(x * , x 1 ), . . . , k(x * , x N )] ⊤ .To get the posterior mean and covariance function 705 we consider ( 12) and ( 13) over an infinite number 706 of potential test points:EQUATION)708where k* ′ = [k(x ′ * , x 1 ), . . . , k(x ′ * , x N )] ⊤ . 709 B.2 Kernel Learning 710The generalization properties of GPs rely on the 711 selection of the appropriate kernel [6, 7] . A common 712 kernel is the Matern kernel, given by 713 15) 714 with d = (x-x ′ ) ⊤ l -2 (x-x ′ ), lengthscale parameter 715 l, outputscale parameter σ f , smoothness parameter 716 ν, Gamma function Γ and a modified Bessel function 717 K ν . As ν → ∞, the Matern kernel approaches the 718 Radial Basis Function (RBF) kernel:EQUATIONEQUATION)720which when used results in smooth infinitely differ-721 entiable functions being sampled.For generalization, it is also important to optimize 723 the GP hyperparameters, such as kernel parameters 724 and noise variance, alongside computing the predic-725 tive distribution as seen in ( 12). This can be done 726 by performing type-II MLE through maximizing the 727 MLL log p(y|X, θ), for hyperparameters θ:EQUATION) 729As y ∼ N (0, K σ ) this integral can be computed as: where the first term is a data fit term, the second 732 term a model fit term and the last a constant. The 733 negative MLL is differentiable with respect to θ, so 734 stochastic gradient descent can be used.EQUATIONUsing black box matrix-matrix multiplication 736 (BBMM) inference [21] , which leverages GPU accel-737 eration, the overall time and space complexity of GP Upper confidence bound (UCB) is an acquisition function defined as:y n+1 = f (x n+1 ) + ϵ n 6: D n+1 = D n ∪ {(x n+1 , y n+1 )}EQUATIONwhere µ n , σ n are the mean and standard deviation outputs as described in (13) . β n is an exploration parameter.Another acquisition function is Thompson sampling. In the context of Bandits, this involves sampling an action-value function q from the GP posterior predictive distribution, and greedily selecting an action according to the sample:EQUATION)The intuition behind Thompson sampling is that exploration is encouraged by maximizing q because the sampled function is within the credible interval (standard deviations around the mean) with high values around the areas with high uncertainty. Maximizing q involves selecting actions where there is potentially high uncertainty on q. At the same time, actions with a high mean value are also likely to be sampled, which promotes exploitation.To address the O(N 2 ) time and space complexity of exact inference, where N = |D|, different approximation approaches can be taken to allow GPs to scale to larger datasets. This is especially of concern in RL, where in theory we have a continually growing dataset. See [7] for a full overview of the available techniques.SVGPs [7, 12, 13, 20, 23, 24] approximate the GP posterior predictive distribution through variational inference. The core idea is to use a set of inducing points Z = (z 1 , . . . , z M ) where M ≪ N , which serve as a sparse approximation of the full dataset.The associated inducing variables are denoted u = [f GP (z 1 ), . . . , f GP (z M )] ⊤ . The variational posterior is defined as:EQUATIONwhere p(f |u, X, Z) is the conditional density of the function values f = f GP (X) given train inputs X, inducing points Z and inducing variables u. p(u) is a Gaussian distribution with mean m and covariance matrix S.Variational inference aims to minimize the Kullback-Leibler (KL) divergence D KL (p ∥ p) between the variational posterior p and the true posterior p. Using the variational posterior in (21) and the lower bound on the marginal likelihood, data and also some vector datasets [24] . DGPs seek 857 to address this issue, while still staying in a Bayesian 858 nonparametric framework.EQUATIONEQUATIONA DGP is a composition of GPs [7, 25] :  
: 
 because they can learn a broad range of functions and extract features from data, as shown in tasks like pixel-based control from Atari images. [3] .A popular choice for interpretable function approximation is a Gaussian Process (GP), a nonparametric kernel-based model for regression or classification. Unlike linear models, they can model complex nonlinear functions [6, 7] and their nonparametric nature allows them to scale in complexity with the dataset size. Additionally, GPs provide uncertainty quantification out of the box, as the predictive variance is directly derived from the model, offering a way to estimate uncertainty in predictions. Engel et al. [8] employed GPs to approximate the value and action-value function, resulting in a variant of the SARSA algorithm for on-policy RL. This approach was extended to the off-policy case by Chowdhary et al. [9] , who proposed a variant of Q-learning using GPs alongside a proof of convergence. Building on these ideas, Grande et al. [10] derived sample complexity results for GPs in RL and introduced Delayed GPQ, a sample-efficient, model-free RL algorithm using GPs, which can obtain an optimal policy in a polynomial number of exploration steps in continuous state spaces. These methods relied on a sparsification method from Csató and Opper [11] to reduce computational complexity. However, they did not make use of more recent inducing point methods initially introduced by Titsias [12] , which offer a more effective and scalable sparsification approach that more accurately approximates the GP model on the full dataset [13] . Kameda and Tanaka [14] applied such inducing point methods to reduce the computational complexity in GP Q-learning.Previous research on GPs in value function approximation mainly explored strategies like upper confidence bound and ϵ-greedy, without considering alternatives like Thompson sampling, which could improve exploration effectiveness. Moreover, extensions to the basic GP model, such as Deep Gaussian Processes (DGPs) or Sparse Variational Gaussian Processes (SVGPs), have not been considered despite their ability to capture complex relationships through their hierarchical structure. Additionally, prior work, such as that of Kameda and Tanaka [14] , required pre-specified inducing points, which limits the algorithm's applicability in environments where selecting inducing points in advance is challenging. In contrast, the proposed algorithm is more gen-The behaviour of the RL agent is modelled by its 125 policy, which is a conditional probability distribution 126 over actions given the state: 127 π(a|s) = P(A t = a|S t = s) ∀s, a ∈ S × A, (1) 128 where S t ∈ S, A t ∈ A are Random Variables 129 (RVs) for the state and action at time t. The goal 130 is to learn an optimal policy π * that maximizes the 131 expected discounted sum of rewards,EQUATION)133where G is the return, R is the reward and γ ∈ (0, 1] 134 is a discount factor. The hypothesis in RL is that 135 by interacting with the environment the agent can 136 learn to make (sub-)optimal decisions. Value-based 137 methods achieve it by learning a value function:EQUATION139 or an action-value function:EQUATION141 and then deriving a policy from it. We want then an approximation v or q that approaches the optimal value function and generalizes well. In the case of a parametric model, we learn parameters w ∈ R d . 
 
 Appendix:  


Decision: The paper presents a study that explores the use of GPs for Q-function approximation, which trades-off interpretability vs. performance. While the reviewers raised concerns about the novelty and significance of the outcomes, the analysis is sound, and the insights are useful for the RL community. Hence, with the expectation that the camera-ready version will place a stronger focus on interpretability to better substantiate the claims, I recommend accepting the paper.

Reviews: review: title: title: Contribution is weak and writing should be formal
paper_summary: paper_summary: Gaussian processes are a non-standard function approximator for reinforcement learning value functions. RL value functions are usually approximated by deep neural networks or kernel regression. The GP is a flexible non-parametric estimator and has been studied in the past in the context of RL function approximation. This paper uses SVGP and Deep GPs for off-policy and on-policy learning. These provide the advantage of uncertainty estimation and interpretability. Numerical comparisons are performed in the CartPole and Lunar Lander environments.
main_review: main_review: Based on the discussion in the committee chat and the discussion with the Area chairs, I am increasing my scores. The paper does not make a significant contribution and the results are weak. However, there is nothing wrong with the method or the evaluation. Hence, I am increasing the scores.
strength_weakness: strength_weakness: There is a definite improvement needed in the contributions and writing before being published. While the experimental results are unfavorable, a clear illustration of the advantages offered by the paper along with more formal explanations of the math and algorithms would have pushed the paper towards a higher score. The idea of using Gaussian processes for function approximation is not new. More strong experimental evidence on other benchmarks could be useful to ascertain further.