\section{Introduction}
\indent State estimation has found significant application in all areas of science and engineering to accelerate technological advancements. The particle filtering (PF) approach has arisen as a crucial means, especially for nonlinear non-Gaussian systems. At the core, PF updates an ensemble of particles {\em via} sequential importance sampling such that their empirical distribution approximates the target distribution of the system's state. However, the PF often suffers from particle degeneracy, whereby a majority of the particle's weights diminish to almost zero, leading to a poor approximation of the target distribution. Even though one can use a large number of particles so that at least some of them have significant weights, the number of particles needed often grows catastrophically with the system's dimension.

A survey of the literature indicates an ongoing search for principled ways to address the above issue. A straightforward, popular method is resampling, which randomly replaces low-weight particles with high-weight ones~\cite{Sarkka:Cambridge:2013}. More sophisticated mechanisms have also been developed, e.g.,~\cite{Chorin:PNAS:2009,Yang:TAC:2013,Raihan:AUTO:2018,Stano:AUTO:2013}. Among them, the method of implicit sampling has inspired much attention. Implicit sampling exploits the idea that fewer particles will be needed as long as they lie in the high-probability regions of the target distribution~\cite{Chorin:PNAS:2009}. By design, it constructs probability distributions assumed for the particles and utilizes them as a reference to select particles from the target distribution's high-probability regions. The resultant implicit PF (IPF) shows effectiveness in keeping the number of particles manageable even for high-dimensional systems~\cite{Chorin:CAMCS:2010,Chorin:Springer:2013}. However, implementing IPF is non-trivial as it entails solving a nonlinear optimization problem and then a nonlinear algebraic equation. While different numerical approaches are used  in~\cite{Chorin:CAMCS:2010,Chorin:Springer:2013,Su:SIAM:2017,Weir:NPG:2013}, the computation is expensive and tedious. This motivates the development of efficient alternative approaches to approximate the IPF.

This paper proposes that the IPF framework can be approximately realized by a bank of nonlinear Kalman filters (KFs). This realization shows more computational tractability while maintaining the IPF's merit to sample from regions of high probability to achieve high estimation performance. Centering around this insight, we first introduce implicit importance sampling. This will provide a basis to understand the IPF and a link to connect the IPF with the KFs. Then, we synthesize the approach of using a bank of nonlinear KFs to approximately execute the IPF. Here, we particularly use the extended and unscented KF (EKF and UKF, respectively), as they are among the most commonly used KFs. Finally, we discuss the connections between the proposed IPF realization with several existing PF methods and validate our approach through simulation.

\section{Implicit Importance Sampling}\label{IIS}
In general, a Bayesian estimation problem involves computation of the conditional expectation of an arbitrary function as follows:
\begin{align}\label{Integral-Bayesian-estimation}
\mathbb{E} \left[g(x) \, | \, y_{1:T}\right] = \int g(x) \p(x | y_{1:T})dx,
\end{align}
where $x \in \mathbb{R}^n$, $g: \mathbb{R}^n \rightarrow \mathbb{R}^m$ is an arbitrary function, and $\p(x | y_{1:T})$ is the posterior probability density function (PDF) of $x$ conditioned on discrete measurements $y_{1:T} = \left\{ y_1, y_2, \ldots, y_T \right\}$~\cite{Sarkka:Cambridge:2013}. The above integral often prohibits  closed-form evaluation for a nonlinear $g$. As a remedy, the Monte Carlo method approximates $\p(x | y_{1:T})$ by a set of particles (samples) and then computes the integral in~\eqref{Integral-Bayesian-estimation} based on the particles. Yet, drawing particles directly from $\p(x | y_{1:T})$ can be challenging.
Importance sampling thus proposes to sample from a different distribution $q(x)$, called importance distribution, and weight the particles accordingly. However, it may require an enormous number of particles to achieve just fair accuracy, because the particles  taken from $q(x)$ may fall into low-probability regions of $\p( x | y_{1:T})$. To ameliorate this problem, a promising approach  is {\em implicit importance sampling}, as introduced below.

 Consider the  task of drawing particles $ x^i$ for $ i=1,\ldots,N $ from high-probability regions of  $\p(x | y_{1:T})$. To this end,  we can  introduce a reference random  vector $\xi$ with a known PDF $p(\xi)$ to  create a map $\xi \rightarrow x^i$ that aligns the target $\p(x^i | y_{1:T})$ with reference $p(\xi)$ density. This notion delivers a two-fold benefit. First, the map will allow to connect high-probability samples of $\xi$ to highly probable $x^i$. Second, one can choose $p(\xi)$ such that it is easy to sample, thus reducing the difficulty of sampling  $x^i$.
Proceeding forward, write $X$ as a shorthand  for $x^i$, and define $F  (X) = - \log \left(\p(X | y_{1:T})  \right)  $ and $s(\xi) =-\log \left ( p(\xi) \right) $. The highest probabilities of $\p(X | y_{1:T})$ and $p(\xi)$ will appear at points where $\min F(X)$ and $\min s(\xi)$ are achieved, respectively. Then, to align $\p(X | y_{1:T})$ with $p(\xi)$ around each other's highest probability point, we can let
\begin{align}\label{map-xi-X}
F  (X) - \min F   (X) = s (\xi) - \min s (\xi),
\end{align}
which implies a map $X(\xi)$. With $X(\xi)$, one can draw a particle from the high-probability regions of $\xi$ and map it to obtain a particle of $X$, which is ensured to lie in the high-probability regions of $\p(X | y_{1:T})$. Here, a trick to pick a high-probable particle for $\xi$ is by using an auxiliary probability distribution that has the same support but concentrates on the high-probability regions of $p(\xi)$.

Next is deciding the sampling weight for $X$. Based on the above, the importance distribution $\pi(\cdot)$ is the PDF of $X(\xi)$ as a mapping of $\xi$ given $p(\xi)$. Assume that $X(\xi)$ is invertible for simplicity of exposition (the result is readily generalizable to the non-invertible case based on~\cite[Corollary 11.3]{Jacod:Springer:2004}). Then, by~\cite[Corollary 11.2]{Jacod:Springer:2004}, it follows that
\begin{align*}
 \pi( X(\xi) )  = \left| J \right|^{-1} \cdot p(\xi),
\end{align*}
where $\pi(X)$ denotes the importance distribution of $X$,   $J = dX(\xi) / d \xi$, and $\left| J\right|$ is the absolute value of J's determinant. The weight of $X$ using~\eqref{map-xi-X} hence is
\begin{align} \label{w_tilde_IIP_compute} \nonumber
\tilde w &= \frac{ \p( X | y_{1:T} ) }{ \pi (X(\xi)) } = \frac{\exp(-F  (X))}{|J|^{-1}\cdot p(\xi)}\\
&= |J | \cdot \exp \left[-\min F  (X)+\min s(\xi ) \right].
\end{align}
For $x^i$, the  normalized weight  is
then given by
\begin{align}\label{w_IIP_compute}
w^i = \frac{\tilde{w}^i}{\sum_{j=1}^{N} \tilde{w}^j} = \frac{|J^i| }{\sum_{j=1}^N |J^j| }.
\end{align}
As a result,  \eqref{Integral-Bayesian-estimation}  is approximated as
\vspace{-5pt}
\begin{align}\label{integral-eval}
\mathbb{E}\left[g(x) \, | \, y_{1:T}\right] \approx \sum_{i=1}^N w^ig(x^i).
\end{align}

In above,~\eqref{map-xi-X}-\eqref{integral-eval} comprise  the method of implicit importance sampling. More details are offered in~\cite{Atkins:MWR:2013,Morzfield:Comms:2015}. While this method promises to use fewer particles, it is computationally challenging to find $\min F(X)$ and then solve \eqref{map-xi-X}, due to the often nonlinear non-convex $F$. To find an alternative way forward, we examine a Gaussian special case here  to gain an insight into tackling this issue.

\noindent\rule{\columnwidth}{1pt}

\noindent{\bf Gaussian Case Study}: Consider the following Gaussian approximation:
\begin{align*}
 p\left(\left[
\begin{matrix}
x \cr y_{1:T}
\end{matrix}
\right] \right) \sim
\mathcal{N}\left(
\left[
\begin{matrix}
\bar m \cr \bar y_{1:T}
\end{matrix}
\right],
\left[
\begin{matrix}
P^x & P^{xy} \cr (P^{xy})^\top & P^y
\end{matrix}
\right]
\right).
\end{align*}
It then follows that
$
 \p( x   |   y_{1:T}) =  \mathcal{N} ( \tilde m, \tilde P^x)$,
where
\vspace{-5pt}
\begin{align*}
\tilde m &=  \bar m  + P^{xy} \left( P^y \right)^{-1} \left( y_{1:T} - \bar y_{1:T} \right),\\
\tilde P^x &= P^x - P^{xy} \left( P^y \right)^{-1} \left( P^{xy}\right)^\top.
\end{align*}
Suppose   $p(\xi) \sim  \mathcal{N} ( 0, I)$ without loss of generality. The map $X(\xi)$ derived from~\eqref{map-xi-X}  is
\vspace{-5pt}
\begin{align*}
X(\xi) = \tilde m +\sqrt{ \tilde P^x }  \xi.
\end{align*}
After a high-probability particle $\xi^i$ is drawn from $p(\xi)$,   it can be mapped by $X(\xi)$ to compute   the particle $x^i$. Its weight  can be determined using~\eqref{w_tilde_IIP_compute}-\eqref{w_IIP_compute}. Specifically,
\begin{align*}
x^i = \tilde m + \sqrt{ \tilde P^x } \xi^i, \ \
w^i = {1 \over N}.
\end{align*}
Based on~\eqref{integral-eval}, we have
\vspace{-10pt}
\begin{align*}
{\mathbb E} \left[ g(x) \, | \, y_{1:T} \right] \approx  {1 \over N} \sum_{i=1}^N g(x^i).
\end{align*}

\vspace{-20pt}

\noindent\rule{\columnwidth}{1pt}

The  above case study indicates that implicit importance sampling is efficiently executable in a Gaussian setting. Extending this insight, we identify that, under the Gaussian approximation, the IPF can be approximately implemented as a bank of KFs.

\section{Approximation of IPF as a Bank of KFs}

Consider the problem of estimating the state of a nonlinear system taking the form
\begin{align}\label{nonlinear-state-space-sys}
\left\{
\begin{aligned}
      x_{k+1} &= f(x_{k}) +w_{k},\\
      y_{k} &= h(x_{k}) +v_{k},\\
\end{aligned}
\right.
\end{align}
where $x_{k} \in \mathbb{R}^{n_x}$ is the state,  $y_{k} \in \mathbb{R}^{n_y}$ is the measurement, and  $w_{k} $ and $v_{k}$ are  zero-mean white noises with covariances $Q$ and $R$, respectively.  For the purpose of state estimation, it is of interest to consider the conditional PDF $\p(x_{0:k} | y_{1:k} )$.  By the Markovian property of  \eqref{nonlinear-state-space-sys},   $\p(x_{0:k} | y_{1:k} )$ satisfies the following recursive relation:
\begin{align*}
\p(x_{0:k}|y_{1:k}) \propto \p(y_{k}|x_{k}) \p(x_{k}|x_{k-1}) \p(x_{0:k-1}|y_{1:k-1}).
\end{align*}
One can assume that   $\p(x_{0:k-1} | y_{1:k-1})$ has been made  available at time $k-1$, as a result of the preceding recursive updates, and that it can be weakly approximated by the empirical distribution of  particles  $ x_{0:k-1}^i$ for $i=1, \ldots, N  $ with   weights $w_{k-1}^i$.  Therefore, to draw   new particles $x_k^i$ at time $k$, we    only need to consider
\begin{align}\label{IPF_recursive_update}
\p(x_{k}^{i}|y_{1:k}) \propto \p(y_{k}|x_{k}^{i}) \p(x_{k}^{i}|x_{k-1}^{i}).
\end{align}
Using $X_k$ as a shorthand for $x_k^i$, we define
\begin{align*}
F_i(X_k)
&= - \log\left( \p(y_{k}| X_k ) \p(X_k |x_{k-1}^{i}) \right),
\end{align*}
To pick $X_k$ in high-probability regions, we  choose  a reference random vector $\xi_k$ with $p(\xi_k)=\exp(-s(\xi_k))$, and following  \eqref{map-xi-X}, let
\begin{align}\label{map-xi-X-PF}
F_i(X_k) - \min F_i(X_k) = s (\xi_k)  -\min s (\xi_k ).
\end{align}
Taking a high-probability sample $\xi_k^i$ from $p(\xi_k)$ and   solving~\eqref{map-xi-X-PF}, we can obtain  a desired particle $x_k^i$, with weight
\begin{align} \label{w_tilde_compute} \nonumber
\tilde {w}_{k}^i  &=   {w}_{k-1}^i   \frac{\p(y_k|X_k ) \p(X_k |x_{k-1}^i)}{\pi( X_k (\xi_k) )}  \\
&=  {w}_{k-1}^i    |J_k^i | \cdot  \exp \left[- \min{F_i(X_k)} + \min{s (\xi_k)} \right] ,
\end{align}
where $J_k^i = d X_k (\xi_k^i) / d \xi_k^i$.
The normalized weight then is
\begin{align}\label{weight-update}
w_k^i =    \frac{\tilde{w}_k^i}{\sum_{j=1}^{N} \tilde{w}_k^j}.
\end{align}
Here,~\eqref{map-xi-X-PF}-\eqref{weight-update} illustrate the particle computation and weighting underlying the IPF. However, a difficulty for the implementation lies in finding $\min F_i(X_k)$ and solving the algebraic equation~\eqref{map-xi-X-PF}. Even though the literature presents several useful methods based on iterative computation~\cite{Chorin:Springer:2013}, their computational complexity is still not-trivial. To address the challenge, we   take inspiration from the Gaussian Case Study in Section~\ref{IIS} to  solve~\eqref{map-xi-X-PF}.

First, let us   introduce a   Gaussian approximation:
\begin{align}\label{Gaussian_assumption_IPF_KF}
p \left(\left. \begin{bmatrix}
x_k^i \\ y_{k}
\end{bmatrix} \, \right| \, x_{k-1}^i\right) \sim \mathcal{N}\left(\begin{bmatrix}
\bar{m}_k^i \\ \bar{y}_k^i
\end{bmatrix}, \begin{bmatrix}
\bar{P}_k ^{i} & \bar{P}_k^{xy,i} \\  \left(\bar{P}_k^{xy,i} \right)^\top & \bar{P}_k ^{y,i}
\end{bmatrix}\right).
\end{align}
By~\eqref{Gaussian_assumption_IPF_KF}, we have
\begin{align}\label{Particle_update}
\p( x_k^i  |  y_k, x_{k-1}^i ) \sim \mathcal{N} \left(  \tilde m_k^i , \tilde P_k^i  \right),
\end{align}
where
\begin{subequations}\label{KF-update}
\begin{align}\label{KF-update-a}
\tilde m_k^i  &= \bar m_k^i +  \bar{P}_k^{xy,i} \left( \bar{P}_k ^{y,i} \right)^{-1} \left( y_k - \bar y_k^i \right)   ,\\ \label{KF-update-b}
\tilde P_k^i &= \bar{P}_k ^{i}  - \bar{P}_k^{xy,i} \left( \bar{P}_k ^{y,i} \right)^{-1} \left(\bar{P}_k^{xy,i} \right)^\top.
\end{align}
\end{subequations}
Assuming $ \p(y_k|x_k^i) = \p(y_k|x_k^i,x_{k-1}^i) $ and inserting it into~\eqref{IPF_recursive_update}, we obtain
\begin{align}\label{IPF_recursive_update_new}
\p(x_{k}^{i}|y_{1:k}) \propto     \p(x_k^i|y_k,x_{k-1}^i) \p(y_k|x_{k-1}^i)
\end{align}
Given~\eqref{IPF_recursive_update_new}, we use~\eqref{map-xi-X-PF}  to align  $\p(x_{k}^{i}|y_{1:k})$ with  the high-probability regions of  a standard Gaussian random vector  $\xi_k$, with $p(\xi_k) \sim \mathcal{N}(0,I)$.  Then,  referring to the Gaussian Case Study in Section~\ref{IIS},
the desired particle $x_k^i$   can be expressed as
\begin{align}\label{particle-compute}
x_k^i = \tilde m_k^i +  \sqrt{\tilde{P}_{k}^{i}}\xi_k^i,
\end{align}

\vspace{35pt}

\noindent\rule{\columnwidth}{1pt}

\vspace{-10pt}

\noindent{\bf Algorithm 1: The {KF-IPF} Algorithm}
\vspace{-5pt}
\begin{itemize}

\item Initialize a set of particles $x_0^i$ for $i=1,\ldots,N$ with   $w_0^i$ and $P_0^i$ at $k=0$
\item Implement the following for $x_k^i$ for $i=1,\ldots,N$   from $k=1, \ldots,T$

\begin{itemize}

\item[-] Compute $\bar m_k^i$ and $\bar P_k^i$ using $x_{k-1}^i$ and $P_{k-1}^i$ {\em via} KF prediction

\item[-] Compute $\bar y_k^i$, $\bar P_k^{y,i}$ and $\bar P_k^{xy,i}$ as in KF prediction

\item[-] Compute $\tilde m_k^i$ and $\tilde P_k^i$  {\em via} KF update in~\eqref{KF-update-a}-\eqref{KF-update-b}

\item[-] Compute $x_k^i$ and $w_k^i$  {\em via}~\eqref{particle-compute}-\eqref{Cov_perturb}

\item[-] Do resampling if necessary

\item[-] Output the state estimate

\item[-] $k \leftarrow k+1$

\end{itemize}

\end{itemize}

\vspace{-15pt}

\noindent\rule{\columnwidth}{1pt}

\noindent where $\xi_k^i$ is a high-probability sample of $p(\xi_k)$.
Based on~\eqref{w_tilde_compute}-\eqref{weight-update}, its weight is computed as
\vspace{-2pt}
 \begin{align}  \label{weight-compute}
{w}_k^i &= \frac {w_{k-1}^i \p(y_k|x_{k-1}^i)  }{\sum_{j=1}^N  w_{k-1}^j  \p(y_k|x_{k-1}^j)}.
 \end{align}

In practice, it is important that $\xi_k^i$ are drawn from the high probability regions of $p(\xi_k)$. An implementation trick to ensure this is to sample from $\mathcal{N}(0, \alpha I )$ with $0< \alpha \ll 1$. We emphasize that $\alpha$ can be considered  as a tunable parameter to dictate  the desired high probability region of $p(\xi_k)\sim\mathcal{N}(0, I)$. It hence must take a small enough value to ensure that only highly probable particles with respect to $\mathcal{N}(0, I)$ are drawn. This is especially true when only a small number of particles are used. The covariance associated to $x_k^i$ is given as
\vspace{-5pt}

\begin{equation}\label{Cov_perturb}
P_k^i = \tilde{P}_k^i .
\end{equation}

Note that the  formulation in~\eqref{Gaussian_assumption_IPF_KF}-\eqref{KF-update}   is equivalent to the well-known KF update procedure. This connection suggests the viability of  approximately implementing the IPF as a bank of parallel KFs applied to individual particles.
Specifically, given $x_{k-1}^i$ at time $k-1$,  one can compute $\bar m_k^i$ and $\bar P_k^i$ {\em via} KF prediction and then compute $\tilde m_k^i$ and $\tilde P_k^i$ {\em via} KF update~\cite{Fang:JAS:2018}; going further, $\tilde m_k^i$ combines with $\xi_k^i$ to generate $x_k^i$. Different nonlinear KFs can be leveraged to enable the   prediction-update   of the particles. Here, we highlight the use of the EKF and UKF to approximate the mean and covariance statistics of the Gaussian approximation in~\eqref{Gaussian_assumption_IPF_KF}. The EKF and UKF have found wide use with proven effectiveness for various  systems, and their detailed equations can be found in~\cite{Fang:JAS:2018}. As extensions of the standard linear KF, both EKF and UKF at their core aim to track the conditional statistics (mean and covariance) of a nonlinear system's state recursively using accumulating measurement data. To achieve this end,   EKF linearizes the system's nonlinear functions, and UKF leverages the unscented transform.

Summarizing the above, Algorithm 1 outlines the KF-IPF algorithm  framework. Based on this framework,
approximate realizations based on the EKF or UKF can be readily derived and named as E-IPF and U-IPF, respectively. The details are omitted here for the sake of space.

The proposed KF-IPF algorithm provides an amenable way to implement the IPF by approximating it as a bank of KFs. The algorithm circumvents the need to numerically find $\min F_i(X_k)$ and solve the nonlinear algebraic equation~\eqref{map-xi-X-PF} with computationally intensive numerical schemes in existing IPF literature. By design, the KF-IPF can also exploit the advantages of a nonlinear KF. For instance, it can achieve derivative-free computation if the UKF is used. These features will enable it to facilitate the computation of the IPF considerably. Further, the algorithm inherits the capability of the IPF in mitigating the issue of particle degeneracy. Also, note that the development of the KF-IPF algorithm uses the Gaussian approximation in~\eqref{Gaussian_assumption_IPF_KF}, which plays a utilitarian role in solving~\eqref{map-xi-X-PF} {\em via} the KF update. In effect, the approximation serves to facilitate computation. Despite it, the KF-IPF algorithm, as an IPF method, is still well applicable to non-Gaussian estimation problems. The simulation study in Section~\ref{num-sim} further highlights this.

The literature also includes PF methods based on banks of EKF and UKF~\cite{Doucet:SC:2000,Merwe:NIPS:2001}, referred to as EPF or UPF, respectively. However, they are different from the proposed KF-IPF algorithm on several aspects. First, the EPF/UPF is designed to use a bank of EKF/UKF to approximate PF's optimal importance distribution, $\p(x_k | x_{k-1}, y_k)$. However, the KF-IPF is based on the notion of {\em implicit importance sampling}, which looks for high-probability regions of the target distribution and then locates particles in these regions. Second, in terms of implementation, in the sampling step~\eqref{particle-compute}, the EPF/UPF draws particles from the Gaussian approximation made by the EKF/UKF, which is analogous to the case where the KF-IPF directly draws $\xi_k^i$ from $\mathcal{N}(0, I)$. However, a significant risk, in this case, is that $x_k^i$ will fall in the low-probability regions when a marginally probable $\xi_k^i$ is taken---this is also likely why the EPF and UPF sometimes fail to perform adequately, as noticed in~\cite{Sarkka:Cambridge:2013}. By contrast, the actual implementation of the KF-IPF requires focusing particles only on the high-probability regions. For instance, to ensure this, we suggested drawing particles $\xi_k^i$ from $\mathcal{N}(0,\alpha I)$ with $0<\alpha \ll 1$ so that the resultant $x_k^i$ is highly probable. Finally, in the particle weighting step~\eqref{weight-compute}, the KF-IPF utilizes $\p(y_k | x_{k-1}^i)$ instead of $\p(y_k | x_{k}^i)$ as in the case of the EPF/UPF. To sum up, the KF-IPF is distinct from the EPF/UPF, while contributing a new insight into the broad notion of using KF banks to implement the PF.

\section{Numerical Simulation}\label{num-sim}

In this section, we assess the performance of the KF-IPF algorithm on the 40-dimensional Lorenz'96 model~\cite{lorenz96}:
\begin{align}\label{Lorenz96}
 \dot x_j &= f(x, F) = \big(x_{j+1} -x_{j-2} \big)x_{j-1} - x_j +F,
\end{align}
where $j = 1,\hdots,n_x = 40$ is the dimension index,  $x_{-1} = x_{n_x-1}$, $x_0 = x_{n_x}$, $x_{n_x+1} = x_1$, $F=5$. The simulation run of~\eqref{Lorenz96} is based on discretization using the explicit fourth-order Runge-Kutta method.  The discret-time system with process noise $w$ is expressed as:
\begin{align*} \label{Lorez96-discrete}
x_{j,k+1} = x_{j,k} + \frac{\Delta t}{6} \big(h_1 + 2h_2 + 2h_3 + h_4\big) + w_{j,k},
\end{align*}
where
\begin{align*}
h_1 &= f\left(x_{j,k}, F\right), \quad   h_2 = f\left(x_{j,k} + \Delta t \cdot h_1/2, F\right),\\
h_3 &= f\left(x_{j,k} + \Delta t\cdot h_2/2,  F\right),  \ h_4 = f\left(x_{j,k} + \Delta t\cdot h_3, F\right),
\end{align*}
$w_{j,k}\sim \mathcal{U}(-0.5, 0.5)$ is added as the i.i.d. process noise,  and $\Delta t = 0.01$ is the discretization time-step. In addition, we consider a nonlinear measurement function that partially observes
every other component of the state:
\begin{align*}
y_{l,k} = x_{2l-1,k}  + \sin(x_{2l-1,k} ) + v_{l,k} ,
\end{align*}
where $l=1,\ldots,n_y = 20$, and $v_l \sim \mathcal{U}(-0.5, 0.5)$. We evaluate the state estimation performance of the E-IPF and U-IPF for the above system by the root mean squared error (RMSE) metric over 50 Monte Carlo (MC) runs. Each run is initialized randomly and with a bias to the true state. Meanwhile, we run the IPF based on numerical iteration (I-IPF) in~\cite{Chorin:ESAIM:2012}, EPF, and UPF for comparison. The PDF of the reference random variable $p(\xi_k)$ for all the considered IPFs is set to $\mathcal{N}(0, \alpha I)$, where $\alpha = 0.05$. Further, the individual Kalman filters use the covariances of $w$ and $v$ as $Q$ and $R$, respectively, in their estimation run.
\begin{figure}[t]\centering
\includegraphics[width=0.49\textwidth, trim={4cm 8.6cm 4.5cm 9cm},clip]{Lorenz_96.pdf}
\caption{RMSE performance of  the EPF/E-IPF (with 1,000 particles) and U-IPF (with only 10 particles) versus UPF/I-IPF (with 100 particles). The shaded regions correspond to the $\pm \sigma$ bounds by the filters over 50 Monte Carlo runs.}
\label{fig:RMSE}
\end{figure}

Fig.~\ref{fig:RMSE} shows the estimation performance of the different filters using different numbers of particles. The observations are summarized as follows. First, the KF-IPF algorithm provides advantageous estimation. Specifically, the E-IPF performs better than the EPF in terms of accuracy, with both using 1,000 particles; further, the U-IPF with only ten particles delivers considerably higher accuracy than the UPF with 100 particles. Second, compared to the I-IPF, the KF-IPF represents a more favorable and promising approach. Even though the I-IPF is understandably more accurate than the E-IPF (since it is based on iterative linearization and solution rather than the first-order linearization involved in the E-IPF), the U-IPF consistently outperforms the I-IPF, and this comes at using just ten particles rather than 100 by the I-IPF. The observations also suggest the need to use a good KF to well execute the KF-IPF algorithm. The simulation was run using MATLAB R2020b on a workstation equipped with a 3.5GHz Intel Core i9-10920X CPU, 128GB of RAM. Table~\ref{table:1} summarizes the average computation time of a MC run for the considered filters. Overall, the computation time increases with the number of particles. The U-IPF and I-IPF are computationally the fastest among the filters. The simulation results indicate the advantages and promise of the KF-IPF algorithm for nonlinear non-Gaussian estimation.
 \begin{table}[t]\caption{Average Computation Time}\label{table:1}
\setlength{\tabcolsep}{4pt}
\centering
        \begin{tabular}{@{}lllllll@{}}
        \toprule
        & EPF  & E-IPF &  UPF & I-IPF & U-IPF\\
	\midrule
	$N$ & 1,000 & 1,000 & 100 & 100 & 10 \\
	\midrule
        Time (s) & 3,462 &3,362 & 652 & 578.9 & 46.4 \\
        \bottomrule
        \end{tabular}
\end{table}
\vspace{-3pt}
\section{Conclusion}
The IPF has emerged as a vibrant tool for dealing with particle degeneracy confronting the PF. However, the present IPF methods are computationally cumbersome, which may limit their application. In this paper, we proposed approximately implementing the IPF framework {\em via} nonlinear KFs for higher computational tractability. We presented implicit importance sampling and then, on this basis, derived the IPF realization based on a bank of nonlinear KFs. We discussed the EKF and UKF realizations of the IPF. Further, we examined how the proposed realization relates to some existing PF methods. Finally, we demonstrated the effectiveness of KF-based IPF through the use of a simulation.

\section*{Acknowledgments}

The authors would like to thank the anonymous reviewers for their constructive comments. I. Askari and H. Fang were partially supported by the U.S. National Science Foundation under Awards CMMI-1763093 and CMMI-1847651.   X. Tu was partially supported by the U.S. National Science Foundation under Award DMS-1723066.

 \balance
\bibliographystyle{elsarticle-num}


\section{Introduction}
The inherently unambiguous nature of parsing expression grammars (PEGs) makes them an attractive choice for modelling structured text such as programming languages and computing data formats, but in practice the superior performance of parsers based on context-free grammars (CFGs) has led to CFGs being more-widely used, despite the difficulty of disambiguating them.
This work is an initial effort toward a unified framework that provides the advantages of both grammar formalisms: it is an algorithm adopted from an efficient, general-purpose CFG parser that supports PEG semantics without discarding support for the unordered choice operator of CFGs in cases where that ambiguity may be desirable.

More specifically, this paper presents a modification of the GLL parser-generator of Scott \& Johnstone\cite{SJ10,SJ16}.
The key contribution is the \emph{FailCRF} data structure, which adds a failure path to Scott \& Johnstone's call-return forest; the addition of a failure path allows the lookahead and ordered choice operations of the PEG formalism to be supported.
The authors have extended Ackerman's GoGLL\cite{Ack19} parser-generator to implement this new algorithm, adding syntactic sugar for common PEG operators and modifying the lexer algorithm to allow the PEG parser to override the usual maximal-munch rule.

\section{Parsing Expression Grammars}
The primary difference between parsing expression grammars and the more familiar context-free grammars is \emph{ordered choice}: PEGs, as a formalism of recursive-descent parsing introduced by Ford \cite{For02}, do not try subsequent alternatives of an alternation if an earlier alternative matches.
The other significant difference between the PEG and CFG formalisms are the PEG \emph{lookahead} expressions, $\nl\alpha$ and $\&\alpha$, which match only if the subexpression $\alpha$ does not (resp. does) match, but consume no input regardless.
These lookahead operators provide the infinite lookahead of the PEG formalism.
The other fundamental PEG operators act much like their CFG equivalents, and are described in Fig.~\ref{expr-fig} as functions over an input string $s$ drawn from some alphabet $\Sigma$ producing either a (matching) suffix of $s$ or the special value $\fail \not\in \Sigma^*$.
In summary, the \emph{string literal} $u$ matches and consumes the string $u$, the \emph{empty expression} $\varepsilon$ always matches without consuming anything, while the \emph{failure expression} $\varnothing$ never matches.
A \emph{nonterminal} $A$ is replaced by the parsing expression $\Rule(A)$ it corresponds to.
The \emph{sequence} expression $\alpha\beta$ matches $\alpha$ followed by $\beta$, while the \emph{ordered choice} expression $\alpha/\beta$ only tries $\beta$ if $\alpha$ does not match.
To differentiate CFG \emph{unordered choice}, it is represented in this paper as $\alpha|\beta$.

\begin{figure}
	\centering
	\begin{equation*}
	\begin{aligned}[c]
	u(s)            & = \begin{cases} v     & s = uv \\
	                                  \fail & \text{otherwise} \end{cases} \\
	\varepsilon(s)  & = s \\
	\varnothing(s)  & = \fail \\
	\nl\alpha(s)    & = \begin{cases} s     & \alpha(s) = \fail \\
	                                  \fail & \text{otherwise} \end{cases} \\
    \&\alpha(s)     & = \begin{cases} s     & \alpha(s) \neq \fail \\
                                      \fail & \text{otherwise} \end{cases}
	\end{aligned}
	~~~~
	\begin{aligned}[c]
	A(s)            & = (\Rule(A))(s) \\
	\alpha\beta(s)  & = \begin{cases} \beta(\alpha(s)) & \alpha(s) \neq \fail \\
	                                  \fail            & \text{otherwise} \end{cases} \\
	\alpha/\beta(s) & = \begin{cases} \alpha(s)  & \alpha(s) \neq \fail \\
	                                  \beta(s)   & \text{otherwise} \end{cases}
	\end{aligned}
	\end{equation*}
	\caption[Expression definitions]{Formal definitions of parsing expressions} \label{expr-fig}
\end{figure}

\section{GLL Parsing}
\emph{Generalized LL} (GLL) parsing, introduced by Scott \& Johnstone\cite{SJ10,SJvB19}, extends the power of LL parsing to all CFGs through use of a \emph{call-return forest} (CRF) to represent the recursive-descent call stack of the LL parsing algorithm.
For efficiency, the CRF is implemented using the \emph{graph-structured stack} (GSS) data structure introduced by Tomita\cite{Tom85} for the GLR parsing algorithm.
The gist of the GLL approach is that each CRF node represents a function call (equivalently, nonterminal invocation) in a recursive-descent parse, and includes an input position, a nonterminal to match, and a grammar slot to return to on completion.
The graph structure of this stack comes from a dynamic de-duplication of CRF nodes which share a nonterminal and input position, changing a stack data structure into a directed acyclic graph (DAG).
The GLL algorithm keeps a queue of CRF nodes which are pending parsing, and handles the nondeterminism of unordered choice by enqueuing a CRF node for each choice.

Scott \etal\cite{SJvB19} introduced \emph{binary subtree representation} (BSR) sets as an output format to represent nonterminal matches in GLL.
The essential insight is that, while the traditional \emph{shared packed parse forest} (SPPF)\cite{Tom85} data structure representing possible parse trees requires significant complication in the parser algorithm to properly store and update edges between parse tree nodes, those edges can be efficiently reconstructed from an indexed set of edgeless parse-tree nodes (the BSR set) with minimal added information.

A BSR element is a 4-tuple containing a \emph{grammar slot} $X ::= \alpha\theta\cdot\beta$, and three input indices $i$, $j$, and $k$, $i \leq j \leq k$.
The BSR element represents a successful match of the nonterminal $X$ up to the end of $\theta$, the single terminal or nonterminal immediately before the dot of the grammar slot; $i$ is the input index where $X$ began to match, $j$ is the index where $\theta$ began to match, and $k$ is the index where $\theta$ finished matching.
Note that if $\beta = \varepsilon$, the BSR node represents a complete match of $X$.
Parse trees can be straightforwardly reconstructed from BSR sets: a predecessor of a BSR element $(X ::= \alpha\theta\cdot\beta, i, j, k)$ is any element $(X ::= \alpha\cdot\theta\beta, i, \ell, j)$, while its child where $\theta$ is some nonterminal $A$ is any element $(A ::= \delta\cdot, j, m, k)$. Successor and parent elements can be defined analogously.

\section{Parsing Expression GLL}
The \emph{Parsing Expression GLL} (PEGLL) algorithm introduced in this paper uses similar data structures and abstractions as GLL for CFGs.
The main loop is outlined in Figure~\ref{main-loop-algo}; it first initializes an empty queue $R$ of slot descriptors to parse, an empty cache $U$ of previously seen slot descriptors, and an empty set $T$ of BSR elements to report.
It then queues the start rule of the grammar at input position 0 for parsing, parses each descriptor, and completes by returning whether or not the start rule matched at position 0.
Note that (unlike CFGs), PEGs match prefixes of their input, so the start rule may only consume the input up to some index $k$; if this behavior is not desired, a match can be returned only if $k$ is the length of the input string.

\begin{figure}
\caption{PEGLL main loop} \label{main-loop-algo}
\begin{algorithmic}
\State $R \gets \emptyset$, $U \gets \emptyset$, $T \gets \emptyset$
\State \Call{addNt}{$S$, 0}
\While{$R \neq \emptyset$}
    \State $(L, c_U, c_I) \gets R$.\Call{remove}{}
    \State $t \gets$ \Call{tokens}{$c_I$} $t' \gets t$
    \Loop
        \Switch{$L$}
            \State $\langle$ generate code for each rule $R$ $\rangle$
        \EndSwitch
    \State nextSlot: \EndLoop
\State nextDesc: \EndWhile
\If{$\exists \alpha, j, k, (S ::= \alpha \cdot, 0, j, k) \in T$}
    \State \Return match at maximal such $k$
\Else
    \State \Return $\fail$;
\EndIf
\end{algorithmic}
\end{figure}

The primary thing the loop in Figure~\ref{main-loop-algo} does is dispatch the current descriptor to the code which executes its parse.
The code for each nonterminal may be generated according to the patterns in Figures~\ref{nt-pattern} and~\ref{empty-nt-pattern}.
This parser-generator assumes any ordered choice expressions are at the very top level of a nonterminal (parenthesized subexpressions are added as syntactic sugar, see Section~\ref{syntax-sec}).
A label is then generated for each alternate, failing over to the next alternate if it does not match, with a synthesized failure alternate at the end of each alternation.

\begin{figure}
\caption[Non-terminal code pattern]{Pattern for a nonterminal $X ::= \tau_1 / \ldots / \tau_p, p \geq 1$} \label{nt-pattern}
\begin{algorithmic}
\Case{$X ::= \cdot \tau_1$}
	\State $\langle$ generate code for $X ::= \cdot \tau_1$ with failure path $X ::= \cdot \tau_2$ $\rangle$
	\State \Call{rtn}{$X$, $c_U$, $c_I$}; \Goto{nextDesc}
\EndCase
\State $\vdots$
\State ~
\Case{$X ::= \cdot \tau_p$}
	\State $\langle$ generate code for $X ::= \cdot \tau_p$ with failure path $X ::= \cdot \varnothing$ $\rangle$
	\State \Call{rtn}{$X$, $c_U$, $c_I$} \Goto{nextDesc}
\EndCase
\Case{$X ::= \cdot \varnothing$}
	\State \Call{rtn}{$X$, $c_U$, $\fail$} \Goto{nextDesc}
\EndCase
\end{algorithmic}
\end{figure}

\begin{figure}
\caption[Empty nonterminal code pattern]{Pattern for a nonterminal $X ::= \varepsilon$} \label{empty-nt-pattern}
\begin{algorithmic}
\Case{$X ::= \cdot \varepsilon$}
	\State $T \gets T \cup {(X ::= \varepsilon \cdot, c_I, c_I, c_I)}$
	\State \Call{rtn}{$X$, $c_U$, $c_I$} \Goto{nextDesc}
\EndCase
\end{algorithmic}
\end{figure}

The \textsc{rtn}$(X, c_U, c_I)$ function in PEGLL code is detailed below in Figure~\ref{call-rtn-code}, but its essential purpose is to modify the ``call stack'' in the CRF graph consistently with a recursive call to the nonterminal $X$ at position $c_U$ returning at position $c_I$.

To parse each of the sequence expressions $\tau_i$ inside a nonterminal alternative, PEGLL repeatedly executes the appropriate code for each expression in the sequence, moving to the failure path if that expression does not work.
The code patterns for the sequence expression are in Figure~\ref{seq-expr-pattern}, while the code patterns for each atomic expression are in Figure~\ref{atom-expr-pattern}.
Terminal matches advance the input index $c_I$ to the right extent $r$ of the token and find the new tokens $t$ at that position, while mid-sequence errors reset $c_I$ and $t$ to their initial values in the descriptor, $c_U$ and $t'$.
Note also that nonterminal calls preserve the failure path from their caller, and in particular that negative-lookahead expressions swap the success and failure results of the call.

\begin{figure}
\caption[Sequence expression code pattern]{Pattern for a sequence $X ::= \cdot x_1 \ldots x_d$ with failure path $L_f$} \label{seq-expr-pattern}
\begin{algorithmic}
\State $r \gets$ \Call{testSelect}{$X ::= \cdot x_1 \ldots x_d$, $t$}
\If{\Call{testSelect}{} failed}
	\State report error
	\State $(L, c_I, t) \gets (L_f, c_U, t')$ \Goto{nextSlot}
\EndIf
\State $\langle$ generate code for $X ::= x_1 \cdot x_2 \ldots x_d$ with failure path $L_f$ $\rangle$
\State $\vdots$
\State $\langle$ repeat for remaining atoms $x_2 \ldots x_d$ in sequence $\rangle$
\end{algorithmic}
\end{figure}

\begin{figure}
\caption[Atomic expression code patterns]{Pattern for an atomic expression $X ::= \alpha \varphi \cdot \beta$ with failure path $L_f$} \label{atom-expr-pattern}
\begin{algorithmic}
\StaticCase{$\varphi = a$}
	\State $R \gets R \cup (X ::= \alpha a \cdot \beta, c_U, c_I, r)$
	\State $c_I \gets r$
	\State $t \gets$ \Call{tokens}{$c_I$}
\EndStaticCase{}
\StaticCase{$\varphi = Y$}
	\State \Call{call}{$X ::= \alpha ~Y \cdot \beta$, $L_f$, $Y$, $c_U$, $c_I$} \Goto{nextDesc}
\EndStaticCase{\textbf{case} $X ::= \alpha ~Y \cdot \beta$:}
\State ~
\StaticCase{$\varphi = \& Y$}
	\State \Call{call}{$X ::= \alpha ~\& Y \cdot \beta$, $L_f$, $Y$, $c_U$, $c_I$} \Goto{nextDesc}
\EndStaticCase{\textbf{case} $X ::= \alpha ~\& Y \cdot \beta$:}
\State ~
\StaticCase{$\varphi = \nl Y$}
	\State \Call{call}{$L_f, X ::= \alpha ~\nl Y \cdot \beta$, $Y$, $c_U$, $c_I$} \Goto{nextDesc}
\EndStaticCase{\textbf{case} $X ::= \alpha ~\nl Y \cdot \beta$:}
\end{algorithmic}
\end{figure}

Together, the helper functions \textsc{call} and \textsc{rtn} simulate the call stack of a recursive-descent PEG parser; full code is in Figure~\ref{call-rtn-code}.
The significant difference between the \emph{FailCRF} of PEGLL and the \emph{CRF} of Scott~\etal \cite{SJvB19} is that edges in the \emph{FailCRF} are labelled as either ``match'' edges or ``fail'' edges, representing successful and unsuccessful return paths, respectively.
The motivation for this modification is handling negative lookahead expressions --- $\nl X$ matches only if $X$ does not, and thus the CRF needs to encode the failure of $X$ as the trigger of a move from slot $Y ::= \alpha \cdot \nl X \beta$ to slot $Y ::= \alpha \nl X \cdot \beta$.
Fail edges are also used to encode the ordered choice rule; while the GLL algorithm attempts to parse all alternates of a nonterminal concurrently, PEGLL attempts one at a time, following its fail edge to the next alternate if that one fails.

In more detail, the \textsc{call}$(L_m, L_f, X, i, j)$ function enqueues parsing of nonterminal $X$ at position $j$, returning to slot $L_m$ on match or $L_f$ on failure, where $L_m$ or $L_f$ began parsing at position $i$.
The \textsc{rtn}$(X, j, h)$ function reports the result of parsing nonterminal $X$ begining at position $j$, where $h$ is either the last consumed position or the special value $\fail$ indicating that the parse did not succeed.
For convenience, this presentation assumes the call-return forest (CRF) is stored in a global cache.
The ``popped cache'' of previously-parsed results is included to allow updating of previous parse results when they are used in a new context and also assumed to be global.

\begin{figure}
\caption{Code for \textsc{call} and \textsc{rtn} functions} \label{call-rtn-code}
\begin{algorithmic}
\Function{call}{slot $L_m$, slot $L_f$, nonterminal $X$, int $i$, int $j$}
	\State $u_m \gets$ CRF node $(L_m, i)$, created if not in cache
	\State $u_f \gets$ CRF node $(L_f, i)$, created if not in cache
	\State $v \gets$ CRF node $(X, j)$, created if not in cache
	\If{$v$ was not previously in cache}
		\State add a match edge from $v$ to $u_m$ and a fail edge from $v$ to $u_f$
		\State \Call{addNt}{$X$, $j$}
	\Else
		\If{there is not a match edge from $v$ to $u_m$}
			\State add a match edge from $v$ to $u_m$
			\ForAll{$(X, j, h), h \neq \fail$ in popped cache}
				\State \Call{addMatch}{$L_m$, $i$, $j$, $h$}
			\EndFor
		\EndIf

		\If{there is not a fail edge from $v$ to $u_f$}
			\State add a fail edge from $v$ to $u_f$
			\ForAll{$(X, j, \fail)$ in popped cache}
				\State \Call{addFail}{$L_f$, $i$, $j$}
			\EndFor
		\EndIf
	\EndIf
\EndFunction
\State ~
\Function{rtn}{nonterminal $X$, int $j$, int $h$}
	\If{$(X, j, h)$ not in popped cache}
		\State add $(X, j, h)$ to popped cache
		\ForAll{children $(L, i)$ of $(X, j)$ in the CRF}
			\If{$h \neq \fail$}
				\State \Call{addMatch}{$L$, $i$, $j$, $h$}
			\Else
				\State \Call{addFail}{$L$, $i$, $j$}
			\EndIf
		\EndFor
	\EndIf
\EndFunction
\end{algorithmic}
\end{figure}

In addition to call-stack management, a PEGLL parser-generator depends on helper functions to fill the descriptor queue $R$, descriptor cache $U$, and the BSR set $T$ representing the parse forest.
The descriptor queue management functions are all in Figure~\ref{add-match-fail}.
\textsc{addMatch} and \textsc{addFail} enqueue the next descriptor after a nonterminal match or failure, accounting for the fact that lookahead expressions consume no input and adding a BSR element for matches.
\textsc{addNt} is a convenience function which enqueues the first alternate of a nonterminal in the descriptor queue.
\textsc{addDesc} checks if a given descriptor has already been processed, adding it to the descriptor queue and the cache of processed descriptors if not.

\begin{figure}
\caption{Code for descriptor queue management functions} \label{add-match-fail}
\begin{algorithmic}
\Function{addMatch}{slot $L$, int $i$, int $j$, int $h$}
	\State $T \gets T \cup {(L, i, j, h)}$
	\If{there is a lookahead expression before the dot in $L$}
		\State \Call{addDesc}{$L$, $i$, $j$}
	\Else
		\State \Call{addDesc}{$L$, $i$, $h$}
	\EndIf
\EndFunction
\State ~
\Function{addFail}{slot $L$, int $i$, int $j$}
	\If{there is a lookahead expression before the dot in $L$}
		\State \Call{addDesc}{$L$, $i$, $j$}
	\Else
		\State \Call{addDesc}{$L$, $i$, $i$}
	\EndIf
\EndFunction
\State ~
\Function{addNt}{nonterminal $X$, int $i$}
	\State $L \gets$ initial slot of first alternate of $X$
	\State \Call{addDesc}{$L$, $i$, $i$}
\EndFunction
\State ~
\Function{addDesc}{slot $L$, int $i$, int $h$}
	\If{$(L,i,h) \not\in U$}
		\State $R \gets R \cup {(L,i,h)}$
		\State $U \gets U \cup {(L,i,h)}$
	\EndIf
\EndFunction
\end{algorithmic}
\end{figure}

\subsection{Lexing}
Traditional maximal-munch lexers interact poorly with PEG parsers; in particular, the recursive-descent structure of a parsing expression grammar may impose context-sensitive priorities between tokens, and the usual scannerless design of a PEG may result in difficult-to-decompose token sets.
Nonetheless, experience \cite{Ack19,Lau19} has shown that a separate lexing pass is a useful performance optimization.
As such, the PEGLL parser-generator uses a regular-expression lexer modified from the traditional maximal-munch approach \cite{Aetc07}.
Rather than returning a single maximal-munch token at a given input position, as in the classical algorithm , the PEGLL lexer returns (from the \textsc{tokens}$(i)$ function) a map of all the tokens which match starting at position $i$ to their greatest right extent, defering the choice of which token to actually match to the PEG rules in the parser.
The \textsc{testSelect} function, described in Figure~\ref{test-select} in the parser is used to check if any tokens in the FIRST set of the expression \cite{Red09} are also present in the token set returned from the lexer.

\begin{figure}
\caption[testSelect]{Code for \textsc{testSelect} function} \label{test-select}
\begin{algorithmic}
\Function{testSelect}{slot $L$, token set $T$, int $c_I$}
	\State $b \gets -1$
	\If{$L$ is nullable} $b \gets c_I$ \EndIf
	\ForAll{$(t,r)$ in $T$}
		\If{$r > b$ \textbf{and} $t \in FIRST(L)$} $b \gets r$ \EndIf
	\EndFor
	\State \Return{$b$}
\EndFunction
\end{algorithmic}
\end{figure}

\subsection{Syntactic Sugar} \label{syntax-sec}

The preceding discussion has outlined how PEG semantics are supported within the framework of the GLR algorithm; however, support for the fundamental PEG operators does not cover how several common PEG idioms are supported.
Much like the EBNF notation for CFGs, PEGs typically support the ``syntactic sugar'' operators defined in Figure~\ref{sugar-fig} for optional and repeated matches.
The $\alpha?$, $\alpha^*$, and $\alpha^+$ operators have their usual semantics of match $\alpha$ zero-or-one, zero-or-more, and one-or-more times, respectively, though note that the ordered-choice semantics of PEGs imply greedy match.
The PEGLL parser supports expressions using these operators by replacing such expressions during AST generation with a fresh nonterminal defined as in Figure~\ref{sugar-fig}.

\begin{figure}
	\centering
	\begin{equation*}
	\begin{aligned}[c]
	\Rule_{\alpha?} & = \alpha / \varepsilon \\
	\Rule_{\alpha^*} & = \alpha \Rule_{\alpha^*} / \varepsilon \\
    \Rule_{\alpha^+} & = \alpha \Rule_{\alpha^*} \\
	\end{aligned}
	\end{equation*}
	\caption[Syntactic sugar definitions]{Each syntactic sugar expression $\varphi$ is replaced by a fresh nonterminal $\Rule_\varphi$ defined as in this table.} \label{sugar-fig}
\end{figure}

\subsection{Unordered Choice}

PEGLL also supports an unordered-choice operator; the main challenge in implementing unordered choice in PEGLL is determining when to trigger the failure path for a nonterminal.
PEGLL introduces the failure paths in the FailCRF data structure to support the sequential nature of ordered choice -- when one alternate fails, it enqueues the next, and when the final alternate fails the algorithm marks the failure of the entire nonterminal.
Unordered choice, by contrast, has concurrent semantics; all alternates can conceptually be executed in parallel (as in Scott \& Johnstone's GLL \cite{SJ10,SJ16}), and the nonterminal only fails if \emph{all} alternates fail.

The approach taken by PEGLL is formalized in Figure~\ref{unordered-choice-algo}, but the essential idea is to duplicate all alternates of an unordered choice into a \emph{pass} alternate and a \emph{fail} alternate.
Matching of an unordered choice begins with the first fail alternate, and proceeds sequentially through the alternates in source order, as in PEGLL's handling of ordered choice.
Fail alternates all enqueue the next fail alternate on failure, or the next pass alternate on success, while pass alternates enqueue the next pass alternate on either success or failure; the final fail alternate marks the failure of the nonterminal on failure, while any alternate that succeeds marks the success of the nonterminal.
The effect of this duplication of alternates is to encode the failure or success of previous unordered alternates in the control-flow of the algorithm, rather than modifying the FailCRF to store this information in the data structure.
This choice was made to avoid a memory-usage penalty in all FailCRF nodes to support unordered choice, but empirical investigation of the relative trade-offs of both approaches would be a fruitful direction for future work.

\begin{figure}
	\caption[Unordered choice code pattern]{Pattern for an unordered-choice nonterminal $X ::= \tau_1 | \ldots | \tau_p, p \geq 1$} \label{unordered-choice-algo}
	\begin{algorithmic}
	\Case{$X_{fail} ::= \cdot \tau_1$}
		\State $\langle$ generate code for $X ::= \cdot \tau_1$ with failure path $X_{fail} ::= \cdot \tau_2$ $\rangle$
		\State \Call{rtn}{$X$, $c_U$, $c_I$};
		\State $(L, c_I, t) \gets (X_{pass} ::= \cdot \tau_2, c_U, t')$ \Goto{nextSlot}
	\EndCase
	\Case{$X_{pass} ::= \cdot \tau_2$}
		\State $\langle$ generate code for $X ::= \cdot \tau_2$ with failure path $X_{pass} ::= \cdot \tau_3$ $\rangle$
		\State \Call{rtn}{$X$, $c_U$, $c_I$};
		\State $(L, c_I, t) \gets (X_{pass} ::= \cdot \tau_3, c_U, t')$ \Goto{nextSlot}
	\EndCase
	\State $\vdots$
	\State ~
	\Case{$X_{fail} ::= \cdot \tau_p$}
		\State $\langle$ generate code for $X ::= \cdot \tau_p$ with failure path $X_{fail} ::= \cdot \varnothing$ $\rangle$
		\State \Call{rtn}{$X$, $c_U$, $c_I$} \Goto{nextDesc}
	\EndCase
	\Case{$X_{fail} ::= \cdot \varnothing$}
		\State \Call{rtn}{$X$, $c_U$, $\fail$} \Goto{nextDesc}
	\EndCase
	\end{algorithmic}
	\end{figure}

\section{Related Work}
There are a number of pre-existing approaches to PEG parsing; Ford~\cite{For02} introduced the PEG formalism and two algorithms, a backtracking \emph{recursive descent} algorithm and a memoized \emph{packrat} algorithm.
In practice, recursive decent is more efficient, but has exponential worst-case runtime, while packrat trades increased space and time in the common case for linear worst-case runtime\cite{Mos17}; most subsequent work has attempted to mitigate one or both of these shortcomings.
Mizushima~\etal~\cite{MMY10} use \emph{cut operators} to reduce the memory usage of packrat parsing, while Kuramitsu~\cite{Kur15} and Redziejowski~\cite{Red07} use heuristic table-trimming mechanisms to similar effect.
Medeiros \& Ierusalimschy~\cite{MI08} and Henglein \& Rasmussen~\cite{HR17} have developed a parsing machine approach and a tabular parsing algorithm, respectively, both of which have some evidence of efficient performance.
Moss~\cite{Mos20} and Garnock-Jones~\etal~\cite{GJWE18} have independently developed derivative parsing algorithms for PEGs, neither of which claim improved performance, but which allow use of PEGs with algorithmic tools based on derivative parsing.
Similarly to this work, Chida \& Kuramitsu \cite{CK17} have extended the PEG formalism to include unordered choice.
Medeiros~\etal \cite{MMI14left} and Hutchinson \cite{Hut20} have developed PEG parsers that accept left-recursive grammars, while Redziejowski \cite{Red21} has developed a grammar re-writing scheme that accomplishes the same end.

Scott \& Johnstone's GLL algorithm \cite{SJ10,SJ16} is based on the GLR algorithm of Tomita\cite{Tom85}.
The PEGLL implementation (based on Ackerman's GoGLL \cite{Ack19}) also uses a modified call-return forest (CRF) for control flow and a binary subtree representation (BSR) for parse results, both as in Scott~\etal \cite{SJvB19}.
Afroozeh \& Izmaylova \cite{AI15} have contributed performance improvements and extensions to the original GLL algorithm.
Parr~\etal \cite{PHF14} have developed the \emph{ALL(*)} parsing algorithm, with similar power to GLL but better reported runtime performance.

\section{Conclusion \& Future Work}

This paper has described PEGLL, a new algorithm for parsing parsing expression grammars, based on Scott \& Johnstone's GLL \cite{SJ10,SJ16}.
The primary contribution is the FailCRF data structure, which adds failure paths to GLL's call-return forest, and code-generation algorithms to implement PEG semantics on top of this data structure.
In addition to supporting the fundamental ordered-choice and lookahead operators of PEG, PEGLL also includes support for repetition operators as syntactic sugar, and also for CFG-style unordered choice using a duplication approach to track whether any of the unordered alternates have matched.

An Apache-licenced implementation of PEGLL is available on GitHub\footnote{\url{https://github.com/bruceiv/pegll}}; this implementation passes the research team's internal correctness tests, but performance testing is left to future work.
In addition to measuring and optimizing the performance of the PEGLL algorithm, a straightforward piece of future work would be to expand the expressivity of PEGLL's operator set; parenthesized subexpressions would be a natural expansion, as would more variants on the repetition operators -- in particular, non-greedy variants could be implemented by using unordered choice instead of ordered choice in the generated nonterminals.

A more ambitious advance would be to support left-recursive rules in PEGLL.
Combined with PEGLL's existing support for unordered choice, allowing left-recursive rules would make the set of languages parsed by PEGLL a superset of both parsing-expression languages and context-free languages, allowing language designers to apply the idioms and capabilities of both formalisms as desired.
Given that GLL already supports left-recursive rules \cite{SJ10} and a number of approaches have been suggested for PEG-parsing as well \cite{MMI14left,Hut20,Red21}, this goal should be achievable.

\bibliographystyle{ieeetr}