Ex-Ante Assessment of Discrimination in Dataset

Data owners face increasing liability for how the use of their data could harm under-priviliged communities. Stakeholders would like to identify the characteristics of data that lead to algorithms being biased against any particular demographic groups, for example, defined by their race, gender, age, and/or religion. Specifically, we are interested in identifying subsets of the feature space where the ground truth response function from features to observed outcomes differs across demographic groups. To this end, we propose FORESEE, a FORESt of decision trEEs algorithm, which generates a score that captures how likely an individual's response varies with sensitive attributes. Empirically, we find that our approach allows us to identify the individuals who are most likely to be misclassified by several classifiers, including Random Forest, Logistic Regression, Support Vector Machine, and k-Nearest Neighbors. The advantage of our approach is that it allows stakeholders to characterize risky samples that may contribute to discrimination, as well as, use the FORESEE to estimate the risk of upcoming samples.

• We provide experimental proofs of how our measure of ex-ante risk is a useful warning sign for the demographic disparity, inequality of odds and opportunities of classifiers trained on a sample of the data.
• We show how model averaging over diverse decision trees provides reliable estimates for our ex-ante risk.

Problem Setting
We frame the task where a set of individuals sampled from X × Y (where x ∈ X describes individual's features and y ∈ Y the ground truth of target variable) can lead into finding g : X → Y used to support decision-making process. We assume that g has a sufficient performance such that the decision-makers consider the function g as part of the process. For example, g can be the classifier for mortgage applications, where upcoming applicants' information are the inputs and low and high risk of mortgage delinquency is the outcome. Additionally, we assume that for each sampled individual we are given their sensitive attributes s ∈ S, which represents protected information forbidden to be used in the classification. Furthermore, sensitive attribute is defined by the union of the different groups, i.e., S = ∪ j S i ∀i ∈ [j], with j the total of different groups and S i is a group of individuals with same sensitive attribute(s) usually represented by a number. Although S is not part of the inputs, they are allowed to be used as the evaluation of discrimination across the groups in the sensitive attributes.
The objective is to determine how likely any individual (x, s, y) ∈ D will lead into discrimination a function g : X → Y learned from D s.t. g show sufficient performance to be used by decision-makers. This objective implies two challenges. First, we have to define what discrimination stands for. And second, we are aimed to determine the discrimination when g is unknown, i.e., under ex-ante notion. We discuss and address these two aspects in the following section under a binary classification task, although this is extensible to other task types.

Risk Scores and Aggregate Misclassification Rates
Discrimination can be understood as the disparity between two groups of individuals. In the mortgage application example, disparities can be found along the process, i.e., if there is explicit evidence of disparity through the use of sensitive attributes as part of the classification, and there is sign of dependency between the outcome and the protected attribute, it can be determined that there exists disparate treatment [21]. One solution is to remove any sensitive attribute, which is known as Fairness Through Unawareness (FTU) [22]. However, FTU is not enough since the correlation of sensitive attribute with other features may still end up with disparities in the outcome, which is known as disparate impact [21]. Based on this fairness notion, Demographic Parity (demP ) measures disparity across groups over the outcome of a function g [21]. Similarly, individual fairness certificates that g maps similar outcomes for similar individuals [23]. Using these notions, we define an ex-ante measure by evaluating disparity on the labeled dataset (instead of the outcome of g) for similar individuals of a given sample. More in detail, we state that fairness risk score is the conditional probability of Y being 1 (from demP ) given X = x with differences in the sensitive attribute S (from individual fairness).
In the following formulation, we first discuss theoretically how the risk can be obtained from infinite population (or at least large enough). However, since this is not true for most of the cases, we then discuss two estimation proposals under finite population in Section 3. Definition 1. For any individual (x, s) ∈ X × S, its fairness risk score r(x) is defined as We claim that Definition 1 is ex-ante since it does not depend on any classifier that will use the data to make decisions. However, this is based only on independence notions of disparate treatments, which can lead into flaws of identification of discrimination. According to Kizilcec and Lee [7], separate notion can solve this issue by including the ground truth. We argue that Definition 1, as ex-ante risk, captures this granular and aggregate fairness properties of classifiers using the data: (i) differences in individual misclassification rates across demographic groups; (ii) differences in aggregate misclassification rates across demographic groups. Moreover, we show that ex-ante risk score is directly related to the differences δ mis (x) in misclassification rate conditional on X = x for an unknown function g, where δ mis is defined as Equation 2 .
Having this we state the Theorem 1, which implies that the ex-ante risk score captures how the misclassification rate of any unaware classifier vary across demographic groups. Consequently, for a given X = x, a large risk score means that individuals in one group will be more likely to be misclassified than in the other, regardless of the classifier. Theorem 1. For any unaware classifier g, x ∈ X , we have: are the characteristic functions of the sets {x ∈ X |g(x) = 1} and {x ∈ X |g(x) = 0}. Moreover, for a deterministic unaware classifier, δ mis (x) = r(x).
Proof. Let η s (x) denote P (Y |X = x, S = s). For any classifier g, we have where I {g(x)=1} is the indicator function equal to one whenever g(x) = 1; and, (a) uses the fact that conditional on X, g(X) is independent of Y . Therefore, If g is deterministic, we have either g(x) = 0 or g(x) = 1, which leads to δ mis (x) = r(x).
Our ex-ante fairness risk score also captures aggregate fairness properties of any classifier g, which allows us to look at the aggregate difference in misclasification rates [15] defined as Equation 3.
Large values for ∆ mis means that the error rate is not equalized across groups. Having this, we state Theorem 2, which implies that if the distribution of features X does not depend on the sensitive attribute, then the total variation between distribution of X conditional on S = s and S = s is zero and low aggregate risk scores implies low ∆ mis (g). Therefore, whenever the distribution shift between X|S = s and X|S = s is small, our measure of risk allows users to identify data that are unlikely to lead to unfair classification, when fairness is measured in terms of difference in misclassification rates (i.e., in terms of separation notions [7]). Theorem 2. For any unaware classifier g, where T V (P s , P ns ) is the total variation between the distribution P s of X conditional on S = s and the distribution P ns of X conditional on S = s.
where (a) uses the fact that c(x) ≤ 1 and η s , η ns ≤ 1, the triangular inequality. The result in Theorem 2 follows from the definition of total variation.

Risk Scores and Sub-population Fairness Metrics
Results in Theorem 2 do not provide us with a lower bound on aggregate fairness metrics, even in the absence of a distributional shift between P s and P ns . In this section, we show a lower bound for misclassification rates at the sub-population level. Specifically, for γ ∈ (0, 1), a sub-population level fairness metric measures the maximum value this metric takes over any subset G of the feature spaces X such that P (G) ≥ γ (e.g. [24,25,26]). Sub-population level definitions of algorithmic fairness are stronger than their aggregate counterparts, since they protect subgroups defined by a complex intersection of many sensitive attributes ( [25,27]) or a structured slicing of the feature space [24]. We define sub-population differences in misclassification rates of a classifier g as where ∆ mis (g|x ∈ G) is the difference in misclassification rates between demographic groups in sub-population G. Theorem 3. Suppose that there exists a sub-population G ⊂ X such that P (G) > γ, P s = P ns , and P (Y = 1|X = x, S = s) > P (Y = 1|X = x, S = s). For any stochastic classifier g : Theorem 3 implies that whenever the response function for a sub-population G, a classifier that predicts Y = 1 on this sub-population will lead to higher misclassification rates for one demographic group and the differences in misclassification rates will be bounded from below by the average risk score. Therefore, by looking at sub-populations with high predicted risk scores, the user can identify sub-populations where differences in sensitive attributes S will correlate with differences in misclassification rates.
Proof. Let 2κ = inf x∈G g(x) − 1/2 > 0. Using derivations from the proof of Theorem 1 and the assumption that P 0 = P 1 , we can show that where (a) uses that η s (x) ≥ η ns (x) on G.

Risk Estimation
To estimate the ex-ante risk score, we propose two methods: (i) a Bayesian-based approach to estimate the probability P (Y = 1|X = x, S = s) followed by the computation of the risk score r(x); (ii) an ensemble approach that averages the risk estimates from diverse decision trees.

Method 1: Risk Score of Average Models
Our first approach consists of first obtaining an estimate of the response function E[Y |X = x, S = s] and E[Y |X = x, S = s] and then, compute the risk score. We follow [28] and use an additive tree model BART ( [20] to obtain a non parametric estimate of the response functions. Given M trees T 1 , T 2 , ..., T M , we compute the risk score as BART uses a prior to regularize the depth of each tree T m . This non-parametric approach has been successful to estimate individual treatment effect following the seminal work of [28].

Method 2: Average of Risk Scores
One limitation of a BART estimate of the ex-ante risk score is that it requires to estimate both response functions BART averages decision trees to obtain a robust Bayesian estimate of each response function.
In this section, we present an algorithm that directly averages the risk score estimates from each decision tree.
Given a collection of trees T 1 , T 2 , ..., T M , we define as L 1 (x), L 2 (x), ..., L M (x) the leafs to which an individual with feature x ∈ X belongs to. For each decision tree T m , we compute its estimate r(x, T m ) of the risk score as the difference between demographic groups in the estimation error rate Y mis in L m (x): Our risk estimate is then Intuitively, decision trees partition the feature space into regions where outcomes are assumed to be constant. Leaf level risk scores characterize how violations of this assumption vary across demographic groups.
To obtain a robust risk estimate, we average a diverse set of trees obtained as followed: (i) we train each tree on a different random sub-sample of the insances and a random subset of the features; (ii) we retain the trees that achieve a minimum performance β ∈ (0, 1).
Step (i) allows to ensemble trees that see different aspect of the data generating process and de-correlates estimation errors.
Step (ii) allows to retain only the trees that captures at least partially the characteristics of the response function E[Y |X = x]. It is worthy to note that we use S in the training. We argue that for a leaf with none members of one group would reflect dependency between Y and S given X = x, hence a high risk of discrimination. Under an optimistic scenario, the expected misclassification rate should be equal across groups in the leaf, while the highest error rate is expected for the not present group under pessimistic one. Since we aim to flag probable discrimination, we argue that the latter is preferable. Note that a better approach might be used, however, according to the results in experiments section, our assumption seems to be sufficient. Finally, more details about the FORESEE algorithm are provided in the appendix.

Experimental Evaluation
We design our experiments to answer the four following research questions.
• RQ1: Do risk estimates from FORESEE and BART average to the ground truth risk in numerical simulations?
• RQ2: Is the ex-ante risk score a good predictor of how unfair classifiers using the data as input will be? • RQ3: Does our measure of risk allow the user to describe the sub-populations that are the most likely to be under-served by classifiers using the data as inputs? • RQ4: How can ex-ante risk score guide the development of fair algorithms?

Data sets
Synthetic Data First, we rely on numerical simulations to test whether FORESEE and BART are on unbiased estimator of the true risk score. We uniformly draw features from [0, 1] 2 and assign a sensitive attribute s with P (S = 1) = 0.5. For S = 1, we assign a label Y = 1. For S = 0, we assign a label Y = 1 with probability equal to 1 − x1+x2 2 . Therefore, for a given X = x, the ground truth risk is equal to r(x) = x1+x2 2 . We draw 5000 samples from this synthetic data, run FORESEE and BART to estimates r F ORESEE and r BART for each sample and compare their values to the ground truth risk r(x). We repeat the protocol 20 times with different seeds to compute mean and standard deviation of the risk estimates.
Real world datasets We use two benchmark data sets in fair machine learning and a new dataset in education data mining. The former two are publicly available at ProPublica webpage and UCI Machine Learning repository.
Adults (a.k.a. Census Income Data Set) consists of 48, 844 individuals described by 14 features that include marital status, education, working hours. The sensitive attribute is the gender to which individuals self-identify to. The outcome is whether an individual's income is larger than 50K [29].
Compas gathers information on 7, 214 individuals in the criminal justice system to assess their recidivism risk. Features include number of felonies, misdemeanors and degree of charges. The sensitive attribute is whether an individual self-identifies as African American. The outcome is whether an individual re-commits a crime within two years [2].
Dropout contains academic records of 4, 706 students in a Latin-American university. Features include grades and participation to academic support programs. The sensitive attribute is the gender to which each student self-identifies to. The outcome reflects drop out within the first two years.
We split each dataset into training and testing sets using a ratio of 70/30.

Risk Scores vs. Downstream Fairness Metrics
To explore how ex-ante risk score is informative of future unfair uses of the data, we train four models -Logistic Regression (LR), Random Forest (RF), k-Nearest Neighbors (KNN), and Support Vector Machine (SVM) -to predict Y from X. We compute accuracy for evaluating models' performance in Compas, and F-1 score in Adults and Dropout. We also compute three standard fairness metrics: demographic disparity δ demP [22]; inequality of opportunity δ opp [30]; and, inequality of odds δ odd [30]. We also classify individuals into High and Low Risk whose risk score are larger and smaller than λ, respectively. However, this threshold is truly context dependent. We leave it to the stakeholders to decide which threshold is appropriate in a specific application.

Risk Score and Mitigation
To show how our ex-ante risk score can guide machine learning developers, we explore two risk-based fairness mitigation strategies. We identify all instances with risk score larger than λ and apply two mitigation approaches: (i) pre-processing and (ii) post-processing.
Preprocessing removes the high-risk instances from the training set and only report classification outcomes on instances with low ex-ante risk score. This is a very relevant approach whenever stakeholder can replace automated decisions by human ones for outcomes with high ex-ante risk score. On the other hand, post-processing consists of adjusting the classification threshold over the classifier's predicted class conditional probabilities, while minimizing one of the aggregate fairness metrics and the error rates [30]. Figure 1 shows the mean and standard deviation of r F ORESEE (x) (left) and r BART (x) (right) against the ground truth risk r(x) for each x in our synthetic dataset. We observe that FORESEE generates a statistically less unbiased estimate of the ground truth risk since the mean of the r F ORESEE (x) across 20 simulations is equal to r(x) for almost all values of r(x) ∈ (0, 1). On the other hand, BART estimates have an upward bias for low values of the true risk score r(x) and a downward bias for large values r(x). This observation confirms our intuition that it is preferable to directly average risk estimates from local approximations instead of first estimating the response functions over the entire dataset and then computing the risk score.  Table 1 shows the classification performance of four models across our three real world datasets. Best classifiers in terms of F1 score/accuracy vary across datasets: Random Forest achieves the highest performance on Adults and Compas, while Logistic Regression on Dropout. Moreover, best classifiers generate significant fairness metrics, where of particular interest is whether our ex-ante risk score anticipates these fairness issues.

Model Adult Dropout
Compas  Table 1: Performances of machine learning models on each data set. Higher F-1 score and accuracy (acc) are better, whereas for Equalized Opportunity (δ opp ), Equalized Odd (δ odd ), and Demographic Parity (δ demP ), lower is better.  Across all datasets there exist instances with high risk scores, which indicates that misclassification rates would vary with changes in sensitive attribute. This is particularly striking for Compas, where whole of the population is estimated to have an ex-ante risk score larger than ≈ 0.4. High risk scores for Compas confirms findings in existing studies (e.g [2]) that this dataset is not well suited for training automated decision systems without significantly harming a demographic group.
Using Figure 2, we define as high risk a sample with risk higher than 0.5 for Adults and Compas, and higher than 0.06 for Dropout. In Table 2, we find that all three fairness metrics are higher for high risk individuals than low risk ones. This is evidence that our ex-ante risk score allows identifying individuals for which a given classifier's aggregate performances would significantly vary across demographic groups. Notice that the differences are not significant for Compas as the other datasets, which is consistent with the distribution of our risk core estimates depicted in Figure 2 and flags once again the issues of using this set for automated training [31].

Fairness Adult Dropout Compas Metric
High | Low High | Low High | Low  Table 2: Fairness metrics in absolute value for groups defined by level of ex-ante fairness risk score. The thresholds between Low and High risk groups are set to 0.5 for Adult and Compas, and 0.06 for Dropout. We truncate figures to three decimal places. Figure 3 shows the differences between the High and Low Risk profiles in the datasets. The high risk group is composed of the top 20% of samples with high risk scores drawn equally from both groups in S. The low risk group is similarly formed with the lowest 20% risk scores. From figure 3 we can characterizes high risks as in the following list (the low profile are the opposite to this description):

RQ3: Profiles for High and Low Risk Instances
• Adult: aged individuals, high number of working hours per week, and married or divorced.
• Dropout: students with low academic performances and higher numbers of academic warnings. • Compas: younger defendants, fewer criminal history, mostly felony as charge degree, and spending less days in jail.
Our risk score provides stakeholders and data owners information on the characteristics of individuals for which predictions would be affected by their sensitive attributes. A useful implication of the profiles in Figure 3 is to use with caution any data mining algorithm that predicts the target variable if the individual' risk estimate is high.

RQ4: Risk Score Guided-Mitigation Strategies.
We test the pre-and post-processing strategies previously described. Table 3 shows that a pre-processing strategy that removes samples at high risk of unfair outcomes is effective at reducing aggregate inequality of odds and opportunities. This result is evidence that our ex-ante risk score is instrumental in identifying the data instances that push the classifier toward unfair outcomes. Moreover, our pre-processing strategy even improves the performance, which are consistent with previous studies that have shown results in favor of pre-processing approaches [32,33]. On the other hand, removing high risk individuals in Compas dataset results in insignificant improvement, which is expected given the found distribution of risk scores. Intuitively, although we remove instances with scores > 0.5, the remaining still has high levels (mostly > 0.4). Consequently, we find that post-processing approaches on the risk sub-sample might be better a approach, where post-processing subjected to demographic parity is the best one. Table 3 also shows that post-processing tailored to high risk group is effective at reducing each fairness metric, although at a cost of performance. Therefore, developers can use the high risk group as a unit test against which they can automatically evaluate the fairness of their algorithm; then, post-process classification outcomes for that group.

Conclusions
Although the fair machine learning literature offers an abundance of metrics to measure the potential harmful effect on individuals of a machine learning model, stakeholders would like to anticipate these fairness pitfalls before investing in model development. We propose an ex-ante risk score to measure whether future performances of any classifier using the data will vary with sensitive attributes. Moreover, we propose a bagging approach (FORESEE) to estimate the risk score from a finite sample. Our experiments demonstrate that our risk score, although measured before model development, anticipates correctly whether a model using the data will under-serve some demographic groups. Moreover, our risk score offers useful guidance for stakeholders to decide whether a data is appropriate for training automated decision systems; for developers to mitigate potentially discriminatory outcomes; and for end-users to decide when to trust the outcome of the algorithm. However, our approach is still computational expensive. Also, although a pessimistic scenario assumption may over-estimate the misclassification rate variation, our results show this seems to be sufficient.


extfieldsofstudy:
Computer Science
Engineering
Mathematics

s2fieldsofstudy:
Computer Science
Business

text:
Multiple Signal Classification Based Joint Communication and Sensing System

Joint communication and sensing (JCS) has become a promising technology for mobile networks because of its higher spectrum and energy efficiency. Up to now, the prevalent fast Fourier transform (FFT)-based sensing method for mobile JCS networks is on-grid based, and the grid interval determines the resolution. Because the mobile network usually has limited consecutive OFDM symbols in a downlink (DL) time slot, the sensing accuracy is restricted by the limited resolution, especially for velocity estimation. In this paper, we propose a multiple signal classification (MUSIC)-based JCS system that can achieve higher sensing accuracy for the angle of arrival, range, and velocity estimation, compared with the traditional FFT-based JCS method. We further propose a JCS channel state information (CSI) enhancement method by leveraging the JCS sensing results. Finally, we derive a theoretical lower bound for sensing mean square error (MSE) by using perturbation analysis. Simulation results show that in terms of the sensing MSE performance, the proposed MUSIC-based JCS outperforms the FFT-based one by more than 20 dB. Moreover, the bit error rate (BER) of communication demodulation using the proposed JCS CSI enhancement method is significantly reduced compared with communication using the originally estimated CSI.

FD capability, as detailed in [10]. Therefore, the self-leakage between arrays is ignored and not considered in the signal model in this paper. One BS array is used for transmitting the DL JCS signal, and the other is used for consistently receiving echoes of the JCS signal. MUE receives the JCS signal to demodulate the communication data, while BS receives the echoes to estimate the AoAs, ranges, and velocities. Moreover, we consider that both BS and MUE receive the superimposed co-channel interference from multiple reflected interference sources (ISs). The MUE is equipped with one UPA for receiving the communication signal. The array sizes of the BS and MUEs are P t × Q t and P r × Q r , respectively. where ϕ k is the azimuth angle, and θ k is the elevation angle. We use A p,q to denote the (p,q)th antenna element, and A 0,0 to represent the reference antenna element. Then, the phase difference between A p,q and A 0,0 is expressed as a p,q (p) = exp −j 2π λ d a (p cos ϕ sin θ+q sin ϕ sin θ) ,

B. UPA Model
where λ = c/f c is the wavelength of the carrier, c is the speed of light in vacuum, and f c is the carrier frequency.
The steering matrix for K far-field signals is then represented as which is a matrix of dimension P Q × K.

C. DL JCS Signal and Channel Model
In this paper, we consider the JCS system using OFDM-based signals. The transmitting signal is s D (t) = where P t is the DL transmit power, d n,m is the mth baseband OFDM symbol of the nth subcarrier, ∆f is the subcarrier interval, T = T s + T g , T s = 1 ∆f is the duration of OFDM symbol, T g is the guard interval, M s and N c are the number of OFDM symbols and subcarriers, respectively, and Rect (t/T ) is the rectangular window function of duration T . When the DL preamble signal for beam alignment and CSI estimation is transmitted, d n,m is replaced by the preamble symbols, denoted byd n,m , which is known and deterministic to both BS and MUE. When the DL data signal is transmitted, d n,m ∈ Θ QAM is a random symbol, where Θ QAM is the constellation of quadrature amplitude modulation (QAM). Note that d n,m is known to BS but unknown to MUE.
Next, we present the JCS channel model. As illustrated in Fig. 1, the DL JCS channel comprises a communication channel and an echo sensing channel.
• The JCS communication channel consists of a line-of-sight (LoS) path and several nonline-of-sight (NLoS) scattering paths.
• The JCS echo sensing channel consists of the echo path from MUE as a scatterer, and the echo paths from other scatterers which may or may not contribute to the communication channel. Since the signals after multiple reflections are much smaller than those with only one reflection, we only consider echoes directly reflected from scatterers.
Then, the JCS sensing echo and communication channels at the nth subcarrier of the mth OFDM symbol are defined as [5], [8] where p T X,l is the AoD of BS's JCS transceivers, a(p T X,l ) ∈ C PtQt×1 is the corresponding transmit steering vectors as given in (2), L is the number of scatterers, l = 0 is for the direct path between BS and MUE, l = 1, · · · , L − 1 is for the reflected paths involved the lth scatterer.
Moreover, i = S and i = C represent the echo sensing and communication channels, respectively; p S RX,l and p C RX,l are the AoAs of BS's echo receiver and the MUE's communication receiver, respectively; and α S,n,m,l and α C,n,m,l are the channel fading for the lth sensing echo path and communication path, respectively. 1) JCS Echo Sensing Channel: When i = S, a(p S RX,l ) ∈ C PtQt×1 is the receive steering vector for the lth echo sensing path, as given in (2). Since the mmWave array is typically small, p S RX,l = p T X,l . Moreover, α S,n,m,l is the fading factor for the lth echo (when l = 0, MUE acts as a scatterer), which is given by where f s,0,1 = 2v r,0,1 λ and τ s,0 = 2d 0,1 c are the echo Doppler frequency shifts and time delay between BS and MUE, with v r,0,1 and d 0,1 being the corresponding radial relative velocity and the distance, respectively; f s,l,1 = 2v r,l,1 λ and τ s,l = 2d l,1 c are the echo Doppler frequency shifts and time delay between BS and the lth scatterer, with v r,l,1 and d l,1 being the corresponding radial relative velocity and distance, respectively. Moreover, b S,l = λ 2 (4π) 3 d l,1 4 β S,l , and β S,l is the random reflection fading factor of the lth scatterer, following the complex Gaussian distribution with zero mean and variance σ 2 Sβ,l . 2) JCS Communication Channel: When i = C, a(p C RX,l ) ∈ C PrQr×1 is the receive steering vector for the lth communication path, as given in (2). Moreover, α C,n,m,l is the fading factor for the lth path, and is expressed as where f c,d,0 = v r,0    and time delay between BS and scatterer, and between the scatterer and MUE of the lth NLoS path, respectively, with v r,l,2 and d l,2 being the radial relative velocity and distance between the lth scatterer and MUE, respectively; b C,0 = λ 2 (4πd 0 ) 2 is the propagation loss of the LoS path, and b C,l = λ 2 (4π) 3 d l,1 2 d l,2 2 β C,l is the path fading factor of the lth NLoS path with β C,l being the scattering factor of the lth scatterer. Here, β C,l is the random reflecting factor of the scatterer in the lth path, which is assumed to follow the complex Gaussian distribution with zero mean and variance σ 2 Cβ,l . Due to the existence of b C,l , the LoS path is much stronger than the NLoS path for mmWave.
Note that H C,n,m is unknown and needs to be estimated by utilizing the DL preambles,d n,m .
The parameters of H S,n,m are unknown, and BS has to estimate the AoA, range and Doppler in H S,n,m . Since BS acts as both the sensing transmitter and receiver, both d n,m andd n,m can be used for DL sensing. Moreover, l = 0 represents a special path, for which the echo time delay and Doppler are twice of those in the communication channel, which is the theoretical basis for the JCS CSI enhancement method to be introduced in Section III-B.

D. JCS Received Signal Model
In this subsection, we present the expressions for DL JCS received signals.

1) DL Communication Received Signal:
The frequency-domain DL communication signal received by MUE at the mth OFDM symbol of the nth subcarrier is expressed as y C,n,m = P t d n,m (w RX ) H H C,n,m w T X + n X C,n,m , where w T X ∈ C PtQt×1 and w RX ∈ C PrQr×1 are the JCS transmit and communication receive beamforming (BF) vectors, respectively; w T X 2 2 = w RX 2 2 = 1. In this paper, the lowcomplexity least-square (LS) method is used to generate w T X and w RX for BF. BS utilizes the known DL preambles, i.e., d n,m =d n,m , to conduct beam alignment with MUE. When beam alignment is completed, w T X = c 1 [a T (p T X,0 )] † and w RX = c 2 [a(p RX,0 )] † , where c 1 and c 2 are both arbitrary complex values with modulus 1, [·] † is the pseudo-inverse operation,p T X,0 ≈ p T X,0 , andp RX,0 ≈ p C RX,0 . Simultaneously, the unknown communication CSI, (w RX ) H H C,n,m w T X , can be estimated by processing the received preambles.
Moreover, n X C,n,m = n C,n,m +ξ C,n,m is the sum of noise and interference, n C,n,m = (w RX ) H n C,n,m and ξ C,n,m = (w RX ) H x C,n,m are transformed noise and interference, the dimensions of n C,n,m and x C,n,m are both P r Q r × 1, n C,n,m is Gaussian noise vectors with each element following , and x C,n,m is the reflected interference signals from other network devices. We assume there are N ic ISs, and the reflected fading for each IS follows a Gaussian distribution.
Since the superimposed one of multiple random OFDM signals is noise-like, the pth element of x C,n,m can be given as where P i,c is the power of incident signal from the ith IS, and β I i,p ∼ CN (0, 1).
. Further, we define the communication SINR (C-SINR) as where h C,n,m = (w RX ) H H C,n,m w T X is the gain of DL communication signal at each antenna element.
2) DL Echo Sensing Received Signal: The echo signal that BS receives for the mth OFDM symbol at the nth subcarrier is given by y S,n,m = P t d n,m H S,n,m w T X +n X S,n,m = P t d n,m L−1 l=0 (α S,n,m,l )χ T X,l a(p S RX,l ) +n X S,n,m , (10) where χ T X,l = a T (p T X,l )w T X represents the gain of the DL JCS transmit BF, n X S,n,m = n S,n,m + x S,n,m is the sum of noise and interference, n S,n,m is the Gaussian noise vector with each element following CN (0, σ 2 N ), x S,n,m is the superimposed interference vector for N is reflected ISs, and the dimensions of n S,n,m and x S,n,m are P t Q t × 1. Similar to x C,n,m , the pth element of x S,n,m can be given as where h S,n,m,l = α S,n,m,l χ T X,l is the gain of DL echo sensing signal at each antenna element.
By defining s n,m,l = √ P t d n,m α S,n,m,l χ T X,l and s n,m = [s n,m,l ]| l=0,1,...,L−1 , (10) can be expressed in the matrix form as where A S,RX = [a(p S RX,l )]| l=0,1,...,L−1 is the steering matrix stacked by steering vectors of L echoes, A S,RX ∈ C PtQt×L , and s n,m ∈ C L×1 . By stacking all the M s OFDM symbols with N c subcarriers, we have

III. DL JCS SIGNAL PROCESSING
In this section, we demonstrate the signal processing for DL JCS sensing and communication, which is shown in Fig. 3. We first present the sensing signal processing scheme, and then elaborate on the JCAS CSI enhancement method.

A. JCS Sensing Signal Processing
In this subsection, we first present the conventional MUSIC method for estimating the 2D AoAs, and then introduce the novel MUSIC-based range and Doppler estimation method.

1) JCS MUSIC 2D Angle Detection:
First, the correlation matrix of Y S is obtained as By applying eigenvalue decomposition to R X , we have where Σ x is the real-value eigenvalue diagonal matrix in descending order, and U x is the orthogonal eigen matrix. Calculate the average of eigenvalues and denote it as m x . Let α t be a preset threshold, which is determined as elaborated in Appendix A. Then, the number of echo paths is determined as the number of eigenvalues no smaller than α t m x , denoted by N x .
Construct U N = U x (:, N x + 1 : P t Q t ) 1 as the noise subspace basis. We then use it to obtain the spatial angular spectrum function as [12] f a (p; where p = (ϕ, θ) is the 2D angle, and a (p) is given in (2). The spatial spectrum is represented The maximum points of S a (p; U N ), i.e., the minimum points of f a (p; U N ) are the estimated AoAs [13]. We first find N x local maximum points of S a (p; U N ) using a grid searching method with relatively large granularity, then we use the Newton descent method to identify the accurate minimum point of f a (p; U N ) by inputting the above local maximum points as initial points for iteration.
2) JCS Range and Doppler Detection: After the AoAs are obtained, through BF at the AoA of interest, the filtered received signal at the nth subcarrier of the mth OFDM symbol can be expressed as where w t,n,m,k =(w RX,S,k ) H n X S,n,m is the transformed noise and interference with zero mean and variance σ W 2 , RX,l,k = (w RX,S,k ) H a(p S RX,l ) is the receive BF gain, and w RX,S,k is the receive BF vector for the kth AoA, k ∈ [0, 1, ..., N x − 1]. Note that RX,k,k is typically larger than RX,l,k (l = k) due to the narrow beam feature of mmWave. By substituting (6) into (18), we obtain (19).
In (19), there are independent complex exponential functions for range and Doppler, i.e., e −j2πn∆f ( r l c ) and e j2πf s,l,1 mT , respectively. Here, we define the range and Doppler steering vectors as respectively. The range and Doppler steering matrices are defined as where A r ∈ C Nc×L , and A f ∈ C Ms×L . Stackȳ S,n,m,k into a matrixȲ S where Ȳ S n,m =ȳ S,n,m,k , then erase the communication where the division is element-wise, andH S ∈ C Nc×Ms .
According to (19),H S can be expressed by A r as where On the other hand, the transpose ofH S , i.e., H S T , can be presented by A f as where The range and Doppler can be estimated via the autocorrelation ofH S and H S T , which are given by respectively. Denote the noise subspaces of R X,r and R X,f as U x,rN and U x,f N , respectively.
Theorem 1. The minimum of U x,rN H a r (r) 2 2 , denoted by r s,l , is linked to the range via r s,l = 2d l,1 . The minimum of U x,f N H a f (f ) 2 2 corresponds to the Doppler value, f s,l,1 .
Proof. The proof is presented in Appendix C.
By applying eigenvalue decomposition to R X,r and R X,f , we have where Σ x,r and Σ x,f are the real-value diagonal matrices of eigenvalues in the descending order, and U x,r and U x,f are the corresponding eigenvector matrices.
We use m x,r to denote the mean value of Σ x,r , and then set the threshold α t,r using the method in Appendix A by replacing Σ x with Σ x,r . The number of targets in the AoA of interest, N x,r , is then determined as the number of eigenvalues no smaller than α t,r m x,r . Then, the noise subspace basis for range estimation is derived as U x,rN = U x,r (:, N x,r + 1 : N c ). Since the number of targets is the same for both Doppler and range estimation, the noise subspace basis for the Doppler estimation can be derived as U x,f N = U x,f (:, N x,r + 1 : M s ).
We use U x,rN and U x,f N to derive the range and Doppler spectrum functions as respectively. The range and Doppler spectra can be given by respectively.
The maximum points of S r (r; U x,rN ) and S f (f ; U x,f N ), i.e., the minimum points of f r (r; U x,rN ) and f f (f ; U x,f N ), are the range and Doppler estimation values, denoted byr s,l andf s,l , respectively. The distance, d l,1 , and radial velocity, v r,0,1 , between BS and the target are given bŷ The minimum of f r (r) and f f (f ) can be identified using a two-step Newton descent method.
We first find the local maximum points of S r (r) and S f (f ) with large-granularity grid searching.
Then, we use the Newton descent method to find the accurate minimum points of f r (r) and f f (f ) using the above local maximum points as the initial points. The iterative expression for the Newton descent method is derived as follows.
By applying the Taylor series decomposition to f r (r) and f f (f ), and taking their first order derivative over r and f , respectively, we obtain ∂f r (r) ∂r and By setting the above first-order derivative to be 0, the iterative descent expression for range and Doppler estimation can be given by respectively. From (29), the first-order and second-order derivatives of f f (f ) and f r (r) are expressed as ∂f r (r) ∂r ∂f f (f ) ∂f =2Re{a (1) where a (1) r (r), and a

B. JCS Communication Signal Processing
By substituting (5) into (8), and taking into consideration that w T X and w RX in (8) where h C,n,m = b C,0 RX,0 χ T X,0 e j2πmTsf c,d,0 e −j2πn∆f τ c,0 is the real communication channel response, and RX,0 = (w RX ) H a(p C RX,0 ) and χ T X,0 = a T (p T X,0 )w T X are the BF transmitting and receiving gains. In the CSI estimation, d n,m =d n,m , and we denote y C,n,m =ȳ C,n,m as the received signal. The CSI estimated with the LS method is expressed as [14] h C,n,m =ȳ where w C,n,m = Output: Filtered sequence [h C,n,m ]| n=0,··· ,Nc−1 .
Step 2: respectively. Based on the prior information obtained by JCS sensing, we propose a Kalman filter-based JCS CSI enhancement method to improve CSI by leveraging the sensing estimation results of JCS.
For the mth OFDM symbol,ĥ C,n,m can be regarded as the observation of h C,n,m as given in (41). Since b C,0 is unchanged for the same OFDM symbol. The state transfer of h C,n,m is given by The Kalman filter algorithm that utilizesΦ = [ĥ C,n,m ]| n=0,··· ,Nc−1 to recursively derive the estimation of Φ = [h C,n,m ]| n=0,··· ,Nc−1 is presented in Algorithm 1, with the details of the Kalman Filter algorithm referenced to [15]. Note that we obtain h C,n,m = e −j2πn∆f τ c,0 h C,0,m from (42), based on which we can further estimate the initial observation variance as After

IV. PERFORMANCE ANALYSIS OF THE JCS PROCESSING
In this section, the analytical MSE results of AoAs, range, Doppler, and location estimation of the proposed MUSIC-based JCS processing are derived using the perturbation method.

A. Analysis of 2D Angle Detection MSE
From (13), the noise term N t can be treated as the perturbation to the useful signal, which is expressed as where Y S,R = A S,RX S is the useful signal. The singular value decomposition of Y S,R can be expressed as where U 0 is the noise subspace basis, and U 0 H Y S,R = 0. Further, we have U 0 H A S,RX = 0.
With noise as perturbation, Y S can be expressed as whereΣ 0 = ∆Σ 0 ,Σ s = Σ s + ∆Σ s , andŨ 0 = U 0 + ∆U 0 . Here,Ũ 0 andŨ s are both orthogonal unitary matrices, andŨ H 0 Y S =∆Σ 0Ṽ H 0 . In the high SINR regime, solving the perturbation problem is equivalent to seeking the optimal ∆U 0 to minimize Ũ H 0 Y S 2 subject to the constraintŨ H 0Ũ 0 = I [13]. By substituting (44) The second-order perturbation (∆U 0 ) H N t and U 0 H Y S,R =0 can be discarded. By using the LS method [13], ∆U 0 can be presented as The MUSIC 2D angle estimation result is distorted by the noise perturbation, which is expressed asp k = p k + ∆p k , where p k is the actual value of AoA. Apply Taylor series decomposition to f a (p k ;Ũ 0 ) in (16), and take the first three terms. Applying first-order derivative to the truncated Taylor series, we have By setting (49) to be 0, we can obtain where H p (p k ;Ũ 0 ) = ∂ 2 fa(p k ;Ũ 0 ) With the perturbation expression, we can obtain From (16), we have where vec(·) is to vectorize a matrix, a p (p) and a (2) p (p) are the first-order and second-order derivatives of a(p) over p, respectively, which can be derived from (2).
Since U 0 H A S,RX = 0, we can obtain We use H p0 to represent H p (p k ; U 0 ). By substituting (55) and (56) into (51) and (52), (50) can be rewritten as Discarding the perturbation terms that are higher than second-order in (57), we can rewrite (57) as where the perturbation expression of ∆G p is derived in Appendix D as By substituting (59) and (56) into (58), we can obtain ∆p k as shown in (60).
The MSE of angle estimation can be expressed as  (25), is rewritten asH whereH S,p = A r S r,s is the useful signal. The singular value decomposition ofH S,p is where U r,s and U r,0 are orthogonal unitary matrices, and U r,0 HH S,p = 0. We further obtain U r,0 H A r = 0.
By treating W tr as a perturbation term,H S can be decomposed as

C. Analysis of Location MSE
The location of the target can be obtained after the AoA, p k = (ϕ k , θ k ), and the range, r k , are detected. The expression for the actual location is given by p loc (r k , ϕ k , θ k ) = (x k , y k , z k ) = (r k sin θ k cos ϕ k , r k sin θ k sin ϕ k , r k cos θ k ).
With the AoA and range estimation perturbation, ∆p k = (∆ϕ k , ∆θ k ) and ∆r k , the location of the target is Comparing (81) with (82) and discarding the second-order perturbation, we can represent the perturbation of x, y, and z axes coordinates as ∆x .
Finally, the location error can be expressed as

D. Cramer-Rao bound of JCS Sensing
We further derive the Cramer-Rao bound (CRB) to characterize the minimum lower bound for sensing. Based on the signal model presented in Section II-D2, the echo signal of the lth target received by the (p, q)th antenna element at the nth subcarrier of the mth OFDM symbol is y p,q S,n,m,l = P t d n,m α S,n,m,l χ T X,l a p,q (p l ) + n p,q S,n,m + x p,q S,n,m , where α S,n,m,l = b S,l e j2πmTs2v s,l /λ e −j2πn∆f 2r s,l /c is given as (6), v s,l = v r,l,1 and r s,l = d l,1 are the radial relative velocity and distance between BS and the lth target, respectively; a p,q (p l ) is given in (1), p l = p T X,l = (ϕ l , θ l ) is the 2D AoA of the lth target; χ T X,l is the transmitting BF gain; n p,q S,n,m and x p,q S,n,m are the noise and interference at the (p, q)th antenna element. Let n X,p,q S,n,m n p,q S,n,m + x p,q S,n,m , then n X,p,q S,n,m is independent and identically distributed, following CN (0, σ 2 W ), where σ 2 W = P IS + σ 2 N . Let ψ = (r s,l , v s,l , ϕ l , θ l ) be the set of estimation parameters. Then, the distribution of y p,q S,n,m,l is p(y; ψ)= 1 πσ 2 W e − y− √ Ptdn,mα S,n,m,l χ T X,l ap,q(p l ) 2 2 /σ 2 Because there are N c M s P t Q t independent symbols used for estimation, the joint distribution of these symbols is (n,m,p,q) y p,q n,m −sn,mα S,n,m,l ap,q( where ρ = ( 1 πσ 2 W ) NcMsPtQt , and s n,m = √ P t d n,m χ T X,l . Note that d n,m is independent and identically distributed with E( d n,m 2 2 ) = 1. According to [16], [17], the CRB of ψ i , ψ i ∈ (r s,l , v s,l , ϕ l , θ l ), is given by  With (89) and (90), the sensing CRBs are derived as (q cos ϕ l sin θ l −p sin ϕ l sin θ l ) 2 , C θ l = λ 2 8π 2 d 2 a γ S NcMs p,q (p cos ϕ l cos θ l +q sin ϕ l cos θ l ) 2 , where γ S is the S-SINR as given in (11).

E. Complexity Analysis and Comparison
In this section, we analyze and compare the complexity of the proposed MUSIC-based JCS method with the conventional FFT-based methods. We consider three schemes: Scheme 1 is the proposed MUSIC-based method; Scheme 2 is the original FFT-based method in [6]; and Scheme 3 is the Code-division OFDM (CD-OFDM) FFT-based method in [8]. It can be seen that Scheme 2 has the lowest complexity. The complexity of Scheme 3 increases due to the additional code-division multiplex processing. The complexity of our proposed MUSIC-based JCS method has the highest complexity to achieve super-resolution detection.

V. NUMERICAL AND SIMULATION RESULTS
In this section, we present extensive simulation results for the proposed MUSIC-based JCS processing method, with comparison to the Schemes 2 and 3 as described in Section IV-E, and verify them against the analytical performance bounds derived in Section IV. We also compare the BER results of communication demodulation for the proposed JCS CSI enhancement method with those in conventional communication systems.

A. System Setup
The system setup largely follows the specification in the 3GPP Vehicles-to-Everything (V2X) applications [18]. The carrier frequency is 63 GHz, the antenna interval, d a , is half of the wavelength, the sizes of antenna arrays of BS and MUE are P t ×Q t = 8×8 and P r ×Q r = 1×1,

B. Sensing Performance
We first demonstrate the sensing spectra of schemes 1, 2, and 3. The normalized range spectrum and radial velocity spectrum are shown in Figs. 4(a) and 4(b), respectively. The S-SINR is γ S,n,m = −20 dB. For range estimation as shown in Fig. 4(a), the peak to sidelobe ratio (PSLR) of scheme 1 is about 26 dB. By contrast, the PSLRs of schemes 2 and 3 are both around 10 dB.
For radial velocity estimation as shown in Fig. 4  the AoA estimation MSE decreases as the receiving signal power increases. As the S-SINR is larger than −27 dB, the AoA estimation MSE is less than 0.5 square degrees. Since the range of azimuth angle, ϕ k , is larger than the elevation angle, θ k , the MSE of ϕ k is larger than θ k at first. With S-SINR becoming large enough, the MSE of ϕ k approaches that of θ k .  achieves better location MSE than scheme 3. The gaps between scheme 1 and scheme 3 are not so large in the high SINR regime. This is because the AoA estimation error dominates the location MSE. More specifically, E{(r s,l − r s,l ) 2 } is smaller than 10 −1 m 2 , while the error of location as shown in (83), (84), and (85) can be much larger than E{(r s,l − r s,l ) 2 }, because they are related to r s,l .

C. Communication Performance
We  , according to (9) and (11) under the assumption γ IN C = γ IN S = 3 dB. Due to the CSI estimation error caused by noise and interference, the BER for case B is significantly larger than that for case A. As C-SINR increases, the BER for case C decreases rapidly and becomes lower than that for case B after C-SINR is larger than 20 dB. This is because the JCS CSI enhancement method exploits the accurate sensing results and improves the estimated CSI. By comparing case D with cases B and C, we can see that the BER for case D is much larger, which indicates that the FFT-based JCS sensing results are not helpful for improving the CSI and for communication. Referring to the sensing MSEs of the MUSIC-based and the FFT-based JCS in Fig. 5(b), we can see that the more accurate the sensing results are, the better CSI enhancement performance is, as the accuracy of range estimation directly determines the accuracy of A in Algorithm 1. This is also the reason that the BER for case C decreases rapidly when the sensing MSE becomes sufficiently low.

VI. CONCLUSION
In this paper, we proposed a novel JCS system that can achieve accurate AoA, range, and velocity estimation based on improved MUSIC algorithms, together with improved communication performance. Compared with the conventional FFT-based sensing method, our proposed [ where P i is the power of the ith incident signal, σ 2 N is the noise power. We define the differential vector of v [v ∆ ] k /(N − (N − 1)/2 ), andv is close to 0. Therefore, the number of detected targets is determined asL where ε is a parameter used to avoid false detection caused by a small error. In the simulation, we set ε = 1.
Therefore, α t is set as α t = [v x ]L/m x . It is a key parameter and has an important impact on sensing accuracy. When α t is too large, the selected noise subspace will include part of the signal subspace, and thus the target may be missed; when α t is too small, the noise subspace is not selected completely, and thus large noise may be taken into the signal space. [v x ] i /(N c − 1). On one hand, since U x,rN is the noise subspace of R X,r , we have


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.