# Modified discrete Laguerre polynomials for efficient computation of
exponentially bounded Matsubara sums

Guanpeng A. Xu and Steven G. Johnson Department of Mathematics, Massachusetts
Institute of Technology, Cambridge, MA<EMAIL_ADDRESS>
###### Abstract

We develop a new type of orthogonal polynomial, the modified discrete Laguerre
(MDL) polynomials, designed to accelerate the computation of bosonic Matsubara
sums in statistical physics. The MDL polynomials lead to a rapidly convergent
Gaussian “quadrature” scheme for Matsubara sums, and more generally for any
sum $F(0)/2+F(h)+F(2h)+\cdots$ of exponentially decaying summands
$F(nh)=f(nh)e^{-nhs}$ where $hs>0$. We demonstrate this technique for
computation of finite-temperature Casimir forces arising from quantum field
theory, where evaluation of the summand $F$ requires expensive electromagnetic
simulations. A key advantage of our scheme, compared to previous methods, is
that the convergence rate is nearly independent of the spacing $h$
(proportional to the thermodynamic temperature). We also prove convergence for
any polynomially decaying $F$.

###### keywords:

Gaussian quadrature, orthogonal polynomials, Matsubara summation

33C45, 65D32, 65B10, 81T55

##  1 Introduction

Our work develops an efficient $N$-point scheme, analogous to Gauss–Laguerre
quadrature, to rapidly evaluate exponentially decaying sums of the form

(1)
$\sum\limits_{n=0}^{\infty}{}^{\prime}\underbrace{f(nh)e^{-nhs}}_{F(nh)}h\approx\sum_{k=1}^{N}w_{k}F(x_{k})\,,$

where $\sum^{\prime}$ denotes that the $n=0$ term is weighted by $1/2$,
exploiting knowledge of an asymptotic exponential decay rate $s>0$ of the
summand $F(x)$ to derive quadrature points $x_{k}$ and weights $w_{k}$ for any
desired order $N$ of accuracy. This calculation arises for bosonic Matsubara
sums in statistical physics [14]—for example, in the computation of finite-
temperature Casimir forces from quantum field theory [14, 21, 8]—in which case
$h$ is proportional to the temperature $T$, evaluation of $F$ can involve
expensive partial-differential equation (PDE) solutions, and the asymptotic
decay rate $s$ can be determined from a Born approximation [15]. Even though
Eq. 1 is exponentially decaying, direct evaluation of this sum requires many
summand evaluations in the common case where $hs$ is small, corresponding to
low temperatures. In the $h\to 0^{+}$ limit, Eq. 1 becomes
$\int_{0}^{\infty}f(x)e^{-sx}dx=\int_{0}^{\infty}f(y/s)e^{-y}dy/s$ and can be
evaluated efficiently by Gauss–Laguerre quadrature [13, 8], so our goal was to
develop a “Gaussian summation” scheme (equivalent to Gaussian quadrature with
a discrete measure [3]) that extends Gauss–Laguerre quadrature to $h>0$ for
this type of sum (Section 2). We demonstrate our “modified discrete Laguerre”
(MDL) approach on a typical Casimir-force computation (Section 4.1) shown in
Fig. 1, where evaluating $F$ requires the solution of Maxwell’s equations [8,
19, 20]. We found that MDL summation requires many fewer points $N$ than naive
summation at low temperatures, and also outperforms a competing Gaussian-
quadrature scheme from the literature on Matsubara sums [16, 11] that did not
exploit knowledge of the asymptotic exponential decay rate. Ours was the only
technique whose performance did not degrade as temperature ($\sim h$)
decreased. We anticipate that similar benefits from MDL summation should apply
to other physical problems involving Matsubara sums, e.g. for density
functional theory where various Matsubara-summation methods have been proposed
[9, 11, 10], and more generally for evaluating any slowly converging
exponentially weighted sums.

Figure 1: We aim to reduce the number of summand evaluations (left) through a
quadrature rule (right). The 4-point quadrature rule attains identical
accuracy to a 20-term partial Matsubara sum at $T=100$ K. We exploit the fact
that the asymptotic decay rate (black dashed line) is known to be $e^{-2xd/c}$
for objects with minimum separation $d$ ($=1\,\mu m$ here).

The key to an efficient summation scheme for Eq. 1, as for all Gaussian-
quadrature methods [13], is the construction of a set of orthogonal
polynomials—in this case for the discrete inner product

(2) $\langle
f,g\rangle_{h,s}^{\prime}=\sum\limits_{n=0}^{\infty}{}^{\prime}f(nh)g(nh)e^{-nhs}h\,.$

We call our resulting polynomials the “modified discrete Laguerre” (MDL)
polynomials: they reduce to the ordinary Laguerre polynomials [13] in the
limit $h\to 0^{+}$ for $s=1$, and they are closely related to the discrete
Laguerre (DL) polynomials derived in Ref. [6], which differ only in the $n=0$
term of the inner product (using $\sum$ instead of $\sum^{\prime}$).
Regardless of whether the inner product is an integral or a discrete sum,
orthogonal polynomials (described by their three-term recurrence, derived in
Section 2) immediately yield the weights $w_{n}$ and points $x_{n}$ in Eq. 1
(e.g. by the Golub–Welsch algorithm [5, 28, 25]) to exactly sum/integrate
polynomials up to degree $2N-1$ and exhibit rapid convergence for a broad
class of smooth functions $f$ [2, 26]. In principle, one could use the
original DL polynomials to evaluate Eq. 1 by separately subtracting $F(0)/2$,
but this would require one additional function evaluation beyond the $x_{k}>0$
quadrature points. Worse, for the Casimir case, $F(0)$ involves a singularity
that must be evaluated by an expensive $x\to 0^{+}$ numerical-extrapolation
procedure [8], whereas MDL summation only evaluates $F$ only at $x>0$ and
implicitly computes the contribution of the limit. In contrast to direct
summation and previous Gaussian quadrature techniques, our MDL scheme
maintains good convergence even at low temperature (small $h$) where the
direct summand decays more and more slowly. Finally, in Appendix B, we adapt
previous work on the convergence of Gaussian quadrature [29] for improper
integrals in order to prove convergence of MDL quadrature for any $f(x)$
decaying faster than $x^{-1-\epsilon}$ for $\epsilon>0$. For fermionic
Matsubara summation ($\sum$ instead of $\sum^{\prime}$ of $F(nh+h/2)$), one
would simply apply an analogous quadrature scheme based on the original DL
polynomials instead of the MDL polynomials.

##  2 MDL polynomials and their three-term recurrence

In this section, we derive certain key properties of the MDL polynomials,
culminating in the three-term recurrence of Theorem 2.6. Given the three-term
recurrence coefficients, standard techniques such as the Golub–Welsch
algorithm [5, 1, 25, 28] immediately yield the quadrature points and weights
of Section 3.1, shown in Fig. 1 and employed in Section 4.1. We begin with
some definitions and preliminary results adapted from Ref. [6], and outline
the steps leading to Theorem 2.6. The details of the inductive proofs are
given in Appendix A.

For physical applications or for comparison to Gauss–Laguerre quadrature, it
is convenient to separate the decay rate $s$ from the spacing $h$ in Eq. 1,
but in deriving the MDL polynomials it is more natural to combine these into a
single parameter $\tau=e^{hs}>1$ and omit the overall $h$ scale factor; we
will recover the original sum by a change of variables in Section 3. In terms
of $\tau$, we define the inner products $\langle\cdot,\cdot\rangle$ and
$\langle\cdot,\cdot\rangle^{\prime}$ of two real-valued functions $f$ and $g$
as follows, analogous to Eq. 1:

(3) $\langle
f,g\rangle(^{\prime})=\sum\limits_{n=0}^{\infty}(^{\prime})\tau^{-n}f(n)g(n)\,,$

where the notation $\sum\limits_{i=0}^{L}{}^{\prime}s_{i}$ represents the
summation $\frac{s_{0}}{2}+\sum\limits_{i=1}^{L}s_{i}$. The discrete Laguerre
(DL) polynomials of Ref. [6] are orthogonal with respect to the
$\langle\cdot,\cdot\rangle$. Here, we introduce “modified” discrete Laguerre
(MDL) polynomials that are orthogonal with respect to
$\langle\cdot,\cdot\rangle^{\prime}$. We denote the DL and MDL polynomials as
$L_{n}(x)$ and $L_{n}^{\prime}(x)$, respectively.

Below, we will make use of the following results proved in Ref. [6]. The DL
polynomials for $n\geq 0$ may be written in the closed form

(4)
$L_{k}(n)=\sum\limits_{i=0}^{k}(-1)^{i}\binom{k}{i}\frac{(1-\frac{1}{\tau})^{i}}{i!}\prod\limits_{j=1}^{i}(n+j).$

Furthermore, the DL polynomials satisfy the following normalization condition:

(5) $L_{n}(0)=\frac{1}{\tau^{n}},\;\;\;\;\langle
L_{n},L_{n}\rangle=\frac{1}{\tau^{n-1}(\tau-1)}.$

We now derive the coefficients of the three-term recurrence for the MDL
polynomials $L_{n}^{\prime}$. We proceed by analyzing the difference between
$L_{n}^{\prime}$ and $L_{n}$ using induction on $n$, arriving at the following
results:

######  Lemma 2.1.

For all nonnegative integers $n$, we have:

$L_{n}^{\prime}=L_{n}+\frac{1}{2}\frac{1}{\tau^{n}}\sum\limits_{i=0}^{n-1}\frac{L_{i}^{\prime}(0)L_{i}^{\prime}}{\langle
L_{i}^{\prime},L_{i}^{\prime}\rangle^{\prime}}.$

######  Proof 2.2.

We apply the Gram–Schmidt procedure [28] to re-orthogonalize $L_{n}$ with
respect to the modified inner product $\langle\cdot,\cdot\rangle^{\prime}$,
using Eq. 5:

$\displaystyle L_{n}^{\prime}$
$\displaystyle=L_{n}-\sum\limits_{i=0}^{n-1}\frac{\langle
L_{n},L_{i}^{\prime}\rangle^{\prime}}{\langle
L_{i}^{\prime},L_{i}^{\prime}\rangle^{\prime}}L_{i}^{\prime}=L_{n}+\frac{1}{2}\sum\limits_{i=0}^{n-1}\frac{L_{n}(0)L_{i}^{\prime}(0)}{\langle
L_{i}^{\prime},L_{i}^{\prime}\rangle^{\prime}}L_{i}^{\prime}$
$\displaystyle=L_{n}+\frac{1}{2}\sum\limits_{i=0}^{n-1}\frac{L_{i}^{\prime}(0)L_{i}^{\prime}}{\langle
L_{i}^{\prime},L_{i}^{\prime}\rangle^{\prime}\tau^{n}}=L_{n}+\frac{1}{2}\frac{1}{\tau^{n}}\sum\limits_{i=0}^{n-1}\frac{L_{i}^{\prime}(0)L_{i}^{\prime}}{\langle
L_{i}^{\prime},L_{i}^{\prime}\rangle^{\prime}}.$

The following statement is proved by induction on $n$ in Appendix A:

######  Lemma 2.3.

For all nonnegative integers $n$, we have:

$L_{n}^{\prime}(0)=\frac{2}{1+\tau^{n}},\;\;\;\;\langle
L_{n}^{\prime},L_{n}^{\prime}\rangle^{\prime}=\frac{1+\tau^{n+1}}{\tau^{n}(1+\tau^{n})(\tau-1)}.$

We now derive the three-term recurrence for the $L_{n}^{\prime}$ polynomials.

######  Theorem 2.4.

For all nonnegative $n$, the MDL polynomials $L_{n}^{\prime}$ satisfy the
recurrence

$L_{n+1}^{\prime}(x)=(-\alpha_{n}x+\beta_{n})L_{n}^{\prime}(x)-\gamma_{n}L_{n-1}^{\prime}(x),$

where

$\displaystyle\alpha_{n}=\frac{\tau-1}{(n+1)\tau}$
$\displaystyle\beta_{n}=\frac{1+\tau^{n}}{1+\tau^{n+1}}+\frac{n}{\tau(n+1)}\frac{1+\tau^{n+1}}{1+\tau^{n}}$
$\displaystyle\gamma_{n}=\frac{n}{n+1}\frac{(1+\tau^{n+1})(1+\tau^{n-1})}{\tau(1+\tau^{n})^{2}}.$

######  Proof 2.5.

See Appendix A.

In order to apply the MDL polynomials to Gaussian-quadrature applications, it
is necessary to define rescaled orthonormal variants
$\hat{L}_{n}^{\prime}=L_{n}^{\prime}/\sqrt{\langle
L_{n}^{\prime},L_{n}^{\prime}\rangle}$. Using Lemma 2.3, we derive the
rescaling

$\hat{L}_{n}^{\prime}=L_{n}^{\prime}\cdot(-1)^{n}\sqrt{\frac{\tau^{n}(1+\tau^{n})(\tau-1)}{1+\tau^{n+1}}}.$

By straightforward algebraic substitution of this rescaling into Theorem 2.4,
we obtain a recurrence for the normalized MDL polynomials (whose coefficients
appear directly in the Jacobi matrix used for Gaussian quadrature in Section
3):

######  Theorem 2.6.

Define two sequences $\\{\hat{\alpha}_{n}\\}_{n},\\{B_{n}\\}_{n}$ as follows:

$\displaystyle\hat{\alpha}_{n}$
$\displaystyle=\frac{(n+1)\tau}{\tau-1}\left(\frac{1+\tau^{n}}{1+\tau^{n+1}}+\frac{n}{\tau(n+1)}\frac{1+\tau^{n+1}}{1+\tau^{n}}\right)$
$\displaystyle\hat{\beta}_{n}\,,$
$\displaystyle=\frac{(n+1)\tau}{\tau-1}\sqrt{\frac{(1+\tau^{n})(1+\tau^{n+2})}{\tau(1+\tau^{n+1})^{2}}}\,.$

Then the normalized polynomials $\hat{L}_{n}^{\prime}$ satisfy the recurrence

$x\hat{L}_{n}^{\prime}(x)=\hat{\beta}_{n}\hat{L}_{n+1}^{\prime}(x)+\hat{\alpha}_{n}\hat{L}_{n}^{\prime}(x)+\hat{\beta}_{n-1}\widehat{L_{n-1}^{\prime}}(x).$

We may see with direct calculation of coefficients that $\hat{\alpha}_{n}$ and
$\hat{\beta}_{n}$ are positive, monotonically increasing, and asymptotically
proportional to $n$ for any fixed $\tau>0$.

##  3 MDL quadrature for Matsubara sums

In order to evaluate the Matsubara sum Eq. 1, we now consider a broader class
of sums, characterized by the generalized (M)DL inner product

(6) $\langle
f,g\rangle_{h,s}(^{\prime})=\sum\limits_{n=0}^{\infty}(^{\prime})\,f(nh)g(nh)e^{-nhs}h.$

Here, we may set $\tau=e^{hs}$; in this case, the generalized (M)DL
polynomials correspond to a simple rescaling with respect to the original
(M)DL polynomials. In particular, we define generalized MDL polynomials for
this inner product via a change of variables:

$\hat{L}_{k}^{(h,s)\prime}(nh)=\hat{L}_{k}^{\prime}(n).$

By substituting this change of variables into Theorem 2.6, we immediately see
that $\hat{L}_{k}^{(h,s)\prime}(nh)$ satisfies an analogous three-term
recurrence with coefficients $h\hat{\alpha}_{n}$ and $h\hat{\beta}_{n}$.

We compute the quadrature points (roots of $\hat{L}_{k}^{(h,s)\prime}$) and
weights from the eigenvalues and eigenvectors of the symmetric tridiagonal
Jacobi matrix (the Golub–Welsch algorithm):

(7) $J_{N}=h\begin{pmatrix}\hat{\alpha}_{0}&\hat{\beta}_{0}&&&&\\\
\hat{\beta}_{0}&\hat{\alpha}_{1}&\hat{\beta}_{1}&&&\\\
&\hat{\beta}_{1}&\ddots&\ddots&&\\\ &&\ddots&\ddots&\hat{\beta}_{N-3}&\\\
&&&\hat{\beta}_{N-3}&\hat{\alpha}_{N-2}&\hat{\beta}_{N-2}\\\
&&&&\hat{\beta}_{N-2}&\hat{\alpha}_{N-1}\end{pmatrix}\,.$

In the following subsections, we examine the properties of these roots and the
resulting quadrature scheme for computing the sum Eq. 7.

###  3.1 Roots and asymptotic properties of MDL polynomials

The roots of $\hat{L}_{N}^{(h,s)\prime}$ correspond to the eigenvalues of the
$N\times N$ Jacobi matrix Eq. 7 [5, 28]. These roots are plotted as a function
of $hs=\ln(\tau)$ in Fig. 3, and exhibit asymptotic behaviors that we briefly
explain here.

Figure 2: The roots of the modified discrete Laguerre (MDL) polynomials for
$n=10$, $s=1$, and $\ln(\tau)=hs\in[0,4]$.

Figure 3: Differences $x_{k}-kh$ between the roots of the MDL polynomials
$x_{k}$ and the Matsubara points $kh$ for the same $n,h,s$.

Due to Proposition 2.3 of Ref. [3], the roots $x_{k}$ satisfy
$kh<x_{k}<(k+1)h$ (for $k=0,1,\ldots$), and in the limit of large $h$ (high
temperatures) the roots $x_{k}$ approach $kh$ from above. Correspondingly, the
quadrature scheme for Eq. 1 reduces to naive summation in this limit.
(Intuitively, in the large-$h$ regime the summand weight $e^{-shn}$ decays so
rapidly that naive summation is nearly optimal.) Quantitatively, the off-
diagonal entries $\hat{\beta}_{n}$ of $J_{N}$ decay exponentially in
$h=\ln(\tau)$, whereas the diagonal entries approach multiples of $h$. In
other words, the Jacobi matrix equals a diagonal matrix plus an exponentially
decaying correction, from which it follows via eigenvalue perturbation theory
[12] that the eigenvalues (roots) $x_{k}$ converge exponentially to the
diagonal entries: $x_{k}\to kh$ as $\tau\to\infty$. In fact, a more careful
analysis shows that the convergence of $x_{k}$ to $kh$ is even faster than
first-order perturbation theory might suggest, due to cancellations between
first-order perturbation theory (which depends on the $kn+O(1/\tau)$ diagonal
elements) and second-order perturbation theory (which depends on the square of
the $O(1/\sqrt{\tau})$ off-diagonal elements). This rapid convergence is
illustrated in Fig. 3, which shows $x_{0}$ empirically converging as
$O(1/\tau^{N})$ for $J_{N}$ with $N=10$. A more precise asymptotic prediction
of $x_{k}-kh$ may be an interesting topic for future work, but in practice the
Gaussian summation scheme is most useful in the $\tau\approx 1$ ($hs\ll 1$)
regime where naive summation is slowly converging.

In the small-$hs$ regime, where naive summation requires a huge number of
terms $\sim 1/hs$, we will demonstrate below that MDL Gaussian summation
retains rapid convergence. As $hs\to 0^{+}$ ($\tau\to 1^{+}$), in fact, we
recover the Laguerre polynomials [13, 24] and the Gauss–Laguerre [22] points
and weights. In this limit, we have

$\displaystyle h\hat{\alpha}_{n}$
$\displaystyle\approx\frac{2n+1}{s}+O(h^{2})$ $\displaystyle h\hat{\beta}_{n}$
$\displaystyle\approx\frac{n+1}{s}+O(h^{2}),$

which tends to the coefficients of the three-term recurrence for the
continuous Laguerre polynomials [24] with a change of variables $x\to sx$ (to
change the weight function from $e^{-x}$ to $e^{-sx}$). The $O(h^{2})$
convergence of the MDL recurrence to the Laguerre recurrence is unsurprising
because the generalized MDL inner product Eq. 6 can be viewed as a trapezoidal
rule approximation for the continuous Laguerre inner product
$\int_{0}^{\infty}f(x)g(x)e^{-sx}dx$, and the trapezoidal rule converges as
$O(h^{2})$ [1].

Figure 4: Comparison of MDL quadrature convergence for the function
$F(x)=\cos{(x)}e^{-1.6x}$ via quadratures generated by
$\tau=e^{-s}=e^{-(1.6+\Delta s)}\in\\{e^{-0.8},e^{-1.6},e^{-3.2}\\}$ ($\Delta
s\in\\{-0.8,0,1.6\\}$). While convergence remains exponential in all cases,
optimal performance occurs at the correct asympototic decay rate $s=1.6$
($\Delta s=0)$.

###  3.2 MDL quadrature convergence

The convergence of Gaussian quadrature is well established for continuous
measures on bounded intervals [26] and on $[a,\infty)$ [29, 2]. Though Ref.
[29] does not directly address discrete measures, in Appendix B we adapt their
arguments to demonstrate the convergence of our MDL quadrature for any
polynomially bounded summand $F$, regardless of the decay rate $s$ used to
derive the quadrature rule.

Nevertheless, in cases where $F(x)$ is asymptotically exponentially decaying
proportional to $e^{-s_{0}x}$, as in many physical applications of Matsubara
summation, selecting the correct decay rate $s=s_{0}$ in the quadrature rule
is useful to obtain the fastest possible convergence. We illustrate this
property for a simple oscillating summand $F(x)=\cos{(x)}e^{-1.6x}$
($s_{0}=1.6$) in Figure 4, summed with $h=1$ ($\sum^{\prime}F(n)$). Here, we
constructed quadrature schemes using three different decay rates $s=0.8$,
$1.6$ ($=s_{0}$), and $3.2$. All three quadrature schemes converged for this
summand, but the scheme assuming the correct decay rate $s=s_{0}=1.6$ (red)
achieved superior convergence compared with either MDL summation with an
“incorrect” $s$ (blue and green) or direct summation (black).

##  4 Application to Casimir Forces

Casimir forces, a generalization of van der Waals forces (between particles)
to macroscopic objects, are an important interaction between electrically
neutral surfaces at sub-$\mu m$ separations and arise as a consequence of
quantum vacuum fluctuations in the electromagnetic fields [14, 21].
Physically, the net force is an integral $\int_{0}^{\infty}d\omega$ over
contributions from fluctuations at all frequencies $\omega$ (or wavelengths
$2\pi c/\omega$ where $c$ is the speed of light), weighted by a Bose–Einstein
distribution $\coth(\hbar\omega/2k_{B}T)$ where $\hbar$ is Planck’s constant,
$k_{B}$ is Boltzmann’s constant, and $T$ is the temperature [14, 8].
Computationally, it is useful to “Wick rotate” this integral (which is wildly
oscillatory on the real-$\omega$ axis) to the imaginary-frequency axis
$\omega=i\xi$, in which case the integral becomes a sum Eq. 1 over the
“Matsubara frequencies” $\xi_{n}=\frac{2\pi k_{B}T}{\hbar}n=hn$ (where $h$ the
spacing parameter of our notation in the previous sections, not Planck’s
constant). This bosonic Matsubara sum also arises in various other contexts of
statistical physics [23]. In the limit of zero temperature ($T\to 0^{+}$),
this sum becomes an integral $\int_{0}^{\infty}F(\xi)d\xi$ over the imaginary-
frequency axis.

The summand $F(\xi)$ in the Casimir force represents the contribution of
random current fluctuations at an frequency $\xi=-i\omega$ everywhere in the
materials [14] and requires the solution of Maxwell’s equations [8]. This
electromagnetic scattering problem has been formulated in a wide variety of
ways in the Casimir-force literature depending upon the geometry in question
[8], but one of the most general formulations uses a boundary-element method
(BEM) [19], which discretizes each surface into a generic mesh (e.g. a
triangular mesh) and formulates electromagnetism as a discretized scattering
problem $Mu=f$, where the matrix $M$ is a Galerkin discretization of an
integral operator based on the Maxwell Green’s function [7]. It turns out that
the Casimir summand/integrand $F$ is given by a simple trace formula involving
the matrix $M$ [19]. Therefore, the evaluation of $F$ requires computational
effort equivalent to that of solving Maxwell’s equations (in order to
factorize the matrix $M$ in BEM, or to obtain a similar information such as a
scattering matrix in other Casimir formulations), and it is desirable to
compute the Matsubara sum with as few $F$ evaluations as possible, especially
at low temperatures where the spacing between Matsubara frequencies is small.

The key to rapid evaluation of the Casimir Matsubara sum is to exploit the
asymptotic decay rate of the summand $F(\xi)$, and this decay turns out to be
very straightforward to obtain analytically. At imaginary frequencies,
Maxwell’s equations yield _exponentially decaying interactions_ , because
spherical waves $e^{-i\omega r/c}/r$ are transformed to decaying solutions
$e^{-\xi r/c}/r$, with an exponential decay rate proportional to $\xi$ [21].
This causes the Casimir summand to decay exponentially. More precisely, using
a Born approximation, one can easily show that the large-$\xi$ interactions
are dominated by a single “round trip” of waves from one object to the other
and back [15]. In consequence, for the force between two objects with minimum
separation $d$, the Casimir summand decays at least as fast as $O(e^{-2\xi
d/c})$ for large $\xi\gg c/d$. This asymptotic decay is illustrated for two
spheres with $d=1\,\mu m$ in Fig. 1.

The $\xi=0$ term in the Casimir Matsubara sum requires special care, because a
singularity arises in most formulations of Maxwell’s equations at zero
frequency, where the electric and magnetic fields decouple [4]. It turns out
that the correct zero-frequency Casimir contribution should be computed as a
limit $\xi\to 0^{+}$ [8]. Computationally, this limit can be obtained by
methods such as Richardson extrapolation [18]. A favorable property of MDL
quadrature Eq. 1 for this application, however, is that the quadrature scheme
_implicitly_ computes the $\xi\to 0^{+}$ limit, since it evaluates the summand
only at strictly positive quadrature points $x_{k}>0$, as illustrated in Fig.
3.

Figure 5: Relative error of the series, quadrature, and Monien approaches at
$T=300,3$ K.

###  4.1 Efficiency for Casimir summation

We compare three methods for approximating the Casimir Matsubara sum in Fig.
5, to evaluate their accuracy as a function of the number $N$ of summand
$F(\xi)$ evaluations. The first approach (blue lines) is naive summation
truncated to $N$ terms. The second method is our MDL quadature scheme of order
$N$ (red lines). The third method (green lines) is an alternative Gaussian-
quadrature scheme proposed in [16], which uses a quadratically decaying weight
function instead of our exponential weight function, so that it does not
exploit knowledge of the asymptotic decay rate of $F$.

Our test geometry consisted of two spheres of radius 1 $\mu$m separated by a
gap of length 0.2 $\mu$m, a typical separation for nanoscale applications of
Casimir forces. For both high ($T=300$ K) and low ($3$ K) temperatures, we
computed the estimated Casimir force between the spheres using the free/open-
source SCUFF-CAS3D BEM solver [19, 20]. In order to evaluate the summand at
the smallest $\xi d/c\ll 1$ point (where our BEM solver encountered numerical
difficulties for all three methods), we used Richardson extrapolation.

All three methods are viable in the high-temperature case $T=300$ K, altbough
naive summation is still noticeably less effective than the quadrature-based
approaches. As expected, however, our method shows superior performance at low
temperatures. At $T=3$ K, the naive summation method requires hundreds of
evaluations to obtain even one digit of accuracy, while the alternative
quadrature method performs only marginally better. Our method, meanwhile,
shows no increase in the number of required function evaluations.

##  5 Conclusion

Overall, our quadrature approach for computing this Matsubara sum is a
significant improvement over naive summation and even compared to Gaussian
quadrature schemes that do not exploit knowledge of the asymptotic decay rate.
This approach could also be applied towards any exponentially decaying
Matsubara sum (e.g. in density functional theory), or more generally to any
similar summation problem. More generally, given the MDL polynomials and an
exponentially decaying function $F(x)$ evaluated at the MDL roots, one could
construct a polynomial approximation that could be used to rapidly compute the
sum or any other functional of $F(nh)$ [27].

There are many possible avenues for future study. A related problem is to
evaluate the Matsubara sum at many different temperatures (different $h$) for
the same $F$, and in this case one would want to construct a _single_
polynomial approximation and re-use that for multiple $h$ values—it would be
interesting to consider what polynomial basis would be optimal for such a
problem. In our work, we used the Golub–Welsch algorithm to obtain Gaussian
quadrature points and weights, but it might also be interesting to explore
nested quadrature rules, such as Gauss–Kronrod or Gauss–Patterson rules [13],
in order to obtain error estimates and enable adaptive determination of the
required quadrature order. It would also be useful to develop tighter
theoretical estimates for the convergence rate of such quadrature schemes for
polynomially bounded, analytic summands on infinite intervals, since previous
such error bounds [26] typically apply only to quadrature on finite intervals.

##  Appendix A Proofs of Results 2.2–2.3

Here, we provide extended proofs of the important results in the text
pertaining to the MDL polynomials:

######  Proof A.1 (Proof of Lemma 2.3).

We proceed via strong induction. It is easy to see that when $n=0$, we have
$L_{0}^{\prime}(0)=1$ and $\langle
L_{0}^{\prime},L_{0}^{\prime}\rangle^{\prime}=\frac{\tau+1}{2\tau-2}=\frac{1+\tau^{1}}{\tau^{0}(1+\tau^{0})(\tau-1)}$.
Now, assume the statement for all $n^{\prime}<n$. By Lemma 2.1, we have

$L_{n}^{\prime}=L_{n}+\frac{1}{2}\frac{1}{\tau^{n}}\sum\limits_{i=0}^{n-1}\frac{L_{i}^{\prime}(0)L_{i}^{\prime}}{\langle
L_{i}^{\prime},L_{i}^{\prime}\rangle}=L_{n}+\frac{1}{\tau^{n}}\sum\limits_{i=0}^{n-1}\frac{\tau^{i}(\tau-1)L_{i}^{\prime}}{(1+\tau^{i+1})},$

which implies

$\displaystyle L_{n}^{\prime}(0)$
$\displaystyle=L_{n}(0)+\frac{1}{\tau^{n}}\sum\limits_{i=0}^{n-1}\frac{\tau^{i}(\tau-1)L_{i}^{\prime}(0)}{(1+\tau^{i+1})}$
$\displaystyle=\frac{1}{\tau^{n}}\left(1+\sum\limits_{i=0}^{n-1}\frac{2\tau^{i}(\tau-1)}{(1+\tau^{i})(1+\tau^{i+1})}\right)$
$\displaystyle=\frac{1}{\tau^{n}}\left(1+2\sum\limits_{i=0}^{n-1}\left(\frac{1}{1+\tau^{i}}-\frac{1}{1+\tau^{i+1}}\right)\right)$
$\displaystyle=\frac{1}{\tau^{n}}\left(1+2\left(\frac{1}{1+\tau^{0}}-\frac{1}{1+\tau^{n}}\right)\right)$
$\displaystyle=\frac{1}{\tau^{n}}\left(\frac{2\tau^{n}}{1+\tau^{n}}\right)=\frac{2}{1+\tau^{n}}.$

Similarly, by the orthogonality of the $L_{i}^{\prime}$, we deduce that

$\displaystyle\langle L_{n}^{\prime},L_{n}^{\prime}\rangle^{\prime}$
$\displaystyle=\langle L_{n}^{\prime},L_{n}\rangle^{\prime}$
$\displaystyle=\langle
L_{n}^{\prime},L_{n}\rangle-\frac{1}{2}L_{n}^{\prime}(0)L_{n}(0)$
$\displaystyle=\langle L_{n},L_{n}\rangle-\frac{1}{\tau^{n}(1+\tau^{n})}$
$\displaystyle=\frac{1}{\tau^{n-1}(\tau-1)}-\frac{1}{\tau^{n}(1+\tau^{n})}=\frac{1+\tau^{n+1}}{\tau^{n}(1+\tau^{n})(\tau-1)},$

as desired.

######  Proof A.2 (Proof of Theorem 2.4).

Note that the leading coefficient of $L_{n}^{\prime}$ is
$\left(\frac{1}{\tau}-1\right)^{n}\frac{1}{n!}$. This shows that
$\alpha_{n}=\frac{\tau-1}{(n+1)\tau}$. Furthermore, since $\langle
L_{n+1}^{\prime},L_{n-1}^{\prime}\rangle^{\prime}=0$, we have

$\displaystyle\gamma_{n}\langle
L_{n-1}^{\prime},L_{n-1}^{\prime}\rangle^{\prime}$
$\displaystyle=\alpha_{n}\langle
xL_{n}^{\prime},L_{n-1}^{\prime}\rangle^{\prime}$
$\displaystyle=\alpha_{n}\langle
L_{n}^{\prime},xL_{n-1}^{\prime}\rangle^{\prime}$
$\displaystyle=\frac{\alpha_{n}}{\alpha_{n-1}}\langle
L_{n}^{\prime},L_{n}^{\prime}\rangle^{\prime},$

implying that

$\displaystyle\gamma_{n}$
$\displaystyle=\frac{\alpha_{n}}{\alpha_{n-1}}\frac{\langle
L_{n}^{\prime},L_{n}^{\prime}\rangle^{\prime}}{\langle
L_{n-1}^{\prime},L_{n-1}^{\prime}\rangle^{\prime}}$
$\displaystyle=\frac{n}{n+1}\frac{(1+\tau^{n+1})(1+\tau^{n-1})}{\tau(1+\tau^{n})^{2}}.$

We then have

$\displaystyle\frac{2}{1+\tau^{n+1}}$ $\displaystyle=L_{n+1}^{\prime}(0)$
$\displaystyle=\beta_{n}L_{n}^{\prime}(0)-\gamma_{n}L_{n-1}^{\prime}(0)$
$\displaystyle=\beta_{n}\frac{2}{1+\tau^{n}}-\frac{2}{1+\tau^{n-1}}\frac{n}{n+1}\frac{(1+\tau^{n+1})(1+\tau^{n-1})}{\tau(1+\tau^{n})^{2}}$
$\displaystyle=\beta_{n}\frac{2}{1+\tau^{n}}-\frac{n}{n+1}\frac{2(1+\tau^{n+1})}{\tau(1+\tau^{n})^{2}},$

which implies

$\displaystyle\beta_{n}=(1+\tau^{n})\left(\frac{1}{1+\tau^{n+1}}+\frac{n}{n+1}\frac{(1+\tau^{n+1})}{\tau(1+\tau^{n})^{2}}\right)=\frac{1+\tau^{n}}{1+\tau^{n+1}}+\frac{n}{\tau(n+1)}\frac{1+\tau^{n+1}}{1+\tau^{n}}.$

##  Appendix B Convergence of MDL Quadrature

Here, we establish the convergence of our MDL quadrature (summation) scheme.
For simplicity, we consider the case of $h=s=1$, and prove convergence for any
summand $F(x)$ satisfying the asymptotic upper bound

$|F(x)|=e^{-x}|f(x)|=O\left(\frac{1}{x^{1+\epsilon}}\right)$

for some $\epsilon>0$. Our proof is an extension of the proof in Ref. [29].
Below, we will let $W(x)$ denote the function taking the value $\frac{1}{2}$
at $x=0$, $e^{-x}$ whenever $x\in\mathbb{Z}^{+}$, and $0$ elsewhere, and let
$w(x)$ denote the corresponding distribution
$w(x)=\sum_{n=0}^{\infty}W(n)\delta(x-n)=\sum_{n=0}^{\infty}{}^{\prime}\delta(x-n)e^{-n}$.

Ref. [29] establishes convergence based on properties of the even moments
$c_{2m}=\int x^{2m}w(x)dx$ of the weight function, and we can employ the same
strategy here. Since the function $x^{n}e^{-x}$ is increasing in $(0,n)$ and
decreasing in $(n,\infty)$, we can bound the even moments $c_{2m}$ for $m>0$
by breaking them into left and right Riemann sums as follows:

$\displaystyle c_{2m}=\sum\limits_{k=1}^{\infty}k^{2m}e^{-k}$
$\displaystyle=\sum\limits_{k=0}^{2m-1}k^{2m}e^{-k}+\sum\limits_{k=2m+1}^{\infty}k^{2m}e^{-k}+\left(\frac{2m}{e}\right)^{2m}$
$\displaystyle\leq\left(\frac{2m}{e}\right)^{2m}+\int\limits_{0}^{\infty}x^{2m}e^{-x}dx$
$\displaystyle=\left(\frac{2m}{e}\right)^{2m}+(2m)!<2(2m!)1^{2m}.$

This bound coincides with the condition given by Ref. [29] (p. 556) with $C=2$
and $R=1$. We then apply identical reasoning to Ref. [29] (p. 546) to show
that, for any non-integer $H$ and quadrature points and weights
$(x_{k},w_{k})$ generated by orthogonal polynomials of increasing order $n$,

$\lim\limits_{n\rightarrow\infty}\sum\limits_{x_{k}\geq
H}w_{k}=\int\limits_{H}^{\infty}w(x)dx=\sum_{k>H}W(k)\,.$

For any nonnegative integer $m$ and $0<\epsilon_{m}<1$, we therefore see that

$\lim\limits_{n\rightarrow\infty}\sum\limits_{x_{k}\in(m-\epsilon_{m},m+\epsilon_{m})}w_{k}=\int\limits_{m-\epsilon_{m}}^{m+\epsilon_{m}}w(x)dx=W(m)>0\,,$

and hence there must be at least one $x_{k}$ in any neighborhood of $m$ for
sufficiently large quadrature order $n$, with combined weight approaching
$W(m)$. When combined with the facts that each quadrature point $x_{k}$ must
be positive [3] and each interval $(m,m+1)$ between nonnegative integers must
contain at most one quadrature point [6, 17], we see that the MDL points and
weights necessarily converge to the nonnegative integers and $W(m)$,
respectively.

Finally, we may conclude using the same methods as Ref. [29] that the MDL
quadrature converges over the infinite interval: For any $\hat{\epsilon}>0$,
we may choose some non-integer $H$ such that

$\int\limits_{H}^{\infty}w(x)f(x)dx<\hat{\epsilon}\,.$

We then see that

$\lim\limits_{n\rightarrow\infty}\sum\limits_{x_{k}<H}w_{k}f(x_{k})=\int\limits_{0}^{H}w(x)f(x)dx$

by our previous result, and the convergence of the remaining integral
$\int\limits_{H}^{\infty}w(x)f(x)dx$ may be handled using the same technique
as Ref. [29] (p. 557). This method makes use of the bounds on the moments
$c_{2m}$ proven previously.

## Acknowledgments

This work was supported in part by the Undergraduate Research Opportunities
program at MIT. We are also grateful to Alan Edelman for helpful discussions.

## References

  * [1] K. E. Atkinson, An Introduction to Numerical Analysis, John Wiley & Sons, 2008.
  * [2] A. Bultheel, C. DÃaz-Mendoza, P. GonzÃĄlez-Vera, and R. Orive, On the convergence of certain Gauss-type quadrature formulas for unbounded intervals, Mathematics of Computation, 69 (2000), pp. 721–747.
  * [3] S. Engblom, Gaussian quadratures with respect to discrete measures, tech. report, Uppsala University, 2006.
  * [4] C. L. Epstein and L. Greengard, Debye sources and the numerical solution of the time harmonic Maxwell equations, Communications on Pure and Applied Mathematics, (2009), pp. NA–NA, https://doi.org/10.1002/cpa.20313, https://doi.org/10.1002/cpa.20313.
  * [5] G. H. Golub and J. H. Welsch, Calculation of Gauss quadrature rules, Mathematics of Computation, 23 (1969), pp. 221–230.
  * [6] M. J. Gottlieb, Concerning some polynomials orthogonal on a finite or enumerable set of points, American Journal of Mathematics, 60 (1938), pp. 453–458.
  * [7] R. F. Harrington, Time-Harmonic Electromagnetic Fields, IEEE Press, 1961\.
  * [8] S. G. Johnson, Numerical methods for computing Casimir interactions, in Casimir Physics, D. Dalvit, P. Milonni, D. Roberts, and F. da Rosa, eds., vol. 834 of Lecture Notes in Physics, Springer, Berlin, June 2011, ch. 6, pp. 175–218.
  * [9] M. Kaltak and G. Kresse, Minimax isometry method: A compressive sensing approach for Matsubara summation in many-body perturbation theory, Physical Review B, 101 (2020), https://doi.org/10.1103/physrevb.101.205145, https://doi.org/10.1103/physrevb.101.205145.
  * [10] A. A. Kananenka, A. R. Welden, T. N. Lan, E. Gull, and D. Zgid, Efficient temperature-dependent Green’s function methods for realistic systems: Using cubic spline interpolation to approximate matsubara green’s functions, Journal of Chemical Theory and Computation, 12 (2016), pp. 2250–2259, https://doi.org/10.1021/acs.jctc.6b00178, https://doi.org/10.1021/acs.jctc.6b00178.
  * [11] C. Karrasch, V. Meden, and K. Schönhammer, Finite-temperature linear conductance from the Matsubara Green’s function without analytic continuation to the real axis, Physical Review B, 82 (2010), p. 125114.
  * [12] T. Kato, Perturbation theory for linear operators, vol. 132, Springer Science & Business Media, 2013.
  * [13] P. K. Kythe and M. R. Schäferkotter, Handbook of Computational Methods for Integration, Chapman and Hall/CRC, 2004.
  * [14] E. M. Lifshitz and L. P. Pitaevskii, Statistical Physics: Part 2, Pergamon, Oxford, 1980.
  * [15] M. F. Maghrebi, Diagrammatic expansion of the Casimir energy in multiple reflections: Theory and applications, Physical Review D, 83 (2011), p. 045004.
  * [16] H. Monien, Gaussian quadrature for sums: A rapidly convergent summation scheme, Mathematics of computation, 79 (2010), pp. 857–869.
  * [17] T. Popoviciu, Sur la distribution des zeros de certains polynomes minimisants, Acad. RP Romîne, Bulletin de la Section Scientifique, 10 (1929), pp. 1–4.
  * [18] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, Cambridge Univ. Press, third ed., 2007.
  * [19] M. H. Reid, J. White, and S. G. Johnson, Fluctuating surface currents: An algorithm for efficient prediction of Casimir interactions among arbitrary materials in arbitrary geometries, Physical Review A, 88 (2013), p. 022514.
  * [20] M. T. H. Reid, “SCUFF-EM” software package.  https://github.com/HomerReid/scuff-em, 2020.
  * [21] A. W. Rodriguez, F. Capasso, and S. G. Johnson, The Casimir effect in microstructured geometries, Nature Photonics, 5 (2011), pp. 211–221.
  * [22] G. Sansone, Orthogonal functions, Interscience, 1959.
  * [23] G. Stefanucci and R. van Leeuwen, Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction, Cambridge University Press, 2013.
  * [24] G. Szegő, Orthogonal polynomials, vol. 23, American Mathematical Society, 1939.
  * [25] A. Townsend, The race for high order Gauss–Legendre quadrature, SIAM News, 48 (2015), pp. 1–3.
  * [26] L. N. Trefethen, Is Gauss quadrature better than Clenshaw–Curtis?, SIAM Review, 50 (2008), pp. 67–87.
  * [27] L. N. Trefethen, Approximation Theory and Approximation Practice, SIAM, 2019.
  * [28] L. N. Trefethen and D. Bau III, Numerical Linear Algebra, vol. 50, SIAM, 1997.
  * [29] J. V. Uspensky, On the convergence of quadrature formulas related to an infinite interval, Transactions of the American Mathematical Society, 30 (1928), pp. 542–559.


# Graph-based Prediction and Planning Policy Network (GP3Net) for scalable
self-driving in dynamic environments using Deep Reinforcement Learning

Jayabrata Chowdhury 1 , Venkataramanan Shivaraman 2 , Suresh Sundaram 1, P B
Sujit 2

###### Abstract

Recent advancements in motion planning for Autonomous Vehicles (AVs) show
great promise in using expert driver behaviors in non-stationary driving
environments. However, learning only through expert drivers needs more
generalizability to recover from domain shifts and near-failure scenarios due
to the dynamic behavior of traffic participants and weather conditions. A deep
Graph-based Prediction and Planning Policy Network (GP3Net) framework is
proposed for non-stationary environments that encodes the interactions between
traffic participants with contextual information and provides a decision for
safe maneuver for AV. A spatio-temporal graph models the interactions between
traffic participants for predicting the future trajectories of those
participants. The predicted trajectories are utilized to generate a future
occupancy map around the AV with uncertainties embedded to anticipate the
evolving non-stationary driving environments. Then the contextual information
and future occupancy maps are input to the policy network of the GP3Net
framework and trained using Proximal Policy Optimization (PPO) algorithm. The
proposed GP3Net performance is evaluated on standard CARLA benchmarking
scenarios with domain shifts of traffic patterns (urban, highway, and mixed).
The results show that the GP3Net outperforms previous state-of-the-art
imitation learning-based planning models for different towns. Further, in
unseen new weather conditions, GP3Net completes the desired route with fewer
traffic infractions. Finally, the results emphasize the advantage of including
the prediction module to enhance safety measures in non-stationary
environments.

## Introduction

Motion planning for AVs still has a long way to go due to the complexity of
real-world driving environments, ensuring the safety and comfort of everyone.
The AV must be able to carefully plan its maneuvers in urban and highway
environments with varying traffic dynamics set by cars, bikers, pedestrians,
and weather conditions. A typical intersection scenario in an urban
environment is shown in Fig.1.

Figure 1: A typical intersection scenario with traffic lights, pedestrians,
and vehicles. The scenario has pedestrians and vehicles moving on the desired
path of AV. The figure has colour coded markers to show the context.

Various traffic participants’ different intentions and final goals (as shown
by arrows in Fig.1) make it challenging for AV to safely and comfortably move
to its destination. Variable speeds of different traffic participants also
make it difficult for AV to drive safely. A research finding on the United
States (Choi 2010) shows that wrong estimation of other traffic participants’
speed is responsible for $8.4\%$ of critical reasons for crashes. Also,
another problem of driving in an unstructured environment is that traffic
pattern varies with the cooperative/noncooperative behaviors of other traffic
participants. By utilizing communication channels as described in (Bazzi et
al. 2021; Cui et al. 2022), it becomes possible to determine the intentions
and objectives of fellow traffic participants. This information simplifies the
process for AVs to make informed decisions regarding their actions. However,
vehicles use different communication protocols, and ensuring reliable real-
time communications between all traffic participants is challenging. Hence,
another way is to model the interactions of traffic participants through time
and predict their future trajectories to understand their behaviors. An AV can
make a safe and efficient decision in a non-stationary environment if it can
model the interactions better and predict the future trajectories of other
surrounding vehicles. Depending on the future trajectories of others’, an AV
can decide on a safe planning maneuver. Hence, there is a need to develop a
motion planning algorithm that can handle non-stationary environments using
the intention understanding from the predicted trajectories of other traffic
participants.

Existing state-of-the-art methods can be broadly classified into rule-based,
imitation learning-based, and reinforcement learning-based works. Recently
(Aksjonov and Kyrki 2021) uses a rule-based decision-making system for
navigating road intersections. However, due to the non-stationary nature of
the environment, rules vary for each situation, and attempting to create
generalized rules across diverse environments suffers from scalability issues.
On the other hand, some models mimicking expert drivers, called Imitation
Learning (IL) based motion planner, has been proposed in (Bansal, Krizhevsky,
and Ogale 2018; Rhinehart, McAllister, and Levine 2020; Codevilla et al. 2018;
Cai et al. 2020, 2021). The driving model in (Bansal, Krizhevsky, and Ogale
2018) learns driving decisions from expert human driver demonstrations. The
work in (Rhinehart, McAllister, and Levine 2020) uses expert driving data
collected from CARLA (Dosovitskiy et al. 2017) simulator to learn expert
behaviors. An end-to-end motion planning algorithm has been learned with IL in
(Cai et al. 2020; Codevilla et al. 2018). However, the IL method suffers from
a data distribution shift from training scenarios to unseen evaluation
scenarios as (Bansal, Krizhevsky, and Ogale 2018) did not perform well in
unknown highly interactive lane change scenarios. In (Cai et al. 2021), a
graph-based attention network has been used to model participant interactions,
and a driver model learns to imitate the expert with IL. In another work (Teng
et al. 2023), a Bird’s Eye View (BEV) representation has been generated from
the front camera view to use with IL for decision-making. The expert bias and
distribution shift problems heavily influence the aforementioned IL-based
works’ performance. Since there are no near-collision scenarios in expert
driving, it will be difficult for IL-based methods to perform in the real
environment.

Another popular approach to overcome the expert-based IL method is
Reinforcement Learning (RL). In RL, the agent learns through their experience
of completing the task. Recently works (Chen, Yuan, and Tomizuka 2019; Chen,
Li, and Tomizuka 2021; Ye et al. 2020; Tang et al. 2022; Alizadeh et al. 2019;
Naveed, Qiao, and Dolan 2021) employed RL for motion planning and maneuver
decision-making for autonomous driving. In (Chen, Yuan, and Tomizuka 2019;
Chen, Li, and Tomizuka 2021), a BEV image-based representation has been used
to represent the AV’s surrounding environment. An RL agent learns to make
maneuver decisions with this input representation. The works in (Ye et al.
2020; Tang et al. 2022) learn lane-changing maneuvers in simulated handcrafted
driving scenarios through policy optimization methods ((Schulman et al. 2017;
Haarnoja et al. 2018)). Recently, a work (Chen et al. 2020) used conditional
Deep Q Network (CDQN) for stable driving decisions in the CARLA environment.
However, this work did not include other vehicles, pedestrians, or traffic
lights to model spatio-temporal interactions with AV. The works mentioned
above need to recognize the intent of surrounding traffic participants. In
(Huang et al. 2019), the intent of surrounding vehicles is calculated
probabilistically to solve motion planning problems. However, the intent is
modeled for a fixed number of vehicles without pedestrians and evaluated in a
handcrafted driving scenario. The intentions of varying numbers of surrounding
vehicles should be inferred to understand how the non-stationary environment
can evolve. Hence, a predictive motion planning algorithm is required to model
the spatio-temporal interactions for trajectory prediction and make a safe and
efficient maneuver decision.

The proposed deep Graph-based Prediction and Planning Policy Network (GP3Net)
framework finds a way to handle distribution shifts, model the interactions to
predict trajectories and plan safe maneuvers. The GP3Net’s benefits are as
follows:

  1. 1.

A deep spatio-temporal graphical model encodes interaction to depict the
behaviors of traffic participants with AV. These modeled interactions are
passed through the trajectory prediction module to provide predicted future
trajectories with associated prediction uncertainties.

  2. 2.

The future mask generation module generates a BEV occupancy map with the
prediction uncertainties, capturing the intentions of other traffic
participants and probable evolution of a non-stationary environment in the
future.

  3. 3.

The AV neural architecture (state-encoder plus policy network) encodes past
BEV masks and predicted future BEV masks together to give the final state.

  4. 4.

The AV neural architecture outputs the values of acceleration and steering
directly. The whole AV neural architecture is updated end-to-end with the RL
algorithm.

  5. 5.

The learned GP3Net driving agent outperforms the previous model DiGNet (Cai et
al. 2021) and HIIL (Teng et al. 2023) in CARLA Leaderboard and No-Crash motion
planning benchmark in different towns and weathers. A Qualitative analysis
highlights the importance of the graph-based trajectory prediction module.

## Related Work

### Trajectory prediction

Trajectory prediction is a crucial part of autonomous driving involving
spatiotemporal interactions between traffic participants. In (Ma et al. 2019;
Xu, Yang, and Du 2020), Long-Short Term Memory (LSTM) network has been used
for modeling. However, these works need to model the influence of one traffic
participant on another. Recently, (Shi et al. 2022) proposes a tree-based
attention method. Also, in (Li and Zhu 2021), spatio-temporal modeling has
been by a graphical model. However, these approaches needed to provide
uncertainties in trajectory prediction. In (Ivanovic and Pavone 2019), a
spatio-temporal dynamic graph models interactions, and LSTM networks model the
state encoding of graph nodes. A Conditional Variational Auto-Encoder (C-VAE)
learns the trajectory distribution. Samples from this probability distribution
are taken for Gaussian Mixture Model (GMM)-based multimodal trajectory
prediction. Trajectron predicts future trajectories with prediction
uncertainties. GP3Net uses Trajectron architecture for trajectory prediction
of surrounding traffic participants.

### End-to-end and semantic-based motion planning

Many works directly use sensors (camera, LIDAR) to understand the surrounding
environment and make decisions. In (Codevilla et al. 2018), a directional
high-level command is used to learn low-level steering and throttle values.
The following work (Codevilla et al. 2019) showed the limitations of IL for
autonomous driving. Recently (Teng et al. 2023) uses a front camera as a
vision sensor and converts it to BEV representation of the driving scene.
However, this work does not model rear-end side of the AV in the BEV. The
rear-end side is crucial for AV because $64.2\%$ of AV crashes happen from the
rear end as reported in (Petrović, Mijailović, and Pešić 2020). GP3Net uses
the semantic BEV representation for deciding safe maneuvers.

Figure 2: Schematic diagram of the GP3Net framework. The diagram shows the
pipeline from input observations O to control outputs $a$ and $\omega$,
highlighting the future occupancy map generation by the trajectory prediction
module, the state-encoder module, the PPO algorithm that optimizes the AV’s
state-encoder and policy network.

## The GP3Net Framework

This section introduces the GP3Net motion planning framework that incorporates
future trajectory predictions generated by a trajectory prediction module into
Deep Reinforcement Learning-based motion planning. The subsections describe
the self-driving task as a Partially Observed Markov Decision Process (POMDP).
The primary input for the GP3Net consists of past and predicted future BEV
occupancy maps containing information about road layout, vehicles,
pedestrians, and traffic lights. The predicted future occupancy maps of
vehicles are generated using a spatio-temporal graph-based trajectory
prediction module based on (Ivanovic and Pavone 2019). Fig. 2 is a flowchart
showing the pipeline between input processing, future BEV occupancy map
generation, and output control actions.

### POMDP Formulation

POMDP is defined for this task as a 6-tuple $(S,A,T,R,O,\gamma)$. Here, $S$ is
a set of possible states in the environment, $A$ is the action space that the
agent/AV can take in the environment, $T$ is the transition probability
$T(s^{\prime}|s,a)$ that represents the probability of AV moving to state
$s^{\prime}$ after taking action $a$ from state $s$. $R$ is the reward
function $R:S\times A\rightarrow\mathbb{R}$ for the AV, $O$ is a set of
observations the AV can get, and $\gamma\in(0,1]$ is the discount factor. The
goal is to learn a policy $\pi$ which can maximize the expected cumulative
reward $\mathbb{E}[\sum_{t=0}^{\infty}\gamma^{t}r_{t}]$, where $r_{t}$ is the
reward obtained in the time step $t$. Following subsections define the
observation space, action space, and reward function for the driving task in
CARLA simulator.

### Input Preprocessing and Observation Space Design

The observation space is $\textbf{O}=[O_{c};O_{p};O_{f};O_{m}]$. This
observation vector O for RL-based AV consists of several components of
particular categories. The observations $O_{c},O_{p}$ and $O_{m}$, except
$O_{f}$, are directly obtained from environment. The observations
$O_{c},O_{p}$ and $O_{f}$ are rasterized, segmented local BEV image/mask
centered around the ego vehicle. The $O_{c}$ contains segmented BEV masks of
the road layout, the mission route, and lane boundaries for contextual
localization at time-step $t$. Observation $O_{p}$ contains past
configurations masks of the surrounding vehicles, pedestrians, traffic lights,
and stop signs for $k$ time-steps in the past, each $4$ time-steps apart.
Observation $O_{m}=[\tau,b,\omega,\dot{\textbf{x}}]$ is a simple vector
containing the AV’s odometric measurements such as throttle $\tau$, braking
$b$, steering $\omega$, and velocity $\dot{\textbf{x}}=(\dot{x},\dot{y})$ at
time-step $t$ as shown in Fig. 2. The vector $O_{f}$ contains the predicted
occupancy maps (added as masks) of surrounding vehicles for $l$ time-steps in
the future, spaced $4$ time-steps apart. These future occupancy maps in
$O_{f}$ are accurately generated by a spatio-temporal graph-based trajectory
prediction module. This trajectory prediction model takes in the observed past
trajectories of surrounding vehicles/pedestrians and models their interactions
as dynamic spatiotemporal graphs to accurately predict future trajectories.
The future occupancy maps are generated using the future trajectories with
respect to our AV.

Graph-based Trajectory Prediction Module: The trajectory prediction module
takes past kinematics information of the surrounding vehicles as input to
generate a distribution of future trajectories for all agents present in a
scene. The most probable future trajectory is then selected for generating the
future maps. This module doesn’t assume a fixed number of vehicles in the
sensor range of the AV. Let $\textbf{x}^{i}_{(t-k,\
t)}=[(x,y)_{t-k},(x,y)_{t-(k-1)},...,(x,y)_{t}]$ are the spatial coordinates
of agent $i$ from timestep $t-k$ to $t$. The objective is first to predict the
future trajectories
$\textbf{x}^{1}_{(t+1,t+T)},\textbf{x}^{2}_{(t+1,t+T)},...,\textbf{x}^{N}_{(t+1,t+T)}$
for $T$ timesteps for all $N$ agents which were present in the scene in that
time window. Formally, the input for the model is
$\textbf{x}=\left[\textbf{x}^{1,..,N}_{1:t_{obs}},\dot{\textbf{x}}^{1,..,N}_{1:t_{obs}},\ddot{\textbf{x}}^{1,..,N}_{1:t_{obs}}\right]$
and the output $\textbf{y}=\dot{\textbf{x}}^{1,..,N}_{(t_{obs}+1,t_{obs}+T)}$.
The module contains a dynamic spatio-temporal graph model, deep generative
CVAE with LSTM-based encoder/decoder, and GMM. This architecture directly
models the velocity rather than the trajectories. When the future velocities
are obtained, they are numerically integrated to produce future trajectories
using single-integrator dynamics, as given below. Here, $\delta t$ is the
difference of consecutive timesteps.

$\textbf{x}^{i}_{t+1}=\textbf{x}^{i}_{t}+\textbf{y}^{i}_{t}.\delta t$ (1)

The dynamic spatiotemporal graph representation effectively models dynamic
interactions between the agents present in the scene. The graphs are created
based on the spatial proximity of the agents. The CVAE framework attempts to
learn the probability distribution function $p(\textbf{y}|\textbf{x})$ as
given below, where $z$ is a discrete latent variable.

$p(\textbf{y}|\textbf{x})=\sum_{z\in
Z}p_{\psi}(\textbf{y}|\textbf{x},z)p_{\phi}(z|\textbf{x})$ (2)

The encoder architecture consists of a node history encoder, a node future
encoder (only during training), and an edge encoder.

$e_{t,node}^{i}=LSTM(e_{t-1,node}^{i},\textbf{x}^{i}_{t};W_{NE})$ (3)

Both graph node history and future encoders have LSTM networks, as shown in
Equ.3 for encoding node data and temporal relations between them.

$f_{t}^{i}=\left[\textbf{x}_{t}^{i};\sum_{j\in
N(i)}\textbf{x}_{t}^{j}\right];a_{t}^{i}=LSTM(a_{t-1}^{i},f_{t}^{i};W_{EE})$
(4)

The edge encoder also uses LSTM network as given in Equ.4 and a modulating
filter to encode the influence of other agents with smooth edge additions and
removals if high-frequency changes occur in the graph. In Equ.4, $[a;b]$ is
concatenation, $N(i):$ set of neighbors of agent $i$. The latent variable $z$
is then sampled and fed to the decoder consisting of an LSTM network and GMMs.
The outputs for the decoder’s LSTM network are the GMMs’ parameters. The GMMs
are then used to produce the future velocities and trajectories. The future
masks containing future occupancy maps for the surrounding vehicles are
generated after obtaining the predicted future trajectories for all the
vehicles inside AV’s sensor range.

Future Mask Generation: After getting the predicted future trajectories
$\textbf{x}^{1,...,N}_{t+1,t+T}$, the cartesian coordinates
$\textbf{x}^{i}_{t+1,t+T}$ are converted to pixel coordinates
$\textbf{p}^{i}_{t+1,t+T}$ centered around the AV for each vehicle/agent $i$.
The existence of a vehicle in the future mask is denoted by a $2D$ symmetric
Gaussian patch as shown in Fig.2. The pixel coordinates p serve as the mean of
the 2D symmetric Gaussian function. The symmetric Gaussian approach helps us
account for the future prediction uncertainties of the vehicles, accounting
for various factors such as sensor noise and occlusion. Equation of future
vehicle patch is given by:

$g(x,y)=exp\left({\frac{((x-x_{o})^{2}+(y-y_{o})^{2}))}{2\sigma^{2}}}\right)$
(5)

In the above equation, $\sigma$ is the standard deviation of the Gaussian,
$x_{o}$ and $y_{o}$ are the mean of the 2D Gaussian. For each timestep $t_{f}$
in future prediction horizon $T$ of the trajectory prediction model, a future
mask is generated using pixel coordinates $\textbf{p}^{1,...,N}_{t_{f}}$ of
all $N$ vehicles in the method mentioned above. Therefore, $O_{f}$ consists of
$l$ future masks. The observations are taken $4$ timesteps apart instead of
$1$, and the inputs to the trajectory prediction module from timesteps $t-k$
to $t$ are passed in $4$ timestep intervals but processed as though inputs are
of consecutive timesteps. This is further explained in the experiments section
of the paper.

Final State Encoder: Once $O_{c}$, $O_{p}$, $O_{f}$, and $O_{m}$ are obtained,
they are passed through a state-encoder (first part of AV Neural Architecture)
which results in a column vector. The architecture for the state encoder, as
shown in Fig. 2, has two parts. The first part processes the BEV image masks
$O_{c},O_{p}$ and $O_{f}$ using a Convolutional Neural Network (CNN) as input
channels and a simple Fully Connected (FC) layer. The measurement vector
$O_{m}$ is encoded using two FC layers and concatenated with the output of the
BEV encoder, passed through an FC layer, finally giving the final state
representation for the RL algorithm.

### Action Space

The CARLA-Gym simulator takes three values for actions: throttle $\tau$, brake
$b$, and steering $\omega$. However, to train the RL algorithm, action space
is defined as a 2D continuous space
$\textbf{a}=[\max\\{\tau,b\\}\cdot\operatorname{sgn}(\tau-b),\ \omega]$, where
sgn(.) is signum function. The first expression in a is acceleration whose
value ranges from -1 (maximum braking) to 1 (maximum throttle), and the
steering value goes from -1 (maximum left steering) to 1 (maximum right
steering). If the acceleration value is positive, it is passed to throttle
$\tau$, and brake value $b$ is zero. If the acceleration value is negative,
the magnitude of the acceleration is passed as brake $b$ value, and throttle
$\tau$ is zero. To select actions from this continuous space, the policy
network $\pi_{\theta}$ is trained to output parameters $\alpha$ and $\beta$
for a Beta distribution $\mathcal{B}(\alpha,\beta)$ for an improved policy
gradient method (Chou, Maturana, and Scherer 2017). The throttle and steering
actions are sampled from this Beta distribution. This work uses separate Beta
distributions for throttle and steering.

### Reward Structure

The reward structure of RL is essential to train agents/AV to perform specific
tasks. The reward is a superposition of different components, as given below.

$R=r_{route}+r_{halt}+r_{vel}+r_{pos}+r_{hd}+r_{act}$ (6)

where the AV is penalized, and the episode is terminated whenever the vehicle
collides with anything, runs a traffic light or a stop sign, routes deviation,
or blocks. $r_{route}=-1$ if route deviation $\delta\geq 3.5m$; $r_{halt}=-1$
if AV’s velocity $v\leq 0.1m/s$; $r_{vel}=1-|v-v_{assigned}|/v_{max}$ where
$v_{max}$ is maximum velocity allowed and $v_{assigned}$ is assigned by the
simulator; $r_{pos}=-0.5x$ where $x$ is the distance from the AV’s center and
the center of the desired route; $r_{hd}=-\Delta_{h}$ absolute difference
between heading of AV and route; $r_{act}=-0.1$ if
$|\omega_{t}-\omega_{t-1}|>0.01$.

### PPO Algorithm

The AV’s policy network $\pi$ (second part of AV neural architecture) learns
using a state-of-the-art policy gradient-based RL algorithm called Proximal
Policy Optimization (PPO) (Schulman et al. 2017). The usual form of loss
function used in policy gradient methods is shown below.

$L_{pg}=\mathbb{E}_{\pi}\left[\hat{A}_{t}\log\pi(a_{t}\,|\,s_{t},\theta)\right]$
(7)

Here, $\pi_{\theta}(a_{t}\,|\,s_{t},\theta)$ is a stochastic policy with
$\theta$ as parameter. $\hat{A}_{t}$ is an estimator of the advantage function
at time $t$, and $\mathbb{E}_{\pi}$ represents the empirical average of a
finite batch of samples.

However, updating the policy network $\pi_{\theta}$ for optimizing the
objective loss function shown in Equ.8 can lead to large policy updates, often
resulting in instability and bad policy updates. A new objective function
called the clipped surrogate objective, which places a constraint on policy
updates, is proposed in PPO paper. This makes the training process much more
stable and reliable when updated over multiple epochs of gradient ascent.
Also, a Generalized Advantage Estimate (GAE) and adding an entropy bonus with
their new clipped surrogate objective function is used, as shown in the
following equation for PPO training of the GP3Net model.

$\displaystyle\mathcal{L}_{ppo}=\mathbb{E}\Big{[}\min(\tilde{r}(\theta)\hat{A},\mbox{clip}(\tilde{r}(\theta),1-\epsilon,1+\epsilon)\hat{A}),$
(8) $\displaystyle-
c_{1}(V_{\phi}(s,a)-V_{targ})-c_{2}H(s,\pi_{\theta})\Big{]}$

Here in Equ.8, $\epsilon,c_{1}\mbox{ and }c_{2}$ are tunable hyper-parameters
and $H(s,\pi_{\theta})$ denotes the entropy bonus. The generalized advantage
estimate makes use of a learned state function $V_{\phi}(s,a)$ for computing
the advantage function $\hat{A}$. The entropy bonus encourages the AV to
explore adequately. Most importantly, during update, the whole AV Neural
Architecture gets updated.

## Experiments

### Training Setup

This subsection describes the trajectory prediction module’s training setup
and the RL-based motion planning training.

#### Trajectory Prediction Module Training

First, trajectory data of surrounding vehicles are collected as observed by
the CARLA simulator’s roaming agent in different towns for training the
trajectory prediction module. The trajectory data contains rectangular
coordinates $\textbf{x}_{t}=(x_{t},y_{t})$ of all vehicles at time-steps
$t,t+4\delta t,t+8\delta t...$ with each time-step spaced $\delta t=0.1$
seconds apart in simulation time. The trajectory prediction module is provided
with 3.2 seconds of observed trajectories as input, divided into eight
discrete time-steps $(t-29\delta t,t-25\delta t,...,t-\delta t)$, and the
model predicts the trajectories up to 2.8 seconds, divided into seven time-
steps $(t+4\delta t,t+8\delta t,...,t+28\delta t)$.

The training starts by initializing the weights and setting the batch size to
$64$ and the learning rate to $0.001$ with a decay rate of $0.9999$ to improve
convergence. Adam optimizer updates the weight parameters of the model.
Trajectory prediction models are typically benchmarked by comparing their
performance against actual or simulated ground truth trajectories to evaluate
accuracy and reliability. Here Mean Squared Error (MSE) metric evaluates the
performance and accuracy of the model. The MSE is the mean $L_{2}$ error
between the ground truth trajectories and predicted trajectories. The model
was trained over 2000 steps and the MSE was evaluated intermittently over a
validation set as shown in (Ivanovic and Pavone 2019).

### RL Simulation Environment Setup

The CARLA simulator is employed as a gym environment for training the RL
algorithm and evaluating proposed RL-based motion planning. The training was
conducted in parallel on six distinct town maps, namely Town 1 to Town 6,
employing a vectorized environment configuration. This training approach on
multiple maps simultaneously enhanced the replay buffer’s diversity and
improved the training process’s efficiency. During the rollout phase, the AV
is provided with a mission route at the beginning of each episode. This route
is generated by randomly selecting a source and destination point and
connecting them using a path search algorithm like $A^{*}$. The AV aims to
successfully navigate the environment map and complete the assigned mission
route. The episode would terminate if any of the following conditions were
met: $1)$ collision with a vehicle, pedestrian, or obstacle in the map layout;
$2)$ deviation from the given route; $3)$ traffic stagnation/halt and $4)$
violation of traffic lights and stop signs. Upon reaching the destination
point, the AV was assigned a new mission route, continuing the rollout phase.
The training and experiments were run in parallel on a Ubuntu 20.04 machine
with 64 GB RAM and used two Nvidia 2080Ti GPU.

#### RL-based Motion Planning Training

The RL-based motion planning model is trained using the CARLA simulator
(Dosovitskiy et al. 2017) as an RL gym environment. The environment setup is
vectorized, and the model trains by running the RL training simulation on
different town maps in parallel. The training process consists of two phases:
(a) the rollout phase and (b) the training phase. In the rollout phase, the AV
keeps interacting with the environment for several episodes, trying to
navigate around the environment and reach its destination. At each step, the
AV collects the tuple $(\textbf{O},\textbf{a},r,\textbf{O}^{{}^{\prime}},d)$
and stores it in its replay buffer. Here, O and $\textbf{O}^{{}^{\prime}}$ are
observations of timestep $t$ and next timestep $t^{{}^{\prime}}$, $a$ is
action, $r$ is reward and $d$ is done/termination condition. The AV mainly
collects experience from various scenarios offered by different town maps. The
AV remains in this phase until it has accumulated sufficient experience in the
replay buffer for moving into the training phase. The training phase mainly
involves the AV updating its policy network $\pi_{\theta}$ based on the
experience it accumulated in the rollout phase. The policy $\pi_{\theta}$ is
updated based on the loss function Equ.8, and the AV again reverts to the
rollout phase. After a few rounds, the AV’s performance is evaluated, and
several metrics are recorded to ensure that the AV is learning correctly and
troubleshooting issues with proper tuning.

## Performance Evaluation

### Evaluation Environment Setup

To ensure a comprehensive evaluation of GP3Net, the performance benchmarking
environment is setup based on the Leaderboard benchmarking suite described in
the literature (Dosovitskiy et al. 2017) and NoCrash environment (Codevilla et
al. 2019). These benchmarking suites offer diverse, complex scenarios across
several towns, encompassing various weather conditions, numerical
configurations, and behavior of vehicles and pedestrians. GP3Net’s performance
is assessed and recorded after conducting several episodes involving different
environment configurations to ensure that the model can generalize well in
complex unseen scenarios.

### Performance Evaluation Metrics

|  | Town 03
---
urban
| Town 04
---
mixed
| Town 05
---
urban
| Town 06
---
highway
| Town 01
---
urban
| Town 02
---
urban
Models | SR $\uparrow$ | DS $\uparrow$ | SR $\uparrow$ | DS $\uparrow$ | SR $\uparrow$ | DS $\uparrow$ | SR $\uparrow$ | DS $\uparrow$ | SR $\uparrow$ | DS $\uparrow$ | SR $\uparrow$ | DS $\uparrow$
CILRS | 43.2 | 0.56 | 16.2 | 0.28 | 38.2 | 0.52 | 35.2 | 0.48 | 42.3 | - | 24.7 | -
MLP | 62.3 | 0.72 | 60.5 | 0.70 | 60.0 | 0.69 | 72.2 | 0.79 | - | - | - | -
GCN (U) | 68.8 | 0.77 | 68.5 | 0.76 | 70.2 | 0.78 | 73.8 | 0.80 | - | - | - | -
GCN (D) | 74.2 | 0.80 | 64.0 | 0.72 | 71.5 | 0.79 | 74.8 | 0.81 | - | - | - | -
GAT | 76.0 | 0.82 | 71.0 | 0.78 | 81.5 | 0.87 | 78.2 | 0.83 | - | - | - | -
DiGNet (CTL) | 72.8 | 0.79 | 67.5 | 0.75 | 69.8 | 0.77 | 73.2 | 0.79 | - | - | - | -
DiGNet (Full) | 80.2 | 0.85 | 67.5 | 0.75 | 69.8 | 0.77 | 78.2 | 0.79 | - | - | - | -
GP3Net (ours) | 82.5 | 0.87 | 75.0 | 0.85 | 82.8 | 0.91 | 82.5 | 0.85 | 92.5 | 0.96 | 93.1 | 0.91
Table 1: A comparative result showing the mean success rate and the driving
score of GP3Net and other SOTA works on different towns.

Specific performance evaluation metrics are employed to measure the
performance of the predictive reinforcement learning-based motion planning
framework GP3Net by the benchmarking suites aforementioned in the previous
subsection. These metrics include the Success Rate (SR) and Driving Score
(DS). The SR is calculated as the percentage of instances where the AV reaches
its intended destination without experiencing any collision. This metric
directly assesses the model’s ability to navigate and complete the designated
mission routes safely and successfully. The DS is a composite metric that
combines the AV’s route completion success with penalties incurred due to
collisions and traffic rule infractions. It is formulated as the product of
the percentage of completed routes and a penalty factor.

### Quantitative Benchmark Results

The GP3Net’s performance was evaluated in different towns and weather
conditions for $100$ episodes with five different seeds. Here Table 1 and
Table 2 display the mean comparative results with other previously proposed
solutions, some of which use a BEV perspective for motion planning similar to
GP3Net. It can be observed in Table 1 that the predictive RL-based motion
planner, GP3Net, outperforms DiGNet (Cai et al. 2021), its variations, and
other works like CILRS (Codevilla et al. 2019) in terms of the success rate
and the driving score in several towns. The mean improvement in the success
rate and driving score is $3.85\%$ and $8\%$, respectively. The standard
deviation is $1.5\%$. GP3Net has a large improvement in Town 4 and Town 6
which consists of both mixed and highway environments, showing capability to
handle well in higher speeds. This generalization to a wide range of scenarios
can be attributed to the GP3Net’s understanding of how a non-stationary
environment can evolve with the trajectory prediction module.

---

---

Figure 3: Plots of (a) reward obtained and (b) vehicle collisions with and without trajectory prediction module. | Weather
---
Condition
| T1
---
(ours)
| T2
---
(ours)
| T1
---
(HIIL)
| T2
---
(HIIL)
Wetnoon | 96 | 95.9 | 84 | 84
ClearSunset | 97.3 | 92 | 88 | 88
ClearNoon | 92 | 94 | 89 | 89
HardRainNoon | 96 | 88 | 80 | 80
SoftRainSunset | 93.3 | 94 | 86 | 86
WetSunset | 98 | 94 | 85 | 85
WetCloudSunset | 98 | 92 | 84 | 84
SoftRainNoon | 88 | 91.4 | 88 | 88
Table 2: Success-rate Benchmarking results on T1: Town 1 and T2: Town 2 in
different weather conditions.

Another experiment was performed to test the resilience of GP3Net to different
weather conditions, whose results are shown in Table 2. The test was conducted
in Towns 1 and 2 under eight weather conditions mentioned in Table 2.
According to Table 2, GP3Net is robust to rainy weather conditions and does a
better job of navigating the environment successfully compared to another
work, HIIL (Teng et al. 2023). The insights about the future trajectories
provided by the proposed trajectory prediction module have positively impacted
the RL-based motion planner’s decision-making during training and evaluation
time, making it safe and cautious of the non-stationary environment. The
following section discusses the effectiveness of our graph-based prediction
module.

### Qualitative Analysis

During training experiments, we plotted quantities such as total reward and
the number of vehicle collisions for monitoring the training process for two
cases: 1) With the trajectory prediction module and 2) without the trajectory
prediction module. Fig.3 shows a significant difference in the training plots.
After the $5$ M steps, both the models have learned correctly not to deviate
from the route or to stop when the lights turn red. However, the models were
yet to learn to avoid collisions with other vehicles. The model with a
trajectory prediction module showed a better capability to avoid collisions as
the training progressed compared to the one which did not use a trajectory
prediction module.

Ablation Study: In addition to the leaderboard suite, GP3Net’s performance was
benchmarked on the NoCrash suite for analyzing the framework’s behavior in
dense environments. Several metrics are utilized, such as end reach %, success
rate %, penalty factor, crash rates, vehicle halts, and traffic lights passed
were recorded to perform this analysis. The quantitative evaluation was done
in Towns 1 and 2 in training and testing weather conditions which are
displayed in Table 3. It can be observed that the success rate of GP3Net in
Town 2, though competitive, is not as high as in Town 1 due to a much denser
population of vehicles and walkers since Town 2 is smaller in size compared to
Town 1. It is to be noted that the NoCrash suite has twice the walker
population than that in the Leaderboard suite.

| Performance
---
Metric
| T1
---
train
| T1
---
new
| T2
---
train
| T2
---
new
End Reach (%) | 100 | 100 | 96 | 90
Succ. Rate (%) | 96 | 98 | 94 | 80
Driving Score | 0.976 | 0.982 | 0.976 | 0.918
Penalty Factor | 0.976 | 0.982 | 0.992 | 0.941
Score Route | 1 | 1 | 0.984 | 0.977
Layout crash | 0 | 0 | 0 | 0.047
Walker crash | 0.039 | 0 | 0 | 0.151
Vehicle crash | 0.018 | 0.039 | 0.083 | 0.167
Vehicle halt | 0 | 0 | 0.313 | 1.556
Lights met | 4.4 | 4.4 | 3.7 | 3.7
Lights passed | 4.38 | 4.36 | 3.7 | 3.7
Table 3: Ablation study on CARLA NoCrash benchmarking suite. T1: Town 1; T2:
Town 2.

## Conclusion

This paper discusses the challenges faced by rule-based and imitation learning
methods that rely on expert demonstrations. It highlights the limited ability
of these methods to recover from domain shifts and near-failure scenarios in
non-stationary environments. To mitigate this problem, this paper presents a
deep Graph-based Prediction and Planning Network (GP3Net) framework that
encodes spatio-temporal dynamic interactions of different traffic participants
using a graphical model to infer their future behaviors and combines the
advantage of exploration in reinforcement learning. GP3Net improves AV
performance compared to the previous graph-encoded imitation learning-based
policy in CARLA benchmarking driving scenarios for urban and highway
environments. The qualitative study also clearly indicates the importance of
the proposed prediction module for safe motion planning in dynamic driving
scenarios. Future work will look at other methods to model interactions with
the dynamics of traffic participants involved and how they can influence the
planning capabilities of AV.

## References

  * Aksjonov and Kyrki (2021) Aksjonov, A.; and Kyrki, V. 2021.  Rule-Based Decision-Making System for Autonomous Vehicles at Intersections with Mixed Traffic Environment.  In _2021 IEEE International Intelligent Transportation Systems Conference (ITSC)_ , 660–666.
  * Alizadeh et al. (2019) Alizadeh, A.; Moghadam, M.; Bicer, Y.; Ure, N. K.; Yavas, U.; and Kurtulus, C. 2019.  Automated Lane Change Decision Making using Deep Reinforcement Learning in Dynamic and Uncertain Highway Environment.  In _2019 IEEE Intelligent Transportation Systems Conference (ITSC)_ , 1399–1404.
  * Bansal, Krizhevsky, and Ogale (2018) Bansal, M.; Krizhevsky, A.; and Ogale, A. 2018.  Chauffeurnet: Learning to drive by imitating the best and synthesizing the worst.  _arXiv preprint arXiv:1812.03079_.
  * Bazzi et al. (2021) Bazzi, A.; Berthet, A. O.; Campolo, C.; Masini, B. M.; Molinaro, A.; and Zanella, A. 2021.  On the Design of Sidelink for Cellular V2X: A Literature Review and Outlook for Future.  _IEEE Access_ , 9: 97953–97980.
  * Cai et al. (2021) Cai, P.; Wang, H.; Sun, Y.; and Liu, M. 2021.  DiGNet: Learning Scalable Self-Driving Policies for Generic Traffic Scenarios with Graph Neural Networks.  In _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 8979–8984.
  * Cai et al. (2020) Cai, P.; Wang, S.; Sun, Y.; and Liu, M. 2020.  Probabilistic End-to-End Vehicle Navigation in Complex Dynamic Environments With Multimodal Sensor Fusion.  _IEEE Robotics and Automation Letters_ , 5(3): 4218–4224.
  * Chen, Li, and Tomizuka (2021) Chen, J.; Li, S. E.; and Tomizuka, M. 2021.  Interpretable end-to-end urban autonomous driving with latent deep reinforcement learning.  _IEEE Transactions on Intelligent Transportation Systems_ , 23(6): 5068–5078.
  * Chen, Yuan, and Tomizuka (2019) Chen, J.; Yuan, B.; and Tomizuka, M. 2019.  Model-free deep reinforcement learning for urban autonomous driving.  In _2019 IEEE Intelligent Transportation Systems Conference (ITSC)_ , 2765–2771. IEEE.
  * Chen et al. (2020) Chen, L.; Hu, X.; Tang, B.; and Cheng, Y. 2020.  Conditional DQN-based motion planning with fuzzy logic for autonomous driving.  _IEEE Transactions on Intelligent Transportation Systems_ , 23(4): 2966–2977.
  * Choi (2010) Choi, E.-H. 2010.  Crash factors in intersection-related crashes: An on-scene perspective.  Technical report.
  * Chou, Maturana, and Scherer (2017) Chou, P.-W.; Maturana, D.; and Scherer, S. 2017.  Improving stochastic policy gradients in continuous control with deep reinforcement learning using the beta distribution.  In _International conference on machine learning_ , 834–843. PMLR.
  * Codevilla et al. (2018) Codevilla, F.; Müller, M.; López, A.; Koltun, V.; and Dosovitskiy, A. 2018.  End-to-End Driving Via Conditional Imitation Learning.  In _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , 4693–4700.
  * Codevilla et al. (2019) Codevilla, F.; Santana, E.; López, A. M.; and Gaidon, A. 2019.  Exploring the limitations of behavior cloning for autonomous driving.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 9329–9338.
  * Cui et al. (2022) Cui, J.; Qiu, H.; Chen, D.; Stone, P.; and Zhu, Y. 2022.  Coopernaut: End-to-End Driving With Cooperative Perception for Networked Vehicles.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 17252–17262.
  * Dosovitskiy et al. (2017) Dosovitskiy, A.; Ros, G.; Codevilla, F.; Lopez, A.; and Koltun, V. 2017.  CARLA: An open urban driving simulator.  In _Conference on robot learning_ , 1–16. PMLR.
  * Haarnoja et al. (2018) Haarnoja, T.; Zhou, A.; Abbeel, P.; and Levine, S. 2018.  Soft Actor-Critic: Off-Policy Maximum Entropy Deep Reinforcement Learning with a Stochastic Actor.  In Dy, J.; and Krause, A., eds., _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , 1861–1870. PMLR.
  * Huang et al. (2019) Huang, X.; Hong, S.; Hofmann, A.; and Williams, B. C. 2019.  Online risk-bounded motion planning for autonomous vehicles in dynamic environments.  In _Proceedings of the International Conference on Automated Planning and Scheduling_ , volume 29, 214–222.
  * Ivanovic and Pavone (2019) Ivanovic, B.; and Pavone, M. 2019.  The trajectron: Probabilistic multi-agent trajectory modeling with dynamic spatiotemporal graphs.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2375–2384.
  * Li and Zhu (2021) Li, M.; and Zhu, Z. 2021.  Spatial-temporal fusion graph neural networks for traffic flow forecasting.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 35, 4189–4196.
  * Ma et al. (2019) Ma, Y.; Zhu, X.; Zhang, S.; Yang, R.; Wang, W.; and Manocha, D. 2019.  TrafficPredict: Trajectory prediction for heterogeneous traffic-agents.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 33, 6120–6127.
  * Naveed, Qiao, and Dolan (2021) Naveed, K. B.; Qiao, Z.; and Dolan, J. M. 2021.  Trajectory Planning for Autonomous Vehicles Using Hierarchical Reinforcement Learning.  In _2021 IEEE International Intelligent Transportation Systems Conference (ITSC)_ , 601–606.
  * Petrović, Mijailović, and Pešić (2020) Petrović, D.; Mijailović, R.; and Pešić, D. 2020.  Traffic accidents with autonomous vehicles: type of collisions, manoeuvres and errors of conventional vehicles’ drivers.  _Transportation research procedia_ , 45: 161–168.
  * Rhinehart, McAllister, and Levine (2020) Rhinehart, N.; McAllister, R.; and Levine, S. 2020.  Deep Imitative Models for Flexible Inference, Planning, and Control.  In _International Conference on Learning Representations_.
  * Schulman et al. (2017) Schulman, J.; Wolski, F.; Dhariwal, P.; Radford, A.; and Klimov, O. 2017.  Proximal policy optimization algorithms.  _arXiv preprint arXiv:1707.06347_.
  * Shi et al. (2022) Shi, L.; Wang, L.; Long, C.; Zhou, S.; Zheng, F.; Zheng, N.; and Hua, G. 2022.  Social interpretable tree for pedestrian trajectory prediction.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 36, 2235–2243.
  * Tang et al. (2022) Tang, X.; Huang, B.; Liu, T.; and Lin, X. 2022.  Highway Decision-Making and Motion Planning for Autonomous Driving via Soft Actor-Critic.  _IEEE Transactions on Vehicular Technology_ , 71(5): 4706–4717.
  * Teng et al. (2023) Teng, S.; Chen, L.; Ai, Y.; Zhou, Y.; Xuanyuan, Z.; and Hu, X. 2023.  Hierarchical Interpretable Imitation Learning for End-to-End Autonomous Driving.  _IEEE Transactions on Intelligent Vehicles_ , 8(1): 673–683.
  * Xu, Yang, and Du (2020) Xu, Y.; Yang, J.; and Du, S. 2020.  CF-LSTM: Cascaded feature-based long short-term networks for predicting pedestrian trajectory.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, 12541–12548.
  * Ye et al. (2020) Ye, F.; Cheng, X.; Wang, P.; Chan, C.-Y.; and Zhang, J. 2020.  Automated lane change strategy using proximal policy optimization-based deep reinforcement learning.  In _2020 IEEE Intelligent Vehicles Symposium (IV)_ , 1746–1752. IEEE.

## Supplementary Materials for GP3Net

The supplementary materials consists of 1) More qualitative analysis 2) More
quantitative results for different towns, 3) Architectural details of the
neural network models used in GP3Net, 4) Hyper-parameter settings for GP3Net
and 5) video recording of trained GP3Net based AV in different towns in
different weathers.

### Driving Videos by GP3Net

There are six videos included in the videos folder. The videos show many
scenarios in different towns (small cities, big cities, highways, and mixed
towns) and weather conditions (rainy and clear). There are urban
intersections, junctions, highway lane merging, and diversions. The GP3Net
performed better in all scenarios.

### Qualitative Analysis

Figure 4: Plot of the layout collisions with and without trajectory prediction
module. Figure 5: Plot of the percentage of cases going outside lane with and
without trajectory prediction module.

This section shows some plots about the importance of the prediction module in
motion planning for AVs. A PPO algorithm is trained without the prediction
module to check its performance. The following plots show the performance in
different metrics defined by CARLA with and without the trajectory prediction
module.

Figure 6: Plot of the score composed with and without trajectory prediction
module.

In Fig.4 and Fig. 5, it can be seen that up to 2M steps, the GP3Net model
without trajectory prediction performed better. However, after 2M training
steps, the GP3Net with trajectory prediction module started performing well.

Performance Metric |  | Town01
---
Urban
| Town02
---
Urban
| Town03
---
Urban
| Town04
---
Mixed
| Town05
---
Urban
| Town06
---
Highway
End Reach (%) | 100 | 95.43 | 96.25 | 93.75 | 97.86 | 95.0
Success Rate (%) | 92.5 | 93.14 | 82.5 | 75.0 | 82.8 | 82.5
Driving Score | 0.96 | 0.91 | 0.87 | 0.85 | 0.91 | 0.91
Penalty Factor | 0.9588 | 0.9840 | 0.9062 | 0.8721 | 0.9223 | 0.914
Score Route | 1.0 | 0.9994 | 0.9866 | 0.9794 | 0.9916 | 0.9832
Layout Crash | 0 | 0 | 0 | 0.0115 | 0 | 0
Walker Crash | 0.0922 | 0 | 0.2141 | 0.0419 | 0.015 | 0.0104
Vehicle Crash | 0 | 0.0508 | 0.0320 | 0.0579 | 0.1291 | 0.0865
Vehicle Halt | 0 | 0.2509 | 0.0394 | 0.0439 | 0.0572 | 0.0296
Lights met | 4.7 | 3.76 | 7.638 | 8.65 | 8.2357 | 6.025
Lights passed | 4.675 | 3.7371 | 7.5375 | 8.5375 | 8.2214 | 5.975
Table 4: Detailed performance of GP3Net on Leaderboard Benchmarking scenarios
in different towns

In Fig.6, the advantage of the trajectory prediction module is visible. After
2M training steps, there is improvement in the metric score composed. This
metric combines infractions and penalty factors and provides the final score
for driving.

### Detailed Quantitative Results

The detailed quantitative results for the CARLA Leaderboard benchmarking
scenarios are given in Table.4.

### Trajectory Prediction Architecture Details

The trajectory prediction module uses a deep generative trajectory prediction
architecture which is has Conditional Variational Auto Encoder (CVAE) to
generate a distribution of potential trajectories. CVAE, a latent variable
model, rely on Gaussian Mixture Models (GMM) to produce multimodal
distribution outputs. The encoder and decoder of the CVAE consists of Long
Short Term Memory (LSTM) networks for modeling temporal data such as
trajectories and system evolution. After representing the nodes and edges as a
vector, they are encoded using a Node Encoder and Edge Encoder respectively,
whose hidden dimensions are shown in Table 5. The encoded nodes and edges are
sent to an attention module and low-pass filter $M$. The probability density
function $p(z|\textbf{x})$ is captured by a Fully Connected Layer. Here, $z$
is the latent variable and x is the input to the model. After sampling $z$
from $p(z|\textbf{x})$, it is passed through the decoder LSTM network which
outputs the parameters for the GMM, from which multimodal trajectories are
sampled. Table 5 shows the Trajectory Prediction Architecture details.

Network Architecture Components | Value
---|---
Batch Size | 16
Learning Rate | 0.001
Minimum Learning Rate | 0.00001
Learning Decay Rate | 0.9999
Edge Addition distance (vehicles) | 45 meters
LSTM Encoder Edge | 8
| LSTM Node History Encoder
---
Hidden Dimensions
32
| LSTM Node Future Encoder
---
Hidden Dimensions
32
FC Layer $p(z|\textbf{x})$ Dimensions | 16
LSTM Decoder Hidden Dimensions | 128
GMM Components | 16
Past Data Length | 8 steps (3.2 sec)
Future Prediction Horizon | 7 steps (2.8 sec)
Table 5: Trajectory Prediction Module Architecture and Hyperparameters Details

### State Encoder Architecture Details

The State Encoder consists of Convolutional Neural Network (CNN) for encoding
the BEV (past and generated future) masks and a Multi-Layered Perceptron to
encode the odometry values of AV. After concatenating the outputs from these
two networks, it is further processed using a Fully-Connected Layer giving the
final state encoding. The details of the State Encoder architecture is
displayed in Table 6.

| State Encoder
---
Layers
| IN
---
Channels/
Dims
| OUT
---
Channels/
Dims
| Filter
---
Size
Conv 1 (BEV) | 21 | 8 | 5
Conv 2 (BEV) | 8 | 16 | 5
Conv 3 (BEV) | 16 | 32 | 5
Conv 4 (BEV) | 32 | 64 | 3
Conv 5 (BEV) | 64 | 128 | 3
Conv 6 (BEV) | 128 | 256 | 3
| Fully Connected
---
(Odometry)
256 | 256 | -
| Fully Connected
---
(Odometry)
256 | 256 | -
Fully Connected | 256 | 512 | -
Fully Connected | 512 | 256 | -
Table 6: State Encoder Architecture

### PPO Policy Network Architecture Details and Hyperparameters

The policy network $\pi$ like in many RL works is a Multi-Layered Perceptron
(MLP). This policy network $\pi$ gives the values of $\alpha$ and $\beta$
parameters of Beta distribution of throttle and steering, resulting in two
$\alpha$s and two $\beta$s. The details for policy network architecture and
hyperparameters used for PPO training is shown in Table 7.

Policy Network Layers | IN dims | OUT dims
---|---|---
Fully Connected | 256 | 256
Fully Connected ($\alpha$) | 256 | 2
Fully Connected ($\beta$) | 256 | 2
PPO Hyperparameters | Values
Clip ratio $\epsilon$ | 0.2
Epochs | 30
Gamma $\gamma$ | 0.99
Target KL | 0.01
GAE $\lambda$ | 0.97
Learning Rate | 0.00003
Table 7: Policy Network Architecture and Hyperparameter Details

### Data collection settings for trajectory prediction module

The training data was collected using the CARLA simulator’s Autopilot/Roaming
Agent in Town 1. For several episodes, this autopilot vehicle navigates from
source to destination, recording the trajectories of all vehicles at each time
step $t$. The trajectories are represented by a sequence of local rectangular
coordinates $\textbf{x}_{t}=(x_{t},y_{t})$ of the agents at time step $t$, as
observed by the roaming agent. The dataset contains the timestamp of the
recorded trajectory data, the index of the agent, and the trajectory
coordinates. The data for vehicles and pedestrians are collected every 4-time
steps (0.4 seconds) and stored in separate datasets, as we train two models
for the former and the latter.