# Understanding Trade-offs in Stellarator Design with Multi-objective
Optimization

David Bindel1  Matt Landreman2  Misha Padidar3<EMAIL_ADDRESS>1Department
of Computer Science, Cornell University, Ithaca, NY 14853, USA 2Institute for
Research in Electronics and Applied Physics, University of Maryland, College
Park MD 20742, USA 3Center for Applied Mathematics, Cornell University,
Ithaca, NY 14853, USA

###### Abstract

In designing stellarators, any design decision ultimately comes with a trade-
off. Improvements in particle confinement, for instance, may increase the
burden on engineers to build more complex coils, and the tightening of
financial constraints may simplify the design and worsen some aspects of
transport. Understanding trade-offs in stellarator designs is critical in
designing high performance devices that satisfy the multitude of physical,
engineering, and financial criteria. In this study we show how multi-objective
optimization (MOO) can be used to investigate trade-offs and develop insight
into the role of design parameters. We discuss the basics of MOO, as well as
practical solution methods for solving MOO problems. We apply these methods to
bring insight into the selection of two common design parameters: the aspect
ratio of an ideal magnetohydrodynamic equilibrium, and the total length of the
electromagnetic coils.

##  1 Problem Description

Design criteria for stellarators stem from many physical phenomena as well as
many disciplines. For instance, designs must consider fusion performance
metrics from multiple models and scales, engineering specifications that
determine the realizability of the design, and financial timelines which set
the pace of research and construction. The collection of these criteria give
rise to highly constrained design optimization problems with many, potentially
competing, objectives. Any design decision ultimately comes with a trade-off.
For instance, improvement in particle confinement may increase the burden on
engineers to build more complex coils, and tightening of financial constraints
may simplify the design and worsen some aspects of transport. When making
design decisions, such as how to weight or balance objectives in an
optimization or how tightly to enforce a constraint, it is critical to
understand how each choice effects other aspects of the design.

To understand trade-offs that appear in stellarator design we turn to the use
of multi-objective optimization (MOO) (Ehrgott, 2005; Miettinen, 2012;
Emmerich & Deutz, 2018). MOO treats problems of the form,

$\displaystyle\min_{\mathbf{x}\in\bm{\Omega}}\mathbf{F}(\mathbf{x}):=(F_{1}(\mathbf{x}),\ldots,F_{{n_{\mathbf{F}}}}(\mathbf{x})),$
(1)

where $\mathbf{x}\in\mathbb{R}^{{n_{\mathbf{x}}}}$ is a vector of design
variables, $\bm{\Omega}\subseteq\mathbb{R}^{{n_{\mathbf{x}}}}$ is a compact
set, and $F_{1},\ldots,F_{{n_{\mathbf{F}}}}$ are, typically differentiable,
objectives. MOO explores the Pareto optimal space of solutions – those which
are neither better nor worse than one another, and provide a distinct trade-
off in the value of objectives – as well as develops an understanding of how a
reduction in one objective may require worsening of another.

In this study we introduce multi-objective optimization in the context of
stellarator design. We discuss the basics of MOO, as well as a practical
solution method for solving MOO problems: the $\bm{\epsilon}$-constraint
method. We also a present continuation method which leverages a local
expansion of the Pareto front to explore it efficiently. We apply these MOO
methods to bring insight into the selection of two common design parameters:
the aspect ratio of an ideal magnetohydrodynamic (MHD) equilibrium, and the
total length of the electromagnetic coils. The aspect ratio has long been
considered to have a trade-off with the degree to which a configuration is
quasi-symmetric. A conjecture by (Garren & Boozer, 1991) states that quasi-
symmetry can only be achieved through second order in the inverse aspect
ratio, and hence that decreasing aspect ratio would result in a worsening of
the degree of quasi-symmetry. Our numerical experiments suggest that this is
indeed the case, but the trade-off is modest. We also explore the relationship
between the total allowable length of the coil pack, and the ability for the
coils to reproduce a target magnetic field. There is a natural trade-off here
since longer more complex coils can reproduce more intricate fields than
shorter coils.

The paper is structured as follows: in Section 2 we review MOO, as well as two
methods for solving MOO problems. Subsequently, in Section 3 we apply the
optimization methods to explore two trade-off problems in stellarator design.
Finally, in Section 4 we look beyond our case studies to the more general use
of MOO for stellarator design and discuss future directions.

We use bold characters, such as $\mathbf{x}$, to denote vectors. We denote a
vector $\mathbf{x}$ with entry $j$ removed as
$\overline{\mathbf{x}}^{j}=(x_{1},\ldots,x_{j-1},x_{j+1},\ldots,x_{n_{\mathbf{x}}})$.
We compare vectors using vector inequalities: if $\mathbf{x}\leq\mathbf{y}$
then $x_{j}\leq y_{j}$ for $j=1,\ldots,{n_{\mathbf{x}}}$.

##  2 Multi-Objective Optimization

Multi-objective optimization problems typically do not have a single solution,
but rather an entire set of solutions. These are points which balance a trade-
off between objectives: improvement in one objective implies a worsening of
the other. The set of solutions to a MOO problem is called the Pareto front.
In this section, we briefly formalize this notion, before introducing
practical solution techniques for MOO problems.

Figure 1: (Left) Visualization of two points on the efficient set (black
lines). (Right) Pareto front is indicated as the thicker black lines on the
edge of the co-domain $\mathbf{F}(\bm{\Omega})$.

The idea of a trade-off can be formalized through the notion of non-dominance.
A point $\mathbf{x}\in\bm{\Omega}$ dominates $\mathbf{y}\in\bm{\Omega}$, with
respect to the objectives $\mathbf{F}$, if
$\mathbf{F}(\mathbf{x})\leq\mathbf{F}(\mathbf{y})$. Informally, $\mathbf{x}$
is at least as good as $\mathbf{y}$ in all of the objectives. If $\mathbf{x}$
is not dominated by any other point in $\bm{\Omega}$, then $\mathbf{x}$ is
called efficient, or Pareto Optimal, and $\mathbf{F}(\mathbf{x})$ is a non-
dominated point. The set of all efficient points forms the efficient set, and
the set of all non-dominated points forms the non-dominated set, which is more
commonly known as the Pareto front. Finding the Pareto front is the goal of
multi-objective optimization. A visualization of an efficient set and Pareto
front are shown in Figure 1. Intuitively, efficient points are those which
cannot be compared by their objective values alone. For instance, $\mathbf{x}$
may have less complex coil shapes than $\mathbf{y}$, but may not reproduce a
target magnetic field as well. Based off of these facts alone, neither
$\mathbf{x}$ nor $\mathbf{y}$ is a dominant configuration, and additional
information would be needed to determine which configuration is preferable.

In the absence of convex objectives, finding globally Pareto efficient
solutions is difficult. Instead, we settle for searching for weakly, locally
efficient points, which is all most MOO algorithms can guarantee in this
setting, given a finite sample size. $\mathbf{x}\in\bm{\Omega}$ is weakly
efficient (Ehrgott, 2005) if there is no $\mathbf{y}\in\bm{\Omega}$ such that
$\mathbf{y}$ strictly dominates $\mathbf{x}$,
$\mathbf{F}(\mathbf{y})<\mathbf{F}(\mathbf{x})$. Moreover,
$\mathbf{x}\in\bm{\Omega}$ is locally efficient (Emmerich & Deutz, 2018) if
there exists a non-empty open ball containing $\mathbf{x}$ such that
$\mathbf{F}(\mathbf{x})$ is non-dominated in the intersection of the ball and
$\bm{\Omega}$. We find that finding weakly, locally efficient solutions is
enough to provide insight into the trade-offs of interest.

MOO problems are most commonly solved by either directly applying MOO
algorithms (Chang & Wild, 2023; Knowles, 2006; Daulton et al., 2022; Wagner et
al., 2010; Deb et al., 2002) or by reformulating the MOO problem as a series
of scalar optimization problems which can be solved with scalar optimization
methods. Algorithmic approaches attempt to find and explore the efficient set
by taking steps to jointly minimize a combination of the objectives. These
methods are a practical option for problems with many objectives, since the
algorithm can weigh how to efficiently explore the high dimensional Pareto
front. Scalarization methods, on the other hand, allow for more user
involvement and incorporation of problem specific information, such as
derivatives, which many MOO algorithms do not use. Scalarization methods can
be formulated and solved efficiently because they rely on scalar optimization
methods, and allow the user to control the exploration of the efficient set. A
particularly popular scalarization approach is the linearization method
Ehrgott (2005). The linearization method finds efficient points by solving the
linearized problem,

$\min_{\mathbf{x}\in\bm{\Omega}}\ \
\sum_{i=1}^{{n_{\mathbf{F}}}}\alpha_{i}F_{i}(\mathbf{x}),$ (2)

where the user-selected weight vector, $\bm{\alpha}$, must be non-negative and
sum to one. The linearization method, however, has two major drawbacks (Das &
Dennis Jr, 1996): 1) if the Pareto front is non-convex, then no weight vector
$\bm{\alpha}$ exists such that the solution to the linearized problem, eq. 2,
lies on the non-convex portion of the Pareto front, 2) a uniformly spaced
selection of weight vectors do not uniformly explore the Pareto front. For
these reasons, it is worthwhile to look beyond the linearization method when
solving non-convex MOO problems.

In the remainder of this section we discuss a scalarization method, the
$\bm{\epsilon}$-constraint method, which overcomes the drawbacks of the
linearization method and is particularly useful for trade-offs in stellarator
design. We also discuss a continuation method for locally exploring the Pareto
front in an efficient manner.

###  2.1 The $\epsilon$-constraint method

Multi-objective optimization problems can be transformed into a series of
scalar optimization problems via scalarization methods. One such scalarization
method, the $\bm{\epsilon}$-constraint method (Ehrgott, 2005), finds (weakly)
Pareto optimal points by minimizing a single objective subject to upper bound
inequality constraints on all other objectives, with an upper bound parameter
$\bm{\epsilon}\in\mathbb{R}^{{n_{\mathbf{F}}}-1}$. To find points on the
Pareto front, the $\bm{\epsilon}$-constraint method solves problems of the
form

$\displaystyle\begin{split}\min_{\mathbf{x}\in\bm{\Omega}}\ \
&F_{j}(\mathbf{x})\\\ \text{s.t.}\ \
&F_{1}(\mathbf{x}),\ldots,F_{j-1}(\mathbf{x}),F_{j+1}(\mathbf{x}),\ldots
F_{m}(\mathbf{x})\leq\bm{\epsilon},\end{split}$ (3)

where the user select the index $j$ and vector $\bm{\epsilon}$. The
$\bm{\epsilon}$-constraint method theoretically guarantees that it can be used
to find any point on the Pareto front: if $\mathbf{x}$ is a (local) solution
to eq. 3, then it is weakly (locally) efficient. Furthermore $\mathbf{x}$ is
efficient if and only if it is a solution to eq. 3 for all choices of
$j\in\\{1,\ldots,{n_{\mathbf{F}}}\\}$ (Ehrgott, 2005).

The $\bm{\epsilon}$-constraint method is particularly relevant to stellarator
optimization since we often have have bounds on the range of the objectives,
which allows us to easily select the upper bound parameter, $\bm{\epsilon}$.
For example, devices often have an aspect ratio between $3$ and $10$, a coil
must be at least as long as the minor circumference of the plasma, and it
should not be longer than a few times that. By selecting $\bm{\epsilon}$ we
can target areas of the Pareto front that we can search. This is more
interpretable than using other scalarization techniques, such as a weighted
sum of objectives, which requires selecting abstract weights that do not have
any clear relation to the Pareto front.

###  2.2 Continuation methods

Continuation methods (Hillermeier et al., 2001; Schütze et al., 2005;
Vasilopoulos et al., 2021; Peitz et al., 2019; Gkaragkounis et al., 2018) in
MOO explore the Pareto front locally around a given efficient point. By using
optimality conditions and the Implicit Function Theorem (Krantz & Parks, 2002)
continuation methods build local models of the Pareto front, which they use to
estimate nearby efficient points. Local exploration of the Pareto front with
continuation methods comes at a relatively low computational cost when
compared to restarting a scalarization solver, such as the
$\bm{\epsilon}$-constraint method, from scratch. However, when combined,
continuation methods and scalarization solvers make a powerful pair: the
scalarization methods find efficient points over distinct parts of the Pareto
front and the low cost continuation method is used to fill in the space
between points. In this section, we discuss a predictor-corrector continuation
method based off of the $\bm{\epsilon}$-constraint method, which first
predicts an estimate of an efficient point using a Taylor expansion, and
subsequently uses the $\bm{\epsilon}$-constraint method to correct the
predictions, see Algorithm 1. While other predictor-corrector continuation
schemes (Schütze et al., 2005; Vasilopoulos et al., 2021; Peitz et al., 2019)
write the expansion in terms of weights, or target points, we expand the
Pareto front in terms of $\bm{\epsilon}$ because it is easily interpreted, and
allows for the $\bm{\epsilon}$-constraint method to be used in the corrector
step.

Input: Efficient point $\mathbf{x}_{0}$, small value
$\Delta\bm{\epsilon}\in\mathbb{R}^{{n_{\mathbf{F}}}-1}$, index
$j\in\\{1,\ldots,{n_{\mathbf{F}}}\\}$.

Result: Approximate efficient points $\\{\mathbf{x}_{k}\\}$.

1 for  _$k=0,1,\ldots$_ do

2 Set $\bm{\epsilon}_{k}=\overline{\mathbf{F}}^{j}(\mathbf{x}_{k})$;

3 Compute the Lagrange multiplier $\lambda_{k}$ from eq. 9 using
$\bm{\epsilon}=\bm{\epsilon}_{k}$;

4 Compute the Hessian approximations for $\mathbf{F}$ at $\mathbf{x}_{k}$;

5 Build the Jacobian $\mathbf{D}(\mathbf{x}_{k},\lambda_{k})$ from the
Hessians approximations and gradients;

6 Solve for $\nabla_{\bm{\epsilon}}\mathbf{x}$ using equation eq. 8;

7 Predict: Approximate $\mathbf{x}(\bm{\epsilon}_{k}+\Delta\bm{\epsilon})$
using eq. 4;

8 Correct: Compute $\mathbf{x}_{k+1}$ with the $\bm{\epsilon}$-constraint
method by solving eq. 3 using
$\bm{\epsilon}=\bm{\epsilon}_{k}+\Delta\bm{\epsilon}$, and starting the
optimization from the predicted point;

9

Algorithm 1 $\bm{\epsilon}$-Constraint Predictor-Corrector Method

The key idea of the predictor step is to notice that the
$\bm{\epsilon}$-constraint method allows us to parametrically describe weakly
efficient points in terms of $\bm{\epsilon}$, i.e. we can write weakly
efficient points as $\mathbf{x}(\bm{\epsilon})$. In fact, as we will show,
under reasonable conditions $\mathbf{x}(\bm{\epsilon})$ is a continuously-
differentiable map. The predictor step leverages this parametric
representation to approximate the efficient set by a Taylor expansion. Given
an efficient point $\mathbf{x}$, which is a solution to eq. 3 with parameter
$\bm{\epsilon}$, the predictor step locally approximates the efficient set by

$\mathbf{x}(\bm{\epsilon}+\Delta\bm{\epsilon})\approx\mathbf{x}(\bm{\epsilon})+\nabla_{\bm{\epsilon}}\mathbf{x}(\bm{\epsilon})^{T}\Delta\bm{\epsilon},$
(4)

where $\Delta\bm{\epsilon}$ is a small change in $\bm{\epsilon}$, and
$\nabla_{\bm{\epsilon}}\mathbf{x}$ is the Jacobian of
$\mathbf{x}(\bm{\epsilon})$ with respect to $\bm{\epsilon}$. A visualization
of this method is shown in Figure 2. In the remainder of this section we will
discuss conditions which determine when this expansion exists and a method of
computing the Jacobian $\nabla_{\bm{\epsilon}}\mathbf{x}$.

Figure 2: Visualization of the predictor-corrector continuation method.

The existence of the expansion, eq. 4, and the computation of the Jacobian
$\nabla_{\bm{\epsilon}}\mathbf{x}$ rely on applying the Implicit Function
Theorem (Krantz & Parks, 2002) to the first order necessary conditions for
optimality of eq. 3, the Karush-Kuhn-Tucker (KKT) conditions (Nocedal &
Wright, 1999). The KKT conditions are necessary conditions for a point
$\mathbf{x}$ to be an optimal solution to an inequality constrained
optimization problem, such as the $\bm{\epsilon}$-constraint problem eq. 3. In
the case of the $\bm{\epsilon}$-constraint method, points satisfying the KKT
conditions for eq. 3 are also weakly, locally efficient.

Suppose that $\mathbf{x}$ is an optimal solution to eq. 3, and
$\mathcal{A}\subseteq\\{1,\ldots,j-1,j+1,\ldots{n_{\mathbf{F}}}\\}$ is the set
of $m$ active constraints, i.e. the set of constraints $i$ for which
$F_{i}(\mathbf{x})=\epsilon_{i}$. The KKT conditions with the Linear
Independence Constraint Qualification (LICQ) (Nocedal & Wright, 1999) state
that there exist $\bm{\lambda}\in\mathbb{R}^{m}$ such that $\bm{\lambda}\geq
0$ and,

$\displaystyle\nabla F_{j}(\mathbf{x})+\sum_{i\in\mathcal{A}}\lambda_{i}\nabla
F_{i}(\mathbf{x})$ $\displaystyle=0,$ (5)
$\displaystyle\lambda_{i}(F_{i}(\mathbf{x})-\bm{\epsilon}_{i})$
$\displaystyle=0\quad\forall i\in\mathcal{A},$ (6) $\displaystyle
F_{i}(\mathbf{x})$ $\displaystyle\leq\bm{\epsilon}\quad\forall i\neq j.$ (7)

Notice that eqs. 5 and 6 show that if a constraint is active and
$\lambda_{i}>0$ then the solution of the $\bm{\epsilon}$-constraint method is
dependent on $\epsilon_{i}$ in the sense that a small change in $\epsilon_{i}$
could change the solution to the system of equations. In this case, it is
reasonable to write $\mathbf{x}$ as a function of
$\bm{\epsilon}_{\mathcal{A}}=(\epsilon_{i})_{i\in\mathcal{A}}$ to stress the
dependence of solutions on the parameters,
$\mathbf{x}(\bm{\epsilon}_{\mathcal{A}})$. The Implicit Function Theorem gives
us the mathematical leeway to do so.

Informally, the Implicit Function Theorem states that if there exists a point
$\mathbf{x}^{*}$, with Lagrange multipliers $\bm{\lambda}^{*}$ and constraint
parameters $\bm{\epsilon}_{\mathcal{A}}^{*}$, which satisfy the KKT
conditions, and if the Jacobian,
$\mathbf{D}(\mathbf{x}^{*},\bm{\lambda}^{*})$, of the left-hand side of eqs. 5
and 6 with respect to $(\mathbf{x},\bm{\lambda})$ is invertible, then
$\mathbf{x}(\bm{\epsilon}_{\mathcal{A}}),\
\bm{\lambda}(\bm{\epsilon}_{\mathcal{A}})$ are continuously differentiable
functions of $\bm{\epsilon}_{\mathcal{A}}$ around
$\bm{\epsilon}_{\mathcal{A}}^{*}$. Taking the form of a continuously
differentiable function, the space of solutions to the
$\bm{\epsilon}$-constraint problem eq. 3 can be approximated with the first
order Taylor expansion, eq. 4. The Jacobian $\nabla_{\bm{\epsilon}}\mathbf{x}$
can be computed by solving the system,

$\mathbf{D}(\mathbf{x}^{*},\bm{\lambda}^{*})\begin{bmatrix}\nabla_{\bm{\epsilon}}\mathbf{x}\\\
\nabla_{\bm{\epsilon}}\bm{\lambda}\end{bmatrix}=\begin{bmatrix}\mathbf{0}\\\
\bm{\lambda}^{*}\end{bmatrix}$ (8)

The system eq. 8 is derived by treating $\mathbf{x}$ and $\bm{\lambda}$ as a
function of $\bm{\epsilon}$ and differentiating eqs. 5 and 6 with respect to
$\bm{\epsilon}$. Importantly, $\mathbf{D}$ depends on second derivatives of
the objective functions. If $\mathbf{D}$ is invertible, then then the Inverse
Function Theorem guarantees existence of the expansion eq. 4, and a unique
solution exists to eq. 8 (Krantz & Parks, 2002).

Invertibility of the Jacobian matrix $\mathbf{D}$ occurs naturally when
$\mathbf{x},\bm{\lambda}$ satisfy the KKT conditions, and $\mathbf{x}$ is a
strict local minimum: $\nabla^{2}F_{j}(\mathbf{x})+\sum_{i\neq
j}\lambda_{i}\nabla^{2}F_{i}(\mathbf{x})$ is positive definite over the
orthogonal subspace to $\\{\nabla F_{i}\\}$ for all $i\in\mathcal{A}$ such
that $\lambda_{i}>0$. On the other hand, the expansion does not exist if there
are no active constraints: if all $F_{i}(\mathbf{x})<\epsilon_{i}$ then a
small variation in any $\bm{\epsilon}_{i}$ will not effect the solution.
Furthermore, if some constraint $i$ is active but the associated Lagrange
multiplier $\lambda_{i}=0$, then a small increase in $\epsilon_{i}$ will not
change the solution to eq. 3, but a small decrease would. In this case,
$\lambda_{i}$ is only sub-differentiable (Clarke, 1990) with respect to
$\epsilon_{i}$. While an expansion can be derived in this case, we ignore it
for simplicity, as this case only appears at the edges of the Pareto front. We
are primarily concerned with the case when all active constraints have stictly
positive Lagrange multipliers, since in this case there is a trade-off to be
made between objectives: reduction in $F_{j}$ comes at the cost of increase in
$\epsilon_{i}$.

Now that we have established when a local expansion of the efficient set
exists, and how to compute it, we can discuss the predictor-corrector method.
The predictor step of the predictor-corrector method takes as input an
efficient point $\mathbf{x}$ and a step size $\Delta\bm{\epsilon}$.
Subsequently, we compute the Lagrange multiplier $\bm{\lambda}$. If some
components of $\bm{\lambda}$ are strictly positive, then we compute
$\mathbf{D}(\mathbf{x},\bm{\lambda})$, solve eq. 8 for
$\nabla_{\bm{\epsilon}}\mathbf{x}$, then predict a weakly efficient point via
the expansion, eq. 4.

For small enough $\Delta\bm{\epsilon}$, the prediction will typically be a
very good estimate of an efficient point. Nonetheless, it is helpful to finely
resolve the point, particularly as this allows us to compute the subsequent
predictor step accurately. The corrector step corrects the prediction by
simply solving the $\bm{\epsilon}$-constraint problem, eq. 3, with new
constraint parameters $\bm{\epsilon}+\Delta\bm{\epsilon}$. This method is
detailed in Algorithm 1.

####  2.2.1 Computation of the Lagrange Multiplier

In practice, we rarely have a point $\mathbf{x}$ that exactly satisfies any
constraint with equality. Typically, an inequality constraint may
approximately be satisfied with equality, i.e.
$\overline{F}^{j}(\mathbf{x})=\bm{\epsilon}-\nu$ for a small value of $\nu>0$.
In this case, theoretically the constraint is not active and so the Lagrange
multipliers are equal to zero. However, the numerics suggest that the
constraint is active and the Lagrange multiplier should not be zero. In
practice, we determine the active set $\mathcal{A}$ as the set of all
constraints with $\nu_{i}<\tau$ for some small tolerance $\tau$. The Lagrange
multiplier can then be computed by solving the KKT conditions for a perturbed
version of eq. 3 where the $\epsilon_{i}$ is shifted to $\epsilon_{i}-\nu$ for
all constraints in $\mathcal{A}$. In this way, the constraints hold with exact
equality, since $\overline{F}^{j}(\mathbf{x})=\bm{\epsilon}-\nu$, and we are
justified in computing the potentially non-zero Lagrange multipliers. The
Lagrange multipliers are then estimated as the solution to

$\min_{\bm{\lambda}\geq 0}\|\nabla
F_{j}(\mathbf{x})+\sum_{i\in\mathcal{A}}\lambda_{i}\nabla
F_{i}(\mathbf{x})\|^{2}.$ (9)

Importantly, this procedure affects the prediction step in that the expansion
is taken around $\bm{\epsilon}+\nu$. However, this is accounted for in step 1
of Algorithm 1.

##  3 Numerical Experiments

In this section we use MOO to bring insight into the selection of two common
stellarator design parameters: the aspect ratio of an ideal MHD equilibrium,
and the total length of the electromagnetic coils.

Our first experiment attempts to determine if there is a trade-off between
achieving precise quasi-symmetry in a stage-one stellarator design (an ideal
MHD equilibrium), and having low aspect ratio. The Garren-Boozer conjecture
(Garren & Boozer, 1991) suggests that exact quasi-symmetry is only possible at
high aspect ratio. In addition, precise quasi-symmetry has been achieved
throughout a volume in high aspect ratio stellarators (Landreman & Paul, 2022;
Wechsung et al., 2022c; Giuliani et al., 2022b; Landreman, 2022). However, it
is unclear at what rate quasi-symmetry decays as the aspect ratio is
increased, or if this trade-off applies to precise quasi-symmetry at all. Our
first experiment answers the following question: To what extent does the
aspect ratio limit the degree to which quasi-symmetry can be achieved
throughout a volume?

Our second experiment considers a trade-off in the stage-two design problem,
where coils are optimized to fit a target magnetic field. A target magnetic
field can be recreated arbitrarily well by coils which have no constraint on
their length. However, when restricted to have a reasonably short length for
engineering purposes, coils may not be able to reproduce a target magnetic
field. In this problem, we aim to understand how reduction in the allowable
coil length worsens the reproduction of the target magnetic field.

###  3.1 Problem 1: the aspect ratio and quasi-symmetry trade-off

We seek a plasma boundary shape, parameterized by $\mathbf{x}$, of an Ideal
MHD equilibrium for a quasi-helical (QH) stellarator configuration that has
minimal aspect ratio $A(\mathbf{x})$ and deviation from quasi-symmetry
$Q_{M,N}(\mathbf{x})$,

$\min_{A(\mathbf{x})\in[A_{l},A_{u}]}(A(\mathbf{x}),Q_{M,N}(\mathbf{x})).$
(10)

Including bound constraints on the aspect ratio $A_{l}=3,A_{u}=10$, restricts
our decision space to a realistic range of configurations. For convenience, we
collect the ${n_{\mathbf{F}}}=2$ objectives into the vector
$\mathbf{F}(\mathbf{x}):=(A(\mathbf{x}),Q_{M,N}(\mathbf{x}))$.

Figure 3: Pareto front of the aspect ratio and quasi-symmetry objectives over
the domain $A_{l}\leq A(\mathbf{x})\leq A_{u}$. The configurations
corresponding to the blue star, orange square and green diamond markers are
plotted in Figure 4. We find that the efficient set undergoes a branch change
as the aspect ratio crosses the large gap from $\approx 6.11$ to $\approx
8.5$.

The violation of quasi-symmetry is defined as the Quasi-symmetry Ratio
Residual objective employed in (Landreman & Paul, 2022). The objective
measures the departure from quasi-symmetry throughout the plasma volume as the
sum of flux surface averages across surfaces $(s_{1},\ldots,s_{n_{s}})$,

$Q_{M,N}(\mathbf{x})=\sum_{s_{j}}\left\langle\left(\frac{1}{B^{3}}[(N-\iota
M)\mathbf{B}\times\nabla B\cdot\nabla\psi-(MG-NI)\mathbf{B}\cdot\nabla
B]\right)^{2}\right\rangle.$ (11)

The helicity parmeters, $M,N$, determine the type of quasi-symmetry which
$Q_{M,N}$ measures: $M=1,N=0$ being quasi-axisymmetry, and $M=1,N=\pm
k{n_{\text{fp}}}$ with non-zero integers $k$ being quasi-helical symmetry
(Landreman & Paul, 2022). For our experiments, the helicity parmeters were set
to $M=1,N=-{n_{\text{fp}}}$ for the ${n_{\text{fp}}}=4$ field period magnetic
field. Discretization of the flux surface average over the poloidal and
toroidal angles results in an objective with sum-of-squares structure,
$Q_{M,N}(\mathbf{x})=\sum_{i=1}^{{n_{Q}}}Q_{M,N,i}^{2}(\mathbf{x})$, where
$Q_{M,N,i}$ measures the violation of quasi-symmetry at a point throughout the
volume. For a given plasma boundary, $\mathbf{x}$, we compute the aspect ratio
using the definition in the Variational Moments Equilibrium code (VMEC)
(Hirshman et al., 1986; Hirshman & Whitson, 1983),

$A(\mathbf{x})=R_{\text{major}}(\mathbf{x})/R_{\text{minor}}(\mathbf{x})\quad
R_{\text{minor}}=\sqrt{\overline{C}_{A}(\mathbf{x})/\pi}\quad
R_{\text{major}(\mathbf{x})}=\frac{V(\mathbf{x})}{2\pi^{2}R_{\text{minor}}^{2}(\mathbf{x})}$
(12)

where $\overline{C}_{A}$ is the average cross-sectional area of the surface
and $V$ is the volume enclosed by the surface.

The decision variables, $\mathbf{x}$, are Fourier amplitudes that describe the
shape of the plasma boundary. The plasma boundary is represented in terms of
the standard cylindrical coordinates $(R,\phi,Z)$ where $R,Z$ are
parameterized as a Fourier series in the poloidal and toroidal angles $\theta$
and $\phi$,

$\begin{split}R(\theta,\phi)&=\sum_{n=0}^{{n_{sm}}}R_{0,n}\cos(-n_{\text{fp}}n\phi)+\sum_{m=1}^{{n_{sm}}}\sum_{n=-{n_{sm}}}^{{n_{sm}}}R_{m,n}\cos(m\theta-
n_{\text{fp}}n\phi),\\\
Z(\theta,\phi)&=\sum_{n=1}^{{n_{sm}}}Z_{0,n}\sin(-n_{\text{fp}}n\phi)+\sum_{m=1}^{{n_{sm}}}\sum_{n=-{n_{sm}}}^{{n_{sm}}}Z_{m,n}\sin(m\theta-
n_{\text{fp}}n\phi).\end{split}$ (13)

The number of modes describing the surface ${n_{sm}}=\\{0,1,2,\ldots\\}$ can
be increased to achieve more intricate boundary representations. Field period
symmetry with ${n_{\text{fp}}}$ periods and stellarator symmetry have been
assumed. The major radius $R_{0,0}$ is held fixed throughout the optimization,
to fix scale of the design. The Fourier amplitudes are collected into the
decision variable via $\mathbf{x}=(R_{0,1},\ldots,Z_{0,0},\ldots)$. The total
number of decision variables satisfies
${n_{\mathbf{x}}}=4{n_{sm}}^{2}+4{n_{sm}}$.

Numerical experiments were performed using SIMSOPT (Landreman et al., 2021) to
handle variables, compute objectives, and interface with VMEC which computed
the Ideal MHD equilibria from the plasma boundary representation. The
$\bm{\epsilon}$-constraint method was used to find points along the Pareto
front where quasi-symmetry was set to be the target function for minimization
while the aspect ratio was constrained by $\epsilon$. The bound constraints on
the aspect ratio made specifying $\epsilon$ straight forward: $\epsilon$ was
set to linearly spaced values between $A_{l}$ and $A_{u}$. The
$\bm{\epsilon}$-constraint problem eq. 3 was solved by reformulating the
constrained optimization problem with a quadratic penalty method (Nocedal &
Wright, 1999). The penalty parameter was increased from an initial value of
$1$ by a factor of $10$ at each iteration, and the quadratic penalty
subproblems were solved by applying a Gauss-Newton optimization routine.
SIMSOPT was used to compute forward difference gradients via MPI-based
concurrent function evaluations. To avoid local minima, the number of Fourer
modes, ${n_{sm}}$, was increased iteratively from $1$ to $5$, reaching
${n_{\mathbf{x}}}=120$ variables. The $\epsilon$-constraint problem was solved
with the penalty method after each increase of ${n_{sm}}$.

Figure 4: (Left column) Contour plots of the magnetic field strength in Boozer
coordinates $(\theta,\zeta)$ on four flux surfaces $s=0,0.25,0.5,1.0$ for
three Pareto optimal configurations denoted by a star (aspect $3.5$), square
(aspect $5.6$) and diamond (aspect $8.8$) in Figure 3. (Right column) Three
dimensional renderings of corresponding Pareto optimal designs. The color of
the 3D-configurations indicates the field strength, red is stronger.

To reduce the computational overhead of using the $\epsilon$-constraint method
we used the predictor-corrector method, introduced in Section 2.2, to explore
the Pareto front between solutions of the $\epsilon$-constraint method. For
the high dimensional decision space of interest, the full Hessians
$\nabla^{2}Q_{M,N}$ and $\nabla^{2}A$ are too expensive to compute with finite
differences. Instead we approximated the Hessian $\nabla^{2}Q_{M,N}$ with a
Gauss-Newton Hessian approximation and use a diagonal second-order central
difference approximation to $\nabla^{2}A$.

Figure 3 shows the Pareto front for the problem eq. 10. Figure 4 shows three
dimensional renderings of the three solutions highlighted by the square, star,
and diamond markers in Figure 3, as well as contour plots of the field
strength in Boozer coordinates. From Figure 3 it is clear that low values of
the quasi-symmetry violation, $Q_{1,-4}(\mathbf{x})<10^{-3}$, can be achieved
at all aspect ratios in our range. There is a slight trend indicating that
quasi-symmetry may be achieved more precisely at higher aspect ratios.
Nonetheless, by viewing the contour plots of the magnetic field strength in
Boozer coordinates in Figure 4, we see that configurations from all parts of
the Pareto front have visibly precise quasi-symmetry.

While it seems that precise quasi-symmetry can be achieved at all aspect
ratios considered, it is not clear how well the reduction in quasi-symmetry
translates to an improvement in particle confinement. After all, achieving
quasi-symmetry is a proxy for the true goal of achieving good confinement. To
this end we computed the fraction of alpha particles lost from $10$ distinct
Pareto optimal configurations, with aspect ratios 3, 3.5, 4, 4.5, 4.9, 5.6,
6.0, 8.5, 8.7, 9. To compute the losses we scaled the Pareto optimal
configurations to the ARIES-CS reactor (Najmabadi et al., 2008) scale ($1.7$m
minor radius and $B_{0,0}(s=0)=5.7$T field strength on axis), and traced
$5000$ particles born on the $s=1/4$ flux surface as well $5000$ particles
born on the $s=1/2$ flux surface until a terminal time of $0.1$ seconds.
Particles were traced according the vacuum guiding center equations in Boozer
coordinates and were deemed lost if they crossed the $s=1$ flux surface. The
loss fractions are shown in table 1. Not a single particle born on the $s=1/4$
flux surface was lost from nine of the ten configurations, and the remaining
configuration, with aspect ratio $5.6$, lost only $0.22\%$ of the alpha
particles. A slightly larger fraction of the alpha particles born on the
$s=1/2$ flux surface were lost, between $0\%-2.6\%$ for each configuration.
The loss fraction of alpha particles born on the $s=1/2$ flux surface
increases slightly as the aspect ratio increases up until $A\approx 6$, then
drops to approximately zero on the right tail of the Pareto front. It seems
that there is no clear trend between the aspect ratio of the Pareto optimal
designs and the loss fractions for particles born on either of the two
surfaces. Thus while there is a slight trade-off between quasi-symmetry and
aspect ratio, the trade-off between confinement and aspect ratio may be more
complex.

|  | Aspect Ratio
---|---|---
|  | 3.0 | 3.5 | 4.0 | 4.5 | 4.9 | 5.6 | 6.0 | 8.5 | 8.7 | 9.0
Surface $s$ | $1/4$ | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0022 | 0.0 | 0.0 | 0.0 | 0.0
$1/2$ | 0.010 | 0.015 | 0.018 | 0.021 | 0.023 | 0.026 | 0.026 | 0.0002 | 0.0002 | 0.0
Table 1: Fraction of alpha particles lost for ten Pareto optimal
configurations from Figure 3. The loss fractions were computed by sampling
$5000$ particles on the $s=1/4$ flux surface or $s=1/2$ flux surface and
evolving their trajectories by the vacuum guiding center equations in Boozer
coordinates until a terminal time of $0.1$ seconds or the particle breached
the $s=1$ surface and was lost.

Spatially, we may expect or hope that the set of desirable configurations
forms a compact region. However, we find the contrary – the efficient set is
not connected. Qualitatively, the large gap in the Pareto front in Figure 3
hints that the efficient set is undergoing a branch change. Quantitatively, we
measure the branch change by evaluating the Lagrange multipliers at edges of
the gap in the Pareto front (roughly $A=6.11,8.5$) in Figure 3; a branch
change is indicated by $\lambda=0$, since this implies that the aspect ratio
can be changed without causing any change in the quasi-symmetry. Indeed we
find that the Lagrange multipliers are approximately zero, and hence that the
efficient set corresponding to the small piece of the Pareto front, with
aspect ratio greater than or equal to $8.5$, is disconnected from the
remainder of the efficient set. While this may be a product of local
optimization or the representation of decision variables, it is interesting to
see nonetheless that the set of desirable configurations is disparate.

All of the Pareto optimal configurations found achieve exceptional levels of
quasi-symmetry and particle confinement. From this it is clear that the aspect
ratio is not limiting the device performance. On the other hand, engineering
criteria on coil shape, fabrication tolerances, and placement tolerances may
exhibit a greater trade-off with particle confinement, since the magnetic
field generated by coils may not align well with the target magnetic field
computed by Ideal MHD. In the next section we consider the trade-off between
coil length, and the ability for coils to reproduce a target magnetic field.

###  3.2 Problem 2: The coil length and quadratic flux trade-off

Given the shape of the last close flux surface, $\mathcal{S}$, the stage-two
stellarator optimization problem seeks to find magnetic coils with magnetic
field $\mathbf{B}$, such that the field is orthogonal to the surface normal,
$\mathbf{B}\cdot\mathbf{n}=0$ (Merkel, 1987; Zhu et al., 2017). Coils with a
tight constraint on their length, will not be able to form intricate shapes
and reduce the normal component of the magnetic field to zero, whereas coils
which are longer can form arbitrary shapes that better minimize
$\mathbf{B}\cdot\mathbf{n}$. Long coils, on the other hand, are undesirable
from a financial and engineering standpoint: longer coils are more expensive
to build, often have higher curvature, and are more difficult to fit into the
space around the device. In this experiment, we seek to understand the trade-
off between coil length and the ability of coils to reproduce a target
magnetic field.

Figure 5: (Left) Pareto front for coil length vs quadratic flux. The $x$-axis
is non-dimensionalized by dividing by the minor circumference of the surface
($2\pi$ times the minor radius). (Right) Three dimensional rendering of the
Pareto optimal coil sets denoted by the blue star (smaller blue coils) and the
orange square (larger orange coils).

We develop insight into this trade-off by solving the multi-objective problem

$\min_{\mathbf{x}\in\bm{\Omega}}(J_{\mathbf{n}},L),$ (14)

where $J_{\mathbf{n}}$ is the quadratic flux,

$J_{\mathbf{n}}(\mathbf{x})=\frac{1}{2}\int_{\mathcal{S}}(\mathbf{B}\cdot\mathbf{n})^{2}ds,$
(15)

$L$ is the total length of the ${n_{C}}$ coils considered,

$L=\sum_{i=1}^{n_{C}}L_{i},$ (16)

and the decision variables $\mathbf{x}$ represent coil shape parameters and
the coil currents. The quadratic flux is a standard metric used to find coils
in the stage-two problem (Merkel, 1987; Zhu et al., 2017; Landreman, 2017;
Singh et al., 2020; Kruger et al., 2021; Glas et al., 2022; Wechsung et al.,
2022b). To solve the bi-objective coil problem, eq. 14, we use the coil
optimization framework implemented in SIMSOPT. The modular coils are described
by a Fourier representation of the Cartesian coordinates with ${n_{cm}}$
modes. The Fourier representation of the $x$-coordinate of the $ith$ coil is,

$x^{i}(t,\mathbf{x})=x_{c,0}^{i}+\sum_{k=1}^{{n_{cm}}}[x_{c,k}^{i}\cos(kt)+x_{s,k}^{i}\sin(kt)],$
(17)

with analogous forms for the $y$ and $z$ coordinates. Each coil is
additionally equipped with a current, though the current of the first coil is
held fixed at $I_{1}=10^{5}$ Amperes throughout the optimization. The number
of Fourier amplitudes per coil is $3(2{n_{cm}}+1)$ making the total number of
design variables ${n_{\mathbf{x}}}={n_{C}}3(2{n_{cm}}+1)+{n_{C}}-1$. The coils
considered satisfy field-period symmetry, as well as stellarator symmetry.
Thus, the quadratic flux is only computed over a half-field-period, the total
coil length is only taken over the coils in a half-field-period, and the
decision variables $\mathbf{x}$ only represent the ${n_{C}}$ coils in a half-
field-period. The plasma boundary shape $\mathcal{S}$ used in this experiment
is that of the ${n_{\text{fp}}}=2$ field-period quasi-axisymmetric(QA)
configuration from Landreman & Paul (2022), which we will refer to as LP-QA.
${n_{cm}}=5$ Fourier modes were used to describe the coil shapes.

The bi-objective problem was solved by applying the $\bm{\epsilon}$-constraint
method, where quadratic flux was the objective and the total coil length was
the constraint. The $\bm{\epsilon}$-constraint problems were solved with a
quadratic penalty method (Nocedal & Wright, 1999). The penalty parameter was
increased from an initial value of $1$ by a factor of $10$ at each iteration,
and the quadratic penalty subproblems were solved by the L-BFGS-B algorithm
(Zhu et al., 1997) with $2000$ iterations or until the norm of the gradient
reached $10^{-6}$.

Figure 5 shows the Pareto front for the bi-objective coil optimization problem
eq. 14 when ${n_{C}}=4$ coils are used per half-field-period. Figure 5 also
shows three dimensional renderings of the two Pareto optimal coils sets,
highlighted as the star and diamond in the plot of the Pareto front. Figure 5
shows that increasing the coil length substantially improves the ability for
coils to reduce the quadratic flux up until the average coil length over minor
circumference of the plasma boundary is about $5.5$, after which we
numerically find no improvement. In addition we find that when the average
coil length over minor circumference of the plasma boundary exceeds roughly
$5$, the coils become exceedingly complex: coils curvature becomes large and
coils begin to pass under one another, competing for space. At these coil
lengths, coil curvature and the coil-to-coil separation should also be
constrained.

Figure 6: (Left) Trade-off between quasi-axisymmetry and coil length for the
Pareto optimal points shown in Figure 5. For each Pareto optimal coil set, the
quasi-axisymmetry metric $Q_{1,0}$ was computed by running VMEC with a
boundary surfaced obtained by fitting a quadratic flux minimizing surface to
the field generated by the coils. (Middle and right) Contour plots of the coil
generated field strength in the Boozer toroidal and poloidal angles
$\zeta,\theta$ on the $s=1$ surface of the LP-QA configuration for the coil
sets of the blue star and orange square shown in Figure 5.

While Figure 5 shows that constraining the total coil length limits the
ability for coils to reproduce a magnetic field, it does not show how the
limits on the total coil length relate to a loss in quasi-symmetry or loss in
confinement. In Figure 6 we measure the extent to which the Pareto optimal
coil configurations from Figure 5 generate quasi-symmetric magnetic fields. To
do so, we fit a quadratic flux minimizing (QFM) surface (Dewar et al., 1994)
to the coil-generated magnetic field $\mathbf{B}$, evaluate VMEC using the QFM
surface as the boundary shape, and compute the quasi-symmetric metric
$Q_{1,0}$. The QFM surface is constrained to have identical volume to the LP-
QA target surface. We find the QFM surface by solving the following problem in
SIMSOPT:

$\displaystyle\min_{S}\frac{\int_{S}(\mathbf{B}\cdot\mathbf{n})^{2}ds}{\int_{S}\|\mathbf{B}\|^{2}ds},$
(18) s.t. $\displaystyle\text{Volume}(S)=\text{Volume}(\text{LP-QA}).$ (19)

Figure 6 shows that the quasi-axisymmetry degrades significantly, by three
orders of magnitude, as the coil length constraints are tightened. Unlike the
longer coil sets, the shorter coil sets generate significant coil ripple,
degrading the quasi-symmetry. While it may not be possible to entirely avoid
the loss of quasi-symmetry when coil constraints are tightened, the loss may
be diminished by using single-stage optimization approaches to design the
coils. Part of the loss of quasi-symmetry is inherent to the two-stage
approach to coil design: the coils are optimized to reproduce a target
magnetic field, but not to generate a quasi-symmetric field. Single-stage
optimization approaches (Giuliani et al., 2022a, b; Jorge et al., 2023) on the
other hand, directly optimize coils for quasi-symmetry, ensuring that they
achieve optimal quasi-symmetry levels over the space of coils satisfying the
coil length constraint. Single stage approaches would hence improve the
unfortunate trade-off between coil length and quasi-symmetry.

##  4 Discussion

Understanding trade-offs in stellarator designs is particularly important in
designing high performance devices that satisfy the multitude of physical,
engineering, and financial criteria. Throughout this study we have shown how
multi-objective optimization can be used to investigate trade-offs and develop
insight into the role of design parameters.

The example problems considered here were biobjective problems, but many
stellarator design problems may have three or more competing objectives. MOO
methods are useful in this context, however the dimension of the Pareto front
grows with the number of objectives which makes visualizing the Pareto front
difficult and thoroughly exploring the Pareto front often intractable. For
this reason, we generally recommend exploring trade-offs between two or three
objectives at a time. For three objectives, scalarization methods (Ehrgott,
2005; Miettinen, 2012) and other gradient-based MOO methods (Fliege & Svaiter,
2000; Désidéri, 2012; Schäffler et al., 2002; Peitz & Dellnitz, 2018;
Gkaragkounis et al., 2018) are an appropriate choice since they allow the use
of derivatives and for the user to dicate which part of the Pareto front is
explored. For problems with more than three objectives, parallel algorithmic
approaches will find promising configurations most efficiently (Chang & Wild,
2023; Knowles, 2006; Daulton et al., 2022; Wagner et al., 2010; Deb et al.,
2002).

Looking beyond the two trade-offs considered here, there are a host of other
trade-offs that appear in stellarator design that deserve attention. There is
a natural trade-off, for instance, between coil complexity and quasi-symmetry
that should be explored using a direct single state optimization method, such
as with the near-axis expansion (Wechsung et al., 2022a). In addition, it is
not well understood how stability criteria trade off with particle confinement
criteria, or how design flexibility for multi-purpose coils trades off with
volume and quasi-symmetry (Lee et al., 2022). Ideally, these problems should
be solved and the Pareto optimal solutions should be tabulated in a way that
practitioners can easily survey a multitude of configurations and analyze
their strengths and weaknesses as a holistic set.

##  5 Data Availability

Code and data can be found at https://doi.org/10.5281/zenodo.7838063.

##  6 Acknowledgements

This work was generously supported by a grant from the Simons Foundation (No.
560651, D.B.).

## References

  * Chang & Wild (2023) Chang, Tyler H & Wild, Stefan M 2023 Parmoo: A python library for parallel multiobjective simulation optimization. Journal of Open Source Software 8 (82), 4468.
  * Clarke (1990) Clarke, Frank H 1990 Optimization and nonsmooth analysis. SIAM.
  * Das & Dennis Jr (1996) Das, Indraneel & Dennis Jr, John E 1996 A closer look at drawbacks of minimizing weighted sums of objectives for pareto set generation in multicriteria optimization problems. Tech. Rep..
  * Daulton et al. (2022) Daulton, Samuel, Eriksson, David, Balandat, Maximilian & Bakshy, Eytan 2022 Multi-objective bayesian optimization over high-dimensional search spaces. In Uncertainty in Artificial Intelligence, pp. 507–517. PMLR.
  * Deb et al. (2002) Deb, Kalyanmoy, Pratap, Amrit, Agarwal, Sameer & Meyarivan, TAMT 2002 A fast and elitist multiobjective genetic algorithm: Nsga-ii. IEEE transactions on evolutionary computation 6 (2), 182–197.
  * Désidéri (2012) Désidéri, Jean-Antoine 2012 Multiple-gradient descent algorithm (mgda) for multiobjective optimization. Comptes Rendus Mathematique 350 (5-6), 313–318.
  * Dewar et al. (1994) Dewar, RL, Hudson, SR & Price, PF 1994 Almost invariant manifolds for divergence-free fields. Physics Letters A 194 (1-2), 49–56.
  * Ehrgott (2005) Ehrgott, Matthias 2005 Multicriteria optimization, , vol. 491. Springer Science & Business Media.
  * Emmerich & Deutz (2018) Emmerich, Michael TM & Deutz, André H 2018 A tutorial on multiobjective optimization: fundamentals and evolutionary methods. Natural computing 17 (3), 585–609.
  * Fliege & Svaiter (2000) Fliege, Jörg & Svaiter, Benar Fux 2000 Steepest descent methods for multicriteria optimization. Mathematical methods of operations research 51, 479–494.
  * Garren & Boozer (1991) Garren, DA & Boozer, Allen H 1991 Existence of quasihelically symmetric stellarators. Physics of Fluids B: Plasma Physics 3, 2822–2834.
  * Giuliani et al. (2022a) Giuliani, Andrew, Wechsung, Florian, Cerfon, Antoine, Stadler, Georg & Landreman, Matt 2022a Single-stage gradient-based stellarator coil design: Optimization for near-axis quasi-symmetry. Journal of Computational Physics 459, 111147.
  * Giuliani et al. (2022b) Giuliani, Andrew, Wechsung, Florian, Landreman, Matt, Stadler, Georg & Cerfon, Antoine 2022b Direct computation of magnetic surfaces in boozer coordinates and coil optimization for quasi-symmetry. arXiv preprint arXiv:2203.03753 .
  * Gkaragkounis et al. (2018) Gkaragkounis, K, Papoutsis-Kiachagias, E, Asouti, V, Giannakoglou, K & Unit, Optimization 2018 Adjoint-based pareto front tracing in aerodynamic shape optimization. In 10th international conference on computational fluid dynamics (ICCFD10), Barcelona, Spain.
  * Glas et al. (2022) Glas, Silke, Padidar, Misha, Kellison, Ariel & Bindel, David 2022 Global stochastic optimization of stellarator coil configurations. Journal of Plasma Physics 88 (2), 905880208\.
  * Hillermeier et al. (2001) Hillermeier, Claus & others 2001 Nonlinear multiobjective optimization: a generalized homotopy approach, , vol. 135\. Springer Science & Business Media.
  * Hirshman et al. (1986) Hirshman, SP, Merkel, P & others 1986 Three-dimensional free boundary calculations using a spectral green’s function method. Computer Physics Communications 43 (1), 143–155.
  * Hirshman & Whitson (1983) Hirshman, Steven P & Whitson, JC 1983 Steepest-descent moment method for three-dimensional magnetohydrodynamic equilibria. The Physics of fluids 26 (12), 3553–3568.
  * Jorge et al. (2023) Jorge, R, Goodman, A, Landreman, M, Rodrigues, J & Wechsung, F 2023 Single-stage stellarator optimization: Combining coils with fixed boundary equilibria. arXiv preprint arXiv:2302.10622 .
  * Knowles (2006) Knowles, Joshua 2006 Parego: A hybrid algorithm with on-line landscape approximation for expensive multiobjective optimization problems. IEEE Transactions on Evolutionary Computation 10 (1), 50–66.
  * Krantz & Parks (2002) Krantz, Steven George & Parks, Harold R 2002 The implicit function theorem: history, theory, and applications. Springer Science & Business Media.
  * Kruger et al. (2021) Kruger, Thomas G, Zhu, C, Bader, A, Anderson, DT & Singh, L 2021 Constrained stellarator coil curvature optimization with focus. Journal of Plasma Physics 87 (2), 175870201\.
  * Landreman (2017) Landreman, Matt 2017 An improved current potential method for fast computation of stellarator coil shapes. Nuclear Fusion 57 (4), 046003.
  * Landreman (2022) Landreman, Matt 2022 Mapping the space of quasisymmetric stellarators using optimized near-axis expansion. Journal of Plasma Physics 88, 905880616.
  * Landreman et al. (2021) Landreman, Matt, Medasani, Bharat, Wechsung, Florian, Giuliani, Andrew, Jorge, Rogerio & Zhu, Caoxiang 2021 Simsopt: A flexible framework for stellarator optimization. Journal of Open Source Software 6 (65), 3525.
  * Landreman & Paul (2022) Landreman, Matt & Paul, Elizabeth 2022 Magnetic fields with precise quasisymmetry for plasma confinement. Physical Review Letters 128 (3), 035001.
  * Lee et al. (2022) Lee, Brandon F, Paul, Elizabeth J, Stadler, Georg & Landreman, Matt 2022 Stellarator coil optimization supporting multiple magnetic configurations. Nuclear Fusion 63 (1), 014002\.
  * Merkel (1987) Merkel, Peter 1987 Solution of stellarator boundary value problems with external currents. Nuclear Fusion 27, 867.
  * Miettinen (2012) Miettinen, Kaisa 2012 Nonlinear multiobjective optimization, , vol. 12. Springer Science & Business Media.
  * Najmabadi et al. (2008) Najmabadi, F, Raffray, AR, Abdel-Khalik, SI, Bromberg, L, Crosatti, L, El-Guebaly, L, Garabedian, PR, Grossman, AA, Henderson, D, Ibrahim, A & others 2008 The aries-cs compact stellarator fusion power plant. Fusion Science and Technology 54 (3), 655–672.
  * Nocedal & Wright (1999) Nocedal, Jorge & Wright, Stephen J 1999 Numerical optimization. Springer.
  * Peitz & Dellnitz (2018) Peitz, Sebastian & Dellnitz, Michael 2018 Gradient-based multiobjective optimization with uncertainties. In NEO 2016: Results of the Numerical and Evolutionary Optimization Workshop NEO 2016 and the NEO Cities 2016 Workshop held on September 20-24, 2016 in Tlalnepantla, Mexico, pp. 159–182. Springer.
  * Peitz et al. (2019) Peitz, Sebastian, Ober-Blöbaum, Sina & Dellnitz, Michael 2019 Multiobjective optimal control methods for the navier-stokes equations using reduced order modeling. Acta Applicandae Mathematicae 161 (1), 171–199.
  * Schäffler et al. (2002) Schäffler, Stefan, Schultz, Reinhart & Weinzierl, Klaus 2002 Stochastic method for the solution of unconstrained vector optimization problems. Journal of Optimization Theory and Applications 114, 209–222.
  * Schütze et al. (2005) Schütze, Oliver, Dell’Aere, Alessandro & Dellnitz, Michael 2005 On continuation methods for the numerical treatment of multi-objective optimization problems. In Dagstuhl Seminar Proceedings. Schloss Dagstuhl-Leibniz-Zentrum für Informatik.
  * Singh et al. (2020) Singh, Luquant, Kruger, TG, Bader, A, Zhu, C, Hudson, SR & Anderson, DT 2020 Optimization of finite-build stellarator coils. Journal of Plasma Physics 86 (4), 905860404.
  * Vasilopoulos et al. (2021) Vasilopoulos, Ilias, Asouti, Varvara G, Giannakoglou, Kyriakos C & Meyer, Marcus 2021 Gradient-based pareto front approximation applied to turbomachinery shape optimization. Engineering with Computers 37 (1), 449–459.
  * Wagner et al. (2010) Wagner, Tobias, Emmerich, Michael, Deutz, André & Ponweiser, Wolfgang 2010 On expected-improvement criteria for model-based multi-objective optimization. In Parallel Problem Solving from Nature, PPSN XI: 11th International Conference, Kraków, Poland, September 11-15, 2010, Proceedings, Part I 11, pp. 718–727. Springer Berlin Heidelberg.
  * Wechsung et al. (2022a) Wechsung, Florian, Giuliani, Andrew, Landreman, Matt, Cerfon, Antoine & Stadler, Georg 2022a Single-stage gradient-based stellarator coil design: stochastic optimization. Nuclear Fusion 62 (7), 076034.
  * Wechsung et al. (2022b) Wechsung, Florian, Giuliani, Andrew, Landreman, Matt, Cerfon, Antoine & Stadler, Georg 2022b Stochastic and a posteriori optimization to mitigate coil manufacturing errors in stellarator design. Plasma Physics and Controlled Fusion 64 (10), 105021.
  * Wechsung et al. (2022c) Wechsung, Florian, Landreman, Matt, Giuliani, Andrew, Cerfon, Antoine & Stadler, Georg 2022c Precise stellarator quasi-symmetry can be achieved with electromagnetic coils. Proceedings of the National Academy of Sciences 119 (13), e2202084119.
  * Zhu et al. (1997) Zhu, Ciyou, Byrd, Richard H, Lu, Peihuang & Nocedal, Jorge 1997 Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization. ACM Transactions on mathematical software (TOMS) 23 (4), 550–560.
  * Zhu et al. (2017) Zhu, Caoxiang, Hudson, Stuart R, Song, Yuntao & Wan, Yuanxi 2017 New method to design stellarator coils without the winding surface. Nuclear Fusion 58 (1), 016008.


# HashMem : PIM-based Hashmap Accelerator

Akhil Shekar  Morteza Baradaran  Sabiha Tajdari  Kevin Skadron

###### Abstract

Hashmaps are widely utilized data structures in many applications to perform a
probe on key-value pairs. However, their performance tends to degrade with the
increase in the dataset size, which leads to expensive off-chip memory
accesses to perform bucket traversals associated with hash collision. In this
work, we propose HashMem, a processing-in-memory (PIM) architecture designed
to perform bucket traversals along the row buffers at the subarray level. Due
to the inherent parallelism achieved with many concurrent subarray accesses
and the massive bandwidth available within DRAM, the execution time related to
bucket traversals is significantly reduced. We have evaluated two versions of
HashMem, performance-optimized and area-optimized, which have a speedup of
$49.1x/17.1x$ and $9.2x/3.2x$ over standard C$++$ map and hyper-optimized
hopscotch map implementations respectively.

###### keywords:

Processing In Memory, HashMaps, In-Situ Computing, DRAM, Memory Systems

##  1 Introduction

As we move further into the digital age, the amount of data being generated
and consumed daily is increasing at an unprecedented rate. One of the popular
data structures for searching large datasets is a hashmap, due to its near-
constant-time lookup performance. For large datasets, the hashmaps usually are
not cache-resident, and hence, any lookup would be an expensive off-chip DRAM
access to read the entire hash chain to perform a single probe.

By supporting hashmap lookup directly in memory, i.e. using a processing-in-
memory (PIM) architecture, applications can avoid costly memory transfers
between the processor and memory, leading to faster and more efficient
lookups. Hashmaps are particularly well-suited to PIM architectures, as they
are able to leverage the parallel processing capabilities of the memory to
perform lookups and hash-chain traversals in constant time.

For a brief background, a hashmap is a data structure that allows efficient
storage and retrieval of key-value pairs. It uses a hash function to compute
an index into an array of buckets or slots from which the desired value can be
found. In other words, a hashmap is an associative array that maps keys to
values.

The computation of hashmaps is memory-bound because it heavily relies on
accessing memory when the dataset cannot fit within the cache. When searching
for a key in a hashmap, the hash function is used to compute the index of the
bucket where the key is stored. The bucket is then accessed in memory to
retrieve the value associated with the key. As the size of the hashmap
increases, the number of memory accesses required to perform a lookup also
increases, making it a memory-bound operation.

Processing-in-memory is a new architecture paradigm that breaks down the
memory wall by integrating processing elements within the memory chips
themselves. This eliminates the need for data to be transferred between memory
and processor, leading to significant improvements in performance and energy
efficiency [9].

In this paper, we propose HashMem, a PIM architecture designed to accelerate
key-value probes on hashmaps. Our architecture comprises two versions: area-
optimized and performance-optimized. Both place processing elements adjacent
to each subarray in the DRAM, but the area-optimized provides one processing
unit per subarray and operates on one value at a time, i.e., element-serial,
bit-parallel; while the performance-optimized provides multiple processing
units per subarray, operating on the entire row at once in an element-parallel
but bit-serial fashion. Thus, in the latter case, the values are laid out in a
column-oriented fashion, so that each row contains a single-bit slice from
thousands of values, achieving high parallelism at the expense of requiring
$b$ steps in order to find a b-bit key.

These organizations are evaluated against standard C$++$ map and hyper-
optimized hopscotch map implementations and found to yield significant
speedups compared to a server-grade CPU: $17.1x/3.2x$ (area-optimized) and
$49.1x/9.2x$ (performance optimized) over standard C$++$ map and hyper-
optimized hopscotch map implementations respectively. The in-situ processing
with PIM also minimizes energy spent on moving data, thus achieving energy
savings in addition to the performance benefit. Quantifying the energy savings
is left for future work.

##  2 HASHMEM ARCHITECTURE

We have designed and implemented a PIM architecture that leverages the lower
access latency and inherent parallelism available within the DRAM structure at
the subarray level to perform hashmap lookups at the subarray interface. Our
key idea behind the design involves mapping an entire hash bucket to a
subarray row within PIM memory. Each subarray row contains between 512-2048
columns, where a column is defined in terms of a multi-bit access length, with
each such column being 4, 8 or 16 bits in length. Internally each of these
subarrays could be broken down into mats for implementation purposes. However,
we consider a subarray to be a set of mats that are activated in parallel for
our understanding of the rest of the paper.

Since an entire hash bucket happens to be mapped to a subarray row, the entire
bucket is activated into the row buffer when the subarray row is accessed. We
designed processing elements (PEs) to sit at the edge of each subarray closer
to the row buffer to perform lookups for necassary keys within the activated
hash bucket. Each PE consists of (i) comparison unit to perform the comparison
operation, (ii) control unit or logic to orchestrate and control the
operations and (iii) output register to hold the resultant value associated
with a matched key during the lookup operation. In our architecture, we
propose two implementations of these comparison units, an area-optimized
version and a performance-optimized version.

Figure 1: Area-optimized HashMem Architecture

###  2.1 Area-optimized version

The area-optimized version of HashMem performs the hash bucket traversal in an
element-serial bit-parallel manner. Fig. 1 demonstrates the architecture for
the area-optimized version. The PE accesses each of the activated hash
bucket’s key-value pairs sequentially from the row buffer and looks for a key
match. Upon a match, the corresponding value would be stored within its output
register to be read out later by the RLU.

###  2.2 Performance-optimized version

As shown in Fig. 3, the performance-optimized version relies on placing many
small comparison units below the row buffer to scan all the keys in parallel
within a single or small number of clock ticks. Compared to the Area-optimized
version, this offers higher performance with lower execution time but suffers
from increased area overhead since many comparison units must be placed and
pitch-mapped along the subarray row buffer. This version resembles the
operation of Content Addressable Memory (CAM) except that this is operating on
a DRAM row buffer at the subarray-level.

Figure 2: RLU mounted at rank-level

###  2.3 Rank-Level Unit (RLU)

The HashMem architecture also involves the usage of a Rank-level Unit (RLU),
which acts as an intermediary orchestrating agent or command processor between
the in-situ PIM processing elements and the host processor (CPU). The job of
the RLU is to : (i) Propagate the key to be searched to the necessary subarray
(ii) Orchestrate probing operations compliant with the DRAM timing parameters
and architecture constraints (iii) Retrieve the output values after the
probing operation is completed from the subarray units and buffer them before
transferring them to the memory controller.

The RLU helps in the overall integration with the rest of the system by
abstracting the PIM operations and interfacing with the Memory Controller (MC)
with special PIM-capable extended DRAM commands. This ensures the host can
support the PIM capabilities with minimal changes to its integrated memory
controllers. The RLU is responsible for communicating with the in-situ
subarray level PIM elements and orchestrating operations amongst them. It is
analogous to the command processor (CP) that exists within GPUs to interface
with its Streaming Multiprocessors (SMs) or Compute Units (CUs) and the PCIe
bus. As shown in figure 2, since the RLU is mounted as a separate chip, the
logic area overhead of the RLU does not affect the memory capacity of the
DIMM, as observed in similar rank-level modifications done with AxDIMM [4].

Figure 3: Performance-optimized HashMem Architecture

###  2.4 Virtualization

The CPU operates on physical addresses in the typical fashion. Hence, in order
to make the PIM memory systems compatible with the current virtualization
scheme, we rely on storing hash buckets at page granularity. With this,
irrespective of the location of page within physical memory, when a hash
bucket is accessed, the corresponding page containing the bucket is activated
and its related subarray processing elements are enabled to perform the lookup
operations. In scenarios where a page is co-located with other pages in the
same row buffer, the page start and end addresses are communicated to the
subarray PEs to access only the necessary address range. Using this scheme,
the traditional virtualisation is still supported without the need for CPU to
concern itself with physical data placement to orchestrate PIM operations.

###  2.5 Probing

Once the initial dataset is populated within the PIM memory, the CPU
communicates the key to be probed in the respective hash bucket. The page
table translation helps in locating the necessary rank (RLU) and subarray row
that holds the hash bucket. The respective RLU receives a compute-capable DRAM
command that informs it of the input key to be probed and the address of the
page to be probed. The RLU, in turn, orchestrates the probing operation by
activating the necessary subarray row and communicating to the in-situ
subarray PE to perform the actual probing operation. The RLU later retrieves
the value from the output register of the same subarray-level PE and passes it
to the Memory Controller (MC) in a cache line format. The cache line can be
padded with additional zeroes if the data being transferred is less than the
size of a cache line. The MC, after receiving the data from the RLU, places it
into the requested CPU’s Last Level Cache (LLC) address. Once the CPU reads
and extracts the value from the cache, the probing operation is complete.

Currently, deletion operations involve putting tombstone values at the place
where a key-value pair is deleted, at the cost of wasted space. This is
similar to software implementation of hashmaps. We aim to further investigate
how to perform efficient deletion operations to reclaim the space back for
further usage while using the hashmap on PIM.

Often times, the load distribution of key-value pairs amongst the hash buckets
is not equal and this might lead to some buckets having too many key-value
pairs and some having too few. We ran a test case scenario where we mapped the
first 350,000 words of a dictionary into a hashmap and measured the length of
each bucket. We observed a significant variance in the lengths as shown in
fig. 4 that demonstrated the under and over-utilisation of buckets.

Under-utilized buckets \- If the page size is N bytes of memory, and the
bucket occupies only P bytes of data (where P < N), then the remaining (N-P)
bytes of the page are wasted and lead to inefficient memory usage. The page
size N is dictated at the boot time without any prior knowledge of the
dataset. The value of P (bucket size) is decided during the runtime and
depends on the input dataset. Hence, there are bound to be certain buckets
which are under-utilized and efforts could be made to map and fit two or more
of such buckets into the same page. This helps in reorganising the memory and
improving its utilization. However, care has to be taken to ensure proper
bookkeeping of the relocated buckets. Also, a strict criterion that the bucket
to be relocated is not split and fragmented across multiple pages is to be
followed.

Over-utilized Buckets \- Some hash buckets may exceed the allocated page size
of N bytes and occupy P bytes (where P > N), resulting in an overflow
situation that needs to accommodate an extra (P-N) bytes. In these scenarios,
an extra page is allocated to accommodate the overflow data, and a bookkeeping
structure is updated to record the presence of hash bucket across two or more
pages which helps while performing a lookup. Essentially, having these extra
pages spread across different channels and ranks helps in probing them in
parallel, thereby improving the performance. This optimization could be
introduced into the Memory Management Unit (MMU) to instruct it to spread
pages containing overly-utilized buckets across different channels evenly to
enable the parallel probing of pages. We have marked this as an avenue for
future work with micro-architectural changes to be investigated to introduce
to support for this optimization strategy.

Alternately, several works such as [8, 12], propose ideal hashing functions to
counter this unequal distribution phenomenon when certain prior knowledge of
the dataset is available to us.

Figure 4:  Length of Hashbuckets after mapping first 350,000 words in a
dictionary

###  2.6 Memory Controller (MC) Changes

The Memory Controller needs to have the capability to differentiate between a
conventional READ/WRITE operation and PIM operation. Specifically, it needs to
understand that a particular page being accessed is a hash bucket that needs
traversal/lookup to occur on the PIM side and retrieve only the value
associated with a specific key. This interaction with the MC should be invoked
by the CPU and abstracted from the programmer and exposed as a simple library
call. Towards this end, it could be suggested as an extension to the ISA and
change on the micro-architectural implementation. Furthermore, we also suggest
that the MC have capabilities to communicate with RLU on the other side of the
memory bus using special physical layer (PHY) commands that are pin-compatible
with the existing DDR standards. The RLU present on the DIMM provides the
necessary abstraction to the memory controllers to understand and parse these
special PIM commands and orchestrate the operations with its related in-situ
subarray-level elements.

##  3 Programming Interfaces

We have provided snippets of code for insertion and lookup in the HashMem PIM-
capable memory below. A pseudo-code is listed that provides the necessary
abstraction for the programmer to utilize this PIM-capable memory using a
high-level programming language. We make use of a bookkeeping structure that
keeps track of the hash buckets and the pages that store them.

###  3.1 Insertion Operation

The overall idea behind the insertion operation is to obtain the hash bucket
(or page/pages the bucket is mapped to) that needs to store the key-value pair
by performing a hash function on the key. However, there are several
constraints, as mentioned in Section 2, that need to be looked after,
especially with regard to bucket utilization and page overflow scenarios.

Initially, in Step-➊, we obtain the page size from the system information
using a library call. This information is usually stored in the operating
system after the boot procedure has initialized the page tables and other
virtualization structures. In ➋, we hash the key and decide the bucket and
page in which to store the key-value pair. Consequently, we need to perform a
check if the bucket is going to overflow in step-➌ while inserting the input-
key-value pair. If the bucket is not going to overflow, then the key-value
pair could be inserted successfully, as in Step-➍. However, if the bucket is
about to overflow, we do pim_malloc() to initialize a fresh page in step-➎. We
update a bookkeeping structure to reflect that a particular hash bucket
extends to a new page. This is to ensure any subsequent lookup operations
probe the new page in addition to the prior existing page/(s). Once the
bookkeeping structure is updated, we go ahead and store the key-value pair on
the new page (step-➏) and provide the necessary return code to reflect the
successful insertion of the key-value pair.

Listing 1: PIM Insertion Operation

⬇

void MapInputKeyValuePairToHashMemPage(keyDataType inputKey, valueDataType
inputValue) {

pageSizeInBytes = getPageSizeFromSystemInfo(); ➊

/* This Page Size is decided during the boot time based on subarray
structures. This information is retrieved using the runtime system information
call. */

numOfKVPairsPerPage = pageSizeInBytes / sizeOfEachKVPair;

// Calculating number of KV pairs per page

destinationPage = getHashValueFromHashingAlgorithm(inputkey); ➋

// Perform Hashing operation to decide which page the key-value pair needs to
be mapped to

currentNumOfKVPairsInPage = getCurrentPageSize(destinationPage);

// Get size current page size in terms of how many number of key-value pairs
it is currently holding

/* Perform a check if the page can accomodate the key-value pair or not */

if(currentSizeOfPage < numOfKVPairsPerPage) ➌ {

storeKVPairIntoPage(destinationPage, inputKey, inputValue); ➍

// Store the KV pair into the destinationPage

}

else{ // if page is already full

pim_page newPage;

ret_code = pim_malloc(newPage); ➎

// allocates a new page and assigns to newPage structure

if(ret_code == PR_ERROR)

return PR_ERROR; // Page allocation failed, function exits

updateBookkeepingStructure(newPage, destinationPage); ➏

/* Attaches and links new page to old page (destination_page) in a Linked List
fashion. This book-keeping structure also informs that these two pages (new
page and old page) need to be probed together when looking for keys. Because
the key can reside either in new page or old page. */

storeKVPairIntoPage(newPage, inputKey, inputValue);

// Store the key-value into the newPage

}

return PR_SUCCESS;

}

###  3.2 Lookup operation

The lookup operation is relatively simple and straightforward, where the user
provides an input key and expects the value associated with it to be returned.
In Step-➊, the hashing function computes the bucket and page to be probed
based on the input key provided. In Step-➋, we make use of a special library
call that performs a hash bucket lookup operation. In the background process,
the library call consults the bookkeeping structure to check the number of
pages to probe and instructs the MMU to perform page-level probing operations
that traverses the entire hash bucket using the subarray PIM processing
elements. These PIM probing operations retrieve the value associated with the
key and provide the necessary return code to reflect the successful lookup
operation.

Listing 2: PIM Probe Operation

⬇

valueDataType probeKey(keyDataType inputKey) {

bucketToProbe = getHashValueFromHashingAlgorithm(inputkey); ➊

// Perform hashing to get the hash bucket to probe

/* Variable to hold value associated with key-value pair */

valueDataType outputValue = NULL;

outputValue = pimProbeBucket(bucketToProbe); ➋

/* This first checks the book keeping structure as to how many pages to probe
and then issues a special pim probe command to the memory controller to
perform probing on PIM-capable memory*/

/* If the value was not found */

if (outputValue == NULL)

return PR_ERROR;

return outputValue;

}

##  4 Evaluation

In this section, we will discuss the performance and area overhead of HashMem
architecture and compare it with traditional CPU-based implementation
baselines. The configuration of the hardware setups utilized in our analysis
are as mentioned in the Table-1.

Table 1: Hardware Configuration Property | Value
---|---
Processor Name | Intel(R) Xeon(R) Silver 4208 CPU (2.1 GHz)
Total Cores / Main Memory | 8 (16 threads) / 512 GB
L1D/L2/L3 Cache Size | 32 KB per core/1 MB per core/11.2 MB shared
HashMem | DDR4_8Gb_x16_3200 (Single Channel)
| 8 banks per rank, 128 subarrays per bank
| 512 rows per subarray
Table 2: Workload Overview Property | Value
---|---
Dataset | Contains 100 million key-value pairs (800 MB)
| 10% i.e. 10 million randomly selected keys probed
Points of Comparison | Standard C++ Map (binary tree)
| C++ unordered_map (hashmap)
| hopscotch map (optmized hashmap)

###  4.1 Benchmarking

Currently, there exists no standard benchmark to test exclusively for the
hashmap performance. Hence, we proposed our microbenchmark to test the hashmap
performances on both the CPU and HashMem. The details of the microbenchmark
and baselines we compare against are provided in the proceeding section. In
order to estimate and model the HashMem performance, we analyzed the timing
data gathered from prior works [14, 6, 1, 7].

####  4.1.1 Microbenchmark

One of the goals of the microbenchmark is for the input dataset to be
sufficiently large enough such that it flows out of the cache and into the
DRAM. This ensures that the cache effects of the CPU are sufficiently
eliminated and the expensive off-chip DRAM accesses are captured while
performing lookups during microbenchmark run. Another goal is to generate
random accesses, i.e., accessing random keys to eliminate any prefetching
based on spatial locality.

Towards this end, we propose a microbenchmark consisting of 100 million key-
value (KV) pairs with a key and value occupying $4$ bytes each. The key and
the value are coded as default uint32_t data type in C$++$. Hence, each KV
pair would occupy $8$ bytes of data, and the overall memory footprint of the
input dataset would be $800$ MB in size, sufficiently large enough to overflow
the L3 cache of most processors used in our testing scenario. Furthermore,
$10\%$ of the keys in the input dataset are probed, which equates to $10$
million keys being searched for in the hashmap. The keys to be probed again
are selected at random and fed into both the CPU and PIM architectures
simultaneously to assess the hashmap probing performance.

Although this benchmark does not account for string values and other types of
data, we envision that they could be pre-processed and dictionary-encoded into
numerical values to be used in HashMem. Performing string value comparison or
any regex operations at the subarray-level units would incur a very high area
overhead and is avoided. Hence we consciously supported probing just numerical
values with HashMem PIM architecture.

###  4.2 Performance of different data structures on CPUs

Figure 5: Probing Times of Different Data Structures

Our initial research indicated that “unordered_map” Standard Template Library
(STL) in C++ has the closest implementation resembling a hashmap. However,
there are several other default C++ data structures that perform similar
operations as a hashmap but are implemented using alternative techniques. One
of the most popular examples is the “map” data structure that is implemented
as a red-black tree, a specialized binary search tree data structure. We
intended compare against both the unordered_map and map to assess the
performance.

Apart from the default C++ libraries, there are various implementations that
are more performance-optimized to yield better throughput in performing key-
value pair lookups. Hopscotch is one example that implements a hopscotch
hashing algorithm [3] to resolve hashing collisions [13]. We found a popularly
used repository online [2] that implemented this hashing mechanism that we
could consider as a state-of-the-art baseline. Interestingly, we discovered
that hopscotch_map is faster than Google’s sparse-hash [10] implementation
during our testing. We intended our evaluations to compare with a mix of
candidates, hence our choice of these three software implementations.

The performance results of the three data structures chosen on CPU are
demonstrated above in Figure-5. Hopscotch is significantly faster than the
other implementations by a wide margin, around $5.3x$ and $3.1x$ compared to
C$++$ map and unordered_map, respectively. The map structure performs the
worst considering that it is implemented as a balanced binary search tree with
many indirect accesses to traverse along the tree both during insertions and
probing. The complexity of a map is $log_{2}^{n}$. Interestingly, we found
that unordered_map rehashes itself when the number of elements to be inserted
exceeds the load factor, i.e., the number of buckets available, rather than
having more nodes attached to extend each bucket size.

###  4.3 Area overhead

To estimate the area overhead of subarray-level PIM processing elements, we
first obtain the area breakdown of the DDR4 chip (DDR4_8Gb_x8_3200) that is
used to build the HashMem. Each subarray contains 512 rows, and there are 128
subarrays per bank. The comparator units of HashMem are implemented in RTL and
the delay, area overheads are evaluated using Synopsys DC Compiler in 14nm. We
use scaling factors from [11] to scale the results to 22nm.

####  4.3.1 Performance optimized

For the performance-optimized version, there is a requirement to pitch-map the
comparator units to fit within the section boundaries of the row buffer
containing the key-value pair segments. This presents significant challenges
related to the evaluation of this version and this is part of the future work
to investigate HashMem further.

####  4.3.2 Area optimized

The area-optimized version does not require significant efforts to pitch-map
as was required with the performance-oriented version. Our estimate revealed
that incorporating 64 additional ALUs, shared by 128 subarrays per bank,
incurred only 5.26% area overhead.

##  5 Results

The performance of hashmap workloads are both dependent on the size and the
distribution of the input dataset. The numbers reported are for the evaluation
of the dataset as detailed in section 4.1.1 describing the microbenchmark.

Figure 6 highlights that our area-optimized HashMem version outperforms the
standard map, unordered map, and hopscotch map implementations, achieving
speedup values of approximately $17.1x$, $5.5x$, and $3.2x$ respectively.

Figure 6: HashMem Speedup Against CPU

As shown in figure 6, our performance-optimized HashMem surpasses the standard
map, unordered map, and hopscotch implementations by even greater factors of
$49.1x$, $15.8x$, and $9.2x$ respectively.

In both scenarios, HashMem outperforms the CPU baseline by a wide margin, even
against the state-of-the-art hopscotch implementation. An interesting aspect
of the experiment is that these results represent a single DRAM channel
competing against a server-class CPU. We could parallelize the lookups across
the independent memory channels and obtain further improvement in the HashMem
performance. These results demonstrate the potential of PIM architectures,
which have tremendous intrinsic parallelism and bandwidth that are not being
harnessed by current Von Neumann-bottlenecked architectures.

##  6 FUTURE WORK

Our future work broadly aims at improving the architecture with a richer set
of evaluations. We aim to look at a wider variety of datasets with different
distribution patterns of key and value strings to assess the hashmap
performance in each setting.

Tiered-Latency DRAM. [5] splits the subarray into low-latency and higher-
latency regions by placing isolation transistors to alter the DRAM READ access
timings. We aim to leverage this work to map our hashmap buckets into the low-
latency region for faster lookups. In a conventional architecture, mapping
several elements to a smaller subset of buckets could degrade the performance.
However, due to parallelism provided by the PIM processing elements and the
lower-latency offered by the Tiered-Latency DRAM, it could improve our
performance significantly. Moreover, the lower the number of buckets available
within a hashmap, the higher the probability of a row hit which further leads
to improved performance and decreased latency.

Channel-level Parallelism. Memory channels are independent and parallel READ /
WRITE operations could be performed on each of the channels separately. Hence,
if there are multiple keys to be probed, these probing operations could be
parallelised amongst different channels to increase the throughput. However,
this type of parallelism could be exploited only if the keys being probed
belong to different channels.

Energy Savings Analysis. Due to the vastly reduced number of expensive off-
chip data access over the memory bus, there are significant energy savings to
be realised. This also leads to reduced instruction overhead on the CPU
related to performing a scan operation on the hash bucket traversals.

Data Types. Currently our evaluations only support 32-bit int operations on
both key and value. The limitation is with regards to the configuration of the
subarray-level PIM elements. We aim to investigate re-configurable sizes for
future evaluations to support other data types too.

Hash Function. We aim to investigate an optimum hashing mechanism that can
evenly distribute the input dataset over several buckets of equal or near
equal length. This is to reduce certain buckets from getting over-utilized and
some from getting under-utilized. There are prior works in this area such as
[8].

Real-world kernels. Many genomics, database and other applications extensively
make use of hashmap structures in their kernels while implementing them. We
aim to test our PIM architecture in these settings and observe the performance
improvements at the application level.

##  7 Conclusions

We proposed HashMem, a subarray-level PIM architecture that accelerates
hashmap lookups by leveraging the existing parallelism available within the
DRAM structure. We have demonstrated that the two variants of HashMem, area
and performance-optimized ones, were able to outperform state-of-art hopscotch
hashmap and default C++ map implementations on CPU by at least $3.2x/9.2x$ and
$17.1x/49.1x$ respectively. The area overhead as observed on the area-
optimized version was found to be $5.26\%$. We have laid out several
directions as part of our future work to further investigate potential of
HashMem architecture.

##  8 Acknowledgements

This work was supported in part by PRISM, one of seven centers in JUMP 2.0, a
Semiconductor Research Corporation (SRC) program sponsored by DARPA. We also
thank the reviewers for their helpful feedback and suggestions.

## References

  * [1] D. T. W. Bruce Jacob, Spencer W. NG, _Memory Systems : Cache, DRAM, Disk_. Burlington, MA: Morgan Kaufmann Publishers, 2008.
  * [2] T. Goetghebuer-Planchon. (2015, Mar.) A c++ implementation of a fast hash map and hash set using hopscotch hashing. [Online]. Available: https://github.com/Tessil/hopscotch-map
  * [3] M. Herlihy, N. Shavit, and M. Tzafrir, “Hopscotch hashing,” in _Proceedings of the 22nd International Symposium on Distributed Computing_ , ser. DISC ’08. Berlin, Heidelberg: Springer-Verlag, 2008, p. 350–364. [Online]. Available: https://doi.org/10.1007/978-3-540-87779-0_24
  * [4] L. Ke, X. Zhang, J. So, J.-G. Lee, S.-H. Kang, S. Lee, S. Han, Y. Cho, J. H. Kim, Y. Kwon, K. Kim, J. Jung, I. Yun, S. J. Park, H. Park, J. Song, J. Cho, K. Sohn, N. S. Kim, and H.-H. S. Lee, “Near-memory processing in action: Accelerating personalized recommendation with axdimm,” _IEEE Micro_ , vol. 42, no. 1, pp. 116–127, 2022.
  * [5] D. Lee, Y. Kim, V. Seshadri, J. Liu, L. Subramanian, and O. Mutlu, “Tiered-latency dram: A low latency and low cost dram architecture,” in _2013 IEEE 19th International Symposium on High Performance Computer Architecture (HPCA)_ , 2013, pp. 615–626.
  * [6] D. Levinthal. Performance analysis guide for intel core i7 processor and intel xeon 5500 processors. [Online]. Available: https://www.intel.com/content/dam/develop/external/us/en/documents/performance-analysis-guide-181827.pdf
  * [7] S. Li, Z. Yang, D. Reddy, A. Srivastava, and B. Jacob, “Dramsim3: A cycle-accurate, thermal-capable dram simulator,” _IEEE Comput. Archit. Lett._ , vol. 19, no. 2, p. 106–109, jul 2020. [Online]. Available: https://doi.org/10.1109/LCA.2020.2973991
  * [8] Y. Lu, B. Prabhakar, and F. Bonomi, “Perfect hashing for network applications,” in _2006 IEEE International Symposium on Information Theory_ , 2006, pp. 2774–2778.
  * [9] O. Mutlu, S. Ghose, J. Gómez-Luna, and R. Ausavarungnirun, “A modern primer on processing in memory,” 2022.
  * [10] Sparse-Hash. (2020) sparse-haspmap. [Online]. Available: https://github.com/sparsehash/sparsehash
  * [11] A. Stillmaker, Z. Xiao, and B. M. Baas, “Toward more accurate scaling estimates of CMOS circuits from 180 nm to 22 nm,” Univ. of California-Davis Tech. Report ECE-VCL-2011-4, 2012.
  * [12] Wikipedia. Open-addressing. [Online]. Available: https://en.wikipedia.org/wiki/Open_addressing
  * [13] Wikipedia. (2015, Mar.) Hopscotch hashing is a scheme in computer programming for resolving hash collisions. [Online]. Available: https://en.wikipedia.org/wiki/Hopscotch_hashing
  * [14] Y. Yuan, M. Alian, Y. Wang, R. Wang, I. Kurakin, C. Tai, and N. S. Kim, “Don’t forget the i/o when allocating your llc,” in _2021 ACM/IEEE 48th Annual International Symposium on Computer Architecture (ISCA)_ , 2021, pp. 112–125.