# Group Convolutional Neural Networks Improve Quantum State Accuracy

Christopher Roth and Allan H. MacDonald  Physics Department, University of
Texas at Austin

###### Abstract

Neural networks are a promising tool for simulating quantum many body systems.
Recently, it has been shown that neural network-based models describe quantum
many body systems more accurately when they are constrained to have the
correct symmetry properties. In this paper, we show how to create maximally
expressive models for quantum states with specific symmetry properties by
drawing on literature from the machine learning community. We implement group
equivariant convolutional networks (G-CNN) cohen2016group , and demonstrate
that performance improvements can be achieved without increasing memory use.
We show that G-CNNs achieve very good accuracy for Heisenberg quantum spin
models in both ordered and spin liquid regimes, and improve the ground state
accuracy on the triangular lattice over other variational Monte-Carlo methods.

##  I Introduction

Because quantum physics problems with interacting degrees of freedom grow
exponentially with the number of degrees of freedom, accurate numerical
solutions are often available only for relatively small system sizes.
Simulations nevertheless play a valuable role in advancing understanding
because they have the important advantages over experiments on real physical
systems that the problem being solved is fully characterized and free from
unintended disorder. Neural network quantum states (NQS) (carleo2017solving ,
carrasquilla2020machine , sharir2020deep , hibat2020recurrent ,
roth2020iterative , nomura2020dirac , liang2021hybrid ,
astrakhantsev2021broken ) have emerged as a competitive tool for understanding
the low temperature properties of quantum many-body physics models. Unlike
traditional variational Monte-Carlo (VMC) methods, such as Gutzwiller
projection, neural networks have the advantage, and also the disadvantage, of
being free from inductive biases about the structure of the solution. They
compensate for the absence of an informed bias by using an enormous number of
parameters. As long as a wide enough model is used, neural networks contain
arbitrarily accurate solutions hornik1989multilayer .

Although reasonably large neural networks are guaranteed to harbor good
solutions in the space of possible parameters, there is no guarantee that
these solutions can be found in a reasonable time. One way to accelerate
training is to constrain the search space by removing non-solutions. This is
especially useful on lattice problems which have highly symmetric low-lying
energy levels. Recent research choo2019two , nomura2020dirac ,
ferrari2019neural , luo2021gauge , luo2020gauge , vieijra2020restricted  has
shown that forcing neural networks to have the correct symmetry tremendously
improves their ability to model the ground state and low-lying excited states
of quantum many-body systems.

Other research in the field of NQS has focused on training deeper models,
specifically convolutional neural networks (CNN) broecker2017machine ,
choo2019two , liang2021hybrid . This is motivated by the field of computer
vision he2016deep , szegedy2017inception , in which deep convolutional neural
networks have excelled at pattern recognition. While both images and quantum
ground states often have translational symmetry, quantum ground states
frequently have additional point group symmetry. In this paper we leverage
literature from the machine learning community cohen2016group  to generalize
convolutional NQS to the full wallpaper group of highly symmetric lattices.

We show that group equivariant convolutional networks (G-CNN) cohen2016group
provide the most expressive perceptron-like model for wavefunctions with
particular symmetry over discrete groups. G-CNNs are constructed using
equivariant convolutions, which are the most complex linear operations that
preserve the structure of a discrete symmetry group cohen2016group . This
approach contrasts with previous research choo2019two , nomura2020dirac ,
ferrari2019neural , which has utilized a symmetry-averaging procedure in which
the model is applied to symmetry-transformations of the input and the output
is averaged. We show that G-CNNs can be mapped to symmetry-averaged linear
models by zeroing out some filters between hidden layers. When we mask a G-CNN
in this manner, i.e. when we constrain it to match previous symmetry-averaging
procedures, the performance degrades.

We use a G-CNN with local translational filters to simulate frustrated
$J_{1}-J_{2}$ Heisenberg models on square and triangular lattices. Both of
these Hamiltonians give rise to rich phase diagrams that include magnetically
ordered states and spin liquid states that do not order but have non-trivial
quantum entanglement. Our model gives accurate results for both triangular and
square lattice phase-diagrams, and compares favorably with state-of-the-art
variational Monte-Carlo (VMC) results on the triangular lattice.

##  II VMC with Neural Network Quantum States

We compute ground state energies for the frustrated $J_{1}-J_{2}$ Heisenberg
model on the square and triangular lattices. The Hamiltonian is given by:

$\displaystyle
H=J_{1}\sum_{i,j\in\langle\rangle}{\bf\sigma}_{i}\cdot{\bf\sigma}_{j}+J_{2}\sum_{i,j\in\langle\langle\rangle\rangle}{\bf\sigma}_{i}\cdot{\bf\sigma}_{j},$
(1)

where $\langle\rangle$ and $\langle\langle\rangle\rangle$ are nearest and
next-nearest neighbor links and $J_{1},J_{2}>0$. We use a neural-network
ansatz that associates a complex number $\psi({\bf\sigma})$ with each
configuration of spins ${\bf\sigma}$. These define the wavefunction:

$\displaystyle|\psi>=\sum_{\sigma}\psi(\sigma)|\sigma>.$ (2)

We optimize our Ansatz using gradient based methods as detailed in Carleo et
al. carleo2017solving .

##  III Symmetric Neural Network Wavefunctions

Symmetry-averaging neural networks has been shown to improve performance in
both machine learning applications cohen2016group , cohen2019gauge ,
dieleman2016exploiting , cohen2018spherical , zaheer2017deep  and quantum
many-body VMC applications choo2018symmetries , nomura2020helping ,
ferrari2019neural . Below, we first discuss symmetry-averaging from the point
of view of equivariance. We then show that that G-CNNs provide a richer model
class, which can be mapped down to symmetry-averaged models by masking
filters.

###  III.1 The Principle of Equivariance

An equivariant function acts on a G-space cohen2016group , pitts2013nominal ,
a set of objects, $S$, that are related by the symmetry transformations of a
group, $G$. In this paper we consider the set of transformations, generated by
lattice translations and the d4(6) point group for square (triangular)
lattices. These are the groups of transformations that leave the lattices
unchanged. An equivariant function is one that transforms a G-set of poses
while preserving the group structure,

$\displaystyle g^{\prime}{\cal F}(\\{{\bf f}\\}_{S})={\cal F}(g\\{{\bf
f}\\}_{S}),$ (3)

where $\\{{\bf f}\\}_{S}$ is a set of features over the set $S$. Here
$g^{\prime}$ and $g$ do not be the same operator, but they must be generators
of isomorphic groups. For example if $g=T_{x,y}$ is a translation on a
periodic 2D square lattice of length L, then
$g^{\prime}=T_{x^{\prime},y^{\prime}}$ must also be a translation that is
cyclic over L in both dimensions.

###  III.2 Symmetry-Averaged Models are Equivariant Models

Symmetry can be imposed nomura2020helping  on a model wavefunction $\psi$ by a
symmetry-averaging procedure that applies the model to all symmetry-
transformations of the input and performs a phase-factor weighted average over
outputs:

$\displaystyle\psi({\bf\sigma})=\sum_{g}\chi_{g}{\tilde{\psi}}(g^{-1}{\bf\sigma}).$
(4)

Here $\chi_{g}$ is the character of symmetry operation $g$, and $g^{-1}$ is
the inverse of g, which means that $gg^{-1}=g^{-1}g=I$. This procedure can
also be thought of as averaging over the output of an equivariant model
$\psi_{g}({\bf\sigma})$ if we define:

$\displaystyle\psi_{g}({\bf\sigma})={\tilde{\psi}}(g^{-1}{\bf\sigma}).$ (5)

This definition recasts a single model ${\tilde{\psi}}$ applied to $N_{G}$
inputs, as an equivariant model $\psi_{g}$ with output dimension $N_{G}$
applied to a single input. Note that a symmetry transformation on the input to
$\psi_{g}$ yields a symmetry operation on the output:

$\displaystyle\psi_{g}(u^{-1}{\bf\sigma})=\tilde{\psi}(g^{-1}u^{-1}{\bf\sigma}).=\psi_{ug}({\bf\sigma}).$
(6)

We see that is still equivariant, as the set of transformations $\\{u\\}$ and
$\\{u^{-1}\\}$ form isomorphic groups. Using this new definition, the
symmetrization procedure expresses the wavefunction as follows:

$\displaystyle\psi({\bf\sigma})=\sum_{g}\chi_{g}\psi_{g}({\bf\sigma}).$ (7)

This still produces a $\psi$ that has the desired character:

$\displaystyle\begin{aligned}
\psi(u{\bf\sigma})=\sum_{g}\chi_{g}\psi_{g}(u{\bf\sigma})=\sum_{g}\chi_{g}\psi_{u^{-1}g}({\bf\sigma})\\\
=\sum_{g}\chi_{ug}\psi_{g}({\bf\sigma})=\sum_{g}\chi_{u}\chi_{g}\psi_{g}({\bf\sigma})=\chi_{u}\psi({\bf\sigma}).\end{aligned}$
(8)

From this discussion we see that symmetry-averaging is not required to
restrict symmetry eigenvalues. Only equivariance is needed. In this paper, we
replace the symmetry averaging procedure by a G-CNN, which consists of stacked
equivarant convolutional layers cohen2016group  interspersed with non-
linearities. The G-CNN architecture allows us to scale up the number of
parameters without using more memory, and improves performance on difficult
optimization problems.

###  III.3 Generalizing Translational Convolutions to Discrete Groups

(a)

(b)

Figure 1: Schematic of an equivariant convolutional model on a triangular
lattice with p6m symmetry. Nearest neighbor convolutions are drawn for
legibility. a) Equivariant convolution connecting the input to the feature
maps. Each feature map uses a set of filters with a specified symmetry
transformation over $d_{6}$. (b) Equivariant convolution connecting feature
maps at adjacent layers. There are 12 sets of filters that connect poses based
on their relative orientations.

Equivariance under the translational group is present in standard
convolutional neural networks (CNNs),

$\displaystyle C^{i}_{x,y}=\sum_{x^{\prime},y^{\prime}<L}{\bf
W}_{x^{\prime}-x,y^{\prime}-y}^{i}\cdot{\bf f}_{x^{\prime},y^{\prime}},$ (9)

that satisfy periodic boundary conditions. The filters ${\bf
W}_{x^{\prime}-x,y^{\prime}-y}^{i}$, look at patterns in the feature map ${\bf
f}_{x^{\prime},y^{\prime}}$ separated from ${x,y}$ by
$\\{x^{\prime}-x,y^{\prime}-y\\}$. CNNs work by alternating convolution
operations with non-linearities to develop increasingly abstract
representations of the input.

The G-CNN generalizes the CNN to act over a discrete group $G$, which may
contain non-commuting operations. For our implementation, we consider the full
wallpaper group, which adds rotation and mirror symmetry on top of
translation. The group convolution operation is written as follows:

$\displaystyle C^{i}_{g}=\sum_{h\in G}{\bf W}_{g^{-1}h}^{i}\cdot{\bf f}_{h}.$
(10)

For the G-CNN, the feature map is defined over the G-space generated by the
wallpaper group, instead of just the translation group. We see that the group
convolution is equivariant, since a group operation on input yields the same
group operation on the output:

$\displaystyle\sum_{h\in G}{\bf W}^{i}_{g^{-1}h}\cdot{\bf f}_{uh}=\sum_{h\in
G}{\bf W}^{i}_{g^{-1}u^{-1}h}\cdot{\bf f}_{h}=C^{i}_{ug}.$ (11)

This linear operation is the building block of our model.

###  III.4 Our Model

(a)

(b)

Figure 2: Filter shapes for our restricted equivariant convolution operation
on the a) triangular lattice (b) square lattice. For both lattices,
interactions are restricted to third nearest neighbors.

We use a model that is similar to a G-CNN cohen2016group , and specified by a
stack of equivariant convolutional layers interspersed with pointwise non-
linearities. The first convolution takes the z-components of the spins,
$\sigma_{\bf x}$, where ${\bf x}$ are the lattice site labels, and outputs a
feature map over the full wallpaper group:

$\displaystyle{\bf f}^{1}_{g}=\Gamma\Big{(}\sum_{{\bf x}}W_{g^{-1}{\bf
x}}^{0}\sigma_{\bf x}\Big{)},$ (12)

where $\Gamma$ is a point-wise non-linearity. This convolution is diagrammed
in Fig. 1(a). The feature-to-feature convolutions are performed as in equation
10:

$\displaystyle{\bf f}^{i+1}_{g}=\Gamma\Big{(}\sum_{h\in G}W_{g^{-1}h}^{i}{\bf
f}^{i}_{h}\Big{)}.$ (13)

The structure of these convolutions is seen in Fig. 1(b). We compute our
wavefunction with character $\\{\chi_{g}\\}$ by phase-weighting over the
exponential our final embedding, $f_{g}^{N}$, which has length one:

$\displaystyle\psi({\bf\sigma})=\sum_{g}\chi_{g^{-1}}\textrm{exp}(f_{g}^{N}).$
(14)

We choose $\Gamma$ to be the SELU nonlinearity klambauer2017self  applied
separately to the real and imaginary parts,

$\displaystyle\Gamma(x)=\textrm{SELU}(\textrm{Re}(x))+i\textrm{SELU}(\textrm{Im}(x)).$
(15)

The SELU nonlineariity moves the distribution of activations in the direction
of zero mean and one variance, which enables stable training of deep networks.
We note that for all models studied in this paper, the ground state is fully
symmetric, i.e. $\chi_{g}=1$ for all g.

While this model may appear complicated, we believe it is actually quite
simple. It is just a multi-layer perceptron, with alternating linear and
pointwise non-linear transformations, in which weights are tied together to
force particular symmetry properties. Training a simple model, with a constant
feature dimension and with neither bells nor whistles, will demonstrate the
power of this approach.

##  IV Results

We begin by showing results from our best performing model, which uses a G-CNN
architecture with local filters. This model achieves state-of-the-art VMC
ground state energies on the triangular lattice, and competitive energies on
the square lattice. We then examine how the perfomance of the G-CNN is
affected by masking off-diagonal filters. As detailed in section VIII.2, this
simplification maps G-CNNs to models that are forced to be equivariant by
symmetry-averaging. The optimization details for both of these experiments are
outlined in section VIII.5.

###  IV.1 Results from our Best Performing Model

In Table 1 we summarize a variety of ground state energy estimates obtained
using different methods for 6$\times$6 finite-size $J_{1},J_{2}$ triangular
lattice Heisenberg models. For this lattice size, the DMRGiqbal2016spin  and
and EDbernu1994exact  estimates agree to five figures and are essentially
exact. For the $J_{2}=0$ case our model achieves $\sim 0.2\%$ accuracy for the
energy, compared to the $\sim 1\%$ accuracy achieved by other variational
methods in the literature. For $J_{2}=J_{1}/8$, in the spin-liquid regime, the
ground state energy error is reduced by a factor of $\sim 3$ compared to
literature VMC estimates. VMC calculations often add Lanczos steps
sorella2001generalized  to improve variational energy accuracies; our ground
state energies surpasses all other VMC methods in accuracy even without this
elaboration.

Figure 3: Performance of the G-CNN on the triangular Heisenberg model ($J_{2}=\frac{1}{8}$) as various off-diagonal filters are set to zero. As described in section VIII.2 this is equivalent to replacing group convolution with symmetry averaging. All models are trained using the Adam optimizer with a learning rate of $3\times 10^{-3}$ | $J_{2}=0$ | $J_{2}=J_{1}/8$
---|---|---
ED bernu1994exact  | -0.5603734 | -0.515564
DMRG iqbal2016spin  | -0.560375 | -0.51557
G-CNN | -0.55922 | -0.51365
NN + Gutzwiller ferrari2019neural  | -0.553 | N/A
VMC + 2LS iqbal2016spin  | N/A | -0.512503
VMC kaneko2014gapless  | -0.55519 | -0.5089
VMC mezzacapo2010ground  | -0.55420 | N/A
VMC iqbal2016spin  | -0.548025 | -0.501788
Table 1: Comparison between finite size $6\times 6$ triangular lattice ground
state energies estimates

We studied the dependence of performance, measured by the accuracy of ground
state energy estimates, on the range of the filter for $J_{2}=J_{1}/8$.
Restricting the connectivity reduces computational overhead, and biases the
model towards learning short range interactions. Changing the filter shape
changes the way in which the model conveys information about the connectivity
of the lattice. We tried a few different symmetric filter shapes and settled
on a version that convolves over third-nearest-neighbors, as shown in Fig. 2.

In order to demonstrate the robustness of our model, we applied the same
architecture to the ordered and spin liquid regimes of the square lattice
$J_{1}-J_{2}$ Heisenberg model. Remarkably, across all four domains, we only
had to adjust a single learning rate. The results on the square lattice are
summarized in Table 2

| G-CNN | CNN + $c_{4}$ | RBM + PP
---|---|---|---
$J_{2}=J_{1}/2$ | $0.17\%$ | $0.39\%$ | $0.022\%$
$J_{2}=0$ | $0.016\%$ | $0.007\%$ | N/A
Max Memory | $16\times N_{g}$ | $10\times N_{g}$ | $576\times N_{g}$
$N_{params}$ | 113360 | 3838 | 22032
Table 2: Errors relative to exact diagnalization for $6\times 6$ square-
lattice finite-size ground state states. The number of variational parameters
and the memory usage are compared for different approaches.

We make a direct comparison to the square lattice CNN model described in Choo
et al. choo2019two  who use a translationally equivariant CNN symmetry-
averaged over $c_{4}$, while we use a model that is equivariant over the full
space group $p4m$. As explained in section VIII.2 of the appendix, their model
can be mapped to a G-CNN over $p4$ with masked off-diagonal rotational
filters. Our model uses slightly more memory, while their model has more
layers and uses an improved optimizer. We see that the G-CNN model
nevertheless has better ground state energy accuracy on the difficult-to-
simulate spin-liquid state. Neither CNN architecture is competetive with the
architecture detailed by Nomura et al. nomura2020dirac , who combine a
restricted Boltzmann machine with a pair product state (RBM+PP). Theirs is a
shallower model that uses far more memory. For the ordered state, both CNN
based methods are very accurate, and the discrepancy may be due to convergence
issues reddi2019convergence  with the adaptive moment (Adam) parameter
optimizer kingma2014adam .

We note that the equivariant model has by far the best parameter-number to
memory ratio, as it adds connectivity along the entire $d_{4}$ symmetry group.
Maximizing the number of parameters in a memory profile kitaev2020reformer ,
fedus2021switch  is important for performance improvement on massively
parallel modern hardware. Even though our model has more parameters than the
RBM + PP model, it is likely much faster to train on a GPU.

###  IV.2 Effect of Off-Diagonal Filters

We document the importance of full connectivity using the $J_{1}-J_{2}$
Heisenberg model on a triangular lattice with $J_{2}=J_{1}/8$. This is a
difficult system to simulate since it gives rise to a spin liquid state with
no magnetic order and non-trivial entanglement. We train four models with
different feature-to-feature connectivities to test the effect of replacing
layer-wise equivariance with symmetry averaging. We do this by applying a
masking procedure, where we set off-diagonal filters to zero. Here, the term
off-diagonal distinguishes between filters that connect poses of the same
orientation in a particular symmetry group and filters that connect different
orientations. As an example, filters in equation 9 with
$\\{x^{\prime},y^{\prime}\\}=\\{x,y\\}$ are translationally diagonal whereas
all others are translationally off-diagonal.

First we train a model with diagonal feature-to-feature interactions along all
groups, corresponding to a fully connected linear model symmetrized over
$p6m$. Second we train models with off-diagonal features in translation and
rotation set to zero. These correspond to models that convolve over some
degrees of freedom, and symmetry average over others. Finally, we compare all
of these with a G-CNN with full connectivity. While all four of these models
use the same amount of memory, they have vastly different numbers of
parameters.

Connectivity | Full | Translational | $d_{6}$ | Diagonal
---|---|---|---|---
$N_{params}$ | 449856 | 38016 | 13056 | 1616
Table 3: Comparison between the number of parameters of different equivariant
architectures

We see in Fig. 3 that adding a richer layer structure systematically improves
the performance. The model with fully-connected layers markedly outperforms
the others, surpassing their performance in $1/10$ of the training time. We
also see in Fig. 3 that the energy estimates are not monotonic functions of
training time. This behavior is related to our use of high learning rates that
enable escapes from local minima. Although higher learning rates worsen the
performance temporarily, we find that the final variational energies are
ultimately superior.

We see that translational filters have a bigger effect on performance than
rotational filters. This could simply be due to the number of parameters in
the models. On a $6\times 6$ lattice, adding translational interactions
increases the number of parameters by a factor of $36$, while adding $d_{6}$
interactions only increases the number of parameters by a factor of $12$.

##  V Discussion and Conclusion

In this paper we use a G-CNN to find the ground wavefunction of the frustrated
Heisenberg model on the square and triangular lattices. We find that the G-CNN
architecture yields accurate ground state energies for both magnetically
ordered states and spin liquids, including the most accurate VMC energies on
the triangular lattice to date.

Among perceptron-like models, G-CNN models are more expressive than symmetry-
averaged models. We explicitly demonstrate how performance degrades as the
model class is constrained. For the wallpaper group $p6m$, we see that it is
important to convolve over both translation operations and the point group.

For our experiments we used the simplest kind of neural network model that
satisfies the symmetry constraints, a multi-layer perceptron with tied
weights. We note that other advances in machine learning, such as depthwise
convolutions chollet2017xception  and attention vaswani2017attention , can
also be implemented with group equivariance romero2020attentive . Implementing
these models could certainly further improve the performance of the quantum
spin model ground state calculations we have discussed, and other quantum
ground state calclations. Although we only studied relatively small systems,
these calculations could be scaled to larger systems by combining
sparsification, transfer learning, and efficient sampling yang2020scalable .

Finally, we note that several other papers have studied the importance of
using equivariant models applied to different symmetry groups. Luo et al.
luo2020gauge  show that gauge equivarant neural networks improve performance
on gauge invariant Hamiltonians luo2021gauge , and Pfau et al pfau2020ab  use
a permutation equivariant model to account for identical electrons. Since many
important physics models have substantial symmetry, applications of
equivariant models are likely to be important for machine-learning methods in
computational physics.

##  VI Code

Code for implementing G-CNNs on NQS will be made publically available through
NetKet carleo2019netket  in the near future. Pytorch code for these
experiments will be uploaded at https://github.com/chrisrothUT/

##  VII Acknowledgements

Chris Roth acknowledges helpful interactions with Mohamed Hibat-Allah, Juan
Carrasquilla and Giuseppe Carleo. This work was supported in part by the Army
Research Office (ARO) Grant $\\#$ W911NF-17-1-0312 (MURI).

##  VIII Appendix

###  VIII.1 Implementation of the Equivariant Convolution

We implement our equivariant convolution exactly as in Cohen et al.
cohen2016group . We learn an index mapping that relates our filters of shape
$[d_{in},d_{out},N_{g}]$ to a group of 2D convolutions with rotated filters as
depicted in 1(b). Then we concatenate the rotation/reflection with the channel
dimension and apply a 2D convolution. For details see the aforementioned
paper.

###  VIII.2 Equivariant Models can be Mapped to other Symmetrized models by
Constraining Filters

Recently several groups nomura2020dirac , choo2019two , ferrari2019neural ,
szabo2020neural , liang2021hybrid  have accurately simulated the frustrated
$J_{1}-J_{2}$ Heisenberg model on a square lattice by constraining the
symmetry eigenvalues. While some of these models used translational
convolutions, none of them convolved over the entire wallpaper group $p4m$ as
done by the G-CNN.

We demonstrate how the G-CNN can be mapped down to other multi-layer linear
models by masking filters. For the sake a brevity, we show the mapping between
the p4m G-CNN and a) a linear model symmetry-averaged over p4m b) a
translationally equivariant model symmetry averaged over $d_{4}$.

We construct a symmetry operation in $p4m$ by first applying a reflection,
followed by a rotation and a translation,

$\displaystyle S_{x,y,\alpha,\beta}=T_{xy}R_{\alpha}R^{m}_{\beta}.$ (16)

Using equation 10, we can write the filter that connects two G-set elements in
neighboring layers.

$\displaystyle
S^{-1}_{x,y,\alpha,\beta}S_{x^{\prime},y^{\prime},\alpha^{\prime},\beta^{\prime}}=R^{m}_{-\beta}R_{-\alpha}T_{x^{\prime}-x,y^{\prime}-y}R_{\alpha^{\prime}}R^{m}_{\beta^{\prime}}.$
(17)

This is a convolution operation with filters transformed by $p4m$ as shown in
Fig. 1(b).

To map to a multi-layer linear model symmetry-averaged over $p4m$, the input-
to-feature filters are kept the same and the feature-to-feature filters are
restricted to be diagonal, $\alpha=\alpha^{\prime}$, $\beta=\beta^{\prime}$,
$x=x^{\prime}$ and $y=y^{\prime}$.

The new feature-to-feature convolution becomes

$\displaystyle{\bf f}^{i+1}_{x,y,\alpha,\beta}=\sigma\Big{(}W^{i}{\bf
f}^{i}_{x,y,\alpha,\beta}\Big{)}.$ (18)

This is equivalent to applying the same linear model to each pose. The
different poses are generated by the symmetry transformed filters shown in
Fig. 1(a). Equivalently, we could generate the poses by symmetry transforming
the input and keeping the filters the same. As a result, we see that the G-CNN
with diagonal feature-to-feature filters is equivalent to a linear model
symmetrized over $p4m$.

We can also map the G-CNN to a translationally equivariant model (CNN)
symmetrized over d4 by constraining $\alpha=\alpha^{\prime}$ and
$\beta=\beta^{\prime}$. With this constraint the feature-to-feature
convolution is

$\displaystyle{\bf
f}^{i+1}_{x,y,\alpha,\beta}=\sum_{x^{\prime},y^{\prime}}\sigma\Big{(}W^{i}_{x^{\prime}-x,y^{\prime}-y,\alpha,\beta}{\bf
f}^{i}_{x^{\prime},y^{\prime},\alpha,\beta}\Big{)}.$ (19)

The feature-to-feature convolution now consists of CNNs applied to each pose
in $d_{4}$ labeled by $\alpha,\beta$ where the filters of each CNN are
transformed according to the group element

$\displaystyle
W^{i}_{x^{\prime}-x,y^{\prime}-y,\alpha,\beta}=R^{m}_{-\beta}R_{-\alpha}W^{i}_{x^{\prime}-x,y^{\prime}-y}R_{\alpha}R^{m}_{\beta}$
(20)

Again, if we change basis and transform the input instead of the filters, this
is a CNN symmetrized over $d_{4}$. Both of these examples show that symmetry-
averaging creates a more restrictive model class than the G-CNN for multi-
layer networks.

###  VIII.3 Triangular Heisenberg on a Rectangular Grid

We implement the triangular Heisenberg model on a rectangular graph by noting
that a triangular graph is isometric to a rectangular graph with coupling
along one of the diagonals. We implement our models on $6\times 6$ tori with
periodic boundary conditions (PBCs) and coupling along one of the diagonals.
This is equivalent to a hexagonal shaped lattice with $36$ sites.

###  VIII.4 Coordinate Mappings

On the square lattice, a $\frac{\pi}{2}$ rotation is defined by the coordinate
mapping $(i,j)\rightarrow(-j,i)$ where the coordinates are defined modulo L.
On the triangular lattice a $\frac{\pi}{3}$ mapping is accomplished by a
coordinate mapping $(i,j)\rightarrow(i-j,i)$. We can arbitrarily choose a
reflection axis, and for simplicity we do the coordinate transform
$(i,j)\rightarrow(j,i)$ for both geometries.

###  VIII.5 Optimization Details

Our best performing model had $4$ layers and $16$ hidden nodes. We used the
Adam optimizer kingma2014adam  with a learning rate of $3\times 10^{-3}$ and
the default parameters. We optimized our learning rate and network
architecture on the triangular lattice with $J_{2}=0.125$ and used the same
hyperparameters for the other systems, except for the square lattice with
$J_{2}$ = 0 where we found that the learning rate was too high, and decreased
it by a factor of 3.

Our optimization procedure was as follows. First, we optimized the phase
structure for $500$ steps by setting $Re(log(\psi))=0$ before outputting
$\psi$. We found that pre-optimizing the phases szabo2020neural , which helps
avoid local-minima, was crucial to successfully finding the ground state. We
then trained both amplitudes and phases for $10^{4}$ steps at batch size $100$
before increasing the batch size to $1000$ and training for an additional
$2\times 10^{3}$ steps.

We tuned the learning rate (on the triangular lattice at $J_{2}=\frac{1}{8})$
by altering it by factors of 3 and modified Nhidden and Nlayers under the
constraint $N_{hidden}\times N_{layers}=64$ We found that a ratio of
$\frac{N_{hidden}}{N_{layers}}=4$ worked best.

## References

  * [1] Taco Cohen and Max Welling.  Group equivariant convolutional networks.  In International conference on machine learning, pages 2990–2999. PMLR, 2016.
  * [2] Giuseppe Carleo and Matthias Troyer.  Solving the quantum many-body problem with artificial neural networks.  Science, 355(6325):602–606, 2017.
  * [3] Juan Carrasquilla.  Machine learning for quantum matter.  Advances in Physics: X, 5(1):1797528, 2020.
  * [4] Or Sharir, Yoav Levine, Noam Wies, Giuseppe Carleo, and Amnon Shashua.  Deep autoregressive models for the efficient variational simulation of many-body quantum systems.  Physical review letters, 124(2):020503, 2020.
  * [5] Mohamed Hibat-Allah, Martin Ganahl, Lauren E Hayward, Roger G Melko, and Juan Carrasquilla.  Recurrent neural network wave functions.  Physical Review Research, 2(2):023358, 2020.
  * [6] Christopher Roth.  Iterative retraining of quantum spin models using recurrent neural networks.  arXiv preprint arXiv:2003.06228, 2020.
  * [7] Yusuke Nomura and Masatoshi Imada.  Dirac-type nodal spin liquid revealed by machine learning.  arXiv preprint arXiv:2005.14142, 2020.
  * [8] Xiao Liang, Shao-Jun Dong, and Lixin He.  Hybrid convolutional neural network and projected entangled pair states wave functions for quantum many-particle states.  Physical Review B, 103(3):035138, 2021.
  * [9] Nikita Astrakhantsev, Tom Westerhout, Apoorv Tiwari, Kenny Choo, Ao Chen, Mark H Fischer, Giuseppe Carleo, and Titus Neupert.  Broken-symmetry ground states of the heisenberg model on the pyrochlore lattice.  arXiv preprint arXiv:2101.08787, 2021.
  * [10] Kurt Hornik, Maxwell Stinchcombe, and Halbert White.  Multilayer feedforward networks are universal approximators.  Neural networks, 2(5):359–366, 1989.
  * [11] Kenny Choo, Titus Neupert, and Giuseppe Carleo.  Two-dimensional frustrated j 1- j 2 model studied with neural network quantum states.  Physical Review B, 100(12):125124, 2019.
  * [12] Francesco Ferrari, Federico Becca, and Juan Carrasquilla.  Neural gutzwiller-projected variational wave functions.  Physical Review B, 100(12):125131, 2019.
  * [13] Di Luo, Zhuo Chen, Kaiwen Hu, Zhizhen Zhao, Vera Mikyoung Hur, and Bryan K Clark.  Gauge invariant autoregressive neural networks for quantum lattice models.  arXiv preprint arXiv:2101.07243, 2021.
  * [14] Di Luo, Giuseppe Carleo, Bryan K Clark, and James Stokes.  Gauge equivariant neural networks for quantum lattice gauge theories.  arXiv preprint arXiv:2012.05232, 2020.
  * [15] Tom Vieijra, Corneel Casert, Jannes Nys, Wesley De Neve, Jutho Haegeman, Jan Ryckebusch, and Frank Verstraete.  Restricted boltzmann machines for quantum states with non-abelian or anyonic symmetries.  Physical review letters, 124(9):097201, 2020.
  * [16] Peter Broecker, Juan Carrasquilla, Roger G Melko, and Simon Trebst.  Machine learning quantum phases of matter beyond the fermion sign problem.  Scientific reports, 7(1):1–10, 2017.
  * [17] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [18] Christian Szegedy, Sergey Ioffe, Vincent Vanhoucke, and Alexander Alemi.  Inception-v4, inception-resnet and the impact of residual connections on learning.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 31, 2017.
  * [19] Taco Cohen, Maurice Weiler, Berkay Kicanaoglu, and Max Welling.  Gauge equivariant convolutional networks and the icosahedral cnn.  In International Conference on Machine Learning, pages 1321–1330. PMLR, 2019.
  * [20] Sander Dieleman, Jeffrey De Fauw, and Koray Kavukcuoglu.  Exploiting cyclic symmetry in convolutional neural networks.  In International conference on machine learning, pages 1889–1898. PMLR, 2016.
  * [21] Taco S Cohen, Mario Geiger, Jonas Köhler, and Max Welling.  Spherical cnns.  arXiv preprint arXiv:1801.10130, 2018.
  * [22] Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Ruslan Salakhutdinov, and Alexander Smola.  Deep sets.  arXiv preprint arXiv:1703.06114, 2017.
  * [23] Kenny Choo, Giuseppe Carleo, Nicolas Regnault, and Titus Neupert.  Symmetries and many-body excitations with neural-network quantum states.  Physical review letters, 121(16):167204, 2018.
  * [24] Yusuke Nomura.  Helping restricted boltzmann machine with quantum-state representation by restoring symmetry.  arXiv preprint arXiv:2009.14777, 2020.
  * [25] Andrew M Pitts.  Nominal sets: Names and symmetry in computer science, volume 57.  Cambridge University Press, 2013.
  * [26] Günter Klambauer, Thomas Unterthiner, Andreas Mayr, and Sepp Hochreiter.  Self-normalizing neural networks.  arXiv preprint arXiv:1706.02515, 2017.
  * [27] Yasir Iqbal, Wen-Jun Hu, Ronny Thomale, Didier Poilblanc, and Federico Becca.  Spin liquid nature in the heisenberg j 1- j 2 triangular antiferromagnet.  Physical Review B, 93(14):144411, 2016.
  * [28] B Bernu, P Lecheminant, C Lhuillier, and L Pierre.  Exact spectra, spin susceptibilities, and order parameter of the quantum heisenberg antiferromagnet on the triangular lattice.  Physical Review B, 50(14):10048, 1994.
  * [29] Sandro Sorella.  Generalized lanczos algorithm for variational quantum monte carlo.  Physical Review B, 64(2):024512, 2001.
  * [30] Ryui Kaneko, Satoshi Morita, and Masatoshi Imada.  Gapless spin-liquid phase in an extended spin 1/2 triangular heisenberg model.  Journal of the Physical Society of Japan, 83(9):093707, 2014.
  * [31] Fabio Mezzacapo and J Ignacio Cirac.  Ground-state properties of the spin-antiferromagnetic heisenberg model on the triangular lattice: a variational study based on entangled-plaquette states.  New Journal of Physics, 12(10):103039, 2010.
  * [32] Sashank J Reddi, Satyen Kale, and Sanjiv Kumar.  On the convergence of adam and beyond.  arXiv preprint arXiv:1904.09237, 2019.
  * [33] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [34] Nikita Kitaev, Łukasz Kaiser, and Anselm Levskaya.  Reformer: The efficient transformer.  arXiv preprint arXiv:2001.04451, 2020.
  * [35] William Fedus, Barret Zoph, and Noam Shazeer.  Switch transformers: Scaling to trillion parameter models with simple and efficient sparsity.  arXiv preprint arXiv:2101.03961, 2021.
  * [36] François Chollet.  Xception: Deep learning with depthwise separable convolutions.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 1251–1258, 2017.
  * [37] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  arXiv preprint arXiv:1706.03762, 2017.
  * [38] David Romero, Erik Bekkers, Jakub Tomczak, and Mark Hoogendoorn.  Attentive group equivariant convolutional networks.  In International Conference on Machine Learning, pages 8188–8199. PMLR, 2020.
  * [39] Li Yang, Wenjun Hu, and Li Li.  Scalable variational monte carlo with graph neural ansatz.  arXiv preprint arXiv:2011.12453, 2020.
  * [40] David Pfau, James S Spencer, Alexander GDG Matthews, and W Matthew C Foulkes.  Ab initio solution of the many-electron schrödinger equation with deep neural networks.  Physical Review Research, 2(3):033429, 2020.
  * [41] Giuseppe Carleo, Kenny Choo, Damian Hofmann, James ET Smith, Tom Westerhout, Fabien Alet, Emily J Davis, Stavros Efthymiou, Ivan Glasser, Sheng-Hsuan Lin, et al.  Netket: A machine learning toolkit for many-body quantum systems.  SoftwareX, 10:100311, 2019.
  * [42] Attila Szabó and Claudio Castelnovo.  Neural network wave functions and the sign problem.  Physical Review Research, 2(3):033075, 2020.


#  Robust Learning-Based Incipient Slip Detection using the PapillArray
Optical Tactile Sensor for Improved Robotic Gripping

Qiang Wang, Student Member, IEEE, Pablo Martinez Ulloa, Robert Burke, David
Cordova Bulens,
and Stephen J. Redmond, Senior Member, IEEE This publication has emanated from
research conducted with the financial support of China Scholarship Council
under grant number 202006540003 and of Science Foundation Ireland under grant
numbers 17$/$FRL$/$4832 and SFI$/$12$/$RC$/$2289$\\_$P2.The authors are with
the School of Electrical and Electronic Engineering, University College
Dublin, Belfield, Dublin, D04 V1W8, Ireland (e-mail:{qiang.wang,
pablo.martinezulloa<EMAIL_ADDRESS>{david.cordovabulens,
stephen.redmond}@ucd.ie).Corresponding author: Stephen J. Redmond.

###### Abstract

The ability to detect slip, particularly incipient slip, enables robotic
systems to take corrective measures to prevent a grasped object from being
dropped. Therefore, slip detection can enhance the overall security of robotic
gripping. However, accurately detecting incipient slip remains a significant
challenge. In this paper, we propose a novel learning-based approach to detect
incipient slip using the PapillArray (Contactile, Australia) tactile sensor.
The resulting model is highly effective in identifying patterns associated
with incipient slip, achieving a detection success rate of 95.6% when tested
with an offline dataset. Furthermore, we introduce several data augmentation
methods to enhance the robustness of our model. When transferring the trained
model to a robotic gripping environment distinct from where the training data
was collected, our model maintained robust performance, with a success rate of
96.8%, providing timely feedback for stabilizing several practical gripping
tasks. Our project website: https://sites.google.com/view/incipient-slip-
detection.

##  I Introduction

###  I-A Background

Autonomous robots have yet to achieve human-like dexterity when performing
gripping tasks, mainly due to a lack of satisfactory tactile perception and
processing abilities. Studies have shown that even humans struggle with simple
gripping tasks in the absence of tactile sensation [1, 2]. The palm of the
human hand contains $\sim$17,000 mechanoreceptors, i.e., specialized nerve
endings that respond to mechanical stimuli such as deformation, pressure, and
displacement [3, 4]. These receptors play a crucial role in sensing and
relaying tactile information to the nervous system [5], allowing humans to
adjust their grip in real-time to account for slipperiness and other factors.
Building on these insights, researchers have designed tactile sensors
replicating part of human hand sensing capabilities and explored slip
detection techniques using these sensors to enhance robotic manipulation
performance [6, 7, 8].

####  I-A1 Types of slip

The two main types of slip are gross slip and incipient slip. Gross slip
refers to the occurrence of slip across the entire contact surface, where the
relative motion between the gripper or tactile sensor and the gripped object
is typically observable at a macro level [9, 10]. On the other hand, incipient
slip refers to the initial stage of slip, when parts of the contact surface
slip while others remain stuck [9, 10]. For example, when an object is held by
elastic fingertips, and an external force is applied to the object in a
direction tangential to the contact surface, some parts of the fingertips will
stretch while others will compress, causing incipient slip at the periphery of
the contact surface while the central part remains stuck. As the applied force
increases, the slip will finally spread across the entire contact surface,
leading to gross slip. Throughout the incipient slip phase, there may not be
any observable relative motion between the object and the finger.

####  I-A2 Slip detection and challenges

Previous studies have proposed techniques to detect gross slip and apply
corrective measures when the slip is detected to prevent objects from dropping
out of the grasp [7, 8, 11]. Detecting gross slip may not always be a wise
strategy, as it occurs when the entire contact has already started slipping.
On the other hand, detecting incipient slip can provide an early warning of an
impending and more dangerous gross slip, allowing corrective measures to be
applied earlier, and increasing the likelihood of maintaining a safe grip.
However, detecting incipient slip is not trivial because it requires the
contact interface of the sensor to possess adequate elasticity, enabling one
part to undergo sufficient and detectable deformation, resulting in slip,
while the other part remains stuck. Furthermore, validating incipient slip can
be challenging since it is not generally associated with macro-level relative
movement between the sensor/finger and the object. To verify the occurrence of
incipient slip, researchers commonly utilize a camera to monitor the contact
surface; by examining the camera images, they can visually confirm the
presence of incipient slip events [12, 13]. However, this method of relying on
cameras may not be feasible in real-world situations, such as when gripping
everyday objects.

Figure 1: (a) Illustration of two versions of the Contactile PapillArray Dev
Kit sensor that were used in this study: the upper shows version 1 (v1) and
the bottom shows v2. (b) Illustration of robotic gripping rig used to deploy
our trained model for online evaluation.

###  I-B Our contribution

Our study presents a new technique for detecting incipient slip using the
PapillArray (Contactile, Australia) tactile sensor. This sensor features a
square array of nine elastic silicone pillars with varying unloaded heights,
promoting different normal forces on the pillars when pressed against a
surface. This design enhances the likelihood of inducing incipient slip on
shorter pillars when a tangential force is applied.

We utilized deep neural networks (NN) to develop our incipient slip detection
algorithm, where we made novel use of the data gathered in a previous study
[13] to construct the dataset for training and evaluating the NN. The primary
objective of the NN was to classify inputs into two distinct categories:
incipient slip and other, functioning as a binary classifier; other refers to
all others states that are not incipient slip, such as gross slip or being
stationary. Furthermore, the tactile data at hand is presented in the form of
a uniformly-sampled time series. Therefore, to effectively capture the serial
nature of the data, we utilize a recurrent neural network (RNN) [14]. The
inclusion of historical data in a NN model has the potential to enhance its
performance in real-time prediction tasks, as it enables the capture of
temporal patterns and dependencies, leading to more robust and accurate
forecasts [14]. We also propose several data augmentation methods designed to
enhance the performance and robustness of our trained model, making it robust
to environmental confounders.

##  II Related work

Similar to the approach we will take in this paper, the approach proposed in
[7] treats slip detection as a classification task; the authors employed a
support vector machine [15] to detect slip using the velocity of embedded pins
on the inner surface of a TacTip camera-based tactile sensor [16]. Labels of
the training data are assigned manually based on the alignment of pin
velocities. In a more recent study [12], the authors modified the TacTip
sensor used in [7] by introducing raised fingerprint-like ridges, decreasing
skin thickness, and increasing pin spacing to reduce mechanical coupling
between ridges and to create the traction differential and facilitating the
shear displacement required for the occurrence of incipient slip. This is
similar to the behavior seen on the human finger pad when sheared against an
object, thus allowing the sensor to experience incipient slip. They used an
external camera to monitor the contact in real-time for data labeling, and
then employed a convolutional neural network [17] as a binary classifier to
detect incipient slip.

The GelSight technology is another camera-based tactile sensing system that
uses an elastic body to establish a contact with an object, with the built-in
camera recording the resulting deformation to obtain tactile data [18]. An
approach was introduced in [19] for detecting incipient slip using the
GelSight sensor. This method determines the degree of incipient slip by
analyzing the inhomogeneity of the displacement field, which is quantified in
terms of entropy. More recently, a more advanced version of the GelSight
technology, called GelSlim, was proposed in [20]; it employed the deviation of
the deformation field from a 2D planar rigid displacement field to determine
slip.

Compared to camera-based tactile sensors, the distributed optical sensor used
in our work, the PapillArray, is less complex in terms of instrumentation[21].
It offers several advantages over other sensor designs, including size,
temporal resolution, and compliance. A heuristic algorithm that employs the
PapillArray tactile sensor to detect incipient slip is proposed in [6]. The
approach is based on the observation that incipient slip happens when some
sensor pillars stop deflecting at the same rate as the contacted object is
moving in the sensor’s frame of reference. Precisely, this approach detects
slip by evaluating the tangential velocity drop with respect to a reference
pillar, which is the pillar under the highest normal force (usually the
center). In the case of rotational movements, with the center of rotation at
the center pillar, the algorithm cannot detect any slip since no movement can
be detected in the center pillar. This heuristic approach is further improved
in [13] to account for rotational slips, detecting the deceleration of each
pillar by comparing it to its own recent maximum velocity, and then it checks
if other pillars are still in motion to confirm that the deceleration
indicates an incipient slip. However, these methods may not be applicable when
dealing with deformable or non-planar surfaces, or when only a subset of the
pillars makes contact with the object. In such cases, establishing a
dependable reference pillar to represent the object’s movement in [6] becomes
challenging; in [13], it is difficult to determine whether the deceleration of
pillars is caused by slip or by the shape of the object’s surface.

In our work, we are motivated to take a learning-based approach in developing
a dedicated incipient slip detection algorithm, where we propose domain
adaptation techniques to enhance the robustness of our trained model, enabling
it to effectively detect incipient slip for more realistic objects and
contacts, overcoming the challenges outlined above.

##  III Materials and methods

###  III-A Hardware

####  III-A1 Contactile sensor

Our study employed the commercial PapillArray sensor from
Contactile111https://contactile.com/, depicted in Fig. 1, which is based on
the concept described in [22]. The sensor outputs the real-time $x-y-z$ force
data experienced by each pillar at a high sampling rate of 1,000 Hz. Our
training data was collected using the Dev Kit v1, while for the online
evaluation of our trained model, we used the Dev Kit v2. Dev Kit v2 and Dev
Kit v1 differ in size and the pillar Shore hardness.

####  III-A2 Robotic gripping rig

Fig. 1 displays the rig used in our study for the gripping task. The rig
features a specialized two-finger gripper (RG2, OnRobot, Germany) with a blue
adapter fixed to one of its fingers. This adapter serves to couple the
Contactile PapillArray Dev Kit v2 sensor to the gripper finger. A white
3D-printed cuboid is used to extend another finger, matching the length of the
finger equipped with the sensor. Moreover, a couple of ArUco markers are
attached to this extended cuboid to track the gripper’s pose. We replaced the
original motor of the RG2 gripper with a stepper motor (MX-28, Dynamixel, US)
to achieve high-frequency interruptible control of the gripper. The modified
gripper was mounted on a six-axis robot arm (UR5e, Universal Robots, Denmark).

###  III-B Data preparation

####  III-B1 Collect slip data and annotate slip events for individual pillars

Our training dataset is sourced from [13]. In brief, the training data was
acquired using a six-degree-of-freedom hexapod robot (H-820, Physik
Instrumente, Germany) with the Contactile PapillArray Dev Kit v1 sensor
mounted on the top. A transplant acrylic plate is fixed above the sensor on a
T-slot frame and a video camera (Logitech Streamcam, Logitech, Switzerland) is
positioned above the acrylic plate to capture videos of the contact between
the sensor and the plate. During the data collection, the hexapod pushes the
sensor vertically against the acrylic place and then moves it laterally to
induce a slip. The horizontal movement could be a translation, a rotation, or
a combination of both. A total of 200 data sequences were collected, covering
a range of compression levels, hexapod movement velocities, and movement
directions. The recorded videos are processed using the Matlab Computer Vision
Toolbox (MathWorks, USA) to track the pillar tip position. The tangential
pillar tip velocity is then used to label the slip state (gross slip or not
gross slip) of individual pillars.

####  III-B2 Collect control data

When the sensor is compressed against a flat surface and moved laterally, the
tangential velocity measured by each pillar will increase at first, as the
sensor starts deforming, before reaching a peak velocity and subsequently
decreasing its speed when the pillar stops deforming (Fig. 2). If a pillar
stops deforming because it is undergoing incipient slip, at least one other
pillar will still be deforming laterally; this is observed by an asynchronous
decrease of the tangential velocity of the nine pillars (Fig. 2 \- Slip).
However, if the object stops moving before any slip occurs, the tangential
velocity magnitude of the nine pillars decreases almost simultaneously (Fig. 2
\- Stop).

Since the stop events display similar temporal feature to slip events, we
collected an additional dataset specifically focusing on stop events,
consisting of a total of 28 data sequences. We label the data points in these
sequences as other. By incorporating this dataset, the NN is less likely to
misclassify between incipient slip and other, thereby improving the accuracy
and reliability of the NN. The data collection process was similar to that of
the slip events, except that the hexapod’s movement was abruptly halted before
any slip occurred. Further details on this process can be found in [13].

Figure 2: Changes in force and yank (rate of change of force) over time
between slip and stop events. The X area (red) and Y area (green) represent
the respective components in the $x$ and $y$ directions. The compression of
the tallest pillar in these two examples is 1.3 mm, and the movement direction
is at a 45∘ angle in the $x-y$ directions, with a speed of 8 mm.s-1.

####  III-B3 Annotate the incipient slip

Based on the definition of incipient slip provided in Section I, we annotate
the incipient slip in the dataset as follows: we consider that incipient slip
has occurred when at least one pillar slips with respect to the contact
surface, while at least one other pillar remains stationary with respect to
the contact surface. In other words, we start annotating incipient slip from
the moment the first slip occurs on any pillar, and this interval continues
until the time when all nine pillars have slipped. The slip label of each
pillar is obtained as described in Section III-B1. It should be noted that
when annotating incipient slip in the rotational data, we only consider the
outer eight pillars. This is because the rotational movement is centered
around the central pillar, which never slips by our definition (remains in the
same location on the contact area), for our data set.

Figure 3: (a) Illustration of how sequence data is input into a model to
generate classifications for each window. The input is represented by the red
rectangle, which undergoes window reshaping, denoted by the dotted rounded
rectangles. $T$ represents the number of data points in one data sequence, and
in our case $T=40$ samples at 1,000 Hz. The transparent circles represent the
hidden states; the initial hidden layer, $h_{0}$, is initialized as a vector
of zeros with the same shape as the hidden state of the GRU cells. The encoder
is shown as grey trapezoids connected to blue GRU cells, with purple circles
representing recurrent computations. The orange trapezoids represent the fully
connected estimator, and the green items represent the output estimations. (b)
Illustration the aggregation process when the ensemble model is deployed,
using $Z=5$ classifiers in the emsemble.

####  III-B4 Refine data sequence

The sensor output exhibits variance due to noise and sporadically produces
glitches that deviate significantly from the mean value, displaying sudden
extreme highs or lows. To address these issues, we apply a median filter with
a window size of 21 samples on the raw sensor signal, which is sampled at
1,000 Hz.

We divided the raw data sequence into non-overlapping windows, with each
window containing 40 samples. This division reduced the data rate to 25 Hz.
This was done for practical limitations in the hardware and software of our
system. More precisely, the maximum refresh rate of our gripper servo is
$\sim$62 Hz, and the computation rate of our classifier is $\sim$40 Hz.
Moreover, it is worth noting that reliable gripping does not necessarily
require a high sampling frequency. Indeed, humans have a reaction time of
approximately 80-120 ms (equivalent to 8.3-12.5 Hz) [23], enabling us to
perform most everyday gripping tasks effectively.

Finally, we only consider the $x-y$ forces on the pillars as input in NN
training, while excluding the $z$ force. During the data collection process,
when the hexapod moves tangentially to induce slip, it remains stationary in
the $z$ direction. As a result, we assume that the $z$ force does not play a
significant role in detecting incipient slip in our case. It should be
acknowledged that in real-world scenarios, the normal force can provide
valuable information for humans to detect slip, and it is likely to vary
appreciably for different gripping objectives. Therefore, another reason for
excluding the $z$ force is to prevent the NN from incorrectly learning that
the $z$ force remains relatively stable during slip events, as occurs in our
data set.

###  III-C Training data augmentation

####  III-C1 Data augmentation by rotational symmetry

During the data collection process, the sensor is placed at the origin of the
world coordinate frame. Its horizontal surface is parallel to the $x-y$ plane
of the world frame of reference, and the side edges align with the $x-y$ axis
directions. Hence we use a rotation transformation to augment the data;
intuitively, it can be understood as rotating the initial position of the
sensor around the $z$ axis by a random angle. For each data point in a
sequence, we perform the following mathematical calculations:

$\begin{bmatrix}F_{x^{\prime}}\\\
F_{y^{\prime}}\end{bmatrix}=\begin{bmatrix}\cos(\theta)&-\sin(\theta)\\\
\sin(\theta)&\cos(\theta)\end{bmatrix}\cdot\begin{bmatrix}F_{x}\\\
F_{y}\end{bmatrix},\theta\in[0,2\pi),$ (1)

where $F_{x}$ and $F_{y}$ represent the force values along the original $x-y$
axis, and $F_{x^{\prime}}$ and $F_{y^{\prime}}$ are the augmented force values
after virtual rotation of the sensor by a randomly sampled angle, $\theta$,
from a uniform distribution of $[0,2\pi)$.

####  III-C2 Advanced data augmentation for domain adaptation

The data used in our study was collected under idealized conditions, where a
hexapod robot was used to compress the sensor against a flat surface and move
laterally in a controled manner. In this setup, the force was nearly
perpendicular or parallel to the contact surface and the movement speed is
nearly constant. However, in real-world robotic gripping, the conditions are
expected to be quite different from this idealized setup, and the performance
of the model trained on such data is expected to be poor. We identify several
issues that may arise when transferring the model trained on idealized data to
real-world gripping scenarios, and we propose a range of advanced data
augmentation methods to address these issues in the following paragraphs.
These methods are designed to generate synthetic data that mimics the real-
world variability of the gripping:

  * •

Issue: The slipping velocity in real-world robotic gripping is not constant,
as it is influenced by various factors such as gravity, friction, and the
shape of the object being gripped. However, during the data collection
process, the hexapod induces slip at a constant velocity. Remedy:  We employ
random sampling to sample a percentage of data points from the raw data
sequence, thereby generating a new data sequence. And we maintain the
frequency of the new sequence at the same rate as the raw sequence (1,000 Hz).
This approach can simulate velocity variations to mimic real-world gripping
scenarios, as it changes the magnitude differences of some temporally adjacent
data points while keeping the time interval unchanged.

  * •

Issue: In some gripping scenarios, a portion of the sensor pillars may not be
in contact with the object. For instance, this can occur when employing
sensors to grip an object with a rounded surface or when gripping an object
smaller than the sensor’s contact area. Remedy: To simulate an unloaded
pillar, we substitute a number of pillar data sequences with zero sequences.
Noise is then added to make the generated sequence resemble a realistic sensor
signal. The noise is derived from a normal distribution with a mean of 0.0 N
and a standard deviation of 0.001 N.

  * •

Issue: Unlike with the hexapod, the force generated by the gripper may not be
perfectly perpendicular to the $x-y$ plane of the sensor frame of reference,
and the force leading to slip may not be perfectly in this plane. This can
occur when the gripped object is not flat or the mechanical linkage of the
gripper flexes when applying force to the object. Remedy: First, we sampled
nine individual pillar sequences from raw sensor sequences with different
sensor compression levels and hexapod movement types, and then combined them
to form a new sensor sequence. Secondly, we scaled (scale factor ranging from
0.2 to 2.0) the magnitude of values for a number of pillar sequences. Lastly,
we randomly permuted the position (by pillar index) of a nine-pillar sequence.
Employing these techniques can encourage the NN capture a broader and more
comprehensive pattern of incipient slip (see Section III-B3), rather than only
learning the limited pattern introduced by the hexapod.

###  III-D Neural networks

The key decision making component of our incipient slip detection approach is
a binary classifier. Initially, we trained a NN capable of estimating the
probability of incipient slip for each time point in a sequence. Next, we set
a threshold to convert the continuous probability into a binary output. To
enhance the accuracy of the classifier, we used an ensemble technique that
trains multiple independent classifiers concurrently and aggregates their
output probabilities to produce the final decision (shown in Fig. 3).

####  III-D1 Architecture

Fig. 3 illustrates the process of inputting a data sequence into the NN and
obtaining the corresponding slip classification. The modified data sequence,
as explained in Section III-B4, is input into an encoder. Subsequently, the
encoder output is passed to a specific type of RNN called a gated recurrent
unit (GRU) [24]. In our approach, we utilize a single layer of GRU for each
propagation step, and we refer to it as a GRU cell. The hidden output from the
GRU cell is generated as a combination of the current input and historical
information. Moreover, an estimator is included that takes the hidden layer
output from the GRU cell and converts it into a probability estimation. The
ground truth label of each window is determined by the label of the last
sample in the window.

####  III-D2 Training

The ensemble model consists of $Z$ ($Z=5$ in our case) independently trained
classifier models. During each training iteration of each classifier model, a
subset comprising a proportion of $\lambda$ sequences ($\lambda=40\%$ in our
case) is randomly sampled with replacement from the entire training set and
used for NN training. The final layer of the estimator utilizes a two-class
softmax activation function, with its outputs interpreted as probabilities for
the occurrence of incipient slip and other. Our chosen loss function is binary
cross-entropy.

####  III-D3 Decision making

We aggregate the output probability from each classifier model in the ensemble
to convert the continuous probability to binary prediction:

$f:=\mathds{1}\left[\frac{\sum_{z=1}^{Z}M_{z}(x=[F_{(n-1)T+1},\cdot\cdot\cdot,F_{nT}])}{Z}>P_{th}\right],$
(2)

where $\mathds{1}[\cdot]$ is an indicator function, $M_{z}$ donates the
$z^{th}$ classifier model in the ensemble, $x$ donates the input vector, and
$P_{th}$ denotes the probability threshold, which is 50% in our work. $Z$
donates the number of classifiers in the ensemble model.

Figure 4: The offline evaluation results. (a) Illustration of two examples of
successfully classified sequences, one corresponding to slip event and the
other to stop event. (b) Illustration of the confusion matrix for the
augmented test set, where TP represent successful detection of incipient slip,
FN represents incipient slip is detected in a sequence of stop event, FP
represents nothing is detected in a sequence of slip event, and TN represents
the absence of detection in stop events, i.e. successful classification. (c)
Illustration of the latency for successful TP cases in predicting incipient
slip. This figure was created with reference to the moment when the first of
the nine pillars slips (the ground truth for the onset of incipient slip).

##  IV Experiments and results

We first explicitly display our method’s high success rate in detecting
incipient slip, including offline and online scenarios. Then, we illustrate
the practical benefits of our approach by showcasing its ability to stabilize
an insecure robotic grasp in a number of practical gripping tasks.

###  IV-A Offline evaluation

The entire dataset is randomly split into two subsets: a training set
($\sim$80% of the entire dataset, comprising 160 data sequences of slip event
and 23 data sequences of stop event) for model training, and a test set
($\sim$20% of the entire dataset, consisting of 40 data sequences of slip
event and 5 data sequences of stop event) for model evaluation. Both subsets
are expanded through the symmetry-based augmentation method described in
Section III-C1, resulting in a five-fold increase in the size of the training
set and test set.

Fig. 4 displays two examples comparing the incipient slip detection results
over slip and stop events. As observed, the algorithm’s confidence in labeling
incipient slip increases rapidly as incipient slip starts and decreases as it
progresses toward gross slip. In comparison, the probability in the stop case
fluctuates slightly but remains well below the threshold.

We define an incipient slip detection as successful if it occurs within a 0.3
second window preceding the true labeled time point of incipient slip (to
accommodate the error of the ground truth) and prior to the occurrence of the
gross slip. For the stop event, a successful estimation is defined as a
classification of the entire sequence as other.

Fig. 4 presents the confusion matrix, displaying the final classification
results over the entire test set; our algorithm achieves an overall success
rate of $\sim$95.6%. The results also demonstrate its effectiveness in
differentiating between the slip and stop events; this indicates that our
algorithm is not simply detecting the changes in the force and yank of the
pillars, as mentioned earlier in Section III-B2.

Our algorithm can effectively detect incipient slip in its early stages. In
Fig. 4, we present the latency between the moment of incipient slip detected
by the algorithm and the ground truth onset of incipient slip. It is evident
that, on average, incipient slip can be detected within 10 ms of its
initiation.

Figure 5: Everyday objects for our online gripping task.

Figure 6: (a) Illustration of a confusion matrix for the classification. (b)
Illustration of our validation of incipient slip in online testing. The $x$
axis in each sub-figure represents the normalized robot translational/angular
displacement, calculated as $D_{norm}=D/D_{th}$, where $D$ is the displacement
before normalizing and $D_{th}$ (2 mm for translation and 2∘ for rotation) is
an assumption that indicates the reasonable minimum distance at which the
gross slip must have happened or is likely to happen. To calculate $D_{norm}$
for the compound Translation + Rotation sequences, we calculate the $D_{norm}$
for both the translation and rotation displacements respectively and then take
their mean value. $D_{norm}$ facilitates the intuitive visualization of when
incipient slip is detected during movement of the robot and also facilitates
comparison across three different types of movement. The blue gradient fill
shows increasing confidence of gross slip, starting at 0.0 when the robot
starts moving and reaching 1.0 at $D_{th}$.

###  IV-B Online evaluation

In the online evaluation stage, we utilized the full data set for training the
final deployed model. Again, to increase the amount of training data, we
applied both symmetry-based (see Section III-C1) and advanced data
augmentation (see Section III-C2) techniques, resulting in a five-fold
increase in data amount (1140 data sequences).

The online evaluation was performed on six everyday objects, depicted in Fig.
5. We include objects of varying surface materials, curvatures, and hardness
to ensure a broad range of conditions are represented in our results.

####  IV-B1 Validating incipient slip detections

We cannot easily validate incipient slip occurrences for everyday objects as
we cannot independently monitor individual pillar contacts. Hence, we choose
to perform the online evaluation based on following well-founded assumptions.
The incipient slip detection is considered successful if it can be detected at
any time-point between the time when the robot’s movement begins ($T_{m}$) and
the time when gross slip occurs ($T_{g}$); the criterion for determining the
occurrence of gross slip has been arbitrarily defined as the occurrence of
relative translational movement greater than 2 mm or relative rotational
movement exceeding 2∘ between the object and the robot’s frame of reference.

To induce a slip, the gripper first grips the object with a constant force.
Then the robot moves the gripper downwards towards a rigid and stationary
table surface, eliciting the slip between the sensor attached to the gripper
tip and the object. In each trial, the gripping force is selected from a range
of 8 N to 30 N. The robot movement can be either translational, rotational or
a combination of translational and rotational. The velocity ($v$) and
acceleration ($a$) of the robot movement have three different levels: low
($v=4\text{ mm.s}^{-1}$, $a=10\text{ mm.s}^{-2}$), medium ($v=10\text{
mm.s}^{-1}$, $a=50\text{ mm.s}^{-2}$), and high ($v=40\text{ mm.s}^{-1}$,
$a=100\text{ mm.s}^{-2}$). All robot movements were performed using the built-
in movel function of the UR script. The tool center position and orientation
are obtained using the built-in getl function of the UR robot. This function
employs forward kinematics calculations based on the read joint angles.

In accordance with the offline evaluation, control trails are also conducted
here for each $v$ and $a$ combination and movement type. The purpose is to
validate that the identified behavior is indeed the incipient slip, rather the
event with similar pattern like the stop event we mentioned above. The control
data involves lifting the robot arm while maintaining a secure grip using a
pre-determined grip force that is sufficient to prevent any slippage. As a
result, when lifting an object, the pillars in contact undergo downward
deformation due to the force of gravity; subsequently, once the object is
securely held by the gripper and remains relatively motionless, these pillars
will remain stationary. Here, for the sake of convenient explanation, we will
also refer to this event as stop, and we label the sequence as other. To
ensure a fair experiment, we add extra weight to lightweight objects to
enhance their downward motion when being lifted, aiming to make the pattern of
the output data sequence more like a slip event. In total, our experiment
consisted of 216 trials, including 162 sequences of slip event (6 objects
$\times$ 3 movements $\times$ 3 forces $\times$ 3 velocity/acceleration
combinations) and 54 sequences of stop event (6 objects $\times$ 3 movements
$\times$ 1 force $\times$ 3 velocity/acceleration combinations).

Fig. 6 illustrates the final validation results. Fig. 6 shows a confusion
matrix, highlighting the high success rate ($\sim$96.8%) of our method in
detecting incipient slip and its ability to differentiate between slip and
stop events. Fig. 6 demonstrates that our algorithm can detect incipient slip
almost immediately upon the initiation of the movement that induces slip, with
a normalized displacement $D_{norm}$ range of 0.2 - 0.4, within which the
incipient slip can be detected (refer to the caption for the definition of
$D_{norm}$). These results provide comprehensive validation of the
effectiveness of our approach in detecting incipient slipping in real-world
gripping tasks.

####  IV-B2 Ablation study

This study aims to showcase the effectiveness of our advanced augmentation
method in bridging the domain gap between the idealized data collected with
the hexapod and more realistic data encountered with the robotic gripper. To
accomplish this, we employed the model training approach described in Section
IV-A. However, instead of splitting the data into separate train and test
sets, we trained the model using the entire dataset here, given the different
objective. Subsequently, we conducted online gripping experiments, as
described in Section IV-B1, using this trained model. Our findings, as
illustrated in Fig. 7, indicate that the model trained without our advanced
augmentation method exhibits a notably high false positive rate in the
subsequent online gripping task when compared to the results shown in Fig. 6
where the model was trained using our advanced augmentation method. In other
words, the model trained without our advanced augmentation is unable to
effectively distinguish patterns between slip and stop events. As a result, it
incorrectly detects incipient slip in many stop events.

Figure 7: The ablation study results are to be compared with those presented
in Fig. 6. The model here is trained without using data augmentation of the
training data set leading to an increased number of misclassifications.

###  IV-C Grasp stabilization after incipient slip detection

This experiment aims to show the benefit of using our incipient slip detection
method in practical gripping tasks. This involve lifting the robot arm while
gripping the object with a pre-determined small force to ensure that slip
occurs. We applied our incipient slip detection method and adjusted the grip
when incipient slip was first detected to prevent the object from slipping
further. In this experiment, we simulate two common scenarios that can trigger
slips. The first involves gripping an object at its center of gravity with
insufficient force and lifting it, causing a translational slip between the
gripper and the object. The second involves gripping an object away from its
center of gravity and lifting it, where rotational slip is likely to occur. We
implemented a simple grip force adaptation that responds to incipient slip
detection as follows: if incipient slip is detected, the robot immediately
stops, and the gripper applies a pre-determined secure force to the object.
The objects used in the experiment are the same as those shown in Fig. 5. The
experiment was conducted 36 times (6 objects $\times$ 2 scenarios (translation
or rotation) $\times$ 3 repetitions). We fix ArUco markers on the objects and
gripper and use Python OpenCV to track the positions and orientations of all.

We report the results in Table I, which demonstrate the quickly and effective
detection of incipient slip using our algorithm. On average, our algorithm can
timely detect incipient slip and prevent the object from slipping when the
relative translation between the object and the gripper reaches 2.5 mm and the
relative rotation reaches 1.9 ∘. Our algorithm showcases its ability to
facilitate timely corrective action, preventing object falls; a demonstration
video can be seen at our project website given in the abstract.

Table I: The translation and angular displacement between the object and the gripper after being lifted in our online gripping and lifting experiment, utilizing our incipient slip detection and grip correction methods. Object | Translation (mm) | Rotation (∘)
---|---|---
Hard paper box | $4.5\pm 0.6$ | $1.9\pm 0.6$
Wooden box | $2.4\pm 0.8$ | $1.5\pm 0.9$
Shampoo bottle | $2.8\pm 1.2$ | $2.1\pm 0.4$
Coffee jar | $2.1\pm 0.5$ | $1.6\pm 0.2$
Pringles can | $1.4\pm 0.4$ | $1.1\pm 0.3$
Thermos flask | $2.1\pm 0.6$ | $1.2\pm 0.4$
Average | $2.5\pm 0.7$ | $1.9\pm 0.5$

##  V Discussion

Our developed algorithm enables the NN to effectively learn the incipient slip
pattern from offline data and demonstrates high accuracy in both offline and
online test sets. Furthermore, our algorithm enhances the security of robotic
gripping.

Compared to previous related works [6, 13], our algorithm offers several
advantages. Firstly, our incipient slip detection algorithm incorporates a
data-driven learning-based approach, minimizing the need for extensive human
involvement in investigating the complex patterns of incipient slip. Secondly,
the improved robustness of our algorithm enables the NN to effectively adapt
to diverse domains with various types of PapillArray sensors and robotic
gripping systems, despite being trained solely on data lacking heterogeneity.
Therefore, our algorithm is more practical and possesses greater potential for
maximizing the utilization of valuable tactile data in real-world scenarios.
Thirdly, our algorithm has the ability to distinguish between incipient slip
and a closely related tactile pattern that we refer to as a stop event.
Notably, previous related work [6, 13, 7, 19] has not adequately considered or
addressed the stop event; however, our investigation has revealed the
importance of including stop events when developing incipient slip detection
algorithms due to their similar patterns but entirely different consequences.

There are limitations to our work that need consideration. Firstly, the
incipient slip detection could be improved by transitioning from a binary
signal to a continuous warning signal. For instance, if incipient slip is
detected in a small portion of the contact surface, the remaining area may
still possess sufficient fraction to prevent significant slippage. In such
cases, the warning level of incipient slip is low and corrective actions may
not be necessary. Conversely, if a significant portion of the contact surface
exhibits incipient slip, the warning level should escalate and it becomes
important to for appropriate corrective actions. Moreover, our current choice
of force adaptation method for reacting to incipient slip falls short when
compared to the state-of-the-art gripping control work [6]. However, it is
important to note that force adjustment is not the primary focus of our
research in this paper, which is focused on improving the incipient slip
detection. In future work, we will develop a more sophisticated force
adaptation technique that incorporates our incipient slip detection method.

##  VI Conclusion

In conclusion, this paper presents an incipient slip detection method that
employs deep learning and several data augmentation techniques to improve the
robustness of the trained NN. Our method is highly effective and reaches the
state-of-art performance, it enable a single pre-trained NN model to be
applied across various domains and tasks. In addition, our method has the
potential to be extended to other approaches that use compliant tactile
sensors.

To train the NN parameters, we use stochastic gradient descent with a momentum
of 0.95 and a learning rate of $10^{-3}$, with a batch size of 512. We also
incorporate a weight decay of $10^{-3}$ using $L_{2}$ regularization during
training. The encoder NN consists of one hidden layer with 1024 units, and the
output dimension is 128. The GRU cell has a hidden layer dimension of 128. The
predictor network comprises two hidden layers with 256 and 128 units,
respectively. To all hidden layers, we apply rectified non-linearity [25] and
batch normalization [26].

We implement our NN using PyTorch (Version 1.12.1, Meta, USA). All our
experiments are conducted on a PC with an Intel 7-10875H CPU and an NVIDIA
2060 GPU. During the online evaluation stage, e utilise ROS [27] to facilitate
communication between various components in our system.

## References

  * [1] R. S. Johansson and G. Westling, “Roles of glabrous skin receptors and sensorimotor memory in automatic control of precision grip when lifting rougher or more slippery objects,” Experimental Brain Research, vol. 56, pp. 550–564, 1984.
  * [2] A.-S. Augurelle, A. M. Smith, T. Lejeune, and J.-L. Thonnard, “Importance of cutaneous feedback in maintaining a secure grip during manipulation of hand-held objects,” Journal of Neurophysiology, vol. 89, no. 2, pp. 665–671, 2003.
  * [3] A. B. Vallbo, R. S. Johansson, et al., “Properties of cutaneous mechanoreceptors in the human hand related to touch sensation,” Hum Neurobiol, vol. 3, no. 1, pp. 3–14, 1984.
  * [4] W. Chen, H. Khamis, I. Birznieks, N. F. Lepora, and S. J. Redmond, “Tactile sensors for friction estimation and incipient slip detection—toward dexterous robotic manipulation: A review,” IEEE Sensors Journal, vol. 18, no. 22, pp. 9049–9064, 2018.
  * [5] B. P. Delhaye, E. Jarocka, A. Barrea, J.-L. Thonnard, B. Edin, and P. Lefevre, “High-resolution imaging of skin deformation shows that afferents from human fingertips signal slip onset,” Elife, vol. 10, p. e64679, 2021.
  * [6] H. Khamis, B. Xia, and S. J. Redmond, “Real-time friction estimation for grip force control,” in 2021 IEEE International Conference on Robotics and Automation (ICRA), pp. 1608–1614, IEEE, 2021.
  * [7] J. W. James, N. Pestell, and N. F. Lepora, “Slip detection with a biomimetic tactile sensor,” IEEE Robotics and Automation Letters, vol. 3, no. 4, pp. 3340–3346, 2018.
  * [8] M. Schöpfer, C. Schürmann, M. Pardowitz, and H. Ritter, “Using a piezo-resistive tactile sensor for detection of incipient slippage,” in ISR 2010 (41st International Symposium on Robotics) and ROBOTIK 2010 (6th German Conference on Robotics), pp. 1–7, VDE, 2010.
  * [9] S. du Bois de Dunilac, D. Córdova Bulens, P. Lefèvre, S. J. Redmond, and B. P. Delhaye, “Biomechanics of the finger pad in response to torsion,” Journal of the Royal Society Interface, vol. 20, no. 201, p. 20220809, 2023\.
  * [10] B. Delhaye, P. Lefevre, and J.-L. Thonnard, “Dynamics of fingertip contact during the onset of tangential slip,” Journal of The Royal Society Interface, vol. 11, no. 100, p. 20140698, 2014.
  * [11] Z. Su, K. Hausman, Y. Chebotar, A. Molchanov, G. E. Loeb, G. S. Sukhatme, and S. Schaal, “Force estimation and slip detection/classification for grip control using a biomimetic tactile sensor,” in 2015 IEEE-RAS 15th International Conference on Humanoid Robots (Humanoids), pp. 297–303, IEEE, 2015\.
  * [12] J. W. James, S. J. Redmond, and N. F. Lepora, “A biomimetic tactile fingerprint induces incipient slip,” in 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 9833–9839, IEEE, 2020\.
  * [13] P. M. Ulloa, D. C. Bulens, and S. J. Redmond, “Incipient slip detection for rectilinear movements using the papillarray tactile sensor,” in 2022 IEEE Sensors, pp. 1–4, 2022.
  * [14] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, “Learning internal representations by error propagation,” tech. rep., California Univ San Diego La Jolla Inst for Cognitive Science, 1985.
  * [15] B. E. Boser, I. M. Guyon, and V. N. Vapnik, “A training algorithm for optimal margin classifiers,” in Proceedings of the Fifth Annual Workshop on Computational Learning Theory, pp. 144–152, 1992.
  * [16] C. Chorley, C. Melhuish, T. Pipe, and J. Rossiter, “Development of a tactile sensor based on biologically inspired edge encoding,” in 2009 International Conference on Advanced Robotics, pp. 1–6, IEEE, 2009.
  * [17] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-based learning applied to document recognition,” Proceedings of the IEEE, vol. 86, no. 11, pp. 2278–2324, 1998.
  * [18] W. Yuan, S. Dong, and E. H. Adelson, “Gelsight: High-resolution robot tactile sensors for estimating geometry and force,” Sensors, vol. 17, no. 12, p. 2762, 2017.
  * [19] W. Yuan, R. Li, M. A. Srinivasan, and E. H. Adelson, “Measurement of shear and slip with a gelsight tactile sensor,” in 2015 IEEE International Conference on Robotics and Automation (ICRA), pp. 304–311, IEEE, 2015.
  * [20] S. Dong, D. Ma, E. Donlon, and A. Rodriguez, “Maintaining grasps within slipping bounds by monitoring incipient slip,” in 2019 International Conference on Robotics and Automation (ICRA), pp. 3818–3824, IEEE, 2019.
  * [21] H. Khamis, B. Xia, and S. J. Redmond, “A novel optical 3d force and displacement sensor–towards instrumenting the papillarray tactile sensor,” Sensors and Actuators A: Physical, vol. 291, pp. 174–187, 2019.
  * [22] H. Khamis, R. I. Albero, M. Salerno, A. S. Idil, A. Loizou, and S. J. Redmond, “Papillarray: An incipient slip sensor for dexterous robotic or prosthetic manipulation–design and prototype validation,” Sensors and Actuators A: Physical, vol. 270, pp. 195–204, 2018.
  * [23] R. S. Johansson and A. B. Vallbo, “Tactile sensibility in the human hand: relative and absolute densities of four types of mechanoreceptive units in glabrous skin.,” The Journal of Physiology, vol. 286, no. 1, pp. 283–300, 1979.
  * [24] J. Chung, C. Gulcehre, K. Cho, and Y. Bengio, “Empirical evaluation of gated recurrent neural networks on sequence modeling,” arXiv preprint arXiv:1412.3555, 2014.
  * [25] X. Glorot, A. Bordes, and Y. Bengio, “Deep sparse rectifier neural networks,” in Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics, pp. 315–323, JMLR Workshop and Conference Proceedings, 2011.
  * [26] S. Ioffe and C. Szegedy, “Batch normalization: Accelerating deep network training by reducing internal covariate shift,” in International Conference on Machine Learning, pp. 448–456, PMLR, 2015.
  * [27] M. Quigley, K. Conley, B. Gerkey, J. Faust, T. Foote, J. Leibs, R. Wheeler, A. Y. Ng, et al., “Ros: an open-source robot operating system,” in ICRA workshop on open source software, vol. 3, p. 5, Kobe, Japan, 2009.