Paper deep dive

Scaling of Gaussian Kolmogorov--Arnold Networks

Amir Noorizadegan, Sifan Wang

Year: 2026Venue: arXiv preprintArea: cs.CEType: PreprintEmbeddings: 86

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 96%

Last extracted: 4/26/2026, 7:41:57 PM

Summary

The paper investigates the role of the Gaussian scale parameter (epsilon) in Gaussian Kolmogorov-Arnold Networks (KANs). It identifies the first layer as the critical bottleneck for scale selection, as the first-layer feature geometry determines the distinguishability of inputs. The authors propose a practical operating interval for epsilon, [1/(G-1), 2/(G-1)], where G is the number of Gaussian centers. This interval is validated across various function-approximation problems (smooth, periodic, sharp, and discontinuous) and a physics-informed Helmholtz problem, demonstrating that proper scaling improves stability, accuracy, and training efficiency.

Entities (7)

B-spline · basis-function · 100%Chebyshev · basis-function · 100%Gaussian Kolmogorov-Arnold Networks · neural-network-architecture · 100%Gaussian scale parameter · hyperparameter · 100%Helmholtz problem · scientific-problem · 100%Kolmogorov-Arnold Networks · neural-network-architecture · 100%First-layer feature matrix · mathematical-object · 90%

Relation Signals (4)

Gaussian scale parameter → hasoptimalinterval → [1/(G-1), 2/(G-1)]

confidence 100% · identify a practical operating interval, ϵ∈[1/(G-1), 2/(G-1)]

Gaussian Kolmogorov-Arnold Networks → isinspiredby → Kolmogorov-Arnold Networks

confidence 100% · The original KAN formulation is built from B-spline bases... Gaussian KANs are part of the development of RBF methods.

Gaussian Kolmogorov-Arnold Networks → uses → Gaussian scale parameter

confidence 100% · The Gaussian scale parameter ϵ is central to the behavior of Gaussian Kolmogorov--Arnold Networks (KANs)

First-layer feature matrix → governs → Gaussian Kolmogorov-Arnold Networks

confidence 90% · scale selection is governed primarily by the first layer... we analyze the first-layer feature matrix

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:The Gaussian scale parameter \(\epsilon\) is central to the behavior of Gaussian Kolmogorov--Arnold Networks (KANs), yet its role in deep edge-based architectures has not been studied systematically. In this paper, we investigate how \(\epsilon\) affects Gaussian KANs through first-layer feature geometry, conditioning, and approximation behavior. Our central observation is that scale selection is governed primarily by the first layer, since it is the only layer constructed directly on the input domain and any loss of distinguishability introduced there cannot be recovered by later layers. From this viewpoint, we analyze the first-layer feature matrix and identify a practical operating interval, \[ \epsilon \in \left[\frac{1}{G-1},\frac{2}{G-1}\right], \] where \(G\) denotes the number of Gaussian centers. For the standard shared-center Gaussian KAN used in current practice, we interpret this interval not as a universal optimality result, but as a stable and effective design rule, and validate it through brute-force sweeps over \(\epsilon\) across function-approximation problems with different collocation densities, grid resolutions, network architectures, and input dimensions, as well as a physics-informed Helmholtz problem. We further show that this range is useful for fixed-scale selection, variable-scale constructions, constrained training of \(\epsilon\), and efficient scale search using early training MSE. Finally, using a matched Chebyshev reference, we show that a properly scaled Gaussian KAN can already be competitive in accuracy relative to another standard KAN basis. In this way, the paper positions scale selection as a practical design principle for Gaussian KANs rather than as an ad hoc hyperparameter choice.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

85,809 characters extracted from source content.

Expand or collapse full text

Scaling of Gaussian Kolmogorov–Arnold Networks Amir Noorizadegan1,∗, Sifan Wang2 1Department of Mathematics, Hong Kong Baptist University, Hong Kong SAR, China 2Institution for Foundations of Data Science, Yale University, New Haven, CT 06520, USA *Corresponding author: amir_noori@hkbu.edu.hk Abstract The Gaussian scale parameter ϵε is central to the behavior of Gaussian Kolmogorov–Arnold Networks (KANs), yet its role in deep edge-based architectures has not been studied systematically. In this paper, we investigate how ϵε affects Gaussian KANs through first-layer feature geometry, conditioning, and approximation behavior. Our central observation is that scale selection is governed primarily by the first layer, since it is the only layer constructed directly on the input domain and any loss of distinguishability introduced there cannot be recovered by later layers. From this viewpoint, we analyze the first-layer feature matrix and identify a practical operating interval, ϵ∈[1G−1,2G−1],ε∈ [ 1G-1, 2G-1 ], where G denotes the number of Gaussian centers. For the standard shared-center Gaussian KAN used in current practice, we interpret this interval not as a universal optimality result, but as a stable and effective design rule, and validate it through brute-force sweeps over ϵε across function-approximation problems with different collocation densities, grid resolutions, network architectures, and input dimensions, as well as a physics-informed Helmholtz problem. We further show that this range is useful for fixed-scale selection, variable-scale constructions, constrained training of ϵε, and efficient scale search using early training MSE. Finally, using a matched Chebyshev reference, we show that a properly scaled Gaussian KAN can already be competitive in accuracy relative to another standard KAN basis. In this way, the paper positions scale selection as a practical design principle for Gaussian KANs rather than as an ad hoc hyperparameter choice. The implementation is available at https://github.com/AmirNoori68/Gaussian-KAN. Keywords: Gaussian KAN; Kolmogorov–Arnold networks; Gaussian shape parameter; conditioning; kernel methods; radial basis functions; physics-informed neural networks. 1 Introduction Kolmogorov–Arnold Networks have recently emerged as a structured alternative to classical multilayer perceptrons (MLPs), replacing fixed node activations with learnable univariate functions placed on edges [1, 2]. This edge-based construction provides a different inductive bias for multivariate approximation and has attracted growing interest in scientific machine learning, regression, and physics-informed modeling [3, 4]. The original KAN formulation is built from B-spline bases [1], and a rapidly expanding literature has explored alternative basis families to improve smoothness, computational efficiency, and approximation quality. Recent KAN variants have employed many different basis functions, including Chebyshev-type constructions [5, 6, 7], Jacobi and rational Jacobi bases [8, 9], Fourier-based activations [10, 11], ReLU-based and adaptive piecewise bases [13, 12, 14], as well as wavelet, finite-basis, and polynomial variants [15, 16, 17]. In parallel, recent work has also focused on improving adaptability and efficiency through geometric or grid refinement strategies [18], grid-adaptive physics-informed KANs [14], sparse-identification frameworks [20], JAX/GPU implementations [6, 14], and improved initialization strategies [19]. Among these basis families, Gaussian radial basis functions (RBFs) are particularly attractive because of their smoothness, locality, and simple analytical form. More broadly, Gaussian KANs are part of the long development of RBF methods. Starting from Hardy’s multiquadric formulation [21] and Kansa’s extension to PDEs [22], RBF methods have developed into a large literature on approximation, interpolation, and scientific computing. In the KAN setting, replacing B-spline activations with Gaussian bases was first proposed in FastKAN [23], where Gaussian functions were introduced as computationally efficient surrogates for spline-based activations. In that work, the scale parameter was fixed, mainly to show that Gaussian bases can approximate spline constructions while reducing implementation complexity and computational cost. Subsequent works extended Gaussian-based KANs in several directions, including residual formulations [24], operator-learning architectures with learnable Gaussian bases [25], improved training efficiency and accuracy [26, 27], and complex-valued Gaussian KANs [35, 36]. Beyond regression and PDE-related applications, Gaussian-based KANs have also shown promising performance in classification, with reduced parameter counts, faster training, and improved stability in large-scale or high-dimensional settings [37]. Taken together, these works show that Gaussian bases are not merely substitutes for splines, but rather constitute a flexible and expressive class of KAN representations. What remains largely unresolved, however, is how to choose the Gaussian scale parameter ϵε in a neural-network setting, or more generally the scale in RBFs that admit such a parameter, since not all RBFs do. In most existing Gaussian-KAN works, ϵε is either fixed heuristically [23, 38], selected by trial and error, or made trainable without a clear understanding of its role [39, 37, 26]. A recent comparative physics-informed study also included Gaussian RBF KANs among several learnable basis families, but treated the Gaussian scale as a fixed preset hyperparameter rather than the subject of a dedicated analysis [38]. In that study, the authors further noted that RBFs may suffer from low accuracy and vanishing gradients when the scaling factor is not chosen appropriately, again highlighting that basis scaling should be handled carefully for stable training. In classical Gaussian RBF interpolation, the scale parameter governs the balance between approximation power and numerical stability, and poor choices can lead either to weak approximation or to severe ill-conditioning [34, 28, 29, 30, 31, 32]. This is closely related to the stability–accuracy trade-off and the uncertainty-principle viewpoint in the RBF literature [33, 34, 40], and the general optimal-scaling problem remains open even in classical meshfree methods. This difficulty in choosing the scale parameter, together with the associated conditioning issues, is one reason why strong-form meshfree RBF methods, despite their potentially much higher accuracy [41], have not replaced the more stable weak-form finite element framework. However, as Larsson and Schaback emphasized, the strong dependence of radial basis function techniques on scaling is a feature, not a bug [28]. Rather than treating scale sensitivity as a defect, it can be viewed as a mechanism that, when understood properly, allows one to control approximation behavior. Noorizadegan et al. [42, 43, 44], inspired by Schaback’s conditioning–accuracy trade-off, showed that working within a numerically safe regime, while pushing the conditioning as high as safely allowed by the implementation and the available reliable digits [45], can lead to good accuracy. One step further in this direction is the flat-limit regime [46, 47]. We do not consider that regime here, and instead keep the analysis within the simple Gaussian-KAN setting used in current practice. The present work does not aim to introduce a new approximation theorem for Gaussian radial basis functions. Rather, it studies a question that is specific to Gaussian KANs: how the Gaussian scale parameter behaves inside a deep edge-based architecture, how this behavior is reflected in first-layer conditioning, and how it can be used to guide practical model design. In this sense, the paper translates classical RBF scaling ideas into a neural setting where the scale parameter affects not only approximation quality, but also representation collapse, trainability, and architectural robustness across layers. Our main contributions are as follows. First, we identify the first layer as the decisive layer for scale selection in Gaussian KANs, both theoretically and experimentally, and show that its conditioning tracks the error behavior of the full network more closely than the deeper layers. Second, we formulate the first Gaussian KAN layer in kernel language and use this viewpoint to define practical conditioning diagnostics directly at the feature-matrix level. Third, for the standard shared-center Gaussian KAN used in current practice, we obtain a simple operating interval for the scale parameter and show that it remains effective across different target functions, collocation densities, grid sizes, architectures, dimensions, and a representative physics-informed PDE problem. Finally, we show that this interval is useful not only for stable fixed-scale selection, but also for variable-scale constructions, constrained scale training, and low-cost scale search using early training MSE. We also use a matched Chebyshev reference to show that a properly scaled Gaussian KAN is already competitive in accuracy relative to another standard KAN basis. The remainder of this paper is organized as follows. Section 2 introduces the motivation from Kolmogorov superposition and presents the Gaussian KAN formulation, together with the experimental setup and benchmark target functions. Section 3 develops the first-layer kernel viewpoint, establishes the first layer as the main representational bottleneck, and supports this role with layer-wise numerical evidence. Section 4 studies first-layer conditioning and derives a practical interval for the Gaussian scale parameter, then examines how this interval behaves under changes in collocation density, number of Gaussian centers, network architecture, and input dimension, and finally discusses shared variable scales, training-MSE-based scale search, and a representative physics-informed Helmholtz problem. Section 5 concludes the paper, and Section 6 discusses limitations and directions for future work. 2 From Kolmogorov Superposition to Gaussian KAN Formulation 2.1 Kolmogorov Superposition as Motivation for KANs Hilbert’s question. At the 1900 International Congress of Mathematicians, Hilbert posed what later became known as Hilbert’s 13th problem [48]. In modern terms, the question asks whether multivariate continuous functions can be represented by finite superpositions of lower-dimensional functions. Hilbert conjectured that such a reduction should fail in general. Kolmogorov superposition theorem. This conjecture was overturned by the work of Kolmogorov and Arnol’d [49, 50, 51]. Theorem 1 (Kolmogorov superposition theorem [51]). Let d≥2d≥ 2 and let C([0,1]d)C([0,1]^d) denote the space of continuous real-valued functions on [0,1]d[0,1]^d. Then there exist continuous functions θq,p:[0,1]→ℝ _q,p:[0,1] , independent of the target function f, such that every f∈C([0,1]d)f∈ C([0,1]^d) admits a representation of the form f(x1,…,xd)=∑q=02dΨq(∑p=1dθq,p(xp)),f(x_1,…,x_d)= _q=0^2d _q\! ( _p=1^d _q,p(x_p) ), (1) where the outer functions Ψq:ℝ→ℝ _q:R are continuous and depend on f. Equation (1) shows that multivariate dependence can be expressed through coordinate-wise univariate transformations, followed by summation and nonlinear recombination. This structural principle is the main connection to KANs. However, KANs should be viewed as architectures inspired by KST rather than exact realizations of the theorem. A more detailed discussion of the relation between KST and KAN can be found in [52]. 2.2 Gaussian KAN Layer Formulation Consider a layer ℓ with input width nℓn_ and output width nℓ+1n_ +1. Given an input vector x(ℓ)∈ℝnℓx^( ) ^n_ , a KAN layer assigns a univariate function ψj,i(ℓ):ℝ→ℝψ^( )_j,i:R to each edge from input coordinate i to output coordinate j. The output x(ℓ+1)∈ℝnℓ+1x^( +1) ^n_ +1 is defined componentwise by xj(ℓ+1)=∑i=1nℓψj,i(ℓ)(xi(ℓ)),j=1,…,nℓ+1.x^( +1)_j= _i=1^n_ ψ^( )_j,i\! (x^( )_i ), j=1,…,n_ +1. (2) Thus, each input coordinate is processed independently through a univariate map, and the resulting contributions are summed to form each output coordinate. This edge-based structure is the defining feature of KANs and the starting point for the Gaussian construction considered here. In the Gaussian KAN studied in this work, every edge function in (2) is expanded in the same finite Gaussian basis. We fix a shared set of centers =c1,…,cG⊂[0,1]C=\c_1,…,c_G\⊂[0,1] (3) and a shared Gaussian scale ϵ>0.ε>0. (4) These quantities are shared across all layers. The restriction ⊂[0,1]C⊂[0,1] is motivated by the first layer, where the original inputs are scaled to this interval. Since deeper-layer coordinates are not constrained to remain in [0,1][0,1], the Gaussian feature map is defined on all of ℝR: φ:ℝ→ℝG,φ(t)=[exp⁡(−(t−c1)2ϵ2)⋮exp⁡(−(t−cG)2ϵ2)]. :R ^G, (t)= bmatrix \! (- (t-c_1)^2ε^2 )\\ \\ \! (- (t-c_G)^2ε^2 ) bmatrix. (5) where G is the number of centers. For each layer ℓ , each edge function is represented as ψj,i(ℓ)(t)=(wj,i(ℓ))⊤φ(t),wj,i(ℓ)∈ℝG.ψ^( )_j,i(t)= (w^( )_j,i ) (t), w^( )_j,i ^G. (6) Hence, the nonlinear behavior of the layer is carried by the Gaussian feature map, while the dependence on the coefficients remains linear. To write the layer compactly, define the stacked feature vector Φ(x(ℓ))=[φ(x1(ℓ))⋮φ(xnℓ(ℓ))]∈ℝnℓG, \! (x^( ) )= bmatrix \! (x^( )_1 )\\ \\ \! (x^( )_n_ ) bmatrix ^n_ G, (7) and the block coefficient matrix W(ℓ)=[(w1,1(ℓ))⊤⋯(w1,nℓ(ℓ))⊤⋮⋱⋮(wnℓ+1,1(ℓ))⊤⋯(wnℓ+1,nℓ(ℓ))⊤]∈ℝnℓ+1×nℓG.W^( )= bmatrix (w^( )_1,1 ) &·s& (w^( )_1,n_ ) \\ & & \\ (w^( )_n_ +1,1 ) &·s& (w^( )_n_ +1,n_ ) bmatrix ^\,n_ +1× n_ G. (8) Then the Gaussian KAN layer takes the form x(ℓ+1)=W(ℓ)Φ(x(ℓ)).x^( +1)=W^( ) \! (x^( ) ). (9) For later use, define the layer operator F(ℓ):ℝnℓ→ℝnℓ+1,F(ℓ)(z)=W(ℓ)Φ(z).F^( ):R^n_ ^n_ +1, F^( )(z)=W^( ) (z). (10) A deep Gaussian KAN with L layers is then written as f=F(L−1)∘F(L−2)∘⋯∘F(0),f:[0,1]n0→ℝnL.f=F^(L-1) F^(L-2) ·s F^(0), f:[0,1]^n_0 ^n_L. (11) Thus, every layer is built from the same mechanism: a Gaussian feature map followed by a linear coefficient map. The key distinction is that the first layer acts directly on the original input coordinates, whereas deeper layers act on learned intermediate representations. For this reason, the first layer is the natural starting point for the kernel analysis in the next section. 3 First-Layer Kernel Structure and Bottleneck Role 3.1 First-Layer Feature Map and Induced Kernel The first layer is the only layer whose Gaussian basis is evaluated directly on the original input domain. For an input x=(x1,…,xd)⊤∈[0,1]dx=(x_1,…,x_d) ∈[0,1]^d, the first-layer feature representation is obtained by stacking the one-dimensional Gaussian features coordinate-wise: Φ(x)=[φ(x1)⋮φ(xd)]∈ℝdG. (x)= bmatrix (x_1)\\ \\ (x_d) bmatrix ^dG. (12) The first hidden representation is therefore x(1)=W(0)Φ(x),x^(1)=W^(0) (x), (13) where W(0)∈ℝn1×dGW^(0) ^n_1× dG. Thus, the first layer is a finite-dimensional feature map followed by a linear transformation. The associated one-dimensional kernel induced by the Gaussian feature map (5) is k(s,t)=φ(s)⊤φ(t),s,t∈[0,1].k(s,t)= (s) (t), s,t∈[0,1]. (14) Accordingly, the first-layer kernel between two inputs x,x′∈[0,1]dx,x ∈[0,1]^d is K0(x,x′)=Φ(x)⊤Φ(x′).K_0(x,x )= (x) (x ). (15) Using the stacked structure in (12), this becomes K0(x,x′)=∑i=1dφ(xi)⊤φ(xi′)=∑i=1dk(xi,xi′).K_0(x,x )= _i=1^d (x_i) (x_i )= _i=1^dk(x_i,x_i ). (16) Hence the first Gaussian KAN layer induces an additive kernel over the input coordinates. This representation is exact and directly reflects the edge-wise structure of KANs in (2): each coordinate is processed separately, and the resulting contributions are summed. For later conditioning analysis, consider a finite sample set =x1,…,xN⊂[0,1]d.X=\x^1,…,x^N\⊂[0,1]^d. (17) For each coordinate i=1,…,di=1,…,d, define the feature block Φ(i)=[φ(xi1)⊤⋮φ(xiN)⊤]∈ℝN×G. ^(i)= bmatrix \! (x_i^1 ) \\ \\ \! (x_i^N ) bmatrix ^N× G. (18) Stacking these blocks horizontally gives the full first-layer feature matrix Φ=[Φ(1)⋯Φ(d)]∈ℝN×dG. = bmatrix ^(1)&·s& ^(d) bmatrix ^N× dG. (19) The empirical first-layer kernel matrix is then K0=ΦΦ⊤∈ℝN×N,K_0= ^N× N, (20) and, by the block structure of Φ , K0=∑i=1dΦ(i)(Φ(i))⊤.K_0= _i=1^d ^(i) ( ^(i) ) . (21) Thus, the empirical first-layer kernel matrix is the sum of the coordinate-wise Gram matrices. This matrix is the starting point for the later analysis of rank, conditioning, and admissible Gaussian scales. 3.2 The First Layer as a Representational Bottleneck Because the first layer acts directly on the original variables, any loss of distinguishability introduced at this stage is inherited by all later layers. Writing the full network as f=T∘G,G(x)=W(0)Φ(x),f=T G, G(x)=W^(0) (x), (22) where T=F(L−1)∘⋯∘F(1)T=F^(L-1) ·s F^(1) denotes the composition of the remaining layers, makes this dependence explicit. Proposition 2 (First-layer bottleneck). Let f be the Gaussian KAN in (22). Then: (i) If two inputs x,x′∈[0,1]dx,x ∈[0,1]^d satisfy G(x)=G(x′),G(x)=G(x ), (23) then f(x)=f(x′).f(x)=f(x ). (24) (i) If two inputs x,x′∈[0,1]dx,x ∈[0,1]^d satisfy Φ(x)=Φ(x′), (x)= (x ), (25) then they induce the same first-layer kernel section, K0(x,z)=K0(x′,z),z∈[0,1]d,K_0(x,z)=K_0(x ,z), z∈[0,1]^d, (26) and therefore f(x)=f(x′)f(x)=f(x ) (27) for every choice of W(0)W^(0) and every downstream map T. (i) On a finite sample set ⊂[0,1]dX⊂[0,1]^d, as ϵ→∞ε→∞, Φ(x)→dGfor every x∈[0,1]d, (x) 1_dG every x∈[0,1]^d, (28) and hence Φ→NdG⊤. 1_N1_dG . (29) Consequently, K0=ΦΦ⊤→dG 1NN⊤,K_0= → dG\,1_N1_N , (30) so the empirical first-layer kernel collapses to a rank-one limit. Moreover, x(1)=W(0)Φ(x)→W(0)dG,x^(1)=W^(0) (x)→ W^(0)1_dG, (31) which means that the first hidden representation becomes constant across the sample set. Proposition 2 formalizes the special role of the first layer. Once two inputs become indistinguishable after the first-layer feature map, no later layer can separate them. In particular, large values of the shared Gaussian scale ϵε force the first-layer features toward a rank-one limit, causing collapse of both the empirical kernel matrix and the representation itself. For this reason, the numerical behavior of the Gaussian KAN is governed first and foremost by the geometry of the first-layer feature matrix. 3.3 Experimental Setup and Benchmark Target Functions Unless stated otherwise, we use a Gaussian KAN with architecture [2,12,12,1],[2,12,12,1], that is, a two-dimensional input, two hidden layers of width 1212, and a one-dimensional output. The Gaussian scale parameter ϵε is swept over 100100 logarithmically spaced values in the interval ϵ∈[5×10−3, 5].ε∈[5× 10^-3,\,5]. The learning rate is fixed at 10−3.10^-3. This dense sweep is used as a brute-force reference, so that the proposed conditioning-based interval can be evaluated against the best empirical choices observed within the tested search range. Training points are sampled on [0,1]2[0,1]^2 using Halton sequences; see Fasshauer [40] and the discussion in Fasshauer and Iske [47]. Halton points are low-discrepancy sequences that provide a more uniform coverage of the domain than purely random samples. Such point sets are widely used in radial basis function meshfree methods. For each experiment, 44 or 55 random seeds are used, accounting for both the random initialization of the coefficients and the point distribution. In the plots, the solid curve shows the geometric mean of the error over the seeds in log scale, while the shaded region indicates the minimum and maximum values across seeds. As the main error measure, we use the validation root mean square error (RMSE), defined by RMSE=(1M∑i=1M(u^(xi,yi)−u(xi,yi))2)1/2,RMSE= ( 1M _i=1^M ( u(x_i,y_i)-u(x_i,y_i) )^2 )^1/2, (32) where (xi,yi)i=1M\(x_i,y_i)\_i=1^M are the validation points, u is the target function, and u u is the network prediction. We consider four benchmark target functions, denoted by F1, F2, F3, and F4. They were chosen to represent smooth localized structure, periodic oscillation, sharp variation, and discontinuity. The functions F1, F2, and F4 are defined on (x,y)∈[0,1]2(x,y)∈[0,1]^2, whereas F3 is defined on (x,y)∈[−1,1]2(x,y)∈[-1,1]^2. • F1: smooth multi-bump function F1(x,y)= F1(x,y)= 0.75exp⁡(−(9x−2)2+(9y−2)24)+0.75exp⁡(−(9x+1)249−9y+110) 75 \! (- (9x-2)^2+(9y-2)^24 )+75 \! (- (9x+1)^249- 9y+110 ) (33) +0.5exp⁡(−(9x−7)2+(9y−3)24)−0.2exp⁡(−((9x−4)2+(9y−7)2)). +5 \! (- (9x-7)^2+(9y-3)^24 )-2 \! (- ((9x-4)^2+(9y-7)^2 ) ). • F2: periodic oscillatory function F2(x,y)=sin⁡(4πx)sin⁡(4πy).F2(x,y)= (4π x) (4π y). (34) • F3: smooth function with a sharp localized peak F3(x,y)=11+103(x2−0.25)2(y2−0.25)2,(x,y)∈[−1,1]2.F3(x,y)= 11+10^3(x^2-0.25)^2(y^2-0.25)^2, (x,y)∈[-1,1]^2. (35) • F4: discontinuous piecewise oscillatory function F4(x,y)=g(x)(1+0.15sin⁡(2πy)),F4(x,y)=g(x) (1+0.15 (2π y) ), (36) where g(x)=5+∑k=14sin⁡(2πkx),x<0.5,cos⁡(20πx),x≥0.5.g(x)= cases5+ _k=1^4 (2π kx),&x<0.5,\\[6.0pt] (20π x),&x≥ 0.5. cases (37) Figure 1 shows the surfaces of these four target functions. Figure 1: Three-dimensional surfaces of the four target functions used in the regression experiments: F1 (smooth multi-bump), F2 (periodic oscillatory), F3 (sharp localized peak), and F4 (discontinuous piecewise oscillatory). These examples span different approximation challenges, including smooth variation, repeated oscillation, strong localization, and discontinuity. 3.4 Numerical Verification of First-Layer Dominance We now provide numerical evidence for Proposition 2. (i) Layer-wise sensitivity of the Gaussian scale parameter. We first study how the validation error changes when the Gaussian scale is varied in only one layer at a time. For the Gaussian KAN used in Figure 2, the architecture is [2,12,12,1][2,12,12,1], with a two-dimensional input, two hidden layers of width 1212, and a one-dimensional output. This gives rise to three successive Gaussian mappings, ϵ1 for 2→12,ϵ2 for 12→12,ϵ3 for 12→1, _1 for 2→ 12, _2 for 12→ 12, _3 for 12→ 1, so three scale parameters are involved, one for each applied layer. Let ϵ=(ϵ1,ϵ2,ϵ3). ε=( _1, _2, _3). (38) We fix the non-varied layers at the reference value 0.10.1 and compare the following schedules: Case 1: (ϵ,0.1,0.1),Case 2: (0.1,ϵ,0.1),Case 3: (0.1,0.1,ϵ),Shared-ϵ case: (ϵ,ϵ,ϵ).Case 1: (ε,0.1,0.1), 2: (0.1,ε,0.1), 3: (0.1,0.1,ε), -ε case: (ε,ε,ε). (39) That is, we vary the scale in only one layer while keeping the others fixed at a good reference value, and compare the resulting response with the fully shared-ϵε case. This allows us to identify the dominant layer, namely the layer whose scale variation has the strongest effect on the overall accuracy. Figure 2 shows the validation RMSE as a function of ϵε for the smooth target F1F1 and the discontinuous target F4F4. In both cases, the first-layer sweep (ϵ,0.1,0.1)(ε,0.1,0.1) tracks the shared-parameter sweep (ϵ,ϵ,ϵ)(ε,ε,ε) much more closely than the sweeps obtained by varying only the deeper layers. This shows that changing the first-layer scale has a much stronger effect on accuracy than changing the scales in the later layers. We also observe that a single good scale can work well across all layers, so there is no clear need to select different scales for different layers. Thus, for this architecture, the dominant sensitivity to the Gaussian scale is carried by the first layer. (i) Conditioning and controlled layer collapse. We next examine the same phenomenon through conditioning and explicit layer collapse. Here the architecture is [2,12,12,12,1],[2,12,12,12,1], so the network contains four successive Gaussian layers. To isolate the effect of representation loss, we force one layer at a time into the large-ϵε collapse regime described in Proposition 2, while keeping the other layers at the fixed value ϵ=0.1ε=0.1. The left panel of Figure 3 shows the resulting validation RMSE for the smooth target F1F1. The no-collapse RMSE corresponds to the case in which no layer is collapsed and all layers use ϵ=0.1ε=0.1. Collapsing an early layer causes the largest degradation, with the most severe effect occurring when the first layer is collapsed. By contrast, collapse in later layers is less harmful, especially in the deepest layers. The right panel of Figure 3 plots the spectral condition number κ(Φℓ)κ( _ ) of each layer as a function of ϵε for the uncollapsed network, that is, without any forced collapse. The first layer is the only one that exhibits a clear low-conditioning regime, whereas the deeper layers remain highly ill-conditioned throughout the sweep. The apparent drop in the condition numbers of Layers 22 to 44 after about ϵ≈0.3ε≈ 0.3 is a finite-precision artifact in float32, not a genuine improvement in conditioning. It occurs because the smallest singular values are no longer reliably resolved at that precision. In higher precision, such as float64, these layers remain highly ill-conditioned and the artificial drop is delayed or disappears. Together, Figures 2 and 3 support Proposition 2 and confirm that the first layer is the primary bottleneck in Gaussian KANs: it is the layer most sensitive to the Gaussian scale, and any loss introduced there has the strongest effect on the final accuracy. Figure 2: Layer-wise sensitivity of the validation RMSE to the Gaussian scale for the architecture 2→12→12→12→ 12→ 12→ 1. The reference value is ϵ0=0.1 _0=0.1, and the four schedules are (ϵ,ϵ0,ϵ0)(ε, _0, _0), (ϵ0,ϵ,ϵ0)( _0,ε, _0), (ϵ0,ϵ0,ϵ)( _0, _0,ε), and (ϵ,ϵ,ϵ)(ε,ε,ε). For both F1F1 and F4F4, varying the first-layer scale (ϵ,ϵ0,ϵ0)(ε, _0, _0) tracks the all-layer response (ϵ,ϵ,ϵ)(ε,ε,ε) much more closely than varying only the deeper layers. This confirms that the first layer carries the dominant sensitivity to the Gaussian scale. Figure 3: Numerical evidence for first-layer dominance for the architecture [2,12,12,12,1][2,12,12,12,1]. Left: validation RMSE after forcing one layer at a time into the large-ϵε collapse regime, while keeping the remaining layers fixed at ϵ=0.1ε=0.1. The no-collapse case corresponds to using ϵ=0.1ε=0.1 in all layers. Collapsing earlier layers, especially the first one, causes the largest increase in error. Right: spectral condition number κ(Φℓ)κ( _ ) of the feature matrix in each layer versus ϵε for the uncollapsed network. The first layer is the only one with a clear low-conditioning regime. 4 First-Layer Conditioning and a Practical Gaussian Scale Interval The goal of this section is to identify a numerically stable operating regime for the first-layer Gaussian feature construction and to translate that regime into a practical rule for choosing the Gaussian scale parameter. Since the empirical first-layer kernel matrix in (20) is generated by the pre-weighting feature matrix Φ∈ℝN×dG, ^N× dG, defined in (19), the relevant object for stability analysis is the feature matrix Φ itself, or equivalently its coordinate-wise blocks Φ(i) ^(i), rather than the learned coefficient matrix applied afterward. 4.1 Conditioning Diagnostics for the First-Layer Feature Matrix To see how the Gram matrix arises, consider one output component and let y∈ℝNy ^N denote the target values on the sample set X. For a fixed first-layer feature matrix Φ , fitting coefficients by least squares gives L(w)=12‖Φw−y‖22,w∈ℝdG.L(w)= 12\| w-y\|_2^2, w ^dG. (40) The corresponding normal equations are Φ⊤Φw=Φ⊤y, \,w= y, (41) so the Gram matrix Φ⊤Φ appears naturally in coefficient recovery. However, column dependence, numerical rank loss, and basis degeneracy originate in Φ before they are inherited and amplified by Φ⊤Φ . Let σ1(Φ)≥σ2(Φ)≥⋯≥σr(Φ)>0 _1( )≥ _2( )≥·s≥ _r( )>0 (42) denote the positive singular values of Φ , where r=rank⁡(Φ)r=rank( ). Define the spectral condition number κ(Φ)=σmax(Φ)σmin(Φ)=σ1(Φ)σr(Φ),κ( )= _ ( ) _ ( )= _1( ) _r( ), (43) when Φ has full column rank, and κ(Φ)=∞κ( )=∞ otherwise. Similarly, κ(Φ⊤Φ)=λmax(Φ⊤Φ)λmin(Φ⊤Φ)κ( )= _ ( ) _ ( ) (44) when Φ⊤Φ is nonsingular, and κ(Φ⊤Φ)=∞κ( )=∞ otherwise. Proposition 3 (Relation between feature-matrix and Gram-matrix conditioning). Let Φ∈ℝN×M ^N× M with M=dGM=dG. Then: (i) If Φ has full column rank, then κ(Φ⊤Φ)=κ(Φ)2.κ( )=κ( )^2. (45) (i) If Φ is rank deficient, then Φ⊤Φ is singular and κ(Φ)=κ(Φ⊤Φ)=∞.κ( )=κ( )=∞. (46) Proof. The eigenvalues of Φ⊤Φ are the squared singular values of Φ , namely λj(Φ⊤Φ)=σj(Φ)2. _j( )= _j( )^2. Hence, if Φ has full column rank, κ(Φ⊤Φ)=σ1(Φ)2σr(Φ)2=κ(Φ)2.κ( )= _1( )^2 _r( )^2=κ( )^2. If Φ is rank deficient, then some singular value of Φ is zero, so Φ⊤Φ has a zero eigenvalue and is singular as well. ∎ Proposition 3 shows that the Gram matrix does not introduce a new instability mechanism; it squares the instability already present in the feature matrix. For this reason, the primary diagnostic for the first-layer Gaussian feature construction is κ(Φ)κ( ), not κ(Φ⊤Φ)κ( ). The numerical quality of the first-layer Gaussian basis is determined by the extreme singular values σmax(Φ)=σ1(Φ),σmin(Φ)=σr(Φ). _ ( )= _1( ), _ ( )= _r( ). (47) To assess numerical rank, we use the standard floating-point tolerance tol=max⁡(N,M)σmax(Φ)εmach,M=dG.tol= (N,M)\, _ ( )\, _mach, M=dG. (48) Accordingly, the first-layer feature matrix is numerically full rank precisely when σmin(Φ)>tol. _ ( )>tol. (49) In float32 arithmetic, a practical additional requirement is that the feature matrix should remain moderately conditioned. Motivated by Proposition 3, we impose κ(Φ)<3×103.κ( )<3× 10^3. (50) This corresponds to κ(Φ⊤Φ)<9×106,κ( )<9× 10^6, which remains within a conservative numerical range for float32-based training. Combining (49) and (50), we define the first-layer diagnostic criteria as σmin(Φ)>tolandκ(Φ)<3×103. _ ( )>tol κ( )<3× 10^3. (51) The first condition detects numerical rank loss, while the second controls near-dependence before rank is lost. Together, they identify the stable operating regime of the first-layer Gaussian feature construction. These diagnostics may be applied either to the full matrix Φ or to the coordinate-wise blocks Φ(i) ^(i), but in all cases the assessment is performed before weighting. 4.2 Derivation of a Practical Gaussian Scale Interval Our goal is to identify a practically useful interval for the Gaussian scale parameter ϵε, of the form ϵ∈[1G−1,ϵcond(N,G)],ε∈ [ 1G-1,\, _cond(N,G) ], (52) where G is the number of Gaussian centers per coordinate. The lower bound 1G−1 1G-1 is the center spacing, while the upper bound ϵcond(N,G) _cond(N,G) is chosen from first-layer conditioning. The reason is that, as ϵε becomes too large, the Gaussian columns become increasingly similar, the first-layer feature matrix becomes ill-conditioned, and, by Proposition 2, the representation moves toward collapse. We then simplify ϵcond(N,G) _cond(N,G) into a rule depending only on G, leading to an explicit interval that is easy to use in practice. The lower endpoint ϵmin=1G−1 _ = 1G-1 (53) is the spacing between adjacent uniformly distributed centers. This is the natural geometric scale below which the Gaussian basis becomes too localized. The main task is therefore to determine a practical upper boundary. To do this, we use the first-layer feature matrix Φ(ϵ) (ε) and define two diagnostic markers. First, let ϵκ=arg⁡minϵ>0⁡|log⁡κ(Φ(ϵ))−log⁡(3×103)| _κ= _ε>0 | κ\! ( (ε) )- \! (3× 10^3 ) | (54) denote the Gaussian scale at which the first-layer condition number is closest to the practical threshold. Second, let ϵrank=infϵ>0:σmin(Φ(ϵ))≤tol(ϵ) _rank= \ε>0: _ \! ( (ε) ) (ε) \ (55) denote the onset of numerical rank loss. Figure 4 compares these two markers with the validation RMSE and the first-layer condition number. Across all four test cases, ϵκ _κ and ϵrank _rank occur very close to one another. Thus, for the first-layer Gaussian feature matrix, the practical onset of ill-conditioning and the onset of numerical rank loss are nearly simultaneous. As an additional accuracy reference, Figure 4 also includes the validation RMSE obtained by a Chebyshev KAN of the same architecture. In that comparison, the Chebyshev degree is chosen as m=G−1,m=G-1, so that each Chebyshev edge expansion has m+1=Gm+1=G basis coefficients, matching the G Gaussian coefficients used per edge in the Gaussian KAN. The Chebyshev result is shown only as a reference accuracy level, not as part of the first-layer conditioning analysis. Figure 4: Empirical localization of the conditioning-based interval. Each panel shows the validation RMSE, the first-layer condition number κ(Φ)κ( ), the lower geometric scale ϵ=1/(G−1)ε=1/(G-1), the conditioning-based marker ϵκ _κ, and the rank-loss marker ϵrank _rank defined by σmin(Φ)≤tol _ ( ) . Across all cases, ϵκ _κ and ϵrank _rank are nearly coincident, and the minimum validation RMSE occurs near this transition. This identifies the upper boundary of the useful scale interval. The horizontal purple line in each panel marks the validation RMSE of the corresponding Chebyshev KAN reference model, with degree m=G−1m=G-1, so that the number of basis coefficients per edge matches the Gaussian case (see paragraph 4.2 for the Chebyshev discussion.) The same figure also shows that the minimum validation RMSE occurs near this transition region. For smaller values of ϵε, the basis is too localized and the error remains relatively large. As ϵε increases toward ϵκ _κ and ϵrank _rank, neighboring Gaussian features overlap sufficiently to improve approximation while the feature matrix remains numerically stable. Beyond this regime, κ(Φ)κ( ) grows rapidly, the numerical full-rank condition fails, and the validation error deteriorates. In most cases, the best-performing region lies between the geometric scale 1/(G−1)1/(G-1) and the conditioning boundary. This motivates defining the upper endpoint in (52) by ϵcond(N,G)=arg⁡minϵ>0⁡|log⁡κ(Φ(ϵ))−log⁡(3×103)|. _cond(N,G)= _ε>0 | κ\! ( (ε) )- \! (3× 10^3 ) |. (56) That is, ϵcond(N,G) _cond(N,G) is the value of ϵε at which the first-layer feature matrix reaches the practical conditioning threshold. Hence the initial conditioning-based interval is ϵ∈[1G−1,ϵcond(N,G)].ε∈ [ 1G-1,\, _cond(N,G) ]. (57) The next question is whether ϵcond(N,G) _cond(N,G) can be simplified. For each pair (N,G)(N,G), we compute ϵcond(N,G) _cond(N,G) by sweeping ϵε and selecting the value whose condition number is closest to the threshold. In the regression study, we use N∈[400,3600]N∈[400,3600] and G∈[6,30]G∈[6,30], with multiple low-discrepancy sample realizations for each configuration. The resulting values are aggregated across seeds by a geometric mean, yielding one representative conditioning scale for each pair (N,G)(N,G). The full set of pairs is then divided into training and validation subsets using a 70%/30%70\%/30\% split. To identify the dominant dependence of ϵcond(N,G) _cond(N,G), we fit the log-linear model log⁡ϵcond=log⁡C−qlog⁡(G−1), _cond= C-q (G-1), (58) or equivalently, ϵcond≈C(G−1)−q. _cond≈ C\,(G-1)^-q. (59) The purpose of this fit is to test whether the conditioning boundary is governed mainly by the center spacing and only weakly by the number of samples. The regression results show that the dominant dependence of ϵcond _cond is on G, while the influence of N is weak once the samples cover [0,1][0,1] sufficiently densely. This is consistent with the fact that the overlap of neighboring Gaussian columns is controlled primarily by the center spacing 1/(G−1)1/(G-1). Fitting (59) on the training set yields ϵcond≈2.8127(G−1)−1.1297, _cond≈ 2.8127\,(G-1)^-1.1297, (60) with validation coefficient of determination R2≈0.997R^2≈ 0.997. Since the fitted exponent is close to 11, this rule can be simplified to ϵcond≈2G−1, _cond≈ 2G-1, (61) which remains highly accurate on the validation set, with R2≈0.985R^2≈ 0.985. To isolate the dependence on G, the validation values of ϵcond _cond are also aggregated over all samples sharing the same number of centers, using a geometric mean across N. The resulting curve follows (61) closely across the full tested range of G, again showing that the principal dependence is on the center spacing rather than on the sample count. We therefore replace the conditioning-based interval (57) by the explicit practical rule ϵ∈[1G−1,2G−1].ε∈ [ 1G-1, 2G-1 ]. (62) The lower endpoint keeps the Gaussian scale comparable to the spacing of adjacent centers, while the upper endpoint places the first-layer feature matrix near the empirically observed conditioning boundary. Figure 5 validates this simplification. The fitted power law gives an excellent approximation to the measured values of ϵcond _cond, and the simplified rule 2/(G−1)2/(G-1) remains highly accurate on held-out data. Thus, little is lost by replacing the fitted law with the simpler expression (61). The conclusion of this section is therefore direct: the useful operating range for ϵε is ϵ∈[1G−1,2G−1],ε∈ [ 1G-1,\, 2G-1 ], and this range is determined primarily by the number of centers G, not by the sample count N. Figure 5: Simplification of the conditioning-based upper bound ϵcond(N,G) _cond(N,G). Left: predicted versus measured values of ϵcond _cond on held-out (N,G)(N,G) pairs. The fitted power law (60) achieves R2≈0.997R^2≈ 0.997, while the simplified rule (61) remains highly accurate with R2≈0.985R^2≈ 0.985. The dashed line denotes exact agreement. Right: geometric-mean aggregation of the measured values over validation samples sharing the same G. The empirical curve is well approximated by (61), showing that the conditioning-based upper bound is governed primarily by the number of centers. As a final remark, one may also define the conditioning-based upper bound for ϵε using classical bounds from meshfree RBF theory, for example through estimates involving both the smallest and largest singular values of the feature matrix [40]. Such bounds are theoretically meaningful, but in the present setting they are too impractical for routine use. Since our aim here is a rule that is both simple and usable in practice, we do not pursue that direction. Chebyshev reference. With the practical interval now identified, we return briefly to Figure 4 to comment on the Chebyshev reference line. The Chebyshev implementation is included only as an additional accuracy benchmark. For a fair comparison, the Chebyshev degree is chosen as m=G−1m=G-1, so that each Chebyshev edge expansion has m+1=Gm+1=G trainable coefficients, matching the G Gaussian basis coefficients per edge. In this reference model, tanh is used between layers to keep the inputs within the admissible range of the Chebyshev basis, and the final layer also uses a Chebyshev expansion rather than a linear head. Under this matched per-edge coefficient budget, Figure 4 shows that a suitably scaled Gaussian KAN can attain validation errors comparable to, and in some cases lower than, the corresponding Chebyshev KAN. This comparison also highlights a practical difference between the two representations. In the Gaussian KAN, conditioning and accuracy can be adjusted directly through the centers and the scale parameter ϵε. In contrast, the Chebyshev KAN is tied more rigidly to the polynomial degree: increasing the degree may improve approximation power, but often at the cost of reduced stability (see Figure 4), whereas lower degrees are usually more stable but may sacrifice accuracy. Various heuristic remedies have therefore been introduced in the Chebyshev literature, including localization of the basis [7], interlayer normalization [53], and architectural modifications such as additional tanh nonlinearities or alternative readout layers [6]. By comparison, the Gaussian KAN offers a cleaner mechanism in which the same scale parameter directly influences both conditioning and accuracy. The role of the present comparison is to show that the Gaussian KAN, when equipped with a properly chosen scale, can already be as accurate as a standard Chebyshev KAN. This supports the main purpose of the present work, namely, to understand how the Gaussian scale parameter ϵε governs both conditioning and accuracy, and how it can be selected in a principled and practically useful way. A definitive comparison between Chebyshev and Gaussian KANs would require a separate study based on carefully tuned models under matched architectures, training strategies, and stabilization techniques. 4.3 Effect of Collocation Density N We now study how the number of collocation points affects approximation in the Gaussian KAN. The goal is to separate two effects: N controls how well the domain is resolved, while ϵε controls the quality and stability of the Gaussian feature representation. We show that increasing N improves accuracy, but only when ϵε is chosen in a suitable range. All experiments in this subsection are performed on the unit square Ω=[0,1]2, =[0,1]^2, (63) using Halton points XN=x1,…,xN⊂ΩX_N=\x^1,…,x^N\⊂ (64) as the training set. Since Halton points are quasi-uniform, their geometric resolution is described by the fill distance hXN,Ω=supx∈Ωminxn∈XN⁡‖x−xn‖2,h_X_N, = _x∈ _x^n∈ X_N\|x-x^n\|_2, (65) which in two dimensions satisfies hXN,Ω∼N−1/2.h_X_N, N^-1/2. (66) Thus, increasing N improves the geometric resolution of the training set. To generate a systematic refinement sequence, we choose N∈k2:k=5,7,9,…,87.N∈\k^2:\;k=5,7,9,…,87\. (67) This gives a monotone progression of collocation densities and makes the scaling in (66) easy to interpret. Throughout this subsection, the grid size is fixed at G=20G=20. Therefore, the practical interval from (62) becomes ϵ∈[119,219]≈[0.0526, 0.1053].ε∈ [ 119, 219 ]≈[0.0526,\,0.1053]. (68) We now examine how this interval behaves as N changes. Error landscape in ϵε for different N. For fixed architecture and grid size, let E(N,ϵ)E(N,ε) (69) denote the validation RMSE obtained with N collocation points and Gaussian scale ϵε. Figure 6 plots E(N,ϵ)E(N,ε) for two representative training sizes, N=961N=961 and N=3249N=3249, with the vertical lines marking the interval in (68). Across all four targets, the same basic structure appears: the error is large for very small ϵε, decreases into a low-error region, and rises again for large ϵε. Thus the dependence on ϵε remains essentially U-shaped as N changes. Increasing N lowers the overall error level, but shifts the useful ϵε-window only mildly. In particular, the error-minimizing region remains inside or close to (68). This shows that the conditioning-based interval remains informative even as the collocation density changes. Best attainable error as a function of N. To isolate the effect of collocation density, define Eopt(N)=minϵ⁡E(N,ϵ),E_opt(N)= _εE(N,ε), (70) that is, the best validation RMSE obtained over the ϵε-sweep for a given N. Figure 7 plots Eopt(N)E_opt(N) against N. The best error decreases monotonically for all four targets, confirming that denser collocation improves approximation quality once the Gaussian scale is chosen well. The observed decay also reflects target regularity. For the smoother targets F1 and F2, the visible pre-asymptotic behavior is approximately aligned with Eopt(N)∼h3∼N−3/2,E_opt(N) h^3 N^-3/2, (71) whereas the less regular targets F3 and F4 are closer to Eopt(N)∼h2∼N−1.E_opt(N) h^2 N^-1. (72) This is consistent with the standard expectation that smoother targets admit faster convergence, while reduced regularity lowers the observed rate; see Chapter 17 of Fasshauer [40]. Interaction between N and ϵε. Figure 8 shows the validation RMSE versus N for the fixed Gaussian scales ϵ∈0.06,0.08,0.1ε∈\0.06,0.08,0.1\. In all cases, the error decreases as N increases, but the decay still depends on the choice of ϵε. Thus, N and ϵε play different roles: N increases geometric resolution, while ϵε controls how effectively the Gaussian feature representation uses that resolution. Figure 6: Validation RMSE as a function of the Gaussian scale ϵε for two representative collocation densities, N=729N=729 and N=3249N=3249, with G=20G=20. The vertical markers indicate ϵ=1/(G−1)ε=1/(G-1) and ϵ=2/(G−1)ε=2/(G-1). Across all four targets, increasing N lowers the error, while the low-error region remains concentrated near the conditioning-based interval. Figure 7: Best validation RMSE Eopt(N)E_opt(N) from (70) as a function of the number of collocation points. For the smoother targets F1 and F2, the visible decay is approximately aligned with the h3h^3 reference, while F3 and F4 are closer to the h2h^2 reference. In all cases, denser Halton sampling improves the best attainable accuracy. Figure 8: Validation RMSE versus N for the fixed Gaussian scales ϵ=0.06,0.08ε=0.06,0.08, together with the near-optimal choice at each N with G=20G=20. 4.4 Effect of the Number of Gaussian Centers G We next study the effect of the grid resolution, that is, the number of shared Gaussian centers in each one-dimensional feature map. For fixed architecture and training size N, let G=c1,…,cG⊂[0,1]C_G=\c_1,…,c_G\⊂[0,1] (73) denote the uniformly distributed centers. The conditioning-based rule in (61) and (62) already shows that G and ϵε cannot be chosen independently. As G increases, the center spacing Δc _c decreases, so the useful scale range must also move to smaller ϵε. Fixed ϵε versus adaptive ϵε. Figure 9 shows E(G,ϵ)E(G,ε) for several fixed values of ϵε, together with a near-optimal choice selected separately for each G. The fixed-ϵε curves depend strongly on the chosen scale: a value that works well for one grid size can be clearly suboptimal for another. By contrast, the near-optimal choice gives a much more stable error curve across the full range of G. Thus, refining the grid alone is not enough; ϵε must be adjusted with the center spacing. Shift of the low-error ϵε-window. Figure 10 compares the full ϵε-sweeps for two representative grid sizes, G=10G=10 and G=32G=32. Both cases show the same U-shaped dependence on ϵε, but the effect of increasing G is different from the effect of increasing N. Increasing N mainly pushes the curve downward, reducing the error level while leaving the useful ϵε-window roughly in the same place. Increasing G, instead, shifts the low-error region to the left, that is, toward smaller ϵε, and also makes the good-ϵε range wider. This is exactly what the proposed range predicts, since that range is tied to the center spacing. At the same time, the minimum error itself changes much less than the location of the minimizer. This is consistent with the matrix viewpoint developed earlier: increasing G mainly changes the geometry of the basis and the position of the stable operating range, rather than acting as a direct accuracy parameter by itself. Figure 9: Validation RMSE versus grid resolution G for several fixed Gaussian scales and for a near-optimal scale chosen separately at each G. The fixed-ϵε curves vary strongly with G, whereas the near-optimal choice yields a much more stable error profile. Figure 10: Validation RMSE as a function of the Gaussian scale ϵε for two representative grid resolutions, G=10G=10 and G=32G=32. As G increases, the low-error region shifts to smaller ϵε and becomes wider. 4.5 Effect of Network Architecture We next study how the network architecture affects the relation between approximation error and the Gaussian scale. Throughout this subsection, the grid resolution is fixed at G=20G=20, so the conditioning-based reference interval is ϵ∈[1/19, 2/19]≈[0.0526, 0.1053]ε∈[1/19,\,2/19]≈[0.0526,\,0.1053]. Figure 11 shows that all tested architectures retain the same U-shaped dependence on ϵε. What changes is the minimum error and the width of the low-error region. Since the first layer is the only layer defined directly on the input domain, it still determines the relevant scale regime. Increasing width or depth does not move the left side of the U-shaped curve, which is the part controlled by the first-layer geometry. Instead, larger architectures mainly enlarge the range of good ϵε from the right. This is different from increasing G, which shifts the useful range to the left. The effect is more visible for the less smooth targets F3 and F4. In the top row, the wider architecture (2,20,20,1)(2,20,20,1) generally attains lower error than (2,4,4,1)(2,4,4,1) and makes the low-error region broader. In the bottom row, the deeper architecture (2,8,8,8,1)(2,8,8,8,1) also lowers the error and broadens the admissible range compared with (2,8,1)(2,8,1), again mainly on the larger-ϵε side. Thus, a scale range that is suitable for a shallower Gaussian representation remains admissible for a deeper Gaussian KAN; deeper layers may enlarge the useful range, but they do not shrink the range already set by the first layer. Figure 11: Validation RMSE versus the Gaussian scale ϵε for different Gaussian KAN architectures at fixed grid size G=20G=20. The two vertical markers indicate ϵ=1/(G−1)ε=1/(G-1) and ϵ=2/(G−1)ε=2/(G-1). 4.6 Effect of Input Dimension We finally investigate how the Gaussian scale behaves as the input dimension increases. For each dimension d, we consider the smooth target fd(x)=exp⁡(1d∑i=1d[sin⁡(πxi)+12xi2]),x∈[0,1]d,f_d(x)= \! ( 1d _i=1^d [ (π x_i)+ 12x_i^2 ] ), x∈[0,1]^d, (74) which remains bounded and of comparable magnitude across dimensions. The experiments in Figure 12 use the settings listed in Table 1. Table 1: Experimental settings used in the input-dimension study. The last column lists the hidden-layer widths, so for example [12][12] denotes one hidden layer with 12 neurons, while [12,12,12][12,12,12] denotes three hidden layers with 12 neurons each. Dimension d Collocation points N Hidden-layer widths 1 30 [12][12] 2 1000 [12,12][12,12] 3 30000 [12,12,12][12,12,12] 4 50000 [12,12,12,12][12,12,12,12] Figure 12 shows that the validation RMSE keeps the same U-shaped dependence on ϵε from 1D to 4D. Very small ϵε leads to overly localized features, while very large ϵε leads to the flat and numerically unfavorable regime. Most importantly, the low-error region remains close to the same two reference lines so increasing the ambient dimension does not introduce a new preferred scale. This is consistent with the first-layer matrix construction in (19). The feature matrix is built by concatenating coordinate-wise Gaussian blocks, Φ=[Φ(1)⋯Φ(d)], = bmatrix ^(1)&·s& ^(d) bmatrix, (75) so increasing d adds more one-dimensional blocks but does not change the spacing mechanism inside each block. The same shared scale ϵε still controls the overlap between neighboring centers, and therefore the same conditioning-based interval remains relevant. The minimum error does increase with dimension, as expected, because the approximation problem becomes harder. However, the shape of the error curve and the location of the useful ϵε-window remain essentially unchanged. Thus, in the present separable Gaussian KAN, the selection of ϵε is still controlled primarily by the first-layer Gaussian geometry rather than by the ambient dimension itself. Figure 12: Validation RMSE versus the Gaussian scale ϵε for input dimensions d=1,2,3,4d=1,2,3,4. The settings are those listed in Table 1. The two vertical markers indicate the reference scales ϵ=1/(G−1)ε=1/(G-1) and ϵ=2/(G−1)ε=2/(G-1). Across dimensions, the same U-shaped dependence on ϵε is observed, and the low-error region remains close to the same conditioning-based interval. 4.7 Gaussian KAN with Shared Centerwise Variable Scales The fixed-scale Gaussian feature map in (5) can be generalized by assigning a different scale to each center while still sharing the same scale vector across all neurons and all layers. Let ℰ=ϵ1,…,ϵG,ϵg>0,g=1,…,G,E=\ _1,…, _G\, _g>0, g=1,…,G, (76) be the vector of centerwise Gaussian scales associated with the shared centers =c1,…,cG⊂[0,1]C=\c_1,…,c_G\⊂[0,1]. Then the Gaussian feature map becomes φℰ(t)=[exp⁡(−(t−c1)2ϵ12)⋮exp⁡(−(t−cG)2ϵG2)]∈ℝG,t∈ℝ. _E(t)= bmatrix \! (- (t-c_1)^2 _1^2 )\\ \\ \! (- (t-c_G)^2 _G^2 ) bmatrix ^G, t . (77) Thus, each center cgc_g has its own width ϵg _g, while the same vector ℰE is reused throughout the whole network. In the experiments below, the variable scales are sampled from a normal distribution on ϵg∈[ϵmin,ϵmax],ϵmax=2G−1, _g∈[ _ , _ ], _ = 2G-1, (78) and then assigned to the shared centers cgg=1G\c_g\_g=1^G, randomly. Hence, this model keeps the same Gaussian KAN architecture, but introduces controlled nonuniformity in the basis widths across centers. Although the centers themselves remain uniformly distributed, allowing the widths to vary makes the effective coverage of the basis nonuniform across the domain. In this sense, the model gains some of the flexibility usually associated with nonuniform or data-adapted center placement, while still retaining the simplicity and interpretability of a fixed uniform center grid. Figure 13 compares this variable-scale model with the fixed-scale choice ϵ=1/(G−1)ε=1/(G-1). As Figure 13 shows, the shared centerwise variable-scale model can be more accurate than the single-ϵε version. In the present experiment, the scales are assigned randomly within the admissible range, so this result should be interpreted only as an initial demonstration rather than an optimized construction. The use of variable shape parameters has a long history in classical RBF methods, where nonuniform scales often improve accuracy over constant-scale formulations [54, 55, 56], and the present results indicate that this advantage can also carry over to the neural-network setting. More generally, one may vary not only the interval but also the sampling law for the scales, for example by using normal, log-normal, log-uniform, gamma, or chi-square distributions [54]. Thus, even with uniformly distributed centers, variable-scale Gaussian KANs can induce a more adaptive and spatially nonuniform representation, making them a flexible and promising extension of the fixed-scale model. Figure 13: Validation RMSE versus epoch for the fixed-scale Gaussian KAN and the Gaussian KAN with shared variable scales. Left: F1F1 with N=961N=961, where the variable model uses range R5R5. Right: F4F4 with N=1681N=1681. The fixed model uses the common scale ϵ=1/(G−1)ε=1/(G-1), while the variable model assigns one scale ϵg _g to each center cgc_g, with the same scale vector shared across all neurons and all layers. 4.8 Training-MSE-Based Scale Search Within the Admissible Interval The first-layer diagnostics do not select a single value of ϵε; they reduce the search to a numerically admissible interval. This is the main practical role of the range: instead of sweeping over all candidate scales, one only needs to search inside the conditioning-controlled region. Figure 14 shows that, within this reduced interval, the training MSE provides a useful criterion for selecting a working scale. For the four representative problems F1F1, F2F2, F3F3, and F4F4, the low-training-error region at 500500 epochs is located in essentially the same part of the ϵε-axis as the low-validation-RMSE region at 1000010000 epochs. The agreement is clearest for the smooth and oscillatory cases F1F1 and F2F2, and remains informative for the more difficult cases F3F3 and F4F4, where the curves are rougher but still identify the same favorable range. The practical implications are straightforward: • The admissible interval sharply reduces the ϵε-search space. • Inside that interval, training MSE is sufficient to localize a good scale. • This localization already appears after a small number of iterations. • Hence one can perform a short sweep over ϵε, choose a candidate from the training-MSE curve, and run full training only for that value. • This is especially useful when the exact solution is unavailable, so validation against the true target cannot be used for scale selection. Therefore, the proposed interval is useful not only for stability control, but also for reducing the cost of scale selection. It turns the search for ϵε into a smaller optimization problem over a numerically justified range, and the early training MSE is usually sufficient to identify the relevant part of that range. Figure 14: Training-MSE-based scale selection inside the admissible interval. Each panel shows the training MSE at 500500 epochs and the validation RMSE at 1000010000 epochs as functions of ϵε. The vertical markers indicate the lower reference scale and the conditioning-based upper boundary of the search interval. Across all four examples, the early training-MSE curve already localizes essentially the same favorable ϵε-region as the final validation RMSE, so the admissible interval can be searched efficiently using only a short preliminary sweep. 4.9 Helmholtz Problem: Physics-Informed Validation To test whether the conditioning-based scale interval also remains useful in a physics-informed setting, we consider the two-dimensional Helmholtz problem −Δu−λu=fin Ω=[0,1]2,- u-λ u=f =[0,1]^2, (79) with homogeneous Dirichlet boundary condition u=0on ∂Ω.u=0 ∂ . (80) We choose the exact solution u(x,y)=sin⁡(a1πx)sin⁡(a2πy),(x,y)∈[0,1]2,u(x,y)= (a_1π x) (a_2π y), (x,y)∈[0,1]^2, (81) which vanishes on the boundary and therefore satisfies (80). Substituting (81) into (79) gives the forcing term f(x,y)=((a12+a22)π2−λ)sin⁡(a1πx)sin⁡(a2πy).f(x,y)= ((a_1^2+a_2^2)π^2-λ ) (a_1π x) (a_2π y). (82) In the experiments below, the PDE residual is weighted by wpde=1,w_pde=1, (83) and the boundary loss is weighted by wbc=100.w_bc=100. (84) The collocation size and Gaussian grid are fixed at (NBC,NPDE)=(800,2000),G=20,(N_BC,N_PDE)=(800,2000), G=20, (85) so the conditioning-based reference interval from (62) becomes ϵ∈[119,219]≈[0.0526, 0.1053].ε∈ [ 119, 219 ]≈[0.0526,\,0.1053]. (86) We study three values of the Helmholtz parameter, λ∈0,10,100,λ∈\0,10,100\, (87) for two exact solutions: the lower-frequency case (a1,a2)=(1,2),(a_1,a_2)=(1,2), (88) and the anisotropic higher-frequency case (a1,a2)=(1,4).(a_1,a_2)=(1,4). (89) Figure 15 shows the validation RMSE as a function of the Gaussian scale ϵε, with the two vertical markers corresponding to the endpoints of (86). In both test cases, the same qualitative behavior seen in the regression experiments is recovered: the error is large for very small ϵε, decreases sharply into a low-error regime, and then increases again as ϵε becomes too large. Thus, the Helmholtz problem exhibits the same basic trade-off between insufficient overlap at small scales and degradation of the Gaussian representation at large scales. Most importantly, the best-performing region remains within or close to the interval (86) for all three values of λ. For λ=0λ=0 and λ=10λ=10, the minimum error is attained clearly inside this interval in both problem settings. For the more challenging case λ=100λ=100, the error level is higher and the low-error region becomes broader and less stable, but it is still centered near the same conditioning-based range. This indicates that the first-layer scale rule derived earlier remains informative even after the loss is replaced by a PDE residual and boundary penalty. The two exact solutions also illustrate a useful distinction. The case (a1,a2)=(1,2)(a_1,a_2)=(1,2) is smoother and easier, so the minimum RMSE is lower and the low-error region is broader. The case (a1,a2)=(1,4)(a_1,a_2)=(1,4) is more oscillatory in the y-direction, and therefore more demanding, which raises the minimum error and makes the curves more sensitive, especially for larger λ. Nevertheless, across both exact solutions and all tested values of λ, the preferred ϵε-window remains within or close to the same conditioning-based interval and stays tied to the same geometric scale set by the shared Gaussian centers. This shows that larger λ and more oscillatory solutions may increase the error and make the curves less stable, but they do not substantially change the preferred ϵε-window. This example is intentionally limited to a single PDE family with a tunable parameter λ. Its role is not to provide an exhaustive PDE study, but to show that the interval suggested by the first-layer conditioning analysis continues to work in a representative physics-informed problem. A broader investigation across other PDEs, boundary conditions, and training regimes can be carried out in future work. Figure 15: Validation RMSE versus the Gaussian scale ϵε for the Helmholtz problem (79)–(80) with (NBC,NPDE)=(800,2000)(N_BC,N_PDE)=(800,2000), G=20G=20, wpde=1w_pde=1, and wbc=100w_bc=100. The two vertical markers indicate ϵ=1/(G−1)ε=1/(G-1) and ϵ=2/(G−1)ε=2/(G-1), i.e. the endpoints of the conditioning-based interval (86). (a) Exact solution (81) with (a1,a2)=(1,2)(a_1,a_2)=(1,2). (b) Exact solution (81) with (a1,a2)=(1,4)(a_1,a_2)=(1,4). In both cases, and for all tested values λ∈0,10,100λ∈\0,10,100\, the best-performing region lies within or close to the same conditioning-based range. 5 Conclusion To the best of our knowledge, this work is the first to study Gaussian KANs from a kernel and conditioning viewpoint, with emphasis on how the Gaussian scale parameter ϵε governs both numerical stability and approximation accuracy. A matched Chebyshev reference was included as an additional benchmark, and the results show that, when ϵε is chosen properly, a Gaussian KAN can achieve accuracy comparable to, and in some cases better than, the corresponding Chebyshev KAN with the same number of trainable coefficients per edge. This comparison is not intended as a definitive optimized benchmark between two basis families; rather, it shows that Gaussian KANs are already highly competitive once the scale parameter is selected appropriately. The results of this study consistently indicate that scale selection in Gaussian KANs is governed primarily by the first layer, since it is the only layer constructed directly on the input domain and any loss of distinguishability introduced there cannot be recovered by later layers. This has been supported both theoretically and experimentally, with the conditioning of the first layer closely reflecting the error behavior of the full network. From this viewpoint, the relevant stability object is the first-layer feature matrix Φ , rather than primarily the Gram matrix Φ⊤Φ . Using numerical full rank together with a practical conditioning threshold, we identified the interval ϵ∈[1G−1,2G−1],ε∈ [ 1G-1,\, 2G-1 ], which provides a simple and effective admissible range for Gaussian KANs. The experiments further show that this interval remains informative across different target functions, collocation densities N, grid sizes G, network architectures, input dimensions, and also in a physics-informed Helmholtz problem. In all cases, the validation error exhibits a clear U-shaped dependence on ϵε, with the best-performing region lying within or close to the proposed interval. An additional practical observation concerns whether the Gaussian scale should be shared across layers. The layer-wise sweeps show that varying the first-layer scale produces an error curve much closer to that of the fully shared-ϵε network than varying the scale only in deeper layers. This further supports the first layer as the dominant layer for scale selection. At the same time, allowing different scales across layers does not lead to a clear improvement in the best attainable accuracy. Although deeper layers often admit a broader range of acceptable scales, using the same well-chosen ϵε in all layers yields nearly the same optimum. Therefore, a shared Gaussian scale across layers is a sensible choice: it keeps the model simpler while maintaining essentially the same accuracy. The numerical studies also clarify the roles of N and G. Increasing the number of collocation points N improves the best attainable accuracy and yields clear convergence with respect to N, provided that ϵε is chosen in a suitable range. Increasing the number of centers G mainly shifts the low-error region toward smaller ϵε and broadens the admissible range, rather than guaranteeing a substantially smaller optimum by itself. In practice, a relatively small G can sometimes attain nearly the same optimal error as a larger one while being computationally much cheaper, which makes the proposed scale rule useful not only for accuracy but also for efficient model design. Beyond fixed-scale selection, the proposed interval also proved useful in more practical settings. A shared centerwise variable-scale construction defined within this range can improve accuracy over a single fixed ϵε. In addition, the training MSE was shown to be an effective proxy for locating a good scale within the admissible interval, often after only a small number of iterations. This greatly reduces the cost of scale search and is especially important in applications where no exact solution is available, so that validation against the true target cannot be used. Overall, the results show that the Gaussian scale parameter in Gaussian KANs should not be viewed as an arbitrary hyperparameter. Rather, it can be selected in a principled way that simultaneously supports stability, accuracy, and efficient computation. In this sense, the conclusion of Larsson and Schaback that the strong dependence of radial basis function techniques on scaling is a feature, not a bug [28] naturally extends to neural architectures based on Gaussian radial basis functions. 6 Limitations and Future Work This work focuses on first-layer conditioning. Deeper layers may still be ill-conditioned, but our results suggest that this does not necessarily harm performance when ϵε is chosen in the proposed range. In later layers, ill-conditioning may occur because the fixed centers are no longer aligned with the transformed coordinates, or because those coordinates become clustered or unevenly distributed. Improving this would likely require adaptive center relocation and a more complex multilevel training strategy, possibly together with updating collocation points during training. Such extensions would substantially complicate the model, while better conditioning would still not necessarily guarantee better accuracy. For this reason, we retain the standard Gaussian KAN architecture used in current practice. Several directions remain open for future work. One natural extension is to study Gaussian KANs beyond regression and the representative PDE setting considered here, for example in classification and operator-learning problems. It would also be interesting to investigate other RBF-based KANs [57, 58, 59], to examine how the present conditioning-based perspective extends beyond the Gaussian case. Another promising direction is to use the proposed interval as a constrained training region for a learnable ϵε, instead of training the scale blindly over an unrestricted domain. This may improve stability, avoid collapse, and reduce tuning cost. Variable-scale Gaussian KANs also deserve further study. Although a shared centerwise variable-scale model has already been introduced here, the preliminary results suggest that it can outperform the fixed-scale version. This opens the possibility of studying not only the admissible range of variable scales, but also different distribution families and sampling laws, such as chi-square distributions, random perturbations with different variances, or other nonuniform choices. Since variable shape parameters have a long history in classical RBF methods and often improve accuracy over constant-scale formulations [54, 55, 56], this direction appears especially promising. Finally, more PDEs should be tested in order to assess how far the present scale-selection rule extends in physics-informed settings. References [1] Z. Liu, Y. Wang, S. Vaidya, F. Ruehle, J. Halverson, M. Soljačić, T. Y. Hou, and M. Tegmark, “KAN: Kolmogorov-Arnold networks,” arXiv preprint arXiv:2404.19756, 2024. [2] Z. Liu, P. Ma, Y. Wang, W. Matusik, and M. Tegmark, “KAN 2.0: Kolmogorov-Arnold networks meet science,” arXiv preprint arXiv:2408.10205, 2024. [Online]. Available: https://github.com/KindXiaoming/pykan [3] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” J. Comput. Phys., vol. 378, p. 686–707, 2019. [4] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, “Physics-informed machine learning,” Nat. Rev. Phys., vol. 3, no. 6, p. 422–440, 2021. [5] S. Sidharth, K. A. R., and A. K. P., “Chebyshev polynomial-based Kolmogorov-Arnold networks: An efficient architecture for nonlinear function approximation,” arXiv preprint arXiv:2405.07200, 2024. [6] N. A. Daryakenari, K. Shukla, and G. E. Karniadakis, “Representation meets optimization: Training PINNs and PIKANs for gray-box discovery in systems pharmacology,” arXiv preprint arXiv:2504.07379, 2025. [7] J. D. Toscano, L.-L. Wang, and G. E. Karniadakis, “KKANs: Kurkova-Kolmogorov-Arnold networks and their learning dynamics,” Neural Netw., vol. 191, p. 107831, 2025. [8] A. A. Aghaei, “fKAN: Fractional Kolmogorov-Arnold networks with trainable Jacobi basis functions,” arXiv preprint arXiv:2406.07456, 2024. [Online]. Available: https://github.com/alirezaafzalaghaei/fKAN [9] A. Kashefi and T. Mukerji, “Kolmogorov–Arnold PointNet: Deep learning for prediction of fluid fields on irregular geometries,” arXiv preprint arXiv:2504.06327, 2025. [Online]. Available: https://github.com/Ali-Stanford/Physics_Informed_KAN_PointNet [10] J. Xu, Z. Chen, J. Li, S. Yang, W. Wang, X. Hu, and E. C. H. Ngai, “FourierKAN-GCF: Fourier Kolmogorov-Arnold network—An effective and efficient feature transformation for graph collaborative filtering,” arXiv preprint arXiv:2406.01034, 2024. [Online]. Available: https://github.com/Jinfeng-Xu/FKAN-GCF [11] J. Zhang, Y. Fan, K. Cai, and K. Wang, “Kolmogorov-Arnold Fourier networks,” arXiv preprint arXiv:2502.06018, 2025. [Online]. Available: https://github.com/kolmogorovArnoldFourierNetwork/KAF [12] C. C. So and S. P. Yung, “Higher-order ReLU-KANs (HRKANs) for solving physics-informed neural networks (PINNs) more accurately, robustly, and faster,” arXiv preprint arXiv:2409.14248, 2024. [13] Q. Qiu, T. Zhu, H. Gong, L. Chen, and H. Ning, “ReLU-KAN: New Kolmogorov-Arnold networks that only need matrix addition, dot multiplication, and ReLU,” arXiv preprint arXiv:2406.02075, 2024. [Online]. Available: https://github.com/quiqi/relu_kan [14] S. Rigas, M. Papachristou, T. Papadopoulos, F. Anagnostopoulos, and G. Alexandridis, “Adaptive training of grid-dependent physics-informed Kolmogorov-Arnold networks,” IEEE Access, vol. 12, p. 176982–176998, 2024. [Online]. Available: https://github.com/srigas/jaxKAN [15] Z. Bozorgasl and H. Chen, “Wav-KAN: Wavelet Kolmogorov-Arnold networks,” arXiv preprint arXiv:2405.12832, 2024. [Online]. Available: https://github.com/zavareh1/Wav-KAN [16] A. A. Howard, B. Jacob, S. H. Murphy, A. Heinlein, and P. Stinis, “Finite basis Kolmogorov-Arnold networks: Domain decomposition for data-driven and physics-informed problems,” arXiv preprint arXiv:2406.19662, 2024. [Online]. Available: https://github.com/pnnl/neuromancer/tree/feature/fbkans/examples/KANs [17] S. T. Seydi, “Exploring the potential of polynomial basis functions in Kolmogorov-Arnold networks: A comparative study of different groups of polynomials,” arXiv preprint arXiv:2406.02583, 2024. [18] J. A. Actor, G. Harper, B. Southworth, and E. C. Cyr, “Leveraging KANs for expedient training of multichannel MLPs via preconditioning and geometric refinement,” arXiv preprint arXiv:2505.18131, 2025. [19] S. Rigas, D. Verma, G. Alexandridis, and Y. Wang, “Initialization schemes for Kolmogorov-Arnold networks: An empirical study,” arXiv preprint arXiv:2509.03417, 2025. https://github.com/srigas/KAN_Initialization_Schemes [20] A. A. Howard, N. Zolman, B. Jacob, S. L. Brunton, and P. Stinis, “SINDy-KANs: Sparse identification of nonlinear dynamics through Kolmogorov-Arnold networks,” 2026. [21] R. L. Hardy, Multiquadric equations of topography and other irregular surfaces, Journal of Geophysical Research (1896-1977), 76(8), 1905–1915, 1971. [22] E. Kansa, Multiquadrics scattered data approximation scheme with applications to computational fluid-dynamics solutions to parabolic, hyperbolic and elliptic partial differential equations, Computers and Mathematics with Applications, 19(8), 147–161, 1990. [23] Z. Li, “Kolmogorov-Arnold networks are radial basis function networks,” arXiv preprint arXiv:2405.06721, 2024. [Online]. Available: https://github.com/ZiyaoLi/fast-kan [24] B. C. Koenig, S. Kim, and S. Deng, “LeanKAN: A parameter-lean Kolmogorov-Arnold network layer with improved memory efficiency and convergence behavior,” arXiv preprint arXiv:2502.17844, 2025. [Online]. Available: https://github.com/DENG-MIT/LeanKAN [25] D. W. Abueidda, P. Pantidis, and M. E. Mobasher, “DeepOKAN: Deep operator network based on Kolmogorov-Arnold networks for mechanics problems,” Computer Methods in Applied Mechanics and Engineering, vol. 436, p. 117699, 2025. GitHub: https://github.com/DiabAbu/Dee [26] S. T. Chiu, S. W. Cheung, U. Braga-Neto, C. S. Lee, and R. P. Li, “Free-RBF-KAN: Kolmogorov-Arnold Networks with Adaptive Radial Basis Functions for Efficient Function Learning,” arXiv preprint arXiv:2601.07760, 2026. [27] A. Delis, “FasterKAN,” GitHub repository, 2024. [Online]. Available: https://github.com/AthanasiosDelis/faster-kan [28] E. Larsson and R. Schaback, “Scaling of radial basis functions,” IMA Journal of Numerical Analysis, vol. 44, no. 2, p. 1130–1152, 2024. [29] R. Cavoretto and A. De Rossi, “An adaptive refinement scheme for radial basis function collocation,” in Proc. Int. Conf. Numer. Comput.: Theory Algorithms, p. 19–26, 2019. [30] R. Cavoretto, S. Lancellotti, and F. Romaniello, “A Bayesian approach for simultaneously radial kernel parameter tuning in the partition of unity method,” in Proc. Int. Conf. Numer. Comput.: Theory Algorithms, p. 215–222, 2023. [31] T. Wenzel and G. Santin, “On the optimal shape parameter for kernel methods: Sharp direct and inverse statements,” arXiv preprint arXiv:2601.14070, 2026. [32] T. Wenzel and A. Iske, “Spectral alignment of kernel matrices and applications,” SIAM Journal on Matrix Analysis and Applications, vol. 47, no. 1, p. 265–281, 2026. [33] R. Schaback, Error estimates and condition numbers for radial basis function interpolation. Advances in Computational Mathematics, 3(3): 251–264, 1995. [34] R. Schaback, “Small errors imply large evaluation instabilities,” Advances in Computational Mathematics, vol. 49, no. 2, p. 25, 2023. [35] M. Wolff, F. Eilers, and X. Jiang, “CVKAN: Complex-valued Kolmogorov-Arnold networks,” arXiv preprint arXiv:2502.02417, 2025. [Online]. Available: https://github.com/M-Wolff/CVKAN [36] R. Che, L. af Klinteberg, and M. Aryapoor, “Improved Complex-Valued Kolmogorov–Arnold Networks with Theoretical Support,” in Proc. 24th EPIA Conf. on Artificial Intelligence (EPIA), Faro, Portugal, Oct. 2025, Part I, p. 439–451. Springer, Heidelberg. [37] Z. Chao, X. Liu, Z. Wu, and X. Li, “RBF-KAN: Radial Basis Function-Kolmogorov-Arnold Networks,” IEEE Internet Things J., 2026. [38] A. Farea, and M.S. Celebi, Learnable activation functions in physics-informed neural networks for solving partial differential equations. Computer Physics Communications, 2025. 315: p. 109753. [39] Z. Wen, Q. Zhang, J. Chen, et al., “Computing-in-memory architecture for Kolmogorov-Arnold networks based on tunable Gaussian-like memory cells,” Nat. Commun., 2026. [40] G.E. Fasshauer, “Meshfree approximation methods with Matlab” (Vol. 6). World Scientific Publishing Company, 2007. [41] J. Li, A. H. D. Cheng, and C. S. Chen, “A comparison of efficiency and error convergence of multiquadric collocation method and finite element method,” Eng. Anal. Bound. Elem., vol. 27, no. 3, p. 251–257, 2003. [42] A. Noorizadegan, C.-S. Chen, D. L. Young, and C. S. Chen, “Effective condition number for the selection of the RBF shape parameter with the fictitious point method,” Applied Numerical Mathematics, vol. 178, p. 280–295, 2022. [43] C.-S. Chen, A. Noorizadegan, D. L. Young, and C. S. Chen, “On the selection of a better radial basis function and its shape parameter in interpolation problems,” Applied Mathematics and Computation, vol. 442, p. 127713, 2023. [44] A. Noorizadegan, A. Naji, T. L. Lee, R. Cavoretto, and D. L. Young, “Bending analysis of quasicrystal plates using adaptive radial basis function method,” J. Comput. Appl. Math., vol. 450, p. 115990, 2024. [45] A. Noorizadegan, and R. Schaback, “Introducing the evaluation condition number: A novel assessment of conditioning in radial basis function methods. Engineering Analysis with Boundary Elements,” 166, p.105827, 2024. [46] E. Larsson and B. Fornberg, “Theoretical and computational aspects of multivariate interpolation with increasingly flat radial basis functions,” Comput. Math. Appl., vol. 49, no. 1, p. 103–130, 2005. [47] G.E. Fasshauer, and M.J. McCourt, “Stable evaluation of Gaussian radial basis function interpolants,” SIAM Journal on Scientific Computing, 34(2), p.A737-A762, 2012. [48] D. Hilbert, “Mathematical problems,” Bull. Amer. Math. Soc., vol. 8, p. 437–479, 1902. [49] A. N. Kolmogorov, “On the representation of continuous functions of several variables as superpositions of continuous functions of a smaller number of variables,” Dokl. Akad. Nauk SSSR, vol. 108, no. 2, p. 179–182, 1956. (In Russian.) [50] V. I. Arnol’d, “On functions of three variables,” Dokl. Akad. Nauk SSSR, vol. 114, p. 679–681, 1957. (In Russian.) [51] A. N. Kolmogorov, “On the representation of continuous functions of many variables by superposition of continuous functions of one variable and addition,” Doklady Akademii Nauk, vol. 114, p. 953–956, 1957. [52] A. Noorizadegan, S. Wang, , L. Ling, and J.P. Dominguez-Morales, A Practitioner’s Guide to Kolmogorov-Arnold Networks. arXiv preprint arXiv:2510.25781, 2025. [53] T. Yu, J. Qiu, J. Yang, and I. Oseledets, “Sinc Kolmogorov-Arnold network and its applications on physics-informed neural networks,” arXiv preprint arXiv:2410.04096, 2024. [Online]. Available: https://github.com/DUCH714/SincKAN [54] J. Wertz, E. J. Kansa, and L. Ling, “The role of the multiquadric shape parameters in solving elliptic partial differential equations,” Computers & Mathematics with Applications, vol. 51, no. 8, p. 1335–1348, 2006. [55] S. N. Chiu, L. Ling, and M. McCourt, “On variable and random shape Gaussian interpolations,” Applied Mathematics and Computation, vol. 377, p. 125159, 2020. [56] G. E. Fasshauer and J. G. Zhang, “On choosing “optimal” shape parameters for RBF approximation,” Numer. Algorithms, vol. 45, no. 1, p. 345–368, 2007. [57] H. Wendland, Scattered Data Approximation, Cambridge, U.K.: Cambridge Univ. Press, 2005. [58] M. S. Floater, and A. Iske, “ Multistep scattered data interpolation using compactly supported radial basis functions,” Journal of Computational and Applied Mathematics, 73(1-2), 65-78, 1996. [59] S. Müller and R. Schaback, “A Newton basis for kernel spaces,” J. Approx. Theory, vol. 161, no. 2, p. 645–655, 2009.

Instant research discovery