TL;DR: Neural networks can easily learn to regress cubic polynomials of type , but how do they do that? Results show that they internally access a variable , which is surprisingly close to the substitution used in the first step of the Cardano method for solving cubic polynomials.
This project was developed as part of Wolfram Summer Research Institute and first published, together with the source code, on Wolfram Community.
How do neural networks regress cubic polynomials? Before answering this question, let's start with a more trivial question: can neural networks (NNs) regress cubic polynomials?
The task is defined as follows: given an input combination compute . Notice that the NN is not tasked to find the roots of that polynomial, which would presumably be a much harder task. We further limit the domain in which the NN is trained to operate for real values .
We choose a simple multi-layer perceptron (MLP) infrastructure that takes a 3-dimensional input vector and passes it through four hidden layers. Each layer has a linear component followed by a ReLU activation. The dimensions of the layers are, respectively, 64, 64, 64 and 15. The final 15-dim vector is passed through a linear layer that outputs a single real number y.
The model is able to achieve lossless regression (corresponding to a loss equal to on validation dataset) in less than 15 minutes of training on my M2 Pro CPU. The figure below shows the NN regression (left) and the expected regression (right) for the polynomial . The accuracy achieved by the NN is such that I had to create two different graphs. If I were to overlay them, the NN regression would completely ablate the underlying curve.
So, how did the NN learn to regress a cubic polynomial? A dummy way to perform regression would be to store a lookup table that associates an output value to every possible combination of inputs (x,a,b).
Let's try to calculate the memory space that such a lookup table would require. Given an integer range length and a step size , the lookup table needs to store combinations in memory. If we consider a step size of 0.001 the lookup table space is equal to units of memory. On the other hand, the actual memory space at disposal of the NN is described by the number of its parameters P. In this particular case, so less than units of memory.
This quick analysis suggests that the model must have developed some abstract and compressed understanding of the polynomial world to successfully perform cubic polynomial regression.
We can also try to use the NN to regress cubic polynomials with coefficients outside its training range. We see that the regression starts to lose some precision. Nevertheless, the result is not completely off range, demonstrating that the model has developed some understanding of the nature of cubic polynomials beyond its training data distribution.
A good question to ask is how we humans learn to regress cubic polynomials. When I look at a polynomial of type , first I identify that, given the highest term , the shape of the function goes from the bottom-left of the graph (as x goes to negative infinity) to the top-right (as x goes to positive infinity), while the is the y-intercept. A more sophisticated observer would also notice that the discriminant of the first derivative is a function of the parameter , so, for the curve is always increasing, while, for the curve will have a sort of S-shape. Lastly, the curve is symmetric over its inflection point located at . All these heuristic techniques are abstractions, that we, computationally bounded humans, have developed to deal with the sheer amount of information required to regress the world around us.
NNs, even though they might have brains bigger than ours, are also computationally bounded. So they must develop some sort of abstraction to survive.
I'm interested to see if and which abstractions they use to perform this task. Are these abstractions similar to the ones that we developed in centuries of calculus or are these completely different? Can we even aspire to interpret them? We notice that the training process for the regression task is completely unsupervised: the NN is not fed with any human-induced hints about how they should or should not approach the problem.
To find an answer to these questions, we execute an MRI scan over the trained NN and extract the features that the NN activates when performing the task. Surprisingly, the NN learned a coordinate that resembles a transformation performed in the first step of Cardano's Method for solving cubic equations. Further manipulation experiments are run to understand how these features influence the behavior of the NN .
What are the techniques used to interpret the behavior of an NN? Some might accept the cliché that NNs are black-box and we cannot really understand their inner workings. If you don't accept this as absolute truth and you are willing to at least try to poke into the brain of an NN there are two main white-box interpretability techniques.
The first technique relies on the analysis of individual neurons and answers the question "what's the concept associated with a specific neuron?". To answer this question, you can feed various inputs to your NN, measure the neurons' activations and try to draw correlations such as "neuron #14 fires for and it is invariant to ". Unfortunately, this often turns out to be an impossible task. Research has shown that neurons are polysemantic. Namely, the model reuses the same neuron for unrelated concepts. This happens because of the training pressure to compress a number of concepts that is inevitably greater than the number of neurons at disposal. In large language models, it was found that the same neuron fires both for "baseball" and "academic citations". The clean 1-to-1 mapping between neurons and concepts that would make interpretability easy simply does not exist in most models. Similar noisy results were gathered in early efforts in trying to interpret the cubic polynomial regression NN.
The second technique is to train a Sparse Autoencoder (SAE) on top of an activation layer of the NN to extract a more fine-grained representation of their concepts. From now on, I am gonna use the term feature, to indicate the NN-equivalent of human concepts. To understand what a SAE is, it is worth first explaining what a standard autoencoder is. A regular autoencoder is an MLP trained to compress and reconstruct an input vector. As an example, it might take a 100-dim vector as input, pass it through the encoder layer to compress it to a 50-dim vector and then feed it through the decoder to reconstruct the original 100-dim output vector. The encoder and the decoder are trained to minimize the reconstruction loss. Of course, the reconstruction cannot be fully accurate and some information is inevitably lost in the compression step.
A SAE is like a standard autoencoder with the difference that it encodes the input vector into a larger dimension vector. If the training goal would simply be to minimize the reconstruction loss, the task would be trivial: the SAE would learn to use the identity matrix as encoder. But that wouldn't tell us anything interesting. Therefore, we add a sparsity penalty to the training loss: the SAE is incentivized to encode the input vector into a sparse intermediate vector, namely a vector that can only have a limited number of non-zero elements.
Going back to our context, we train a SAE on top of the last 15-dim activation layer of our NN. More concretely, we train (using an external python library) a 64-dim SAE with TopK=4. This means that it would be trained to learn 64 features and each encoded activation will be expressed in, at most, 4 features. To train the SAE, we sample 200k activations from the last activation layer of our NN for random input combinations . The learned parameters are and (used for encoding) and (used for decoding). More specifically, given an activation vector a, the sparse feature activation vector is defined as and the reconstructed activation vector is defined as .
To assess the quality of a trained SAE, we measure the coefficient of the reconstruction on held-out activation vectors. The coefficient of the trained SAE is 0.995, indicating a close-to-lossless reconstruction.
But why even bother training a SAE? It turns out that, while the activation vector is dense and difficult to interpret, the vector breaks down the polysemantic neurons into monosemantic components. This is possible thanks to the two key mechanisms that regulate the SAE: dimension expansion and sparsity. The expanded dimension provides enough space for every single feature to have its own dedicated "slot". The sparsity requirement forces disentanglement.
So, how can SAE help us understand what's going on inside an NN? This is not so straight-forward.
The previous figure tells us that, for the given set of input , layer-4 activates four features: feature #1, #4, #14 and #18. And the second and fourth features are the ones that activate the most.
So this, by itself, doesn't help us to understand how an NN thinks. The SAE does not provide any human-interpretable label to interpret a feature. In fact, the SAE is trained in a completely unsupervised fashion and there's no training pressure that privileges interpretability. So, how do we make sense out of these features?
Let's look at how features are interpreted in the context of LLMs. Assume that we have trained a SAE on top of the activations of layer . To understand what feature corresponding to the index represents, the practice is to simply look at the inputs that maximally activate that feature and give a gut reaction on their interpretability. For example, Anthropic trained SAEs on top of Claude Sonnet and found a feature that activates on contexts (text or images) related to the Golden Gate Bridge. Much progress has been made since then and now we have complete SAE mapping of LLMs such as Gemma-4 and plausible labels for each of these features.
The gist is that, in order to understand the meaning of a feature, it is necessary to identify the least common denominator across all the scenarios in which that specific feature fires. So the first step is to come up with a hypothetical scenario in which we expect a certain feature to fire.
Based on the way we humans interpolate cubic polynomials, I made the hypothesis that the NN would recognize landmark points of a polynomial such as roots, critical points (roots of the derivative) and inflection point (where the concavity of the curve changes, or root of the second order derivative).
To understand whether an NN recognizes landmark points, I draw, for a polynomial parametrized by , the activation pattern over for all the 64 features. I also draw vertical dotted lines in proximity of the landmark points: in particular, red indicates the roots, blue indicates the critical points and green the inflection point. The figure shows, for the polynomial parametrized by and , the activation pattern of the 64 features and, right below, the actual shape of the polynomial.
The picture is pretty complicated to understand. We can try to make it easier by isolating the activation patterns of individual features such as feature #12 or #18.
Looking at the activation pattern of feature #12 I see somewhat a pattern of bumps in proximity to critical points, and a fall in proximity of the second root. Let's analyze the activation pattern of the same feature for a different polynomial.
Looking at the activation pattern of the same feature over various cubic polynomials, the correlation seems more difficult to infer, although we can still see the dip in proximity of the second root.
After having spent too much time staring at these plots, I tried to approach feature discovery in a more structured and scalable way.
For example, let's say that we want to measure the correlation between feature #12 and the roots for the curve . To do that we calculate the Pearson correlation score between the activation pattern of that feature (see figure above) and an artificially made curve with gaussian-like bumps in proximity of the roots. The result is a correlation score of , which makes sense looking at the previous diagram: the feature does indeed dip in proximity of the root, but it is not localized to that region only. Now let's try to analyze the behavior of feature #12 over a broader range of cubic polynomials by sliding over a range of combinations . The heat-map should be read as follows: the pixel indicates the magnitude of the correlation between the activation profile of feature 12 and the roots of the cubic polynomial parametrized by .
Let's try to find out other "interesting" features. What are the features that have the highest correlation with the roots?
We find some interesting symmetric pattern over the axis for features #4, #12, #18 and #47. We call these features "wings in the clouds" (thanks Paula for the suggestion!). We next analyze features that highly correlate with critical points. Black pixels indicate polynomials for which no critical points exist.
Feature #57 shows a positive correlation with the critical points . In particular, it seems to fire in proximity for cubic polynomials with , and it is invariant to .
Features #4, #8 and #12, identified before as highly correlated with roots, record high correlation to critical points (and inflection point too). These are sort of "generalist" features and therefore not very interesting to us.
Lastly, we spot various features (#61, #48, #56, #55, #14, #51) that show correlation (between 0.4 and 0.2) with inflection points for very localized portions of polynomials.
At this point I keep staring at these graphs to figure out what they really mean. The "wings in the clouds" feature looks bizarre to me, but I don't really know how to interpret it. One possibility is that my initial hypotheses were simply wrong and our NN had formed abstractions to interpolate polynomials that are intractable by humans and therefore cannot be discovered by any human-generated hypothesis.
So, what is the tool at my disposal that is closer to the NN alien mind? AIs are amazing at pattern matching at scale. So I provided Claude Code with a large set of SAE vector activations for thousands of combinations and tasked it with finding "interesting" patterns across the activations.
So, what did Claude discover? First of all, let's have a look at the lens that Claude used to find patterns. The approach used in the previous section is highly supervised: it starts with a hypothesis "features might track landmark points" and tries to find a correlation.
Claude, instead, adopted a more exploratory lens. It started by fixing b and plotting the raw activation value itself directly over the plane. Without any a-priori hypothesis to test. This approach allows the raw geometry of the activation to emerge. For , here's the activation profile of the 64 features. Notice that the pixels of the heat-map have a different meaning: each pixel measures the activation of a feature for a given combination , where is fixed.
What immediately emerges is the diagonal band behavior shown by features #7, #16, #18, #22, #33, #47, #53, #54 and #64. There's clearly a pattern here. We can also see that most of the features are all blue so mostly inactive. We notice that these features are mostly b-dependent and are likely to fire for different values of b. We now provide a representation of the diagonal band features only.
How to decipher these heat-maps? It seems that each of those nine features activates for a combination of x and a. In other words, the feature does not seem to care separately about the value of x or about the value of a. They only care about a specific linear combination of them that can be expressed, generically, as
What's even more impressive is that the behavior of these diagonal band features is invariant to the parameter . Let's reproduce the same plot by fixing . We see an identical activation pattern.
Ok, so we have identified nine diagonal band features that activate for a specific linear combination and are invariant to . They also appear to cover the whole spectrum of the possible cubic polynomials in the search space for .
The next step to better describe these features is to identify the slope of the diagonal band. In particular, we run the following test for each feature:
We run this test for each of the nine diagonal band features and for candidate values in the range .
The table above displays, for each feature, the best (and the related score). We see that the median of the best candidate across the nine features is 0.3. This suggests that we can somehow approximate the linear combination that dictates the activations of these features as . This variable is surprisingly close to the substitution operated in the first step of the Cardano method.
The Cardano method is a popular trick to solve cubic polynomials of type (yeah, I know that this is not the common convention for naming coefficients but I wanted to be consistent with the notation used throughout the essay).
The method was introduced to the world by the publication of Ars Magna by Gerolamo Cardano in 1545, although the attribution of the invention is quite disputed.
The first step of the method is to perform a variable replacement inside the cubic polynomial in order to remove the quadratic term. Let's do it step by step.
What you are left with, after the variable substitution, is a depressed cubic polynomial: a cubic polynomial that is missing its quadratic term. Which makes it easier to solve. You can follow along with this article to understand the second step of the Cardano method and find the roots.
We can also explain this substitution trick from a geometrical point of view. The transformation is a horizontal shift of the graph along the x-axis. More precisely, this transformation moves the graph such that its inflection point, located at the x-coordinate sits precisely on the y-axis.
After the transformation, the curve is symmetric over the y-axis which makes it somehow easier to reason about.
Ok, so let's go back to our NN and try to summarize the discovery. We have an NN that is trained to regress cubic polynomials given an input combination to compute . We observe that our NN, at its fourth and last layer, developed nine different features that activate for a specific linear combination .
Similarly, if a high-school student is presented with the task of solving a cubic polynomial , they will start by using the Cardano method and apply the substitution .
So what can we claim on the basis of these two observations? The NN has discovered a new coordinate and it activates it to regress a cubic polynomial. Similarly, humans seem to call the same concept any time they need to find the solution for a polynomial of the same type. This is an amusing finding!
There are a few nuances here. The slope of the linear transformation is not quite the same. For the NN the median slope is , while for the Cardano method it is . Additionally, the tasks are not exactly the same: the Cardano formula has been developed to solve a cubic polynomial, while the NN is only required to regress it. Lastly, we identified the existence of a feature, but have not tested yet how this is being used.
We registered the occurrence of multiple diagonal band features. It appears that the NN calls these features when it comes to compute a new coordinate , as part of the broader task of regressing cubic polynomials. We also observed, surprisingly, that a similar coordinate is computed as part of the first step of the Cardano method to solve cubic polynomials in order to remove the quadratic term and obtain a depressed cubic polynomial, which is easier to solve.
So the real question now is: how are these features even being used ? What is the specific role of these features? How do they influence the downstream behavior of the model?
Let's take a step back and describe how one can perform intervention over an NN to verify if their interpretation of a feature is correct. Starting from the SAE vector , it is possible to artificially manipulate it by pushing up (or down) some specific feature in order to affect the reconstructed activation . The figure displays a legit reconstruction of the last activation layer (top) and a manipulated reconstruction (bottom) after clamping features #7 and #18 by a factor of 1.5x and suppressing feature #44.
To measure the impact of the feature clamping, the forward pass is run on a model in which the real activation is replaced by the activation derived after manipulation.
A good example of model manipulation comes from Anthropic. The researchers identified, via SAE, a specific feature in Sonnet that looks like firing in context related to the Golden Gate Bridge. In order to verify their interpretation, they tuned the activation of that feature artificially high and observed the behavior of the model. The model started mentioning the Golden Gate Bridge in every response, demonstrating that the feature was labelled correctly.
In the case of my experiment, it is somewhat harder to predict what is going to happen when manipulating the diagonal band features. The only thing that I can think about is that the diagonal band features are associated with a sort of geometrical transformation. Let's see what happens if we try to regress a given polynomial () while manipulating the diagonal band features via various factors .
We can try to do the same with another cubic polynomial .
It seems like the diagonal band features are somewhat associated with a multiplicative effect or vertical stretch of the curve and the magnitude of the transformation is proportional to the multiplicative factor applied to manipulate such features.
We notice that the error is not a consequence of the feature manipulation but it is attributable to the "standard" SAE reconstruction error. You can see that, where the multiplicative factor is 1 (no manipulation was applied), the error is still present.
Does this confirm or refute our Cardano-based interpretation of the diagonal-band features? Honestly, I don't know. If the feature directly mapped to a horizontal shift transformation by , that would have been a confirmation of the interpretation. What we observe is that these features cause a vertical stretch of the curve. This result leaves the door open for further testing.
One further experiment that I plan to run is to train a SAE on top of the middle layers of the NN, which are often associated with abstract thinking.
It's finally time to re-abuse the word abstraction. Abstractions are compressed (and lossy) representations of reality. Such compression is necessary to let a computationally bounded observer deal with the over complexity of the world. In the introduction I hinted at the abstractions (or heuristics, if you prefer) that we humans use when it comes to regressing polynomials: by staring at an equation of type we can derive information such as the asymptotic behavior of the curve and its shape.
For the same task, we considered the point of view of an NN. The bounded space they have at their disposal, measured by their parameters, implies that they must have developed some compressed way to perform the task. Or, in other words, they must have developed some sort of abstract representation. Throughout the essay we referred to these abstract representations as features and we were able to extract them using SAEs. Surprisingly, the NN seems to have developed at least one abstraction that is similar to one that humans use when performing the same task.
This finding makes me wonder how inevitable are the abstractions that we humans developed. Spelled differently, how wide is the ruliology slice of abstractions that can be used to regress cubic polynomials? My feeling is that it is very large, and the occurrence of a shared abstraction between us and the NN is purely coincidental. I'm not even sure that we have access to the same space of abstractions. If that's the case, any effort to interpret an NN will always be incomplete.
While fully interpreting the abstractions developed by an NN might be an overreaching task, I still feel that NNs can be used as a playground to study how abstractions emerge. As an example, one can study how abstractions emerge as a function of time. For NNs, time is measured by training rounds. Extracting the learned features at different training epochs (and analyzing their impact over the output accuracy) can tell us a lot about how abstractions emerge. Without necessarily having to understand them. If I had to guess, I suspect that NNs' first instinct is to operate at the lowest possible level of abstraction (which corresponds to the highest level of complexity) and soon discover that this is too much to handle given their space boundedness. Training will then push the learning in two directions:
Analogously, it would be cool to study how abstractions emerge as a function of space (or computational boundedness) of the observer. In the case of an NN, space is determined by the number of parameters at its disposal. What I highly suspect is that bigger models should have more room to operate at a lower abstraction level and therefore feel less pressured to move up in the abstraction ladder. One early result that I didn't share in this essay already makes a point. Consider the following 7x7 grid of cubic polynomials.
For each polynomial, I trained an NN of size 15x15x15x15 to predict its regression. Compared to the NN trained in this essay, each of the 49 trained NNs only takes x as input. In other words, each model has to regress a single cubic polynomial. Given a specific NN (in the case of the following figure, the chosen model is the one approximating the cubic polynomial parametrized by and ), the SAE feature analysis revealed that the features are highly correlated with the x input range. As if the model "simply" learned to associate a behavior based on the magnitude of the input. Since each model had to deal with a single polynomial, there is not enough training pressure to develop any higher-level representation. For example, ones related to the geometry of the curve. As if the models were lazy.
Another interesting question is to ask how the emergence of abstractions is path-dependent. How does the discovery of abstraction X impact the chance of abstraction Y to emerge? Can you jump from abstraction level to level or is it necessary to pass through level ?
The long term goal is to develop a model that describes how abstractions emerge, using NNs as a proxy for the human brain.