#  $1$-Bit SubTHz RIS with Planar Tightly Coupled Dipoles: Beam Shaping and
Prototypes

Xianjun Ma1, Yonggang Zhou1, Qi Luo2, Yihan Ma2, Kyriakos Stylianopoulos3, and
George C. Alexandropoulos3  1 Nanjing University of Aeronautics and
Astronautics, Nanjing, China  2 University of Hertfordshire, Hatfield, UK  3
Department of Informatics and Telecommunications, National and Kapodistrian
University of Athens, Greece  emails<EMAIL_ADDRESS>{q.luo2,
<EMAIL_ADDRESS>{kstylianop<EMAIL_ADDRESS>
###### Abstract

In this paper, a proof-of-concept study of a $1$-bit wideband reconfigurable
intelligent surface (RIS) comprising planar tightly coupled dipoles (PTCD) is
presented. The developed RIS operates at subTHz frequencies and a $3$-dB gain
bandwidth of $27.4\%$ with the center frequency at $102$ GHz is shown to be
obtainable via full-wave electromagnetic simulations. The binary phase shift
offered by each RIS unit element is enabled by changing the polarization of
the reflected wave by $180^{\circ}$. The proposed PTCD-based RIS has a planar
configuration with one dielectric layer bonded to a ground plane, and hence,
it can be fabricated by using cost-effective printed circuit board (PCB)
technology. We analytically calculate the response of the entire designed RIS
and showcase that a good agreement between that result and equivalent full-
wave simulations is obtained. To efficiently compute the $1$-bit RIS response
for different pointing directions, thus, designing a directive beam codebook,
we devise a fast approximate beamforming optimization approach, which is
compared with time-consuming full-wave simulations. Finally, to prove our
concept, we present several passive prototypes with frozen beams for the
proposed $1$-bit wideband RIS.

###### Index Terms:

Reconfigurable intelligent surface, tightly coupled dipole, THz, beam shaping,
wideband.

##  I Introduction

Reconfigurable intelligent surfaces (RISs) have been identified as one of the
key physical layer technologies for the upcoming sixth generation (6G) of
wireless communications [1, 2, 3]. Such metasurfaces are capable of reflecting
their incident waves to desired angles, enriching specular channels with
scattering, and enabling over-the-air analogue computations [4, 5], paving the
way for the envisioned smart wireless environments [6, 7]. Another highly
expected feature of 6G is the promising utilization of the THz frequency band
that encapsulates extensive unallocated bandwidth, which can be used to
support rate- and localization/sensing-demanding applications, while
exhibiting increased confidentiality and anti-interference capabilities [8].

The combination of the latter two 6G technologies is lately attracting
significant research [9, 10], mainly due to the fact that reconfigurable
metasurfaces can contribute in extending the coverage of THz links. RIS unit
cells operating at THz were recently presented in [11, 12, 13], where the
tunable characteristics of their reflective beams were achieved through a
tunable dielectric substrate equipped with voltage control. In [14], the beam
tuning characteristics were dynamically controlled based on a liquid crystal
unit structure, indicating that this design method of reconfigurable
reflection arrays is a promising one for THz RISs. A broadband tightly coupled
reflectarray structure has been also reported for THz RIS unit elements in
[15, 16]. According to the relevant literature, the tightly coupled dipole
technology has the advantage of being ultra wideband for antenna application,
which makes tightly coupled antenna elements another candidate design solution
for THz RISs. In fact, remarkable progress has been made in the design of
tightly coupled RIS unit elements for those frequencies, and $1$-bit and
$2$-bit wideband and ultra-wideband reflectarrays have been designed [17, 18,
19].

In this paper, we capitalize on our recent ultra-wideband subTHz RIS [20],
which was based on tightly coupled dipoles that are vertically placed on a
ground plane, and present a novel efficient RIS unit element design. Although
the elements in the design of [20] can be fabricated by using conventional
printed circuit board (PCB) technology, the assembly of each individual RIS
element is challenging and the mechanical robustness of the overall design is
poor. We have thus improved that design by modifying the RIS unit cell to a
planar structure. This design can be easily fabricated via multi-layer PCB
technology with significantly reduced complexity and cost for the assembly
process. Compared with the design in [20], the bandwidth of the presented
subTHZ RIS design is reduced, but the resulting planar tightly coupled dipoles
(PTCD)-based RIS still possesses a wideband response. In particular, the
simulated $3$-dB gain bandwidth is $27.4\%$ with the center frequency at $102$
GHz.

The remainder of the paper is organized as follows. Section II details the
design of the proposed RIS unit cell, including full-wave simulations for the
entire RIS. Section III presents the beamforming modeling and algorithmic
approach to calculate the response of the developed $1$-bit wideband RIS to
arbitrary directions, while Section IV concludes the paper.

##  II The Proposed PTCD-based subTHz RIS Design

In this section, we present the configuration of the planar tightly coupled
dipole, which is used as the unit element of the proposed subTHz RIS.
Simulated radiation patterns via electromagnetic (EM) software of the entire
RIS are also included and discussed.

###  II-A $1$-Bit Unit Cell Design

Fig. 1 illustrates the configuration of the proposed PTCD structure which
constitutes the unit cell of the designed subTHz RIS. The unit cell is
designed with a printed butterfly dipole element placed on a RO3003 substrate
with a thickness of $0.25$ mm. Below the RO3003, there is a RO4450F substrate
with a thickness of $0.2$ mm, which is also backed by a ground plane on the
other side. Additionally, another dielectric layer, FR4 with a thickness of
$1.5$ mm, is introduced below the metallic ground to enhance the overall
structure’s mechanical strength. As a result, the total thickness of the
reflecting unit cell is $2$ mm, including the thickness of the copper layers.
By controlling the ON/OFF states of the RF switch placed at the same layer of
the printed dipole, the reflection phase of the unit cell can be tuned with
$180^{\circ}$ phase difference. As a proof-of-concept study, at this stage, we
model the RF switch as an ideal one and do not consider the equivalent circuit
model of the switch, which is left for the extension of this work.

Figure 1: The illustration in 3D of the designed RIS unit cell, which is a
planar tightly coupled butterfly dipole printed on the top layers of the
substrate. The RIS element is linearly polarized in the $x$-axis. An RF switch
is placed in the middle of the dipole to rotate the polarization of the
impinging wave.

The current distribution of the dipole for different states of the RF switch
is depicted in Figs. 2 and 3. As shown, the surface current on the planar
dipole was rotated by $180^{o}$, implying that the polarization of the
reflected wave is changed by $180^{o}$. In this simulation, a linear polarized
normal incident wave along the $x$ direction is employed to illuminate the
dipole antenna with the same linear polarization. Figs. 4 and 5 demonstrate
the corresponding reflection phase and amplitude of the unit cell with binary
phase states, considering the normal incident angle. As shown in Fig. 4, the
RIS unit cell shows low insertion loss. It is noted that the actual insertion
loss of the RIS unit cell will be higher when considering the equivalent
circuit model of the subTHz switch. From Fig. 5, it can be seen that the unit
cell exhibits a wideband phase response with stable $180^{o}$ phase
differences from the frequency $88$ to $116$ GHz for the two different states
of the RF switch.

Figure 2: The current distribution of the designed $1$-bit RIS unit cell when
the RF switch is ON. Figure 3: The current distribution of the designed
$1$-bit RIS unit cell when the RF switch is OFF. Figure 4: The reflection
amplitude of the designed $1$-bit RIS unit cell when the RF switch is ON and
OFF as a function of the operating frequency in the $[88.5,116.7]$ GHz range.
Figure 5: The reflection phase of the designed $1$-bit RIS unit cell when the
RF switch is ON and OFF as a function of the operating frequency in the
$[88.5,116.7]$ GHz range.

###  II-B RIS Design and Simulated Response

We use the previously presented unit cell to synthesize an RIS aperture
comprising $31\times 31$ elements in CST Microwave Studio. Considering the
limitation of the testing facility, we choose to place the source of the
incident wave relatively close to the RIS, and in particular, we set
$f/D=0.5$, where $f$ is the focal distance and $D$ is the diameter of the
aperture of the reflecting surface. The simulation model for the designed RIS
is illustrated in Fig. 6.

Figure 6: The simulation model of the wideband subTHz RIS consisting of
$31\times 31$ tightly-coupled dipole elements.

In Fig. 7, the simulated radiation pattern of the designed RIS at the
$H$-plane is demonstrated when the reflection beam is configured to
$22.5^{\circ}$. It is noted that despite the phase error resulting from the
$1$-bit phase quantization, the beam pointing error is small and the sidelobe
level (SLL) is larger than $10$ dB. We have also conducted full-wave EM
simulations to extract the active radiation pattern of the unit cell and
analytically calculated the response of the entire RIS. Figs. 10 and 11 depict
the results of the pattern calculated by the active method and that of
equivalent EM simulations, when the RIS is reconfigured to $0^{\circ}$ and
$22.5^{\circ}$ beam scanning, respectively. The approximated beams will be
discussed next in Section III. As observed, the analytical results agree well
with the simulated results for the main beam, while some differences are
present in the SLLs. The disagreement is mainly caused due to the fact that,
during the analytical calculation, we ignored the change of the mutual
couplings at the edge elements and the response of the RIS unit cell with
arbitrary incident angles. According to the simulation results, it was found
that when the RIS reflects the incident wave to foresight, the $3$-dB gain
bandwidth is $27.4\%$ with the center frequency at $102$ GHz.

Figure 7: The simulated radiation of the designed wideband subTHz RIS when the
beam is configured to point to the $22.5^{\circ}$ direction.

To validate the proposed wideband subTHz RIS design with $1$-bit PTCD-based
unit cells and evaluate the fabrication accuracy, several passive prototypes
have been fabricated. These prototypes were designed to provide reflection
beams pointing at the fixed angles $0^{\circ}$, $22.5^{\circ}$, and
$45^{\circ}$. Fig. 8 shows a photo of one of the fabricated prototypes. Fig. 9
shows the comparison between the measured and simulated radiation pattern with
the beam steered to $22.5^{\circ}$ degree. As shown, good agreement at the
main lobe is obtained.

Figure 8: Photo of the fabricated $31\times 31$ subTHZ RIS with $1$-bit unit
cells. Figure 9: The measured and simulated radiation of the designed wideband
subTHz RIS when the beam is configured to point to the $22.5^{\circ}$
direction.

##  III RIS Beam Shaping Optimization

In this section, we capitalize on the RIS beamforming optimization approach
recently developed in [21] to devise reflection beams for the proposed
wideband PTCD-based $1$-bit RIS design to any desired direction. This approach
can be quite efficient for efficiently designing static as well as dynamic RIS
beam codebooks for communications [22, 23, 24], localization/sensing [10, 25],
as well integrated sensing and communications [26].

To deploy the approach in [21], we consider a system model comprising a
single-antenna transmitter (TX) placed in front of the center of the broadside
of the designed RIS as in Section II-B, and a candidate receiver (RX) antenna
at an intended position $\mathbf{p}$. Considering the pure line of sight (LoS)
reflected link from the TX to the RX, and given an arbitrary RIS configuration
vector $\boldsymbol{\omega}^{T}$ $\in\mathbb{C}^{M\times 1}$ for an RIS of $M$
unit cells, the corresponding beam pattern toward $\mathbf{p}$ has the form:

$G(\mathbf{p})=\mathbf{\boldsymbol{\omega}}^{\rm{T}}\mathbf{b}(\mathbf{p},\mathbf{p}_{\rm
TX}),$ (1)

where $\boldsymbol{\omega}$ is an $M$-element column vector including the RIS
beam pattern and $\rm{T}$ represents vector transposition, $\mathbf{p}_{\rm
TX}$ is the TX position,
$\mathbf{b}(\mathbf{p},\mathbf{p}_{\rm{TX}})=\mathbf{a}(\mathbf{p})\odot\mathbf{a}(\mathbf{p}_{\rm{TX}})$,
with $\odot$ denoting the element-wise product, and
$[\mathbf{a}(\mathbf{p})]_{m}=\exp\left(-j\frac{2\pi}{\lambda}\left(\|\mathbf{p}-\mathbf{p}_{m}\|-\|\mathbf{p}-\mathbf{p}_{\rm{RIS}}\|\right)\right)$
denotes the beam-focusing vector with $\lambda$ being the wavelength,
$\mathbf{p}_{m}$ is the position of the $m$-th ($m=1,2,\dots M$) RIS unit
cell, and $\mathbf{p}_{\rm{RIS}}$ represents the RIS center location. This
formulation holds both for near- and far-field scenarios, however, for our
specific setup, we choose the RX coordinates at a fixed distance well beyond
the Fraunhofer limit to investigate beam patterns that depend only on the
targeted RX’s angles. The positions of the TX, $\mathbf{p}_{\rm{TX}}$, and
that of the RIS, $\mathbf{p}_{\rm{RIS}}$, are considered fixed. According to
[21, Sec. 3.2], the RIS beam pattern $\boldsymbol{\omega}$ in (1) can be
optimized to realize any desired reflection beam for any given RX position,
disregarding the environmental factors of path loss, transmission power, and
RX noise. To this end, the beam pattern for a directional beam can be
represented by:

$G(\mathbf{p})=\left(\mathbf{b}^{\ast}\left(\mathbf{p},\mathbf{p}_{\rm
TX}\right)\right)^{\rm{T}}\mathbf{b}\left(\mathbf{p},\mathbf{p}_{\rm
TX}\right),$ (2)

where $\mathbf{b}^{\ast}\left(\mathbf{p},\mathbf{p}_{\rm
TX}\right)=\boldsymbol{\omega}^{T}$ is the conjugate of
$\mathbf{b}\left(\mathbf{p},\mathbf{p}_{\rm TX}\right)$.

Following this system modeling approach, the obtained phase shifts of the RIS
configuration fall under the continuous regime in $[0,2\pi]$. Then, using the
derivations presented in [21], a discretized beam pattern
$\hat{G}(\mathbf{p})$ (i.e., a discrete approximation of $G(\mathbf{p})$) can
be numerically computed, whose states are guaranteed to belong in an arbitrary
discrete set $\mathcal{K}$. For the $1$-bit subTHz RIS design presented in the
previous Section II-B, the phase shift of each RIX unit cell can take values
from the discrete set $\mathcal{K}=\\{0,\pi\\}$, while their reflection
amplitudes can be safely assumed to be unitary. The RIS beam-shaping
methodology involves iteratively taking appropriately scaled gradient steps
over a modified version of the unconstrained objective, before projecting onto
the feasible set of solutions. The complete formulation and algorithmic
approach for RIS beam shaping is given in Algorithm $3$ of [21].

Figure 10: Full-wave EM simulated, analytical, and approximated beam patterns
when the RIS is configured to reflect the incident wave to the $0^{\circ}$
direction. Figure 11: Full-wave EM simulated, analytical, and approximated
beam patterns when the RIS is configured to reflect the incident wave to the
$22.5^{\circ}$ direction.

In the previously discussed Figs. 10 and 11, patterns for the approximated
beams designed using [21]’s Algorithm $3$ are included, and compared with
equivalent analytically computed ones as well as beam patterns obtained by
full-wave EM simulations of the fabricated unit cells. As it can be seen, the
gain of each beam pattern is normalized in the dB scale, i.e., the main lobe’s
gain is $0$ dB. It can be observed that the approximated beam patterns can
achieve effective focusing at the desired pointing angles since the beam
widths at these angles almost coincide among all three cases. The SLL levels
of the approximated beams are between $6.5$ and $7$ dB, and the overall beam
profiles are flatter and with SLLs placed more closely to the main lobe, as
compared to full-wave EM simulations and theoretical calculations for the
fabricated responses. All in all, the results in these two figures indicate
that, while arbitrary discrete configurations can be produced to perform beam
steering at desired angles, more accurate approximation approaches and
algorithms are required to reduce reflections at the SLLs. In future work, we
plan to extend the approach in [21] to consider tunable tradeoffs between the
gains of the main lobe and those of the SLLs.

##  IV Conclusion

A wideband RIS design operating at the subTHz frequency band was presented in
this paper. The proposed design comprised PTCD-based unit cells, each being
integrated with an RF switch for realizing a $1$-bit phase shift. The
conducted full-wave EM simulations showcased that the designed RIS has
excellent beam scanning ability. Directional beams in different direction
angles were realized by modifying the $1$-bit phase distribution of the RIS
structure.

To deal with the increased complexity of extremely large RIS/array full-wave
EM simulations, we presented a fast approximate beam shaping approach that,
based on a projection gradient descent technique, approximates ideal
continuous reflection patterns to obtain beam patterns of arbitrary large
realistic RISs with unit cells of effectively discrete phase shifts. It was
demonstrated that the proposed beam shaping approximation technique exhibits
sufficient agreement with equivalent full-wave EM simulations in the main
lobes while requiring further refinement to match the SLL behavior. The latter
constitutes our future research direction which will pave the way for
efficient designs of static as well as dynamic beam codebooks for extremely
large RIS in subTHz and THz bands.

## Acknowledgments

This work was supported by the UK Royal Society with grant number ${\rm
IES}\\!\ \\!{\rm R}2\\!\ \\!212064$ and the Smart Networks and Services Joint
Undertaking (SNS JU) under the European Union’s Horizon Europe research and
innovation programme under Grant Agreement No 101097101, including top-up
funding by UK Research and Innovation (UKRI) under the UK government’s Horizon
Europe funding guarantee.

## References

  * [1] “Reconfigurable intelligent surfaces (RIS); Use cases, deployment scenarios and requirements,” ETSI GR RIS 001, ETSI, Apr. 2023.
  * [2] G. C. Alexandropoulos, G. Lerosey, M. Debbah, and M. Fink, “Reconfigurable intelligent surfaces and metamaterials: The potential of wave propagation control for 6G wireless communications,” _IEEE ComSoc TCCN Newslett._ , vol. 6, no. 1, pp. 25–37, Jun. 2020.
  * [3] M. Jian, G. C. Alexandropoulos, E. Basar, C. Huang, R. Liu, Y. Liu, and C. Yuen, “Reconfigurable intelligent surfaces for wireless communications: Overview of hardware designs, channel models, and estimation techniques,” _Intell. Converged Netw._ , vol. 3, no. 1, pp. 1–32, Mar. 2022.
  * [4] Q. Li, M. Wen, S. Wang, G. C. Alexandropoulos, and Y.-C. Wu, “Space shift keying with reconfigurable intelligent surfaces: Phase configuration designs and performance analysis,” _IEEE Open J. Commun. Society_ , vol. 2, pp. 322–333, Feb. 2021.
  * [5] G. C. Alexandropoulos, N. Shlezinger, and P. del Hougne, “Reconfigurable intelligent surfaces for rich scattering wireless communications: Recent experiments, challenges, and opportunities,” _IEEE Commun. Mag._ , vol. 59, no. 6, pp. 28–34, Jun. 2021.
  * [6] E. Calvanese Strinati, G. C. Alexandropoulos, H. Wymeersch, B. Denis, V. Sciancalepore, R. D’Errico, A. Clemente, D.-T. Phan-Huy, E. De Carvalho, and P. Popovski, “Reconfigurable, intelligent, and sustainable wireless environments for 6G smart connectivity,” _IEEE Commun. Mag._ , vol. 59, no. 10, pp. 99–105, Oct. 2021.
  * [7] G. C. Alexandropoulos, D.-T. Phan-Huy, K. D. Katsanos, M. Crozzoli, H. Wymeersch, P. Popovski, P. Ratajczak, Y. Bénédic, M.-H. Hamon, S. Herraiz Gonzalez, P. Mursia, M. Rossanese, V. Sciancalepore, J.-B. Gros, S. Terranova, G. Gradoni, P. Di Lorenzo, M. Rahal, B. Denis, R. D’Errico, A. Clemente, and E. Calvanese Strinati, “RIS-enabled smart wireless environments: Deployment scenarios, network architecture, bandwidth and area of influence,” _EURASIP J. Wireless Commun. Netw._ , to appear, 2023\.
  * [8] C. Han, Y. Wang, Y. Li, Y. Chen, N. A. Abbasi, T. Kürner, and A. F. Molisch, “Terahertz wireless channels: A holistic survey on measurement, modeling, and analysis,” _IEEE Communications Surveys & Tutorials_, vol. 24, no. 3, pp. 1670–1707, 2022.
  * [9] “TERahertz ReconfigurAble METAsurfaces for ultra-high rate wireless communications (TERRAMETA) SNS JU Project,” https://terrameta-project.eu/.
  * [10] K. Keykhosravi, B. Denis, G. C. Alexandropoulos, Z. S. He, A. Albanese, V. Sciancalepore, and H. Wymeersch, “Leveraging RIS-enabled smart signal propagation for solving infeasible localization problems,” _IEEE Veh. Technol. Mag._ , vol. 18, no. 2, pp. 20–28, Jun. 2023.
  * [11] B. A. Çelebi, I. Aras, K. E. Özçelik, K. Tekbiyik, G. Kurt, O. Özdemır, and A. R. Ekti, “Accurate modelling of reconfigurable intelligent surfaces in THz band,” in _Proc. 29th Signal Process. Commun. Appl. Conf._ , Istanbul, Turkey, Jun. 2021.
  * [12] K. Rasilainen, T. D. Phan, M. Berg, A. Pärssinen, and P. J. Soh, “Hardware aspects of Sub-THz antennas and reconfigurable intelligent surfaces for 6G communicationss,” _IEEE J. Sel. Areas Commun._ , vol. 41, no. 8, pp. 2530–2546, Aug. 2023.
  * [13] S. Saab, S. A. Surra, G. Xu, and R. W. Heath, “An inverted $L$-shaped multi-layer reconfigurable intelligent surface for THz communications,” in _Proc. Asilomar Conf. Signals, Sys., Comp._ , Pacific Grove, USA, Jun. 2022\.
  * [14] D. Jun, Y. Youn, C. Lee, M. Hwang, W. Hong, and C. B. Chae, “Real-time implementation of semi-active reconfigurable intelligent surfaces for mmWave and Sub-THz systems,” in _Proc. IEEE Consumer Commun. Netw. Conf._ , Las Vegas, USA, Jan. 2023.
  * [15] W. Li, H. Tu, Y. He, L. Zhang, S.-W. Wong, and S. Gao, “A novel wideband tightly coupled dual-polarized reflectarray antenna,” _IEEE Trans. Antennas Propag._ , vol. 71, no. 6, pp. 5422–5427, Jun. 2023.
  * [16] J. Wang, Y. Zhou, S. Gao, and Q. Luo, “An efficiency-improved tightly coupled dipole reflectarray antenna using variant-coupling-capacitance method,” _IEEE Access_ , vol. 8, pp. 37 314–37 320, Feb. 2020.
  * [17] H. Kamoda, T. Iwasaki, J. Tsumochi, T. Kuki, and O. Hashimoto, “60-GHz electronically reconfigurable large reflectarray using single-bit phase shifters,” _IEEE Trans. Antennas Prop._ , vol. 59, no. 7, pp. 2524–2531, Jul. 2011.
  * [18] J. Han, L. Li, G. Liu, Z. Wu, and Y. Shi, “A wideband 1 bit 12×12 reconfigurable beam-scanning reflectarray: Design, fabrication, and measurement,” _IEEE Antennas Wireless Propag. Lett._ , vol. 18, no. 6, pp. 1268–1272, Jun. 2019.
  * [19] Q. Luo, S. Gao, W. Li, X. Yang, and G. Wen, “Ultra-wideband and multiband reflectarrays for intelligent multi-functional platforms,” in _Proc. European Conf. Antennas Propag._ , Krakow, Poland, Mar. 2019, pp. 1–5.
  * [20] Q. Luo, G. C. Alexandropoulos, Y. Zhou, and X. Ma, “Ultra-wideband metasurface at SubTHz: Hardware design and reflection optimization,” in _Proc. European Conf. Antennas Propag._ , Florence, Italy, Mar. 2020, pp. 1–5.
  * [21] M. Rahal, B. Denis, K. Keykhosravi, M. F. Keskin, B. Uguen, G. C. Alexandropoulos, and H. Wymeersch, “Performance of RIS-aided nearfield localization under beams approximation from real hardware characterization,” _EURASIP J. Wireless Commun. Netw._ , to appear, 2023.
  * [22] V. Jamali, G. C. Alexandropoulos, R. Schober, and H. V. Poor, “Low-to-zero-overhead IRS reconfiguration: Decoupling illumination and channel estimation,” _IEEE Commun. Lett._ , vol. 16, no. 4, pp. 932–936, Apr. 2022.
  * [23] G. C. Alexandropoulos, V. Jamali, R. Schober, and H. V. Poor, “Near-field hierarchical beam management for RIS-enabled millimeter wave multi-antenna systems,” in _Proc. Sensor Array Multich. Signal Process Works._ , Trondheim, Norway, Jun. 2022.
  * [24] J. An, C. Xu, Q. Wu, D. W. K. Ng, M. Di Renzo, C. Yuen, and L. Hanzo, “Codebook-based solutions for reconfigurable intelligent surfaces and their open challenges,” _IEEE Wireless Commun._ , to appear, 2023.
  * [25] G. C. Alexandropoulos, I. Vinieratou, and H. Wymeersch, “Localization via multiple reconfigurable intelligent surfaces equipped with single receive RF chains,” _IEEE Wireless Commun. Lett._ , vol. 11, no. 5, pp. 1072–1076, May 2022.
  * [26] S. P. Chepuri, N. Shlezinger, F. Liu, G. C. Alexandropoulos, S. Buzzi, and Y. C. Eldar, “Integrated sensing and communications with reconfigurable intelligent surfaces,” _IEEE Signal Process. Mag._ , vol. 40, no. 6, pp. 41–62, Sep. 2023.


# Study of Firecracker MicroVM

Madhur Jain  Khoury College of Computer Sciences
Northeastern University
Boston, MA
<EMAIL_ADDRESS>
###### Abstract

Firecracker is a virtualization technology that makes use of Kernel Virtual
Machine (KVM). Firecracker belongs to a new virtualization class named the
micro-virtual machines (MicroVMs). Using Firecracker, we can launch
lightweight MicroVMs in non-virtualized environments in a fraction of a
second, at the same time offering the security and workload isolation provided
by traditional VMs and also the resource efficiency that comes along with
containers [1]. Firecracker aims to provide a slimmed-down MicroVM, comprised
of approximately 50K lines of code in Rust and with a reduced attack surface
for guest VMs. This report will examine the internals of Firecracker and
understand why Firecracker is the next big thing going forward in
virtualization and cloud computing.

###### Index Terms:

firecracker, microvm, Rust, VMM, QEMU, KVM

##  I Introduction

Firecracker is a new open source Virtual Machine Monitor (VMM) developed by
AWS for serverless workloads. It is a virtualization technology which makes
use of KVM, meaning it can only be run on KVM-supported and enabled hosts and
a corresponding Linux Kernel v4.14+. With the recommended Linux kernel guest
configuration, Firecracker claims to offer a memory overhead of less than 5MB
per container, boots to application code within 125ms, and allows for the
creation of up to 150 MicroVMs per second. The number of Firecracker MicroVMs
running simultaneously on a single host is only limited by the availability of
hardware resources. [1]

Firecracker provides security and resource isolation by running the
Firecracker userpsace process inside a jailer process. The jailer sets up
system resources that require elevated permissions (e.g., cgroup, chroot),
drops privileges, and then exec()s into the Firecracker binary, which then
runs as an unprivileged process. Past this point, Firecracker can only access
resources to which a privileged third-party grants access (e.g., by copying a
file into the chroot, or passing a file descriptor). [13]

Seccomp filters limit the system calls that the Firecracker process can use.
There are 3 possible levels of seccomp filtering, configurable by passing a
command line argument to the jailer: 0 (disabled), 1 (whitelists a set of
trusted system calls by their identifiers) and 2 (whitelists a set of trusted
system calls with trusted parameter values), the latter being the most
restrictive and the recommended one. The filters are loaded in the Firecracker
process, immediately before the execution of the untrusted guest code starts.
[13]

Firecracker was developed to handle serverless workloads and has been running
in production for AWS Lambda and Fargate since 2018. Serverless workloads
require isolation and security, and at the same time, have container benefits
of faster boot time. A lot of research has been done with respect to the cold-
start and warm-start of the VM instances for serverless workloads. [17] The
idea behind developing Firecracker was to make use of the KVM module loaded in
the Linux kernel and get rid of the legacy devices that other virtualization
technologies like Xen and VMWare offer. This way, Firecracker can create a VMM
with a smaller memory footprint and also provide improved performance.

Section 2 explains the high-level architecture of the Firecracker Micro VM,
Section 3 dives deep into the boot sequence of the Firecracker, Section 4
describes the device model emulation provided, and Section 5 shares some light
on the conclusions garnered through this study.

##  II Architecture of Firecracker

KVM is an enabler of hardware extensions provided by vendors such as Intel and
AMD with their virtualization extensions such as SVM and VMX. These extensions
allow KVM to directly execute the guest code on the host CPU. There are three
sets of ioctls that make up the KVM API and are issued to control the various
aspects of the virtual machine. The three classes that the iocltls belongs to
are [7] -

  * •

System IOCTLs: These query and set global attributes, which affect the whole
KVM subsystem. In addition, a system ioctl is used to create virtual machines.

  * •

VM IOCTLs: These query and set attributes that affect an entire virtual
machine—for example, memory layout. In addition, a VM ioctl is used to create
virtual CPUs (vCPUs). VM ioctls are run from the same userspace process
(address space) that was used to create the VM.

  * •

vCPU IOCTLs: These query and set attributes that control the operation of a
single virtual CPU. They run vCPU ioctls from the same thread that was used to
create the vCPU.

Figure 1: Firecracker Design - Threads

QEMU is an open source machine emulator and a virtualization technique to run
KVM-enabled Virtual Machines. Similar to QEMU, Firecracker uses a
multithreaded event-driven architecture. Each Firecracker process runs one and
only one MicroVM. The process consists of three main threads: API, VMM and
vCPU. The API server thread is an HTTP server that is created to accept
requests and performs actions on those requests. The API server thread
structure comprises of a socket connection to listen to requests on the port,
epoll_fd to listen to events on the socket port, and a hashmap of connections.
The hashmap consists of token and connection pairs corresponding to the file
descriptor of the stream—in this case, the socket.

Traditionally, QEMU has made use of select or poll system calls to maintain an
event loop of file descriptors on which to listen for new events. [16] select
or poll system calls requires a list of all open file descriptors maintained
in the VMM structure and then it would go through each of the file descriptors
to determine which FDs have new events. This would take up O(N) time where N
is the number of file descriptors to listen for new events. [14]

Firecracker takes the epoll approach where the host kernel maintains a list of
file descriptors for VM process and notifies the VMM whenever there is a new
event that occurs in any of the file descriptors. This is called as a Edge
Triggered mechanism (”pull”), whereas the select/poll was a Level Triggered
mechanism (”push”). The epoll_fd structure created by Firecracker has the
‘close_on_exec‘ flag set, which means if a process forks the Firecracker
process, the file descriptors would not be shared.

The API server exposes a number of requests as REST APIs which can be used to
configure the MicroVM. Once the MicroVM has been started, i.e., once it
receives the ”InstanceStart” command, the API server thread will just block on
the epoll file descriptor until new events come in. Firecracker creates a
channel (see Rust channels) to enable communication between the API Server
thread and the VMM thread. Rust channels are similar to Unix pipes for
comparison.

The VMM server thread manages the entire MicroVM. Once the VMM server thread
is created, it runs an event loop which takes the parsed request one by one
from the API server thread and dispatches it to the appropriate handlers. The
handlers are defined according to the dispatch table set by the event loop.
For now, Firecracker supports the following handlers - Exit, Stdin,
DeviceHandler, VMMActionRequest, WriteMetrics. The dispatch table is managed
by the epoll fd. The dispatch table maintains a map of file descriptors that
are to be monitored, and the kind of events to be monitored for. When the vCPU
thread creation request is received by the VMM thread, the VMM spawns the
required number of vCPU by using the KVM vCPU ioctls. A vCPU thread is created
for each virtual CPU. These vCPU threads are nothing but POSIX threads created
by KVM. To run guest code, the vCPUs execute an ioctl with KVM_RUN as its
argument.

Software in the root mode is the hypervisor. Hypervisor or the VMM forms a new
plane that runs in root mode while the VMs run in non-root mode. KVM uses
virtualization extensions to provide these different modes on the host CPUs.
In the case of Intel CPUs, VT-x is the CPU virtualization and VT-d is the IO
virtualization. For vCPUs, VT-x provides two modes of guest code execution:
root and non-root. Whenever a VM attempts to execute an instruction that is
prohibited by the non-root model, vCPU immediately switches to a root mode in
a trap-like way. This is called a VMEXIT. Hypervisor deals with the reason for
VMEXIT and then executes VMRESUME to re-enter non-root mode for that VM
instance. This interaction between root and non-root is the essence of
hardware virtualization.

##  III Firecracker Boot Sequence

Traditional PCs boot Linux with a BIOS and a bootloader. The primary
responsibilities of the BIOS includes booting the CPU in real mode and
performing a Power on Self Test (POST) setup before loading the bootloader.
BIOS determines the candidates for boot devices, and once a bootable device is
found, the bootloader is loaded into RAM and executed. Different systems have
different stages of the bootloader to be executed. LILO, for example, has a
two stage bootloader while GRUB contains a 3 stage bootloader. Multiple stages
of a bootloader are used because of the system limitations of some of the
older devices that were used to boot Linux.

Linux kernels actually do not require a BIOS and a bootloader. Instead,
Firecracker uses what is known as Linux Boot Protocol. [15] There are multiple
versions of the Linux Boot Protocol standard that exist. Firecracker follows
the 64-bit Linux Boot Protocol Standard. Thus, Firecracker can directly load
the Linux kernel and mount the corresponding root file system.

The Linux kernel is an uncompressed bzImage (big compressed image, usually
larger than 512KB). The bzImage format consists of a real mode kernel code and
a protected mode kernel code. Instead of booting into the entry point defined
by the real mode, Firecracker directly boots into the 64-bit entry point
located at 0x200 in the protected mode kernel code. Firecracker loads the
uncompressed Linux kernel as well as the init process, thereby saving
approximately 20 to 30ms of the time taken to decompress the kernel.

                  ~                        ~
                  |  Protected-mode kernel |
          100000  +------------------------+
                  |  I/O memory hole       |
          0A0000  +------------------------+
                  |  Reserved for BIOS     |      Leave as much as possible unused
                  ~                        ~
                  |  Command line          |      (Can also be below the X+10000 mark)
          X+10000 +------------------------+
                  |  Stack/heap            |      For use by the kernel real-mode code.
          X+08000 +------------------------+
                  |  Kernel setup          |      The kernel real-mode code.
                  |  Kernel boot sector    |      The kernel legacy boot sector.
          X       +------------------------+
                  |  Boot loader           |      <- Boot sector entry point 0000:7C00
          001000  +------------------------+
                  |  Reserved for MBR/BIOS |
          000800  +------------------------+
                  |  Typically used by MBR |
          000600  +------------------------+
                  |  BIOS use only         |
          000000  +------------------------+

Figure 2: Linux BzImage Memory Layout

Linux kernel also contains another component, namely the initramfs [5]. There
are four primary reasons to have an initramfs in the LFS environment: loading
the rootfs from a network, loading it from an LVM logical volume, having an
encrypted rootfs where a password is required, or for the convenience of
specifying the rootfs as a LABEL or UUID. Anything else usually means that the
kernel was not configured properly.

Since Firecracker doesn’t need any of the above stated reasons for loading the
initramfs before mounting the root file system, it is recommended to avoid
loading the initramfs at boot time, thereby further reducing the overall boot
time and the memory footprint for the kernel. So, if no initramfs is
configured externally, then at boot time, Firecracker replaces the initramfs
with a default empty, 134 byte initramfs.

##  IV Firecracker Device Model

Until this section, we have talked about the similar architectures and
execution flow for Firecracker and QEMU. So, what is different between QEMU
and Firecracker? One of the main differences is with the device emulations.
There are only 5 Device emulations available in Firecracker: network, block
devices, sockets, serial console and minimal keyboard controller, as shown in
Figure 3. Firecracker does not provide support for device emulations like USB,
GPU and 9P filesystem in order to provide increased security compared to other
virtualization technologies like QEMU. On the other hand, QEMU has most device
model emulations available in the VMM. The careful reader will notice that
Firecracker does not the use the vhost implementation in the host kernel that
provides more efficient IO performance without doing VMEXITS.

An open specification for emulating device models in virtualization has been
developed, named Virtio. Virtio is defined as a straightforward, efficient and
a standard mechanism to allow guest OS to talk to the virtual device driver in
a similar way the host OS would call the actual hardware device driver. It
takes advantage of the fact that the guest can share memory with the host for
IO.

The general flow for the virtio specification [8] includes a front-end driver
representing the virtual device in the guest, and the corresponding device
being exposed by the hypervisor or the VMM. A transport layer enables
communication between the host and the guest. For the transport layers, Virtio
employs a ring-buffer virtqueue structure. A virtuqueue is a queue of guest-
allocated buffers that the host interacts with either by reading or writing to
them. Each device can have zero or more virtqueues. A back-end driver present
in the host kernel completes the communication flow, to which the virtqueue is
connected.

Firecracker device model architecture using Virtio is shown in Figure 3. The
following list provides a description of the devices available within
Firecracker:

  * •

virtio-net: implementation for the network driver (tun/tap devices)

  * •

virtio-blk: implementation for the block devices

  * •

virtio-vsock: implementation for VM sockets providing N:1 serial
communications

  * •

serial console: implementation for the legacy console devices for serial
communication - terminal

  * •

keyboard controller: implementation for the keyboard device, though only one
function is implemented - Ctrl+Alt+Del to reboot/shutdown the system.

Figure 3: Firecracker Device Model

##  V Scope for Improvements

Even with all the excellent features providing near-native performance of the
guest code using KVM, as well as faster boot times and lower memory footprint
of the VMs due to fewer support for available device emulations, there still
exist some areas for improvements that can make Firecracker more suitable for
general use cases and not just for serverless workloads.

  * •

Support for virtio-fs: virto-fs is the interface to provide efficient sharing
between the host and the guest filesystem avoiding context switches (VMEXITS)
thereby providing more performance. virtio-fs is an upgrade on the existing
virtio-9p interface for the same purpose. Though more research is required for
security purposes before including it as part of Firecracker.

  * •

Increased IO Performance: The results of the tests performed between the
Firecracker, QEMU and Cloud Hypervisor show limitations in Firecracker’s
virtio implementation and serial execution. [1] [2]

  * •

Larger number of device emulations: Currently, Firecracker can emulate only 10
devices, since each device gets its own IRQ. [10]

  * •

Support for attaching devices at runtime: Firecracker only allows specifying
the devices at booting time. Devices can only be attached when the MicroVM is
shut down.

  * •

Hotplugging Support: For any workload, it is beneficial to allow guest
memory/CPU hotplugging within a VM at runtime in order to avoid interference
to the workload. Firecracker oversubscribes the allocated memory required for
the guest, but there is no way to expand the allocated memory for the guest.

  * •

Memory Ballooning Support: At present, Firecracker does not have any support
for reclaiming unused memory from the guest, since no communication is present
between the host and the guest. This, along with the hotplugging feature would
make it very easy to dynamically add/remove memory/CPU at runtime thereby
providing elasticity to the MicroVM. [9]

##  VI Conclusion/Thoughts

This paper reviews the implementation of a minimalist and modular VMM in the
form of Firecracker MicroVM. It also identifies the process of how Firecracker
provides resource isolation and security through the use of seccomp filters
and jailer process and provides faster boot times and lower memory footprint
due to KVM and minimal device model emulation.

One other thing to note is that Firecracker embodies the modular design in the
development of the hypervisor. The modular design approach has also led to the
development of community-driven high-quality rust-vmm crates which provide us
with the core modules required for the implementation of a hypervisor [11].
rust-vmm [12] is a community approach initiated by AWS. Amazon along with
Intel, Redhat and Google, are trying to provide a platform to build a
hypervisor from scratch by only consuming the modules required from the rust-
vmm crates. This approach also enables the development of a plug-n-play
architecture in hypervisors, which we haven’t seen so far.

## References

  * [1] https://aws.amazon.com/blogs/aws/firecracker-lightweight-virtualization-for-serverless-computing
  * [2] Alexandru Agache, Marc Brooker, Andreea Florescu, Alexandra Iordache, Anthony Liguori, Rolf Neugebauer, Phil Piwonka, and Diana-Maria Popa https://www.usenix.org/system/files/nsdi20-paper-agache.pdf Firecracker: Lightweight Virtualization for Serverless Applications
  * [3] https://unixism.net/2019/10/how-aws-firecracker-works-a-deep-dive/
  * [4] https://developer.ibm.com/technologies/linux/articles/l-linuxboot/
  * [5] http://www.linuxfromscratch.org/blfs/view/cvs/postlfs/initramfs.html
  * [6] Swift Birth and Quick Death: Enabling Fast ParallelGuest Boot and Destruction in the Xen Hypervisor https://ssrg.ece.vt.edu/papers/vee_2017.pdf
  * [7] Mastering KVM Virtualization by Prasad Mukhedkar; Humble Devassy Chirammal; Anil Vettathu https://learning.oreilly.com/library/view/mastering-kvm-virtualization/9781784399054/ch02s04.html
  * [8] Deep Dive into virtio and vhost-net driver https://www.redhat.com/en/blog/deep-dive-virtio-networking-and-vhost-net
  * [9] Memory Ballooning Support https://github.com/firecracker-microvm/firecracker/issues/1571
  * [10] Enable IRQ Sharing to remove current limit of 10 devices https://github.com/firecracker-microvm/firecracker/issues/1268
  * [11] Building virtualization stack for the future https://opensource.com/article/19/3/rust-virtual-machine
  * [12] rust-vmm https://github.com/rust-vmm
  * [13] Firecracker Design Doc https://github.com/firecracker-microvm/firecracker/blob/master/docs/design.md
  * [14] Understand epoll and its madness https://medium.com/@copyconstruct/the-method-to-epolls-madness-d9d2d6378642
  * [15] Linux 64-bit Boot Protocol https://www.kernel.org/doc/html/latest/x86/boot.html#id1
  * [16] QEMU Internals http://blog.vmsplice.net/2011/03/qemu-internals-overall-architecture-and.html
  * [17] Cold Start/Warm Start with AWS Lambda https://blog.octo.com/en/cold-start-warm-start-with-aws-lambda/


# Structural Neural Additive Models:
Enhanced Interpretable Machine Learning

Mattias Luber  Anton Thielmann  Benjamin Säfken

###### Abstract

Deep neural networks (DNNs) have shown exceptional performances in a wide
range of tasks and have become the go-to method for problems requiring high-
level predictive power. There has been extensive research on how DNNs arrive
at their decisions, however, the inherently uninterpretable networks remain up
to this day mostly unobservable "black boxes". In recent years, the field has
seen a push towards interpretable neural networks, such as the visually
interpretable Neural Additive Models (NAMs). We propose a further step into
the direction of intelligibility beyond the mere visualization of feature
effects and propose Structural Neural Additive Models (SNAMs). A modeling
framework that combines classical and clearly interpretable statistical
methods with the predictive power of neural applications. Our experiments
validate the predictive performances of SNAMs. The proposed framework performs
comparable to state-of-the-art fully connected DNNs and we show that SNAMs can
even outperform NAMs while remaining inherently more interpretable.

Machine Learning, ICML

##  1 Introduction

In the recent years deep neural networks (DNNs) turned out to be extremely
effective in modelling difficult high dimensional problems and in recognizing
complex patterns in unstructured data. While the reasons for their performance
is yet not fully understood, it is undoubted that DNNs form the current state-
of-the-art-approach in many disciplines such as computer vision, machine
translation and speech recognition (Sejnowski, 2020). Their effectiveness
however, comes at a price. Due to their "black box" nature, DNNs lack
interpretability, which limits their application in data sensitive domains
like health care, finance or insurance (Adadi & Berrada, 2018). To address
this limitations, several approaches have been developed based on post-hoc
investigations of the model predictions. Common examples here are LIME
(Ribeiro et al., 2016) or layer wise relevance propagation (LRP) (Bach et al.,
2015). Although those approaches might be able to indicate how individual
predictions are generated, they do not provide a global and complete picture
of the underlying decision making process. Neural Additive Models (NAMs)
(Agarwal et al., 2021b) were recently proposed as a class of Neural Networks,
that impose an additivity constraint on the input data and thus allow to
directly derive the feature-wise contribution onto the generated predictions
as a function of the input domain. While this indeed yields an exact
representation of the decision making process, NAMs are nevertheless highly
complex functions that are characterized by hundreds of thousands of
parameters and thus fail to address additional dimensions of interpretability
(Murdoch et al., 2019). To this end, we propose the use of Structural Neural
Additive Models (SNAMs) as a way to achieve the same (and even better)
predictive performance with a fraction of the parameters required, while
providing intelligibility beyond mere visualizations. The contributions of
SNAMs can be summarized as follows:

  * •

SNAMs enable direct interpretation of the estimated parameters, as well as a
measure of model uncertainty based on Bayesian confidence intervals.

  * •

Parameter sparsity is achieved through spline-based activation layers that
mimic the behaviour of cubic regression splines and their ability to model
nonlinear functions with a minimum of parameters.

  * •

We propose fast and computationally efficient neural splines based on the
Silverman kernel.

  * •

Due to the extreme parameter sparsity, intelligibility beyond visualization is
possible and quantification of parameter uncertainty and hence the
construction of confidence bands is achieved.

  * •

Learnable knots allow SNAMs to model highly non-linear functional
relationships turning the spline based approach into a highly flexible model
class, despite the extreme parameter sparsity.

  * •

A special hallmark of the spline based activation layers is, that they
directly generalize to higher dimensions and thereby form a native way to
process and interpret spatial data. The effectiveness of this property is
demonstrated on the Carlifornia Housing dataset along with an evaluation of
the full interpretability of SNAMs.

  * •

The predictive performance of SNAMs is evaluated with a benchmark study. We
analyze multiple datasets and compare the results to multiple neural- as well
as non-neural models.

The paper is structured as follows: First, we shortly summarize generalized
additive models (Hastie & Tibshirani, 1986) and traditional smoothers.
Subsequently, we reformulate the basis expansion with cubic regression splines
as a neural activation layer and then propose a fast and computationally
efficient approximation based on the Silverman kernel. Additionally, we
introduce learnable knots and thus adapt neural splines to account for jagged
shape functions. Section 3 demonstrates the interpretability of SNAMs via the
construction of the Fisher Information Matrix. In section 4 we conduct the
benchmark study and compare SNAMs to other state-of-the-art models. Section 5
concludes.

##  2 Structural Neural Additive Models

Generalized Additive Models (GAM) were developed by Hastie and Tibshirani
(Hastie & Tibshirani, 1986) to model flexible nonlinear effects in regression
and classification problems as an extension of the more restrictive
Generalized Linear Models (GLM) (Nelder & Wedderburn, 1972). In its
fundamental form, a GAM can be expressed as

$\eta=\beta_{0}+\sum_{i}^{p}f_{i}(x_{i})+\sum_{i,j}^{p}f_{ij}(x_{i},x_{j})+\epsilon,\mu=g(\eta),$
(1)

where $\eta$ is used to model the expectation of the target variable $y$ given
some features $x\in\mathbb{R}^{(n,p)}$ via link-function $g(\eta)=\mu\equiv
E(y|x)$ and $y\in\mathbb{R}^{(n,1)}$ is further assumed to follow a
distribution of the exponential family.

Figure 1: Visualization of the cubic regression splines. The left hand side
displays the basis functions and the right hand side shows how the function
approximation is done by scaling the bases appropriately

For the construction of the smoothers, $f_{i}(x_{i})$, a basis expansion with
nonlinear bases is applied and the final function is represented as their
linear combination:

$f(x_{i})=\sum_{j}^{k}B_{j}(x_{i})\beta_{j}.$ (2)

The location of the basis functions, $B_{j}$, is parametrized by an ordered
set of knots $\kappa_{1}<\kappa_{2}<...<\kappa_{k}$, where the positioning is
traditionally chosen to be either uniformly distributed over the input domain
or located at the positions of the quantiles. The later one allows the model
to be more flexible in regions with more data and is more restrictive in
sparse areas.

To control the smoothers “wigglyness” and hence overfitting their squared
second derivative can be penalized. This allows a trade-off between fit to the
data and good generalization. However, due to their smoothness the splines are
by design already robust against overfitting and since they are taking their
local neighborhood into account, they also yield good results on respectively
small datasets (Wood, 2017). Parameter estimation is hence done by penalized
likelihood maximization, which aims to solve the following optimization
problem:

$\underset{\theta}{\text{arg
min}}-\frac{1}{n}l(\theta,x,y)+\sum_{j}^{k}\lambda_{j}\int
f^{\prime\prime}(x)^{2}dx,$ (3)

with $\theta=\left[\beta_{1},\dots,\beta_{k}\right]$ with penalized iterative
reweighted least squares (PIRLS). $l(\theta,x,y)$ defines the log-likelihood
and $\lambda_{j}$ determines the impact of the penalization of the
“wigglyness”. Further hyperparameters are the positions of the knots $\kappa$,
as well as the number of basis functions. Theoretically a higher number of
bases leads to more flexible functions. In practice, the dependence on the
number of bases is only moderately big, as it is compensated by the
penalization (Wood, 2017).

Figure 2: The figure sketches how the different types of approximations could
be implemented as neural networks.

###  2.1 Basis Functions and Activation Functions

To model explainable terms in neural networks with smooth functions, we
propose spline units according to equation 2. Hence, a unit consists of a
basis expansion layer in connection with an ordinary linear layer to form the
weighted sum. For simplicity, we introduce the model framework based on the
cubic spline unit, where the basis expansion layer is constructed around the
$\mathbf{G}$ matrix similarly as the cubic expansion in (Wood, 2017):

$\mathbf{G}=\begin{pmatrix}\mathbf{0}\\\ \mathbf{\mathbf{B}^{-1}}\mathbf{D}\\\
\mathbf{0}\end{pmatrix},$ (4)

with $h_{j}=\kappa_{j+1}-\kappa_{j}$,

$\displaystyle\mathbf{B}=\begin{pmatrix}\frac{1}{2}(h_{1}+h_{2})&\frac{1}{6}h_{2}&0&\dots&0\\\
\frac{1}{6}h_{2}&\frac{1}{3}(h_{2}+h_{3})&\frac{1}{6}h_{3}&\dots&0\\\
0&\ddots&\ddots&&\vdots\\\ \vdots&&&\ddots&\frac{1}{6}h_{k-2}\\\
0&\dots&0&\frac{1}{6}h_{k-2}&\frac{1}{3}(h_{k-1}+h_{k}-1)\end{pmatrix}$

and

$\displaystyle\mathbf{D}=\begin{pmatrix}\frac{1}{h_{2}}&-\frac{1}{h_{1}}-\frac{1}{h_{2}}&\frac{1}{h_{2}}&0&\dots&0\\\
0&\frac{1}{h_{2}}&-\frac{1}{h_{2}}-\frac{1}{h_{3}}&\frac{1}{h_{3}}&&\vdots\\\
\vdots&&\ddots&\ddots&\ddots&\\\
0&\dots&0&\frac{1}{h_{k-2}}&-\frac{1}{h_{k-2}}-\frac{1}{h_{k-1}}&\frac{1}{h_{k-1}}\end{pmatrix}.$

Based on that, the expansion is a function
$B:\mathbb{R}\rightarrow\mathbb{R}^{k}$ with the explicit formulation:

$\begin{split}B(x)=(c^{-}_{j}\mathbf{G}_{j,\bullet}+a^{-}_{j}e_{j}+c^{+}_{j}\mathbf{G}_{j+1,\bullet}+a_{j}^{+}e_{j+1})\end{split}$
(5) $\displaystyle
c^{-}_{j}=\frac{(\kappa_{j+1}-x)^{3}}{h_{j}}-\frac{h_{j}(\kappa_{j+1}-x)}{6}$
$\displaystyle a_{j}^{-}=\frac{(k_{j+1}-x)}{h_{j}}$ $\displaystyle
c^{+}_{j}=\frac{(x-\kappa_{j})^{3}}{h_{j}}-\frac{h_{j}(x-\kappa_{j})}{6}$
$\displaystyle a_{j}^{+}=\frac{(x-k_{j})}{h_{j}}$ $\displaystyle
h_{j}=\kappa_{j+1}-\kappa_{j}$ $\displaystyle j\text{ s.t.
}\kappa_{j}<x<\kappa_{j+1}$

for values $x\in[\kappa_{1},\kappa_{k}]$, whereby the extrapolation is handled
separately. Furthermore, $e_{j}$ represents the $j$-th standard basis.
Eventually $B(x)$ is fully defined by the knot positions and thus it can be
implemented straightforward as a parametrized layer with fixed knots. In case
of the penalized splines the layer is further equipped with the smoothing
matrix $\mathbf{S}=\mathbf{D}^{T}\mathbf{B}^{-1}\mathbf{D}$ which can be used
to determine the penalty.

####  2.1.1 Fast Spline Units

The cubic splines represent a coherent way to achieve nonlinearity in a model
and come with convenient properties such as a guaranteed smoothness of the
resulting function, a good fit and implicit penalization. In context of deep
learning it could be viewed as a special activation layer, but in contrast to
the established activations such as ReLU or Sigmoid, the gradient calculation
is more expensive. This section generalizes the idea of basis expansion to
ordinary activation functions and hence shows the generalizability of the
presented approach. Fast approximations of the cubic spline layer based on the
Silverman kernel (Silverman, 1984) are derived. Since polynomial splines can
be represented as a truncated power series of degree $d$, (2) can be expressed
as a polynomial part and a truncated polynomial part, with
$\kappa_{1}<\dots<\kappa_{k}$ representing a sequence of knots (Fahrmeir et
al., 2013):

$\begin{split}f(x)=\beta_{0}+\beta_{1}x+\dots+\beta_{d}x^{d}+\beta_{d+1}(x-\kappa_{2})_{+}^{d}+\dots+\\\
\beta_{k+d-2}(x-\kappa_{k-1})_{+}^{d}\end{split}$ (6)

with

$(x-\kappa)_{+}^{d}=\begin{cases}(x-\kappa)^{d},&\text{if }x\geq\kappa\\\
0,&\text{otherwise}.\end{cases}$ (7)

Especially for the truncated part the idea is the same as for the cubic
splines in a way that the input domain is populated with shifted versions of
the basis function.

Figure 3: Basis expansion based on different types of activation functions
along with their approximations of an exemplary function.

The essential difference is that truncated polynomials are used instead of the
bell-shaped cubic bases. The approximation is then again performed by
searching for an appropriate scaling, such that the bases sum up roughly to
the function of interest (see figure 3). For $d=1$ the truncated polynomials
are actually equivalent to shifted ReLU functions and thus the exact same
behaviour could be accomplished by shifting the data rather than the function
and passing it through an ordinary ReLU. From that point it is easy to see how
the idea of basis expansion can be generalized to various types of
activations, which is also visualized in figure 3. To get closer to the cubic
splines the Silverman kernel is used as an activation function.

$K(x)=\frac{1}{2}e^{\frac{-|x|}{\sqrt{2}}}\text{sin}\left(\frac{|x|}{\sqrt{2}}+\frac{\pi}{4}\right)$
(8)

It is considered to be the kernel method that is equivalent to cubic smoothing
splines (Silverman, 1984) and thus has the same characteristic bell shape as
the cubic spline bases. One difference is that the width of the basis is not
determined by some knot positions and has to be handled by adjusting the
bandwidth of the kernel, or respectively scaling the data. More precisely, the
data is passed through the kernel as $K(\frac{x-d}{\sigma})$ with $d$ as the
displacement and $\sigma$ as the scaling. Other properties that are lost by
approximating the cubic spline unit with Silverman kernels, are the implicit
control of the smoothness and the linear extrapolation, but on the other hand
the calculation of the gradients becomes easier and furthermore the position
of each kernel as well as its bandwidth can easily be set to be learnable
parameters.

####  2.1.2 Learnable Knots

Modeling jagged shape functions and accurately representing jumps is an
advantage of deep neural networks because of their flexibility. Multi-layer
perceptrons can accurately represent jumps using an appropriate number of
neurons and appropriate activation functions. In contrast, (Caruana et al.,
2015) found that GAMs can overregulate and do not accurately represent details
in data. Flexible choice of node locations, however, in addition to an
appropriate number of knots, can easily overcome this problem (Mohanty &
Fahnestock, 2021; Gálvez & Iglesias, 2011). Since the expansion can be broken
down to a couple of simple arithmetic operations and matrix products in
respect to $\kappa$, backpropagation can be used to make the knot positions
trainable parameters. The only part that requires extra care is the assumption
that the knots are an ordered sequence $\kappa_{1}<...<\kappa_{k}$. Thus,
before each forward pass the knot parameters have to be sorted appropriately
and the $\mathbf{F}$ matrix has to be recalculated accordingly. In equation 5
it can be seen that the construction of $\mathbf{F}$ involves the division by
$h_{j}=\kappa_{j+1}-\kappa_{j}$, which can lead to numerical issues when two
neighbouring knots are getting too close and can be circumvented with a knot
distance penalty $g(\kappa)=\sum_{i}^{k-1}\frac{1}{h_{j}}$.

Figure 4: Turning the knot positions into learnable model parameters, enables
the SNAMs to fit jagged functions. The black marks on the x-axis indicate the
knot positions after training.

###  2.2 Identifiability

Without further restrictions, the identifiability is not fulfilled for GAMs
and hence not fulfilled for SNAMs. One could easily shift an additive constant
from one component $f_{i}(x)$ to another $f_{j}(x)$ without affecting the
predictive behaviour of the model at all. To prevent this, additional
constraints are imposed which ensure that

$\sum_{i}^{n}f(x_{1})=\dots=\sum_{i}^{n}f(x_{k})=0$ (9)

and all the additive constants are collected in a global intercept. Those
constraints do not affect the capabilities of the model, but ensure the
desired uniqueness and full interpretability of its components. One way of
implementing such a constraint is by subtracting the column means of the basis
expansion. To make the presented units identifiable an intermediate layer
between the expansion and the summation is implemented which tracks the
running column means of the expansion and subtracts it during each forward
pass.

##  3 Enhancing interpretability with SNAMs

One major advantage of SNAMs is that they offer a toolbox for improving
intelligibility of the model while retaining the high-accuracy of state-of-
the-art DL models. While NAMs turn black-box DL models into glass-box models
mainly by visualization, SNAMs can be considered as sand-box models. This
highlights the fact that SNAMs can be used for a post-hoc analysis widening
the understanding of the underlying model’s behaviour and therefore the data-
generating process. These improvements can be connected to different
properties described in this section.

###  3.1 Low dimensional parameter space

Standard DL models and also recently proposed NAMs are notoriously
overparameterized. Owning their flexibility to the high dimensionality, they
often consist of hundreds of thousands or even millions of parameters. On the
contrary SNAMs considerably reduce the number of parameters to the amount of a
few hundred or so.

The stability of the fitting process is largely improved with less multiple
and more pronounced minima of the loss landscape, leading to more reliable
parameter estimates. Furthermore the parameters are intelligible as the
weights scale local basis function, see formulation (2), while the biases
define the knot positions. This allows the local behaviour of each single
feature to be investigated and for instance formally specified by a
derivative. This well-behaved parameterization is precisely the reason for the
openness of SNAMs to methods from the statistical community such as
calculating confidence bands, hypothesis testing or using AIC / BIC based
model selection.

SNAMs maintain the well calibrated estimates between over- and underfitting
for low to moderate sized datasets. This supersedes the need to distinguish
between “too small” and “large enough” data DL models are exposed to.

###  3.2 Understanding spatial effects

Although limiting the amount of parameters SNAMs can model interactions in
high-dimensional feature spaces if need be. This is especially useful for
spatial data. While in the NAM framework for the California Housing dataset
longitude and latitude are included separately (Agarwal et al., 2021a) both
can be included in an interaction term by SNAMs.
Two (or more) dimensional smoothers can be modelled as tensor product splines
by evaluating two univariate smoothers along each dimension and estimating a
coefficient for the combinations of their Cartesian product

$f_{ij}(x_{i},x_{j})=\sum_{l,r}^{p,q}B_{l}(x_{i})B_{r}(x_{j})\beta_{l,r}.$
(10)

The learned shape function $f_{ij}(x_{i},x_{j})$ can be visualized as a 2d
heat map for the California Housing data, see figure (10). The effect size is
given by the color with yellow indicting a strong positive effect and blue
indicting a strong negative effect. This gives human interpretable effect
strengths and provides a sense of confidence with regard to the $p\times q$
learned parameters.

Figure 5: The results of the SplineNAM on the CA Housing dataset displaying a
heatmap of the spatial effect.

###  3.3 Intelligibility beyond visualization

Intelligibility goes far beyond visualization. While with NAMs effects of
single features can be visualized, no further statements on the underlying
data generating process can be made. However with SNAMs further
intelligibility is available. This is due to the low dimensional parameter
space giving meaning to the parameters. This allows for the quantification of
the parameter uncertainty. As a result meaningful confidence bands can be
calculated and questions regarding the presence or absence of feature effects
can be answered withing the hypothesis testing framework. The key quantity to
asses the uncertainty is the Fisher information.

Fisher Information: The Fisher information matrix is one of the core concepts
in the maximum-likelihood framework. It is defined as the covariance of the
gradient of the log-likelihood. Hence the Fisher information is defined as

$\mathbf{F}(\theta)=E\left(\frac{\delta
l(\theta,x,y)}{\delta\theta}\frac{\delta
l(\theta,x,y)}{\delta\theta}^{T}\right).$ (11)

Furthermore under certain regularity conditions the Fisher information can be
expressed as second partial derivatives the log-likelihood, i.e. the Hessian
(Ly et al., 2017)

$\mathbf{H(\theta)}=-E\left(\frac{\delta^{2}l(\theta,x,y)}{\delta\theta\delta\theta^{T}}\right).$
(12)

Confidence Bands: In context of parameter estimation the Fisher information is
used together with the asymptotic distribution of

$\mathbf{\hat{\beta}}\overset{asy}{\sim}N\left(\beta,\mathbf{F}^{-1}(\beta)\right)$
(13)

to asses the uncertainty of the parameters $\hat{\beta}$. In the classic GLMs
$\alpha$-level confidence intervals are in general constructed as

$\left[\mathbf{\hat{\beta}}_{i}\pm
z_{1-\frac{\alpha}{2}}\sqrt{\mathbf{F^{-1}(\mathbf{\hat{\beta}})}_{ii}}\right]$
(14)

(Hastie et al., 2009) with $z_{1-\frac{\alpha}{2}}$ as the
$(1-\frac{\alpha}{2})$ quantile of the standard normal distribution. Even
though theoretically the confidence bands for the GAMs could be derived in the
same pointwise manner for each of the coefficients, in practice this leads to
poor coverage results since it neglects the structure between basis functions
and thus, the construction of Bayesian credibility bands is preferred (Wood,
2017). Hence, new coefficients are sampled from the posterior

$\beta^{*}\sim
N\left(\mathbf{\hat{\beta}},(\mathbf{F(\mathbf{\hat{\beta}}})+\lambda\mathbf{S})^{-1}\right)$
(15)

and the bands are represented by the quantiles of the resulting functions

$f_{i}^{*}(x)=\sum_{j}^{k}B_{j}(x)\beta^{*}_{j}$ (16)

with smoothing matrix $\mathbf{S}$ as in (2). From a deep learning perspective
the Fisher information is especially interesting in context of optimization,
as it can capture the local geometry of the loss landscape (Karakida et al.,
2019) and thus helps to improve the gradient steps. Since the true Fisher
information is not always easy to derive, but the gradients are calculated
anyways during the optimization, the so-called empirical Fisher information is
used instead to estimate $\mathbf{F(\mathbf{\hat{\beta}}})$ as an averaged
outer product of the (stochastic) gradients (Martens, 2014; Soen & Sun, 2021)

$\mathbf{\tilde{F}}(\theta)=\frac{1}{N}\sum_{i}^{N}\left[\Delta_{\theta}l(\theta,x,y)\Delta_{\theta}l(\theta,x,y)^{T}\right]$
(17)

While the empirical Fisher is easy to derive, in practice it is barely used
directly in deep learning since the storage capacity grows quadratic with the
number of parameters. This makes it infeasible for ordinary neural networks
and instead diagonal- or factorized approximations are used to cope with the
complexity, but those approaches are still under active research (Martens &
Grosse, 2015; Frantar et al., 2021). Here, the extreme parameter sparsity of
the SplineNAMs pays off since $\tilde{\mathbf{F}}$ is actually feasible for
reasonable settings.111(Kunstner et al., 2019) provides a critical discussion
of the use of the empirical Fisher and even though he admits that it can
improve the optimization he also points out, that it is only a consistent
estimate for the Fisher information close to the minimum, with a sufficient
amount of data and for a correctly specified model. Especially the first two
conditions justify the use of the empirical Fisher to construct confidence
intervals for the SplineNAMs as it is only evaluated after the fitting
procedure and respectively few parameters are involved.

##  4 Benchmarking the Accuracy of SNAMs

Table 1: Results of the benchmark study: We orientate on the benchmark study conducted in Agarwal et al. (2021b) and compare SNAMs to state-of-the-art neural as well as non-neural approaches. We report the means and standard deviations of 5-fold cross validation on multiple datasets. We report the results on widely used and publicly available datasets. A detailed description on the datasets as well as the preprocessing steps can be found in the appendix A.1.  Model | No. Params |  CA Housing (RMSE)$\downarrow$ |  Insurance (RMSE) $\downarrow$ |  Credit (AUC) $\uparrow$ |  FICO (AUC) $\uparrow$
---|---|---|---|---|---
MLP | $>90,000$ | 0.433 $\pm$ (0.009) | 0.425 $\pm$ (0.031) | 0.942 $\pm$ (0.012) | 0.795 $\pm$ (0.006)
XGBoost | - | 0.403 $\pm$ (0.008) | 0.427 $\pm$ (0.025) | 0.961 $\pm$ (0.022) | 0.714 $\pm$ (0.008)
EBM | - | 0.441 $\pm$ (0.011) | 0.384 $\pm$ (0.025) | 0.954 $\pm$ (0.018) | 0.727 $\pm$ (0.011)
NAM | $>500,000$ | 0.483 $\pm$ (0.011) | 0.497 $\pm$ (0.028) | 0.904 $\pm$ (0.012) | 0.799 $\pm$ (0.008)
NAM(ExU) | $>500,000$ | 0.510 $\pm$ (0.011) | 0.499 $\pm$ (0.029) | 0.913 $\pm$ (0.013) | 0.790 $\pm$ (0.008)
SNAM | $<2.400$ | 0.396 $\pm$ (0.009) | 0.487 $\pm$ (0.039) | 0.948 $\pm$ (0.011) | 0.782 $\pm$ (0.011)
FAST SNAM | $<9,600$ | 0.390 $\pm$ (0.016) | 0.482 $\pm$ (0.049) | 0.968 $\pm$ (0.005) | 0.788 $\pm$ (0.008)

In this section we asses the predictive performance of SNAMs in comparison to
several state-of-the-art models including neural as well as non-neural
approaches. The following four datasets were chosen for this purpose to
reflect real-world uses cases where model interpretability is a key factor:

  * •

California Housing: A widely used regression dataset for predicting home
prices in California from the U.S. census in 1990.

  * •

Insurance: A regression dataset for predicting billed medical expenses
containing 1338 observations.

  * •

Credit: A fraud detection dataset with a binary response and heavily
unbalanced class distributions.

  * •

FICO: A classification dataset for modelling credit score predictions.

A detailed description of the datasets and the contained features can be found
in the Appendix. The data pre-processing was adapted from the benchmark study
of Gorishniy (Gorishniy et al., 2021), who evaluated various DNNs on tabular
data. In particular a quantile scaling was applied to all features and the
target variable was centered to zero mean and standardized to unit variance.
For The regression problems (California Housing, Insurance) the averaged RMSE
over a 5-fold cross validation is reported in table 1 and for the
classification problems (FICO, Credit) the AUC is taken respectively.

We used the following benchmarks to asses the performance of the presented
method:

  * •

Multilayer Perceptron (MLP): Unrestricted fully connected deep neural network.
The loss function is either defined as the mean squared error (regression) or
binary cross entropy (logistic regression).

  * •

Gradient Boosted Trees (XGBoost): Decision tree based gradient boosting. We
use the implementation provided by (Chen & Guestrin, 2016).

  * •

Explainable Boosting Machines (EBMs): State-of-the-art Generalized Additive
Models leveraging shallow boosted trees (Lou et al., 2013).

  * •

Neural Additive Models (NAMs): Linear combination of deep neural networks that
allow visual interpretability. Each feature is modeled by a different
subnetwork. We implement two NAM versions, one leveraging exp-centered (ExU)
hidden units (Agarwal et al., 2021b) and one using ReLU activation functions
throughout the hidden layers.

The results show that the SNAMs are all throughout able to reach at least the
same performance range as traditional NAMs and even outperform them on
California Housing, FICO and Insurance. The insurance dataset has respectively
few observations(N=X). In this case the improvement of the SNAMs might be
attributed to their parameter sparsity and consequently a enhanced small-
sample behaviour. However, in particular noteworthy is the substatial
performance gain that the SNAMs achieved on the California Housing dataset.
This is likely due to the SNAM’s ability to jointly process the geographical
coordinates as a spatial effect and highlight the importance of integrating
structured data appropriately into DNNs.

###  4.1 Interpretability

While achieving remarkable results on the benchmark datasets, SNAMs remain to
be fully explainable glass-box-models and not only allow to investigate the
exact decision making process, but yield the identifiable effect of each
feature as shown in figure 6. The spatial features result in a 2D geographical
heatmap, see figure 5. That indicates which areas have a respectively positive
or negative effect on the house price. Usually, those estimations could be
interpreted directly, but due to the quantile scaling of the features, this
interpretation has slight restrictions in terms of the scale of the x-axis.
Nevertheless, valuable conclusions can be derived such as that the median
income, as well as the average rooms have an increasing positive effect,
whereas for the average number of household members it is in reverse. The
house age indicates that in certain timespans more valuable houses were built
and the average number of bedrooms seems to not have a big effect, except for
large numbers, which decreases the house value.

Figure 6: The results of the SNAM on the CA Housing dataset showing the
estimated nonlinear effects of the features with Bayesian confidence
intervals.

##  5 Related Work

Achieving interpretability in deep neural networks is inherently difficult.
While several approaches have successfully created interpretability at the
sample level (Shapley, 1953; Sundararajan & Najmi, 2020; Baehrens et al.,
2010; Ribeiro et al., 2016), the interpretability of classical statistical
generalized additive models have yet to be matched in a neural framework. The
idea of transferring the advantages of generalized additive models to neural
approaches was already introduced by Potts (1999). However, while being
remarkably parameter sparse, the framework did not use backpropagation and was
therefore ultimately unsuccessful in implementing the GAM structure into a
neural framework. de Waal & du Toit (2007) translate GAMs into a neural
framework. More recently, Agarwal et al. (2021b) introduced NAMs, a more
flexible approach than the GANNs introduced by de Waal & du Toit (2007) and
exploiting the advances made in the field of deep learning. Similar to SNAMs,
NAMs learn the non-linear additive functions with seperate networks for each
feature. Using large numbers of hidden units for each feature and exploiting
flexible exponentially centered activation functions, they successfully
integrate the predictive power of neural networks into an additive framework.
However, while the linear combination of neural sub-networks allows for good
visual interpretation of the results, any interpretability beyond plotting the
model predictions on a feature level is lost in their black-box sub-networks.
Additionally, the predictive power comes at the cost of an enormous number of
parameters and a huge hyperparameter space for optimization. Building upon
NAMs, several extensions have been proposed. Chang et al. (2021) introduced
NODE-GAM, extending the works of Lou et al. (2012) a differentiable oblivious
decision trees based model designed for high-risk domains. Yang et al. (2021)
expand NAMs to adjust for pairwise interaction effects. Rügamer et al. (2020)
extends statistical structural additive models by neural networks constructing
an orthogonalization cell. Thielmann et al. (2023) extend the framework to
account for arbitrarily many distributional parameters.

##  6 Conclusion

In this work we present Structural Neural Additive Models as an intelligible
extension of the recently proposed Neural Additive Models. We demonstrate the
integration of basis expansion layers into neural networks as a way to
flexibly model nonlinear functions with extremely few parameters. Furthermore,
we derive computationally efficient neural splines based on the Silverman
kernel. In a benchmark study we show that this approach is able to reach the
same predictive performance as other state-of-the-art models, while remaining
inherently more interpretable. In conclusion, Structural Neural Additive
Models (SNAMs) offer a promising architecture for interpretable neural
networks, combining classical statistical methods with the predictive power of
neural architectures.

## References

  * Adadi & Berrada (2018) Adadi, A. and Berrada, M.  Peeking Inside the Black-Box: A Survey on Explainable Artificial Intelligence (XAI).  _IEEE Access_ , 6:52138–52160, 2018.  ISSN 2169-3536.  doi: 10.1109/ACCESS.2018.2870052.  URL https://ieeexplore.ieee.org/document/8466590/.
  * Agarwal et al. (2021a) Agarwal, R., Melnick, L., Frosst, N., Zhang, X., Lengerich, B., Caruana, R., and Hinton, G. E.  Neural additive models: Interpretable machine learning with neural nets.  _Advances in Neural Information Processing Systems_ , 34:4699–4711, 2021a.
  * Agarwal et al. (2021b) Agarwal, R., Melnick, L., Frosst, N., Zhang, X., Lengerich, B., Caruana, R., and Hinton, G. E.  Neural additive models: Interpretable machine learning with neural nets.  _Advances in Neural Information Processing Systems_ , 34, 2021b.
  * Bach et al. (2015) Bach, S., Binder, A., Montavon, G., Klauschen, F., Müller, K.-R., and Samek, W.  On pixel-wise explanations for non-linear classifier decisions by layer-wise relevance propagation.  _PloS one_ , 10(7):e0130140, 2015.
  * Baehrens et al. (2010) Baehrens, D., Schroeter, T., Harmeling, S., Kawanabe, M., Hansen, K., and Müller, K.-R.  How to explain individual classification decisions.  _The Journal of Machine Learning Research_ , 11:1803–1831, 2010.
  * Caruana et al. (2015) Caruana, R., Lou, Y., Gehrke, J., Koch, P., Sturm, M., and Elhadad, N.  Intelligible models for healthcare: Predicting pneumonia risk and hospital 30-day readmission.  In _Proceedings of the 21th ACM SIGKDD international conference on knowledge discovery and data mining_ , pp. 1721–1730, 2015.
  * Chang et al. (2021) Chang, C.-H., Caruana, R., and Goldenberg, A.  Node-gam: Neural generalized additive model for interpretable deep learning.  _arXiv preprint arXiv:2106.01613_ , 2021.
  * Chen & Guestrin (2016) Chen, T. and Guestrin, C.  XGBoost: A scalable tree boosting system.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’16, pp. 785–794, New York, NY, USA, 2016. ACM.  ISBN 978-1-4503-4232-2.  doi: 10.1145/2939672.2939785.  URL http://doi.acm.org/10.1145/2939672.2939785.
  * Dal Pozzolo (2015) Dal Pozzolo, A.  Adaptive machine learning for credit card fraud detection.  2015\.
  * de Waal & du Toit (2007) de Waal, D. A. and du Toit, J. V.  Generalized additive models from a neural network perspective.  In _Seventh IEEE International Conference on Data Mining Workshops (ICDMW 2007)_ , pp. 265–270. IEEE, 2007.
  * Fahrmeir et al. (2013) Fahrmeir, L., Kneib, T., Lang, S., and Marx, B.  _Regression_.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2013.  ISBN 978-3-642-34332-2 978-3-642-34333-9.  doi: 10.1007/978-3-642-34333-9.
  * FICO (2018) FICO.  Fico explainable machine learning challenge.  2018\.
  * Frantar et al. (2021) Frantar, E., Kurtic, E., and Alistarh, D.  M-FAC: Efficient Matrix-Free Approximations of Second-Order Information.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Gálvez & Iglesias (2011) Gálvez, A. and Iglesias, A.  Efficient particle swarm optimization approach for data fitting with free knot b-splines.  _Computer-Aided Design_ , 43(12):1683–1692, 2011\.
  * Gorishniy et al. (2021) Gorishniy, Y., Rubachev, I., Khrulkov, V., and Babenko, A.  Revisiting deep learning models for tabular data.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Hastie & Tibshirani (1986) Hastie, T. and Tibshirani, R.  Generalized Additive Models.  _Statistical Science_ , 1(3), August 1986.  ISSN 0883-4237.  doi: 10.1214/ss/1177013604.
  * Hastie et al. (2009) Hastie, T., Tibshirani, R., Friedman, J. H., and Friedman, J. H.  _The elements of statistical learning: data mining, inference, and prediction_ , volume 2.  Springer, 2009.
  * Karakida et al. (2019) Karakida, R., Akaho, S., and Amari, S.-i.  Universal statistics of fisher information in deep neural networks: Mean field approach.  In Chaudhuri, K. and Sugiyama, M. (eds.), _Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics_ , volume 89 of _Proceedings of Machine Learning Research_ , pp. 1032–1041. PMLR, 16–18 Apr 2019.  URL https://proceedings.mlr.press/v89/karakida19a.html.
  * Kingma & Ba (2014) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kunstner et al. (2019) Kunstner, F., Hennig, P., and Balles, L.  Limitations of the empirical fisher approximation for natural gradient descent.  _Advances in neural information processing systems_ , 32, 2019.
  * Lantz (2019) Lantz, B.  _Machine learning with R: expert techniques for predictive modeling_.  Packt publishing ltd, 2019.
  * Lou et al. (2012) Lou, Y., Caruana, R., and Gehrke, J.  Intelligible models for classification and regression.  In _Proceedings of the 18th ACM SIGKDD international conference on Knowledge discovery and data mining - KDD ’12_ , pp. 150, Beijing, China, 2012. ACM Press.  ISBN 978-1-4503-1462-6.  doi: 10.1145/2339530.2339556.  URL http://dl.acm.org/citation.cfm?doid=2339530.2339556.
  * Lou et al. (2013) Lou, Y., Caruana, R., Gehrke, J., and Hooker, G.  Accurate intelligible models with pairwise interactions.  In _Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining_ , pp. 623–631, 2013.
  * Ly et al. (2017) Ly, A., Marsman, M., Verhagen, J., Grasman, R. P., and Wagenmakers, E.-J.  A tutorial on Fisher information.  _Journal of Mathematical Psychology_ , 80:40–55, 2017.  Publisher: Elsevier.
  * Martens (2014) Martens, J.  New insights and perspectives on the natural gradient method.  _arXiv preprint arXiv:1412.1193_ , 2014.
  * Martens & Grosse (2015) Martens, J. and Grosse, R.  Optimizing neural networks with kronecker-factored approximate curvature.  In _International conference on machine learning_ , pp. 2408–2417. PMLR, 2015.
  * Mohanty & Fahnestock (2021) Mohanty, S. D. and Fahnestock, E.  Adaptive spline fitting with particle swarm optimization.  _Computational Statistics_ , 36(1):155–191, 2021\.
  * Murdoch et al. (2019) Murdoch, W. J., Singh, C., Kumbier, K., Abbasi-Asl, R., and Yu, B.  Definitions, methods, and applications in interpretable machine learning.  _Proceedings of the National Academy of Sciences_ , 116(44):22071–22080, 2019.
  * Nelder & Wedderburn (1972) Nelder, J. A. and Wedderburn, R. W.  Generalized linear models.  _Journal of the Royal Statistical Society: Series A (General)_ , 135(3):370–384, 1972.  Publisher: Wiley Online Library.
  * Pace & Barry (1997) Pace, R. K. and Barry, R.  Sparse spatial autoregressions.  _Statistics & Probability Letters_, 33(3):291–297, 1997.  Publisher: Elsevier.
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., and Duchesnay, E.  Scikit-learn: Machine Learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011\.
  * Potts (1999) Potts, W. J.  Generalized additive neural networks.  In _Proceedings of the fifth ACM SIGKDD international conference on Knowledge discovery and data mining_ , pp. 194–200, 1999.
  * Ribeiro et al. (2016) Ribeiro, M. T., Singh, S., and Guestrin, C.  "Why Should I Trust You?": Explaining the Predictions of Any Classifier.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’16, pp. 1135–1144, New York, NY, USA, 2016. Association for Computing Machinery.  ISBN 978-1-4503-4232-2.  doi: 10.1145/2939672.2939778.  URL https://doi.org/10.1145/2939672.2939778.  event-place: San Francisco, California, USA.
  * Rügamer et al. (2020) Rügamer, D., Kolb, C., and Klein, N.  Semi-Structured Deep Distributional Regression: Combining Structured Additive Models and Deep Learning.  _arXiv preprint arXiv:2002.05777_ , 2020.
  * Sejnowski (2020) Sejnowski, T. J.  The unreasonable effectiveness of deep learning in artificial intelligence.  _Proceedings of the National Academy of Sciences_ , 117(48):30033–30038, 2020.  Publisher: National Acad Sciences.
  * Shapley (1953) Shapley, L.  Quota solutions op n-person games1.  _Edited by Emil Artin and Marston Morse_ , pp. 343, 1953.
  * Silverman (1984) Silverman, B. W.  Spline smoothing: the equivalent variable kernel method.  _The annals of Statistics_ , pp. 898–916, 1984.  Publisher: JSTOR.
  * Soen & Sun (2021) Soen, A. and Sun, K.  On the Variance of the Fisher Information for Deep Learning.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Sundararajan & Najmi (2020) Sundararajan, M. and Najmi, A.  The many shapley values for model explanation.  In _International conference on machine learning_ , pp. 9269–9278. PMLR, 2020.
  * Thielmann et al. (2023) Thielmann, A., Kruse, R.-M., Kneib, T., and Säfken, B.  Neural additive models for location scale and shape: A framework for interpretable neural regression beyond the mean.  _arXiv preprint arXiv:2301.11862_ , 2023.
  * Wood (2017) Wood, S. N.  _Generalized additive models: an introduction with R_.  Chapman & Hall/CRC texts in statistical science. CRC Press/Taylor & Francis Group, Boca Raton, second edition edition, 2017.  ISBN 978-1-4987-2833-1.
  * Yang et al. (2021) Yang, Z., Zhang, A., and Sudjianto, A.  Gami-net: An explainable neural network based on generalized additive models with structured interactions.  _Pattern Recognition_ , 120:108192, 2021.

##  Appendix A Benchmarking

###  A.1 Data sets & Preprocessing

##### California Housing

The California Housing (CA) dataset (Pace & Barry, 1997) contains the house
prices for areas in California from the U.S. census in 1990. The data was
obtained from sklearn (Pedregosa et al., 2011) and besides the logarithmic
median house price of the blockwise areas as the target variable it contains
the eight predictors from the following list. Of particular interest is, that
the CA dataset contains the longitude and latitude as spatial features.
Considering that those features are more meaningful when they are interpreted
mutually, they can serve as a good example for a usecase of tensor product
splines in SplineNAMs. Besides that, the remaining features are modelled
univariately with cubic spline units and for the sake of simplicity a
sophisticated hyperparameter tuning is omitted. Instead the number of bases
and feature-wise penalties are chosen heuristically by trying out a small
number of configurations on a single training set.

  * •

Targets: Median House Price

  * •

Features:

    1. 1.

MedInc: Median income of the area

    2. 2.

HouseAge: Median house age of the area

    3. 3.

AveRooms: Area-wise average number of rooms per household

    4. 4.

AveBedrms: Area-wise average number of bedrooms per household

    5. 5.

Population:Total number of people living in the area

    6. 6.

AveOccup: Area-wise average number of people within a household

    7. 7.

Latitude: Latitude of the center of the area

    8. 8.

Longitude: Longitude of the center of the area

##### Insurance

The Insurance dataset is another regression type dataset for predicting billed
medical expenses (Lantz, 2019). It is a relatively small dataset with 1338
observations. Besides the target variable, charges, it contains 6 variables.

  * •

Targets: Charges

  * •

Features:

    1. 1.

Age: The patient’s age

    2. 2.

Sex: Gender of the patient

    3. 3.

BMI: The body-mass index of the patient

    4. 4.

Children: The number of children the patient has

    5. 5.

Smoker: Whether the patient is a smoker or not

    6. 6.

Region: The region where the patient lives, categorized into 4 sections,
southeast, southwest, northwest and northeast

##### Credit

The credit card dataset (Dal Pozzolo, 2015) is a dataset used for fraud
detection and hence highly unbalanced. It is comprised of 108754 observations
and only 0.172% of the observations are fraudulent. It is also used as a
benchmark dataset by Agarwal et al. (2021b). The dataset is comprised of 30
numerical input variables in addition to the target variable which is a binary
identifier whether there was financial fraud or not. The meaning for most of
the variables is not revealed due to confidentiality reasons (Dal Pozzolo,
2015). All variables have already been transformed by principal component
analysis. The few variables which are revealed are listed below.

  * •

Targets: Class $\rightarrow$ Fraud or not

  * •

Features:

    1. 1.

Time: The seconds elapsed between each transaction and the first transaction
in the dataset

    2. 2.

Amount: The transaction amount

##### FICO

As a second dataset for a classification benchmark we use the FICO dataset
(FICO, 2018). We use the Risk Performance as the target variable and all other
variables as features. A detailed description of the features and their
meaning is available at FICO (2018). The dataset is comprised of 10459
observations.

##### Preprocessing

We implement the same preprocessing for all four dataset and closely follow
Gorishniy et al. (2021) in their preprocessing steps. Hence all numerical
variables are scaled between -1 and 1, all categorical features are one-hot
encoded and we implement quantile smoothing. We use 5-fold cross-validation
and report mean results as well as the standard deviations. For
reproducability we use the sklearn (Pedregosa et al., 2011) Kfold function
with a random state of 101 and shuffle equal to true for all datasets.

###  A.2 Hyperparameter Tuning

We implement light hyperparameter tuning over all models and closely follow
Agarwal et al. (2021b) reported model structures. Hence, we implement the
Multilayer Perceptrons with 10 layers and 100 neurons each. Each layer uses
the rectified linear unit (ReLU) activation functions and is followed by a
0.25 dropout layer. Similar to Agarwal et al. (2021b) we implement an early
stopping with a patience of 100 epochs and find no overfitting for all models.
Additionally, we reduce the learning rate by a factor of 0.995 once a plateau
is reached and start with a learning rate of 1e-04 using the ADAM optimizer
(Kingma & Ba, 2014). All models are trained with an initial number of 1000
epochs, however, over all neural models, all datasets and all folds every
model stopped training before reaching the maximum number of epochs. We
iterate over different batch sizes and find that, similar to Agarwal et al.
(2021b), 1024 is the best batch size for the California Housing dataset and
the FICO dataset. For the creditcard dataset we find a batchsize of 512 to
perform best and for the Insurance dataset, which is the smallest of the four
datasets, we find a smaller batchsize of 64 to perform best. We implement the
Neural Additive Models (Agarwal et al., 2021b) identically to the
implementations described in Agarwal et al. (2021b). Additionally, we
experiment with the described ExU activation functions but find that simple
ReLU activation functions outperform the ExU functions for 3 out of 4
datasets. For XGBoost we use the following hyperparameters: number of
estimators=1000, maximum depth=25, eta=0.1, subsample=1.0, colsample by
tree=0.95. Interestingly, XGBoost performs extremely well on 3 out of 4
datasets but is clearly outperformed by all neural approaches for the FICO
dataset. For EBM we use the following hyperparameters: max leaves=10, learning
rate=0.005, early stopping rounds=100. EBM perform very similar to XGBoost but
outperform all models for the small Insurance dataset.