# Multiway Ensemble Kalman Filter

Yu Wang
University of Michigan
<EMAIL_ADDRESS>
&Alfred Hero
University of Michigan
<EMAIL_ADDRESS>

###### Abstract

In this work, we study the emergence of sparsity and multiway structures in
second-order statistical characterizations of dynamical processes governed by
partial differential equations (PDEs). We consider several state-of-the-art
multiway covariance and inverse covariance (precision) matrix estimators and
examine their pros and cons in terms of accuracy and interpretability in the
context of physics-driven forecasting when incorporated into the ensemble
Kalman filter (EnKF). In particular, we show that multiway data generated from
the Poisson and the convection-diffusion types of PDEs can be accurately
tracked via EnKF when integrated with appropriate covariance and precision
matrix estimators.

##  1 Introduction

There has recently been a resurgence of interest in integrating machine
learning with physics-based modeling. Much of the recent work has focused on
black-box models such as deep neural networks [20, 14, 30, 24, 16, 25].
However, seeking shallower models that capture mechanism in a physically
interpretable manner has been a recurring theme in both machine learning and
physics [28]. In this paper, we introduce a high-dimensional statistical
approach that naturally integrates physics and machine learning through
Kronecker-structured Gaussian graphical models. The learned representation can
then be incorporated into a high dimensional predictive model using the
ensemble Kalman filtering framework.

Multiway covariance/precision models. High-dimensional multiway/tensor-variate
data arise naturally in physical sciences. For example, weather satellites
measure spatio-temporal climate variables such as temperature, wind velocity,
sea level, pressure, etc. Due to the non-homogeneous nature of these data,
second-order information that encodes (conditional) dependency structure
within the data is of interest. Assuming the data are drawn from a tensor
normal distribution, a straightforward way to estimate this structure is to
vectorize the tensor and estimate the underlying Gaussian graphical model
associated with the vector. However, such an approach ignores the tensor
structure and requires estimating a rather high dimensional precision matrix,
often with insufficient sample size. In many scientific applications the
sample size can be as small as one when only a single tensor-valued
measurement is available.

To address sample complexity challenges in learning second-order
representations for multiway (tensor) data, sparsity is often imposed on the
covariance $\mathbf{\Sigma}$ or the inverse covariance $\mathbf{\Omega}$. Such
approaches include the sparse Kronecker product (KP) or Kronecker sum (KS)
decomposition of $\mathbf{\Sigma}$ or $\mathbf{\Omega}$. Statistical models
and corresponding learning algorithms can be derived using generative models
or matrix approximations. The former include: KGlasso/Tlasso [23, 15] for
estimating $\mathbf{\Omega}=\mathbf{A}\otimes\mathbf{B}$, using a
representation $\mathbf{A}\mathbf{X}\mathbf{B}=\mathbf{Z}$ for data
$\mathbf{X}$ when $\mathbf{Z}$ is white noise. Another generative model is
SyGlasso/SG-PALM [27, 26] that models the precision matrix as
$\mathbf{\Omega}=(\mathbf{A}\oplus\mathbf{B})^{2}$, which corresponds to
assuming the data $\mathbf{X}$ obeys a Sylvester equation
$\mathbf{X}\mathbf{A}+\mathbf{B}\mathbf{X}=\mathbf{Z}$. Matrix approximation
methods include: KPCA [22, 8] that approximates the covariance matrix as
$\mathbf{\Sigma}=\sum_{i=1}^{r}\mathbf{A}_{i}\otimes\mathbf{B}_{i}$. Another
matrix approximation method is the TeraLasso [9] that models the precision
matrix as $\mathbf{\Omega}=\mathbf{A}\oplus\mathbf{B}$. TeraLasso is
equivalent to approximation of the conditional dependency graph (encoded by
the precision matrix) with a Cartesian product of smaller graphs 111Note that
Tlasso, TeraLasso, Syglasso/SG-PALM are generalizable to precision matrices of
the form $\bigotimes_{k=1}^{K}\mathbf{A}_{k}$,
$\bigoplus_{k=1}^{K}\mathbf{A}_{k}$, and
$(\bigoplus_{k=1}^{K}\mathbf{A}_{k})^{2}$, respectively, for $K\geq 2$..

Multiway second-order characterization of dynamic processes. Physical systems
often exhibits sparsity and low-rank structures in their covariance or inverse
covariance matrix. This is due to the fact that many physical systems are
governed by differential equations, which are characterized by sparse
differential operators. For instance, Wang and Hero [26] showed that for
multivariate data generated by the Poisson equation, the discretized data has
an inverse covariance matrix equal to a squared Kronecker sum of smaller
sparse matrices. In related work, Lindgren et al. [13] elucidated a link
between certain classes of Gaussian Fields (GF) and Gaussian Markov Random
Fields (GMRF) via stochastic partial differential equations, and showed that
efficient learning algorithms can be developed using the fact that GMRFs have
sparse precision matrices.

For multiway/tensor-variate Gaussian data, the aforementioned multiway
(inverse) covariance estimators have been shown to be statistically consistent
in high-dimensional regimes when sample sizes ($N$) are much less than the
dimensionality ($d$) of the covariates. An important question is whether these
Kronecker structures can be integrated into the Kalman filter for tracking the
states of a physical system that generates multiway data.

##  2 Numerical experiments: ensemble Kalman filtering

The Kalman filter is a well-known technique to track the states of a linear
system over time, and many variants have been proposed to deal with non-linear
systems, such as the extended and ensemble Kalman filters. The ensemble Kalman
filter (EnKF) is particularly effective when the dynamical system is
complicated and its gradient is infeasible to calculate, which is often the
case in physical systems [5, 1]. However, such systems are often high-
dimensional and the EnKF operates in the regime where the number of ensemble
members, $N$, is much less than the size of the state, $d$, suggesting that
sparse inverse covariance models will be especially attractive. Hou et al.
[11] introduced a sparsity-penalized EnKF, which uses an estimator of the
forecasting covariance whose inverse is sparsity regularized. Here we propose
incorporating the multiway covariance / inverse covariance models discussed in
Section 1 into the EnKF of Hou et al. [11].

For motivation, we consider the Poisson equation, an elliptical PDE that
governs many physical processes including electromagnetic induction, heat
transfer, and convection [2]. On a rectangular region
$\Omega=(0,d_{1})\times(0,d_{2})$ in the 2D Cartesian plane, the Poisson
equation with homogeneous Dirichlet boundary condition is expressed as

$\displaystyle\mathcal{D}u=(\partial^{2}_{x}+\partial^{2}_{y})u$
$\displaystyle=f\quad\text{in }\Omega,$ (1) $\displaystyle u$
$\displaystyle=0\quad\text{on }\partial\Omega$

where $f:\Omega\to{\mathbb{R}}$ is the given source function and
$u:\Omega\to{\mathbb{R}}$ is the unknown. Using the finite difference method
with a square mesh grid with unit spacing, the unknown and the source can be
expressed as $d_{1}$-by-$d_{2}$ matrices, ${\mathbf{U}}$ and ${\mathbf{F}}$,
respectively, that are related to each other via

$\displaystyle U_{i+1,j}+U_{i-1,j}+U_{i,j+1}+U_{i,j-1}-4U_{i,j}=F_{i,j}$ (2)

for any interior grid point $(i,j)$. Defining $n$-by-$n$ square matrix

$\mathbf{A}_{n}=\begin{bmatrix}2&-1&&\\\ -1&2&\ddots&\\\ &\ddots&\ddots&-1\\\
&&-1&2\end{bmatrix},$

the relation (2) can be expressed as the (vectorized) Sylvester equation with
$K=2$:

$({\mathbf{A}}_{d_{1}}\oplus{\mathbf{A}}_{d_{2}}){\mathbf{u}}={\mathbf{f}},$
(3)

where ${\mathbf{u}}=\operatorname{vec}({\mathbf{U}})$,
${\mathbf{f}}=\operatorname{vec}({\mathbf{F}})$. Note that $\mathbf{A}$ is
tridiagonal. In the case where $\mathbf{f}$ is white noise with variance
$\sigma^{2}$, the inverse covariance matrix of ${\mathbf{u}}$ has the form
${}\mathrm{cov}^{-1}({\mathbf{u}})=\sigma^{-2}({\mathbf{A}}_{d_{1}}\oplus{\mathbf{A}}_{d_{2}})^{T}({\mathbf{A}}_{d_{1}}\oplus{\mathbf{A}}_{d_{2}})$
and hence sparse.

Here, we discuss two ways to extend the spatial Poisson equation described
above to incorporate temporal dynamics, and illustrate how multiway (inverse)
covariance models can be used to track spatio-temporal systems. The first
extension, which we call the Poisson-AR(1) process, imposes an autoregressive
temporal model of order 1 on the source function $f$ in the Poisson equation
(1). Specifically, we say a sequence of discretized spatial observations
$\\{{\mathbf{U}}^{k}\in{\mathbb{R}}^{d_{1}\times d_{2}}\\}_{k}$ indexed by
time step $k=1,\cdots,T$ is from a Poisson-AR(1) process if

$\displaystyle(\mathbf{A}_{d_{1}}\oplus\mathbf{A}_{d_{2}})\operatorname{vec}(\mathbf{U}^{k})=\operatorname{vec}(\mathbf{Z}^{k}),$
(4)
$\displaystyle\operatorname{vec}(\mathbf{Z}^{k})=a\operatorname{vec}(\mathbf{Z}^{k-1})+\operatorname{vec}(\mathbf{W}^{k}),\quad|a|<1,$
(5)

where $\mathbf{Z}^{0}\sim\mathcal{N}(\mathbf{0},\sigma^{2}_{z}\mathbf{I})$ and
$\\{\mathbf{W}^{k}\in{\mathbb{R}}^{d_{1}\times d_{2}}\\}_{k}$ is
spatiotemporal white noise, i.e.,
$W_{i,j}^{k}\sim\mathcal{N}(0,\sigma^{2}_{w})$, i.i.d.

The second time-varying extension of the Poisson PDE model (1) is the
convection-diffusion process [2]

$\frac{\partial u}{\partial
t}=\theta\sum_{i=1}^{2}\frac{\partial^{2}u}{\partial
x_{i}^{2}}-\epsilon\sum_{i=1}^{2}\frac{\partial u}{\partial x_{i}}.$ (6)

Here, $\theta>0$ is the diffusivity; and $\epsilon\in{\mathbb{R}}$ is the
convection velocity of the quantity along each coordinate. Note that for
simplicity of discussion here, we assume these coefficients do not change with
space and time (see, Stocker [19], for example, for a detailed discussion).
These equations are closely related to the Navier-Stokes equation commonly
used in stochastic modeling for weather and climate prediction [2, 19].
Coupled with Maxwell’s equations, these equations can be used to model
magneto-hydrodynamics [17], which characterize solar activities including
flares.

A solution of Equation (6) can be approximated similarly as in the Poisson
equation case, through a finite difference approach. Denote the discrete
spatial samples of $u(\mathbf{x},t)$ at time $t_{k}$ as a matrix
$\mathbf{U}^{k}\in{\mathbb{R}}^{d_{1}\times d_{2}}$. We obtain a discretized
update propagating $u(\mathbf{x},t)$ in space and time, which locally
satisfies

$\displaystyle\frac{U_{i,j}^{k}-U_{i,j}^{k-1}}{\Delta t}=$ $\displaystyle\
\theta\left(\frac{U_{i+1,j}^{k}+U_{i-1,j}^{k}+U_{i,j+1}^{k}+U_{i,j-1}^{k}-4U_{i,j}^{k}}{h^{2}}\right)$
(7)
$\displaystyle-\epsilon\left(\frac{U_{i+1,j}^{k}-U_{i-1,j}^{k}+U_{i,j+1}^{k}-U_{i,j-1}^{k}}{2h}\right),$

where $\Delta t=t_{k+1}-t_{k}$ is the time step and $h$ is the mesh step
(spatial grid spacing). For $\Delta t=1$ and $h=1$, $\mathbf{U}^{k}$ can be
shown to obey the Sylvester matrix update equation [21]
$\mathbf{A}_{\theta,\epsilon}\mathbf{U}^{k}+\mathbf{U}^{k}\mathbf{B}_{\theta,\epsilon}^{T}=\mathbf{U}^{k-1}$,
or equivalently,

$\displaystyle(\mathbf{B}_{\theta,\epsilon}\oplus\mathbf{A}_{\theta,\epsilon})\operatorname{vec}({\mathbf{U}}^{k})=\operatorname{vec}({\mathbf{U}}^{k-1}),$
(8)

where $\mathbf{A}_{\theta,\epsilon}\in{\mathbb{R}}^{d_{1}\times d_{1}}$ and
$\mathbf{B}_{\theta,\epsilon}\in{\mathbb{R}}^{d_{2}\times d_{2}}$ are
symmetric tridiagonal matrices whose entries depend on $\theta,\epsilon$, and
$h$ [7]. Note that, although $\\{\mathbf{U}^{k}\\}_{k}$ updates according to a
sparse Kronecker sum, the inverse covariance of the marginal distribution that
smooths over time is not sparse.

We turn Equation (8) into a state space model for the observations
$\mathbf{X}$ by adding i.i.d. state noise
$\mathbf{W}^{k}\sim\mathcal{N}(0,\sigma^{2}_{w}\mathbf{I})$, i.i.d.
observation noise $\mathbf{V}^{k}\sim\mathcal{N}(0,\sigma^{2}_{v}\mathbf{I})$
and a state measurement matrix $\mathbf{H}$:

$\displaystyle\operatorname{vec}(\mathbf{X}^{k})$
$\displaystyle=\mathbf{H}\operatorname{vec}(\mathbf{U}^{k})+\operatorname{vec}(\mathbf{V}^{k}),$
(9) $\displaystyle\operatorname{vec}(\mathbf{U}^{k})$
$\displaystyle=(\mathbf{B}_{\theta,\epsilon}\oplus\mathbf{A}_{\theta,\epsilon})^{-1}\operatorname{vec}(\mathbf{U}^{k-1})+\operatorname{vec}(\mathbf{W}^{k}).$
(10)

We model the observed process $\mathbf{X}$ as an incomplete noisy version of
the convection-diffusion state $\mathbf{U}_{t}$ obeying the Sylvester matrix
update equation above discretization the convection-diffusion becomes assume a
linear Gaussian state-space model for the observed process $\mathbf{X}_{t}$
governed by convection-diffusion dynamics:

$\displaystyle\mathbf{A}_{\theta,\epsilon}\mathbf{U}_{t}+\mathbf{U}_{t}\mathbf{A}_{\theta,\epsilon}=\mathbf{U}_{t-1},$
$\displaystyle\mathbf{X}_{t}=\mathbf{U}_{t}+\mathbf{V}_{t},$

where $\mathbf{V}_{t}\sim\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I})$ is an
i.i.d. Gaussian white noise. Note that in general one might not fully observe
the state, leading to a partially observed measurement process

$\operatorname{vec}(\mathbf{X}_{t})=\mathbf{H}\operatorname{vec}(\mathbf{U}_{t})+\operatorname{vec}(\mathbf{V}_{t}),$

where $\mathbf{H}$ is a measurement matrix that can incorporate effects such
as unobserved, masked, or superposed states. It is also possible to take into
account the error in the dynamic process, i.e.,

$\operatorname{vec}(\mathbf{U}_{t})=(\mathbf{A}_{\theta,\epsilon}\oplus\mathbf{A}_{\theta,\epsilon})^{-1}\operatorname{vec}(\mathbf{U}_{t-1})+\operatorname{vec}(\mathbf{W}_{t},)$

where $\operatorname{vec}(\mathbf{W}_{t})$ is assumed to be white noise. Under
the case of perfect observation ($\mathbf{H}$ is the identity matrix and
$\mathbf{W}_{t}=\mathbf{0}$), Note that although the state variable evolves
via a Sylvester equation, similar to the Poisson equation case, the state
(inverse) covariance matrix at time step $t_{k}$ admits different structures.
Specifically, the state precision matrix
$\mathbf{\Omega}^{k}={}\mathrm{cov}^{-1}(\operatorname{vec}({\mathbf{U}^{k}}))\in{\mathbb{R}}^{d_{1}d_{2}\times
d_{1}d_{2}}$ evolves as
$\mathbf{\Omega}^{k}=(\mathbf{B}_{\theta,\epsilon}\oplus\mathbf{A}_{\theta,\epsilon})\mathbf{\Omega}^{k-1}(\mathbf{B}_{\theta,\epsilon}\oplus\mathbf{A}_{\theta,\epsilon})+\sigma^{-2}_{w}\mathbf{I}$
(see Katzfuss et al. [12], for example). This matrix is not necessarily sparse
for finite $k$ but, assuming that the eigenvalues of the matrix
$\mathbf{B}_{\theta,\epsilon}\oplus\mathbf{A}_{\theta,\epsilon}$ are in
$(-1,1)$, the limiting precision matrix
$\mathbf{\Omega}^{\infty}=\lim_{k\rightarrow\infty}\mathbf{\Omega}^{k}$ is
$\mathbf{\Omega}^{\infty}=(\mathbf{B}_{\theta,\epsilon}\oplus\mathbf{A}_{\theta,\epsilon})\mathbf{\Omega}^{\infty}(\mathbf{B}_{\theta,\epsilon}\oplus\mathbf{A}_{\theta,\epsilon})+\sigma^{-2}\mathbf{I}$.
The $\mathbf{\Omega}^{\infty}$ matrix is sparse because
$\mathbf{A}_{\theta,\epsilon}$ and $\mathbf{B}_{\theta,\epsilon}$ are both
tridiagonal.

To illustrate, we consider a 2D spatio-temporal process of dimension $64\times
64$ where only half of the entries are observed, which leads to a measurement
matrix $\mathbf{H}\in\\{0,1\\}^{2048\times 4069}$. We generated the true
states and the corresponding observations according to Convection-Diffusion
and Poisson dynamics for $T=20$ time steps. Several realizations of the true
state variables are shown in Figure 1. At each time step, we generated an
ensemble of size $N=25$ and estimated the state covariance / inverse
covariance using several sparse (multiway) inverse covariance estimation
methods, including Glasso [6], KPCA [8], KGlasso [23], TeraLasso [9], SG-PALM
[26]. Figure 2 shows evolution of the computed root mean squared errors
(RMSEs) for the estimated states across all ensemble members. In Figure 3 we
show the true and estimated (inverse) covariance matrices obtained for the
last time step. The Poisson process is time-invariant, and at each time step
the EnKF involves estimation of a sparse Kronecker sum squared inverse
covariance matrix. Hence, the SG-PALM method operates under the correct model
assumption in this situation. On the other hand, the inverse covariance
structure under the convection-diffusion dynamics model is dense due to the
smoothing effect of Kalman filtering and the nature of the temporal dynamics.
But, its steady-state covariance has low-dimensional structures. The KPCA in
this case was able to approximate this structure, as illustrated in Figure 3.

Figure 1: 2D convection-diffusion (bottom) and Poisson (top) state variables at three different time stamps. Note that there is temporal correlation exists in the convection-diffusion states while the Poisson states are temporally independent. Noise-free states |  Noisy states
---|---
Convection-Diffusion |
Poisson |
Figure 2: RMSEs of the estimated states via EnKF over $20$ time steps using
different (inverse) covariance estimators. RMSEs over all ensemble members are
shown here with the mean highlighted using solid lines. Here, each state is of
dimension $64\times 64$ and is generated via either a convection-diffusion
(top row) or Poisson equation (bottom row). The best performers in terms of
mean RMSE over all ensemble members are KPCA for convection-diffusion and SG-
PALM for Poisson. Figure 3: Covariance/precision structures for Poisson and
Convection-Diffusion dynamics and their estimates. Here, white/blank entries
indicate zeros in the (inverse) covariance matrix. For Poisson dynamics the
Sylvester graphical model matches the true structure of the precision matrix.
For convection-diffusion dynamics the covariance instead of the precision
matrix is structured and sparse.

Remark. The proposed multiway EnKF is able to track systems governed by
elliptic (e.g., Poisson) parabolic, and hyperbolic (e.g., convection-
diffusion) PDEs. There are indeed other important PDEs / dynamical systems
that cannot be modelled by these types of equations. Furthermore, the
Sylvester matrix equations arise when the finite-difference discretization is
performed on a rectangular grid. The relations (2) and (7) might not hold for
finite-difference on, for example, spherical coordinates.

##  3 Future work

### Applications.

Spatiotemporal PDEs are prominent techniques for modeling real-world physical
systems. One such system arises in space physics, where solar flares and
coronal mass ejections are associated with rapid changes in filed connectivity
and are powered by partial dissipation of electrical currents in the solar
atmosphere [18]. The nonlinear force-free filed model is often used to
describe the solar coronal magnetic field [4, 29] and can be derived from the
convection-diffusion process described in this work. Additionally, global maps
of the solar photospheric magnetic flux are fundamental drivers for
simulations of the corona and solar wind. However, observations of the solar
photosphere are only made intermittently over approximately half of the solar
surface. Hickmann et al. [10] introduced the Air Force Data Assimilative
Photospheric Flux Transport model that uses localized ensemble Kalman
filtering to adjust a set of photospheric simulations to agree with the
available observations. In future work, we plan to incorporate our proposed
multiway EnKF framework for tracking these solar physical systems.

### Interpretability of factorization-based multiway model.

While the Kronecker products expansion used in KPCA captures dense structures
in the covariance matrix of data generated from more complex spatio-temporal
physical processes, it lacks physical interpretability. In contrast to the
case of Sylvester graphical model and Poisson processes, it is not obvious
whether the sum of Kronecker products structure corresponds to any true
physical models. Recent work in quantum informatics [3] has demonstrated a
link between estimation of the density matrix for entangled quantum states and
the structured tensor approximation via
$\sum_{i=1}^{r}\mathbf{A}_{i}\otimes\mathbf{B}_{i}$. Further characterizing
these connections and extending them to study its connections with certain
classes of discretized PDEs would be an interesting future direction.

## Acknowledgments and Disclosure of Funding

The authors thank Zeyu Sun at the University of Michigan for many helpful
discussions of potential applications of the work on solar physics. The work
is partially supported by NASA DRIVE Science Center grant 80NSSC20K0600 and by
the National Nuclear Security Administration within the US Department of
Energy under grant DE-NA0003921.

## References

  * Burgers et al. [1998] Gerrit Burgers, Peter Jan van Leeuwen, and Geir Evensen.  Analysis scheme in the ensemble kalman filter.  _Monthly weather review_ , 126(6):1719–1724, 1998\.
  * Chandrasekhar [1943] Subrahmanyan Chandrasekhar.  Stochastic problems in physics and astronomy.  _Reviews of modern physics_ , 15(1):1, 1943.
  * Chu and Lin [2021] Moody T Chu and Matthew M Lin.  Nonlinear power-like and svd-like iterative schemes with applications to entangled bipartite rank-1 approximation.  _SIAM Journal on Scientific Computing_ , (0):S448–S474, 2021.
  * DeRosa et al. [2015] ML DeRosa, MS Wheatland, KD Leka, Graham Barnes, Tahar Amari, Aurélien Canou, SA Gilchrist, JK Thalmann, Gherardo Valori, Thomas Wiegelmann, et al.  The influence of spatial resolution on nonlinear force-free modeling.  _The Astrophysical Journal_ , 811(2):107, 2015\.
  * Evensen [1994] Geir Evensen.  Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics.  _Journal of Geophysical Research: Oceans_ , 99(C5):10143–10162, 1994.
  * Friedman et al. [2008] Jerome Friedman, Trevor Hastie, and Robert Tibshirani.  Sparse inverse covariance estimation with the graphical lasso.  _Biostatistics_ , 9(3):432–441, 2008.
  * Grasedyck [2004] Lars Grasedyck.  Existence and computation of low kronecker-rank approximations for large linear systems of tensor product structure.  _Computing_ , 72(3-4):247–265, 2004.
  * Greenewald and Hero [2015] Kristjan Greenewald and Alfred O Hero.  Robust kronecker product pca for spatio-temporal covariance estimation.  _IEEE Transactions on Signal Processing_ , 63(23):6368–6378, 2015.
  * Greenewald et al. [2019] Kristjan Greenewald, Shuheng Zhou, and Alfred Hero III.  Tensor graphical lasso (teralasso).  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 81(5):901–931, 2019.
  * Hickmann et al. [2015] Kyle S Hickmann, Humberto C Godinez, Carl J Henney, and C Nick Arge.  Data assimilation in the adapt photospheric flux transport model.  _Solar Physics_ , 290(4):1105–1118, 2015.
  * Hou et al. [2021] Elizabeth Hou, Earl Lawrence, and Alfred O Hero.  Penalized ensemble kalman filters for high dimensional non-linear systems.  _PloS one_ , 16(3):e0248046, 2021.
  * Katzfuss et al. [2016] Matthias Katzfuss, Jonathan R Stroud, and Christopher K Wikle.  Understanding the ensemble kalman filter.  _The American Statistician_ , 70(4):350–357, 2016\.
  * Lindgren et al. [2011] Finn Lindgren, Håvard Rue, and Johan Lindström.  An explicit link between gaussian fields and gaussian markov random fields: the stochastic partial differential equation approach.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 73(4):423–498, 2011.
  * Long et al. [2018] Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong.  Pde-net: Learning pdes from data.  In _International Conference on Machine Learning_ , pages 3208–3216. PMLR, 2018.
  * Lyu et al. [2019] Xiang Lyu, Will Wei Sun, Zhaoran Wang, Han Liu, Jian Yang, and Guang Cheng.  Tensor graphical model: Non-convex optimization and statistical inference.  _IEEE transactions on pattern analysis and machine intelligence_ , 42(8):2024–2037, 2019.
  * Reichstein et al. [2019] Markus Reichstein, Gustau Camps-Valls, Bjorn Stevens, Martin Jung, Joachim Denzler, Nuno Carvalhais, et al.  Deep learning and process understanding for data-driven earth system science.  _Nature_ , 566(7743):195–204, 2019.
  * Roberts [2006] B Roberts.  Slow magnetohydrodynamic waves in the solar atmosphere.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 364(1839):447–460, 2006.
  * Schrijver et al. [2008] CJ Schrijver, ML DeRosa, Thomas Metcalf, Graham Barnes, Bruce Lites, T Tarbell, James McTiernan, Gherardo Valori, Thomas Wiegelmann, MS Wheatland, et al.  Nonlinear force-free field modeling of a solar active region around the time of a major flare and coronal mass ejection.  _The Astrophysical Journal_ , 675(2):1637, 2008\.
  * Stocker [2011] Thomas Stocker.  _Introduction to climate modelling_.  Springer Science & Business Media, 2011.
  * Takeishi et al. [2017] Naoya Takeishi, Yoshinobu Kawahara, and Takehisa Yairi.  Learning koopman invariant subspaces for dynamic mode decomposition.  _arXiv preprint arXiv:1710.04340_ , 2017.
  * Thomas [2013] James William Thomas.  _Numerical partial differential equations: finite difference methods_ , volume 22.  Springer Science & Business Media, 2013.
  * Tsiligkaridis and Hero [2013] Theodoros Tsiligkaridis and Alfred O Hero.  Covariance estimation in high dimensions via kronecker product expansions.  _IEEE Transactions on Signal Processing_ , 61(21):5347–5360, 2013.
  * Tsiligkaridis et al. [2013] Theodoros Tsiligkaridis, Alfred O Hero III, and Shuheng Zhou.  On convergence of kronecker graphical lasso algorithms.  _IEEE transactions on signal processing_ , 61(7):1743–1755, 2013.
  * Vlachas et al. [2018] Pantelis R Vlachas, Wonmin Byeon, Zhong Y Wan, Themistoklis P Sapsis, and Petros Koumoutsakos.  Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 474(2213):20170844, 2018.
  * Wang et al. [2020a] Rui Wang, Karthik Kashinath, Mustafa Mustafa, Adrian Albert, and Rose Yu.  Towards physics-informed deep learning for turbulent flow prediction.  In _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, pages 1457–1466, 2020a.
  * Wang and Hero [2021] Yu Wang and Alfred Hero.  Sg-palm: a fast physically interpretable tensor graphical model.  _International Conference on Machine Learning (ICML), arXiv preprint arXiv:2105.12271_ , 2021.
  * Wang et al. [2020b] Yu Wang, Byoungwook Jang, and Alfred Hero.  The sylvester graphical lasso (syglasso).  In _International Conference on Artificial Intelligence and Statistics_ , pages 1943–1953. PMLR, 2020b.
  * Weinan et al. [2020] E Weinan, Jiequn Han, and Linfeng Zhang.  Integrating machine learning with physics-based modeling.  _arXiv_ , 2020.
  * Wheatland and Gilchrist [2013] MS Wheatland and SA Gilchrist.  The state of nonlinear force-free magnetic field extrapolation.  In _Journal of Physics: Conference Series_ , volume 440, page 012037\. IOP Publishing, 2013.
  * Zhang et al. [2018] Linfeng Zhang, Jiequn Han, Han Wang, Wissam Saidi, Roberto Car, and Weinan E.  End-to-end symmetry preserving inter-atomic potential energy model for finite and extended systems.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 31. Curran Associates, Inc., 2018.  URL https://proceedings.neurips.cc/paper/2018/file/e2ad76f2326fbc6b56a45a56c59fafdb-Paper.pdf.

## Checklist

  1. 1.

For all authors…

    1. (a)

Do the main claims made in the abstract and introduction accurately reflect
the paper’s contributions and scope? [Yes]

    2. (b)

Did you describe the limitations of your work? [Yes] Both at the end of the
Section 2 and the future work discussed in Section 3.

    3. (c)

Did you discuss any potential negative societal impacts of your work? [N/A]

    4. (d)

Have you read the ethics review guidelines and ensured that your paper
conforms to them? [Yes]

  2. 2.

If you are including theoretical results…

    1. (a)

Did you state the full set of assumptions of all theoretical results? [N/A]

    2. (b)

Did you include complete proofs of all theoretical results? [N/A]

  3. 3.

If you ran experiments…

    1. (a)

Did you include the code, data, and instructions needed to reproduce the main
experimental results (either in the supplemental material or as a URL)? [Yes]
The code is included as part of a software package, available at:
https://github.com/ywa136/TensorGraphicalModels.jl

    2. (b)

Did you specify all the training details (e.g., data splits, hyperparameters,
how they were chosen)? [Yes]

    3. (c)

Did you report error bars (e.g., with respect to the random seed after running
experiments multiple times)? [Yes]

    4. (d)

Did you include the total amount of compute and the type of resources used
(e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] All numerical
experiments were run on a laptop equipped with Intel Core i5 Quad-Core 1.4 GHz
with a 16 GB RAM.

  4. 4.

If you are using existing assets (e.g., code, data, models) or
curating/releasing new assets…

    1. (a)

If your work uses existing assets, did you cite the creators? [N/A]

    2. (b)

Did you mention the license of the assets? [N/A]

    3. (c)

Did you include any new assets either in the supplemental material or as a
URL? [N/A]

    4. (d)

Did you discuss whether and how consent was obtained from people whose data
you’re using/curating? [N/A]

    5. (e)

Did you discuss whether the data you are using/curating contains personally
identifiable information or offensive content? [N/A]

  5. 5.

If you used crowdsourcing or conducted research with human subjects…

    1. (a)

Did you include the full text of instructions given to participants and
screenshots, if applicable? [N/A]

    2. (b)

Did you describe any potential participant risks, with links to Institutional
Review Board (IRB) approvals, if applicable? [N/A]

    3. (c)

Did you include the estimated hourly wage paid to participants and the total
amount spent on participant compensation? [N/A]


# Data Augmentation for Automated Essay Scoring using Transformer Models

Kshitij Gupta  Department of Electrical and Electronics Engineering
BITS Pilani, Pilani Campus
Pilani, India
<EMAIL_ADDRESS>
###### Abstract

Automated essay scoring is one of the most important problem in Natural
Language Processing. It has been explored for a number of years, and it
remains partially solved. In addition to its economic and educational
usefulness, it presents research problems. Transfer learning has proved to be
beneficial in NLP. Data augmentation techniques have also helped build state-
of-the-art models for automated essay scoring. Many works in the past have
attempted to solve this problem by using RNNs, LSTMs, etc. This work examines
the transformer models like BERT, RoBERTa, etc. We empirically demonstrate the
effectiveness of transformer models and data augmentation for automated essay
grading across many topics using a single model.

###### Index Terms:

Automated System, Transformers, BERT

##  I Introduction

As a result of the COVID-19 pandemic, online schooling system became
necessary. From elementary schools to colleges, almost all educational
institutions have adopted the online education system. The majority of
automated evaluations are accessible for multiple-choice questions, but
evaluating short and essay type responses remains unsolved since, unlike
multiple-choice questions, there is no one correct solution for these kind of
questions. It is an essential education-related application that employs NLP
and machine learning methodologies. It is difficult to evaluate essays using
basic computer languages and methods such as pattern matching and language
processing.

Among the most important pedagogical uses of NLP is automated essay scoring
(AES), the technique of using a system to score short and essay type questions
without manual assistance. Initiated by Page’s [1966] groundbreaking work on
the Project Essay Grader system, this area of study has seen continuous
activity ever since. The bulk of AES research has been on holistic scoring,
which provides a quantitative summary of an essay’s quality in a single
number. At least two factors contribute to this concentration of effort. To
begin with, learning-based holistic scoring systems may make use of publically
accessible corpora that have been manually annotated with holistic scores.
Second, there is a market for holistic scoring algorithms because they may
streamline the arduous process of manually evaluating the millions of essays
for tests like GRE, IELTS, SAT.

Past research on automated essay grading has included training models for
essays for which training data is available and those models are topic
specific. This model is trained on all the topics thus could be used for
assessment of essays of all those topics without training model specific for
each topic. This would be useful in the scenario where we did not have enough
data to train a model that is specific to a particular topic, but we still
needed to evaluate essays on that topic. Therefore, in order to assess them,
We may utilize a model that has been trained on essays on a variety of topics
and a tiny amount of data on the topic for which we need to develop a model,
which will then be fine-tuned using the limited data available on the subject
being assessed.

This paper is organized as follows: In Section II, we explore pertinent prior
research on automated essay scoring; in Section III, we cover experimental
setup; and in Section IV, we describe our methodology for augmenting essay
data. In Section V, we give the results and analysis of the automated essay
grading model. Section VI comprises of conclusion and future work for
Automated Essay Scoring.

##  II Related Works

Project Essay Grader (PEG) by [1] started the research on Automated essay
scoring. Shermis (2001)[2] improved the PEG system by incorporating the
grammatical features as well in the evaluation. Around the turn of century,
great majority of essay scoring systems used conventional methods like latent
semantic analysis by Foltz (1999)[3], as pattern matching and statistical
analysis like Bayesian Essay Test Scoring System by [4]. These systems employ
natural language processing (NLP) approaches that concentrate on grammar,
content to determine an essay’s score.

Multiple studies studied AES systems, from the earliest to the most recent.
Blood (2011)[6] reviewed the PEG literature from 1984 to 2010, it has
discussed just broad features of AES systems, such as ethical considerations
and system performance. However, they have not addressed the implementation
aspect, nor has a comparison research been conducted, nor have the real
problems of AES systems been highlighted.

After 2014, Automated grading systems like as those by [5] and others,
employed deep learning approaches to induce syntactic and semantic
characteristics, producing greater outcomes than previous systems. Burrows
(2015)[7] reviewed on aspects, including datasets, NLP approaches, model
construction, model grading, model assessment, and model efficacy. Ke
(2019)[8], Hussein (2019)[9] and Klebanov (2020)[10] offered overview of the
AES system.

Ramesh (2019)[19] offered us a comprehensive summary of all accessible
datasets and the machine learning algorithms used to grade essays. Recently
Park (2022)[23] tried using GAN’s for evaluation of essays. Our motivation of
using Deep-learning models for automated essay scoring was due to recent study
by Elijah (2020)[20] and Wang (2022)[25], it showed us that using BERT based
models it certainly improves the accuracy of the Automated essay scoring.

Jong (2022)[24] demonstrated the effectiveness of data augmentation techniques
on automated essay scoring. Ludwig (2021)[21], Ormerod (2021)[22] and Sethi
(2022)[26] examined the transformer based models on automated essay scoring,
these recent findings prompted us to test our data augmentation strategy on
transformer-based models for automated essay grading.

##  III Experimental Setup

For the purpose of this study, we will constrain the study and experiments to
ASPA1 dataset. The Automated Student Assessment
Prize(ASAP1)111https://www.kaggle.com/c/asap-aes corpus was released as part
of a Kaggle competition in 2012. This corpus is used for holistic scoring of
essays and consists of essays from eight topics with around 3000 essays on
each of the topic. Each topic had a different scoring method, so, we
normalized of each essay score from 0 to 10 so that we could train all the
data together.

We run our experiments using the BERT, RoBERTa, ALBERT, DistilBERT, XLM-
RoBERTa, implementation available in the Simple Transformers library provided
by Thilina Rajapakse222https://simpletransformers.ai/. We run our experiment
on Google Colab. For baseline trials, we consider the essays which are not
supplemented with our data augmentation technique and trained on models based
on Transformers.

During training, we additionally fine-tuned our parameters by adjusting the
parameters like learning rate, weight decay rate, etc. For the purpose of
evaluating our model, we will use the accuracy score measure given by the
Scikit library333https://scikit-learn.org/stable/index.html. Since we are
approaching our issue as a multi-label classification, the accuracy metric is
often used to evaluate multi-label classifications.

Figure 1: Transformers Architecture

##  IV Methodology

###  IV-A Large Pre-trained models

The models which we used for training are based on the Transformers(Fig. 1)
architecture introduced by [11]. The Transformer’s architecture follows an
encoder-decoder structure.

Given below is the brief description of each of the model which we are using
for training:

####  IV-A1 BERT

The Bidirectional Encoder Representations (BERT) introduced by [12] is a deep
learning model in which every output element is linked to every input element
and the weightings between them are dynamically determined depending on their
relationship. BERT is pre-trained on two tasks: Masked Language Modeling and
Next Sentence Prediction.

####  IV-A2 RoBERTa

Robustly optimized BERT Pre-training Approach (RoBERTa) introduced by [13]
builds upon BERT’s language masking method, in which the system learns to
anticipate purposely masked bits of text inside unannotated language samples.
RoBERTa changes critical hyperparameters in BERT, such as eliminating BERT’s
next-sentence pretraining target and training with much bigger mini-batches
and learning rates. This enables RoBERTa to outperform BERT at the masked
language modeling goal and improves the performance of subsequent tasks.

####  IV-A3 ALBERT

ALBERT, introduced by [14] is an encoder-decoder based model with self-
attention at the encoder and attention to encoder outputs at the decoder end.
It is a modified version of BERT and it stands for ”A Lite BERT”. It builds on
parameter sharing, embedding factorization, and sentence order
prediction(SOP).

####  IV-A4 DistilBERT

DistilBERT, introduced by [15] aims to optimize the training by reducing the
size of BERT and increasing the speed of BERT—all while trying to retain as
much performance as possible. Specifically, DistilBERT is smaller than the
original BERT-base model, is faster than it, and retains its functionality.

####  IV-A5 XLM-RoBERTa

XLM-RoBERTa, Unsupervised Cross-lingual Representation Learning at Scale,
introduced by [16] is a scaled cross-lingual sentence encoder. It is trained
on 2.5 TB of data across 100 languages filtered from Common Crawl. XLM-RoBERTa
achieves state-of-the-art results on multiple cross-lingual benchmarks.

Figure 2: Augmentation Technique

###  IV-B Data Augmentation

Researchers have attempted to use several RNNs and LSTMs as training models
for automated essay scoring. However, the fundamental disadvantage of such
models is that they are topic-specific, and we want to construct an automated
essay scoring system that can perform well not just on subjects for which we
have an abundance of data but also on subjects for which we have a limited
amount of data. Now, we want to augment the essay so that it can accurately
assess a essay on a different topic for which we have a very small amount of
data.

When training a model, we add each essay with its topic at certain intervals.
We are including essay topics after an interval because essays are lengthy,
and if we train only on essays, the trained model will be topic-specific. In
order to construct a more robust model, we append the essay’s subject to each
essay so that the training model may learn the relationship between the
essay’s topic and the essay itself. This is so that it accurately grades
essays on a different topic, which we are using majorly for fine-tuning since
we have very little data available to us and we can not train topic-specific
models using this less amount of data.

We cannot add a complete topic to an essay because sometimes essay topics are
quite large. Hence, we summarize the topics of essays using the summarization
pipeline provided by [17] implementation of BART, which was introduced by
[18]. Now, the second issue is, after how many lines should we put the subject
for optimum precision? We conducted extensive data trials by inserting them at
different places.

After a comprehensive investigation of several transformer-based models and
essay subject insertions, we determined that the tenth place is optimal for
inserting the topic. It implies that after every tenth line, a summary of the
subject is added to the essay. It may be due to the fact that the model
struggles with lengthy text classification; thus, we keep this in mind by
inserting topics at regular intervals. Fig. 2 shows the complete data
augmentation of essay. Now that the topic has been inserted into the essay,
the data is modified for training transformer models. The next section
provides a study of modified essays.

TABLE I: Sklearn Accuracy Score of Mutli-Label Classification of all the four
topics in the testing set

TABLE II: Topic 1 Transformer Models |  Unaugmented Data |  Augmented Data
---|---|---
LSTM |  30.8% |  38.3%
XLM-RoBERTa |  47.4% |  59.2%
RoBERTa |  51.1% |  60.4%
ALBERT |  48.7% |  59.7%
DistilBERT |  50.3% |  58.1%
BERT |  50.7% |  60.2%

TABLE III: Topic 2 Transformer Models |  Unaugmented Data |  Augmented Data
---|---|---
LSTM |  29.4% |  37.6%
XLM-RoBERTa |  47.8% |  58.6%
RoBERTa |  50.9% |  60.8%
ALBERT |  49.6% |  60.2%
DistilBERT |  49.8% |  58.8%
BERT |  50.2% |  61.3%

TABLE IV: Topic 3 Transformer Models |  Unaugmented Data |  Augmented Data
---|---|---
LSTM |  30.3% |  38.1%
XLM-RoBERTa |  46.9% |  58.9%
RoBERTa |  51.3% |  61.3%
ALBERT |  49.2% |  59.9%
DistilBERT |  49.4% |  58.9%
BERT |  49.3% |  61.5%

TABLE V: Topic 4 Transformer Models |  Unaugmented Data |  Augmented Data
---|---|---
LSTM |  30.3% |  38.9%
XLM-RoBERTa |  46.9% |  58.4%
RoBERTa |  50.5% |  61.4%
ALBERT |  47.6% |  59.3%
DistilBERT |  51.1% |  58.2%
BERT |  50.3% |  61.7%

##  V Results and Analysis

The ASAP1 dataset contains around 17K essays on eight topics. We are using
that data for pre-training by augmenting those essays using our technique. For
fine tuning For research and testing purposes, we used this dataset. This
dataset consists of 1241 essays on four subjects. We fine-tuned and tested our
models on each subject individually. We used around two-thirds of the above
mentioned dataset for training and fine-tuning, and the remaining one-third
for testing.

We followed a very simple yet state-of-the-art modeling technique for multi-
label classification using transformer models. We bucketed scores into each
interval class, resulting in 11 buckets. These 11 classes correspond to a
score from 0 to 10. Our methodology assigns each essay to a particular
category. If an essay is categorized by my model as being in Bucket 6, then it
receives a score of 5. Models like BERT, RoBERTa, ALBERT, DistilBERT, and XLM-
RoBERT444we used base models for our experiment were used to teach augmented
essays how to recognize multiple labels.

We approached this issue in the same manner as sentiment analysis, which
utilizes classification algorithms and yields extremely positive results. We
have discovered that using this data augmentation training method contributes
to an improvement in accuracy as referenced in Table V. From these results, it
is quite evident that BERT and RoBERTa outperform other models, although by a
small margin. The analysis of these findings demonstrates that utilizing
Transformer-based models is considerably superior than using LSTMs.

When used on top of huge pre-trained classification models, our data
augmentation strategies significantly enhance the performance of automated
essay grading. The accuracy of all those pre-trained models increases after
applying our augmentation technique. We believe this performance is because
after mixing a summary of the topic with each essay, it encourages the
internal representation of each essay to align with the topic, so that when we
test it on a essay with a different topic after fine tuning, it checks for the
alignment between the topic and the essay as a result of training and fine-
tuning and grades it accordingly.

##  VI Conclusion and Future work

In this paper, we automated essay grading using transformer based models and
an augmentation approach. The conclusion are as follows:

  * •

Pre-trained transformer-based models BERT, RoBERTa, ALBERT, DistilBERT, and
XLM-RoBERTa are very proficient at Automated Essay Scoring.

  * •

Data augmentation approaches might further enhance its performance for
analyzing lengthier resources such as essays to attain accurate essay scores.

In the future, we’ll come up with a better way to include elements that are
relevant to the topic instead of just a summary, so that training with this
data will lead to a more accurate model for automatically grading essays.

## References

  * [1] Ajay HB, Tillett PI, Page EB (1973) Analysis of essays by computer (AEC-II) (No. 8-0102). Washington, DC: U.S. Department of Health, Education, and Welfare, Office of Education, National Center for Educational Research and Development.
  * [2] Shermis MD, Mzumara HR, Olson J, Harrington S (2001) On-line grading of student essays: PEG goes on the World Wide Web. Assess Eval High Educ 26(3):247–259.
  * [3] Foltz PW, Laham D, Landauer TK (1999) The Intelligent Essay Assessor: Applications to Educational Technology. Interactive Multimedia Electronic Journal of Computer-Enhanced Learning, 1, 2, http://imej.wfu.edu/articles/1999/2/04/index.asp.
  * [4] Rudner, L. M., Liang, T. (2002). Automated essay scoring using Bayes’ theorem. The Journal of Technology, Learning and Assessment, 1(2).
  * [5] Dong F, Zhang Y, Yang J (2017) Attention-based recurrent convolutional neural network for automatic essay scoring. In: Proceedings of the 21st Conference on Computational Natural Language Learning (CoNLL 2017) p 153–162.
  * [6] Blood, I. (2011). Automated essay scoring: a literature review. Studies in Applied Linguistics and TESOL, 11(2).
  * [7] Burrows S, Gurevych I, Stein B (2015) The eras and trends of automatic short answer grading. Int J Artif Intell Educ 25:60–117. https://doi.org/10.1007/s40593-014-0026-8.
  * [8] Ke Z, Ng V (2019) “Automated essay scoring: a survey of the state of the art.” IJCAI
  * [9] Hussein, M. A., Hassan, H., Nassef, M. (2019). Automated language essay scoring systems: A literature review. PeerJ Computer Science, 5, e208.
  * [10] Klebanov, B. B., Madnani, N. (2020, July). Automated evaluation of writing–50 years and counting. In Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics (pp. 7796–7810).
  * [11] Vaswani, Ashish, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Łukasz Kaiser, and Illia Polosukhin. ”Attention is all you need.” Advances in neural information processing systems 30 (2017).
  * [12] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019. BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding. In Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), pages 4171–4186, Minneapolis, Minnesota. Association for Computational Linguistics.
  * [13] Liu, Yinhan, Ott, Myle, Goyal, Naman, Du, Jingfei, Joshi, Mandar, Chen, Danqi, Levy, Omer, Lewis, Mike, Zettlemoyer, Luke and Stoyanov, Veselin RoBERTa: A Robustly Optimized BERT Pretraining Approach. (2019). , cite arxiv:1907.11692.
  * [14] Lan, Zhenzhong, Chen, Mingda, Goodman, Sebastian, Gimpel, Kevin, Sharma, Piyush and Soricut, Radu. ”ALBERT: A Lite BERT for Self-supervised Learning of Language Representations..” Paper presented at the meeting of the ICLR, 2020.
  * [15] Sanh, Victor, Lysandre Debut, Julien Chaumond, and Thomas Wolf. ”DistilBERT, a distilled version of BERT: smaller, faster, cheaper and lighter.” arXiv preprint arXiv:1910.01108 (2019).
  * [16] Conneau, Alexis, Kartikay Khandelwal, Naman Goyal, Vishrav Chaudhary, Guillaume Wenzek, Francisco Guzmán, Edouard Grave, Myle Ott, Luke Zettlemoyer, and Veselin Stoyanov. ”Unsupervised cross-lingual representation learning at scale.” arXiv preprint arXiv:1911.02116 (2019).
  * [17] Wolf, Thomas, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac et al. ”Huggingface’s transformers: State-of-the-art natural language processing.” arXiv preprint arXiv:1910.03771 (2019)
  * [18] Lewis, Mike, Yinhan Liu, Naman Goyal, Marjan Ghazvininejad, Abdelrahman Mohamed, Omer Levy, Ves Stoyanov, and Luke Zettlemoyer. ”Bart: Denoising sequence-to-sequence pre-training for natural language generation, translation, and comprehension.” arXiv preprint arXiv:1910.13461 (2019).
  * [19] Ramesh, Dadi, and Suresh Kumar Sanampudi. ”An automated essay scoring systems: a systematic literature review.” Artificial Intelligence Review (2021): 1-33.
  * [20] Elijah Mayfield and Alan W Black. 2020. Should You Fine-Tune BERT for Automated Essay Scoring?. In Proceedings of the Fifteenth Workshop on Innovative Use of NLP for Building Educational Applications, pages 151–162, Seattle, WA, USA → Online. Association for Computational Linguistics.
  * [21] Ludwig, Sabrina, et al. ”Automated essay scoring using transformer models.” Psych 3.4 (2021): 897-915.
  * [22] Ormerod, Christopher M., Akanksha Malhotra, and Amir Jafari. ”Automated essay scoring using efficient transformer-based language models.” arXiv preprint arXiv:2102.13136 (2021).
  * [23] Park, Yo-Han, et al. ”EssayGAN: Essay Data Augmentation Based on Generative Adversarial Networks for Automated Essay Scoring.” Applied Sciences 12.12 (2022): 5803.
  * [24] Jong, You-Jin, Yong-Jin Kim, and Ok-Chol Ri. ”Improving Performance of Automated Essay Scoring by using back-translation essays and adjusted scores.” Mathematical Problems in Engineering 2022 (2022).
  * [25] Wang, Yongjie, et al. ”On the Use of BERT for Automated Essay Scoring: Joint Learning of Multi-Scale Essay Representation.” arXiv preprint arXiv:2205.03835 (2022).
  * [26] Sethi, Angad, and Kavinder Singh. ”Natural Language Processing based Automated Essay Scoring with Parameter-Efficient Transformer Approach.” 2022 6th International Conference on Computing Methodologies and Communication (ICCMC). IEEE, 2022.