Lookup multivariate Kolmogorov-Arnold Networks
Authors: Sergey Pozdnyakov
Philippe Schwaller
Venue: ICLR 2026 Poster (2026)
Keywords: KAN
inference efficiency
CUDA kernels
Primary area: other topics in machine learning (i.e., none of the above)
TL;DR: We propose a fully connected layer that decouples inference efficiency from the number of trainable parameters and empirically find it to be Pareto optimal across a wide range of macro-architectural backbones.

Abstract: High-dimensional linear mappings, or linear layers, dominate both the parameter count and the computational cost of most modern deep-learning models. We introduce a general-purpose drop-in replacement, lookup multivariate Kolmogorov-Arnold Networks (lmKANs), which deliver a substantially better trade-off between capacity and inference cost. Our construction expresses a general high-dimensional mapping through trainable low-dimensional multivariate functions. These functions can carry dozens or hundreds of trainable parameters each, and yet it takes only a few multiplications to compute them because they are implemented as spline lookup tables. Empirically, lmKANs reduce inference FLOPs by up to 6.0× while matching the flexibility of MLPs in general high-dimensional function approximation. In another feedforward fully connected benchmark, on the tabular-like dataset of randomly displaced methane configurations, lmKANs enable more than 10× higher H100 throughput at equal accuracy. Within the framework of Convolutional Neural Networks, lmKAN-based CNNs cut inference FLOPs at matched accuracy by 1.6–2.1× and by 1.7× on the CIFAR-10 and ImageNet-1k datasets, respectively. Our code, including dedicated CUDA kernels, is available online at https://github.com/schwallergroup/lmkan.

Published as a conference paper at ICLR 2026
LOOKUP MULTIVARIATE KOLmOGOROY-ARNOI. NETWORKS
Sergey Pozdnyakov' sergey. sergey.pozdnyakove pozdnyakov@epfl.ch
Philippe Schwaller
philippe. wtapppessccaallleepl...oo ch
École Polytochnique Fédéralc de Lausanne (EPFL)
National Centre of Competence in Research (NCCR) Catalysis
A ABSTRACT
High-dimensional linear mappings, or lincar layers, dominate both the param- eter count and the computational cost of most modern deep-lcarning models. We introduce a general-purpe purpose drop-in replacement, lookup multivariate Kolmogorov-Arnold Networks (ImKANs), which deliver substantially bet- ter trade-off between capacity and inference cost. Our construction expresses general high- ligensiinn mapping through trainable low- owwdiiensioona mul- tivariate functions. These functions can carry dozens or hundreds trainable parameters each, and yet it takes only few multiplications to compute them because they are implemented as spline lookup tables. Empirically, ImKANs reduce inference FLOPs by up 6.0x while matching the flexibility of MLPs in general high-dimensional function approximation. In another feed- forward fully connected benchmark, on the tabular-like dataset of randomly displaced methane configurations, ImKANs enable more than 10x higher H100 throughput at equal accuracy. Within the framework Convolutional Neural Networks, ImK AN-based CNNs cut inference FTOOPs at matched accuracy by 1.6-2.1x and by 1.7x on the CIFAR-10 and ImageNet-1k 1k datasets, respectively. Our code, including dedicated CUDA kernels, is available online at https:/ / tigghh...b.... com/ http//hhawalllrrrppaeeeo /mmkan.
General function approximation
ImKAN
60-1
Text RASE/SD
3x 10
Mean Squared
10
1x-1
30-4
10a
10
10-0
10-7
Inference FLOPs
H100 seconds/mleccll
llarge batch sizw)
CIFAR-10 10
ImageNet
80
Test Accurcc (5)
55
70
2.1%
1.7%
65
45
Vallaaaton
55
MLP CAN
40
ImKAN CNN
ImKAN CNN
10a
100
10f
10
20'
Inference FLOPs
Inference FLOPs
Figure 1: Performance summary. Soe Soc. for details.
INTRODUCTION
With sufficient amount of training data, the capabilities of deep-learning models systemati- cally improve with the number of trainable parameters (Zhai et al., 2022; Kaplan et al., 2020). However, deploying very large models is challenging because of the associated inference cost.

---

Published as a conference paper at ICLR 2026
In most models, high-dimensional linear mappings dominate both the parameter count and the computational cost. Standard multilayer perceptrons (MLPs) alternate linear layers with activations, and sometimes with a few other layers (loffe & Szegedy, 2015; Hinton et al., 2012). Tf is the width of a layer, the parameter count and inference cost of these linear mappings scale O(N2), whereas most other layers scale only O(N). The same observation holds for many other architectures. Transformers (Vaswani et al., 2017), when applied very long sequences, are one of the few notable exceptions because the cost of attention grows quadratically with the number of tokens. Even in that case, however, the cost of the linear mappings remains substantial, not mention the potential use of fast approximations of allention (Choromanski et al., 2020).
The computational cost of a linear layer proportional the number of its parameters: at inference, each parameter induces one multiplication per input, where the input is whole input object the case of MI.Ps, a token the case of Recurrent Neural Networks (Elman, 1990) and Transformers (Vaswani et al., 2017), node or edge in the case of Graph Neural Networks (Zhou et al., 2020), a patch of an image in the case of Convolutional Neural Networks (LeCun et al., 2002; Krizhevsky et al., 2012), and similarly for other architectures.
Spline lookup tables make it possible to do better than that. Consider, for example, a one-dimensional serrewise-linea function f(x) on the interval from 0 to 1 with uniform grid. On each interval, it is given as f(z) = afi] x + b(t), where is the interval index. With G intervals, the function has 2G parameters, out of which G G+1 are independent once continuity at the internal grid points is enforced. Yet the computational cost of evaluating such function at any given point is (1), not depending on G. The computational pipeline involves first determining the current grid interval as = [x * GJ, and then evaluating only one linear function as f(x) ali] x + b[i].
Golmogorov-Arnoke Networks (KANs) (Liu et al., 2024), designed as a general alternative to MLPs, are natural hosts for spline lookup tables as they construct general high- dimensional mapping through collection of trainable univariate functions. This design enables the introduction of numerous trainable parameters within each inner function without commensurate increase in inference cost.
The main contributions of this work are the following:
We propose lookup multivariate Kolmogorov- Networks (ImKANs) that are built upon multivariate din-einnsional functions instead of the univariate ones that standard KANs employ. We find that the multivariate version is better suited for a rich parametrization of inner functions in practice. When empirically comparing the 2D version of ImKANs against 1) ID FastKAN (Li, 2024b) and 2) the ID version of ImKANs, we find that the 2D ImKANs are more accurate and easier to train in both cases.
We implement the inner functions as spline lookup tables. Ignoring a non-asymptotic ((N) term, the required inference FLOPs are only 2x those of a lincar layer of the same shape, while the number of trainable parameters can be dozens or hundreds of times higher.
We provide custom CUDA kernels that enable efficient inference of ImKANs modern GPUs When using the 8 8 tile size, on the H100 GPU, our implementation enables up 88x faster inference compared linear layer with the same number of trainable parameters.
We empirically compare ImKANs and MLPs across diverse datasets, scales, and backbones, using varied experimental setups obtain a comprchensive view of performance. Across these conditions, ImKANs are consistently Pareto-optimal with respect to inference FLOPs. The performance of ImKANs is summarized in Fig,
RELATED WORK
The Kolmogorov-Arnol Representation Theorem (KART) (Kolmogorov; Arnold, 2009)
states that continuous function f: 10m -> IR be represented
2

---

Published as a conference paper at ICLR 2026
2n+1
f(r1...., In) =
Σ )am(op))
(1)
where 6op: : :0.1] R. and 0. R RR are contimuous univariate functions. There has been a long debate (Girosi & Poggio, 1989; Schmidt--Hieber, 2021) on the usefulness of this theorem for machine learning because of the general non-smoothness and wild behavior of the inner functions. Nevertheless, it inspired the construction of Kolmogorov-Armole Networks (Hecht- Nielsen, 1987; Igelnik & Parikh, 2003), whose layers are defined as B/ Ep fup(zp), where fap are trainable functions. Liu et al. (2024) introduced the modern version, which suggests stacking an arbitrarily large number of KAN layers and using an arbitrarily large number of neurons, similarly to MLPs. While Liu et al. (2024) illustrated strong performance of KANs, many test cases involve ground-truth functions with known, reasonably smooth KART or KART-like (matching KANs with more than one hidden layer and a larger number of neurons) closed-form representations. Subsequent works such as Yang & Wang (2024), Kundu et al. (2024), and Kashefi (2025) further reinforced the efficiency of KANs. On the contrary, Yu et al. (2024) found that KANs can fall short compared to MLPs for some tasks.
The idea of lookup-based ((1) computations of KAN univariate functions is sometimes briefly mentioned but rarely implemented in practice (Somvanshi al., 2024; Ji et al., 2024), likely because of challenges associated with an efficient GPU implementation, Surprisingly, most of the research goes in the somewhat opposite direction. B-splines, piecewise polynomial basis functions used in the original KAN paper, have compact support and thus are well suited for 011 inference. Subsequent works often replace them with dense basis functions, such as Chebyshev polynomials (SS et al., 2024) Fourier harmonics (Xu et al., 2024). The case of FastKAN (Li, 2024b), which replaces sparse B-splines with looking dense Gaussian radial basis functions exclusively for the sake of optimization, is especially notable.
few works, such as Moradzadeh et al. (2024) and Huang el al. (2025), implement the lookup idea. Moradzadeh et al. (2024), however, predict B-spline coefficients using MLP for the given grid points, which, thus, are not fully independent of each other. Additionally, it benchmarks inference cost against a naive implementation of KANs and not against MLPs, thus not claiming that their construction is Pareto optimal compared to MLPs, which are still seate-of-the-arr for general tasks. Huang al. (2025) achieves remarkable efficiency on co-design using
a small-scale problem from the original KAN paper by algorithm-hardwaar the TSMC 22 nm RRAM-ACIM chip. Polucktov & Polar (2025) and Polar & Polucktov (2021) employ piecewise linear parametrization suitable for O(1) inference but do not focus on inference efficiency. While Shtoff (2025) provides a pure PyTorch implementation of splined parametric curves and, in particular, the KAN architectures derived from them, it does not benefit from low-level CUDA optimizations as the present work.
In this work, we provide CUDA kernels for efficient inference and benchmark the introduced models against MLPs on general tasks where KART representations are not known in closed form and where there is no reason to believe that they are smoother than in other cases.
LOOKUP MULTIVARIATE KOLMOGOROV
NETWORKS
At first glance, given that the inference cost. of spline lookup tables does not depend on the number of parameters, very expressive univariate functions with tens of thousands of trainable parameters each are an ideal match for the Kolmogorov- Arnold Representation Theorem. However, KANs rarely use more than a few dozen parameters per function in practice. The difference between univariate function parametrized by tens of thousands of parameters and just a few dozen the capability of the former to parametrize a very high frequency band. On the one hand, this expressivity is necessary for closely approximating the 'wild behavior' of KART inner functions, but on the other, it raises concerns about training stability and generalization. On the contrary, multivariate functions can 'accommodate" significantly larger number of parameters without spilling expressive power to exceedingly high frequency bands. For instance, four-dimensional function with just 10 grid points
3

---

Published as a conference paper at ICLR 2026
along each dimension has roughly the same number of trainable parameters univariate one with ~ 10 grid intervals.
A layer of a multivariate version of Kolmogorov-Ariol Networks with dimension d defines the output as:
Nin/d-1
Lap Edpt Idp+d-1),
p=0
where fop are trainable d-dimensional functions and Na is the input dimensionality (assumed to be divisible by d). We implement CUDA kernels for the one- and two-dimensional cases. The motivation behind this choice is detailed in Section 3.3. An example of such layer in the 2D case is depicted in Fig. 2. Similar to KANs, these layers do not need additional activations in between and can be stacked arbitrarily, substituting linear mapping-activation pairs in MLP-based backbones.
In Sec. 4.4, empirically compare the two-dimensional version of ImKANs with both one-dimensional ImKAN and me-dimensional FastKAN, The outcomes of these numerical experiments reinforce the intuitive considerations given here and suggest that multidimen- sional building blocks can indeed be more effective hosts for a large number of parameters in a practical setup.
Additionally, it is worth noting that, if necessary, multivariate functions fap can always fall back sums of univariate ones, which would make the whole ImKAN fall back to standard KAN. Thus, the Kolmogorov-Arool Representation Theorem is applicable also to our construction.
3/
Soal
far (r
JT1
X/1
fix (c.))
M/
f
(rx(...)
33
Figure 2: Schematic representation of 2D ImKAN layer with 4 inputs and outputs. This layer defines outputs as: yo = foo (Zo, 21) + for (I2, 1z), V/1 = fio 21) + f12 rs), and y/2 f20(00. 21) S21 (x2, za). The functions f..(:.) are to be trained during fitting.
3.1 FUNCTION PARAMETRIZATION
During training, activations of neurons can evolve arbitrarily, making the use of grids defined on bounded regions challenging. Therefore, we designed an unbounded grid which is still regular enough to allow ((1) computations.
Sigma grid The oue-dimensional sigma grid, which illustrated in the left panel of Fig. is generated by any sigmoid-like function (z). If the desired number of grid intervals is G, then the grid points are given as the intersections of G equispaced percentile levels with (xr) Such construction spans the entire real axis. The grid has the finest. resolution near the origin, and becomes progressively coarser as | increases. For given v, the index of the corresponding grid interval can be computed as [a(x)G], which makes such grid suitable for O(1) computations. For multivariate functions, we apply the same one-dimensional grid independently to cach coordinate.
To balance occupancy across intervals, we precede each ImKAN layer with a batch normal- ization layer without affine parameters to keep the activations in controlled range. The choice of the (2x) function was primarily motivated by computational efficiency. The exact definition and other details available in Appendix B.
4

---

Published as a conference paper at ICLR 2026
1.25
1.0
(
1.00
lines
Grid points
0.75
Grid points
0.50
0.25
-4
-3
-2
Figure 3: (Left) Construction of the sigma grid; (Right) Second-order B-spline.
Compact basis functions We use B-splines of second order built on top of the described grid as basis functions to parametrize the ImKAN inner functions. A one-dimensional sccond-order B-spline centered around grid point i is given in the right panel of Fig. 3. takes non-zero values on only the two adjacent grid intervals around the center grid point. II there are G intervals, then we G 1 basis functions: G 1 such B-splines centered around each inner grid point, and two linear functions on the leftmost. and rightmost infinite intervals. Appendix B provides the definition of the edge-interval B- splines used in this work, along with other details. A two-dimensional B-spline is defined as Ba (I1,
Bi (z1) Brg (I2). All the functions defining 2D ImKAN layer are parametrized as f(11, 22) En,,, Piç Bi (21,x2) where Piú are trainable coefficients.
With such construction, there (G + 1)2 independent parameters for each 2D function, parametrized functions are bilincar on each 2D grid interval, all the functions are continuous for arbitrary Piçéz and all but the edge coefficients Piréa have a simple interpretation the value of the function on the corresponding grid point.
For any given point (2), there only four non-zero B-splines; thus, one needs to evaluate only four terms to compute ((x), which forms the basis of O(1) computations. The full algorithm compute a standalone two-dimensional function given in Appendix B.1.
3.2 COMPUTATIONAL COST
Overall, the functional form of ImKANs involves computations of many low-dimensional functions with exactly the same arguments (those in the same column, see Fig. 2). When doing it is possible to reuse many of the intermediate values, such as indices of the grid intervals and B-splines. These intermediate values can be computed once for each pair of inputs, and then utilized to compute the value of given function with just four multiply-add operations for the 2D case. Given that the total mumber of 2D functions in ImKAN layer is [N /2) Nout, the total number of required multiply-add operations for the dominant, O(N2), part is 4[Nmm/2 Nout = 2N.Noott just 2x that of a linear layer of the same shape. Algorithm 3, given in Appendix C, pinpoints this estimation. Following the common practice (He et al., 2016), we estimate FLOPs as the number of fused multiply-adde of the main asymptotic term for both MLPs and ImKANs.
It is worth noting that the non-asymptotic O(N) term is not an additional cost relative to MLPs; replaces the per-unit bias additions and activation evaluations that ImKANs do not require - operations that can be quite expensive when the activations are transcendental functions such as tanh.
3.3 PERSPECTIVE DIMENSIONS AND SPLINE ORDERS
All the constructions introduced so far straightforwardly generalize to higheer-dimensional case. Evaluation single d-dimensional function parametrized by d-dimensional second- order B-splines would take 2d multiply-adds, which stem from the summation of B-spline contributions residing on all the corners the corresponding d-dimensional hypercube. Given that the number of such d-dimensional functions would be d times smaller than the number

---

Published as a conference paper at ICLR 2026
of weights of linear layer of the same shape, the inference FLOPs of the d-dimensional ImKAN layer would be /d of that the linear layer with the same shape.
These slowdown factors identical, 2x, for to grow for higher dimensions. Thus, we chose the one- and wwo-dimensional versions for the implementation of CUDA kernels the 2D extension of the standard univariate KAN comes essentially for free.
and two-dimensional ImKANS, and start
If B-splines of order & are employed for paramotrization instead of the second-order ones described so far, the inference cost becomes k- /d of that of an MLP. For more details, see the B-spline definition at De Boor & De Boor (1978), and how they used in KANs (Liu et al., 2024). Increasing the B-spline order brings a few benefits, but they likely do not justify the associated increase in the computational cost.
The classical theorem about B-splines (De Boor, 1968) indicates that while the order of B- splines affects the convergence rate, any spline order is sufficient to approximate any function arbitrarily closely by increasing the resolution of the grid'. Given that the computational cost of spline look-up tables does not depend on the number of grid points, the grid resolution can be arbitrarily increased without any computational overhead at inference.
On top of that, increasing the B-spline order introduces additional smoothness to the model. Functions parametrized by B-splines of order k belong to C*-2, but in general not to Ck-11 The smoothness of second-order B-splines employed in this work matches that of ReLU (Glorot et al., 2011), one of the most popular and successful activation functions, which also continuous, but not contimuuully differentiable. Thus, it is questionable if the additional smoothness available through higher orders & is necessary and would justify the associated computational overhead.
3.4 CUDA KERNELS
We implement CUDA kernels for efficient inference of 1D and 2D ImKANs on modern GPUs. Our kernels use the classic shared memory tiling used in GEMM (Volkov & Demmel, 2008). In the following, all benchmarks run in full-precision float32.
In the 2D case, with a 16 x 16 tile, our implementation is ~ 8x slower than a dense linear layer with the same shape on an H100-SXM, irrespective of the grid resolution. The slowdown exceeds the 2x FLOPs-based estimate from Sec. 3.2, primarily because of the less coherent memory-accesss pattern algorithm 3. Additionally, dense matrix multiplication has been the cornerstone of many computational pipelines and, thus, has enjoyed decades of thorough optimization.
Finite shared memory capacity limits the number of grid intervals to G < 20 on H100. At this limit, an ImKAN layer holds (20 + 1)3/2 ss 220x more parameters than the linear baseline with the same shape. This means that an ImKAN layer delivers [20 + 1)2/2) x 27.5x faster inference compared to linear layer with the same number of trainable parameters.
Reducing the tile to x 8 raises the slowdown to ~ 9.5x but lets us increase to 40, yielding 88.5x better per-parameter efficiency.
4 RESULTS
We have demonstrated so far that ImKANs can have significantly better inference cost per trainable parameter compared to linear layers in terms of both FLOPs and wall-clock time modern GPUs. The question is, however, whether this nominal efficiency translates to real-life performance. Do ImKANs indeod represent a better tradcoff between performance and inference cost?
In this section, we empirically compare the efficiency of ImKANs and MLPs across the following settings: (i) approximating general high-dinensional functions, (ii) tabular-like dataset of randomly displaced methane configurations, and (iii) within CNN frameworks
Thhis applicable, though, functions bounded domain.

---

Published as a conference paper at ICLR 2026
evaluated on CIFAR-10 and ImageNet. Across all experiments, use identical macro- architectural backbones for ImKANs and MLPs. Overall, to obtain a comprehensive picture very large scale the
of the performance, we prioritized the diversity of the setups architectural complexity of a particular backbone. We found that ImKANs are consistently Pareto-optimal, with the largest gains the methane dataset. Finally, we compare ImKANs with FastKANs.
inference
4.1 GENERAL PUNCTION APPROXIMATION
Our first experiment is set to measure the crude flexibility of ImKANs in approximating general high-dimensional functions, which we model by large teacher MLPs with fixed random weights. We define ground-truth TR32 Rd function hidden layers, each with 1024 neurons, and hyperbolic tangent activations.
ML.P with 32 input neurons, 10
We fit both MLP and ImKAN students to approximate this ground-truth function and compare their performance. We use the same fully connected backbone for both types of models with two hidden layers and varying hidden dimensions and grid resolutions G (for more details, see Appendix G.2).
In our first experiment, we fixed the hidden dimension at 256 to investigate the performance of ID and 2D ImKANs as a function of grid resolution G. The results, presented in Fig. 4 (left panel), show that performance does not improve indefinitely with higher resolution; instcad, it follows U-shaped curve. However, there is a notable distinction between the one- and two-dimensional cases: the 2D ImKAN saturates a much larger number of trainable parameters per function and achieves significantly higher accuracy. Notably, the 2D lnKAN with optimal grid resolution performs comparably to an MLP that roughly 16 times larger (with 4x larger hidden dimension). Furthermore, the second and third panels confirm that inference cost in terms of both FLOPs and H100 SXM wall-clock time remains independent of the grid resolution G.
Maan Squueed Ereer
0.34
ImKAN 2D
InKK4N 1D
Inforesoo sze TLOPs
20
InKAN baseliee
ImKAN 20
InKAN
InKAN 2D
Vittr dr.
0.52
HiP beseline
HLP baseine
0.30
0.40
10-1
0.46
Number paremeters frear funation
Number garamelers tsmer functiion (G ImICAN 27 InikAN 20
Number of panameters per inter fameliom
ImKAN + 117
In1AN 2D
12 1) for ImKKN
Figure 4: Performance of 1D and 2D ImKANS depending on the grid resolution.
In our second experiment, we fixed the grid resolution its optimal value and varied the hidden dimensions for both MLPs and 2D ImKANs. As expected, the final Mean Squared Error (MSE) decreases monotonically with the hidden dimension, reflecting the increased flexibility of the models. However, larger hidden dimensions inevitably incur higher inference costs; the key question, therefore, is which architecture offers a better tradcoff between accuracy and efficiency. Our findings are summarized in Fig. 5. The left panel plots the final MSE against the hidden dimension, while the middle and right panels depict the Pareto fronts for MSE versus inference FLOPs and H100 wall-clock time, respectively.
In order to justify the claims that ImKANs are indeed more efficient, we converge the baseline MLP-based models very tightly here, and in all similar experiments in this manuscript. For the MLP baseline, we have two lines: one with a full training budget and one with only half of it. The close similarity between them demonstrates very tight convergence.
Fig. 5 clearly indicates that ImKANs are significantly more FLOP-efficient at the same accuracy level, up to 6x, for the largest dimensions. Furthermore, ImKANs also appeared to be H100 wall-clock time optimal for all the scales, with a speedup factor of about 1,8x for the largest hidden dimension.

---

Published as a conference paper at ICLR 2026
MSE vs Hidden dimension
MSE vs Inference FLOPs
MSE vs Inverse H100 SXM throughput
KCr 122 MLP Im0
K2
1/2
Mean Squared Errer
MLP
MiP
kmW
6 xx
5 xx10-1
512
1024
10a
10a
10a
10
10-3 size)
Hidden dimension
Figure ImKAN vs MLP for general function approximation. The 'MLP 1/2" line corresponds to the outcome of the fitting procodure with only half of the training steps compared to the 'MLP' onc.
Appendix G.2 contains analogous experiment for an IR + 32 function, ablation studies, and other details.
4.2 RANDOMLY DISPLACED METHANE CONFIGURATIONS
Having demonstrated that ImKANs are Pareto-optimal when ap- proximating a general function, we proceed benchmark their efficiency on real data. We chose the tabular-like dataset of ran- domly displaced methane configurations for this comparison, it is particularly suitable for this purpose (sce Appendix G.4). The dataset consists of multiple off-equilibrium methane configurations, as illustrated in Fig, 6. The larget is given by the corresponding quantum-mechancen
energy (Turney al., 2012; Kohn Sham, 1965). Hydrogen atoms are placed around the carbon atom ran- domly, varying from instance to instance, which leads to different target energies.
We encode the geometry of each methanc molecule by the Cartesian components of displacement vectors from the carbon atom to all the hydrogen atoms. Therefore, the input dimension of both the the
Figure A methane configuration
ImKAN and the ML.P networks is 12. Otherwise, backbone with two hidden layers as in the previous experiment.
Relative RMSE vs Hidden dimension
Relative RMSE FLOPs
Relotive RMSE vs inverse H100 SXM throughput 10"
10a
10
The SAEE stD
ImCAN
InKAN
ImKAN
10-1
10-1
78.00
12, Bx
302
30
10
124
13
10
10-1
Inforence
Midden dimmension
Herre beitch siee
Figure 7: ImKAN vs MLP the dataset of randomly displaced methane configurations. On the vertical axis, we plot the relative Root Mean Squared Error, which given as test RMSE normalized by the standard deviation of the target in the dataset. The "MLP 1/2" and 'MLP 1/4" curves correspond to outcomes of fitting procedures with half and quarter of the training budget, respectively.
Our findings are given Fig. 7. Similar to the previous experiment, ImKANs were found to be Pareto-optimal in terms of both formal FLOOPs and H100 wall-clock time. Moreover, the
8

---

Published as a conference paper at ICLR 2026
relative speedup is higher for this dataset: up to 78.0x reduced FLOPs, and 12.9x faster H100 inference at matched accuracy.
In Appendix G.4 we provide ablation studies, analogous experiments using representations other than the Cartesian components the displacement vectors, and further details.
4.3 IMKAN-BASEI CONVOLUTIONAL NEURAL NETWORKS
In the introduction, we brielly outlined that high-dimensionall linear mappings are the primary building blocks most architectures, not only leedforward fully connected neural networks. Convolutional Neural Networks (CNNs) are exception.
standard two-dimensional convolution with kernel size & x parametrized by a linear mapping Ci, + RCo,, where Cin and Cout are the numbers of input and output channels, respectively. Since Kolmogorov-Aroold layers can be used as general substitute for high- dimensional linear mappings, one can construct a KAN-based convolut ional neural network well suited for image processing, as was done, e.g., in Bodner et al. (2024). In this section, we compare the performance of ImKAN- and MLP-based CNNs on the CIFAR-10 10 (Krizhevsky et al., 2009) and ImageNet (Deng et al., 2009) datasets.
For CIFAR-10, our backbone architecture consists of five 2 x 2 convolutions, each with stride 2, and two fully connected layers at the end. We use identical backbones for ImKAN and MLP CNNs, and exactly the same pool of augmentations (see more details in Appendix G.5). Similarly to the previous experiments, we vary the width of the neural networks, which the number of filters in the case of convolutions, and the hidden dimension in the case of fully connected layers.
For ImageNet, we downsample the images to 81 81 pixels. Our backbone consists of four convolutional layers each with a 3 x kernel and a stride of 3 followed by two fully connected layers. In contrast to the CIFAR-10 experiment, we progressively increase the number of filters as the spatial resolution of the image decreases through the neural network.
Accuracy Hidden dimension
Accuracy Inference FLOPS
10
Test Accuracy (16)
CIFAR 10
70
70
60
CNN
ENN
MLP CMN
InKON ChN
50
ImKAN CNN
Width
Inference FLOPS
Accuracy Base width
Accuracy Inference FLOPs
Vel Top-5 Accuracy (16)
175
55
ImageNet
50
1.7
45
ImKIAN
ImKAN C3N
MLP CNN
40
MLP CNN
MUP CNN 1/2
YUP CNN 1/2
35
16
Base width
Inference FLOPs
Figure 8: Comparison the performance of standard MLP-based CNNs and ImKAN-based CNNs the CIFAR-10 and ImageNet datasets. The 'MLP CNN 1/2" line corresponds to the outcome of the fitting procedure with only half of the training steps compared to the "MLP CNN' one.
Our findings are illustrated in Fig. Similarly to previous experiments, LmKKNNbbee based CNNs were found to be more FLOPs efficient at the same accuracy level. The observed reduction in FLOPs 1.6-2.1x for CIFAR-10, and 1.7x for ImageNet. We have not yet implemented dedicated CUDA kernels for efficient inference of ImKAN-basee bassedd convolutions. Any type of convolution can be cast to a fully connected layer by the corresponding memory manipulations, which employed during fitting.
9

---

Published as a conference paper at ICLR 2026
4.4 COMPARISON WITH FASTKAN
In this section, we compare the performance 2D ImKAN, ID ImKAN, and FastKAN models. We the training script² for the CIFAR-10 dataset. available in the FastKAN GitHub repository (Li, 2024a) as the basis for the comparison. This script fits fully connected FastKAN model with one hidden layer. We provide several modifications, particularly enabling the same augmentation pipeline we used in the previous section (Appendix G.6).
We fix the hidden dimensionality 256 as in the original script for all the models, vary grid resolutions, and analyze the accuracy depending on the number of trainable parameters in each inner function. The result given Fig.
Kn (23)
oporss 32A moorss
66
VHCAA (12)
TeaSSAN (20)
44-
Figure 9: Comparison of the one- and two-dimensional ImKAN and FastKAN models within the fully connected framework on the CIFAR-10 dataset for fitting budgets of one and ten thousand epochs. The number of parameters per inner function is (G + 1)² for 2D ImKAN, (G + 1) for 1D ImKAN, and (num grids 1) for FastKAN.
The first observation is that the performance of both FastKAN and 1D ImKAN models degrades for excessively fine grid resolutions. For the training budget of one thousand epochs, the final models become even less accurate than the ML.P baseline. For ten thousand epochs, the effect is less pronounced, but still takes place. This degradation if present at all, for the 2D ImKAN models- tooyls--heyy-eoo they are easier to fit even for a very rich parametrization of each inner function. Another distinct feature of Fig. 9 that 2D ImKANs achieve notably better accuracy compared to both 1D ImKANs and FastKANs, even when the latter have a very rich paramctrization of inner functions.
much less severe,
To sum up, these findings reinforce the inluitive considerations given in Sec. 3 and suggest that building blocks of multivariate trainable functions are indeed more effective.
SUMMARY
Higl-dinensional linear mappings are the cornerstone of modern deep learning, dominating both the parameter count and the computational cost in most models. We introduce lookup multivariate Kolmogorov-Arnold Networks (ImKANs) that offer substantially better capacity-inference inference cost ratio. Across all experiments, ImKANs were found to be inference FLOPs Pareto-optimal. The efficiency gains are task-dependent: for general high-dimensional function approximation, modelled as a distillation from a large ground-truth teacher MLP with random weights, ImKANs achieved up to 6x fewer FLOPs at matched accuracy. On randomly displaced methane configurations efficiency improved by to 78x. Within convolutional networks, the gains were smaller but still significant: 1.6-2.1x on CIFAR-10 and 1.7x on ImageNet.
Our CUDA kernels compete directly with highly optimized dense matrix multiplications the backbone of many numerical pipelines for decades. Even so, the gains were sufficient to H100 wall-clock time for both the general high-dimensional function approximation and the methane datasct, achieving speedup of more than an order of magnitude in the latter case.
make ImKANs Parcto-optimal
10

---

Published as a conference paper at ICLR 2026
6 REPRODUCIBILITY STATEMENT
The associated supplementary material contains not only the implementation of ImKANs, but also training scripts for the experiments presented in this manuscript, along with the corresponding settings as a collection of YAML files. These asscts enable reproducibility of the numerical experiments presented in this paper.
ETHICS STATEMENT
This work adheres to the ethical standards expected at ICLR 2026. To the best of our knowledge, the proposed method is not expected to have adverse effects.
REFERENCES
Alice EA Allen, Geneviève Dusson, Christoph Ortner, and Gábor Csányi. Atomic permuta- tionally invariant polynomials for fitting molecular force fields. Machine Learning: Science and Technology, 2(2):025017, 2021.
Vladimir Arnold. On functions of three variables. Collected Works: Representations of Functions, Celestial Mechanics and KAM Theory, 1957-1965, pp. 5-8, 2009.
Albert P Bartók, Mike C Payne, Risi Kondor, and Gábor Csányi Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons. Physical review letters, 104(13) 2010.
Jörg Behler and Michele Parrinello. Generalized neural-network representation of high-
dimensional potential-cnergy surfaces. Physical review letters, 98(14):146401, 2007.
Alexander Dylan Bodner, Antonio Santiago Tepsich, Jack Natan Spolski, and Santiago Pourteau. Convolutional kolmogorov-arrold networks. arXiv preprint arXiv:2406.1315. 2024.
Krzysalof Choromanski, Valerii Likhosherstov, David Dohan, Xingyou Song, Andreea Gane, Tamas Sarlos, Peter Hawkins, Jared Davis, Afroz Mohiuddin, Lukasz Kaiser, et al. Re- thinking attention with performers. arXiv preprint arXi::2009. 14794, 2020.
Ekin D Cubuk, Barret Zoph, Jonathon Shlens, and Quoe V Le. Randaugment: Practical au- lomated data augmentation with reduced search space. In Proceedings of the IEEE/CVP conference on compuler vision and pallern recognilion workshops, pp. 702-703) 2020.
George Cybenko. Approximation by superpositions of a sigmoidal function. Mathematics of control, signals and systems, 2(4):303 -14, 1989.
Carl De Boor. On uniform approximation by splines. Approz. Theory, 1(1):219 235, 1968.
Carl De Boor and Carl De Boor. A practical guide splines, volume 27. springer New York, 1978.
Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei. Imagenet: A large-scale hierarchical image database. 2009 IEEE conference on computer vision and pattern recognition, pp. 248 255. Ieee, 2009.
Harm Derksen and Gregor Kemper. Compulational invoriont theory. Springer, 2015.
Jeffrey L Elman. Finding structure in time. Cognitive science, 14(2):179 211, 1990.
Federico Girosi and Tomaso Poggio. Representation properties of networks: Kolmogorov's theorem is irrelevant. Neural Computation, (4):4465-46, 1989.
Xavier Glorot, Antoine Bordes, and Yoshua Bengio. Deep sparse rectifier neural networks. In Proceedings of the fourteenth international conference on artificial intelligence and statistics, pp. 315 323..23.442 JMLR Workshop and Conference Proceedings, 2011.
11

---

Published as a conference paper at ICLR 2026
David Harrison Jr and Daniel L Rubinfeld. Hedonic housing prices and the demand for clean air. Journal of environmental economics and management, 5(1):81-102, 1978.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pallern recognilioo,, pp. 770 778, 2016.
Robert Hecht-Niclsen. Kolmogorov's mapping neural network existence theorem. In Proceed- ings of the international conference on Neural Networks, volume 3, pp. 11-14. IEEE press New York, NY, USA, 1987.
Geoffrey E Hinton, Nitish Srivastava, Alex Krizhevsky, llya Sutskever, and Ruslan R Salakhutdinov. Improving neural networks by preventing co-adaptation of feature detectors. arXiv preprint arXiv:1207.00 1207.0580, 2012.
Wei-Hsing Huang, Jianwei Jía, Yuyao Kong, Faaiq Waqar, Tai-Hao Wen, Meng-Fan Chang, and Shimeng Yu. Hardware acceleration of kolmogorov-arnold network (kan) for lightweight edge inference. In Proceedings of the 30th Asia and South Pacific Design Automation 2025.
Conference, pp.
Boris Igelnik and Neel Parikh. Kolmogorov's spline network. IEEE transactions neural nehworks, 14(4):725-73) 733, 2003.
Sergey Ioffe and Christian Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. International conference on machine learning, pp. 448-456. pmlr, 2015.
Tianrui Ji, Yuntian Hou, and Di Zhang. A comprehensive survey on kolmogorov arnold networks (kan). arxiv preprint orXiv:2407.11 11075, 2024.
Kaggle. Titanic machine learning from disaster. https://www.kagee wwt.://wwwwwag..oooe/.//.. com/ competitions/titane Accessed 2025-08-19.
Jared Kaplan, Sam McCandlish, Tom Henighan, Tom B Brown, Benjamin Chess, Rewon Child, Scott Gray, Alec Radford, Jelsrey Wu, and Dario Amodei. Scaling laws for neural language models. arXiv preprinl or Xiv:2001.08300 2020.
Ali Kashefi. Pointnet with kan versus pointnet with mlp for 3d classification and segmentation of point sets. Computers & Graphics, pp. 104319, 2025.
Walter Kohn and Lu Jeu Sham. Self-consistent equations including exchange and correlation effects. Physical review, 140(4A):A1133, 1965.
Andreí Kolmogorov. On the representation of continuous functions of several variables by superpositions of continuous functions of a smaller number of variables.
Alex Krizhevsky, Geoffrey Hinton, et al. Learning multiple layers of features from tiny images. 2009.
Alex Krizhevsky, Πya Sutskever, and Geoffrey E Hinton. Imagenet classification with deep convolutional neural networks. Advances in neural information processing systems, 25, 2012.
Akash Kundu, Aritra Sarkar, and Abhishek Sadhu, Kanqas: Kolmogorov-arrold network for quantum architecture scarch. EPJ Quantum Technology, 11(1):76, 2024.
Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEFE, 86(11):2278-23)) 2002.
Ziyao Li. fast-kan: Fastkan very fast implementation of kolmogorov arnold networks (kan). https://github. com//iyyoliifastttaaaaaaoaaa 2024a. GitHuu repository.
Ziyao Li. Kolmogorov-arnold networks radial basis function networks. 2024b.
12

---

Published as a conference paper at ICLR 2026
Ziming Liu, Yixuan Wang, Sachin Vaidya, Fabian Ruchle, James Halverson, Marin Soljačié, Thomas Y Hou, and Max Tegmark. Kan: Kolmogorov-arnole networks. arXiv preprint ar 2rXXt:2404.19956 2024.
llya Loshchilov and Frank Hutter. Sgdr: Stochastic gradient descent with warm restarts. Xive preprint or arXiv: 1608.600.0000 160800383388, 2016.
Alireza Moradzadch, Lukasz Wawrzyniak, Miles Macklin, and Sace G Paliwal. Ukan: Unbound kolmogorov-arnold network accompanied with accelerated library. arXiv preprint ar08::2408. 1 1200, 2024.
Andrew Polar and Michael Poluektov. A deep machine learning algorithm for construction of the kolmogorov arnold representation. Engineering Applications of Artificial Intelligence, 99:104137, 2021.
Michael Poluektov and Andrew Polar. Construction of the kolmogorov-arnold networks using the newton-kaczmarz method. Machine Learning, 114(8):185, 2025.
Sergey Pozdnyakov and Michele Ceriotti. Smooth, exact rotational symmetrization for deep learning on point clouds. Advances in Neural Information Processing Systems, 36: 79469-7999 19469, 2023.
Johannes Schmidt-Hiebe Hieber. The kolmogorov arnold representation theorem revisited. Newral nelworks, 137:119 126, 2021.
Alex Shtoff. torchcurves: Differentiable parametric curves in pytorch, 2025. URL https: //github. h/w.ttub.ccmm//aeeeeeeeeeeeeeeet htmm//llexstf/toocccaeeee
Shriyank Somvanshi, Syed Aaqib Javed, Md Monzurul Islam, Diwas Pandit, and Subasish Das. A survey on kolmogorov-arnold network. ACM Computing Surveys, 2024.
Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, llya Sutskever, and Ruslan Salakhutdi- nov. Dropout: a simple way to prevent neural networks from overfitting. The journal of machine learning research, 15(1):1929-1958 2014.
Sidharth SS, Keerthana AR, Anas KP, et al. Chebyshev polynomial- kolnogorov-arnold networks: An efficient architecture for nonlinear function approximation. ar preprint arxiv:2405.00 07200,, 2024.
Justin M Turney, Andrew C Simmoneut, Robert M Parrish, Edward G Hohenstein, Francesco A Evangelista, Justin T Fermann, Benjamin J Mintz, Lori A Burns, Jeremiah Wilke, Micah L Abrams, et al. Psi4: an open-source ab initio electronic structure program. Wiley Interdisciplinary Reviews: Computational Molecular Science, 2(4):556- 565, 2012.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. Advances in neural information processing systems, 30, 2017.
Vasily Volkov and James W Demmel. Benchmarking gpus to tune dense linear algebra. In SC'08: Proceedings of the 2008 ACM/IEEE conference on Supercomputing, pp. 1-11. IEEE, 2008.
Stephen R Xie, Matthias Rupp, and Richard G Hennig. Ultra-fast interpretable machine- learning potentials. upj Compulational Materials, 9(1):162, 2023.
Jinfeng Xu, Zheyu Chen, Jinze Li, Shuo Yang, Wei Wang, Xiping Hu, and Edith C-H Ngai. Fourierkan-gef: Fourier kolmogorov-arnold network-an effective and efficient fcature transformation for graph collaborative filtering, arXiv preprint ar arXiv:2406.0100 01044, 2024.
Xingyi Yang and Xinchao Wang. Kolmogorov-armold transformer. arXiv:2/09. 10594, 2024.
arXiv preprint
Runpeng Yu, Weihao Yu, and Xinchao Wang. Kan or mlp: A fairer comparison. arXiv preprint 16767 2024.
13

---

Published as a conference paper at ICLR 2026
Sangdoo Yun, Dongyoon Han, Seong Joon Oh, Sanghyuk Chun, Junsuk Choe, and Youngjoon Yoo. Cutmix: Regularization strategy to train strong classifiers with localizable features. In Proceedings of the IEEE/CVF international conference computer vision, pp. 6023-6032, 2019.
Xiaohua Zhai, Alexander Kolesnikov, Neil Houlsby, and Lucas Beyer. Scaling vision trans- formers. In Proceedings of the IEEE/CVE conference on compuler vision and patlern recognition, pp. 12104 12104-12113, 2022.
Hongyi Zhang, Moustapha Cisse, Yann N Dauphin, and David Lopez-Paz. mixup: Beyond empirical risk minimization. orXiv preprint arXio: 1rx.::1710.0941,, 2017.
Yaolong Zhang, Junfan Xia, and Bin Jiang. Physically motivated recursively embedded atom neural networks: incorporating local completeness and nonlocality. Physical Review Letters, 127(15):156002, 2021.
Jie Zhou, Ganqu Cui, Shengding Hu, Zhengyan Zhang, Cheng Yang, Zhiyuan Liu, Lifeng Wang, Changeheng Li, and Maosong Sun. Graph neural networks: A review of methods and applications. Al open, 1:57-81, 81, 2020.
A LLM USAGE STATEMENT
We used LLMSs for (i) polishing the text of the manuscript; (ii) finding related work and locating specific details in the corresponding papers; (iii) as a coding assistant during implementation of our codebase.
B MORE DETAILS ON FUNCTION PARAMETRIZATION
The (x) function Section 3.1 of the main text describes the construction of the sigma grid used in the parametrization of all the functions ImKAN consists of. involves any sigmoid-like function (x), and, as mentioned in the main text, choice, which is given in Eq. 3 and illustrated in Fig. 10, was motivated by computational efficioncy.
1.0
0.8
0.6
e
x 0,
(x))
(3)
0.4
0.5 -2 0.
0.2
-4
-2
Figure 10: The (zx) function defined in
Eq. 3.
This construction cheap to compute because the computational pipeline consists of a single exponential call and few arithmetic operations, elaborated in algorithm 1.
14

---

Published as a conference paper at ICLR 2026
Algorithm 1 Evaluation of (z) with single exponential call
Input: x E R
Output: o(x)
1: 4- exp [x])
compute expensive exponential only once
2: if then
3:
(zz) 1- 0.5
4: else
5:
((:) 0.5t
6: end if
7: return (x))
Edge cases Section 3.1 the main text introduced sigma grids and the corresponding basis of second-order B-splines. For a grid with G intervals and G grid points, there are G+1 1 basis functions, out of which G - I are given by second-order B-splines centered around all grid points, as illustrated in the right panel of Fig. 3 of the main text. The other two are given as linear functions on the leftmost and rightmost infinite intervals.
The second-order B-splines given in the right panel of Fig. 3 are defined linearly increase from 0 to 1 from the left grid point to the central one and then linearly decrease from 1 to 0 from the central grid point to the right grid point. This construction is well-defined for all the inner grid points but requires additional definitions for the B-splines centered around the leftmost and rightmost grid points, as these do not have left and right neighboring grid points, respectively.
In order define these edge B-splines, we introduce 'ghost' left and right grid points. The position of the left ghost point is given g (g[1] g(0),, where ( and g[1] are the positions leftmost and second leftmost grid points respectively. The right ghost point is defined similarly. With such notation of additional grid points, can define edge B-splines similarly to all the others.
Finally, there two linear basis functions the leftmost and rightmost infinite intervals. We define them to be 0-s in the leftmost and rightmost grid points, linearly increasing to 1-8 the left and right ghost. points and continuing to left and right infinities with the same slope, respectively.
Direct o-grid vs. uniform grid after pre-normalization by (xx) natural question is how the proposed sigma grid construction differs from simpler alternative that first maps the input via x' (()) E (0, 1) and then fits a heeeewse-llnear function g(x') on uniform grid over 1].
In general, the resulting functions in the original z-domain, f(x) and f((x) g()())), are not identical. Most importantly, f'(x) is no longer piecewise linear within each x-grid interval. The most notable discrepancy appears the tails: f'(x) has horizontal (constant) asymptotes when r + too, whereas f(z)-by construction on the exhibits lincar asymptotes.
Horizontal asymptotes entail vanishing gradients in the tails, behavior widely implicated training difficulties and one reason ReLU activations often outperform tanh (Glorot et al., 2011). For these reasons, adopt the direct o-grid parametrization in practice.
B.1 STANDALLONE TWO-DIMENSIONAL FUNCTION COMPUTATION
Algorithun provides a recipe to compute a standalone two-dinensional function given our parametrization scheme. See discussion in Sec. 3.1 of the main text.
C COMPUTATIONAL COST
Algorithm 3 represents a computational pipeline to compute the forward pass an entire 2D ImKAN layer.
15

---

Published as a conference paper at ICLR 2026
Algorithm 2 (1) evaluation of a standalone 2D ImKAN function. Red lines indicate computations that can be reused when computing many 2D iunctions for the same arguments, while green lines indicate computations that have to be repeated for each 2D function.
Input:
scalars 21, 22 R; grid points G; parameters P [[C+1,0+1]] P(á, i2] stores the function value on the (i1, 22) grid point. Output: output y IR
1: function Eval2D (11, Handling of edge-interval cases described above is omitted for brevity.
2: i - G
3: t - G
9 +1 31
g +1 32
4:
x22) +
g + - 9 1
9 +
1
31 1
22
5:
+
g +1|
9 +1
9 12
g +1]+ T1
22 9 i2
6:
82)
g +1
g 11
9 +
g i2
371
02
9
i2
7:
Buttiitl ②2)
g +
11
9 1/
+
1|
9
8:
+ 0
9:
y += Bitlla 11
(2 22
P 11 i2
10:
P 11 i2
y+= += Ba (1 22
11, 12+1
12:
+ Bi
p[ +1
13:
return y
14:
end function
As the algorithm shows, the preamble part contributes only to an asymptotically insignificant O(N) term, where N is the input (Nm) output (Nout) dimension. Given that the total number of 2D functions multiply-add operations for the dominant, O(N2), part is 4[N,m/2] Nout = 2N, Nout, just 2x that of a linear layer of the same shape.
an ImKAN layer is
the total number of required
D HESSIAN REGULARIZATION
Direct fitting of splined functions with fine grid imposes additional challenges related to generalization. The problem is illustrated in Fig. 1la for the simple case of fitting standalone one-dimensional function parametrized by second-order B-splines uniform grid on [0, 1] with G = 40 intervals, f(x) = Σ, (z).
Ground truth is an exact parabola, and the training illustrated on the plot. Since the number of grid points is large, the model has enough flexibility to reproduce the training set exactly, but generalization is quite off. With such fine grid, only a few B-splines, marked as bold on the bottom panel of Fig. 11a, take non-zero values for the training points. Thus, only the coefficients p; associated with these active B-splines receive non-zero gradient during training, while all the others have incentive to evolve from the random values assigned at initialization. Standard 1.2 regularization not much better, it simply pushes all non-active coefficients p: to zero, which also results in non-meaningful approximation after training.
consists of 5 points, which are
When dealing with similar problem, Xie et al. (2023) employed off-diagonal regularization based on a finite-difference scheme for the second derivative. One can put regularization terms X E, p for first derivative, A (p: 2p P242)2 for second, and on. Such regularization schemes result in meaningful approximations after training, as illustrated in Fig. 11b.
16

---

Published as a conference paper at ICLR 2026
Algorithm 3 Forward pass of 2D ImKAN layer.
Input: imput vector RN=,, parameter tensor P IR G+I, G+1, Mujaa, Nout
P( i2, imput inder, output inder] is the value of fimput fimput_inder inden, outpue indez at the
i1= 22-th grid point.
Output: output vector y E R Nom
y+0
2: for inpul indez = 0 to Nu/2 - do
3: 12, Ba (T1, [2), B +1 (1 x2); B i2 (21, B1 12 (2)
Prcamble x imputt index), x in2-innut__mmm inder 1])
4:
for output_ing pepput_inder---ee = 0 to Nout do
5: yloutput_ inder + (z P i i2, input_ inder, output_ index]
6:
oulput index + (z
P
1i+ 12, inpul inder, oulpul index
7: y output index
+ Bu (2 T2)
P 11 12 inpul inder, oulpul inder Plitl, igt1, input indez, output (nde__ndder_
8: yloutput index
B (1
9: end for
10: end for
11: function Preamble(z), x2) Handling of edge-interval cases omitted for brevity.
See Appendix for more details. Input:
æ1, Ø2i precomputed grid points 9;
precomputed inverse areas Ainv (1, 12] = [g + +1 - g( (2 + Output: indices 82 and B-splines B, is (x1, x2), Bi (21, Bi, i2 (x), Big (x,
g(v2))
(x1,
12:
i - G
13:
is + -+ G(12)
14:
(1, x2) t + -
Aiy (1, 22
15:
B2t+1i x2) + (z) [0 22) iny 12
16:
i1 22)
+
(12 iny i2
17:
B1i+1/++100 (2)
(21 9 9 Anv 82
18:
return i1, i2, Bi s (x1, B.) B1 Bi (x1, B (211)
end function
Problem
Solution
1.0
1.0
Ground truth regularisation regularisation
Ground truth (D) regularization
regularization
0.8
"(x) regularization
0.6
Rx)
0.4
0.0
0.2
0.0
0.2
0.6
0.8
10
X
X
Figure 11: Generalization pitfall.
We implemented CUDA kernels for 2D ImKANs, which express general high-dimensional mapping in terms of building blocks two-dimensional functions. For 2D functions, we off-diagonal regularization based on the squared Frobenius norm of the Hessian,
17

---

Published as a conference paper at ICLR 2026
rotationally invariant measure of curvature in any direction, which is not to be confused with the Laplacian. Furthermore, our finite-difeerence schemes take into account that the grids, introduced in Sec. 3.1, are not uniform. The exact formulas given below.
The Hessian of a function is zero if and only the function linear. Therefore, the use of very strong Hossian-based regularization leads to lincarization of the trained functions, enforcing them to converge to [()) = azı + bxg + This makes the whole ImKAN equivalent an MLP of the same shape, modulo training dynamics. In other words, the Hessian regularization coefficient À can be used to smoothly adjust the ImKAN behavior between fully unconstrained ImKAN and MLP extremes. This observation that ImKANs with heavy Hessian regularization match non-regularized MLPs suggests that one should use a combination of the proposed regularization scheme with standard ones, such as 1.2 or dropout (Srivastava et al., 2014), for the best results.
D.1 DETAILED FORMULATION OF THE HESSIAN REGULARIZATION
We use the following fiffteddiferrenes approximation for the second derivative with respect to 21:
f Ou7 (21
f(a + hy. 22) (he + h1) (x1, 12)) + hy s(x1 - he, x2))
(4)
he h, (he h,)
where h4 is the spacing between left and central grid points, while h, is the spacing between central and right grid points. The corresponding expression in terms the coefficients pi.j is given as:
2 hh,ii415 - (h; + h141) Pi.j + Ri41 Pi-
Day
(5)
h, h441 (he hx41).
For the second derivative with respect to 22 we use an analogous expression:
2 2 f(r1, 12 + h2) (hs + hw) f(11, x2) - hu f(x1, x2 hb)
f
(6)
Ou3 (21,22)
hg Nt (ha + lw))
where by and hB are upper and bottom spacings, respectively.
2( Pigti
(hy + h541) + h941 Pis-2)
Drn.zaii
(7)
hy hot1 (h, + hzt1)
For the mixed derivative our finiteddifferenees scheme the following:
ff
2
Da Owr (21,22)
(8)
f(: + by 32 hp) f(r hr, x2 -h)) f(x1 hu, 32 hp) + f(11 hi, x2 hb)
(h, + h1) (hpp + hb)
Denaeis Ps+1, Pit1.111 - 11-1 + Pi-1,5-1
(9)
(h, + hit1) (by hj+1)
The final regularization term is the following:
H D2 D...a + 2D 21 axdss + D2 D
(10)
In order to compute the total regularization term for the whole model, we average Ha, across all the grid points and 2) sum these values across all the 2D functions within all the ImKAN AN layers in the model.
18

---

Published as a conference paper at ICLR 2026
E CUDA KERNELS
The performance of our CUDA kornels with 16x16 tile size in the limit of large dimensions sunmarized in Fig. 12. All the ImKAN curves computed with the largest number of grid intervals G = 20 available for the 16x16 tile size. We compare the inference efficiency of ImKAN and linear layers. On the left panel, normalize time by the shape of the layers. Fig. 12 illustrates clear convergence of these normalized times to the same value for all the dimensions. In the limit of large batch sizes, the forward pass of an ImKAN layer is -8x slower compared to linear layer with the same shape. At the same time, an ImKAN AN layer contains a significantly larger number of parameters than a linear layer of the same shape. Thus, inference time per parameter is significantly better for ImKAN layers, about 27 times, as illustrated on the right panel Fig. 12.
- No Rox 2045
No Nor 2048
Mue
10>41
N. Now 1324 Linear
1024
Uime #perameteed
10-13
30-11
20
25
2
2.2
2.3
2
Bazch
Butch $20
Figure 12: The performance of our CUDA kernels on the H100 SXM GPU in comparison with the lincar layer in the limit of large dimensions. Left panel time normalized by shape. Right panel time normalized by the number of parameters.
Fig. 13 is an analogous illustration but for small dimensions 16 and 32. Our CUDA kernels are better adjusted for such small dimensions, and thus, relative performance compared to linear layers is even higher in this case.
InkaN
No
32
10>2
13-11
10-13
200
214
255
231
Botch
Figure 13: The performance of our CUDA kernels on the H100 SXM GPU in comparison with the linear layer for small dimensions. Left panel - time normalized by shape. Right panel time normalized by the number of parameters.
Finally, Fig. 14 illustrates the inference efficiency depending the number of grid intervals G, which control the mumber of parameters. The time indeed does not depend on G the large batch size limit.
19

---

Published as a conference paper at ICLR 2026
2 10
Uime (batch Site seconds
M
==
10:12
16
4x
3x
24
2
2
2
2
Figure 14: Inference efficiency of an ImKAN layer depending on the number of grid intervals G.
F PRECONDITIONING AND FITTING SCHEME
The first thing attempted upon implementing the CUDA kernels was to fit a model with the highest grid resolution, G 40, supported for the 8 tile on the H100 GPU. this setup, each 2D function had many 41 = 1681 trainable parameters. We found that the training was unstable, so we designed a preconditioning and multi-stage fitting pipeline to stabilize it. We employed this pipeline consistently for all our experiments.
The subsequent evidence revealed that ImKANs (similarly to KANs, as Sec. 4.4 illustrates) progressively harder to fit as grid resolution increases. In other words, our very first experiment was the most challenging one. At more moderate grid resolutions, preconditioning measures can likely be simplified, if not omitted altogether. Specifically, think that fitting scheme omitting additional preconditioning terms, but preserving the Hessian regularization decay phase, which is described in the following, could be effective. With that, below is the description of the current pipeline.
F.1 PRECONDITIONING
We precondition ImKAN layers by adding linear terms into the overall functional form. We use one of the following:
= ImKAN())) + ReLU(Linear(())
= ImKAN(2)) + Linear(ReLU(z))
(1la)
(11b)
where the ImKAN weight, 7, is first set to 0 and then is gradually increased in our multistaged fitting procedure described later. In the case of ReLU-last preconditioning of Eq. 1la, insert ReLU into all the layers except the last one; for ReLU-first preconditioning of Eq. 11b, we insert ReLU into all the layers except the first one. Therefore, at the beginning. when the ImKAN weight is zero, the model is equivalent to pure MI.P-based one for both types of preconditioning.
A merit of ReLU-first preconditioning is that during inference the whole Eq. 11b can be absorbed into a single ImKAN layer whenever the number of grid intervals is even, that is, when the origin is one of the grid points, more details in Appendix G.1. Thus, this type of preconditioning does not increase inference cost in any way. This is an advantage over the original KAN preconditioning scheme (Liu et al., 2024), which requires the additional computation of the computationally expensive transcendental function SiLU for each edge at inference.
20

---

Published as a conference paper at ICLR 2026
Because of the possibility of such an absorption, the total inference FLOPs of ReLU-first preconditioned ImKAN layer are 2x those of a linear layer of the same shape, while for the ReLU-last preconditioning, the slowdown factor is 3x, taking into account the linear branch.
F.2 FITTING PROCEDURE
Our fitting scheme consists of several phases:
Phase I- pure MLP: is set to 0, so the whole architecture is operating in pure MLP mode. This part is typically very short.
ImKAN: is linearly (over time) increased from 0 to 0.3. After that, there is some part with the constant 7 = 0.3. At this phase, we use very strong (= with very high coefficient X) Hessian regularization introduced in Appendix D. For all the subsequent phases, ImKAN weight y is fixed at 0.3. The pipeline is also stable increasing y to 1.0, but in a few, though non-systematic, experiments, we found that using 0.3 value leads to slightly better final accuracy.
Phase II turning
Phase III Hessian regularization decay: At this phase, we gradually decay the strength of the Hessian regularization À from the initial very high value to the target value this regularization intended to be utilized in this fitting procedure and to nearly zero otherwise.
Phase IV Main ImKAN fitting part: In this final phase, we keep Hessian regularization constant at the value reached in the previous phase. The model is fit with the given learning rate schedule. In the experiments this work, we use a constant learning rate for the most part of this phase, and step exponential learning rate decay at the end.
An example of the described fitting procedure for one of the training runs we did for numerical experiments described in Sec. 4.1 of the main text and in the Appendix G.2 is given in Fig. 15.
F.3 COMPARISON OF THE PRECONDITIONING SCHEEMES
We fitted ImKAN models with both types of preconditioning for the CIFAR-10 dataset and for general function approximation. The results are given in Fig. 16 and Fig. 17. Overall, the ReLU-last type of preconditioning appeared to lead to slightly more accurate models, but this small gain accuracy does not justify additional computational cost.
When designing some of our experiments we did not know this yet. Therefore, some of them use the ReLU-last type of preconditioning.
We use the ReLU-last type of preconditioning for Figures 21, 22, 23, and 24. We use the ReLU-first type of preconditioning for Figures 5, 8, 9, 26, and 27.
In other words, the performance of ImKANs on the methane datasets can likely be further improved by switching from the ReL U-lat type of preconditioning to the ReLU-first one. However, since the observed gains in efficiency are already more than an order of magnitude in terms the H100 wall-clock time, we left this for future work.
G EXPERIMENTS
G.1 GENERAL DETAILS ABOUT THE BENCHMARKING PROTOCOLS
Within the scope of this work, we primarily focus on the saturated throughput in the limit of large batch sizes. Thus, we benchmark all the models for progressively large batch sizes until reaching saturation. All the models are benchmarked with 10 warm-up dry runs, and 20 timed runs, Overall, we tried to optimize each model as much as possible while staying within the limits of full precision float32 data type.
MLPs employed in this work consist three types of layers lincar ones, ReLU activations, and batch normalizations. At inference, batch normalizations simply perform clementwise linear transformation and thus can be absorbed into the weights of linear layers. We perform this operation manually and, on top of that, compile the model with the torch,
torch tensorrt
21

---

Published as a conference paper at ICLR 2026
5.8
Pure Loss (MSE, ESAI
1.6
Lose
5x10-1
4.8
0.00
0.00
regulaiiatiinn
10-4
10-4
10
10-1
10-2
Leaming
10-4
Training Steps
The multi-stage fitting procedure we use. Total loss indicates the full loss, including the Hessian regularization term. Pure loss is only the MSE part. For clarity, plot exponential moving averages losses. Note that the horizontal scale is logarithmic. it is linear, the first couple of phases are hard to discern they very short. The fourth phase takes most of the training budget. This training run corresponds to the unregularized ImKAN, where Hessian regularization is turned only at the beginning of the fitting to ensure stability. At the end of phase III, it reaches a nearly zero value, which, in this case, is 10 10-20 20
Figure 15:
backend (with tf32 disabled to ensure full precision float32). We use the same compilation strategy for FastKANs.
For ImKANs, when using ReLU-first preconditioning (see more details in Appendix F), we absorb the entire expression in Fq. 11b into the weights of the ImKAN layer. This modification (x)) +
requires updating the ImKAN 2D functions (xz)) Our construction allows for doing this absorption exactly whenever the origin is of the grid points, which, in turn, is the case when the number of grid intervals G even. On top of that, batch normalizations are absorbed similarly to MLPs. We do not compile ImKAN models.
f(x1,22) + 7f(1),12)
For ImKAN models with the ReLU-last preconditioning we absorb only the ImKAN weight
y-
22

---

Published as a conference paper at ICLR 2026
Accuracy vs Hidden dimensian
Accuracy vs Inference FLOPs
75
Text Accuracy (80
70
60
MLP CNN 1M 1/2
MLP CNN 1/4
55
MLP CNN
MLP CNN
ImKAN CNN ReLU--ast IMKAN RelU-ffrst
InKAN CNN ReLU-Iatt
ImKAN CNN ReLU-frstt
Hidden dimension
Inference FLOPs
Figure 16: Comparison of the preconditioning schemes when fitting ImKANs on the CIFAR- 10 dataset.
MSE vs Hidden dimension
MSE vs Inference FLOPs
MSE vs Inverse H100 SXM throughput
MiP 1/2
VLP 1/2
MLF 3/
MiP
MLP
HLF
IrIANN Redl
IMKN ReLU IKKN RetU-frst
InCAN Rest
IYIAN Rell-Fist
InCAN Rell-ffr
IR 32 R
60-10-1
MSE
50-10-1
8x
IYKAN
ReLu
Rott
IYIGAN ReLL-Fist
Rest
Retlifrst
32 32
1
MSE
60-10-1
128
235
1024
10
10'
10'
10
10
Inference ELOPs
H100 SXM secondsinput (lagg batch si2))
Hidden dimension
Figure 17: Comparison of the preconditioning schemes when fitting ImKANs within general function approximation setup.
G.2 GENERAL FUNCTION APPROXIMATION
Both MLPs and ImKANs use batch normalizations. We set affine=True for MLPs as it is the standard choice, and af iff for ImKANs in accordance with sigma grids introduced in Sec. 3.1. MLPs use ReLU activations, while ImKANs do not require any additional activation functions. We use = 12 for all the ImKAN models, as this was the optimal value found in the ablation study described below. Pseudocode for both models is available in Fig. 18. Both students have two hidden layers, which is one more than both Cybenko (Cybenko, 1989) (the one for MLPs) and Kolmogorov-Arus Arnolld universal approximation theorems require. This setup, however, is more realistic, as MLPs with exactly one hidden layer are rarely used in practice.
23

---

Published as a conference paper at ICLR 2026
MLP student
ImKAN student
Linear (input_ dim -> hidden, dim) Batchhorrmdd (wwdoeendddddddddddddaaaeeaaaaaa affine=True) ReLU()
InKANLayer (input_dd dim + hidden dim) BatchNormid (hidden_dddddddddddd affine=False)
Linear (hidden, dim hidden_dimg dim)
1mKANLayer (hidden_dinp dim + hidden dim) aarchlorrmddhiiddee.c dim, af affime=alsse)
BatchNormid (hidden dim, affine=True) ReLU() Linear (hidden dim output_das (imm
ImKANLayer (hidden_dis + output autput_dim)
Figure 18: Pseudo code for MLP and ImKAN students.
The performance of ImKANs when approximating RR3 -> R function generated similarly as the R232 IR_ described the main text is given in Fig. 19.
MSE Hidden dimension
MSE Inference FLOPs
MSE Inverse H100 SXM throughout
1z MKAN
MLP 1/2
Squared Erer
VLP
1x
6x:10
4.2%
128
250
517
1021
10
10'
10*
10
10
5XM seconds/input (large
Hidden dimension
Figure 19: ImKAN MLP for general, R + 32 function approximation. The 'MLP 1/2' line corresponds to the outcome of the fitting procedure with only half of the training steps compared to the 'MLP" one.
There is a trend that the relative performance of ImKANS improves with the scale. It is clearly seen on the MSE vs FLOPs panel. On the MSE vs H100 wall-clock time panel, it is first. masked by the non-homogeneous efficiency of the code, but next still reveals itself for the largest hidden dimensions.
We investigate the effect of the chosen number of grid intervals G on the resulting ImKAN accuracy when approximating the R + R!! function with hidden dim = 256. The result is given in Fig. 21.
Mean Squared Error
0.4700
0.4675
0.4650
0.4625
0.4600
34
8
12
16
20
26
32
40
Figure 20: Final MSE vs G for the hidden dim-256 ImKAN model.
24

---

Published as a conference paper at ICLR 2026
G.3 PLOT WITH SHARED AXIS FOR THE NUMBER OF PARAMETERS PER EACH INNER FUNCTION
0.50
Intolt
InfiN 30
MLP D1P33350 k1
0.52
0.50
0.48
0.48
In6AN
Intrenec 1.00
1032 209
10"
Mamber of amazara inner Nanczion
16 12 Inkck 23
Figure 21: The same plot as Fig. 4 the main text, but in vertical format and with the shared axis.
G.4 МETHANE
Tabular datasets are the natural choice for feedforward fully connected neural networks. Popular tabular datasets, such the Titanic (Kaggle) housing prices (Harrison Jr & Rubinfeld, 1978), however, are not particularly convenient for this purpose. First, they are typically stochastic in nature for instance, while it is possible to improve guess on the survival based on the data available for the Titanic dataset, it is impossible to say for sure. Thus, even an arbitrarily large model fitted on arbitrarily many data points would have non-zero limitation on the accuracy. In other words, the performance of model translates into an error metric not so directly, making the comparisons between different models less illustrative. Second, these dalasets are typically relatively small, making it challenging to sweep across wide range of model scales obtain a comprehensive picture of performance.
Machine learning models fit on such datasets belong to the class of so-called machine learning interatomic potentials (Behler & Parrinello, 2007; Bartók et al., 2010). This dataset is sufficiently large for the comparisons, containing more than seven million configurations. Additionally, this dataset deterministic the geometry of the corresponding methane configuration completely determines the target (formally, there can be stochastic term due to the lack of complete convergence of ab initio computations for the quantum- quantum-mechanich. mechanical energy, but it negligible in practice).
25

---

Published as a conference paper at ICLR 2026
Table 1: Summary of methane representations
Rotational symmetry
Permutational symmetry
Label
Features
Cartesian Components
Augmentations
Augmentations
12
Distances
Features
Augmentations
Cartesian Components Polynomials
Augmentations
Features
34
Distances Polynomials
Features
Features
31
The target, the potential energy of the system, is invariant with respect to rotations and permutations identical atoms Therefore, there are several viable representations the methane molecules depending on how these symmetries are addressed:
Cartesian Components: The simplest representation just a collection all the Cartesian components of all displacement vectors from the carbon atom to all the hydrogen atoms. Since each methane molecule contains 4 hydrogen atoms, the total number of displacement vectors is and the total number of components is 12. When using this representation, simply concatenate all these components together and feed them to a fully connected MLP or ImKAN whose input dimension 12. This representation is not invariant with respect to both rotations and permutations; thus, we use the corresponding augmentations during training. We randomly permute hydrogen atoms and rotate each molecule whenever we sample a minibatch from the training subset for each step of stochastic gradient descent.
Distances: Another possible representation a collection of all the interatomic distances between all the atoms. Since the total number of atoms is 5, the number of all the interatomic distances is 5 4/2 = 10. Therefore, the input dimension of fully connected networks applied to this representation is 10. This representation is invariant with respect to rotations but not with respect to permutations. During training, we use only permutational augmentations.
FFormally there is an additional symmetry, inversion, but the corresponding group contains only two elements, thus it does not make much sense to treat it separately. We unite it with the group of rotations, and in the following, by rotation we mean proper improper rotation.
26

---

Published as a conference paper at ICLR 2026
127
MLP IMMA
MLP
MLP
KMN
RAN
10-1
1.7
124
10.
10a
13"
MUP
NLP1/2
NLP
MLP
MIP
MIP
mwW
1021.5K
71.0x
100
100
10-
10:4
10"7
MLP 2/4
MLP
MLP 1/2
MEP
MLP 1/2
MLP
HLP
MLP
Companent
4820
2x:0
19:1
4.14
1.7
200
1.a
12a
307
15-
20-4
1017
MLP 1-4
MIP
MLP 1/4
10"1
S0LP 3/2
YLP 3/2
SULP 3/2
S0L
VLP
SIL
IMKAN
reymomia
Rop.
Rop.
4x20-
3x20
2x10
100
1.24
Inforrnce
PO00 sucoTeda melcccul butch slas!
Friddon dimetaioo
Figure 22: ImKAN vs MLP on the dataset of randomly displaced methane configurations. "InKAN Reg." curve corresponds to ImKAN fitted with Hessian regularization introduced in Appendix D. On the vertical axis, we plot the relative Root Mean Squared Error, which is given test RMSE normalized by the standard deviation of the target in the dataset. The 'MLP 1/2" and "MLP 1/4" curves correspond to outcomes of fitting procedures with half and quarter of the training budget, respectively.
Cartesian Components Polynomials: We compute power sum symmetric polynomials on top of the Cartesian components of the displacement vectors: Pa..a,,a, 1i-1- x, y, 2 2 +a, +a, < 4. The total number of such synmetric polynomials 34 (excluding trivial Po.o.o). This representation is invariant with respect to permutations but not with respect to rotations. Thus, during training we use only rotational augmentations.
for nen-tegative integer ag + ak,
Distances Polynomials: The final representation is a collection of non-trivial symmetric polynomials on top the interatomic distances, constructed similarly to Allen et al. (2021). The total number of such polynomials is 31, and their exact formulas are given in Ap- pendix G.4. This representation invariant with respect to both rotations and permutations. Thus, we do not any augmentations during training for this representation.
The described representations are summarized Table 1. We systematically evaluate all four possible combinations of how the rotational and permutational symmetries can be
27

---

Published as a conference paper at ICLR 2026
incorporated into the fitting pipeline. Within the Distances Polynomials representation, the methane dataset tabular in the classical sense - it is a table with about 7.7 million rows and 31 columns. For other representations, the dataset is tabular-like given the available augmentation strategies. We randomly split the data into 7000000, 300000, and 432488 train, validation, and test molecules, respectively.
For cach ropresentation, we fit the same families of MLP and ImKAN models as in the previous section. The result is given in Fig. 22. For this dataset, we G = 28, the optimal value we found in ablation studies. Similarly to the previous experiment, we demonst rate tight convergence of the baseline MLP models by providing three lines corresponding to full, half, and quarter of the training budget, respectively. Overall, when compared to domain-specific architectures, typically given by GNNs (Zhang al., 2021) and/or transformers (Pozdnyakov & Ceriotti, 2023), the introduced feedforward fully connected models occupy a non-overlapping part of the Pareto frontier they are less accurate, but also orders of magnitude faster.
The figure illustrates that ImKANs consistently outperform MLPs across all modalities. Furthermore, the performance improvement is much larger compared to our previous ex- periment. At the same accuracy level, ImKANs require up to many dozen times (or even more for the Distances modality) less inference FLOPs, which results in more than 10x improvement of the inference H100 wall-clock time.
Furthermore, Fig. 22 provides carly indications that ImKANs sometimes can be more accurate in the limit of large scale, that is, to have better generalizabillty. The second row of the figure, corresponding to the Distances modality, rllustrates that the rate of improvement of MLP models becomes very slow, and it is questionable if this family of models would ever surpass the accuracy achieved by ImKAN models at any scale.
On the other hand, depending on the nature of the data, raw ImKANs, without the Hessian regularization we proposed in Appendix D, can be more prone to overfitting. This happens for the Distances Polynomials modality as the last row of Fig. 22 illustrates. This modality incorporates all the symmetrics into the representation and does not involve any sort of augmentations. Therefore, it is likely that the generalization problem we outlined in Appendix D takes place for this fitting setup. As the green line the fourth row Fig. 22 illustrates, the Hessian regularization is sulficient to overcome the overfitting. Properly regularized ImKANs were found to outperform the MLPs and be Pareto-optimal from the point of view of both inference FLOPs and inference H100 wall-clock time.
G.4.1 ABLATIONS
The ablation study on the effect of Hessian regularization given in Fig. 23, and the effect of the number of grid intervals G in Fig. 24.
Figure 23: Effecc of the strength of Hessian regularization on the validation error when fitting lmKAN with hidden dim 256 on the methane dataset using the Distances Polynomials representation.
28

---

Published as a conference paper at ICLR 2026
0.253
9922 510
0.30
0.355
Figure 24: Effect of the number of grid intervals G the validation error when fitting ImKAN with hidden dim 128 on the methane dataset using the Cartesian Components representation.
29

---

Published as a conference paper at ICLR 2026
G.4.2 THE Distances Polynomials REPRESENTATION
was mentioned that the Distances Polynomials representation is given by non-trivial invariant polynomials computed on top of interatomic distances. These polynomials are constant with respect to changing the order of identical hydrogen atorss.
P1 = 26 +
+
27
48 + 2g
210
P2 Iι
P3 of + 20 +
27
23
xig 110
P4 = T5 + #537 2647 + 55 + Tg + 2520 + 172gg + #g
Tg + + 10 Ø9
Ps = 01 + + + 23 + 21 + F4 1228 + F3
F22g +33000 14
Pg 12
P7 + ,
Ps = 0 +
262 + + zZøE
+ + ry + rzaI
+ +o assi 19
Pg وι1د +u + + 2z2
22 2s0
P10 = 25 + +6070
P1 2 + +2 IzJeTss I2 + 24 + +
1 +20
P12 = + +2zs +zass
12.0 +2
P13 22 +1 12431s -24409 1334310
P44 = z + +
P1s = o + 2 + 2
P1 03
N62 + agas + agast t +z xz + rsunt
agan + 2 + Ig + zy
3
P17 mas + +2 -x x mzng + xa
22ag + + 13
Pis 2 aage
+ + x 2i 22 + 2zg
22252 + 2z + +2 22 + z2 + و + 12
+ t + ca + n + 13 + 1
P19 + 12 m + uaal + 11a
+ + + + 2, 24 13
+ +
+
P20 2/ -d
2 23 2
12 2500
(12)
The exact form of these polynomials given in Eq. 12 and Eq. 13, where Tt, D2. Z correspond to interatomic distances CH1, CH2, CH3, CH4, H H2, Hi Ha, H H,H,, H₂ HsH
30

---

Published as a conference paper at ICLR 2026
H2H4, and HyH, respectively. The presented polynomials form invariant generators (Derksen & Kemper, 2015) of the group corresponding to arbitrary permutations of the hydrogen atoms. Therefore, the Distances Polynomials representation preserves all the information about the initial CH4 molecule.
P21 310223 2y zyazat + + 1g3
P22 ota500 ufu507 + +āosas 22 + xaayag + x2
z6 + 1 + +
P23 uios +2 + + 2/07 ojes + rzgss
22209 + +
P24 x -
P25 23 + +
P26 mn tr + m2ng + maggh+
122g +30
P2t z 1 + +
+ + 22 + rgrjast
F2 +2 123 + + 22rro + x20s2ÿ+
x2nsa
Is + 14002/0
P2s = +2 +2 +2
+2
P29 01222 = +a rgr
P30 u/ +
2. 2nd + rggu
P31 x/ +
+ 2irre7 aseaus + x2032s+
13244 + 222 1 22
(13)
G.5 LMKAN-BASED CONVOLUTIONAL NEURAL NETWORKS
For CIFAR-10, our backbone architecture consists of five x 2 convolutions, each with stride and two fully connected layers at the end. Since the resolution of CIFAR-10 images is 32 x 32, where 32 = 25, five 2 2 convolutions with stride 2 transform the spatial dimensions of an image exactly to
All the layers use the same width (= number of filters in case of convolutions, and hidden dimension in case of fully connected layers), which vary for both families of the models. In other aspects, the models are similar those we employed in previous sections and without for ImKAN-CNNs; MLP-CNNs use ReLU activations, while ImKAN-CNNs do not require additional activation layers.
use batch normalizations with affine transforms for
The dataset comes with pre-defined full training and test subsots. We split the full training subset into training and validation parts in a 90% 110% ratio. Our augmentation pipeline consists established techniques, such as RandAngment Cubuk et al. (2020), MixUp Zhang et al. (2017), CutMix Yun et al. (2019), and few others.
31

---

Published as a conference paper at ICLR 2026
CIFAR-10 CNN backbone
ImageNet CNN backbone
# Only convolutional and fully connected layers shown
# Only convolutional and fully connected layers shown
# [32, 32, 3] + [16, 16, width] Conv2D (3 + width, kernel size 2, stride 2)
[81, 81, 3] - [27, 27, base base_vidta width] Conv2D (3 + base width, kernel_stie size 3, stride = 3)
[16, 16, width] + [8, 8, width] (width + width, kernel_size 2, stride = 2)
[27, 27, base_width, width]
+ [9, 9, 3*
Conv2D (base_widtde width - 3*base width, kernel size 3, stride 3)
[8, 8, width] + [4, 4, width] (withD + width, kernel_size 2, stride - 2)
[9, 9, 3*base_widee + [3, 3, 9* base base_widee width]
Conv2D(3*base (3+base width + 9*base width, kernel size 3, stride 3)
[2, 2, width] CondhD(width + width, kernel_sit size 2, stride = 2)
[4, 4, width]
# [3, 3, g*base_widthl + [1, 1, 27* base Conv2D width 27+base __idth, bernel_size 3, stride 3)
base_width] width]
[2, 2, width]
[1, 1, width]
Conv2D(width >> width, @ernel_size
2, stride 2)
(27 width + 27+ base width) FullyConnectedd274aaae width 1000)
FullyConnected (widhh width) (ully (width + 10)
Figure 25: CIFAR-10 and ImageNet CNN backbones. MLP-based CNNs additionally have ReLU activations and batch normalizations with enabled affine transforms. ImK AA-baed CNNs do not require additional activations and use batch normalizations without affine transforms as suggested by our sigma grids described in Sec. 3.1.
We use the following pool of augmentations for the CIFAR-10 dataset:
CIFAR-10 augmentation pipeline
MEAN (0.4914, 0.4822, 0.4465) STD (0.2470, 0424.5, 0.2616)
nn. Sequential
T. (3ndonCrrop322,, padding=4), TandoonHorrzzontalllio
T. ColorJitter 0.3, 0.3, 0.05),
T. RandAugment (2, 7),
T. RandomErasing
oc.le=(0.055 0.2), ratio= 3.3)), Normalize (MEAN, STD),
)
On top of these, we use MixUp (Zhang al., 2017) (a= 0.2) and CutMix (Yun al., 2019) (B 1.0) augmentations, both with 50% probability.
When fitting the families of convolutional neural networks described in the main text, we use the above pool of augmentations consistently for MLP-based and ImKAN-based CNNs.
For ImageNet, the standard data preparation pipeline introduced by (lexNet(Krizhevsik et al., 2012) involves first resizing image to 256 pixcls along the smallest dimension, then performing a random crop of 224 x 224 pixels during training, and a center crop of 224 x 224 pixels during validation.
We mimic this procedure by first resizing the image to 81 256/224 93 pixels across the smallest dimension, and then performing random center crops of 81 81 pixels.
32

---

Published as a conference paper at ICLR 2026
86
84
Validation
2.5
5.0
7.5
10.0
12.5
15.0
17.5
20.0
Figure 26: Accuracy of the ImKAN-basel CNNs on the CIFAR-10 datasct depending on the grid resolution G.
88
Train Validation
84
Accuracy
80
10-1
10-4
10-²
10a
regularization strength,
Figure 27: Accuracy of the ImKAN-based CNNs on the CIFAR-10 dataset depending on the strength of Hessian regularization.
Next, we the following augmentation pipeline:
ImageNet augmentation pipeline
(
RandonfforizznntaFF1i..))
T. (obl:riitterrbrightteeeet 4, contrast=0.44 4, saturation=0.4. 44, 1))=0.1),
T. Randdugment(()..
RandonErasing((per.... (p=0.25 25, scale=(0.00 02, 0. 0.33), ratio=(0.3, 3.3), value=0), Normoorllll ize (mean= [0. 485, 456, 0.406], std= 229, 0.224, 0.225]),
)
On top of these, we use MixUp and CutMix.
G.6 COMPARISON WITH FASTKAN
As was already mentioned in the main text, we use the training script for the CIFAR- 10 datasct available in the FastKAN GitHub repository (Li, 2024a) the basis for the comparison of ImKAN and FastKAN. However, we provide several modifications to the pipeline.
The original script implements the fitting procedure fully connected FastKAN model the CIFAR-10 dataset. without augmentations. The model has only one hidden layer
33

---

Published as a conference paper at ICLR 2026
with 256 neurons. Without augmentations, it overfits the data quickly. Thus, the very short fitting procedure in the original script is sufficient.
We extend the script by the CIFAR-10 dataset. We observed that because of augmentations, one has to the model for longer time, so the training budget substantially increased. The data split properly into train, validation, and test subsets, while the original script employed only a train-validation split. We use the cosine (without restarts) learning rate scheduler (Loshchilov & Hutter, 2016) instead of the exponential decay one in the original script. Finally, the normalization was performed with real values of the mean and standard deviation for the CIFAR-10 dataset, instead of the dummy 0.5 values of the original script.
augmentation pipeline
used in Sec. 4.3 for the
These modifications significantly improve the performance of FastKAN models (x= 54 55% validation accuracy in the original script), see Fig. 9. Furthermore, even the MLP bascline of the same shape yiclds better accuracy compared the performance of the FastKAN model in the original script.
34


Reviews:
[{"review_id": "ITy4WnwMR3", "reviewer": "Reviewer_UnT9", "summary": "This paper introduces lmKANs, a lookup multivariate KAN that replaces the high dimensionality of the univariate functions employed in KANs with multivariate low-dimensional functions, where the inner functions are spline lookup tables. The paper delivers custom CUDA kernels for accelerated inference on lmKANs, conferring 88x faster inference than standard KANs. The paper then compares this implementation with standard MLPs across a range of tasks and demonstrates that lmKANs outperform in terms of accuracy and performance.", "questions": "Do you have statistics around training time when training lmKANs to convergence versus an equivalent MLP? Figure 18 starts to get at this result with the MLP 1/2 trace, but I would love full training stats.\n\nHow sensitive/brittle is training as a function of HPO over G? Is it possible to provide some sensitivity analysis of lmKAN performance to G? I'm trying to understand how much increasing the parameter capacity of KAN's improves expressivity.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "The paper present an important contribution to KAN development, making KANs easier to train and more scalable. The method presented is novel and is a significant improvement to address hurdles when training KANs.\n\nThe paper does a good job benchmarking lmKANs against MLPs. I particularly appreciate the MLP 1/2 traces, giving a sense of the compute benefits afforded by lmKANs.\n\nOverall, the approach is interesting and novel, and provides a way to increase the expressivity of KAN training.", "weaknesses": "In the main paper, the authors focus far more on the computational benefit afforded by the accelerations than the theoretical foundation of how this model. However, the supplementary information dives into the theory well and motivates lmKANs, all of which is missing from the main manuscript. This unfortunately detracts from the impact of the main paper. I also struggled a bit to connect the dots between the motivation of the paper and the results presented. Adding dozens or hundreds of more parameters should significantly increase expressivity, but for some plots (e.g. Figure 4) the difference in MSE is on the order of 1e-2, which is surprising given the motivation.\n\nThe paper also claims that e.g. Moradzadeh et al do not compare against MLPs, which is untrue, taking a look at the manuscript. This oversight detracts from the novelty of the proposed method, given that other similar approaches to also result in KANs that match performance or outperform MLPs.\n\nWhile the idea of using a lookup table to convert internal functions to O(1) complexity is interesting, the authors don't explore the added burden of HPO to find the optimal grid size to match or exceed MLP performance.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "R9tSGm39RT", "reviewer": "Reviewer_HStq", "summary": "This paper proposes lookup multivariate Kolmogorov-Arnold Networks (lmKANs), a novel architecture that replaces high-dimensional linear mappings with trainable low-dimensional multivariate functions (primarily d=2) implemented via spline lookup tables.\n\nlmKANs split high-dimensional inputs into low-dimensional subgroups (e.g., 32D input → 16×2D subgroups), each processed by a multivariate spline function.\nSpline lookup tables leverage the local support of B-splines to achieve O(1) inference per subgroup. This design reduces inference FLOPs by up to 6× for general function approximation.\n\nCustom CUDA kernels further accelerate computation, achieving up to 88× faster inference on H100 GPUs. Across tasks, lmKANs consistently outperform MLPs and FastKANs in both accuracy and efficiency.", "questions": "Appendix D mentions that second-order B-splines match ReLU’s smoothness. Did the authors try higher-order splines?\n\nAlso, Hessian regularization is only mentioned in the appendix. Was it applied to all experiments in the main paper? If so, how does it affect the results?\n\nDid the author tried other options for mapping unbounded activation to bounded regions? For example, Erf-based function, such that we may assume the activation follows a normal distribution.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The idea is original and overall design logics is very well organized. By grouping high-dimensional inputs into low-dimensional subgroups, lmKANs avoid the exponential complexity of high-dimensional splines. The look-up table take advantage of this The Sigmoid-like σ-grid very elegant way to adapts to unbound activations into bounded region.\n\nExperiments are solid: across four task types, the paper tests different model scales, training budgets, and data formats to show lmKANs’ robustness.\nEngineering details, including CUDA implementation, are clear too.\n\nAlso, the paper does a good job of explaining a relatively complex idea in a clear and accessible way.", "weaknesses": "Limited high-dimensional multivariate validation\nThe paper focuses exclusively on 2D lmKANs. It would provide values to validate for other d=4/8. Even on very small experiments this would help.\n\nCNN architecture and  implementation\nThe experiments on CIFAR uses “For CIFAR-10, our backbone architecture consists of five 2 × 2 convolutions, each with stride 2” and for ImageNet “For ImageNet, we downsample the images to a resolution of 81 × 81 pixels. Our backbone consists of four convolutional layers with the 3 × 3 kernel size and stride 3 and two fully connected layers” These are not typical architectures for these datasets. Does this imply lmKANs cannot integrate with standard models like ResNet or AlexNet?\n\nAdditionally, Section 4.3 mentions that lmKAN-CNNs “cast convolutions to fully connected layers via memory manipulations.” This seems like a temporary workaround rather than a native multivariate spline convolution.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "JIlqMYSHL3", "reviewer": "Reviewer_63pm", "summary": "The paper introduces lookup multivariate Kolmogorov-Arnold Networks (lmKANs). The central idea is to replace the 1D univariate activation functions used in standard KANs with 2D multivariate functions. The key component of this architecture is the implementation of these 2D functions as $\\mathcal{O}(1)$ spline lookup tables. The authors posit that this approach \"decouples\" the model's parameter capacity which could be increased by grid point from its inference cost. This, in theory, allows lmKANs to possess a significantly larger parameter count than an MLP with only a $\\sim 2\\times$ increase in theoretical FLOPs. The paper presents experiments on general function approximation, the methane molecule dataset, and image classification to demonstrate that lmKANs achieve a superior trade-off between accuracy and inference cost compared to MLPs.", "questions": "- Can the authors please provide the critical experiment described in Weakness #2? Specifically, for the \"General Function Approximation\" task (from Fig 4), please provide a plot with a fixed hidden dimension (e.g., $N=256$) while sweeping $G$ (e.g., $G \\in \\{4, 8, 16, 32\\}$). This plot must show (a) Test MSE vs. $G$ and (b) Inference FLOPs vs. $G$. This is essential to validate the paper's core premise.\n- Why did the authors choose to compare against the RBF-based FastKAN to justify the 1D vs. 2D claim, rather than the much more appropriate baseline of a 1D B-spline KAN using the same lookup table method? Without this direct comparison, the claim that 2D functions are inherently more stable is unsupported.\n- The appendix figures (e.g., Fig 19) show that performance degrades as $G$ becomes very large. This seems to contradict the motivation that multivariate functions can \"accommodate\" a large number of parameters without issue. Does this imply that the \"free parameter\" benefit is limited in practice? How do the authors reconcile this with their primary motivation?", "limitations": "", "rating": 2, "confidence": 5, "soundness": 1, "presentation": 2, "contribution": 2, "strengths": "- The idea of using low-dimensional multivariate functions (specifically 2D) implemented as $\\mathcal{O}(1)$ spline lookup tables is novel.\n- The authors have clearly put significant effort into the implementation, including custom CUDA kernels and a sophisticated multi-stage training pipeline with Hessian regularization.", "weaknesses": "- The paper's primary thesis is that lmKANs are superior because parameters can be added \"for free\" in terms of FLOPs by increasing the grid resolution $G$. This is the claimed advantage over MLPs, where adding parameters requires increasing width $N$ at an $\\mathcal{O}(N^2)$ FLOPs cost. However, the main experimental comparisons (Figs 4, 6, 7) do not test this thesis. These experiments fix $G$ and vary $N$, just like an MLP. These results show, at best, that lmKAN's $\\mathcal{O}(N^2)$ scaling factor is better than an MLP's, but they provide zero evidence for the central claim about the benefit of scaling $G$.\n- The most important experiment is missing from the paper. To validate the core claim, the authors must provide an analysis where network width $N$ is fixed, while grid resolution $G$ is varied (e.g., $G \\in \\{4, 8, 16, 32\\}$). This experiment must plot two curves on the same y-axis: (1) model performance (e.g., MSE) and (2) inference cost (FLOPs and/or wall-clock time). The paper's hypothesis predicts that performance will improve while inference cost remains flat. Without this plot, the main premise is entirely unsubstantiated.\n- The only analysis of $G$ is relegated to the appendix (e.g., Fig 19, 22, 24). These plots are incomplete, showing only \"Performance vs. $G$\". They critically omit the corresponding \"Inference Cost vs. $G$\" data, which is essential to validating the trade-off. Furthermore, these plots (e.g., Fig 19) show a U-shaped curve where performance worsens after $G=12$, which contradicts the paper's narrative that multivariate functions can \"accommodate\" large parameters without issue. This is not adequately discussed.\n- The paper motivates the move from 1D to 2D functions by claiming 1D functions with dense grids are \"unstable\". To support this, it compares its 2D B-spline lmKAN to FastKAN (Fig 8), which uses 1D Gaussian RBFs. This is an invalid. The degradation in FastKAN could be a property of RBFs, not a fundamental flaw of 1D spline functions. The correct baseline would be a 1D B-spline KAN using the same $\\mathcal{O}(1)$ lookup technology. This baseline is missing.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "eIIxAuQzb9", "reviewer": "Reviewer_jcL7", "summary": "This paper proposes to increase KANs' inference efficiency via lookup tables. Their construction expresses a general high-dimensional mapping through trainable low-dimensional multivariate functions. Although these functions may carry a large number of trainable parameters, they sparsely activate during inference time; as a result, it can be implemented as a sparse lookup table. They show the effectiveness of the proposed method on a wide range of datasets, outperforming MLPs and FastKANs.", "questions": "* The default setup is $d=2$, how about $d=1$ or $d=3$ or larger? What's the philosophy for choosing $d=2$?\n* How about training efficiency? When the batch is larger than 1, different samples can activate different basis functions.\n* Section 4.3, it should be \"IMKAN\", not \"LMKAN\"", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The paper is well written and the presentation is good\n* The proposed sparse lookup table method is simple yet effective\n* Experiments are extensive and convincing", "weaknesses": "* The application range of the method should be discussed, i.e., when do we expect this method to be less of an advantage?\n* Although the proposed method can increase inference efficiency, it does not necessarily help training when a larger than 1 batch is being used.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)