# KAN-SAs: Efficient Acceleration of Kolmogorov-Arnold Networks on Systolic Arrays

Sohaib Errabii, Olivier Sentieys, Marcello Traiola

► To cite this version:

Sohaib Errabii, Olivier Sentieys, Marcello Traiola. KAN-SAs: Efficient Acceleration of Kolmogorov-Arnold Networks on Systolic Arrays. IEEE/ACM Design, Automation & Test in Europe Conference (DATE) 2026, Apr 2026, Verona, Italy. hal-05372642

HAL Id: hal-05372642

https://inria.hal.science/hal-05372642v1

Submitted on 19 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# KAN-SAs: Efficient Acceleration of Kolmogorov-Arnold Networks on Systolic Arrays

Sohaib Errabii, Olivier Sentieys, Marcello Traiola

University of Rennes, CNRS, Inria, IRISA, Rennes, France

**Abstract**—Kolmogorov-Arnold Networks (KANs) have garnered significant attention for their promise of improved parameter efficiency and explainability compared to traditional Deep Neural Networks (DNNs). KANs' key innovation lies in the use of learnable non-linear activation functions, which are parametrized as splines. Splines are expressed as a linear combination of basis functions (B-splines). B-splines prove particularly challenging to accelerate due to their recursive definition. Systolic Array (SA)-based architectures have shown great promise as DNN accelerators thanks to their energy efficiency and low latency. However, their suitability and efficiency in accelerating KANs have never been assessed. Thus, in this work, we explore the use of SA architecture to accelerate the KAN inference. We show that, while SAs can be used to accelerate part of the KAN inference, their utilization can be reduced to 30%. Hence, we propose KAN-SAs, a novel SA-based accelerator that leverages intrinsic properties of B-splines to enable efficient KAN inference. By including a non-recursive B-spline implementation and leveraging the intrinsic KAN sparsity, KAN-SAs enhances conventional SAs, enabling efficient KAN inference, in addition to conventional DNNs. KAN-SAs achieves up to 100% SA utilization and up to 50% clock cycles reduction compared to conventional SAs of equivalent area, as shown by hardware synthesis results on a 28nm FD-SOI technology. We also evaluate different configurations of the accelerator on various KAN applications, confirming the improved efficiency of KAN inference provided by KAN-SAs<fnref n="1" />.

**Index Terms**—Kolmogorov-Arnold Networks, hardware accelerator, systolic array

## I. INTRODUCTION

The Kolmogorov-Arnold Network (KAN) is a neural network architecture [1] that has garnered interest and been adopted in various applications, such as time series analysis [2], recommender systems [3], and medical image segmentation [4]. The main interest of KANs lies in their improved parameter efficiency and explainability compared to conventional Deep Neural Networks (DNNs). This is enabled by replacing the conventional scalar weights with learnable spline-based activation functions, which are parameterized in a basis function (B-spline). However, this leads to an increase in computational complexity proportional to the size of the basis. Indeed, to compute KAN inference, instead of a single scalar multiply, each function in the basis must first be evaluated at the input, and then a linear combination of all the basis functions is performed. Moreover, the B-spline functions are evaluated recursively through the Cox-de Boor formula (See Eq. 3), which makes their acceleration challenging.

Recent artificial intelligence (AI) accelerators increasingly rely on spatial architectures, which have become the de facto standard for AI acceleration. Such architectures efficiently execute general

matrix multiplication (GEMM), the core operation underlying many AI workloads. Prominent examples include Google's Tensor Processing Unit (TPU) [5] and NVIDIA's Tensor Cores, first introduced in the Volta GPU microarchitecture [6]. The efficiency of spatial architectures stems from their ability to maximize data reuse and minimize data movement, which dominates the energy cost [7]. Building on this foundation, numerous works have proposed further optimizations tailored to specific workloads. For instance, Eyeriss [8] introduced dataflow optimizations to improve the efficiency of convolutional neural networks. Other efforts, such as SCNN [9], explored computation on compressed weights and activations to exploit zero weights stemming from model pruning and zero activations that occur in ReLU-based networks. The unique design of KANs limits the effectiveness of existing solutions in accelerating KAN inference compared to other AI applications. Specifically, the recursive nature of B-spline function computation creates a considerable bottleneck, making it challenging to take full advantage of the efficient GEMM executions in spatial architectures such as modern GPUs or Systolic Arrays (SAs) [10], which are the foundation of modern TPUs.

Therefore, research efforts have been focusing on new approaches to accelerate KANs (details in Sec. II-B). Among recent studies, compute-in-memory (CIM) approaches have been proposed [11], [12]. In such studies, different ways to approximate the learned non-linear KAN functions or their basis are utilized, such as piece-wise linear (PWL) approximation. However, no insights are offered into how spatial (non-CIM) architectures may be optimized for KANs. A recent approach, ArKANe [13], focuses on the B-spline evaluation bottleneck. It proposes an efficient dataflow acceleration method for the Cox-de Boor recursive formula, achieving a considerable speedup compared to CPU and GPU implementations. While all these efforts considerably improve knowledge of KAN acceleration, unfortunately, the current literature lacks solutions to efficiently accelerate KANs on spatial architectures, and particularly on SAs. To address this gap, this paper analyzes and utilizes the KAN properties to enhance SAs, enabling efficient end-to-end inference acceleration of KAN while maintaining the generality of the accelerator for non-KAN DNN workloads.

As shown in Section III, once the B-splines have been evaluated, their linear combination is nothing more than a GEMM operation, which can be accelerated on SAs. Regarding the B-spline computation, we observe that a direct floating-point implementation of the recursive evaluation of B-splines is quite costly. For inference-only acceleration, an efficient tabulation strategy of B-splines is possible, thanks to their properties [12]. Furthermore, B-spline computation results in an $N:M$ sparsity pattern in the values processed by the SA, thereby leading to low Processing Element

<footnote n="1">¹The source code will be publicly released upon acceptance.</footnote>
(PE) utilization. Designing a PE of the SA that can handle inputs with a KAN-specific $N:M$ sparsity pattern is needed to improve the overall efficiency of the SA [14]. To the best of our knowledge, there is currently no SA-based accelerator that, in addition to standard DNN workloads, can accelerate KANs by utilizing non-recursive B-spline computation and achieving high PE utilization. In summary, this work makes the following contributions:

* We show how a KAN layer can be transformed into a GEMM formulation for execution on a systolic array, and analyze the causes of the inefficiencies that lead to low throughput and poor PE utilization.

* Building upon our analysis, we include in KAN-SAs the needed architectural modifications to handle such inefficiencies. These consist of a PE (i) enhanced with a non-recursive B-spline unit and (ii) providing efficient handling of the $N:M$ sparsity pattern generated by B-splines.

* We synthesize KAN-SAs on a 28nm FD-SOI technology and evaluate its improvements compared to a conventional systolic array and to the state-of-the-art approach for B-spline acceleration. Our results show that the proposed KAN-SAs achieves 39.9% average improvement in PE utilization and 50% clock cycle reduction on average compared to conventional SAs. Moreover, the B-spline evaluation approach utilized in KAN-SAs achieves more than 72× improvement compared with the state-of-the-art approach.

## II. BACKGROUND

In this section, we present the core operation of a KAN layer and the implications for execution on a systolic array.

### A. KAN Layer

The main difference between a layer in a Multi-Layer Perceptron (MLP) and a KAN layer is illustrated in Fig. 1. It essentially involves replacing the weights at the connections with learnable activation functions $\phi_i$.

In general, a KAN layer as proposed by [1] can be expressed as follows

$$ \text{KANLayer}(x) = \sum w_i \phi_i(x) + w_b b(x) \quad (1) $$

The second term of the equation, which is not represented in Fig. 1 acts as a sort of bias. It is a typical MLP layer with some fixed non-linear activation applied before the dot product. It is typically a SiLU but we replace it with a ReLU. In the first term, the $w_i$ scales suggested by [1] offer better control of the magnitude of the activation functions $\phi_i$; however, at inference time, they can be absorbed in the functions. While this paper focuses on the first term, the proposed accelerator is capable of executing conventional MLP workloads and therefore any KAN layer that follows (1).

In an MLP layer, a matrix-vector multiplication between the inputs and the weights is followed by an activation using fixed functions such as ReLU. In the KAN layer, the learned functions are evaluated on the inputs, followed by the sum of all activations (see Fig. 1(b)).

These activation functions are learned through parameterization in some function basis ($\phi(x) = \sum_i c_i B_i(x)$). That is, the learnable parameters of the KAN layer are the coefficients $c_i$ for the linear combination in this basis. Since the parameterization is a linear combination, we can view the KAN layer as a matrix multiply, as

<figure><img src="image_2.png" /><figcaption>Fig. 1: A single layer of dimensions [2, 1] for MLP (left) and KAN (right).</figcaption></figure>

shown in Fig. 1(c). First, all the basis functions are evaluated at the inputs to obtain an intermediate matrix $B$ of dimensions ($M, N_b \times K$), with ($M, K$) being the dimensions of the input matrix, and $N_b$ the size of the function basis. The matrix multiplication is then performed with the coefficient matrix of dimensions ($N_b \times K, N$), with $N$ being the output dimension of the layer (in Fig. 1, $N = 1, K = 2$).

Therefore, the second part of the KAN layer can be efficiently accelerated by conventional hardware platforms, such as GPUs and SAs. However, the first part, i.e., the evaluation of the basis functions, poses a significant computational challenge. Indeed, the basis first proposed by [1] and also adopted by many KAN applications [2], [4], [15], [16], is the B-spline basis. This basis is defined by the following Cox-de Boor [17] recursion formula:

$$ B_{i,P}(x) = \frac{x - t_i}{t_{i+P} - t_i} B_{i,P-1}(x) + \frac{t_{i+P+1} - x}{t_{i+P+1} - t_{i+1}} B_{i+1,P-1}(x) \quad (2) $$

with

$$ B_{i,0}(x) = \begin{cases} 1 & \text{if } t_i \le x < t_{i+1} \\ 0 & \text{otherwise.} \end{cases} \quad (3) $$

Where the points $t_i$ are referred to as knot sequence of the B-spline. $P$ the spline degree and $G$ the grid size of discretization of the input domain of the layer. These are hyperparameters of the KAN layer. Importantly, the recursive formulation of Eq. 3 introduces dependencies between computational stages, making it not efficiently parallelizable on GPUs.

Fig. 2 illustrates the case of a uniform grid ($t_{k+1} - t_k = \Delta$) with $G = 3$ and $P = 3$. As shown in the figure, the grid must be extended to account for the contributions of B-splines that are non-zero within the input domain, $P$ extra intervals are added both before and after the input domain, leading to $G+2P$ total intervals, and $N_b = G + P$ B-spline functions.

### B. Related Work

Regarding the hardware acceleration of KANs, a few approaches have been proposed in the literature. As for the acceleration of B-spline evaluation, ArKANe [13] has proposed a dataflow acceleration strategy, based on an efficient unroll of the recursive Cox-de Boor formula from Eq. 2. While this can help during training, for inference it incurs significant latency overhead (i.e., several cycles for each B-spline function and input), as it still relies on the recursive computing of B-splines. Moreover, ArKANe approach requires floating-point MAC units and is tightly coupled
<figure><img src="image_3.png" /><figcaption>Fig. 2: B-spline basis functions of degree $P = 3$ on a uniform grid of size $G = 3$ for the input domain $[t_3, t_6]$ that is extended on both ends to account for the contributions of B-spline functions that are non-zero within this input domain.</figcaption></figure>

to AMD-Xilinx FPGA boards, as it relies on Xilinx's recent spatial computing architecture (the Adaptive Intelligent Engine, AIE), for deployment on AMD-Xilinx Versal ASoCs. Finally, ArKANe does not target the acceleration of the full KAN, but only of B-splines. For an inference accelerator, a tabulation-based strategy is much more attractive for its low resource cost and minimal latency. Moreover, it can be adapted for both floating-point and integer-only inference [18].

A tabulation-based B-spline implementation was proposed in [12] for a CIM chip based on resistive random access memory (RRAM). Unfortunately, besides not targeting systolic arrays, that work does not address the structured sparsity generated by the B-splines. Moreover, that strategy introduces constraints specific to the accelerated software application. For example, grid points were restricted to be powers of two, limiting the acceleration of more generic cases. In this work, instead, we aim to implement a more generic systolic-array-based accelerator. The only assumption we make is that of a uniform grid, which does not limit the generality of the B-spline unit. Indeed, as demonstrated by [1], it is possible to fine-grain the grid without retraining, using least squares to compute the new coefficients. This enables the approximation of non-uniform grids through finer uniform grids.

Another CIM approach based on CMOS technology was proposed by [11]. This strategy consists of approximating the learned non-linear KAN activations $\phi_i$ of Eq. 1 through piece-wise linear (PWL) functions. In contrast, our approach focuses on accelerating the evaluation of $\phi_i(x) = \sum_i c_i B_i(x)$ on a systolic array architecture through efficient tabular evaluation of B-splines $B_i(x)$ and efficient GEMM operations.

### III. GENERAL PRINCIPLE AND BASIS FUNCTION UNIT OF KAN-SAS

In this work, we target a weight-stationary systolic array for GEMM acceleration. In our design depicted in Fig. 3, the weights (i.e., the B-spline coefficients) are pre-loaded into the Processing Elements (PEs) and remain stationary while the input activations are propagated horizontally to other PEs and the resulting partial sums flow vertically towards an accumulator memory. This minimizes data movement by reusing weights and activations across processing elements, thereby enabling high energy efficiency. The memory hierarchy and additional architectural components needed for a practical DNN accelerator are beyond the scope of this paper.

#### A. General Architecture of KAN-SAS

Following the formulation presented in Section II-A, the only requirement for enabling hardware acceleration of KANs through a GEMM systolic array is the generation of the intermediate matrix $B$ of dimensions $(BS, (G + P)K)$, with $BS$ the batch size. $B$ represents the B-spline activations, as shown in Fig. 1(c). To minimize the data movement required for KAN inference, we aim to enable on-the-fly generation of this intermediate matrix $B$ close to the systolic array. This can be implemented by adding a basis function unit (the BSpline block in Fig. 3) that directly streams its values $B_i(x)$ into the systolic array. The rest of this section focuses on this basis function unit, whereas Section IV will detail the PEs of the systolic array.

<figure><img src="image_4.png" /><figcaption>Fig. 3: GEMM weight stationary (WS) systolic array with an additional B-spline unit for on-the-fly B-spline activation computation</figcaption></figure>

#### B. Basis Function Unit

Using the Cox-de Boor from Eq. 2 is not really feasible, as computing a single $P = 3$ function would require 20 multipliers. Expressing the B-spline in the canonical basis $1, x, \dots, x^3$ is still expensive, considering we need to evaluate the $G + P$ functions. Moreover, as mentioned previously, the B-spline functions exhibit several properties enabling an efficient Look Up Table (LUT)-based implementation.

1) Tabulation Strategy: A B-spline is invariant under a translation and scaling transformation of its knot sequence [19]. i.e.,

$$B_{t,P}(x) = B_{\alpha t+\beta,P}(\alpha x + \beta), \quad \alpha, \beta \in \mathbb{R}, \quad \alpha \neq 0$$

where $t = (t_0, \dots, t_{G+2P+1})$. If the grid is uniform, then $t_{k+1} - t_k = \Delta$. Applying the previous property with $\alpha = \frac{1}{\Delta}$ and $\beta = -\frac{t_0}{\Delta}$, we can map our B-spline defined on an arbitrary knot sequence $t_i$ to one defined on integer knots $0, \dots, G + 2P + 1$, also called cardinal B-spline, such as

$$B_{t,P}(x) = B_{[0, \dots, G+2P+1], P} \left( \frac{x - t_0}{\Delta} \right).$$

Furthermore, with the translation-invariance applied to $B_{k,P}$ we can further write

$$B_{t_k,P}(x) = B_{k,P}\left(\frac{x - t_0}{\Delta}\right) = B_{0,P}\left(\frac{x - t_0}{\Delta} - k\right). \quad (4)$$

This makes the B-spline independent of the KAN layer grid points $t_i$, since we only need to tabulate the function $B_{0,P}$, thereby enabling a ROM-based implementation of the B-spline unit. Of course, the unit must still perform the alignment required for $x$ shown in Eq. 4.

More importantly, the key property of B-splines for our purposes is their local support, which is given by the extreme knots used in
<figure><img src="image_5.png" /><figcaption>Fig. 4: Illustration of efficient tabular storing of half a cubic B-spline and computation of the others by aligning and shifting the input $x$. Only the values of $B_{0,3}$ corresponding to the black sample points on the horizontal axis quantizing the interval $[0, 2]$ are stored. For a given input $x$, all the B-spline values can be obtained by retrieving $B_{0,3}$ stored values at the aligned input $x_a$ and by shifting the evaluation point. E.g., $B_{t_{k-2,3}}(x) = B_{0,3}(2 - x_a)$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5: Simplified version of the B-Spline unit. The address is inverted, and the corresponding values are reverse-packed to obtain the remaining non-zero activations. The example values 0, 32 at the output correspond to $x_{addr} = 0$. The unit also returns the values at $\sim x_{addr} = 3$ in reverse. i.e. 127, 32.</figcaption></figure>

their definition of Eq. 3, $B_{t_k,P}(x) = 0$, $x \notin [t_k, t_{k+P+1}]$. This is also shown clearly in Fig. 2, and it implies that for any input $x$, there is at most $P + 1$ non-zero B-spline activations.

To summarize, the B-spline unit must first handle the alignment required by Eq. 4, which enables the use of stored $B_{0,P}$ values. And then, relying on translation-invariance, it must infer the value of all $P + 1$ non-zero B-splines. Moreover, to efficiently tabulate $B_{0,P}$, we must leverage the fact that the cardinal B-spline is symmetric with respect to the midpoint of the support (for $B_{0,P}$, that is $[0, \dots, P+1]$), $\frac{P+1}{2}$ [19]. Therefore, we only need to store half the B-spline corresponding to the interval $[0, \frac{P+1}{2}]$. For instance, for $P = 3$, we only need to sample points from the interval $[0, 2]$.

The example in Fig. 4 helps clarify the approach. For an input $x$, only $P + 1 = 4$ non-zero B-splines exist: $B_{t_{k-3,3}}(x)$, $B_{t_{k-2,3}}(x)$, $B_{t_{k-1,3}}(x)$, $B_{t_k,3}(x)$. The aligned $x_a$ input shown in the figure ensures that these values correspond respectively to $B_{-3,P}(x_a)$, $B_{-2,3}(x_a)$, $B_{-1,3}(x_a)$, $B_{0,3}(x_a)$.

2) Implementation of the BSpline Unit: Fig. 5 shows the table storing the B-spline values of $B_{0,P}$ required to infer all the non-zero B-spline activations. The Align unit implements the previously mentioned alignment in Eq. 4 and the Compare unit performs an interval search to compute $k$, necessary for both the alignment as well as the systolic array, as will be shown in Sec. IV. The LUT address for retrieving the stored value for $x$ is the quantized aligned input $x_a \in [0, 1]$ to the integer range $[0, 255]$.

The unit must obtain $x_{addr}$ using integer arithmetic based on its quantized inputs $x_q$ and $t_q$. Using affine integer quantization scheme with Eq. 4 we obtain the LUT address computed in Align unit as

$$x_{addr} = \text{clip}((G + 2P)(x_q - t_{q0}) - 255 \times k, 0, 255). \quad (5)$$

Since we use an address based on $x_a \in [0, 1]$ and we need to store values in $[0, 2]$. We simply store the value corresponding to $x_a + 1$ in the same address as $x_a$. This value corresponds to $B_{t_{k-1,3}}(x)$ as shown in Fig. 4, hence the two values per row in Fig. 5. For the other values $B_{t_{k-3,3}}(x)$, $B_{t_{k-2,3}}(x)$, which correspond to $B_{0,3}(1 - x_a)$ and $B_{0,3}(1 - x_a + 1)$, they are stored in the address

corresponding to $1 - x_a$, i.e., $255 - x_{addr}$. This is implemented in the inversion unit $\sim$ in Fig. 5.

## IV. LEVERAGING B-SPLINES PROPERTIES IN KAN-SAS PROCESSING ELEMENTS

In this section, we present the analysis of the B-spline properties presented previously, and how they enable the efficient acceleration of KAN networks on SA-based accelerators. In particular, we show how to provide efficient handling of the $N:M$ sparsity pattern generated by B-splines.

### A. On the Inherent Sparsity of B-splines

As illustrated in Fig. 2 and explained in Section III-B, B-splines have local support, i.e., the B-spline $B_{t_k,P}$ is only non-zero within $[t_k, t_{k+P+1}]$. Therefore, if $x \in [t_k, t_{k+1}]$, then only the $P + 1$ functions $B_{t_{k-P}}(x), \dots, B_{t_k}(x)$ are non-zero among all $G + P$ functions. Moreover, when $x$ is within the grid extension ($k < P$ or $k > G + P - 1$), there are even fewer non-zero B-splines. For example, with $G = 10, P = 3$, each input would at most contribute 4 non-zero B-spline activations among 13, which leads to at most $\frac{4}{13} \approx 30\%$ PE utilization (i.e., computations involving non-zero B-spline activations) in the SA (as shown later in Fig. 8).

To address this inefficiency, the B-spline unit must output only the $P+1$ contiguous non-zero activations along with the integer index $k$ specifying their positions among all the $G+P$ B-splines. This nice property of B-splines enables the use of an $N:M$ structured sparsity-aware PE (where $N = P+1, M = G+P$) to avoid useless multiplications with zero. This means that while the classical scalar PE in Fig. 3 performs $psum_i + c_i B_i(x)$, a KAN-optimized PE has to be vectorized to perform $psum_i + \sum_{i=0}^{P} c_{k-i} B_{k-i}(x)$, i.e. only involving $B_{k-i}(x)$, $i \in [0, P]$ non-zero B-spline values. This $N:M$ structured pattern, also known as density bound block (DBB), has been previously enforced in conventional DNNs static weights [20]–[22] through ad-hoc pruning. Conversely, in the KAN case, this $N:M$ sparsity appears dynamically in the B-spline activations, and it is guaranteed by their local support property.

### B. Designing the Processing Elements of KAN-SAS

A conventional systolic array uses a scalar processing element performing a multiply-accumulate. To harness the benefit of the intrinsic B-spline sparsity, we design an $N:M$ sparsity-aware vector processing element performing a multiply-accumulate between the $N$ non-zero B-spline values and their corresponding coefficients among all the $M$ coefficients.
<figure><img src="image_7.png" /><figcaption>Fig. 6: Proposed KAN-SAs PE architecture showing how the B-spline unit interacts with the *N*:M vector PE. Different from the scalar PE SA in Fig. 3, each B-spline unit streams its output to a single row of *N*:M PEs.</figcaption></figure>

This is shown in the Fig. 6, where each B-spline unit sends the *P* + 1 non-zero activations to its corresponding row in the array along with in the index *k*₀ as a control signal for the multiplexer selecting the right *N* = *P* + 1 coefficients among *M* = *G* + *P*. Therefore, while the *N* multiplications occur in parallel, the additional *M*-to-*N* multiplexer and multi-operands adder for summing the partial products do increase the delay of the critical path. Similar to the activations, the indices *k* are propagated and reused horizontally along the array.

## V. EVALUATION

In this section, we evaluate the proposed KAN-SAs architecture against the conventional SA setup. In the conventional SA, we assume B-spline units feeding a systolic array with scalar PEs, not able to handle the intrinsic *N*:M B-spline sparsity. We synthesized the conventional SA solution and the proposed one with Synopsys Design Compiler v2025.6, targeting the ST 28nm FD-SOI PDK [23]. The integer-only implementation, quantized as proposed by [18], was validated with a software baseline. We observe less than 1% in accuracy drops for all the models (e.g., MNIST-KAN drops from 96.58% to 96.0%).

### A. PE Comparison

In this section, we compare the hardware synthesis results for the scalar conventional PE and the *N*:M PE, as provided in Table I. Post-synthesis power estimation was carried out using activity-based analysis at a frequency of 500MHz. We first study the delay for different sparsity parameters *N* and *M*. As mentioned in Section IV-B, there is an increase in the delay of the critical path due to both the *M*-to-*N* multiplexer and the *N* + 1-operands adder. For instance, when increasing *N* from 2:6 to 4:6 in Table I, the increase in delay is due to the adder since the multiplexing occurs in parallel. Similarly, when increasing *M*, the adder is not affected, but the multiplexer's delay increases. As a consequence of the added logic, as shown in the table, the power of *N*:M PEs also increases with respect to 1:1. The switching activity was extracted from simulation traces for a typical KAN workload. The coefficients are loaded in the PE and then reused for several cycles with different activations. However, the normalized energy reported in the table shows the advantage of the *N*:M PE, as a KAN workload running on an *N*:M PE array takes (*G* + *P*) times fewer cycles than a scalar PE, with *N* = *P* + 1, *M* = *G* + *P*. The energy is estimated based on the reported power in the table and the number of cycles a KAN workload requires to run on a

<table><caption>TABLE I: ST28nm FD-SOI post-synthesis delay and power for 8-bit inputs and 32-bit output PE at a target frequency of 500 MHz and estimated normalized energy. Columns refer to the sparsity pattern *N*:M; 1:1 represents the scalar PE.</caption><thead><tr><th>N:M</th><th>1:1</th><th>1:2</th><th>2:4</th><th>2:6</th><th>4:6</th><th>4:8</th></tr></thead><tbody><tr><th>Delay (ns)</th><td>1.02</td><td>1.05</td><td>1.15</td><td>1.19</td><td>1.28</td><td>1.31</td></tr><tr><th>Power (mW)</th><td>0.35</td><td>0.40</td><td>0.62</td><td>0.77</td><td>0.98</td><td>1.12</td></tr><tr><th>Normalized Energy</th><td>1.00</td><td>0.57</td><td>0.44</td><td>0.37</td><td>0.47</td><td>0.40</td></tr></tbody></table>

scalar PE by multiplying the power in Table I by the number of cycles needed to run a typical KAN workload. The 1:1 PE takes (*G* + *P*) times more cycles than a PE handling *N*:M sparsity. E.g., for *P* = *G* = 3 (sparsity 4:6), the 1:1 PE will need 6× more cycles than 4:6.

### B. B-spline Acceleration Comparison with Previous Work

As already mentioned, ArKANe [13] proposes an acceleration of the recursive implementation of B-splines. A direct comparison is not straightforward, as ArKANe B-spline implementation evaluates the B-splines using floating-point arithmetic across multiple cycles for the objective of KAN training acceleration. Moreover, the reported results correspond to a mapping on Xilinx FPGA over multiple AIE tiles. However, to propose a fair comparison, we estimate the possible area consumption for ArKANe proposal and compare it to KAN-SAs B-spline unit. ArKANe wavefront algorithm requires *P* + 1 processing elements and has a latency of (*P* + 1) × PE<sub>latency</sub> cycles to evaluate a specific B-spline *B*<sub>*i*,*P*</sub>. Then, thanks to pipelining, it requires *G* + *P* − 1 cycles for all *G* + *P* activations. Hence, we can compute the ArKANe number of cycles for *M* inputs as (*P* + 1)PE<sub>latency</sub> + *G* + *P* − 1 + *M*. We consider an existing single-precision FMA implementation, FPMax [24], as a reference for the ArKANe FP32 multiply-accumulate (FMA) operation in each PE. In FPMax, a single-precision FMA circuit has PE<sub>latency</sub> = 4 and is estimated to occupy approximately 0.0081mm² of standard cell area. By contrast, our tabulation-based B-spline unit occupies 450µm², and requires at most a single cycle to retrieve the values of all *G* + *P* B-splines for a certain input. Therefore, in the same estimated area for ArKANe, i.e., 4 × 0.0081mm² we can fit 72 B-spline units to feed 72 rows of the systolic array. Tabulating the B-splines can offer a minimum of 72× speedup for high values of *M* over the recursive floating-point implementation.

### C. KAN Applications Benchmark

We perform a design space exploration of different configurations of the proposed architecture on a representative set of KAN applications collected from prior work and reported in Table II. Each of the applications contributes a certain number of KAN workloads, i.e., matrix multiplications where the left matrix is B-spline activations *B* mentioned in Section II-A and the right matrix is the coefficients. For instance, the application Catch22-KAN relies on a single KAN layer [22, X] where *X* is the number of classes for one of the UCR datasets [25]. This layer implies a matrix *B* of dimensions (*BS*, 22 × (*G* + *P*)). Some of these workloads also include the MLP bias term mentioned in Section II-A, which is included in the evaluation. Some studies explore different values of *G* and *P*. However, we limit our evaluation to workloads with
<table><caption>TABLE II: Collected KAN workloads from prior work. X in CF-KAN takes values from [2810, 34395, 6969]. In Catch22-KAN, X is the number of classes in UCR time series dataset, and we consider it smaller than 60 in our experiments.</caption><thead><tr><th>Application</th><th>Layers</th><th>G</th><th>P</th></tr></thead><tbody><tr><td>5G-STARDUST [2]</td><td>[168, 40, 40, 40, 24]</td><td>5</td><td>3</td></tr><tr><td>Catch22-KAN [26]</td><td>[22, X]</td><td>3</td><td>3</td></tr><tr><td>CF-KAN [3]</td><td>[X, 512, X]</td><td>2</td><td>3</td></tr><tr><td>U-KAN [4]</td><td>[512, 1024, 512], [512, 512]</td><td>5</td><td>3</td></tr><tr><td>GKAN [15]</td><td>[200, 16, 7], [100, 20, 7]</td><td>2,3</td><td>1,2,3</td></tr><tr><td>Prefetcher [27]</td><td>[5, 64, 128]</td><td>4</td><td>3</td></tr><tr><td>MNIST-KAN [28]</td><td>[784, 64, 10]</td><td>10</td><td>3</td></tr><tr><td>ResKAN18 [29]</td><td>20 ConvKAN layers</td><td>3</td><td>3</td></tr></tbody></table>

$P \le 3$. The applications MNIST-KAN and ResKAN18 are models we implemented and trained on MNIST [30] and CIFAR10 [31], respectively. MNIST-KAN is a two-layer KAN [784, 64, 10], and the ResKAN18 is a ResNet18 architecture, where the scalar filter weights in convolutions are replaced by learnable splines (referred to as ConvKAN or KANConv in prior works [16], [32]). The average results reported in Figs. 7a and 7b consider different dimensions *R* rows and *C* columns for the array. In the graphs, for reference, we mark some points that correspond to square SAs (e.g., 2x2, 4x4, etc.). The other parameters are fixed as int8 multiplication and int32 accumulation, *G* = 5 and *P* = 3. Hence, results are averaged over all collected workloads except MNIST-KAN, as it requires *G* = 10. Since the main objective of the study is to evaluate the proposed improvements, we focus solely on B-spline sparsity without considering other dynamic sources of sparsity, such as zero coefficients (i.e., weights) or activations. The workloads are tiled for running on the weight-stationary SA.

In Fig. 7a, we plot, for both conventional SA and KAN-SAs, the average PE utilization across all applications vs. the area obtained post-synthesis, for various sizes of the SA ($R \times C$). The figure

<figure><img src="image_8.png" /><figcaption>Fig. 7: Average PE utilization (a) and runtime in clock cycles (b) across all applications vs. the area obtained post-synthesis for both conventional SA and KAN-SAs, for various sizes of the SA ($R \times C$).</figcaption></figure>

clearly shows that the proposed solution consistently improves the average PE utilization of the SA for different design parameters. However, it can also be noticed that, while KAN-SAs consistently provides PE utilization > 65%, there is still some effect preventing 100% utilization. The PE utilization is fundamentally impacted by two distinct effects. The first and most important one is addressed by KAN-SAs, which is the sparsity generated by the B-splines. The second is imperfect tiling, i.e., some workloads' dimensions are not

a multiple of the SA dimensions, and is intrinsic to the workload being executed. To give more details, in Fig. 8, we show the per-application average PE utilization for specific configurations of the arrays with similar areas, i.e., $0.47mm^2$ for KAN-SAs $16 \times 16$ and $0.50mm^2$ for $32 \times 32$ scalar PE SA. Since KANs are still in their

<figure><img src="image_9.png" /><figcaption>Fig. 8: PE Utilization(%) for KAN workloads.</figcaption></figure>

early stages, most applications have small dimensions. This leads to a noticeable imperfect tiling effect in these applications. Also, the values of *G* impact the utilization. When *G* has high values, two effects are present: (i) the overall number of parameters also increases, and (ii) the sparsity grows, as shown in Sec. IV. The first effect intrinsically reduces the imperfect tiling, and the second reduces the utilization in conventional SA. For instance, MNIST-KAN with *G* = 10 has low utilization (30%) when executed on SA because of the high *G* value. Conversely, on KAN-SAs, it has a high utilization rate (99.25%) because the high *G* value leads to a large number of parameters, thereby reducing the impact of imperfect tiling. As can be seen from Table II, ResKAN18, CF-KAN, U-KAN, and Prefetcher applications are large compared to the others, which reduces the imperfect tiling impact, even for low *G* values. Additionally, low *G* values result in higher PE utilization, also for conventional SA. The absolute utilization improvement of KAN-SAs over conventional SA is 39.9% on average over all the applications, with a maximum of 69.3% for MNIST-KAN.

Finally, Fig. 7b shows, for both conventional SA and KAN-SAs, the average runtime (in clock cycles) across all applications vs. the area. KAN-SAs consistently provides reduced runtime by $2 \times$ compared to conventional SA at the same area and reduced area by $2 \times$ at the same runtime. The improvement stems from the fact that while the scalar PEs in the conventional SA can load weight tiles of size $R \times C$, the $N \times M$ PEs in KAN-SAs load ($R \times M, C$) tiles of B-spline activations and ($R \times N, C$) tiles of non-KAN workloads (i.e., the MLP term in Eq. 1), leading to faster execution.

## VI. CONCLUSION

In this paper, we presented KAN-SAs, the first systolic-array-based efficient hardware accelerator for KANs. KAN-SAs efficiency is enabled by a thorough analysis of the inefficiencies that the B-spline functions inherently introduce, i.e., recursive nature and sparsity. KAN-SAs addresses them through an enhanced Processing Element including a non-recursive B-spline unit and efficient handling of B-spline sparsity. Results show that, compared to conventional SAs, KAN-SAs achieves 39.9% average PE utilization improvement, 50% average clock cycle reduction, and $72 \times$ speedup in B-spline evaluation compared to the state-of-the-art approach.
ACKNOWLEDGMENT

This work was supported by the French National Research Agency (ANR) through the RADYAL project ANR-23-IAS3-0002.

REFERENCES

[1] Z. Liu, Y. Wang, S. Vaidya, F. Ruehle, J. Halverson et al., "Kan: Kolmogorov-arnold networks," 2024, arXiv:2404.19756 [cs.LG].

[2] C. J. Vaca-Rubio, L. Blanco, R. Pereira, and M. Caus, "Kolmogorov-Arnold Networks (KANs) for Time Series Analysis," Sep. 2024, arXiv:2405.08790 [esssp].

[3] J.-D. Park, K.-M. Kim, and W.-Y. Shin, "Cf-kan: Kolmogorov-arnold network-based collaborative filtering to mitigate catastrophic forgetting in recommender systems," 2024, arXiv:2409.05878 [cs.IR].

[4] C. Li, X. Liu, W. Li, C. Wang, H. Liu et al., "U-KAN Makes Strong Backbone for Medical Image Segmentation and Generation," Aug. 2024, arXiv:2406.02918 [essxiv].

[5] N. P. Jouppi, C. Young, N. Patil, D. Patterson, G. Agrawal et al., "In-Datacenter Performance Analysis of a Tensor Processing Unit," in ACM/IEEE 44th Annual International Symposium on Computer Architecture (ISCA), 2017, pp. 1-12.

[6] "Nvidia tesla v100 gpu architecture: The world's most advanced data center gpu," Whitepaper, https://images.nvidia.com/content/volta-architecture/pdf/volta-architecture-whitepaper.pdf, 2017.

[7] M. Horowitz, "1.1 computing's energy problem (and what we can do about it)," in 2014 IEEE International Solid-State Circuits Conference Digest of Technical Papers (ISSCC), 2014, pp. 10-14.

[8] Y.-H. Chen, T. Krishna, J. S. Emer, and V. Sze, "Eyeriss: An Energy-Efficient Reconfigurable Accelerator for Deep Convolutional Neural Networks," IEEE Journal of Solid-State Circuits, vol. 52, no. 1, pp. 127-138, Jan. 2017.

[9] A. Parashar, M. Rhu, A. Mukkara, A. Puglielli, R. Venkatesan et al., "Scnn: An accelerator for compressed-sparse convolutional neural networks," 2017, arXiv:1708.0448 [cs.NE].

[10] H. T. Kung, "Why systolic architectures?" Computer, vol. 15, no. 1, pp. 37-46, Jan. 1982.

[11] C. Sudarshan, P. Manea, and J. P. Strachan, "A Kolmogorov-Arnold Compute-in-Memory (KA-CIM) Hardware Accelerator with High Energy Efficiency and Flexibility," Jan. 2025, preprint 10.21203/rs.3.rs-5804189/v1.

[12] W.-H. Huang, J. Jia, Y. Kong, F. Waqar, T.-H. Wen et al., "Hardware Acceleration of Kolmogorov-Arnold Network (KAN) for Lightweight Edge Inference," in 30th Asia and South Pacific Design Automation Conference. ACM, Jan. 2025, pp. 693-699.

[13] Y. Wu and M. T. Arafin, "Arkane: Accelerating kolmogorov-arnold networks on reconfigurable spatial architectures," IEEE Embedded Systems Letters, pp. 1-1, 2025.

[14] X. He, S. Pal, A. Amarnath, S. Feng, D.-H. Park et al., "Sparse-TPU: adapting systolic arrays for sparse matrices," in 34th ACM International Conference on Supercomputing, Jun. 2020, pp. 1-12.

[15] M. Kiamari, M. Kiamari, and B. Krishnamachari, "GKAN: Graph Kolmogorov-Arnold Networks," Jun. 2024, arXiv:2406.06470 [cs].

[16] A. D. Bodner, A. S. Tepsich, J. N. Spolski, and S. Pourteau, "Convolutional Kolmogorov-Arnold Networks," Mar. 2025, arXiv:2406.13155 [cs].

[17] C. De Boor, "On calculating with b-splines," Journal of Approximation theory, vol. 6, no. 1, pp. 50-62, 1972.

[18] B. Jacob, S. Kligys, B. Chen, M. Zhu, M. Tang et al., "Quantization and training of neural networks for efficient integer-arithmetic-only inference," 2017, arXiv:1712.05877 [cs.LG].

[19] T. Lyche, C. Manni, and H. Speleers, Foundations of Spline Theory: B-Splines, Spline Approximation, and Hierarchical Refinement. Springer International Publishing, 2018, pp. 1-76.

[20] Z.-G. Liu, P. N. Whatmough, Y. Zhu, and M. Mattina, "S2ta: Exploiting structured sparsity for energy-efficient mobile cnn acceleration," 2022, arXiv:2107.07983 [cs.AR].

[21] H.-J. Kang, "Accelerator-aware pruning for convolutional neural networks," IEEE Transactions on Circuits and Systems for Video Technology, vol. 30, no. 7, pp. 2093-2103, 2020.

[22] M. Zhu, T. Zhang, Z. Gu, and Y. Xie, "Sparse tensor core: Algorithm and hardware co-design for vector-wise sparse neural networks on modern gpus," in 52nd Annual IEEE/ACM International Symposium on Microarchitecture (MICRO), ser. MICRO-52, 2019, p. 359-371.

[23] A. Cathelin, "Fully Depleted Silicon on Insulator Devices CMOS: The 28-nm Node Is the Perfect Technology for Analog, RF, mmW, and Mixed-Signal System-on-Chip Integration," IEEE Solid-State Circuits Magazine, vol. 9, no. 4, pp. 18-26, 2017.

[24] J. Pu, S. Galal, X. Yang, O. Shacham, and M. Horowitz, "FPMax: a 106GFLOPS/W at 217GFLOPS/mm2 Single-Precision FPU, and a 43.7GFLOPS/W at 74.6GFLOPS/mm2 Double-Precision FPU, in 28nm UTBB FDSOI," 2016, arXiv:1606.07852 [cs.AR].

[25] H. A. Dau, A. Bagnall, K. Kamgar, C.-C. M. Yeh, Y. Zhu et al., "The ucr time series archive," IEEE/CAA Journal of Automatica Sinica, vol. 6, no. 6, pp. 1293-1305, 2019.

[26] A. Ismail-Fawaz, M. Devanne, S. Berretti, J. Weber, and G. Forestier, "Feature-based time series classification with kolmogorov-arnold networks," https://github.com/MSD-IRIMAS/Simple-KAN-4-Time-Series, 2024.

[27] D. Kulkarni, B. Bhammar, H. Thaker, P. Dhobi, R. P. Gohil et al., "A case for kolmogorov-arnold networks in prefetching: Towards low-latency, generalizable ml-based prefetchers," 2025, arXiv:2504.09074 [cs.AR].

[28] "An efficient implementation of kolmogorov-arnold network," https://github.com/Blealtan/efficient-kan.

[29] V. Starostin, "Convolutional kan layer," https://github.com/StarostinV/convkan, 2024.

[30] L. Deng, "The mnist database of handwritten digit images for machine learning research [best of the web]," IEEE Signal Processing Magazine, vol. 29, no. 6, pp. 141-142, 2012.

[31] A. Krizhevsky, "Learning multiple layers of features from tiny images," 2009.

[32] I. Drokin, "Kolmogorov-arnold convolutions: Design principles and empirical studies," 2024, arXiv:2407.01092 [cs.AR].