MatrixKAN:
Parallelized Kolmogorov-Arnold Network

In addition to the features already noted, KANs differ from MLPs in their theoretical foundation. Unlike MLPs, whose theoretical basis as a function modeling framework derives from the universal approximation theorem, KANs rely on the Kolmogorov-Arnold representation theorem, which provides that every multivariate continuous function $f$ can be represented as a finite composition of continuous, univariate functions and the binary operation of addition:

where $\phi_{q,p}:[0,1]\rightarrow\mathbb{R}$ and $\Phi_{q}:\mathbb{R}\rightarrow\mathbb{R}$ (Kolmogorov, 1957; Liu et al., 2024a). This traditional formulation of the Kolmogorov-Arnold representation theorem is analogous to a small neural network with (i) $n$ input neurons, (ii) a single hidden layer with $2n+1$ neurons, and (iii) a single output neuron. $\phi_{q,p}$ represent the edges between input neurons and neurons in the hidden layer, and $\Phi_{q}$ represent edges between neurons in the hidden layer and the output neuron.

To generalize this approach for modeling more complex systems and to facilitate training of this using traditional neural network regression strategies, KANs expand this two-layer network to networks of arbitrary width and depth and parameterize all edge functions using B-splines. The shape of a KAN is represented as $[n_{0},n_{1},...,n_{L}]$, where $n_{i}$ is the number of nodes in the $i^{th}$ layer. The activation value for a given node is represented as $x_{l,i}$ and the B-spline activation function between two nodes in adjacent layers is represented as $\Phi_{l,j,i}$, where $l$ denotes the layer, $i$ denotes the neuron within the given layer and $j$ denotes the connected neuron in the subsequent layer (Liu et al., 2024a):

The output of a given B-spline activation function is represented as $\tilde{x}\equiv\Phi_{l,j,i}(x_{l,i})$. For any given node, the activation value is equal to the sum of the outputs of the B-spline activation functions on each incoming network edge (Liu et al., 2024a):

where $j=1,...,n_{l+1}$.

To facilitate diversity in the learned activation function on each edge of a KAN, KANs implement activation functions using B-splines. B-splines are piecewise polynomial curves whose name derives from the elastic beams (i.e., splines) used by draftsmen to construct sweeping curves in ship design (Prautzsch et al., 2002). They are defined as affine combinations of control points $c_{i}$ and associated basis functions $B_{i}(x)$ (Liu et al., 2024a):

B-splines can be used to model a wide variety of shapes due to their ability to represent functions of arbitrary degree. A curve $spline(x)$ is called a B-spline of degree $k$ with knots $t_{o},...,t_{m}$, where $t_{i}\leq t_{i+1}$ and $t_{i}<t_{i+k+1}$ for all $i$, which delineate the intervals over which the B-spline is defined (see Figures 1 and 2), and control points $c_{o},...,c_{i}$, which control the local shape of the B-spline. Due to their control over the shape of the B-spline, KANs parameterize B-splines via their control points. Expressions for the basis functions $B_{i,k}(x)$ of a given B-spline are derived via the Cox-De Boor recursion algorithm, examples of which are shown in Figure 2:

As Equations 4 and 5 demonstrate, the value of a basis function $B_{i,k}(x)$ of degree $k$ is calculated recursively by reference to the position of the input value $x$ within the knot vector and the output of certain basis functions of degree $k-1$, and the final output of the B-spline is determined by multiplying all applicable basis functions by their associated control points $c_{i}$.

Given that basis functions of degree $k$ are defined by reference to basis functions of degree $k-1$, it follows that calculation of a degree $k$ basis function requires $k$ levels of sequential, recursive calls. This sequential dependency in the Cox-De Boor recursion algorithm prevents parallelization, resulting in a proportional increase in computation time as B-spline degree increases, even with memorization or dynamic programming techniques. Our testing shows that B-spline calculations constitute a significant portion of total execution time for KANs, ranging from 18% for degree-2 B-splines to over 50% for degree-30 B-splines, making B-spline computational efficiency a fundamental concern in improving the computational efficiency of KANs as a whole. While low-degree B-splines enable KANs to outperform MLPs (Liu et al., 2024a; Peng et al., 2024; Hu et al., 2024), our tests have demonstrated that high-degree B-splines are necessary for optimal performance in modeling certain functions, with optimal convergence occurring at degrees of 30 or higher. To bridge the gap between the need for high-degree B-splines and their high computational cost, we need a way to parallelize the recursive B-spline calculations for efficient execution on a parallel processor.

To address the computational cost of KANs, a number of approaches have been taken. For certain applications, despite the general additional computational overhead of KANs as compared to MLPs with similar parameter counts, KANs have been shown to produce superior performance with lower parameter counts, resulting in overall computational efficiency, as is the case with the convolutional KAN and U-KAN (Bodner et al., 2024; Li et al., 2024). However, these results involve hybrid KAN architectures and are domain-specific, suggesting the need for more general optimization techniques. To achieve more general computational efficiency, many modified KAN architectures have been proposed that replace underlying B-splines with different basis functions, such as radial basis functions (FastKAN / FasterKAN) (Li, 2024), trigonometric functions (e.g., Larctan-SKAN) and left-shift softplus functions (LSS-SKAN) (Chen & Zhang, 2024b, a), and rational functions (rKAN) (Aghaei, 2024). While each of these approaches has shown speedups over the B-spline-based KAN architecture, each basis function type has its own strengths and weaknesses in function approximation, and none identified here exhibit the local support feature of B-splines that allows for its granular function modeling capabilities. In certain research, adding wavelet functions to B-spline-based KANs (Wav-KAN) (Bozorgasl & Chen, 2024) or implementing neuron grouping / weight sharing to minimize the overall parameter count of B-spline-based KANs (Free Knots KAN) (Zheng et al., 2025) has yielded efficiency gains, but in all cases these speedups are only by a constant factor and none resolves the underlying computational complexity of the recursive B-spline calculations that is addressed by our proposed approach. That said, each is complementary to our approach and could be used in combination to achieve further improved computational efficiency.

Given the above theoretical analysis, we now present the procedure for precalculating the basis matrix $\Psi^{k+1}(i)$ of order $k+1$ (i.e., of degree $k$) and executing all B-spline calculations in a MatrixKAN network in a fully parallelized manner:

Step 1: Calculate the basis matrix $\Psi^{k}$. At model initialization, calculate the basis matrix applicable to B-splines of the degree, $k$, specified for the model. This results in a $k$ x $k$ tensor.

Step 2: Calculate applicable knot intervals and power bases $\begin{bmatrix}1&u&...&u^{k}\end{bmatrix}$, where $u=\frac{t-t_{i+k}}{t_{i+k+1}-t_{i+k}},t_{i+k}\leq t<t_{i+k+1}$ and $t_{o},...,t_{m}$ represents the knot vector⁴⁴4Equation 9 is described in terms of B-spline order, while the formula here is described in terms of B-spline degree.. Thus, for each B-spline, we must first calculate $u$, which requires determining the applicable knot interval for a given input value. This is achieved by comparing the input values to the knot values in each B-spline knot vector, returning two tensors—one representing the lower knot in the applicable grid interval for each value $t_{i+k}$ and one representing the upper knot in the applicable grid interval for each value $t_{i+k+1}$. These tensors are then used, along with the tensor containing the input values $t$, to calculate $u$, per the formula above. Once $u$ is calculated, the power bases tensors is then populated with the appropriate powers of $u$.

Step 3: Calculate applicable basis function outputs. Next, we calculate the basis function outputs for each B-spline based on Equation 9 by multiplying the basis matrix tensor calculated in Step 1 with the power bases tensor calculated in Step 2. To ensure that only applicable basis function outputs are calculated, we utilize the applicable knot intervals calculated in Step 2 to zero portions of the basis matrix tensor representing inapplicable basis functions. We then execute the final matrix multiplication, which results in a tensor containing the basis function outputs for all B-splines across all input values⁵⁵5The basis matrix tensor may be expanded to match the dimension of the power bases tensor, which has additional dimensions corresponding to the number of input values and number of dimensions of each input value..

Step 4: Calculate B-spline outputs. At this point, the output of Step 3 has completed the parallelized calculation of all equivalent recursive calculations required by the Cox-De Boor recursion algorithm and is equivalent to the basis function output tensor produced by KAN, so we can employ the original KAN algorithm from here forward to complete the B-spline output calculations.

All of the above operations are tensor operations implemented in PyTorch, and unlike the original KAN implementation, the recursive operations required to calculate basis function outputs have been replaced by matrix multiplication operations, which can be effectively parallelized.

In this section, we analyze and compare the computational efficiency of the algorithms implemented in KAN and MatrixKAN. For purposes of this discussion, we assume a network that mirrors what is presented in (Liu et al., 2024a), where all $L$ layers have equal width $n_{0}=n_{1}=...=n_{L}=N$ with each B-spline of degree $k$ on $G$ intervals.

Following the evaluation methodology presented in (Liu et al., 2024a), we conducted experiments on (i) the dataset generated by the hellokan.ipynp script (for the efficiency tests) and (ii) a random subset of the Feynman equations generated by the feynman.py script (for the performance tests), each of which is provided by (Liu et al., 2024b). For each test, we generated a training set containing 1000 training samples and a test set containing 1000 validation samples. After each training epoch, each model was validated against the test set.

For the efficiency tests, our aim was to compare the computation time of KAN⁶⁶6For all testing, we utilized the latest MultKAN implementation of KAN provided by (Liu et al., 2024b). and MatrixKAN relative to grid size, B-spline degree, and dataset size. We trained models using one NVIDIA 1080ti-8MB GPU for 20 steps, and to mirror the metrics in (Liu et al., 2024a), we measured computational efficiency in seconds per step, averaged across all steps. With this general setup, we ran three tests comparing–(1) training time vs. grid size, (2) training time vs. spline degree, and (3) training time vs. dataset size.

For performance tests, our aim was to demonstrate (i) the validity of MatrixKAN, by demonstrating model performance consistent with that of KAN and (ii) the superior performance of models using high-degree B-splines (i.e., of degree 4 or higher) relative to models using low-degree B-splines (i.e. of degree 3 or less) for certain modeling tasks. We trained KAN and MatrixKAN models in tandem using the shapes prescribed in (Liu et al., 2024a) Table 2 Feynman dataset. To ensure comparability between KAN and MatrixKAN models, models of each architecture were trained using identical training and test sample sets and identical random seed values. Training was performed using the Adam optimizer with a learning rate of 0.001. We trained models of each architecture for 1000 steps using the following values for the underlying B-spline degree: 2, 4, 6, 8, 10, 20, 30.

To compare the computational efficiency of MatrixKAN to KAN, we trained models of shape [2,2,1] and [2,5,1] using the LBFGS optimizer, each model initialized with a grid size of 3, a grid eps of 1, a seed value of 42, and grid update enabled.

For our performance tests, we compared the performance of MatrixKAN to KAN by plotting the loss level of each model against B-spline degree. The results of select tests are shown in Figures 7, 7, and 7.

Two things can be observed. First, MatrixKAN and KAN performed identically modeling all Feynman equations (i.e., completely overlapped curves). Although MatrixKAN and KAN perform B-spline calculations in different ways, the calculations are identical. Thus, when initialized with identical parameters and trained on identical datasets, models from each architecture should produce identical outputs and follow the same path of gradient descent to identical levels of convergence, which is demonstrated here. Second, the results demonstrate the superior performance of high-degree B-splines (i.e. of degree 4 or higher) in modeling various functions. Of the nine Feynman functions modeled, all but two models achieved optimal convergence at B-spline degree 4 or higher, with the majority of models converging at B-spline degree 6 or higher with loss level improvements of up to 27.0%. Further, Figure 7 demonstrates that for some modeling tasks model performance generally improves as a function of increasing B-spline degree and that B-splines of higher degree than tested here may yield improved results.