# Finite Sums and Products involving Special Functions

Robert Reynolds  Department of Mathematics and Statistics, York University,
Toronto, ON, Canada, M3J1P3<EMAIL_ADDRESS>
###### Abstract.

Various product and sum relationships are established using special functions,
specifically involving Special functions. These relationships are derived from
formulas inspired by the finite sum that incorporates the Hurwitz-Lerch zeta
function.

###### Key words and phrases:

Finite sum, finite product, Special function, contour integral, Nielsen, Gauss

###### 2020 Mathematics Subject Classification:

Primary 30E20, 33-01, 33-03, 33-04

##  1\. Introduction

Finite sums and products involving Special functions refers to mathematical
expressions that involve the summation or multiplication of a finite number of
terms, where the terms themselves are represented by special functions. These
special functions are typically mathematical functions that have specific
properties or applications in various branches of mathematics and physics. The
finite sums and products involving special functions can arise in the study of
series, sequences, combinatorics, and other areas where these specialized
mathematical functions are utilized.

The finite sum of the Hurwitz-Lerch zeta function finds applications in
various areas of mathematics and physics. A few examples are; Number theory
section (27.17) in [1], where the finite sum of the Hurwitz-Lerch zeta
function is closely related to the theory of special values of $L$-functions
and Dirichlet series. It has been used in the study of modular forms, elliptic
curves, and other topics in number theory. Another example is in Analytic
number theory [2], where the finite sum of the Hurwitz-Lerch zeta function
appears in the analysis of certain arithmetic functions, such as the divisor
function and the Riemann zeta function. It can be used to study the
distribution of prime numbers and investigate properties of arithmetic
sequences. In Statistical physics, see section (24.18) in [1], the finite sum
of the Hurwitz-Lerch zeta function arises in the study of thermodynamic
properties of physical systems, particularly in the context of lattice models
and quantum field theory. It is used to analyze partition functions,
correlation functions, and critical phenomena. In Quantum mechanics, see
section (36.14) in [1], the finite sum of the Hurwitz-Lerch zeta function is
encountered in the calculation of energy levels and wave functions of quantum
mechanical systems with certain potentials or boundary conditions. It provides
a mathematical tool for solving differential equations in quantum mechanics.
In the area of Special functions and mathematical physics, [3, 1, 4, 5], the
finite sum of the Hurwitz-Lerch zeta function is part of a broader class of
special functions used in mathematical physics. It appears in the context of
integral transforms, differential equations, and solutions to boundary value
problems. These applications demonstrate the versatility and utility of the
finite sum of the Lerch function in various branches of mathematics and
physics.

The book authored by Prudnikov et al. [6] provides a comprehensive compilation
of both indefinite and definite integrals, as well as finite and infinite sums
and products involving elementary and special functions. The book offers an
extensive inventory of these mathematical expressions, catering to a wide
range of applications and research areas. In this current study, we aim to
build upon previous research that focused on the finite sum of special
functions. To achieve this, we employ the contour integral method [7],
specifically applied to equation (4.4.7.12) mentioned in [6]. Consequently, we
obtain the contour integral representation as a result given by;

(1.1) $\frac{1}{2\pi
i}\int_{C}\sum_{p=0}^{n}\left(-2^{-p}a^{w}w^{-k-1}\tan\left(2^{-p-1}(m+w)\right)\sec\left(2^{-p}(m+w)\right)\right)dw\\\
\\\ =\frac{1}{2\pi
i}\int_{C}\left(2^{-n}a^{w}w^{-k-1}\left(\csc\left(2^{-n}(m+w)\right)-2^{n+1}\csc(2(m+w))\right)\right)dw$

where $a,m,k\in\mathbb{C},Re(m+w)>0,n\in\mathbb{Z^{+}}$. Using equation (1.1)
the main Theorem to be derived and evaluated is given by

(1.2)
$\sum_{p=0}^{n}2^{-p}e^{im2^{-p}}\left(\left(i2^{-p}\right)^{k}e^{im2^{-p}}\Phi\left(-e^{i2^{1-p}m},-k,1-i2^{p-1}\log(a)\right)\right.\\\
\left.-\left(i2^{-p-1}\right)^{k}\Phi\left(-e^{i2^{-p}m},-k,1-i2^{p}\log(a)\right)\right)\\\
\\\
=i\left(i2^{-n}\right)^{k+1}e^{im2^{-n}}\Phi\left(e^{i2^{1-n}m},-k,\frac{1}{2}\left(1-i2^{n}\log(a)\right)\right)\\\
+i^{k}2^{k+1}e^{2im}\Phi\left(e^{4im},-k,\frac{1}{2}-\frac{1}{4}i\log(a)\right)$

The expression provided involves the variables $k$, $a$, $m$, which can be any
complex numbers, and $n$, which represents any positive integer. This
expression is subsequently utilized to obtain specific instances using
trigonometric functions. The derivations employed in this process are based on
the approach presented in our previous work referenced as [7]. This approach
incorporates a variant of the generalized Cauchy’s integral formula, as given
by;

(1.3) $\frac{y^{k}}{\Gamma(k+1)}=\frac{1}{2\pi
i}\int_{C}\frac{e^{wy}}{w^{k+1}}dw,$

In this context, we consider variables $y$ and $w$ belonging to the set of
complex numbers ($\mathbb{C}$), while $C$ represents a general open contour in
the complex plane. The bilinear concomitant, as described in the reference
[7], vanishes at the endpoints of the contour. This approach involves
utilizing a specific form of equation (1.3), multiplying both sides by a
function, and subsequently summing the resulting finite terms on both sides.
As a result, a finite sum expressed in terms of a contour integral is
obtained. Additionally, by multiplying both sides of equation (1.3) by another
function and summing infinitely on both sides, the contour integrals in both
equations become identical.

##  2\. The Hurwitz-Lerch Zeta Function

In our analysis, we employ equation (1.11.3) from the reference [8], where
$\Phi(z,s,v)$ denotes the Hurwitz-Lerch zeta function. This function serves as
a generalization of both the Hurwitz zeta function $\zeta(s,v)$ and the
Polylogarithm function $\text{Li}_{n}(z)$. The Lerch function can be expressed
using a series representation, which is given by.

(2.1) $\Phi(z,s,v)=\sum_{n=0}^{\infty}(v+n)^{-s}z^{n}$

where $|z|<1,v\neq 0,-1,-2,-3,..,$ and is continued analytically by its
integral representation given by

(2.2)
$\Phi(z,s,v)=\frac{1}{\Gamma(s)}\int_{0}^{\infty}\frac{t^{s-1}e^{-vt}}{1-ze^{-t}}dt=\frac{1}{\Gamma(s)}\int_{0}^{\infty}\frac{t^{s-1}e^{-(v-1)t}}{e^{t}-z}dt$

where $Re(v)>0$, and either $|z|\leq 1,z\neq 1,Re(s)>0$, or $z=1,Re(s)>1$.

##  3\. Contour Integral Representation for the Finite Sum of the Hurwitz-
Lerch Zeta Functions

In this section we derive the contour integral representations of the left-
hand side and right-hand side of equation (1.1) in terms of the Hurwtiz-Lerch
zeta and trigonometric functions.

###  3.1. Derivation of the left-hand side first contour integral

We employ the methodology described in reference [7]. By utilizing equation
(1.3), our initial step involves substituting $\log(a)+i2^{-p}(y+1)$ and
subsequently multiplying both sides of the equation by
$-i2^{1-p}(-1)^{y}e^{im2^{-p}(y+1)}$. We then proceed to perform finite and
infinite summations over the ranges $p\in[0,n]$ and $y\in[0,\infty)$
respectively. Finally, we simplify the resulting expression in terms of the
Hurwitz-Lerch Zeta function, yielding

(3.1)
$-\sum_{p=0}^{n}\frac{i2^{k-p+1}\left(i2^{-p-1}\right)^{k}e^{im2^{-p}}\Phi\left(-e^{i2^{-p}m},-k,1-i2^{p}\log(a)\right)}{\Gamma(k+1)}\\\
\ =-\frac{1}{2\pi
i}\sum_{y=0}^{\infty}\sum_{p=0}^{n}\int_{C}(-1)^{y}a^{w}w^{-k-1}e^{i2^{-p}(y+1)(m+w)}dw\\\
\\\ =-\frac{1}{2\pi
i}\int_{C}\sum_{p=0}^{n}\sum_{y=0}^{\infty}(-1)^{y}a^{w}w^{-k-1}e^{i2^{-p}(y+1)(m+w)}dw\\\
\\\ =\frac{1}{2\pi
i}\int_{C}\sum_{p=0}^{n}\left(2^{-p}a^{w}w^{-k-1}\tan\left(2^{-p-1}(m+w)\right)-i2^{-p}a^{w}w^{-k-1}\right)dw$

from equation (1.232.1) in [9] where $Re(w+m)>0$ and $Im\left(m+w\right)>0$ in
order for the sums to converge. We apply Tonelli’s theorem for multiple sums,
see page 177 in [10] as the summands are of bounded measure over the space
$\mathbb{C}\times[0,n]\times[0,\infty)$.

###  3.2. Derivation of the additional contour integral

Using equation (1.3) we replace $y$ with $\log(a)$ and multiply both sides by
$-i2^{-p}$ take the finite sum over $p\in[0,n]$ and simplify to get;

(3.2) $-\frac{1}{2\pi
i}\sum_{p=0}^{n}\frac{i2^{-p}\log^{k}(a)}{\Gamma(k+1)}=-\frac{1}{2\pi
i}\int_{C}\sum_{p=0}^{n}i2^{-p}a^{w}w^{-k-1}dw$

###  3.3. Derivation of the left-hand side second contour integral

We employ the methodology described in reference [7]. By utilizing equation
(1.3), our initial step involves substituting $\log(a)+i2^{1-p}(y+1)$ and
subsequently multiplying both sides of the equation by
$-i2^{1-p}(-1)^{y}e^{im2^{1-p}(y+1)}$. We then proceed to perform finite and
infinite summations over the ranges $p\in[0,n]$ and $y\in[0,\infty)$
respectively. Finally, we simplify the resulting expression in terms of the
Hurwitz-Lerch Zeta function, yielding

(3.3)
$\sum_{p=0}^{n}\frac{i2^{k-p+1}\left(i2^{-p}\right)^{k}e^{im2^{1-p}}\Phi\left(-e^{i2^{1-p}m},-k,1-i2^{p-1}\log(a)\right)}{\Gamma(k+1)}\\\
\ =-\frac{1}{2\pi
i}\sum_{y=0}^{\infty}\sum_{p=0}^{n}\int_{C}(-1)^{y}a^{w}w^{-k-1}e^{i2^{1-p}(y+1)(m+w)}dw\\\
\\\ =-\frac{1}{2\pi
i}\int_{C}\sum_{p=0}^{n}\sum_{y=0}^{\infty}(-1)^{y}a^{w}w^{-k-1}e^{i2^{1-p}(y+1)(m+w)}dw\\\
\\\ =\frac{1}{2\pi
i}\int_{C}\sum_{p=0}^{n}\left(-2^{-p}a^{w}w^{-k-1}\tan\left(2^{-p}(m+w)\right)+i2^{-p}a^{w}w^{-k-1}\right)dw$

from equation (1.232.1) in [9] where $Re(w+m)>0$ and $Im\left(m+w\right)>0$ in
order for the sums to converge. We apply Tonelli’s theorem for multiple sums,
see page 177 in [10] as the summands are of bounded measure over the space
$\mathbb{C}\times[0,n]\times[0,\infty)$.

###  3.4. Derivation of the additional contour integral

Using equation (1.3) we replace $y$ with $\log(a)$ and multiply both sides by
$i2^{-p}$ take the finite sum over $p\in[0,n]$ and simplify to get;

(3.4) $\frac{1}{2\pi
i}\sum_{p=0}^{n}\frac{i2^{-p}\log^{k}(a)}{\Gamma(k+1)}=\frac{1}{2\pi
i}\int_{C}\sum_{p=0}^{n}i2^{-p}a^{w}w^{-k-1}dw$

###  3.5. Derivation of the right-hand side first contour integral

We use the method in [7]. Using equation (1.3) we first replace
$\log(a)+2i(2y+1)$ and multiply both sides by $4ie^{2im(2y+1)}$ then take the
finite and infinite sum over $y\in[0,\infty)$ and simplify in terms of the
Hurwitz-Lerch Zeta function to get

(3.5)
$\frac{(4i)^{k+1}e^{2im}\Phi\left(e^{4im},-k,\frac{1}{2}-\frac{1}{4}i\log(a)\right)}{\Gamma(k+1)}\\\
\ =\frac{1}{2\pi i}\sum_{y=0}^{\infty}\int_{C}4ia^{w}e^{2i(2y+1)(m+w)}dw\\\
\\\ =-\frac{1}{2\pi i}\int_{C}\sum_{y=0}^{\infty}4ia^{w}e^{2i(2y+1)(m+w)}dw\\\
\\\ =-\frac{1}{2\pi i}\int_{C}2a^{w}w^{-k-1}\csc(2(m+w))dw$

from equation (1.232.3) in [9] where $Re(w+m)>0$ and $Im\left(m+w\right)>0$ in
order for the sums to converge. We apply Tonelli’s theorem for multiple sums,
see page 177 in [10] as the summands are of bounded measure over the space
$\mathbb{C}\times[0,\infty)$.

###  3.6. Derivation of the right-hand side second contour integral

We use the method in [7]. Using equation (1.3) we first replace
$\log(a)+i2^{-n}(2y+1)$ and multiply both sides by
$-i2^{1-n}e^{im2^{-n}(2y+1)}$ then take the finite and infinite sum over
$y\in[0,\infty)$ and simplify in terms of the Hurwitz-Lerch Zeta function to
get

(3.6)
$-\frac{i2^{k-n+1}\left(i2^{-n}\right)^{k}e^{im2^{-n}}\Phi\left(e^{i2^{1-n}m},-k,2^{n-1}\left(2^{-n}-i\log(a)\right)\right)}{\Gamma(k+1)}\\\
\ =-\frac{1}{2\pi
i}\sum_{y=0}^{\infty}\int_{C}i2^{1-n}a^{w}e^{i2^{-n}(2y+1)(m+w)}dw\\\ \\\
=-\frac{1}{2\pi
i}\int_{C}\sum_{y=0}^{\infty}i2^{1-n}a^{w}e^{i2^{-n}(2y+1)(m+w)}dw\\\ \\\
=-\frac{1}{2\pi i}\int_{C}2^{-n}a^{w}w^{-k-1}\csc\left(2^{-n}(m+w)\right)dw$

from equation (1.232.3) in [9] where $Re(w+m)>0$ and $Im\left(m+w\right)>0$ in
order for the sums to converge. We apply Tonelli’s theorem for multiple sums,
see page 177 in [10] as the summands are of bounded measure over the space
$\mathbb{C}\times[0,\infty)$.

##  4\. Hurwitz-Lerch zeta function identity

In this section we will derive the finite sum of Hurwitz-Lerch Zeta functions
in terms of the Hurwitz-Lerch Zeta function along with its composite
functions.

######  Theorem 4.1.

For all $k,a,m\in\mathbb{C},n\in\mathbb{Z^{+}}$ then,

(4.1)
$\sum_{p=0}^{n}2^{-p}e^{im2^{-p}}\left(\left(i2^{-p}\right)^{k}e^{im2^{-p}}\Phi\left(-e^{i2^{1-p}m},-k,1-i2^{p-1}\log(a)\right)\right.\\\
\left.-\left(i2^{-p-1}\right)^{k}\Phi\left(-e^{i2^{-p}m},-k,1-i2^{p}\log(a)\right)\right)\\\
=i\left(i2^{-n}\right)^{k+1}e^{im2^{-n}}\Phi\left(e^{i2^{1-n}m},-k,\frac{1}{2}\left(1-i2^{n}\log(a)\right)\right)\\\
+i^{k}2^{k+1}e^{2im}\Phi\left(e^{4im},-k,\frac{1}{2}-\frac{1}{4}i\log(a)\right)$

###### Proof.

With respect to equation (1.1) and observing the addition of the right-hand
sides of relations (3.1), (3.2), (3.3) and (3.4), and the addition of
relations (3.5) and (3.6) are identical; hence, the left-hand sides of the
same are identical too. Simplifying with the Gamma function yields the desired
conclusion. ∎

######  Example 4.2.

The degenerate case.

(4.2)
$\sum_{p=0}^{n}2^{-p-1}\tan\left(m2^{-p-1}\right)\sec\left(m2^{-p}\right)=\csc(2m)-2^{-n-1}\csc\left(m2^{-n}\right)$

###### Proof.

Use equation (4.1) and set $k=0$ and simplify using entry (2) in Table below
(64:12:7) in [4]. Note equation (4.4.7.12) in [6] is in error. ∎

######  Example 4.3.

The product of the ratio of cosine functions in terms of the ratio of tangent
functions.

(4.3)
$\prod_{p=0}^{n}\frac{\cos\left(2^{-p}m\right)}{\cos\left(2^{-p}r\right)}\left(\frac{\cos\left(2^{-1-p}r\right)}{\cos\left(2^{-1-p}m\right)}\right)^{2}=\frac{\tan\left(2^{-1-n}m\right)\tan(r)}{\tan(m)\tan\left(2^{-1-n}r\right)}$

###### Proof.

Use equation (4.1) and form a second equation by replacing $m\to r$ take the
difference of both these equations then set $k=-1,a=1$ and simplify using
entry (3) of Section (64:12) in [4]. ∎

######  Example 4.4.

A functional equation for the Hurwitz-Lerch zeta function.

(4.4)
$\Phi(z,s,a)=8^{-s}\left(4^{s}z\Phi\left(-z^{2},s,\frac{a+1}{2}\right)+4^{s}\Phi\left(z^{2},s,\frac{a}{2}\right)\right.\\\
\left.-2z^{3}\left(2^{s}\Phi\left(-z^{4},s,\frac{a+3}{4}\right)-2\Phi\left(z^{8},s,\frac{a+3}{8}\right)\right)\right)$

###### Proof.

Use equation (4.1) and set $n=1,m=2\log(-z)/i,k=-s,a=e^{(a-1)i/2}$ and
simplify.. ∎

######  Example 4.5.

The product of the ratio of cosine functions in terms of the ratio of tangent
functions.

(4.5)
$\prod_{p=0}^{n}\frac{\cos^{3}\left(2^{-1-p}x\right)}{\cos^{2}\left(2^{-2-p}x\right)\cos\left(2^{-p}x\right)}\\\
=\frac{\tan(x)\tan\left(2^{-2-n}x\right)}{\tan\left(\frac{x}{2}\right)\tan\left(2^{-1-n}x\right)}$

###### Proof.

Use equation (4.1) and set $k=1,a=1,m=x$ and simplify using the method in
section (8.1) in [11]. ∎

######  Example 4.6.

The product of the exponential function and ratio of cosine functions in terms
of the ratio of tangent functions.

(4.6)
$\prod_{p=0}^{n}\cos^{2}\left(2^{-p-2}x\right)\cos\left(2^{-p}x\right)\sec^{3}\left(2^{-p-1}x\right)\\\
\exp\left(-2^{1-p}\left(\cos\left(2^{-p-1}x\right)+\cos\left(3\
2^{-p-1}x\right)-3\cos\left(2^{-p}x\right)+1\right)\csc\left(2^{1-p}x\right)\right)\\\
=\tan\left(\frac{x}{2}\right)\cot(x)\tan\left(2^{-n-1}x\right)\cot\left(2^{-n-2}x\right)\\\
\exp\left(2^{-n}\left(\csc\left(2^{-n}x\right)-\csc\left(2^{-n-1}x\right)\right)+\tan\left(\frac{x}{2}\right)-\tan(x)+\cot\left(\frac{x}{2}\right)-\cot(x)\right)$

###### Proof.

Use equation (4.1) and set $k=1,a=e,m=x$ and simplify using the method in
section (8.1) in [11]. ∎

######  Example 4.7.

Sum of the log-gamma function.

(4.7)
$\sum_{p=0}^{n}2^{-p}\left(2\text{log$\Gamma$}\left(-i2^{p-2}\log(a)\right)-2\text{log$\Gamma$}\left(-i2^{p-1}\log(a)\right)\right.\\\
\left.-2\text{log$\Gamma$}\left(\frac{1}{4}\left(-i2^{p}\log(a)-2\right)\right)+2\text{log$\Gamma$}\left(\frac{1}{2}\left(-i2^{p}\log(a)-1\right)\right)+\log\left(\frac{2\left(2^{p}\log(a)-i\right)^{2}}{\left(2^{p}\log(a)-2i\right)^{2}}\right)\right)\\\
=2^{-n-1}\left(-2^{n}\left(8\text{log$\Gamma$}\left(-\frac{1}{4}i\log(a)-\frac{1}{2}\right)+\log(a)\left(2i\log\left(i2^{-n}\right)+\pi-2i\log(2)\right)\right.\right.\\\
\left.\left.+8\log(-2-i\log(a))-4\log(32\pi)\right)+4\text{log$\Gamma$}\left(\frac{1}{2}\left(-i2^{n}\log(a)-1\right)\right)\right.\\\
\left.+4\log\left(-1-i2^{n}\log(a)\right)-2\log(\pi)-6\log(2)\right)$

###### Proof.

Use equation (4.1) and set $m=0$, then simplify using equation (25.14.2) in
[1]. Next take the first partial derivative with respect to $k$ and set $k=0$
and simplify using equation (25.11.18) in [1]. ∎

######  Example 4.8.

Sum of the log-gamma function alternate form.

(4.8)
$\sum_{p=0}^{n}2^{-p}\left(2\text{log$\Gamma$}\left(2^{p-2}a\right)-2\text{log$\Gamma$}\left(2^{p-1}a\right)-2\text{log$\Gamma$}\left(\frac{1}{4}\left(2^{p}a-2\right)\right)\right.\\\
\left.+2\text{log$\Gamma$}\left(\frac{1}{2}\left(2^{p}a-1\right)\right)+\log\left(\frac{2\left(a2^{p}-1\right)^{2}}{\left(a2^{p}-2\right)^{2}}\right)\right)\\\
=2^{-n}\left(2\text{log$\Gamma$}\left(\frac{1}{2}\left(2^{n}a-1\right)\right)+2\log\left(\frac{a2^{n}-1}{2\sqrt{2\pi}}\right)\right)-4\text{log$\Gamma$}\left(\frac{a-2}{4}\right)\\\
+a\log\left(2^{-n-1}\right)+4\log\left(\frac{4\sqrt{2\pi}}{a-2}\right)$

###### Proof.

Use equation (4.1) and set $m=0,a=e^{ai}$ then simplify using equation
(25.14.2) in [1] Next take the first partial derivative with respect to $k$
and set $k=0$ and simplify using equation (25.11.18) in [1]. ∎

######  Example 4.9.

Sum involving the digamma function.

(4.9)
$\sum_{p=0}^{n}\left(\frac{4}{a2^{p}\left(a2^{p}-3\right)+2}-\psi^{(0)}\left(2^{p-2}a\right)+2\psi^{(0)}\left(2^{p-1}a\right)+\psi^{(0)}\left(\frac{1}{4}\left(2^{p}a-2\right)\right)\right.\\\
\left.-2\psi^{(0)}\left(\frac{1}{2}\left(2^{p}a-1\right)\right)\right)\\\
=-2\left(\frac{2}{a2^{n}-1}+\psi^{(0)}\left(\frac{1}{2}\left(2^{n}a-1\right)\right)-\frac{4}{a-2}-\psi^{(0)}\left(\frac{a-2}{4}\right)+\log\left(2^{-n-1}\right)\right)$

###### Proof.

Use equation (4.8) and take the first partial derivative with respect to $a$
and simplify using equation (5.2.2) in [1]. A similar form is given in
equation (2.1) in [12]. ∎

######  Example 4.10.

Log-gamma transformation.

(4.10)
$\log\left(\Gamma\left(\frac{a}{4}\right)\right)+\log\left(\Gamma\left(\frac{a-2}{4}\right)\sqrt{\frac{\Gamma\left(\frac{a-1}{2}\right)}{\Gamma\left(\frac{a}{2}\right)\Gamma(a)}}\right)=2\log\left(\frac{\pi^{3/8}2^{2-\frac{a}{2}}\sqrt[4]{a+\frac{1}{a-1}-3}}{a-2}\right)$

###### Proof.

Use equation (4.9) and set $n=1$ and simplify. Similar form is given by
equation (8.1) in [12]. ∎

######  Example 4.11.

Extended Nielsen product form.

(4.11)
$\prod_{p=1}^{n}\left(\frac{2^{-2^{p-1}x}\Gamma\left(\frac{1}{2}\left(2^{p}x+1\right)\right)}{\Gamma\left(\frac{1}{4}\left(2^{p}x+2\right)\right)^{2}}\right)^{2^{-p}}=\frac{2^{-\frac{nx}{2}-2^{-n}}\left(2^{n}x-1\right)^{2^{-n}}e^{2^{-n}\text{log$\Gamma$}\left(\frac{1}{2}\left(2^{n}x-1\right)\right)}}{\Gamma\left(\frac{x+1}{2}\right)}$

###### Proof.

Use equation (4.8) and take the exponential of both sides and simplify. Note
the exponential of the $\log\Gamma(x)$ on the right-hand side can be
simplified for real $x$. Similar forms are given by equations (6.9) and
(38.12) in [12]. ∎

######  Example 4.12.

The infinite limiting case.

(4.12)
$\prod_{p=1}^{\infty}\left(\frac{2^{-2^{p-1}x}\Gamma\left(\frac{1}{2}\left(2^{p}x+1\right)\right)}{\Gamma\left(\frac{1}{4}\left(2^{p}x+2\right)\right)^{2}}\right)^{2^{-p}}=\frac{(2e)^{-x/2}x^{x/2}}{\Gamma\left(\frac{x+1}{2}\right)}$

###### Proof.

Use equation (4.11) and take the limit of the right-hand side as $n\to\infty$
and simplify. This expression holds for $0<Re(x)<1,-1<Im(x)<1$ where $x$ is
small. Similar forms are in equation (90.1.10) in [13] and equation (8.323) in
[9]. ∎

######  Example 4.13.

Finite sum involving the generalized Stieltjes constant.

(4.13)
$\sum_{p=0}^{n}\left(\log\left(i2^{1-p}\right)\left(H_{2^{p-2}a}-H_{\frac{1}{4}\left(2^{p}a-2\right)}\right)+2\log\left(i2^{-p}\right)\left(H_{\frac{1}{2}\left(2^{p}a-1\right)}-H_{2^{p-1}a}\right)\right.\\\
\left.-\gamma_{1}\left(2^{p-2}a+1\right)+2\gamma_{1}\left(2^{p-1}a+1\right)-2\gamma_{1}\left(\frac{1}{2}\left(2^{p}a+1\right)\right)+\gamma_{1}\left(\frac{1}{4}\left(2^{p}a+2\right)\right)\right)\\\
=\frac{1}{4}\left(-8\gamma_{1}\left(\frac{1}{2}\left(2^{n}a+1\right)\right)+8\log\left(i2^{-n}\right)\psi^{(0)}\left(\frac{1}{2}\left(2^{n}a+1\right)\right)\right.\\\
\left.+8\gamma_{1}\left(\frac{a+2}{4}\right)+(-8\log(2)-4i\pi)\psi^{(0)}\left(\frac{a+2}{4}\right)+4\log^{2}\left(i2^{-n}\right)+(\pi-2i\log(2))^{2}\right)$

###### Proof.

Use equation (4.1) and set $m=0$ then simplify using equation (25.14.2) in
[1]. Next take the first partial derivative with respect to $k$ and take the
limit as $k\to-1$ and simplify using equation (25.6.12) in [1]. Similar is
given equation (1.1) in [14]. ∎

######  Example 4.14.

Finite sum involving the polylogarithm function.

(4.14)
$\sum_{p=0}^{n}2^{-p}\left(\left(2^{-p-1}\right)^{k}\text{Li}_{-k}\left(-e^{i2^{-p}m}\right)-\left(2^{-p}\right)^{k}\text{Li}_{-k}\left(-e^{i2^{1-p}m}\right)\right)\\\
=2^{k+1}e^{2im}\Phi\left(e^{4im},-k,\frac{1}{2}\right)-\left(2^{-n}\right)^{k+1}e^{im2^{-n}}\Phi\left(e^{i2^{1-n}m},-k,\frac{1}{2}\right)$

###### Proof.

Use equation (4.1) set $a=1$ and simplify using equation (25.14.3 ) in [1].
Similar form is given by equation (3.3.47) in [15]. ∎

######  Example 4.15.

Finite product of the exponential of trigonometric functions.

(4.15)
$\prod_{p=0}^{n}\exp\left(4^{1-p}\left(\sec^{2}\left(2^{-p-2}x\right)-3\sec^{2}\left(2^{-p-1}x\right)+2\sec^{2}\left(2^{-p}x\right)\right)\right)\\\
=\exp\left(2^{1-2n}\left(-\csc^{2}\left(2^{-n-2}x\right)+\csc^{2}\left(2^{-n-1}x\right)+\sec^{2}\left(2^{-n-2}x\right)\right.\right.\\\
\left.\left.-\sec^{2}\left(2^{-n-1}x\right)\right)+32\cot(x)\csc(x)-32\cot(2x)\csc(2x)\right)$

###### Proof.

Use equation (4.1) and set $k=2,a=1,m=x$ and simplify using the method in
section (8.1) in [11]. ∎

##  5\. Conclusion

This paper introduces a technique that enables the derivation of a finite sum
identity incorporating the Hurwitz-Lerch zeta function, along with intriguing
sums and products involving Special functions. The obtained results were
numerically validated using Mathematica by Wolfram, for parameters in the
integrals encompassing real, imaginary, and complex values.

## References

  * [1] Olver, F.W.J.; Lozier, D.W.; Boisvert, R.F.; Clark, C.W. (Eds.) _NIST Digital Library of Mathematical Functions_ ; U.S. Department of Commerce, National Institute of Standards and Technology: Washington, DC, USA; Cambridge University Press: Cambridge, UK, 2010; With 1 CD-ROM (Windows, Macintosh and UNIX). MR 2723248 (2012a:33001).
  * [2] H. Rademacher (1973) _Topics in Analytic Number Theory_. Springer-Verlag, New York.
  * [3] Bell, W. W, _Special Functions for Scientists and Engineers_. London ; Princeton, N.J., [etc.], Van Nostrand, 1968.
  * [4] Oldham, K.B.; Myland, J.C.; Spanier, J. _An Atlas of Functions: With Equator, the Atlas Function Calculator_ , 2nd ed.; Springer: New York, NY, USA, 2009.
  * [5] Andrews, G., Askey, R., & Roy, R. (1999). _Special Functions (Encyclopedia of Mathematics and its Applications)_. Cambridge: Cambridge University Press. doi:10.1017/CBO9781107325937
  * [6] Prudnikov, A. P. , Brychkov, Yu. A. and Marichev, O. I., _Integrals and Series: Elementary Functions_ , Vol. 1. Gordon & Breach Science Publishers, New York, (1986a).
  * [7] Reynolds, R.; Stauffer, A. A Method for Evaluating Definite Integrals in Terms of Special Functions with Examples. _Int. Math. Forum_ 2020, _15_ , 235–244, doi:10.12988/imf.2020.91272
  * [8] Erdéyli, A.; Magnus, W.; Oberhettinger, F.; Tricomi, F.G. _Higher Transcendental Functions_ ; McGraw-Hill Book Company, Inc.: New York, NY, USA; Toronto, ON, Canada; London, UK, 1953; Volume I.
  * [9] Gradshteyn, I.S.; Ryzhik, I.M. _Tables of Integrals, Series and Products_ , 6th ed.; Academic Press: Cambridge, MA, USA, 2000\.
  * [10] Gelca, Răzvan., Andreescu, Titu., _Putnam and Beyond_. Germany: Springer International Publishing, 2017.
  * [11] Reynolds, R., & Stauffer, A. (2022). , _A Note on the Infinite Sum of the Lerch function_. European Journal of Pure and Applied Mathematics, 15(1), 158–168. https://doi.org/10.29020/nybg.ejpam.v15i1.4137
  * [12] Nielsen, Niels. _”Handbuch der Theorie der Gammafunktion.”_ , Leipzig : Teubner, 1906,
  * [13] E. R. Hansen (1975) _A Table of Series and Products_. Prentice-Hall, Englewood Cliffs, N.J.
  * [14] Mark W. Coffey. _”Series representations for the Stieltjes constants.”_ , Rocky Mountain J. Math. 44 (2) 443 - 477, 2014. https://doi.org/10.1216/RMJ-2014-44-2-443
  * [15] Choi, Junesang., Srivastava, Hari Mohan. _Zeta and Q-Zeta Functions and Associated Series and Integrals_. Netherlands: Elsevier Science, 2011.


# Probabilistic Trajectory Prediction with Structural Constraints

Weiming Zhi1,∗, Lionel Ott2, Fabio Ramos1,3 Correspondence to: W. Zhi,
<EMAIL_ADDRESS>School of Computer Science, the University of
Sydney, Australia2 Autonomous Systems Lab, ETH Zurich, Zurich, Switzerland3
NVIDIA, USA

###### Abstract

This work addresses the problem of predicting the motion trajectories of
dynamic objects in the environment. Recent advances in predicting motion
patterns often rely on machine learning techniques to extrapolate motion
patterns from observed trajectories, with no mechanism to directly incorporate
known rules. We propose a novel framework, which combines probabilistic
learning and constrained trajectory optimisation. The learning component of
our framework provides a distribution over future motion trajectories
conditioned on observed past coordinates. This distribution is then used as a
prior to a constrained optimisation problem which enforces chance constraints
on the trajectory distribution. This results in constraint-compliant
trajectory distributions which closely resemble the prior. In particular, we
focus our investigation on collision constraints, such that extrapolated
future trajectory distributions conform to the environment structure. We
empirically demonstrate on real-world and simulated datasets the ability of
our framework to learn complex probabilistic motion trajectories for motion
data, while directly enforcing constraints to improve generalisability,
producing more robust and higher quality trajectory distributions.

##  I Introduction

Autonomous robots, such as service robots and self-driving vehicles, are often
required to coexist in environments with other dynamic objects. For robots to
safely navigate and plan in an anticipatory manner in dynamic environments,
accurate motion predictions for other nearby objects are needed. Many recent
developments [1, 2, 3] in motion prediction have relied on learning based
approaches, whereby future positions, or probabilities over future positions,
are learned from motion trajectories of previously observed dynamic objects.
Learning-based methods to motion prediction typically extract patterns from
the training data. These patterns are then used to extrapolate future
positions based on new observations without explicitly considering the
feasibility of the prediction.

However, there are often rules for motion trajectories that are known a
priori, but are not explicitly enforced during extrapolation by learning-based
models. In particular, we address the problem of incorporating constraints
that arise from the environment’s structure, such as obstacles, into motion
predictions. In most environments, motion trajectories are not highly
unstructured, and depend closely on the environment occupancy structure.

To address the shortcomings of purely learning-based methods, we view motion
prediction as a combined probabilistic trajectory learning and constrained
trajectory optimisation problem. Constrained optimisation methods allow for
constraints to be imposed directly onto the distribution of trajectories
without incorporating the collected data. We propose a novel framework that
incorporates structural constraints into probabilistic motion prediction
problems. Our framework leverages the ability of learning based models to
extract patterns from data, as well as the ability of constrained trajectory
optimisation approaches to explicitly specify constraints, giving better
generalisability and higher quality extrapolations. Specifically, our
contributions consist of: (1) a method to learn a distribution over
trajectories from observations that is amenable to being optimised; (2) an
approach to enforce constraints on distributions of trajectories, particularly
collision chance constraints when provided an occupancy map.

Figure 1: A simple example of enforcing predictions to be consistent with
environment: Provided an observed trajectory segment (black), a purely
learning method may predict collision-prone distributions of trajectories
(blue). Our framework formulates an optimisation problem with chance
constraints arising from environmental structure, as given by an occupancy
map. We solve to obtain predictions which are compliant with the environment,
while remaining similar to the learned future.

##  II Related Work

Motion Trajectory Prediction. Motion prediction is a problem central to robot
autonomy. The most simple methods of motion prediction are physics-based
models, such as constant velocity/acceleration models [4]. More complex
physics-based models often fuse multiple dynamic models, or select a best fit
model from a set [5]. Recent developments in machine learning have led to
learning-based methods that learn motion patterns [1, 6] from a dataset of
observed trajectories. In particular, long-short term memory (LSTM) and
generative adversarial network (GAN)-based models [1, 7, 8, 2] have gained
popularity. Methods of learning entire trajectories, similar to our learning
component, but with relaxed dependencies, is outlined in [3, 9]. Further
attempts to incorporate static obstacles in learning attempt to do so as a
part of the learning pipeline [10] and attempts to extrapolate behaviour.
Approaches of this kind cannot guarantee that valid solutions satisfy
specified chance constraints, and there have not been approaches that directly
incorporate occupancy maps, which are typically built by robots from depth
sensor data.
Trajectory Optimisation. Trajectory optimisation has been successfully applied
in motion planning, by directly expressing obstacles as constraints to solve
for collision-free paths between start and goal points. Popular methods
include _Trajopt_ [11], _CHOMP_ [12], _STOMP_ [13], and _GPMP2_ [14]. Trajopt
utilises sequential quadratic programming to find locally optimal trajectories
compliant with constraints. CHOMP exploits obstacle gradients for efficient
optimisation. Further work on gradient-based trajectory optimisation for
planning can be found in [15, 16]. Although trajectory optimisation is a
popular choice to obtain collision-free trajectories for motion-planning,
there has yet to be extensive work in utilising trajectory optimisation for
motion prediction. The optimisation required in this work is also different in
that we optimise to obtain distributions of trajectories, rather than a single
trajectory as required in motion planning.
Combined Learning and Optimisation. There has been previous explorations of
combining learning and trajectory optimisation. [17] introduce a framework
which combines a known analytical cost function with black box functions in
reinforcement learning for manipulation. [18] introduces CLAMP, a framework
that combines learning from demonstration with optimisation based motion
planning. Theoretical aspects of using learning in conjunction with
optimisation, in a Bayesian framework known as posterior regularisation has
also been investigated in [19]. All of these methods aim to add more structure
to be encoded in learning problems. Similarly our proposed framework allows
for encoding of structure, but specifically for the trajectory prediction
problem, where we are required to do learning and optimisation on
distributions of trajectories with respect to environment occupancy.

##  III Continuous-time Motion Trajectories

Our framework operates over continuous motion trajectories, capable of being
queried at arbitrary resolution, without forward simulating. In this section,
we briefly introduce the continuous-time motion trajectory representation used
in this work. We restrict ourselves to the 2-D case.

Similar to methods described in [16, 15, 3], we represent motion trajectories
as functions $\bm{\xi}:[0,T]\rightarrow\mathbb{R}^{2}$. $T$ represents the
prediction time horizon. To we construct our trajectories as dot products
between weight vectors and _reproducing kernels_. We limit our discussion to
squared exponential radial basis function (RBF) kernels, though methods in
this study generalise to various kernels with minimal modifications. Using
RBFs result in smooth trajectories, and have shown good empirical convergence
properties for functional trajectory optimisation [15]. A trajectory,
$\bm{\xi}(t)$, can be expressed as:

$\displaystyle\mathbf{\xi}(t)\\!=\begin{bmatrix}\mathbf{w}^{x^{\top}}\bm{\phi}(t)\\\
\mathbf{w}^{y^{\top}}\bm{\phi}(t)\end{bmatrix},$
$\displaystyle\mathbf{\phi}(t)\\!=\\!\left[\phi(t,t^{\prime}_{1}),\dots,\phi(t,t^{\prime}_{M})\right]^{\top},$
(1)

where the basis functions are given as,
$\phi(t,t^{\prime})\\!=\\!\exp(-\gamma||t-t^{\prime}||^{2})$. The coordinates
of a dynamic object at time, $t\in[0,T]$, is given as $x(t)$ and $y(t)$,
$\mathbf{w}^{x}$, $\mathbf{w}^{y}\in\mathbb{R}^{M}$ are weight parameters,
$\bm{\phi}(t)\in\mathbb{R}^{M}$, and
$\bm{t^{\prime}}=[t^{\prime}_{1},t^{\prime}_{2},\ldots,t^{\prime}_{M}]^{\top}$
are $M$ pre-define time points at which RBF kernels are centered. $\gamma$ is
a length-scale hyper-parameter which determines the impact a time point
neighbouring times on the trajectory.

Motion data typically comes in the form of discrete sequences of timestamped
coordinates $\\{t_{i},x_{i},y_{i}\\},i=0,\dots,N$. To construct a best-fit
continuous trajectory of discrete sequences from $N$ timestamped coordinates,
we solve the ridge regression problem:

$\displaystyle\min\limits_{\mathbf{w}^{x},\mathbf{w}^{y}}\sum\limits_{n=1}^{N}\big{\\{}(x_{n}-{\mathbf{w}^{x}}^{\top}\bm{\phi}(t_{n}))^{2}$
$\displaystyle+(y_{n}-{\mathbf{w}^{y}}^{\top}\bm{\phi}(t_{n}))^{2}$
$\displaystyle+\ldots\lambda(||\mathbf{w}^{x}||^{2}+||\mathbf{w}^{y}||^{2})\big{\\}},$
(2)

where $\lambda$ is a regularisation hyper-parameter, with the set of $M$
uniform times $\bm{t}^{\prime}$ selected _a priori_ , as centers where the
kernels are fixed. Evaluating the closed form ridge regression solution
results in weight vectors $\mathbf{w}^{x}$ and $\mathbf{w}^{x}$ which forms
the trajectory
$\bm{\xi}(t)=[{\mathbf{w}^{x}}^{\top}\bm{\phi}(t),{\mathbf{w}^{y}}^{\top}\bm{\phi}(t)]$.
We denote the stacked vector of $\mathbf{w}^{x}$ and $\mathbf{w}^{y}$ as
$\mathbf{w}$. As future motion trajectories are inherently uncertain, the
motion prediction problem often requires predicting a distribution over future
trajectories a dynamic object may follow. We denote distributions of
trajectories as $\bm{\Xi}$, with the vector of probabilistic weight
parameters, $\bm{\omega}$, and the parameters of the distribution over
$\bm{\omega}$ is denoted as $\Psi$. $\bm{\Xi}(t)$ gives the distribution of
coordinate position at time $t$. Every sample function of $\bm{\Xi}$ is a
trajectory $\bm{\xi}$, and samples of vector of random variables $\bm{\omega}$
give $\mathbf{w}$. Denoting the parameter distributions associated with $x$
and $y$ as $\bm{\omega}_{x}$ and $\bm{\omega}_{y}$, distributions of
trajectories can be constructed as
$\bm{\Xi}=[{\bm{\omega}^{x}}^{\top}\bm{\phi}(t),{\bm{\omega}^{y}}^{\top}\bm{\phi}(t)]$.
In practice, to obtain $\bm{\xi}$ we can also sample $\mathbf{w}$ from the
distribution $\bm{\omega}$, and take dot products with $\bm{\phi}$.

Training data Learning Model Trajectory weight distribution
$\bm{\omega}_{\Psi_{0}}$ Occupancy map Construct cost operator $\mathcal{C}$
Optimisation to obtain $\bm{\omega}_{\Psi}^{*}$, via eq. 4 Constraint-
compliant trajectory distribution, $\bm{\Xi}_{\Psi}^{*}$ Sample possible
future trajectories $\bm{\xi}$ Kernel features $\bm{\phi}$ Train Predict
$\mathcal{C}$ $\bm{\omega}_{\Psi_{0}}$

Figure 2: Overview of our proposed framework. We combine learning and
optimisation to generalise data and enforce constraints.

##  IV Motion Prediction with Environmental Constraints

###  IV-A Problem Overview

Figure 3: (Left) Examples from real-world dataset [20]. Historical
trajectories in blue, and future in red. (Right) Provided a new unseen
trajectory (blue), extrapolate futures (red). Learned future trajectories may
collide with obstacles. Our framework tackles this by directly enforcing
constraints on trajectory distributions, such that the trajectories are
collision free.

Our goal is to reason about the possible future trajectories of dynamic
objects, based on the immediate past trajectories of the objects, along with
occupancy information from the environment. In this work, we limit our
discussion to trajectories in 2-D, which usually arise from the tracking of
humans and vehicles in video data. We are assumed to have a dataset containing
coordinates of an object up to a time, along with its coordinates thereafter,
and an occupancy map of the environment. Examples from a real-world dataset
are shown in fig. 3, as well as a prediction given by a learning-based model.
Given the observed historical coordinates, we aim to extrapolate the
distribution of trajectories representing likely future motion thereafter. We
require the distribution of trajectories compliant with the environment
structure indicated in the occupancy map, such that the probability of
hypothesising infeasible trajectories to be constrained to a small value.

###  IV-B Framework Overview

An overview of our framework is shown in fig. 2. We learn from our training
data to predict a vector of random variables of weights
$\bm{\omega}_{\Psi_{0}}$ with distribution parameters $\Psi_{0}$, and produce
a prior distribution of trajectories, $\bm{\Xi}_{0}$, as described in eq. 1.
Our aim is to find an alternative random vector $\bm{\omega}_{\Psi}^{*}$,
defined by optimised parameters $\Psi^{*}$, which satisfies enforced
constraints, while remaining close to the learned distribution
$\bm{\omega}_{\Psi_{0}}$. The optimisation objective uses
$\bm{\omega}_{\Psi_{0}}$, while constraints are provided through the obstacle
cost operator, $\mathcal{C}$, which requires an occupancy map of the
environment. The desired random vector $\bm{\omega}_{\Psi}^{*}$ defines a
constraint-compliant distribution of trajectories, $\bm{\Xi}_{\Psi}^{*}$,
provided pre-defined RBF features $\bm{\phi}$, by eq. 1. Individual
constraint-compliant future trajectories, $\bm{\xi}$, can then be realised
from $\bm{\Xi}_{\Psi}^{*}$.

To find the constrain-compliant parameters $\Psi^{*}$, we formulate the
optimisation problem in a similar manner to Bayesian posterior regularisation
methods [19]. The desired parameters $\Psi^{*}$ is the solution to the
constrained optimisation problem:

$\displaystyle\min_{\Psi}\quad\mathrm{D}_{KL}(\bm{\omega}_{\Psi}\lVert\bm{\omega}_{\Psi_{0}})$
(3) $\displaystyle\textrm{s.t.}\quad\mathcal{C}(\bm{\Xi}_{\Psi})-\epsilon\leq
0$ (4)

where $\mathrm{D}_{KL}$ is the Kullback-Leibler (KL) divergence [21], and
$\mathrm{D}_{KL}(\bm{\omega}_{\Psi}\lVert\bm{\omega}_{\Psi_{0}})$ gives us a
definition of “close-ness” between our predicted $\bm{\omega}_{\Psi_{0}}$ and
optimised $\bm{\omega}_{\Psi}$ distributions. $\mathcal{C}(\bm{\Xi}_{\Psi})$
is the obstacle cost operator, which maps distributions of trajectories
$\bm{\Xi}_{\Psi}$ to a scalar cost value. The cost is interpreted as the time-
averaged probability of generating a trajectory that collides with obstacles.
We designate $\epsilon$ to be the allowed limit of the obstacle cost.
Intuitively, our optimisation problem returns the closest distribution of
trajectories, as defined by the KL divergence term, to that given by our
predictive model, subject to the obstacle collision constraint. By separating
the learning and optimisation, we disentangle the complexities of the learning
model and the structural constraints. Our learned $\bm{\omega}_{\Psi_{0}}$
also provides us with a “warm-start” initial solution to our optimisation. In
contrast to learning-based approaches, the constraints in our optimisation
problem is enforced in valid solutions. The main challenges to solving this
trajectory optimisation problem lies in how to obtain the vector of random
variables $\bm{\omega}_{\Psi_{0}}$, which defines our prior trajectory
distribution $\bm{\Xi}_{0}$, as well as how to derive the obstacle cost
operator $\mathcal{C}$ required for optimisation. These are addressed in the
following subsections.

###  IV-C Learning Distributions of Trajectories

In this subsection, we outline a method of learning a mapping from observed
historical trajectories to prior distributions of future trajectories,
parameterised by $\bm{\omega}_{\Psi_{0}}$. The resulting prior distributions
of trajectories are to be used directly in the optimisation problem eq. 4. We
model the distribution of trajectories, $\bm{\Xi}$, as mixtures of multi-
output stochastic processes. To capture the correlation between x and y
coordinates, the distribution over weight parameters are assumed to be
mixtures of matrix normal distributions [22]. A matrix normal distribution is
the generalisation of the normal distribution to matrix-valued random
variables, where each sample from a matrix normal distribution gives a matrix.
We work with the matrix of random variables defined by vertically stacking the
parameters of coordinate dimensions, denoted as
$\mathcal{W}=[\bm{\omega}^{x},\bm{\omega}^{y}]\in\mathbb{R}^{M\times 2}$:

$\displaystyle\mathcal{W}=\begin{bmatrix}{w}^{x}_{1}&{w}^{y}_{1}\\\
\smash{\vdots}&\smash{\vdots}\\\ {w}^{x}_{M}&{w}^{y}_{M}\end{bmatrix},$
$\displaystyle p(\mathcal{W})$
$\displaystyle=\sum_{r=1}^{R}\alpha_{r}\mathrm{MN}\Big{(}\mathcal{M}_{r},\mathcal{U}_{r},\mathcal{V}_{r}\Big{)}.$
(5)

We assume that there are $R$ mixture components, similar to a Gaussian mixture
model [21], with each component being a matrix normal distribution, and
associated weight components $\alpha_{r}$. Matrix normal distributions have
location matrices
$\mathcal{M}_{r}=[\bm{\mu}_{r}^{x},\bm{\mu}_{r}^{y}]\in\mathbb{R}^{M\times
2}$, and positive-definite scale matrices
$\mathcal{U}_{r}\in\mathbb{R}^{M\times M}$,
$\mathcal{V}_{r}\in\mathbb{R}^{2\times 2}$ as parameters [23].

We utilise a neural network to learn a mapping from a fixed-length vectorised
coordinate sequences, $\bm{\varphi}$, that encodes observed partial
trajectories, to component weights $\bm{\alpha}$, along with location matrices
$\mathcal{M}_{r}$ the lower triangular $\mathcal{V}_{r}^{\frac{1}{2}}$ and
diagonal $\mathcal{U}_{r}^{\frac{1}{2}}$ matrices. $\mathcal{V}$ is positive-
definite and will be constructed by taking the product of two lower triangular
matrices. The dependence between different times on the trajectory are largely
captured the RBF features in the time domain $\bm{\phi}(t)$, as described in
section III. Hence, it is sufficient for $\mathcal{U}$ to be a positive
definite diagonal matrix. We obtain the input vector $\bm{\varphi}$ by taking
a 10 time-step window of the latest observed coordinates, and vectorising to
obtain a vector of 20 dimensions. The neural network used is relatively light-
weight, consisting of 4 Relu-activated dense layers with number of neurons:
$(15\times M\times R)\rightarrow(5\times M\times R)\rightarrow(5\times M\times
R)\rightarrow(R\times[3M+3])$. In this work, our main goal is to outline a
representation of future trajectories outputted by the network, such that the
result can be directly optimised. Our method does not preclude alternative
encoding schemes of conditioned history trajectories to be inputted into the
network, such as features from LSTM networks [24] or trajectory features
outlined in [3].

Provided a dataset of history and future trajectories, given by $N$ pairs of
timestamped coordinate sequences. After flattening the history trajectory
sequence to obtain $\bm{\varphi}$, and applying eq. 2 on trajectory futures to
obtain weight matrices, we have pairs denoted as
$\\{\bm{\varphi}_{n},\mathcal{W}_{n}\\}_{n=1}^{N}$. We learn a mapping from
$\bm{\varphi}_{n}$ to parameters
$\\{\alpha_{r},\mathcal{M}_{r},\mathcal{V}_{r},\mathcal{U}_{r}\\}_{r=1}^{R}$
with the fully-connected neural network, with the negative log-likelihood of
the weight matrices distributions as the loss function:

$\displaystyle\mathcal{L}$ $\displaystyle=-\sum_{n=1}^{N}\log$
$\displaystyle\dots\sum_{r=1}^{R}\alpha_{r}\frac{\exp\\{-\frac{1}{2}\mathrm{tr}[{\mathcal{V}_{r}}^{-1}\\!(\mathcal{W}_{n}-\mathcal{M}_{r})^{\top}\\!\mathcal{U}_{r}(\mathcal{W}_{n}-\mathcal{M}_{r})]\\}}{(2\pi)^{M}|\mathcal{V}_{r}||\mathcal{U}_{r}|^{\frac{M}{2}}}.$
(6)

and by applying the activation functions:

$\displaystyle\mathrm{diag}[(\mathcal{U}_{r})^{\frac{1}{2}}]$
$\displaystyle=\exp(z^{\mathrm{diag}(\mathcal{U}_{r})}),$ (7)
$\displaystyle\mathrm{LowerTrig}[(\mathcal{V}_{r})^{\frac{1}{2}}]$
$\displaystyle=z^{\mathrm{LowerTrig}(\mathcal{V}_{r})},$ (8)
$\displaystyle\mathrm{diag}[(\mathcal{V}_{r})^{\frac{1}{2}}]$
$\displaystyle=\exp(z^{\mathrm{diag}(\mathcal{V}_{r})}),$ (9)
$\displaystyle\alpha_{r}$ $\displaystyle=\mathrm{softmax}(z^{\alpha_{r}}),$
(10)

where
$z^{\alpha_{r}},z^{\mathrm{LowerTrig}(\mathcal{V}_{r})},z^{\mathrm{diag}(\mathcal{V}_{r})},z^{\mathrm{diag}(\mathcal{U}_{r})}$
are neural network outputs, and $\mathrm{diag}(\cdot)$,
$\mathrm{LowerTrig}(\cdot)$ indicate the matrix diagonal, and strictly lower
triangular matrix elements respectively. The $\exp$ activation function
guarantees the validity of the scale matrices, by enforcing
$z^{\mathrm{diag}(\mathcal{V}_{r})}>0,z^{\mathrm{diag}(\mathcal{U}_{r})}>0$.
We also enforce $\sum_{r=1}^{R}\alpha_{r}=1$ using the $\mathrm{softmax}$
function. After we train the neural network, provided a vectorised historical
trajectory, $\bm{\varphi}$, the neural network gives the corresponding random
vector of weights of our prior trajectory distribution,
$\bm{\omega}_{\Psi_{0}}=\mathrm{vec}(\mathcal{W})$, where
$\mathrm{vec}(\cdot)$ is the vectorise operator. The prior distributions of
trajectories is then the dot product of the trajectory parameters and the RBF
features, $\bm{\phi}(t)$, as described in eq. 1.

###  IV-D Constraints on Distributions of Trajectories

After obtaining the weights, $\bm{\omega}_{\Psi_{0}}$, of a prior distribution
of trajectories via learning, our main challenge is to design a cost operator,
$\mathcal{C}(\bm{\Xi})$, to encode our desired constraints. Specifically we
wish to design $\mathcal{C}(\bm{\Xi})$ to account for collisions of
distributions of trajectories, provided an occupancy map.
$\mathcal{C}(\bm{\Xi})$ maps a distribution of smooth trajectories,
$\bm{\Xi}$, to a real-valued scalar cost, which quantifies how collision-prone
a given $\bm{\Xi}$ is. $\mathcal{C}$ differs from the cost operator in
trajectory optimisation motion planning problems, such as those in [12, 15,
16], in that the input is a distribution over trajectories rather than a
single trajectory. Our optimisation can not only optimise over trajectory
mean, but also over trajectory variances.

Like trajectory optimisation methods [12], we assume the cost operator takes
the form of average costs over prediction horizon $T$, where the cost at a
given point in time, $\mathcal{C}(\bm{\Xi}_{\Psi}(t))$, is the probability of
a collision at a given time. $\bm{\Xi}_{\Psi}$ contains weight variables with
a distribution parameterised by $\Psi$. We can then interpret
$\mathcal{C}(\bm{\Xi}_{\Psi})$ as the average probability of collision over a
time horizon, defined as:

$\displaystyle\mathcal{C}(\bm{\Xi}_{\Psi})=$
$\displaystyle\frac{1}{T}\int\displaylimits_{0}^{T}\mathcal{C}(\bm{\Xi}_{\Psi}(t))\mathrm{d}t$
(11) $\displaystyle=$
$\displaystyle\frac{1}{T}\int\displaylimits_{0}^{T}\int\displaylimits_{\mathbb{R}^{2}}p(m=1|\mathbf{x})p(\mathbf{x}|\Psi,t)\mathrm{d}\mathbf{x}\mathrm{d}t.$
(12)

The probability of collision at some time and location is computed by taking
the product of the probability of the trajectory reaching the coordinate at
the given time and the probability that the coordinate is occupied. We are
provided with an occupancy map and the probability a coordinate of interest,
$\mathbf{x}\in\mathbb{R}^{2}$, being occupied is denoted as
$p(m=1|\mathbf{x})$. To evaluate the cost, we could marginalise over the all
trajectory parameters, but this may be difficult due to the high
dimensionality of $\bm{\omega}$, instead work in the 2-D world space, where
the environmental constraints are defined. As $\bm{\Xi}_{\Psi}$ is constructed
by taking the dot product of weight vectors and feature vectors, trajectory
coordinate distribution in world space at time $t$ is:

$\displaystyle p(\bm{\Xi}_{\Psi}(t))$
$\displaystyle\\!=\\!\sum_{r=1}^{R}\\!\alpha_{r}\overbrace{{\mathrm{MN}(\mathcal{M}_{r}^{\top}\bm{\phi}(t),\mathcal{V}_{r},\bm{\phi}(t)^{\top}\mathcal{U}_{r}\bm{\phi}(t)\\!)}}^{\text{matrix
normal distribution}}$ (13)
$\displaystyle\\!=\\!\sum_{r=1}^{R}\\!\alpha_{r}\overbrace{\mathcal{N}(\underbrace{\mathcal{M}_{r}^{\top}\bm{\phi}(t)}_{=:\bm{\mu}_{r}},\underbrace{\bm{\phi}(t)^{\top}\mathcal{U}_{r}\bm{\phi}(t)\mathcal{V}_{r}}_{=:\Sigma_{r}}).}^{\text{multivariate
normal}}$ (14)

By the equivalences between matrix normal and multivariate normal
distributions [22], we have a mixture of two-dimensional multivariate normal
distributions in world-space,
$\bm{\Xi}_{\Psi}(t)\sim\sum_{r=1}^{R}\alpha_{r}\mathcal{N}(\bm{\mu}_{r},\Sigma_{r})$.
This is the probability distribution over trajectory coordinates at time $t$,
$p(\mathbf{x}|\Psi,t)$. The mixture means and covariances of each component
are defined as ${\bm{\mu}_{r}:=\mathcal{M}_{r}^{\top}\bm{\phi}(t)},$ and
$\Sigma_{r}:={\bm{\phi}(t)^{\top}\mathcal{U}_{r}\bm{\phi}(t)\mathcal{V}_{r}}.$
Substituting eq. 14 into the inner integral of eq. 12, gives the cost operator
equation at given time:

$\displaystyle\mathcal{C}(\bm{\Xi}_{\Psi}(t))=\sum_{r=1}^{R}\alpha_{r}$
$\displaystyle\ldots\int\displaylimits_{\mathbb{R}^{2}}\Big{\\{}\frac{\exp[-\frac{1}{2}(\mathbf{x}-\bm{\bm{\mu}_{r}})^{\top}\Sigma_{r}^{-1}\\!(\mathbf{x}-\bm{\mu}_{r})]}{(2\pi)|\Sigma_{r}|^{\frac{1}{2}}}\Big{\\}}\underbrace{p(m=1|\mathbf{x})}_{\text{From
map}}\mathrm{d}\mathbf{x}$ (15)

Note that $\bm{\mu}_{r}$ and $\Sigma_{r}$ are dependent on
$\mathcal{M}_{r},\mathcal{U}_{r},\mathcal{V}_{r}$. During the optimisation, we
assume that each mode can be optimised independently and that the importance
of each mode, represented by $\alpha$, is accurately captured by the learning
model and does not change during trajectory optimisation. We then define the
set of parameters to optimise as
$\Psi=\big{\\{}\mathcal{M}_{r},\mathcal{U}_{r},\mathcal{V}_{r}\\}_{r=1}^{R}$.
The integral defined in eq. 15 may be intractable to evaluate analytically, so
we aim to find an approximation to it. The coordinate distribution at a given
time in world space, $\bm{\Xi}_{\Psi}(t)$, is a mixture of bi-variate Gaussian
distributions. We use Gauss-Hermite quadrature, suitable for numerically
integrating exponential-class functions, to approximate eq. 15. Quadrature
methods approximate an integral by taking the weighted sum of integrand values
sampled at computed points, called abscissae. We introduce the change of
variables, $\bm{z_{r}}=\frac{1}{\sqrt{2}}L_{r}^{-1}(\mathbf{x}-\bm{\mu}_{r})$,
where the Cholesky factorisation gives $\Sigma_{r}=L_{r}L_{r}^{T}$. Using
Gauss-Hermite quadrature schemes, we have:

$\displaystyle\mathcal{C}(\bm{\Xi}_{\Psi}(t))$
$\displaystyle\approx\sum_{r=1}^{R}\alpha_{r}\pi^{-1}\sum_{i=1}^{I}\sum_{j=1}^{J}\beta^{i}\beta^{j}p(m=1|\sqrt{2}L\bm{z}+\bm{\mu}_{r}),$
(16)

where we find the values of $\bm{z}^{i,j}=[z^{i},z^{j}]^{\top}$ as abscissa
obtained from roots of Hermite polynomials, and $\beta^{i}$, $\beta^{j}$ are
the associated quadrature weights. Details of abscissae and weight
calculations, as well as a review on multi-dimensional Gauss-Hermite
quadrature, can be found in [25]. An example of abscissae points when
estimating a multi-variate Gaussian is shown in fig. 4. Next, we evaluate the
outer integral in eq. 12 using rectangular numerical integration for
${\mathcal{C}(\bm{\Xi}_{\Psi})=\int_{0}^{T}\mathcal{C}(\bm{\Xi}_{\Psi}(t))}\mathrm{d}t$.
With $\mathcal{C}(\bm{\Xi}_{\Psi})$ defined, we are in a position to solve the
combined learning and constrained optimisation problem, defined earlier in eq.
4, using a non-linear optimisation solver, such as the sequential least
squares optimiser SLSQP [26]. Analytical gradients with respect to occupancy
is provided when using a continuous differential occupancy map such as those
introduced in [27, 28]. By utilising the equivalence between matrix normal
distributions and multivariate Gaussians [23], the KL divergence term between
$\bm{\omega}_{\Psi}$ and $\bm{\omega}_{\Psi_{0}}$ components can be computed
[21]. After obtaining the solution of the optimisation we have the optimised
distribution of $\bm{\omega}_{\Psi}^{*}$. By taking dot products between
$\bm{\omega}_{\Psi}^{*}$ and $\bm{\phi}$, we can construct
$\bm{\Xi}_{\Psi}^{*}$ the constraint-compliant trajectory distribution of
futures. Individual realised trajectories $\bm{\xi}$ can then be generated
from $\bm{\Xi}_{\Psi}$.

Although we have focused on collision constraints with respect to environment
structure, the optimisation allows for constraints on velocity or
acceleration. The continuous nature of our trajectories allows time
derivatives of trajectories to be obtained. Specifically, the _derivative
reproducing property_ , ensures the derivatives of smooth kernels also
correspond to kernels [29]. We can reason about the $n^{th}$-order derivatives
of displacement by replacing $\bm{\phi}(t)$ with
$\frac{\partial^{n}\bm{\phi}(t)}{\partial t^{n}}$.

[capbesideposition=right,top,capbesidewidth=0.3]figure

Figure 4: The abscissae (black) for estimating the distribution in world space
with the Gauss-Hermite quadrature scheme

##  V Experimental Evaluation

We empirically evaluate (1) the ability of our learning component to learn
probabilistic representations of future trajectories to provide good prior
trajectories distributions; (2) the ability of the proposed framework to
enforce collision constraints, and its effect on trajectory quality.

###  V-A Datasets and Metrics

We run our experiments on one simulated dataset, and two real-world datasets,
including:

  * •

Simulated dataset (Simulated): This contains simulated trajectories on a
floor-plan;

  * •

THÖR dataset [20] (Indoor): This dataset contains real human trajectories
walking in a room with three large obstacles in the center of the room.

  * •

Lankershim dataset (Traffic): This dataset contains real vehicle data along a
long boulevard. We use a particular busy intersection between $x=[-80,80]$ and
$y=[250,500]$.

Throughout each of our experiments, we ensure that the observed trajectory
conditioned on contains 10 timesteps. The prediction horizons, $T$, for the
simulated, indoor and traffic datasets are 15, 10, and 20 timesteps
respectively. We split long trajectories, such that we have 200 pairs of
partial trajectories in the simulated dataset, while the indoor dataset
contained 871 pairs of partial trajectories, and the traffic dataset contained
6175 pairs. Trajectories that were shorter than the length of the prediction
horizon were filtered out. A train-to-test ratio of 80:20 was used for each
dataset. For our learning model, we set the number of mixture components
$R=2$, hyper-parameters $M=8,10,11$, and $\gamma=0.05,0.1,0.2$ for the
simulated, indoors, and traffic datasets respectively. Maps for the simulated
and indoor datasets are available, and we construct a map for the traffic
dataset by considering the road positions. Our occupancy map is represented as
a continuous Hilbert Map [27], giving smooth occupancy probabilities over
coordinates, and occupancy gradients.

Metrics used in the quantitative evaluation include: (i) Average displacement
error (ADE): The average euclidean distance error between the expectation of
the nearest distribution component of trajectories with ground truth
trajectory; (iii) Final displacement error (FDE): The euclidean distance error
at the end of time horizon between the expectation of the nearest distribution
component of trajectories with ground truth trajectory. (iii) Average
Likelihood (AL): To take into account the uncertain nature of motion predict,
for probabilistic models, we record the likelihood of drawing the trajectory
averaged over the timesteps. Note that likelihoods are difficult to interpret
and not comparable over different datasets. Lower ADE and FDE, and higher AL
indicate better performance. For the evaluation of collisions, we also record
the percentage of constraint-violating trajectory distributions in test set.

###  V-B Learning Continuous-Time Trajectory Distributions as Priors

The trajectory distribution learning component of our framework provides a
prior distribution for the optimisation. We examine whether the learning
component of our framework is able to learn distributions of future
trajectories from data, comparing the performance with several benchmark
models. The benchmark models are: (1) Kernel trajectory maps (KTMs) [3], (2)
Constant velocity (CV) model, and (3) Naive neural-network (NN) model, where
we learn a mapping between past trajectories and future trajectories by
minimising the mean squared error loss. The neural network comprised of Relu-
activated fully-connect layers with number of neurons:
$(560)\rightarrow(180)\rightarrow(180)\rightarrow(N_{t})$, where $N_{t}$ gives
the trajectory time horizon.

We see that the learning component of our framework is able to learn high
quality distributions of trajectories. Figure 5 shows examples of extrapolated
trajectory distributions in red, and the observed trajectory in blue. We see
that the learned distribution is relatively close to the ground truth green
trajectory, and is able to capture the uncertain nature of motion prediction.
Although a non-negligible number of sampled trajectories collide with the
environment. Table I summarises the performance results of our method against
our benchmarks. We see that our learning model is comparable to the method
proposed in [3], while outperforming a constant velocity and naive neural
network model. While the KTM method performs slightly better on the traffic
dataset, it is purely learning-based and lacks mechanisms to enforce
constraints. We note that our method encodes the observed trajectory by simply
flattening the input coordinates, as described in section IV-C. More
sophisticated methods of encoding the observed trajectory (such as with LSTM
networks [24]), independent of our outlined output of trajectory
distributions, could improve the performance of the learning model.

|  | ADE | FDE | AL
---|---|---|---|---
Simulated | Ours | 1.29 | 2.27 | 0.12
KTM [3] | 1.44 | 2.41 | 0.11
CV | 4.27 | 11.66 | -
NN-naive | 2.11 | 6.07 | -
Indoors | Ours | 0.94 | 1.32 | 3.32
KTM | 0.65 | 1.26 | 2.88
CV | 2.08 | 4.66 | -
NN-naive | 1.84 | 3.35 | -
Traffic | Ours | 4.99 | 7.24 | 0.04
KTM | 4.98 | 7.03 | 0.043
CV | 5.14 | 10.38 | -
NN-naive | 21.1 | 39.5 | -
TABLE I: We evaluate the quality of our learned trajectory prior with
benchmark models. Lower ADE and FDE, and higher AL indicate better
performance.

Figure 5: Quantitative evaluations of learned distributions of trajectories
(red), Conditioned observations (blue), ground truth (green). We see that our
learning model gives relatively good priors, but many sampled trajectories
collide with obstacles. The prior distributions are further enforced to be
collision-free (See fig. 6)

Figure 6: Examples from each dataset of the optimisation component of the
framework applying constraints. The top row shows abscissae points of the main
component, with notable gradients with respect to occupancy shown as arrows
(teal) at the abscissae. The bottom row shows the same abscissae points after
optimisation. We see that the constrained distribution subtly conforms better
to structures in the environment, while still closely resembling the learned
trajectory distribution. The left, middle and right sub-figures are examples
from the simulated, indoors, and traffic datasets respectively.

###  V-C Enforcing Occupancy Compliance via Trajectory Optimisation

After obtaining a prior from our learned model, we apply constraints on the
prediction via trajectory optimisation. We will evaluate whether collision
constraints can be enforced, and their effect on trajectory quality. We define
the constraint such that the time-averaged probability of collision to be less
than or equal to 0.05, $\mathcal{C}-0.05\leq 0$. The SLSQP optimiser was able
to solve the optimisation quickly, under half a second, with a python
implementation. We obtain trajectory distribution priors from our learned
model, and select all of the observations where the learned prior trajectory
is constraint-violating. We evaluate the performance of the trajectory priors,
and compare them with the trajectory distribution predictions after enforcing
constraints. The quantitative results are listed in table II. We see after
solving with structural constraints, all predictions are feasible.
Additionally, the trajectory distribution quality improves after enforcing
constraints for all the datasets considered.

|  | ADE | FDE | AL | C.V.P
---|---|---|---|---|---
Simulated | Unoptimised | 1.48 | 2.87 | 0.11 | 40%
| Optimised | 1.28 | 2.33 | 0.17 | 0%
Indoors | Unoptimised | 1.62 | 2.41 | 3.42 | 15.4%
| Optimised | 1.44 | 2.38 | 3.3 | 0%
Traffic | Unoptimised | 6.56 | 11.81 | 0.034 | 5.7%
| Optimised | 5.14 | 8.95 | 0.043 | 0%
TABLE II: We evaluate the quality of our optimised trajectory distribution
with collision constraints, relative to the prior. We see that providing
constraints improves prediction quality across all dataset. The percentage of
trajectory distributions violating constraint, $\mathcal{C}-0.05\leq 0$,
relative to the entire test set is also given, as Constraint Violation
Percentage (C.V.P). After optimisation all trajectory distributions are
constraint-compliant.

Qualitative results are shown in fig. 6. We see the abscissae points of the
largest trajectory distribution component before (top) and after (bottom) the
trajectory optimisation of examples in each dataset. Our framework is able to
ensure that the distributions comply with environmental structure. Notably, we
see that the trajectory optimisation not only adjusts the mean of the
trajectory distribution, but also optimises the variances, giving more robust
trajectory distributions. This effect is most prominent in the left sub-
figure, where the predicted distribution of trajectories in a simulated
hallway recovers a much tighter variance which follows the environment
structure. We also observe instances collision-avoidance with objects in open
space (demonstrated in middle sub-figure), as well as increased adherence to
road structure (right sub-figure).

##  VI Conclusion

We introduced a novel framework to learn motion predictions from observations
while conforming to constraints, such as obstacles. Our framework consists of
a learning component which learns a distribution over future trajectories.
This distribution is used as a prior in a trajectory optimisation step which
enforces chance constraints on the trajectory distribution. We empirically
demonstrate that our framework can learn complex trajectory distributions and
enforce compliance with environment structures via optimisation. This results
in reduced variance and the avoidance of obstacles by the predicted
trajectories, leading to improved prediction quality.

## References

  * [1] A. Alahi, K. Goel, V. Ramanathan, A. Robicquet, F.-F. Li, and S. Savarese, “Social lstm: Human trajectory prediction in crowded spaces,” IEEE Conference on Computer Vision and Pattern Recognition, 2016.
  * [2] A. Gupta, J. E. Johnson, L. Fei-Fei, S. Savarese, and A. Alahi, “Social gan: Socially acceptable trajectories with generative adversarial networks,” IEEE Conference on Computer Vision and Pattern Recognition, 2018.
  * [3] W. Zhi, L. Ott, and F. Ramos, “Kernel trajectory maps for multi-modal probabilistic motion prediction,” Conference on Robot Learning (CoRL), 2019\.
  * [4] R. Schubert, E. Richter, and G. Wanielik, “Comparison and evaluation of advanced motion models for vehicle tracking,” International Conference on Information Fusion, 2008.
  * [5] J. F. P. Kooij, F. Flohr, E. A. I. Pool, and D. Gavrila, “Context-based path prediction for targets with switching dynamics,” International Journal of Computer Vision, 2018.
  * [6] W. Zhi, R. Senanayake, L. Ott, and F. Ramos, “Spatiotemporal learning of directional uncertainty in urban environments with kernel recurrent mixture density networks,” IEEE Robotics and Automation Letters, 2019.
  * [7] P. Zhang, W. Ouyang, P. Zhang, J. Xue, and N. Zheng, “Sr-lstm: State refinement for lstm towards pedestrian trajectory prediction,” IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2019.
  * [8] J. Amirian, J.-B. Hayet, and J. Pettré, “Social ways: Learning multi-modal distributions of pedestrian trajectories with gans,” in CVPR Workshops, 2019.
  * [9] W. Zhi, T. Lai, L. Ott, G. Francis, and F. Ramos, “Octnet: Trajectory generation in new environments from past experiences,” 2019.
  * [10] A. Sadeghian, V. Kosaraju, A. Sadeghian, N. Hirose, H. Rezatofighi, and S. Savarese, “Sophie: An attentive gan for predicting paths compliant to social and physical constraints,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2019.
  * [11] J. Schulman, J. Ho, A. X. Lee, I. Awwal, H. Bradlow, and P. Abbeel, “Finding locally optimal, collision-free trajectories with sequential convex optimization,” in Robotics: Science and Systems, 2013.
  * [12] M. Zucker, N. D. Ratliff, A. D. Dragan, M. Pivtoraiko, M. Klingensmith, C. M. Dellin, J. A. Bagnell, and S. S. Srinivasa, “Chomp: Covariant hamiltonian optimization for motion planning,” The International Journal of Robotics Research, 2013.
  * [13] M. Kalakrishnan, S. Chitta, E. Theodorou, P. Pastor, and S. Schaal, “Stomp: Stochastic trajectory optimization for motion planning,” IEEE International Conference on Robotics and Automation, 2011.
  * [14] J. Dong, M. Mukadam, F. Dellaert, and B. Boots, “Motion planning as probabilistic inference using gaussian processes and factor graphs,” in Robotics: Science and Systems, 2016.
  * [15] Z. Marinho, A. D. Dragan, A. Byravan, B. Boots, S. Srinivasa, and G. J. Gordon, “Functional gradient motion planning in reproducing kernel hilbert spaces,” Robotics: Science and Systems, 2016.
  * [16] G. Francis, L. Ott, and F. Ramos, “Fast stochastic functional path planning in occupancy maps,” International Conference on Robotics and Automation (ICRA), 2019.
  * [17] P. Englert and M. Toussaint, “Learning manipulation skills from a single demonstration,” International Journal of Robotics Research, 2018.
  * [18] M. A. Rana, M. Mukadam, S. R. Ahmadzadeh, S. Chernova, and B. Boots, “Towards robust skill generalization: Unifying learning from demonstration and motion planning,” in Conference on Robot Learning, 2017.
  * [19] K. Ganchev, J. a. Graça, J. Gillenwater, and B. Taskar, “Posterior regularization for structured latent variable models,” J. Mach. Learn. Res., 2010.
  * [20] A. Rudenko, T. P. Kucner, C. S. Swaminathan, R. T. Chadalavada, K. O. Arras, and A. J. Lilienthal, “Thör: Human-robot navigation data collection and accurate motion trajectories dataset,” IEEE Robotics and Automation Letters, 2020.
  * [21] C. M. Bishop, Pattern Recognition and Machine Learning (Information Science and Statistics).  2006\.
  * [22] A. P. Dawid, “Some matrix-variate distribution theory: Notational considerations and a bayesian application,” Biometrika, 1981.
  * [23] P. Dutilleul, “The mle algorithm for the matrix normal distribution,” Journal of Statistical Computation and Simulation, 1999.
  * [24] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural Comput., 1997.
  * [25] Q. Liu and D. A. Pierce, “A note on gauss-hermite quadrature,” Biometrika, vol. 81, 1994.
  * [26] D. Kraft, “A software package for sequential quadratic programming. tech,” tech. rep., DLR German Aerospace Center — Institute for Flight Mechanics, 1988\.
  * [27] F. Ramos and L. Ott, “Hilbert maps: Scalable continuous occupancy mapping with stochastic gradient descent,” The International Journal of Robotics Research, 2016.
  * [28] W. Zhi, L. Ott, R. Senanayake, and F. Ramos, “Continuous occupancy map fusion with fast bayesian hilbert maps,” in 2019 International Conference on Robotics and Automation (ICRA), 2019.
  * [29] D.-X. Zhou, “Derivative reproducing properties for kernel methods in learning theory,” Journal of Computational and Applied Mathematics, 2008.