Bayesian policy gradient and actor-critic algorithms

In plain words

A learning agent improves by nudging its policy, the rule that picks its actions, in the direction that earns more reward. The usual way to find that direction is to run many trial episodes and average what happened. The average is noisy, so it takes a lot of episodes before the direction is reliable.

This paper treats the unknown direction as a smooth random function. It starts from a prior guess about how that function behaves, looks at the episodes it has seen, and computes a best estimate together with a measure of how unsure it still is. The prior is chosen from the policy's own scoring function, so the math works out in closed form and the estimate costs little extra. The same idea is applied to a critic, the part that judges how good each action is, so the method can learn from single steps when the problem has the Markov property. The authors report that their estimates need fewer episodes than plain averaging on small test problems, and that the critic version learns faster on several control tasks.

TL;DR

Bayesian policy gradient (BPG) places a Gaussian process prior on the integrand of the policy gradient and applies Bayesian quadrature, which gives a posterior mean, a posterior covariance and a natural-gradient variant in closed form when the kernel is a Fisher kernel §4 tab. 1. Bayesian actor-critic (BAC) uses a GPTD critic with a state-action Fisher kernel and the same closed-form gradient moments, with single transitions as the observation unit §7. Author-reported: BPG's estimates need fewer samples than Monte-Carlo on a bandit and a linear quadratic regulator (LQR) §6.1 tab. 2 fig. 1; BAC reaches better policies than Monte-Carlo policy gradient on a random walk, mountain car and ship steering §8 fig. 7 fig. 8 fig. 9.

Problem

Policy gradient methods estimate the gradient of the expected return from sampled trajectories with Monte-Carlo averages. The estimator is unbiased but has high variance, so accurate estimates need many samples and convergence is slow §1. The authors also argue that Monte-Carlo estimation violates the likelihood principle, and that it does not make full use of the data at hand §1 §3.

Assumptions: the policy must be continuously differentiable in its parameters (Assumption 1) §2. The return may be discounted, with gamma in [0, 1] §2. Episodes are finite in the LQR and random-walk experiments, and the authors state that BQ-based trajectory models degrade when trajectory lengths vary §8.

Cost of the Bayesian method: the closed forms need a specific kernel. The authors note that this restriction can introduce approximation error that is hard to quantify, and that they did not see a significant error in their experiments §4 §7.2. BPG treats whole trajectories as the observation unit, so it cannot use the Markov property §1 §4. BAC restores that property but needs a state kernel and a Fisher information matrix estimate §7.2.

Method

Background: Bayesian quadrature §3

To estimate an integral rho = integral of f(x) g(x) dx, model f as a Gaussian process with prior mean f_bar and kernel k. Condition on noisy samples y_i of f at points x_i. The posterior of rho is Gaussian, with mean rho_0 + b^T C (y − f_bar) and variance b_0 − b^T C b, where C = (K + Sigma)^{-1}, b = integral of k(x) g(x) dx and b_0 = integral of integral of k(x, x') g g' §3. Samples may come from any distribution, not only from g, which is the stated advantage over Monte-Carlo §3. For the formulas to be usable, the three integrals must be analytic; the authors call this the requirement for the kernel and the factorisation to be chosen so that they can be solved §3. Section 3.1 extends the formulas to vector-valued integrals and states the simplifying assumption that each component is independent with a shared kernel §3.1 [prop. 1] [prop. 2].

Gradient as a Bayesian integral §4

The paper uses the score-function form of the gradient §2 eq. 4 eq. 5: the gradient of the expected return is the integral over paths xi of R(xi) times the score u(xi) = grad log Pr(xi; theta) times Pr(xi; theta). The score is a sum over time steps of grad log mu(a_t | x_t; theta) eq. 5. The Fisher information matrix is G(theta) = E[u u^T] §4.1. The BPG posterior is taken over the gradient of the expected return, with the Bayesian expected return treated as random §4.

Two ways of splitting the integrand give two models §4 tab. 1.

Model 1 (vector-valued GP, §4.1). The deterministic factor is g = Pr(xi; theta) and the GP factor is f = R_bar(xi) grad log Pr(xi; theta). Measurements are y(xi) = R(xi) u(xi). The kernel is the quadratic Fisher kernel k(xi, xi') = (1 + u(xi)^T G^{-1} u(xi'))^2 §4.1. Proposition 3 gives the closed forms b_i = 1 + u_i^T G^{-1} u_i and b_0 = 1 + n, where n is the number of policy parameters §4.1 [prop. 3]. The posterior is E[grad eta | D] = Y C b and Cov = (b_0 − b^T C b) I, where Y stacks the measurement vectors as rows §4.1 tab. 1.

Model 2 (scalar-valued GP, §4.2). The deterministic factor is g = grad Pr(xi; theta) and the GP factor is f = R_bar(xi), with measurement y(xi) = R(xi). The kernel is the Fisher kernel k(xi, xi') = u^T G^{-1} u' §4.2. Proposition 4 gives B = U, the matrix whose columns are the score vectors, and B_0 = G §4.2 [prop. 4]. The posterior is E[grad eta | D] = B C y and Cov = G − B C B^T §4.2 tab. 1.

Model 1 uses one kernel and one noise covariance for all components, so the covariance comes out as a scalar times the identity §4.1 §3.1. Model 2 is simpler: one prior mean, one kernel and one noise matrix §4.2.

Algorithm 1: BPG evaluation alg. 1

Input: policy parameters theta, sample size M, and the model. Output: posterior mean and covariance of the gradient alg. 1. The algorithm samples M paths and updates the kernel matrix K path by path, using the new kernel column and diagonal entry §4.3. The noise matrix Sigma is added before the inverse C = (K + Sigma)^{-1} is formed alg. 1, line 12.

# BPG_Eval  [alg. 1]
# Input:  theta (n,), sample size M, model in {1, 2}, Fisher G (n, n)
# Output: posterior mean (n,) and covariance (n, n) of the gradient

def bpg_eval(theta, M, model, fisher):
    # Sample M paths under the current policy  [alg. 1, lines 3-4]
    paths = [sample_path(theta) for _ in range(M)]
    # Score vector u(xi_i), a sum over steps of grad log mu  [§2] [eq. 5]
    scores = [score(path, theta) for path in paths]      # (M, n)
    # Return R(xi_i), the sum of the rewards along the path  [alg. 1]
    returns = [total_reward(path) for path in paths]     # (M,)
    fisher_inv = inverse(fisher)                         # (n, n)
    # Gram matrix of the kernel of this model  [§4.1] [§4.2]
    gram = kernel_matrix(scores, fisher_inv, model)      # (M, M)
    # Noise covariance Sigma; the text gives it only for noisy rewards  [§6.2.1]
    C = inverse(gram + noise_matrix(scores, returns))    # (M, M)
    if model == 1:
        # Measurement R(xi_i) u(xi_i), one row per path  [§4.1]
        Y = returns[:, None] * scores                    # (M, n)
        # b_i = 1 + u_i^T G^{-1} u_i  [prop. 3]
        b = 1 + sum(scores * (scores @ fisher_inv), axis=1)  # (M,)
        mean = Y.T @ C @ b                               # (n,)
        # Scalar b_0 = n + 1 times the identity  [prop. 3] [tab. 1]
        cov = (n + 1 - b @ C @ b) * identity(n)          # (n, n)
    else:
        # Score matrix B = U and B_0 = G  [prop. 4]
        B = scores.T                                     # (n, M)
        mean = B @ C @ returns                           # (n,)
        cov = fisher - B @ C @ B.T                       # (n, n)
    return mean, cov

The Fisher matrix G is assumed known in this algorithm. In practice it must be estimated §4.3. The paper gives two options: a Monte-Carlo estimate from the sampled paths, G_hat_M = (1/M) sum u u^T (unbiased) with a recursive or Sherman-Morrison update of its inverse, or a model-based estimate when the transition probabilities P are modelled and fitted by maximum likelihood §4.3.

Sparsification §4.4

A new path is added to a dictionary only if k(xi, xi) − k_t^T K_t^{-1} k_t exceeds a threshold tau §4.4. With the sparse dictionary matrices, Proposition 5 gives the posterior moments as below, where A is a matrix that maps each sample to its dictionary element (or to a projection of it) [prop. 5] §4.4. Sigma^{-1} appears in every form, so the sparse formulas require a non-singular Sigma [prop. 5].

E[\nabla\eta_B | D] = Y \Sigma^{-1} A (\tilde{K} A^T \Sigma^{-1} A + I)^{-1} \tilde{b} [prop. 5, Model 1]

Cov = [(n+1) − \tilde{b}^T A^T \Sigma^{-1} A (\tilde{K} A^T \Sigma^{-1} A + I)^{-1} \tilde{b}] I [prop. 5, Model 1]

E[\nabla\eta_B | D] = \tilde{B} (A^T \Sigma^{-1} A \tilde{K} + I)^{-1} A^T \Sigma^{-1} y [prop. 5, Model 2]

Cov = G − \tilde{B} (A^T \Sigma^{-1} A \tilde{K} + I)^{-1} A^T \Sigma^{-1} A \tilde{B}^T [prop. 5, Model 2]

The paper says sparsification makes the algorithms faster and more memory efficient, and may reduce accuracy §4.4. It also says that sparsification helps numerically when the kernel matrix is singular or nearly singular §4.4.

Algorithm 2: BPG policy update alg. 2

Repeat N times: evaluate the posterior mean of the gradient with Algorithm 1 at the current parameters, then step theta_{j+1} = theta_j + beta_j * Delta. The conventional gradient uses Delta = posterior mean. The natural-gradient variant multiplies the mean by G(theta_j)^{-1} alg. 2. The posterior covariance is not used in the update in this algorithm §6.2.2.

# BPG  [alg. 2]
# Input:  theta_0 (n,), learning rates beta_j, number of updates N, sample M
def bpg(theta_0, betas, N, M, natural):
    theta = theta_0
    for j in range(N):
        # Posterior mean of the gradient from Algorithm 1  [alg. 1]
        mean, cov = bpg_eval(theta, M, model=1, fisher=fisher_of(theta))
        if natural:
            # Natural gradient: precondition by the inverse Fisher  [alg. 2]
            step = inverse(fisher_of(theta)) @ mean    # (n,)
        else:
            step = mean                                 # (n,)
        theta = theta + betas[j] * step
    return theta

Partially observable case §5

The same models and algorithms apply to POMDPs with observations in place of states, because the score of a path has the same form §5. The paper notes that the gradient can also be used with conjugate-gradient and line-search methods, as in GPOMDP §5. These options are not used in the experiments.

Bayesian actor-critic §7

BAC uses the state-action pair as the observation unit. The action-value function Q gets a GP prior with kernel k(z, z') = k_x(x, x') + k_F(z, z'), where k_F is the Fisher kernel u(z)^T G^{-1} u(z') on state-action pairs and k_x is any state kernel §7.2. The critic is GPTD, which models each reward as r(z_i) = Q(z_i) − gamma Q(z_{i+1}) + N §7.1. Its posterior mean and covariance are Q_hat = k_t(z)^T alpha and S_hat = k(z, z') − k_t(z)^T C k_t(z') §7.1.

With Proposition 6, the gradient posterior is E = U_t alpha_t and Cov = G − U_t C_t U_t^T, where U_t stacks the score vectors of the visited pairs and G is the Fisher matrix §7.2 [prop. 6]. The Fisher matrix is estimated from data: for a stationary state distribution, the paper uses G_hat_t = U_t U_t^T / (t + 1) §7.2. For the discounted state distribution, which the paper uses, it gives Algorithm 3 (G-EST) §7.2 alg. 3.

# GPTD posterior for one run of observed transitions  [§7.1]
# K: Gram matrix of the visited state-action pairs, (t+1, t+1)
# rewards: observed rewards, (t,)
# sigma: noise scale; the value used in the experiments is not stated
def gptd(K, rewards, gamma, sigma):
    # H is (t, t+1): ones on the diagonal, -gamma just above it  [§7.1]
    H = bidiagonal_matrix(t, gamma)                      # (t, t+1)
    S = inverse(H @ K @ H.T + sigma**2 * H @ H.T)        # (t, t)
    alpha = H.T @ S @ rewards                            # (t+1,)
    C = H.T @ S @ H                                      # (t+1, t+1)
    return alpha, C

# BAC  [alg. 4]
# Input:  theta (n,), M episodes per update, epsilon, natural flag
def bac(theta, M, epsilon, natural):
    while True:
        # Run GPTD over M episodes  [alg. 4] [§7.1]
        alpha, C = run_gptd(theta, M)
        # Scores of the visited pairs, U = [u(z_0), ..., u(z_t)]  [prop. 6]
        U = score_matrix(theta)                          # (n, t+1)
        # Fisher estimate from Eq. Ghat or from G-EST  [alg. 3] [§7.2]
        fisher = estimate_fisher(theta)                  # (n, n)
        if natural:
            delta = inverse(fisher) @ U @ alpha          # (n,)
        else:
            delta = U @ alpha                            # (n,)
        theta = theta + beta * delta
        # Stop when the update is small  [alg. 4]
        if norm(delta) < epsilon:
            return theta

The sparse BAC uses a dictionary of state-action pairs and the sparse posterior of Proposition 7, with the threshold tau on the same squared distance as in §4.4 §7.3 [prop. 7]. The ship-steering experiments use online sparsification §8.3. The pseudocode above shows the dense form.

Hyperparameters

Values are as stated in the paper. "not stated" means the text does not give it.

Learning rates, as printed:

tab. 3 — initial learning rates for the LQR policy-gradient algorithms, one per gradient component.

tab. 4 — learning rates for random-walk policy learning.

tab. 5 — learning-rate parameters for mountain car; beta_c = inf means a fixed rate.

tab. 6 — learning rates for ship steering.

Results

Author-reported. Setup: simulated tasks with the parameters in Hyperparameters; each number is the paper's own. Apart from tab. 2 and the LQR optimum, the numbers appear only in figures fig. 1–fig. 9, which this card does not transcribe.

Bandit gradient estimates §6.1 tab. 2: exact gradient, Monte-Carlo (MC) and Bayesian quadrature (BQ) estimates, mean ± standard deviation over 10^4 runs.

The text says BQ's variance is about an order of magnitude below MC's at M = 10 and about six orders below at M = 100 §6.1.

LQR gradient estimates §6.2.1 fig. 1 fig. 2: mean squared error and mean absolute angular error versus M, for Model 1 and Model 2 with sparsification. The text says BQ beats MC at these sample sizes. With noisy rewards (sigma_r^2 = 1), the error increases slightly and the picture stays similar §6.2.1 fig. 2.

LQR policy optimisation §6.2.2 fig. 3 fig. 4 fig. 5: BPG and BPNG (natural) beat MCPG for M = 10, 20 and 40 and at M = 5 only MCPG is better §6.2.2 fig. 3. Using ML instead of MC for the Fisher matrix still beats MCPG fig. 4. BPG-var, which scales the step by the posterior covariance, beats BPG and MCPG at all four sample sizes fig. 5. Reported optimum: lambda* ≈ −0.92, sigma* = 0.001, eta_B = 0.3067, in the cost convention §6.2.2.

Random walk §8.1 fig. 6 fig. 7: BAC's gradient estimates are more accurate and have lower variance than BQ's and MC's §8.1 fig. 6. In policy learning, the text says BAC learns clearly faster than the other algorithms §8.1 fig. 7. The text also says BPG improves with M and is similar to BAC at M = 100, but the experiments use only M = 1, 25, 50 and 75 §8.1 (see Pitfalls).

Mountain car §8.2 fig. 8: BAC is better and has lower variance than MCPG. BAC reaches a good policy with M = 5 and changes little as M grows. MCPG catches up only at M = 40 §8.2.

Ship steering §8.3 fig. 9: BAC converges to a higher success ratio than MCPG for all M. For M = 5, MCPG is slightly better early on §8.3.

Paper vs code

No official code found. The paper's footnote points to a project page, https://sequel.lille.inria.fr/Software/BAC §8, which could not be reached (two connection attempts were refused). GitHub search for "bayesian actor-critic" with the author name, and for "bayesian policy gradient", found no repository under the authors or their lab. The third-party repositories found were not used.

Pitfalls

• Sign convention for the LQR. The paper's rewards are costs and the LQR is minimised, but the update rules in alg. 2 and alg. 4 are written as ascent (theta + beta * delta). The authors explain the cost convention in a footnote §6.2 but do not show the sign change in the update. An implementer must choose one convention and apply it throughout.
• Noise matrix Sigma is not stated for noise-free rewards. The noise covariances are given only for the noisy-reward LQR experiment §6.2.1. For Model 2 the kernel is u^T G^{-1} u', so the Gram matrix has rank at most n (the number of policy parameters). For M > n it is singular unless Sigma or a nugget is added. This is an inference from the kernel form §4.2. The sparse formulas need Sigma^{-1} [prop. 5], so the noise-free sparse runs need a non-zero Sigma or nugget that the paper does not give.
• Fisher matrix estimates break the closed forms. Proposition 3 holds exactly only if G is the true expectation E[u u^T] app. A. With an estimated G, the b and b_0 identities are approximations, and the covariance formula is no longer exact. The paper estimates G in its experiments §4.3 §7.2.
• Model 1 covariance is isotropic. One kernel and one noise matrix are shared across components, so the posterior covariance is a scalar times I §4.1. Cross-component correlations are ignored, and the authors say this may use the data sub-optimally §3.1.
• Two notions of the sample. BPG and BAC use different observation units: whole trajectories and single transitions §1 §7. Switching between them changes the score and Fisher definitions §7.2. Do not mix them.
• BPG is not used when trajectory lengths vary. The authors say BQ's performance degrades when trajectories differ in length, so they used BAC only for mountain car and ship steering §8 remark.
• Gaussian arguments are ambiguous. The LQR gives x_0 ~ N(0.3, 0.001) and the noise as N(0, 0.01), while the reward-noise covariance is written sigma_r^2 §6.2. The paper does not say whether the second argument is a variance or a standard deviation. The reward-noise notation suggests variance, but this is not stated.
• BPNG step size. The BPNG learning rate is a fixed rate multiplied by the determinant of the Fisher matrix §6.2.2. This is unusual, and the paper gives no reason for it.
• Mountain car grid. The text gives the Gaussian centres on the grid {0, 0.25, 0.5, 1} × {0, 0.25, 0.5, 1}, with spacing 0.5 at the top end. This may be a typo for 0.75, but the paper does not say §8.2. The numeric values of the three actions (forward, reverse, zero) are also not stated §8.2.
• Ship steering naming. The text calls the map from the Gaussian sample to the action a sigmoid, but the formula is an arctan §8.3. The symbol T is used both for the lag (5 steps) and for the number of updates (3000 or 1000) §8.3.
• Ship steering state scale. The state is in metres and degrees with sigma_k = 1. The paper does not normalise the state, and sigma_k = 1 is set without comment §8.3.
• G-EST is hard to read. Algorithm 3 uses termination flags and a 1/(1 − gamma) weight on the final segment alg. 3. The paper gives the restart interpretation in words §7.2, but not a derivation of the normalisation. Check the estimator against the stated restart distribution before relying on it.
• Stationary-state estimate in an episodic task. Eq. Ghat (1/(t+1) U U^T) assumes a stationary state distribution §7.2. The random-walk and mountain-car tasks are episodic with absorbing states, so this is an approximation there. The paper does not say how large the error is.
• Figures are not transcribed. Most results appear only in figures §6.2 §8. Numerical comparisons beyond tab. 2 require reading the figures.
• BPG results do not match the stated sample sizes. The text says BPG is similar to BAC at M = 100 §8.1, but the random-walk policy-learning runs use only M = 1, 25, 50 and 75 §8.1 tab. 4. Treat the M = 100 claim as unverified.
• Source template. The TeX preamble uses the JMLR 2016 template (volume 17). The conference versions are the 2006 (BPG) and 2007 (BAC) papers §1. This paper extends them with the BQ treatment of vector integrals, the POMDP extension, the proofs and two new domains §1.

Lineage

• Builds on: Bayesian quadrature (O'Hagan 1991) §3; Fisher kernels §4.1 §4.2; GPTD, from the Engel works the paper cites §7.1; the GPOMDP setting §5. The 2006 (BPG) and 2007 (BAC) conference papers are the earlier forms of this work §1.
• Compared against: Monte-Carlo policy gradient (MCPG, Algorithm 1 of the GPOMDP paper, cited as Baxter01IP) §8.
• Related: a 2011 paper (Vien11HM) that uses BQ to set the learning-rate schedule from the Hessian §6.2.2.
• Superseded by: not known from the content of this paper.