# Heterogeneous Computing on Mixed Unstructured Grids with PyFR

F. D. Witherden111Corresponding author; e-mail
<EMAIL_ADDRESS>B. C. Vermeire, P. E. Vincent

Department of Aeronautics, Imperial College London, SW7 2AZ

###### Abstract

PyFR is an open-source high-order accurate computational fluid dynamics solver
for mixed unstructured grids that can target a range of hardware platforms
from a single codebase. In this paper we demonstrate the ability of PyFR to
perform high-order accurate unsteady simulations of flow on mixed unstructured
grids using heterogeneous multi-node hardware. Specifically, after
benchmarking single-node performance for various platforms, PyFR v0.2.2 is
used to undertake simulations of unsteady flow over a circular cylinder at
Reynolds number $3\,900$ using a mixed unstructured grid of prismatic and
tetrahedral elements on a desktop workstation containing an Intel Xeon E5-2697
v2 CPU, an NVIDIA Tesla K40c GPU, and an AMD FirePro W9100 GPU. Both the
performance and accuracy of PyFR are assessed. PyFR v0.2.2 is freely available
under a 3-Clause New Style BSD license (see www.pyfr.org).

##  1 Introduction

There is an increasing desire amongst industrial practitioners of
computational fluid dynamics (CFD) to perform scale-resolving simulations of
unsteady compressible flows within the vicinity of complex geometries.
However, current generation industry-standard CFD software—predominantly based
on first- or second-order accurate Reynolds Averaged Navier-Stokes (RANS)
technology—is not well suited to the task. Henceforth, over the past decade
there has been significant interest in the application of high-order methods
for mixed unstructured grids to such problems. Popular examples of high-order
schemes for mixed unstructured grids include the discontinuous Galerkin (DG)
method, first introduced by Reed and Hill [1], and the spectral difference
(SD) methods originally proposed under the moniker ‘staggered-gird Chebyshev
multidomain methods’ by Kopriva and Kolias in 1996 [2] and later popularised
by Sun et al. [3]. In 2007 Huynh [4] proposed the flux reconstruction (FR)
approach; a unifying framework for high-order schemes on unstructured grids
that incorporates both the nodal DG schemes of [5] and, at least for a linear
flux function, any SD scheme. In addition to offering high-order accuracy on
unstructured mixed grids, FR schemes are also compact in space, and thus when
combined with explicit time marching offer a significant degree of element
locality. This locality makes such schemes extremely good candidates for
acceleration using either the vector units of modern CPUs or graphics
processing units (GPUs) [6, 7, 8]. There exist a variety of approaches for
writing accelerated codes. These include directive based methodologies such as
OpenMP 4.0 and OpenACC, and language frameworks such as OpenCL and CUDA.
Unfortunately, at the time of writing, there exists no single approach which
is _performance portable_ across all major hardware platforms. Codes desiring
cross-platform portability must therefore incorporate support for multiple
approaches. Further, there is also a growing interest from the scientific
community in _heterogeneous computing_ whereby multiple platforms are employed
simultaneously to solve a problem. The promise of heterogeneous computing is
improved resource utilisation on systems with a mix of hardware. Such systems
are becoming increasingly common.

PyFR is a high-order FR code for solving the Euler and compressible Navier-
Stokes equations on mixed unstructured grids [8]. Written in the Python
programming language PyFR incorporates backends for C/OpenMP, CUDA, and
OpenCL. It is therefore capable of running on conventional CPUs, as well as
GPUs from both NVIDIA and AMD, as well as heterogeneous mixtures of the
aforementioned. The objective of this paper is to demonstrate the ability of
PyFR to perform high-order accurate unsteady simulations of flow on mixed
unstructured meshes using a heterogeneous hardware platform—demonstrating the
concept of ‘heterogeneous computing from a homogeneous codebase’.
Specifically, we will undertake simulations of unsteady flow over a cylinder
at Reynolds number $3\,900$ using a mixed unstructured grid of prismatic and
tetrahedral elements using a desktop workstation containing an Intel Xeon
E5-2697 v2 CPU, an NVIDIA Tesla K40c GPU, and an AMD FirePro W9100 GPU. At the
time of writing these represent the high-end workstation offerings from the
three major hardware vendors. All are designed for workstations with support
for error-correcting code (ECC) memory and double precision floating point
arithmetic.

The paper is structured as follows. In section 2 we provide a brief overview
of the PyFR codebase. In section section 3 we provide details of the cylinder
test case. In section 4 single node performance is discussed. In section 5
multi-node heterogeneous performance is discussed, and finally in section 6
conclusions are drawn.

##  2 PyFR

###  2.1 Overview

For a detailed overview of PyFR the reader is referred to Witherden et al.
[8]. Key functionality of PyFR v0.2.2 is summarised in Table 1. We note that
PyFR achieves platform portability via use of an innovative ‘backend’
infrastructure.

Table 1: Key functionality of PyFR v0.2.2. Dimensions | 2D, 3D
---|---
Elements | Triangles, Quadrilaterals, Hexahedra,
| Tetrahedra, Prisms
Spatial orders | Arbitrary
Time steppers | Euler, RK4, RK45
Precisions | Single, Double
Platforms | CPUs, NVIDIA GPUs, AMD GPUs
Inter-Node Communication | MPI
Governing Systems | Euler, Compressible Navier-Stokes

The majority of operations within an FR time-step can be cast as matrix
multiplications of the form

$\bm{\mathsf{C}}\leftarrow
c_{1}\bm{\mathsf{A}}\bm{\mathsf{B}}+c_{2}\bm{\mathsf{C}},$ (1)

where $c_{1,2}\in\mathbb{R}$ are scalar constants, $\bm{\mathsf{A}}$ is a
constant operator matrix, and $\bm{\mathsf{B}}$ and $\bm{\mathsf{C}}$ are row-
major state matrices. Within the taxonomy proposed by Goto et al. [9] the
multiplications fall into the block-by-panel (GEBP) category. The specific
dimensions of the operator matrices are a function of both the polynomial
order $\wp$ used to represent the solution in each element of the domain, and
the element type. A breakdown of these dimensions for the volume-to-surface
interpolation operator matrix $\bm{\mathsf{M}}^{0}$ can be found in Table 2.
In PyFR platform portability of the matrix multiply operations is achieved by
deferring to the GEMM family of subroutines provided by a Basic Linear Algebra
Subprograms (BLAS) library for the target platform.

Table 2: Dimensions of the volume-to-surface interpolation operator matrix $\bm{\mathsf{M}}^{0}$ at orders $\wp=1,2,3,4$ for tetrahedral, prismatic, and hexahedral element types. | $\dim{\bm{\mathsf{M}}^{0}}$
---|---
Type | $\wp=1$ | $\wp=2$ | $\wp=3$ | $\wp=4$
Tet | $4\times 12$ | $10\times 24$ | $20\times 40$ | $\hphantom{0}35\times 60$
Pri | $6\times 18$ | $18\times 39$ | $40\times 68$ | $\hphantom{0}75\times 105$
Hex | $8\times 24$ | $27\times 54$ | $64\times 96$ | $125\times 150$

All other operations involved in an FR time-step are point-wise, concerning
themselves solely with data at an individual solution/flux point. In PyFR
platform portability of these point-wise operations is achieved via use of a
bespoke domain specific language based on the Mako templating engine [10].
Mako specifications of point-wise operations are converted into backend-
specific low-level code for the target platform at runtime, which is then
compiled, linked and loaded into PyFR.

###  2.2 C/OpenMP backend

The C/OpenMP backend can be used to target CPUs from a range of vendors. The
BLAS implementation employed by PyFR for the C/OpenMP backend must be
specified by the user at runtime. Both single- and multi-threaded libraries
are supported. When a single-threaded library is specified PyFR will perform
its own parallelisation. Given a state matrix $\bm{\mathsf{B}}$ of dimension
$(K,N)$ the decomposition algorithm splits $\bm{\mathsf{B}}$ into $N_{t}$
slices of dimension $(K,N/N_{t})$ where $N_{t}$ is the number of OpenMP
threads. Each thread then multiplies its slice through by $\bm{\mathsf{A}}$ to
yield the corresponding slice of $\bm{\mathsf{C}}$. A visualisation of this
approach is shown in Figure 1. For the block-by-panel multiplications
encountered in FR this strategy has been found to consistently outperform
those employed by multi-threaded BLAS libraries.

$\underbrace{\left[\vphantom{\int}\;\;\hphantom{\ldots}\;\;\middle|\;\;\hphantom{\ldots}\;\;\middle|\;\;\ldots\;\;\middle|\;\;\hphantom{\ldots}\;\;\middle|\;\;\hphantom{\ldots}\;\;\right]}_{\bm{\mathsf{C}}}=\underbrace{\left[\vphantom{\int}\;\;\hphantom{\ldots}\;\;\right]}_{\bm{\mathsf{A}}}\underbrace{\left[\vphantom{\int}\;\;\hphantom{\ldots}\;\;\middle|\;\;\hphantom{\ldots}\;\;\middle|\;\;\ldots\;\;\middle|\;\;\hphantom{\ldots}\;\;\middle|\;\;\hphantom{\ldots}\;\;\right]}_{\bm{\mathsf{B}}}$
Figure 1: How matrix multiplications are parallelised in the C/OpenMP backend
of PyFR.

###  2.3 CUDA backend

The CUDA backend can be used to target NVIDIA GPUs of compute capability 2.0
or later. PyCUDA [11] is used to invoke the CUDA API from Python. Matrix-
multiplications are performed using the cuBLAS library which ships as part of
the CUDA distribution. The cuBLAS library is exclusively column-major.
Nevertheless it is possible to directly multiply two row-major matrices
together by utilising the fact that

$\bm{\mathsf{C}}=\bm{\mathsf{A}}\bm{\mathsf{B}}\implies\bm{\mathsf{C}}^{T}=(\bm{\mathsf{A}}\bm{\mathsf{B}})^{T}=\bm{\mathsf{B}}^{T}\bm{\mathsf{A}}^{T},$
(2)

and observing the effect of passing a row-major matrix to a column-major
subroutine is to implicitly transpose it.

###  2.4 OpenCL backend

The OpenCL backend can be used to target CPUs and GPUs from a range of
vendors. The PyOpenCL package [11] is used to interface OpenCL with Python.
OpenCL natively supports runtime code generation. BLAS support is provided via
the open source clBLAS library, which is primarily developed and supported by
AMD. For GPU devices clBLAS utilises auto-tuning in order to effectively
target a wide range of architectures. Performance is heavily dependent on the
underlying OpenCL implementation.

###  2.5 Distributed memory systems

To scale out across multiple nodes PyFR utilises the Message Passing Interface
(MPI). By construction the data exchanged between MPI ranks is independent of
the backend being used. A natural consequence of this is that different MPI
ranks can transparently utilise different backends—hence enabling PyFR to run
on heterogenous multi-node systems.

##  3 Cylinder Test Case

###  3.1 Overview

Flow over a circular cylinder has been the focus of numerous previous
experimental and numerical studies. Characteristics of the flow are known to
be highly dependent on the Reynolds number $Re$, defined as

$Re=\frac{u_{\infty}D}{\nu},$ (3)

where $u_{\infty}$ is the free-stream fluid speed, $D$ is the cylinder
diameter, and $\nu$ is the fluid kinematic viscosity. Roshko [12] identified a
stable range between $Re=40$ and $150$ that is characterised by the shedding
of regular laminar vortices, as well as a transitional range between $Re=150$
and $300$, and a turbulent range beyond $Re=300$. These results were
subsequently confirmed by Bloor [13], who identified a similar set of regimes.
Later, Williamson [14] identified two modes of transition from two dimensional
to three dimensional flow. The first, known as Mode-A instability, occurs at
$Re\approx 190$ and the second, known as Mode-B instability, occurs at
$Re\approx 260$. The turbulent range beyond $Re=300$ can be further sub-
classified into the shear-layer transition, critical, and supercritical
regimes as discussed in the review by Williamson [15].

In the present study we consider flow over a circular cylinder at $Re=3\,900$,
and an effectively incompressible Mach number of $0.2$. This case sits in the
shear-layer transition regime identified by Williamson [15], and contains
several complex flow features, including separated shear layers, turbulent
transition, and a fully turbulent wake.

###  3.2 Domain

In the present study we use a computational domain with dimensions
$[-9D,25D]$; $[-9D,9D]$; and $[0,\pi D]$ in the stream-, cross-, and span-wise
directions, respectively. The cylinder is centred at $(0,0,0)$. The span-wise
extent was chosen based on the results of Norberg [18], who found no
significant influence on statistical data when the span-wise dimension was
doubled from $\pi D$ to $2\pi D$. Indeed, a span of $\pi D$ has been used in
the majority of previous numerical studies [16, 19, 20], including the recent
DNS study of Lehmkuhl et al. [21]. The stream-wise and cross-wise dimensions
are comparable to the experimental and numerical values used by Parnaudeau et
al. [22], whose results will be directly compared with ours. The overall
domain dimensions are also comparable to those used for DNS studies by
Lehmkuhl et al. [21]. The domain is periodic in the span-wise direction, with
a no-slip isothermal wall boundary condition applied at the surface of the
cylinder, and Riemann invariant boundary conditions applied at the far-field.

###  3.3 Mesh

The domain was meshed in two ways. The first mesh consisted of entirely
hexahedral elements, and the second a mixture of prismatic elements in the
near wall boundary layer region, and tetrahedral elements in the wake and far-
field. Both meshes employed quadratically curved elements, and were designed
to fully resolve the near wall boundary layer region when $\wp=4$.
Specifically, the maximum skin friction coefficient was estimated _a priori_
as $C_{f}\approx 0.075$ based on the LES results of Breuer [19]. The height of
the first element was then specified such that when $\wp=4$ the first solution
point from the wall sits at $y+\approx 1$, where non-dimensional wall units
are calculated in the usual fashion as $y+=u_{\tau}y/\nu$ with
$u_{\tau}=\sqrt{C_{f}/2}u_{\infty}$.

The hexahedral mesh had 200 elements in the circumferential direction, and 10
elements in the span-wise direction, which when $\wp=4$ achieves span-wise
resolution comparable to that used in previous studies; as discussed by Breuer
[19] and the references contained therein. The prism/tetrahedral mesh has 120
elements in the circumferential direction, and 20 elements in the span-wise
direction, these numbers being chosen to help reduce face aspect ratios at the
edges of the prismatic layer; which facilitates transition to the fully
unstructured tetrahedral elements in the far-field. In total the hexahedral
mesh contained $118\,070$ elements, and the prism/tetrahedral mesh contained
$79\,344$ prismatic elements and $227\,298$ tetrahedral elements. Both meshes
are shown in Figure 2.

(a) Hexahedral, far-field.

(b) Prism/tetrahedral, far-field.

(c) Hexahedral, wake.

(d) Prism/tetrahedral, wake.

Figure 2: Cutaways through the two meshes.

###  3.4 Methodology

The compressible Navier-Stokes equations, with constant viscosity, were solved
on each of the two meshes shown in Figure 2. A DG scheme was used for the
spatial discretisation, a Rusanov Riemann solver was used to calculate the
inviscid fluxes at element interfaces, and the explicit RK45[2R+] scheme of
Carpenter and Kennedy [23] was used to advance the solution in time. No sub-
grid model was employed, hence the approach should be considered ILES/DNS, as
opposed to classical LES. Values of $\wp$ from one to four were used with each
mesh. The approximate memory requirements of PyFR for these simulations can
are detailed in Table 3. The total number of floating point operations
required by an RK45[2R+] time-step for each simulation are detailed in Figure
3. When running with $\wp=1$ both meshes require ${\sim}1.5\times 10^{10}$
floating point operations per time-step. This number can be seen to increase
rapidly with the polynomial order.

Table 3: Approximate memory requirements of PyFR for the two cylinder meshes. | Device memory / GiB
---|---
Mesh | $\wp=1$ | $\wp=2$ | $\wp=3$ | $\wp=4$
Hex | $0.8$ | $2.0$ | $4.1$ | $7.2$
Pri/tet | $1.1$ | $2.6$ | $4.7$ | $7.7$
Figure 3: Computational effort required for the $118\,070$ element hexahedral
mesh and the mixed mesh with $79\,344$ prismatic elements and $227\,298$
tetrahedral elements.

##  4 Single-Node Performance

###  4.1 Overview

In this section we will analyse performance of PyFR on an Intel Xeon E5-2697
v2 CPU, an NVIDIA Tesla K40c GPU, and an AMD FirePro W9100 GPU. These are the
individual platforms used to construct the multi-node heterogeneous system
employed in section 5.

###  4.2 Hardware specifications

Various attributes of the E5-2697, K40c, and W9100 are detailed in Table 4.
The theoretical peaks for double precision arithmetic and memory bandwidth
were obtained from vendor specifications. However, in practice it is usually
only possible to obtain these theoretical peak values using specially crafted
code sequences. Such sequences are almost always platform specific and seldom
perform useful computations. Consequently, we also calculate and tabulate
_reference_ peaks. Reference peaks for double precision arithmetic are defined
here as the maximum number of giga floating point operations per second
(GFLOP/s) obtained while multiplying two large double precision row-major
matrices using DGEMM from an appropriate BLAS library. Reference peaks for
memory bandwidth are defined here as the rate, in gigabytes per second (GB/s),
that data can be copied between two one gigabyte buffers. Reference peaks for
the E5-2697 were obtained using DGEMM from the Intel Math Kernel Libraries
(MKL) version 11.1.2, and with the E5-2697 paired with four DDR3-1600 DIMMs
(with Turbo Boost enabled). Reference peaks for the K40c were obtained using
DGEMM from cuBLAS as shipped with CUDA 6, with GPU boost disabled and ECC
enabled. Reference peaks for the W9100 were obtained using DGEMM from clBLAS
v2.0 with version 1411.4 of the AMD APP OpenCL runtime. ECC settings for the
W9100 were left unchanged.

We note that on the K40c ECC is implemented in software, and when enabled
error-correction data is stored in global memory. A consequence of this is
that when ECC is enabled there is a reduction in available memory and memory
bandwidth. This partially accounts for the discrepancy observed between the
theoretical and reference memory bandwidths for the K40c. We also note that
for the K40c and the E5-2697, reference peaks for double precision arithmetic
are in excess of 80% of their theoretical values. However, for the W9100 the
reference peak for double precision arithmetic is only 34% of its theoretical
value. This value is not significantly improved via the auto-tuning utility
that ships with clBLAS. It is hoped that this figure will improve with future
releases of clBLAS.

Table 4: Baseline attributes of the three hardware platforms. For the NVIDIA Tesla K40c GPU Boost was left disabled and ECC was enabled. The Intel Xeon E5-2697 v2 was paired with four DDR3-1600 DIMMs with Turbo Boost enabled. | Platform
---|---
| K40c | W9100 | E5-2697
Arithmetic / GFLOP/s |  |  |
theoretical peak | 1430 | 2620 | 280
reference peak | 1192 | 890 | 231
Memory Bandwidth / GB/s |  |  |
theoretical peak | 288 | 320 | 51.2
reference peak | 190 | 261 | 37.1
Thermal Design Power / W | 235 | 275 | 130
Memory / GB | 12 | 16 |
Clock / MHz | 745 | 930 | 3000
Transistors / Billion | 7.1 | 6.2 | 4.3

Finally, as an aside, we note that the number of ‘cores’ available on each
platform have been deliberately omitted from Table 4. It is our contention
that the term is both ill-defined and routinely subject to abuse in the
literature. For example, the E5-2697 is presented by Intel as having 12 cores,
whereas the K40c is described by NVIDIA as having 2880 ‘CUDA cores’. However,
whereas the cores in the E5-2697 can be considered linearly independent those
in the K40c can not. The rough equivalent of a CPU core in NVIDIA parlance is
a ‘streaming multiprocessor’, or SMX, of which the K40c has 15\. Additionally,
the E5-2697 has support for two-way simultaneous multithreading—referred to by
Intel as Hyper-Threading—permitting two threads to execute on each core. At
any one instant it is therefore possible to have up to 24 independent threads
resident on a single E5-2697. The AMD equivalent of a CUDA core is a ‘stream
processor’ of which the W9100 has 2816. This is not to be confused with the
aforementioned streaming multiprocessor of NVIDIA; for which the AMD
equivalent is a ‘Compute Unit’. Practically, both CUDA cores and stream
processors are closer to the individual vector lanes of a traditional CPU
core. Given this minefield of confusing nomenclature we have instead opted to
simply state the peak floating point capabilities of the hardware.

###  4.3 Performance

By measuring the wall clock time required for PyFR to take a defined number of
time-steps, and utilising the operation counts per time-step detailed in
Figure 3, one can calculated the sustained performance of PyFR in GFLOP/s when
running with the meshes detailed in subsection 3.3 and $\wp=1,2,3,4$.

Sustained performance of PyFR for the various hardware platforms is shown in
Figure 4. From the figure it is clear that the computational efficiency of
PyFR increases with the polynomial order. This is consistent with higher order
simulations having an increased compute intensity per degree of freedom. This
additional intensity results in larger operator matrices that are better
suited to the tiling schemes employed by BLAS libraries. The OpenCL
implementation shipped by NVIDIA as part of CUDA only supports the use of
32-bit memory pointers. As such a single context is limited to $4\,\text{GiB}$
of memory, _cf._ Table 3. It was therefore not possible to perform the third
and fourth order simulations for either of the two meshes using the OpenCL
backend with the K40c.

Figure 4: Sustained performance of PyFR in GFLOP/s for the various pieces of
hardware. The backend used by PyFR is given in parentheses. For the OpenCL
backend the initial of the vendor is suffixed. As the NVIDIA OpenCL platform
is limited to $4\,\text{GiB}$ of memory no results are available for
$\wp=3,4$.

The Intel and AMD implementations of OpenCL, when used in conjunction with
clBLAS, are only competitive with the C/OpenMP backend when $\wp=1$ for the
hexahedral mesh, and $\wp=1,2$ for the prism/tetrahedral mesh. This is also
the case when comparing performance between the CUDA backend and the NVIDIA
OpenCL backend on the K40c. Prior analysis by Witherden et al. [8] suggests
that at these orders a reasonable proportion of the wall clock time will be
spent in the bandwidth-bound point-wise kernels as opposed to DGEMM. On
account of being bandwidth-bound such kernels do not extensively test the
optimisation capabilities of the compiler. By the time $\wp=4$ both
implementations of OpenCL on the E5-2697 are delivering between one third and
one quarter of the performance of the native backends. This highlights the
difficulties associated with writing a performant DGEMM routine and hence
severely impacts the performance portability of the OpenCL backend. This
result lends credence to our opening assertion that there are currently no
programming methodologies that are performance portable across a range of
platforms. Further, it also justifies the approach that has been adopted by
PyFR.

Performance of the K40c culminates at $647\,\text{GFLOP/s}$ for the $\wp=4$
hexahedral mesh. This represents some $45\%$ of the theoretical peak and
$54\%$ of the reference peak. By comparison the E5-2697 obtains
$132\,\text{GFLOP/s}$ for the same simulation equating to $47\%$ and $57\%$ of
the theoretical and reference peaks, respectively. Performance does improve
slightly to $140\,\text{GFLOP/s}$ for the $\wp=4$ prism/tetrahedral mesh,
however. On this same mesh at $\wp=4$ the W9100 can be seen to sustain
$657\,\text{GFLOP/s}$ of throughput. Although, in absolute terms, this
observation represents the highest sustained rate of throughput it corresponds
to just $25\%$ of the theoretical peak. However, working in terms of
realisable peaks, we find PyFR obtaining some $74\%$ of the reference value.

##  5 Multi-Node Heterogeneous Performance and Accuracy

###  5.1 Overview

Having determined the performance characteristics of PyFR on various
individual platforms, we will now investigate the ability of PyFR to undertake
simulations on a multi-node heterogeneous system containing an Intel Xeon
E5-2697 v2 CPU, an NVIDIA Tesla K40c GPU, and an AMD FirePro W9100 GPU.

###  5.2 Mesh partitioning

In order to distribute a simulation across the nodes of the heterogeneous
system it is first necessary to partition the mesh. High quality partitions
can be readily obtained using a graph partitioning package such as METIS [24]
or SCOTCH [25].

When partitioning a mixed element mesh for a _homogeneous_ cluster it is
necessary to suitably weight each element type according to its computational
cost. This cost depends both upon the platform on which PyFR is running and
the order at which the simulation is being performed. In principle it is
possible to measure this cost; however in practice the following set of
weights have been found to give satisfactory results across most polynomial
orders and platforms

$\text{hex}:\text{pri}:\text{tet}=3:2:1,$

where larger numbers indicate a greater computational cost. One subtlety that
arises here, is that from a graph partitioning standpoint there is no penalty
associated with placing a sole vertex (element) of a given weight inside of a
partition. Computationally, however, there is a very real penalty incurred
from having just a single element of a certain type inside of the partition.
It is therefore desirable to avoid mesh partitions where any one partition
contains less than around a thousand elements of a given type. An exception is
when a partition contains no elements of such a type—in which case zero
overheads are incurred.

When partitioning a mesh with one type of element for a _heterogeneous_
cluster it is necessary to weight the partition sizes in line with the
performance characteristics of the hardware on each node. However, in the case
of a mixed element mesh on a heterogeneous cluster the weight of an element is
no longer static but rather depends on the partition that it is placed in—a
significantly richer problem. Solving such a problem is currently beyond the
capabilities of most graph partitioning packages. Accordingly, mixed element
meshes that are partitioned for heterogeneous clusters often exhibit inferior
load balancing than those partitioned for homogeneous systems. Moreover, for
consistent performance it is necessary to dedicate a CPU core to each
accelerator in the system. The amount of useful computation that can be
performed by the host CPU is therefore reduced in accordance with this.

Given the single-node performance numbers of Figure 4 it comports to pair the
E5-2697 with the C/OpenMP backend, the K40c with the CUDA backend, and the
W9100 with the OpenCL backend, in order to achieve optimal performance.
Employing the results of Figure 4, in conjunction with some light
experimentation, a set of partitioning weights were obtained and are tabulated
in Table 5.

Table 5: Partition weights for the multi-node heterogeneous simulation. | $\text{E5-2697}:\text{W9100}:\text{K40c}$
---|---
Mesh | $\wp=1$ | $\wp=2$ | $\wp=3$ | $\wp=4$
Hex | $3{:}27{:}23$ | $3{:}27{:}24$ | $4{:}24{:}26$ | $4{:}24{:}28$
Pri/tet | $5{:}33{:}17$ | $5{:}33{:}17$ | $5{:}30{:}20$ | $5{:}27{:}23$

###  5.3 Performance

Sustained performance of PyFR on the multi-node heterogeneous system for each
of the meshes detailed in subsection 3.3 with $\wp=1,2,3,4$ is shown in Figure
5. Under the assumptions of perfect partitioning and scaling one would expect
the sustained performance of the heterogeneous simulation to be equivalent to
the sum of the E5-2697 (C/OpenMP), K40c (CUDA), and W9100 (OpenCL A) bars in
Figure 4. However, for reasons outlined in the preceding paragraphs these
assumptions are unlikely to hold. Some of the available FLOP/s can therefore
be considered as ‘lost’. For the hexahedral mesh the fraction of lost FLOP/s
varies from $22.5\%$ when $\wp=1$ to $8.7\%$ in the case of $\wp=4$. With the
exception of $\wp=1$ the fraction of lost FLOP/s are a few percent higher for
the mixed mesh. This is understandable given the additional complexities
associated with mixed mesh partitioning and can likely be improved upon by
switching to order-dependent element weighting factors.

Figure 5: Sustained performance of PyFR on the multi-node heterogeneous system
for each mesh with $\wp=1,2,3,4$. Lost FLOP/s represent the difference between
the achieved FLOP/s and the sum of the E5-2697 (C/OpenMP), K40c (CUDA), and
W9100 (OpenCL A) bars in Figure 4.

###  5.4 Accuracy

In this section we present instantaneous and time-span-averaged (henceforth
referred to as averaged) results obtained using the multi-node heterogeneous
system and the mixed unstructured prism/tetrahedral mesh with $\wp=1,2,3,4$.
Following the approach of Breuer [19] all averaged results were obtained over
$100D/u_{\infty}$ in time, and the full span of the domain. The $\wp=1$
simulation was initialised with spatially constant free-stream initial
conditions, and run for a lead in time of ${\sim}50D/u_{\infty}$ before time-
averaging began. Subsequent simulations, at $\wp>1$, were initialised with the
final solution state from the previous simulation at $\wp-1$, and time-
averaging began immediately. For all $\wp$ this approach led to averaged
results exhibiting U-shape or Mode-L characteristics (following the
terminology of Ma et al. [16] and Lehmkuhl et al. [21] respectively). Hence,
all averaged results were compared with the experimental results of Norberg et
al. [18] and Parnaudeau et al. [17], which also exhibited U-shape/Mode-L
characteristics, as well as Mode-L results from the recent DNS study of
Lehmkuhl et al. [21].

Instantaneous surfaces of iso-density are shown in Figure 6 for each
polynomial order. We observe that the $\wp=1$ simulation captures
predominantly large-scale structures in the turbulent wake behind the
cylinder. As $\wp$ is increased we are able to capture a greater number of
small scale turbulent structures. This is due to an increase in the total
number of degrees of freedom in the domain, as well as a relative decrease in
discretisation errors for the higher-order schemes. We observe laminar flow at
the leading edge of the cylinder for all test cases, transition near the
separation points, and finally turbulent flow in the wake region. These are
the characteristic features of the shear-layer transition regime as described
by Williamson [15].

(a) $\wp=1$.

(b) $\wp=2$.

(c) $\wp=3$.

(d) $\wp=4$.

Figure 6: Instantaneous surfaces of iso-density coloured by velocity
magnitude.

Plots of averaged pressure coefficient $\overline{C_{p}}$ on the surface of
the cylinder are shown in Figure 7. The results are shown alongside the
experimental results of Norberg et al. [18] (extracted from Kravchenko and
Moin [20]), and the Mode-L DNS results of Lehmkuhl et al. [21]. At $\wp=1$ we
observe a large negative pressure coefficient near the top and bottom of the
cylinder, which includes the location of maximum skin friction coefficient
[19]. However, with increasing $\wp$ the results tend towards those of the
previous studies.

The averaged pressure coefficient at the base of the cylinder
$\overline{C_{pb}}$, and the averaged separation angle $\theta_{s}$ measured
from the leading stagnation point are tabulated in Table 6, along with
experimental data from Norberg et al. [18] at a similar $Re=4020$ (extracted
from Kravchenko and Moin [20]), experimental data from Parnaudeau et al. [17],
and Mode-L DNS data from Lehmkuhl et al. [21]. Once again, with increasing
$\wp$ the results tend towards those of the previous studies.

Figure 7: Averaged pressure coefficient for $\wp=1,2,3,4$ compared with the experimental results of Norberg et al. [18] (from Kravchenko and Moin [20]) and DNS results of Lehmkuhl et al. [21]. Table 6: Comparison of quantitative values with experimental and DNS results. | $-\overline{C_{pb}}$ | $\theta_{s}/{}^{\circ}$
---|---|---
PyFR pri/tet $\wp=1$ | $0.752$ | $92$
PyFR pri/tet $\wp=2$ | $0.796$ | $89$
PyFR pri/tet $\wp=3$ | $0.808$ | $88$
PyFR pri/tet $\wp=4$ | $0.881$ | $88$
Parnaudeau et al. [17] |  | $88$
Lehmkuhl et al. [21] | $0.877$ | $87.8$
Norberg et al. [18, 20] | $0.880$ |

Plots of averaged stream-wise velocity at $x/D=1.06$, $1.54$, and $2.02$ are
shown in Figure 8. The results are shown alongside the experimental results of
Parnaudeau et al. [17] and the Mode-L DNS results of Lehmkuhl et al. [21].
With $\wp=1$ the profiles deviate significantly from the previous studies. On
increasing the order to $\wp=2$ the results are improved. We observe a U-shape
profile at $x/D=1.06$, with strong gradients near the mixing layer between the
wake and the free stream. The $\wp=4$ results agree well those of the previous
studies. The only exception is the reduced magnitude of the averaged velocity
deficit near the centre of the wake at $x/D=1.54$.

Plots of averaged cross-wise velocity at $x/D=1.06$, $1.54$, and $2.02$ are
shown in Figure 9. The results are also shown alongside the experimental
results of Parnaudeau et al. [17] and the Mode-L DNS results of Lehmkuhl et
al. [21]. The $\wp=4$ results agree well with those of the previous studies.

Figure 8: Time-span-average stream-wise velocity for $\wp=1,2,3,4$ compared
with the experimental results of Parnaudeau et al. [17] and DNS results of
Lehmkuhl et al. [21]. Figure 9: Time-span-average cross-stream velocity for
$\wp=1,2,3,4$ compared with the experimental results of Parnaudeau et al. [17]
and DNS results of Lehmkuhl et al. [21].

##  6 Conclusions

In this paper we have demonstrated the ability of PyFR to perform high-order
accurate unsteady simulations of flow on mixed unstructured grids using
_heterogeneous_ multi-node hardware. Specifically, after benchmarking single-
node performance for various platforms, PyFR v0.2.2 was used to undertake
simulations of unsteady flow over a circular cylinder at Reynolds number
$3\,900$ using a mixed unstructured grid of prismatic and tetrahedral elements
on a desktop workstation containing an Intel Xeon E5-2697 v2 CPU, an NVIDIA
Tesla K40c GPU, and an AMD FirePro W9100 GPU. Results demonstrate that PyFR
achieves performance portability across various hardware platforms. In
particular, the ability of PyFR to target individual platforms with their
‘native’ language leads to significantly enhanced performance _cf._ targeting
each platform with OpenCL alone. PyFR was also found to be performant on the
heterogeneous multi-node system achieving a significant fraction of the
available FLOP/s. Further, the numerical results obtained using a mixed
unstructured grid of prismatic and tetrahedral elements were found to be in
good agreement with previous experimental and numerical data.

## Acknowledgements

The authors would like to thank the Engineering and Physical Sciences Research
Council for their support via a Doctoral Training Grant and an Early Career
Fellowship (EP/K027379/1). The authors would also like to thank AMD, Intel,
and NVIDIA for hardware donations.

## References

  * [1] WH Reed and TR Hill.  Triangular mesh methods for the neutron transport equation.  Technical Report LA-UR-73-479, Los Alamos Scientific Laboratory, 1973.
  * [2] David A Kopriva and John H Kolias.  A conservative staggered-grid Chebyshev multidomain method for compressible flows.  Journal of computational physics, 125(1):244–261, 1996.
  * [3] Yuzhi Sun, Zhi Jian Wang, and Yen Liu.  High-order multidomain spectral difference method for the Navier-Stokes equations on unstructured hexahedral grids.  Communications in Computational Physics, 2(2):310–333, 2007.
  * [4] HT Huynh.  A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods.  AIAA paper, 4079:2007, 2007.
  * [5] Jan S Hesthaven and Tim Warburton.  Nodal discontinuous Galerkin methods: algorithms, analysis, and applications, volume 54.  Springer Verlag New York, 2008.
  * [6] Andreas Klöckner, Tim Warburton, Jeff Bridge, and Jan S Hesthaven.  Nodal discontinuous Galerkin methods on graphics processors.  Journal of Computational Physics, 228(21):7863–7882, 2009.
  * [7] Patrice Castonguay, David M Williams, Peter E Vincent, Manuel Lopez, and Antony Jameson.  On the development of a high-order, multi-gpu enabled, compressible viscous flow solver for mixed unstructured grids.  AIAA paper, 3229:2011, 2011.
  * [8] F.D. Witherden, A.M. Farrington, and P.E. Vincent.  PyFR: An open source framework for solving advection–diffusion type problems on streaming architectures using the flux reconstruction approach.  Computer Physics Communications, 185(11):3028–3040, 2014.
  * [9] Kazushige Goto and Robert A Geijn.  Anatomy of high-performance matrix multiplication.  ACM Transactions on Mathematical Software (TOMS), 34(3):12, 2008\.
  * [10] Michael Bayer.  Mako: Templates for python, 2013.
  * [11] Andreas Klöckner, Nicolas Pinto, Yunsup Lee, Bryan Catanzaro, Paul Ivanov, and Ahmed Fasih.  Pycuda and pyopencl: A scripting-based approach to gpu run-time code generation.  Parallel Comput., 38(3):157–174, 2012.
  * [12] Anatol Roshko.  On the development of turbulent wakes from vortex streets.  Technical Report No. NACA TR 1191, California Institute of Technology, 1953.
  * [13] MS Bloor.  The transition to turbulence in the wake of a circular cylinder.  Journal of Fluid Mechanics, 19:290–304, 1964.
  * [14] CHK Williamson.  The existence of two stages in the transition to three dimensionality of a cylinder wake.  Physics of Fluids, 31:3165–3168, 1988.
  * [15] CHK Williamson.  Vortex dynamics in the cylinder wake.  Annual Review of Fluid Mechanics, 28:477–539, 1996.
  * [16] X Ma, GS Karamanos, and GE Karniadakis.  Dynamics and low-dimensionality of a turbulent near wake.  Journal of Fluid Mechanics, 310:29–65, 2000.
  * [17] P Parnaudeau, J Carlier, D Heitz, and E Lamballais.  Experimental and numerical studies of the flow over a circular cylinder at reynolds number 3900.  Journal of Fluid Mechanics, 310:29–65, 2000.
  * [18] C Norberg.  Ldv measurements in the near wake of a circular cylinder.  International Journal for Numerical Methods in Fluids, 28(9):1281–1302, 1998.
  * [19] M Breuer.  Large eddy simulation of the subcritical flow past a circular cylinder.  International Journal for Numerical Methods in Fluids, 28(9):1281–1302, 1998.
  * [20] AG Kravchenko and P Moin.  Numerical studies of flow over a circular cylinder at $re_{d}=3900$.  Physics of Fluids, 12:403–417, 2000.
  * [21] O Lehmkuhl, I Rodriguez, R Borrell, and A Oliva.  Low-frequency unsteadiness in the vortex formation region of a circular cylinder.  Physics of Fluids, 25(8):3165–3168, 2013.
  * [22] Philippe Parnaudeau, Johan Carlier, Dominique Heitz, and Eric Lamballais.  Experimental and numerical studies of the flow over a circular cylinder at reynolds number 3900.  Physics of Fluids, 20(8), 2008.
  * [23] Christopher A Kennedy, Mark H Carpenter, and R Michael Lewis.  Low-storage, explicit runge–kutta schemes for the compressible navier–stokes equations.  Applied numerical mathematics, 35(3):177–219, 2000.
  * [24] George Karypis and Vipin Kumar.  A fast and high quality multilevel scheme for partitioning irregular graphs.  SIAM Journal on Scientific Computing, 20(1):359–392, 1998.
  * [25] François Pellegrini and Jean Roman.  Scotch: A software package for static mapping by dual recursive bipartitioning of process and architecture graphs.  In High-Performance Computing and Networking, pages 493–498. Springer, 1996.


# Hybrid Reinforcement Learning Framework for Mixed-Variable Problems

Haoyan Zhai
<EMAIL_ADDRESS>
School of Mathematics
Georgia Institute of Technology &Qianli Hu
<EMAIL_ADDRESS>
School of Mathematics
Georgia Institute of Technology &Jiangning Chen
<EMAIL_ADDRESS>
UKG AI
Ultimate Kronos Group

###### Abstract

Optimization problems characterized by both discrete and continuous variables
are common across various disciplines, presenting unique challenges due to
their complex solution landscapes and the difficulty of navigating mixed-
variable spaces effectively. To Address these challenges, we introduce a
hybrid Reinforcement Learning (RL) framework that synergizes RL for discrete
variable selection with Bayesian Optimization for continuous variable
adjustment. This framework stands out by its strategic integration of RL and
continuous optimization techniques, enabling it to dynamically adapt to the
problem’s mixed-variable nature. By employing RL for exploring discrete
decision spaces and Bayesian Optimization to refine continuous parameters, our
approach not only demonstrates flexibility but also enhances optimization
performance. Our experiments on synthetic functions and real-world machine
learning hyperparameter tuning tasks reveal that our method consistently
outperforms traditional RL, random search, and standalone Bayesian
optimization in terms of effectiveness and efficiency.

##  1 Introduction

Optimization is a cornerstone of modern technological advancements,
underpinning critical applications in machine learning (Sun et al., 2019),
operation research (Rardin & Rardin, 1998), finance (Cornuejols & Tutuncu,
2005), and engineering (Deb, 2012). As the problem complexity in these domains
escalates - along with the exponential growth in the data sizes - the
limitations of traditional optimization methods become increasingly apparent.
This is particularly true for challenges that require tuning discrete and
continuous variables simultaneously. While a plethora of optimization
techniques have been developed to address either discrete or continuous
variables effectively, seamlessly optimizing across a space imbued with both
remains a formidable challenge.

Historically, researchers have treated continuous and discrete optimization
tasks separately. On the one hand, Reinforcement Learning (RL) has carved its
niche as a powerful tool for discrete optimization problems, lauded for its
ability to learn and adapt strategies through direct interaction with the
environment, see Neunert et al. (2020); Mazyavkina et al. (2021); Poupart et
al. (2006). Parallelly, continuous optimization techniques, notably Bayesian
Optimization (BO), Gradient Descent, and Evolutionary Algorithms, have
established their efficacy in navigating the subtleties of continuous
parameter space (Andréasson et al., 2020; Zhan et al., 2022; Omidvar et al.,
2015; Lauer et al., 2011; Kistler & Franz, 2003).

Despite these advancements, the quest for a unified approach that adeptly
handles the intricacies of mixed-variable optimization tasks remains. These
standalone applications of either discrete or continuous optimization methods
tend to falter when faced with the multifaceted nature of real-world problems,
underscoring a distinct gap in the current optimization toolbox (Liao et al.,
2013; Dimopoulos, 2007; Lin et al., 2018).

Figure 1: Overview of our hybrid Reinforcement Learning system. Each RL
iteration is augmented by an n-step Bayesian Optimization for continuous
variables until the stop criterion is met.

To bridge this divide, we introduce a flexible hybrid Reinforcement Learning
framework (depicted in Figure 1) that melds the strategic acumen of RL with
the precision of continuous optimization techniques. This innovative approach
not only navigates the complex terrain of mixed-variable problems with
dexterity but also sheds some light on a potential paradigm in optimization,
where the synergistic potential of integrating discrete and continuous
methodologies is fully realized (Wang & Zheng, 2001; Piccoli, 1999; Cavazos et
al., 2006; Kelner et al., 2008).

Our contributions through this work are manifold. Firstly, we introduce a
novel Hybrid Framework that aims to integrate Reinforcement Learning (RL) with
a diverse set of continuous optimization techniques, which provides a
comprehensive strategy for effectively addressing mixed-variable optimization
challenges. Secondly, our contributions are substantiated through empirical
validation, involving rigorous experimentation on both synthetic functions and
practical scenarios. The results consistently demonstrate the superior
performance of our Hybrid Framework when compared to traditional RL methods,
random search approaches, and standalone continuous optimization methods.
Lastly, we highlight the adaptability and generalization of our framework
across various continuous optimization strategies, showcasing its versatility
and applicability across a wide range of optimization landscapes.

This paper is structured to guide the reader through the evolution of our
framework and its validation. Section 2 provides a review of related work,
setting the stage for our methodology. Section 3 delves into the flexible
hybrid Reinforcement Learning framework, detailing its components and
operational dynamics. Section 4 outlines our experimental design and the
ensuing results, highlighting the effectiveness and versatility of the
approach. Section 5 concludes the paper, summarizing our contributions and
envisioning the path forward.

##  2 Related Work

The exploration of optimization strategies within artificial intelligence has
led to significant advancements across various domains. Reinforcement Learning
(RL) and continuous optimization techniques stand out as pivotal components in
this evolutionary trajectory. This section reviews the pertinent literature,
focusing on the development and application of RL in optimization, the role of
continuous optimization methods, and the inception of hybrid frameworks that
seek to leverage the strengths of both approaches.

Reinforcement Learning has emerged from its roots in dynamic programming to
become a key methodology for solving decision-making problems characterized by
uncertainty and complexity. The foundational work by Sutton & Barto (2018)
provides a comprehensive overview of RL principles, outlining how agents learn
to make decisions through interactions with their environment to maximize some
notion of cumulative reward.

In the context of optimization, RL has been particularly impactful in discrete
spaces, where the goal is to find optimal sequences or configurations from a
finite set of possibilities. Applications range from combinatorial
optimization (Mazyavkina et al., 2021) to more nuanced domains like network
configuration (Gao et al., 2019; Mammeri, 2019) and inventory management
(Giannoccaro & Pontrandolfo, 2002; Sultana et al., 2020). These studies
underscore RL’s versatility and capacity to adapt strategies based on
feedback, a trait that is invaluable for navigating the complex landscapes of
optimization problems.

Parallel to the advancements in RL, continuous optimization techniques have
been refined and adapted to meet the challenges of tuning continuous
parameters. Techniques such as Bayesian Optimization (BO), Gradient Descent,
and Evolutionary Algorithms have been at the forefront of this endeavor.
Bayesian Optimization, in particular, has gained prominence for its efficiency
in hyperparameter tuning of machine learning models (Wu et al., 2019),
leveraging probabilistic models to guide the search process. Gradient Descent
remains a staple in training neural networks by iteratively minimizing loss
functions (Bengio, 2000), while Evolutionary Algorithms offer a population-
based search strategy, drawing inspiration from biological evolution to
optimize complex functions (Young et al., 2015; Tani et al., 2021).

Recognizing the complementary strengths of RL and continuous optimization
techniques, recent research has begun to explore hybrid frameworks that
integrate these methodologies. Such frameworks aim to capitalize on RL’s
ability to make strategic decisions in discrete spaces, while employing
continuous optimization methods for fine-tuning continuous variables. This
synergy has the potential to offer more robust and efficient solutions to
mixed-variable optimization problems, a frontier that is increasingly relevant
in areas such as automated machine learning (Liu et al., 2020) and dynamic
resource allocation (Bagci & Tekalp, 2018).

Initial attempts at such integration have shown promising results, suggesting
that the combined approach could outperform the application of either
methodology in isolation (Lin et al., 2018; Wang et al., 2017; 2021). These
hybrid models not only offer a nuanced understanding of the optimization
landscape but also open new avenues for research in both theoretical
development and practical applications.

Despite these advancements, there remains a notable gap in the literature
concerning efficient, scalable, and adaptable frameworks that can effectively
navigate the complexities of mixed-variable optimization landscapes. The
development of a flexible hybrid Reinforcement Learning framework, as proposed
in this paper, aims to address this gap by offering a novel integration of RL
and continuous optimization techniques, tailored to the demands of
contemporary optimization challenges.

##  3 Methodology

###  3.1 Problem Definition

This study focuses on optimization problems that feature a combination of
discrete and continuous decision variables, with the aim of maximizing an
objective function. Let $\mathcal{A}$ represent the space of discrete
variables and $\mathcal{X}$ the space of continuous variables. The objective
is to find the optimal pair $(a^{*},x^{*})$ that maximizes the objective
function $f:\mathcal{A}\times\mathcal{X}\rightarrow\mathbb{R}$, which can be
formally expressed as:

$\max_{a\in\mathcal{A},x\in\mathcal{X}}f(a,x)$

###  3.2 Hybrid Reinforcement Learning Framework

The flexible hybrid Reinforcement Learning framework addresses mixed-variable
optimization problems by bifurcating the optimization process: it employs
Reinforcement Learning (RL) for discrete variable selection and a specific
continuous optimization technique for continuous variable adjustment. Herein,
we employ the Gradient Bandit model with Softmax action selection for RL, and
Bayesian Optimization for updating the reward function.

####  3.2.1 Reinforcement Learning for Discrete Variable Optimization

For the discrete variable optimization component, we utilize the Gradient
Bandit agent, which is designed to select actions that optimize cumulative
rewards over time. The action selection is governed by the Softmax policy,
defined as follows:

$Pr\\{A_{t}=a\\}\frac{\exp(H_{t}(a))}{\sum_{b=1}^{n}\exp(H_{t}(b))}=\pi_{t}(a),$
(1)

where $A_{t}$ represents the action taken at time $t$, $\pi_{t}(a)$ the
probability of taking action $a$ at time $t$, and $H_{t}(a)$ the learned
numerical preference for each action $a$. The preferences are updated upon
selecting action $A_{t}$ and observing reward $R_{t}$, using the formulae:

$\displaystyle
H_{t+1}(A_{t})=H_{t}(A_{t})+\alpha(R_{t}-\bar{R_{t}})(1-\pi_{t}(A_{t})),$
$\displaystyle and$ (2) $\displaystyle
H_{t+1}(a)=H_{t}(a)-\alpha(R_{t}-\bar{R_{t}})\pi_{t}(a),$
$\displaystyle\forall a\neq A_{t}$ (3)

where $\alpha>0$ is a step-size parameter, and $\bar{R}_{t}\in\mathbb{R}$
denotes the average of all rewards up to and including time $t$ (Sutton &
Barto, 1999). This mechanism ensures an exploration-exploitation balance in
action selection, enabling the agent to prioritize actions based on their
reward potential. Through interactions with the objective function, the RL
agent dynamically refines its policy based on reward feedback, enhancing its
decision-making accuracy over time.

####  3.2.2 Continuous Optimization Techniques

A crucial aspect of the Gradient Bandit’s effectiveness is defining the reward
function $R_{t}$. For this, we employ Bayesian Optimization. After setting the
discrete variables $A_{t}$, we conduct several iterations of Bayesian
Optimization on the problem: $\max_{x\in\mathcal{X}}f(A_{t},x)$, retaining the
solver status for $A_{t}$. If $A_{t}$ has been previously encountered, we
reload its continuous problem status and continue optimization; otherwise, we
initiate a new problem. The outcome of Bayesian Optimization, denoted as
$\\{f(A_{t},X_{1}),f(A_{t},X_{2}),\ldots,f(A_{t},X_{s})\\}$, informs the
definition of the reward:

$R_{t}(A_{t})=\max_{s}f(A_{t},X_{s}).$ (4)

The whole algorithm is given as Alg.1. While this methodology employs the
Gradient Bandit and Bayesian Optimization for demonstration, it is important
to note that the framework is flexible and can incorporate other RL methods
and continuous optimization algorithms to suit different problem contexts and
requirements.

Input:

$f(a,x)$,

Bayesian Optimization Iteration Steps $n$,

$H_{0}(a)=0$,

$Cache[a]=\\{\\}$, $\forall a\in\mathcal{A}$,

1 for  _$t=0,1,2,3\cdots$_ do

2 Calculate $\pi_{t}(a)$ $\forall a\in\mathcal{A}$.

3 Choose an action $A_{t}\in\mathcal{A}$ according to the probability
distribution $\pi_{t}$.

4 if  _$Cache[A_{t}]$ is empty_ then

5 Initialize a Bayesian Optimizer $opt$.

6 else

7 $opt=Cache[A_{t}]$.

8 end if

9 Iterate $opt$ for $n$ steps, and get the reward $R_{t}=\max_{s\leq
t}f(A_{t},X_{sn})$, where $X_{sn}$ is the iteration results in $opt$.

10 Save the current Bayesian Optimization status $Cache[A_{t}]=opt$.

11 Update $H_{t}$ according to 3.

12 if  _$(X_{t},R_{t}(X_{t}))$ appears $m$ times in the last $T$ iterations_
then

13 Break

14 end if

15

16 end for

Algorithm 1 Hybrid RL Optimizer

##  4 Experiments

###  4.1 Experimental Setup

The experiments were conducted on a computing environment equipped with AMD
Ryzen Threadripper 3970X 32-Core Processor with 256GB RAM. The framework was
implemented in Python, with reinforcement learning components based on
Galbraith (2023), and continuous optimization techniques, including Bayesian
Optimization, Gradient Descent, and Evolutionary Algorithms, were utilized as
per the framework’s design.

###  4.2 Datasets and Scenarios

To assess the framework’s performance, we use a combination of synthetic
functions and real-world scenarios. In the experiments, we use Random Search,
Gradient Bandit, Bayesian Optimization, and our Hybrid method to perform the
optimization problems. Since Gradient Bandit only takes discrete variables, we
discretize continuous variables in their corresponding domains. For Bayesian
Optimization, we choose variable values in the range of the minimum and the
maximum of the variable, and round it to the nearest discrete value in the
domains.

####  4.2.1 Synthetic Functions

Three Benchmark functions are employed to evaluate the framework’s capability
to find global optima in complex landscapes. The objective function is simply
the function values, and we report the absolute difference between the
calculated function values and the known global maximum in Section 4.3

  1. 1.

Shekel Function: The first one is the Shekel function, which is a benchmark
test function in global optimization problems. The function is given as below:

$f_{1}(x_{1},x_{2},x_{3},x_{4})=\sum_{i=1}^{4}\left(c_{i}+\sum_{j=1}^{4}(x_{j}-a_{ji})^{2}\right)^{-1},$
(5)

where

$C=(c_{1},\cdots,c_{10})=\frac{1}{10}(1,2,2,4,4,6,3,7,5,5),$

and

$A=(a_{ji})=\left(\begin{array}[]{cccccccccc}4&1&8&6&3&2&5&8&6&7\\\
4&1&8&6&7&9&3&1&2&3.6\\\ 4&1&8&6&3&2&5&8&6&7\\\ 4&1&8&6&7&9&3&1&2&3.6\\\
\end{array}\right).$

To make it a function containing both discrete and continuous variables, we
set the first two variables to be integers. This function has 10 local
maximums and the global maximum is around $10.536283726219603$ at $(4,4,4,4)$,
and the searching region is
$\\{0,1,2,3,4,5,6,7,8,9,10\\}\times\\{0,1,2,3,4,5,6,7,8,9,10\\}\times[0,10]\times[0,10]$.

  2. 2.

Composition Function: The second function is a composition of Rastrigin
(global minimum $f(0,0)=0$), Ackley (global minimum $f(0,0)=0$), and Sphere
(global minimum $f(0,0)=0$) functions:

$\displaystyle\text{Rastrigin}(\mathbf{x})=10n+\sum_{i=1}^{n}\left[x_{i}^{2}-10\cos(2\pi
x_{i})\right],\ \ \ \text{Sphere}(x,y)=x^{2}+y^{2}$ (6)
$\displaystyle\text{Ackley}(x,y)=-20\exp\left[-0.2\sqrt{0.5(x^{2}+y^{2})}\,\right]-\exp\left[0.5\left(\cos
2\pi x+\cos 2\pi y\right)\right]+e+20$ (7)

We define our function as:

$f_{2}(u,x,y)=\left\\{\begin{array}[]{ll}-\text{Rastrigin}(x,y)&u=0\\\
-\text{Ackley}(x,y)+10&u=1\\\ -\text{Sphere}(x,y)+20&u=2\\\
\end{array}\right.$ (9)

Therefore, the global maximum is $f(2,0,0)=20$. We also define the variable
$x$ to be discrete and the search area to be
$\\{0,1,2\\}\times\\{-1,0,1,2,3\\}\times[-5,5]$.

  3. 3.

Sine Permutation Function: At last, we define a two-variable function as below

$g(x,y)=\frac{x\sin\left((-x+7)^{2}\pi/(2(y-4)^{2}+1)\right)}{((x-5)^{2}+1)}.$

Then, we define $u,v,w\in\\{1,4,7,10,13\\}$, $x\in[0.5,8]$, and $y\in[0.1,5]$.
Define the permutations $Perm_{u}(x)=(7,1,13,10,4)$, in which the value of
this function is the integer next to the given integer $x$ in the array, i.e.
$Perm_{u}(1)=13$. Following the same logic, also define
$Perm_{v}(x)=(13,1,4,7,10)$, and $Perm_{w}(x)=(7,4,10,1,13)$. Now we can
define our objective function as

$f_{3}(u,v,w,x,y)=g(Perm_{u}(u)+Perm_{v}(v)+Perm_{w}(w)+x,y).$ (10)

####  4.2.2 Real-World Scenario

Referred to the Walmart dataset in Kaggle (H, 2022), we predict the upcoming
weekly sales with the tree model XGBoost. Seeking suitable hyperparameters to
minimize the mean absolute error (MAE) of the regression problem, we tune
$max\\_bin$, $n\\_estimators$, and $max\\_depth$, which are discrete
parameters, incorporating continuous parameters $colsample\\_bytree$,
$learning\\_rate$, $reg\\_alpha$, and $reg\\_lambda$. To prevent over-fitting
and under-fitting, we split the whole dataset as a training set and a
validation set. Fixing a set of hyperparameters, we train the model on the
training dataset and then calculate the MAE of the validation set. As our
framework is maximizing the objective and to avoid number overflow, we set the
objective function of this problem as $-\text{MAE}/1e5$.

###  4.3 Results

Figure 2: The results for synthetic functions. Each graph is a trajectory of
the optimization. The x-axis is the number of iterations, while the y-axis is
the absolute gap between the searched maximum and the known global maximum.
The red dots are the result of one step of the experiment, and the blue curve
is the rolling average of the adjacent 50 points. Specifically in hybrid RL,
we need to define the Bayesian Opt steps to perform in each reinforcement
learning step $n$: for Shekel function $n=3$, for Composition Function $n=3$,
for Sine Permutation Function $n=2$. Figure 3: The results for machine
learning hyperparameter tuning problem. Each curve is the rolling average of
the adjacent 50 points. The x-axis is the number of iterations, while the
y-axis is the reward. The reward is defined as the objective function value at
each step for Random Search, Bayesian Opt, and RL, while it is defined as
Equation 4. In Hybrid RL, the Bayesian Opt step $n=2$. Figure 4: The statistic
result for all the problems repeated 5 times for synthetic functions and 10
times for the machine learning problem with different random seeds. We can see
that the hybrid RL gives the lowest std in all problems and achieves the
optimal values in most cases.

We have conducted four sets of optimization experiments detailed in section
4.2.1, covering the Shekel function, a composition function, the Sine
Permutation function, and a machine-learning-based hyper-parameter tuning
task). Each experiment is evaluated by the same four optimization methods:
random search, Bayesian Optimization, vanilla Reinforcement Learning, and our
hybrid method. The results demonstrate that our flexible hybrid Reinforcement
Learning framework significantly outperforms the comparative methods. More
specifically:

Convergence Assurance: As shown in Figure 2, our hybrid RL method demonstrates
a high level of convergence certainty. The experiments support the assertion
that the hybrid method has amalgamated the benefits of reinforcement learning
and conventional optimization methods tailored for continuous variables such
as Bayesian Optimization.

Convergence Efficiency: Our hybrid method, enhanced with periodic
Reinforcement Learning checks, outpaces Bayesian Optimization in computational
speed despite requiring more iterations, as Bayesian Optimization’s per-
iteration time increases with expanded search spaces, evident in Figure 3.
This results in our hybrid approach achieving faster convergence within
practical time frames, even beyond 600 rounds where Bayesian Optimization
becomes untenably slow. Compared to vanilla Reinforcement Learning, which
necessitates segmenting continuous spaces into discrete intervals—thus
expanding the action space and computational load—our hybrid method maintains
the action space’s continuity. This unique balance allows it to converge in
fewer iterations than pure RL, despite the individual iterations taking
slightly longer, showcasing an optimal trade-off between iteration count and
time efficiency for complex optimization tasks.

Low Standard Deviation Between Rounds of Experiments: Compared to conventional
methods, our hybrid approach demonstrates remarkable consistency, reflected in
its consistently low standard deviations across the four experiments. Figure 4
illustrates that although our hybrid model lags slightly behind in the ML
hyper-parameter tuning task, the statistical analysis of standard deviations
suggests a notably stable performance. Furthermore, in the remaining three
tests, the hybrid model consistently identifies the global optimum.

##  5 Conclusion and Future Research

This paper presented a flexible hybrid Reinforcement Learning Framework that
innovatively combines reinforcement learning with continuous optimization
techniques such as Bayesian Optimization to tackle mixed-variable optimization
problems. Through our experimentation, we have demonstrated its efficacy and
efficiency over traditional approaches. The framework’s flexibility and
adaptability suggest its potential applicability across a broad spectrum of
optimization challenges. We believe that this work not only advances the
fusion of continuous optimization techniques and reinforcement learning but
also opens new doors for solving complex problems in various domains. Future
research will aim to expand the framework’s capabilities and its application
in more diverse and large-scale scenarios, to theoretically substantiate the
efficacy of our algorithm, and to explore more sophisticated reinforcement
learning and continuous optimization algorithms to further improve the
framework’s effectiveness.

## References

  * Andréasson et al. (2020) Niclas Andréasson, Anton Evgrafov, and Michael Patriksson.  _An introduction to continuous optimization: foundations and fundamental algorithms_.  Courier Dover Publications, 2020.
  * Bagci & Tekalp (2018) K Tolga Bagci and A Murat Tekalp.  Dynamic resource allocation by batch optimization for value-added video services over sdn.  _IEEE Transactions on Multimedia_ , 20(11):3084–3096, 2018.
  * Bengio (2000) Yoshua Bengio.  Gradient-based optimization of hyperparameters.  _Neural computation_ , 12(8):1889–1900, 2000.
  * Cavazos et al. (2006) John Cavazos, J Eliot B Moss, and Michael FP O’Boyle.  Hybrid optimizations: Which optimization algorithm to use?  In _Compiler Construction: 15th International Conference, CC 2006, Held as Part of the Joint European Conferences on Theory and Practice of Software, ETAPS 2006, Vienna, Austria, March 30-31, 2006. Proceedings 15_ , pp. 124–138. Springer, 2006.
  * Cornuejols & Tutuncu (2005) Gerard Cornuejols and Reha Tutuncu.  Optimization methods in finance.  _Camegie Mellon University, Pittsburgh_ , 2005.
  * Deb (2012) Kalyanmoy Deb.  _Optimization for engineering design: Algorithms and examples_.  PHI Learning Pvt. Ltd., 2012.
  * Dimopoulos (2007) George G Dimopoulos.  Mixed-variable engineering optimization based on evolutionary and social metaphors.  _Computer methods in applied mechanics and engineering_ , 196(4-6):803–817, 2007.
  * Galbraith (2023) Byron Galbraith.  Bandits: A collection of bandit algorithms.  https://github.com/bgalbraith/bandits.git, 2023.
  * Gao et al. (2019) Yuanqi Gao, Jie Shi, Wei Wang, and Nanpeng Yu.  Dynamic distribution network reconfiguration using reinforcement learning.  In _2019 IEEE International Conference on Communications, Control, and Computing Technologies for Smart Grids (SmartGridComm)_ , pp. 1–7. IEEE, 2019.
  * Giannoccaro & Pontrandolfo (2002) Ilaria Giannoccaro and Pierpaolo Pontrandolfo.  Inventory management in supply chains: a reinforcement learning approach.  _International Journal of Production Economics_ , 78(2):153–161, 2002.
  * H (2022) Yasser H.  Walmart dataset, 2022.  URL https://www.kaggle.com/datasets/yasserh/walmart-dataset.
  * Kelner et al. (2008) Vincent Kelner, Florin Capitanescu, Olivier Léonard, and Louis Wehenkel.  A hybrid optimization technique coupling an evolutionary and a local search algorithm.  _Journal of Computational and Applied Mathematics_ , 215(2):448–456, 2008.
  * Kistler & Franz (2003) Thomas Kistler and Michael Franz.  Continuous program optimization: A case study.  _ACM Transactions on Programming Languages and Systems (TOPLAS)_ , 25(4):500–548, 2003.
  * Lauer et al. (2011) Fabien Lauer, Gérard Bloch, and René Vidal.  A continuous optimization framework for hybrid system identification.  _Automatica_ , 47(3):608–613, 2011.
  * Liao et al. (2013) Tianjun Liao, Krzysztof Socha, Marco A Montes de Oca, Thomas Stützle, and Marco Dorigo.  Ant colony optimization for mixed-variable optimization problems.  _IEEE Transactions on evolutionary computation_ , 18(4):503–518, 2013.
  * Lin et al. (2018) Ying Lin, Yu Liu, Wei-Neng Chen, and Jun Zhang.  A hybrid differential evolution algorithm for mixed-variable optimization problems.  _Information Sciences_ , 466:170–188, 2018.
  * Liu et al. (2020) Sijia Liu, Parikshit Ram, Deepak Vijaykeerthy, Djallel Bouneffouf, Gregory Bramble, Horst Samulowitz, Dakuo Wang, Andrew Conn, and Alexander Gray.  An admm based framework for automl pipeline configuration.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pp. 4892–4899, 2020.
  * Mammeri (2019) Zoubir Mammeri.  Reinforcement learning based routing in networks: Review and classification of approaches.  _Ieee Access_ , 7:55916–55950, 2019.
  * Mazyavkina et al. (2021) Nina Mazyavkina, Sergey Sviridov, Sergei Ivanov, and Evgeny Burnaev.  Reinforcement learning for combinatorial optimization: A survey.  _Computers & Operations Research_, 134:105400, 2021.
  * Neunert et al. (2020) Michael Neunert, Abbas Abdolmaleki, Markus Wulfmeier, Thomas Lampe, Tobias Springenberg, Roland Hafner, Francesco Romano, Jonas Buchli, Nicolas Heess, and Martin Riedmiller.  Continuous-discrete reinforcement learning for hybrid control in robotics.  In _Conference on Robot Learning_ , pp. 735–751. PMLR, 2020.
  * Omidvar et al. (2015) Mohammad Nabi Omidvar, Xiaodong Li, and Ke Tang.  Designing benchmark problems for large-scale continuous optimization.  _Information Sciences_ , 316:419–436, 2015.
  * Piccoli (1999) Benedetto Piccoli.  Necessary conditions for hybrid optimization.  In _Proceedings of the 38th IEEE Conference on Decision and Control (Cat. No. 99CH36304)_ , volume 1, pp. 410–415. IEEE, 1999.
  * Poupart et al. (2006) Pascal Poupart, Nikos Vlassis, Jesse Hoey, and Kevin Regan.  An analytic solution to discrete bayesian reinforcement learning.  In _Proceedings of the 23rd international conference on Machine learning_ , pp. 697–704, 2006.
  * Rardin & Rardin (1998) Ronald L Rardin and Ronald L Rardin.  _Optimization in operations research_ , volume 166.  Prentice Hall Upper Saddle River, NJ, 1998.
  * Sultana et al. (2020) Nazneen N Sultana, Hardik Meisheri, Vinita Baniwal, Somjit Nath, Balaraman Ravindran, and Harshad Khadilkar.  Reinforcement learning for multi-product multi-node inventory management in supply chains.  _arXiv preprint arXiv:2006.04037_ , 2020.
  * Sun et al. (2019) Shiliang Sun, Zehui Cao, Han Zhu, and Jing Zhao.  A survey of optimization methods from a machine learning perspective.  _IEEE transactions on cybernetics_ , 50(8):3668–3681, 2019.
  * Sutton & Barto (1999) Richard S Sutton and Andrew G Barto.  Reinforcement learning: An introduction.  _Robotica_ , 17(2):229–235, 1999.
  * Sutton & Barto (2018) Richard S Sutton and Andrew G Barto.  _Reinforcement learning: An introduction_.  MIT press, 2018.
  * Tani et al. (2021) Laurits Tani, Diana Rand, Christian Veelken, and Mario Kadastik.  Evolutionary algorithms for hyperparameter optimization in machine learning for application in high energy physics.  _The European Physical Journal C_ , 81:1–9, 2021.
  * Wang et al. (2021) Feng Wang, Heng Zhang, and Aimin Zhou.  A particle swarm optimization algorithm for mixed-variable optimization problems.  _Swarm and Evolutionary Computation_ , 60:100808, 2021.
  * Wang & Zheng (2001) Ling Wang and Da-Zhong Zheng.  An effective hybrid optimization strategy for job-shop scheduling problems.  _Computers & Operations Research_, 28(6):585–596, 2001.
  * Wang et al. (2017) Ling Wang, Ji Pei, Muhammad Ilyas Menhas, Jiaxing Pi, Minrui Fei, and Panos M Pardalos.  A hybrid-coded human learning optimization for mixed-variable optimization problems.  _Knowledge-Based Systems_ , 127:114–125, 2017.
  * Wu et al. (2019) Jia Wu, Xiu-Yun Chen, Hao Zhang, Li-Dong Xiong, Hang Lei, and Si-Hao Deng.  Hyperparameter optimization for machine learning models based on bayesian optimization.  _Journal of Electronic Science and Technology_ , 17(1):26–40, 2019.
  * Young et al. (2015) Steven R Young, Derek C Rose, Thomas P Karnowski, Seung-Hwan Lim, and Robert M Patton.  Optimizing deep learning hyper-parameters through an evolutionary algorithm.  In _Proceedings of the workshop on machine learning in high-performance computing environments_ , pp. 1–5, 2015.
  * Zhan et al. (2022) Zhi-Hui Zhan, Lin Shi, Kay Chen Tan, and Jun Zhang.  A survey on evolutionary computation for complex continuous optimization.  _Artificial Intelligence Review_ , pp. 1–52, 2022.