# Cold hybrid electrical-optical ion trap

Jin-Ming Cui<EMAIL_ADDRESS>CAS Key Laboratory of Quantum Information,
University of Science and Technology of China, Hefei 230026, China  CAS Center
for Excellence in Quantum Information and Quantum Physics, University of
Science and Technology of China, Hefei 230026, China  Hefei National
Laboratory, University of Science and Technology of China, Hefei 230088, China
Shi-Jia Sun  CAS Key Laboratory of Quantum Information, University of Science
and Technology of China, Hefei 230026, China  CAS Center for Excellence in
Quantum Information and Quantum Physics, University of Science and Technology
of China, Hefei 230026, China  Xi-Wang Luo<EMAIL_ADDRESS>CAS Key
Laboratory of Quantum Information, University of Science and Technology of
China, Hefei 230026, China  CAS Center for Excellence in Quantum Information
and Quantum Physics, University of Science and Technology of China, Hefei
230026, China  Yun-Feng Huang  CAS Key Laboratory of Quantum Information,
University of Science and Technology of China, Hefei 230026, China  CAS Center
for Excellence in Quantum Information and Quantum Physics, University of
Science and Technology of China, Hefei 230026, China  Hefei National
Laboratory, University of Science and Technology of China, Hefei 230088, China
Chuan-Feng Li<EMAIL_ADDRESS>CAS Key Laboratory of Quantum Information,
University of Science and Technology of China, Hefei 230026, China  CAS Center
for Excellence in Quantum Information and Quantum Physics, University of
Science and Technology of China, Hefei 230026, China  Hefei National
Laboratory, University of Science and Technology of China, Hefei 230088, China
Guang-Can Guo  CAS Key Laboratory of Quantum Information, University of
Science and Technology of China, Hefei 230026, China  CAS Center for
Excellence in Quantum Information and Quantum Physics, University of Science
and Technology of China, Hefei 230026, China

###### Abstract

Advances in research such as quantum information and quantum chemistry require
subtle methods for trapping particles (including ions, neutral atoms,
molecules, etc.). Here we propose a hybrid ion trapping method by combining a
Paul trap with optical tweezers. The trap combines the advances of the deep-
potential feature for the Paul trap and the micromotion-free feature for the
optical dipole trap. By modulating the optical-dipole trap synchronously with
the radio frequency voltage of the Paul trap, the alternating electrical force
in the trap center is fully counteracted, and the micromotion temperature of a
cold trapped ion can reach the order of nK while the trap depth is beyond
300K. These features will enable cold collisions between an ion and an atom in
the $s$-wave regime and stably trap the produced molecular ion in the cold
hybrid system. This will provide a unique platform for probing the
interactions between the ions and the surrounding neutral particles and enable
the investigation of new reaction pathways and reaction products in the cold
regime.

_Introduction.—_ Cold chemistry near absolute zero reveals quantum process in
the reaction process, which attracts interests in both chemical and physical
fields Heazlewood and Softley (2021); He  _et al._ (2020); Hu  _et al._
(2019); Zhang  _et al._ (2020); Hood  _et al._ (2020); Liu and Ni (2022); Yang
_et al._ (2022); Su  _et al._ (2022); Cao  _et al._ (2023); Cheuk  _et al._
(2020); Xie  _et al._ (2020); de Jongh  _et al._ (2020). In the last decade,
an emerging area of cold hybrid ion-atom systems has been developed as a new
platform for fundamental research, which is dedicated to controlling the cold
collision and quantum chemistry of ions and atoms Härter and Denschlag (2014);
Tomza  _et al._ (2019); Zuber  _et al._ (2022); Zou  _et al._ (2023);
Weckesser  _et al._ (2021); Hirzler  _et al._ (2022); Katz  _et al._ (2022).
Besides, binding the neutral atom and ion and then stably trapping the
produced molecular ion are prerequisites for exploring the quantum features of
the molecular ion. To investigate collision physics in the quantum regime of
the cold hybrid ion-atom system, reaching the single partial wave limit (i.e.,
the $s$-wave limit) is an important step Tomza  _et al._ (2019). Although the
trapped ion can be cooled close to the quantum collision regime using buffer
atomic gases with a large ion-atom mass ratio Feldker  _et al._ (2020), the
micromotion of the Paul trap (PT) heats the system away from a stable bound
state Cetina  _et al._ (2012); Pinkas  _et al._ (2022); Hirzler  _et al._
(2023), which hinders progress in this area. Pioneering experiments using
optical dipole traps (ODTs) to capture ions were demonstrated Schneider  _et
al._ (2010); Huber  _et al._ (2014); Schmidt  _et al._ (2018, 2020), which is
a potential solution to the micromotion heating problem in the cold hybrid
ion-atom system. However, because the dipole potential is shallow, it is quite
challenging to keep the ion during the reaction and stably trap the produced
molecular ion.

In this work, we propose a novel hybrid electrical-optical ion trap that
combines the advantages of the deep-potential feature of a PT and the
micromotion-free feature of an ODT. The key idea is to counteract the
alternating electrical force around the center of a PT that is responsible for
the micromotion problem, which is achieved by introducing optical-dipole traps
that are modulated synchronously with the radio frequency (RF) voltage of the
PT. This hybrid electrical-optical trap corresponds to a deep PT with a
micromotion-free dip at the center, and this potential structure provides an
ideal platform for studying low-energy physics (such as cold collisions and
cold chemistry). Since the center dip provides a micromotion-free area for the
trapped ion to interact with a neutral atom, the produced molecular ion can be
stably kept in a deep background potential profile provide by the PT. To
analyze the performance of the trap, we consider counteracting $99\%$ of the
alternating force which is experimentally reasonable. The intrinsic
micromotion (iMM) energy can be suppressed by four orders of magnitude. Under
current voltage compensation techniques in trapped ion fields Gloger  _et al._
(2015), the DC stray field can be suppressed below 1 V/m, which leads to the
ultimate cold ion excess micromotion (eMM) temperature in the hybrid trap
below 1 nK. Under such cold ion temperatures, many ion-atom $s$-wave
scattering process can be investigated. Moreover, by simulating a cold
collision process between a $\mathrm{Yb^{+}}$ ion and a Rb atom, we find that
the hybrid trap effectively solves the micromotion heating problem, which
dissociates the ion-atom bound state in the cold collision process. The
proposed hybrid trap offers a simple yet powerful platform to study ion-atom
collisions in the quantum regime, which paves the way for exploring novel
reactions and manipulating the produced molecular ion in such cold hybrid
systems.

_Scheme and Results.—_ Our scheme employs modulating optical tweezers to
counteract the alternating force caused by the RF electrical field. If we
consider a hybrid trap consisting of a PT and $N$ modulating ODTs, the total
potential of the trap is

$\displaystyle U(\mathbf{r},t)=$ $\displaystyle
U_{E}(\mathbf{r})+\widetilde{U}_{E}(\mathbf{r})\cos(\Omega t)$
$\displaystyle+\sum_{i}^{N}U_{D,i}(\mathbf{r})\left[1+\eta_{i}\cos(\Omega
t)\right],$ (1)

where $U_{E}(\mathbf{r})$ and $\widetilde{U}_{E}(\mathbf{r})$ are the
potentials generated by the DC field and RF field of the PT respectively,
$U_{D,i}(\mathbf{r})$ is the potential of $i$’th ODT, and $\eta_{i}$ is the
modulation depth of $i$’th ODT with $\left|\eta_{i}\right|\leq 1$. Supposing
the centers of all traps overlap each other for simplicity, the potentials
around the trap center can be expanded as harmonic forms as

$\displaystyle U_{E}(\mathbf{r})$
$\displaystyle=k_{0}^{x}x^{2}+k_{0}^{y}y^{2}+k_{0}^{z}z^{2}+C_{E},$
$\displaystyle\widetilde{U}_{E}(\mathbf{r})$
$\displaystyle=k_{1}^{x}x^{2}+k_{1}^{y}y^{2}+k_{1}^{z}z^{2}+\widetilde{C}_{E},$
$\displaystyle U_{D,i}(\mathbf{r})$
$\displaystyle=K_{i}^{x}x^{2}+K_{i}^{y}y^{2}+K_{i}^{z}z^{2}+C_{D,i},$

where $C_{E}$, $\widetilde{C}_{E}$ and $C_{D,i}$ are constants.
$U_{E}(\mathbf{r})$ and $\widetilde{U}_{E}(\mathbf{r})$ generated by
electrical fields have to fulfill the Laplace equation $\Delta
U(\mathbf{r})=0$, which leads to restrictions in the geometric factors,

$\displaystyle k_{0}^{x}+k_{0}^{y}+k_{0}^{z}$ $\displaystyle=0,$ (2)
$\displaystyle k_{1}^{x}+k_{1}^{y}+k_{1}^{z}$ $\displaystyle=0,$ (3)

while $U_{D,i}(\mathbf{r})$ is not restricted by this condition. The
alternating force strength of the hybrid trap is
$\mathbf{\widetilde{F}(\mathbf{r})}=\nabla\left(\widetilde{U}_{E}+\sum\eta_{i}U_{D,i}\right)$,
so eliminating the force to zero leads to the condition

$k_{1}^{\gamma}+\sum_{i}^{N}\eta_{i}K_{i}^{\gamma}=0,$ (4)

for all $\gamma\in\left\\{x,y,z\right\\}$, which can be achieved by setting
the modulation depths and geometric shapes of ODTs.

Figure 1:  Scheme of a hybrid electrical-optical ion trap. (a) The trap
consists of a Paul trap (e.g., a four-rod trap or a blade trap) and two
modulating optical dipole traps (ODTs) with different oriented astigmatic
spots. (b) The phase relation between the RF voltage of the Paul trap and the
optical intensities of the two ODTs. (c, d) two feasible realizations of the
hybrid trap.

Supposing we use astigmatic Gaussian spots to generate the dipole traps
required. The intensity profile of a focused astigmatic Gaussian beam along
axis $c$ in coordinate $(\hat{a},\hat{b},\hat{c})$ can be written as

$I(\mathbf{r})=I_{0}\frac{w_{0a}}{w_{a}(c)}\frac{w_{0b}}{w_{b}(c)}\exp\left\\{-2\left[\frac{a^{2}}{w_{a}(c)^{2}}+\frac{b^{2}}{w_{b}(c)^{2}}\right]\right\\},$
(5)

where $w_{a/b}(c)=w_{0,a/b}\sqrt{1+(c/z_{Ra/b})^{2}}$ are the beam waists for
$a$ and $b$ axes, i.e., short and long axes respectively, $z_{Ra/b}=\pi
w_{0,a/b}^{2}/\lambda$ are the Rayleigh ranges and $\lambda$ is the
wavelength. As the potential of ODT is proportional its intensity
$U_{D,i}(\mathbf{r})\propto I_{D,i}(\mathbf{r})$, by expanding it near the
zero point, the potential coefficients are

$\displaystyle K_{i}^{a}$
$\displaystyle=\frac{2U_{0i}}{w_{0,ai}^{2}},K_{i}^{b}=\frac{2U_{0i}}{w_{0,bi}^{2}},$
$\displaystyle K_{i}^{c}$
$\displaystyle=\frac{U_{0i}\lambda^{2}}{2\pi^{2}}\left(\frac{1}{w_{0,ai}^{4}}+\frac{1}{w_{0,bi}^{4}}\right),$
(6)

where $U_{0i}=\frac{\pi
c^{2}\Gamma}{2\omega_{0}^{3}}(\frac{2}{\Delta_{2}}+\frac{1}{\Delta_{1}})I_{0i}$
Grimm  _et al._ (2000). Finding a proper set of
$\left\\{\eta_{i}U_{0i},\,w_{0,ai},\,w_{0,bi}\right\\}$ to satisfy Eq. 4 and
Eq. 6 will eliminate the alternation force at the trap center .

Figure 1(c,d) shows two feasible configurations to demonstrate the scheme in a
symmetric linear PT with four rods or blades Berkeland  _et al._ (1998);
Drewsen and Brøner (2000). For simplicity, we take Fig. 1(c) for example,
where two optical tweezers with focused astigmatic Gaussian spots are applied
on the trapped ion with similar beam waists, i. e., $w_{0,a1}=w_{0,a2}=w_{a}$,
$w_{0,b1}=w_{0,b2}=w_{b}$, $U_{01}=U_{02}=U_{0}$, and the two trap long axes
are perpendicular to each other, i.e.
$\hat{a}_{1}\perp\hat{a}_{2}\parallel\hat{y}$,
$\hat{b}_{1}\perp\hat{b}_{2}\parallel\hat{x}$. The symmetric linear PT
condition requires $k_{1}^{x}=-k_{1}^{y}=k_{1}$, $k_{1}^{z}=0$. Suppose
$\eta_{1}=-\eta_{2}$=1, which indicates modulations between two ODTs have a
$\pi$ phase difference (shown in Fig. 1(b)), then the configurations of the
ODTs can be

$\displaystyle k_{1}$ $\displaystyle=k_{a}-k_{b},$ $\displaystyle K_{1}^{x}$
$\displaystyle=K_{2}^{y}=k_{a},$ $\displaystyle K_{1}^{y}$
$\displaystyle=K_{2}^{x}=k_{b},$ (7)

where $k_{a,b}=2U_{0}/w_{a,b}^{2}$.

Figure 2:  Pseudo potential of a hybrid trap for $\mathrm{{}^{171}Yb^{+}}$.
(a) The total pseudo potential of the trap, indicating a deep potential. (b,c)
The constant part and the alternating part of the pseudo potential near the
trap center, respectively. The alternating pseudo potential is flat near the
center, which indicates the alternating force is zero. The simulation
condition is $\Omega=2\pi\times 2$ MHz, $V_{RF}=10$ V, $\kappa=0.6$,
$V_{DC}=1$ V, $R=0.5$ mm, $Z_{0}=1$ mm, $\omega_{x}=\omega_{y}=2\pi\times 180$
kHz, $\omega_{z}=2\pi\times 131$ kHz for the Paul trap, laser wavelength at
554.6 nm, $w_{a}=w_{b}/3=478$ nm can be generated through an objective with a
numerical aperture (NA) of 0.6, $U_{0}/h=-622$ MHz for optical dipole traps.

The pseudo potential of the hybrid trap can be decomposed into a direct
potential (DP) and an alternating potential (AP)

$\varPsi(\mathbf{r})=\Phi_{0}(\mathbf{r})+\frac{1}{4m\Omega^{2}}\left\|\nabla\Phi_{1}(\mathbf{r})\right\|^{2},$
(8)

where $\Phi_{0}(\mathbf{r})=U_{E}(\mathbf{r})+\sum_{i}^{N}U_{D,i}(\mathbf{r})$
and
$\Phi_{1}(\mathbf{r})=\widetilde{U}_{E}(\mathbf{r})+\sum_{i}^{N}\eta_{i}U_{D,i}(\mathbf{r})$
are the direct and alternating potentials respectively. For a symmetric linear
PT Berkeland  _et al._ (1998), $k_{0}^{x}=k_{0}^{y}=-2k_{0}^{z}=-2\kappa
QV_{DC}/Z_{0}^{2}$, $k_{1}^{x}=-k_{1}^{y}=QV_{RF}/2R^{2}$, $k_{1}^{z}=0$,
where $\kappa$ is a geometric factor of the trap, $Q$ is the charge of the ion
, $V_{DC}$ and $V_{RF}$ are the DC and RF voltages applied to the trap, $R$ is
the perpendicular distance from the trap axis to the trap electrodes, $Z_{0}$
is the distance from the trap center to the endcap electrodes. The pseudo
potential can be calculated by substituting Eq. 5-7 into Eq. 8. When the
alternating force caused by the RF electrical field is completely compensated
by the modulating ODTs, which we call the zero alternating trapping condition
(ZATC), we can derive
$U_{0}=QV_{RF}w_{a}^{2}w_{b}^{2}/4R^{2}(w_{a}^{2}-w_{b}^{2})$. To make the
condition realistic, a low RF voltage of PT and small spots of ODTs should be
used to minimize the optical power of the ODT. Consequently, a lower RF
frequency $\Omega$ is preferred to make a deeper pseudo potential under the
low RF voltage. A numerical simulation result is show in Fig. 2. In a large
range, it acts as a PT with deep potential as Fig. 2(a). At the center of the
trap, it acts as a tight ODT as Fig. 2(b). The alternating pseudo potential in
Fig. 2(c) is flat, which indicates the alternating force is zero. With the
simulation result, the required laser power is 0.7 W at 554.6 nm, if we use an
objective with NA=0.6 to generate the required ODTs.

Figure 3:  Ion motion simulation for trapping temperature in a stray
electrical field (SEF). (a) Excess micromotion (eMM) temperature of a trapped
ion along with different SEFs. The hybrid trap provides a lower temperature
than a PT under the same SEF. When the hybrid trap works at zero alternating
trapping condition (ZATC), the ion’s eMM temperature reaches nK under a SEF of
1 V/m. The direction of the SEF for simulation is (1,1,1). (b) Ion’s eMM
temperature simulation under experimental imperfections, $\delta x$:
misalignment between Paul trap and the ODT, $\Delta q$: uncompensated
alternating force factor. The figure is calculated under a SEF of 1 V/m along
the $x$ axis. Hybrid trap and Paul trap parameters are the same as those in
Fig. 2.

Next, we will discuss the cold trapped ion temperature limit in the hybrid
trap. Supposing the ion oscillates near $\mathbf{u}_{0}$ in the hybrid trap
under a stray electrical field $\mathbf{E}$, and has been cooled to a low
temperature $\left|\mathbf{u}-\mathbf{u}_{0}\right|\ll w_{a}$, then the motion
equation can be approximated as a Mathieu equation near $\mathbf{u}_{0}$

$\ddot{u_{l}}+\left[a_{l}(\mathbf{u}_{0})+2q_{l}(\mathbf{u}_{0})\cos\left(\Omega
t\right)\right]\frac{\Omega^{2}}{4}u_{l}=\frac{Q\cdot E_{l}}{m},$ (9)

where
$a_{l}(\mathbf{u})=4\hat{e_{l}}\cdot\nabla\Phi_{0}(\mathbf{u})/m\Omega^{2}u_{l}$
and
$q_{l}(\mathbf{u})=2\hat{e_{l}}\cdot\nabla\Phi_{1}(\mathbf{u})/m\Omega^{2}u_{l}$.
The equilibrium position $\mathbf{u}_{0}$ can be numerically solved by
iterating

$u_{0,l}=\frac{4Q\cdot
E_{l}}{m\Omega^{2}\left[a_{l}(\mathbf{u})+q_{l}^{2}(\mathbf{u})/2\right]}.$
(10)

The final temperature caused by excess micromotion is

$T_{\mathrm{eMM}}=\frac{m\Omega^{2}}{16k_{B}}\sum_{l}q_{l}^{2}(\mathbf{u}_{0})u_{0,l}^{2}.$
(11)

Using the above numerical method, we can study the ion eMM temperature under a
large electrical field. Fig. 3(a) shows the comparison of trapped ion eMM
temperature in three different traps under a stray electrical field along the
(1,1,1) direction. The stray field will shift the trap center of the PT part
of the hybrid trap while the ODT is not affected, so the total pseudo
potential is a double well under a small stray field. The ion stays in the
tight ODT under a small stray field, and the ion’s eMM temperature will be low
even without modulating the ODTs. By increasing the stray field, the ion will
eventually be pushed from the ODT to PT, so the ion’s eMM temperature of the
hybrid trap is close to that of the PT, when the stray field is larger than 5
V/m.

Table 1:  ODT laser wavelengths of hybrid ion traps for different ion-atom systems. Ion | Atom |  $\lambda_{\mathrm{ODT}}$ (nm) |  $P$ (W) |  $E_{k}$/$k_{B}$ (nK)
---|---|---|---|---
$\mathrm{Yb^{+}}$ | Rb | 420.6 | 0.24 | 44.7
$\mathrm{Yb^{+}}$ | Rb | 422.3 | 0.25 | 44.7
$\mathrm{Yb^{+}}$ | Rb | 787.4 | 1.48 | 44.7
$\mathrm{Yb^{+}}$ | Yb | 554.63 | 0.70 | 44.7
$\mathrm{Ca^{+}}$ | Na | 589.46 | 0.48 | 1370.7
$\mathrm{Ca^{+}}$ | Li | 670.97 | 1.09 | 10601.5
$\mathrm{Sr^{+}}$ | Rb | 787.41 | 0.94 | 77.8

In order to further insight the ion eMM temperature relation with experiment
parameters (e. g., misalignment between PT and the ODT $\delta x$, stray
electrical field $E$, residual alternating force (RAF) parameter $\Delta
q=4(k_{1}+k_{a}-k_{b})/m\Omega^{2}$ ), we study one dimensional motion for
simplification. Considering the condition of small misaligned distance $\delta
x\ll w_{0}$, the approximated ion eMM temperature can be derived as 111See
Supplemental Material at [URL] for details of the ion dynamics in the hybrid
trap potential (Sec. I), details of ion eMM temperature with experiment
imperfections (Sec.II), ion-atom Collision simulation (Sec. III), and the work
down by RF field on the ion in hybrid trap (Sec. IV).

$T_{\mathrm{eMM}}\approx\frac{m\Omega^{2}}{16k_{B}}\left(\Delta
q+\frac{(\delta x-x_{0}){}^{2}}{x_{b}^{2}}\right)^{2}x_{0}^{2},$ (12)

where $x_{0}\thickapprox-Q\cdot E/(k_{0}+2k_{a}+2k_{b})$ is the ion position
offset pushed by the stray electrical field, and
$x_{b}^{2}=m\Omega^{2}/8(k_{b}/w_{b}^{2}-k_{a}/w_{a}^{2})$. Fig. 3(b) shows
that the ion eMM temperature can reach 100 nK level even with a misalignment
about 50 nm and 5% of RAF, under 1 V/m stray field along the $x$ axis. Several
ion imaging methods with a position accuracy below 10 nm have been
demonstrated Qian  _et al._ (2021); Blums  _et al._ (2018); Wong-Campos  _et
al._ (2016), which can be used to compensate the misalignment in future. The
ion iMM temperature and energy is proportional to ${\Delta q}^{2}$ Note (1),
which can also be significantly reduced by four orders of magnitude with $1\%$
RAF.

Figure 4:  Numerical simulations of a ion-atom (Rb and $\mathrm{Yb^{+}}$)
collision in a hybrid trap (a) and a Paul trap (b). For comparison, the hybrid
trap has the same electric fields with configurations as Fig. 2. Rb
approaching an $\mathrm{Yb^{+}}$ initially stationed in trap center, the first
hard-core collision time is at $t=0$. In Paul trap, the ion is heated by the
RF field and bounces out quickly, whereas in hybrid trap, the ion-atom can be
weakly bounded. The harmonic frequencies of two traps are
$\omega_{\mathrm{Hybrid}}=2\pi\times 0.84$ MHz and
$\omega_{\mathrm{Paul}}=2\pi\times 0.18$ MHz, leading to the corresponding
characteristic lengths $R_{\mathrm{Hybrid}}=33.4$ nm and
$R_{\mathrm{Paul}}=55.9$ nm, respectively. The simulation method is modified
from Ref. Cetina  _et al._ , 2012.

Finally, we will investigate the collision process between an ion and an atom
in the hybrid trap. To perform this study, we note that the laser wavelength
of the ODT should be specified for a given composition of ion and atom
elements. The alternating ODT should not affect the ground state of the atom,
otherwise, the alternating ODT will drive atom motion when the atom
approaching the ion. Tab. 1 lists several laser wavelengths for some
composition of ion-atom species, and the laser power with hybrid trap setting
as Fig. 2. We take Rb and $\mathrm{Yb^{+}}$ in this investigation following
the method presented in Ref. Cetina  _et al._ (2012), with details in Note
(1). When an ion moves in a PT without any disturbance, the RF field does zero
work on the ion for an RF period. However, in the cold hybrid ion-atom system,
the polarized atom will exert an attraction force on the ion (by attraction
potential $V(r)=-C_{4}/r^{4}$ ) and thus disturb the conserved cycle of energy
transfer, and finally the ion-atom system gains net energy from the RF field,
which increases the temperature of the system. A characteristic length scale
$R=\sqrt[6]{2C_{4}/m_{i}\omega^{2}}$ is preferred to be used (when the
interaction potential $V(r)$ is equal in magnitude to the trap harmonic
potential, where $\omega=\Omega\sqrt{a+q^{2}/2}/2$ is the secular frequency of
the ion near the trap center). In a PT (Fig. 4(b)), the ion gains enough
energy from the RF field when the atom approaches, and after one collision,
the atom escapes. In the hybrid trap (Fig. 4(a)), since the alternating force
induced by the RF field near the trap center is compensated, the ion-atom
system nearly undergoes an energy-conservative process and forms a weakly
bound ion-atom pair, keeping a low temperature. In the hybrid trap, the
corresponding characteristic energy scale is no longer restricted by the RF
heating but determined by the atomic gas temperature, which enables access to
the regimes of the $s$-wave limit or quantized ion motion in the cold hybrid
system. The work done by the RAF during a collision depends on the phase
$\phi$ at the time of the close-range collision, which can be written as Note
(1)

$W=\frac{1}{2}\mu\Omega^{2}R^{2}\sin\phi\cdot\left(\Delta
q\eta_{1}+\frac{R^{4}}{r_{0}^{2}x_{b}^{2}}\eta_{2}\right)$

where $\mu=m_{i}m_{a}/(m_{i}+m_{a})$ is the reduced mass,
$\eta_{1}=2/3\intop_{0}^{\Omega
t_{1}}(A-B\tau^{1/3})\cdot\tau^{-2/3}\sin\tau\cdot d\tau$ and
$\eta_{2}=2/3\int_{0}^{\Omega
t_{1}}(A-B\tau^{1/3})^{3}\cdot\tau^{-2/3}\sin\tau\cdot d\tau$ are integrates
with parameters $A=r_{c}r_{0}/R^{2}$, $B=\mu r_{0}^{2}/m_{i}R^{2}$, the ion’s
displacement from the trap center $r_{c}$ for the collision position and a
characteristic length $r_{0}=\sqrt[6]{18C/\mu\Omega{}^{2}}$. Fig. 4(a) shows
$r_{c}$ is distributed in $0\sim 20$ nm for multi hard-core collsions,
indicationg that the max work done by RAF for one collision is on the order of
tens of nK$\cdot k_{B}$.

_Conclusion and discussion.—_ In summary, we proposed a novel hybrid trap
consisting of a PT and modulated ODTs that takes advantage of the deep-
potential feature of the PT while eliminating the micromotion problem in the
cold hybrid ion-atom system. The trap can suppress the micromotion temperature
to the order of nK, offering a new route to reach the quantum collision in the
ion-atom mixture. Our proposed schemes are quite feasible since an ion trap
system with high optical accessibility has been demonstrated in He  _et al._
(2021). Compared with previous optical trapping methods for ions where the PT
is switched off Schneider  _et al._ (2010); Huber  _et al._ (2014); Schmidt
_et al._ (2018, 2020), our scheme can stably trap elementary and molecular
ions all through the cold ion-atom collision process due to the always-on PT,
which provides an ideal platform for studying cold ion-atom collision and
reaction and investigating the properties of the produced molecular ion.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China
(Grants No. 11774335, No. 11734015, and No. 12204455), the Key Research
Program of Frontier Sciences, CAS (Grant No. QYZDY-SSW-SLH003), the Science
Foundation of the CAS (Grant No. ZDRW-XH-2019-1), the Fundamental Research
Funds for the Central Universities (Grants No. WK2470000027, and No.
WK2470000028), Innovation Program for Quantum Science and Technology (Grants
No. 2021ZD0301604, and No. 2021ZD0301200). J.-M. Cui and S.-J. Sun contributed
equally to this work.

## References

  * Heazlewood and Softley (2021) B. R. Heazlewood and T. P. Softley, Nature Reviews Chemistry 5, 125 (2021).
  * He  _et al._ (2020) X. He, K. Wang, J. Zhuang, P. Xu, X. Gao, R. Guo, C. Sheng, M. Liu, J. Wang, J. Li, G. V. Shlyapnikov, and M. Zhan, Science 370, 331 (2020).
  * Hu  _et al._ (2019) M.-G. Hu, Y. Liu, D. D. Grimes, Y.-W. Lin, A. H. Gheorghe, R. Vexiau, N. Bouloufa-Maafa, O. Dulieu, T. Rosenband, and K.-K. Ni, Science 366, 1111 (2019).
  * Zhang  _et al._ (2020) J. T. Zhang, Y. Yu, W. B. Cairncross, K. Wang, L. R. B. Picard, J. D. Hood, Y.-W. Lin, J. M. Hutson, and K.-K. Ni, Physical Review Letters 124 (2020), 10.1103/physrevlett.124.253401.
  * Hood  _et al._ (2020) J. D. Hood, Y. Yu, Y. W. Lin, J. T. Zhang, K. Wang, L. R. Liu, B. Gao, and K. K. Ni, Physical Review Research 2 (2020), 10.1103/physrevresearch.2.023108.
  * Liu and Ni (2022) Y. Liu and K.-K. Ni, Annual Review of Physical Chemistry 73, 73 (2022).
  * Yang  _et al._ (2022) H. Yang, X.-Y. Wang, Z. Su, J. Cao, D.-C. Zhang, J. Rui, B. Zhao, C.-L. Bai, and J.-W. Pan, Nature 602, 229 (2022).
  * Su  _et al._ (2022) Z. Su, H. Yang, J. Cao, X.-Y. Wang, J. Rui, B. Zhao, and J.-W. Pan, Physical Review Letters 129 (2022), 10.1103/physrevlett.129.033401.
  * Cao  _et al._ (2023) J. Cao, H. Yang, Z. Su, X.-Y. Wang, J. Rui, B. Zhao, and J.-W. Pan, Physical Review A 107 (2023), 10.1103/physreva.107.013307.
  * Cheuk  _et al._ (2020) L. W. Cheuk, L. Anderegg, Y. Bao, S. Burchesky, S. S. Yu, W. Ketterle, K.-K. Ni, and J. M. Doyle, Physical Review Letters 125 (2020), 10.1103/physrevlett.125.043401.
  * Xie  _et al._ (2020) Y. Xie, H. Zhao, Y. Wang, Y. Huang, T. Wang, X. Xu, C. Xiao, Z. Sun, D. H. Zhang, and X. Yang, Science 368, 767 (2020).
  * de Jongh  _et al._ (2020) T. de Jongh, M. Besemer, Q. Shuai, T. Karman, A. van der Avoird, G. C. Groenenboom, and S. Y. T. van de Meerakker, Science 368, 626 (2020).
  * Härter and Denschlag (2014) A. Härter and J. H. Denschlag, Contemporary Physics 55, 33 (2014).
  * Tomza  _et al._ (2019) M. Tomza, K. Jachymski, R. Gerritsma, A. Negretti, T. Calarco, Z. Idziaszek, and P. S. Julienne, Reviews of Modern Physics 91, 035001 (2019).
  * Zuber  _et al._ (2022) N. Zuber, V. S. V. Anasuri, M. Berngruber, Y.-Q. Zou, F. Meinert, R. Löw, and T. Pfau, Nature 605, 453 (2022).
  * Zou  _et al._ (2023) Y.-Q. Zou, M. Berngruber, V. S. V. Anasuri, N. Zuber, F. Meinert, R. Löw, and T. Pfau, Physical Review Letters 130, 023002 (2023).
  * Weckesser  _et al._ (2021) P. Weckesser, F. Thielemann, D. Wiater, A. Wojciechowska, L. Karpa, K. Jachymski, M. Tomza, T. Walker, and T. Schaetz, Nature 600, 429 (2021).
  * Hirzler  _et al._ (2022) H. Hirzler, R. S. Lous, E. Trimby, J. Pérez-Ríos, A. Safavi-Naini, and R. Gerritsma, Physical Review Letters 128, 103401 (2022).
  * Katz  _et al._ (2022) O. Katz, M. Pinkas, N. Akerman, and R. Ozeri, Nature Physics 18, 533 (2022).
  * Feldker  _et al._ (2020) T. Feldker, H. Fürst, H. Hirzler, N. V. Ewald, M. Mazzanti, D. Wiater, M. Tomza, and R. Gerritsma, Nature Physics 16, 413 (2020).
  * Cetina  _et al._ (2012) M. Cetina, A. T. Grier, and V. Vuletić, Physical Review Letters 109, 253201 (2012).
  * Pinkas  _et al._ (2022) M. Pinkas, O. Katz, J. Wengrowicz, N. Akerman, and R. Ozeri, “Observation of trap-assisted formation of atom-ion bound states,”  (2022).
  * Hirzler  _et al._ (2023) H. Hirzler, E. Trimby, R. Gerritsma, A. Safavi-Naini, and J. Pérez-Ríos, Physical Review Letters 130, 143003 (2023).
  * Schneider  _et al._ (2010) C. Schneider, M. Enderlein, T. Huber, and T. Schaetz, Nature Photonics 4, 772 (2010).
  * Huber  _et al._ (2014) T. Huber, A. Lambrecht, J. Schmidt, L. Karpa, and T. Schaetz, Nature Communications 5, 5587 (2014).
  * Schmidt  _et al._ (2018) J. Schmidt, A. Lambrecht, P. Weckesser, M. Debatin, L. Karpa, and T. Schaetz, Physical Review X 8, 021028 (2018).
  * Schmidt  _et al._ (2020) J. Schmidt, P. Weckesser, F. Thielemann, T. Schaetz, and L. Karpa, Physical Review Letters 124, 053402 (2020).
  * Gloger  _et al._ (2015) T. F. Gloger, P. Kaufmann, D. Kaufmann, M. T. Baig, T. Collath, M. Johanning, and C. Wunderlich, Physical Review A 92, 043421 (2015).
  * Grimm  _et al._ (2000) R. Grimm, M. Weidemüller, and Y. B. Ovchinnikov, in  _Advances In Atomic, Molecular, and Optical Physics_, Vol. 42, edited by B. Bederson and H. Walther (Elsevier, 2000) pp. 95–170.
  * Berkeland  _et al._ (1998) D. J. Berkeland, J. D. Miller, J. C. Bergquist, W. M. Itano, and D. J. Wineland, Journal of Applied Physics 83, 5025 (1998).
  * Drewsen and Brøner (2000) M. Drewsen and A. Brøner, Physical Review A 62, 045401 (2000).
  * Note (1) See Supplemental Material at [URL] for details of the ion dynamics in the hybrid trap potential (Sec. I), details of ion eMM temperature with experiment imperfections (Sec.II), ion-atom Collision simulation (Sec. III), and the work down by RF field on the ion in hybrid trap (Sec. IV).
  * Qian  _et al._ (2021) Z.-H. Qian, J.-M. Cui, X.-W. Luo, Y.-X. Zheng, Y.-F. Huang, M.-Z. Ai, R. He, C.-F. Li, and G.-C. Guo, Physical Review Letters 127, 263603 (2021).
  * Blums  _et al._ (2018) V. Blums, M. Piotrowski, M. I. Hussain, B. G. Norton, S. C. Connell, S. Gensemer, M. Lobino, and E. W. Streed, Science Advances 4, eaao4453 (2018).
  * Wong-Campos  _et al._ (2016) J. D. Wong-Campos, K. G. Johnson, B. Neyenhuis, J. Mizrahi, and C. Monroe, Nature Photonics 10, 606 (2016).
  * He  _et al._ (2021) R. He, J.-M. Cui, R.-R. Li, Z.-H. Qian, Y. Chen, M.-Z. Ai, Y.-F. Huang, C.-F. Li, and G.-C. Guo, Review of Scientific Instruments 92, 073201 (2021).


# SIMAP: A simplicial-map layer for neural networks

Rocio Gonzalez-Diaz, Miguel A. Gutiérrez-Naranjo,
Eduardo Paluzo-Hidalgo  R. Gonzalez-Diaz is with the Department of Applied
Math I, School of Engineering, University of Sevilla, Seville, Spain.M.A.
Gutiérrez-Naranjo is with the Department of Computer Science and Artificial
Intelligence, School of Engineering, University of Sevilla, Seville, Spain.E.
Paluzo-Hidalgo is with the Department of Quantitative Methods, Universidad
Loyola Andalucía, Campus Sevilla, Dos Hermanas, Seville, Spain.
Authors are listed alphabetically.

###### Abstract

In this paper, we present SIMAP, a novel layer integrated into deep learning
models, aimed at enhancing the interpretability of the output. The SIMAP layer
is an enhanced version of Simplicial-Map Neural Networks (SMNNs), an
explainable neural network based on support sets and simplicial maps
(functions used in topology to transform shapes while preserving their
structural connectivity). The novelty of the methodology proposed in this
paper is two-fold: Firstly, SIMAP layers work in combination with other deep
learning architectures as an interpretable layer substituting classic dense
final layers. Secondly, unlike SMNNs, the support set is based on a fixed
maximal simplex, the barycentric subdivision being efficiently computed with a
matrix-based multiplication algorithm.

##  1 Introduction

In the current landscape of Artificial Intelligence, Deep Learning (DL) models
have evolved to possess enormous architectures characterized by a multitude of
parameters. These models, equipped with substantial computational capacity,
have found widespread applications in many real-world domains.

In response to the increasing complexity of DL architectures, it is imperative
to develop methods that improve the interpretability and explainability of the
models. In this way, addressing the opacity inherent in current DL
architectures is not only a technical challenge but also a key step toward
fostering trust and understanding in the deployment of artificial intelligence
systems. This paradigm shift emphasizes the importance of developing models
that not only deliver accurate predictions but also provide clear and
understandable reasoning for their decisions.

Interpretable layers can significantly contribute to reliable AI systems by
improving understanding to better assess whether the AI system is functioning
as intended. This term is typically used to describe a layer in a neural
network where the operations and transformations it performs are
understandable to humans. In this way, interpretability focuses on the
transparency of the process.

In the literature, several approaches can be found to this aim. One of the
first attempts to develop interpretable layers is [1], where each filter in
the convolutional layer is activated by a certain object part of a certain
category and remains inactivated on images of other categories. A recent
survey that gives an explanation of a wide range of interpretable neural
networks is [2]. In that paper, the authors divide interpretable neural
network methods into two primary categories: model decomposition neural
networks that involve the fusion of a conventional model-based method’s
interpretability with the neural network’s learning capability for a clearer
understanding of the neural network’s operations; and semantic interpretable
neural networks that derive their interpretability from a human-centric
perspective, using visualization and semantic understanding.

According to this taxonomy, the approach developed in this paper can be
considered as a model-decomposition neural network since we specifically
construct a neural network layer, the so-called SIMAP layer, based on a
simplicial map, which is a concept of algebraic topology that links two
triangulated spaces in such a way that incidence relations are preserved.
Furthermore, SIMAP layers can be considered as an evolution of Simplicial Map
Neural Networks (SMNNs), first defined in [3]. SMNNs are a type of feedforward
neural network that specifically uses simplicial maps. This makes them useful
as interpretable machine learning models. As already demonstrated in [4],
SMNNs are explainable neural network models, as they can rationalize the
outputs. Specifically, an SMNN provides a justification based on similarities
and dissimilarities of the data instance to be explained with the instances of
the training dataset that correspond to the vertices of the simplex that
contains it, after considering the dataset as a point cloud embedded in a
metric space. The drawbacks of the SMNN approach are mainly twofold. First,
the arbitrariness of the selection of a small subset of the training dataset
needed to compute the Delaunay triangulation and, second, the computation of
the Delaunay triangulation itself. Let us recall that the computational
complexity of the Delaunay triangulation increases significantly in higher
dimensions. In particular, for higher dimensions, the construction of a
Delaunay complex becomes challenging, even for datasets containing more than a
few hundred points. This phenomenon is associated with the curse of
dimensionality (see [5]).

The SIMAP layers overcome the mentioned drawbacks by first training them using
the barycentric coordinates of the input data with respect to a specific
simplex surrounding the dataset. In doing so, we avoid the need to extract a
small subset of the input data set and compute the Delaunay triangulation. In
addition, we demonstrate that the capacity of a SIMAP layer increases with
successive barycentric subdivisions of the simplex. We also prove that the
barycentric coordinates of the input data after the subdivision are obtained
just by matrix multiplications. In this way, the vertices of the simplex that
contain an input point are no longer part of the training set, improving, at
the same time, the interpretability of the model, as the entire process
becomes transparent and easily understandable to humans.

Another approach that establishes a topology-based layer is [6] where the
authors introduce a topology layer that computes persistent homology (a tool
that captures how topological features change over an increasing sequence of
complexes). In contrast to the SIMAP layer, the layer defined in [6] needs to
have on hand an ordering of the points of the data set to induce filtration,
i.e. an increasing sequence of triangulations. A different approach is the
recent field of topological deep learning (TDL), where the key idea is that
input data can have a rich structure, including graphs, hypergraphs or cell
complexes [7]. As explained in [7, page 60], “the architecture of a TDL
network can make it difficult to interpret the learnt representations and
understand how the network is making predictions”. On the contrary, so far,
SIMAP layers only support datasets embedded in a Euclidean space
$\mathbb{R}^{n}$ but, as we will see later, they are clearly interpretable.

This paper is organized as follows. In Section 2, we introduce the background
and some basic definitions. Section 3 is the paper’s main section and is
devoted to introducing the novel simplicial map layer based on barycentric
coordinates and barycentric subdivisions. In Section 4, some experiments are
performed to evaluate the effectiveness and reliability of the proposed
method. We end the paper with a section devoted to conclusions and future
work.

##  2 Background

In the following, we recall concepts such as simplicial complexes, simplicial
maps, barycentric coordinates, and barycentric subdivisions, which will be
later combined to improve the interpretability of deep neural network models.

###  2.1 Combinatorial topology

Let $V$ be a finite nonempty set whose elements are called vertices. A
simplicial complex $K$ with vertex set $V$, is a collection of nonempty
subsets of $V$ satisfying that:

  * •

for each vertex $v$ in $V$, the set $\\{v\\}$ is in $K$.

  * •

If $\tau$ is in $K$ and $\sigma\subset\tau$, then $\sigma$ is in $K$.

Every non-empty subset in $K$ is called a simplex; its dimension is given by
its cardinality minus one. Given two simplices $\sigma$ and $\tau$ in $K$, we
say that $\sigma$ is a face of $\tau$ whenever $\sigma\subset\tau$. The
simplices in $K$ that are not the face of other simplices of $K$ are called
maximal simplices.

Let us consider a simplicial complex $K$ with a set of vertices
$V\subset\mathbb{R}^{n}$ such that any simplex in $K$ is a set of (affinely
independent) points. Let $\sigma=\\{v^{0},v^{1},\dots,v^{d}\\}$ be a
$d$-simplex of $K$. Then, the geometric realization of $\sigma$, denoted by
$|\sigma|$, is defined by the convex hull of the vertices of $\sigma$:

$\begin{split}|\sigma|=\Big{\\{}\mbox{$x=\sum_{j=0}^{d}b_{j}(x)v^{j}$}\ \mid\
\mbox{$\sum_{j=0}^{d}b_{j}(x)=1$ and \hskip 34.14322pt }\\\ b_{j}(x)\geq 0\
\forall j\in\\{0,1,\dots,d\\}\Big{\\}}\,.\end{split}$

The vector $b(x)=(b_{0}(x),b_{1}(x),\dots,b_{d}(x))$ is called the barycentric
coordinates of $x$ with respect to $\sigma$. Furthermore, the geometric
realization or polytope of $K$, denoted by $|K|$, is the union of the
geometric realization of the simplices of $K$.

The barycentric subdivision of a simplicial complex $K$, denoted by
$\operatorname{Sd}K$, is a simplicial complex whose vertex set is the set of
the barycenters of all the simplices of $K$ (see Example 2.1). Recall that the
barycenter of a simplex  $\sigma=\\{v^{0},\dots,v^{d}\\}$ is the point

$\mbox{$\operatorname{bar}\sigma=\frac{1}{d+1}\sum_{i=0}^{d}v^{i}$}.$

Besides, for each dimension $i\geq 0$, $\mu$ is an $i$-simplex of
$\operatorname{Sd}K$ if and only if its vertices can be written as an ordered
set $\\{w^{0},w^{1},\dots,w^{i}\\}$ such that
$w^{k}=\operatorname{bar}\sigma_{k}$ for $k\in\\{0,1,\dots,i\\}$ and

$\sigma_{0}\subset\sigma_{1}\subset\cdots\subset\sigma_{i}\in K\ .$

The barycentric subdivision can be iterated as many times as needed. If $K$ is
subdivided $n$ times, then the corresponding simplicial complex is denoted by
$\operatorname{Sd}^{n}K$.

Let us consider now a simplicial complex $K$ with vertex set
$V=\\{v^{1},v^{2},\dots,v^{\beta}\\}\subset\mathbb{R}^{n}$. Consider also a
simplex $\sigma=\\{v^{i}\\}_{i\in I}$ of $K$, with
$I\subset\\{1,2,\dots,\beta\\}$, and a point
$x\in|\sigma|\subset\mathbb{R}^{n}$. Let $b(x)$ be the barycentric coordinates
of $x$ with respect to $\sigma$. Then, the barycentric coordinates of $x$ with
respect to the simplicial complex $K$ are given by the vector

$\xi(x)=(\xi_{1}(x),\xi_{2}(x),\dots,\xi_{\beta}(x))$

where $\xi_{i}=0$ if $i\not\in I$ and $(\xi_{i}(x))_{i\in I}=b(x)$.

Figure 1: The barycentric subdivision of the simplicial complex described in
Example 2.1.

###### Example 2.1

Let $K$ be a simplicial complex with a set of vertices
$\big{\\{}v^{0}=(0,0),v^{1}=(2,0),v^{2}=(0,2)\big{\\}}$. Then,
$K=\big{\\{}\\{v^{0}\\},\\{v^{1}\\},\\{v^{2}\\},\\{v^{0},v^{1}\\},\\{v^{1},v^{2}\\},$
$\\{v^{0},v^{1},v^{2}\\}\big{\\}}$ and $\operatorname{Sd}K$ is composed of the
following simplices (see Figure 1):

  * •

Vertices: $\\{v^{0},v^{1},v^{2},v^{01},v^{02},v^{12},v^{012}\\}$;

  * •

$1$-simplices:
$\big{\\{}\\{v^{0},v^{01}\\},\\{v^{0},v^{02}\\},\\{v^{0},v^{012}\\},\\\ \hskip
61.17325pt\\{v^{1},v^{01}\\},\\{v^{1},v^{12}\\},\\{v^{1},v^{012}\\},\\\ \hskip
61.17325pt\\{v^{2},v^{02}\\},\\{v^{2},v^{12}\\},\\{v^{2},v^{012}\\},\\\ \hskip
61.17325pt\\{v^{01},v^{012}\\},\\{v^{02},v^{012}\\},\\{v^{12},v^{012}\\}\big{\\}}$;

  * •

$2$-simplices:
$\big{\\{}\\{v^{0},v^{01},v^{012}\\},\\{v^{0},v^{02},v^{012}\\},\\\ \hskip
61.17325pt\\{v^{1},v^{01},v^{012}\\},\\{v^{1},v^{12},v^{012}\\},\\\ \hskip
61.17325pt\\{v^{2},v^{02},v^{012}\\},\\{v^{2},v^{12},v^{012}\\}\big{\\}}$.

Let us now consider a point $x=(\frac{1}{2},\frac{1}{2})$. The barycentric
coordinates of $x$ with respect to $K$ are
$\xi_{K}(x)=(\frac{1}{2},\frac{1}{4},\frac{1}{4})$. The barycentric
coordinates of $x$ with respect to the simplex $\\{v^{0},v^{01},v^{012}\\}$
are $(\frac{1}{4},0,\frac{3}{4})$. Hence,
$\xi_{\operatorname{Sd}K}(x)=(\frac{1}{4},0,0,0,0,0,\frac{3}{4})$.

Given two simplicial complexes $K$ and $L$ with vertex set $V$ and $W$,
respectively. A vertex map $\varphi^{(0)}:V\rightarrow W$ is a correspondence
between vertices such that for each simplex $\sigma=\\{v^{i}\\}_{i\in I}$ of
$K$, the set obtained from
$\big{\\{}\varphi^{\scriptscriptstyle(0)}(v^{i})\big{\\}}_{i\in I}$ after
removing duplicated vertices is a simplex of $L$. A vertex map induces a
continuous function called a simplicial map defined for all $x\in|K|$ as
follows:

$\mbox{$\varphi(x)=\sum_{i\in I}b_{i}(x)\varphi^{(0)}(v^{i})$ such that
$x\in|\sigma|$ and $\sigma\in K$}.$

Figure 2: All possible dichotomies of a dataset of size $3$ in
$\mathbb{R}^{2}$ shattered by a line.

###  2.2 Neural networks for classification

A dataset  $\mathcal{D}=(X,\lambda)$ of a classification problem is a pair
composed of a set of points $X\subset\mathbb{R}^{d}$ of size $N$ and a
function $\lambda$ from $X$ to $\Lambda=\\{0,1,\dots,k\\}$ where $k+1$ is the
number of classes. We denote by $y_{i}$ the image by $\lambda$ of $x_{i}\in X$
for all $i\in\\{1,2,\dots,N\\}$. The set of points $\\{x\ |x\in
X,\lambda(x)=c\\}$ with $c\in\Lambda$ is called the class $c$. The
classification task consists of finding a function
$\mathcal{N}:\mathbb{R}^{d}\rightarrow\Lambda$ that approximates $\lambda$.

In this paper, the function $\mathcal{N}$ is a neural network defined between
$\mathbb{R}^{d}$ and $\mathbb{R}^{k}$ and composed of $m+1$ functions,

$\mathcal{N}=f_{m+1}\circ f_{m}\circ\cdots f_{1}\,,$

where the integer $m>0$ is the number of hidden layers and, for
$i\in\\{1,2,\dots,m+1\\}$, the function $f_{i}:X_{i-1}\rightarrow X_{i}$ is
defined as $f_{i}(y)=\phi_{i}(W^{(i)};y_{i};b_{i})$ where
$X_{0}=\mathbb{R}^{d}$, $X_{m+1}=\mathbb{R}^{k}$,
$X_{i}\subseteq\mathbb{R}^{d_{i}}$, with $d_{0}=d$ and $d_{m+1}=k$, $W^{(i)}$
is a matrix $d_{i-1}\times d_{i}$, $b_{i}$ is a $d_{i}$-dimensional vector,
and $\phi_{i}:\mathbb{R}^{d_{i-1}}\to\mathbb{R}^{d_{i}}$ is a continuous
function. The matrices $\\{W^{(i)}\\}$ are usually updated using an
optimization algorithm such as stochastic gradient descent and are called the
weight matrices of $\mathcal{N}$.

From learning theory, we will use the concept of the Vapnik-Chervonenkis
dimension (VC-dimension) to quantify the capacity of a neural network (see
[8]). Given a point cloud $X$, a dichotomy of $X$ is one of the possible ways
to label $X$ in the binary classification context. Then we say that a model
shatters $X$ if, for any possible dichotomy of $X$, we can find the
architecture of the model that correctly classifies it. The VC-dimension of a
model is the maximum size of a dataset that can be found to be shattered by
it. A good introduction to the VC-dimension of neural networks is [9] where it
is shown that the VC-dimension of perceptrons with $n+1$ entries (including
the bias term) and hyperplanes in $\mathbb{R}^{n}$ is $n+1$. For example,
consider the two-dimensional case illustrated in Figure 2. All possible
dichotomies of a dataset of size $3$ in $\mathbb{R}^{2}$ can be shattered by a
hyperplane. However, no data set of size $4$ in $\mathbb{R}^{2}$ can be
shattered by a hyperplane. Therefore, the VC dimension of a hyperplane in
$\mathbb{R}^{2}$ is $3$.

##  3 Simplicial-map layer

In this section, we introduce the main novelty of this paper to add
interpretability to deep neural networks: the simplicial-map layer, also
called the SIMAP layer,

Unlike the previous definitions of simplicial-map neural networks [10], the
SIMAP layers are not based on Delaunay triangulation, but on barycentric
subdivisions, so we do not need to establish the support set. Furthermore, we
prove that we do not need to explicitly calculate the successive barycentric
subdivisions, since we can deduce the barycentric coordinates of a point with
respect to a subdivision from its barycentric coordinates with respect to the
simplicial complex that existed before the subdivision. Next, let us see step
by step how to compute the SIMAP layers for a dataset $D=(X,\lambda)$ of a
classification problem.

###  3.1 Convex polytope and barycentric coordinates

First, we consider that $X\subset\mathbb{R}^{n}$ lies inside the unit
hypercube $\mathcal{H}$ whose corners are the $2^{n}$ points with coordinates
$0$ or $1$. If that is not the case, we can always do an affine transformation
to achieve it.

Then, we compute an $n$-simplex $\sigma$ whose geometric realization (convex
polytope) contains $\mathcal{H}$ and, therefore, the point cloud $X$. The
following lemma provides the coordinates of the vertices of $\sigma$.

###### Lemma 3.1

Let $\mathcal{H}$ be the $n$-dimensional hypercube in $\mathbb{R}^{n}$ whose
corners are the $2^{n}$ points with coordinates $0$ or $1$. Then, the vertices
of an $n$-simplex $\sigma$ satisfying $\mathcal{H}\subset|\sigma|$ are the
rows of the following $(n+1)\times n$ matrix

$S=\begin{pmatrix}0&0&\cdots&0&0\\\ n&0&\cdots&0&0\\\ 0&n&\cdots&0&0\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&\cdots&n&0\\\ 0&0&\cdots&0&n\\\
\end{pmatrix}.$

Proof.  The corner $(0,0,\dots,0)$ of $\mathcal{H}$ is a vertex of $\sigma$
(first row of $S$). The corner $c=(c_{1},c_{2},\dots,c_{n})$ of $\mathcal{H}$
with $c_{i}=1$ for $i\in I\subset\\{1,2,\dots,n\\}$ lies in the face of
$\sigma=\\{v^{i}\\}_{i\in I\cup\\{0\\}}$ whose vertices $v^{i}$ are the
$(i+1)$-th rows of $S$ for $i\in I\cup\\{0\\}$. Specifically,
$c=\frac{1}{n}\sum_{i\in I\cup\\{0\\}}v^{i}$. The corner $(1,1,\dots,1)$ of
$\mathcal{H}$ is the barycenter of the face of $\sigma$ whose vertices are the
$(i+1)$-th rows of $S$ for $i\in\\{1,2,\dots,n\\}$. Since $|\sigma|$ is the
convex hull of the vertices of $\sigma$ and the corners of $\mathcal{H}$ are
in $|\sigma|$ then $\mathcal{H}\subset|\sigma|$. $\Box$

Now, let $\sigma$ be the simplex defined in Lemma 3.1, let $x\in|\sigma|$ and
let $b(x)$ be the barycentric coordinates of $x$ with respect to $\sigma$. Let
$T=({\bf 1}\ |\ S)$ be the matrix consisting of $S$ with an additional column
whose entries are all $1$. Taking into account that, by definition of
barycentric coordinates, $b(x)\ T=(1\ |\ x)$, the barycentric coordinates of
$x$ with respect to $\sigma$ can be easily computed, as the following result
shows.

###### Lemma 3.2

Let $\sigma$ be the $n$-simplex defined in Lemma 3.1 and let $x=(x_{1},\dots$,
$x_{n})\in|\sigma|$. Then, the barycentric coordinates
$b(x)=(b_{0}(x),\dots,b_{n}(x))$ of $x$ with respect to $\sigma$ can be
obtained by multiplying the vector $(1,x_{1},\dots.x_{n})$ by the matrix $M$,
i.e.,

$(b_{0}(x),\dots,b_{n}(x))=(1,x_{1},\dots.x_{n})\,M$

where $M$ is the $(n+1)\times(n+1)$ matrix

$M=\begin{pmatrix}1&0&0&\cdots&0&0\\\ -\frac{1}{n}&\frac{1}{n}&0&\cdots&0&0\\\
-\frac{1}{n}&0&\frac{1}{n}&\cdots&0&0\\\ \vdots&\vdots&&\ddots&&\vdots\\\
-\frac{1}{n}&0&0&\cdots&\frac{1}{n}&0\\\
-\frac{1}{n}&0&0&\cdots&0&\frac{1}{n}\\\ \end{pmatrix}$

Proof.  Using that

$(b_{0}(x),\dots,b_{n}(x))\,T=(1,x_{1},\dots,x_{n}),$

and since $T\cdot M$ is the identity matrix then:

$\begin{split}(b_{0}(x),\dots,b_{n}(x))&=(b_{0}(x),\dots,b_{n}(x))\,T\cdot
M\\\ &=(1,x_{1},\dots,x_{n})\,M.\end{split}$

$\Box$

###  3.2 Barycentric coordinates after a barycentric subdivision

Once we know how to compute the barycentric coordinates of a point
$x\in|\sigma|$, the next step is to find how to calculate the barycentric
coordinates of $x$ with respect to a subdivision of $\sigma$.

Let us consider an $n$-simplex $\sigma=\langle v^{0},\dots,v^{n}\rangle$ in
$\mathbb{R}^{n}$. Then, each ordering $(i_{0},\dots,i_{n})$ of the indices
$(0,\dots,n)$ represents one of the maximal simplices obtained by the
barycentric subdivision of $\sigma$. This way, the set of vertices of the
$n$-simplex $\sigma^{i_{0}\cdots i_{n}}$, represented by the ordering
$(i_{0},\dots,i_{n})$, are

$\\{v^{i_{0}},v^{i_{0}i_{1}},\dots,v^{i_{0}\dots i_{n}}\\}$

where $v^{i_{0}\dots i_{k}}=\operatorname{bar}\langle
v^{i_{0}},\dots,v^{i_{k}}\rangle$ for all $k\in\\{0,\dots,n\\}$, and
$\\{v^{i_{0}},\dots,v^{i_{k}}\\}$ is a face of $\sigma$.

###### Example 3.3

Following Example 2.1, the simplex $\\{v^{1},v^{12},v^{012}\\}$ will be
identified with the ordering $(1,2,0)$ and will be denoted by $\sigma^{120}$.

###### Lemma 3.4

Let $\sigma^{i_{0}\cdots i_{n}}$ be a maximal simplex of the barycentric
subdivision of $\sigma$ represented by an ordering $(i_{0},\dots,i_{n})$. Let
$x\in\mathbb{R}^{n}$ such that its barycentric coordinates with respect to
$\sigma$ are $(b_{0}(x),\dots,b_{n}(x))$. Then, if $x\in|\sigma^{i_{0}\dots
i_{n}}|$, the barycentric coordinates $(b^{1}_{0}(x),\dots,b^{1}_{n}(x))$ of
$x$ with respect to $\sigma^{i_{0}\dots i_{n}}$ satisfies that:

$(b^{1}_{0}(x),\dots,b^{1}_{n}(x))=(b_{i_{0}}(x),\dots,b_{i_{n}}(x))\,P$

where $P$ is the $(n+1)\times(n+1)$ matrix

$P=\begin{pmatrix}1&0&0&\cdots&0&0&0\\\ -1&2&0&\cdots&0&0&0\\\
0&-2&3&\cdots&0&0&0\\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&\cdots&n-1&0&0\\\ 0&0&0&\cdots&-(n-1)&n&0\\\
0&0&0&\cdots&0&-n&n+1\end{pmatrix}.$

Proof.  By definition, $x=\sum_{k=0}^{n}b^{1}_{k}(x)v^{i_{0}\dots i_{k}}$.
Since

$v^{i_{0}\dots i_{k}}=\operatorname{bar}\langle
v^{i_{0}},\dots,v^{i_{k}}\rangle=\frac{1}{k+1}\sum_{j=0}^{k}v^{i_{j}}$

then

$\begin{split}x&\mbox{$=b^{1}_{0}(x)\,v^{i_{0}}+b^{1}_{1}(x)\,\frac{1}{2}(v^{i_{0}}+v^{i_{1}})+\dots$}\\\
&\mbox{$\;\;\;+b^{1}_{n}(x)\,\frac{1}{n+1}(v^{i_{0}}+\dots+v^{i_{n}})$}\\\
&\mbox{$=\big{(}b^{1}_{0}(x)+\frac{1}{2}b^{1}_{1}(x)+\cdots+\frac{1}{n+1}b^{1}_{n}(x)\big{)}v^{i_{0}}+\cdots$}\\\
&\mbox{$\;\;\;+\frac{1}{n+1}b^{1}_{n}(x)v^{i_{n}}$.}\end{split}$

As a matrix equation,

$(b^{1}_{0}(x),\dots,b^{1}_{n}(x))\,Q=(b_{i_{0}}(x),\dots,b_{i_{n}}(x)),$

where

$Q=\begin{pmatrix}1&0&0&\cdots&0&0&0\\\
\frac{1}{2}&\frac{1}{2}&0&\cdots&0&0&0\\\
\frac{1}{3}&\frac{1}{3}&\frac{1}{3}&\cdots&0&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
\frac{1}{n}&\frac{1}{n}&\frac{1}{n}&\cdots&\frac{1}{n}&\frac{1}{n}&0\\\
\frac{1}{n+1}&\frac{1}{n+1}&\frac{1}{n+1}&\cdots&\frac{1}{n+1}&\frac{1}{n+1}&\frac{1}{n+1}\end{pmatrix}.$

Therefore,

$(b^{1}_{0}(x),\dots,b^{1}_{n}(x))\,Q\cdot
P=(b_{i_{0}}(x),\dots,b_{i_{n}}(x))\,P$

concluding the proof since $Q\cdot P$ is the identity matrix. $\Box$

We end this subsection with the following result that is useful to easily
detect the simplex $\sigma^{i_{0}\cdots i_{n}}\in\operatorname{Sd}\sigma$
where the given point $x$ is located.

###### Lemma 3.5

Let $\sigma$ be the simplex defined in Lemma 3.1 and let $x$ be a point in
$|\sigma|$. The point $x$ lies in $|\sigma^{i_{0}\cdots i_{n}}|$ if and only
if the ordering $(i_{0},\dots,i_{n})$ makes that the barycentric coordinates
$(b_{i_{0}}(x),\dots,b_{i_{n}}(x))$ are arranged in a non-increasing manner,
that is, $b_{i_{0}}(x)\geq\dots\geq b_{i_{n}}(x)$.

Proof.  Recall that the point $x$ lies in $|\sigma^{i_{0}\cdots i_{n}}|$ if
its barycentric coordinates $(b^{1}_{0}(x),\dots,b^{1}_{n}(x))$ are
nonnegative. Applying Lemma 3.4, we have $b^{1}_{n}(x)=(n+1)b_{n}(x)$, which
is nonnegative since $x\in|\sigma|$. Furthermore,
$b^{1}_{j}(x)=(j+1)(b_{j}(x)-b_{j+1}(x))$ for any $j\in\\{0,1,\dots,n-1\\}$,
which is nonnegative if and only if $b_{j}(x)\geq b_{j+1}(x)$, concluding the
proof. $\Box$

###  3.3 Training

This is the core of the SIMAP layer. In the first step, without considering
any subdivision, it is simply a perceptron with $n+1$ inputs and $k$ outputs
if the given dataset is a point cloud in $\mathbb{R}^{n}$ labelled with $k$
classes. After successive subdivisions, this part of the SIMAP layer can be
seen as the union of $m$ perceptrons with $n+1$ inputs and $k$ outputs since
the idea of this part of the SIMAP layer is that a point is classified
according to the maximal simplex of the subdivision on which it lies.

Let $\\{v^{0},\dots,v^{n}\\}$ be the vertices of the $n$-simplex $\sigma$
defined in Lemma 3.1. Assume that our input dataset $D=(X,\lambda)$ is a
finite set of points in $|\sigma|$, together with a set of $k+1$ labels
$\Lambda=\\{0,\dots,k\\}$ such that each $x\in X$ is labelled with a label
$\lambda(x)$ taken from $\Lambda$. A one-hot encoding representation is:

$L=\big{\\{}(0,\stackrel{{\scriptstyle
j}}{{\dots}},0,1,0,\stackrel{{\scriptstyle
k-j}}{{\dots}},0):\;j\in\\{0,\dots,k\\}\big{\\}}\,,$

where the one-hot vector $(0,\stackrel{{\scriptstyle
j}}{{\dots}},0,1,0,\stackrel{{\scriptstyle k-j}}{{\dots}},0)$ encodes the
$j$-th label of $\Lambda$ for $j\in\\{0,\dots,k\\}$ and
$\ell(x)=(\ell_{0}(x),\dots,\ell_{k}(x))$ encodes the label $\lambda(x)$ for
any $x\in X$.

Observe that if we define a vertex map $\varphi^{\scriptscriptstyle(0)}:V\to
L$ mapping each vertex $v\in V$ to a point in $\mathbb{R}^{k+1}$. Then, for
any $x\in|\sigma|$ with barycentric coordinates
$b(x)=(b_{0}(x),\dots,b_{n}(x))$, we have that
$\varphi(x)=\sum_{i=1}^{n}b_{i}(x)\varphi^{(0)}(v^{i})$ is a probability
distribution satisfying that the $i$-th coordinate of $\varphi(x)$ indicates
the probability that $x$ belongs to the class $i\in\Lambda$.

In general, the values of $\varphi^{\scriptscriptstyle(0)}(v^{i})$ are unknown
and the procedure explained can not be applied. Therefore, the idea for the
next step of our algorithm is to define a multiclass perceptron denoted by
$\mathcal{N}$ with an input layer with $n+1$ neurons that computes the
probability that $x\in|\sigma|$ is labelled with a label of $\Lambda$ using
the formula:

$\mathcal{N}(x)=(\mathcal{N}_{0}(x),\dots,\mathcal{N}_{k}(x))=\operatorname{softmax}\big{(}b(x)\cdot\Omega\big{)}$

where $\Omega\in\mathcal{M}_{(n+1)\times(k+1)}$ is a weight matrix and
$b(x)\in\mathbb{R}^{n+1}$ are the barycentric coordinates of $x$ with respect
to $\sigma$. The training procedure has the aim of learning the values of
$\varphi^{\scriptscriptstyle(0)}(v^{i})$ that minimizes the error:

$\mbox{$\mathcal{L}(x)=-\sum_{h=1}^{k}\ell_{h}(x)\log(\mathcal{N}_{h}(x))$}\,,$

for $x\in D$. As proven in [10],

$\frac{\partial\mathcal{L}(x)}{\partial
p_{j}^{t}}=(\mathcal{N}_{j}(x)-\ell_{j}(x))b_{t}(x).$

Then, using, for example, gradient descent, we have to update the entries
$p_{j}^{t}$ of $\Omega$, as follows:

$p_{j}^{t}:=p_{j}^{t}-\eta(\mathcal{N}_{j}(x)-\ell_{j}(x))b_{t}(x)$.

Observe that the $i$-th row of $\Omega$ can be thought of as the value of
$\varphi^{(0)}(v^{i})$ for $i\in\\{0,\dots,n\\}$.

Since the capacity of the perceptron $\mathcal{N}$ is limited, the idea is to
train another perceptron $\mathcal{N}^{1}$ from it, emulating a barycentric
subdivision. This new perceptron $\mathcal{N}^{1}$ is initialized as

$\mathcal{N}^{1}(x)=\operatorname{softmax}\big{(}\xi_{\operatorname{Sd}\sigma}(x)\cdot\Omega^{1}\big{)}$

where $\Omega^{1}$ is a $(2^{n+1}-1)\times(k+1)$ matrix whose rows are in one-
to-one correspondence with the value of $\varphi$ applied to the vertices of
$\operatorname{Sd}\sigma$ and $\xi_{\operatorname{Sd}\sigma}(x)$ is computed
as explained in subsection 3.2. Specifically, by construction, we have

$\xi_{\operatorname{Sd}\sigma}(x)\cdot\Omega^{1}=\xi_{\operatorname{Sd}\sigma}(x)\cdot
Q\cdot\Omega=b(x)\cdot\Omega\,.$

So, initially, $\mathcal{N}^{1}(x)=\mathcal{N}(x)$ for any $x\in|\sigma|$.
Then, the weights of $\Omega^{1}$ are updated using gradient descent as
explained above. The process can be iterated until a given error is reached.
We will denote by $VC(\mathcal{N})$ to the VC dimension of the neural network
$\mathcal{N}$. We have the following result.

###### Theorem 3.6

Let $\\{v^{0},\dots,v^{n}\\}\subset\mathbb{R}^{n}$ be the vertices of the
$n$-simplex $\sigma$ defined in Lemma 3.1. Let $\mathcal{N}$ be the neural
network defined from $\sigma$ and $\mathcal{N}^{k}$ the neural network defined
from $\operatorname{Sd}^{k}\sigma$ with $k\geq 0$ as explained above. Then,

$VC(\mathcal{N}^{k})=((n+1)!)^{k}\cdot(n+1)\,.$

Proof.  For $k=0$, the neural network $\mathcal{N}$ defined from $\sigma$ is a
perceptron with $n+1$ entries (the $n+1$ barycentric coordinates with respect
to $\sigma$) and no bias whose VC-dimension is $n+1$. For $k>0$ and
$x\in|\operatorname{Sd}^{k}\sigma|$, there exists a maximal simplex $\mu$ in
$\operatorname{Sd}^{k}\sigma$ such that $x\in|\mu|$. Then,
$\xi_{\operatorname{Sd}^{k}\sigma}(x)$ has at most $n+1$ non-null coordinates.
Since $\operatorname{Sd}^{k}\sigma$ has $((n+1)!)^{k}$ maximal simplices, the
neural network $\mathcal{N}^{k}$ defined from $\operatorname{Sd}^{k}\sigma$
acts as $((n+1)!)^{k}$ independent perceptions. Then, the VC-dimension of
$\mathcal{N}^{k}$ is $((n+1)!)^{k}\cdot(n+1)$. $\Box$

Let us remark that, by definition, there exists no simplicial map that maps a
simplex to another simplex of higher dimension.

###### Lemma 3.7

Let $D=(X,\lambda)$ be a dataset with $X\subset\mathbb{R}^{n}$. Let $\sigma$
be the $n$-simplex defined in Lemma 3.1. Let $\mathcal{N}^{0}$ be the neural
network defined from $\sigma$ and $\mathcal{N}^{k}$ the neural network defined
from $\operatorname{Sd}^{k}\sigma$ with $k>0$ as explained above. Fixing
$k\geq 0$, let $\mu\in\operatorname{Sd}^{k}\sigma$ such that
$Y=\\{x^{0},\dots,x^{n+1}\\}\subseteq X$ is a set of points in $|\mu|$
satisfying that no two points of $Y$ have the same label. Then,
$\mathcal{N}^{k}$ cannot correctly classify $X$.

Proof.  The dimension of $\mu$ is $n+1$. Then, for all $x\in Y$,
$\xi_{\operatorname{Sd}^{k}\sigma}(x)$ is possibly not null for fixed $n+1$
coordinates. Therefore, the output of
$\mathcal{N}^{k}(x)=\operatorname{softmax}(\xi_{\operatorname{Sd}^{k}\sigma}(x)\cdot\Omega^{k})$
is again possibly not null for fixed $n+1$ coordinates, so it can cannot
correctly classify $Y$ (nor $X$) since the number of classes for $Y$ (and
then, $X$) is at least $n+2$. $\Box$

Summing up, given a dataset $D=(X,\lambda)$ with $X\subseteq\mathbb{R}^{n}$
and $k+1$ classes, the SIMAP layer $\mathcal{N}^{k}$ is a neural network with
no hidden layer and acting as $\big{(}(n+1)!\big{)}^{k}$ independent
perceptrons with at most $n+1$ activated neurons in the input layer and $k+1$
neurons in the output layer. A pseudocode for SIMAP layers can be found in
Algorithm 1.

Figure 3: Dataset of Example 3.8. On the left, the input data for the binary
classification is shown. On the right, its translation into the simplex
$\sigma$ of Lemma 3.1. Algorithm 1 Pseudocode to compute SIMAP layers

Input: A dataset $D=(X,\lambda)$ with a set of labels $\Lambda$.

An integer $\ell$ (number of subdivisions).

Output: A SIMAP layer that generalizes $D$.

if $X$ does not lie in the hypercube $\mathcal{H}$ then

Transform the points of $X$ so that they lie in $\mathcal{H}$ (see Subsection
3.1)

end if

$k=0$

Compute the barycentric coordinates $b(x)$ of the all points $x\in X$ (see
Subsection 3.1).

Train the perceptron
$\mathcal{N}(\;)=\operatorname{softmax}(b(\;)\cdot\Omega)$ (as explained in
Subsection 3.3).

while $k\leq\ell$ do

Compute the barycentric coordinates $\xi(x)$ with respect to
$\operatorname{Sd}^{k}\sigma$ of the all points $x\in X$ (as explained in
Subsection 3.2).

Train the perceptron
$\mathcal{N}^{k}(\;)=\operatorname{softmax}(\xi(\;)\cdot\Omega)$ (as explained
in Subsection 3.3).

end while

Below, we provide a step-by-step example following Algorithm 1 showing how to
compute a SIMAP layer for a given dataset $D$.

Figure 4: Training curves of Example 3.8. On the top: the loss function is
shown. On the bottom: the accuracy values for the different epochs are shown.

Figure 5: From top to bottom, the decision boundaries for Example 3.8 for the
different models with respect to $\sigma$, $\operatorname{Sd}\sigma$, and
$\operatorname{Sd}^{2}\sigma$, respectively.

###### Example 3.8

Let us consider a dataset $D=(X,\lambda)$, with $X\subset\mathbb{R}^{2}$,
inspired by the classic XOR-dataset (see Figure 3). Specifically,
$X=\big{\\{}v^{0}=(1,1),v^{1}=\left(\frac{3}{2},\frac{3}{2}\right),v^{2}=\left(\frac{5}{2},\frac{5}{2}\right),v^{3}=(3,3),v^{4}=(1,3),v^{5}=\left(\frac{3}{2},\frac{5}{2}\right),v^{6}=\left(\frac{5}{2},\frac{3}{2}\right),v^{7}=(3,1)\big{\\}}$.
The point cloud $X$ is labelled by $\lambda(v^{i})=0$ if $i<4$ and
$\lambda(v^{i})=1$ in any other case. The steps to calculate a SIMAP layer
that classifies $X$ are the following:

  1. 1.

Dataset preparation: the points of $X$ are centered inside the simplex
$\sigma$ of Lemma 3.1 by subtracting the mean and then rescaling the
coordinates of the points between $0$ and $1$ (between $0$ and half of the
dimension for the general case). Hence, we obtain the following matrix whose
rows are the coordinates of the points of the transformed point cloud $X$:

$V=\begin{pmatrix}0&0\\\ 0.25&0.25\\\ 0.75&0.75\\\ 1&1\\\ 0&1\\\ 0.25&0.75\\\
0.75&0.25\\\ 1&0\end{pmatrix}.$

  2. 2.

Barycentric coordinates computation: we compute the barycentric coordinates of
each point of $X$ with respect to the simplex $\sigma$. To do so, we have to
multiply the matrix $V$ by the matrix $M$. obtained using the formula from
Lemma 3.2:

$\begin{pmatrix}1&0&0\\\ 1&0.25&0.25\\\ 1&0.75&0.75\\\ 1&1&1\\\ 1&0&1\\\
1&0.25&0.75\\\ 1&0.75&0.25\\\ 1&1&0\end{pmatrix}\cdot\begin{pmatrix}1&0&0\\\
-\frac{1}{2}&\frac{1}{2}&0\\\ -\frac{1}{2}&0&\frac{1}{2}\end{pmatrix}.$

  3. 3.

SIMAP layer training: A first SIMAP layer can be trained using as input the
matrix obtained in the previous step. After training, we obtained the
following weight matrix:

$W_{0}=\begin{pmatrix}0.46&0.48\\\ -0.15&-0.34\\\ 0.79&0.93\end{pmatrix}$

  4. 4.

Transfer learning: to increase the capacity of the SIMAP layer, we apply
barycentric subdivisions and inherit the previous weight matrix. Let us denote
by $B_{0}$ the $25\times 7$ matrix whose entries are the barycentric
coordinates of the vertices of $\operatorname{Sd}\sigma$ with respect to
$\sigma$ in the correct row ordering, i.e.
$B_{0}=\xi_{\operatorname{Sd}\sigma}(x)\cdot\Omega^{1}$. Then, the init weight
matrix $W_{1}$ of the new SIMAP layer is the $25\times 2$ matrix:

$W_{1}=B_{0}\cdot W_{0}\,.$

The steps 3 and 4 can be iterated as many times as needed. In our case, we
applied the methodology for two barycentric subdivisions and obtained the
curves for the training accuracy and loss function shown in Figure 4.

In figure 5 we can see the decision boundaries for the different models with
respect to $\sigma$, $\operatorname{Sd}\sigma$, and
$\operatorname{Sd}^{2}\sigma$. As expected, the complexity of the decision
boundary increases when the barycentric subdivision process is iterated.

##  4 Experiments

This section presents various experiments. In the first set of experiments,
synthetic datasets were used. In the second set of experiments, the proposed
methodology is applied as a final layer of a convolutional neural network to
classify the MNIST dataset. In all experiments, we used the Adam [11] training
algorithm.

###  4.1 Experiments with synthetic datasets

In this set of experiments, we used synthetic datasets composed of 500 points
for binary classification from [12] scikit-learn implementation for different
numbers of features (2,3,4 and 5). The datasets were split into training and
test sets with a proportion of $80\%$ and $20\%$, respectively. The steps
followed are the same as in Example 3.8. The SIMAP layers were trained for
$1000$ epochs and barycentric subdivisions were applied two times. The results
are shown in Table 1 where the values are the mean of $5$ repetitions. As we
can see from the results, when applying barycentric subdivisions, the SIMAP
layer quickly reaches overfitting, being, in general, enough to apply just one
barycentric subdivision for this dataset.

Table 1: Loss and accuracy values on training and test set for the experiment in Section 4.1.  n | no. subdivisions | Loss | Accuracy
---|---|---|---
Train | Test | Train | Test
2 | 0 | 0.32 | 0.33 | 0.9 | 0.87
1 | 0.26 | 0.27 | 0.91 | 0.91
2 | 0.22 | 0.26 | 0.91 | 0.89
3 | 0 | 0.44 | 0.4 | 0.82 | 0.84
1 | 0.39 | 0.37 | 0.83 | 0.81
2 | 0.26 | 0.48 | 0.87 | 0.83
4 | 0 | 0.38 | 0.36 | 0.87 | 0.87
1 | 0.29 | 0.29 | 0.88 | 0.85
2 | 0.14 | 1.04 | 0.95 | 0.83
5 | 0 | 0.34 | 0.31 | 0.87 | 0.87
1 | 0.24 | 0.29 | 0.92 | 0.87
2 | 0.03 | 0.99 | 0.92 | 0.76

###  4.2 Convolutional network in conjunction with a SIMAP layer

In this experiment, we show the case where a convolutional layer is combined
with a SIMAP layer. We used the MNIST dataset.

  1. 1.

Dataset: The MNIST dataset is composed of $60000$ training grayscale images
and $10000$ test images of size $28\times 28$.

  2. 2.

Convolutional network training: A convolutional neural network is trained for
$10$ epochs using the training set. The architecture used is the following:

Conv2D $\displaystyle(26,26,28)$ MaxPooling2D $\displaystyle(13,13,28)$ Conv2D
$\displaystyle(1,11,64)$ MaxPooling2D $\displaystyle(5,5,64)$ Conv2D
$\displaystyle(1,1,4)$ Flatten $\displaystyle(4)$ Dense $\displaystyle(64)$
Dense $\displaystyle(10)$

This architecture reached an accuracy of $0.98$ and a loss value of $0.057$ on
the test set.

  3. 3.

Dataset preparation: The output of the flatten layer is used as a four-
dimensional point cloud. Next, it is scaled and translated within the simplex
$\sigma$ defined in Lemma 3.1. We denote by $X$ the transformed point cloud.

  4. 4.

Computation of the barycentric coordinates: The barycentric coordinates of the
points of $X$ with respect to the simplex $\sigma$, as well as their
barycentric coordinates with respect to the first and second barycentric
subdivisions ($\operatorname{Sd}\sigma$, and $\operatorname{Sd}^{2}\sigma$),
are calculated following the methodology presented in Lemma 3.2.

  5. 5.

SIMAP layer training: Sequentially, three SIMAP layers were trained, with
weights transferred from  a SIMAP layer to the next as explained in Section
3.3. The accuracy and loss values reached in the test set are shown in Table
2.

Table 2: Loss and accuracy values of the SIMAP layer with the different barycentric subdivision iterations and the Convolutional Neural Network without the SIMAP layer. | Number of subdivisions | Loss | Accuracy
---|---|---|---
| 0 | 0.39 | 0.97
| 1 | 0.06 | 0.98
SIMAP layer | 2 | 0.05 | 0.98
CNN |  | 0.06 | 0.98

Based on the results obtained, we can conclude that achieving good performance
does not require many barycentric subdivisions. Furthermore, the performance
of convolutional neural networks is not adversely affected by the application
of the SIMAP layer.

##  5 Conclusions and future work

Blind optimization methods have shown to be an efficient method for finding an
accurate set of weights to solve real-world problems with neural networks.
This ability to solve problems anyhow has been the main target during the
first years of the development of Deep Learning. Nevertheless, the wide use of
Artificial Intelligence (AI) methods in many domestic areas is turning the
perception of the social use of AI. The development of AI methods that solve
problems anyhow is no longer socially accepted. In this way, in the next
years, the interpretability of the model, i.e., the ability to interpret the
decision-making protocols in AI in a human-readable way, will be an
inexcusable requirement in the development of AI. This paper is a contribution
to this research line, and, to the best of our knowledge, no other layer based
on simplicial map methods for improving the interpretability of neural
networks has been presented to establish a comparison.

Specifically, the main contributions of this work are the SIMAP layer
definition and an efficient matrix-based algorithm for their implementation.
The benefits of this new layer are their interpretability and their
possibility to be applied with deep learning models such as convolutional
neural networks. In the future, we want to exploit their explainability
abilities, as well as develop specific barycentric subdivisions to reduce the
number of needed simplices and, therefore, to reduce the complexity.

## Code availability

The code for the experiments is available in
https://github.com/Cimagroup/SIMAP-layer.

## Acknowledgments

The work was supported in part by the European Union HORIZON-
CL4-2021-HUMAN-01-01 under grant agreement 101070028 (REXASI-PRO) and by
TED2021-129438B-I00 / AEI/10.13039/501100011033 / Unión Europea
NextGenerationEU/PRTR.

## References

  * [1] Q. Zhang, Y. N. Wu, and S.-C. Zhu, “Interpretable convolutional neural networks,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2018. [Online]. Available: https://doi.org/10.3389/frai.2023.974295
  * [2] Z. Liu and F. Xu, “Interpretable neural networks: Principles and applications,” _Frontiers in Artificial Intelligence_ , vol. 6, p. 974295, 2023.
  * [3] E. Paluzo-Hidalgo, R. Gonzalez-Diaz, and M. A. Gutiérrez-Naranjo, “Two-hidden-layer feed-forward networks are universal approximators: A constructive approach,” _Neural Networks_ , vol. 131, pp. 29–36, 2020. [Online]. Available: https://doi.org/10.1016/j.neunet.2020.07.021
  * [4] E. Paluzo-Hidalgo, M. A. Gutiérrez-Naranjo, and R. Gonzalez-Diaz, “Explainability in simplicial map neural networks,” _arXiv_ , 2023. [Online]. Available: arXiv:2306.00010v1
  * [5] T. H. Chang, L. T. Watson, T. C. Lux, A. R. Butt, K. W. Cameron, and Y. Hong, “Algorithm 1012: Delaunaysparse: Interpolation via a sparse subset of the delaunay triangulation in medium to high dimensions,” _ACM Transactions on Mathematical Software_ , vol. 46, no. 4, pp. 1–20, 2020. [Online]. Available: https://doi.org/10.1145/3422818
  * [6] R. B. Gabrielsson, B. J. Nelson, A. Dwaraknath, and P. Skraba, “A topology layer for machine learning,” in _Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics_ , ser. Proceedings of Machine Learning Research, S. Chiappa and R. Calandra, Eds., vol. 108. PMLR, 2020, pp. 1553–1563. [Online]. Available: https://proceedings.mlr.press/v108/gabrielsson20a.html
  * [7] M. Hajij, G. Zamzmi, T. Papamarkou, N. Miolane, A. Guzmán-Sáenz, K. N. Ramamurthy, T. Birdal, T. K. Dey, S. Mukherjee, S. N. Samaga, N. Livesay, R. Walters, P. Rosen, and M. T. Schaub, “Topological deep learning: Going beyond graph data,” 2023. [Online]. Available: https://arxiv.org/abs/2206.00606
  * [8] V. N. Vapnik and A. Y. Chervonenkis, “On the uniform convergence of relative frequencies of events to their probabilities,” _Theory of Probability and its Applications_ , vol. 16, no. 2, pp. 264–280, 1971. [Online]. Available: http://link.aip.org/link/?TPR/16/264/1
  * [9] E. D. Sontag, “Vc-dimension of neural networks,” in _Neural Networks and Machine Learning_ , 2018, p. 69–95. [Online]. Available: http://www.sontaglab.org/FTPDIR/vc-expo.pdf
  * [10] E. Paluzo-Hidalgo, M. A. Gutiérrez-Naranjo, and R. Gonzalez-Diaz, “Explainability in simplicial map neural networks. under review,” 2023.
  * [11] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , Y. Bengio and Y. LeCun, Eds., 2015. [Online]. Available: http://arxiv.org/abs/1412.6980
  * [12] I. M. Guyon, “Design of experiments for the nips 2003 variable selection benchmark,” 2003. [Online]. Available: https://archive.ics.uci.edu/ml/machine-learning-databases/dorothea/Dataset.pdf