extfieldsofstudy:
Computer Science
Medicine

text:
Extracting automata from neural networks using active learning

Deep learning is one of the most advanced forms of machine learning. Most modern deep learning models are based on an artificial neural network, and benchmarking studies reveal that neural networks have produced results comparable to and in some cases superior to human experts. However, the generated neural networks are typically regarded as incomprehensible black-box models, which not only limits their applications, but also hinders testing and verifying. In this paper, we present an active learning framework to extract automata from neural network classifiers, which can help users to understand the classifiers. In more detail, we use Angluin’s L* algorithm as a learner and the neural network under learning as an oracle, employing abstraction interpretation of the neural network for answering membership and equivalence queries. Our abstraction consists of value, symbol and word abstractions. The factors that may affect the abstraction are also discussed in the paper. We have implemented our approach in a prototype. To evaluate it, we have performed the prototype on a MNIST classifier and have identified that the abstraction with interval number 2 and block size 1 × 28 offers the best performance in terms of F1 score. We also have compared our extracted DFA against the DFAs learned via the passive learning algorithms provided in LearnLib and the experimental results show that our DFA gives a better performance on the MNIST dataset.

INTRODUCTION
Deep learning is one of the most advanced forms of machine learning, which has been applied to various fields, including computer vision, speech recognition, natural language processing, audio recognition, social network filtering, machine translation, bioinformatics, drug design and board game programs (Schmidhuber, 2014;Lecun, Bengio & Hinton, 2015). Most modern deep learning models are based on an artificial neural network, such as deep neural networks (DNN), deep belief networks, convolutional neural networks (CNN) and recurrent neural networks (RNN). Benchmarking studies reveal that neural networks have produced results comparable to and in some cases superior to human experts.
However, the generated neural networks are typically regarded as incomprehensible black-box models. They are in practice unlikely to generalise exactly to the concept being trained, and what they eventually learn actually is unclear (Omlin & Giles, 2000). The opaqueness of neural networks not only limits their applications, but also hinders is found, it may be not that the hypothesis is incorrect, but rather that the abstract model is not precise enough and needs to be refined.
Finally, we have implemented our approach in Java, wherein we use the library LearnLib (Howar et al., 2012) to implement the active learning framework. To evaluate our approach, we conducted a series of experiments on a classifier for the MNIST dataset, a large database of handwritten digits that is commonly used for training various image processing systems. We first test the measures of the MNIST classifier, namely, safety, conflict, the size of alphabet and the length of words, under the abstractions with different interval numbers (i.e., the number of partitioning) and block sizes, and have identified some suitable abstractions. Secondly, we conduct some experiments to learn DFAs from the MNIST classifier with the suggested abstractions. The results shows that the abstraction with interval number 2 and block size 1 × 28 offers the best performance in terms of F1 score. At last, we also conduct the experiments to compare our resulted DFA against the DFAs learned via the passive learning algorithms (see "Passive Learning") provided in LearnLib and the MNIST classifier itself. Although worse than the classifier, our DFA gives a better performance than the other DFAs in our experiments. Nevertheless, there are still some limitations for our approach.
In summary, our contributions are as follows: We have proposed an MAT framework to extract automata from neural networks, employing abstraction interpretation of the neural network for answering membership and equivalence queries.
We have conducted several experiments on a MNIST classifier, and the experimental results show that our approach is viable, and the resulted DFA has a better performance in terms of F1 score than the DFAs learned via the passive learning algorithms provided in LearnLib on the MNIST dataset.
The remainder of this paper is organised as follows. "Preliminary" gives the preliminaries of DFA and active learning. "Approach" describes our approach, followed by the experimental results in "Experiments". "Limitations" discusses some limitations of our approach. "Related Work" presents the related work, followed by some concluding remarks in "Conclusion".

PRELIMINARY
In this section, we present the notion about neural networks, deterministic finite automata, active learning and passive learning.

Neural networks as functions
Neural network models can be viewed as mathematical models defining a possible nonlinear function N : X ! Y, where X is the input and Y is the output. In more detail, N can be defined as a composition of other (layer) functions, which can further be decomposed into other functions. This can be conveniently represented as a network structure, with arrows depicting the dependencies between functions. In this paper we assume that the framework of neural networks is unknown and only network-acceptors are considered. So the function we consider here is the one representing the whole network N : X ! Bool, where X is a multi-dimensional array. We say an input data X is positive if NðXÞ ¼ true, and otherwise negative.

Deterministic finite automata
Definition 2.1 A deterministic finite automaton (DFA) is a 5-tuple (Q, σ, δ, q 0 , F), where Q is a finite set of states, σ is a finite set of input symbols and is called the alphabet, d : Q Â AE ! Q is the transition function, q 0 ∈ Q is the starting state, F Q is the set of accepting states.
A word or string over an alphabet σ is a finite sequence of symbols from σ. The length of a word is the number of symbols it contains. Note that a word can be empty: the empty word, denoted as ε, has length 0 and contains no symbols.

Definition 2.2
Let M ¼ ðQ; AE; d; q 0 ; FÞ be a DFA and w = a 1 a 2 … a n be a word of length n over σ. The automaton M accepts the word w if and only if there exists a sequence of states r 0 ,r 1 , …, r n with the following conditions: The set of words recognised by a DFA M, called the language of M, is the following set: LðMÞ ¼ fw 2 Σ Ã jw is accepted by Mg Active learning framework Angluin (1987) proposed the first active learning algorithm, the L* algorithm, to learn finite automata from a MAT in 1987, and today all the most efficient learning algorithms that are being used follow Angluin's approach. In the following, we briefly introduce the MAT framework.
Finite automata can be learned precisely from a MAT, that is, an oracle capable of answering the so-called membership and equivalence queries: membership queries: the learner asks whether a given word is accepted by the automaton or not, and the teacher answers with the result. equivalence queries: the learner asks whether a given hypothesis automaton H is equal to the automaton model M held by the teacher. The teacher answers yes if this is the case. Otherwise she answers no and supplies a word, the so-called counterexample, on which the hypothesis automaton H and the automaton model M disagree.
The MAT framework is shown in Fig. 1. Initially, the learner knows the static interface of the SUL, that is, the sets of input (i.e., multi-dimensional array for neural networks) and output (i.e., yes or no for recognizers). Then the learner starts to ask a sequence of membership queries (MQs) and receives the corresponding responses from the teacher. After a "sufficient" number of queries, the learner builds a hypothesis H from the obtained information, and then sends an equivalence query (EQ). If the teacher answers yes, then the hypothesis H is returned. Otherwise, the learner refines the information with the returned counterexample, and continues on querying.

Passive learning
Different from active learning, passive learning constructs automata from sets of examples directly. Many approaches in grammatical inference can be described as passive learning. In the paper, we consider the polynomial-time RPNI algorithms provided in the library LearnLib (Howar et al., 2012). Oncina & Garca (1992) proposed the Regular Positive and Negative Inference (RPNI) algorithm for DFA learning. RPNI starts with a prefix tree acceptor, a tree-like DFA built from the learning examples by taking all the prefixes of the examples as states, and then greedily creates clusters of states (by merging) in order to come up with an automaton that is always consistent with the examples. Two heuristic strategies can be employed in state merging: Evidence Driven State Merging (EDSM) (Cicchello & Kremer, 2002) and Minimum Description Length (MDL) (Adriaans & Jacobs, 2006).

APPROACH
In this section, we present an active learning framework to extract automata from neural network classifiers. Our framework is shown in Fig. 2, which is a classic MAT framework with an abstraction 2 . In a nutshell, we make an abstraction between the learner and the SUL. When the membership queries are sent, the abstraction maps the abstract words into the concrete ones (i.e., the up arrow in Fig. 2), which are then fed into the neural network under learning. When the equivalence queries are sent, the abstraction does the opposite (i.e., the down arrow in Fig. 2) and from the abstract words an abstract representation is built for checking. In the following, we explain how to define an abstraction for neural networks and how to instantiate the active learning framework on neural networks. Abstraction For simplicity, we focus on network-acceptors, that is, there are only two outputs for the SUL, and thus we do not need to abstract them. In other words, only the inputs need to be abstracted. Generally, the inputs of neural network classifiers are always multidimensional arrays. As mentioned in "Introduction", we aim to abstract an input as a word, rather than a symbol. So the aim of the abstraction is to convert multi-dimensional arrays into words and vice-versa. Just like the serialization of multi-dimensional arrays, a naive and simple solution to the abstraction is to convert the input multi-dimensional array into a 1-dimensional array in row (or column) major order, and then concatenate the string representation of each value in the converted array in order, yielding a word. For example, Fig. 3 shows an array with size 4 × 4, on which the simple abstraction is applied. However, there are two problems for this solution: (1) the size of alphabet may be too large, even infinite.
(2) the length of the abstracted word may be too long. Both of them can make the automata learning too time-consuming.
So for practicality, we propose a three-layer abstraction, which consists of: value abstraction: each value in an input array is mapped into an integer via partitioning, which helps reduce the size of alphabet; symbol abstraction: a block of multi-dimensional integer array is abstracted as a symbol, which enables us to reduce the length of word; word abstraction: the whole input array is encoded into a word, wherein value abstraction and symbol abstraction are applied.

Value abstraction
In order to reduce the size of alphabet, inspired by Omlin & Giles (1996) work, we first split the values of the input space into n (equal) intervals, and map each interval into an integer, that is, the index of the corresponding intervals. Formally, let the input space be I and it be split into n intervals I 0 , …, I n − 1 . Then a value abstraction function a v : I ! f0; . . . ; n À 1g is defined as follows: This value abstraction function maps a concrete value in the input to an abstract integer. Figure 4 shows an example of value abstraction applying on the array given in Fig. 3, where the input space I is [0,1] and I is split into I 0 = [0,0.5) and I 1 = [0.5,1].
While for concretization, an abstract integer is mapped to its corresponding interval, that is, the value concretization function g v : f0; . . . ; n À 1g ! 2 I is defined: Both a v and β v can be extended on sets of elements in a natural way: It is easy to get that a v (γ v (i)) = i for each integer i, and d ∈ γ v (a v (d)) for any given value d. Therefore, (a v , γ v ) forms a Galois connection (Nielson, Nielson & Hankin, 1999). While in practice, especially when we query a word, we are unable to test all the values in the interval for each integer. For that, we randomly select at most k v (which can be dependent on the intervals) values to represent the corresponding interval. That is to say, we define a weak value concretization function: g prac v ðiÞ ¼ fd j j d j 2 I i and 0 j < k v g Obviously, the larger k v is, the closer g prac v ðiÞ is to γ v (i) . So concerning Galois connections, the larger k v , the better. But only Galois connections are not enough here. We also need to consider the safety of neural networks (Huang et al., 2017), that is, a vibration of values should not flap the outputs, since different values may be abstracted into an identity integer. In fact, the composition of the functions a v and γ v can be viewed as a kind of manipulations (Huang et al., 2017). We say a k-value manipulation vm k with respect to a v and γ v is a function such that for any input array in . Intuitively, k-value manipulation replaces (at most) k values of the input array by some values, which share the same intervals with the corresponding original values. Figure 5 shows an example of 4-value manipulation applying on the input array in Fig. 4, where the input space I is [0,1] and it is split into I 0 = [0,0.5) and I 1 = [0.5,1]. And we say a network N is safe with respect to this value manipulation vm k if for every input array in

Nðvm k ðinÞÞ ¼ NðinÞ
That is to say, performing this value manipulation should not result in a different classification. This requires that every interval I i should be as small as possible or the number n of intervals should be as large as possible. However, the safety cannot easily be preserved in practice, unless the abstraction is an identity function or the network is robust enough. So instead, we use a weak notation called σ-safety: we say a network N is σ-safe with respect to k-value manipulation vm k under a given input set D if

Symbol abstraction
After the value abstraction, each integer can be used as a symbol. But this could yield words that are too long to learn the model. So for scalability, we add a symbol abstraction, which abstracts input arrays into symbols by blocks. For simplicity, in this paper we consider 2-dimensional array with size iRow × iCol. We say a slice of an input array starting from the index (ri, ci) to the index (ri + oRow − 1, ci + oCol − 1) is a block, and the size of the block is oRow × oCol. A natural way to abstract blocks into symbols is to map the blocks into one dimension in row (or column) major order and then encode the one dimension into a base-n number (or a string consisting of the integers in the one dimension). We denote this mapping as a B s . Figure 6 gives an example of a B s that are performed on the array that are obtained by the value abstraction shown in Fig. 4, where the number n of intervals is 2 and the size of blocks is 2 × 2. Moreover, by decoding the base-n number (or the string), it is easy to obtain the inverse mapping g B s . It is clear ða B s ; g B s Þ forms a Galois connection. But a drawback of this solution is that the size of alphabet is n oRow × oCol , which could be too large in practice. For example, the size of alphabet of the example shown in Figure 6 is 16.
In this paper we use an alternative way to represent a block as its sum. In more detail, we define a symbol abstraction function a S s that maps integer blocks of size oRow ×oCol into the sum of the integers in blocks: It is easy to compute the set of the possible sums of blocks, that is, {0,…, oRow×oCol× (n − 1)}. So the size of alphabet is oRow×oCol× (n − 1) + 1. Compared to the natural solution, the size of alphabet is quite smaller (from n oRow × oCol reduced to oRow×oCol× (n − 1) + 1). Take the input array in Fig. 6 for example. Under this sum abstraction, its abstraction is shown in Fig. 7 and the size of alphabet is 5.
But it is pity that this mapping is not bijective. So in order to form Galois connections, similar to value abstraction, we define the inverse mapping from symbols (i.e. sums) to sets consisting of blocks of size oRow × oCol whose sum is exactly the symbol: Likewise, these two functions can be lifted to sets of elements in a natural way. It is easy to get that a s (γ s (sum)) = sum for each symbol sum, and b ∈ γ s (a s (b)) for any given block b of size oRow × oCol. Therefore, (a s , γ s ) forms a Galois connection. While for practicality, similar to value abstraction γ v , we use a weak sum concretization function: That is, we select at most k s (which can be dependent on the block size and n) blocks to represent the corresponding sum. Clearly, the larger k s is, the closer g Sprac s ðsumÞ is to g S s ðsumÞ. So concerning Galois connections, the larger k s , the better. In addition, the composition of a S s and g S s forms a manipulation as well, and the network N should be safe with respect to this manipulation. Formally, we say a k-shift manipulation sm k with respect to a S s and g S s is a function such that for any input array in where b i is a block of size oRow ×oCol belonging to in and b 0 i 2 g S s ða S s ðb i ÞÞ. And a network N is safe with respect to sm k if for every input array in Nðsm k ðinÞÞ ¼ NðinÞ In a word, k-shift manipulation replaces (at most) k blocks of the input array by some blocks, which share the same sums with the corresponding original blocks; and performing this k-shift manipulation should not result in a different classification. Similar to value manipulation, the safety cannot easily be preserved in practice. So we use a weak notation called σ-safety: a network N is σ-safe with respect to sm k under a given input set D if jfin 2 DjNðsm k ðinÞÞ 6 ¼ NðinÞgj jDj r where |D| ≥ 1. Considering the safety, the distance between two blocks of the same sum or the size of block should be as small as possible. Figure 8 shows an example of 1-shift manipulation applying on the input array in Fig. 7, where n is 2 and the size of blocks is 2 × 2. Finally, let us consider the alphabet size. As discussed above, the alphabet size for the abstraction function a S s (a s resp.) is linear (exponential resp.) in the block size and the number of intervals n. So concerning the alphabet size, the smaller the block size and the interval number n, the better.

Word abstraction
Finally, we split the input array into blocks, and map them into a sequence of symbols (i.e., a word) in row (or column) major order. Algorithm 1 shows the detail of the word abstraction function a w . The algorithm first invokes the value abstraction a v to map the values in the input array into integers (Line 1). Then it slides over the integer array block by block (Lines 2-10) and maps each block into a symbol by the symbol abstraction a S s (Line 6). Note that here we use a narrow slide on the input array, that is, the blocks to be abstracted are fully contained in the input array. One can use the wide slide with zeropadding as well. Just like the convolution operation of CNN, one can further set the stride sizes for each dimension. In addition, one can further encode the sequence of symbols into a final word in a more compact format, such as run-length encoding (RLE).
As mentioned above, the symbol abstraction aims to reduce the length of words. According to the word abstraction, we have that the larger the block size, the shorter the word. Figure 9 shows our abstraction applying on the input array given in Fig. 3, where the number n of intervals is 2 and the size of blocks is 2 × 2. Compared to the simple abstraction shown in Fig. 3, our abstraction yields words with smaller size of alphabet and shorter length.
The word concretization function γ w , which is shown in Algorithm 2, does the opposite: it maps a sequence of symbols (i.e., a word) into a sequence of sets of blocks (Lines 5-7), and combines them into a set of arrays (Lines 8−10). Note that we require the length of Algorithm 1 Word abstraction function α w (in).

Input: an input in
Output: a word w 7: ci = ci + oCol 8: end while 9: ri = ri + oRow 10: end while 11: return w word to be concretized should conform to the size of input (Lines 2−4). One can release this length condition by zero-padding or discarding the superfluous symbols. But this may break the Galois connections.
Theoretically, if (a v , γ v ) and ða S s ; g S s Þ form Galois connections, then so does (a w , γ w ). While for practicality, we use , and the number of data in g prac w ðwÞ depends on the word w as well as k v and k s . In particular, in our implementation we collect sets of inputs (including values and blocks) that are mapped into an identity word from existing data and then select inputs from the corresponding set.
Finally, considering safety, if (a v , γ v ) and ða S s ; g S s Þ does not cause the flapping, then neither does (a w , γ w ). But the safety of a network focus on the vibrations of local parts of Algorithm 2 Word Concretization Function γ(w).

Input: a word w
Output: a set matrix_set of arrays inputs (Huang et al., 2017). To evaluate the whole inputs, we use another notation conflict, that is, inputs of different classifications should not abstracted into an identity word. Formally, we say a network N is non-conflict with respect to a w and γ w if for every input array in and for every array in′ ∈ γ w (a w (in))

Nðin 0 Þ ¼ NðinÞ
In other words, the abstraction itself should not be over-approximated. We say a word w is conflict, if there exist two inputs of different classifications that are abstracted into it. So to avoid over-approximation, the number of the conflict words caused by the abstraction should be as few as possible. Similar to the safety, the non-conflict cannot easily be preserved in practice. For that, we evaluate the conflicts words under a given dataset. In detail, we say a network N is σ-conflict with respect to a w under a given set D if jfin j 9in 0 :NðinÞ 6 ¼ Nðin 0 Þ^a w ðinÞ ¼ a w ðin 0 Þgj jDj r where |D| ≥ 1 and in, in′∈ D.
To sum up, to obtain a suitable abstraction (e.g., scalable, safe and non-conflict), one needs to take into account the number n of intervals, the block size oRow and oCol, and the other factors.

Active learning
In this section, we present how to instantiate the active learning framework on neural networks, in particular the membership and equivalence queries.

Membership query
Membership queries can be answered by the neural networks via the word concretization function. In our abstraction, we map a word into a set of data. As mentioned above, the abstraction may flap the results or yield some conflict words, that is, the classifications of different data in the set of an identity word may not be the same. To address this, we count the numbers of different classifications of the data in the set and take the classification which gets the most votes as the result for the word.
Given a network N and a word concretization function a w , we say a word w is positive if jfin 2 g w ðwÞjNðinÞ ¼ truegj ! jfin 2 g w ðwÞjNðinÞ ¼ falsegj and otherwise negative. Intuitively, a word is positive (negative resp.) if there are more positive (negative resp.) input arrays that are abstracted into it than the negative (positive resp.) ones.
Algorithm 3 gives the procedure of membership query checking, where N denote the neural network under learning. Firstly, the algorithm concretises the word w that is being queried into a set matrix_set of possible data, using the word concretization function γ w (Line 1). If the set matrix_set is null, that is, the length of word w does not conform to the size of input data, then the algorithm returns false immediately. Otherwise, the algorithm feeds each data into the neural network N under learning and counts the numbers of different classifications (Lines 5-12). Finally, it returns the classification that gets the most votes (Line 13).

Equivalence query
As there is no finite interpretation for neural networks (Weiss, Goldberg & Yahav, 2018), equivalence queries are more challenging than membership queries. To address this, similar to Weiss, Goldberg & Yahav (2018)'s work, we use an abstract representation of the neural network under learning. But different from Weiss, Goldberg & Yahav (2018)'s work, we start with the automaton that is learned passively via the RPNI algorithm (Oncina & Garca, 1992) from some test queries, which are selected from the training dataset. Then we perform the equivalence query against this abstract model. As discussed in Weiss, Goldberg & Yahav (2018), when a counterexample is found, it may be not that the hypothesis is incorrect, but rather that the abstract model is not precise enough (i.e., different behaviors from the neural network under learning) and needs to be refined.
The procedure 3 of equivalence query checking is given in Algorithm 4. Firstly, the algorithm tries to find a word that can separate the hypothesis H and the abstract model M (Line 3). If such a word does not exist, then it returns null (Lines 4−6), which means the equivalence query is yes. Assume a word w is found. Then it checks whether this word is a true counterexample, that is, the classifications of the abstract model and the neural network under learning are the same (Line 7). If it is in that case, then it returns this word as a counterexample to the learner (Line 8). Otherwise, it refines the abstract model with this word (Line 10): it adds the counterexample into the positive set or the After that, the algorithm continues on the equivalence query against this refined model.

EXPERIMENTS
We have implemented our approach in a prototype in Java, wherein we use the library LearnLib (Howar et al., 2012) to implement the MAT learning framework and the RPNI algorithm. Moreover, to find the true counterexamples faster, we use the Wpmethod test (Fujiwara et al., 1991) 4 in the equivalence query between the hypothesis and the abstract models. To evaluate our approach, we conduct a series of experiments on a classifier for the MNIST dataset, a large database of handwritten digits that is commonly used for training various image processing systems. Firstly, we conduct experiments to see the measures of the MNIST classifier, namely, σ-safety, σ-conflict, the size of alphabet and the length of words, under the abstractions with different interval numbers and block sizes. Secondly, we present the experiments to learn DFAs from the MNIST classifier under different selected abstractions. Thirdly, we also conduct experiments to compare the resulted DFAs against the DFAs learned via the passive learning algorithms provided in LearnLib and the MNIST classifier itself. The experiments were conducted on a workstation with Intel Processor i7-7820HQ (2.90GHz) and 32GB memory.

MNIST classifier
The MNIST classifier under learning is a binary classification version of MnistClassifier from the tutorial examples of DeepLearning4J (https://github.com/deeplearning4j/dl4jexamples), which recognises the number 1. It is built on a convolution neural network, which consists of six layers, namely, a convolution layer, a pooling layer, another convolution layer, another pooling layer, a dense layer and an output layer. The training

Abstraction experiments
As discussed in "Abstraction", the interval number n and the block size affect the definition of the abstraction, especially the safety and the conflict of the neural network under learning, the size of alphabet and the length of words. For that, we present in this section some experiments to see these measures of the abstractions with different interval numbers and block sizes.

Safety
The first measure to test is the safety. For that, we present some experiments to test the flapping of the MNIST classifier on some selected inputs from the training set (i) via performing k-value manipulations vm k with different interval numbers and (ii) via performing k-shift manipulations sm k with different block sizes. First, in the experiments about k-value manipulation vm k , for a given interval number n, we randomly select k values from a selected input, and replace each selected value by a random value which shares the same interval with the corresponding selected value. Then we fed the resulted data into the MNIST classifier and see whether the classifications are flapped. We select 59,838 inputs in total from the training set, which are classified correctly by the MNIST classifier. Table 1 shows the results, where Flaps denotes the number of inputs whose results are flapped by the manipulation, and Ratio denotes the percentage of the number of flapped input to the total number of selected inputs.
From the results we can see that, the number of flapped inputs increases as the number k of selected values increases, since the larger the number k, the larger the vibration for the inputs. In contrast, as the number of intervals increase, the number of flapped inputs decreases, which indicates that the larger the interval number, the better. This conforms to the discussion in "Abstraction". Moreover, the results also show that the MNIST classifier is about 0.053%-safety, with respect to the 100-value manipulation vm 100 with the interval number 2. And the 100-value manipulation means 12.76% (100/784) of an input has been modified, such that we believe 100 is enough for local vibration. Therefore, we suggest to set the interval number n as 2. Next, in the experiments about k-shift manipulation sm k , for a given block size oRow × oCol, we randomly select k blocks from a selected input, and rearrange the values in each selected block. Then we feed the resulted data into the MNIST classifier and see whether the classifications are flapped. Similarly, we select the 59,838 inputs that are classified correctly by the MNIST classifier from the training set. For simplicity and scalability, we consider the block sizes whose row sizes or column sizes are 28. The results are given in Table 2, where the notations are the same as the ones in Table 1.
First, the results show that, as the size of block increases, the number of flapped inputs increases, which conforms to the discussion in "Abstraction". The results also show that the number of flapped inputs increases as the number k of selected blocks increases. This is because that, the larger the number k, the larger the vibration for the inputs. Moreover, we found that the MNIST classifier is more safe under k-shift manipulation built on rows than the one on columns. The reason may be that the digit number of 1 is more regular in row order than in column order. Finally, assume the size allowed for local vibration is about 100. All the σ-safeties of the MNIST classifier with respect to the k-shift manipulation with the block size 1 × 28, 2 × 28, 4 × 28, 28 × 1, 28 × 2 or 28 × 4 are smaller than 1.8%. In particular, the MNIST classifier is about 0.055%-safety, with respect to the 4-shift manipulation vm 4 with the block size 1 × 28.

Conflict
The second measure to test is the non-conflict, which indicates whether the abstraction with the given block size is over-approximated. In other words, we would like to conduct experiments to test how many conflict words that are generated by the abstractions with different block sizes under the training set. For that, we perform the abstractions with different block sizes on some selected inputs from the training set, and do a statistic analysis on the abstracted words with respect to their classifications, where we take the interval number n as the suggested one 2. The test inputs that are selected from the training set is 59,840 in total, with 6,700 positive inputs and 53,140 negative ones.
The statistic results are given in Table 3, where TW denotes the total number of words, PW (NW resp.) denotes the number of positive (negative resp.) words, CPW (CNW resp.) denotes the number of positive (negative resp.) words that have both positive and negative inputs and CPD (CND resp.) denotes the number of positive (negative resp.) inputs that are abstracted into a negative (positive resp.) word.
The results show that as the block size increases, the number of abstracted words decreases, which conforms to the discussion in "Abstraction". Thus it could be easier to extract the automaton for a larger block size. For example, when taking the whole input as a symbol, there are 239 words in total. But both the number of conflict words and the number of conflict data increase as the block size increases, which indicates that an abstraction with a larger block size is prone to be an over-approximation. In particular, when taking the whole input as a symbol, 56.194% of the positive inputs are abstracted into negative words and 64.865% of the positive words are conflict. Moreover, from the results we can also see that all the σ-conflicts for the MNIST classifier with respect to the abstractions with the block size 1 × 28, 2 × 28, 28 × 1 or 28 × 2 are smaller than 0.015% (8/59,840). And the abstractions with the block size 2 × 28 and 1 × 28 perform best, yielding none conflict data nor words.

Word complexities
Finally, we also conduct experiments to see the size of alphabets and the length of words, which are dubbed as word complexities. Table 4 shows the word complexities under different abstractions with different block sizes, where Size denotes the size of alphabet, dSize denotes the number of symbols occurring in the selected inputs and Length denotes the length of words.
The results show that the larger the block size, the larger the alphabet size and the shorter the word length, which conforms to the discussion in "Abstraction". Moreover, we found that the products of the alphabet size and the word length are almost the same. So for scalability, any block size seems fine. But if considering the practical alphabet (i.e., symbols occurring in the inputs), the larger block size could be better.
To sum up, based on the experiments above, we suggest to use for the MNIST classifier the abstractions with the interval number n = 2 and the block size 1 × 28, 2 × 28, 28 × 1 or 28 × 2.

Automata learning
In this section, we present the experiments to learn DFAs from the MNIST classifier under the suggested abstractions.
To quantitatively validate the models, we use the following performance measures. Accuracy is the most intuitive performance measure and it is simply a ratio of correctly predicted observation to the total observations. Precision is the ratio of correctly predicted positive observations to the total predicted positive observations, and Recall is the ratio of correctly predicted positive observations to all observations in actual class. F1 score is the weighted average of Precision and Recall, that is, (2 · Precision · Recall)/(Precision + Recall). Moreover, there are two kinds of observations in our experiments, namely, the input arrays and the abstracted words (i.e., the abstractions of the input arrays). So we compute the measures above with respect to both kinds of observations. Intuitively, the higher the measures above, the better the model.
The results show that all the learned DFAs perform well on the testing dataset, with the F1 score more than 50% 5 . In particular, the DFAs learned under the abstraction with block size 1 × 28 performs best. The results also show that DFA learned via the abstraction with a smaller block size can obtain a higher F1 score. In detail, the F1 score of DFA learned via the abstraction with block size 1 × 28 (28 × 1 resp.) is higher than the one with block size 2 × 28 or 28 × 2 (28 × 2 resp.). This is because that a smaller block size can generate a more preciser abstraction, which conforms to the discussion in "Abstraction". Moreover, from the results, we can see that the DFAs learned under the abstractions in rows perform better than the ones under the abstractions in columns in terms of all performance measures with respect to both words and data. For example, the F1 score of DFA learned via the abstraction with block size 1 × 28 is higher than the one with block size 28 × 1. The reason may be that the digit number of 1 is more regular in row order than in column order. In addition, we also perform the abstraction mapping a whole input as a symbol as does in Weiss, Goldberg & Yahav (2018)'s work. Due to this abstraction is over-approximated, the extracted DFA gets the worst performance in the input data layer, although it has a better performance than the other models on the word layer, especially the Precision.

Learning complexities
During the experiments, we also count the learning times in seconds needed by the resulted DFAs and the sizes of the resulted DFAs. The results are given in Table 6, where aTime, xTime and iTime respectively denote the average time, the maximum time and the minimum time needs by the resulted DFAs, and aState, xState and iState respectively denote the average number, the maximum number and the minimum number of states of the resulted DFAs.
From the results, we can see that learning the DFA via a smaller block size needs more time. As discussed before, an input array can be abstracted into a longer word under the abstraction with a smaller block size, which thus requires more time to proceed. Concerning the size of learned DFA, learning via a smaller block size can yield a larger DFA. For example, the number of the states of the DFA learned under the abstraction with block size 28 × 1 is the largest one among the results. Similar to the learning time, the reason is that an abstraction with a smaller block size yields longer words, which could enlarge the learned DFA. In addition, both the learning times and the sizes of the learned DFAs under the abstractions in rows are larger than the ones under the abstractions in columns. The reason may be that more blocks are abstracted into 0 under the abstractions in columns than the abstractions in rows.

Learned automata
Finally, we convert the learned DFAs into the format used in JFLAP (Rodger & Finley, 2006), which enables us to view the DFAs and convert DFA into regular expressions step by step. For simplicity, we consider a DFA learned via the abstraction with block size 14 × 28, which is given in Fig. 10. From this figure, we can see that (i) the learned DFA has a very clear structure: a starting node, an intermediate layer with several nodes, an accepting node and a trap node; and (ii) the learned DFA accepts words with length 2.
Let us see a DFA learned via the abstraction with block size 7 × 28, which is shown in Fig. 11. Compared with the one in Fig. 10, this DFA has a more complex structure. But we still can identity some hierarchical structures in it. It is pity that we are not able to convert this DFA into a regular expression via JFLAP due to a runtime error.
In addition, we also present a DFA learned via the abstraction with block size 1 × 28, which has 949 states and is given in Fig. 12. It is a little complex to understand, so we can only identity a rough hierarchical structure. We believe that one can understand this DFA more if he gets a closer look on it. In addition, the learned DFAs can help to generate test cases to test the networks, which are left as a future work.

Comparison
To further evaluate the resulted DFA, we compare it against the DFAs learned via the passive learning algorithms provided in LearnLib, namely, the RPNI algorithm, the RPNI- EDSM algorithm and the RPNI-MDL algorithm, and the MNIST classifier itself. In these experiments, we perform our abstraction on the training data and then learn a DFA via each passive learning algorithms provided in LearnLib, wherein the abstraction we used is the one with interval number 2 and block size 1 × 28, all the arrays are selected for , and all the positive arrays and only one ninth of the negative ones are selected for the RPNI-EDSM algorithm (to avoid memory overflow and time-consuming). Next, we evaluate all the models with the testing dataset. The results are given in Table 7, where the notations are the same as the ones of Table 5.
Compared to the RPNI one, our DFA performs better in all the performance measures, since the abstract model we use is the DFA learned via the RPNI algorithm from some inputs in the training dataset and is refined with respect to the classifier continually during learning. While compared to the RPNI-MDL and RPNI-EDSM ones, our DFA has a better Accuracy, Precision and F1 score, but a worse Recall. This is because that these two DFAs take all the positive inputs in the training set into account such that it can recognise more positive inputs in the testing dataset, while only part of positive inputs are selected for our abstracted model. The results also show that our DFA is still worse than the classifier. There are several reasons for this. The first one is that we have set some bounds (e.g., the refining time for the abstract representation) in our implementation for the learning procedure for efficiency and to avoid memory overflow. The second one is that the Wp-method test used in our experiments may miss some true counterexamples. The third one is that the abstraction may be over-approximated to yield too many conflict words. Nevertheless, our approach still needs to be improved.

LIMITATIONS
Although our approach works for the MNIST classifier, there are still some limitations. Firstly, to figure out a suitable abstraction for the neural network under learning is not an easy task. As shown in Biggio et al. (2013), Szegedy et al. (2013) and Huang et al. (2017), several DNN, including highly trained and smooth networks optimised for vision tasks, are unstable with respect to so called adversarial perturbations. Hence, some neural networks may be too sensitive to the abstraction manipulation to find a reasonable interval number. Even if a reasonable interval number were found, one need to make a compromise between the abstraction and the scalability to find a block size. Moreover, whether a turing machine can simulate a natural neural network is an open question (Zenil & Quiroz, 2006). So in some sense, we cannot define an abstraction without the conflict or the flapping.
Secondly, the scalability is another problem. Generally, the size of inputs of neural networks is in thousands. For such a neural network, either the alphabet may be too large (if a large block size is taken) or the word may be too long (if a small block size is taken) for us to extract the automaton. Taking the MNIST classifier for example, it may last several hours for some abstractions (see the abstraction with block size 1 × 28 in Table 6) to extract the automaton, even several days. There are two possible reasons for this issue in our implementation: (i) we use the implementation of RPNI algorithm from LearnLib, which does not support incremental learning proposed in Dupont (1996) and can be improved with it; and (ii) there are too many queries for the non-accepting words with invalid lengths. Thirdly, our approach is dependent on the dataset. In "Experiments", we selected the interval number and the block size via an analysis on the training dataset. Different datasets may derive different abstractions. To make things worse, it may be the case that an abstraction is suitable for the training dataset, but unsuited for some other testing dataset. Moreover, our abstract model is built from some existing testing data. Different data yields different abstract models, which could affect the results, such as the learning time and the learned DFA.
Fourthly, the implementation of equivalence query is a practical problem requiring attention. One may think that a precise equivalence check can be performed in polynomial time on the hypothesis automaton and the abstract automaton. However, the precise equivalence check could return too many false counterexamples such that it takes too much time for the learning. This is because the precise equivalence check is prone to generate a short and false counterexample that is invalid with respect to the abstraction. Indeed, we have tried this precise equivalence check, but we only succeed on the abstraction of 28 × 28 in 1 h. So considering the efficiency, we use the Wp-method test in the equivalence check, which enables us to find the counterexamples whose lengths are in a given range.
Fifthly, although our approach is black-box, the structures of neural networks may affect the performances of the networks themselves, so as the performances of the learned DFAs. We have performed our approach with the abstraction 2 × 28 on the MNIST classifiers whose hidden layer numbers range in {1, 2, 5, 10}, which are binary classification versions from the tutorial examples of DeepLearning4J. The accuracies of all the classifiers are above 99%. Figure 13 shows the F1 scores with respect to words and input data of Figure 13 The F1 scores with respect to words and input data via the abstraction 2 × 28.
Full-size  DOI: 10.7717/peerj-cs.436/ fig-13 the learned DFA. We can see that the F1 scores are quite close to each other. And the DFA learned from the larger network is not necessary to get the best performance in terms of F1 score. Therefore, we present the experimental results on only one MNIST classifier in "Experiments".

RELATED WORK
In this section, we review some related work. Existing work on DFA extraction from neural networks targets RNNs, which was extensively explored in Jacobsson (2005) and Wang et al. (2017). Omlin & Giles (1996) proposed a global partitioning of the network state space according to q equal intervals along every dimension, and then exploring the network transitions in the partitioned space. Our value abstraction adopts this partitioning, but we work on the input space, instead of the state space. Cechin, Simon & Stertz (2003) presented a approach to extract DFA using k-means and fuzzy clustering. The key idea is to classify a large sample set of reachable network state using k-means. Hou & Zhou (2018) proposed another approach to extract DFA from RNN using two clustering algorithms, namely LISOR-k and LISOR-x, on hidden states. There are several other work that adopted cluster analysis on state space, including k-means clustering (Zeng, Goodman & Smyth, 1993;Frasconi et al., 1996;Gori et al., 1998;Cohen et al., 2017), hierarchical clustering (Sanfeliu & Alquezar, 1994) and selforganizing maps (Tiňo & Šajda, 1995). These approaches have to access the state-vectors, while our approach is a black-box one.
Recently, Weiss, Goldberg & Yahav (2018) adopted active learning to extract automata from RNN. Our work is inspired by and similar to this, but different in the follows: (1) we target general neural network, not only RNN; (2) we consider an input is a word, rather than a symbol; (3) we use a DFA that is inferred from some training data as an abstract model for equivalent queries.

CONCLUSION
In this work, we have proposed a MAT framework to extract automata from neural networks, employing abstraction interpretation of the neural networks for answering membership and equivalence queries. We have implemented our approach in a prototype and have carried out some interesting experiments on a MNIST classifier. Through experiments, we have found that the DFA extracted from the MNIST classifier under the abstraction with the interval number 2 and the block size 1 × 28 performs the best. In the experiments, that our resulted DFA has a better performance than the DFAs learned via the passive algorithms provided in LearnLib on the MNIST dataset.
As for future work, we may consider a better encoding such as RLE to improve the approach. We can improve the RPNI algorithm with incremental learning to reduce the learning time. We can also perform experiments on other neural network classifiers. Other models to be extracted from neural network are under consideration.


Second order Kalman filtering channel estimation and machine learning methods for spectrum sensing in cognitive radio networks

We address the problem of spectrum sensing in decentralized cognitive radio networks using a parametric machine learning method. In particular, to mitigate sensing performance degradation due to the mobility of the secondary users (SUs) in the presence of scatterers, we propose and investigate a classifier that uses a pilot based second order Kalman filter tracker for estimating the slowly varying channel gain between the primary user (PU) transmitter and the mobile SUs. Using the energy measurements at SU terminals as feature vectors, the algorithm is initialized by a K-means clustering algorithm with two centroids corresponding to the active and inactive status of PU transmitter. Under mobility, the centroid corresponding to the active PU status is adapted according to the estimates of the channels given by the Kalman filter and an adaptive K-means clustering technique is used to make classification decisions on the PU activity. Furthermore, to address the possibility that the SU receiver might experience location dependent co-channel interference, we have proposed a quadratic polynomial regression algorithm for estimating the noise plus interference power in the presence of mobility which can be used for adapting the centroid corresponding to inactive PU status. Simulation results demonstrate the efficacy of the proposed algorithm.

Introduction
The enormous growth of wireless communications has created an ever increasing demand for frequency spectrum. Consequently, pressure is being mounted on the service providers to meet the end users' requirement for higher bandwidth. The existing policies governing spectrum management are however, based on static spectrum allocation for a specific technology and service [1]. According to this static spectrum allocation policy, the systems are bound to operate within the assigned frequency band which often results in under-utilization and scarcity of the limited wireless spectrum [2,3]. The wastage is owing to the fact that when there is no transmission, the assigned frequency bands remain idle. Several studies and measurements have been conducted around the world regarding spectrum utilization. The results point to the fact that spectrum utilization is concentrated only on certain portions of the allocated spectrum. Besides, a significant portion of the allocated spectrum remains unutilized while a considerably large portion is used occasionally [4]. A paradigm shift from static spectrum allocation currently being practiced to dynamic spectrum access is required [5]. Hence, the problems of under-utilization and spectrum scarcity have received a lot of attention in the global research community. Novel techniques for optimising spectrum usage are also being proposed to provide the urgently required panacea for addressing the problem of the unprecedented demand for spectrum access [6]. Mitola et. al. in [7] and [8] proposed the concept of cognitive radio (CR) as an effective solution to the problem of under-utilization and scarcity of the wireless spectrum. A CR is an intelligent radio device that can reason and adapt itself to the changes in activities within its radio frequency (RF) environment. It can autonomously detect which communication channels are in use and which are not, quickly migrate into the vacant channel while avoiding the occupied ones. Upon finding a vacant band to use in an opportunistic manner, it is required that a CR quickly vacates such band immediately the original assignee returns. The CR technology would enable wireless devices to efficiently utilize the spectrum in a dynamic manner. The users in the CR network can be categorized into two, namely; primary users (PUs) and secondary users (SUs). The PUs are the licensed users that have the right to use a specific portion of the spectrum, and the SUs are the unlicensed users that use the spectrum allocated to the PUs in an opportunistic manner, without causing harmful interference to the licensed users [5,9,10]. Spectrum sensing is the process through which CR devices are able to detect unused frequency bands in their RF neighbourhood. Thus, it is core to the successful implementation of the CR technology [9,11].
Over the years, a lot of techniques have been studied and proposed for achieving spectrum sensing. Most popular methods include energy detection, matched filter and cyclostationary feature detection. All these techniques involve the detection of PU activities solely by utilising observations made by the SU. It is evident from literature that the energy detection technique, which is a non-coherent detection method is the most widely adopted. This is because it is very simple to implement compared to other techniques and does not require prior knowledge of the PU signal. However, it performs poorly in shadowing and fading environment especially when PU operates at low signal-to-noise ratio (SNR) [12][13][14]. The matched filter based detector compares the PU's signal with a set of received signals in order to determine the presence or absence of the PU. Although, it is a very effective method for performing spectrum sensing, it requires that the SU has prior knowledge of the patterns in the PU's signal [9]. Cyclostationary detectors take advantage of the certain statistical properties that are inherent in the PU's signal which vary cyclically with time. These features also have a periodic statistics and spectral correlation that cannot be found in any interference signal or stationary noise. Thus, the method is effective for detecting the presence of the PU transmissions. As in the case of the matched filter, though, the scheme requires prior knowledge of the PU's signal patterns at the SU. Besides, the method would fail if the PU alters any of these properties without notifying the SU. Furthermore, to mitigate the problem of shadowing, cooperative sensing has been proposed whereby multiple SUs are deployed at different geographical locations to perform local sensing of the PU signal [15]. The scheme exploits the spatial diversity of the multiple nodes by combining the local sensing information to make an overall decision about PU activities [16,17].
It is noteworthy that CR devices need to be equipped with intelligence in order to be able to perform reasoning functionalities [5,18]. To this end, various machine learning (ML) schemes have been proposed to effectively address the problem of spectrum sensing. For example, in [19], the support vector machine, K-nearest neighbours, Kmeans clustering and Gaussian mixture model methods are proposed while the unsupervised scheme based on the variational Bayesian learning technique is presented in [17]. A comprehensive review of parametric and nonparametric learning methods previously proposed is provided in [20]. Of late, the deep learning methods based signal classification techniques have been considered in [21][22][23][24] and [25].

Motivations and related works
The capability of parametric ML methods to effectively address the spectrum sensing problem is demonstrated in [19]. However, the performance of these proposed schemes would be severely degraded in a built up area due to the slowly fading channel occasioned by the mobility of SUs in the presence of scatterers [17,26].
A consideration of existing literature shows that until now, there is dearth of study considering the deployment of these machine learning schemes in realistic scenarios involving the mobility of PU or SUs. In [27], the compromising effect of multipath fading and shadowing on the detection performance of spectrum sensing schemes was studied. In particular, overhead problem caused by the cooperation among SUs received attention and a reinforcement learning-based cooperative sensing scheme is proposed to address it. Also, in [28], the K-means clustering based learning method is proposed for enhancing the performance of cooperative spectrum sensing in generalized j À l fading channels. Furthermore, concern about nodes placement in ML based cooperative sensing schemes was addressed in [29].
This study considers scenarios involving mobile SUs operating under flat fading channel condition and investigates the deployment of parametric ML methods for spectrum sensing under these conditions. To address the attendant sensing performance degradation, a real time channel tracking and decision boundary updating strategy is developed.

Contributions
In this paper, extending our work in [26], we consider a CR network of mobile SUs and investigate the K-means clustering technique 1 in a non-cooperative sensing scenario. The K-means clustering algorithm like other similar machine learning algorithms for pattern classification lends itself readily to decision plane tracking by updating the clusters centroids. However, different from [26] where the approximation of the time-varying channel was based on the first order Kalman filter combined with a first-order auto-regressive model, we propose a pilot based secondorder (order-2) Kalman filter channel estimation technique for tracking the changes in the fading channel, and demonstrate the advantage within the context of spectrum sensing to improve the performance of the PU detection system. In addition, unlike the work in [26], no explicit cooperation is assumed between the PU and SUs. The second order Kalman filter yields a better approximation of the time varying channel [31,32]. Furthermore, in comparison with its higher order counterparts, the order-2 Kalman filter scheme requires estimation of fewer parameters. The main contributions of this paper are as follows: -We have addressed the problem of spectrum sensing by mobile secondary devices under Rayleigh fading channel condition and modeled the task of monitoring the ensuing channel as a real time channel gain tracking problem. A pilot based second order Kalman filter algorithm has been developed for addressing it. -Using the knowledge of the Kalman channel gain estimates, we have proposed a completely unsupervised (no cooperation from PU) K-means clustering technique for adapting the centroids according to channel variations and for making classification decisions on the fly. -As a beneficial mechanism for estimating the location dependent noise plus interference power in the presence of secondary users mobility, we have proposed a model for estimating the noise plus interference power that can be used to adapt the centroid corresponding to the inactive PU.
We quantified the performance of the tracker in terms of mean square error and tracking accuracy whereas the performance of the enhanced spectrum sensing scheme is evaluated in terms of average detection probability, probability of false alarm and receiver operating characteristics.
The remainder of the paper is organized as follows. In Sect. 2 we provide the system model and assumptions. The process for the feature vectors realization, clustering algorithm and pilot based second order Kalman filter fading channel tracker are described in Sect. 3. Section 4 explains the quadratic polynomial regression model for the estimation of noise plus interference power. Numerical results are discussed in Sect. 5 followed by conclusion in Sect. 6.

System model and assumptions
We consider a CR network consisting of a PU transmitter (PU TX ), PU receivers (P-Rx) and a decentralized system of SUs. The SU system has multiple clusters, each consisting of a secondary base station (SBS) and J mobile SU receivers as illustrated in Fig. 1.
The PU transmitter is assumed to alternate between the ON (active) and OFF (inactive) states and its transmission power is fixed. The secondary networks share the same frequency channel with the primary system and the SUs in each cluster are assumed to be cooperating to jointly sense the PU's activities. The essence of the joint sensing is to overcome the limitations caused by shadowing and to take advantage of the spatial diversity of the SUs. It is assumed that PU TX periodically broadcasts to the PU receivers pilot signals which capture the channel state information of the paths traversed by the PU signals [30]. The attributes of the pilot signals are assumed to be known to the SUs so that the SUs can estimate the slowly fading channel between the PU TX and the SUs during the training period. However, we assume there is no explicit cooperation between PU TX and SU network, hence spectrum sensing is performed in an unsupervised manner.
Our proposed learning process is accomplished in two stages. During the first stage, the SUs within a particular cluster are assumed relatively stationary. They estimate the energy received from the transmitted PU signal and collectively report these to their serving SBS. The SBS uses the energy feature vectors to perform K-means clustering and to obtain centroids corresponding to ON and OFF status of PU TX . Since there are only two clusters, it is reasonable to assume that the cluster with its centre closer to the origin corresponds to OFF status while the other cluster corresponds to the ON status. There is no need for the labels of PU transmission and hence, there is no explicit cooperation between PU and SU networks.
In the second stage, it is assumed that SUs make random walks and the channels between the PU and the SUs are treated as time-variant as the SUs move around in the presence of scatterers within the coverage area of the PU. Due to this mobility, as the underlying channel varies, the energy that is supposed to be received during the PU active status varies. Hence, the cluster centroid corresponding to ON status of PU TX is updated continuously with the aid of the proposed Kalman filter channel tracking algorithm, that is, we estimate the predicted channel coefficients of each moving SUs and compute the estimate of the received energy and update the K-means cluster.
If the instantaneous PU signal received by the j 2 J SU is represented by x j ðnÞ, the goal of the spectrum sensing is to decide between two hypotheses, namely; the PU is absent which corresponds to the hypothesis H 0 , and the PU is present which corresponds to the hypothesis H 1 . This can be expressed as where s(n) is the signal transmitted by the PU, assumed to be an independent and identically distributed (i.i.d) process with zero mean and variance, E½jsðnÞj 2 ¼ r 2 s . The parameter c j ðnÞ is the non-frequency-selective channel coefficient between the PU and the j th SU which is assumed to be Rayleigh distributed. However, we estimate this parameter at every block interval using pilot signals. Hence, for the purpose of channel estimation, at every block interval, c j is treated as deterministic, but it is assumed to vary according to a Rayleigh distribution between the blocks. Finally, g j ðnÞ is an i.i.d circularly symmetric complex zero-mean additive white Gaussian noise (AWGN) with variance, E½jg j ðnÞj 2 ¼ r 2 g and E½: denotes the expectation operator.
3 Features realization, clustering algorithm and fading channel tracker

Energy vectors realization for SUs
During the learning interval, it is assumed that the PU alternates between H 0 and H 1 . So, the normalized energy estimated at j th SU node is given as where N s is the total number of samples taken during the observation window.
The energy feature, / j is a random variable whose probability density function (PDF) follows a central Chisquare distribution under H 0 and a non-central Chi-square distribution conditional on the channel gain under H 1 , both with 2N s degrees of freedom [33][34][35]. Furthermore, when N s is large enough (say, N s ' 250) [35], under H 0 this PDF can be approximated as Gaussian via the central limit theorem with mean, l 0 = r 2 g and variance, 17]. Similarly, the distribution of / j under H 1 and with large N s can be approximated as Gaussian with mean, l 1 = jc j j 2 r 2 s þ r 2 g and variance, When Kalman gain is used for estimating energy, we estimate the energy as where P Tu is the PU transmission power. For constructing energy clusters, Z energy samples, / j are collected and reported to the serving SBS by each of the cooperating SUs. The set of energy vectors obtained at the SBS from all SUs during the clustering period is represented by S ¼ f/ 1 ; :::; / Z g where / z 2 R M , z 2 Z and M J is the corresponding number of collaborating SUs. We employ Kmeans clustering that clusters the set S of energy vectors / z to two classes H 1 and H 0 . The centroid that is closer to the origin is deemed to be the centroid for H 0 and the other one

The K-means clustering algorithm
The K-means clustering algorithm partitions a set of data into K disjoint clusters. In our case, K is 2 and the algorithm aims to minimize the measure between the centroids of the clusters and the given set of data by iteratively assigning the data into the cluster with the closest centroid. The end result of the K-means algorithm is a set of clusters that are tightly packed and well-separated for every group.
Let S ¼ f/ z g Z z¼1 2 fH 1 ; H 0 g be the set of energy feature vectors obtained using (2). Also, let C k denote the set of training energy vectors that belong to cluster k which has a centroid h k defined as the mean of all the training energy vector in C k : S ¼ fC k g K k¼1 . According to [36], the Kmeans algorithm aims to obtain K clusters which minimize the within-cluster sum of squares as where card implies the cardinality function. The clustering algorithm in (4) is employed at the SBS to derive the required clusters' centroids. A generalized pseudo-code for implementing (4) is provided in Algorithm 1. The cluster centroid corresponding to H 0 can be identified as where h 0 is a zero vector, 0 2 R M representing the origin of the clusters' plane. It should be noted however, that when mobility of SUs is considered, the centroids obtained under H 1 must be continuously updated for reliable detection of PU's activities. In the Sub-Sect. 3.3, we describe a Kalman filter based scheme for tracking the gradual changes in PU-SU channel which would facilitate ''on the fly'' updating of the centroid-based decision plane for our spectrum sensing system.

Fading channel model and estimation
We consider the mobility of SUs in heavily built-up areas resulting in the presence of scatterers in their RF operating environment. Our aim is to estimate the slowly varying channel gain between the PU transmitter and SUs which can be modeled as flat fading Rayleigh channel [31,37] using the second order Kalman filter tracker. Let the discrete-time observation at the mobile SU terminal be given by where n is the symbol time index, bðnÞ is a known modulated pilot signal periodically broadcasted by PU TX , cðnÞ is the complex Gaussian channel coefficient whose variance is r 2 c , gðnÞ is the zero mean AWGN at the receiver with variance r 2 g . The noisy measurement of the channel gain from the pilot at the SU isĉ.
Referring to (1) and ignoring the user index j, the magnitude of the instantaneous channel gain, a = jcðnÞj is a Rayleigh distributed random variable. The autocorrelation of the fading channel coefficient R a ½l is defined for lag l as where J 0 ðÁÞ is the zeroth order Bessel function of the first kind, r 2 a = r 2 c , f d is the Doppler frequency and T s is the symbol interval. For Kalman filter based channel estimation, in addition to the noisy observationĉ, an approximate state-space representation of the evolution of the channel gain is required. For the first order Kalman tracker approach described in [26], the state-space model is described by where 0\d 1 \1 and nðnÞ is the complex AWGN with variance, r 2 n = ð1 À d 2 1 Þr 2 a . However, in the current work, we propose to use the second order autoregressive model (AR 2 ) being a more suitable approximation to the fading channel in the estimation problem being considered [37,40].
The AR 2 model with Gaussian assumption which approximates the flat fading Rayleigh channel can be expressed as a difference equation of the form where the model coefficients, d 1 and d 2 are the first and second order autoregressive coefficients respectively, whose values are obtained by imposing the correlation matching criterion as [37] ð10Þ which may be further re-expressed as and The remaining parameter nðnÞ in (9) is the AWGN with variance, r 2 n and n remains the symbol index. According to the Kalman filter framework, the tracking equations can be separated into two parts. These are the prediction equations that implement the prediction stage by means of appropriately chosen models and the correction equations for correcting the errors and inaccuracies arising in the estimation process with the aid of the measurement received [43].
The prediction stage: In general, the prediction equation that models a system's process evolution takes the form of the first order linear difference equation expressed as [26] aðnÞ ¼ Aaðn À 1Þ þ BuðnÞ þ wðn À 1Þ ð 13Þ where the associated parameters can be initialized with A = d 1 , B = 0 and u(n) = 0 while wðn À 1Þ is the process white noise. The measurement equation for correcting errors in the prediction is also described by an expression of the form where H is the measurement model which maps the true state space into the observed space and v(n) is a random variable representing the measurement noise. In this work, we adopted (6) as a form of (14) to derive the required multiple sequential measurements while (9) was used as the prediction equation that models the dynamics of the slowly fading channel. Furthermore, the state prediction error covariance is accounted for by the projection equation where P is the covariance of the state prediction error and the matrix A is defined by while ð:Þ T denotes the transpose operation. Also, the parameter Q in (15) is the prediction process noise covariance given by The correction stage: To correct the errors associated with the prediction of the channel gain process evolution using the measurements received from the PU's pilot, we first compute the Kalman gain as where the parameter R represents the measurement noise covariance whose value depends on the peculiarity of the PU-SU operating environment. For the second order approximation to Kalman filter, the matrix H is defined as [41] The optimal estimate of the PU-SU channel gain is then obtained by combining the Kalman gain in (18) with the measurement and prediction Eqs. (6) and (9) yielding the joint a posteriori estimate as From (20), the desired estimate of the magnitude of the channel gain is extracted aŝ It should be noted that a AR 1 ðnÞ and a AR 2 ðnÞ in (20) are the initial rough estimates obtained from the prediction process before correction is done by using the measurement update and I in (21) is the identity matrix. After each cycle of projection in time and measurement update, the process is repeated such that the a posteriori estimates become the initial values for the next cycle. Once the optimal estimate of the PU-SU channel gain has been obtained, the corresponding energy samples will be derived using (3).

Quadratic polynomial regression algorithm for energy prediction
As highlighted in Sect. 2, using the Kalman flter channel estimate, the energy at each SU receiver can be estimated using (3). In order to apply this, we need to have a reliable estimate of the noise power at the receiver, which can be estimated in the absence of any PU transmission. The noise power may include the co-channel interference due to any other transmitters in different cells. In the presence of mobility, the power of noise plus co-channel interference may vary due to various spatial locations. Estimation of this noise power would become non-trivial if the PU is also active. Hence, we propose a method of estimating the noise power during mobility of users and the PU transmission as below. Let / j;x and a j;x = jc j j denote the estimate of the energy and channel gain at various instants x for SU receiver j.
Hence for each SU j, we may relate the energy and channel gain as follows: where c 0;j is the noise plus interference power and c 2;j is the PU transmission power, both are the model parameters we seek to determine, e j;x is an unobserved random error with zero mean conditioned on variableâ j;x . Considering each SU receiver has N number of energychannel estimate samples corresponding to various time instants, the regression model in (23) can be expressed in terms of a response energy vector /, a design matrix A of estimated channel gain, a model parameter vector c, and random errors vector e aŝ   which may be written compactly as Using the ordinary least square estimation method [42], the vector of estimated polynomial regression coefficients is then obtained aŝ Hence, the first element of c j would provide the estimate of the noise plus interference power at the j th SU receiver. Once the estimates of the energy from all participating SUs are obtained from the estimates of the channel gains and the noise power as in Eq. (3), the centroid for H 1 cluster is updated as where b is an adaptation parameter, 0\b\1. In Algorithm 2, the proposed second-order Kalman filter based channel estimation technique applied within the context of spectrum sensing problem is presented.

Numerical results
For the Jake's fading model, the velocity of each SU is set to 6 km/hr. We generated 30,000 symbols received at the central frequency, f c of 2 GHz. Without loss of generality, we assume that the average power P avg: of the fading process is unity, but the proposed method would be applicable to any other value of P avg . The duration of the PU symbol is set to 33.33 ls transmitted at the central frequency, f c . The SUs are assumed to be equipped with omnidirectional antennas while the scatterers are assumed to cause the signals to be received at uniformly distributed angles of arrival. Furthermore, the PU is considered to alternate between the active and inactive states so that the number of clusters, K is set to 2. For the study, each energy sample is computed using the channel estimate which was based on varied pilot length, L. We first considered the case where the PU is cooperating with the SUs and assumed that the pilot is available all the time, only to compare the performance of the order-1 Kalman and order-2 filters. However, a more realistic case where there is no PU cooperation is considered for the rest of the simulations.
In Fig. 2, by 1000 Monte Carlo repetitions, we show the comparison between the estimation capabilities of the first order and the second order Kalman filter in terms of the mean square error (MSE) estimate of the fading over a range of PU's operating SNR from -8 dB to 30 dB and length of pilot symbols for channel estimation, L, equals 10. The PU was considered to be operating at 100%. It can be seen that for the proposed order-2 Kalman tracker, the MSE dropped appreciably from 0.164 to 1.08Â10 À4 as SNR is increased from -8 dB to 30 dB while a drop from 0.382 to 1.8Â10 À4 is observed for the order-1 tracker. Similarly, a shift from the order-1 to order-2 filter indicates about 57.89% reduction in estimation error at SNR of -8 dB, that is, the MSE dropped from 0.382 to 0.164 whereas at the SNR of 30 dB, a decrease of about 40.13 % is seen. These results clearly indicate that the proposed order-2 tracker is a better estimator compared to the order-1 tracker at all cases of SNR considered.
The tracking capability of the second order Kalman filter is presented in Fig. 3. Here, the tracker's performance is evaluated in terms of tracking accuracy over the duration of 6000 symbols at the PU's operating SNRs of -10 dB and 10 dB. The cases considered are both when the PU is operating at 100 percent duty cycle (the PU is active with pilot transmissions all the time) and when operating at 50 percent duty cycle (PU is active only 50% of the time alternating between ON and OFF) and when there is no cooperation with the secondary network. As observed throughout the observation window, the order-2 filter shows a good tracking performance and closely matches up with the true channel gain in both scenarios.
In Fig. 4, we show the comparison between the order-1 and order-2 filter based scheme in terms of their performance for detecting the presence of the PU under the assumption that there is a cooperation between the PU and SUs and the pilot is available all the time. It can be seen here that the order-2 based sensing scheme performs better that the order-1 scheme especially as the SNR is degraded. For example, although both schemes maintained zero average false alarm probability (Pfa ave ) throughout the SNR range considered, an average detection probability (Pd ave ) performance gain of about 4 % is offered by the order-2 scheme as the SNR is reduced from 16 dB to -4 dB. We opine that the improvement offered by the order-2 Kalman filter based scheme over the order-1 counterpart is likely as a result of its capability to capture the dynamics of the varying channel evolution better and the consequent decrease in the channel estimation errors as previously highlighted.
In Figs. 5 and 6, we focus on the performance of the two filters in a more realistic scenario where the PU and the SU networks do not cooperate and the PU operates at 50 % duty cycle. For the study, we set the pilot length, L to 10 and quantified the performance of both schemes in terms of   Pd ave and Pfa ave . In Fig. 5, about 3% detection performance improvement is seen to be offered by the order-2 scheme at the SNR of -4 dB where the Pd ave rises from about 0.71 to 0.74 with the use of order-2 based sensing scheme instead of the order-1 based scheme. Similar improvement in Pd ave can be observed within the SNR range of -4 dB to 6 dB. Similarly, the Pfa ave is reduced from 0.164 to 0.149 (about 9%) at -4 dB in favor of the order-2 based scheme over the order-1 scheme. This improvement in Pfa ave performance can also be observed throughout the SNR range of -4 dB to 6 dB.
We show the performance of the proposed order-2 filter based channel tracking spectrum sensing scheme under mobile SUs scenario and compare with the performance of the scheme where channel tracking is not used, that is, nontracking mode (NTM) in Fig. 7. The PU's operation is maintained at 50% duty cycle and L is kept at 10. As seen, the channel tracking based scheme exhibits superior performance in terms of Pd ave compared to the NTM system throughout the entire SNR range considered. It is especially noteworthy that as the SNR is improved from -10 dB to 20 dB, while the proposed channel tracking based scheme provides improvement in both Pd ave and Pfa ave , there is no noticeable performance gain in the NTM scheme. Furthermore, due to fading, the NTM scheme is characterized by significant missed detection. This indicates that the NTM scheme fails to cater to the challenges of PU detection in spectrum sensing for cognitive radio system when the SUs are mobile within the network and the PU-SU channel experiences aging/fading. It should be noted that for the NTM scheme, the centroids remain the same throughout the testing period, however, the true channel has been varied according to the Jakes fading model. Hence, the Pd ave and Pfa ave performance would be influenced by the initial channel values of the fading profile.
In Fig. 8, we study the effect of the variation in the duty cycle of the PU on the performance of the proposed channel tracking sensing scheme and the NTM scheme over the SNR range of -10dB to 20dB. It can be observed that as the duty cycle is increased from 25% to 75%, the performance of the proposed order-2 scheme improves as measured by Pd ave and Pfa ave . For example, at the SNR of -10dB, the proposed scheme offers an   14. An overall improvement in the performance of the channel tracking based scheme can also be observed as the SNR is improved from -10dB to 20dB. On the other hand, although the NTM scheme maintains zero Pfa ave throughout the range of SNR considered, its capability to detect the presence of the PU is reduced as the duty cycle of the PU is increased. The reduction in the NTM's capability for PU's detection is noticed across all SNR where Pd ave is maintained at approximately 0.26 and 0.75 for 75% to 25% duty cycle respectively. Fig. 9 shows the effect of varying the pilot length used for channel estimation on the performance of the proposed order-2 filter tracking based sensing scheme. It can be observed that as the pilot length is increased from 10 symbols to 30 symbols at 50% duty cycle operation of the PU, the system's performance is improved. For example, at the SNR of -10 dB, Pd ave is raised from 0.62 to 0.72 (about 16% gain) while the Pfa ave is reduced from 0.29 to 0.17 (about 41% reduction). It is noteworthy that at the SNR of 0 dB, the proposed scheme achieves zero Pfa ave and 0.96 Pd ave when the number of symbols in the pilot length is 30 indicating the capability of the proposed scheme to meet the performance requirement of CRs in terms of detection accuracy.
The receiver operating characteristics performance of the proposed second order Kalman filter channel gain estimation based sensing scheme is shown in Fig. 10. We considered 10 symbols and 20 symbols for the pilot length while the PU operating SNR of -5 dB, 0 dB and 5 dB were used for the investigation. It can be seen that the system's performance is enhanced both when the SNR improves and when the pilot length is increased. For example, at 5 dB, given that the operational duty cycle of the PU is kept at 50%, a Pd ave of 0.94 is achieved at the Pfa ave of 0.1. Furthermore, given the same PU operational duty cycle and Pfa ave of 0.1, an increase in Pd ave from 0.53 to 0.63, 0.78 to 0.9, and 0.94 to 0.98 are obtained at the SNR of -5 dB, 0 dB and 5 dB respectively as symbol size of the pilot used for the channel estimation is raised from 10 symbols to 20 symbols. These results further lend credence to the robustness of the proposed second order Kalman filter channel estimation method for spectrum sensing in CR networks of mobile SUs.
Finally, to conclude our investigations, we consider the implementation complexity of the proposed Kalman filter channel estimation based spectrum sensing scheme. The evaluation was done on an Intel core i5 À 5200U dual processors with CPU clocking speed of 2.20GHz and using MATLAB software version R2016b. Our comparison is based on the computation of average learning and prediction durations for implementing the first and second order channel gain trackers. For the same parameters as used in the previous simulations, the computation time averaged over 10,000 Monte Carlo runs for the first and he second order trackers are 18.97 ls and 22.64 ls respectively. As expected, the second order channel gain tracker has more computational complexity as compared to the first order channel gain tracker, however, this is compensated by the improvement in sensing performance as seen in the previous simulations.

Conclusion
In this paper, we studied the problem of spectrum sensing in cognitive radio networks of mobile secondary users. In particular, we considered the performance degradation that occurs when mobile secondary users are cooperatively sensing primary user's activities under flat fading channel conditions in the presence of scatterers. We presented a second order Kalman filter technique for real time estimation of the slowly fading channel gain between the primary user and secondary users. The estimated channel gain was then used to update the centroids of the K-means clustering based spectrum sensing scheme ''on the fly''. The performance of the proposed channel estimation technique was evaluated using mean square error estimation and tracking accuracy metrics while the performance of the ensuing spectrum sensing algorithm was quantified in terms of probability of detection, probability of false alarm and receiver operating characteristics. Simulation results show that the proposed techniques offer significant advantage for spectrum sensing in cognitive radio networks of mobile secondary users in comparison with previously proposed alternatives.
Acknowledgements Partial results of this paper were presented at the 10th Int. Conf. on Cognitive Radio Oriented Wireless Networks (CROWNCOM), Doha, Qatar, Apr. 2015 [26].

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.