# Fortran High-Level Synthesis: Reducing the barriers to accelerating HPC
codes on FPGAs

Gabriel Rodriguez-Canal § EPCC, The University of Edinburgh  Nick Brown  EPCC,
The University of Edinburgh  Tim Dykes  HPC/AI EMEA Research Lab, Hewlett
Packard Enterprise  Jessica R. Jones  HPC/AI EMEA Research Lab, Hewlett
Packard Enterprise  Utz-Uwe Haus  HPC/AI EMEA Research Lab, Hewlett Packard
Enterprise

###### Abstract

In recent years the use of FPGAs to accelerate scientific applications has
grown, with numerous applications demonstrating the benefit of FPGAs for high
performance workloads. However, whilst High Level Synthesis (HLS) has
significantly lowered the barrier to entry in programming FPGAs by enabling
programmers to use C++, a major challenge is that most often these codes are
not originally written in C++. Instead, Fortran is the lingua franca of
scientific computing and-so it requires a complex and time consuming initial
step to convert into C++ even before considering the FPGA.

In this paper we describe work enabling Fortran for AMD Xilinx FPGAs by
connecting the LLVM Flang front end to AMD Xilinx’s LLVM back end. This
enables programmers to use Fortran as a first-class language for programming
FPGAs, and as we demonstrate enjoy all the tuning and optimisation
opportunities that HLS C++ provides. Furthermore, we demonstrate that certain
language features of Fortran make it especially beneficial for programming
FPGAs compared to C++. The result of this work is a lowering of the barrier to
entry in using FPGAs for scientific computing, enabling programmers to
leverage their existing codebase and language of choice on the FPGA directly.

###### Index Terms:

FPGAs, Fortran, High Level Synthesis, HPC

©2023 IEEE. Personal use of this material is permitted. Permission from IEEE
must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works.

§§footnotetext: Corresponding author<EMAIL_ADDRESS>
##  I Introduction

First introduced in the late 1950s, Fortran has been the lingua franca of
scientific computing for over 60 years [1]. Whilst many alternative
programming languages have come and gone, it has regained its popularity for
writing high performance codes. Indeed, over 80% of the applications running
on ARCHER2, a 750,000 core Cray EX which is the UK national supercomputer, are
written in Fortran. This ubiquity has meant that there are many applications,
both actively maintained and legacy, that use the language ranging from
weather forecasting models used by the world’s leading meteorology
organisations [2] [3] to simulation codes used for modelling state of the art
aircraft and jet engines [4] [5].

As the HPC community turns to exploring accelerators other than GPUs, for
potential performance and energy usage advantages, Field Programmable Gate
Arrays (FPGAs) have significant potential. Whilst FPGAs are yet to gain wide
acceptance in HPC, there are a number of testbeds and smaller HPC systems that
contain FPGAs. Numerous studies report that the ability to tailor the logic
and select a bespoke on-chip memory configuration can be beneficial compared
to general purpose architectures, especially for those codes which are memory
bound [6] [7]. However, a major blocker is that currently, when using High
Level Synthesis (HLS), one must write their codes in C or C++. Consequently
there is a significant initial overhead for many HPC developers where they
must first convert their Fortran code into C++, before getting anywhere near
an FPGA. This not only requires expertise in another language, but furthermore
is a time consuming and error prone process due to the differences between
Fortran and C++, for instance in array index ordering (row major in C verses
column major in Fortran), and default array start indices.

This paper describes work bringing Fortran programming to FPGAs by integrating
Flang with AMD Xilinx’s LLVM-based HLS back end. The objective is to lower the
barrier to entry in porting existing HPC codes to FPGAs, and ultimately
increase the accessibility of FPGAs to the HPC community. The rest of the
paper is structured as follows, Section II describes the background to this
work in more detail and surveys previous efforts bringing Fortran to FPGAs.
Section III reports the hardware and software configurations used in this
work, before Section IV then explores our Fortran HLS solution. A performance
comparison of our work against the existing C/C++ HLS approach is undertaken
in Section V, highlighting benefits of Fortran for writing HLS code. Section
VI then draws conclusions and describes further work.

The contributions of this paper are:

  1. 1.

Describing an approach to enabling direct support for Fortran on AMD Xilinx
FPGAs through the Xilinx Vitis ecosystem.

  2. 2.

A case study of leveraging AMD Xilinx’s open source LLVM front end components
to integrate new programming languages and tools into Vitis HLS, without
changes needed to the standard Flang front end or Xilinx back end.

  3. 3.

A comparison of performance between Fortran HLS and C/C++, demonstrating how
Fortran specific language constructs, such as dynamic N-dimensional arrays,
can provide an advantage over the existing HLS C/C++ flow.

##  II Background and related work

A common approach to programming FPGAs is to use C or C++ via HLS, with both
AMD Xilinx and Intel providing their own HLS solutions. Several tools have
been developed [8] [9] which aim to undertake automatic code transformations,
often via bespoke MLIR dialects [10], in order to enable non-experts to
program FPGAs effectively. There are mixed results when using these tools, and
code must be in C or C++ to start with, but it demonstrates the ability of
third parties to develop their own transformations and then integrate with
existing tools.

DACE [11] is an alternative approach which is more focused on HPC, where
developers write code in Python and this is then transpiled to a variety of
target architectures including FPGAs (C/C++ HLS for AMD Xilinx and OpenCL for
Intel FPGAs). However, the limitation is that one must first undertake
conversion of their code to integrate with DACE before running on FPGAs,
similar to the disadvantages associated with having to first convert Fortran
to C/C++ for HLS.

TyTra [12] is a compiler which synthesises Fortran to HDL and is intended as a
platform for undertaking auto-optimisation steps to automatically convert the
programmer’s Von Neumann-based code into an efficient dataflow representation.
Whilst, at the time of writing TyTra is not actively maintained, it is an
interesting comparison point. The major disadvantages with TyTra is that one
must rely on their own bespoke compiler resulting in risk (whether the
compiler will continue to be maintained in the future), potentially poorer
future performance (as they will not automatically benefit from the
improvements made by vendors to their tool chains) and a lack of integration
with the wider ecosystem (such as the Vitis profiling tools).

###  II-A AMD Xilinx’s LLVM-based HLS tooling

AMD Xilinx’s HLS synthesis tool is based upon LLVM and, as illustrated in
Figure 1, this comprises two parts; a front end built upon Clang [13] which
compiles a programmer’s C/C++ code to LLVM’s Intermediate Representation
format known as LLVM-IR, and a back end which synthesises LLVM-IR to the
underlying HDL level for their FPGAs. Whilst much of Clang in Vitis HLS has
remained unchanged, AMD Xilinx have added support for their own bespoke HLS
pramas which correspond to additional constructs in the generated LLVM-IR.

Figure 1: Default workflow of AMD Xilinx Vitis v++ Figure 2: Overall workflow
of our Fortran, f++, tool where Fortran is converted into the _xpirbc_ format
by our tool and this is then provided to the AMD Xilinx back end which
generates the HDL.

AMD Xilinx have released their HLS Clang front end via the Illinois/NCSA open
source licence, and this is kept up to date with each release of the Vitis
suite, although their back end remains proprietary. Based upon these
components it is possible to generate LLVM-IR and provide this to the AMD
Xilinx back end for it to synthesise to HDL.

###  II-B LLVM Flang

Flang [14] is the open source Fortran compiler provided as part of the LLVM
compiler suite. Based around an active and growing open source community, this
compiler was developed from the ground-up around five years ago and replaced
the previous Flang compiler in the main LLVM repository around two years ago.
Built around MLIR, a framework in LLVM for defining bespoke IR dialects and
mixing these together, Flang currently supports the Fortran 2003 standard with
several features from later standards, such as Fortran 2008 and Fortran 2018,
also supported and many others in development.

Ultimately, Flang is a front end for the Fortran language and after language
specific lexing, parsing, and optimisation, generates LLVM-IR. This is then
passed to a specific LLVM back end and built into the final executable, for
instance targeting CPUs or GPUs.

##  III Experimental setup

An AMD Xilinx U280 FPGA is used for the FPGA experiments conducted in this
paper, at the default clock frequency of 300MHz. This contains an FPGA chip
with 1.08 million LUTs, 4.5MB of on-chip BRAM, 30MB of on-chip UltraRAM, and
9024 DSP slices. This PCIe card also contains 8GB of High Bandwidth Memory
(HBM2) and 32GB of DDR DRAM. All bitstreams are built for the U280 using
Xilinx’s Vitis framework version 2021.2 which at the time of writing is the
latest compatible version for the U280. We use Flang version 16 throughout.
Reported runtime is the execution time of the HLS kernel on the FPGA and
excludes the overhead of data transfer to or from the device. Kernel execution
time is gathered via profiling information on the OpenCL event.

For the CPU comparison runs, which are intended to provide context around the
optimisations performed by the Clang and Flang front ends respectively, were
undertaken on a single core of a Xeon Platinum (Cascade Lake) 8260M CPUs
running at 2.40GHz. All reported performance times are based on wall clock
time and numbers are averaged over 10 runs.

##  IV Fortran FPGA programming via Flang HLS

The behaviour of AMD Xilinx’s Vitis v++, the tool that synthesises C/C++ code
by calling HLS, can be altered to provide AMD Xilinx’s LLVM-based back end
externally generated LLVM-IR as long as the file has the extension .xpirbc. We
leverage this approach to enable Fortran on Xilinx’s FPGAs by connecting LLVM
Flang to the AMD Xilinx back end via Flang’s generated LLVM-IR, with the major
objective being to leave LLVM Flang and the AMD Xilinx back end unchanged.

There are numerous challenges that must be overcome for AMD Xilinx’s HLS back
end to properly interpret the LLVM-IR. This is because Flang must conform to
certain requirements when generating LLVM-IR based upon what is expected by
the AMD Xilinx back end. In summary, these are:

  1. 1.

The generated LLVM-IR must conform to LLVM version 7.

  2. 2.

HLS pragmas must be lowered into AMD Xilinx’s specific LLVM-IR metadata or
directives.

  3. 3.

HLS Streams are defined and managed with custom IR primitives generated by AMD
Xilinx’s custom Clang front end.

The challenge is that none of these requirements are directly met by Flang.
For instance, at the time of writing the latest version of LLVM is version 16,
in fact version 7 was from 2018 and is not supported by the latest Flang.
Similarly, Flang has no internal knowledge of HLS streams or HLS annotations
and as a result, the IR that is generated is not compatible with the back end
unless these are explicitly added.

To address these issues, our approach places a layer between Flang and AMD
Xilinx’s back end which, as illustrated in Figure 2, performs passes and
regular expressions on the Flang generated LLVM-IR to transform it into a form
that is acceptable to the back end. These are all wrapped in our f++ tool,
which launches separate components of the Fortran HLS workflow.

###  IV-A Lowering the generated LLVM-IR to version 7

The differences between the latest versions of LLVM-IR and version 7 will
either result in the LLVM-IR generated by Flang being rejected by the AMD
Xilinx HLS back end, or accepted and subsequent erroneous results. The later
can be especially challenging, for instance if the LLVM-IR contains the byval
attribute for memory pointers, which is used by the latest LLVM version, then
HLS will succeed and generate HDL, but after running the synthesis, place, and
route stages the resulting bitstream will be incompatible with the Xilinx
Runtime library, XRT. Whilst it might seen an obvious approach for Xilinx to
update their LLVM back end to accept the latest LLVM-IR, this is non-trivial
because of the differences between LLVM-IR version 7 and the current version.
An example of this is that, as of version 16, LLVM-IR only supports opaque
pointers at the IR level which are incompatible with HLS. Furthermore, the
rapid pace of change in and around LLVM means that it is not always possible
for mature products, such as AMD Xilinx’s HLS, built upon this technology to
track the latest bleeding-edge version.

Therefore, we developed a pass that will downgrade the LLVM-IR to version 7.
This was initially based upon the prototype developed by the Khronos team [15]
as part of their integration with SYCL. However, their work was built for
Clang and only partially amends the generated LLVM-IR. For example, our pass
also transforms procedure names into the format expected by the back
end.Furthermore, to enable compatibility with Flang our pass also transforms
procedure argument names, which cannot begin with a number, but are generated
as such by Flang. Finally, we strip the LLVM-IR of all metadata which does not
conform to LLVM-IR version 7.

###  IV-B HLS pragmas

Our research prototype supports a subset of the most commonly used HLS
pragmas, and it is trivial to add others following the same approach we have
developed. It was decided to concentrate on those pragmas which are
fundamental for controlling performance, namely the _pipeline_ , _unroll_ ,
_array_partition_ , and _dataflow_ pragmas, as well as the _interface_ pragma
which enables programmers to define the data interfaces of the generated HLS
IP block. Crucially, all pragmas are handled in a similar way, and supporting
a pragma involves the following steps:

  1. 1.

Preprocessing to replace the pragma by a call to a placeholder Fortran
subroutine with no arguments. This will lead to the Flang compiler generating
a placeholder function call in the LLVM-IR, and we followed this approach
because it does not introduce extra instructions before or after the call. Any
arguments passed to the pragma are encoded as part of the subroutine name and
this approach has been adopted because, if we were to pass them as arguments,
then Flang will generate corresponding values in the IR, which makes it more
difficult to then process and can in some circumstances result in bottlenecks
in the generated HDL.

  2. 2.

Run our bespoke LLVM pass which operates on the Flang generated LLVM-IR. Each
placeholder subroutine call is identified and any arguments are decoded and
extracted. The corresponding LLVM-IR constructs for the specific pragma that
this corresponds to are inserted into the IR as appropriate.

  3. 3.

Remove placeholder subroutine calls from the LLVM-IR

We followed this approach because it means that we can undertake pragma
replacement via regular expressions in a preprocessing phase, rather than
having to modify Flang itself to be aware of the pragmas. This would have
required changes to Flang which would then diverge it from the main tree.

An important aspect to highlight, and a further challenge, is how the
contention of instructions operate in loop trees in LLVM-IR, as this is needed
to correctly correlate the _pipeline_ and _unroll_ pragmas with the loops that
they are operating over. In LLVM, given a loop $L$ and any of its sub loops
$S$ and the instruction $I$, $I\in S\implies I\in L$. Therefore, our pass must
identify the deepest sub loop in the loop tree that contains $I$ using depth
first search with backtracking. Otherwise, for nested loops, if the
corresponding LLVM-IR was generated at the location of the place holder
subroutine call, then the parent loop would be qualified in every case
instead, altering the semantics of the program. Figure 3 illustrates the
equivalence of a series of nested _for_ loops and their corresponding loop
tree, where we undertake this search for the annotated loop.

Figure 3: Nested for loops and their loop tree, traversed with DFS with
backtracking. Pointed by the red arrow, the loop containing instruction I.
Black nodes contain the instruction.

###  IV-C HLS streams and polymorphism

In AMD Xilinx’s HLS, the C++ template mechanism is leveraged to implement the
HLS streams library and this uses polymorphism to enable many different
datatypes to be streamed. Whilst polymorphism is available from Fortran 2008
onwards, we decided to keep our tool backwards-compatible with Fortran 90
because that is the most common Fortran version used by HPC codes.
Consequently, it was undesirable to encode any features that required later
versions of Fortran. Therefore, we replicate the polymorphic behaviour with a
library of preprocessor macros and preprocessing of the Fortran code using
regular expressions.

⬇

1proto_hls_stream(integer)

2

3subroutine f(a, b)

4 use hls_stream

5

6 integer, dimension(100) :: a, b

7 type(HLSStream) :: s

8

9 set_hls_stream_type(s, integer)

10

11 do i = 1,100

12 call hls_write(a(i), s)

13 end do

14

15 do i = 1,100

16 b(i) = hls_read(s)

17 end do

18end subroutine

Listing 1: Example Fortran program using our Fortran HLS Stream macro library.

Listing 1 sketches an example of Fortran code using our macro-based HLS Stream
library. The stream functions provided in Vitis HLS (read, write, empty and
full) are also provided in Fortran HLS. The module is instantiated in line 1
for the types that will be used with streams in the code, just an integer in
this example, with twenty Fortran types and precisions supported overall by
the macro. Listing 2 illustrates the macro expansion which will replace line 1
in Listing 1 by this preprocessor generated code.

The generated _hls_stream_ module contains the user derived type HLSStream
which enables Flang to check the correct use of the library at compile-time.
For the types provided by the user in the macro function proto_hls_stream, the
corresponding set of HLS stream subroutines are generated for each of these
types inside the module. For instance, the _hls_lowered_write_integer_
subroutine at line 8 in Listing 2 has been generated to accept an integer
because that is the type provided to proto_hls_stream.

Similarly to the template instantiation mechanism in C/C++ HLS for typing
streams, the macro at line 9 of Listing 1 is equivalent to hls::stream<int> s.
This macro is expanded to call set_depth_integer(s%data_integer, 0) and has
two purposes. Firstly, the type of the stream is added to a dictionary that
keeps track of the stream types for each subroutine in the code. Secondly, the
function call in the IR generated by Flang is replaced by the corresponding
AMD Xilinx LLVM-IR primitive @llvm.fpga.set.stream.depth. This primitive is
fundamental for AMD Xilinx’s HLS back end to correctly handle streams, because
otherwise it will attempt to implement these as Master AXI-4 interfaces.

⬇

1module hls_stream

2 type(HLSStream)

3 end type

4

5 type(HLSStream_integer)

6 integer :: data_integer

7 end type

8 subroutine hls_lowered_write_integer(

9 input_data, stream_data)

10 integer, value :: input_data

11 integer :: stream_data

12 end subroutine

13 …

14end module

Listing 2: Expansion of proto_hls_stream macro. Only the generated
hls_lowered_write_integer subroutine shown.

This approach of representing stream polymorphic types as a Fortran derived
type, similar to a _typedef_ in C, with one field of the derived type being
the type of the stream, is beneficial because we found that it ensures Flang
generates similar LLVM-IR to Clang when accessing the member of a structure.
This makes it more straightforward to locate the argument(s) of the AMD
Xilinx’s LLVM-IR primitive and/or the result of invoking such functions.

##  V Evaluation

TABLE I: TeaLeaf mini-app kernels in HLS using C and Fortran (unoptimised) for the FPGA, also comparing against corresponding CPU performance. |  | Alveo U280 FPGA | Xeon Platinum CPU
---|---|---|---
Benchmark | Problem size |  HLS C Runtime(s) |  HLS Fortran Runtime(s) |  C Runtime(s) |  Fortran Runtime(s)
calc_2norm | 480000 | 0.056 | 0.038 | 0.002 | 0.089
calc_residual | 480000 | 0.196 | 0.215 | 0.008 | 0.087
cg_calc_p | 480000 | 0.315 | 0.319 | 0.002 | 0.003
cg_calc_ur | 480000 | 0.781 | 0.947 | 0.015 | 0.034
cg_calc_w | 480000 | 0.198 | 0.201 | 0.009 | 0.006
cg_calc_w_norxy | 480000 | 0.317 | 0.299 | 0.008 | 0.005
cg_init | 480000 | 0.214 | 0.214 | 0.007 | 0.007
cheby_init | 480000 | 0.850 | 0.845 | 0.025 | 0.093
cheby_iterate | 480000 | 0.628 | 1.072 | 0.016 | 0.083
common_init | 480000 | 0.711 | 0.631 | 0.025 | 0.093
field_summary | 480000 | 0.105 | 0.054 | 0.008 | 0.095
finalise | 480000 | 0.003 | 0.055 |  $\approx$0.000 | 0.107
generate_chunk | 480000 | 0.369 | 0.898 | 0.008 | 0.096
initialise_chunk | 480000 | 0.116 | 0.143 |  $\approx$0.000 | 0.048
jacobi_solve | 480000 | 0.192 | 0.230 | 0.018 | 0.001
ppcg_calc_rrn | 39677401 | 3.871 | 3.866 | 0.184 | 0.203
ppcg_calc_zrnorm | 480000 | 0.059 | 0.059 | 0.002 | 0.086
ppcg_init | 480000 | 0.343 | 0.509 | 0.009 | 0.109
ppcg_init_sd | 480000 | 0.227 | 0.188 | 0.005 | 0.069
ppcg_inner | 480000 | 0.754 | 0.727 |  $\approx$0.000 | 0.106
ppcg_inner_norxy | 480000 | 0.748 | 0.704 |  $\approx$0.000 | 0.098
ppcg_pupdate | 480000 | 0.043 | 0.044 | 0.002 | 0.084
ppcg_store_r | 480000 | 0.043 | 0.044 | 0.002 | 0.086
ppcg_update_z | 480000 | 0.043 | 0.044 | 0.002 | 0.099
set_field | 480000 | 0.043 | 0.044 | 0.002 | 0.087
tea_block_init | 480000 | 0.616 | 0.615 | 0.005 | 0.020
tea_block_solve | 480000 | 0.098 | 0.267 |  $\approx$0.000 | 0.026
tea_diag_init | 480000 | 0.083 | 0.083 | 0.002 | 0.019
tea_diag_solve | 480000 | 0.083 | 0.099 | 0.002 | 0.019
update_halo | 480000 | 0.788 | 0.889 | 0.011 | 0.059
update_halo_cell | 480000 | 0.001 | 0.001 |  $\approx$0.000 | 0.008
update_in_halo_bt | 480000 | 0.016 | 0.016 |  $\approx$0.000 | 0.067
update_in_halo_cell_bt | 480000 | 0.002 | 0.002 |  $\approx$0.000 | 0.004
update_in_halo_cell_lr | 480000 | 0.049 | 0.080 | 0.001 | 0.006
update_in_halo_lr | 480000 | 0.874 | 0.898 | 0.011 | 0.067

We evaluate the work presented in this paper with TeaLeaf [16], an HPC mini-
app that solves the linear heat conduction equation using 5 point stencils
with implicit solvers. In total, this suite comprises 35 kernels which
corresponds to 2587 lines of code ranging from Conjugate Gradient iteration to
the Chebyshev method. Table I presents performance results for the baseline
implementation of these kernels on an AMD Xilinx Alveo U280 FPGA and a single
core of a Xeon Platinum CPU. It should be highlighted that the standard
deviation across the 10 averaged runtimes is negligible, only affecting the
4th decimal figure of the results, and hence this has been omitted for
brevity. It is not our intention in this paper to study the performance
differences between architectures for these benchmarks and Table I reports a
direct translation from the CPU to FPGA without any FPGA specific
optimisations being undertaken because, in our first experiment, we are
interested in the performance that our HLS Fortran tool can provide _out of
the box_. Indeed, the comparison between FPGA and CPU performance is provided
to highlight where there are performance differences resulting from the
activities of the front ends of Clang and Flang, for instance undertaking
language specific code-level optimisations.

Several high level observations can be made from the runtime results reported
in Table I. Firstly, it can be seen that there are performance differences
between the HLS C and HLS Fortran kernels, but in the main these differences
are fairly small. Furthermore, there is no clear pattern around whether C or
Fortran performs best with HLS. For instance, the _field_summary_ Fortran
kernel outperforms C, whereas the opposite is true for the _generate_chunk_
kernel. Table II reports resource usage for these TeaLeaf kernels when using
HLS C and HLS Fortran. The numbers in Table II are expressed as the percentage
of overall resources that are used on the Alveo U280, and it can be seen that
whilst the BRAM resource usage between HLS C and HLS Fortran is fairly
consistent, for the LUTs, FFs, and DSP slices it is more varied. In general
terms, HLS Fortran tends to require more resources, especially LUTs, compared
to HLS C but there are exceptions to this such as with the _cg_calc_p_ kernel.

| HLS C % usage | HLS Fortran % usage
---|---|---
Kernel | BRAM | LUT | FF | DSP | BRAM | LUT | FF | DSP
calc_2norm | 0.10 | 0.32 | 0.14 | 0.19 | 0.10 | 0.43 | 0.16 | 0.16
calc_residual | 0.10 | 0.50 | 0.30 | 0.32 | 0.10 | 0.64 | 0.35 | 0.29
cg_calc_p | 0.10 | 0.52 | 0.22 | 0.32 | 0.10 | 0.49 | 0.22 | 0.27
cg_calc_ur | 0.10 | 2.24 | 1.20 | 1.30 | 0.10 | 2.88 | 1.23 | 1.00
cg_calc_w | 0.10 | 0.69 | 0.36 | 0.42 | 0.10 | 0.80 | 0.39 | 0.32
cg_calc_w_norxy | 0.10 | 0.60 | 0.31 | 0.32 | 0.10 | 0.75 | 0.37 | 0.32
cg_init | 0.10 | 1.82 | 1.09 | 1.26 | 0.10 | 2.19 | 0.95 | 0.81
cheby_init | 0.10 | 2.15 | 1.31 | 1.43 | 0.10 | 2.50 | 1.19 | 0.98
cheby_iterate | 0.10 | 2.31 | 1.36 | 1.46 | 0.10 | 2.85 | 1.31 | 1.10
common_init | 0.10 | 2.29 | 1.42 | 1.13 | 0.10 | 2.82 | 1.53 | 0.79
field_summary | 0.10 | 0.44 | 0.24 | 0.22 | 0.10 | 0.56 | 0.26 | 0.27
finalise | 0.10 | 0.22 | 0.11 | 0.10 | 0.10 | 0.30 | 0.13 | 0.03
generate_chunk | 0.10 | 1.05 | 0.54 | 0.32 | 0.10 | 1.63 | 0.74 | 0.70
initialise_chunk | 0.10 | 0.91 | 0.33 | 0.12 | 0.10 | 1.07 | 0.46 | 0.22
jacobi_solve | 0.10 | 0.63 | 0.37 | 0.35 | 0.10 | 0.85 | 0.44 | 0.32
ppcg_calc_rrn | 0.10 | 0.37 | 0.17 | 0.19 | 0.10 | 0.45 | 0.20 | 0.23
ppcg_calc_zrnorm | 0.10 | 0.35 | 0.15 | 0.19 | 0.10 | 0.46 | 0.18 | 0.16
ppcg_init | 0.10 | 1.95 | 1.15 | 1.33 | 0.10 | 2.39 | 1.04 | 0.84
ppcg_init_sd | 0.10 | 0.54 | 0.26 | 0.22 | 0.10 | 0.68 | 0.30 | 0.20
ppcg_inner | 0.10 | 2.88 | 1.61 | 1.66 | 0.10 | 3.62 | 1.70 | 1.46
ppcg_inner_norxy | 0.10 | 2.55 | 1.43 | 1.46 | 0.10 | 3.16 | 1.43 | 1.14
ppcg_pupdate | 0.10 | 0.23 | 0.09 | 0.07 | 0.10 | 0.31 | 0.11 | 0.03
ppcg_store_r | 0.10 | 0.23 | 0.09 | 0.07 | 0.10 | 0.31 | 0.11 | 0.03
ppcg_update_z | 0.10 | 0.23 | 0.09 | 0.07 | 0.10 | 0.31 | 0.11 | 0.03
set_field | 0.10 | 0.23 | 0.09 | 0.07 | 0.10 | 0.31 | 0.11 | 0.03
tea_block_init | 0.10 | 0.78 | 0.35 | 0.48 | 0.10 | 0.89 | 0.35 | 0.35
tea_block_solve | 0.10 | 1.23 | 0.82 | 0.93 | 0.10 | 1.30 | 0.65 | 0.61
tea_diag_init | 0.10 | 0.26 | 0.38 | 0.07 | 0.10 | 0.36 | 0.42 | 0.07
tea_diag_solve | 0.10 | 0.24 | 0.12 | 0.16 | 0.10 | 0.39 | 0.16 | 0.12
update_halo | 0.20 | 1.95 | 0.44 | 0.20 | 0.74 | 2.76 | 0.65 | 0.31
update_halo_cell | 0.20 | 1.35 | 0.32 | 0.20 | 0.20 | 1.40 | 0.40 | 0.31
update_in_halo_bt | 0.74 | 1.48 | 0.34 | 0.13 | 0.74 | 1.67 | 0.42 | 0.17
update_in_halo_cell_bt | 0.10 | 0.40 | 0.15 | 0.13 | 0.10 | 0.56 | 0.23 | 0.17
update_in_halo_cell_lr | 0.10 | 0.39 | 0.18 | 0.13 | 0.10 | 0.50 | 0.22 | 0.13
update_in_halo_lr | 0.74 | 1.52 | 0.35 | 0.13 | 0.74 | 1.63 | 0.38 | 0.13
TABLE II: Percentage of resource utilisation on the Alveo U280 FPGA for
TeaLeaf mini-app kernels written in HLS using C or Fortran.

⬇

1int x_dim = x_max - x_min;

2int y_dim = y_max - y_min;

3int w_x = x_max - x_min +

4 2 * halo_exchange_depth + 1;

5for(int k=0; k <= y_dim; k++) {

6 for(int j=0; j <= x_dim; j++) {

7 w[(k + halo_exchange_depth) * w_x +

8 (j + halo_exchange_depth)] = …

9 }

10}

Listing 3: Access to a flattened 2D array with manually calculated offsets in
C for _cg_calc_w_norxy_ kernel.

Whilst the same algorithms from the TeaLeaf suite are being compiled using HLS
C and HLS Fortran, the performance and resource usage differences between
these flows reported in Tables I and II demonstrate that there are
implications of whether one uses C or Fortran to drive HLS. This is
interesting because, based upon the fact that thus far we are using Von-
Neumann based algorithms, then one would expect the unoptimised nature of
these codes to dominate and performance to be very similar or the same between
HLS C and HLS Fortran. Indeed our objective in this work has been for HLS
Fortran to match the performance of HLS C. However, clearly the results are
more nuanced and we focus on two kernels to explore these differences in more
detail. Starting with the _cg_calc_w_norxy_ kernel which undertakes Conjugate
Gradient, as shown in Table I, this performs better on the FPGA in Fortran
compared to C. The reason for this is because of the explicit calculation of
offsets for accessing the flattened 2-dimensional array, as shown in Listing
3. Vitis imposes a limitation that external arrays have a maximum of one
dimension, in C this must be handled manually by the programmer but with
Fortran it is possible to use N-dimensional arrays. This is because dynamic
arrays require pointer to pointers in C and are not supported by Vitis HLS
[17]. Arrays in Fortran are higher-level than in C, and-so during compilation
Flang is able to calculate the offsets and flatten the dimensions of the
arrays statically, rather than requiring this be undertaken at runtime,
yielding an IR compatible with the aforementioned restriction. The
programmer’s Fortran code is illustrated in Listing 4, and this demonstrates
an important advantage of using Fortran compared to C/C++ for programming
FPGAs, where the additional expressiveness of Fortran, in this case via
N-dimensional arrays, means that the compiler is able to undertake such offset
calculations statically which are low level, time consuming, and can be error
prone.

⬇

1REAL(KIND=8),

2DIMENSION(x_min-halo_exchange_depth:

3 x_max+halo_exchange_depth,

4 y_min-halo_exchange_depth:

5 y_max+halo_exchange_depth) :: w

6

7do k=y_min,y_max

8 do j=x_min,x_max

9 w(j,k) = …

10 end do

11end do

Listing 4: Access to a 2D array in Fortran for _cg_calc_w_norxy_ kernel.

The _ppcg_calc_rrn_ kernel, which undertakes a reduction and whose runtime is
reported in Table I, is also slightly faster using our Fortran-based approach
compared to C on the FPGA. It is instructive to compare the runtime between C
and Fortran on the FPGA against that of the CPU to understand whether
performance differences come from these respective front ends. For instance,
it can be seen that for the _cg_calc_w_norxy_ kernel the Fortran runtime is
lower than the C runtime for both the FPGA and CPU, whereas this is not the
case for the _ppcg_calc_rrn_ kernel where Fortran is slightly faster than C on
the FPGA but slower on the CPU. It can be seen that there is no clear pattern,
and just because one language is faster than the other on the FPGA then it
does not mean that the same will be true on the CPU, or vice versa.
Consequently, it can be deduced that it is the combination of both the front
end and back end which is determining runtime performance and differences
between languages are at the individual kernel level.

Results reported in Table I are unoptimised for the FPGA, and whilst this was
an important baseline to understand the performance differences delivered for
such codes, it is well known that Von-Neumann based codes must be transformed
to their dataflow counterparts for best performance on an FPGA [18]. A key
question was therefore whether support in our tool was rich enough to enable
these same optimisations to be undertaken at the code level in HLS via
Fortran, as is currently possible in C or C++. We focus on optimising the
_cg_calc_w_norxy_ and _ppcg_calc_rrn_ kernels, with the results of optimising
these kernels reported in Table III. We optimised the _cg_calc_w_norxy_
Conjugate Gradient stencil kernel using shift registers. This technique
replaces the strided memory accesses by sequential accesses to a shift
register, where every element can be accessed in one cycle due to the initial
prefetching of part of the array and a sliding window, where the access to new
elements from memory is pipelined via the shift register operation [6]. It can
be seen that this optimisation reduces the runtime of both the C and Fortran
versions compared to their unoptimised counterparts.

TABLE III: Optimised kernels in HLS C and HLS Fortran for the FPGA, and also a synthetic streaming kernel. Problem size of 480000 for _cg_calc_w_norxy_ , 39677401 for _ppcg_calc_rrn_ , and 32 million elements for Synthetic Streaming. | Alveo U280 FPGA
---|---
Benchmark |  HLS C Runtime(s) |  HLS Fortran Runtime(s)
cg_calc_w_norxy | 0.200 | 0.222
ppcg_calc_rrn (partial sums) | 0.630 | 0.631
ppcg_calc_rrn (dataflow) | 0.820 | 0.535
Synthetic Streaming | 0.131 | 0.132

The bottleneck with the _ppcg_calc_rrn_ reduction kernel is in the floating
point adder, which requires 7 cycles to complete a sum, therefore leading to
an initiation interval of 7 because results from consecutive sum operations
are required. We ameliorate this with partial sums, where Vitis pipelines the
adder to undertake seven independent sum operations every 7 cycles instead
[19], and then at the end adds these seven partial sums together to calculate
the overall result.

There are two versions of the _ppcg_calc_rrn_ kernel reported in Table III,
_partial sums_ uses the partial sum technique only and it can be seen that
both Fortran and C deliver very similar performance on the FPGA, which is
significantly faster than the unoptimised version in Table I. The second
version of the _ppcg_calc_rrn_ kernel, _dataflow_ , uses the _dataflow_ pragma
and HLS streams to pipeline both the compute and interaction with external
memory to load and store data.

This second optimised version of the _ppcg_calc_rrn_ kernel, using dataflow
regions and HLS streams, yielded surprising results, as it can be seen from
Table III that the Fortran version is 53.27% faster than the C HLS code. Based
upon the report generated by Vitis HLS, we observed that Vitis is unable to
pipeline the reduction loop when using C, whereas it is able to do so in
Fortran. Because of the failure of Vitis to pipeline the loop in C, the
performance of the dataflow optimisation does not improve with respect to the
partial sums. Whilst our kernels all ran at 300Mhz on the FPGA, it should also
be noted that the maximum operating frequency reported by Vitis HLS is also
different, up to 371.47 MHz possible with Fortran and 357.68 MHz for C. It is
interesting to observe that small differences in the LLVM-IR resulting from
different language front ends, can have a significant impact.

Finally, we developed a synthetic example, Synthetic Streaming in Table III,
which sums two vectors from external memory as data is streamed in. This
kernel relies on three dataflow regions connected via HLS streams; loading
data from external memory, undertaking the element sum, and writing results
back to external memory. Compared to the _ppcg_calc_rrn_ kernel, these loops
are much simpler and-so can be pipelined by Vitis when using either the C and
Fortran flows, and both languages deliver similar performance.

TABLE IV: Percentage total resource usage in the optimised kernels and synthetic streaming in HLS C and HLS Fortran. Percentage overhead of Fortran also reported in brackets, in green the cases where Fortran leads to less resource utilisation. | HLS kernel usage C/Fortran (% total resources)
---|---
Benchmark | LUTs | FFs | BRAM | DSP
cg_calc_w _norxy | 6.13 / 5.76 | 6.32 / 5.13 | 3.17 / 3.17 | 0.74 / 0.74
| (-6.04%) | (-18.83%) | (0.00%) | (0.00%)
ppcg_calc_rrn (partial sums) | 0.65 / 0.60 | 0.34 / 0.30 | 0.40 / 0.40 | 0.19 / 0.22
| (-7.69%) | (-11.76%) | (0.00%) | (15.79%)
ppcg_calc_rrn (dataflow) | 0.64 / 0.78 | 0.34 / 0.40 | 0.40 / 0.40 | 0.19 / 0.39
| (21.88%) | (17.65%) | (0.00%) | (105.26%)
Synthetic Streaming | 0.29 / 0.22 | 0.10 / 0.08 | 0.15 / 0.10 | 0.00 / 0.00
| (-24.13%) | (-20.00%) | (-33.33%) | -

Another important consideration is the resource utilisation that results from
these tools when undertaking code optimisation, and whether one approach is
more frugal than another. Table IV reports resource utilisation for the four
optimised benchmarks that were explored in Table III and these can be compared
against the resource utilisation for unoptimised kernels reported in Table II.
For each benchmark and resource, the percentage of overall resource usage is
reported for the C and then the Fortran code, with these separated by a slash.
Underneath, in brackets, the percentage overhead of Fortran compared to C is
reported, where green is used to highlight situations where Fortran resulted
in reduced resource usage.

It can be seen in Table IV that for both _cg_calc_w_norxy_ and _ppcg_calc_rrn_
with partial sums kernels, Fortran requires fewer LUTs and FFs than C, with
BRAM and DSP being comparable. With the synthetic streaming benchmark, where
performance was comparable between the C and Fortran versions, it can be seen
that Fortran also uses fewer resources.

Fortran requires resources than C for the _ppcg_calc_rrn_ with dataflow
kernel, especially for DSP slices, however the hardware that is generated is
fundamentally different between these versions due to the lack of pipelining
on the reduction loop in C which led to lower performance. Therefore, overall,
these utilisation results demonstrate that Fortran is competitive against C
for resource utilisation when optimising the kernel and, in some cases, such
as synthetic streaming, can result in fewer resources with no loss in
performance.

##  VI Conclusions and further work

The work described in this paper enables programmers to leverage Fortran when
programming AMD Xilinx FPGAs, which is especially important for HPC workloads
as many of these are written in Fortran. Whilst Fortran programmers still need
to transform their Von-Neumann algorithms to a dataflow representation via
algorithmic restructuring and apply appropriate HLS pragmas, crucially our
tool avoids the time consuming and error prone preliminary step of converting
code into C/C++. This, combined with the fact that code remains in one single
language which HPC developers are familiar with, is beneficial in encouraging
HPC developers to explore accelerating their Fortran codes using FPGAs. The
tool that we have described connects LLVM’s Flang with the AMD Xilinx HLS back
end, without any modification required to either Flang or the back end. As a
result, the programmer is, apart from our conversion utility, using standard
tools that are maintained by the wider community. There is no divergence from
LLVM Flang or the AMD Xilinx HLS back end, and programmers benefit from
enhancements made both in Flang by the community and in the HLS back end by
AMD Xilinx. Furthermore, wider features of the Vitis ecosystem, such as
profiling, are also made available to the Fortran programmer.

We have compared performance of our approach against the existing C/C++ Vitis
HLS tooling. It was observed that, in the main, performance between these
tools is fairly comparable, although there are some differences for instance
with Fortran where the programmer can leverage higher level features, such as
N-dimensional arrays, that improves expressiveness and enables the compiler to
undertake additional optimisations. The major benefit of our approach is that
we unlock the ability for large HPC applications that are written in Fortran
to be run on FPGAs, and as we demonstrated these pragmas currently supported
by our tool enables these codes to be optimised for FPGAs, with additional
pragmas and constructs being trivial to add to the flow if required.

An objective has been for our HLS Fortran flow to match the performance of HLS
C so that HPC programmers are able to accelerate their Fortran codes on FPGAs
without any reduction in performance compared to leveraging C/C++. However, it
can be observed that some HLS Fortran kernels execute slightly slower than
their HLS C versions and consequently, for further work, optimisation at the
LLVM-IR level would be beneficial. For instance, the _cg_calc_w_norxy_ Fortran
kernel performs slightly worse than its C counterpart and we believe that this
is because of the use of long (64-bit) integers for array indicies in the IR
generated by Flang. Currently these require integer conversion from the
programmer’s 32-bit indicies to 64-bit in the IR, and this could be optimised
by our tool to avoid such conversions. It would also be interesting to apply
our code to programming the AI engines that are present on the Versal,
potentially providing a unified programming interface between the PL and AI
engines. Because AMD Xilinx have already developed MLIR dialects for their AI
engine compiler, and internally Flang uses MLIR via the FIR dialect, it should
be possible to connect these tools and leverage our HLS Fortran flow.
Furthermore, it would be interesting to explore generalising our tool, for
instance to support other programming languages, such as Python, Rust or
Julia, to be integrated with the AMD Xilinx back end.

We intend to make our tool open source and available at
_https://fpga.epcc.ed.ac.uk/community/fortran.html_ .

## Acknowledgment

The authors would like to thank HPE who funded this work via the EMEA Research
Lab internship programme, the ExCALIBUR xDSL project, and the ExCALIBUR H&ES
FPGA testbed and AMD Xilinx HACC program for access to compute resource used
in this work.

## References

  * [1] W. Brainerd, “The importance of Fortran in the 21st century,” _Journal of Modern Applied Statistical Methods_ , vol. 2, no. 1, p. 3, 2003.
  * [2] P. E. Strazdins, M. Kahn, J. Henrichs, T. Pugh, and M. Rezny, “Profiling methodology and performance tuning of the met office unified model for weather and climate simulations,” in _2011 IEEE International Symposium on Parallel and Distributed Processing Workshops and Phd Forum_. IEEE, 2011, pp. 1322–1331.
  * [3] N. Brown, M. Weiland, A. Hill, B. Shipway, C. Maynard, T. Allen, and M. Rezny, “A highly scalable Met Office NERC Cloud model,” _arXiv preprint arXiv:2009.12849_ , 2020.
  * [4] B. Yang, X. Fang, L. Zhang, F. Zhuang, M. Bi, C. Chen, G. Li, and X. Wang, “Applicability of empirical models of isentropic efficiency and mass flow rate of dynamic compressors to jet engines,” _Progress in Aerospace Sciences_ , vol. 106, pp. 32–42, 2019.
  * [5] O. Lehmkuhl, A. Lozano-Durán, and I. Rodriguez, “Active flow control for external aerodynamics: from micro air vehicles to a full aircraft in stall,” in _Journal of Physics: Conference Series_ , vol. 1522, no. 1. IOP Publishing, 2020, p. 012017\.
  * [6] N. Brown, “Accelerating advection for atmospheric modelling on Xilinx and Intel FPGAs,” in _2021 IEEE International Conference on Cluster Computing (CLUSTER)_. IEEE, 2021, pp. 767–774.
  * [7] M. A. Dávila-Guzmán, R. G. Tejero, M. Villarroya-Gaudó, and D. S. Gracia, “An Analytical Model of Memory-Bound Applications Compiled with High Level Synthesis,” in _2020 IEEE 28th Annual International Symposium on Field-Programmable Custom Computing Machines (FCCM)_. IEEE, 2020, pp. 218–218.
  * [8] H. Ye, C. Hao, J. Cheng, H. Jeong, J. Huang, S. Neuendorffer, and D. Chen, “ScaleHLS: A New Scalable High-Level Synthesis Framework on Multi-Level Intermediate Representation,” in _2022 IEEE International Symposium on High-Performance Computer Architecture (HPCA)_. IEEE, 2022, pp. 741–755.
  * [9] S. Curzel, S. Jovic, M. Fiorito, A. Tumeo, F. Ferrandi _et al._ , “Higher-Level Synthesis: experimenting with MLIR polyhedral representations for accelerator design,” in _IMPACT 2022: 12th International Workshop on Polyhedral Compilation Techniques_ , 2022, pp. 1–10.
  * [10] C. Lattner, M. Amini, U. Bondhugula, A. Cohen, A. Davis, J. Pienaar, R. Riddle, T. Shpeisman, N. Vasilache, and O. Zinenko, “MLIR: Scaling compiler infrastructure for domain specific computation,” in _2021 IEEE/ACM International Symposium on Code Generation and Optimization (CGO)_. IEEE, 2021, pp. 2–14.
  * [11] A. N. Ziogas, T. Schneider, T. Ben-Nun, A. Calotoiu, T. De Matteis, J. de Fine Licht, L. Lavarini, and T. Hoefler, “Productivity, portability, performance: data-centric Python,” in _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ , 2021, pp. 1–13.
  * [12] S. W. Nabi and W. Vanderbauwhede, “Automatic pipelining and vectorization of scientific code for FPGAs,” _International Journal of Reconfigurable Computing_ , vol. 2019, pp. 1–12, 2019.
  * [13] C. Lattner, “LLVM and Clang: Next generation compiler technology,” in _The BSD conference_ , vol. 5, 2008, pp. 1–20.
  * [14] L. J. Kedward, B. Aradi, O. Čertík, M. Curcic, S. Ehlert, P. Engel, R. Goswami, M. Hirsch, A. Lozada-Blanco, V. Magnin _et al._ , “The state of Fortran,” _Computing in Science & Engineering_, vol. 24, no. 2, pp. 63–72, 2022.
  * [15] L. Forget, G. Harnisch, R. Keryell, and F. de Dinechin, “A single-source C++ 20 HLS flow for function evaluation on FPGA and beyond.” in _International Symposium on Highly-Efficient Accelerators and Reconfigurable Technologies_ , 2022, pp. 51–58.
  * [16] S. McIntosh-Smith, M. Martineau, T. Deakin, G. Pawelczak, W. Gaudin, P. Garrett, W. Liu, R. Smedley-Stevenson, and D. Beckingsale, “TeaLeaf: A mini-application to enable design-space explorations for iterative sparse linear solvers,” in _2017 IEEE International Conference on Cluster Computing (CLUSTER)_. IEEE, 2017, pp. 842–849.
  * [17] A. XILINX, _Vitis High-Level Synthesis User Guide-UG1399 (v2022.2)_. June, 2022.
  * [18] N. Brown, “Porting incompressible flow matrix assembly to FPGAs for accelerating HPC engineering simulations,” in _2021 IEEE/ACM International Workshop on Heterogeneous High-performance Reconfigurable Computing (H2RC)_. IEEE, 2021, pp. 9–20.
  * [19] J. Hrica, “Floating-point design with Vivado HLS,” _Xilinx Application Note_ , 2012.


# Predicting Long-Term Human Behaviors in Discrete Representations via
Physics-Guided Diffusion

Zhitian Zhang1, Anjian Li2, Angelica Lim1, Mo Chen1 1 Simon Fraser University
{zhitianz, angelica<EMAIL_ADDRESS>Princeton University
<EMAIL_ADDRESS>
###### Abstract

Long-term human trajectory prediction is a challenging yet critical task in
robotics and autonomous systems. Prior work that studied how to predict
accurate short-term human trajectories with only unimodal features often
failed in long-term prediction. Reinforcement learning provides a good
solution for learning human long-term behaviors but can suffer from challenges
in data efficiency and optimization. In this work, we propose a long-term
human trajectory forecasting framework that leverages a guided diffusion model
to generate diverse long-term human behaviors in a high-level latent action
space, obtained via a hierarchical action quantization scheme using a VQ-VAE
to discretize continuous trajectories and the available context. The latent
actions are predicted by our guided diffusion model, which uses physics-
inspired guidance at test time to constrain generated multimodal action
distributions. Specifically, we use reachability analysis during the reverse
denoising process to guide the diffusion steps toward physically feasible
latent actions. We evaluate our framework on two publicly available human
trajectory forecasting datasets: SFU-Store-Nav and JRDB, and extensive
experimental results show that our framework achieves superior performance in
long-term human trajectory forecasting.

##  I Introduction

The ability to anticipate future behavior is a remarkable yet common skill for
humans. We can navigate naturally through complex environments and perform
long-term planning. Having an agent that is capable of predicting realistic
and plausible human behaviors is essential to human-robot interaction and
autonomous systems [2, 3]. Such prediction ability could not only benefit the
decision-making of an autonomous system, but also generate high-quality data
that mimics human behaviors. However, creating such an agent has been a
challenging problem in the domain of machine learning and robotics.

Reinforcement learning (RL) has been widely studied in behavior learning, but
can also be problematic due to reward constraints [4, 5], data limitations
[6], high costs of sample complexity [7, 8] and relatively short-horizon
behaviors [9, 10]. Imitation learning (IL) is a type of offline RL [11], where
an agent is trained to mimic the actions from an offline human demonstration
dataset. IL has shown promising performance in autonomous driving [12],
robotics [13], and game-playing [14] on imitating human behaviors. However, IL
faces challenges in producing diverse distributions and long-term predictions
with only mean squared error (MSE) as the optimization objective. In this
work, we seek to model a diverse yet realistic human agent’s behavior
distribution for long-term trajectory prediction. While an agent often faces
an uncountable set of states in the continuous trajectory space, our key
insight is that by tokenizing the continuous trajectory space, generative
models can better capture the multimodality of the distribution while avoiding
the curse of dimensionality. We design a hierarchical action quantization
(HAQ) scheme based on vector quantized variational autoencoder (VQ-VAE) [15]
to learn the mapping between the continuous trajectory space and discrete
latent action space. We then apply imitation learning in the discrete latent
action space.

There are several recent advances in generating human behaviors with discrete
representations [16, 17, 18, 19, 20], mostly based on autoregressive models
using transformers [21]. These autoregressive models have two major drawbacks
in the imitation learning setting. First, the learned policy with a
distribution shift at test time could cause the agent to drift away from
optimal states [22]. Second, the autoregressive model only predicts one token
at a time, making it slow in generating long sequences [23]. Recently,
diffusion models have been studied to imitate human behaviors [24, 25] or
robot actions [26], showing a strong capability of modeling complex action
distributions and a stable training process.

In this paper, we use a denoising diffusion probabilistic model (DDPM) [27]
with the analog bit approach [23] to achieve accurate long-term human
trajectory prediction. Our proposed diffusion model takes states (motion
history, environmental information) as the conditional input and samples
future discrete latent actions, which are learned from our hierarchical VQ-
VAEs. We hypothesize that diffusion models as expressive generative models,
with an efficient discrete representation of continuous trajectory space, can
generate plausible human behaviors in long-horizon. To the best of our
knowledge, this is the first work to generate discrete latent human actions
with diffusion models. Furthermore, we propose reachability guidance to
improve the physical feasibility of sampled trajectories. Our reachability
guidance does not require the training of an extra classifier in the pipeline
which could be cumbersome [28]. We demonstrate that with the proposed physics-
guided diffusion policy, our approach significantly improves the performance
in long-term human trajectory prediction tasks.

To summarize, our contribution is threefold. (1) We propose a simple yet
efficient approach using a hierarchical VQ-VAE to quantize a continuous human
trajectory space to a discrete set of actions. Our discrete representations
capture the multimodality of behaviors and enable the use of IL methods to
generate long-term realistic human behaviors. (2) We present a novel paradigm
where we model future human behavior distributions in discrete action space
using a denoising diffusion probabilistic model (DDPM)-based behavior
generator. (3) We introduce reachability guidance: an intuitive physics-
inspired guidance that incorporates the laws of physics and safety into the
denoising process, allowing our model to generate physically feasible human
behaviors (Figure 1).

##  II Related Work

###  II-A Off Policy Learning

Off-policy learning has proven to be effective in trajectory generation.
Offline imitation learning (IL) [29, 30, 5, 31] has emerged as a popular
method for learning human actions from pre-collected datasets. When the reward
information is available, offline reinforcement learning (RL) [11, 32, 33, 34]
has also been explored for learning human actions. However, IL, as a form of
supervised learning, often assumes a unimodal action distribution [35] that
rarely holds in real-world datasets. To address this issue, positive unlabeled
learning [36] is used to focus learning on the expert part of the dataset,
thereby reducing the modality in the learning process, and behaviour
transformers [19] use transformers with k-means clustering to model multiple
modes of human behaviors. In addition, most existing IL methods suffer from
covariant shift when applied to online trajectory prediction, limiting the
effectiveness in long-term prediction [37, 3, 38, 39]. Online RL has been
integrated with IL to address distribution shift through a closed-loop
learning fashion [40] but it requires extensive online training time.

###  II-B Discretizing continuous action space

Although real-world human actions occur in a continuous space, directly
learning continuous actions from a limited pre-collected dataset often results
in unsatisfactory performance [17]. Therefore, action space discretization has
been investigated while avoiding the exponential growth of the state
dimension. One common approach treats each action dimension as independent
[41, 42, 43, 44], while an alternative approach is sequential discretization
through learned causal dependence [45, 46, 47]. Leveraging human demonstration
for action quantization has also been proven effective in generating a
discrete set of reasonable actions [16]. Recently, vector quantized
variational autoencoders (VQ-VAE) have been adopted for off-policy learning
[41], resulting in a small discrete action space while improving the learning
performance. Furthermore, action discretization plays an important role in
hierarchical learning methods to help identify sub-goals and facilitate action
generation over long horizons [48].

###  II-C Diffusion Models

Diffusion models [49, 50] are a class of generative models that are designed
to sample from complex distributions through a reverse stochastic differential
equation (SDE) process. Later on, the development of denoising diffusion
probabilistic models (DDPM) [27] and denoising diffusion implicit models
(DDIM) [51] showed remarkable efficacy in image generation tasks. To enhance
the conditional input in the generation process, classifier guidance [52] and
classifier-free guidance [28] are introduced. Moreover, diffusion models also
demonstrate their exceptional ability to sample behaviors or trajectories from
multimodal distributions in the field of robotics [53, 54, 26, 55, 24]. With
the attention mechanism first proposed in DDPM [27], diffusion models also
show superior performance in modeling sequential correlation of the data.

##  III Background

###  III-A Problem Formulation

Our objective is to learn a robust policy that could generate long-term
plausible human behaviors given the current and previous states in the last
$H$ time steps. In this paper, we are particularly interested in human
navigational behaviors. We aim to generate a long-term future trajectory for
the target agent over the next $T$ time steps. Formally, this means learning
the future trajectory conditional distribution $p(\mathbb{X}|X,C)$, where
$\mathbb{X}=\\{x^{t+1},y^{t+1},\cdots,x^{t+T},y^{t+T}\\}$ represents the
future trajectory, $X=\\{x^{t-H},y^{t-H},\cdots,x^{t},y^{t}\\}$ the past
trajectory including the current state, and $C=\\{c^{t-H},\cdots,c^{t}\\}$ any
other possible contextual observation features such as information related to
the scene or the other nearby agents. We define the state of the agent at a
given time $t$ to be $s^{t}=(x^{t},y^{t},c^{t})$, with $(x,y)$ representing
position. A detailed time horizon is shown in Figure 2.

As mentioned earlier, we incorporate physics priors during the learning
process to encourage feasible human behavior generation. To facilitate the use
of the proposed control-theoretic methods, we use a dynamically extended
Dubins Car to approximately model the dynamics of a human:

$\begin{gathered}\dot{x}=v\cos\theta,\enspace\dot{y}=v\sin\theta\\\
\dot{v}=u_{1},\enspace\dot{\theta}=u_{2}\end{gathered}$ (1)

where $v$ is the speed of movement and $\theta$ is the orientation of the
human. The state of the above system $(x,y,v,\theta)$ represents the augmented
human state (ignoring the context $c$) used only for the purposes of
incorporating physics priors. We define $\hat{u}=(u_{1},u_{2})$ as the control
variables of our system, where $u_{1}$ is the acceleration, $u_{2}$ is the
angular velocity.111Note that since Equation 1 is differentially flat, all
other variables can be easily computed given $(x^{t},y^{t})$.

###  III-B Imitation Learning (IL)

Imitation learning constructs an optimal policy by mimicking a set of expert
demonstrations, without knowing any reward function. Assume the agents have
access to a dataset of expert demonstrations $\mathcal{D}=$
$\left\\{s^{0},\hat{u}^{0},\cdots,s^{T},\hat{u}^{T}\right\\}$ produced by the
expert policy $\pi_{\beta}$, the goal is to learn a policy $\pi$ that imitates
$\pi_{\beta}$. The simplest approach is through behavioral cloning (BC) where
the policy is trained by maximizing the log-likelihood of actions,
$\mathbb{E}_{s,\hat{u}\sim\mathcal{D}}[\log\pi(\hat{u}|s)]$.

###  III-C Vector quantized variational autoencoder

In this work, we use the VQ-VAE [15] to tokenize continuous human
trajectories. It comprises an encoder that maps observations onto a sequence
of discrete latent variables and a decoder that reconstructs the observations
from the discrete latent variables. The encoder $E$ outputs a continuous
embedding $E(\mathbf{s})$ from the input space
$\mathbf{s}=\\{s^{t-T_{vq}},\cdots,s^{t}\\}$ and discretization is done by
finding the index of the nearest prototype vector in the codebook
$\mathbf{e}_{j}$ to the encoder embedding $E(\mathbf{s})$ based on the
distance, where $j\in\\{1,\cdots,J\\}$. The process can be described as:

$\displaystyle\operatorname{Quantize}(E(\mathbf{s}))=\mathbf{e}_{i}\enspace\text{where}\enspace
i=\underset{j}{\arg\min}\left\|E(\mathbf{s})-\mathbf{e}_{j}\right\|$ (2)

Let $\mathbf{z}$ denote the final quantized latent code for input
$\mathbf{s}$, i.e., $\mathbf{z}=\mathbf{e}_{i}$; the decoder $D$ will
reconstruct the original input $\mathbf{s}$ based on latent code
$\mathbf{z}$.222In practice, a collection of vectors are quantized and decoded
in parallel. In addition to the reconstruction objective, a codebook loss and
a commitment loss are added to move codebook vectors closer to the encoder
embeddings. The overall optimization objective of VQ-VAE can be described as:

$\displaystyle\mathcal{L}_{vq}=\|\mathbf{s}-D(\mathbf{z})\|_{2}^{2}+\|sg[E(\mathbf{s})]-\mathbf{z}\|_{2}^{2}+\beta\|sg[\mathbf{z}]-E(\mathbf{s})\|_{2}^{2}$
(3)

where $sg$ stands for the stop gradient operator and $\beta$ is a hyper-
parameter for the commitment loss. The parameters of the encoder and decoder
are both optimized by this objective.

Figure 2: Illustration of past horizon $H$ and future horizon $T$ in
continuous trajectory space (Top). A discrete action $\mathbf{a}^{\tau}$
tokenizes all states in a period of $T_{vq}$ in continuous trajectory space
(Bottom). Figure 3: Overview of our framework. Hierarchical Action
Quantization (HAQ) encoder learns a discrete representation of human
behaviors. Our diffusion policy generates 6 discrete future actions
conditioned on past observations. During each reverse denoising process,
reachability guidance is used to enforce some physical constraints. The final
output is a long-term future human trajectory reconstructed from discrete
future actions using the HAQ decoder.

##  IV Method

Our goal is to generate diverse and feasible long-term human behaviors. The
overall framework consists of three parts: Hierarchical Action Quantization,
Action Diffusion Policy, and Reachability Guidance, which are illustrated in
Figure 3. We consider a goal-oriented navigational problem where the human
agent’s objective is to navigate in the environment to reach a terminal state.
In our proposed framework, our Action Diffusion Policy models a multimodal and
high-level future action distribution $p(A|S)$ conditioned on low-level past
observations $S=\\{X,C\\}$, where $A$ represents the predicted discrete latent
action sequence in the future. Our Hierarchical Action Quantization maps each
discrete latent action $A$ back to the continuous trajectory $\mathbb{X}$. To
incorporate physics information into the reverse diffusion process, we
introduce the Reachability Guidance. We now explain the details of each
component.

###  IV-A Hierarchical Action Quantization

Naturally, human actions can be represented with discrete representations.
Previous works in human trajectory prediction [3, 56, 57] mostly represent
human trajectories in continuous space. However, we aim to learn a context-
aware discrete representation of human high-level actions by training a VQ-VAE
model.

In contrast to the vanilla VQ-VAE, in this work we propose to formulate a
hierarchical action quantization (HAQ) structure to produce the discrete
latent code, inspired by [58]. The motivation behind this is that we want to
capture the multimodality of human actions by modeling the context information
$c^{t}$ separately from $\\{x^{t},y^{t}\\}$. In order to achieve this, we use
a two-level hierarchical structure. The hierarchy contains a top level encoder
$E_{top}$ that learns the context information
$C_{vq}=\\{c^{t-T_{vq}},\cdots,c^{t}\\}$, i.e. body pose; a bottom level
encoder $E_{top}$ that learns the trajectory information
$X_{vq}=\\{x^{t-T_{vq}},y^{t-T_{vq}},\cdots,x^{t},y^{t}\\}$, and a HAQ decoder
$D$ that reconstructs the original input $X_{vq}$. The HAQ encoding process of
the network can be written as:

$\begin{gathered}\mathbf{h}_{\text{top}}=E_{\text{top}}(C_{vq}),\enspace\mathbf{a}_{\text{top}}=\operatorname{Quantize}(\mathbf{h}_{\text{top}})\\\
\mathbf{h}_{\text{bot}}=E_{\text{bot}}(X_{vq},\mathbf{z}_{\text{top}}),\enspace\mathbf{a}_{\text{bot}}=\operatorname{Quantize}(\mathbf{h}_{\text{bot}})\end{gathered}$
(4)

where $\mathbf{h}$ is the continuous latent variables obtained from the
encoders, $\operatorname{Quantize}$ process is defined as in Equation 2,
$\mathbf{a}_{bot}$ is our final discrete action representation, and later we
refer it as $\mathbf{a}$ for brevity. The overall network is optimized by the
objective defined in Equation 3 for both top and bottom codebooks, with the
reconstruction objective only applying to input $X_{vq}$. After training our
VQ-VAE, we obtain a semantic-rich codebook that is conditioned on contextual
information and we can represent a sequence of low-level continuous states
with a discrete high-level action (see LABEL:fig:overview(Right)).

###  IV-B Action Diffusion Policy

In this section, we introduce our diffusion policy under the imitation
learning formulation. With a learned HAQ encoder, we can represent future
continuous trajectory space $\mathbb{X}$ in discrete space:
$A=\\{\mathbf{a}^{1},\mathbf{a}^{2},\cdots,\mathbf{a}^{\tau}\\}$, where
$\tau=T/T_{vq}$ (See Figure 2). Given previous continuous states $S$, we can
represent a behavior cloning policy as $\pi(A|S)=p(A|S)$ (subsection III-B).
We apply a conditional diffusion model to the discrete latent action space, to
model $p(A|S)$. The behavior cloning policy $\pi(A|S)$ can be then optimized
by the diffusion objective, which aims to sample $A$ from the same
distribution as $D$ [25].

In this work, we adopt denoising diffusion probabilistic models (DDPMs) [27]
as our behavior cloning policy. And to allow our diffusion policy to generate
discrete actions, we use the analog bits approach from the bit diffusion model
[23]. We start with a short introduction to diffusion models. Starting from a
noisy discrete action sequence $A_{k}$, where $A_{K}\sim$
$\mathcal{N}(\mathbf{0},\mathbf{I})$, a sequence of $A_{K-1},\cdots,A_{0}$ is
predicted through $K$ iterations of denoising steps, each with a decreasing
level of noise until the “clean” output $A_{0}$ is formed. During training,
which is also called forward diffusion process, noisy input is generated as
$A_{k}=\sqrt{\bar{\alpha}_{k}}A+\sqrt{1-\bar{\alpha}_{k}}\epsilon$, for some
variance schedule $\bar{\alpha}_{k}$, random noise
$\epsilon\sim\mathcal{N}(\mathbf{0},\mathbf{I})$. A denoising network
$\mathcal{G}_{\theta}$ is trained to predict the noise that was added to the
input, conditioned on some past states $S$, by minimizing the following
objective:

$\displaystyle\mathcal{L}_{\text{DDPM},\theta}:=\mathbb{E}_{S,A,k,\epsilon}\left[\left\|\mathcal{G}_{\theta}\left(S,A_{k},k\right)-\epsilon\right\|_{2}^{2}\right]$
(5)

During the reverse diffusion process, which often refers to the sampling time,
with the variance schedule parameters $\alpha$ and $\sigma$, the “cleaner”
input is generated as follows:

$\displaystyle
A_{k-1}=\frac{1}{\sqrt{\alpha_{k}}}\left(A_{k}-\frac{1-\alpha_{k}}{\sqrt{1-\bar{\alpha}_{k}}}\mathcal{G}_{\theta}\left(S,A,k\right)\right)+\sigma_{k}\epsilon$
(6)

Discrete action sequence prediction. While DDPMs are often used in continuous
space for image generation, we adopt the analog bits [23] method to generate
discrete action sequences with the same continuous diffusion models. The
analog bits approach is twofold: First, during training, discrete data are
represented by bits and then cast into real numbers, which can be directly
modeled by the DDPMs. We denote this process as int2bit, where in our case a
discrete action $\mathbf{a}^{\tau}$ from a codebook of size $J$ (subsection
III-C) can be represented using $n=[\log_{2}I]$ bits, as $\\{0,1\\}^{n}$. Then
during sampling, we draw samples following the same procedures in DDPMs, and
apply a simple thresholding before decoding back into discrete data. We denote
this process as bit2int. The whole process can be described with the
following:

$\displaystyle A_{0}=\textit{int2bit}(\text{A}),\text{forward diffusion}$ (7)
$\displaystyle\text{A}=\textit{bit2int}(A_{k-1}),\text{reverse diffusion}$ (8)

###  IV-C Reachability Guidance

Given that a diffusion policy does not inherently have knowledge on
fundamental physics laws, it is highly possible to generate human behaviors
that are infeasible and clearly disobey the physics world. In this work, we
apply a physical constraint during the diffusion process, which we call
reachability guidance. The motivation is simple: every step of our diffusion
process produces an intermediate sequence of discrete actions $A_{k-1}$, and
by applying the reachability guidance, we could guide the diffusion process
towards generating samples that are physically feasible. The complete
procedure is shown below:

Algorithm 1 Reachability guided diffusion sampling

1:Past states $S$, denoising network $\mathcal{G}_{\theta}$, denoising
timestep $k$.

2:Denoised action $A_{0}$.

3:$A_{K}\sim$ $\mathcal{N}(\mathbf{0},\mathbf{I})$

4:for $k$ from $K$ to 0 do

5: $\mu,\Sigma\leftarrow\mathcal{G}\left(A_{k},k,S\right)$

6: $A^{\prime}_{k}\sim\mathcal{N}(\mu,\Sigma)$ $\triangleright$ Sample
$A^{\prime}_{k}$ without guidance

7: Compute $p_{\gamma}(r|A^{\prime}_{k})$ $\triangleright$ Reachability
guidance

8: $\tilde{\mu}\leftarrow\mu+s\Sigma\nabla_{A_{k}^{\prime}}\log
p_{\gamma}(r|A_{k}^{\prime})$ $\triangleright$ Compute the new mean

9: $A_{k-1}\sim\mathcal{N}(\tilde{\mu},\Sigma)$ $\triangleright$ Sample
$A_{k-1}$ with guidance

10:end for

11:return $A_{0}$

Reachability. Hamilton-Jacobi (HJ) reachability analysis is a formal method
that can verify the performance and safety of a dynamic system [59]. Given the
assumed dynamics of the human in Equation 1, we can compute a Backward
Reachable Set (BRS) based on the discrete action $\mathbf{a}^{\tau}$. The BRS
represents the set of states such that the trajectories that start from this
set can reach the given discrete action $\mathbf{a}^{\tau}$ (See
LABEL:fig:overview(Left)). To calculate the BRS, we first decode the discrete
action $\mathbf{a}^{\tau}$ using HAQ decoder $D$ to obtain the states
$\\{s^{t-T_{vq}},\cdots,s^{t}\\}$ where $(x^{t},y^{t},\theta^{t},v^{t})$ can
be derived from $s^{t}$, then the process can be written as:

$\displaystyle\mathcal{S}=\textbf{BRS}(D(\mathbf{a}^{\tau}))$ (9)

where $\mathcal{S}$ represents the possible positions that the agent could be
located in order to feasibly reach the starting states in $\mathbf{a}^{\tau}$,
and BRS represents the mathematical calculation of the backward reachable set.
For the sake of brevity, we will defer to [59] for readers interested in the
details of reachability. Intuitively, the BRS should cover all the states from
the previous discrete action: $D(\mathbf{a}^{\tau-1})\subseteq\mathcal{S}$
Thus, to determine whether an action is physically feasible, we can formulate
a classification problem where:

$p(\mathbf{a}^{\tau-1})=\begin{cases}1,&\quad
D(\mathbf{a}^{\tau-1})\subseteq\mathcal{S}\\\
0,&\quad\text{otherwise}\end{cases}$ (10)

Then the probability that a sequence of actions
$A=\\{\mathbf{a}^{1},\mathbf{a}^{2},\cdots,\mathbf{a}^{\tau}\\}$ is physically
feasible can be easily calculated as:

$\displaystyle
p_{\gamma}(r|A)=\frac{1}{\tau-1}\sum^{\tau-1}_{\tau^{\prime}=1}p(\mathbf{a}^{\tau^{\prime}})$
(11)

More generally, we can write $p_{\gamma}$333$\gamma$ does not represent
network parameters here. in term of any intermediate latent $A_{k}$,
$p_{\gamma}(r|A_{k})$. For brevity, we leave out the quantize operation as
defined in Equation 7 in the reachability analysis process.

Guidance. Classifier guidance is a useful technique for improving diffusion
models [52, 49, 50]. The process involves training an additional classifier
with class labels on noisy inputs and using a classifier gradient to guide the
diffusion sampling process. Similarly, one can also use the gradient of
reachability “classification” to guide the diffusion process. We believe that
this is a simple approach to force some physical constraint to the network,
and only applied during the diffusion sampling time. Here we show a brief
derivation on how to modify an unconditional reverse diffusion process
$p_{\theta}\left(A_{k}|A_{k+1}\right)$ to condition on the reachability
analysis result $p_{\gamma}\left(r|A_{k}\right)$, where

$\displaystyle
p_{\theta,\gamma}\left(A_{k}|A_{k+1},r\right)=Zp_{\theta}\left(A_{k}|A_{k+1}\right)p_{\gamma}\left(r|A_{k}\right)$
(12)

And $Z$ is a normalizing constant. Recall that:

$\displaystyle p_{\theta}\left(A_{k}|A_{k+1}\right)$
$\displaystyle=\mathcal{N}(\mu,\Sigma)$ (13) $\displaystyle\log
p_{\theta}\left(A_{k}|A_{k+1}\right)$
$\displaystyle=-\frac{1}{2}\left(A_{k}-\mu\right)^{T}\Sigma^{-1}\left(A_{k}-\mu\right)+C$
(14)

We can approximate $\log p_{\gamma}\left(r|A_{k}\right)$ using a Taylor
expansion around $A_{k}=\mu$:

$\displaystyle\log p_{\gamma}\left(r|A_{k}\right)$
$\displaystyle\left.\approx\log
p_{\gamma}\left(r|A_{k}\right)\right|_{A_{k}=\mu}$
$\displaystyle+\left.\left(A_{k}-\mu\right)\nabla_{A_{k}}\log
p_{\gamma}\left(r|A_{k}\right)\right|_{A_{k}=\mu}$ (15)

Letting $g=\left.\nabla_{A_{k}}\log
p_{\gamma}\left(r|A_{k}\right)\right|_{A_{k}=\mu}$, we can derive an
approximation of desired distribution:

$\displaystyle\log\left(p_{\theta}\left(A_{k}|A_{k+1}\right)p_{\gamma}\left(r|A_{k}\right)\right)$
$\displaystyle\approx$
$\displaystyle-\frac{1}{2}\left(A_{k}-\mu\right)^{T}\Sigma^{-1}\left(A_{k}-\mu\right)$
$\displaystyle+\left(A_{k}-\mu\right)g$ (16)

And finally:

$\displaystyle\log p_{\theta,\gamma}\left(A_{k}|A_{k+1},r\right)=\log
p(z),z\sim\mathcal{N}(\mu+\Sigma g,\Sigma)$ (17)

This derivation suggests that the conditional distribution can be approximated
by shifting the mean of unconditional distribution by $\Sigma g$. Following
[52], a scale factor $s$ is also added to the gradient calculation.

Dataset | SSN | JRDB
---|---|---
Methods | ADE/FDE $\downarrow$ | MModality $\uparrow$ | Goal Rate $\uparrow$ | ADE/FDE $\downarrow$ | MModality $\uparrow$
LSTM[60] | 1.19/2.01 | - | 0.39 | 3.71/4.78 | -
LSTM-CNN [60] | 1.03/1.97 | - | 0.45 | 3.51/4.61 | -
CVAE [57] | 0.79/1.14 | 0.19 | 0.67 | 2.12/2.77 | 0.41
Transformer [56] | 0.76/1.32 | 0.27 | 0.71 | 1.89/2.85 | 0.49
Diffusion-BC [24] | 0.84/1.21 | 0.24 | 0.59 | 2.33/3.13 | 0.38
Transformer-discrete | 0.73/1.16 | 0.21 | 0.76 | 1.91/2.81 | 0.33
Ours (MLP) |  0.72/1.19 | 0.23 | 0.75 | 1.77/2.93 | 0.36
Ours (no guidance) | 0.73/1.05 | 0.27 | 0.87 | 1.63/2.61 | 0.43
Ours | 0.68/0.97 | 0.25 | 0.88 | 1.49/2.53 | 0.39
TABLE I: Quantitative comparison of our method and baselines with T = 30. Our model achieves the best or second-best performances on all datasets. The first-best is highlighted by bold, and the second-best is highlighted by underline. | SSN | JRDB
---|---|---
Method | ADE/FDE $\downarrow$ | ADE/FDE $\downarrow$
Transformer + VQVAE | 0.79/1.39 | 1.93/2.83
Ours + VQVAE | 0.77/1.20 | 1.72/2.74
Diffusion-BC | 0.84/1.21 | 2.33/3.13
Diffusion-BC + guidance | 0.75/1.17 | 2.01/2.89
Ours (Full) | 0.68/0.97 | 1.49/2.53
TABLE II: Ablation studies of quantization choice and reachability guidance on SSN and JRDB datasets. | SSN
---|---
# of BRS | ADE/FDE $\downarrow$ | MModality $\uparrow$
1 (v = 1.5) | 0.73/1.03 | 0.27
2 (v = 1, 1.5) | 0.69/1.01 | 0.25
3 (v = 0.5, 1, 1.5) | 0.68/0.97 | 0.25
TABLE III: Ablation study of different number of reachable sets on SSN
dataset.

(a)

(b)

Figure 4: Comparison between different forecasting timesteps (T = 10, 20 ,30)
on SSN and JRDB datasets in terms of ADE. Our model is better at generating
long-term future trajectories, while still maintaining comparable performance
in short-term prediction.

##  V Experiments

We evaluate our methods on human trajectory forecasting on two publicly
available datasets: SFU-Store-Nav (SSN) [61] and JRDB [62] datasets. Both
datasets consist of real human trajectories with associated visual
information. We sample both datasets at 3Hz and split the dataset into
training, validation, and testing sets with proportions of 80%, 5%, and 15%
respectively. Following prior works in human trajectory forecasting [3, 63,
60], we use two error metrics to evaluate the generated human trajectories:

  * •

Average Displacement Error (ADE): Mean $l_{2}$ distance between ground truth
trajectories and generated trajectories.

  * •

Final Displacement Error (FDE): $l_{2}$ distance between ground truth
trajectories and generated trajectories at the end of time horizon $T$.

Since the models generate multimodal output, we report the minimum ADE and FDE
from 20 randomly generated samples instead. In addition to these two commonly
used metrics, we introduced two more metrics that are relatively new to the
domain of human trajectory forecasting:

  * •

Multimodality [18]: Mean $l_{2}$ distance between $N$ pairs of generated
trajectories under same input condition. We set $N=20$.

  * •

Goal Rate ($M$): Proportions of generated trajectories which reach at least
$M$ number of goals for the entire sequence. We let $M=1$. (Only applicable to
SSN dataset.)

Note that we perform all evaluations on low-level representations. We always
decode the generated discrete actions into continuous human trajectories and
run evaluations.

###  V-A Implementation Details

For our action quantization, the codebook size is set to 256 $\times$ 128,
where the number of discrete actions is 256, and the dimension of each action
token is 128. We set $\beta=1$ and $T_{vq}=5$ for training. We use multilayer
perceptrons (MLPs) as our HAQ encoders and decoder. The VQ-VAE network is
optimized using AdamW [64] with a learning rate of 1e-6 and batch size of 128.
For both SSN and JRDB datasets, visual images are extracted into 2D body pose
features with dimensions of 50. Our discrete diffusion model has two variants:
diffusion-MLP and diffusion-Transformer [24]. The diffusion-MLP is an MLP
model with 3 hidden layers. And the diffusion-Transformer is a standard
transformer [21] model with 2 encoder blocks and multi-head attentions. We use
10 diffusion steps and square cosine noise scheduler [65]. The diffusion
network is optimized with AdamW [64] with a learning rate of 1e-4. We
calculate a set of backward reachable sets (BRS) based on different maximum
traveling speed assumptions. We assume a maximum turn rate of 1 rad/s, and a
maximum acceleration of 0.5 $\text{m/s}^{2}$. The maximum speed is set to be
{0.5, 1, 1.5} m/s for SSN dataset, and {1, 2, 3} m/s for JRDB dataset.
Reachable set calculation is done with the helperOC Matlab toolbox
444https://github.com/HJReachability/helperOC. Each learned discrete action
represents 5 timesteps (1.5s). For evaluation, the observation length is set
as 10 timesteps (3s) and we are predicting 6 future discrete actions, which is
30 timesteps (9s).

###  V-B Results

We seek to answer the following questions in our evaluation.

How does our model compare to prior work in long-term human trajectory
forecasting? In Table I, we compare our method to a large number of commonly
used baselines. Our methods achieve the best ADE and FDE while maintaining a
good level of diversity on both datasets. To further investigate our method’s
ability in long-term forecasting, we evaluate it across different forecasting
horizons. In Figure 4, we evaluate from $T=10$ (short-term) to $T=30$ (long-
term) and observe that as the forecasting horizon increases, our model
significantly outperforms all other baselines. This result demonstrates that
our model has superior performance in long-term human trajectory forecasting
while still maintaining comparable performance in the short-term.

How does diffusion policy compare to the autoregressive model? We are curious
how a diffusion model would perform against the autoregressive model when the
inputs are discrete actions. Transformer-discrete is a variant of Transformer
[56] which is trained with the same quantized actions as our method. As can be
seen in Table I, the variant Transformer-discrete outperforms the original
model Transformer, suggesting the advantage of using discrete actions. Both
our model and MLP-variant outperform the Transformer-discrete with lower
ADE/FDE and higher Multimodality score, which indicates that the use of the
discrete diffusion model is more effective in modeling a multimodal action
distribution in discrete space.

How effective is the reachability guidance to the performance of our model? To
understand how reachability guidance influences the performance of our model,
we implement a variant of our model which does not use reachability guidance.
In Table I, both ADE and FDE increase with the reachability guidance, but the
multimodality score decreases. This suggests that there is a trade-off between
physical feasibility and multimodality. We hypothesize that more diverse
multimodal actions could lead to some states that are not physically feasible;
our reachability guidance places a constraint on those states.

###  V-C Ablation Study

We further perform an extensive ablation study to understand the contributions
of individual components to our model performance. In Table II, we investigate
(1) the effect of our proposed hierarchical action quantization and (2) the
performance of reachability guidance with the continuous diffusion model.
First, we replace the hierarchical action quantization component with a simple
VQ-VAE on Transformer-discrete and our method. We could see a significant
decrease in performance on both datasets, indicating the importance of our
hierarchical action quantization scheme. Then we integrate our reachability
guidance with the continuous diffusion baseline: Diffusion-BC, resulting in a
12% performance increase. This highlights the flexibility and importance of
the reachability guidance. Finally, our second ablation study Table III
investigates the effect of the backward reachable set. As one would expect,
calculating a more restrictive condition for the reachable sets can result in
a slight improvement in terms of ADE and FDE with a little sacrifice in
diversity. When there is only one BRS calculated with a very relaxed
assumption, the model has a similar performance as the one without
reachability guidance. Thus, a more carefully designed set of dynamic
assumptions could maximize the performance of the reachability guidance.

##  VI CONCLUSIONS

In this work, we present our discrete diffusion framework that generates
future human behaviors with physics-inspired guidance. Our goal is for robots
to be able to imitate realistic and diverse human behaviors in the long term.
To achieve this, we propose reachability guidance to enforce physical
constraints during the diffusion process in discrete action space. The
proposed reachability guidance can be used on any diffusion model without
retraining. Experimental results on human trajectory forecasting datasets
demonstrate the superior performance of our framework. The limitations of our
framework include: the dynamics of humans, and assumptions for reachable set
calculation. Recent developments in the diffusion models [66, 67] have shown
promising results incorporating physics in the sampling process, we hope the
proposed framework and reachability guidance could open up new directions for
future work. We believe our framework is robust to other tasks, e.g.,
autonomous driving, motion generation, etc, which are yet to be explored.

## References

  * [1] B. Biro, Z. Zhang, M. Chen, and A. Lim, “The sfu-store-nav 3d virtual human platform for human-aware robotics,” in _ICRA, Workshop on Long-term Human Motion Prediction_ , 2021.
  * [2] R. Mirsky, I. Carlucho, A. Rahman, E. Fosong, W. Macke, M. Sridharan, P. Stone, and S. V. Albrecht, “A survey of ad hoc teamwork: Definitions, methods, and open problems,” in _EUMAS_ , 2022.
  * [3] T. Salzmann, B. Ivanovic, P. Chakravarty, and M. Pavone, “Trajectron++: Dynamically-feasible trajectory forecasting with heterogeneous data,” in _ECCV 2020_ , 2020.
  * [4] X. Pan, Y. You, Z. Wang, and C. Lu, “Virtual to real reinforcement learning for autonomous driving,” _arXiv preprint arXiv:1704.03952_ , 2017.
  * [5] Z. Zhang, A. Liniger, D. Dai, F. Yu, and L. Van Gool, “End-to-end urban driving by imitating a reinforcement learning coach,” in _CVPR_ , 2021.
  * [6] D. Yarats, R. Fergus, A. Lazaric, and L. Pinto, “Mastering visual continuous control: Improved data-augmented reinforcement learning,” _arXiv preprint arXiv:2107.09645_ , 2021.
  * [7] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski _et al._ , “Human-level control through deep reinforcement learning,” _nature_ , 2015.
  * [8] Y. Duan, X. Chen, R. Houthooft, J. Schulman, and P. Abbeel, “Benchmarking deep reinforcement learning for continuous control,” in _ICML_ , 2016.
  * [9] S. Gu, E. Holly, T. Lillicrap, and S. Levine, “Deep reinforcement learning for robotic manipulation with asynchronous off-policy updates,” in _ICRA_ , 2017.
  * [10] V. Kumar, E. Todorov, and S. Levine, “Optimal control with learned local models: Application to dexterous manipulation,” in _ICRA_ , 2016.
  * [11] S. Levine, A. Kumar, G. Tucker, and J. Fu, “Offline reinforcement learning: Tutorial, review, and perspectives on open problems,” _arXiv preprint arXiv:2005.01643_ , 2020.
  * [12] D. Chen and P. Krähenbühl, “Learning from all vehicles,” in _CVPR_ , 2022.
  * [13] P. Florence, C. Lynch, A. Zeng, O. A. Ramirez, A. Wahid, L. Downs, A. Wong, J. Lee, I. Mordatch, and J. Tompson, “Implicit behavioral cloning,” in _CoRL_ , 2022.
  * [14] T. Pearce and J. Zhu, “Counter-strike deathmatch with large-scale behavioural cloning,” in _CoG_ , 2022.
  * [15] A. Van Den Oord, O. Vinyals _et al._ , “Neural discrete representation learning,” _Advances in neural information processing systems_ , 2017.
  * [16] R. Dadashi, L. Hussenot, D. Vincent, S. Girgin, A. Raichuk, M. Geist, and O. Pietquin, “Continuous control with action quantization from demonstrations,” _arXiv preprint arXiv:2110.10149_ , 2021.
  * [17] J. Luo, P. Dong, J. Wu, A. Kumar, X. Geng, and S. Levine, “Action-quantized offline reinforcement learning for robotic skill learning,” in _CoRL_ , 2023.
  * [18] J. Zhang, Y. Zhang, X. Cun, S. Huang, Y. Zhang, H. Zhao, H. Lu, and X. Shen, “T2m-gpt: Generating human motion from textual descriptions with discrete representations,” _arXiv preprint arXiv:2301.06052_ , 2023.
  * [19] N. M. Shafiullah, Z. Cui, A. A. Altanzaya, and L. Pinto, “Behavior transformers: Cloning $k$ modes with one stone,” _NeurIPS_ , 2022.
  * [20] T. Lucas, F. Baradel, P. Weinzaepfel, and G. Rogez, “Posegpt: Quantization-based 3d human motion generation and forecasting,” in _ECCV_ , 2022.
  * [21] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” _NeurIPS_ , 2017.
  * [22] S. Reddy, A. D. Dragan, and S. Levine, “Sqil: Imitation learning via reinforcement learning with sparse rewards,” _arXiv preprint arXiv:1905.11108_ , 2019.
  * [23] T. Chen, R. Zhang, and G. Hinton, “Analog bits: Generating discrete data using diffusion models with self-conditioning,” _arXiv preprint arXiv:2208.04202_ , 2022.
  * [24] T. Pearce, T. Rashid, A. Kanervisto, D. Bignell, M. Sun, R. Georgescu, S. V. Macua, S. Z. Tan, I. Momennejad, K. Hofmann _et al._ , “Imitating human behaviour with diffusion models,” _arXiv preprint arXiv:2301.10677_ , 2023.
  * [25] Z. Wang, J. J. Hunt, and M. Zhou, “Diffusion policies as an expressive policy class for offline reinforcement learning,” _arXiv preprint arXiv:2208.06193_ , 2022.
  * [26] C. Chi, S. Feng, Y. Du, Z. Xu, E. Cousineau, B. Burchfiel, and S. Song, “Diffusion policy: Visuomotor policy learning via action diffusion,” _arXiv preprint arXiv:2303.04137_ , 2023.
  * [27] J. Ho, A. Jain, and P. Abbeel, “Denoising diffusion probabilistic models,” _NeurIPS_ , 2020.
  * [28] J. Ho and T. Salimans, “Classifier-free diffusion guidance,” _arXiv preprint arXiv:2207.12598_ , 2022.
  * [29] M. Bojarski, D. Del Testa, D. Dworakowski, B. Firner, B. Flepp, P. Goyal, L. D. Jackel, M. Monfort, U. Muller, J. Zhang _et al._ , “End to end learning for self-driving cars,” _arXiv preprint arXiv:1604.07316_ , 2016.
  * [30] F. Codevilla, M. Müller, A. López, V. Koltun, and A. Dosovitskiy, “End-to-end driving via conditional imitation learning,” in _ICRA_ , 2018.
  * [31] N. Nayakanti, R. Al-Rfou, A. Zhou, K. Goel, K. S. Refaat, and B. Sapp, “Wayformer: Motion forecasting via simple & efficient attention networks,” in _ICRA_ , 2023.
  * [32] J. Fu, A. Kumar, O. Nachum, G. Tucker, and S. Levine, “D4rl: Datasets for deep data-driven reinforcement learning,” _arXiv preprint arXiv:2004.07219_ , 2020.
  * [33] A. Kumar, A. Zhou, G. Tucker, and S. Levine, “Conservative q-learning for offline reinforcement learning,” _NeurIPS_ , 2020.
  * [34] Y. Wu, G. Tucker, and O. Nachum, “Behavior regularized offline reinforcement learning,” _arXiv preprint arXiv:1911.11361_ , 2019.
  * [35] S. Legine, “Supervised learning of behaviors,” https://rail.eecs.berkeley.edu/deeprlcourse-fa21/static/slides/lec-2.pdf.
  * [36] Q. Wang, R. McCarthy, D. C. Bulens, K. McGuinness, N. E. O’Connor, F. R. Sanchez, N. Gürtler, F. Widmaier, and S. J. Redmond, “Improving behavioural cloning with positive unlabeled learning,” in _CoRL_ , 2023.
  * [37] Y. Chai, B. Sapp, M. Bansal, and D. Anguelov, “Multipath: Multiple probabilistic anchor trajectory hypotheses for behavior prediction,” _arXiv preprint arXiv:1910.05449_ , 2019.
  * [38] N. Rhinehart, R. McAllister, K. Kitani, and S. Levine, “Precog: Prediction conditioned on goals in visual multi-agent settings,” in _CVPR_ , 2019.
  * [39] J. Ngiam, B. Caine, V. Vasudevan, Z. Zhang, H.-T. L. Chiang, J. Ling, R. Roelofs, A. Bewley, C. Liu, A. Venugopal _et al._ , “Scene transformer: A unified multi-task model for behavior prediction and planning,” _arXiv preprint arXiv:2106.08417_ , 2021.
  * [40] Y. Lu, J. Fu, G. Tucker, X. Pan, E. Bronstein, R. Roelofs, B. Sapp, B. White, A. Faust, S. Whiteson _et al._ , “Imitation is not enough: Robustifying imitation with reinforcement learning for challenging driving scenarios,” in _IROS_ , 2023.
  * [41] A. Tavakoli, F. Pardo, and P. Kormushev, “Action branching architectures for deep reinforcement learning,” in _AAAI_ , 2018.
  * [42] Y. Tang and S. Agrawal, “Discretizing continuous action space for on-policy optimization,” in _AAAI_ , 2020.
  * [43] Y. Chebotar, Q. Vuong, K. Hausman, F. Xia, Y. Lu, A. Irpan, A. Kumar, T. Yu, A. Herzog, K. Pertsch _et al._ , “Q-transformer: Scalable offline reinforcement learning via autoregressive q-functions,” in _CoRL_ , 2023.
  * [44] N. Vieillard, M. Andrychowicz, A. Raichuk, O. Pietquin, and M. Geist, “Implicitly regularized rl with implicit q-values,” _arXiv preprint arXiv:2108.07041_ , 2021.
  * [45] L. Metz, J. Ibarz, N. Jaitly, and J. Davidson, “Discrete sequential prediction of continuous actions for deep rl,” _arXiv preprint arXiv:1705.05035_ , 2017.
  * [46] C. Tessler, G. Tennenholtz, and S. Mannor, “Distributional policy optimization: An alternative approach for continuous control,” _NeurIPS_ , 2019.
  * [47] A. Tavakoli, M. Fatemi, and P. Kormushev, “Learning to represent action values as a hypergraph on the action vertices,” _arXiv preprint arXiv:2010.14680_ , 2020.
  * [48] K. Kujanpää, J. Pajarinen, and A. Ilin, “Hierarchical imitation learning with vector quantized models,” _arXiv preprint arXiv:2301.12962_ , 2023.
  * [49] J. Sohl-Dickstein, E. Weiss, N. Maheswaranathan, and S. Ganguli, “Deep unsupervised learning using nonequilibrium thermodynamics,” in _ICML_ , 2015.
  * [50] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole, “Score-based generative modeling through stochastic differential equations,” _arXiv preprint arXiv:2011.13456_ , 2020.
  * [51] J. Song, C. Meng, and S. Ermon, “Denoising diffusion implicit models,” _arXiv preprint arXiv:2010.02502_ , 2020.
  * [52] P. Dhariwal and A. Nichol, “Diffusion models beat gans on image synthesis,” _NeurIPS_ , 2021.
  * [53] A. Ajay, Y. Du, A. Gupta, J. Tenenbaum, T. Jaakkola, and P. Agrawal, “Is conditional generative modeling all you need for decision-making?” _arXiv preprint arXiv:2211.15657_ , 2022.
  * [54] M. Janner, Y. Du, J. B. Tenenbaum, and S. Levine, “Planning with diffusion for flexible behavior synthesis,” _arXiv preprint arXiv:2205.09991_ , 2022.
  * [55] A. Sridhar, D. Shah, C. Glossop, and S. Levine, “Nomad: Goal masked diffusion policies for navigation and exploration,” _arXiv preprint arXiv:2310.07896_ , 2023.
  * [56] Y. Yuan, X. Weng, Y. Ou, and K. M. Kitani, “Agentformer: Agent-aware transformers for socio-temporal multi-agent forecasting,” in _CVPR_ , 2021.
  * [57] K. Mangalam, H. Girase, S. Agarwal, K.-H. Lee, E. Adeli, J. Malik, and A. Gaidon, “It is not the journey but the destination: Endpoint conditioned trajectory prediction,” in _ECCV_ , 2020.
  * [58] A. Razavi, A. Van den Oord, and O. Vinyals, “Generating diverse high-fidelity images with vq-vae-2,” _Advances in neural information processing systems_ , 2019.
  * [59] S. Bansal, M. Chen, S. Herbert, and C. J. Tomlin, “Hamilton-jacobi reachability: A brief overview and recent advances,” in _CDC_ , 2017.
  * [60] Z. Zhang, J. Rhim, A. Lim, and M. Chen, “A multimodal and hybrid framework for human navigational intent inference,” in _IROS_ , 2021.
  * [61] Z. Zhang, J. Rhim, M. TaherAhmadi, K. Yang, A. Lim, and M. Chen, “Sfu-store-nav: A multimodal dataset for indoor human navigation,” _Data in Brief_ , 2020.
  * [62] E. Vendrow, D. T. Le, J. Cai, and H. Rezatofighi, “Jrdb-pose: A large-scale dataset for multi-person pose estimation and tracking,” in _CVPR_ , 2023.
  * [63] A. Gupta, J. Johnson, L. Fei-Fei, S. Savarese, and A. Alahi, “Social gan: Socially acceptable trajectories with generative adversarial networks,” in _CVPR_ , 2018.
  * [64] I. Loshchilov and F. Hutter, “Decoupled weight decay regularization,” _arXiv preprint arXiv:1711.05101_ , 2017.
  * [65] A. Q. Nichol and P. Dhariwal, “Improved denoising diffusion probabilistic models,” in _ICML_ , 2021.
  * [66] D. Rempe, Z. Luo, X. Bin Peng, Y. Yuan, K. Kitani, K. Kreis, S. Fidler, and O. Litany, “Trace and pace: Controllable pedestrian animation via guided trajectory diffusion,” in _CVPR_ , 2023.
  * [67] Y. Yuan, J. Song, U. Iqbal, A. Vahdat, and J. Kautz, “Physdiff: Physics-guided human motion diffusion model,” in _CVPR_ , 2023.