# Fully-fused Multi-Layer Perceptrons on Intel Data Center GPUs

Kai Yuan2*, Christoph Bauinger2*, Xiangyi Zhang2*,
Pascal Baehr2, Matthias Kirchhart2, Darius Dabert3, Adrien Tousnakhoff3,
Pierre Boudier2 and Michael Paulitsch2  2 Intel Corporation  3 École
Polytechnique

###### Abstract

This paper presents a SYCL implementation of Multi-Layer Perceptrons (MLPs),
which targets and is optimized for the Intel Data Center GPU Max 1550. To
increase the performance, our implementation minimizes the slow global memory
accesses by maximizing the data reuse within the general register file and the
shared local memory by fusing the operations in each layer of the MLP. We show
with a simple roofline model that this results in a significant increase in
the arithmetic intensity, leading to improved performance, especially for
inference. We compare our approach to a similar CUDA implementation for MLPs
and show that our implementation on the Intel Data Center GPU outperforms the
CUDA implementation on Nvidia’s H100 GPU by a factor up to 2.84 in inference
and 1.75 in training. The paper also showcases the efficiency of our SYCL
implementation in three significant areas: Image Compression, Neural Radiance
Fields, and Physics-Informed Machine Learning. In all cases, our
implementation outperforms the off-the-shelf Intel Extension for PyTorch
(IPEX) implementation on the same Intel GPU by up to a factor of 30 and the
CUDA PyTorch version on Nvidia’s H100 GPU by up to a factor 19. The code can
be found at https://github.com/intel/tiny-dpcpp-nn.

###### Index Terms:

Machine Learning, Performance Optimization, SYCL, Intel Data Center GPU Max
1550

**footnotetext: Equal contribution

##  I Introduction

Multi-Layer Perceptrons (MLPs) [1] play a vital role in today’s Machine
Learning (ML) and Artificial Intelligence (AI) landscape next to the prevalent
Transformer architecture [2] and Convolutional Neural Networks [3], which are
mostly used in the Natural Language Processing and Computer Vision domains.
MLPs are used as the main Neural Network architecture for several Machine
Learning applications, such as the representation of the solution operator of
partial differential equations [4], the density or color function in Neural
Radiance Fields (NeRFs) objects [5], and replacing classical ray-tracing with
Neural Ray Tracing [6] (see Section II for details). In contrast to the
aforementioned architectures, MLPs are characterized by their fully connected
layers, where the neuron of every layer is connected to every neuron in the
previous and next layer. A key distinguishing factor is that in MLPs, each
neuron’s output is independent of its neighbors in the same layer, making it
suitable for fully-fusing operations as described in this work.

The present contribution focuses on the efficient implementation on Intel GPUs
of “narrow” MLPs, which consist of an arbitrary number of layers (depth), and
a small and constant number of neurons per layer (width). These narrow MLPs
are of particular interest since i) they are universal approximators [7], ii)
they have several relevant use cases (see Sec. II), and iii) their theoretical
peak performance is severely limited by the small width of the layers. More
precisely, their small width results in a reduced arithmetic intensity of the
requisite matrix-multiplications in each layer for the training and, in
particular, the inference. Thus, in a “classic” implementation of MLPs, where
necessary operations in each layer are performed in separate compute kernels,
the performance of narrow MLPs is typically severely bound by the memory
bandwidth of the global memory or the last level cache.

To alleviate the issues arising from the low arithmetic intensity and the
memory bandwidth of the global memory, a common strategy [6] is the fusion of
the layers into a single kernel to keep relevant data in faster memories,
i.e., the register file, shared memory or faster caches. This approach, termed
“fully-fused MLPs”, has been implemented for Nvidia GPUs utilizing Nvidia’s
CUDA language [8].

In this contribution we focus on a SYCL implementation for Intel GPUs of
fully-fused MLPs with arbitrary depth and fixed layer width of
$2^{i},i\in\\{4,...,7\\}$ neurons in each layer. Note that, as indicated
above, fixed widths are not limiting the expressive power of the Neural
Network as fixed-width MLPs are still universal approximators, i.e., they can
approximate any continuous function to any desired accuracy as proven by the
Universal Approximation Theory for width-bounded networks [7]. Indeed, in
practise MLPs rarely exceed the maximum width of $2^{7}$ elements supported by
this work as networks tend to be deeper to gain more expressiveness rather
than wide (see Section II).

Our implementation of fully-fused MLPs is based on Intel’s joint_matrix SYCL
extension [9] to utilize the XMX hardware [10] in Intel’s Data Center GPU Max
1550 [11], which is the device targeted with our optimized implementation.

Our method is especially well-suited to optimize the training and inference
performance for models that require large data throughput with batch sizes
$2^{i},15<i\in\mathbb{N}$, since those sizes maximize the occupancy of the
device. As shown in Section IV-A, our SYCL implementation on Intel hardware
has improved performance over an equivalent CUDA implementation (tiny-cuda-nn
[8]) for MLPs with width 64 by a factor up to 2.84 in inference and 1.75 in
training.

Furthermore, we argue with a roofline analysis (see Section III-C) that our
approach to fully-fused MLPs is especially well-suited for the acceleration of
the inference and that it significantly increases the arithmetic intensity,
and thus the theoretical peak performance, compared to the approach shown in
[6] by reducing the accesses to global memory.

To further show the performance improvements and potential applications for
our implementation, we demonstrate the performance on a regression benchmark
and the following three applications: Image Compression, Neural Radiance
Fields (NeRFs), and Physics-Informed Machine Learning (see Section IV). In all
cases, our implementation outperforms the off-the-shelf Intel Extension for
PyTorch (IPEX) [12] implementation on the same Intel GPU by up to a factor of
30 and the CUDA PyTorch version on Nvidia’s H100 GPU [13] by up to a factor
19.

Summarised, the contributions of this paper are:

  * 1.

The first SYCL implementation of fully-fused MLPs applied on Intel GPU’s that
support XMX instructions and open-sourced repository of the implementation.

  * 2.

A roofline model of our implementation and comparison to the roofline of the
fully-fused implementation [6]. We argue an improvement of the arithmetic
intensity of up to a factor 2.15.

  * 3.

Demonstrating higher performance on four example applications: regression
benchmark, image compression, Neural Radiance Fields, and Physics-Informed
Neural Networks. Demonstrating an improvement of up to 1.75x and 2.84x
performance increase for training and inference respectively over another
fully-fused implementation, and a performance increase of a factor up to 30
over off-the-shelf PyTorch implementations.

The following sections are structured as follows. First, we outline the
applications of MLPs and their impact. Next, the fully-fused MLPs are
described. Finally, we demonstrate our results on four example applications,
and conclude this paper.

##  II Applications of Multi-Layer Perceptrons

To show the importance of MLPs in today’s AI landscape, this section reviews
the applications of MLP variants and implementations in multiple areas. Please
note, that due to the fast-paced adoption of AI and ML in various fields and
industries, this list is non-exhaustive and should rather provide the reader
with a qualitative intuition about the importance of accelerating MLPs and
their impact to the community. For a more comprehensive overview, please refer
to [14, 15, 16].

Further, to provide a better understanding of what performance gains can be
achieved in key areas, we highlight and provide more details on three
prominent applications in the Appendix -A--C: Neural Rendering, Model
Compression, and Partial Differential Equations. The performance and results
generated with our implementation for these applications are shown in Section
IV.

Figure 1: Taxonomy categorizing various applications of MLP. Each category has
multiple subcategories that illustrate the specific tasks and domains that MLP
can address.

Most MLP applications can be categorized into the following four categories:

  * •

Recognition: using MLPs to identify objects, or patterns.

  * •

Representation learning: aiming to extract meaningful patterns from raw data
to create representations that are easier to understand and process by MLPs.

  * •

Reinforcement learning: utilizing MLPs to learn from interactions with the
environment to make better decisions.

  * •

Regression: predicting values based on input data.

The taxonomy in Figure 1 categorizes contemporary MLP applications, which
demonstrates the diversity and versatility of MLP in solving various real-
world problems. We present these categories in the following sections in some
detail.

###  II-A Recognition

In vision, MLPs are applied to image classification [17, 18, 19, 20, 21],
objects detection and semantic segmentation [22, 23, 24, 25, 26], video
analysis [27, 28, 29], image processing and generation [30, 31, 32, 33]. Other
notable works in NLP include sentiment analysis and opinion classification
[34, 35, 36], malware classification [37, 38, 39], adversarially robust [40]
and multilingual transition [41] which leverages MLPs for language
understanding and generation.

###  II-B Representation Learning

MLPs may solve partial differential equations through Physics-Informed Neural
Networks (PINNs) [42, 43, 44, 4, 45] and neural differential operators [46,
47, 48], which approximate the solution to differential equations using neural
networks.

In the field of neural graphics, MLPs have been developed and integrated into
neural rendering [5, 49, 50, 51, 52, 53], 3D reconstruction [54, 55, 56],
super-resolution [57] and depth-estimation [58, 59]. In this work, we are
interested in exploiting the potential of MLPs in the above three fields and
will show their results in the Appendix -A--C for details.

In the field of Compression, MLPs aim to reduce the size of data or models by
removing redundant or irrelevant parameters and operations. MLPs are utilized
to perform data compression by learning compact representations of data that
can be reconstructed with minimal distortion[60, 61]. In addition, MLPs can be
utilized for model compression by learning to approximate the function of
large or complex models with smaller or simpler models [62, 63].

###  II-C Reinforcement Learning

In the field of robotics, MLPs are integral for various tasks, including
navigation [64, 65, 66, 67], autonomous-driving [68], and manipulator [69]. In
system research, MLPs have been employed for structured control nets [69, 70],
recommender systems [71, 72], pure-feedback systems [73, 74] and Gaming [75,
76]. which utilize MLP-based deep learning techniques to discover system
design and control.

###  II-D Regression

Most existing works using MLPs on regression are in the fields of social
media, biochemistry, and Energy.

In social media, MLPs play a pivotal role in sentiment analysis [77, 78, 79,
80], financial analysis [81, 82] and fraud detection[83, 84], for analysis the
social problems. Furthermore, in biochemistry research, MLPs are able to model
complex biology and chemistry to explore the molecular basis of life in radial
basis function [85] and bioactivity estimation [86]. In the field of energy
and utilities, MLPs work well on load forecasting [87], energy consumption
[88], and machine optimization [89] to deal with some biological issues and
improve the effect of activities.

##  III Fully-fused MLPs on Intel GPUs

In the present section, we introduce our approach to the implementation of
fully-fused MLPs on Intel GPUs. Section III-A we describe the inference and
training algorithms. Section III-B then describes our SYCL joint_matrix
implementation for the Intel Data Center GPU Max 1550 in detail.

###  III-A Inference and Training Operations

For the sake of conciseness and simplicity, we describe our implementation of
fully-fused MLPs for the special case where the width of each layer, including
input and output, is 64. Cases where the input and output size differs may be
reduced to this special case by utilizing, e.g., encodings for the input data
or appropriate padding of the data and the weights.

The matrix operations during training and inference are described in what
follows. Let $M\in\mathbb{N}$ denote the batch size and let $N=K=64$ denote
the layer width. Let $\sigma:\mathbb{R}\to\mathbb{R}$ denote an activation
function [90]. We use the notation $\sigma(A)$, for a matrix $A$, to indicate
the element-wise application of the activation function to the matrix $A$. Let
$2\leq\mathrm{nlayers}\in\mathbb{N}$ denote the number of layers in the MLP,
which includes the input and the output layer. I.e., it denotes the number of
matrices as in Alg. 1. Further, let $A_{i}\in\mathbb{R}^{M\times K}$,
$B_{i}\in\mathbb{R}^{K\times N}$, and $C_{i}\in\mathbb{R}^{M\times N}$,
$D_{i}\in\mathbb{R}^{M\times N}$, $G_{i}\in\mathbb{R}^{K\times N}$ denote
matrices for each $i=1,\ldots,\mathrm{nlayers}$. The transpose of any matrix
$A$ is given by $A^{T}$. An element of a matrix $A$ with the row coordinate
$r$ and the column coordinate $c$ is given by $A(r,c)\in\mathbb{R}$. The
matrices $B_{i}$ denote weights of each layer $i$, $1\leq
i\leq\mathrm{nlayers}-1$.

A pseudo code of the inference is shown in Algorithm 1, which consists of
repeated matrix multiplications $A_{i}B_{i}$ and applications of the
activation function $\sigma$.

Algorithm 1 Inference

$\mathrm{nlayers}$, $\sigma$, $\mathrm{Input}$,
$B_{1},\ldots,B_{\mathrm{nlayers}-1}$

Initialize: $A_{1}=\mathrm{Input}$

for $i\leftarrow 1$ to $\mathrm{nlayers}-1$ do

$C_{i}\leftarrow A_{i}B_{i}$

$A_{i+1}\leftarrow\sigma(C_{i})$

end for

Return $A_{\mathrm{nlayers}}$

Our implementation of the inference step is based on several observations.
First, each weight matrix $B_{i}$ fits in the shared local memory (SLM) [91]
of the targeted Intel Data Center GPU Max 1550. The loads of the $B_{i}$
matrices from global memory, i.e., high bandwidth memory (HBM), can therefore
be minimized by pre-fetching them one after the other into SLM. Secondly, up
to eight rows of the matrices $A_{i}$ and $C_{i}$ (for a single $i$) fit into
the general register file (GRF) for all network widths which are 64 or less.
For network widths between 64 and 128, it still holds but requires large GRF
mode [92]. This is relevant since the targeted XMX hardware may be applied to
sub-matrices with up to eight rows. However, as analyzed in the following Sec.
III-C, the highest arithmetic intensity is achieved when using sub-matrices
consisting of exactly eight rows. Lastly, only $A_{i}$ and $C_{i}$ for a
single $i$ are required at the same time and only the last matrix
$A_{\mathrm{nlayers}}$ is returned. The matrices associated with other layers
$j$, $1\leq j\leq\mathrm{nlayers}-1$ are discarded during inference thus
minimizing the amount of memory accesses required, consequently increasing the
arithmetic intensity and performance.

The pseudo code for the training of the network is shown in Algorithm 2. In
particular, we focus on the forward pass and the subsequent backward pass
after a loss is calculated. We do not consider the optimization step in this
pseudo code. In contrast to the inference, the training takes as input two
activation functions. One for the forward pass and one for the backward pass
denoted by $\sigma_{f}$ and $\sigma_{b}$, respectively. Note, that the
activation function $\sigma_{b}$ is the derivative of $\sigma_{f}$ resulting
from the chain rule during back-propagation [93]. In addition, the loss
computation requires target values given as an $M\times K$ array of real
values, named $\mathrm{Target}$. The outputs of the training algorithm are
$\mathrm{nlayers}-1$ matrices $G_{1},\ldots,G_{\mathrm{nlayers}-1}$. The
matrices $G_{i}$ represent the gradients of the weight matrices $B_{i}$. The
gradients are required for the subsequent optimization step, which is
typically some type of gradient descent algorithm and is a real square matrix
of size $K\times N$ (note that $K=N$ as we have fixed layer widths of
$2^{i},i=4,...,7$) .

The forward pass in Algorithm 2 is nearly identical to the inference. The only
difference is that all matrices $A_{i}$, $i=1,\ldots,\mathrm{nlayers}$, have
to be stored since they are required for the backward pass. The loss
calculation in Algorithm 2 is an exemplary L2 loss. Note, that the loss
calculation is explicitly mentioned as
$\mathrm{Loss}(\mathrm{row},\mathrm{col})$ for illustration purposes but not
stored as the derivative of the loss is used for the back-propagation step
(see step i) in the following paragraph).

The backward pass consists of three steps: i) propagating the gradient
backward to the next layer, ii) calculate the gradient of the loss with
respect to the weighted inputs, and iii) calculate the gradient of the loss
with respect to the weights. As follows from the chain rule for back-
propagation, i) in each layer $i$ the gradient propagation is realized as
matrix multiplications between the matrix $D_{i+1}$ and the transpose of the
weights matrix $B_{i}$, ii) gradients with respect to weighted inputs are
calculated as activation of the previous layer, and iii) gradients with
respect to the weights are a matrix multiplication between the transpose of
the forward outputs $A_{i}^{T}$ and the activated output of the backward pass.
For a derivation of the backward pass, please refer to Chapter 6.5 in [93].

Our implementation of the training computes the forward pass, the loss and the
backward pass within the same kernel to maximize the re-use of cached data and
registers. Only the product $A_{i-1}^{T}D_{i}$ is not part of our training
kernel, although the time required for it is included in all of the
performance data presented in the following Section IV. We found that it is in
general preferable to launch a separate kernel for the products
$A_{i-1}^{T}D_{i}$ due to the differing dimensions of the matrix
multiplications. In detail, a possible version which fuses the matrix
multiplications $A_{i-1}^{T}D_{i}$ into the training and which reuses the data
already in the GRF would perform a block outer product of a block row in the
matrix $A_{i}$ with a block row at the same position in the matrix $D_{i}$ to
compute a partial sum of the whole $K\times N$ output matrix. A subsequent
reduction over all work-items would then be necessary for the final result. At
the time of writing, we were not able to find an efficient implementation of
the above approach or any other approach to the fusing of the product
$A_{i-1}^{T}D_{i}$. We instead found that performing the operations
$A_{i-1}^{T}D_{i}$ with the highly optimized matrix-multiplication routines in
Intel’s oneMKL, in a separate kernel after our training implementation
finishes, delivers the best results.

Algorithm 2 Training

$\mathrm{nlayers}$, $\sigma_{f}$, $\sigma_{b}$, $\mathrm{Input}$,
$\mathrm{Target}$, $B_{1},\ldots,B_{\mathrm{nlayers}-1}$

Initialize: $A_{1}=\mathrm{Input}$

for $i\leftarrow 1$ to $\mathrm{nlayers}-1$ do $\triangleright$ Forward Pass

$C_{i}\leftarrow A_{i}B_{i}$

$A_{i+1}\leftarrow\sigma_{f}(C_{i})$

end for

for $\mathrm{row}\leftarrow 1,\ldots,M$ do $\triangleright$ Loss Calculation

for $\mathrm{col}\leftarrow 1,\ldots,K$ do

$\mathrm{Diff}(\mathrm{row},\mathrm{col})=A_{\mathrm{nlayers}}(\mathrm{row},\mathrm{col})-\mathrm{Target}(\mathrm{row},\mathrm{col})$

$\mathrm{Loss}(\mathrm{row},\mathrm{col})=\mathrm{Diff}(\mathrm{row},\mathrm{col})^{2}/(MK)$

$D_{\mathrm{nlayers}}(\mathrm{row},\mathrm{col})=2\mathrm{Diff}(\mathrm{row},\mathrm{col})/(MK)$

end for

end for

$D_{nlayers}\leftarrow\sigma_{b}(A_{nlayers},D_{nlayers})$

$G_{nlayers-1}\leftarrow A_{nlayers-1}^{T}D_{nlayers}$

for $i\leftarrow\mathrm{nlayers-1}$ to $2$ do $\triangleright$ Backward Pass

$D_{i}\leftarrow D_{i+1}B_{i}^{T}$

$D_{i}\leftarrow\sigma_{b}(A_{i},D_{i})$

$G_{i-1}\leftarrow A_{i-1}^{T}D_{i}$

end for

Return $G_{1},\ldots,G_{\mathrm{nlayers-1}}$

###  III-B SYCL joint_matrix Implementation of MLPs

Our SYCL implementation is based on Intel’s joint_matrix extension [9] to
write code which can be executed on Intel’s XMX hardware [94] for the matrix
multiplications. Efficient utilization of the XMX hardware is highly relevant
since the theoretical multiply-add (MAD) peak performance of an Intel Data
Center GPU Max 1550 on the targeted bfloat16 [95] (bf16) data type is
approximately 838 tera floating point operations per second (Tflops/s), which
is sixteen times the theoretical peak single-precision MAD throughput ($\sim
52$ Tflops/s) when utilizing the vector engines [94].

The Intel SYCL extension provides a joint_matrix object, which represents a
matrix of a small and fixed size distributed across a SYCL sub-group [9], and
several related functions to facilitate computations with these joint_matrix
objects. In particular, we utilize the function joint_matrix_mad, which takes
three joint matrices, say $A$, $B$, and $C$, and returns the joint_matrix
resulting from the operation $AB+C$, which is performed on the XMX hardware.

The joint_matrix_mad function is only available for specific matrix sizes. In
detail, the matrix $A$ has to be of size $\mathrm{TM}\times\mathrm{TK}$ where
$\mathrm{TM}\in\\{1,\ldots,8\\}$ can be chosen arbitrarily and $\mathrm{TK}$
depends on the device and the data type [96]. For the targeted bfloat16 data
type, $\mathrm{TK}=16$ is required. The matrix $B$, in turn, has to be of size
$\mathrm{TK}\times\mathrm{TN}$, where $\mathrm{TN}$ is device dependent. For
the Intel Data Center GPU Max the value $\mathrm{TN}=16$ is required. The
matrix $C$ as well as the output matrix of the joint_matrix_mad function have
to be of size $\mathrm{TM}\times\mathrm{TN}$. It is important to note that
joint_matrix_mad performs the accumulation in single precision when utilizing
bfloat16 for the inputs $A$ and $B$. The matrix $C$ and the output of the
function are therefore of type float.

In the following detailed description, we focus on the inference (cf. Alg. 1)
for the special case where $K=N=64$. We launch the kernel with
$M/\mathrm{TM}\times 16$ work-items, a sub-group size of $16$ and a work-group
size of $1024$ work-items, i.e., the maximum possible work-group size to
minimize the loads from HBM as discussed in the following Section III-C. Each
sub-group with id $\mathrm{sg\\_id}$,
$\mathrm{sg\\_id}=1,\ldots,M/\mathrm{TM}$ loads a unique block-row from the
input into the register file starting at the row
$\mathrm{TM}\times\mathrm{sg\\_id}$ consisting of $\mathrm{TM}$ rows and $K$
columns. Each sub-group stores this block-row as four joint_matrix objects,
each of size $\mathrm{TM}\times\mathrm{TK}$. In each layer $i$, each work-
group (consisting of 64 sub-groups) loads the weight matrix $B_{i}$ jointly
from HBM into SLM. At the end of the algorithm, each sub-group stores its
$\mathrm{TM}$ rows of the matrix $A_{\mathrm{nlayers}}$ to HBM.

Figure 2 illustrates the loads to the register file (green) per sub-group and
the joint load to SLM (blue) per work-group for a single layer $i$. Each sub-
group then computes the product of its block-row with the weights matrix
$B_{i}$ and keeps the resulting four joint_matrix objects ($K=64$,
$\mathrm{TK}=16$) in the registers. The weight sub-matrices sized
$\mathrm{TK}\times\mathrm{TN}$, which are required to facilitate the
joint_matrix_mad function, are loaded on-the-fly from SLM as needed using the
joint_matrix_load [97] function. The weights are stored in a packed format
[98] to increase the performance of the load.

After the matrix-matrix product, each sub-group applies the activation
function $\sigma$ to its resulting block-row and utilizes the output of the
activation function in the next iteration as input to the subsequent matrix-
matrix product. Thus, it keeps the data in the register file and avoids
accesses to HBM or caches. This idea is similar to [6]. However, our approach
keeps the weights matrix in SLM and the input $A_{i}$ and output $C_{i}$ in
the registers. Our algorithm requires two local synchronizations per layer due
to the utilization of the SLM for the weights. The first synchronization
ensures that the copy of the weight matrix into SLM is finished before
accessing it. The second synchronization is to ensure that every work-item in
the work-group finished accessing the weight matrix in SLM before starting the
copy the weight matrix for the next layer.

The training algorithm 2 is similar to the inference described above and
follows the same structure. Although, in contrast to the inference, each layer
$i$ in the training has to store the matrices $A_{i}$ and $D_{i}$, which may
impact the performance. As indicated above, the final matrix-matrix
multiplications, $A_{i}^{T}D_{i}$, are performed by oneMKL outside of our
training kernel since no fusing is possible for these operations.

Figure 2: Illustration of our implementation. On the left hand-side, all
multiple layers are sketched. Each layer is parallelized along the batch size
(i.e., the $M$ dimension). The right hand-side sketches a single layer. The
data of a single sub-group is colored. Green indicates data in the register
file. Blue indicates data in the SLM. Each sub-group performs a matrix-matrix
multiplication of size
$(\mathrm{TM}\times\mathrm{K})\times(\mathrm{K}\times\mathrm{N})$ utilizing
joint_matrix objects of size $\mathrm{TM}\times\mathrm{TK}$,
$\mathrm{TK}\times\mathrm{TN}$, and $\mathrm{TM}\times\mathrm{TN}$.

To conclude this section, we discuss the limitations of our approach and
propose mitigation strategies.

First, to achieve an occupancy of 100% on the Intel Data Center GPU Max 1550,
8192 sub-groups have to be launched. Thus, problems with a batch size of less
than $8192\times 8=65536=2^{16}$ ($\mathrm{TM}=8$) do not fully occupy the
device with our strategy and show reduced performance. This restriction is
alleviated by choosing a smaller value for $\mathrm{TM}$, which increases the
number of launched sub-groups. As discussed in Section III-C this reduces the
arithmetic intensity and may lead to a higher occupancy at the cost of a lower
arithmetic intensity.

Secondly, due to the fixed dimensions of the inputs to the joint_matrix_mad
functions, the batch size $M$ has to be a multiple of $\mathrm{TM}$ and the
network width $N=K$ has to be a multiple of $\mathrm{TN}$ and $\mathrm{TK}$.
This limits the possible width of the network to multiples of 16 on the Intel
Data Center GPU Max 1550. These limitations can be removed by simply padding
the batch size to the next multiple of $\mathrm{TM}$ and the network width to
the next multiple of 16. Another alternative would be to use two-dimensional
load functions and hardware-supported zero-padding of out-of-bounds accesses
for these two-dimensional loads. Although, at the time of writing, support of
two-dimensional loads is currently lacking in Intel’s SYCL implementation and
may required to use intrinsics and built-in functionalities.

Finally, our algorithm is based on storing whole block-rows in the register
file. The maximum possible width of the network is therefore determined by the
size of the register file and given as 128 elements when utilizing the large
register mode [92]. In fact, assuming that the sizes of the data types of the
inputs, the weights, and the output of each layer are given as, say, $s_{A}$,
$s_{B}$, and $s_{C}$ bytes, respectively, we allocate $\mathrm{TM}\times
K\times s_{A}+\mathrm{TM}\times N\times
s_{C}+\mathrm{TK}\times\mathrm{TN}\times s_{B}$ bytes per sub-group in the
register file at the same time to facilitate the joint_matrix multiplications
(3584 bytes for $N=K=64$, $\mathrm{TM}=8$, $\mathrm{TK}=\mathrm{TN}=16$,
$s_{A}=s_{B}=2$, and $s_{C}=4$.) While the default size of the register file
of the Intel Data Center GPU Max 1550 is 8192 bytes per sub-group, i.e., more
than twice the size required for the joint_matrix multiplications, we
empirically found that a width larger than 64 elements does not fit in the GRF
and thus results in register spills. Although, the register allocation depends
on the driver version and may change in the future. Since the large register
mode doubles the GRF size, MLPs with a width of 128 elements fit. The register
pressure may be reduced by choosing a smaller value for $\mathrm{TM}$ at the
cost of performance (see following Section III-C for the relation between
$\mathrm{TM}$ and the arithmetic intensity).

###  III-C Roofline Analysis

Similar to [6], our goal is to minimize the accesses to the relatively slow
HBM. To analyze how well our implementation achieves this goal, we compare our
implementation to the CUDA implementation in [6] utilizing a roofline model
[99] with bandwidth limits. We compute and compare the arithmetic intensities
of the algorithms for $\mathrm{TM}=8$, $K=N=64$ and the bfloat16 data type
(i.e., 2 byte per value in memory). A higher arithmetic intensity indicates
that more floating point operations per byte can be performed and thus shows
that the performance is less limited by the memory bandwidth. Note, that a
higher arithmetic intensity does not necessarily translate to a higher
throughput since other factors (e.g. latency, frequency throttling, etc.) may
be limiting the performance. Further note, that we are not considering the
floating point operations performed by the activation function in the
following analysis.

We use 1024 work-items per work-group and a sub-group size of 16 work-items.
Thus 64 sub-groups (1024/16) constitute one work-group. The number of loads of
the weights matrix from HBM is inversely proportional to the number of sub-
groups in a work-group (see Section III-B). As each sub-group computes
$\mathrm{TM}$ rows of the $M$ rows of the input, we launch $M/\mathrm{TM}$
sub-groups or $M/\mathrm{TM}/64=M/\mathrm{TM}\times 16/1024$ work-groups. Each
work-group has its own SLM and thus each work-group needs to load the weights
matrix from HBM into its SLM. Note, that in our particular case the weights
matrix may be cached in L2. We neglect the effects of the L2 cache in our
roofline analysis, as it would depend on too many parameters, including $N$,
$K$, $M$, $\mathrm{TM}$, $\mathrm{nlayers}$, and datatypes.

As a consequence, for a single layer $i$, $i=1,\ldots,\mathrm{nlayers}-1$,
under the above mentioned assumption that none of the data is cached in the L2
cache, our algorithm loads $M/\mathrm{TM}\times 16/1024\times 64\times
64\times 2$ bytes from HBM to facilitate a matrix-matrix product of $2\times
64\times 64\times M$ flops for an arithmetic intensity of
$64\times\mathrm{TM}$ flops per byte loaded from HBM. For our choice of
$\mathrm{TM}=8$ our algorithm thus achieves an arithmetic intensity of 512
flops per byte loaded from HBM for the matrix-matrix product in every hidden
layer, i.e., all layers except the input and output layers (where additional
loads of the input from HBM or stores of the output to HBM are required).
Together with a HBM bandwidth of approximately 2 TB/s (theoretical peak is
$\sim 3.2$ TB/s, 2 TB/s is a realistic bandwidth for complex workloads), this
shows that the arithmetic intensity is high enough such that the available HBM
bandwidth is not a limiting factor and that each layer is compute bound.

Analogously, the tiny-cuda-nn implementation presented in [6] loads the weight
matrix for each layer $i$, $i=1,\ldots,\mathrm{nlayers}-1$, $M/128$ times from
HBM (again discarding that the values may be cached in the L2 cache). This
results in an arithmetic intensity of 128 flops per byte loaded from HBM. The
Nvidia H100 PCIe GPU delivers a theoretical peak HBM bandwidth of
approximately 2 TB/s [13]. The theoretical peak throughput of the tiny-cuda-nn
inference algorithm is thus reduced to 256 Tflops/s. Each layer in the CUDA
implementation is therefore bound by the HBM bandwidth.

Extending the above model to consider the whole inference algorithm, including
the loads of the input data and stores of the output, we get the following
arithmetic intensities (also known as operational intensities; abbreviated as
OI to prevent confusion with AI – Artificial Intelligence) in dependence of
the number of layers:

$\displaystyle\mathrm{OI}_{\text{SYCL}}=\frac{512(\mathrm{nlayers}-1)}{16+(\mathrm{nlayers}-1)}$
,
$\displaystyle\mathrm{OI}_{\text{CUDA}}=\frac{128(\mathrm{nlayers}-1)}{4+(\mathrm{nlayers}-1)}.$

For example, for a number of $\mathrm{nlayers}=6$, the arithmetic intensity of
the SYCL implementation is 121.9 flops per byte, limiting the theoretical peak
performance to approximately 243.8 Tflops/s. The arithmetic intensity for the
CUDA implementation is 71.1 flops per byte, resulting in a theoretical peak
performance of approximately 142.2 Tflops/s.

Further extending the above model to the forward pass in the training step,
where each $A_{i}$ is additionally stored in HBM, we get an arithmetic
intensity of $512(\mathrm{nlayers}-1)/(9\mathrm{nlayers}+7)$ flops per byte
for the SYCL implementation ($\sim 49$ flops per byte for
$\mathrm{nlayers}=6$; $\sim 98$ Tflops/s) and
$128(\mathrm{nlayers}-1)/(3\mathrm{nlayers}+1)$ for the CUDA implementation
($\sim 33.7$ flops per byte for $\mathrm{nlayers}=6$; $\sim 67.4$ Tflops/s)
thus further reducing the theoretical peak performance for both codes.

Finally, considering the forward pass and the backward pass (loading the
input, storing each $A_{i}$ and $D_{i}$ once, loading $B_{i}$ and $B_{i}^{T}$
according to the work-group size and the batch size, loading $A_{i}$ and
$D_{i}$ once for the final calculation of the $G_{i}$) results in theoretical
peak performance as displayed in Fig. 3.

Similar considerations as above for the accesses to SLM lead to an arithmetic
intensity of 7.877 flops per byte accessed in SLM for the SYCL implementation
and 25.6 flops per byte for the CUDA code. Thus, based on the theoretical peak
SLM bandwidth of the targeted device, the SLM accesses limit the performance
of our SYCL code to 393 Tflops/s (instead of 838 Tflops/s). This is never an
issue for the training since the HBM bandwidth imposes a lower limit on the
performance than the SLM bandwidth. Although, the SLM bandwidth limits the
performance for the inference to at most to 393 Tflops/s when
$\mathrm{nlayers}\geq 12$.

To emphasize the advantages of fused MLPs compared to non-fused MLPs, we
compare the above arithmetic intensities to the arithmetic intensities of non-
fused implementations in what follows. A non-fused implementation of the
inference algorithm 1 would require for each layer a load of the matrices
$A_{i}$ and $B_{i}$ with a subsequent store of $A_{i+1}$ to facilitate a
single matrix-matrix product. Thus, the arithmetic intensity is bounded by 32
flops per byte loaded (independent of the number of layers) for the inference,
which is up to 1/16-th of the arithmetic intensity of a fused implementation.
Similar considerations for the training show that the fused training only
increases the arithmetic intensity by up to a factor 2 at most (in the limit
$\mathrm{nlayers}\to\infty$).

To summarize, Figure 3 shows the theoretical peak performance based on the
above roofline analysis for the CUDA implementation on Nvidia’s H100 GPU and
the SYCL implementation on Intel’s Data Center GPU Max 1550. It shows that our
algorithm improves the theoretical peak performance of the CUDA implementation
[6] by increasing the arithmetic intensity for the HBM accesses and that
decreasing the arithmetic intensity for SLM accesses does not have a negative
impact. The performance which may be achieved in practice is discussed in the
following Sec. IV.

Figure 3: Comparison of the theoretical peak performance based on the the
roofline analysis in Section III-C.

##  IV Results

In this section, we demonstrate how our fully-fused MLP implementation
increases the performance of four commonly used AI tasks: non-linear function
approximation, image compression, Neural Radiance Fields (NeRFs), and solving
differential equations with PINNs (see Section II). For every comparison, we
compare our SYCL implementation on an Intel Data Center GPU Max 1550 with the
CUDA implementation [8] on a Nvidia H100 GPU and PyTorch using both Intel
Extension for PyTorch (IPEX) [12] and CUDA backend. The results of the
comparison can be seen in Table I.

TABLE I: Training and Inference times for our implementation (SYCL), the fully-fused implementation from [6] (CUDA) and PyTorch using both IPEX and CUDA. The numbers in the brackets indicate the time relative to our implementation. |  | Training [s] $\downarrow$ | Inference [s] $\downarrow$
---|---|---|---
Benchmark | SYCL | 0.58 | 0.083
CUDA | 0.79 (1.36) | 0.236 (2.84)
PyTorch (IPEX) | 4.57 (7.88) | 2.178 (26.24)
PyTorch (CUDA) | 4.32 (7.45) | 1.543 (18.55)
Image compr. | SYCL | 9.20 | 2.041
CUDA | 16.10 (1.75) | 3.514 (1.72)
PyTorch (IPEX) | 41.75 (4.54) | 12.84 (6.29)
PyTorch (CUDA) | 36.66 (3.98) | 15.46 (7.57)
NeRF | SYCL | 1.93 | 0.302
CUDA | 2.05 (1.06) | 0.477 (1.58)
PyTorch (IPEX) | 15.09 (7.82) | 9.050 (29.97)
PyTorch (CUDA) | 10.07 (5.21) | 3.841 (12.72)
PINNs | SYCL | 0.55 | 0.088
CUDA | 0.60 (1.09) | 0.108 (1.22)
PyTorch (IPEX) | 1.92 (3.49) | 0.788 (8.95)
PyTorch (CUDA) | 1.85 (3.36) | 0.646 (7.34)

To ensure a fair comparison of the forward and backward performances and to
disregard differences in the implementation of the loss functions and
optimiser, we only measure the time of the code where the fusion of the
operators is relevant for, i.e., forward and backward calculation (cf. Algs.
1, 2). The benchmark protocol is identical for both implementations: every
forward and backward pass is called for the same amount of episodes, the same
batch sizes, MLP architecture, and ReLU activation functions and a linear
output activation function.

The CUDA implementation [8] is tested on a dual-socket system with two Intel
Xeon Platinum 8480+ CPUs, 512GB DRAM and Nvidia H100 GPUs with Ubuntu 22.04
LTS as operating system. The code was compiled with CUDA version 12.3.52 (nvcc
version built on Sept. 8 19:17:24PDT 2023). The driver version is 535.129.03.

The SYCL code, on the other hand, is tested on a dual-socket system with two
Intel Xeon Platinum 8480+ CPUs, 512GB DRAM and Intel Data Center GPU Max 1550
GPUs with Ubuntu 22.04.3 LTS as operating system. We utilized the mpiicpx
compiler (Intel oneAPI DPC++/C++ Compiler 2023.2.0 (2023.2.0.20230721))
included in the oneAPI 2023.2 toolkit and an unreleased engineering
driver***agama-ci-devel 682.16. To scale our implementation to both Xe Stacks
[94] in an Intel GPU, we use Intel MPI version 2021.10.0.

###  IV-A Non-linear Function Approximation

The goal of non-linear function approximation is to train a Neural Network
$y=\hat{f}_{\theta}(x)$ to approximate a given non-linear function
$f:\mathbb{R}^{K}\rightarrow\mathbb{R}^{N}$. During the training the residual
loss between target function $f$ and approximation $\hat{f}_{\theta}$ is
typically chosen as the mean square error, which is depicted in what follows
for a batch size $M$:

$L(f,\hat{f}_{\theta},x_{1},\ldots,x_{M})=\sum_{i=1}^{M}\|f(x_{i})-\hat{f}_{\theta}(x_{i})\|^{2}.$
(1)

In this section, the performance of our SYCL implementation is compared to the
CUDA implementation [8] for various batch sizes $M$. To determine relevant
batch sizes, we investigated the distribution of batch sizes used in MLP-based
methods [100, 101, 51, 52, 50, 102, 103, 104, 104, 49, 105, 106, 46, 107, 108,
109, 45, 110, 111, 4, 112, 58, 113, 114, 115, 116, 117, 118, 119, 119, 120,
121, 53, 42, 43, 47, 122, 123, 124, 125, 126, 127] in the fields of NeRFs,
Neural Compression, and PINNs mentioned in Appendix -A\- -C. The results of
this investigation are summarized in Table II. The distribution appears to be
approximately normal, with the majority of batch sizes falling between
$2^{16}$ and $2^{22}$ suggesting that many MLP-based methods in these fields
use large batch sizes, where our fully-fused implementation reaches maximal
occupancy and performance (see Fig. 4).

TABLE II: Batch Sizes of commonly used MLP applications Batch Size Interval | $[2^{5},2^{11})$ | [$2^{11},2^{15}]$ | $[2^{16},2^{22}]$ | $(2^{22},2^{24}]$
---|---|---|---|---
# Occurences | 6 | 17 | 22 | 1

For NeRFs, the batch size is the number of rays per batch and corresponding
samples per ray during ray tracing. We the formula in [128] to calculate the
batch size: number of rays per batch $\times$ number of samples per ray for
the full method. For Neural Compression, and Partial Differential Equations
(PDEs), since many methods in these representation learning use full-batch
gradient descent as the training strategy, we adopt the dataset size as the
full batch size for both the training and inference steps.

Based on this data, batch sizes larger than $2^{22}$ elements are not commonly
used and thus not considered in this test. Small batch sizes, in turn, are
also removed from the test since the occupancy decreases proportional to the
batch size for batch sizes smaller than $2^{16}$ (cf. Sec. III-B). The
performance is thus limited by the occupancy and fusing the operations does
not provide any benefit.

In our test we learn an MLP that approximates a non-linear function
$f:\mathbb{R}^{64}\rightarrow\mathbb{R}^{64}$. We measure the throughput for
our SYCL implementation on a single Xe Stack and both Xe Stacks of an Intel
Data Center GPU Max 1550 and compare the results to the throughput measured
with the CUDA implementation [8] on a Nvidia H100 GPU. For both
implementations, we use a network width of $64$ with input and output width of
$64$ elements. As indicated above, the batch size $M$ varies from $2^{11}$ to
$2^{22}$, and we run the benchmark for $\mathrm{niter}$ iterations calculated
as $\mathrm{niter}=\text{max}(1000\cdot\frac{2^{18}}{M},250)$. Further, we
choose four hidden layers, i.e., $\mathrm{nlayers}=6$.

Figure 4: Performance of the inference (top) and training (bottom) of our SYCL
implementation on a single tile of the Intel Data Center GPU Max 1550 (dark
blue) and two tiles of the same Intel GPU (light blue) compared to the CUDA
code on a Nvidia H100 GPU (green) for $\mathrm{nlayers}=6$. The occupancy of
the SYCL code on the two tiles of the Intel Data Center GPU is given as the
blue line. The x-axis denotes the batch size (i.e., $M$) from $2^{11}$ inputs
to $2^{22}$.

The resulting performance graphs for training and inference can be seen in
Figure 4 and the relative performance gain of our method compared to CUDA can
be seen in Table III. The figure shows that the performance in the training
case (bottom graph of Fig. 4) is similar for both implementations. This
coincides well with the roofline predictions from Section III-C. It is also
interesting to note that the peak performance of approximately 70 Tflops/s is
close to the roofline prediction of $\sim 87$ Tflops/s thus indicating that
the problem is indeed bound by the HBM bandwidth.

TABLE III: Relative Performance Gain of SYCL compared to CUDA for Function Approximation | Mean | Min | Max
---|---|---|---
Training throughput | 1.01 | 0.74 | 1.37
Inference throughput | 1.34 | 1.03 | 1.68

The inference results (top graph of Fig. 4), in contrast, show that the SYCL
implementation outperforms the CUDA implementation. For the smaller batch
sizes even when utilizing only half of the Intel device. This coincides again
well with the roofline analysis in Section III-C. Although, the measured peak
performance of approximately 140 Tflop/s falls short of the predicted $\sim
244$ Tflops/s. Our performance investigations have shown that this is due to
scoreboard ID stalls occurring when loading the weights matrices from SLM thus
indicating that the performance of our problem is not bound by the memory
bandwidth but by the SLM memory latency. Investigation how to eliminate this
latency are ongoing.

To fully show demonstrate the performance of our approach, we measure the time
required to calculate 1000 iterations for non-linear function approximation
using an MLP of 11 hidden layers and network width of 64. The batch size is
chosen as $M=2^{17}$. The results can be seen in the Benchmark row of Table I.
Our SYCL implementation reaches performance gains of up to 7.88 times faster,
and up to 26.24 times faster for training and inference respectively.

As mentioned in Section III-B, the Intel GPU is not fully occupied for batch
sizes smaller than $2^{16}$ resulting in reduced performance. Note that
maximum performances for the SYCL implementation are not achieved for the
largest batch sizes, but for those batch sizes where the occupancy is maximal
while the problem is sufficiently small to fit in the L2 cache.

###  IV-B Image Compression

For Image Compression, the MLP $\hat{f}_{\theta}$ learns the color
distribution of an image represented as function
$f:\mathbb{R}^{2}\rightarrow\mathbb{R}^{1}$ that maps from 2-D coordinates
(pixels) to corresponding color. We use a network width of $64$. As NNs poorly
learn high-frequent details of the target function $f$, encodings are
typically used to mitigate this issue. Encodings
$g:\mathbb{R}^{2}\rightarrow\mathbb{R}^{N}$ pre-process the 2D pixel input by
directly mapping them into a spectrum of high-frequent signals, which are then
fed into the MLP. The Multiresolution hash encoding
$g:\mathbb{R}^{2}\rightarrow\mathbb{R}^{32}$ [128] is used, which can be seen
as the input to the MLP. The final function approximator is thus
$h(x)=\hat{f}(g(x))$, with Multiresolution Hash Encoding $g$ and MLP
$\hat{f}:\mathbb{R}^{32}\rightarrow\mathbb{R}^{1}$.

The learning progress can be seen in Figure 5. We used a network width of $64$
with input $32$, output width $1$, and batch size $K=2304\times 3072=7077888$
for 1000 iterations. After 10 steps the colors are not entirely correct, but
the silhouette is recognisable. Then, over the next 900 steps, the colors and
finer details are becoming gradually clearer until the training progress
converges after 1000 steps. The whole training takes $9.2$s for the SYCL
implementation and is 1.75 faster than the CUDA implementation and 4 times
faster than PyTorch. The image can be reconstructed (inference) within $2.41$s
for SYCL (see row Image compr. in Table I).

For comparison, storing the image at a resolution of $2304\times 3072$ with
greyscale values costs 2.3MB for a JPEG image, while storing 12288 half
precision weights costs $25$ kb.

Figure 5: Training progress of Image Compression. The training converges after
1000 steps. The training and inference process is performed for both the SYCL
and CUDA implementation and the visualised progress if implemented as per [8].

###  IV-C Neural Radiance Fields (NeRF)

For NeRFs, the goal is to learn a Radiance Field (see Appendix -A)
$f:\mathbb{R}^{5}\rightarrow\mathbb{R}^{4}$. Similar to Image Compression,
NeRFs require an encoding as input to the MLP for better representations of
the high-frequent details of the object. Multiresolution hash encoding is used
for this task resulting in a 32-dimensional feature vector as input to the
MLP. A network width of $64$ is chosen with input width $32$, output width
$4$, and a batch size $K=1048576$. For details about the implementation of the
NeRF algorithm, please see [129].

The training of the MLP without encoding and volume rendering is conducted
within $1.93$s for training and $0.302$s inference in our SYCL implementation.
Compared to the CUDA implementation, this is 1.06 and 1.58 times faster for
training and inference respectively. The generated images can be seen in
Figure 6.

Figure 6: CLNeRF results on 5 scenes for breville, kitchen, spa, community,
and living room. Each row represents a different scene and each column shows
the scene rendering results from different views using the MLP as NeRF.

###  IV-D Physics-Informed Neural Networks (PINN)

Figure 7: Visualisation of learned solution to Navier Stokes Equation. Mean
Absolute Percentage Error: 0.1%. Results generated using DeepXDE [130].

To demonstrate how our method can be applied to solve Partial Differential
Equations (PDE), we leverage the PINNs framework [123] to train an MLP that
represents the solution of the following two-dimensional Navier-Stokes
Equation:

$\displaystyle
u_{t}+\lambda_{1}(uu_{x}+vu_{y})+p_{x}-\lambda_{2}(u_{xx}+u_{yy})$
$\displaystyle=0,$ (2) $\displaystyle
v_{t}+\lambda_{1}(uv_{x}+vv_{y})+p_{y}-\lambda_{2}(v_{xx}+v_{yy})$
$\displaystyle=0,$ (3) $\displaystyle u_{x}+v_{y}$ $\displaystyle=0,$ (4)

with velocity fields $u(t,x,y)$ and $v(t,x,y)$ for the $x$ and $y$ components
respectively, pressure $p(t,x,y)$, and unknown parameters
$\lambda=(\lambda_{1},\lambda_{2})$. For more details of the problem, please
refer to Section 4.1.1 in [123].

We used DeepXDE [130] to find the PINNs $f(t,x,y)$, $g(t,x,y)$ that represent
the solution of the Navier stokes equations (2). For details, see Appendix -C.

The dataset and thus batch size consists of $M=2^{17}$ data points. For 1000
iterations, the training time takes $0.55$s for SYCL, $0.60$s for CUDA
yielding a performance gain of up to 3.49 compared to PyTorch. The evaluation
of the PDE solution on $2^{17}$ collocation points requires $0.088$, which
results in an improvement of up to 8.95x.

##  V Conclusion and Future Work

In this paper, we presented a SYCL implementation of fully-fused Multi-Layer
Perceptrons (MLPs) on Intel Data Center GPU Max. Our approach focuses on
maximizing data reuse within the general register file and the shared local
memory, minimizing the need for slower global memory accesses. This results in
a significant increase in the arithmetic intensity, which leads to improved
performance, especially for inference.

Our implementation outperforms an equivalent CUDA implementation for MLPs with
width 64 by a factor of up to 2.84 in inference and 1.75 in training in our
tests, demonstrating the effectiveness of our approach, and outperforms the
PyTorch implementation by up to a factor of 30. We further showcased the
efficiency of our implementation in three significant areas: Image
Compression, Neural Radiance Fields (NeRF), and Physics-Informed Machine
Learning. Across all these domains, our approach demonstrated substantial
improvements, achieving factors up to 30 times when compared to conventional
PyTorch implementations and up to 2.84 times over highly optimized CUDA
implementations.

In the future, we aim to further optimize our implementation. A strong focus
will be a more efficient usage of registers which may enable further
prefetching of the weights matrices, thus reducing the stalls. In addition, we
may be able to reduce the utilization of SLM and enable the loads of multiple
weights matrices in the SLM, which would reduce the number of necessary
barriers. Another important aspect will be to increase the occupancy for small
batch sizes by either increasing the arithmetic intensity for smaller values
of $\mathrm{TM}$ or by partitioning the work along the $K$ dimension, which
requires the usage of atomic operations. Lastly, an important optimization
would be the fusion of the final matrix products $A_{i}^{T}D_{i}$ (cf. Alg. 2)
into the backward pass. It could potentially increase the performance of the
training by up to a factor 2, based on the roofline model.

Our implementation is open-sourced to allow for wider usage and contributions
from the community. The code can be found at https://github.com/intel/tiny-
dpcpp-nn.

In addition to further performance optimization, we also plan to explore the
use of Intel’s ESIMD SYCL extension for our implementation and compare it to
our existing SYCL implementation. Our past experience with ESIMD showed that
it enables a better control over the register usage and exposes cache controls
for load and store operations.

In future work, we aim to generalize our library to various data types and
larger network widths. In addition, we plan to implement an optimized version
for Intel’s Arc GPUs.

## Acknowledgment

We’d like to thank Jia Xu and Jing Xu for their help and advise during the
development of tiny-dpcpp-nn.

## Disclaimers

Performance varies by use, configuration and other factors. Learn more on the
Performance Index site.

Performance results are based on testing as of dates shown and may not reflect
all publicly available updates. No product or component can be absolutely
secure.

Your costs and results may vary.

Intel technologies may require enabled hardware, software or service
activation.

©Intel Corporation. Intel, the Intel logo, and other Intel marks are
trademarks of Intel Corporation or its subsidiaries. Other names and brands
may be claimed as the property of others.

## References

  * [1] B. Yegnanarayana, _Artificial neural networks_. PHI Learning Pvt. Ltd., 2009.
  * [2] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [3] Z. Li, F. Liu, W. Yang, S. Peng, and J. Zhou, “A survey of convolutional neural networks: analysis, applications, and prospects,” _IEEE transactions on neural networks and learning systems_ , 2021.
  * [4] S. Cuomo, V. S. Di Cola, F. Giampaolo, G. Rozza, M. Raissi, and F. Piccialli, “Scientific machine learning through physics-informed neural networks: Where we are and what’s next,” _Journal of Scientific Computing_ , vol. 92, no. 88, pp. 1–35, 2022.
  * [5] B. Mildenhall, P. P. Srinivasan, M. Tancik, J. T. Barron, R. Ramamoorthi, and R. Ng, “Nerf: Representing scenes as neural radiance fields for view synthesis,” _Communications of the ACM_ , vol. 65, no. 1, pp. 99–106, 2021\.
  * [6] T. Müller, F. Rousselle, J. Novák, and A. Keller, “Real-time neural radiance caching for path tracing,” _arXiv preprint arXiv:2106.12372_ , 2021\.
  * [7] S. Park, C. Yun, J. Lee, and J. Shin, “Minimum width for universal approximation,” _arXiv preprint arXiv:2006.08859_ , 2020.
  * [8] T. Müller, “tiny-cuda-nn,” 4 2021. [Online]. Available: https://github.com/NVlabs/tiny-cuda-nn
  * [9] Intel Corporation. (2023) Programming Intel® XMX using SYCL: Joint Matrix Multiplication. Intel oneAPI Optimization Guide for GPU. [Online]. Available: https://www.intel.com/content/www/us/en/docs/oneapi/optimization-guide-gpu/2024-0/programming-intel-xmx-using-sycl-joint-matrix.html
  * [10] J. Peddie, “The sixth era gpus: Ray tracing and mesh shaders,” in _The History of the GPU-New Developments_. Springer, 2023, pp. 323–360.
  * [11] H. Jiang, “Intel’s ponte vecchio gpu: Architecture, systems & software,” in _2022 IEEE Hot Chips 34 Symposium (HCS)_. IEEE Computer Society, 2022, pp. 1–29.
  * [12] Intel Corporation, “Intel extension for pytorch,” https://github.com/intel/intel-extension-for-pytorch, 2023, gitHub repository.
  * [13] Nvidia Corporation. (2023) H100 Tensor Core GPU. [Online]. Available: {https://www.nvidia.com/en-us/data-center/h100/}
  * [14] R. Liu, Y. Li, L. Tao, D. Liang, and H.-T. Zheng, “Are we ready for a new paradigm shift? a survey on visual deep mlp,” _Patterns_ , vol. 3, no. 7, 2022.
  * [15] W. H. Delashmit, M. T. Manry _et al._ , “Recent developments in multilayer perceptron neural networks,” in _Proceedings of the seventh annual memphis area engineering and science conference, MAESC_ , 2005, pp. 1–15.
  * [16] M.-H. Guo, Z.-N. Liu, T.-J. Mu, D. Liang, R. R. Martin, and S.-M. Hu, “Can attention enable mlps to catch up with cnns?” _Computational Visual Media_ , vol. 7, pp. 283–288, 2021.
  * [17] Q. Hou, Z. Jiang, L. Yuan, M.-M. Cheng, S. Yan, and J. Feng, “Vision permutator: A permutable mlp-like architecture for visual recognition,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 45, no. 1, pp. 1328–1334, 2021.
  * [18] X. Ding, C. Xia, X. Zhang, X. Chu, J. Han, and G. Ding, “Repmlp: Re-parameterizing convolutions into fully-connected layers for image recognition,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022.
  * [19] Y. Rao, W. Zhao, Z. Zhu, J. Lu, and J. Zhou, “Global filter networks for image classification,” in _Advances in Neural Information Processing Systems_ , vol. 34, 2021.
  * [20] H. Touvron, P. Bojanowski, M. Caron, M. Cord, A. El-Nouby, E. Grave, G. Izacard, A. Joulin, G. Synnaeve, J. Verbeek _et al._ , “Resmlp: Feedforward networks for image classification with data-efficient training,” _arXiv preprint arXiv:2105.03404_ , 2021.
  * [21] H. Liu, Z. Dai, D. R. So, and Q. V. Le, “Pay attention to mlps,” _arXiv preprint arXiv:2105.08050_ , 2021.
  * [22] S. Chen, E. Xie, C. Ge, R. Chen, D. Liang, and P. Luo, “Cyclemlp: A mlp-like architecture for dense prediction,” in _Proceedings of the International Conference on Learning Representations_ , 2022.
  * [23] J. Lahoud and B. Ghanem, “2d-driven 3d object detection in rgb-d images,” in _Proceedings of the IEEE International Conference on Computer Vision_ , 2017, pp. 4622–4630.
  * [24] S. Pan, C.-W. Chang, T. Wang, J. Wynne, M. Hu, Y. Lei, T. Liu, P. Patel, J. Roper, and X. Yang, “Unext: Mlp-based rapid medical image segmentation network,” in _Medical Image Computing and Computer Assisted Intervention – MICCAI 2022_. Springer, 2022, pp. 23–33.
  * [25] L. An, L. Wang, and Y. Li, “Abdomen ct multi-organ segmentation using token-based mlp-mixer,” _Medical Physics_ , 2021.
  * [26] H.-P. Lai, T.-T. Tran, and V.-T. Pham, “Axial attention mlp-mixer: A new architecture for image segmentation,” in _2022 IEEE Ninth International Conference on Communications and Electronics (ICCE)_. IEEE, 2022, pp. 381–386.
  * [27] Z. Qiu, T. Yao, C.-W. Ngo, and T. Mei, “Mlp-3d: A mlp-like 3d architecture with grouped time mixing,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 3062–3072.
  * [28] L. An, L. Wang, and Y. Li, “Morphmlp: An efficient mlp-like backbone for spatial-temporal representation learning,” _Sensors_ , vol. 22, no. 18, p. 7024, 2021.
  * [29] J. Xia, M. Zhuge, T. Geng, S. Fan, Y. Wei, Z. He, and F. Zheng, “Skating-mixer: Long-term sport audio-visual modeling with mlps,” _arXiv preprint arXiv:2203.03990_ , 2022.
  * [30] Z. Tu, H. Talebi, H. Zhang, F. Yang, P. Milanfar, A. Bovik, and Y. Li, “Maxim: Multi-axis mlp for image processing,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 5769–5780.
  * [31] N. Sharma, R. Sharma, and N. Jindal, “New approaches to determine age and gender in image processing techniques using multilayer perceptron neural network,” _International Journal of Computer Applications_ , vol. 164, no. 1, pp. 1–5, 2017.
  * [32] G. Cazenavette and M. Ladron De Guevara, “Mixergan: An mlp-based architecture for unpaired image-to-image translation,” _arXiv preprint arXiv:2105.14110_ , 2021.
  * [33] Y. Mansour, K. Lin, and R. Heckel, “Image-to-image mlp-mixer for image reconstruction,” _arXiv preprint arXiv:2202.02018_ , 2022.
  * [34] Z. Al-Makhadmeh and A. Tolba, “Improving sentiment analysis in arabic and english languages by using multi-layer perceptron model (mlp),” in _2020 IEEE 7th International Conference on Data Science and Advanced Analytics (DSAA)_. IEEE, 2020, pp. 1328–1334.
  * [35] J. Sultana, N. Sultana, K. Yadav, and F. AlFayez, “Prediction of sentiment analysis on educational data based on deep learning approach,” in _2018 21st Saudi Computer Society National Computer Conference (NCC)_. IEEE, 2018, pp. 1–6.
  * [36] Z. Al-Makhadmeh and A. Tolba, “Automatic hate speech detection using killer natural language processing optimizing ensemble deep learning approach,” _Computing_ , vol. 102, pp. 501–522, 2019.
  * [37] T. K. Tran and H. Sato, “Nlp-based approaches for malware classification from api sequences,” in _2017 21st Asia Pacific Symposium on Intelligent and Evolutionary Systems (IES)_. IEEE, 2017, pp. 101–105.
  * [38] D. Dang, F. Di Troia, and M. Stamp, “Malware classification using long short-term memory models,” _arXiv preprint arXiv:2103.02746_ , 2021.
  * [39] Y. Qiao, W. Zhang, X. Du, and M. Guizani, “Malware classification based on multilayer perception and word2vec for iot security,” _ACM Transactions on Internet Technology (TOIT)_ , vol. 22, no. 1, pp. 1–22, 2021.
  * [40] M. Raman, P. Maini, J. Z. Kolter, Z. C. Lipton, and D. Pruthi, “Model-tuning via prompts makes nlp models adversarially robust,” _arXiv preprint arXiv:2303.07320_ , 2023.
  * [41] F. Fusco, D. Pascual, and P. Staar, “pnlp-mixer: an efficient all-mlp architecture for language,” _arXiv preprint arXiv:2202.04350_ , 2022.
  * [42] T. Alkhalifah and X. Huang, “Physics informed neural learning of wavefields using gabor basis functions.”
  * [43] M. Takamoto, F. Alesiani, and M. Niepert, “Learning neural pde solvers with parameter-guided channel attention,” in _Proceedings of the 38th International Conference on Machine Learning_ , 2023, pp. 9801–9811.
  * [44] H. Eivazi, M. Tahani, P. Schlatter, and R. Vinuesa, “Physics-informed neural networks for solving reynolds-averaged navier–stokes equations,” _Physics of Fluids_ , vol. 34, no. 7, p. 075117, 2022.
  * [45] S. Wang, S. Sankaran, and P. Perdikaris, “Respecting causality is all you need for training physics-informed neural networks,” _Journal of Scientific Computing_ , vol. 92, no. 88, pp. 1–35, 2022.
  * [46] L. Lu, P. Jin, and G. E. Karniadakis, “Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators,” _Journal of Computational Physics_ , vol. 404, p. 109108, 2020.
  * [47] L. He, Y. Chen, Z. Shen, Y. Yang, and Z. Lin, “Neural epdos: Spatially adaptive equivariant partial differential operator based networks,” in _Proceedings of the 38th International Conference on Machine Learning_ , 2023, pp. 9801–9811.
  * [48] G. Zhang, B. Patuwo, and M. Y. Hu, “Multilayer perceptrons and radial basis function neural network methods for the solution of differential equations: a survey,” _Neural Networks_ , vol. 15, no. 2, pp. 241–258, 2002.
  * [49] S. J. Garbin, M. Kowalski, M. Johnson, J. Shotton, and J. P. C. Valentin, “Fastnerf: High-fidelity neural rendering at 200fps,” _arXiv preprint arXiv:2103.10380_ , 2021.
  * [50] J. T. Barron, B. Mildenhall, M. Tancik, P. Hedman, R. Martin-Brualla, and P. P. Srinivasan, “Mip-nerf: A multiscale representation for anti-aliasing neural radiance fields,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2021, pp. 10 456–10 465.
  * [51] C. Reiser, S. Peng, Y. Liao, and A. Geiger, “Kilonerf: Speeding up neural radiance fields with thousands of tiny mlps,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2021, pp. 10 466–10 475.
  * [52] C. Sun, M. Sun, and H.-T. Chen, “Direct voxel grid optimization: Super-fast convergence for radiance fields reconstruction,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 10 304–10 313.
  * [53] K. Park, U. Sinha, J. T. Barron, S. Bouaziz, D. B. Goldman, S. M. Seitz, and R. Martin-Brualla, “Nerfies: Deformable neural radiance fields,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2021, pp. 5865–5874.
  * [54] Y. Cao, G. Chen, K. Han, W. Yang, and K.-Y. K. Wong, “Jiff: Jointly-aligned implicit face function for high quality single view clothed human reconstruction–supplementary material–,” 2022.
  * [55] Y. Zhu, X. Xiao, W. Wu, and Y. Guo, “3d reconstruction of deformable linear objects based on cylindrical fitting,” _Signal, Image and Video Processing_ , vol. 17, pp. 2617–2625, 2023.
  * [56] A. Božič, P. Palafox, M. Zollhöfer, J. Thies, A. Dai, and M. Nießner, “Neural deformation graphs for globally-consistent non-rigid reconstruction,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2021, pp. 1040–1049.
  * [57] C. Peyrard, F. Mamalet, and C. Garcia, “A comparison between multi-layer perceptrons and convolutional neural networks for text image super-resolution,” in _2014 22nd International Conference on Pattern Recognition_. IEEE, 2014, pp. 4350–4355.
  * [58] S. Lee, J. Kim, and S. Lee, “Estimating gaze depth using multi-layer perceptron,” in _2017 International Symposium on Ubiquitous Virtual Reality (ISUVR)_. IEEE, 2017, pp. 1–4.
  * [59] P. Chopade and P. Kulkarni, “Single image super-resolution based on modified interpolation method using mlp and dwt,” in _2019 3rd International Conference on Trends in Electronics and Informatics (ICOEI)_. IEEE, 2019, pp. 1–5.
  * [60] M. Juszczyk, “Application of pca-based data compression in the ann-supported conceptual cost estimation of residential buildings,” in _AIP Conference Proceedings_ , vol. 1738, no. 1. AIP Publishing LLC, 2016, p. 200007.
  * [61] Y. Strümpler, J. Postels, R. Yang, L. Van Gool, and F. Tombari, “Implicit neural representations for image compression,” _arXiv preprint arXiv:2112.04267_ , 2021.
  * [62] L. Zhao, Z. Dong, and K. Keutzer, “Analysis of quantization on mlp-based vision models,” in _Proceedings of the International Conference on Learning Representations_ , 2022.
  * [63] C. Yin, B. Acun, X. Liu, and C.-J. Wu, “Tt-rec: Tensor train compression for deep learning recommendation models,” in _Proceedings of the Conference on Machine Learning and Systems_ , 2021.
  * [64] H. Liu, W. Chen, S. Li, and J. Wang, “Path-following control of underactuated ships using actor-critic reinforcement learning with mlp neural networks,” in _2016 Sixth International Conference on Information Science and Technology (ICIST)_. IEEE, 2016, pp. 1–6.
  * [65] S. M. J. Jalali, S. Ahmadian, A. Khosravi, S. Mirjalili, M. R. Mahmoudi, and S. Nahavandi, “Neuroevolution-based autonomous robot navigation: a comparative study,” _Cognitive Systems Research_ , vol. 62, pp. 35–43, 2020\.
  * [66] B. Ko, H.-J. Choi, C. Hong, J.-H. Kim, O. C. Kwon, and C. D. Yoo, “Neural network-based autonomous navigation for a homecare mobile robot,” in _2017 IEEE International Conference on Big Data and Smart Computing (BigComp)_. IEEE, 2017, pp. 403–406.
  * [67] N. H. Singh and K. Thongam, “Mobile robot navigation using mlp-bp approaches in dynamic environments,” _Arabian Journal for Science and Engineering_ , vol. 43, no. 12, pp. 8013–8028, 2018.
  * [68] Y. Zhang and A. Srinivasa, “Autonomous rl: Autonomous vehicle obstacle avoidance in a dynamic environment using mlp-sarsa reinforcement learning,” in _2019 IEEE 5th International Conference on Mechatronics System and Robots (ICMSR)_. IEEE, 2019, pp. 1–6.
  * [69] Y. Song and W. Sun, “Pc-mlp: Model-based reinforcement learning with policy cover guided exploration,” in _Proceedings of the 38th International Conference on Machine Learning_ , 2021, pp. 9801–9811.
  * [70] M. Srouji, J. Zhang, and R. Salakhutdinov, “Structured control nets for deep reinforcement learning,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 4742–4751.
  * [71] X. Chen, L. Yao, J. McAuley, G. Zhou, and X. Wang, “Deep reinforcement learning in recommender systems: A survey and new perspectives,” _Knowledge-Based Systems_ , vol. 264, p. 110335, 2023.
  * [72] P. Zhao, K. Xiao, Y. Zhang, K. Bian, and W. Yan, “Ameir: Automatic behavior modeling, interaction exploration and mlp investigation in the recommender system,” _arXiv preprint arXiv:2006.05933_ , 2020.
  * [73] J. Liu, X.-M. Zhang, and W. Wang, “Mlp technique based reinforcement learning control of discrete pure-feedback systems,” _Journal of the Franklin Institute_ , vol. 356, no. 7, pp. 3824–3840, 2019.
  * [74] W. Bai, T. Li, and S. Tong, “Nn reinforcement learning adaptive control for a class of nonstrict-feedback discrete-time systems,” _IEEE Transactions on Cybernetics_ , vol. 50, no. 11, pp. 4573–4584, 2020.
  * [75] J. Bjorck, C. P. Gomes, and K. Q. Weinberger, “Towards deeper deep reinforcement learning,” _arXiv preprint arXiv:2106.01151_ , 2021.
  * [76] M. Wagenaar, “Learning to play the game of hearts using reinforcement learning and a multi-layer perceptron,” Ph.D. dissertation, Faculty of Science and Engineering, 2017.
  * [77] A. Carosia, G. P. Coelho, and A. Silva, “Analyzing the brazilian financial market through portuguese sentiment analysis in social media,” _Applied Artificial Intelligence_ , vol. 34, no. 1, pp. 1–19, 2020.
  * [78] S. Bairavel and M. Krishnamurthy, “Novel ogbee-based feature selection and feature-level fusion with mlp neural network for social media multimodal sentiment analysis,” _Soft Computing_ , vol. 24, pp. 18 431–18 445, 2020\.
  * [79] H. Sun, H. Wang, J. Liu, Y.-W. Chen, and L. Lin, “Cubemlp: An mlp-based model for multimodal sentiment analysis and depression estimation,” in _Proceedings of the 30th ACM International Conference on Multimedia_ , 2022, pp. 3722–3729.
  * [80] Y. Ramdhani, H. Mustofa, S. Topiq, D. P. Alamsyah, S. Setiawan, and L. Susanti, “Sentiment analysis twitter based lexicon and multilayer perceptron algorithm,” in _2022 10th International Conference on Cyber and IT Service Management (CITSM)_. IEEE, 2022, pp. 1–6.
  * [81] M. Z. Abedin, C. Guotai, F.-E. Moula, A. S. Azad, and M. S. U. Khan, “Topological applications of multilayer perceptrons and support vector machines in financial decision support systems,” _International Journal of Finance & Economics_, vol. 24, no. 1, pp. 474–507, 2019.
  * [82] V.-E. Neagoe, A.-D. Ciotec, and G.-S. Cucu, “Deep convolutional neural networks versus multilayer perceptron for financial prediction,” in _2018 International Conference on Communications (COMM)_. IEEE, 2018, pp. 201–206.
  * [83] D. Pehlivanli, S. Eken, and E. Ayan, “Detection of fraud risks in retailing sector using mlp and svm techniques,” in _Turkish Journal of Electrical Engineering & Computer Sciences_, vol. 27, no. 5. Tübitak, 2019, pp. 3657–3669.
  * [84] R. Makuyana, “Fraud detection in e-transactions using deep neural networks-a case of financial institutions in zimbabwe,” _International Journal of Scientific Research in Computer Science, Engineering and Information Technology_ , vol. 6, no. 9, pp. 1–7, 2020.
  * [85] F. Cheshmberah, H. Fathizad, G. Parad, and S. Shojaeifar, “Comparison of rbf and mlp neural network performance and regression analysis to estimate carbon sequestration,” _International Journal of Environmental Science and Technology_ , vol. 17, no. 8, pp. 3891–3900, 2020.
  * [86] Y. Liu, C. Li, X. Shi, and W. Wang, “Mlp-based regression prediction model for compound bioactivity,” _Frontiers in Bioengineering and Biotechnology_ , vol. 4, pp. 1–10, 2016.
  * [87] J. Park and S. Jo, “Approximate bayesian mlp regularization for regression in the presence of noise,” _Neural Networks_ , vol. 83, pp. 75–85, 2016.
  * [88] M. Taki, A. Rohani, F. Soheili-Fard, and A. Abdeshahi, “Assessment of energy consumption and modeling of output energy for wheat production by neural network (mlp and rbf) and gaussian process regression (gpr) models,” _Journal of cleaner production_ , vol. 172, pp. 3028–3041, 2018.
  * [89] M. Esmaeili, M. Osanloo, F. Rashidinejad, A. A. Bazzazi, and M. Taji, “Multiple regression, ann and anfis models for prediction of backbreak in the open pit blasting,” _Engineering with Computers_ , vol. 30, no. 4, pp. 549–558, 2014.
  * [90] S. Sharma, S. Sharma, and A. Athaiya, “Activation functions in neural networks,” _Towards Data Sci_ , vol. 6, no. 12, pp. 310–316, 2017.
  * [91] Intel Corporation. (2023) Shared Local Memory. [Online]. Available: {https://www.intel.com/content/www/us/en/docs/oneapi/optimization-guide-gpu/2024-0/shared-local-memory.html}
  * [92] ——. (2023) GRF mode. [Online]. Available: {https://www.intel.com/content/www/us/en/docs/oneapi/optimization-guide-gpu/2024-0/small-register-mode-vs-large-register-mode.html}
  * [93] I. Goodfellow, Y. Bengio, and A. Courville, _Deep learning_. MIT press, 2016.
  * [94] Intel Corporation. (2024) Xe GPU Architecture. Intel oneAPI Optimization Guide for GPU. [Online]. Available: https://www.intel.com/content/www/us/en/docs/oneapi/optimization-guide-gpu/2024-0/intel-xe-gpu-architecture.html
  * [95] S. Wang and P. Kanwar. (2019) BFloat16: The secret to high performance on Cloud TPUs. Article on bfloat16 data type. [Online]. Available: https://cloud.google.com/blog/products/ai-machine-learning/bfloat16-the-secret-to-high-performance-on-cloud-tpus
  * [96] Intel Corporation. (2023) joint_matrix extension. [Online]. Available: {https://github.com/intel/llvm/blob/sycl/sycl/doc/extensions/experimental/sycl_ext_matrix/sycl_ext_oneapi_matrix.asciidoc}
  * [97] ——. (2023) DPC++ Documentation joint_matrix_load. [Online]. Available: {https://intel.github.io/llvm-docs/doxygen/namespacesycl_1_1__V1_1_1ext_1_1oneapi_1_1experimental_1_1matrix.html#a525506bc79a9d1f675555150e7e97435}
  * [98] ——. (2023) joint_matrix extension. [Online]. Available: {https://github.com/intel/llvm/blob/sycl/sycl/doc/extensions/experimental/sycl_ext_matrix/sycl_ext_intel_matrix.asciidoc}
  * [99] G. Ofenbeck, R. Steinmann, V. Caparros, D. G. Spampinato, and M. Püschel, “Applying the roofline model,” in _2014 IEEE International Symposium on Performance Analysis of Systems and Software (ISPASS)_ , 2014, pp. 76–85.
  * [100] M. Tancik, B. Mildenhall, J. T. Barron, R. Martin-Brualla, N. Radwan, and P. P. Srinivasan, “Block-nerf: Scalable large scene neural view synthesis,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 1070–1080.
  * [101] Q. Wu, D. Bauer, Y. Chen, and K.-L. Ma, “Hyperinr: A fast and predictive hypernetwork for implicit neural representations via knowledge distillation,” _arXiv preprint arXiv:2304.00000_ , 2023.
  * [102] S. Hadadan, S. Chen, and M. Zwicker, “Neural radiosity,” _arXiv preprint arXiv:2105.12319_ , 2021.
  * [103] S. J. Garbin, M. Kowalski, M. Johnson, J. Shotton, and J. Valentin, “Stochastic texture filtering,” _arXiv preprint arXiv:2305.05810_ , 2023\.
  * [104] A. Yu, S. Fridovich-Keil, M. Tancik, Q. Chen, B. Recht, and A. Kanazawa, “Tensorf: Tensorial radiance fields,” _arXiv preprint arXiv:2203.10492_ , 2022.
  * [105] S. Fridovich-Keil, A. Yu, M. Tancik, Q. Chen, B. Recht, and A. Kanazawa, “Plenoxels: Radiance fields without neural networks,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 10 811–10 820.
  * [106] A. Yu, R. Li, M. Tancik, H. Li, R. Ng, and A. Kanazawa, “Plenoctrees for real-time rendering of neural radiance fields,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2021, pp. 10 486–10 495.
  * [107] J. Wynn and D. Turmukhambetov, “Diffusionerf: Regularizing neural radiance fields with denoising diffusion models,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2023, pp. 4180–4189.
  * [108] C. Uchytil and D. Storti, “A function-based approach to interactive high-precision volumetric design and fabrication,” _ACM Transactions on Graphics_ , vol. 43, no. 1, p. 3, 2023.
  * [109] D. Kim, M. Lee, and K. Museth, “Neuralvdb: High-resolution sparse volume representation using hierarchical neural networks,” _arXiv preprint arXiv:2208.04448_ , 2022.
  * [110] P. Wang, L. Liu, Y. Liu, C. Theobalt, T. Komura, and W. Wang, “Neus: Learning neural implicit surfaces by volume rendering for multi-view reconstruction,” _arXiv preprint arXiv:2106.10272_ , 2021.
  * [111] N. Raghavan, Y. Xiao, K.-E. Lin, T. Sun, S. Bi, Z. Xu, T.-M. Li, and R. Ramamoorthi, “Neural free-viewpoint relighting for glossy indirect illumination,” in _Computer Graphics Forum_ , vol. 42, no. 4. Wiley Online Library, 2023, p. e14885.
  * [112] S. Devkota and S. Pattanaik, “Efficient neural representation of volumetric data using coordinate-based networks.” in _Computer Graphics Forum_. Wiley Online Library, 2023, p. e14955.
  * [113] X. Sun, Z.-F. Gao, Z.-Y. Lu, J. Li, and Y. Yan, “A model compression method with matrix product operators for speech enhancement,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 28, pp. 2837–2847, 2020.
  * [114] V. Saragadam, J. Tan, G. Balakrishnan, R. G. Baraniuk, and A. Veeraraghavan, “Miner: Multiscale implicit neural representation,” in _European Conference on Computer Vision_. Springer, 2022, pp. 318–333.
  * [115] Y. Zhu, X. Xiao, W. Wu, and Y. Guo, “3d reconstruction of deformable linear objects based on cylindrical fitting,” _Signal, Image and Video Processing_ , vol. 17, no. 5, pp. 2617–2625, 2023.
  * [116] Y. Mao, Y. Wang, C. Wu, C. Zhang, Y. Wang, Y. Yang, Q. Zhang, Y. Tong, and J. Bai, “Ladabert: Lightweight adaptation of bert through hybrid model compression,” _arXiv preprint arXiv:2004.04124_ , 2020.
  * [117] D. Rebain, W. Jiang, S. Yazdani, K. Li, K. M. Yi, and A. Tagliasacchi, “Derf: Decomposed radiance fields,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2021, pp. 14 153–14 161.
  * [118] D. B. Lindell, J. N. Martel, and G. Wetzstein, “Autoint: Automatic integration for fast neural volume rendering,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2021, pp. 14 556–14 565.
  * [119] T. Neff, P. Stadlbauer, M. Parger, A. Kurz, J. H. Mueller, C. R. A. Chaitanya, A. Kaplanyan, and M. Steinberger, “Donerf: Towards real-time rendering of compact neural radiance fields using depth oracle networks,” in _Computer Graphics Forum_ , vol. 40, no. 4. Wiley Online Library, 2021, pp. 45–59.
  * [120] V. Sitzmann, S. Rezchikov, B. Freeman, J. Tenenbaum, and F. Durand, “Light field networks: Neural scene representations with single-evaluation rendering,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 19 313–19 325, 2021.
  * [121] H. Yu, J. Julin, Z. A. Milacski, K. Niinuma, and L. A. Jeni, “Dylin: Making light field networks dynamic,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2023, pp. 12 397–12 406.
  * [122] J. Cho, S. Nam, H. Yang, S.-B. Yun, Y. Hong, and E. Park, “Separable physics-informed neural networks,” _arXiv preprint arXiv:2306.15969_ , 2023\.
  * [123] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” _Journal of Computational physics_ , vol. 378, pp. 686–707, 2019.
  * [124] S. Sankaran, H. Wang, L. F. Guilhoto, and P. Perdikaris, “On the impact of larger batch size in the training of physics informed neural networks,” in _The Symbiosis of Deep Learning and Differential Equations II_ , 2022.
  * [125] R. Sharma and V. Shankar, “Accelerated training of physics-informed neural networks (pinns) using meshless discretizations,” _Advances in Neural Information Processing Systems_ , vol. 35, pp. 1034–1046, 2022.
  * [126] S. Wang, Y. Teng, and P. Perdikaris, “Understanding and mitigating gradient flow pathologies in physics-informed neural networks,” _SIAM Journal on Scientific Computing_ , vol. 43, no. 5, pp. A3055–A3081, 2021.
  * [127] V. Biesek and P. H. d. A. Konzen, “Burgers’ pinns with implicit euler transfer learning,” _arXiv preprint arXiv:2310.15343_ , 2023.
  * [128] T. Müller, A. Evans, C. Schied, and A. Keller, “Instant neural graphics primitives with a multiresolution hash encoding,” _ACM Transactions on Graphics (ToG)_ , vol. 41, no. 4, pp. 1–15, 2022.
  * [129] Z. Cai and M. Müller, “Clnerf: Continual learning meets nerf,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2023, pp. 23 185–23 194.
  * [130] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis, “Deepxde: A deep learning library for solving differential equations,” _SIAM review_ , vol. 63, no. 1, pp. 208–228, 2021.

###  -A Neural Radiance Fields

Recently, MLPs have been instrumental in the emerging Neural Rendering (NeRFs)
field [5]. NeRF uses MLPs to represent the scene as a continuous function
$f_{\theta}(x,y,z,\theta,\phi)=(r,g,b,\sigma)$, which is parameterized by a
MLP network that takes a 3D point $(x,y,z)$ and a viewing direction
$(\theta,\phi)$ as input and outputs a 4D vector $(r,g,b,\sigma)$ where
$(r,g,b)$ is the color and $\sigma$ is the density.

To render an image from a given camera pose, we use ray casting to sample
points along each ray and compute the accumulated color and density using
volume rendering by the following function:

$\displaystyle\mathbf{{C(r)}}=\int^{t_{far}}_{t_{near}}T(t)\sigma(r(t))c(r(t),d)dt$
(5)

where $C(r)$ is the final color of the ray $r$,
$T(t)=\exp(-\int_{t_{near}}^{t_{far}}\sigma(t^{\prime})dt^{\prime})$ is check
the occlusion through the ray traveling from the near plane $t_{near}$ to the
far plane $t_{far}$ and $r(t)=o+td$ is the ray equation with origin $o$ and
direction $d$.

To train the MLP, NeRF uses a combination of coarse and fine rendering losses.
The reconstruction loss measures the difference between the rendered color and
the ground-truth color for each pixel in the training images. The total loss
is defined as:

$L=\sum_{r\in
R}\|\hat{C}_{c}(r)-C_{gt}(r)\|_{2}^{2}+\|\hat{C}_{f}(r)-C_{gt}(r)\|_{2}^{2}$
(6)

where $r$ is the ray that passes through the pixel, $C_{gt}(r)$ is the ground
truth color of the pixel. By minimizing this loss, the MLP learns to represent
the scene as a neural radiance field that can synthesize novel views from any
camera pose.

###  -B Neural Compression

There has also been an increased focus on using MLPs for Neural Compression,
where MLPs aim to learn non-linear mapping as a function approximation to
compress data and preserve the most relevant information. To this end, the MLP
can be used as the encoder and/or the decoder in the compression pipeline. The
encoder transforms the original data into a compressed representation, and the
decoder aims to reconstruct the data. MLPs can be optimized to minimize the
difference between the original and the reconstructed image. The loss function
for neural data compression can be defined as:

$L=\lambda D(x,\hat{x})+R(z)$ (7)

where $x$ is the original image, $\hat{x}$ is the reconstructed image, $z$ is
the representation, $D$ is a distortion measure such as mean squared error,
$R$ is a rate measure such as entropy, and $\lambda$ is a trade-off parameter
that balances the distortion and the rate.

###  -C Partial Differential Equations (PDEs)

Partial Differential Equations (PDEs) have attracted much attention recently.
However, solving PDEs analytically is often impossible or impractical, and
numerical methods can be costly or inaccurate.

MLPs can be used to solve PDEs in a data-driven and physics-informed way. One
of the approaches is to use Physics-Informed Neural Networks (PINNs), which
leverage MLPs to learn and approximate the solutions to complex physical
systems. Given the underlying PDE and initial, boundary conditions embedded in
a loss function, a coordinate-based neural network is trained to approximate
the desired solution. PINNs take the form of:

$u(x,t)=f_{\theta}(x,t)$ (8)

where $u(x,t)$ is the unknown solution, $f_{\theta}(x,t)$ is an MLP with
parameters $\theta$, $x$ is the spatial variable, and $t$ is the temporal
variable. The MLP is trained to satisfy the boundary and initial conditions of
the PDE, as well as the PDE itself. The loss function can be defined as:

$L=\sum_{i=1}^{N}\|u(x_{i},t_{i})-f_{\theta}(x_{i},t_{i})\|^{2}+\sum_{j=1}^{M}\|F(x_{j},t_{j})\|^{2}$
(9)

where $N$ is the number of data points, $M$ is the number of collocation
points, $u(x_{i},t_{i})$ is the observed or prescribed value of the solution
at $(x_{i},t_{i})$, $f_{\theta}(x_{i},t_{i})$ is the predicted value of the
solution at $(x_{i},t_{i})$, and $F(x_{j},t_{j})$is the residual of the PDE at
$(x_{j},t_{j})$, which is computed by applying automatic differentiation to
the MLP. By minimizing this loss, the MLP learns to approximate the solution
of the PDE that is consistent with the data and the physics.