# Globally Stable Neural Imitation Policies

Amin Abyaneh1, Mariana Sosa Guzmán2, and Hsiu-Chin Lin3 1Amin Abyaneh is with
the Department of Electrical and Computer Engineering, McGill University,
Montreal, Canada<EMAIL_ADDRESS>Sosa Guzmán is with the
Department of Electrical and Computer Engineering, Universidad Veracruzana,
Veracruz, Mexico., 3Hsiu-Chin Lin is with the School of Computer Science and
the Department of Electrical and Computer Engineering, McGill University,
Montreal, Canada<EMAIL_ADDRESS>work is sponsored by NSERC
Discovery Grant and Mitacs Globalink Research Internship.

###### Abstract

Imitation learning mitigates the resource-intensive nature of learning
policies from scratch by mimicking expert behavior. While existing methods can
accurately replicate expert demonstrations, they often exhibit
unpredictability in unexplored regions of the state space, thereby raising
major safety concerns when facing perturbations. We propose SNDS, an imitation
learning approach aimed at efficient training of scalable neural policies
while formally ensuring global stability. SNDS leverages a neural architecture
that enables the joint training of the policy and its associated Lyapunov
candidate to ensure global stability throughout the learning process. We
validate our approach through extensive simulations and deploy the trained
policies on a real-world manipulator arm. The results confirm SNDS’s ability
to address instability, accuracy, and computational intensity challenges
highlighted in the literature, positioning it as a promising solution for
scalable and stable policy learning in complex environments.

##  I Introduction

Imitation learning (IL) is a pivotal notion attempting to overcome the diverse
safety and complexity challenges of policy learning by emulating expert
behavior [1, 2]. IL facilitates the training of intricate motion policies
without resorting to an exhaustive search in the robot’s state space [3, 4].
Nonetheless, only a handful of IL methods offer formal guarantees pertaining
to the stability of the resulting policies. Global stability assumes a
critical role when deploying the policy in stochastic environments prone to
external perturbations. Notably, it assures that the policy can recover
effectively and _predictably_ to a predetermined target, even in uncharted
regions of the state space that lie beyond the scope of expert demonstrations.

Figure 1: Overview of the proposed stable neural policy learning method.
Policy learning (top) optimizes a Lyapunov-stable neural policy over the
expert demonstration data. The optimized policy is then deployed (bottom) to
plan globally stable trajectories resistant to unpredictable perturbations.

Prior research in safe IL is focused on autonomous and time-invariant
dynamical systems to model and train motion policies [4]. These policies are
optimized across a set of expert trajectories, yielding a proper action
(velocity) contingent on the current state (position). The optimization
process ensures global stability through established theoretical frameworks,
such as Lyapunov or contraction theories [5, 6, 7, 8, 9, 10, 11].
Nevertheless, reproducing complex trajectories in high-dimensional state
spaces with the previous methods can be impractical due to computational
complexity, non-convex optimization, and sample inefficiency. To address
computational cost, previous research [7, 9, 11] tends to restrict the class
of functions of the contraction metric or the Lyapunov candidate, which in
turn limits the imitation performance of the optimized policy. Most notably, a
common limitation of these methods appears in the absence of neural policy
representation, which are recognized for their scalability, efficient
gradient-based optimization, and domain transfer capabilities [12].

Recently developed stable neural policies [13, 14] mainly employ
diffeomorphism and invertible mappings to transform a simple stable policy
into a highly nonlinear one. Yet, the training stage requires multiple
demonstrations, leads to policies with quasi-stability, and hinders
incremental learning for slightly different expert trajectories. Neural IL
methods, on the other hand, offer flexible and precise policies [15, 16, 17,
18], but lack the required safety and reliability outside the region covered
by expert demonstrations and are often data hungry. Similar limitations
persist in inverse reinforcement learning [19, 20], where computational
complexity is fueled by having reinforcement learning in an inner loop [15].

To tackle these challenges, we represent policies with Stable Neural Dynamical
Systems (SNDS), where we jointly train a neural dynamical system policy,
alongside a secondary, convex-by-design network that guarantees global
stability. Subsequently, a joint optimization process trains the neural policy
by minimizing a novel differentiable loss aimed at strengthening the alignment
between policy rollouts and expert demonstrations. SNDS, therefore, benefits
from an expressive and stable neural representation, allowing for safe and
scalable approximation of the underlying dynamical system. We present an
overview of our framework in Figure 1, and outline our key contributions
below.

  * •

Designing SNDS—a stable-by-design neural representation for nonlinear
dynamical systems as motion policies based on expert demonstrations.

  * •

Providing formal stability analysis founded on Lyapunov theory and convex
neural networks.

  * •

Formulating a differentiable trajectory alignment loss function, inspired by
forward Euler’s method.

  * •

Empirical evaluation of SNDS’s effectiveness in higher state space dimensions
for complex trajectories, both in simulation and real-world scenarios.

##  II Background

###  II-A Preliminaries

Dynamical system.  A dynamical system (DS) describes the evolution of a state,
$\mathbf{{x}}$, in the ambient space $\mathcal{X}\subset\mathbb{R}^{n}$, over
time [21]. We consider an autonomous, time-invariant DS modeled with a first-
order ordinary differential equation, $\dot{\mathbf{{x}}}=f(\mathbf{{x}})$,
where the system yields the time derivative for each $\mathbf{{x}}$, without
any control input and independent of time.

Lyapunov stability theory.  $f(\mathbf{{x}})$ is globally asymptotically
stable (GAS) at an equilibrium, $\mathbf{{x}}^{*}$, if for any initial state,
the system approaches $\mathbf{{x}}^{*}$ as time progresses towards infinity.
The Lyapunov stability theory is a widely used tool to analyze the GAS
property of dynamical systems. According to the theory, a system exhibits GAS
if there exists a positive-definite function
$\psi(\mathbf{{x}}):\mathcal{X}\xrightarrow{}\mathbb{R}$, referred to as the
Lyapunov potential function (LPF), satisfying $\dot{\psi}(\mathbf{{x}})<0$ for
all $\mathbf{{x}}\neq\mathbf{{x}}^{*}$, and $\dot{\psi}(\mathbf{{x}}^{*})=0$
at the equilibrium.

Input convex neural networks (ICNN).  ICNN’s special architecture enables
universal approximation of convex functions while maintaining convexity during
training [22, 23]. An $L$-layer ICNN, $\hat{\psi}$, is formulated as follows
for $l\in[L]$:

$\displaystyle\hat{\psi}(\mathbf{{x}})=z_{L},\quad
z_{l+1}=\sigma_{l}(U_{l}x+W_{l}z_{l}+b_{l}),$

where $z_{l}$ denotes the layer output, and $W_{l}$ and $b_{l}$ are real-value
weights. In contrast to feedforward networks, ICNN design exploits $U_{l}$,
that are positive weights connecting the input directly to each layer, and
$\sigma_{l}$ can be any convex, non-decreasing activation function.

###  II-B Problem statement

Consider a policy that functions within a state space denoted by
$\mathcal{X}\subset\mathbb{R}^{n}$. The state space may correspond to a
robot’s task ($\mathit{T}$) or configuration ($\mathit{C}$) space. The policy
outputs an action in $\mathcal{A}\subset\mathbb{R}^{n}$, such as velocity or
torque, which determines the change in the state over time. We denote the
state variable by
$\mathbf{{x}}\triangleq[{x}_{1};{x}_{2};\ldots;{x}_{n}]^{T}\in\mathcal{X}$,
and assume that the yielded action indicates the time derivative of the
current state, denoted by $\dot{\mathbf{{x}}}=\frac{\partial x}{\partial
t}\in\mathcal{A}$. In this context, our primary objective is to learn a
globally stable imitation policy to map $\mathbf{{x}}$ to
$\dot{\mathbf{{x}}}$, provided a dataset of expert’s state-action pairs.

Let $N_{d}\in\mathbb{N}$ be the number of trajectories demonstrated by the
expert. Each trajectory contains $N_{s}\in\mathbb{N}$ state-actions. The
dataset of expert trajectories stacks all state-action pairs,

$\displaystyle\mathbf{D}\triangleq\Bigl{\\{}\big{(}\mathbf{{x}}^{d}[s],\;\dot{\mathbf{{x}}}^{d}[s]\big{)}\;\big{|}\;d\in[N_{d}],\;s\in[N_{s}]\Bigr{\\}},$
(1)

where $(\mathbf{{x}}^{d}[s],\;\dot{\mathbf{{x}}}^{d}[s])$ is the dataset entry
corresponding to the $s$-th sample of the $d$-th trajectory. The dataset
$\mathbf{D}$ holds $N_{t}=N_{d}N_{s}$ samples. We assume that the trajectories
share a common target endpoint, $\mathbf{{x}}^{*}\in\mathcal{X}$, and have
zero velocity at the target, i.e., $\mathbf{{x}}^{d}[N_{s}]=\mathbf{{x}}^{*}$
and $\dot{\mathbf{{x}}}^{d}[N_{s}]=\mathbf{0},\;\forall d\in[N_{d}]$.

The mapping between states and actions can be modeled with a time-invariant
autonomous DS [5], denoted by:

$\displaystyle\dot{\mathbf{{x}}}=f(\mathbf{{x}})+\epsilon=\phi(\mathbf{{x}}),\;\;f,\phi:\mathcal{X}\xrightarrow{}\mathcal{A}.$
(2)

In Equation 2, the function $f$ corresponds to an ordinary differential
equation describing the true underlying DS. The term
$\epsilon\in\mathbb{R}^{n}$ accounts for the effect of measurement noise
present in the expert’s demonstrations, which is incorporated into the
estimated DS. Our subsequent objective is to acquire a noise-free
approximation of $f(\mathbf{{x}})$, denoted as $\phi(\mathbf{{x}})$. This
estimated function, $\phi(\mathbf{{x}})$, can be perceived as a policy that
maps states to actions in real-time, thereby guiding the robot with the right
set of actions to imitate the demonstrated trajectories within the state
space.

Alternatively, a trajectory can be generated through forward Euler’s method.
Imagine a robot situated at $\mathbf{{x}}[0]$, the policy generates an action
$\dot{\mathbf{{x}}}[0]=\phi(\mathbf{{x}}[0])$, and the next state is
subsequently calculated through $\mathbf{{x}}[1]=\mathbf{{x}}[0]+\Delta
t\phi(\mathbf{{x}}[0])$, where $\Delta t$ determines the granularity of the
discretization.

##  III Methodology

We model the policy, $\phi(\mathbf{{x}})$, and the corresponding LPF,
$\psi(\mathbf{{x}})$, with two neural networks in Section III-A. This allows
us to accurately imitate an expert’s behavior, while providing formal GAS
analysis in Section III-B. Lastly, we introduce a differentiable loss function
in Section III-C to improve both sample efficiency and imitation accuracy of
SNDS.

###  III-A Dynamical system policy formulation

Let $\hat{\phi}(\mathbf{{x}}):\mathcal{X}\rightarrow\mathcal{A}$ be a standard
and unrestricted feedforward neural network which models the unstable policy.
Although the model can capture intricate expert demonstrations, training
$\hat{\phi}(\mathbf{{x}})$ solely on $\mathbf{D}$ results in an unstable
policy which cannot predictably recover even from mild perturbations, as
depicted in Figure 2. Hence, in this section, we enforce the trained policy to
have the GAS property by satisfying the conditions of Lyapunov stability
theory in Section II-A.

The first step is to define a positive-definite and differentiable LPF. As
mentioned in Section II-A, ICNNs can approximate any convex function. To
ensure that the LPF is positive definite, we employ a well-studied technique
[17, 18] to define it as,

$\displaystyle\psi(\mathbf{{x}})=\big{[}\hat{\psi}(\mathbf{{x}})-\hat{\psi}(\mathbf{{x}}^{*})\big{]}+\delta||\mathbf{{x}}-\mathbf{{x}}^{*}||_{2}^{2},$
(3)

where $\psi(\mathbf{{x}}):\mathcal{X}\rightarrow\mathbb{R}$ is still a convex
function. The first term ensures $\psi(\mathbf{{x}}^{*})=0$, while the
addition of the state’s ${l}_{2}$-norm with a negligible $\delta>0$, along
with convexity of ICNN guarantee that
$\psi(\mathbf{{x}})>0,\;\forall\;\mathbf{{x}}\neq\mathbf{{x}}^{*},\;\mathbf{{x}}\in\mathcal{X}$.
Hence, the LPF with this architecture, $\psi(\mathbf{{x}})$, satisfies the
positive-definite condition required by Lyapunov stability theory.

Next, to satisfy the negative derivative condition of Lyapunov theory, i.e.,
$\frac{d\psi(\mathbf{{x}})}{dt}=\nabla\psi(\mathbf{{x}})^{T}\hat{\phi}(\mathbf{{x}})<0$,
we modify the projection expression in [23] to enforce GAS while training
$\hat{\phi}(\mathbf{{x}})$. The policy $\phi(\mathbf{{x}})$ resulting from
projecting $\hat{\phi}(\mathbf{{x}})$ onto the half space,
$\nabla\psi(\mathbf{{x}})^{T}\hat{\phi}(\mathbf{{x}})<0$, is formulated as
follows,

$\displaystyle\phi(\mathbf{{x}})$
$\displaystyle=\begin{cases}\hat{\phi}(\mathbf{{x}})-\nabla\psi(\mathbf{{x}})\frac{\nabla\psi(\mathbf{{x}})^{T}\hat{\phi}(\mathbf{{x}})}{\lVert\nabla\psi(\mathbf{{x}})\rVert_{2}^{2}}&\nabla\psi(\mathbf{{x}})^{T}\hat{\phi}(\mathbf{{x}})>=0\\\
\hat{\phi}(\mathbf{{x}})&\nabla\psi(\mathbf{{x}})^{T}\hat{\phi}(\mathbf{{x}})<0\end{cases}$
(4)
$\displaystyle=\hat{\phi}(\mathbf{{x}})-\nabla\psi(\mathbf{{x}})\frac{\text{ReLU}(\nabla\psi(\mathbf{{x}})^{T}\hat{\phi}(\mathbf{{x}}))}{\lVert\nabla\psi(\mathbf{{x}})\rVert_{2}^{2}}.$

Figure 2: An example of unstable (left) vs. a stable (right) policies
optimized on expert’s data from the handwriting dataset [24]. While both
policy rollouts can reproduce the expert motion, an unstable policy cannot
recover from perturbation that push the robot to unknown state space regions.

In Equation 4, if the network output $\hat{\phi}(\mathbf{{x}})$ fails to meet
the Lyapunov condition, i.e.,
$\nabla\psi(\mathbf{{x}})^{T}\hat{\phi}(\mathbf{{x}})>=0$, the output is
projected so that $\phi(\mathbf{{x}})$ always fulfills the condition. The
process is simplified with a ReLU activation function. Note that the original
formulation in [23] enforces exponential stability, which is too restrictive
for our problem. Hence, we further relax the projection to guarantee GAS.

Figure 3 portrays an example of trained $\phi(\mathbf{{x}})$ and
$\psi(\mathbf{{x}})$, and the resulting global stability. The differentiable
and globally stable $\phi(\mathbf{{x}})$, can be trained using efficient
gradient-based methods within its parameters, as long as
$\hat{\phi}(\mathbf{{x}})$ and $\psi(\mathbf{{x}})$ are defined using
automatic differentiation tools.

###  III-B Global asymptotic stability guarantees

To establish that $\phi(\mathbf{{x}})$ defines a function with GAS property,
we use proof by contradiction to show that regardless of initial state,
$\mathbf{{x}}[0]$, every trajectory converges to the target, that is,
$\lim_{k\rightarrow\infty}\mathbf{{x}}[k]=\mathbf{{x}}^{*}$.

######  Proposition 1.

The dynamical system, $\phi(\mathbf{{x}})$, in Equation 4 is globally
asymptotically stable with the Lyapunov function, $\psi(\mathbf{{x}})$,
defined in Equation 3, and any two arbitrary networks,
$\hat{\phi}(\mathbf{{x}})$ and $\hat{\psi}(\mathbf{{x}})$, with bounded real-
valued weights.

Figure 3: Joint training of globally stable neural policies (left) and the
corresponding LPF (right) optimized on expert’s data from the handwriting
dataset [24]. The contours for the LPF (green) illustrate both positive-
definiteness and convexity of the trained function.

We leverage proof by contradiction for this proposition. Let
$\\{\mathbf{{x}}[0],\ldots,\mathbf{{x}}[k]\\}$ be a trajectory not converging
to the target, meaning that
$\lim_{k\rightarrow\infty}\mathbf{{x}}[k]\neq\mathbf{{x}}^{*}$. Since
$\psi(\mathbf{{x}})$ is decreasing, as formulated in Equation 4, and
$\psi(\mathbf{{x}})\geq 0$, as defined in Equation 3, there must exist a value
$\mu\in\mathbb{R}^{+}$ such that
$\lim_{k\rightarrow\infty}\psi(\mathbf{{x}}[k])=\mu$. We reason that
$\mu\leq\psi(\mathbf{{x}}[i])\leq\psi(\mathbf{{x}}[0])$ for every $0\leq i\leq
k$.

Now consider the set $\mathbf{S}$ of all $\mathbf{{x}}[i],\forall 0\leq i\leq
k$. This set is _compact_. More precisely, for any open cover 111An open cover
for a set is a collection of open sets whose union contains the original set.
of that set, there exists a finite subcover, which means that a finite number
of the open sets from the cover are sufficient to cover the original set.
Hence, $\dot{\psi}(\mathbf{{x}})$ takes its supremum,
$\sup_{\mathbf{S}}\dot{\psi}(\mathbf{{x}})=-\alpha$ and
$\alpha\in\mathbb{R}^{+}$ over this set. This conclusion is justified based on
the projection introduced in Equation 4. The core integration in the proof of
Lyapunov theorem,

$\displaystyle\psi(\mathbf{{x}}[T])=\psi(\mathbf{{x}}[0])+\int_{0}^{T}\dot{\psi}(\mathbf{{x}}[t])dt\leq\psi(\mathbf{{x}}[0])-\alpha
T,$

indicates a _contradiction_ when,
$T=\frac{\psi(\mathbf{{x}}[0])}{\alpha}+\beta$ and $\beta\in\mathbb{R}^{+}$:

$\displaystyle\psi(\mathbf{{x}}[T])\leq\psi(\mathbf{{x}}[0])-\alpha(\frac{\psi(\mathbf{{x}}[0])}{\alpha}+\beta)=-\alpha\beta<0,$

The above equation contradicts the fact that $\psi(\mathbf{{x}})>0$. Hence,
all trajectories must converge to the desired target, $\mathbf{{x}}^{*}$.
$\square$

######  Remark 2.

Proposition 1 only requires $\psi(\mathbf{{x}})$ to be convex. The policy
network, $\phi(\mathbf{{x}})$, remains fairly unrestricted to capture
intricate expert demonstrations.

Thus, the proposed SNDS formally establishes global stability using an
efficient gradient-based search for the LPF within the class of convex
functions.

###  III-C SRVF training loss

When constructing an optimization loss, the literature often opts for the
conventional Mean Squared Error (MSE) [5, 13, 3] to gauge the dissimilarity
between the policy’s output and the actual velocity labels over each batch of
data. Consequently, the MSE loss is solely focused on the precision of
estimated velocity. While reproducing exact velocities is essential, even
slight inaccuracies at one instance can lead to accumulated error in the
reproduced trajectories. The accumulation happens as a result of using forward
Euler’s method that relies on previous estimates of $\mathbf{{x}}$ to yield
the next ones. Therefore, a more effective loss function should promote the
accuracy of trajectories generated through forward Euler’s method during
policy rollouts.

Employing forward Euler’s method to generate a full trajectory makes the
gradient-based optimization intractable as a result of repetitive
integrations. To mitigate this effect, we propose a limited horizon
combination of the Square-Root Velocity Functions (SRVF) [25] and traditional
MSE to design the loss function. SRVF methods involve normalizing the curve’s
velocity, with a focus on changes in curve shape. Assuming a trajectory
$\mathbf{{x}}^{d}_{k}=\\{\mathbf{{x}}^{d}[0],\ldots,\mathbf{{x}}^{d}[k]\\}$,
we shape the features by the definition:

$\displaystyle\texttt{SRVF}(\mathbf{{x}}^{d}_{k})=\sqrt{\frac{\mathbf{{x}}^{d}_{k}[s+\Delta
t]-\mathbf{{x}}^{d}_{k}[s]}{\|\mathbf{{x}}^{d}_{k}[s+\Delta
t]-\mathbf{{x}}^{d}_{k}[s]\|}},\quad\Delta t=1.$ (5)

Given the state-action pairs,
$\mathbf{{x}}^{d}[s],\;\dot{\mathbf{{x}}}^{d}[s]\sim\texttt{SRVF}(\mathbf{D})$,
and the differentiable dynamics $\phi(\mathbf{{x}})$, we define the loss
$\mathcal{L}$ as,

$\displaystyle\begin{aligned}
&\mathcal{L}(\mathbf{{x}};\;\theta)=\;\gamma_{0}\mathop{\mathbb{E}}\Big{[}\big{(}\phi(\mathbf{{x}}^{d}[s])-\dot{\mathbf{{x}}}^{d}[s]\big{)}^{2}\Big{]}\;+\\\
\sum_{i\in[N_{w}]}\gamma_{i}\mathop{\mathbb{E}}&\Big{[}\big{(}\mathbf{{x}}^{d}[s+i-1]+\phi(\mathbf{{x}}^{d}[s+i-1])\Delta
t-\mathbf{{x}}^{d}[s+i]\big{)}^{2}\Big{]},\end{aligned}$ (6)

where $\gamma_{i}$ are discount factors decreasing as the loss horizon
expands. The intuition behind Equation 6 is rather straightforward:
$\mathcal{L}$ is a differentiable function encapsulating both the current
label mismatch, and deviation from the generated trajectory for a fixed
horizon of $N_{w}$ consecutive samples. As long as the horizon is limited,
stochastic optimizers, such as ADAM [26], can optimize the loss on expert’s
training data.

##  IV Experiments

We utilize two sets of motion planning data in our experiments. Our primary
dataset is sourced from the well-known LASA Handwriting dataset [24], which
contains records of handwritten trajectories on a tablet. The second dataset
consists of more complex motions gathered by [7] in the same way. Our real-
world experiments are only based on the latter due to its complex nature.

###  IV-A Evaluation

Figure 4: Comparing the reproduction accuracy (top) and computational cost
(bottom) of SNDS against baseline methods using the MSE and DTW metrics
introduced in Equation 7 and Equation 8, respectively. The accuracy of policy
actions and the learning time against other _stable_ methods remain comparably
lower for SNDS across most designated trajectories.

We compare our approach against three existing baselines: Behavioral Cloning
(BC) [27, 28], Linear Parameter-Varying Dynamical System (LPV-DS) [7], and
Stable Dynamical System Learning Using Euclideanizing Flows (SDS-EF) [13].
Among these methods, only BC lacks formal global asymptotic stability
guarantees. Nevertheless, comparing SNDS against BC enables us to gauge SNDS’s
accuracy and global stability against an unrestricted neural policy.

To assess the effectiveness of SNDS in comparison with the named baselines, we
carry out the learning process for each method on designated handwriting
dataset motions. For each motion, we randomly split the demonstrated dataset
in the dataset into train ($0.8$) and test ($0.2$) sets. The policy learning
stage is carried out on the training data, and the optimal policy is used to
generate policy rollouts, $\mathbf{{x}}^{g}$, using the Euler’s method.
Policies are subsequently evaluated by calculating Mean-Squared Error (MSE),

$MSE(\dot{\mathbf{{x}}}^{g},\dot{\mathbf{{x}}}^{d})=\frac{\sum_{d=1}^{N_{d}^{test}}\sum_{s=1}^{N_{s}}(\phi(\mathbf{{x}}^{d}[s])-\dot{\mathbf{{x}}}^{d}[s])^{2}}{2N_{d}^{test}N_{s}},$
(7)

and Dynamic Time Warping (DTW),

$DTW(\mathbf{{x}}^{g},\mathbf{{x}}^{d})=\min_{p\in\mathcal{P}({\mathbf{{x}}}^{g},{\mathbf{{x}}}^{d})}\left(\sum_{(i,j)\in\pi}\mathit{Euc}(\mathbf{{x}}^{g}[i],\mathbf{{x}}^{d}[j])^{q}\right)^{\frac{1}{q}},$
(8)

for handwriting dataset and experiments in higher state space dimensions to
better capture discrepancies between policy rollouts and expert
demonstrations. Within Equation 8, $p$ represents an alignment path,
$\mathcal{P}$ is the set of all admissible paths, $\mathit{Euc}$ is Euclidean
distance [29], and $q=2$. Notice that to calculate DTW, we need to generate an
entire trajectory, while with MSE, only the generated actions are compared
against true labels. We conduct the training and assessment procedures
repeatedly with 10 random seeds, and present the mean and standard deviation
of the results.

Figure 5: Policy rollouts for SNDS and other baselines. Notice the highlighted
inaccuracies or instabilities (yellow) for other methods. The acquired
policies are optimized using the N-shaped data of the handwriting dataset.

###  IV-B Handwriting dataset policies

We apply SNDS and selected baselines to each set of demonstrations as
explained in Section IV-A, and compare the reproduction accuracy and
computation time 222Note that we restrict our computational resources to a
single Core-i7 Gen8 CPU, despite that SNDS and SDS-EF can leverage parallel
computing. across various motions in the handwriting dataset. Figure 4
presents a numerical comparison of both reproduction accuracy and computation
time between SNDS and the selected baselines. The results indicate an
acceptable level of reproduction accuracy while providing GAS certificates. In
rare cases, SNDS shows a higher error compared to baselines, but outperforms
others for the majority of elected motions, and maintains a reasonable
computation time.

Figure 6: Policy rollouts when training on a single Sine-shaped demonstration
of the handwriting data. Inaccuracies or instabilities are highlighted. SNDS
performs better owing to the expressive architecture and customized loss.

To closely examine the GAS property of each method, we plot the acquired
policies using streamlines and generated two simulated rollouts with Euler’s
method in Figure 5. It becomes apparent that SNDS, and LPV-DS provide
predictability through guaranteeing GAS, while other approaches, such as BC
and SDS-EF fail to render similar certificates. On a general note, even though
we observe unstable policies trained with SDS-EF, instabilities occur only for
a subset of motions in the handwriting dataset, and tend to intensify in
regions further away from expert’s demonstrations.

To further showcase the sample inefficiency associated with Gaussian mixture
model and diffeomorphism-based methods, such as LPV-DS and SDS-EF,
respectively, we repeat the experiments with the training set reduced to only
a single demonstration. The learned policies in this scenario are illustrated
in Figure 6, indicating the sample efficiency of SNDS when compared to LPV-DS
and SDS-EF in the presence of limited demonstration samples.

###  IV-C SE(3) policy training

To illustrate the applicability of our method in high-dimensional state
spaces, such as SE(3), we consider both the position and orientation of the
robot and redefine the state variable, $\mathbf{{x}}$, accordingly. We pick
the _snake motion_ —an overly complex and nonlinear motion presented in [7]
and 2 different motions in the handwriting dataset, namely _Sine-_ and _G-
shaped_ motions, to consolidate our higher dimensional performances. We use
synthetic demonstrations for orientation, as the demonstrations only comprise
translation and linear velocity. After training SNDS on this data, the policy
generates trajectories in SE(3) that guide the end-effector toward the desired
pose.

In Figure 7, we compare the performance of our method against the baselines
trained with the Messy Snake data in SE(3). We utilize the DTW metric to
evaluate the policy rollouts against the original trajectories. Next, we
deploy the policy trained with SNDS on the simulated arm in PyBullet. Figure 8
offers an overview of policy deployment in simulation trained on the Messy
Snake data.

Figure 7: DTW between policy rollouts and expert’s data for the complex Messy
Snake demonstrations in SE(3) (left). SNDS can better extend to higher state
space dimensions compared to former baselines. The trained SNDS policy over
the expert demonstrations is depicted on the right.

Finally, we deploy the learned policy directly on a similar real-world
manipulator platform. Recordings of experiments conducted with the Kinova
Jaco2 arm are part of the supplementary video. Upon viewing the footage, it is
evident that the robot reproduces the same motion as the simulation. The
smooth sim-to-real transfer is made possible by the presence of Lyapunov
stability conditions. The learned policy is robust to model error, since the
Lyapunov condition ensures that the robot consistently progresses toward the
target.

##  V Discussion

The experiments demonstrate that employing SNDS to model a policy results in
improved predictability and safety when recovering from states far beyond
those covered in the original demonstrations. SNDS proves to be more effective
by utilizing expressive and scalable neural architecture in comparison to the
state of the art. Consequently, SNDS can acquire highly nonlinear stable
policies in larger state spaces with reduced computational overhead. Our
approach also reduces the restrictive assumptions regarding the class of
Lyapunov functions to convexity, while the adopted projection technique
ensures a smooth training process on a variety of expert data, as evident by
the experiments.

While SNDS offers global convergence and predictability, there remains a
concern regarding physically impossible trajectories for robots to replicate,
particularly when considering manipulators and their joint or torque limits.
Future enhancements could address this by incorporating safe regions using
control barrier functions. Moreover, it is worth noting that modeling the
Lyapunov candidate as a strictly convex function is not an essential
requirement. Future research could explore invertible transformations to ease
this restriction. Additionally, integrating obstacle avoidance and joint
constraints into our formulation are promising paths to explore beyond the
current scope. More ambitious extensions might delve into the implications of
utilizing SNDS policies in reinforcement learning or training and deploying
stable policies on legged robots, especially for gait control.

##  VI Conclusion

We outlined a training process to learn expressive neural policies through
imitating expert demonstrations. We effectively enforce global stability
throughout the training process by adhering to the conditions of Lyapunov
stability theory. This ensures that the resulting policy reliably converges to
a predetermined target, regardless of initial conditions, velocity and time
variations, or unexpected perturbations encountered during deployment. Our
theoretical findings are accompanied by simulation and real-world benchmarking
against best-performing methods in the field.

Figure 8: Policy deployment in simulation (top) and real-world Kinova Jaco2
arm (bottom) scenarios. The policy is trained on a single Messy Snake
demonstration. The simulation rollouts (red) are close to the expert’s
demonstration (blue), and are similar to the rollouts in Figure 7.

## Appendix

Network architecture and hyperparameters. Our approach employs a 4-hidden
layer feed-forward neural network for policy, with Leaky ReLU activations and
layer sizes of 2-256-256-128-128-2. For the Lyapunov function, we utilize an
ICNN architecture [22] with 3 hidden layers of 2-128-128-128-1 nodes along
with Leaky ReLU and Soft plus activations as required by the original design.
We set $\gamma_{i}=\frac{1}{2^{i}}$ and $N_{w}=2$ to limit the cost function’s
horizon. Other parameters, e.g., learning rate and batch size, are selected
and explained in our code repository. We use grid-search to pick the
hyperparameters if computationally feasible.

Datasets specifications. The Handwriting Dataset [24] contains 30 sets of 2D
handwriting motions recorded from a Tablet-PC via user input. There are seven
demonstrations per motion, starting from slightly different initial positions
but ending at (0, 0). Each demonstration encompasses position (2 × 1000) as
learning features and velocity (2 × 1000) as labels. We pick the following
representative set of motion for our experiments: Sine, P, Worm, G, N, DBLine,
C, Angle. Note that the Messy Snake motion is not a part of this dataset. The
Messy Snake data, gathered in [7], comprises three demonstrations with various
sample sizes per demonstration, totaling to (2 × 1592), and we sample a single
demonstration for policy learning and deployment.

Computational resources. We utilize a single machine with NVIDIA GeForce RTX
4050 GPU, Intel Core i7-13620H CPU, and 32 GB DDR4 RAM.

Codebase and reproduction. Consult README.md in our GitHub repository
github.com/aminabyaneh/stable-imitation-policy to access SNDS experiments.

## Acknowledgment

We express our heartfelt gratitude to Anya Forestell for her invaluable
hardware support during our experiments.

## References

  * [1] S. Schaal, “Is imitation learning the route to humanoid robots?” _Trends in cognitive sciences_ , vol. 3, no. 6, pp. 233–242, 1999.
  * [2] A. Hussein, M. M. Gaber, E. Elyan, and C. Jayne, “Imitation learning: A survey of learning methods,” _ACM Computing Surveys (CSUR)_ , vol. 50, no. 2, pp. 1–35, 2017.
  * [3] N. Figueroa and A. Billard, “Locally active globally stable dynamical systems: Theory, learning, and experiments,” _The International Journal of Robotics Research_ , vol. 41, no. 3, pp. 312–347, 2022.
  * [4] M. Hersch, F. Guenter, S. Calinon, and A. Billard, “Dynamical system modulation for robot learning via kinesthetic demonstrations,” _IEEE Transactions on Robotics_ , vol. 24, no. 6, pp. 1463–1467, 2008.
  * [5] S. M. Khansari-Zadeh and A. Billard, “Learning stable nonlinear dynamical systems with gaussian mixture models,” _IEEE Transactions on Robotics_ , vol. 27, no. 5, pp. 943–957, 2011.
  * [6] S. M. Khansari-Zadeh and A. Billard, “Learning control Lyapunov function to ensure stability of dynamical system-based robot reaching motions,” _Robotics and Autonomous Systems_ , vol. 62, no. 6, pp. 752–765, 2014.
  * [7] N. Figueroa and A. Billard, “A physically-consistent bayesian non-parametric mixture model for dynamical system learning,” in _2nd Annual Conference on Robot Learning_ , 2018, pp. 927–946.
  * [8] K. Neumann and J. J. Steil, “Learning robot motions with stable dynamical systems under diffeomorphic transformations,” _Robotics and Autonomous Systems_ , vol. 70, pp. 1–15, 2015.
  * [9] H. Ravichandar, I. Salehi, and A. Dani, “Learning partially contracting dynamical systems from demonstrations,” in _Conference on Robot Learning_. PMLR, 2017, pp. 369–378.
  * [10] F. Khadivar, I. Lauzana, and A. Billard, “Learning dynamical systems with bifurcations,” _Robotics and Autonomous Systems_ , vol. 136, p. 103700, 2021\.
  * [11] A. Abyaneh and H.-C. Lin, “Learning Lyapunov-stable polynomial dynamical systems through imitation,” in _7th Annual Conference on Robot Learning_ , 2023.
  * [12] Y.-C. Chang, N. Roohi, and S. Gao, “Neural Lyapunov control,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * [13] M. A. Rana, A. Li, D. Fox, B. Boots, F. Ramos, and N. Ratliff, “Euclideanizing flows: Diffeomorphic reduction for learning stable dynamical systems,” in _Learning for Dynamics and Control_. PMLR, 2020, pp. 630–639.
  * [14] J. Zhang, H. B. Mohammadi, and L. Rozo, “Learning riemannian stable dynamical systems via diffeomorphisms,” in _6th Annual Conference on Robot Learning_ , 2022.
  * [15] J. Ho and S. Ermon, “Generative adversarial imitation learning,” _Advances in neural information processing systems_ , vol. 29, p. 4572–4580, 2016.
  * [16] J. Fu, K. Luo, and S. Levine, “Learning robust rewards with adversarial inverse reinforcement learning,” in _International Conference on Learning Representations_ , 2018.
  * [17] H. Dai, B. Landry, L. Yang, M. Pavone, and R. Tedrake, “Lyapunov-stable neural-network control,” in _Robotics: Science and Systems_ , 2021.
  * [18] A. Coulombe and H.-C. Lin, “Generating stable and collision-free policies through Lyapunov function learning,” in _International Conference on Robotics and Automation_ , 2023, pp. 3037–3043.
  * [19] P. Abbeel and A. Y. Ng, “Apprenticeship learning via inverse reinforcement learning,” in _International conference on Machine learning_ , 2004, p. 1.
  * [20] B. D. Ziebart, A. L. Maas, J. A. Bagnell, A. K. Dey, _et al._ , “Maximum entropy inverse reinforcement learning.” in _Association for the Advancement of Artificial Intelligence_ , vol. 8, 2008, pp. 1433–1438.
  * [21] R. L. Devaney, _An introduction to chaotic dynamical systems_. CRC press, 2021.
  * [22] B. Amos, L. Xu, and J. Z. Kolter, “Input convex neural networks,” in _International Conference on Machine Learning_. PMLR, 2017, pp. 146–155.
  * [23] J. Z. Kolter and G. Manek, “Learning stable deep dynamics models,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * [24] S. M. Khansari-Zadeh and A. Billard, “LASA Handwriting Dataset,” https://cs.stanford.edu/people/khansari/download.html#SEDS_reference, 2011\.
  * [25] M. Bruveris, “Optimal reparametrizations in the square root velocity framework,” _SIAM Journal on Mathematical Analysis_ , vol. 48, no. 6, pp. 4335–4354, 2016.
  * [26] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _International Conference on Learning Representations_ , 2015.
  * [27] D. A. Pomerleau, “Alvinn: An autonomous land vehicle in a neural network,” in _Advances in neural information processing systems_ , vol. 1, 1988, pp. 305–313.
  * [28] F. Torabi, G. Warnell, and P. Stone, “Behavioral cloning from observation,” in _Proceedings of the 27th International Joint Conference on Artificial Intelligence_ , 2018, pp. 4950–4957.
  * [29] E. Keogh and C. A. Ratanamahatana, “Exact indexing of dynamic time warping,” _Knowledge and information systems_ , vol. 7, pp. 358–386, 2005.


# Solving the Funarg Problem with Static Types

Caleb Helbling<EMAIL_ADDRESS>0000-0002-6118-3061 Purdue University610
Purdue MallWest LafayetteIndianaUSA47907 and  Fırat Aksoy
<EMAIL_ADDRESS>0000-0003-3109-9544 Istanbul
UniversityIstanbulTurkey

(2021)

###### Abstract.

The difficulty associated with storing closures in a stack-based environment
is known as the funarg problem. The funarg problem was first identified with
the development of Lisp in the 1970s and hasn’t received much attention since
then. The modern solution taken by most languages is to allocate closures on
the heap, or to apply static analysis to determine when closures can be stack
allocated. This is not a problem for most computing systems as there is an
abundance of memory. However, embedded systems often have limited memory
resources where heap allocation may cause memory fragmentation. We present a
simple extension to the prenex fragment of System F that allows closures to be
stack-allocated. We demonstrate a concrete implementation of this system in
the Juniper functional reactive programming language, which is designed to run
on extremely resource limited Arduino devices. We also discuss other solutions
present in other programming languages that solve the funarg problem but
haven’t been formally discussed in the literature.

funarg, closure, stack, heap, memory, functional, embedded

††copyright: acmcopyright††journalyear: 2021††doi:
10.1145/1122445.1122456††conference: 33rd Symposium on Implementation and
Application of Functional Languages; September 01–03, 2021; Online††booktitle:
IFL ’21: ACM Symposium on Implementation and Application of Functional
Languages, September 01–03, 2021, Online††price: 15.00††isbn:
978-1-4503-XXXX-X/18/06††ccs: Software and its engineering Procedures,
functions and subroutines††ccs: Software and its engineering Functional
languages

##  1\. Introduction

The funarg problem refers to the difficulty associated with compilation of
first-class functions in programming languages for stack-based environments.
The problem arises in the bodies of nested functions where the nested function
refers to identifiers defined in the parent function’s lexical scope which are
not available in it’s own lexical scope. The standard solution is to allocate
closures on the heap. There are two variants of this problem: upwards and
downwards. The upwards funarg problem refers to the management of the stack
when a function returns another function (that may potentially have a non-
empty closure), whereas the downwards funarg problem refers to the management
of the stack when a function is passed to another function (that may also have
a non-empty closure).

Solving the upwards funarg problem is considered to be much more difficult
than the downwards funarg problem. The downwards funarg problem is easily
solved by copying the closure’s lexical scope (captured variables) to the top
of the program’s stack. For the upwards funarg problem, an issue arises with
deallocation. A returning function may capture local variables, making the
size of the closure dynamic. Due to this behavior, the closure itself cannot
be copied down the stack without causing stack corruption. We must keep in
mind that, as soon as the function returns, any local variables will be lost
when they are popped off the stack.

In 1968, Weizenbaum (Weizenbaum, 1968) made an analysis of function closures
using lambda calculus, and showed that function closures in lambda calculus
cannot have a stack based data structure and must in fact consist of a tree.
Moses (Moses, 1970) further analyzed the problem for the Algol programming
language. Early attempts to solve the funarg problem (Sandewall, 1971) appear
to be hampered by a number of factors, including dynamic scoping, lack of
static types and mutable variables.

Despite the titles of these early papers, we contend that the funarg problem
(in particular, the upwards funarg problem) has not been adequately solved. In
our approach, we solve the upwards funarg problem by statically determining a
closure’s size at compile time. This allows the closure from the upwards
funarg problem to be copied into the calling function, as if it were any other
statically sized piece of data. Since the type of a closure’s lexical scope
can be viewed as a structural record type, we can make our language
polymorphic over closures and their lexical scopes. This allows higher order
functions such as map (downwards funarg problem) and compose (upwards and
downwards funarg problem) to be written very concisely and naturally.

Analyzing closures from a type theoretic perspective is not a new exercise
(Minamide et al., 1996; Morrisett and Harper, 1998; Morrisett et al., 1999;
Bowman and Ahmed, 2018). These papers prefer to treat closures as explicit
environments, whose type is determined by an existential. This is not helpful
to solving the funarg problem, since existentials are typically treated as
boxed values and must be allocated on the heap. Therefore our contribution is
unique in the sense that we seek to solve the funarg problem rather than
analyzing closures from a purely type theoretic environmental perspective.

###  1.1. Contributions

We describe the following contributions:

  * •

We demonstrate how the prenex fragment of System F can be extended to enable
stack allocated closures.

  * •

We discuss a concrete implementation of this approach in the Juniper
programming language.

  * •

We analyze the limitations of the extension and identify the connection
between types and closures.

  * •

We make an analysis of other approaches found in existing programming
languages, but haven’t been discussed formally in the literature.

##  2\. Extending System F

In this section, we restrict and extend System F to enable stack allocated
closures. We restrict System F to it’s prenex fragment which separates regular
System F types into type schemes and types. This is a necessary restriction to
prevent the emergence of existential types in the type system which prevents
determining the sizes of types at compile-time. The primary addition is that
of the closure types, which are defined as function types with their lexical
scopes attached to them. The lexical scopes $\delta$ are either type variables
or records of captured variables and their respective types. These lexical
scope types are are considered equivalent if both the field names and the
types of those fields are identical. These additions make it possible for
closures to be polymorphic and enable stack allocation of the closure since
the size of the closure is known statically at compile-time. Table 1 defines
the grammar of our language.

Expressions
---
$e:=x\mid\lambda x:A.e\mid e_{1}\,e_{2}\mid\Lambda\alpha.e\mid e\,A$
Type Schemes
$\sigma:=A\mid\forall{a}.\sigma$
Types
$A,B,C:=\alpha\mid A-C\rightarrow B\mid\delta$
Lexical Scope Types
$\delta:=\\{x_{0}:A_{0},...,x_{n}:A_{n}\\}$
Values
$v:=\left(\lambda x.e\right)\left[\Delta\right]\mid\Lambda.e$
Contexts
$\Gamma::=\emptyset\mid\Gamma,\alpha\mid\Gamma,x:A$
Environments
$\Delta:=\\{x_{0}\mapsto v_{0},...,x_{n}\mapsto v_{n}\\}$
Table 1. Grammar

###  2.1. Big-step semantics

Stack allocation of closures require the lambda values and their lexical
scopes to be propagated and used in pairs. For this reason, we combine the
lambda values with their environments into a single value representation. When
a lambda is applied, a flat representation of the lambda’s lexical scope is
constructed and used as the environment along with the lambda’s argument for
further evaluation. Figure 1 presents the big-step semantics for our language.

$\frac{v=\Delta(x)}{\Delta\vdash x\Downarrow v}\hskip 5.0pt[\text{\small EVAL-
VAR}]$

$\frac{}{\Delta\vdash\lambda x:A.e\Downarrow(\lambda
x.e)[\\{y\mapsto\Delta(y)\mid y\in\textbf{fv}(\lambda x:A.e)\\}]}$

$[\text{\small EVAL-ABS}]$

$\frac{\Delta_{1}\vdash e_{1}\Downarrow(\lambda
x.e_{3})[\Delta_{2}]\quad\;\Delta_{1}\vdash e_{2}\Downarrow
v^{\prime}\quad\;\Delta_{2},x\mapsto v^{\prime}\vdash e_{3}\Downarrow
v}{\Delta_{1}\vdash e_{1}\,e_{2}\Downarrow v}$

$[\text{\small EVAL-APP}]$

$\frac{}{\Delta\vdash\Lambda\alpha.e\Downarrow\Lambda.e}\hskip
5.0pt[\text{\small EVAL-TABS}]$

$\frac{\Delta\vdash e_{1}\Downarrow\Lambda.e_{2}\qquad\Delta\vdash
e_{2}\Downarrow v}{\Delta\vdash e_{1}\,A\Downarrow v}\hskip 5.0pt[\text{\small
EVAL-TAPP}]$

Figure 1. Big step semantics

###  2.2. Typing

Figure 2 presents the typing rules for our language. There are several changes
made to the type rules for the prenex fragment of System F. Of primary
interest is the T-ABS rule which also constructs the type of the lexical
scope. The construction of the lexical scope itself is performed by T-DELTA,
which finds all free variables within the lambda, looks up their type in
$\Gamma$ and returns the corresponding lexical scope type. The rule T-TABS is
also changed to restrict type abstractions to contain no free variables
themselves which simplifies the type system and assists in compilation to a
stack based environment.

$\frac{x:A\in\Gamma}{\Gamma\vdash x:A}\hskip 5.0pt[\text{\small T-VAR}]$

$\frac{\Gamma,x:A\vdash
e:B\qquad\Gamma\vdash\lambda{x:A}.e\rightsquigarrow\delta}{\Gamma\vdash\lambda{x:A}.e:A-\delta\rightarrow
B}\hskip 5.0pt[\text{\small T-ABS}]$

$\frac{\Gamma\vdash e_{1}:A-\delta\rightarrow B\qquad\Gamma\vdash
e_{2}:A}{\Gamma\vdash e_{1}\,e_{2}:B}\hskip 5.0pt[\text{\small T-APP}]$

$\frac{\Gamma,\alpha\vdash
e:\sigma\qquad\textbf{fv}(e)=\emptyset}{\Gamma\vdash\Lambda{\alpha}.e:\forall{\alpha}.\sigma}\hskip
5.0pt[\text{\small T-TABS}]$

$\frac{\Gamma\vdash e:\forall{\alpha}.\sigma}{\Gamma\vdash
e\,A:\sigma[\alpha:=A]}\hskip 5.0pt[\text{\small T-TAPP}]$

$\frac{y_{i}\in\mathbf{fv}(\lambda x:A.e)\qquad\Gamma\vdash
y_{0}:C_{0},...,\Gamma\vdash y_{n}:C_{n}}{\Gamma\vdash\lambda
x:A.e\rightsquigarrow\\{y_{0}:C_{0},...,y_{n}:C_{n}\\}}\hskip
5.0pt[\text{\small T-DELTA}]$

Figure 2. Typing rules.

###  2.3. Examples

In this section, we will look at two examples that demonstrate that our
solution for the funarg problem works for accurately typing two of the most
commonly used functions in functional programming: map and compose. For the
map example, we will assume that the type system is enriched with a list type
constructor. The type of these functions are considerably more noisy than the
standard System F types. Fortunately, most of the noise can be elided through
type inference.

$map:\forall\alpha.\forall\beta.\forall\delta.(\alpha-\delta\rightarrow\beta)-\\{\\}\rightarrow\alpha\,\mathbf{list}-\\{f:\alpha-\delta\rightarrow\beta\\}\rightarrow\beta\,\mathbf{list}$

For the compose example we give the definition of compose as well as its type
in our System F extension.

$compose:=\Lambda\alpha.\Lambda\beta.\Lambda\gamma.\Lambda\delta_{1}.\Lambda\delta_{2}.\\\
\lambda f:(\beta-\delta_{1}\rightarrow\gamma).\lambda
g:(\alpha-\delta_{2}\rightarrow\beta).\lambda x:\alpha.f\,(g\,x)$

$compose:\forall\alpha.\forall\beta.\forall\gamma.\forall\delta_{1}.\forall\delta_{2}.\\\
(\beta-\delta_{1}\rightarrow\gamma)-\\{\\}\rightarrow(\alpha-\delta_{2}\rightarrow\beta)-\\{f:\beta-\delta_{1}\rightarrow\gamma\\}\rightarrow\alpha\\\
-\\{f:\beta-\delta_{1}\rightarrow\gamma,g:\alpha-\delta_{2}\rightarrow\beta\\}\rightarrow\gamma$

##  3\. Limitations and Extensions

The most obvious limitation of the extension presented in this paper is that
all branches in a program that return a closure must have the exact same
lexical scope. One would think that this would be a problem since in most
languages, branches of conditionals may return lambdas with different lexical
scopes, as long as the types are equivalent. A close analysis of the cond
function in System F when combined with Church-encoded Boolean values true and
false reveals a surprising result, this problem simply does not appear in the
extension described in this paper.

To further elucidate our point, we define cond operator as a function which
takes in a true branch and a false branch and a Church-encoded Boolean value
(as described in equations 3, 4, 5, and 6) and returns the respective branch
in relation to it’s Boolean input.

(1) $\text{cond}:=\Lambda\alpha.\Lambda\beta.\Lambda\gamma.\Lambda\delta\\\
\lambda t:\alpha.\lambda f:\beta.\lambda
c:(\alpha-\\{\\}\rightarrow\beta-\delta\rightarrow\gamma).c\,t\,f$

(2) $\text{cond}:\forall\alpha.\forall\beta.\forall\gamma.\forall\delta.\\\
\alpha-\\{\\}\rightarrow\beta-\\{t:\alpha\\}\rightarrow(\alpha-\\{\\}\rightarrow\beta-\delta\rightarrow\gamma)-\\{t:\alpha,f:\beta\\}\rightarrow\gamma$

(3) $\text{true}:=\Lambda\alpha.\Lambda\beta.\lambda t:\alpha.\lambda
f:\beta.t$

(4)
$\text{true}:\forall\alpha.\forall\beta.\alpha-\\{\\}\rightarrow\beta-\\{t:\alpha\\}\rightarrow\alpha$

(5) $\text{false}:=\Lambda\alpha.\Lambda\beta.\lambda t:\alpha.\lambda
f:\beta.f$

(6)
$\text{false}:\forall\alpha.\forall\beta.\alpha-\\{\\}\rightarrow\beta-\\{\\}\rightarrow\beta$

Applying the cond operator with true and false in equations 7 and 8 reveals
why Church-encoding does not suffer from the branching problem. In System F,
the branch chosen by cond is already present in the type system by selecting
the type for $\gamma$. Fundamentally, this occurs because true and false have
different type signatures, that also determine the branch taken.

(7)
$((\text{cond}\,\alpha\,\beta\,\alpha\,\\{t:a\\})\,t\,f\,(\text{true}\,\alpha\,\beta\,)):\alpha$

(8)
$((\text{cond}\,\alpha\,\beta\,\beta\,\\{\\})\,t\,f\,(\text{false}\,\alpha\,\beta\,)):\beta$

However, in a language with sum types and pattern matching (as Boolean types
are in most languages), the branching problem with closures will be present.
Whether or not this is a major issue is debatable, however, we can say for
certain that this removes some generality from our approach. An interesting
approach to recovering this generality is a paper on open closures by Scherer
and Hoffmann (Scherer and Hoffmann, 2013), however their approach only applies
to simply typed lambda calculus and not for System F.

The second limitation appears under the presence of mutually recursive
functions. When compiling top level functions, we can consider their closure
to be empty since we can refer to them using a simple function pointer. For
inner functions that are often mutually recursive, the situation is more
complicated. The lexical scope of all functions declared in a let-rec block
can be constructed by taking the union of all of the lexical scopes of all the
functions declared in the block. This constructs a common lexical scope that
can be passed to all the related functions. Therefore, all functions declared
in the let-rec block will have identical lexical scope types.

When one function in the let-rec block calls another function in the same let-
rec block, it simply passes the lexical scope that it was passed. Within the
lexical scope type itself, we place all the functions within the let-rec block
with the non-closure function type $A\rightarrow B$. This indicates that this
function should be passed the lexical scope that this function is contained
within.

Extending the system presented in this paper to the impredicative System F at
first glance seems feasible. We believe that the funarg problem is orthogonal
issue to that of first-class polymorphism and stack allocated existential
types. However we have yet to make a full analysis of how the lexical
environments presented here interact with System F. One possible issue that
may crop up is subtle interactions with existential types, which can be
encoded in System F (Girard et al., 1989).

##  4\. Implementation

We have implemented the closure system given here in Juniper (Helbling and
Guyer, 2016), a functional reactive programming language for the Arduino. The
Juniper language transpiles all of its code to C++. The analysis of C++ in
section 5 reveals the issues facing the use of lambdas as a transpile target.
To work around the issues presented in that section, we implement our own
function class. This class is constructed by passing in a lexical scope, along
with a function pointer. In C++, lambda types that do not capture anything can
be casted to type safe function pointers. We provide our own stack allocated
closure support by passing in the lexical scope as the first argument to the
underlying lambda function. Within this lambda, the compiler inserts variable
assignments which extracts the contents of the lexical scope into local
variables. Lambdas that do not capture anything are not passed a lexical
scope.

⬇

1template <typename ClosureType, typename Result, typename …Args>

2class function;

3

4template <typename Result, typename …Args>

5class function<void, Result(Args…)> {

6private:

7 Result(*F)(Args…);

8

9public:

10 function(Result(*f)(Args…)) : F(f) {}

11

12 Result operator()(Args… args) {

13 return F(args…);

14 }

15};

16

17template <typename ClosureType, typename Result, typename …Args>

18class function<ClosureType, Result(Args…)> {

19private:

20 ClosureType Closure;

21 Result(*F)(ClosureType&, Args…);

22

23public:

24 function(ClosureType closure, Result(*f)(ClosureType&, Args…)) :
Closure(closure), F(f) {}

25

26 Result operator()(Args… args) {

27 return F(Closure, args…);

28 }

29};

A lexical scope struct is defined for every combination of environment names
discovered in the program. The types of the variables are made generic, so
that these lexical scope structs can be declared before any of the other types
in the program. Below is an example of the lexical scope struct generated by
the Juniper compiler for the compose function. These structs are created and
passed to the function constructor the moment the lambda declared. This
ensures that the closure contains a snapshot of the current environment.
Importantly, if a mutation of any local variables occurs after a closure has
been constructed, the closure will not be updated (since it essentially takes
a snapshot of the local environment the moment it is created). For purely
functional languages this is not an issue since variables cannot be mutated.

⬇

1// Compose lexical scope

2template<typename T1,typename T2>

3struct closuret_0 {

4 T1 f;

5 T2 g;

6

7 closuret_0(T1 init_f, T2 init_g) :

8 f(init_f), g(init_g) {}

9};

##  5\. Existing Approaches

In this section we will analyze a number of different existing languages and
how they solve the funarg problem. To our knowledge the system presented in
this paper is the only language that will allow the type of a closure to be
fully written out with text (ie, it does not generate an implicit type or
class).

In C++, each declaration of a lambda implicitly declares a new class (ISO,
2017) and therefore type in the system. Since these classes are auto-
generated, it is impossible to write their types using text, and programmers
must rely on the auto keyword. These C++ lambdas are indeed stack allocated
since their underlying representation is a class with a operator() defined.
The problems begin however as soon as you start passing around or returning
lambdas. Higher order functions that consume lambdas must be made generic over
all of these implicit lambda classes. This makes it is impossible to add any
sort of constraints on the types of the arguments of the lambda, or the return
type of the lambda. Returning a lambda also becomes very annoying, and the
auto keyword must be used as the return type. Once again, we see that it is
impossible to constrain the return type to anything in particular. The
solution used by most C++ programmers is to make use of the standard library
std::function, which turns the stack allocated closure into a heap allocated
entity. Perhaps in the future, clever use of C++ concepts could add extra type
safety to the use of these implicit lambda classes.

To our knowledge, the only other language that has similar capabilities to
what is presented in this paper is Rust (Matsakis and Klock II, 2014). Just as
in C++, the Rust system operates by generating a unique type for every
closure. However Rust goes a step further with its trait system. Like in C++,
accepting a closure as a parameter to another function requires the use of a
single generic parameter. Unlike C++, the type of this parameter can be
constrained with traits. Additionally, returning a closure requires the
textual use of one of these traits.

##  6\. Future Work

In the future we would like to further extend our type system for enriched
languages. We are interested in how adding conditionals and sum types to the
language will force changes in the handling of the lexical scopes. We also
would like to give a more thorough description of how to handle the let-rec
case.

We would also like to extend this system for the impredicative System F
without restrictions as much as possible, which we predict to potentially
force additional changes in the representation of the lexical scopes.

##  7\. Conclusion

In this paper we have described an extension to the prenex fragment of System
F that solves the funarg problem. Our solution solves both the downwards and
the upwards funarg problem and allows closures to be completely stack
allocated. Functional languages can take advantage of this new system to
optimize memory usage, and this extension is particularly useful on embedded
systems. We also discuss extensions to our language and identify the
possibility to extend this work to solve the conditional problem as well as
stack allocation for first-class polymorphism. We discussed how existing
languages solve the funarg problem, and conclude that our approach is the
first so far that allows writing closure type signatures without auto-
generated types or classes.

## References

  * (1)
  * Bowman and Ahmed (2018) William J Bowman and Amal Ahmed. 2018.  Typed closure conversion for the Calculus of Constructions.  (Aug. 2018).  arXiv:1808.04006 [cs.PL]
  * Girard et al. (1989) Jean-Yves Girard, Paul Taylor, and Yves Lafont. 1989\.  _Proofs and types_. Vol. 7.  Cambridge university press Cambridge.
  * Helbling and Guyer (2016) Caleb Helbling and Samuel Z Guyer. 2016.  Juniper: a functional reactive programming language for the Arduino. In _Proceedings of the 4th International Workshop on Functional Art, Music, Modelling, and Design_. 8–16.
  * ISO (2017) ISO. 2017.  _ISO/IEC 14882:2017 Information technology — Programming languages — C++_ (fifth ed.).  1605 pages.  https://www.iso.org/standard/68564.html
  * Matsakis and Klock II (2014) Nicholas D Matsakis and Felix S Klock II. 2014.  The rust language. In _ACM SIGAda Ada Letters_ , Vol. 34. ACM, 103–104.
  * Minamide et al. (1996) Yasuhiko Minamide, Greg Morrisett, and Robert Harper. 1996\.  Typed closure conversion. In _Proceedings of the 23rd ACM SIGPLAN-SIGACT symposium on Principles of programming languages - POPL ’96_ (St. Petersburg Beach, Florida, United States). ACM Press, New York, New York, USA.
  * Morrisett and Harper (1998) Greg Morrisett and Robert Harper. 1998.  Typed closure conversion for recursively-defined functions.  _Electron. Notes Theor. Comput. Sci._ 10 (1998), 230–241.
  * Morrisett et al. (1999) Greg Morrisett, David Walker, Karl Crary, and Neal Glew. 1999\.  From system F to typed assembly language.  _ACM Trans. Program. Lang. Syst._ 21, 3 (May 1999), 527–568.
  * Moses (1970) Joel Moses. 1970\.  The function of FUNCTION in LISP or why the FUNARG problem should be called the environment problem.  _SIGSAM Bull._ 15 (July 1970), 13–27.
  * Sandewall (1971) Erik Sandewall. 1971\.  A proposed solution to the FUNARG problem.  _SIGSAM Bull._ 17 (Jan. 1971), 29–42.
  * Scherer and Hoffmann (2013) Gabriel Scherer and Jan Hoffmann. 2013.  Tracking Data-Flow with Open Closure Types. In _Logic for Programming, Artificial Intelligence, and Reasoning_. Springer Berlin Heidelberg, 710–726.
  * Weizenbaum (1968) Joseph Weizenbaum. 1968\.  The funarg problem explained.  _unpublished memorandum, MIT_ (1968).