Kolmogorov Arnold networks (KAN) for aerodynamic prediction: a comparison with MLPs and GNNs Miguel Jaraiz1 , Fermin Gutierrez1 , Pablo Yeste1 , Miguel Sánchez-Domı́nguez1 , Eusebio Valero1,3 , Gonzalo Rubio1,3 , and Lucas Lacasa2
arXiv:2606.27126v1 [cs.LG] 25 Jun 2026
1
ETSIAE-UPM-School of Aeronautics, Universidad Politécnica de Madrid, Plaza Cardenal Cisneros 3, E-28040 Madrid, Spain 2 Institute for Cross-Disciplinary Physics and Complex Systems (IFISC, CSIC-UIB), 07122 Palma de Mallorca (Spain) 3 Center for Computational Simulation, Universidad Politécnica de Madrid, Campus de Montegancedo, Boadilla del Monte, 28660 Madrid, Spain June 26, 2026 Abstract Kolmogorov Arnold networks (KAN) have recently been introduced as a (deep) neural network architecture whose trainable parameters adapt the activation functions, instead of the coefficients of the affine transformations at the core of traditional architectures such as deep multilayer perceptrons (MLPs). This architecture builds on the Kolmogorov-Arnold theorem, which endows it with universal approximation properties. While the advent of KANs has been received with excitement, there is a current debate about the possible KAN supremacy over deep multilayer perceptrons (MLPs) for classic fields such as symbolic regression, generic-purpose machine learning, natural language processing or computer vision. Here we assess the performance of KANs –and its nuanced comparison against MLPs and graph neural networks (GNNs)– in the realm of fluid dynamics surrogate modelling. To that aim, we consider the task of predicting the surface pressure distribution over subsonic and transonic airfoils, a canonical task in aerodynamics. Our results show that KAN models show good performance in predicting the whole pressure coefficients and is able to interpolate across Mach numbers and angles of attack, however its performance is comparable –marginally inferior– to a suitably trained MLP, where best performance is achieved by a GNN at the expense or requiring lengthier training. While the optimal KAN model have typically much lower complexity than MLP and GNN –hence resulting in faster training–, we find that KANs suffer from training instabilities, and their performance is highly dependent on a proper hyperparameter optimisation.
1
Introduction
The simulation of physical processes across a wide range of disciplines boils down to numerically solving partial differential equations (PDE) over a specific domain geometry and a range of physical parameters. A paradigmatic case is Computational Fluid Dynamics (CFD) [1] –a fundamental simulation tool e.g. for aircraft design– which proceeds to integrate flow equations of motion (i.e. Navier-Stokes equations or adequate approximations such as Reynolds-Averaged Navier-Stokes (RANS) or Large Eddy Simulations (LES)). CFD is often seen as the gold standard for the numerical analysis of fluid-dynamical problems as it provides high-fidelity solutions of complex flow phenomena. However, CFD-based methods are computationally expensive, both in terms of simulation time and required computational resources. For example, when estimating aerodynamic loads, each unique flight condition within an aircraft’s operational envelope requires a separate CFD simulation, leading to an immense computational burden. Surrogate models [2] provide an efficient alternative, allowing for rapid approximation of CFD results across varying conditions without the need for costly repeated simulations. Traditional surrogate modeling has been approached through data-driven methods such as proper orthogonal decomposition (POD) 1
combined with interpolation techniques [3–5] or dynamic mode decomposition (DMD) [2]. These models are relatively straightforward to construct and perform well for predictable linear behavior, but struggle to accurately approximate intrinsically nonlinear behavior. The advent of powerful machine learning methods, particularly deep supervided learning (DL) [6] has opened in recent years a new approach to construct surrogate models [2]. In the context of fluid-dynamical surrogate models, several studies have demonstrated the potential of machine learning for reduced-order modeling [7–18], including not only academically-oriented models, but also realistic, industrial-level 3D geometries, where unstructured grids with millions of points are necessary to accurately capture aircraft configurations. A paradigmatic problem in aerodynamic CFD with crucial industrial interest is the prediction of aerodynamic coefficients over the surface of an aircraft, e.g. the prediction of pressure coefficients over an airfoil for specific physical parameters (Reynolds or Mach number) and flight conditions (e.g. angle of attack). Comparative studies highlight various techniques in surrogate regression models for aerodynamic coefficient prediction [19]. For instance, in [20] a coordinate-based MLP –the default DL architecture for supervised learning regression tasks– was proposed for the pointwise prediction of such pressure coefficients. Graph neural networks (GNN) have also been proposed as more sophisticated supervised learning architectures for this task [21], and more recently generative AI (flow matching) has also been repurposed for this task [22]. Within supervised learning, observe that both MLPs and GNNs rely on universal approximation theorems [23] to approximate nonlinear functions in terms of an sequential composition of affine combinations and nonlinear activation functions. One of the latest advancements in this field is the inception of so-called Kolmogorov–Arnold Networks (KANs) [24]. At odds with previous architectures, KANs rely on a different representation theorem (Kolmogorov-Arnold representation theorem [25, 26]) and were introduced as an alternative to MLPs where, unlike MLPs, KANs feature learnable activation functions on edges, promising greater flexibility and interpretability. Since their emergence, KANs have received wide popularity and shown great promise in a number of applications, both in science and industry [27]. The supremacy of KANs over MLPs is, however, currently a matter of debate [28]. For instance, recent works [29] argue that MLP might actually remain superior than KANs over a large set of tasks including machine learning, computer vision, audio processing, or natural language processing. In this context, it is unclear whether KANs generally shows higher performance –as compared to MLPs or related architectures– for the design of surrogate models of physical processes [30]. In this paper, we build a KAN-based surrogate model for a vanilla aerodynamic task: the prediction of pressure coefficients over a boundary airfoil, and compare its performance in interpolating across physical parameters and flight conditions (Mach number, angle of attack) against two well-established techniques: MLPs [20] and GNNs [21]. We build on [21], replicate their results on the performance of both MLPs and GNNs on the same learning task and then compare these against the KAN-based surrogate model. Our results indicate that while all three models achieve good results (where GNN reaches the best performance), the KAN-based model ranks third, achieving comparable yet marginally inferior performance with respect to its MLP-based counterpart, and suffers from training instabilities and poorer generalisation error in the more difficult regions where the pressure coefficients show stronger gradients. At the same time, training of KAN model converges substantially faster –i.e. needing fewer training epochs– than the MLP and GNN, arguably because the complexity of optimal MLP and GNN architectures is substantially larger (in terms of number of trainable parameters) than the optimal KAN architecture. This effect has recently been observed in a KAN vs MLP comparison within reinforcement learning [31]. Now, by training smaller MLPs –with a complexity comparable or even smaller than the KAN model–, we show that we can still reach generalisation performance on par with KANs. Despite the recent popularity of KANs, these results overall challenge their acclaimed supremacy over more traditional deep learning architectures. The rest of the paper is organised as follows. In Section 2, we introduce the entire methodology: the learning task, the data, and details on all three surrogate models. Section 3 reports our results, including a comprehensive hyperparameter optimisation analysis –highlighting the strong dependency of KAN’s performance on a proper choice of hyperparameters, even on the random seed [32]–, the comparison of the performance for optimal MLPs, GNNs and KANs in terms of learning curve convergence, pointwise test error metrics and preservation of correlation structures. In this section we also discuss how MLP and KAN of comparable complexities compare in terms of performance. Finally, in Section 4 we conclude.
2
2
Methodology
2.1
Problem setup and data preparation
The learning task considered here is identical across all data-driven models: given a vector of operating conditions, predict the distribution of pressure coefficients over an airfoil via a RANS-based CFD simulation. To that aim, the airfoil geometry is first placed inside a computational domain and discretised with an unstructured mesh that resolves the boundary layer near the surface and extends to a far-field boundary representing the freestream conditions. The RANS equations (together with a turbulence model) are then discretised on this mesh using a finite-volume method and solved iteratively until the velocity and pressure fields converge. Once the solution is obtained, the static pressure values at the nodes/faces that lie on the airfoil boundary are extracted from the computed flow field. These pressures are then converted into pressure coefficients by subtracting the freestream pressure and normalizing by the freestream dynamic pressure. Mapping these values along the airfoil surface as a function of the chordwise coordinate produces the pressure-coefficient distribution Cp (x) on the airfoil.
Figure 1: Data split. The database available for this study was kindly shared by DLR [21], but we nonetheless provide details on how this was constructed for the sake of self-containedness. As the surface geometry, we use the supercritical NLR7301 airfoil introduced in Hines and Bekemeyer [21]. High-fidelity RANS-based CFD computations were performed with the DLR flow solver TAU [33], employing the Spalart-Allmaras turbulence model. The computational domain is discretised with an unstructured mesh of 29, 441 nodes, where a total of n = 597 nodes are located on the airfoil boundary. The Reynolds number is fixed in all cases at Re = 1.7 × 106 and two independent operating parameters are varied across simulations: the Mach number M ∈ [0.3, 0.75] and the angle of attack α ∈ [−3◦ , 5◦ ]. The dataset S [21] comprises a total of m = 98 different CFD simulations of the RANS equations obtained by varying (M, α) so as to adequately cover the M − α plane. Accordingly, S = {(zk , yk )}m k=1 , where the input vector of the k-th datapoint zk = (M, α)k determines the operating condition, and its output vector yk = [Cp (x1 ), Cp (x2 ), . . . , Cp (xn )]k determines the pressure coefficients Cp at each of the i = 1, 2, . . . , n airfoil’s boundary nodes. Following [21], S is then adequately split into training, validation and test sets S = Straining ∪ Svalidation ∪ Stest , with mtraining = 39, mvalidation = 20, and mtest = 39, see Figure 1 for an illustration.
3
2.2
Surrogate models
Three different surrogate models –a MLP, a GNN and a KAN, see below– are all initially trained on Straining . Both the MLP and KAN-based surrogate models are here considered as local predictors, i.e. at inference both these models are given a pair (M, α) and a specific node location xi as inputs, and these models yield the prediction of the pressure coefficient at that specific location MLP KAN Cˆp = f MLP (M, α, xi ) ∈ R; Cˆp = f KAN (M, α, xi ) ∈ R, ∀i = 1, 2, . . . , n.
In other words, both MLP and KAN will take a triplet (z, x) as the input vector. On the other hand, the GNN-based surrogate model is a global predictor, such that, at inference, this model only requires specification of the pair (M, α), and the model predicts the pressure coefficient at all the n spatial locations at once: GNN Cˆp = f GNN (M, α);
GNN Cˆp = [Cˆp (x1 ), Cˆp (x2 ), . . . , Cˆp (xn )] ∈ Rn .
For convenience, in what follows we will use z ≡ (M, α) or z ≡ (M, α, x) depending on the context. Below we provide a reminder on the basic properties of the MLP, GNN and KAN architectures. In order to compare our baseline results on MLP and GNN with previous works, our selection of the MLP and GNN specification matches the one in [21]. MLP – Multi-Layer Perceptrons (MLPs) are the foundation of many deep supervised learning models, thanks to their expressivity as universal approximators. A multilayer perceptron (MLP) with L hidden layers defines a nonlinear mapping from an input vector (in our case z = (M, α, x) ∈ R3 , where x defines a spatial node location in the airfoil boundary) to an output vector (in our case ŷ = Cˆp (x) ∈ R) through a sequence of affine transformations followed by element-wise nonlinear activation functions. Formally, denoting the hidden-layer activations by h(l) , the MLP is a nonlinear function that can be written recursively as h(l) = σ W(l) h(l−1) + b(l) for l = 1, . . . , L, with h(0) = z, where W(l) and b(l) are the weight matrices and bias vectors, and σ(·) is a scalar nonlinear activation function, which in our case is the exponential linear unit ( w, w > 0, ELU(w) = ew − 1, w ≤ 0. The optimal number of layers L is an hyperparameter found after an hyperparameter optimisation, as is explained in the next subsection. The optimal assignment of network parameters θ = {W(l) , b(l) }L+1 l=1 are learned by minimizing a loss function L(θ) over the training dataset Straining via an Adam optimiser. In our case, we use a mean-square error (MSE) type, where L(θ) =
1 mtraining
mtraining
X
k=1
n 1X ˆ 2 [Cp (xi , zk ) − Cp (xi , zk )] . n i=1
(1)
where Cˆp (xi , zk ), Cp (xi , zk ) are respectively the predicted pressure coefficient for operating conditions zk evaluated at the spatial node xi , and its ground true value (the one provided by the RANS simulation) at each location of the airfoil boundary. GNN – A Graph Neural Network (GNN) [34, 35] extends the multilayer perceptron (MLP) paradigm to input data that are naturally defined on (or benefit from) graphs, instead of vectors. Accordingly, while an MLP applies a sequence of dense affine transformations to a fixed input vector, a GNN operates (0) on a graph G = (V, E) with node features hi , and updates each node representation by aggregating information from its neighbors. At layer l, the representation of node i is typically updated as (l) (l−1) hi = σ W (l) AGG {hj : j ∈ N (i)} , where N (i) denotes the neighbors of node i, AGG is a permutation-invariant aggregation operator (e.g., sum, mean, or max), W (l) is a learnable weight matrix, and σ is a nonlinear activation function. In contrast to an MLP, which treats inputs as independent features, the GNN explicitly propagates 4
information along graph edges, allowing the learned representation of each node to depend on the features of its local neighborhood. When applied to an airfoil, a Graph Neural Network (GNN) operates on a graph whose nodes correspond to discrete points along the airfoil surface, typically ordered along the contour so that the graph forms a closed chain. Each node carries local features such as its coordinates or other geometric descriptors of the airfoil. Edges connect neighboring surface nodes, reflecting the physical adjacency of points along the airfoil. During each GNN layer, the feature vector at a node is updated by aggregating information from its neighboring nodes and combining it with its own features through a learnable transformation and nonlinearity. This structure is particularly suitable for modeling aerodynamic quantities such as pressure coefficients, which are not independent at each surface point: due to the continuity and smoothness of the flow field along the surface, the pressure at a given location is strongly correlated with the pressure at nearby locations. By propagating information along the graph, the GNN can naturally learn and exploit these local correlations while still capturing the influence of the global airfoil geometry. The loss function matches Eq. 1, and training is also done via Adam optimisation. Unlike MLP, for GNN some additional geometric information is needed as input features such as surface normals and face normals, see [21] for details. KAN – Both MLPs and GNNs rely on the their expressivity as universal approximators, as per the universal approximation theorem [23]. In turn, the Kolmogorov-Arnold representation theorem [25] states that any continuous multivariate function can be represented as a finite sum of compositions of univariate functions. Of course, the challenge lies in finding the precise shape of these univariate functions. Very recently, Liu and collaborators introduced Kolmogorov–Arnold Networks (KANs) [24, 36] along with a practical way to learn those univariate functions from data (albeit chosen within a finite basis of candidate functions –initially a spline basis–, hence representation is not, in general, exact). Accordingly, KANs build a nonlinear mapping by replacing the (learnable) scalar weights of a standard MLP with learnable (parametric) univariate functions on the edges of the neural network. For our input vector z ∈ R3 , each layer computes a sum of transformed inputs, where the transformations are parameterised one–dimensional functions. Denoting the activations at layer l by h(l) , a KAN layer can be written as X (l) (l−1) (l) hj = ϕij hi , l = 1, . . . , L, i (l)
with h(0) = z = (M, α, x), where ϕij (·) are learnable univariate functions (represented using spline bases in the seminal work [36], but here we use Chebyshev polynomials as the orthonormal basis). The network output is obtained from the final layer in the same additive form. The parameters defining these functions are learned again by minimizing a MSE-type loss similar to the case of the MLP over Straining . Each of the three models have an additional set of parameters which are not trained via loss minimization (e.g. the number of layers, some parameters of the optimiser, etc). For reproducibility, in the case of the MLP and the GNN surrogates, these hyperparameters are initially chosen following the hyperparameter optimisation in [21], and are depicted in Table 1 and Table 2 for the MLP and the GNN, respectively. The hyperparameter optimisation of the KAN model is detailed in Sec. 3.1. All three models are trained on the training set Straining . For hyperparameter optimisation, we train models with different configurations of hyperparameters in Straining and select as the best configuration the one for which the model shows smaller MSE in the validation set Svalidation . The generalisation performance of the best models is then assessed in the test set Stest , see below for details. Hyperparameter
Value
Initial learning rate LR decay factor Batch size Dimension of hidden layers Number of hidden layers
3.45 × 10−4 0.9954 415 256 11
Trainable parameters
725249
Table 1: Hyperparameters used for the MLP model, taken from [21].
5
Hyperparameter
Value
Initial learning rate LR decay factor Parameter batch size Node batch size Dimension of hidden layers Latent dimension Encoder hidden layers GNN layers M LPMt hidden layers M LPUt hidden layers Decoder hidden layers
6.50 × 10−4 0.9954 15 32 256 16 6 1 2 2 1
Trainable parameters
299285
Table 2: Hyperparameters used for the GNN model, taken from [21].
Error quantification – Once the three models are properly trained, their generalisation peformance is initially tested on Stest using a variety of standard global quantifiers [6, 37]. Instead of a MSE-type error (Eq. 1), we use the mean absolute error (MAE) and the root mean squared error (RMSE), where the latter one is essentially the square root of Eq. 1 (computed on the test set), whereas the former is like Eq. 1 but using the absolute value of the residuals, instead of the square value. Additionally, from the population of all absolute values of the residuals {|Cˆp (xi , zk ) − Cp (xi , zk )|} we extract their 95 and 99 percentiles, as indicators of extreme errors. Finally, we also compute the coefficient of determination R2 , computed from a scatter plot between the ground true pressure coefficient and the prediction, for all spatial locations and operating conditions in the test set. In addition, observe that spatial correlations are expected to emerge among the set of surface pressure coefficients, due to physical constraints and simple continuity arguments. Accordingly, those correlations should also be captured by the surrogate model, even if it is a local regressor. To that aim, for each model we construct two 597 × 597 correlation matrices Ctrue and Cpred . The ij-entry of these matrices denote, respectively, the (properly normalised) correlation between the pressure coefficient at spatial location xi and the one at spatial location xj (computed over all the operating conditions zk in the test set) for the ground truth and the model predictions, i.e. m
Ctrue = ij
test 1 X [Cp (xi , zk ) − µi ][Cp (xj , zk ) − µj ] , mtest σi σj
(2)
k=1
where µX and σX are the average and standard deviation of all the pressure coefficients at spatial location X (over the ensemble of operating conditions considered in Stest ). Finally, to quantify the discrepancy between Cpred and Ctrue for each model we use the Frobenius relative error [38]: εF =
∥Cpred − Ctrue ∥F . ∥Ctrue ∥F
(3)
Lower values of εF indicate a closer match to the true correlation structure.
3
Results
3.1
Hyperparameter optimisation and impact of RNG seed: KAN vs MLP
At odds with the MLP and the GNN architectures [21], KAN-based surrogate models have not, to the best of our knowledge, been tested before in the context of this learning task. Accordingly, for this model we need to perform an adequate hyperparameter optimisation in Svalidation . In addition, this section also explores the possible dependence on the performance of both KAN and MLP models with respect to the seed of the pseudo random number generator (RNG), which impacts (i) the weight initialization of 6
the networks, and (ii) the selection of mini-batches. Such potential dependence is seldom explored in the literature and, as we shall see, its impact for the KAN performance is (unexpectedly) substantially larger than for the MLP.
3.1.1
Hyperparameter optimisation for the KAN model
The set of hyperparameters include the number of layers and number of neurons per layer, the learning rate decay factor in the Adam optimiser, or the degree of the Chebyshev polynomials in the trainable activation functions, among others (see Table 3 for details). After selecting the initial range for each of the hyperparameters, we perform two complementary analysis. In the first one, we sample at random (quasi-Montecarlo) a total of 2472 different hyperparameter configurations. For each configuration, we train the KAN-based model on the training set Straining , and compute its MSE-based loss on the validation set Svalidation . Among the initial 2472 different hyperparameter configurations, the KAN-based model only converged in 312 cases, and the rest were discarded mostly due to exploding gradients. This feature highlight the intrinsic instability of the KAN architecture with respect to hyperparameters [28]. Results for the remaining 312 trials is reported in Figure 2. Observe that the distribution of MSE losses is long-tailed (note that the frequency histogram is in linear-log scales), highlighting that whereas a large number of configurations have small validation MSE loss, a non-negligible amount of configurations substantially deviate and can attain much larger errors, i.e. the dependency of the KAN-based model’s performance on the selection of its hyperparameters is strong, with many poorly-performing outliers emerging if the hyperparameters are not carefully chosen. We then investigate whether tendencies emerge between increasing/decreasing values of the hyperparameters and performance in the validation set, by computing Spearman correlation coefficients (and their statistical significance), see Figure 3b. The correlation between selected pairs of hyperparameters is generally low, although some for some pairs such correlation appear to be statistically significant (gray entries in the p-value matrix). The last row of these matrices indicate the dependence on the loss of different hyperparameters, showing moderate yet significant tendency for the learning rate (larger learning rates lead to larger loss error) and the number of layers (a larger number moderately correlate with a smaller error). These dependencies are overall weak and non-informative, see for instance table 3 where the optimal number of layers is found to be small. Hyperparameter
Minimum value
Maximum value
Optimal value
Number of hidden layers Neurons per layer Batch size (2p ) Learning rate decay factor (γ) Degree of Chebyshev polynomials Dropout probability Number of epochs Initial learning rate
1 10 p=4 0.9 4 10−5 100 10−5
6 400 p = 10 0.999 8 10−1 800 10−2
2 25 p=7 9.5071 × 10−1 4 1.7878 × 10−5 474 1.3024 × 10−3
Trainable parameters
6875
Table 3: Range of the different hyperparameters and their optimal configuration obtained for the KAN model via Optuna.
The second analysis repeats this whole procedure albeit performing a non-uniform sampling of the hyperparameters with the Optuna Python library [39], concretely using a Tree-Structured Parzen Estimator (TPE). A total of 311 converged trials (hyperparameter configurations) are found, and the frequency histogram of resulting MSE losses in the validation set for this hyperparameter sampling is reported in Figure 4. Compared to Figure 2, we observe that this non-uniform sampling concentrates the measure in the region of small MSE losses, although the long tail is not completely removed. We use this sampling to select the hyperparameter configuration that minimises the MSE on the validation set. Interestintly, the optimal configuration (the left end of the histogram) itself is not an outlier among the distribution
7
of configurations for which the model is properly trained. This suggests that the performance of this optimal KAN model is reasonably robust against small changes in hyperparameters. The optimal configuration is depicted in Table 3.
Figure 2: Performance of trained KAN-based models (in terms of MSE loss in the validation set) for a quasi-Montecarlo sampling of hyperparameters. (Left) The left panel reports the MSE loss for the different trials of hyperparameters (the black dots depict how the smallest detected MSE changes as trials progress). Observe that we initially check 2472 trials, but only 312 produced converged resuls (blue dots). (Right) Linear-log frequency histogram of MSE loss (in the validation set) for the 312 sets of hyperparmaters. Observe that the histogram has a long right tail, suggesting that the performance of the KAN-based model is highly dependent on its hyperparameters.
(a) Spearman matrix.
(b) p-values matrix.
Figure 3: Spearman correlation analysis between all pairs of hyperparameters (and the validation loss) computed over the 312 hyperparameter configurations (those surviving from the initial 2472 trials).
3.1.2
Impact of the RNG seed
As an additional analysis, we explore –both for the KAN and MLP-based models– the impact that different RNG seeds might have on the model’s performance in the validation set. Observe that changing the seed of the pseudo-random number generator changes (i) the model’s initialisation, and (ii) the order in which the data batches enter into the optimisation loop within each epoch. Different seeds are sampled uniformly at random and without replacement (a total of 468 seeds are considered for the KAN-based model, and 280 seeds for the MLP). Results are shown in panels (a) and (b) of Figure 5. While the MLP model is relatively robust against changing the RNG seed as expected [32], we find that in the 8
Figure 4: Same as Figure 2 but applied to a selection of hyperparameters from Optuna [39], that uses a Tree-Structured Parzen Estimator (TPE). The tail of large MSE loss in the frequency histogram (right panel) has less mass, but is not completely removed. KAN model there are some RNG seeds that induce the model to perform poorly at validation set. This reinforces the evidence presented in Figs. 2 and 4 about the strong dependence of KAN-based model’s performance on hyperparameters.
(a) KAN (468 different seed trials).
(b) MLP (334 different seed trials).
Figure 5: Linear-log frequency histograms of MSE loss (in the validation set) for the KAN-based model (panel a) and the MLP model (panel b) obtained when uniformly samply different RNG seeds. The total number of different sampled seeds is 468 for the KAN model and 334 for the MLP model. The impact of choosing different RNG seeds has a substantially larger impact for the KAN-based model, where we find seed instances that lead to outlier results in terms of very large MSE loss, whereas in the MLP model performance is relatively robust against seed variation.
3.2
Learning curves
Once the hyperparameters of the models have been chosen, a final training run is carried out using the best configurations. The corresponding learning curves (semi-log plots of MSE-based loss as a function of the training progress, calculated both in the training –solid lines– and validation sets –dashed lines– ) for the three models are presented in Figure 6. Both the MLP and GNN models are trained over 1000 epochs, whereas the KAN model requires smaller number of epochs (hyperparameter optimisation yields roughly below 500 epochs). On the training set, both MLP and GNN exhibit comparable learning behavior (with the GNN achieving lower training losses), and continue to decrease the training loss even after 1000 epochs. Comparing these against the curves on the validation set, which converge after around 500 epochs, suggest that for longer epochs both MLP and GNN are overfitting to the training set, with no impact on validation. MLP and GNN attain comparable loss performance on the validation set, stabilising at a MSE loss around 10−3 . The KAN model in turn seems to learn much faster, with both training and validation curves stabilising after about 50 epochs –i.e. only 10% of the number of epochs needed for both MLP and GNN–, achieving a MSE loss in the validation set of 8 × 10−4 , i.e. slightly lower than the MLP and GNN models.
9
Figure 6: Training and validation MSE loss curves as a function of the number of epochs for the MLP, KAN, and GNN models after hyperparameter optimization. DLR [21] MSE MAE RMSE AEq95 AEq99 R2 εF
this work
MLP
GNN
MLP
GNN
KAN
0.0068 0.0335 0.0225 0.1164 0.9954
0.0061 0.2876(∗) 0.0202 0.1176 0.9966
9.7 × 10−4 0.0085 0.0311 0.0305 0.1359 0.9960 0.1158
8.4 × 10−4 0.0065 0.0290 0.0204 0.1195 0.9965 0.1142
1.2 × 10−3 0.0137 0.0346 0.0566 0.1667 0.9951 0.1457
Table 4: Pointwise error metrics on the test set for all three models. Results are grouped by source and model type. We believe [21] has a typo in (*), and possibly the RMSE of the GNN is 0.02876 instead of 0.2876.
3.3
Generalisation performance across models: global metrics
The global performance of the three models on the test set is summarised in Table 4, which also includes the reference values reported by Hines and Bekemeyer [21] for MLP and GNN. This comparison shows that the implementations developed here for both MLP and GNN yield metrics that are broadly consistent with those of the reference study –with GNN showing superior performance over MLP–. The performance of all three surrogate models is arguably very good (e.g. reaching R2 ≥ 0.995 in the three cases). This is a commendable result, reminding that the data split is performed only with respect to the operating conditions (M, α): in this sense we are assessing the ability of the surrogate models to interpolate –and in some cases, also extrapolate– over operating conditions, rather than over spatial locations on the airfoil’s surface. Additionally, observe that the MSE metric on the test set is comparable (only slightly smaller) than the one on the validation set (Figure 6), overall certifying that the training and validation processes have been done correctly and that the data split in Figure 1 is adequate. Now, when comparing the KAN model with the established MLP and GNN baselines, a small performance gap appears. The KAN exhibits slightly higher values across all error metrics, indicating a marginally inferior predictive accuracy. Recall, however, that the KAN model requires only about 10% of the number of training epochs required by both MLP and GNN, hence KAN’s comparatively inferior performance needs to be traded-off with its apparent improvement on training efficiency. We shall come back to this
10
point later.
3.4
Generalisation performance across models: additional metrics
In this section we provide a more detailed analysis of the generalisation performance of the three models. In the three panels of Figure 7, we plot scatter plots of the true pressure coefficient Cp (xi , zk ) vs its prediction Ĉp (xi , zk ), for all location points xi and operating conditions zk = (M, α)k in the test set (the coefficients of determination R2 are depicted in Table 4). We color-code each point (Cp (xi , zk ), Ĉp (xi , zk )) in terms of its residue value Cp (xi , zk ) − Ĉp (xi , zk ). For instance, points above the diagonal correspond to cases where Ĉp (xi , zk ) > Cp (xi , zk ). This kind of color-coded scatter-plots allows us to assess in what regions of the output space the model is more or less conservative, something important in safetycritical applications [40]. Comparing across models (Panels a-c), the MLP and GNN models exhibit their largest deviations in the central region of the pressure coefficient distribution, particularly for values in the range of approximately −2 to 0. These values are typically associated with the suction side (extrados) of the airfoil, where the flow acceleration leads to negative pressure coefficients. The presence of larger deviations in this range suggests that part of the underlying physics governing the flow in this region is not being adequately captured by the MLP and GNN models. In turn, the KAN shows a slightly broader dispersion across the entire range of Cp , indicating difficulties in accurately capturing the pressure coefficient over all regions of the domain. This behavior is consistent with the global performance metrics, which are slightly inferior compared to those of the MLP and GNN models..
(a) MLP.
(b) GNN.
(c) KAN.
Figure 7: True versus predicted surface pressure coefficients for the three models. The color scale indicates the residual magnitude.
To further understand where prediction errors are located within the airfoil, the pressure coefficient distributions for four operating conditions are presented in Figure 8, following the selection criteria of Hines and Bekemeyer [21]. The first two correspond to cases in the low subsonic regime, whereas the last two represent transonic conditions. These cases are chosen such that, for each regime, one case lies inside the convex hull of the training set (in the M − α plane, hence interpolating operating conditions), and the other lies outside the convex hull, i.e. effectively extrapolating operating conditions. In the subsonic interpolation case (Figure 8a), all three models are able to reproduce the true pressure coefficient distribution with good accuracy, although the KAN slightly underestimates the suction peak near the leading edge. In contrast, for the subsonic extrapolation case (Figure 8b), the KAN exhibits more pronounced deviations throughout the distribution compared to the MLP and GNN, indicating a slightly weaker extrapolation capability. For the transonic interpolation case in Figure 8c, the three models successfully capture the shock that appears on the suction side of the airfoil, but the KAN overestimates the suction peak and exhibits a less accurate reconstruction of the pressure plateau downstream of the shock. Finally, the transonic extrapolation in Figure 8d represents the most challenging scenario: none of the models accurately predicts the precise shape of the shock in the airfoil. The MLP and GNN tend to shift it downstream, but they represent the pressure jump with a more realistic gradient. In contrast, although the KAN predicts the geometric location on the suction side where the shock begins more accurately, it tends to excessively smooth the transition, leading to a more diffuse pressure gradient.
11
(a) Subsonic interpolation case: M = 0.4441, α = 3.123 deg.
(b) Subsonic extrapolation: M = 0.3176, α = 4.901 deg.
(c) Transonic interpolation: M = 0.6304, α = 1.642 deg.
(d) Transonic extrapolation: M = 0.6972, α = 4.605 deg.
Figure 8: Surface pressure distribution at four different operating conditions of the test set. Each panel additionally shows an inset with a zoom of a region of interest.
To round-off the analysis, Figure 9 compares the test-set correlation matrices (Eq. 2) for the ground truth and for the three surrogate models, whereas its Frobenius relative error εF (Eq. 3) is depicted in the last line of Table 4. All three models preserve the correlation structure, capturing the dominant structures and the transition zones. In comparison, the GNN model is the one that, overall (according to εF ) better preserves the correlation structure among pressure coefficients over the airfoil. This is somewhat expected since, being a global predictor, this model explicitly uses such spatial information during training. In this sense, the performance of both the MLP and the KAN is commendable, given that these are trained as local regressors where no inductive bias related to the spatial continuity of Cp is explicitly incorporated.
12
Figure 9: Correlation matrices for the test set: (top left) reference, (top right) MLP, (bottom left) KAN, and (bottom right) GNN.
3.5
MLP vs KAN with comparable complexity
In previous sections we have found that, while both MLP and KAN-based surrogate models achieve good performance on the learning task studied here, the performance of MLP is slightly yet systematically superior over the KAN for a variety of global and local metrics. At the same time, we also found that the optimal KAN model was much smaller than the optimal MLP one given in [21] (in terms of number of trainable parameters), and its training was accordingly significantly faster. Thus, it remains difficult to conclude the supremacy of one against the other, since a direct comparison of the two optimised models is not really fair as their model complexities are very different (in particular, the number of hidden layers of the optimal MLP is much larger than for the KAN model). To be able to make a fair comparison, here we depart from the optimal MLP configuration given in [21] and compare the performance of MLP vs KAN for smaller MLP models (in terms of number of hidden layers and/or number of neurons per layer), so that the overall complexity of both models (in terms of total number of trainable parameters) is made comparable. Results are shown in Figure 10, and certify that when we reduce the complexity of the MLP model –and thus reducing its training requirements–, its performance (measured in the validation set) is similar to the KAN model although typically a bit worse. Observe that the best MLP are the ones that preserve a large layer depth (11 hidden layers). In other words, it seems that the expressivity of the MLP comes from its layer depth, something which is in contrast to the KAN model, that based on the hyperparameter optimisation operates with only two hidden layers. This analysis indicates that while MLP might seem more suitable than KAN if we only consider performance, if we are restricted to low-complexity models (i.e. due to restrictions in computational resources), a KAN model might be preferred, provided that we carefully optimise its hyperparameters.
4
Conclusions
All three deep learning surrogate models constructed in this work –a multilayer perceptron, a graph neural network and a Kolmogorov-Arnold network– showed very good performance in the prediction of pressure coefficients alongside the boundary of an airfoil, and were able to accurately interpolate across
13
Figure 10: Comparison of the KAN model against a MLP of comparable complexity. (Left panel): Heatmap of the validation set MSE of a MLP model as a function of its model’s complexity (number of hidden layers and layer width). The cases exceeding the performance level of the KAN are highlighted in pink, with this reference level indicated in the color bar by a red line. (Right panel): Validation set MSE of a MLP model, as a function of its model complexity (total number of trainable parameters). For reference, the complexity of the optimal KAN model and its validation MSE is highlighted (dashed lines). different Mach numbers and angles of attack. Among the three models, typically the GNN showed the best results, followed by the MLP. The KAN model was marginally inferior to the MLP, hence providing relevant insights on the debate around the KAN supremacy over MLP [29]. Incidentally, the training of KAN model proved to be significantly more challenging than that of the other architectures considered in this study. In particular, the optimization process exhibited notable instabilities during training, frequently leading to extremely large gradients. This behavior complicated the learning process and made the optimization of hyperparameters substantially more challenging. As a consequence, achieving stable convergence requires careful tuning and increases the computational effort associated with training these models. Regarding predictive performance, as mentioned previously all three models showed notably low generalisation (test) error, but KANs exhibited slightly inferior results. Regions with strong pressure gradients along the airfoil surface were more difficult to accurately predict for all cases. One positive property of KANs over both MLPs and GNNs is that its optimal model complexity (in terms of number of training parameters) is orders of magnitude smaller than the MLP and GNN ones, and thus the number of epochs needed to achieve good performance is systematically smaller in the KAN case –what results in faster training– when comparing optimal model configurations. When comparing MLP and KAN models of similar number of trainable parameters, the KAN architecture showed to be competitive and typically at least as good as the MLP. Taken together, these results provide a nuanced view on the MLP vs KAN debate in the context of building fluid-dynamics-based surrogate models, with no clear supremacy. Classical deep learning architectures such as MLP or GNN show marginally superior or marginally inferior performance as compared to KAN when considering different aspects such as training stability and efficiency, overall performance, or performance restricted to low-complexity models. Further work should clarify whether these insights hold beyond aerodynamic cases within fluid-dynamics-based surrogate modeling. Code – All stages of the pipeline followed throughout this work, including training, hyperparameter optimization, and evaluation of the three surrogate models has been executed on the FLEXO cluster at the Universidad Politécnica de Madrid (UPM), which includes a computing node equipped with 4 AMD Instinct MI210 GPUs and 80 Intel(R) Xeon(R) Silver 4416+ CPU cores. A single GPU and 10 CPU cores were used for each model, allowing parallel data loading and preprocessing. Code details are available at https://github.com/ArnauMiro/pyLowOrder.git [41], where the implementation of the model architectures is provided. Acknowledgments — We thank Derrick Hines and Philipp Bekemeyer (DLR) for sharing their database 14
[21] and for insightful discussions. The authors acknowledge funding from project TIFON (PLEC2023010251) funded by MCIN/AEI/10.13039/501100011033, Spain. The authors acknowledge funding from the European Union (project HERFUSE) under GA No 101140567. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or Clean Aviation Joint Undertaking. Neither the European Union nor Clean Aviation JU can be held responsible for them. LL acknowledges partial support from project CSxAI (PID2024-157526NB-I00) funded by MICIU/AEI/10.13039/501100011033/FEDER, UE, project Maria de Maeztu (CEX2021-001164-M) funded by the MICIU/AEI/10.13039/501100011033, and from the European Commission Chips Joint Undertaking project No. 101194363 (NEHIL). GR acknowledges partial financial support received by the Grant DeepCFD (Project No. PID2022-137899OB-I00) funded by MICIU/AEI/10.13039/501100011033 and by ERDF, EU. Finally, all authors gratefully acknowledge the Universidad Politécnica de Madrid for providing computing resources on Magerit Supercomputer.
References [1]
Jiri Blazek. Computational fluid dynamics: principles and applications. Butterworth-Heinemann, 2015.
[2]
Steven L Brunton and J Nathan Kutz. Data-driven science and engineering: Machine learning, dynamical systems, and control. Cambridge University Press, 2022.
[3]
Emiliano Iuliano and Domenico Quagliarella. Proper orthogonal decomposition, surrogate modelling and evolutionary optimization in aerodynamic design. Computers & Fluids 84 (2013), pp. 327– 350.
[4]
Thomas Franz, Ralf Zimmermann, Stefan Görtz, and Niklas Karcher. Interpolation-based reducedorder modelling for steady transonic flows via manifold learning. International Journal of Computational Fluid Dynamics 28.3-4 (2014), pp. 106–121.
[5]
Marco Fossati. Evaluation of aerodynamic loads via reduced-order methodology. AIAA Journal 53.8 (2015), pp. 2389–2405.
[6]
Ian Goodfellow, Yoshua Bengio, and Aaron Courville. Deep Learning. Cambridge, MA: MIT Press, 2016. url: https://www.deeplearningbook.org.
[7]
Steven L Brunton, Bernd R Noack, and Petros Koumoutsakos. Machine learning for fluid mechanics. Annual review of fluid mechanics 52.1 (2020), pp. 477–508.
[8]
Soledad Le Clainche, Esteban Ferrer, Sam Gibson, Elisabeth Cross, Alessandro Parente, and Ricardo Vinuesa. Improving aircraft performance using machine learning: A review. Aerospace Science and Technology 138 (2023), p. 108354.
[9]
Jiaqing Kou and Weiwei Zhang. Layered reduced-order models for nonlinear aerodynamics and aeroelasticity. Journal of Fluids and Structures 68 (2017), pp. 174–193.
[10]
Xiaowei Jin, Peng Cheng, Wen-Li Chen, and Hui Li. Prediction model of velocity field around circular cylinder over various Reynolds numbers by fusion convolutional neural networks based on pressure on the cylinder. Physics of Fluids 30.4 (2018).
[11]
Saakaar Bhatnagar, Yaser Afshar, Shaowu Pan, Karthik Duraisamy, and Shailendra Kaushik. Prediction of aerodynamic flow fields using convolutional neural networks. Computational Mechanics 64 (2019), pp. 525–545.
[12]
Haizhou Wu, Xuejun Liu, Wei An, Songcan Chen, and Hongqiang Lyu. A deep learning approach for efficiently and accurately evaluating the flow field of supercritical airfoils. Computers & Fluids 198 (2020), p. 104393.
[13]
Nils Thuerey, Konstantin Weißenow, Lukas Prantl, and Xiangyu Hu. Deep learning methods for Reynolds-averaged Navier–Stokes simulations of airfoil flows. AIAA Journal 58.1 (2020), pp. 25– 36.
[14]
Vladyslav Rozov and Christian Breitsamter. Data-driven prediction of unsteady pressure distributions based on deep learning. Journal of Fluids and Structures 104 (2021), p. 103316.
[15]
Daniel Soler, Oscar Mariño, David Huergo, Martı́n de Frutos, and Esteban Ferrer. Reinforcement learning to maximize wind turbine energy generation. Expert Systems with Applications 249 (2024), p. 123502. 15
[16]
David Ramos, Lucas Lacasa, Eusebio Valero, and Gonzalo Rubio. Transfer learning-enhanced deep reinforcement learning for aerodynamic airfoil optimization subject to structural constraints. Physics of Fluids 37.8 (2025).
[17]
Jiaqing Kou et al. Aeroacoustic airfoil shape optimization enhanced by autoencoders. Expert Systems with Applications 217 (2023), p. 119513.
[18] Ángel Ladrón et al. A certifiable machine learning-based pipeline to predict fatigue life of aircraft structures. Engineering Failure Analysis (2025), p. 110334. [19]
Esther Andrés-Pérez and Carlos Paulete-Periáñez. On the application of surrogate regression models for aerodynamic coefficient prediction. Complex & Intelligent Systems 7.4 (2021), pp. 1991– 2021.
[20]
Christian Sabater, Philipp Stürmer, and Philipp Bekemeyer. Fast predictions of aircraft aerodynamics using deep-learning techniques. AIAA Journal 60.9 (2022), pp. 5249–5261.
[21]
Derrick Hines and Philipp Bekemeyer. Graph neural networks for the prediction of aircraft surface pressure distributions. Aerospace Science and Technology 137 (2023), p. 108268.
[22]
David Ramos, Lucas Lacasa, Fermı́n Gutiérrez, Eusebio Valero, and Gonzalo Rubio. FluidFlow: a flow-matching generative model for fluid dynamics surrogates on unstructured meshes. arXiv preprint arXiv:2604.08586 (2026).
[23]
Kurt Hornik. Approximation capabilities of multilayer feedforward networks. Neural networks 4.2 (1991), pp. 251–257.
[24]
Ziming Liu et al. Kan: Kolmogorov-Arnold Networks. arXiv preprint arXiv:2404.19756 (2024).
[25]
Andre Nikolaevich Kolmogorov. On the representation of continuous functions of several variables by superpositions of continuous functions of a smaller number of variables. American Mathematical Society, 1961.
[26]
Vugar E Ismailov. Addressing common misinterpretations of KART and UAT in neural network literature. Neural Networks (2025), p. 108361.
[27]
Shriyank Somvanshi, Syed Aaqib Javed, Md Monzurul Islam, Diwas Pandit, and Subasish Das. A survey on kolmogorov-arnold network. ACM Computing Surveys 58.2 (2025), pp. 1–35.
[28]
Khemraj Shukla, Juan Diego Toscano, Zhicheng Wang, Zongren Zou, and George Em Karniadakis. A comprehensive and FAIR comparison between MLP and KAN representations for differential equations and operator networks. Computer Methods in Applied Mechanics and Engineering 431 (2024), p. 117290.
[29]
Runpeng Yu, Weihao Yu, and Xinchao Wang. KAN or MLP: A fairer comparison. arXiv preprint arXiv:2407.16674 (2024).
[30]
Ali Kashefi. Kolmogorov-Arnold PointNet: Deep learning for prediction of fluid fields on irregular geometries. arXiv preprint arXiv:2408.02950 (2024).
[31]
Victor A Kich, Jair A Bottega, Raul Steinmetz, Ricardo B Grando, Ayano Yorozu, and Akihisa Ohya. “Kolmogorov-Arnold Networks for Online Reinforcement Learning”. 2024 24th International Conference on Control, Automation and Systems (ICCAS). IEEE. 2024, pp. 958–963.
[32]
David Picard. Torch. manual seed (3407) is all you need: On the influence of random seeds in deep learning architectures for computer vision. arXiv preprint arXiv:2109.08203 (2021).
[33]
Dieter Schwamborn, Thomas Gerhold, and Ralf Heinrich. “The DLR TAU-code: recent applications in research and industry”. ECCOMAS CFD 2006: Proceedings of the European Conference on Computational Fluid Dynamics, Egmond aan Zee, The Netherlands, September 5-8, 2006. Delft University of Technology; European Community on Computational Methods . . . 2006.
[34]
Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini. The graph neural network model. IEEE transactions on neural networks 20.1 (2008), pp. 61–80.
[35]
Gabriele Corso, Hannes Stark, Stefanie Jegelka, Tommi Jaakkola, and Regina Barzilay. Graph neural networks. Nature Reviews Methods Primers 4.1 (2024), p. 17.
[36]
Ziming Liu. pykan. Apr. 2024. url: https://github.com/KindXiaoming/pykan/tree/master.
[37]
Max Kuhn, Kjell Johnson, et al. Applied predictive modeling. Vol. 26. Springer, 2013.
[38]
Gene H Golub and Charles F Van Loan. Matrix computations. JHU press, 2013. 16
[39]
Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama. “Optuna: A Next-generation Hyperparameter Optimization Framework”. Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. 2019.
[40]
Lucas Lacasa et al. Towards certification: A complete statistical validation pipeline for supervised learning in industry. Expert Systems with Applications 277 (2025), p. 127169.
[41]
Benet Eiximeno, Beka Begiashvili, Arnau Miro, Eusebio Valero, and Oriol Lehmkuhl. pyLOM: Low order modelling in Python. 2022. url: https://github.com/ArnauMiro/UPM_BSC_LowOrder.
17