# A Review of the Unknown Input Observer with Examples

Sam Nazari  Department of Electrical & Computer Engineering, Northeastern
University, Boston, MA,02210<EMAIL_ADDRESS>
###### Abstract.

We consider the Unknown Input Observer and focus on three examples found in
the literature.

###### Key words and phrases:

Unknown Input Observer, Disturbance Estimation

##  1\. INTRODUCTION

This memo considers the design of observers for a class of linear dynamic
systems in which system uncertainty can be modeled as an additive unknown
disturbance term in the dynamic equation:

(1.1) $\displaystyle\dot{x}(t)$ $\displaystyle=Ax(t)+Bu(t)+Ed(t)$
$\displaystyle y(t)$ $\displaystyle=Cx(t)$

where $x(t)\in\mathbb{R}^{n}$ is the state vector, $y(t)\in\mathbb{R}^{m}$ is
the output vector, $u(t)\in\mathbb{R}^{r}$ is the known input vector and
$d(t)\in\mathbb{R}^{q}$ is the unknown input or disturbance vector. If the
unknown input distribution matrix, $E$, is not full column rank, then the rank
decomposition:

$\displaystyle Ed(t)=E_{1}E_{2}d(t)$

can be applied where, $E_{1}$, is a full column rank matrix and $E_{2}d(t)$
can be considered as a new unknown input. The term, $Ed(t)$, may be used to
model additive disturbances or model uncertainties that are not known a
priori. Typical uncertainties are noise, plant nonlinearities, model reduction
errors, parametric uncertainties and difficult to model interconnection terms
in large scale systems. While not considered in this memo, a disturbance term
may also appear in the output equation:

$\displaystyle y(t)=Cx(t)+E_{y}d(t)$

Additionally, this memo does not consider the $Du(t)$ term in the system
output equations. For the sake of brevity and without loss of generality, we
will not treat systems with this term. Note that the development of the
theoretical material in this memo closely follows Chapter 3 of [1].

##  2\. Unknown Input Observer Theory

The problem of designing observers for a linear system with both known and
unknown inputs has been studied since the 1970’s [1], [2], [4]. The motivation
for studying this problem is that in practice, for a variety of reasons, many
plants are modeled with disturbance terms as shown in Eq 1.1. When modeled in
such a way, traditional observers with the Luenberger structure become
difficult to implement since they use all of the input signals in order to
compute an estimate of the state vector. This limitation renders them
ineffective for many practical applications. Contrastingly, the observer
structure considered in this memo assumes no a priori knowledge about the
unknown inputs, rendering it more amenable for use in fault detection
applications. Many investigators have addressed this problem with varying
success and we refer the reader to Section 3.2 of [1] for a broad literature
review. The class of observers examined in this memo is based on the Unknown
Input Observer (UIO) scheme proposed by Chen, Patton and Zhang.

######  Definition 2.1.

An observer is defined as an unknown input observer [3, 2] for the system
described by Eq 1.1, if its state estimation error vector, $e(t)$, defined as:

(2.1) $\displaystyle e(t)=x(t)-\hat{x}(t)$

approaches zero asymptotically regardless of the presence of unknown inputs,
$d(t)$, in the system. Furthermore, the structure for a full order UIO is
given by the dynamic system:

(2.2) $\displaystyle\dot{z}(t)$ $\displaystyle=Fz(t)+TBu(t)+Ky(t)$
$\displaystyle\hat{x}(t)$ $\displaystyle=z(t)+Hy(t)$

Figure 1. UIO Structure where $K=K_{1}+K_{2}$.

where $\hat{x}\in\mathbb{R}^{n}$ is the state estimate, $z\in\mathbb{R}^{n}$
is the state of the full-order dynamic observer, and $F,T,K,H$ are matrices to
be designed for achieving unknown input decoupling. A block diagram of the UIO
is shown in Figure 1. From the block diagram, one can see that the UIO is
essentially a dynamic system that decouples the state estimation dynamics from
the disturbance term in the original system. By expanding $\dot{e}(t)$ one
obtains:

$\displaystyle\dot{e}(t)$ $\displaystyle=(A-HCA-K_{1}C)e(t)+\left[F-(A-HCA-
K_{1}C)\right]z(t)$ $\displaystyle+[K_{2}-(A-HCA-K_{1}C)]y(t)$
$\displaystyle+[T-(I-HC)]Bu(t)+(HC-I)Ed(t)$

and it is easy to see that in order to make the estimation error a function of
$Fe(t)$:

(2.3) $\displaystyle\dot{e}(t)=Fe(t)$

the equations:

(2.4) $\displaystyle 0$ $\displaystyle=(HC-I)E$ (2.5) $\displaystyle T$
$\displaystyle=I-HC$ (2.6) $\displaystyle F$ $\displaystyle=A-HCA-K_{1}C$
(2.7) $\displaystyle K_{2}$ $\displaystyle=FH$

must hold. If all the eigenvalues in $F$ are stable, then $e(t)$ will approach
zero asymptotically. Note that Eq 2.3 is not a function of $E$ or $d$.
Therefore, the estimation error approaches zero independently of the
disturbance terms, achieving the desired decoupling of the state estimate from
the unknown disturbance inputs.

######  Theorem 2.2.

The necessary and sufficient conditions for the system in Eq 2.2 to be an UIO
for the system defined in 1.1 are:

  1. (1)

$rank(CE)=rank(E)$

  2. (2)

$(C,A_{1})$ is a detectable pair, where $A_{1}=A-E[(CE)^{T}CE]^{-1}(CE)^{T}CA$

###### Proof.

Please see [1] ∎

According to [4], the disturbance estimate can now be obtained, if needed, by
the relation:

(2.8) $\displaystyle\hat{d}=(CE)^{\dagger}[\dot{\hat{y}}-CA\hat{x}-CBu]$

Next we provide a step by step design procedure for the UIO that will be
utilized in the examples that ensue.

###  2.1. Design Procedure for UIOs

We will now describe a procedure that can lead to effective UIO design. It can
be shown that Eqn 2.4 is solvable if and only if the condition
$rank(CE)=rank(E)$ is met [1]. Therefore, the rank condition
$rank(CE)=rank(E)$ must first be checked. If this condition is not met, then a
UIO does not exist. If the condition is met, then we may compute the following
observer matrices:

$\displaystyle H$ $\displaystyle=E[(CE)^{T}CE]^{-1}(CE)^{T}$ $\displaystyle T$
$\displaystyle=I-HC$ $\displaystyle A_{1}$ $\displaystyle=TA$

Next, the observability of the pair $(C,A_{1})$ must be checked. As shown in
Eq 2.3, an important step in designing a UIO is to stabilize $F=A_{1}-K_{1}C$.
If the pair $(C,A_{1})$ is detectable, then this can be achieved by using pole
placement to choose the appropriate gain matrix $K_{1}$. It can be shown that
the observability of the pair $(C,A_{1})$ is equivalent to the observability
of the pair $(C,A)$ [3]. Therefore, if the pair $(C,A_{1})$ is observable or
at least detectable then the gain, $K_{1}$, in Eq 2.6 can be obtained for the
UIO. Hence, if $(C,A_{1})$ is observable then a UIO exists and $K_{1}$ should
be computed using pole placement.

On the other hand, if the pair $(C,A_{1})$ is not observable, then one must
construct a transformation matrix, $P$, for the observable canonical
decomposition of the system:

$\displaystyle PA_{1}P^{-1}$
$\displaystyle=\left[\begin{array}[]{cc}A_{11}&0\\\
A_{12}&A_{22}\end{array}\right]\hskip
14.22636ptA_{11}\in\mathbb{R}^{n_{1}\times n_{1}}$ $\displaystyle CP^{-1}$
$\displaystyle=\left[C^{*}\hskip 5.69054pt0\right]\hskip
48.36967ptC^{*}\in\mathbb{R}^{m\times n_{1}}$

where $n_{1}$ is the rank of the observability pair $(C,A_{1})$, and
$(C^{*},A_{11})$ is observable. If any of the eigenvalues of $A_{22}$ are
unstable, then a UIO does not exist. If the eigenvalues of $A_{22}$ are
stable, then $(C,A_{1})$ is detectable and we may select $n_{1}$ desirable
eigenvalues and assign them to $A_{11}-K_{p}^{1}C^{*}$ using pole placement to
obtain $K_{p}^{1}$. Next compute

$\displaystyle K_{1}=P^{-1}K_{p}=P^{-1}[(K_{p}^{1})^{T}\hskip
5.69054pt(K_{p}^{2})^{T}]^{T}$

where $K_{p}^{2}$ can be any $(n-n1)\times m$ matrix, because it does not
affect the eigenvalues of $F$. Finally, compute the observer matrices:

$\displaystyle F$ $\displaystyle=A_{1}-K_{1}C$ $\displaystyle K$
$\displaystyle=K_{1}+K_{2}=K_{1}+FH.$

In order to clarify the design procedure we consider some numerical examples
next.

##  3\. Numerical Examples

###  3.1. Third order system with no inputs

The first example considered is taken from page 76 of [1]. There, the
following parameter matrices are used:

(3.1) $\displaystyle A=\left[\begin{array}[]{rrr}-1&1&0\\\ -1&0&0\\\
0&-1&-1\end{array}\right],\hskip 5.69054ptC=\left[\begin{array}[]{rrr}1&0&0\\\
0&0&1\end{array}\right],\hskip
5.69054ptE=\left[\begin{array}[]{rrr}-1&0&0\end{array}\right]$

First we check that $rank(CE)=rank(E)=1$:

    rank(C*E)
    rank(E)

    ans =       ans =

          1             1

The rank is equal to $1$ as needed. The observer matrices are computed next:

    H = E*inv((C*E)’*(C*E))*(C*E)’
    T = eye(3)-H*C
    A1 = T*A

    H =                 T =                 A1 =

         1     0           0    0   0            0   0   0
         0     0           0    1   0           -1   0   0
         0     0           0    0   1            0  -1  -1

Next, the observability of the pair $(A_{1},C)$ must be checked. If it is full
rank ($rank(obsv(A_{1},C))=3$), then we may apply pole placement:

    rank(obsv(A1,C))

    ans =

         3

Now we proceed to apply pole placement in order to place the observer poles at
$[-2,-10,-5]$.

    K1 = place(A’,C’,[-2,-10,-5])’

    K1 =

        6.7231   -3.7816
        7.5581  -11.3458
       -3.0543    8.2769

Finally, the observer $F$ and $K$ matrices are computed:

    F = A1-K1*C
    K = K1 + F*H

    F =                                  K =

       -6.7231         0    3.7816              0    -3.7816
       -8.5581         0   11.3458             -1.0  -11.3458
        3.0543   -1.0000   -9.2769              0     8.2769

Figure 2. Simulink block diagram of the dynamics subsystem. Figure 3. Simulink
block diagram of the UIO.

Let us examine the Simulation results with a disturbance occurring at $t=10$
seconds. This is showin in Figure 4. The initial conditions of the dynamic
system as:

    x1_0 = 100;
    x2_0 = -100;
    x3_0 = 1;

In Figure 4, one is able to see the state trajectories in the solid lines. The
dotted ’x’ represent the estimate from the UIO. For the given disturbance and
the chosen observer pole locations, the observer seems to converge (remove all
estimation error) within five seconds. Figures 2 and 3 show the block diagram
in detail. Note that we have chosen not to estimate the disturbance itself.
All the files for this simulation can be found in the

    C:\Users\sqn4594\Documents\FDIR

folder. Specifically, the file

     bookEx.m

and the Simulink model

     exOne.xsl

contain everything.

Figure 4. State response and UIO estimate for example

###  3.2. Car Suspension System

The next example is a quarter car model. The system matrices are:

(3.2) $\displaystyle A=\left[\begin{array}[]{rrrr}0&1&0&0\\\
\frac{-ks}{mb}&\frac{-bs}{mb}&\frac{ks}{mb}&\frac{bs}{mb}\\\ 0&0&0&1\\\
\frac{ks}{mw}&\frac{bs}{mw}&\frac{-ks-
kt}{mw}&\frac{-bs}{mw}\end{array}\right]\hskip
5.69054ptB=\left[\begin{array}[]{rr}0&0\\\ 0&\frac{1000}{mb}\\\ 0&0\\\
\frac{kt}{mw}&\frac{-1000}{mw}\end{array}\right]\hskip
5.69054ptC=\left[\begin{array}[]{rrrr}1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
0&0&0&1\end{array}\right]\hskip 5.69054ptE=\left[\begin{array}[]{r}1\\\ 1\\\
1\\\ 1\end{array}\right]$

where $m_{b}$, represents the car body and the mass, $m_{w}$, represents the
wheel assembly. The spring ,$k_{s}$, and damper, $b_{s}$, represent the
passive spring and shock absorber placed between the car body and the wheel
assembly. The spring, $k_{t}$, models the compressibility of the pneumatic
tire. The variables $x_{b}$, $x_{w}$ and $r$ are the car body travel, the
wheel travel, and the road disturbance, respectively. The force, $f_{s}$,
which is applied between the body and the wheel assembly, is controlled by
feedback. We assume that the disturbance force may enter at any of the state
variables, not just the road disturbance. In MATLAB, we assign the following
values to the parameters:

    %----System Physical Constants----%
    mb = 300;    % kg
    mw = 60;     % kg
    bs = 1000;   % N/m/s
    ks = 16000 ; % N/m
    kt = 190000; % N/m

Therefore, the system matrices become:

    % --------System Matrices(It is a SIMO system)---------%
    A = [ 0 1 0 0; [-ks -bs ks bs]/mb ; ...
          0 0 0 1; [ks bs -ks-kt -bs]/mw];
    B = [0 0; 0 10000/mb ; 0 0; [kt -10000]/mw];
    C = eye(4);
    E=[1;1;1;1];

Figure 5. Free Body Diagram of the Quarter Car Suspension model
(idealization).

    A =                                             B =

       1.0e+03 *                                        1.0e+03 *

             0    0.0010         0         0                  0         0
       -0.0533   -0.0033    0.0533    0.0033                  0         0.0333
             0         0         0    0.0010                  0         0
        0.2667    0.0167   -3.4333   -0.0167                  3.1667   -0.1667

    C =                                             E =

         1     0     0     0                              1
         0     1     0     0                              1
         0     0     1     0                              1
         0     0     0     1                              1

Now we can obtain an optimal gain for control of the suspension system by
using LQR below. We will not dwell too much on this part of the design work
since our main focus is the UIO design.

    %---------LQR design----------%
    Q=[.25 0 0 0;0 4 0 0;0 0 1 0;0 0 0 4];
    R=50*eye(2);
    [K1 l s]=lqr(A,B,Q,R);

In order to design the UIO, we first check to see if the $rank(CE)=rank(E)=1$

    rank(C*E)
    rank(E)

    ans =       ans =

         1              1

Figure 6. Root leve Simulink Block diagram architecture of the car suspension
system in example two.

Next we compute the observer matrices

    H = E*inv((C*E)’*(C*E))*(C*E)’
    T = eye(4)-H*C
    A1 = T*A

    f=zeros(4,4)
    v=[-4 -4 -4 -4]
    F=diag(v)

    k1=inv(C)*(A-F-H*C*A)
    k2=F*H
    k=zeros
    k=k1+k2

Figure 7. Quarter car suspension system dynamics implementation in Simulink.
Figure 8. Simulink implementation of the Unknown Input Observer for the
quarter car suspension model.

    H =                                         T =

        0.2500    0.2500    0.2500    0.2500         0.7500   -0.2500   -0.2500   -0.2500
        0.2500    0.2500    0.2500    0.2500        -0.2500    0.7500   -0.2500   -0.2500
        0.2500    0.2500    0.2500    0.2500        -0.2500   -0.2500    0.7500   -0.2500
        0.2500    0.2500    0.2500    0.2500        -0.2500   -0.2500   -0.2500    0.7500

    A1 =                                            f =

       1.0e+03 *

       -0.0533   -0.0026    0.8450    0.0031               0    0   0   0
       -0.1067   -0.0069    0.8983    0.0064               0    0   0   0
       -0.0533   -0.0036    0.8450    0.0041               0    0   0   0
        0.2133    0.0131   -2.5883   -0.0136               0    0   0   0

    v =                                             F = eye(-4)

        -4    -4    -4    -4

    k1 =                                            k2 = 4X4 matrix ’-1’

       1.0e+03 *

       -0.0493   -0.0026    0.8450    0.0031
       -0.1067   -0.0029    0.8983    0.0064
       -0.0533   -0.0036    0.8490    0.0041
        0.2133    0.0131   -2.5883   -0.0096

    k =

       1.0e+03 *

       -0.0503   -0.0036    0.8440    0.0021
       -0.1077   -0.0039    0.8973    0.0054
       -0.0543   -0.0046    0.8480    0.0031
        0.2123    0.0121   -2.5893   -0.0106

Next we form $A1$ and check observability

    rank(obsv(A1,C))

    ans =

         4

Since the observability matrix has full rank, we can go on and simulate the
system to obtain the dynamic estimator performance. The initial conditions are
set to be:

    X_0 = [-1;10;3;5]

    X_0 =

        -1
        10
         3
         5

Figure 9. UIO results for the car suspension model.

Figure 9 shows the state estimation error for the UIO used in the car
suspension model. The observer seems to converge within 1.5 seconds. Figures 8
shows the implementation of the UIO in Simulink. All the files for this
simulation can be found in the FDIR  folder. Specifically, the file

     carsus.m

and the Simulink model

     car_sus.xsl

contain everything.

###  3.3. Coupled Mass-Spring System

Figure 10. The free body diagram for the two mass and spring system.

Consider the coupled spring mass system shown in Figure 10. The coupled mass
spring system is a well studied mechanical system in control theory. The input
to this system is the sinusoidal motion of the end of the rightmost spring,
and the output is the position of each mass, $q_{1}$ and $q_{2}$. The
disturbance input is also shown in the figure as $e$. In state space form, the
equations of motion are:

(3.3) $\displaystyle x=\left[\begin{array}[]{r}q_{1}\\\ q_{2}\\\
\dot{q}_{1}\\\ \dot{q}_{2}\\\ \end{array}\right],\hskip
5.69054ptA=\left[\begin{array}[]{rrrr}0&0&1&0\\\ 0&0&0&1\\\
\frac{-2k}{m}&\frac{k}{m}&\frac{-c}{m}&0\\\
\frac{k}{m}&\frac{-2k}{m}&0&\frac{-c}{m}\end{array}\right],\hskip
5.69054ptB=\left[\begin{array}[]{c}0\\\ 0\\\ 0\\\
\frac{k}{m}\end{array}\right],\hskip
5.69054ptC=\left[\begin{array}[]{rrrr}1&0&0&0\\\ 0&1&0&0\\\
\end{array}\right],\hskip 5.69054ptE=\left[\begin{array}[]{r}0\\\ 0\\\ 0\\\
1\end{array}\right]$

In MATLAB, we assign the following values to the parameters:

    %----------System Parameters---------%
    k   = 1   %
    m   = 1   % kg
    c   = 1 %

    k =         m =         c =

         1          1               1

to obtain the following state space formulation:

    %----------State Space Formulation---------%
    A   = [0 0 1 0;
           0 0 0 1;
           -2*k/m k/m -c/m 0;
           k/m -2*k/m 0 -c/m]

    B   = [0;0;0;k/m]
    C   = [1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1]
    D   = 0
    E   = [0;0;0;1]

    A =                         B =             C =                  D =        E =

         0  0   1   0               0                1  0   0   0        0              0
         0  0   0   1               0                0  1   0   0                       0
        -2  1  -1   0               0                0  0   1   0                       0
         1 -2   0  -1               1                0  0   0   1                       1

In order to control the placement of the masses in the system, we utilized
LQR. The following MATLAB code accomplishes this objective:

Figure 11. The root level Simulink block diagram for the coupled mass spring
dynamic system.

    %----------LQR Control---------%
    Q=[100 0 0 0;
        0 100 0 0;
        0 0 100 0;
        0 0 0 100]
    R=1
    [Klqr ll ss]=lqr(A,B,Q,R)

    Q =                     R =       Klqr =

       100  0    0    0         1           -1.3620   10.4615    2.0165   10.0419
         0  100  0    0
         0  0    100  0
         0  0    0    100

Now we are ready to begin the UIO design. The first step is to check that
$rank(CE)=rank(E)$:

    %----------Step 1: rank(CE) = rank(E) ?= 1---------%
    rank(C*E)
    rank(E)

    ans =     ans =

         1          1

Next we are ready to compute the first set of observer matrices:

    %----------Step 2: Compute observer matrices---------%
    H = E*inv((C*E)’*(C*E))*(C*E)’
    T = eye(4)-H*C
    A1 = T*A

    H =                 T =                     A1 =

         0  0   0   0       1   0   0   0               0   0   1   0
         0  0   0   0       0   1   0   0               0   0   0   1
         0  0   0   0       0   0   1   0              -2   1  -1   0
         0  0   0   1       0   0   0   0               0   0   0   0

Figure 12. The mass spring system dynamics implemented in linearized form in
Simulink.

Now that $A_{1}$ has been obtained, we must check the observability of the
pair $(C,A_{1})$:

    %----------Step 3:Check (C,A1) rank---------%
    rank(obsv(A1,C))

    ans =

         4

Since the pair $(C,A_{1})$ is full rank, it is observable. Therefore, we can
use pole placement for the observer gain $K_{1}$:

    %----------Step 4: Choose observer poles---------%
    K1 = place(A’,C’,[-2,-10,-5,-3])’

    K1 =

         2     0     1     0
         0    10     0     1
        -2     1     4     0
         1    -2     0     2

Figure 13. The mass spring system dynamics implemented in linearized form in
Simulink.

and to conclude the observer design we compute the $F$ and $K$ matrices:

    %----------Step 5: Finish observer design---------%
    F = A1-K1*C
    K = K1 + F*H

    F =                     K =

        -2  0   0   0           2   0   1   0
         0 -10  0   0           0   10  0   1
         0  0  -5   0          -2   1   4   0
        -1  2   0  -2           1  -2   0   0

The simulation results are considered next. The top level Simulink block
diagram is shown in Figure 11. The mass spring system dynamics are implemented
as shown in Figure 12. The UIO implementation is shown in Figure 13. The state
estimation error is shown in Figure 14. The following initial conditions were
used for the simulation.

    x1_0       = -1;
    x2_0       = -5;
    x3_0       = 1;
    x4_0       = 1;

Figure 14. The state estimation error for the mass spring system using the
UIO.

##  4\. Conclusion

This memo begins by motivating the UIO paradigm for state estimation. A main
advantage of this approach is to decouple the plant disturbance inputs from
the state estimation process. The UIO design procedure was reviewed and
illustrated with three common examples found in control theory. The estimation
error for the coupled mass spring system, the car suspension system and the
linear system with no inputs was considered. In each example, it was shown
that this error asymptotically approaches zero by utilization of a UIO.

## References

  * [1] Ron J. Patton Jie Chen, _Robust model-based fault diagnosis for dynamic systems_ , Springer, New York, 1999.
  * [2] R. Clark R. Patton, P. Frank, _Fault diagnosis in dynamic systems: Theory and applications_ , Prentice Hall, New York, 1989.
  * [3] Thoralf A. Schwarz, _Uncertainty Analysis of a Fault Detection and Isolation Scheme for Multi-Agent Systems_ , Master’s thesis, KTH, Sweden, 2012\.
  * [4] Bahram Shafai, _Proportional integral observer in robust control, fault detection and decentralized control of dynamic systems_ , To appear in book 21 (2015), 1087–1091.


# Probabilistic Symmetry for Improved Trajectory Forecasting

Sophia Sun  Robin Walters  Jinxi Li  Rose Yu

###### Abstract

Trajectory prediction is a core AI problem with broad applications in robotics
and autonomous driving. While most existing works focus on deterministic
prediction, producing probabilistic forecasts to quantify prediction
uncertainty is critical for downstream decision-making tasks such as risk
assessment, motion planning, and safety guarantees. We introduce a new metric,
mean regional score (MRS), to evaluate the quality of probabilistic trajectory
forecasts. We propose a novel probabilistic trajectory prediction model,
Probabilistic Equivariant Continuous COnvolution (PECCO) and show that
leveraging symmetry, specifically rotation equivariance, can improve the
predictions’ accuracy as well as coverage. On both vehicle and pedestrian
datasets, PECCO shows state-of-the-art prediction performance and improved
calibration compared to baselines.

Machine Learning, ICML

##  1 Introduction

Predicting the future trajectory of multiple agents is a critical task with
applications in autonomous driving (Chang et al., 2019), single-cell RNA
sequencing (Tong et al., 2020), and hurricane tracking (Hall and Jewson,
2007). The problem is difficult due to the strong dependency among the agents
where the number of interactions grows quadratically with the number of
agents; moreover, agent movements are influenced by the environmental features
such as road boundaries on the map. In practice, it is impossible to exactly
measure all dependencies among the agents and the influence from the
environment. Such a partially observed setting introduces a significant amount
of uncertainty.

Probabilistic modeling of trajectory dynamics provides a principled framework
to represent uncertainty. However, it remains elusive how to quantify and
evaluate the probabilistic forecasts. Consider the example of autonomous
driving. Although along a single-lane road with little traffic it may be
relatively certain drivers will continue at their current speed, it is
difficult to know with certainty which way a vehicle approaching an
intersection will turn purely from the past trajectories. By looking at map
information such as turn lanes, historical traffic flow data, and inferred
preferences, it is possible to place a wider or narrower probability
distribution over the potential future trajectories. To succeed, such a
probabilistic model must have good coverage and be well-calibrated; that is,
it must cover all the likely future scenarios without being so broad and
uncertain as to be useless.

(a) Original scene

(b) Rotated scene

Figure 1: Prediction on the same scene rotated by 90 degrees. PECCO is
consistent in trajectory and uncertainty prediction, whereas the non-
equivariant model fails.

Trajectory prediction has been investigated extensively in deep learning
(Liang et al., 2020b; Gao et al., 2020; Wang et al., 2018), mostly focusing on
deterministic point estimates. There have been some efforts in making
probabilistic forecasts for trajectories. For example, Tang and Salakhutdinov
(2019) introduce latent variables to the encoder-decoder framework, Lee et al.
(2017); Salzmann et al. (2020) propose conditional variational autoencoders
(CVAE), but little analysis of the quality of such probabilistic predictions
has been done.

Probabilistic symmetry plays an important role in the presence of data
scarcity (Bloem-Reddy and Teh, 2020). Our key insight is to exploit symmetry
to estimate high-dimensional conditional distributions with limited data. We
assume the predicted probability distribution is rotation and translation
equivariant. That is, if the input data is transformed, the probability
distribution will also be likewise transformed. In Figure 1, we see the same
car approaching an intersection from either the south or east. The scenes are
related by a $\pi/2$ rotation. As the absolute compass direction is not
particularly meaningful for local trajectory prediction, the model should thus
output the same probability distributions over future trajectories for the car
coming from the east as that coming from the north, but rotated by $\pi/2$.

In this paper, we propose a Probabilistic Equivariant Continuous COnvolutional
model, PECCO, an equivariant probabilistic trajectory prediction model. PECCO
uses continuous convolution (Wang et al., 2018) and enforces rotational
equivariance as in (Walters et al., 2021), achieving more accurate and
calibrated probablistic forecasts. Continuous convolution generalizes discrete
convolution and is useful for modeling interactions in many particle systems.
Rotational equivariance not only allows our model to produce physically
consistent predictions, the multiplicative nature of equivariance allows us to
model a density space with a smaller sample size. For each sample which an
equivariant model is trained on, an equivariant model learns as if it were
trained on all transformations of that sample by the symmetry group (Wang et
al., 2021).

Other methods for enforcing equivariance include data augmentation and
normalization. Data augmentation adds rotated versions of the samples to the
training dataset and allows the model to learn rotational equivariance.
However, this slows training drastically, requires greater model capacity, and
rarely achieves the level of equivariance or accuracy as equivariant neural
networks (Salzmann et al., 2020). To apply data normalization, we could rotate
the scene so the agent is pointing up in each case as in (Gao et al., 2020).
However, it is impossible to rotate the scene for multiple agents without a
canonical reference frame. PECCO allows the weights to be relative to the
local orientation of each agent without the need to rotate the scene
repeatedly.

Our main contributions are as follows:

  * •

We incorporate symmetry into probabilistic modeling of high-dimensional
sequence data.

  * •

We design an equivariant neural network, PECCO, to forecast multi-agent
trajectories with uncertainty.

  * •

We introduce a new metric, mean regional score (MRS), to evaluate the quality
of probabilistic forecasts.

  * •

We demonstrate that PECCO improves both accuracy and calibration of
probabilistic forecasts for two real-world benchmark datasets.

##  2 Related Work

#### Trajectory prediction.

Multi-agent trajectory forecasting has been extensively studied, approaches
ranged from Kalman filters (Kalman, 1960) to non-linear Gaussian Process
Regression models (Jordan, 1998). However, these methods either rely on strong
assumptions of the dynamics, or do not explicitly model multi-agent
interactions. We refer readers to (Rudenko et al., 2020) for a comprehensive
survey of such methods. Advancements in deep learning have allowed flexible
modeling of trajectory dynamics (Alahi et al., 2016; Lee et al., 2017; Wang et
al., 2018; Deo and Trivedi, 2018; Sadeghian et al., 2019; Liang et al., 2020a,
b; Walters et al., 2021; Gao et al., 2020; Roddenberry et al., 2021), but they
mostly focused on point estimation without uncertainty.

Recent methods have shifted to predicting distributions of future
trajectories, capturing the inherent uncertainty. There are two main classes
of models: (1) explicitly via exact likelihood (Tang and Salakhutdinov, 2019;
Chai et al., 2019; Gu et al., 2021) and via variational inference (Sohn et
al., 2015; Ivanovic and Pavone, 2019; Salzmann et al., 2020; Lee et al.,
2017), or (2) implicitly with Generative Adversarial Networks (GANs) (Gupta et
al., 2018; Liu et al., 2019). Our work falls into the first category where we
model the distributions parametrically. A parametric model allows us to
evaluate the likelihood of any future trajectories, which are useful for
downstream planning tasks (Chai et al., 2019; Schwarting et al., 2018).

Despite the development in probabilistic modeling, there is no standard metric
for quantifying the prediction uncertainty. Negative log likelihood often
overfits the distribution (Guo et al., 2017), and best-of-n-sample results do
not evaluate the full distribution (Ivanovic and Pavone, 2021). We argue that
probabilistic forecasts should accurately reflect the uncertainty in the model
predictions. We propose a nonparametric metric extending mean interval score
(Gneiting and Raftery, 2007). Furthermore, by leveraging the internal
symmetries of trajectories, our model produces more accurate and calibrated
probabilistic forecasts.

#### Equivariant Deep Learning.

Geometric deep learning that leverages invariance and symmetries has found
wide applications in areas ranging from image recognition (Bao and Song, 2019;
Worrall et al., 2017; Weiler and Cesa, 2019) to reinforcement learning (van
der Pol et al., 2020). Equivariant neural networks are studied for modeling
dynamics as well - Fuchs et al. (2020) use $\mathrm{SE}(3)$-equivariant
transformers to predict trajectories for a small number of particles as a
regression task, and (Walters et al., 2021) proposed a $\mathrm{S0}(2)$
equivariant continuous convolution for traffic trajectory prediction. All
methods mentioned above are deterministic. (Köhler et al., 2020) and (Satorras
et al., 2021) studies equivariant normalizing flows for modeling symmetric
densities, however their domains focus on generative modeling rather than
forecasting, and therefore differs from our work significantly. To our
knowledge, no previous work has studied equivariant neural networks for
probabilistic forecasting.

#### Uncertainty Quantification (UQ).

Uncertainty quantification is critical for risk assessment in safety-critical
domains. Properly quantified uncertainties can be used to generate
probabilistic constraints and whereby generate more robust planning and
control strategies (Ostafew et al., 2016; Bujarbaruah et al., 2019). With the
increasing use of deep learning in forecasting tasks, many works have UQ for
neural networks (Luo et al., 2021; Wu et al., 2021; Guo et al., 2017).
(Stankevičiūtė et al., 2021) proposes a conformal prediction algorithm for 1D
RNN forecasters. (Luo et al., 2021) explores using conformal methods to
produce collision warnings with safety guarantees. However, these works focus
only on classification or 1-dimensional forecasts. To our knowledge, PECCO is
the first paper to study quality of UQ for multivariate spatiotemporal
forecasts.

In our experiments, we show that the uncertainty distributions produced by
deep generative models for trajectory forecasting are often not calibrated.
Hence, better UQ methods are needed. In this work, we expand metrics from
quantile and conformal prediction to the multivariate setting, and introduce a
method to evaluate forecasted uncertainties. We also present a model design
that produces uncertainty distributions with better calibration.

##  3 Background

We give a short background on symmetry and equivariance and their
probabilistic extension.

#### Symmetry and Equivariance.

A _symmetry group_ $G$ is a set together with a composition operation
$\circ\colon G\times G\to G$ which is associative and has an identity and
inverses.

The group $G$ can transform a vector space $V$ by specifying a
_representation_ which is a mapping $\rho\colon G\to\mathrm{GL}_{n}(V)$
sending each element of the group $G$ to an invertible $n\times n$ matrix such
that $\rho(g_{1}\circ g_{2})=\rho(g_{1})\rho(g_{2})$.

Given a function $f\colon X\to Y$ such that $G$ has representations $\rho_{X}$
and $\rho_{Y}$ acting on $X$ and $Y$ respectively, we say $f$ is
$G$-_equivariant_ if for all $x\in X$ and $g\in G$, we have
$\rho_{Y}(g)f(x)=f(\rho_{X}(g)x)$. That is, a transformation of the input of
$x$ induces a corresponding transformation of the output. _Invariance_ for the
function $f$ is a special case in which $\rho_{Y}(g)y=y$.

####  $\mathbf{SO(2)}$ Rotational Symmetry.

In this paper, we consider equivariance to 2D planar rotations,
$\mathrm{SO}(2)=\\{\mathrm{Rot}_{\theta}:0\leq\theta<2\pi\\}$ which has
composition
$\mathrm{Rot}_{\theta_{1}}\circ\mathrm{Rot}_{\theta_{2}}=\mathrm{Rot}_{(\theta_{1}+\theta_{2})\mathrm{mod}2\pi}$.
As a Lie group, $\mathrm{SO}(2)$ also has an underlying topological structure
given by the manifold
$S^{1}\cong\\{{\mathbf{x}}\in\mathbb{R}^{2}:\|{\mathbf{x}}\|=1\\}$. Here, we
consider two specific representations of $\mathrm{SO}(2)$. The standard
representation $\rho_{1}$ on $\mathbb{R}^{2}$ is by $2\times 2$ rotation
matrices. The regular representation $\rho_{\mathrm{reg}}$ on
$L^{2}(\mathrm{SO}(2))=\\{\varphi:\mathrm{SO}(2)\to\mathbb{R}:|\varphi|^{2}\text{
is integrable}\\}$ is
$\rho_{\mathrm{reg}}(\mathrm{Rot}_{\phi})(\varphi)=\varphi\circ\mathrm{Rot}_{-\phi}$.
Here $L^{2}$ is a Hilbert space.

####  $\mathbf{SO(2)}$ Equivariant Continuous Convolution

Continuous convolution is the generalization of the discrete convolution for
point clouds. The feature vector
$\mathbf{f}^{(i)}\in\mathbb{R}^{c_{\mathrm{in}}}$ of particle $i$ forms a
vector field $\mathbf{f}$, and the kernel of the convolution
$K:\mathbb{R}^{2}\mapsto\mathbb{R}^{c_{\mathrm{out}}\times c_{\mathrm{in}}}$
forms a matrix field: for each point $\mathbf{x}\in\mathbb{R}^{2}$,
$K(\mathbf{x})$ is a $c_{\mathrm{out}}\times c_{\mathrm{in}}$ matrix. The
continuous convolution is then be defined by

$\mathbf{g}^{(i)}=\sum_{j}a(\|\mathbf{x}^{(j)}-\mathbf{x}^{(i)}\|)K(\mathbf{x}^{(j)}-\mathbf{x}^{(i)})\cdot\mathbf{f}^{(j)}.$

By (Weiler and Cesa, 2019), this is $\mathrm{SO}(2)$-equivariant if

$K(gv)=\rho_{\mathrm{out}}(g)K(v)\rho_{\mathrm{in}}(g^{-1}).$

ECCO (Walters et al., 2021) defines the convolution kernel $K$ in polar
coordinates $K(\theta,r)$. Let $\mathbb{R}^{c_{\mathrm{in}}}$ and
$\mathbb{R}^{c_{\mathrm{out}}}$ be $\mathrm{SO}(2)$-representations
$\rho_{\mathrm{in}}$ and $\rho_{\mathrm{out}}$ respectively, then the
convolution kernel satisfies the equivariance condition as

$K(\theta+\phi,r)=\rho_{\mathrm{out}}(\mathrm{Rot}_{\theta})K(\phi,r)\rho_{\mathrm{in}}(\mathrm{Rot}_{\theta}^{-1}),$

making the continuous convolution $\mathrm{SO}(2)$-equivariant.

#### Conditional Coverage.

Given an input variable $X$ and output variable $Y$, a prediction interval
with significance (error rate) $\alpha$ is valid if it satisfies

$\displaystyle P(Y\in\hat{C}^{\alpha}|X=x)\geq 1-\alpha$ (1)

which means the ground truth output ${Y}$ falls into the prediction interval
$\hat{C}^{\alpha}$ with at least $1-\alpha$ probability, conditioned on the
input. The left hand side of the inequality is known as conditional coverage,
which we will refer to as simply coverage in this paper. Theoretically, it is
difficult to ensure Equation (1) with statistical methods without further
assumptions on data (Barber et al., 2020). Nevertheless, we show that our
method satisfies Equation 1 in experiments.

##  4 Score Function for Probabilistic Forecast

We consider a probabilistic approach to trajectory forecasting. Ideally, we
want the predicted distribution to be both sharp and valid under the
definition in Equation 1. To quantify the uncertainty, we review the classic
score function Mean Interval Score (MIS) and derive a 2 dimensional extension,
Mean Regional Score (MRS).

#### Mean Interval Score.

Mean Interval Score (MIS) is a proper scoring rule (Gneiting and Raftery,
2007) for interval forecasts that rewards narrower confidence intervals and
encourages coverage. Specifically, let $Y\sim{P}_{Y}$ be a one-dimensional
random variable, and an upper bound $u$, a lower bound $l$ be its estimated
$\frac{\alpha}{2}$ and $(1-\frac{\alpha}{2})$ quantiles respectively, MIS is
defined using samples $y_{i}\sim{P}_{Y}$.

$\displaystyle\mathrm{MIS}(u,l;Y)=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}[(u-l)+\frac{2}{\alpha}(y_{i}-u)\mathbb{1}\\{y_{i}>u\\}$
$\displaystyle+$
$\displaystyle\frac{2}{\alpha}(l-y_{i})\mathbb{1}\\{y_{i}<l\\}]$ (2)

Where the three terms in Equation 2 accounts for the width of the interval,
and how much the sample point exceeds the upper bound and the lower bound,
respectively. MIS has the advantage of being easy to compute and does not
require the model to be parametric (Gneiting and Raftery, 2007). It can also
be used as an objective function for probabilistic forecasting (Wu et al.,
2021).

#### Mean Regional Score.

For n-D distributions, we introduce a new metric called the Mean Region Score
(MRS), which generalizes MIS to higher dimensions. Let $Z\sim{P}_{Z}$ be an
n-dimensional continuous random variable. First, we need to generalize
$\alpha$-quantile and define the confidence region. Let the region bordering
$z$ as $R(z)$, the confidence region of level $\alpha$

$\displaystyle R_{P}(\alpha)=\inf\\{R(z)|\int_{R(z)}P(z)dz\geq\alpha\\}$ (3)

is the smallest region whose probability exceeds level $\alpha$.

When the distribution is Gaussian $P_{Z}=\mathcal{N}(\mu,\Sigma)$ with mean
$\mu\in\mathbb{R}^{n}$ and the covariance matrix $\Sigma\in\mathbb{R}^{2\times
2}$ that is positive definite. The $(1-\alpha)$ confidence region of an n-D
multivariate normal distribution $P=\mathcal{N}(\mu,\Sigma)$ is an ellipsoid
that can be written as (Slotani, 1964):

$R(1-\alpha)=\\{z|(z-\mu)^{T}\Sigma^{-1}(z-\mu)\leq\chi_{n}^{2}(1-\alpha)\\}.$

where $\chi_{n}^{2}$ is the chi-squared distribution with n degrees of
freedom. Then, for a given sample $z_{i}\in\mathbb{R}^{n}$, we can draw an
ellipsoid whose edge coincides $z_{i}$ defined as

$R(z_{i})=\\{z|(z-\mu)^{T}\Sigma^{-1}(z-\mu)\leq c^{\prime}\\}$

where $c^{\prime}=(z_{i}-\mu)^{T}\Sigma^{-1}(z_{i}-\mu)$.

The MRS score of the $(1-\alpha)$ prediction interval thus can be evaluated
as:

$\displaystyle\mathrm{MRS}(R;Z)$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}[|R(1-\alpha)|+$
$\displaystyle\frac{1}{\alpha}|R(z_{i})\backslash
R(1-\alpha)|\mathbb{1}\\{z_{i}\in R^{C}(1-\alpha)\\}]$ (4)

where the first term corresponds to the area of the confidence region and the
second term measures how much each data $z_{i}$ deviates from the region.
$|\cdot|$ measures the size of the set. It is also easy to prove (see Appendix
A.4) that MRS is a proper scoring rule (Gneiting and Raftery, 2007), which
means MRS is optimized if and only if the distribution $P$ recovers that of
$Z$.

Figure 2 illustrates MRS for both unimodal and multimodal distributions. Note
that when the distribution is single-modal, the confidence interval in
Equation 4 is both equally tailed and the shortest confidence interval. For
multi-modal distributions, we solve for the smallest confidence region to
cover probable events. In appendix C.1 we present a method to estimate MRS for
arbitrary 2-D distributions, where we divide the domain space into small grids
to estimate density, and obtain $R(1-\alpha)$ and $R(z_{i})$ by assembling
high density regions.

Figure 2: MRS illustration for unimodal and multimodal distributions. Given a
sample $z_{i}$, MRS calculates the area of blue region plus that of the the
orange region scaled by $1/\alpha$.

##  5 PECCO

Figure 3: Overview of PECCO’s model architecture. Agent trajectories
consisting of velocities and position uncertainties are encoded along with map
information by equivariant continuous convolution and fully connected layers.
The model outputs $v_{j+1}$ and $\Sigma_{v,j+1}$ for all agents, which we use
to calculate position uncertainty $\Sigma_{x,j+1}$ via dynamics. The model
takes in the forecast and predicts autoregressively.

###  5.1 Problem Setup

Given past trajectories of $n$ agents over $t$ timesteps
${x}^{(i)}_{j}\in\mathbb{R}^{2}$ where $1\leq i\leq n$ and $1\leq j\leq t$,
and the environmental context information $\mathbf{e}$ including marker
positions of map lane boundaries, we model the probability distribution of
agents’ positions over $k$ future time steps as
$p_{\theta}({x}_{t+1:t+k}|{x}_{1:t},\mathbf{e})$, with
${x}_{j}=({x}^{(1)}_{j},\cdots,{x}^{(n)}_{j})$ being the positions of all
agents. We introduce PECCO, a deep learning model that leverages rotational
equivariance to produce probabilistic forecasts.

The high-level architecture of our model is illustrated in Figure 3. PECCO
takes as input the positions of all agents ${x}_{1:t}$ in the past, and the
covariance matrix $\Sigma_{x,j}$ at time $j$. It outputs the probability
distribution of each agents’ velocity as a 2-D Gaussian
${\mathcal{N}}(\mu_{v,j+1}^{(i)}$, $\Sigma_{v,j+1}^{(i)})$ for the next time
step. The velocity distribution is then integrated into a position
distribution $\mathcal{N}$($\mu_{x,j+1}^{(i)}$, $\Sigma_{x,j+1}^{(i)}$)
through dynamics integration. PECCO predicts the future $k$ timesteps
autoregressively. The output distributions are guaranteed to be equivariant
under rotation.

###  5.2 Probabilistic Symmetry through Equivariant Neural Networks

Rotational equivariance effectively reduces the dimension of the data space by
placing different samples in the same equivalence class. This improves data
coverage for better probabilistic modeling.

Intuitively, real-world trajectory dynamics have intrinsic symmetry. That is,
if the past trajectories and environmental data such as the lane and map are
rotated, then the probability of a rotated future trajectory will be equally
likely. We can model the probability density function $p_{\theta}$ as an
invariant function of its inputs,

$p_{\theta}({x}_{t+1:t+k}|{x}_{1:t},\mathbf{e})=p_{\theta}(g{x}_{t+1:t+k}|g{x}_{1:t},g\mathbf{e}).$
(5)

Here, past and future positions $x^{(i)}_{j}\in\mathbb{R}^{2}$ are transformed
according the standard representation $\rho_{1}$.

In order to implement Equation 5, we assume future positions follow a
multivariate normal distribution ${x}^{(i)}_{j}\sim
N(\mu^{(i)}_{j},\Sigma^{(i)}_{j})$. This is a common assumption in trajectory
forecasting literature (Rudenko et al., 2020) and provides a parametric closed
form for computing confidence regions and the score function.

We aim to construct an equivariant neural network $f_{\theta}$ which outputs
the parameters $\mu^{(i)}_{j}$ and $\Sigma^{(i)}_{j}$ autoregressively, taking
as input probability distributions over positions of all agents in the past
$k$ timesteps

$\mu_{j+1},\Sigma_{j+1}=f_{\theta}(\mu_{j-t+1:j},\Sigma_{j-t+1:j},\mathbf{e}).$

where $\mu_{j}=(\mu_{j}^{(1)},\ldots,\mu_{j}^{(n)})$ and
$\Sigma_{j}=(\Sigma_{j}^{(1)},\ldots,\Sigma_{j}^{(n)})$ and $\mathbf{e}$
denotes environmental information.

In this case, the equivariance of $f_{\theta}$ leads to the desired invariance
of $p_{\theta}$. This may be seen as a partial evaluation or currying of the
conditional probability density function which has the effect of transforming
invariance to equivariance. The following proposition relates equivariant
networks with probabilistic symmetry. See Appendix A for a proof.

######  Proposition 1.

If the one-step probabilistic forecasting model $f_{\theta}$ is
$G$-equivariant, then the probability distribution
$p_{\theta}(x_{t+1:t+k}|x_{1:t},\mathbf{e})$ is invariant as in Equation 5.

#### Equivariant Conditional Normal Distributions.

In order to enforce $\mathrm{SO}(2)$-equivariance for $f_{\theta}$, the
following proposition describes how the mean and covariance matrix for a 2-D
Gaussian transforms under a rotation of $\mathbb{R}^{2}$.

######  Proposition 2.

Given multivariate normal distribution $\mathcal{N}(\mu,\Sigma)$ over
$\mathbb{R}^{2}$ with probability density function $p_{\mu,\Sigma}$ and
$g\in\mathrm{SO}(2)$, then $\mathcal{N}(g\mu,g\Sigma g^{T})$ is also a
multivariate normal distribution and $p_{g\mu,g\Sigma
g^{T}}(v)=p_{\mu,\Sigma}(g^{-1}v)$ for all $v\in\mathbb{R}^{2}$.

To guarantee the covariance matrix of $f_{\theta}$ is positive-definite and
symmetric, i.e. $\Sigma\in\mathrm{PosDefSym}_{2}(\mathbb{R})$, we make use of
the following fact.

######  Proposition 3.

The map

$\displaystyle\varphi\colon\mathrm{GL}_{2}(\mathbb{R})$
$\displaystyle\to\mathrm{PosDefSym}_{2}(\mathbb{R})$ $\displaystyle M$
$\displaystyle\mapsto MM^{T}$

is surjective and equivariant. For $g\in\mathrm{SO}(2)$,

$\displaystyle\varphi(gM)$ $\displaystyle=g\varphi(M)g^{T}.$ (6)

Moreover, $\varphi$ admits a one-sided inverse which is also equivariant,

$\displaystyle\psi\colon\mathrm{PosDefSym}_{2}(\mathbb{R})$
$\displaystyle\to\mathrm{GL}_{2}(\mathbb{R})$ (7) $\displaystyle\Sigma$
$\displaystyle\mapsto Q\Lambda^{\frac{1}{2}}$ (8)

where $Q\Lambda Q^{T}$ is the eigendecomposition of $\Sigma$ and $Q$ is
orthogonal. That is, $\varphi(\psi(\Sigma))=\Sigma$.

As a consequence of Proposition 2 and Proposition 3, we can ensure that
$f_{\theta}$ outputs a normal distribution which transforms correctly under
rotations by setting

$\displaystyle\mu_{j+1},M_{j+1}$
$\displaystyle=\tilde{f}_{\theta}(\mu_{j-k:j},M_{j-k:j},\mathbf{e})$ (9)
$\displaystyle\Sigma_{j+1}$ $\displaystyle=M_{j+1}M_{j+1}^{T}$

and by constraining $\tilde{f}_{\theta}$ to be equivariant with respect to the
action in Equation 6. In this case, $\mathrm{SO}(2)$ acts by transforming the
the columns of $M$ independently as vectors in $\mathbb{R}^{2}$. Thus the data
$(\mu_{j}^{(i)},\Sigma_{j}^{(i)})$ for each agent and time step is comprised
of 3 copies of the standard representation $\rho_{1}$ as defined in Section 3.
Given this $\mathrm{SO}(2)$-action, we can enforce
$\mathrm{SO}(2)$-equivariance in the neural network $\tilde{f}_{\theta}$.
Details of implementation of this layer is provided in Appendix B.

###  5.3 Dynamics Integration (dyna)

Instead of predicting the position directly, PECCO outputs a Gaussian
distribution over the _velocity_ as $\mathcal{N}(\mu_{v,j},\Sigma_{v,j})$.
More specifically, it predicts $(\mu_{v,j},M_{v,j})$ at each time step and the
covariance matrix $\Sigma_{v,j}$ is calculated as in Proposition 3. However,
we want to obtain the uncertainty over position as $\Sigma_{x,j}$ and perform
autoregressive forecasting. We leverage dynamics integration to propagate the
uncertainty from velocity to position.

Assuming that all agents in the system can be approximated as linear discrete
time dynamics ${x}_{j+1}={x}_{j}+\Delta t\cdot{v}_{j}$, we can obtain the
uncertainty of predicted position $\Sigma_{x,j+1}$ by

$\Sigma_{x,j+1}=\Sigma_{x,j}+(\Delta t)^{2}\,\Sigma_{v,j}+2\Delta
t\cdot\text{cov}(x_{j},v_{j}).$

As in previous works (Salzmann et al., 2020), we make the assumption that the
cross covariance matrix $\mathrm{cov}(x_{j},v_{j})$ is zero for simplicity of
implementation. A consequence of such setup is that it enforces that the
uncertainty grows monotonically over time. During training, the gradients are
calculated after the entire trajectory is predicted.

##  6 Experiments

Model | minADE${}_{6}\downarrow$ | minFDE${}_{6}\downarrow$ | NLL$\downarrow$ | MRS $\downarrow$ | Cov@1s(%) | Cov@2s | Cov@3s
---|---|---|---|---|---|---|---
Argoverse
LSTM-NLL | 1.64 $\pm$ .05 | 4.17 $\pm$ .10 | 3.07 $\pm$ .08 | 19.1 $\pm$ 1.5 | 8.8 $\pm$ 0.7 | 8.5 $\pm$ 0.7 | 7.0 $\pm$ 0.8
LSTM-MRS | 1.86 $\pm$ .04 | 4.71 $\pm$ .08 | 12096\. $\pm$ 305. | 7.4 $\pm$ 0.3 | 99.4 $\pm$ 0.2 | 99.8 $\pm$ 0.1 | 99.8 $\pm$ 0.1
LSTM-NLL-r | 1.61 $\pm$ .02 | 4.15 $\pm$ .08 |  2.48 $\pm$ .03 | 20.3 $\pm$ 2.9 | 10.1 $\pm$ 1.5 | 10.5 $\pm$ 1.0 | 9.8 $\pm$ 1.9
LSTM-MRS-r | 1.66 $\pm$ .04 | 4.22 $\pm$ .08 | 37.26 $\pm$ 2.31 |  7.3 $\pm$ 0.2 |  86.4 $\pm$ 5.1 |  88.5 $\pm$ 4.7 |  91.6 $\pm$ 5.4
CtsConv-NLL | 1.74 $\pm$ .03 | 4.43 $\pm$ .06 | 49.13 $\pm$ 2.16 | 1409\. $\pm$ 277. | 6.3 $\pm$ 2.2 | 0.02 $\pm$ .01 | 0.01 $\pm$ .01
CtsConv-MRS | 1.76 $\pm$ .02 | 4.45 $\pm$ .05 | 16.80 $\pm$ .51 | 114\. $\pm$ 10. | 65.4 $\pm$ 3.5 | 11.4 $\pm$ 1.2 | 0.21 $\pm$ .10
CtsConv-NLL-aug | 1.66 $\pm$ .02 | 4.23 $\pm$ .06 | 11.81 $\pm$ .01 | 857\. $\pm$ 126. | 22.3 $\pm$ 2.1 | 1.7 $\pm$ 0.5 | 0.02 $\pm$ .01
Trajectron++ | 1.83 $\pm$ .02 | 3.85 $\pm$ .07 | 3.15 $\pm$ .17 | 55.2 $\pm$ 2.5 | 45.5 $\pm$ 5.3 | 37.6 $\pm$ 3.2 | 34.9 $\pm$ 2.5
MFP | 1.53 $\pm$ .04 | 3.77 $\pm$ .06 | 3.56 $\pm$ .02 | 12.5 $\pm$ 1.2 | 81.2 $\pm$ 5.1 | 53.0 $\pm$ 4.9 | 21.3 $\pm$ 6.8
PECCO |  1.39 $\pm$ .02 |  3.41 $\pm$ .03 | 4.26 $\pm$ 0.1 | 8.36 $\pm$ 0.16 | 74.9 $\pm$ 0.2 | 93.6 $\pm$ 0.8 | 94.5 $\pm$ 0.9
TrajNet++
LSTM-NLL-r | 0.85 $\pm$ .02 | 1.64 $\pm$ .03 | 2.78 $\pm$ .02 | 4.5 $\pm$ 0.6 | 21.9 $\pm$ 4.3 | 23.2 $\pm$ 4.2 | 23.7$\pm$ 3.9
LSTM-MRS-r | 0.94 $\pm$ .01 | 1.84 $\pm$ .02 | 123.7 $\pm$ 5.0 |  3.3 $\pm$ 0.1 | 95.7 $\pm$ 5.2 |  94.3 $\pm$ 5.6 | 96.5 $\pm$ 3.4
CtsCov-NLL | 1.08 $\pm$ .02 | 2.36 $\pm$ .09 | 5.33 $\pm$ .08 | 33.1 $\pm$ 3.2 | 73.7 $\pm$ 20.6 | 20.7 $\pm$ 5.2 | 15.2 $\pm$ 6.7
CtsConv-MRS | 0.96 $\pm$ .02 | 2.28 $\pm$ .07 | 5.97 $\pm$ .20 | 8.2 $\pm$ 0.9 |  93.4 $\pm$ 5.2 | 64.5 $\pm$ 5.2 | 46.1 $\pm$ 4.8
CtsCov-NLL-aug | 0.92 $\pm$ .01 | 1.76 $\pm$ .03 | 6.74 $\pm$ .11 | 20.0 $\pm$ 1.0 | 62.1 $\pm$ 3.3 | 36.3 $\pm$ 4.9 | 34.1 $\pm$ 5.8
Trajectron++ | 1.14 $\pm$ .03 | 2.31 $\pm$ .05 | 2.83 $\pm$ .12 | 45.1$\pm$ 6.4 | 50.2 $\pm$ 2.2 | 45.8 $\pm$ 3.5 | 32.9 $\pm$ 3.5
MFP | 0.85 $\pm$ .02 | 1.70 $\pm$ .04 |  2.20 $\pm$ .04 | 15.5 $\pm$ 0.6 | 79.1 $\pm$ 4.3 | 32.5 $\pm$ 3.1 | 22.8 $\pm$ 3.2
PECCO |  0.72 $\pm$ .02 |  1.41 $\pm$ .05 | 2.37 $\pm$ .04 | 6.1 $\pm$ 0.1 | 80.8 $\pm$ 4.5 | 81.9 $\pm$ 2.3 |  84.5 $\pm$ 3.0
Table 1:  Performance comparison on the benchmark Argoverse and TrajNet++
dataset. Bold-faced numbers indicate the best performance. Cov@1s(%) is the
percentile coverage at the 1 second prediction; prediction is more calibrated
if closer to 90%. We can see PECCO achieves superior regression accuracy while
maintaining good coverage.

We show that our model achieves state of the art prediction accuracy and
predicts a more calibrated confidence region compared to baseline models in
two trajectory prediction scenarios - autonomous vehicle and pedestrian.

Baselines. We compare our model PECCO against the following baselines:

  * •

LSTM: An encoder-decoder LSTM model for deterministic prediction optimizing
mean square error. We also trained a version with cannonicalized data LSTM-r.

  * •

CtsConv (Ummenhofer et al., 2019): Continuous convolution over point cloud
data for trajectory prediction, a non-equivariant counterpart of ECCO.
CtsConv-aug is trained with an additional data augmentation step where we
randomly rotate the scenes.

  * •

ECCO (Walters et al., 2021): State-of-the-art equivariant model for
deterministic trajectory prediction.

  * •

Multiple Futures Prediction (MFP) (Tang and Salakhutdinov, 2019): A encoder-
decoder model for multimodal probabilistic forecasts.

  * •

Trajectron++ (Salzmann et al., 2020): State-of-the-art probabilistic
trajectory prediction model with conditional VAE architecture.

For the deterministic models (LSTM, CtsConv), we modify them by training with
Negative Log-Likelihood (NLL) and Mean Regional Score (MRS) as loss functions,
the results for the deterministic baselines are in Appendix D.2.

Evaluation Metrics. We include the following metrics to evaluate both point
estimate performance, as well as the quality of the predicted distributions.

  * •

Minimum Average Displacement Error (minADE6): average $l_{2}$ displacement
error over $k$ steps between the prediction and ground truth. We report ADE
for deterministic models and the minimum ADE over 6 samples for probabilistic
models.

  * •

Minimum Final Displacement Error (minFDE6): $l_{2}$ displacement error between
the predicted and the ground truth final position. We report FDE for
deterministic models and the minFDE over 6 samples for probabilistic models.

  * •

Negative Log Likelihood (NLL): Mean NLL of the ground truth trajectory under
the output distribution.

  * •

Mean Regional Score (MRS): Score function for confidence region prediction
introduced in Equation 4. We report MRS for $\alpha=0.1$, the 90% cofidence
region.

  * •

Coverage: The empirical estimate of probability of the true value lying in the
predicted interval, defined in equation 1. We report the coverage of 90%
quantile of the predicted Gaussian, i.e. $\alpha=0.1$. The prediction is more
calibrated if the coverage is closer to 90%.

Datasets. The Argoverse autonomous vehicle motion forecasting (Chang et al.,
2019) is a widely used vehicle trajectory prediction benchmark. The task is to
predict 3 second trajectories based on all vehicle history in the past 2
seconds recorded at 10Hz. TrajNet++ (Sadeghian et al., 2018) is a popular
pedestrian trajectory benchmark with a focus on agent-agent scenarios. The
task is to predict 12 time steps for all agents given 9. For data and training
details see D.1.

We report the performance on the official validation set for Argoverse, and on
a 10% test split for TrajNet++. PECCO predicts all agents’ trajectories
jointly. Following literature, we report on metrics for the agent vehicle
(Argoverse) and primary pedestrian (TrajNet++).

Prediction Accuracy. Table 1 show the prediction performance on two benchmark
datasets. PECCO achieves better regression accuracy compared to non-
equivariant baseline in terms of minADE and minFDE, with a notable 9%
improvement in minADE over the the best performing baseline, MFP. PECCO’s
improved probability coverage allows for more diverse sampling and hence can
produce trajectories closer to ground truth. For TrajNet++, PECCO outperforms
baseline models in both prediction accuracy and calibration.

(a) LSTM-MRS-r

(b) CtsConv-NLL

(c) CtsConv-MRS

(d) PECCO

Figure 4: Comparison of uncertainty predicted at a lane change. Red trajectories are the agent of interest in the same scenario. Note how the LSTM predicted uncertainty explodes after a few time steps, and PECCO is able to model possibility of both staying and lane change.  Model | minADE${}_{6}\downarrow$ | minFDE${}_{6}\downarrow$ | MRS $\downarrow$ | Cov@1s(%) | Cov@2s | Cov@3s
---|---|---|---|---|---|---
Conformal LSTM | 2.45 $\pm$ 0.09 | 4.68 $\pm$ 0.15 | 198.1 $\pm$ 12.0 | 90.1 $\pm$ 0.1 | 92.7 $\pm$ 0.1 | 92.8 $\pm$ 0.1
Conformal ECCO | 1.96 $\pm$ 0.06 | 4.32 $\pm$ 0.10 | 220.9 $\pm$ 8.1 |  90.0 $\pm$ 0.1 |  90.1 $\pm$ 0.1 |  90.0 $\pm$ 0.1
PECCO |  1.39 $\pm$ .02 |  3.41 $\pm$ 0.03 |  8.52 $\pm$ 0.16 | 74.9 $\pm$ 0.2 | 93.6 $\pm$ 0.8 | 92.5 $\pm$ 0.9
Table 2: Comparison with conformal prediction methods on Argoverse dataset.
PECCO produces a parametric distribution with a tighter confidence region
(small MRS), whereby achieving better regression accuracy while maintaining
competitive coverage.

Uncertainty Evaluation. PECCO’s uncertainty estimates are able to achieve
consistently good coverage compared to other methods whose coverage
deteriorates over the prediction horizon. By comparing LSTM-NLL and LSTM-MIS
with their augmented counterparts LSTM-NLL-aug and LSTM-MRS-aug, we can see
that data augmentation through rotation improves both accuracy and
calibration. Models trained with MRS loss (LSTM-MRS, LSTM-MRS-aug, CtsConv-
MRS) can consistently achieve high coverage, but their predictions suffer from
the wide distribution as reflected by their worse NLL and minADE/FDE numbers.
On the other hand, even though trained using NLL loss, PECCO is able to
maintain good calibration while achieving improved accuracy by leveraging
equivariance.

Figure 4 visualizes a typical situation to illustrate this difference. We plot
the predicted distribution at 10, 20, and 30 timesteps of prediction (1
timestep is 0.1 seconds). We can see the probable region predicted by LSTM-
MRS-r explodes at timesteps 20 and 30, whereas CtsConv-NLL and CtsConv-MIS
tend to be overconfident in their predictions. PECCO is able to predict a
Gaussian that covers the cases for both staying in the lane and changing to
the right lane.

###  6.1 Model Analysis

Comparison with data augmentation and cannonicalization. Data augmentation and
cannonicalization are popular methods to implicitly exploit symmetry in
trajectory data. In Table 3 we compare PECCO to its non-equivariant
counterparts with the addition of data augmentation and cannonicalization.
Equivariance dramatically improves the prediction performance in all aspects,
especially calibration.

| minA/FDE6 | NLL | MRS | Cov(%)
---|---|---|---|---
CtsConv | 1.74 / 4.43 | 49.13 | 1498.4 | 2.2
+cannon | 1.66 / 4.28 | 47.46 | 1246.9 | 4.5
+aug | 1.67 / 4.23 | 11.81 | 857.0 | 17.1
equivariant |  1.39 / 3.41 | 4.26 | 8.52 | 87.5
Table 3: Data Augmentation Comparison on the Argoverse Dataset

Dynamics Integration Ablation. Dynamic integration (dyna) introduced in
Section 5.3 enforces the uncertainty to grow monotonically over time. As an
ablative study, Table 4 show that PECCO with dynamic integration has much
better calibration compared to if without; for the Argoverse dataset, the
improved calibration also informs better performance on minADE/FDE.

Argoverse | minA/FDE6 | NLL | MRS | Cov(%)
---|---|---|---|---
no-dyna | 1.52 / 3.76 | 9.72 | 39.1 | 38.6
dyna |  1.39 / 3.41 | 4.26 | 8.36 | 87.5
Pedestrian | minA/FDE6 | NLL | MRS | Cov(%)
no-dyna |  0.72 / 2.12 | 4.71 | 24.6 | 39.6
dyna | 0.73 / 1.98 | 2.37 | 6.12 | 83.7
Table 4: Dynamics Integration (dyna) Ablation: using dyna encourages the
uncertainty to grow over time and improves coverage.

Comparison with Conformal Prediction. Conformal prediction is an alternative
uncertainty quantification method with theoretical guarantees. It does not
impose any assumptions about the underlying distribution. We implemented two
conformal prediction baselines based on (Stankevičiūtė et al., 2021) (Appendix
C.2). We use half of the dataset for training and calibration respectively,
and show performance in Table 2. Conformal methods achieve guaranteed $\geq
90\%$ coverage, but suffer in prediction accuracy due to less training data.
The predicted regions are also much larger (higher MRS), which are less
desirable in downstream tasks.

##  7 Conclusion

We propose Probabilistic Equivariant Continuous Convolution (PECCO), a novel
method for improving uncertainty quantification for trajectory forecasting. We
design an equivariant neural network under which the predicted distributions
transform correspondingly as inputs are transformed. We introduce the MRS
metric to bring attention to the calibration of probabilistic forecasts. By
leveraging equivariance, PECCO produces more accurate and calibrated
probabilistic forecasts compared to existing methods. A future direction of
this work is to expand PECCO to produce multimodal forecasts. In this paper,
we made the simplifying assumption that the uncertainties are 2-D normal
distributions. By modifying PECCO to output a mixture of Gaussians, we can
describe more complex trajectory distributions.

## References

  * Alahi et al. [2016] Alexandre Alahi, Kratarth Goel, Vignesh Ramanathan, Alexandre Robicquet, Li Fei-Fei, and Silvio Savarese.  Social lstm: Human trajectory prediction in crowded spaces.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016.
  * Bao and Song [2019] Erkao Bao and Linqi Song.  Equivariant neural networks and equivarification.  _arXiv preprint arXiv:1906.07172_ , 2019.
  * Barber et al. [2020] Rina Foygel Barber, Emmanuel J. Candès, Aaditya Ramdas, and Ryan J. Tibshirani.  The limits of distribution-free conditional predictive inference, 2020\.
  * Bloem-Reddy and Teh [2020] Benjamin Bloem-Reddy and Yee Whye Teh.  Probabilistic symmetries and invariant neural networks.  _Journal of Machine Learning Research_ , 21(90), 2020.
  * Bujarbaruah et al. [2019] Monimoy Bujarbaruah, Xiaojing Zhang, Marko Tanaskovic, and Francesco Borrelli.  Adaptive mpc under time varying uncertainty: Robust and stochastic.  _arXiv preprint arXiv:1909.13473_ , 2019.
  * Chai et al. [2019] Yuning Chai, Benjamin Sapp, Mayank Bansal, and Dragomir Anguelov.  Multipath: Multiple probabilistic anchor trajectory hypotheses for behavior prediction.  _arXiv preprint arXiv:1910.05449_ , 2019.
  * Chang et al. [2019] Ming-Fang Chang, John W Lambert, Patsorn Sangkloy, Jagjeet Singh, Slawomir Bak, Andrew Hartnett, De Wang, Peter Carr, Simon Lucey, Deva Ramanan, and James Hays.  Argoverse: 3d tracking and forecasting with rich maps.  In _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2019.
  * Deo and Trivedi [2018] Nachiket Deo and Mohan M Trivedi.  Convolutional social pooling for vehicle trajectory prediction.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_ , 2018.
  * Fuchs et al. [2020] Fabian B Fuchs, Daniel E Worrall, Volker Fischer, and Max Welling.  $\mathrm{SE}(3)$-transformers: 3D roto-translation equivariant attention networks.  _arXiv preprint arXiv:2006.10503_ , 2020.
  * Gao et al. [2020] Jiyang Gao, Chen Sun, Hang Zhao, Yi Shen, Dragomir Anguelov, Congcong Li, and Cordelia Schmid.  Vectornet: Encoding hd maps and agent dynamics from vectorized representation.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020.
  * Gneiting and Raftery [2007] Tilmann Gneiting and Adrian E Raftery.  Strictly proper scoring rules, prediction, and estimation.  _Journal of the American statistical Association_ , 102(477), 2007.
  * Gu et al. [2021] Junru Gu, Chen Sun, and Hang Zhao.  Densetnt: End-to-end trajectory prediction from dense goal sets.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2021.
  * Guo et al. [2017] Chuan Guo, Geoff Pleiss, Yu Sun, and Kilian Q Weinberger.  On calibration of modern neural networks.  In _International Conference on Machine Learning_. PMLR, 2017.
  * Gupta et al. [2018] Agrim Gupta, Justin Johnson, Li Fei-Fei, Silvio Savarese, and Alexandre Alahi.  Social gan: Socially acceptable trajectories with generative adversarial networks.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2018.
  * Hall and Jewson [2007] Timothy M Hall and Stephen Jewson.  Statistical modelling of north atlantic tropical cyclone tracks.  _Tellus A: Dynamic Meteorology and Oceanography_ , 59(4), 2007.
  * Ivanovic and Pavone [2019] Boris Ivanovic and Marco Pavone.  The trajectron: Probabilistic multi-agent trajectory modeling with dynamic spatiotemporal graphs.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2019.
  * Ivanovic and Pavone [2021] Boris Ivanovic and Marco Pavone.  Rethinking trajectory forecasting evaluation.  _arXiv preprint arXiv:2107.10297_ , 2021.
  * Jordan [1998] Michael I. Jordan, editor.  _Prediction with Gaussian Processes: From Linear Regression to Linear Prediction and Beyond_.  The MIT Press, 1998.
  * Kalman [1960] Rudolph Emil Kalman.  A new approach to linear filtering and prediction problems.  _Journal of Basic Engineering_ , 82(1), 1960.
  * Köhler et al. [2020] Jonas Köhler, Leon Klein, and Frank Noé.  Equivariant flows: exact likelihood generative learning for symmetric densities.  In _International Conference on Machine Learning_. PMLR, 2020.
  * Lee et al. [2017] Namhoon Lee, Wongun Choi, Paul Vernaza, Christopher B Choy, Philip HS Torr, and Manmohan Chandraker.  Desire: Distant future prediction in dynamic scenes with interacting agents.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2017.
  * Liang et al. [2020a] Ming Liang, Bin Yang, Rui Hu, Yun Chen, Renjie Liao, Song Feng, and Raquel Urtasun.  Learning lane graph representations for motion forecasting.  In _European Conference on Computer Vision_. Springer, 2020a.
  * Liang et al. [2020b] Ming Liang, Bin Yang, Rui Hu, Yun Chen, Renjie Liao, Song Feng, and Raquel Urtasun.  Learning lane graph representations for motion forecasting.  In _European Conference on Computer Vision_. Springer, 2020b.
  * Liu et al. [2019] Yukai Liu, Rose Yu, Stephan Zheng, Eric Zhan, and Yisong Yue.  Naomi: Non-autoregressive multiresolution sequence imputation.  _Advances in Neural Information Processing Systems_ , 2019.
  * Luo et al. [2021] Rachel Luo, Shengjia Zhao, Jonathan Kuck, Boris Ivanovic, Silvio Savarese, Edward Schmerling, and Marco Pavone.  Sample-efficient safety assurances using conformal prediction.  _arXiv preprint arXiv:2109.14082_ , 2021.
  * Ostafew et al. [2016] Chris J Ostafew, Angela P Schoellig, Timothy D Barfoot, and Jack Collier.  Learning-based nonlinear model predictive control to improve vision-based mobile robot path tracking.  _Journal of Field Robotics_ , 33(1), 2016.
  * Roddenberry et al. [2021] T Mitchell Roddenberry, Nicholas Glaze, and Santiago Segarra.  Principled simplicial neural networks for trajectory prediction.  In _International Conference on Machine Learning_. PMLR, 2021.
  * Rudenko et al. [2020] Andrey Rudenko, Luigi Palmieri, Michael Herman, Kris M Kitani, Dariu M Gavrila, and Kai O Arras.  Human motion trajectory prediction: a survey.  _The International Journal of Robotics Research_ , 39(8), 2020.
  * Sadeghian et al. [2018] Amir Sadeghian, Vineet Kosaraju, Agrim Gupta, Silvio Savarese, and Alexandre Alahi.  Trajnet: Towards a benchmark for human trajectory prediction.  _arXiv preprint_ , 2018.
  * Sadeghian et al. [2019] Amir Sadeghian, Vineet Kosaraju, Ali Sadeghian, Noriaki Hirose, Hamid Rezatofighi, and Silvio 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.
  * Salzmann et al. [2020] Tim Salzmann, Boris Ivanovic, Punarjay Chakravarty, and Marco Pavone.  Trajectron++: Dynamically-feasible trajectory forecasting with heterogeneous data.  In _ECCV_. Springer, 2020.
  * Satorras et al. [2021] Victor Garcia Satorras, Emiel Hoogeboom, Fabian Bernd Fuchs, Ingmar Posner, and Max Welling.  E (n) equivariant normalizing flows.  In _Thirty-Fifth Conference on Neural Information Processing Systems_ , 2021.
  * Schwarting et al. [2018] Wilko Schwarting, Javier Alonso-Mora, and Daniela Rus.  Planning and decision-making for autonomous vehicles.  _Annual Review of Control, Robotics, and Autonomous Systems_ , 1, 2018\.
  * Shafer and Vovk [2008] Glenn Shafer and Vladimir Vovk.  A tutorial on conformal prediction.  _Journal of Machine Learning Research_ , 9(3), 2008.
  * Slotani [1964] Mlnoru Slotani.  Tolerance regions for a multivariate normal population.  _Annals of the Institute of Statistical Mathematics_ , 16(1), 1964.
  * Sohn et al. [2015] Kihyuk Sohn, Honglak Lee, and Xinchen Yan.  Learning structured output representation using deep conditional generative models.  _Advances in neural information processing systems_ , 28, 2015.
  * Stankevičiūtė et al. [2021] Kamilė Stankevičiūtė, Ahmed Alaa, and Mihaela van der Schaar.  Conformal time-series forecasting.  In _Advances in Neural Information Processing Systems_ , 2021.
  * Tang and Salakhutdinov [2019] Yichuan Charlie Tang and Ruslan Salakhutdinov.  Multiple futures prediction.  _Advances in Neural Information Processing Systems_ , 2019.
  * Tong et al. [2020] Alexander Tong, Jessie Huang, Guy Wolf, David Van Dijk, and Smita Krishnaswamy.  Trajectorynet: A dynamic optimal transport network for modeling cellular dynamics.  In _International Conference on Machine Learning_. PMLR, 2020.
  * Ummenhofer et al. [2019] Benjamin Ummenhofer, Lukas Prantl, Nils Thuerey, and Vladlen Koltun.  Lagrangian fluid simulation with continuous convolutions.  In _International Conference on Learning Representations_ , 2019.
  * van der Pol et al. [2020] Elise van der Pol, Daniel Worrall, Herke van Hoof, Frans Oliehoek, and Max Welling.  Mdp homomorphic networks: Group symmetries in reinforcement learning.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Walters et al. [2021] Robin Walters, Jinxi Li, and Rose Yu.  Trajectory prediction using equivariant continuous convolution.  _International Conference on Learning Representations_ , 2021.
  * Wang et al. [2021] Rui Wang, Robin Walters, and Rose Yu.  Incorporating symmetry into deep dynamics models for improved generalization.  _International Conference on Learning Representations_ , 2021.
  * Wang et al. [2018] Shenlong Wang, Simon Suo, Wei-Chiu Ma, Andrei Pokrovsky, and Raquel Urtasun.  Deep parametric continuous convolutional neural networks.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2018.
  * Weiler and Cesa [2019] Maurice Weiler and Gabriele Cesa.  General E(2)-equivariant steerable CNNs.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2019.
  * Worrall et al. [2017] Daniel E Worrall, Stephan J Garbin, Daniyar Turmukhambetov, and Gabriel J Brostow.  Harmonic networks: Deep translation and rotation equivariance.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2017.
  * Wu et al. [2021] Dongxia Wu, Liyao Gao, Xinyue Xiong, Matteo Chinazzi, Alessandro Vespignani, Yi-An Ma, and Rose Yu.  Quantifying uncertainty in deep spatiotemporal forecasting.  _ACM SIGKDD Conference on Knowledge Discovery and Data Mining (KDD)_ , 2021.

##  Appendix A Proofs

###  A.1 Proofs for Equivariance of Gaussians [Proposition 2]

In this section, we prove that the Gaussian probability density function is
$\mathrm{SO}(2)$-equivariant.

######  Proposition.

Given multivariate normal distribution $\mathcal{N}(\mu,\Sigma)$ over
$\mathbb{R}^{2}$ with probability density function $p_{\mu,\Sigma}$ and
$g\in\mathrm{SO}(2)$, then $p_{g\mu,g\Sigma g^{T}}(v)=p_{\mu,\Sigma}(g^{-1}v)$
for all $v\in\mathbb{R}^{2}$.

###### Proof.

$\displaystyle P_{g\mu,g\Sigma g^{T}}(v)$
$\displaystyle=\frac{1}{2\pi\det{(g\Sigma
g^{T})}}\exp{(-\frac{1}{2}(v-g\mu)^{T}(g\Sigma g^{T})^{-1}(v-g\mu))}$
$\displaystyle=\frac{1}{2\pi\det{(\Sigma})}\exp{(-\frac{1}{2}(g^{-1}v-\mu)^{T}g^{T}g^{-T}\Sigma^{-1}g^{-1}g(g^{-1}v-\mu))}$
$\displaystyle=\frac{1}{2\pi\det{(\Sigma})}\exp{(-\frac{1}{2}(g^{-1}v-\mu)^{T}\Sigma^{-1}(g^{-1}v-\mu))}$
$\displaystyle=P_{\mu,\Sigma}(g^{-1}v)$

∎

######  Proposition.

Given a normal distribution $\mathcal{N}(\mu,\Sigma)$, its group-actioned
$\mathcal{N}(g\mu,g\Sigma g^{T})$ is also a valid normal distribution.

###### Proof.

We prove the proposition by proving that $g\Sigma g^{T}$ is a symmetric
positive definite matrix.

1) (Symmetry) $(g\Sigma g^{T})^{T}=g^{TT}\Sigma^{T}g^{T}=g\Sigma g^{T}$.

2) (Positive Definite) Let $v\in\mathbb{R}^{2}$ and $v\neq 0$. Let $w=g^{T}v$,
note that $w\neq 0$. Then $v^{T}g\Sigma g^{T}v=w^{T}\Sigma w>0$, which implies
$\Sigma$ is a positive definite matrix. ∎

###  A.2 Equivariance of construct for positive definite and symmetric
$\Sigma$ [Proposition 3]

Note that every positive definite symmetric matrix can be written as
$\Sigma=MM^{T}$ for some $M\in GL_{2}(\mathbb{R})$, so we can represent
$\Sigma$ by a map $\varphi$ where

$\displaystyle\varphi\colon\mathrm{GL}_{2}(\mathbb{R})$
$\displaystyle\to\mathrm{PosDefSym}_{2}(\mathbb{R})$ $\displaystyle M$
$\displaystyle\mapsto MM^{T}$

######  Proposition.

$\varphi$ is $\mathrm{SO}(2)$-equivariant.

###### Proof.

$\varphi(gMg^{T})=gMg^{T}(gMg^{T})^{T}=gMg^{T}g^{TT}M^{T}g^{T}=gMM^{T}g^{t}=g\varphi(M)g^{T}$
∎

We defined the one-sided inverse of $\varphi$ as $\psi$:

$\displaystyle\psi\colon\mathrm{PosDefSym}_{2}(\mathbb{R})$
$\displaystyle\to\mathrm{GL}_{2}(\mathbb{R})$ $\displaystyle\Sigma$
$\displaystyle\mapsto Q\lambda^{\frac{1}{2}}$

where $Q\lambda Q^{T}$ is the eigendecomposition of $\Sigma$ and $Q$ is
orthogonal. That is, $\varphi(\psi(\Sigma))=\Sigma$.

######  Proposition.

$\psi$ is $\mathrm{SO}(2)$-equivariant.

###### Proof.

For symmetric positive definite matrix $\Sigma$, its eigenvector $v_{i}$ and
eigenvalue $\lambda_{i}$ follows the eigenequation

$\Sigma v_{i}=\lambda_{i}v_{i}.$

For $\forall g\in\mathrm{SO}(2)$,

$g\Sigma v_{i}=\lambda_{i}gv_{i}\implies g\Sigma
g^{-1}gv_{i}=\lambda_{i}gv_{i}\implies(g\Sigma
g^{-1})(gv_{i})=\lambda_{i}(gv_{i}),$

so $gv_{i}$ is an eigenvector to $g\Sigma g^{-1}$, with corresponding
eigenvalue $\lambda_{i}$.

Therefore,

$\psi(g\Sigma)=gQ\Lambda^{\frac{1}{2}}=g\psi(\Sigma).$

∎

###  A.3 Proof of Proposition 1: The Equivariant Model Gives an Invariant
Conditional Distribution.

We now prove Proposition 1.

######  Proposition (Proposition 1).

If the one-step probabilistic forecasting model $f_{\theta}$ is
$G$-equivariant, then the conditional probability distribution
$p_{\theta}(\mathbf{x}_{t+1:t+k}|\mathbf{x}_{1:t},\mathbf{e})$ is invariant as
in Equation 5.

###### Proof.

Since $\mathbf{x}_{1:t}$ are known, we define
$\Sigma_{1:t}=\epsilon\mathrm{Id}$ for $\epsilon$ small approximating a delta
distribution at $\mu_{1:t}=\mathbf{x}_{1:t}$. We then evaluate
$\mu_{j+1},\Sigma_{j+1}=f_{\theta}(\mu_{j-t+1:j},\Sigma_{j-t+1:j},\mathbf{e})$
autoregressively for $j=1,\ldots,k$. Finally define

$\displaystyle p_{\theta}(\mathbf{x}_{t+1:t+k}|\mathbf{x}_{1:t},\mathbf{e})$
$\displaystyle=\prod_{j=1}^{k}p_{\theta}(\mathbf{x}_{t+j}|\mathbf{x}_{j:t+j-1},\mathbf{e})$
(10)
$\displaystyle=\prod_{j=1}^{k}p_{N(\mu_{t+j},\Sigma_{t+j})}(\mathbf{x}_{t+j}).$

For $g\in\mathrm{SO}(2)$, we want to show

$p_{\theta}(g\mathbf{x}_{t+1:t+k}|g\mathbf{x}_{1:t},g\mathbf{e})=p_{\theta}(\mathbf{x}_{t+1:t+k}|\mathbf{x}_{1:t},\mathbf{e}).$

Thus for $g\mathbf{x}_{1:t},g\mathbf{e}$ we initialize
$g\mu_{1:t},g\Sigma_{1:t}g^{T}$ since
$g\Sigma_{1:t}g^{T}=gg^{T}\Sigma_{1:t}=\Sigma_{1:t}$. Applying $f_{\theta}$
repeatedly and invoking equivariance, we obtain

$\displaystyle
g\mu_{t+j},g\Sigma_{t+j}g^{T}=f_{\theta}(g\mu_{j:t+j-1},g\Sigma_{j:t+j-1}g^{T},g\mathbf{e}).$
(11)

for $j=1,\ldots,k$. Then by Equation 10,

$p_{\theta}(g\mathbf{x}_{t+1:t+k}|g\mathbf{x}_{1:t},g\mathbf{e})=\prod_{j=1}^{k}p_{N(g\mu_{t+j},g\Sigma_{t+j}g^{T})}(g\mathbf{x}_{t+j}).$

Then by Proposition 2 this is

$\prod_{j=1}^{k}p_{N(\mu_{t+j},\Sigma_{t+j})}(g^{-1}g\mathbf{x}_{t+j}),$

which cancels $g^{-1}g=1$, and applying Equation 10 again gives
$p_{\theta}(\mathbf{x}_{t+1:t+k}|\mathbf{x}_{1:t}\mathbf{e})$. ∎

###  A.4 Proof that MRS is a strictly proper scoring rule

Scoring rules are summary measures to evaluate probabilistic forecasts; they
take in the forecasted distribution and the event or value that materializes
and assign a numerical score. As defined in [Gneiting and Raftery, 2007], a
scoring rule is proper if it is maximized when the forecaster recovers the
ground truth distribution. It is strictly proper if the maximum is unique.

######  Proposition.

MRS is a strictly proper score.

###### Proof.

For arbitrary function $h$, we first define the scoring rule

$S(R;z)=-|R|+\frac{1}{\alpha}(|R|-|R(z)|)\mathbb{1}\\{z\in R^{C}\\}+h(z),$

which means if the forecaster quotes region $R$ at the level $\alpha\in(0,1)$
and $z$ materializes, then the score $S(R;z)$ will be rewarded, where $R$ and
$R(z)$ are defined as ellipsoid parametric by $\mu$ and $\Sigma$. Then the
expected score under the probability measure $P\in\mathcal{P}$ which $z$
follows is defined as

$S(R;P)=\int S(R;z)\,dP(z).$

For $P\in\mathcal{P}$, let $R_{*}$ denote the unique true $P$-region at level
$\alpha$. We say that scoring rule $S$ is strictly proper if

$S(R_{*};P)\geq S(R;P),$

for all region $R$ and for all probability measures $P\in\mathcal{P}$, where
equality holds if and only if $R=R_{*}$.

We identify $P$ with the associated distribution function so that
$P(R_{*}^{C})=\alpha$. If $R_{*}\subset R$, then

$\displaystyle S(R_{*};P)-S(R;P)$
$\displaystyle=\int-|R_{*}|\,dP(z)+\frac{1}{\alpha}\int_{R_{*}^{C}}(|R_{*}|-|R(z)|)\,dP(z)$
$\displaystyle+\int|R|\,dP(z)-\frac{1}{\alpha}\int_{R^{C}}(|R|-|R(z)|)\,dP(z)$
$\displaystyle=-|R_{*}|+|R|+\frac{1}{\alpha}|R_{*}|P(R_{*}^{C})-\frac{1}{\alpha}|R|P(R^{C})-\frac{1}{\alpha}\int_{R_{*}^{C}\setminus
R^{C}}|R(z)|\,dP(z)$
$\displaystyle\geq|R|-\frac{1}{\alpha}|R|P(R^{C})-\frac{1}{\alpha}|R|(P(R_{*}^{C})-P(R^{C}))$
$\displaystyle=0,$

as it is supposed to be. If $R\subset R_{*}$, then an analogous argument
applies.

Putting $h(z)=0$ and reversing the sign of the scoring rule, yields the
negatively oriented regional score. Finally, the mean score of this oriented
regional score is our MRS, which is then proved to be proper. ∎

##  Appendix B Equivariant Matrix Output Layer Construction

Let

$\rho_{1}(\theta)=\begin{pmatrix}\cos{\theta}&-\sin{\theta}\\\
\sin{\theta}&\cos{\theta}\end{pmatrix}\qquad\qquad M=\begin{pmatrix}a&b\\\
c&d\end{pmatrix}$

where $a,b,c,d$ are trainable parameters. Consider input $f\colon
S^{1}\to\mathbb{R}$, then

$F(f)=\int_{\theta=0}^{2\pi}f(\theta)\rho_{1}(\theta)M\rho_{1}(-\theta)d\theta.$

Define $\rho_{reg}(\phi)(f)(\theta)=f(\theta-\phi)$.

######  Proposition.

The convolution $F$ is equivariant, i.e.,

$F(\rho_{reg}(\phi)f)=\rho_{1}(\phi)F(f)\rho_{1}(-\phi).$

###### Proof.

We compute

$\displaystyle F(\rho_{reg}(\phi)f)$
$\displaystyle=\int_{\theta=0}^{2\pi}\rho_{reg}(\phi)(f)(\theta)\rho_{1}(\theta)M\rho_{1}(-\theta)d\theta$
$\displaystyle=\int_{\theta=0}^{2\pi}f(\theta-\phi)\rho_{1}(\theta)M\rho_{1}(-\theta)d\theta$
$\displaystyle=\int_{u=-\phi}^{2\pi-\phi}f(u)\rho_{1}(u+\phi)M\rho_{1}(-u-\phi)du\quad(\text{substituting
}u=\theta-\phi)$
$\displaystyle=\rho_{1}(\phi)\cdot\left(\int_{u=-\phi}^{2\pi-\phi}f(u)\rho_{1}(u)M\rho_{1}(-u)du\right)\cdot\rho_{1}(-\phi)$
$\displaystyle=\rho_{1}(\phi)F(f)\rho_{1}(-\phi).$

The last substitution follows from the fact the integrand is periodic of
period $2\pi$. Since the integral is over the whole circle, this equal to
taking the limits of integration to be $0$ and $2\pi$. ∎

#### Discritization.

If the function $f$ is discretized, such that $f_{i}=f(i2\pi/n)$. Then the
mapping $F$ is given by

$F(f)=\sum_{i=0}^{n-1}f_{i}\rho_{1}(i2\pi/n)M\rho_{1}(-i2\pi/n).$

#### Compatibilitiy with matrix exponential

For a $n\times n$ matrix $M$ define

$\exp(M)=\sum_{k=0}^{\infty}M^{k}/k!.$

Matrix exponentiation is equivariant with respect to conjugation. That is for
$g\in\mathrm{GL}_{n}(\mathbb{R})$,

$\exp(gMg^{-1})=\sum_{k=0}^{\infty}(gMg^{-1})^{k}/k!.$

Now note that

$(gMg^{-1})^{k}=(gMg^{-1})(gMg^{-1})\ldots(gMg^{-1})=gMM\ldots
Mg^{-1}=gM^{k}g^{-1}.$

So then

$\sum_{k=0}^{\infty}(gMg^{-1})^{k}/k!=\sum_{k=0}^{\infty}gM^{k}g^{-1}/k!=g\left(\sum_{k=0}^{\infty}M^{k}/k!\right)g^{-1}=g\exp(M)g^{-1}.$

#### Structured Network Output

So the last layers of the neural network will be $\exp\circ F$ which outputs
an invertible $2\times 2$ matrix and satisfies

$(\exp\circ F)(\rho_{reg}(\phi)f)=\rho_{1}(\phi)(\exp\circ
F)(f)\rho_{1}(-\phi).$

##  Appendix C Algorithms

###  C.1 Empirical method for estimating MRS for general 2D distribution

We present an empirical method to estimate MRS given any 2D distribution. We
first discretize the 2D space that is the domain of the probability
distribution into $n$ grid cells, indexing them as $g_{i}$ where
$i\in\\{1,\ldots,j\\}$. Then we can numerically estimate the density of each
grid cell $p(g_{i})$ for $i=1,\ldots,n$. The confidence region can be
estimated by selecting the grid cells with the highest likelihoods until they
sum to $1-\alpha$.

$R(1-\alpha)=\bigcup_{i\in K}g_{i},\quad K=\text{argmin}_{K}:\Sigma_{i\in
K}p(g_{i})\geq 1-\alpha$

where $K$ is the superlevel set of grid cells $g_{i}$ whose empirical
probability density is at least $1-\alpha$.

Similarly, we define the estimated region

$R(z_{i})=\bigcup g_{i}\quad\text{ for }p(g_{i})\geq p(z_{i})$

as the superlevel set of data point $z_{i}$. This way, we can calculate MRS
for any 2D using the estimated $R(1-\alpha)$ and $R(z_{i})$ in the MRS
formula.

$\mathrm{MRS}(R;Z)=\frac{1}{N}\sum_{i=1}^{N}[|R(1-\alpha)|+\frac{1}{\alpha}|R(z_{i})\backslash
R(1-\alpha)|\mathbb{1}\\{z_{i}\in R^{C}(1-\alpha)\\}]$

Figure 5 illustrates the idea, this method is used to estimate the MRS scores
for multimodal methods Trajectron++ and MFP in our experiments.

Figure 5: Illustration for empirical MRS estimation, where we calculate the
probability density at each grid, and select the highest density grids to
construct $R(z_{i})$ and $R(1-\alpha)$. This way, we may approximate MRS
numerically.

###  C.2 Conformal algorithm for 2-D Time Series

We extend the inductive conformal prediction (ICP) algorithm [Shafer and Vovk,
2008] to conformally estimate a confidence region for 2D forecasts. Given an
exchangeable set of trajectory observations
$\mathcal{D}=\\{(\mathbf{x}_{1:t}^{(i)},\mathbf{e}^{(i)}),\mathbf{x}_{t+1:t+k}^{(i)}\\}_{i=1}^{l}$
and a new sample $(\mathbf{x}_{1:t}^{(l+1)},\mathbf{e}^{(l+1)})$, the ICP
algorithm returns $k$ prediction
intervals,$[\Gamma^{\alpha}_{1},\ldots,\Gamma^{\alpha}_{k}]$, one for each
timestep, such that:

$\mathbb{P}[\>\forall
j\in\\{1,\ldots,k\\},\>x_{t+j}\in\Gamma^{\alpha}_{j}(\mathbf{x}^{(l+1)})\>]\geq
1-\alpha$

for any underlying predictive model. This inequality is the validity
condition. [Stankevičiūtė et al., 2021] proves the 1-D case for time-series
forecasting validity; we leave proving it for the 2-D case for future work.

ICP requires an underlying forecasting model and a choice of nonconformity
score. In our case, LSTM and ECCO are both models that take as input
$(\mathbf{x}_{1:t},\mathbf{e})$ and outputs $\mathbf{x}_{t+1:t+k}$. We select
root mean square error (RMSE) as our nonconformity score as it is commonly
used in the 1-d setting and naturally extends to high dimensions. We split the
training set into the proper training set and a calibration set of equal size:
$\mathcal{D}=\mathcal{D}_{train}\bigcup\mathcal{D}_{cal}$. We train our model
$M$ on $\mathcal{D}_{train}$ and obtain the critical nonconformity scores
$\hat{\gamma_{1}},\ldots,\hat{\gamma_{k}}$ as in algorithm 1 in [Stankevičiūtė
et al., 2021]. Then, for input $(\mathbf{x}_{1:t},\mathbf{e})$ and
$\hat{x}_{t+1:t+k}=M(\mathbf{x}_{1:t},\mathbf{e})$ we can construct a set of
confidence intervals

$\Gamma^{\alpha}_{j}(\mathbf{x}^{(l+1)}_{(1:t)})=\\{x\in\mathbb{R}^{2}\;|\;RMSE(x,\hat{x}_{j})\leq\hat{\gamma}_{j}\\}\;\text{
for }j\in\\{1,\ldots,k\\}$

We illustrate some example scenes with confidence regions provided by
conformal-ECCO in figure 6. Note that the regions for timestep $j$ are given
by $\hat{\gamma}_{j}$, hence are the same for every scenario. With our choice
of RMSE as the nonconformity score, the region is a circle; they appear
elliptical due to scale compression.

Figure 6: Scene illustrations for confidence regions given by conformal-ECCO.
Note that in order to achieve 90% coverage, the region is larger than needed,
especially in straight-lane cases like left two scenarios. Due to the lack of
underlying assumption about the distribution, the model isn’t able to adjust
the shape of uncertainty according to data either.

##  Appendix D Experiment Details

###  D.1 Implementation Details

#### Dataset and Preprocessing

The Argoverse autonomous vehicle dataset contains 205,942 samples, consisting
of diverse driving scenarios from Miami and Pittsburgh. We split 90/10 into a
training set and validation set of size 185,348 and 20,594 respectively. The
official validation set of size 39,472 is used for testing and reporting
performance. We preprocess the scenes to filter out incomplete trajectories
and cap the number of vehicles modeled to 60. If there are less than 60 cars
in the scenario, we insert dummy cars into them to achieve consistent car
numbers. For map information, we only include center lanes with lane
directions as features. Similar to vehicles, we introduce dummy lane nodes
into each scene to make lane numbers consistently equal to 650.

We use the latest release of the TrajNet++ dataset (Update 4.0) for our
pedestrian experiments. TrajNet++ is a compiled set of pedestrian trajectories
captured in both indoor and outdoor locations such as in universities, hotels,
Zara, and train stations. The sample in this dataset is 21 timestamps long,
and the goal is to predict the 2D spatial positions for each pedestrian in the
future 12 timestamps given the first 9. The pedestrian dataset contains
240,896 samples, which we split 80/10/10 into train, validation, and test
sets. Similar to Argoverse, we filter out incomplete trajectories in
processing and either cap or insert dummy pedestrians such that each scene has
60 agents. No map information was used in the pedestrian dataset.

We include our code for preprocessing, model implementation, and training in
the supplementary materials.

#### Hyperparameters and Training Details

We trained the PECCO model with 4 equivariant continuous convolution layers of
hidden size of (8, 16, 16, 16) respectively for Argoverse, and (4, 8, 16, 16)
for the Trajnet++ pedestrian dataset. Our models are all trained with Adam
optimizer with a base learning rate $r=0.001$, and linear learning rate
scheduler set to $\gamma=0.95$. For Argoverse task, we set the CtsConv radius
to be 40, and for the pedestrian task we set it to be 6. PECCO’s Arogoverse
model has 129k parameters. For comparison, Trajectron++ has 127k and MFP has
67K.

All our models without map information are trained for 10K iterations with
batch size 32 with learning rate updating every 150 iterations. Most of our
experiments are performed on a server with 4 RTX 2080 Ti GPUs, and it takes
around 9-12 hours to finish training. We run each experiment 3 times with
different random initialization and data order. The numbers reported in tables
1 2, and 5 are the mean and standard deviation of those 3 runs.

###  D.2 Deterministic Baseline Results

We present numbers for deterministic baseline models for Argoverse and
Trajnet++. By sampling, we achieve better minADE/minFDE performance with
probabilistic models, but they serve as valuable baselines to illustrate the
difficulty of the task.

| Argoverse | Trajnet++
---|---|---
Model | ADE$\downarrow$ | FDE$\downarrow$ | ADE$\downarrow$ | FDE$\downarrow$
Constant Velocity | 2.77 | 6.16 | 1.21 | 2.37
Nearest Neighbor | 3.52 | 7.85 | 1.25 | 2.61
LSTM | 1.97 $\pm$ .05 | 4.98 $\pm$ .31 | 1.01 $\pm$ .02 | 1.98 $\pm$ .08
CtsConv | 1.87 $\pm$ .06 | 4.43 $\pm$ .28 | 1.35 $\pm$ .05 | 2.97 $\pm$ .16
ECCO | 1.68 $\pm$ .04 | 3.98 $\pm$ .19 | 0.94 $\pm$ .01 | 2.05 $\pm$ .03
Table 5:  Deterministic baselines on the Argoverse and TrajNet++ dataset

###  D.3 Visualizations

We include 3 more scenarios for qualitatively comparison between models. They
provide insights consistent with Figure 4: LSTM-MRS-r tends to overestimate
uncertainty while CtsConv models tend to underestimate collapsing into one
mode. PECCO, on the other hand, forecasts a tighter probability that covers
the possibilities well.

(a) LSTM-MRS-r

(b) CtsConv-NLL

(c) CtsConv-MRS

(d) PECCO

Figure 7: Another lane change example. Red trajectories is the agent of
interest in the same scenario. CtsConv-NLL and CtsConv-MRS captured two
different modes of futures, whereas PECCO is able to cover both.

(a) LSTM-MRS-r

(b) CtsCov-NLL

(c) CtsConv-MRS

(d) PECCO

(e) LSTM-MRS-r

(f) CtsConv-NLL

(g) CtsConv-MRS

(h) PECCO

Figure 8: Comparison of uncertainty predicted at 2 crossroads. Red
trajectories is the agent of interest in the same scenario. Note how the LSTM
predicted uncertainty explodes at the final time step, and PECCO is able to
model the possibility of turns.