Building and Visualising Kolmogorov–Arnold Networks

Building and Visualising Kolmogorov–Arnold Networks | Towards Data Science

Machine Learning

Building and Visualising Kolmogorov–Arnold Networks

Exploring what happens when neural networks learn functions instead of weights.

Sam Everest

October 28, 2026

27 min read

Prefer us on Google

Share

Visualisation of a Chebyshev-based Kolmogorov–Arnold Network (ChebyKAN) modelling the Lorenz system, showing network connections and plots of learned functions.

Build a Chebyshev-based Kolmogorov–Arnold Network (ChebyKAN) and explore how it models the Lorenz system through interactive, reproducible experiments.

Two years ago , discussions about Kolmogorov-Arnold networks started to overtake my news feeds. These networks replace the scalar weights of a conventional neural network with learnable functions, this sounds like a relatively small change, but it raises an interesting possibility: instead of only asking if a network makes accurate predictions, can we actually explore the functions it learns, and perhaps turn the trained network back into an equation?

I wanted to understand KANs by building them myself and exploring what they could do. In the article I will work through the increasingly flexible edge functions in LabLayer (an early concept tool I’m building for prototyping and visualising networks) that I’ve been playing with. My first attempts are simple quadratics, B-splines moving to Chebyshev polynomials, bringing up to my latest experiments on recovering the Lorenz Equations from synthetic data.

What is a KAN?

Proposed in 2024, KANs, or Kolmogorov–Arnold Networks, are neural-network architectures inspired by the Kolmogorov–Arnold representation theorem (Liu et al. 2025). In simplified terms, the theorem states that a continuous multivariate function on a bounded domain can be represented through compositions and sums of continuous single-variable, or univariate, functions.

For our purposes, the important idea is that a complex function of several variables can be constructed from a collection of simpler one-dimensional functions. KANs take inspiration from this by replacing the scalar weights used in conventional neural networks with learnable univariate functions on the connections between nodes. This gives the network a way to build complex multivariate relationships through repeated composition and summation of simpler functions.

In practice the core idea of a KAN is conceptually quite simple. Broadly speaking, the main change from a conventional feed-forward neural network is to replace each scalar edge weight with a learnable univariate function. In a standard network, an incoming connection typically contributes a term of the form w j i x i w_{ji}x_i w j i ​ x i ​ , after which the node sums its inputs and applies an activation function. In a KAN, the edge instead applies a learned transformation ϕ j i ( x i ) \phi_{ji}(x_i) ϕ j i ​ ( x i ​ ) , and the node sums the transformed edge outputs. Mechanically, this can be viewed as two changes: first, moving the nonlinear transformation from the node to the individual edges; and second, parameterising those edge functions using B-splines or another suitable function family.

From Weights to Learnable Edge Functions

Before discussing splines my first step was to implement training on edges or more precisely training a parameterized function on edges as opposed to the traditional single value weight w. In a normal network we learn (or initiate) values for w w w and b b b , we apply them in a linear configuration w x + b wx+b w x + b , we sum ∑ w x + b \sum{wx+b} ∑ w x + b , then we apply a non-linear transformation known as the activation function such as ReLU or Tanh which we will call σ. This might be represented as:

y

=

σ

(

∑

i

w

x

+

b

)

y = \sigma\left( \sum_{i}w_{i}x_{i} + b \right)

​

→

j

\boxed{x_{i}\,\overset{\, w_{ji}\,}{\rightarrow}\sum_{i}\,\overset{\,\sigma\,}{\rightarrow}y_{j}}

In order to train on edges we need to make one change we need to apply the non-linearity to each input prior to the summation as opposed to the node output.

y = \sum_{i}\sigma\left( w_{i}x_{i} + b \right)

Our new order of operations for each input x i x_i x i ​ becomes:

ϕ

\boxed{x_{i}\,\overset{\,\phi_{ji}\left( x_{i} \right)\,}{\rightarrow}\sum_{i}\, \rightarrow \, y_{j}}

Here the nonlinear transformation is applied independently to each incoming edge before the results are summed. This is not yet a full KAN formation, for this we need to replace the scalar weight w and the fixed activation function with a learnable univariate edge function we might call ϕ \phi ϕ .

y = \sum_{i}\phi_{ij}\left( x_{i} \right) + b_{j}

ij

ϕ \phi ϕ connects input i i i to node j j j . This potentially allows for a wide range of transformations to be applied to each input. However, unlike a conventional network where each edge may be represented by a single scalar weight, every edge must now evaluate its own function ϕ i j \phi_{ij} ϕ ij ​ . The computational cost therefore depends both on the number of edges and on the complexity of the chosen function ϕ \phi ϕ , including the number of trainable parameters it contains. As a result, both forward evaluation and training can become considerably more expensive than in a conventional feed-forward network. Despite this fairly significant drawback, the flexibility gained from using learnable edge functions is interesting enough that it is worth exploring what advantages they can offer.

Building the first edge function

By this point I had progressed far enough to want to test it out.

For this first experiment I created a dataset based on the following function:

2

1

−

0.25

0.7

2\left( x_{1} - 0.25 \right)^{2} - \left( x_{2} - 0.7 \right)^{2} + 1

The dataset is deliberately simple, with only two inputs and a single output. The target function is a simple quadratic function of the two inputs, which allows us to easily visualise the learned transformations, and to provide an easy target to get Lablayer functional.

An early implementation of a synthetic dataset generator allowed me to generate the dataset.

Synthetic data

While my first priority was to get a simple edge function working, I also wanted to be able to visualise the learned transformation and really get a good understanding of what was happening. We have two opposing quadratic terms in the target function, so I created a simple edge function that is also quadratic, with two trainable parameters. The function is defined as:

p

,

0

y_{j} = \sum_{i}p_{ji,0}\left( p_{ji,1} + x_{i} \right)^{2} + b_{j}

with p i j , n p_{ij,n} p ij , n ​ representing our trainable parameters, n is just the parameter identifier.

At the moment these are created as node types using a basic equation language influenced by tex. This part of Lablayer is still under development, but the idea is that I can quickly define architecture mathematically and the system itself generates the corresponding code to implement things such as node types. Even in this simple example the equation started to feel overly complex, so I added ‘levels’ allowing me to define sub-expressions that can be used in the main node expression. In this case I defined the main expression as:

\sum_{i}\phi_{ji}\left( x_{i} \right) + b_{j}

If the summation is placed prior to the φ then the visualisation knows to draw edges as opposed to static weights.

I defined the sub-expression as:

[

]

p[ 0]\left( p[ 1] + x_{i} \right)^{2}

where p [ 0 ] p[0] p [ 0 ] and p [ 1 ] p[1] p [ 1 ] are the trainable parameters.

In the spirit of being able to control and specify all properties, I also defined the gradient of the edge function, which is used to compute the backpropagation updates. The gradient is defined by partial derivatives as:

∂

/

\partial y/\partial x_{i} = 2p[ 0]\left( p[ 1] + x \right)

\partial y/\partial p[ 0] = \left( p[ 1] + x \right)^{2}

\partial y/\partial p[ 1] = 2p[ 0]\left( p[ 1] + x \right)

\partial y/\partial b = 1

The partial derivative input boxes are created automatically from the parameters used in the function definition. This gives the system a great deal of flexibility: I can define almost any function I like, specify its gradient, and LabLayer will automatically generate the corresponding implementation. There is still considerable work to do to make this process more user friendly, particularly in supporting all of the required mathematical building blocks and allowing users to choose between explicitly defined gradients and automatic differentiation.

To demonstrate these differences in practice, I first trained a simple three-node ReLU network on the same dataset. The architecture is deliberately neither optimised nor particularly reliable, but it usually produces a reasonable approximation within a few training runs. In this conventional network, each connection transforms its input using the single scalar weight described above. The same is currently true of the output connections, as LabLayer does not yet have the UI options to configure the output layer. These scalar weights are straightforward to visualise: at any point in training, an edge can simply be represented by its current numerical value. LabLayer’s weight plots therefore show how each weight changes over the course of training. This provides a useful baseline for comparison with a KAN, where an edge can no longer be adequately described by a single number and instead needs to be visualised as the function it has learned.

A three node ReLU network trained on the synthetic data.

I've included the experiment links, all parameters and reference to the data should be encoded so you can explore and modify what I did. This tool is part hobby, part concept so the links may break as it changes.

Lablayer

Defining a Quadratic Edge function in an early but functional prototype of the node builder I used to run this experiment

Learned Edge Functions

We can easily switch the layer over to our quadratic function and now see the differences. The above image shows the learned edge functions for each of the two inputs after training for 1000 epochs, the learned edge functions are visible and plotted out. An interesting observation in this case is that we only need one node: each input is transformed independently.

Learned edge functions for the two inputs.

Even with this deliberately restricted quadratic edge function, the learned edge plots begin to reveal structure in the data. In this example the target contains a x 2 2 x_2^2 x 2 2 ​ term, and certainly the transformations are showing this. Rather than inspecting only scalar weights we are now able to inspect the transformations on each connection and see what they are learning.

Comparison between the learned model and the exact target function.

The surface view shows the learned transformations on the left and the target function on the right. We have a small learning rate as to avoid explosive gradients. We can see that the learned model and the ‘true’ model are very similar.

Error surface for the learned model, plotted as predicted output minus the exact target function.

The error across the surface only shows small deviations at the edge of the domain, curiously the error itself is quadratic in appearance, which is inline with what you would expect where the differences in co-efficient become largest at the edges of the domain due to the polynomial term.

Interestingly it is fairly straightforward to collapse the entire architecture into a polynomial. Once trained there is really no reason to use it as a feedforward network, instead we can simply evaluate the polynomial directly. The following figure shows the learned polynomial, which is also an excellent approximation of the target function.

The collapsed architecture summary

0.991

1.393

2.053

0.994

0.613

y = - 0.991x_{2}^{2} + 1.393x_{2} + 2.053x_{1}^{2} - 0.994x_{1} + 0.613

as opposed to the target function:

1.4

0.635

y = - x_{2}^{2} + 1.4x_{2} + 2x_{1}^{2} - x_{1} + 0.635

This example is intentionally simple: each edge has only two trainable parameters and is restricted to a quadratic transformation. The natural question is what happens if the function associated with each edge is itself made flexible enough to learn a much wider family of shapes. This is the idea that takes us to the next and final step of Kolmogorov-Arnold Networks.

Making the edges flexible with B-splines

So far, I’ve managed to implement a learnable function on an edge, the function has been tailor made for the simple problem I set myself. Indeed not every problem will be a set of contrasting quadratics, so we might want to capture more possibilities per edge.

Enter basis splines, or B-spline, is a piecewise polynomial function defined over a set of knots. The degree of the polynomial and the number and spacing of the knots determine the flexibility of the resulting spline representation.

With a sufficiently dense set of knots, a B-spline basis can approximate a continuous function to arbitrary accuracy. This property made splines suitable for the original KAN proposal in which they were used as a weighted combination of basis functions in order to construct a learnable univariate function on each edge of the network. It’s important to note at this step, splines are not the only possibility, and we will check out another further on.

In the preceding example, a single quadratic function was used to transform the input. A B-spline follows a similar principle, but instead of applying one polynomial globally, different polynomial expressions are used over different regions of the input domain. This piecewise structure is conceptually similar to ReLU, using a conditional statement to separate the input domains:

R

e

L

U

{

≤

ReLU(x) = \begin{cases} x, & x > 0 \\ 0, & x \leq 0 \end{cases}

LU

For our testing we are going to use the following cubic B-spline basis function:

4

6

3

<

otherwise

∣

c

h

B_{r}\left( x_{i} \right) = \begin{cases} \frac{4 - 6u_{ir}^{2} + 3u_{ir}^{3}}{6} & 0 \leq u_{ir} < 1 \\ \frac{\left( 2 - u_{ir} \right)^{3}}{6} & 1 \leq u_{ir} < 2 \\ 0 & \text{otherwise} \end{cases} ,\quad u_{ir} = \left| \frac{x_{i} - c_{r}}{h} \right|

⎩

⎨

⎧

where c r c_r c r ​ denotes the centre of basis function r r r , and h h h controls its width. Each basis function therefore has local support: it only contributes when x i x_i x i ​ lies within a limited region around c r c_r c r ​ .

The spline basis functions are combined to form the learnable edge function:

a

\phi = \sum_{r}a_{jir}B_{r}\left( x_{i} \right)

where aⱼᵢᵣ controls the contribution, or amplitude, of basis function r r r on the edge connecting input i i i to node j j j .

The output of node j j j is then

The coefficients a j i r a_{jir} a j i r ​ are learned during training. Depending on the implementation, the spline centres cᵣ and width h may either be fixed by the spline grid or treated as additional learnable parameters.

Increasing the number of spline basis functions increases the flexibility of the edge function, since more local regions of the input domain can be modelled independently. However, this also increases the number of parameters that must be learned.

Returning to our simple two-input dataset, we can now replace the quadratic edge function with a B-spline edge function. The following figure shows our single cubic B-spline basis function prior to training.

Example of a cubic spline basis function.

Defining the KAN basis.

There is still work to be done to make the node definition more user friendly, but the basic functionality is there. The above shows how the B-spline basis functions are defined, and how they are combined to form the edge function. Using the level system, the user can define the node expression in terms of the edge function and then also define the edge function and required derivatives for backpropagation. The system will then automatically generate the corresponding code to implement the node type.

I explicitly defined the derivatives used for backpropagation, although these can also be calculated using lower-degree splines. The spline pieces are constructed so that they join smoothly, with matching derivatives where neighboring pieces meet. This makes them somewhat better behaved than a piecewise function such as ReLU, which contains a sharp change in gradient at its breakpoint.

Currently the definitions actually input are:

This is the main node expression, which sums the edge functions and adds a bias. The preceding sum to the phi helps us later identify that the edges are plottable functions. The edge function is defined as:

\sum_{r}a_{jir}B_{r}\left( x_{i} \right)

This is our ϕ \phi ϕ function, which sums the B-spline basis functions weighted by their corresponding coefficients. The B-spline basis function is defined as:

\begin{cases} \frac{4 - 6u_{ir}^{2} + 3u_{ir}^{3}}{6} & 0 \leq u_{ir} < 1 \\ \frac{\left( 2 - u_{ir} \right)^{3}}{6} & 1 \leq u_{ir} < 2 \\ 0 & \text{otherwise} \end{cases} ,\quad u_{ir} = \left| \frac{x_{i} - c_{r}}{h} \right|

Finally the B-spline input needs some simplification.

The next step is to train the network on the dataset and visualise the learned edge functions.

Training using B-Splines on the dataset.

Learned KAN edge functions. Each panel shows the functions learned for inputs x 1 x_1 x 1 ​ and x 2 x_2 x 2 ​ by each node in the hidden layer.

Comparison between the final KAN model response and the target function over the two-dimensional input domain.

We can see that the learned edge functions are now capturing much more complex information. In this case it is actually far too complex for the problem. As we learn a function for each input, with each input only being a simple quadratic we can reduce the number of nodes to one and see if the network can learn the target function.

Learned edge function.

The learned edge functions for the two inputs.

Comparison between the learned KAN model response and the target function.

The learned output is quite close to the target, however inspecting the learned edge functions we can see that they still differ significantly from the domains target function. This is somewhat to be expected as our initial edge function was selected on the basis that it could effectively model the target function. The B-spline edge function is far more flexible, and can learn a wide range of functions, but it is not guaranteed to learn the target function exactly.

We now experience another issue.

Attempting to collapse the KAN equations.

If we try to collapse the network into a single polynomial, we encounter a new difficulty. The learned edge functions are constructed from piecewise-polynomial B-spline basis functions, so their exact algebraic representation is also piecewise. Although the combined spline may approximate a simple global function extremely well, recovering the exact function learned by the network requires a different polynomial expression over each region of the input domain. Even in this single-node example, the result must be represented by multiple local polynomial expressions. As the number of inputs, knots and layers increases, the symbolic representation can therefore become substantially more complicated.

The original KAN paper (Liu et al. 2025) discusses in depth methods for extracting symbolic representations from learned edge functions. However, these methods generally involve fitting or “snapping” the numerically learned function to a candidate symbolic form, followed by further optimisation. In this sense, the symbolic representation is not necessarily the exact function learned during training, but a simplified analytical approximation to it. If KANs are to be used as tools for interpretability and equation recovery, there is therefore considerable value in architectures whose learned functions can be converted directly into explicit symbolic expressions with little or no post-hoc modification.

This issue has been very much acknowledged. While B-Splines have provided a flexible entry, they are definitely not required. If we want to learn a global function from the architecture we can use a different function family. Enter Chebyshev polynomials (SS et al. 2024):

Chebyshev polynomials: a different approach

Chebyshev polynomials are a sequence of orthogonal polynomials that can be defined recursively. They have several properties that make them useful for function approximation, including their ability to minimize the maximum error between the polynomial and the target function over a given interval.

They can be calculated using the following recurrence relation:

T

T_{0}(x) = 1,\quad\quad T_{1}(x) = x

n

.

T_{n + 1}(x) = 2xT_{n}(x) - T_{n - 1}(x).

Essentially the first few are

T_{0}(x) = 1

T_{1}(x) = x

T_{2}(x) = 2x^{2} - 1

T_{3}(x) = 4x^{3} - 3x

Chebyshev polynomials on [ − 1 , 1 ] [-1,1] [ − 1 , 1 ] .

8

T_{4}(x) = 8x^{4} - 8x^{2} + 1

Chebyshev polynomials are well behaved over the interval [ − 1 , 1 ] [-1,1] [ − 1 , 1 ] . Outside this region their magnitude can grow rapidly, particularly for higher-order polynomials. For this reason, inputs are normally scaled to [ − 1 , 1 ] [-1,1] [ − 1 , 1 ] before being passed to the Chebyshev basis.

We can construct a flexible function by taking a linear combination of the Chebyshev polynomials. An edge function of degree n can therefore be written as:

k

\phi(x) = \sum_{k = 0}^{n}c_{k}T_{k}(x)

or for degree 3:

f

7

5

f(x) = \frac{7}{2}T_{0}(x) - \frac{5}{4}T_{1}(x) + 3T_{2}(x) - \frac{3}{4}T_{3}(x) + \frac{1}{2}T_{4}(x)

Chebyshev KAN node definition.

In Lablayer we can enter this similar to the B-spline example, there is still extensive work to be done on the node builder, but the basic functionality is there. We define level 1 as:

By now this should be familiar, it is the main node expression that sums the edge functions and adds a bias. The edge function is defined as:

\sum_{r}a_{jir}T_{r}\left( x_{i} \right)

Level 2 now includes T and r opening level 3 to define these functions in which we will place the Chebyshev polynomials. Level 3 is defined as:

\begin{cases} 1 & r = 0 \\ x_{i} & r = 1 \\ 2x_{i}^{2} - 1 & r = 2 \\ 4x_{i}^{3} - 3x_{i} & r = 3 \end{cases}

At the moment we have an order selector for r however this will likely be removed as has become somewhat redundant. Using this node the network will iterate over the order of the polynomial and sum the results.

Finally we have the derivatives for backpropagation, which are defined as:

12

\begin{cases} 0 & r = 0 \\ 1 & r = 1 \\ 4x_{i} & r = 2 \\ 12x_{i}^{2} - 3 & r = 3 \end{cases}

At present, I have implemented the polynomial explicitly in a case format within LabLayer, although I intend to revisit this approach in future. Explicit derivative specification will remain part of the final system because one of the core design principles of LabLayer is to allow arbitrary functions to be defined while keeping the assumptions used during training visible and explicit. The aim is to make it easy for me to prototype unusual architectures without hiding too much of the underlying mathematics. I also plan to add automatic differentiation as an optional alternative. For most use cases this should save me considerable time, while still allowing for me to manually specify derivatives where greater control or transparency is useful.

Training this on our existing problem:

Chebyshev KAN network.

Learned Chebyshev edge functions for each hidden node. Each panel shows the transformations learned for inputs x 1 x_1 x 1 ​ and x 2 x_2 x 2 ​ .

Comparison between the final ChebyKAN model response and the target function.

0.001

0.002

1.001

1.401

2.000

1.000

y = 0.001x_{1}^{3} - 0.002x_{2}^{3} - 1.001x_{2}^{2} + 1.401x_{2} + 2.000x_{1}^{2} - 1.000x_{1} + 0.635

Chebyshev KAN collapsed.

Show us a fantastic result!! x 3 x^3 x 3 terms are minimised and the learned function is very close to the target function! This is great news, as we have on our very early and simplistic tests managed to create a more general architecture that in itself has become the function we are trying to learn.

Can it recover a physical system?

Emboldened by how smoothly things had been going so far, but not yet brave enough to dive into something complex, I decided to try an equation describing a real physical system, I chose a simple equation of motion many of us might have seen in school:

s

s = ut + \frac{1}{2}at^{2}

For this experiment, the dataset is generated using fixed values of u and a, giving:

1.5

2x + 1.5x^{2}

over the domain:

∈

x \in [ 0,5]

This creates a very straightforward one-in, one-out dataset.

The same ChebyKAN architecture is used, but the input is normalised to the interval [−1, 1]. This is important because Chebyshev polynomials are numerically well behaved over this domain. Without this normalisation, higher-order terms can grow rapidly and the learned parameters may diverge.

Training on the equation of motion.

Learned Chebyshev edge functions for the five hidden nodes, evaluated over the original input domain.

Comparison between the reconstructed ChebyKAN output and the target equation of motion. They are very close.

After training, the collapsed network gives:

2.88

×

10

9

1.500

1.02

y = 2.88 \times 10^{- 9}x^{3} + 1.500x^{2} + 2.000x + 1.02 \times 10^{- 8}

The cubic and constant terms are effectively zero, leaving:

^

≈

\widehat{y} \approx 1.5x^{2} + 2x

which recovers the generating equation almost exactly.

Testing on the Lorenz Equations

Feeling much braver with the successes so far I decided to try and recover the Lorenz equations. The equations were developed to model atmospheric convection and are a classic example of a system that can exhibit chaotic behaviour. The equations themselves are quite straightforward but combined they create quite a complex trajectory.

˙

28

z

\begin{matrix} \dot{x} & = 10y - 10x \\ \dot{y} & = 28x - xz - y \\ \dot{z} & = xy - \frac{8}{3}z \end{matrix}

using the standard parameters:

ρ

β

\sigma = 10,\quad\quad\rho = 28,\quad\quad\beta = \frac{8}{3}.

For this experiment I created the data externally and uploaded it to the tool. The current function generator is intended primarily to allow me to do rapid creation of simple test datasets and does not yet support coupled multi-equation systems such as the Lorenz equations. In practice, I would typically be wanting to upload my own prepared datasets. I was able quickly run the existing single layer ChebyKAN architecture.

Training my architecture on the Lorenz equations.

Which produced the following collapsed expressions.

7.55

1.28

7.62

0.005

6.32

7.05

0.058

9.991

9.996

1.117

{\widehat{y}}_{1} = 7.55 \times 10^{- 5}z^{3} + 1.28 \times 10^{- 5}y^{3} - 7.62 \times 10^{- 6}x^{3} - 0.005z^{2} + 6.32 \times 10^{- 4}y^{2} + 7.05 \times 10^{- 4}x^{2} + 0.058z + 9.991y - 9.996x - 1.117

0.012

0.079

0.108

7.70

0.006

2.416

6.572

0.925

14.351

{\widehat{y}}_{2} = - 0.001z^{3} + 0.012y^{3} - 0.079x^{3} + 0.108z^{2} - 7.70 \times 10^{- 4}y^{2} + 0.006x^{2} - 2.416z + 6.572y - 0.925x + 14.351

2.75

5.34

0.142

0.235

1.052

3.846

0.104

0.064

1.282

{\widehat{y}}_{3} = - 0.005z^{3} - 2.75 \times 10^{- 4}y^{3} + 5.34 \times 10^{- 4}x^{3} + 0.142z^{2} + 0.235y^{2} + 1.052x^{2} - 3.846z + 0.104y - 0.064x + 1.282

First Lorenz equation

True equation:

\dot{x} = 10y - 10x

Dominant learned terms:

\widehat{\dot{x}} \approx 9.991y - 9.996x - 1.117

The two principal coefficients are recovered very closely. The remaining constant offset is comparatively small, and the higher-order terms in the full expression are negligible.

Second Lorenz Equation

\dot{y} = 28x - xz - y

\widehat{\dot{y}} \approx - 0.079x^{3} - 2.416z + 6.572y - 0.925x + 14.351

For the first time the equation is possibly unsurprisingly not recovered at all.

Third Lorenz Equation

\dot{z} = xy - \frac{8}{3}z

\widehat{\dot{z}} \approx 1.052x^{2} + 0.235y^{2} + 0.142z^{2} - 3.846z + 1.282

Again, the governing equation is not recovered.

This result is consistent with a limitation of the current single-layer architecture. Each input is transformed independently and the resulting edge functions are summed, giving a model of the form:

f(x,y,z) = \phi_{x}(x) + \phi_{y}(y) + \phi_{z}(z)

Such an architecture can represent additive relationships, but it cannot directly construct interaction terms such as:

xz

or

xy

This explains why the first Lorenz equation, which is purely additive, can be recovered closely, while the second and third equations are not.

The next step was to move towards an ideal function-reclamation network, so that I could at least demonstrate that recovery of the underlying equations was possible. To make interaction terms such as x y xy x y and x z xz x z mechanically representable, I added an additional layer to the network. For example, if the first layer produces h = x + y h=x+y h = x + y , applying the second-order Chebyshev polynomial T 2 ( h ) = 2 h 2 − 1 T₂(h) = 2h² − 1 T 2 ​ ( h ) = 2 h 2 − 1 produces 2 x 2 + 4 x y + 2 y 2 − 1 2x² + 4xy + 2y² − 1 2 x 2 + 4 x y + 2 y 2 − 1 . The additional composition therefore allows cross terms to emerge naturally. However, this introduced new training difficulties. In practice, the deeper architecture required an additional normalisation layer to prevent gradients from exploding, and the learning rate had to be reduced substantially to 0.00001 to keep training stable.

y_{j} = \sum_{k}\phi_{jk}^{(2)}\left( \sum_{i}\phi_{ki}^{(1)}\left( x_{i} \right) + b_{k} \right) + b_{j}

x x x slice view of the ChebyKAN

y y y slice view of the ChebyKAN

z z z slice view of the ChebyKAN

Initially these looked extremely promising.

After expanding the trained ChebyKAN into an explicit polynomial, the full expression contains a large number of higher-order terms. To make the learned structure easier to inspect, the dominant terms can be retained.

For this example, retaining terms with coefficient magnitude |c| ≥ 0.1 gives:

21.272

5.604

0.734

0.697

0.304

2.268

0.183

1.011

3.022

\begin{matrix} \widehat{\dot{x}} \approx & 21.272x - 5.604y + 0.734z - 0.697x^{2} + 0.304xy \\ & - 2.268xz - 0.183z^{2} + 1.011yz + 3.022 \end{matrix}

26.880

1.911

2.404

1.179

0.842

0.495

0.143

0.159

0.557

1.366

\begin{matrix} \widehat{\dot{y}} \approx & 26.880x - 1.911y + 2.404z - 1.179x^{2} + 0.842xy \\ & - 0.495xz - 0.143y^{2} + 0.159yz - 0.557z^{2} - 1.366 \end{matrix}

2.769

1.343

7.373

2.112

0.317

0.488

0.526

0.128

16.174

\begin{matrix} \widehat{\dot{z}} \approx & - 2.769x - 1.343y - 7.373z + 2.112x^{2} + 0.317xz \\ & + 0.488yz + 0.526z^{2} - 0.128x^{2}z + 16.174 \end{matrix}

The true Lorenz equations, using the standard parameters σ = 10, ρ = 28, and β = 8/3, are:

\dot{y} = 28x - y - xz

The dominant learned coefficients can therefore be compared directly with the known governing coefficients.

Table 1: Comparison of the physically relevant coefficients recovered from the expanded ChebyKAN with the true Lorenz-system coefficients.

Lorenz term

True coefficient

Learned coefficient

Relative error

Comment

ẋ:x

-10.000

313.0%

Wrong sign

ẋ:y

10.000

-5.604

156.0%

ẏ:x

28.000

4.0%

Very good

ẏ:y

-1.000

-1.911

91.1%

Correct sign, poor magnitude

ẏ:xz

-0.495

50.5%

Correct sign, approximately half magnitude

ż:xy

0.075

92.5%

Poor recovery

ż:z

-2.667

-7.373

176.4%

After repeating the experiment several times, it became clear that the trainer was primarily learning relationships sufficient to reproduce the sampled Lorenz attractor. The variables along a Lorenz trajectory are strongly correlated and occupy only a restricted region of the full state space. Consequently, multiple different polynomial expressions can produce similar predictions over the observed trajectory, making recovery of the original governing equation poorly constrained. I had in effect predicted the trajectory but had an incorrect underlying relationship.

To test whether this was a limitation of the architecture or of the sampled data, I repeated the experiment using independently sampled points throughout the state space. Armed with the power of prior knowledge, along with the necessity of needing to do 100 runs, I also limited the Chebyshev polynomials to x², and reduced the layers to three nodes each.

Table 2: Comparison between the true Lorenz coefficients, the mean learned coefficients across all runs, and the five runs with the lowest overall coefficient error.

Run / Summary

Mean error

True Lorenz

—

Mean ± SD

-10.293 ± 0.646

9.757 ± 0.381

26.075 ± 5.681

0.250 ± 1.318

-0.940 ± 0.225

0.970 ± 0.020

-3.902 ± 1.682

Run 25

-9.984

9.903

27.393

-1.037

-0.986

0.992

-2.629

1.52%

Run 66

-9.986

9.926

27.660

-0.930

-0.962

0.995

-2.675

1.95%

Run 75

-9.990

9.927

27.753

-0.914

-1.012

0.998

-2.771

2.24%

Run 80

-10.101

9.889

28.366

-0.860

-0.978

0.997

-2.725

3.16%

Run 2

-10.276

10.013

27.473

-0.868

-0.991

0.993

-2.833

3.68%

The results are striking: the underlying Lorenz equations are beginning to emerge directly from the collapsed network. Several individual runs recover all seven physically relevant coefficients to within only a few percent of their true values.

Distributions of recovered coefficients across 100 independently initialised ChebyKAN models. The dashed line represents the true Lorenz coefficient, while the dotted line represents the median recovered coefficient.

Prediction is not the same as equation recovery

By this point in my journey I had only scratched the surface of what might be possible, but I had started to gain a much better understanding of KANs. My experiments had also led me rather naturally down the much more specific rabbit hole of equation recovery.

The challenge is that, when approximating real physical systems, there may be many mathematical relationships that reproduce the observed data reasonably well. A single prediction loss only measures how closely the model matches the supplied outputs; it does not necessarily tell us whether the model has recovered the relationship that generated them.

The way the data are sampled matters too. If the inputs only occupy a restricted region of the possible state space, relationships between variables in that region can provide alternative ways of reproducing the same outputs. A network may therefore learn the structure of the sampled data rather than the underlying governing equation.

To explore this, I plotted the final prediction loss from each run against the accuracy with which it recovered the seven relevant Lorenz coefficients.

For each run I calculated the mean relative error across the seven physically relevant coefficients,

E_{c} = \frac{1}{7}\sum_{i = 1}^{7}\left| \frac{{\widehat{c}}_{i} - c_{i}}{c_{i}} \right|,

and report coefficient recovery accuracy as

100

max

⁡

A_{c} = 100\max\left( 0,1 - E_{c} \right).

Final training loss versus recovery accuracy of the seven Lorenz coefficients across repeated ChebyKAN runs. Loss is shown on a logarithmic scale; the five runs with the lowest overall coefficient error are labelled.

The result supports an important distinction: minimising prediction loss does not necessarily result in the best recovery of the underlying governing equation. Lower-loss models generally recovered the Lorenz coefficients more accurately, but the relationship was far from exact. Several models achieved very low prediction loss while recovering comparatively poor coefficients, while some higher-loss models reproduced the governing coefficients more faithfully.

This distinction becomes particularly important when neural networks are used as scientific models. A network may reproduce the observed outputs extremely well while relying on a substantially different internal mathematical relationship. In other words,

good predictive performance

⇏

correct governing relationship

\text{good predictive performance} \nRightarrow \text{correct governing relationship}.

Additional constraints can of course be introduced to encourage the network towards physically meaningful solutions. Physics-informed neural networks (PINNs), for example, incorporate known physical relationships into the optimisation process. But what interested me here was how far the network could get without being explicitly forced towards the correct equation.

And, despite the failures, there had already been some encouraging results.

There is an important caveat, though: so far I had been armed with a considerable amount of prior knowledge. I knew the form of the Lorenz equations, I knew which coefficients I was looking for, and I could design the network and experiments around them. Recovering an unknown physical law from data is a much harder problem.

What I learned building this

Exploring KANs so far I was surprised by how easy they actually are to implement, the hardest part by a long margin was the B-splines. I implemented them in my own prototyping tool that I wrote from scratch using no ML libraries. I’ve included links if you wish to try it with the warning that it is really only a toy project that may or may not work for you, you are welcome to contact me if you want to discuss this further. There is a pytorch implementation of KANs, in the references (Liu 2024) along with a ChebyKAN library you might try (SynodicMonth 2024).

My exploration so far led me to an interesting alternative use of neural networks: rather than treating the trained network as the final predictive model, the network can instead be used as a mechanism for discovering a function that can subsequently be inspected, simplified, validated, and reused independently of the network.

This leads naturally into the question of interpretability. Terms such as ‘interpretable’, ‘black box’, and ‘explainable’ are not defined by a clean mathematical framework, and can vary in their use (Molnar 2025). For the purpose of this work, preceding B-spline example illustrates the most common usage. Even for a deliberately simple problem, the network can become quickly explosively complex algebraically. Although it may still be possible to obtain a symbolic description using curve fitting or symbolic regression, this introduces an additional approximation step: the symbolic expression becomes an approximation of the function already approximated by the trained network.

By contrast, the Chebyshev network allows for a direct algebraic representation of the transformation performed by the trained model. Because the edge functions are expressed as finite combinations of known polynomial basis functions, the network can be collapsed exactly, apart from normal floating-point precision, into a global polynomial expression. In our synthetic examples we are able to recover the initial ground truth although certainly as the interactions increase we do trade off the generality . Far more exploration can be done here.

Having access to this explicit function gives us considerably more information than prediction accuracy alone. We can inspect how changes in each input affect the output, determine which terms dominate the learned relationship, and compare the learned structure with existing domain knowledge. For example, if theory suggests that a relationship should be sinusoidal, we can inspect whether the model has actually learned a sinusoidal dependency or whether it has instead found some alternative function that achieves an acceptable loss without representing the expected physical behavior.

Taken further, this suggests that the architecture itself could be interrogated as a source of scientific information. If a functional form is partially known but one or more coefficients are uncertain, the trained model may provide estimates of those parameters directly. Likewise, terms that consistently collapse towards zero may indicate that parts of the architecture are redundant, potentially allowing the model to be simplified and its computational cost reduced. In this sense, the aim is not only to explain a prediction after it has been made, but to investigate whether the structure learned by the network can itself be converted into an interpretable mathematical model.

###### Note: All images are by the author.

References

Liu, Ziming. 2024. “Pykan: Kolmogorov-Arnold Networks.” https://github.com/KindXiaoming/pykan .

Liu, Ziming, Yixuan Wang, Sachin Vaidya, Fabian Ruehle, James Halverson, Marin Soljačić, Thomas Y. Hou, and Max Tegmark. 2025. “KAN: Kolmogorov-Arnold Networks.” https://arxiv.org/abs/2404.19756 .

Molnar, Christoph. 2025. Interpretable Machine Learning: A Guide for Making Black Box Models Explainable . 3rd ed. https://christophm.github.io/interpretable-ml-book .

SS, Sidharth, Keerthana AR, Gokul R, and Anas KP. 2024. “Chebyshev Polynomial-Based Kolmogorov-Arnold Networks: An Efficient Architecture for Nonlinear Function Approximation.” https://arxiv.org/abs/2405.07200 .

SynodicMonth. 2024. “ChebyKAN: Kolmogorov-Arnold Networks Using Chebyshev Polynomials.” https://github.com/SynodicMonth/ChebyKAN .