Scalable Neural Network Kernels
We introduce the concept of scalable neural network kernels (SNNKs), the
replacements of regular feedforward layers (FFLs), capable of approximating the
latter, but with favorable computational properties. SNNKs effectively
disentangle the inputs from the parameters of the neural network in the FFL,
only to connect them in the final computation via the dot-product kernel. They
are also strictly more expressive, as allowing to model complicated
relationships beyond the functions of the dot-products of parameter-input
vectors. We also introduce the neural network bundling process that applies
SNNKs to compactify deep neural network architectures, resulting in additional
compression gains. In its extreme version, it leads to the fully bundled
network whose optimal parameters can be expressed via explicit formulae for
several loss functions (e.g. mean squared error), opening a possibility to
bypass backpropagation. As a by-product of our analysis, we introduce the
mechanism of the universal random features (or URFs), applied to instantiate
several SNNK variants, and interesting on its own in the context of scalable
kernel methods. We provide rigorous theoretical analysis of all these concepts
as well as an extensive empirical evaluation, ranging from point-wise kernel
estimation to Transformers' fine-tuning with novel adapter layers inspired by
SNNKs. Our mechanism provides up to 5x reduction in the number of trainable
parameters, while maintaining competitive accuracy.
Arijit Sehanobish, Krzysztof Choromanski, Yunfan Zhao, Avinava Dubey, Valerii Likhosherstov
# Scalable Neural Network Kernels

###### Abstract

We introduce the concept of _scalable neural network kernels_ (SNNKs), the replacements of regular _feedforward layers_ (FFLs), capable of approximating the latter, but with favorable computational properties. SNNKs effectively disentangle the inputs from the parameters of the neural network in the FFL, only to connect them in the final computation via the dot-product kernel. They are also strictly more expressive, as allowing to model complicated relationships beyond the functions of the dot-products of parameter-input vectors. We also introduce the _neural network bundling process_ that applies SNNKs to compactify deep neural network architectures, resulting in additional compression gains. In its extreme version, it leads to the fully bundled network whose optimal parameters can be expressed via explicit formulae for several loss functions (e.g. mean squared error), opening a possibility to bypass backpropagation. As a by-product of our analysis, we introduce the mechanism of the _universal random features_ (or URFs), applied to instantiate several SNNK variants, and interesting on its own in the context of scalable kernel methods. We provide rigorous theoretical analysis of all these concepts as well as an extensive empirical evaluation, ranging from point-wise kernel estimation to Transformers' fine-tuning with novel adapter layers inspired by SNNKs. Our mechanism provides up to 5x reduction in the number of trainable parameters, while maintaining competitive accuracy.

## 1 Introduction

Consider a kernel function: \(\mathrm{K}:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}\), taking as input two feature vectors encoding latent embeddings of their corresponding objects and returning their similarity. Kernel methods are among the most theoretically principled approaches to statistical machine learning (ML) and have proven effective in numerous real-world problems (Scholkopf & Smola, 2002; Kontorovich et al., 2008). Despite their theoretical guarantees and applicability in a rich spectrum of ML settings, the main drawback of these techniques is a high computational complexity, at least quadratic in the size \(N\) of the training dataset. For example, the kernel regression has time complexity \(\mathcal{O}(N^{3})\).

To address this issue, Rahimi & Recht (2007) proposed to construct a random feature (RF) map \(\Phi:\mathbb{R}^{d}\to\mathbb{R}^{m}\) that transforms an input point \(\mathbf{z}\) to a finite-dimensional feature vector \(\Phi(\mathbf{z})\in\mathbb{R}^{m}\) such that: \(\mathrm{K}(\mathbf{x},\mathbf{y})=\mathbb{E}[\Phi(\mathbf{x})^{\top}\Phi( \mathbf{y})]\) (effectively approximately linearizing kernel function). Approximating general kernels \(\mathrm{K}(\mathbf{x},\mathbf{y})\) via linear (dot-product) kernels \(\mathrm{K}(\mathbf{x},\mathbf{y})\approx\widehat{\mathbf{x}}^{\top}\widehat{ \mathbf{y}}\) for \(\widehat{\mathbf{z}}=\Phi(\mathbf{z})\) drastically changes computational complexity landscape, which is now dominated by the number \(m\) of random features, thus providing computational gains if \(m\ll N\). Since their seminal work, there has been a variety of works proposing random features for a broad range of kernels like Gaussian, Matern (Choromanski et al., 2018) and polynomial (Kar & Karnick, 2012; Wacker et al., 2022; 1).

In the meantime, with the advances in: the optimization algorithms for deep ML architectures and the accelerators' hardware, neural network (NN) models (Goodfellow et al., 2016; Schmidhuber, 2014; LeCun et al., 2015) have become predominant in machine learning. The _feedforward layer
(FFL) is the core computational module of NNs and is of the following form:

\[\mathbf{x}\to f(\mathbf{W}\mathbf{x}+\mathbf{b}) \tag{1}\]

for \(\mathbf{x}\in\mathbb{R}^{d},\mathbf{W}\in\mathbb{R}^{l\times d},\mathbf{b}\in \mathbb{R}^{l}\) (_bias_) and an _activation function_\(f:\mathbb{R}\to\mathbb{R}\) (applied element-wise). The expressiveness of deep NNs, far surpassing standard kernel methods, comes from stacking together several FFLs, each encoding **non-linear** mapping with **learnable \(\mathbf{W},\mathbf{b}\)**.

In this work, we draw a deep connection between scalable kernel methods and neural networks. We reinterpret the FFL as outputting the expected vector of dot-products of: **(1)** the latent embeddings of the input \(\mathbf{x}\) and **(2)** the parameters: \(\mathbf{W}\), b of the FFL, effectively disentangling input from model's parameters in the computational graph, only to connect them in the final computation via the dot-product kernel. To be more specific, we think about the FFL as the following transformation:

\[\begin{cases}\overline{\mathrm{K}}_{f}(\mathbf{x},(\mathbf{W},\mathbf{b})) \stackrel{{\mathrm{def}}}{{=}}\left(\mathrm{K}_{f}(\mathbf{x},( \mathbf{w}^{0},b_{0})),...,\mathrm{K}_{f}(\mathbf{x},(\mathbf{w}^{l-1},b_{l- 1}))\right)^{\top},\\ \mathrm{K}_{f}(\mathbf{x},(\mathbf{w},b))\stackrel{{\mathrm{def} }}{{=}}\mathbb{E}[\Phi_{f}(\mathbf{x})^{\top}\Psi_{f}(\mathbf{w},b)],\end{cases} \tag{2}\]

where mappings: \(\Phi_{f}:\mathbb{R}^{d}\to\mathbb{R}^{m},\Psi_{f}:\mathbb{R}^{d}\times\mathbb{ R}\to\mathbb{R}^{m}\) satisfy: \(f(\mathbf{w}^{\top}\mathbf{x}+b)=\mathbb{E}[\Phi_{f}(\mathbf{x})^{\top}\Psi_{f }(\mathbf{w},b)]\) and \(\mathbf{w}^{0},...\mathbf{w}^{l-1}\) are the transposed rows of \(\mathbf{W}\). Then, in the instantiation of the layer the expectations are dropped out. Rewriting an FFL in terms of two towers: one corresponding to the input and one to its learnable parameters has several advantages:

1. **network compression:** in the above formulation, instead of transforming layer parameters with \(\Psi_{f}\), one can directly learn vectors \(\Psi_{f}(\mathbf{w}^{i},b_{i})\) for \(i=0,...,l-1\). Then the number of trainable parameters becomes \(O(lm)\) rather than \(O(ld)\) and for \(m\ll d\) the layer effectively has a reduced number of parameters.
2. **computational savings:** if RFs can be constructed in time \(o(dl)\) per point and \(m\ll d\), the overall time complexity \(o(dl)\) of the FFL (given pre-computed embeddings \(\Psi_{f}(\mathbf{w}^{i},b_{i})\)) is **sub-quadratic** in layers' dimensionalities,
3. **deep NN bundling process:** a two-tower representation can be used iteratively to compactify multiple FFLs of NNs, the process we refer to as _neural network bundling_ (Sec. 3.3); this also leads to the computational gains.
4. **deep NNs as scalable kernels:** the extreme version of the bundling procedure, involving all the layers, provides a two-tower factorization of the entire deep NN with several potential practical and theoretical implications (Sec. 3.3). In particular, it leads to an explicit formula for the optimal parameters of the fully-bundled network under several loss objectives (e.g. mean squared loss), opening a possibility to bypass backpropagation.

Figure 1: Pictorial representation of different NN layers discussed in the paper. Pink arrays represent NN weight matrices and grey ones, Gaussian projections matrices applied in SNNKs. Nonlinear transformations applied in mappings \(\Phi\) and \(\Psi\) are symbolically represented as functions \(g\) and \(h\) respectively. **Upper left:** Regular FFL with activation \(f\). **Upper right:** SNNK applied to a single FFL. **Bottom:** Bundling process using SNNKs and applied two a deep neural network module.
In order to find mappings: \(\Phi_{f}\), \(\Psi_{f}\) from Eq. 2, we develop a new bounded random feature map mechanism, called _universal random features_ (or URFs) that leads to the unbiased estimation of \(f(\mathbf{w}^{\top}\mathbf{x}+b)\) as long as \(f\) has a well-defined Fourier Transform (FT), either in the classical Riemannian or distributional sense. To derive URFs, we combine Fourier analysis techniques with recent methods for softmax-kernel estimation from Likhosherstov et al. (2022).

**Note:** We do not put any additional assumptions regarding \(f\), in particular \(f\) is not required to be differentiable. Furthermore, function \(\mathrm{K}_{f}\)**does not need** to be positive semi-definite. This is critical for applications in neural networks, where the activation function \(f\) usually does not correspond to a positive semi-definite kernel.

To summarize, our main contributions in this paper are as follows:

* We introduce the _scalable neural network kernel_ module (SNNK) as a replacement of a traditional FFL (Sec. 3), providing the disentanglement of the network's input and its parameter-set before final dot-product computation, as given in Eq. 2 (see also: Fig. 1).
* We accompany SNNKs with our universal random features mechanism (URFs) to efficiently: (1) construct mappings \(\Phi_{f}\) and \(\Psi_{f}\) from Eq. 2 and consequently: (2) implement SNNKs (Sec. 3.1). We provide explicit formulae for URFs for trigonometric maps. Those produce SNNK-based replacements of the SIREN networks from Sitzmann et al. (2020).
* We propose new NN-layers corresponding to the specific SNNK instantiation, called \(\mathrm{ReLU}\)-\(\mathrm{SNNK}\) (Sec. 3.2), that we found particularly effective in downstream applications (see: Sec. 4.3.2). We show that they are related to the class of the _arc-cosine kernels_Cho & Saul (2011). We also demonstrate using them that SNNKs are **strictly more expressive** than regular FFLs, as allowing to compute the functions of the inputs and parameters that cannot be defined as point-wise transformed vectors of their dot-products.
* We introduce the neural network compactification process, that we refer to as _neural network bundling_, leveraging SNNKs (see: Sec. 3.3 and Fig. 1).
* We provide an exhaustive empirical evaluation of SNNKs, from point-wise kernel estimation to the adapter-based Transformers' fine-tuning, providing about 5x reduction of the number of trainable parameters (Sec. 4).

## 2 Related Work

The literature on random features is vast, yet most of the works focus on approximating positive definite kernels. The results on dimensionality reduction and the so-called _Johnson-Lindenstrauss Transform_ (or JLT) (Dasgupta & Gupta, 2003; Dasgupta et al., 2010; Ailon & Liberty, 2013) for the dot-product kernel marked the birth of the subject as as an archetype mechanism that Rahimi & Recht (2007) extended from linear to non-linear shift-invariant kernels. A substantial effort was made to further improve the accuracy of RF-methods by entangling projections used to construct RFs (Choromanski et al., 2017; Yu et al., 2016; Choromanski et al., 2018; Rowland et al., 2018).

For certain classes of functions \(f\), RF-mechanisms leading to the linearization of \(\mathrm{K}_{f}\) have been already developed. In addition to the rich recent literature on the approximation techniques for the softmax-kernel \(\mathrm{K}_{\mathrm{exp}}(\mathbf{x},\mathbf{y})=\exp(\mathbf{x}^{\top} \mathbf{y})\)(Likhosherstov et al., 2022; 2023; Choromanski et al., 2021), algorithms for analytic \(f\) with positive coefficients of their Taylor series expansion were given (Kar & Karnick, 2012). Other RF-methods assume that kernel inputs are taken from the unit-sphere (Scetbon & Harchaoui, 2021; Han et al., 2022). Both assumptions are unrealistic for the neural network applications as far as inputs \(\mathbf{x}\) are concerned (interestingly, the latter one would be however more justifiable for the parameter-tower as long as bounded-norm weight matrices are considered, e.g. _orthogonal neural networks_(Helfrich et al., 2018). We would like to emphasize that our two-tower mechanism, effectively leading to the linearization of the FFLs from Eq. 2, can in principle work with various RF-algorithms, and not only our proposed URFs.

The kernels applied in connection to neural networks have been widely studied (Bengio & Lecun, 2007). Such kernels are generally constructed using dot-products of outputs of the shallow neural networks with various non-linearities like ReLU (Cho & Saul, 2009; Bresler & Nagaraj, 2020) and tanh (Williams, 1996) or the gradients of the network like the NTK kernel (Jacot et al., 2020). Most
of the work in linearizing NNs via kernels have been done in the case of a \(2\)-layer network where

\[J(\mathbf{x};\mathbf{\theta})=\sum_{i=1}^{N}a_{i}f(\mathbf{x}^{\top}\mathbf{w}^{i}), \ \mathbf{\theta}=(a_{1},...,a_{N};\mathbf{w}^{1},...,\mathbf{w}^{N})\in\mathbb{R}^{N(d+1)} \tag{3}\]

It is assumed that \(\mathbf{w}^{i}\) and \(f\) (non-linear activation) are fixed and scalars and \(a_{i}\) are trainable. Under various assumptions, one can write a compact linearized form of this neural network (Cho & Saul, 2009; 2011; Ghorbani et al., 2020). Moreover, in the above setting, \(J(\mathbf{x};\mathbf{\theta})\) corresponds to the first-order Taylor expansion of \(J\) with respect to the top-layer weights \(a_{i}\) which was first explored by (Neal, 1996). Even though our setting is fundamentally different, as our goal is to linearize single layers to disentangle the weights and the inputs, we build on the above intuition to create our SNNK-layers (see also: discussion in Appendix A). Note that NTK-based analysis, as leveraging Taylor-based linearization of the NN, is valid only for the mature stage of training/finetuning when weights do not change much and thus such a linearization is accurate (Malladi et al., 2022). SNNKs do not rely on this assumption. Furthermore, SNNKs can be used also in the context of non-positive definite (Ong et al., 2004) and asymmetric (He et al., 2023) kernels since mappings \(\Phi\) and \(\Psi\) in principle are different (on expectation they can produce both symmetric and asymmetric functions).

Arc-cosine kernels were studied in the context of deep NNs before (Cho & Saul, 2009). However, in (Cho & Saul, 2009), the weights are still entangled with the FFL-input, as the initial latent representations of the inputs (for random parameters) are interpreted as RFs for the arc-cosine kernel.

## 3 Scalable Neural Network Kernels (SNNKs)

The scalable neural network kernel (SNNK) computational module is defined as follows:

\[\begin{cases}\overline{\mathrm{SNNK}}_{f}(\mathbf{x},(\mathbf{W},\mathbf{b}) )\stackrel{{\mathrm{def}}}{{=}}\left(\mathrm{SNNK}_{f}(\mathbf{x },(\mathbf{w}^{0},b_{0})),...,\mathrm{SNNK}_{f}(\mathbf{x},(\mathbf{w}^{l-1}, b_{l-1}))\right)^{\top},\\ \mathrm{SNNK}_{f}(\mathbf{x},(\mathbf{w},b))\stackrel{{\mathrm{ def}}}{{=}}\Phi_{f}(\mathbf{x})^{\top}\Psi_{f}(\mathbf{w},b),\end{cases} \tag{4}\]

for, \(\mathbf{x}\in\mathbb{R}^{d}\), some mappings: \(\Phi_{f}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{m}\), \(\Psi_{f}:\mathbb{R}^{d}\times\mathbb{R}\rightarrow\mathbb{R}^{m}\) and transposed rows of \(\mathbf{W}\in\mathbb{R}^{l\times d}\): \(\mathbf{w}^{0},...\mathbf{w}^{l-1}\). As we show in Sec. 3.1, functions \(\Phi_{f},\Psi_{f}\) can be constructed in such a way that the SSNK module approximates a particular FFL, i.e.: \(\overline{\mathrm{SNNK}}_{f}(\mathbf{x},(\mathbf{W},\mathbf{b}))\approx f( \mathbf{W}\mathbf{x}+\mathbf{b})\), but mechanisms that do not imitate known FFLs are also of interest (see: Sec. 3.2).

**Time complexity:** If we denote by \(t_{m}(d)\) time complexity for constructing an embedding: \(\Phi_{f}(\mathbf{x})\), then time complexity for constructing \(\overline{\mathrm{SNNK}}_{f}(\mathbf{x},(\mathbf{W},\mathbf{b}))\) (given the pre-computed \(\Psi_{f}(\mathbf{w}^{i},b_{i})\) for \(i=0,...,l-1\)) is: \(T_{m,l}(d)=ml+t_{m}(d)\). In Sec. 3.1 we show an algorithms for constructing URFs in time \(t_{m}(d)=O(md)\) and thus computational gains are provided as compared to the regular FFL (with time complexity \(O(ld)\)) as long as \(m\ll\min(l,d)\).

**FFL compression:** As already mentioned in Sec. 1, the key observation is that in the setting, where the layer is learned (and thus \(\mathbf{w}^{0},...,\mathbf{w}^{l-1}\) are learnable), mapping \(\Psi_{f}\)**does not even need to be applied**, since vectors \(\mathbf{\omega}^{\mathrm{def}}\equiv\Psi_{f}(\mathbf{w}^{j},b)\) for \(j=0,...,l-1\) can be interpreted as unstructured learnable vectors. Thus the number of trainable parameters of the SNNK layer is \(O(ml)\), instead of \(O(dl)\) and consequently, the FFL is effectively compressed if \(m\ll d\).

### Universal Random Features (URFs)

In this section, we show how to construct embeddings \(\Phi_{f}(\mathbf{x})\) and \(\Psi_{f}(\mathbf{w},b)\). We denote by \(\mathrm{FT}_{f}\) the _Fourier Transform_ of \(f\), where \(i\in\mathbb{C}\) satisfies: \(i^{2}=-1\):

\[\mathrm{FT}_{f}(\xi)=\int_{\mathbb{R}}f(z)\exp(-2\pi i\xi z)dz \tag{5}\]

If the integral does not exist in the classical Riemannian sense, we use its distributional interpretation. We rewrite \(\mathrm{FT}_{f}\) as: \(\mathrm{FT}_{f}=\mathrm{FT}_{f}^{\mathrm{re},+}-\mathrm{FT}_{f}^{\mathrm{re},- }+i\mathrm{FT}_{f}^{\mathrm{im},+}-i\mathrm{FT}_{f}^{\mathrm{im},-}\), where: \(\mathrm{FT}_{f}^{\mathrm{re},+},\mathrm{FT}_{f}^{\mathrm{re},-},\mathrm{FT}_{ f}^{\mathrm{im},+},\mathrm{FT}_{f}^{\mathrm{im},-}:\mathbb{R}\rightarrow\mathbb{R}_{\geq 0}\). Without loss of generality, we will assume that all four functions are not identically zero.

Let us denote by \(\overline{\mathcal{P}}_{0},\overline{\mathcal{P}}_{1},\overline{\mathcal{P}}_{2 },\overline{\mathcal{P}}_{3}\) some probabilistic distribution on \(\mathbb{R}\) (e.g. Gaussian) and by \(\overline{p}_{0},\overline{p}_{1},\overline{p}_{2},\overline{p}_{3}:\mathbb{R} \rightarrow\mathbb{R}_{\geq 0}\) their corresponding density functions. Furthermore, denote by
\(\mathcal{P}_{0},\mathcal{P}_{1},\mathcal{P}_{2},\mathcal{P}_{3}\) probabilistic distributions of densities: \(p_{0},p_{1},p_{2},p_{3}:\mathbb{R}\rightarrow\mathbb{R}_{\geq 0}\) proportional to: \(\mathrm{FT}_{f}^{\mathrm{re,+}},\mathrm{FT}_{f}^{\mathrm{re,-}},\mathrm{FT}_{f}^ {\mathrm{im,+}},\mathrm{FT}_{f}^{\mathrm{im,-}}\) respectively. We can then write:

\[\begin{split} f(z)=\int_{\mathbb{R}}\mathrm{FT}_{f}(\xi)\exp(2 \pi i\xi z)d\xi&=\sum_{j=0}^{3}c_{j}\int_{\mathbb{R}}\frac{p_{j}( \xi)}{\bar{p}_{j}(\xi)}\exp(2\pi i\xi z)\overline{p}_{j}(\xi)d\xi\\ &=\sum_{j=0}^{3}c_{j}\mathbb{E}_{\xi\sim\overline{p}_{j}}\left[ \frac{p_{j}(\xi)}{\bar{p}_{j}(\xi)}\exp(2\pi i\xi z)\right],\end{split} \tag{6}\]

where: \(c_{0}=\int_{\mathbb{R}}\mathrm{FT}_{f}^{\mathrm{re,+}}(\tau)d\tau,c_{1}=-\int_ {\mathbb{R}}\mathrm{FT}_{f}^{\mathrm{re,-}}(\tau)d\tau,c_{2}=i\int_{\mathbb{R }}\mathrm{FT}_{f}^{\mathrm{im,+}}(\tau)d\tau\), and furthermore \(c_{3}=-i\int_{\mathbb{R}}\mathrm{FT}_{f}^{\mathrm{im,-}}(\tau)d\tau\). For \(\mathbf{x},\mathbf{w}\in\mathbb{R}^{d},b\in\mathbb{R}\), let us denote:

\[\widehat{f}_{j}(\mathbf{x},\mathbf{w},b)=c_{j}\mathbb{E}_{\xi\sim\overline{p} _{j}}\left[\frac{p_{j}(\xi)}{\bar{p}_{j}(\xi)}\exp\left(2\pi i\xi(\mathbf{x}^{ \top}\mathbf{w}+b)\right)\right]=c_{j}\mathbb{E}_{\xi\sim\overline{\mathcal{P} _{j}}}[S_{j}(\xi,b)\exp(\widehat{\mathbf{x}}^{\top}(\xi)\widehat{\mathbf{w}}( \xi))] \tag{7}\]

for \(S_{j}(\xi,b)=\frac{p_{j}(\xi)}{\bar{p}_{j}(\xi)}\exp(2\pi i\xi b)\), \(\widehat{\mathbf{x}}(\xi)=\rho(\xi)\mathbf{x}\), \(\widehat{\mathbf{w}}(\xi)=\eta(\xi)\mathbf{w}\), where \(\rho(\xi),\eta(\xi)\in\mathbb{C}\) satisfy: \(\rho(\xi)\eta(\xi)=2\pi i\xi\). Inside the expectation in Eq. 7, we recognize the softmax-kernel value \(\mathrm{K}_{\exp}(\widehat{\mathbf{x}}(\xi),\widehat{\mathbf{w}}(\xi))=\exp( \widehat{\mathbf{x}}^{\top}(\xi)\widehat{\mathbf{w}}(\xi))\). We thus disentangle \(\widehat{\mathbf{x}}(\xi)\) from \(\widehat{\mathbf{w}}\) there, by applying softmax-kernel linearization mechanism from Likhosherstov et al. (2022): \(\exp(\widehat{\mathbf{x}}^{\top}(\xi)\widehat{\mathbf{w}}(\xi))=\mathbb{E}_ {\mathbf{g}\sim\mathcal{N}(0,\mathbf{I}_{d})}[\Lambda_{\mathbf{g}}(\widehat{ \mathbf{x}})\Lambda_{\mathbf{g}}(\widehat{\mathbf{w}})]\), where \(\Lambda_{\mathbf{g}}:\mathbb{R}^{d}\rightarrow\mathbb{R}\) is defined as follows for \(A\leq 0\):

\[\Lambda_{\mathbf{g}}(\mathbf{z})=(1-4A)^{\frac{d}{4}}\exp(A\|\mathbf{g}\|_{2}^ {2}+\sqrt{1-4A}\mathbf{g}^{\top}\mathbf{z}-\frac{\|\mathbf{z}\|_{2}^{2}}{2}) \tag{8}\]

Thus \(\widehat{f}_{j}(\mathbf{x},\mathbf{w},b)=\mathbb{E}_{(\xi,\mathbf{g})\sim \overline{\mathcal{P}}_{j}\otimes\mathcal{N}(0,\mathbf{I}_{d})}[\Gamma_{ \mathbf{g},\xi}^{1}(\mathbf{x})\Gamma_{\mathbf{g},\xi}^{2}(\mathbf{w},b)]\) for \(\Gamma_{\mathbf{g},\xi}^{1}(\mathbf{x}),\Gamma_{\mathbf{g},\xi}^{2}(\mathbf{ w},b)\) given as:

\[\Gamma_{\mathbf{g},\xi}^{1}(\mathbf{x})=\Lambda_{\mathbf{g}}(\rho(\xi)\mathbf{ x}),\ \Gamma_{\mathbf{g},\xi}^{2}(\mathbf{w},b)=c_{j}S_{j}(\xi,b)\Lambda_{\mathbf{g}}(\eta( \xi)\mathbf{w}) \tag{9}\]

That observation directly leads to the RF mechanism for the estimation of \(\widehat{f}_{j}(\mathbf{x},\mathbf{w},b)\). We can rewrite: \(\widehat{f}_{j}(\mathbf{x},\mathbf{w},b)=\mathbb{E}[\Phi^{j}(\mathbf{x})^{ \top}\Psi^{j}(\mathbf{w},b)]\) for \((\xi_{1},\mathbf{g}_{1}),...,(\xi_{m},\mathbf{g}_{m})\sim\overline{\mathcal{P} }_{j}\otimes\mathcal{N}(0,\mathbf{I}_{d})\) and:

\[\begin{split}\Phi^{j}(\mathbf{x})=\frac{1}{\sqrt{m}}(\Gamma_{ \mathbf{g}_{1},\xi_{1}}^{1}(\mathbf{x}),...,\Gamma_{\mathbf{g}_{m},\xi_{m}}^{1 }(\mathbf{x}))^{\top},\\ \Psi^{j}(\mathbf{w},b)=\frac{1}{\sqrt{m}}(\Gamma_{\mathbf{g}_{1}, \xi_{1}}^{2}(\mathbf{w},b),...,\Gamma_{\mathbf{g}_{m},\xi_{m}}^{2}(\mathbf{w},b ))^{\top}\end{split} \tag{10}\]

Several strategies can be used to construct samples \((\xi_{1},\mathbf{g}_{1}),...,(\xi_{m},\mathbf{g}_{m})\), e.g. iid sampling or block-iid sampling with a fixed \(\xi\) used within a block, but constructed independently for different blocks. In the experiments, we also choose: \(\rho(\xi)=2\pi i\xi\) and \(\eta(\xi)\)=1.

The case of discrete \(\overline{\mathcal{P}}_{j}\) with finite number of atoms:Assume that \((\xi^{1},...,\xi^{K})\) is a sequence of atoms with the corresponding positive probabilities: \((p_{1},...,p_{K})\). Then one can also construct \(K\) pairs of RF-vectors \((\Phi^{j}(\mathbf{x};k),\Psi^{j}(\mathbf{w},b;k))_{k=1}^{K}\), each obtained by replacing \(\overline{\mathcal{P}}_{j}\) with a distribution corresponding to a deterministic constant \(p_{i}\) and get \(\Phi^{j}(\mathbf{x}),\Psi^{j}(\mathbf{w},b)\) by concatenating vectors from \((\Phi^{j}(\mathbf{x};k))_{k=1}^{K}\) and \((\Psi^{j}(\mathbf{w},b;k))_{k=1}^{K}\) respectively. This strategy is effective if \(K\) is small.

Note that: \(f(\mathbf{x}^{\top}\mathbf{w}+b)=\sum_{j=0}^{3}\widehat{f}_{j}(\mathbf{x}, \mathbf{w},b)\) and thus \(\Phi_{f}(\mathbf{x})\) and \(\Psi_{f}(\mathbf{w},b)\) can be defined as:

\[\Phi_{f}(\mathbf{x})=\mathrm{concat}\left(\left(\Phi^{j}(\mathbf{x})\right)_{j=0 }^{3}\right),\ \Psi_{f}(\mathbf{w},b)=\mathrm{concat}\left(\left(\Psi^{j}(\mathbf{w},b) \right)_{j=0}^{3}\right) \tag{11}\]

for the vector concatenation operation \(\mathrm{concat}\), completing the description of the URF mechanism.

**Remark 3.1** (boundedness): _We observe that for upper-bounded \(\|\mathbf{x}\|_{2},\|\mathbf{w}\|_{2},|b|\), the entries of \(\Phi_{f}(\mathbf{x})\) and \(\Psi_{f}(\mathbf{w},b)\) are also upper-bounded as long as \(A<0\). This follows directly from the formula for \(\Lambda_{\mathbf{g}}(\mathbf{z})\) in Eq. 8._
Trigonometric activation functions:Let us assume now that \(f(z)=\sin(z)\) or \(f(z)=\cos(z)\). Note that even though none of them has a Fourier Transform in the classical Riemannian sense, both have trivial Fourier Transforms in the broader distributional sense. To see that, we can rewrite both activations as: \(\sin(z)=\frac{exp(iz)-\exp(-iz)}{2i}\) and \(\cos(z)=\frac{exp(iz)+\exp(-iz)}{2}\). Therefore the corresponding distributions used in the URF derivations above become binary distributions over \(\{-\frac{1}{2\pi},\frac{1}{2\pi}\}\). This observation has interesting practical consequences, since it leads to the conceptually simple linearization of the FFLs applied in SIREN networks (see: Sec. 4.2).

### Beyond regular FFLs: the curious case of the ReLU-SNNK layer

We also propose another SNNK layer which is not directly inspired by any known FFL, but turns out to work very well in practice (see: Sec. 4.3.2). In this case, the mappings \(\Phi\) and \(\Psi\) are defined as: \(\Phi(\mathbf{x})=\mathrm{ReLU}(\frac{1}{\sqrt{l}}\mathbf{G}\mathbf{x})\), \(\Psi(\mathbf{w},b)=\mathrm{ReLU}(\frac{1}{\sqrt{l}}\mathbf{G}\mathbf{w})\) for the Gaussian matrix: \(\mathbf{G}\in\mathbb{R}^{l\times d}\) with entries sampled independently at random from \(\mathcal{N}(0,1)\). One can ask a question what kernel does this pair of maps correspond to. It turns out that the answer is particularly elegant.

**Theorem 3.2** (arc-cosine kernels; Cho & Saul (2011)): _The nth-order arc-cosine kernel \(\mathrm{K}_{n}:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}\) is defined as: \(\mathrm{K}_{n}(\mathbf{x},\mathbf{y})=\frac{1}{n}\|\mathbf{x}\|_{n}^{2}\| \mathbf{y}\|_{2}^{2}J_{n}(\alpha_{\mathbf{x},\mathbf{y}})\), where \(\alpha_{\mathbf{x},\mathbf{y}}\in[0,\pi]\) stands for an angle between \(\mathbf{x}\) and \(\mathbf{y}\) and \(J(\theta)\stackrel{{\mathrm{def}}}{{=}}(-1)^{n}(\sin(\theta))^{n +1}\frac{\partial^{\pi}}{\partial\theta^{\pi}}\left(\frac{\pi-\theta}{\sin( \theta)}\right)\). Then, \(\mathrm{K}_{n}\) can be linearized as: \(\mathrm{K}_{n}(\mathbf{x},\mathbf{y})=2\mathbb{E}[\Gamma_{n}(\mathbf{x})^{\top }\Gamma_{n}(\mathbf{y})]\) for \(\Gamma_{n}(\mathbf{v})\stackrel{{\mathrm{def}}}{{=}}\mathrm{ReLU }((\mathbf{v}^{\top}\boldsymbol{\omega})^{n})\) and \(\boldsymbol{\omega}\thicksim\mathcal{N}(0,\mathbf{I}_{d})\)._

We conclude that our proposed \(\mathrm{ReLU}\)-\(\mathrm{SNNK}\) layer is a scalable version of the FFL defined as: \(\mathbf{x},\mathbf{W}\in\mathbb{R}^{l\times d},\mathbf{b}\rightarrow(\frac{1} {2}\mathrm{K}_{1}(\mathbf{w}^{1},\mathbf{x}),...,\frac{1}{2}\mathrm{K}_{1}( \mathbf{w}^{l},\mathbf{x}))^{\top}\).

**Remark 3.3**: _The \(\mathrm{ReLU}\)-\(\mathrm{SNNK}\) layer is not a regular FFL since the values of its output dimensions cannot be re-written as \(f(\mathbf{x}^{\top}\mathbf{w}^{i}+b_{i})\) for some \(f:\mathbb{R}\rightarrow\mathbb{R}\) (interestingly, after \(\Gamma\)-base pre-processing, it can be still interpreted as a dot-product kernel). This shows that the SSNK mechanism is capable of modeling relationships beyond those of regular FFLs._

### Bundling neural networks with SNNKs

We are ready to propose the neural network _bundling process_, relying on the SSNK-primitives. Consider the following deep NN module with input \(\mathbf{x}=\mathbf{x}_{0}\in\mathbb{R}^{d_{0}}\) and output \(\mathbf{y}=\mathbf{x}_{L}\in\mathbb{R}^{d_{L}}\):

\[\begin{cases}\mathbf{x}_{i+1}=f_{i+1}(\mathbf{W}_{i}\mathbf{x}_{i}+\mathbf{b} _{i})\text{; }i=0,...,L-1,\\ \mathbf{x}_{0}=\mathbf{x}\end{cases} \tag{12}\]

for: (1) matrices \(\mathbf{W}_{i}\in\mathbb{R}^{d_{i+1}\times d_{i}}\), (2) bias vectors: \(\mathbf{b}_{i}\in\mathbb{R}^{d_{i+1}}\), and (3) activations: \(f_{i}:\mathbb{R}\rightarrow\mathbb{R}\).

To understand how the bundling process works, we start by replacing first FFL in Eq. 12 with its SNNK analogue. We obtain the following computational block:

\[\begin{cases}\widehat{\mathbf{x}}_{i+1}=\widehat{f}_{i+1}(\widehat{\mathbf{W} }_{i}\widehat{\mathbf{x}}_{i}+\widehat{\mathbf{b}}_{i})\text{ for }i=0,...,L-2,\\ \widehat{\mathbf{x}}_{0}=\Phi_{f_{1}}(\mathbf{x}_{0}),\\ \widehat{\mathbf{W}}_{0}=\mathbf{W}_{1}\Psi_{f_{1}}(\mathbf{W}_{0},\mathbf{b }_{0})\text{; }\widehat{\mathbf{W}}_{i}=\mathbf{W}_{i+1}\text{ for }i=1,...,L-2,\\ \widehat{f}_{i+1}=f_{i+2},\ \widehat{\mathbf{b}}_{i}=b_{i+1}\text{ for }i=0,...,L-2\end{cases} \tag{13}\]

In the system of equations above, \(\Psi_{f}(\mathbf{W}_{0},\mathbf{b}_{0})\) is a matrix with transposed rows of the form: \(\Psi_{f}(\mathbf{W}_{0}^{j},\mathbf{b}_{0}^{j})\), where \(\mathbf{W}_{0}^{j}\) for \(j=0,...,d_{1}-1\) are the transposed rows of \(\mathbf{W}_{0}\) and \(\mathbf{b}_{0}=(\mathbf{b}_{0}^{0},...,\mathbf{b}_{0}^{d_{1}-1})^{\top}\). We have thus successfully replaced a module of \(L\) feedforward layers with a module of \((L-1)\) feedforward layers. By continuing this procedure, we can ultimately get rid of the FFLs and obtain an estimator \(\overline{\mathbf{y}}\) of \(\mathbf{y}\), given as: \(\overline{\mathbf{y}}=\overline{\mathbf{W}\mathbf{x}}\), where

\[\begin{cases}\overline{\mathbf{x}}=\Phi_{f_{L}}\left(\Phi_{f_{L-1}}(...\Phi_{ f_{1}}(\mathbf{x}_{0})...)\right)\\ \overline{\mathbf{W}}=\Psi_{f_{L}}\left(\mathbf{W}_{L-1}\Psi_{f_{L-1}}(... \mathbf{W}_{2}\Psi_{f_{2}}(\mathbf{W}_{1}\Psi_{f_{1}}(\mathbf{W}_{0},\mathbf{b }_{0}),\mathbf{b}_{1})...,),\mathbf{b}_{L})\right)\in\mathbb{R}^{d_{L}\times m} \end{cases} \tag{14}\]

This has several important consequences. In inference, replacing matrices \(\mathbf{W}_{0},...,\mathbf{W}_{L-1}\) with one matrix \(\overline{\mathbf{W}}\) is a effectively a compression scheme (that does not necessarily need to be applied to
all the layers, but a particular consecutive set of layers of interest). If we apply bundling process to the entire deep neural network, we effectively provide its two-tower factorization with input disentangled from the parameters. In training, we can treat \(\overline{\mathbf{W}}\) as an unstructured parameter matrix and directly learn it (see results in Appendix H.3, table 5). Since the output \(\overline{\mathbf{y}}\) is now modeled as an action of the unstructured learnable matrix \(\overline{\mathbf{W}}\) on the _pre-processed_ input \(\overline{\mathbf{x}}\), for several loss functions there exists an explicit formula for the optimal \(\overline{\mathbf{W}}\). This is the case in particular for the standard regression loss (see discussion in Appendix H.3). If bundling is applied to a particular module, backpropagation through it is not necessary since there exists an explicit formula for the corresponding Jacobian.

## 4 Experiments

We present an extensive empirical evaluation on SNNK on a wide range of experiments. More details on each of the experiments can be found in the Appendix E.

### Pointwise Kernel Estimation

As a warm-up, we test the accuracy of the applied RF-mechanisms on synthetic data. We take \(d=2000\) and \(l=1\). We consider: **(a)** a SIREN-FFL with the activation function \(f(u)=\sin(u)\) and bias \(b=0.5\), **(b)** an arc-cosine-FFL from Sec. 3.2. The entries of the weight vectors \(\mathbf{w}\) and the inputs to the layers are taken independently from \(\frac{1}{\sqrt{d}}\text{Unif}(0,1)\). We report the mean relative error of the NN output (by averaging over \(s=500\) instantiations of the RF-mechanism) made by the RF-based estimator as well as the empirical standard deviation as a function of the number of random projections. This setup corresponds to quantifying the accuracy of the kernel estimator pointwise. The results are presented in Fig. 2 (g and h). Our SNNK provided an accurate approximation with a much smaller number of random projections than the dimensionality \(d\) of the input vectors.

Figure 2: Architecture for **(a)** SNNK layer (see Section A), **(b)** SNNK-Adpt layer **(c)** image fitting (SIREN), MNIST and UCI experiments, **(d)** SNNK-QPNN model, **(e)** SNNK-inspired Adapter-ViT layer, **(f)** SNNK-inspired Adapter-BERT layer. **(g,h)**: The relative error (obtained by averaging over \(s=500\) instantiations of the RF-mechanism) made by the RF-based estimator on the particular entry of the output of the. **(g)** SIREN-FFL and **(h)** are-cosine-FFL as a function of the number of random projections \(p\) (see: Sec. 4.1). The maximum \(p\) for (g) is larger than for (h), as (g) in theory produces larger variance per random projection. The corresponding standard deviations are negligible: **(g)**\(5\cdot 10^{-8}\),\(10^{-12}\), \(5\cdot 10^{-8}\), \(10^{-8}\),\(10^{-12}\), \(2.5\cdot 10^{-9}\), \(10^{-12}\), \(5\cdot 10^{-9}\), \(10^{-12}\), \(10^{-12}\), \(10^{-10}\), \(10^{-12}\), **(h)**\(10^{-12}\), \(3\cdot 10^{-8}\),\(3\cdot 10^{-8}\),\(2\cdot 10^{-8}\), \(10^{-12}\), \(5\cdot 10^{-9}\).
### Toy Experiments

SNNKs are versatile and can be used as a drop-in replacement for FFLs in a wide variety of NNs like the SIREN network Sitzmann et al. (2020), QPNN - a physics-inspired Neural Network (PINN) to solve the Hamiltonian for quantum physical systems (Sehanobish et al., 2021) and a simple multi-layer perceptron (MLP) for classification on MNIST (LeCun and Cortes, 2010). We use the sine activation variant for the first two experiments and the ReLU variant for MNIST. \(32\) random features are used for the solution of the 2-body problem and MNIST and 64 random features for the image fitting problem. We match the performance of the baseline NNs on the 2-body and the image fitting problem (see figure 3) and outperform the baseline on MNIST (Figure 9), while incurring lower training costs. For additional details regarding these experiments, see Appendix E.1.

### Finetuning Experiments

In this subsection, we show how SNNKs can be used for parameter efficient finetuning. For experiments on text, we use the GLUE benchmark consisting of 8 different natural language understanding tasks (Wang et al., 2018). For vision tasks, we use CiFAR-10 and CiFAR-100 (Krizhevsky et al., 2009). BERT-base (Devlin et al., 2019) is used as a backbone for text experiments and ViT (Kolesnikov et al., 2021) for image experiments. Our code is built on top of Transformers (Wolf et al., 2020) and adapter Transformer library (Pfeiffer et al., 2020). Detailed comparisons with various baselines can be found in Appendix I and additional experiments can be found in Appendix H.

#### 4.3.1 Linearizing the Pooler Layer in Transformers

For text classification tasks, a SNNK layer can be used as a drop-in replacement for the pooler layer which is a linear layer with a tanh activation. For these set of experiments, the base model is frozen and only the pooler and the classifier weights are tuned. We get computational gains as the number of random features employed by SNNK is smaller than that of the hidden size of the Transformers. More details are presented in Appendix E.2.

On GLUE dev set, our SNNK-linearized pooler models outperform the baselines on 5 out of 8 tasks (Table 1 (top half)). Additional results can be found in Appendix H.

In this setting, the linearized pooler weights can be merged with the classifier weights to create a weight matrix of size (# number of random features \(\times\) number of classes) and then one can simply store the newly merged layer instead of separately storing the trained classifier and pooler layers. This dramatically _reduces_ the storage from 18.92 Megabit to only.**02** Megabit leading to a compres

Figure 4: Comparison of trainable parameters between various layers/modules and the drop in replacement NNK layers. Results for CiFar-10 and CiFar-100 for SNNK-Adapter models.

Figure 3: (1) **Left column** : Injecting SNNK in a PINN network to approximate the potential energy of the 2-body system. Top to bottom : Ground truth potential, Learned potential by QPNN (Sehanobish et al., 2021) and QPNN-SNNK. QPNN-SNNK can learn the potential function perfectly even using less trainable parameters than the baseline QPNN. (2) **Rightmost three column** : Siren network on the first row, fitting not only the image, but also the gradients. SNNK on the bottom row produces an accurate approximation of the above.
sion factor of \(\mathbf{1/1000}\). More details are presented in Appendix C. Ablation studies on the number of random parameters for this experimental setting are presented in Appendix G.

#### 4.3.2 SNNK-inspired Adapter Layers

Adapters in Transformers were first introduced in (Houlsby et al., 2019) and there has been a lot of work designing different architectures (Pfeiffer et al., 2020; Karimi Mahabadi et al., 2021; Moosavi et al., 2022) and unifying various paradigms (Moosavi et al., 2022; He et al., 2022). Adapters are bottleneck MLPs which are (generally) added twice to each Transformer layer. In our work, we replace each adapter block by a single SNNK layer (Figure 2 (e) and (f)) using only \(\mathbf{16}\) random features resulting in a big drop of training parameters (see Figure 4). Figure 4 (b) shows the architecture of SNNK-inspired adapter layers. Additional details are presented in Appendix B.

As is customary for adapter experiments, base model is frozen and only the adapters and classifier are tuned. Table 1 (bottom half) shows our results on using SNNK layers in place of adapters on the GLUE dev set. We outperform the baseline on \(5\) out of \(8\) datasets while employing only \(\mathbf{1/3}\) of the training parameters. On MNLI, it is noted in (Houlsby et al., 2019), that using smaller adapter size causes worse performance and performance boost can be achieved by increasing the size of the adapter (\(256\) is used in their case). Similar to this observation, we note that we can improve performance and match the baselines on large datasets (ex. MNLI, QNLI) as we increase the number of random features (see Figure 5). Our method also produces competitive performance on image datasets like Cifar-10 and Cifar-100 (see Figure 4 (right 2 figures)). Detailed comparisons with SOTA parameter efficient finetuning methods can be found in Table 6 (vision tasks) and in Table 7 (GLUE tasks).

Moreover, we note that our methods are completely orthogonal to techniques such as gating mechanism in (Mao et al., 2022) or algorithms relying on dropping suitable adapter layers (Moosavi et al., 2022; Ruckle et al., 2021). Thus it can be easily combined with them.

### Experiments on UCI datasets

We have conducted experiments with a variety of real-world datasets found in the UCI Machine Learning Repository (UCI MLR).1. We trained a three-layer MLP model as baseline (see Appendix Sec. F.5 for details). We varied the output of the middle-layer to train MLPs with different sizes. For

Figure 5: Ablation with different number of random features for the ReLU-SNNK-adapter experiments on the GLUE dev set. \(AA\) is the reported adaptable adapter numbers in Moosavi et al. (2022).

\begin{table}
\begin{tabular}{l c c c c c c c c} \hline Dataset & RTL & MRPC & QUNLI & QQP & SST2 & MNLI & STSB & COLA \\ \hline End-baseline (J.e. et al., 2019) & 52.5 & 51.3 & 74.5 & **72.0** & 51.9 & **56.4** & 75.1 & 20.4 \\ \hline End-baseline (SNNK-based (ours)) & \(\mathbf{61.36\pm 1.15}\) & \(\mathbf{82.07\pm 1.07}\) & \(\mathbf{76.52\pm 0.27}\) & \(\mathbf{76.01\pm 0.17}\) & \(\mathbf{65.21\pm 0.34}\) & \(\mathbf{56.02\pm 0.27}\) & \(\mathbf{76.82\pm 0.27}\) & \(\mathbf{26.81\pm 0.30}\) \\ \hline End-baseline (J.e. et al., 2019) & \(\mathbf{61.39\pm 1.27}\) & \(\mathbf{51.90\pm 1.06}\) & \(\mathbf{90.30\pm 0.25}\) & \(\mathbf{88.09\pm 0.16}\) & \(\mathbf{91.31\pm 0.31}\) & \(\mathbf{52.09\pm 0.47}\) & \(\mathbf{83.52\pm 0.17}\) & \(\mathbf{51.41\pm 1.28}\) \\ \hline End-baseline (J.e. et al., 2019) & \(\mathbf{60.68\pm 1.24}\) & \(\mathbf{91.26\pm 1.39}\) & \(\mathbf{90.44\pm 0.16}\) & \(\mathbf{58.52\pm 0.23}\) & \(\mathbf{92.31\pm 0.27}\) & \(\mathbf{52.06\pm 0.17}\) & \(\mathbf{88.81\pm 0.14}\) & \(\mathbf{58.21\pm 0.63}\) \\ \hline \end{tabular}
\end{table}
Table 1: SNNK experiments on GLUE benchmarks. MCC score is reported for CoLA, F1 score is reported for MRPC and QQP, Spearman correlation is reported for STSB. Accuracy scores are reported for the other tasks. All results are obtained by averaging over 5 seeds.
our method, we replace the middle-layer with SNNK (Figure 2 (c)). SNNK matches or outperforms the baseline while using only a fraction of the training parameters (Figure 6).

## 5 Conclusion

We present scalable neural network kernels (SNNK), a novel efficient NN computational model that can be used to replace regular feedforwards layers in MLPs, where inputs and parameters are disentangled and connected only in the final computation via a dot-product kernel. We introduce a general mechanism of the universal random features (URFs) to instantiate SNNKs, show that SNNKs are capable of encoding subtle relationships between parameter- and input-vector beyond functions of their dot-products and finally, explain how they lead to the compactification of the NN stack via the so-called bundling process. We complement our theoretical findings with the exhaustive empirical analysis, from pointwise kernel estimation to training Transformers with adapters.

## Ethics Statement

This paper focuses mostly on the algorithmic properties of the techniques linearizing kernels associated with feedfoward layers' calculations for the computational gains. The experiments with adapter-based fine-tuning of Transformers are presented to illustrate the main concepts. It should be noted though that Transformers should be used cautiously given their considerable computational footprint (improved though while adapters are applied) and the corresponding carbon footprint.

## Reproducibility Statement

Hyperparameters to reproduce each experiment is detailed in section F. All the code to reproduce the experiments will be provided on acceptance of this preprint.

## Author Contributions

AS and KC led the project. AS ran several empirical studies on GLUE, CiFar-10 and CiFar-100 datasets and proposed several strategies for efficiently using SNNK-layers within Transformer models. KC proposed FFL linearization-schemes, URFs, the bundling mechanism and implemented all linearization-schemes. YZ ran empirical studies on GLUE, CiFar-10 and CiFar-100 datasets. AD implemented and ran all UCI experiments, helped with GLUE/image experiments, proposed strategy for efficiently using SNNK-layers and created all figures in experiments. VL proposed an idea to linearize FFLs by disentangling inputs from weights. All authors contributed to the writing of the manuscript.
