Exploring Hidden Semantics in Neural Networks with Symbolic Regression

Many recent studies focus on developing mechanisms to explain the black-box behaviors of neural networks (NNs). However, little work has been done to extract the potential hidden semantics (mathematical representation) of a neural network. A succinct and explicit mathematical representation of a NN model could improve the understanding and interpretation of its behaviors. To address this need, we propose a novel symbolic regression method for neural works (called SRNet) to discover the mathematical expressions of a NN. SRNet creates a Cartesian genetic programming (NNCGP) to represent the hidden semantics of a single layer in a NN. It then leverages a multi-chromosome NNCGP to represent hidden semantics of all layers of the NN. The method uses a (1+$\lambda$) evolutionary strategy (called MNNCGP-ES) to extract the final mathematical expressions of all layers in the NN. Experiments on 12 symbolic regression benchmarks and 5 classification benchmarks show that SRNet not only can reveal the complex relationships between each layer of a NN but also can extract the mathematical representation of the whole NN. Compared with LIME and MAPLE, SRNet has higher interpolation accuracy and trends to approximate the real model on the practical dataset.


INTRODUCTION
Neural networks (NNs) have been successfully applied in many problems, such as CNN for object recognition [26], RNN for time series analysis [11], and Bert for natural language processing (NLP) [5]. However, neural networks are often seen as black-boxes because their input-output (IO) relationships are difficult for a human to understand [12]. Sometimes, it is almost impossible to interpret the NN behaviours when models make unexpected predictions on some datasets, such as adversarial examples [31] and white noise * Corresponding author.
images [22]. Some of the recent work has been centred on researching and explaining the black box behaviours, as summarized in several survey papers [1,29].
In this paper, we aim to develop a new symbolic regression-based method to explore hidden semantics in NNs. Here, a typical hidden semantics interpretation method refers to explaining a NN with an explicit math function. If a function ( ( ) ) explains a single hidden layer ℎ( ( ) ) on a certain input-output ( ) , ( ) or a small range of inputs-outputs, it is called local explanation [35], such as LIME [27] and MAPLE [24]. If a function ( ) is able to explain a hidden layer ℎ( ) on the whole dataset, it is called global explanation, such as Visualization method [3,17] and Net2vec [10]. Although these methods show some degree of success in extracting hidden semantics from NN, they have the following three limitations. (1) Local explanation methods can give a mathematical expression, such as a linear model and a decision tree, for each ( ( ) , ( ) ). However, they cannot obtain a general expression for the whole dataset. Although most global explanation methods can visualize NNs on the whole dataset, they cannot give a mathematical expression for explaining the dataset. (2) These local or global methods often leverage pre-defined interpretable models to explain the hidden semantics of NNs. For example, LIME can use linear models, decision trees, or falling rule lists as interpretable models. However, these pre-defined models may not capture hidden semantics in some situations because the real characteristics of a NN model are often unknown, and applying a predefined model to explain the networks may be inappropriate. (3) These methods cannot generate a mathematical expression that can represent the hidden semantics of all layers in a NN.
To overcome these limitations, this paper leverages the symbolic regression (SR) method to explain a NN. In SR, for a given dataset { , }, the algorithm can find a symbolic function ( ) = ′ that minimizes the distance between and ′ in the mathematical expression space. SR has great flexibility in generating mathematical expressions; hence, it does not need a predefined model to capture the relationships in the dataset. However, classical SR methods, such as GP [14,30], GEP [8], and linear GP [2], usually handle the symbolic function ( ) with a single output ′ , i.e., ′ is a number and not a vector. They cannot represent the relationship ( ℎ −1 + ) of each layer in a NN because each layer's output is a vector, matrix, or tensor. Therefore, when these GPs explain NN [6,9], they can only give a mathematical expression to show the semantics of the whole NN with a single output value, not each layer in the NN. Although Cartesian Genetic Programming (CGP) [21] supports multiple outputs, CGP cannot represent semantics in a NN. Because each CGP output corresponds to a hidden node, and it cannot provide a general model that represents hidden semantics in the layer. To obtain a general model, we assume that the relationship between input and output in a layer (or a NN) has the mathematical expression format (ℎ −1 ) + , where and may be a number, vector, matrix, or tensor, and is a mathematical function that represents the hidden semantics of the layer. Based on the above assumption, this paper proposes a novel SR method (called SRNet) to mine hidden semantics of all layers in a NN simultaneously, as shown in Figure 1. SRNet is an evolutionary computing algorithm. In each evolution, SRNet first leverages the Cartesian Genetic Programming (CGP) [19,20] to find each layer's mathematical function (ℎ −1 ). It then uses the Newton-Raphson method [28] (or L-BFGS method [16]) for few (or many) variables to obtain and so that ℎ = (ℎ −1 ) + approximates the output ℎ of the layer in a NN. At the end of the evolution, SRNet will capture hidden semantics of all layers in a NN when ℎ ≈ ℎ (including ≈ ). The main contributions in the paper are summarized as follows: • The paper proposes a new method called SRNet to explain hidden semantics of all layers in a NN. SRNet generates a mathematical expression in the format of (ℎ −1 ) + that can be used to explain a NN.
• To speed up SRNet, we create a multi-chromosome CGP [33] evolutionary strategy embedded in the Newton-Raphson method. • Experiments show that the proposed SRNet can capture hidden semantics in NN in 12 SR benchmarks and 5 classification benchmarks. Compared with LIME and MAPLE, SRNet has higher interpolation accuracy and trends to approximate the real model on the practical dataset.
The remainder of this paper is organized as follows. In Section 2, we introduce the background knowledge about Cartesian Genetic Programming. Then, we propose SRNet to explore hidden semantics in NNs in Section 3. Section 4 and 5 report the experimental results. We conclude the paper in Section 6.

CARTESIAN GENETIC PROGRAMMING
CGP is a directed acyclic graph-based genetic programming algorithm for addressing the SR problem [21]. In CGP, the graph consists of a two-dimensional grid of computational nodes, as shown in Figure 2. These nodes are classified into three categories: input, output, and function. The input (or output) nodes represent the input (or output) values (or ), which is encoded into an integer, such as 0 encoded by "1" and encoded by "4". The function nodes are computational expressions. Each function node has three parts, input, computational expression, and output, encoded by a series of integers. As each node has only one output, the output code is used to index the node. For example, a function node "+" is regarded as the code "⟨0012⟩". In this code, the first integer 0 is the code of the function "+". The two middle integers "0" and "1" represent two inputs of the function "+", and they are also the outputs of two previous nodes. The last integer "2" is the index of the node. In the example CGP diagram, there is no edge between any two nodes in the same column. Two nodes at different columns can be linked if one's input code equals the other's output code. A genotype is used to represent a CGP, as shown in Figure 2. The genotype contains two categories of nodes: the functional nodes and the output nodes. As the genotype has multiple output nodes, the genotype can describe multiple computational expressions. For example, the genotype in Figure 2 generates three mathematical expressions, = − 1 , = 2 0 1 + 2 1 , and 2 0 + 1 . CGP usually leverages the (1 + ) evolutionary strategy [13,25] to find the best fitted mathematical expression. In each evolution, (1 + ) EA utilizes mutation to generate offsprings. For CGP, the mutation randomly chooses a gene location and changes the allele at the location to another valid random value. A valid value is from the function look-up table if a computational expression gene is chosen for mutation. If an input gene is chosen, a valid value is from the output set of its previous nodes. The mutation does not change the output gene. For example, in Figure 2, a mutation changes the input gene "1" of the node "7" to "3". Then, the output becomes "0".
The multi-chromosome Cartesian genetic programming (MCGP) [33] encodes multiple chromosomes into a single genotype. Each chromosome code is similar to the genotype of CGP. So, MCGP can provide a solution to a large problem by dividing it into many smaller sub-problems. For simultaneously exploring hidden semantics of all layers in a NN, MCGP encodes each layer semantics as a chromosome into a genotype. Thus, a chromosome represents the semantics of one layer, and the genotype represents all NN layers' semantics. After MCGP uses the (1 + ) multi-chromosome evolutionary strategy [34] to acquire a best-fitted individual, it also obtains these semantics.

SRNET
This section proposes the SRNet, a method based on MCGP [33] to simultaneously explain the hidden semantics of layers in a NN. As shown in Figure 1, SRNet can find a group of ℎ = (ℎ −1 ) + that approximates each NN layer output ℎ , i.e., To find the group {ℎ } quickly, SRNet needs to address the following problems: 1) how to encode these ℎ s, and 2) how to find functions to explain hidden semantics in a NN. Section 3.1 shows our solution to the first problem, and section 3.2 provides a multi-NNCGP evolutionary strategy to address the second problem.

SRNet Encoding
The hidden semantics ℎ of each layer in a NN are represented as where (ℎ −1 ) is a mathematical expression that represents the general semantics in the layer , and and are a weight vector, matrix, or a tensor, as shown in Figure 3. To capture a NN layer's semantics ℎ , we define that the neural network CGP (NNCGP) consists of three components, including general semantic model, constant, and operator. The general semantic model is used to generate (ℎ −1 ). It has three parts, inputs, × functions and an output. 2) The constant component includes the weight vector and the bias vector , where | | = | | = |ℎ |.
SRNet uses multiple NNCGPs (called MNNCGP) to encode genotype. Since a NNCGP can represent one layer's semantic, multiple NNCGPs can be used to capture hidden semantics of all layers in a NN. These NNCGPs are regarded as chromosomes that constitute a genotype. In the genotype, each NNCGP 's output ℎ is the input of its next NNCGP +1 .

Evolution Strategy
SRNet leverages a multi-NNCGP evolutionary strategy embedded by the Newton-Raphson method (called MNNCGP-ES) to find the best-fitted genotype that represents hidden semantics of all layers in a NN. MNNCGP-ES is similar to the (1 + ) multi-chromosome evolutionary strategy [34]. MNNCGP-ES includes the following operations: mutation, fitness evaluation, and selection. Mutation, with a certain probability, change each allele in MMCGP to another valid random value [21].

Fitness Evaluation.
To evaluate a genotype encoded by MN-NCGP, the fitness function is defined as the following equation, where L is the mean squared error (MSE) of each middle layer. L is an error function of the output layer. It is a cross-entropy loss in the classification task, while it is MSE in the regression task. ℎ (ℎ ) is the output of the th layer in a NN (NNCGP ). and are the outputs of the NN and SRNet, respectively. Equation 2 indicates that obtaining ℎ needs two computations, as shown in Figure 3. One is (ℎ −1 ) that generates an output by the CGP code. The other is the parameter computation that obtains the constant vectors, and , by the Newton-Raphson method that performs an update operation according to the following equation.
where is or , ( ) is Hessian, ∇ ( ) is the gradient of the loss function ℎ − ( (ℎ −1 ) + ). The Hessian is difficult to obtain if it is a high-dimensional matrix (i.e., many neurons). Therefore, to solve this problem, the Limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm (L-BFGS) [16] is used for limited memory and time-saving.

Selection.
Selection aims at generating a 'super' individual from the population. The 'super' individual consists of a set of super chromosomes that have the best fitness at each position from all individuals, i.e., L (ℎ , ℎ ), as shown in Figure 4. Since these chromosomes constitute an input-output sequence where each chromosome output is the input of its next chromosome, they need to be selected in the order of their positions in the population. After evaluating the finesses of the chromosomes at a certain position ( 0 ) of the population, the selection picks up the chromosome that has the best fitness score as the super chromosome ( 0 ). Then, it evaluates chromosome at the next position( 1 ) with 0 as input. Moreover, it obtains a chromosome with the best fitness as the next super chromosome ( 1 ). Repeating the above evaluation policies results in a group of selected chromosomes { 0 , 1 , ...}. These selected chromosomes form a super individual.
However, for high-dimensional problems, evaluating the finesses of the chromosomes is very time-consuming when using the Newton-Raphson method. So, L-BFGS is used to replace the Newton-Raphson method to compute the two weight vectors ( and ) of all individuals. L-BFGS can speed up the fitness evaluation.

MNNCGP-ES.
The MNNCGP-ES pseudocode is listed in Algorithm 1. MNNCGP-ES combines the fitness evaluation and the selection method mentioned before, using the (1 + ) evolution strategy to evolve generation-by-generation to obtain the optimal individual.

EXPERIMENTS
To validate the SRNet's ability to explain hidden semantics in the neural network (NN), we tested the SRNet on the built NNs of 12 symbolic regression benchmarks as well as 5 classification benchmarks, listed in Table 1. Moreover, the sample sizes and feature sizes of the 17 benchmarks are illustrated in Figure 5.

Alias
Function/Name Training Dataset MLP 0   Table 2 lists the parameters of the three algorithms, LIME, MAPLE, and SRNet. All parameters of the three algorithms are fixed on all benchmarks.

Regression Task
To validate the SRNet's ability to explore hidden semantics of NN in the regression task, we chose 12 symbolic regression benchmarks 0 − 5 and 0 − 5. The benchmarks 0 − 5 were chosen from the commonly used SR Benchmarks [18], while 0 − 5 are from physical laws [32]. We generated 12 datasets (called true datasets) according to these benchmarks. Each of these datasets has different sample sizes (see 'Training Dataset' in Table 1). For example, for the 1 problem in Table 1, we randomly sampled 200 and values in the range of [−1, 1], respectively. Moreover, combining these values can generate 200 samples for the 1 dataset. For each of the 12 datasets, we randomly took 80% samples from it as the training dataset, and the other as the test dataset. Then, we built 12 Multi-Layer Perception neural networks (MLP) with a sigmoid activation function using the training datasets. The training parameters are listed in Table 1. For example, for 1, we created an MLP with two hidden layers where each layer had three hidden nodes. The MLP was trained by the SGD optimization method with a learning rate of 0.01.
After each MLP was trained, we collected each NN layer's input and output data of the 12 MLPs as the NN explanation datasets. We then ran the MNNCGP-ES, LIME, and MAPLE 30 times on each NN explanation dataset.

Classification Task
To validate the SRNet's ability to explore hidden semantics of NN in the classification task, we chose 5 classification benchmarks named 0 − 4 from the PMLB [23] with the different number of samples and features (see Figure 5). The column 'Training Dataset' on the rows "P0-P4" indicates the number of samples, as shown in Table  1. Training 5 MLPs is similar to the regression task except for the function "softmax" that replaces their output functions. After training these classification MLPs, we need to compare SRNet with their decision boundaries, not their outputs. Because the MLP's outputs on the training datasets are sparse and do not fully represent the NN's classification ability, as shown in Figure  6(a). Using these outputs to train SRNet may result in wrong results. To obtain the decision boundaries of the trained MLP, we leverage a uniform sample method around its decision boundary (called USDB). USDB first evaluates the range of the training dataset. It then randomly samples points in the range. It finally selects points with the shortest distance to the decision boundary according to Equation 5 [7,15].
, where is a sample, is the decision boundary of a NN, and is the number of sample classifications. is the probability that the belongs to the th classification, which is the NN output owing to its activation function "softmax".
After each classification MLP was trained, we utilized USDB to generate samples around the decision boundary of the MLP. We then fed these samples into the MLP and collected each NN layer's input and output as the MLP explanation dataset. We finally ran the MNNCGP-ES, LIME, and MAPLE 30 times on the explanation dataset.

RESULT AND ANALYSIS 5.1 Regression Task
In the following regression tasks we only show results on the six benchmarks, 0 − 3, 0, and 1. The regression results on all benchmarks are in the appendix.  mathematical expressions that approximate the hidden semantics of each layer in NNs according to Equation 3, as shown in Table  3. It also indicates the feasibility to use the combination of these mathematical expressions to represent the whole semantics of each MLP, as shown in Table 4.   Not all ranges of fitness scores (light blue areas) become smaller as the MNNCGP-ES runs, such as 1. However, the low bounds of these lines always become smaller and trend to be zero at the later stage. It means that the more episodes the MNNCGP-ES runs, the more likely MNNCGP-ES is able to find the mathematical expressions that can be used to explain the hidden semantics of a NN. The slow decrease of fitness curves also indicates the need to run MNNCGP-ES with sufficient times to obtain the best-fitted mathematical expressions.

Semantics Evaluation.
The results show that the proposed SRNet method can acquire a fitted mathematical expression to explain hidden semantics of each layer, as shown in Table 3. Each of these mathematical expressions represents the general semantics in the expression (ℎ −1 ) + (Equation 2). The number below each is its fitness. For example, for 1 , MNNCGP-ES finds the general semantics −(ℎ 0 ) 0 + (ℎ 0 ) 2 at the second hidden layer ℎ 1 . The     fitness of × (−(ℎ 0 ) 0 + (ℎ 0 ) 2 ) + is 1.39 − 03. All fitness values in the hidden layers are less than 0.1. It means that captured by MNNCGP-ES can represent (approximate) the semantics of each hidden layer in a NN. Figure 8 shows the details of mathematical expressions that MNNCGP-ES finds. The expressions approximate the output of each layer in a NN. To evaluate each layer, we input 9 random values into the mathematical expression and the hidden nodes in the NN, respectively. We then obtained 9 heat maps to show the difference between the output of the mathematical expression and that of the NN hidden nodes. For example, the sub-figure " 1 − ℎ1" represents the 1 outputs of SRNet vs NN in the layer ℎ1 with 9 random input values. In the first heat map in " 1−ℎ1", [0.6, 0.6, 0.3] are the outputs of three hidden nodes in the NN layer ℎ1 with the first input value, while [0.9, 0.9, 0.3] are the output of the mathematical expression in the SRNet layer ℎ1. So, Figure 8 explains why a mathematical expression has low fitness, and another has high fitness. Comparing 1−ℎ1 with 0−ℎ1 in Figure 8, the outputs between NN and SRNet in 0−ℎ1 are closer than the output between them in 1−ℎ1. Thus, the fitness "7.54 − 05" of "0.88 − (ℎ 0 ) 0 " is less than "1.39 − 03" of " −(ℎ 0 ) 0 + (ℎ 0 ) 2 " in Table 3.
For the output layers in the NNs, the last column in Table  3 lists the mathematical expressions to present the outputs. Their fitness scores of the output layer are better than the scores of the previous layers. There are two formulas of 5 and 5 whose fitness scores are greater than 1. The reason causing the higher scores at the output layers is that, in a NN, the network structure of the outer layer is different from the hidden layers. The computation function on the output layer is = +1 ℎ + +1 , while that on the hidden layer is ℎ = ℎ −1 + . Since is a number and ℎ is a vector, = +1 ℎ + +1 is a multivariate linear equation. However, = +1 +1 (ℎ )+ +1 on the output layer in SRNet is one-variable linear equation because +1 and +1 are two numbers, not vectors. Although MNNCGP-ES can find a fitted mathematical function +1 to represent the NN layers, the one-variable linear equation is still hard to approximate the multi-variable linear equation in the layers. Table 4 lists the mathematical expressions that represent the whole NN semantics for the regression tasks. Each of final expressions is obtained by combining the mathematical expressions in different layers shown in Table 3. The complexity and length of the mathematical expressions could increase substantially, such as the mathematical expression in 3, due to the combination of expressions on every layer. If there are more layers in a NN, the length of the mathematical expression could be longer. Although the math representation generated by SRNet could be lengthy, it provides a straightforward expression to show all layers' hidden semantics of the whole NN.

Performance Comparison.
To evaluate the SRNet performance on regression tasks, we ran and compared SRNet, LIME [27], and MAPLE [24] on an interpolation dataset and an extrapolation dataset. The interpolation dataset consists of the NN input-output values. In contrast, the extrapolation dataset consists of the data sampled directly from the original symbolic expression in the column "Function" in Table 1. The interpolation domain is the same as the range of the training dataset shown in Table 1. The size of the extrapolation domain is five times that of the interpolation domain. Figure 9 illustrates the curves (or distribution points) of the true dataset, as well as the results of NN(MLP)s, SRNet, LIME, and MAPLE, on different symbolic regression benchmarks. For the high dimension datasets, it is not easy to visualize the curves. So, the curves are projected into samples on multiple planes. The curves between two vertical blue lines represent interpolated results, while those outside the two lines are the extrapolated results. The interpolated results show that SRNet can find the mathematical expressions close to MLP on most of these benchmarks. SRNet can find smoother results than LIME, while it can find results closer to MLP compared with MAPLE. The extrapolated results show that LIME can find the model closest to MLP because LIME is the local explanation method that generates a model for local (several) samples. For a extrapolate dataset, LIME divides the datasets into many groups of local samples and generates a model for each group of samples. Therefore, it needs many local models to explain the extrapolated dataset and cannot provide a general model to describe the whole dataset. Unlike LIME, SRNet finds the mathematical expression that represents the whole dataset. In addition, it can trend towards the true dataset generated by the symbolic regression benchmark. Therefore, SRNet has better extrapolation results in our evaluation [4]. The ability of SRNet extrapolation could help SRNet judge why a NN fails to predict some test data. Under the assumption that the mathematical expression found by SRNet can represent the real model on the practical dataset. For example, on the benchmark 1, if a NN is unable to predict certain output, the prediction result could be compared with the output of −0.01 1 + 1.69 sin ( 0 ) + 0.0021 found by SRNet. The difference between the NN output and SRNet can be used to analyze the gap between the NN and the real model.

Classification Task
In the following classification tasks we only show the results on the two benchmarks, 0, and 1. The classification results on all classification benchmarks are included in the appendix. Figure 10 shows the fitness convergences curves of the SRNet for the classification tasks on the two benchmarks, 0 and 2. SRNet converges rapidly, especially within about 100th generations. As MNNCGP-ES runs L-BFGS to obtain weight vectors ( and ) of all individuals ( (ℎ −1 ) + ) every 50 generations, MNNCGP-ES only runs L-BFGS two times, and it can converge to an accurate result. The reason is that SRNet does not need a whole dataset to be trained, but a thousand samples around the decision boundary of a NN. For example, on the benchmark 0 , although it has 48842 points, these points are only used to train a classification NN. After training the NN, USDB ( Section 4.2) randomly samples 1000 points around the classification NN. Then, the 1000 points are used to train SRNet. So, the process of training SRNet in the classification task is fast.  Table 5 lists a mathematical expression of each NN layer obtained by SRNet for the two classification tasks, 0 (adult) and 2 (agaricus_lepiota). Moreover, Table 6 shows the final mathematical expressions obtained by SRNet. For the adult dataset, the two mathematical expressions, 0 and 1, represent the prediction probability of SRNet for class 0 (adult makes over $50K a year) and class 1 (adult makes below $50k a year), respectively. Moreover, 0 and 1 indicate that the classification results only depends on the three features, 6 (relationship), 7 (race), and 9 (capital-gain). In addition, they provide a mathematical explanation of how the trained NN classifies data. For the agaricus_lepiota dataset, SRNet gives a simple linear model as the explanation of the classification NN, as shown in 2 in Table 6. According to 0 and 1, the classification NN mainly focused on the four features, 10 (stalk-root), 20 (population), 8 (gill-color) and 9 (stalk-shape). In this case, SRNet degenerates into a simple linear model. So, when SRNet explains a NN on the classification task, it not only shows how the NN computes on the dataset, but also represents which features (variables) the NN focuses on.   Table 7 lists their average (+ std) prediction accuracy on the train and test dataset. Interestingly, SRNet is better than LIME and MAPLE on most classification tasks and only fails on the 'agaricus_lepiota' ( 2) task. The fail reason is that the dataset 'agaricus_lepiota' has 22 input features, making    1000 points sampled by USDB very sparse in the high-dimensional space. However, SRNet still shows better or competitive accuracy in the other four classification tasks than LIME and MAPLE. Although LIME and MAPLE have good performance in explaining local classification samples, their decision boundaries cannot approximate the MLP's decision boundary. The reason is that they only leverage a linear model to explain several local samples, as shown in red lines and purple lines in Figure 11. So, once the MLP's decision boundary is a complex curve at some local samples, the line generated by LIME or MAPLE cannot approximate it. In contrast, SRNet can approximate the complex decision boundary of MLP on all samples. Therefore, SRNet is more suitable for explaining NN on the classification task than LIME and MAPLE.

CONCLUSION
This paper proposes a new evolutionary algorithm called SRNet to address the NN's black box problem. SRNet leverages MNNCGP-ES to find the mathematical expressions that can represent each NN layer's hidden semantics. The combination of every layer's expression represents the whole NN. Compared with the models found by LIME and MAPLE, the mathematical expression provided by SRNet is closer to the NNs in the interpolated domain. The experiment also shows the SRNet models trend to approximate the real data model that used trains NN. The close alignment of SRNet with the real model and its explicit mathematical expression can be used to facilitate the explanation of NN prediction behaviours, such as regression and classification.