Anomaly detection for blueberry data using sparse autoencoder-support vector machine

Low-dimensional feature extraction

Given a d-dimension input sample X^d = [x₁, x₂, x_n]^d, X^d ∈ ℜ^D×N and n ≥1, d>1 is the dimensionality of the input sample, in order to achieve anomaly detection, firstly, the sparse autoencoder is used to capture the low-dimensional features F^l ={f₁, f₂,…, f_m}^l of X^d, l < d, and m ≥1 is the number of low-dimensional features.

The proposed sparse autoencoder consists of two hidden layers, an input layer and an output layer, as shown in Fig. 3, of which the first hidden layer is denoted as $H_1=\left[h_1^{\left(1\right)},h_1^{\left(2\right)},\dots ,h_1^{\left(m1\right)}\right]$, and the second hidden layer is denoted as, $H_2=\left[h_2^{\left(1\right)},h_2^{\left(2\right)},\dots ,h_2^{\left(m2\right)}\right]$, where H₁ ∈ ℜ^D×m1, H₂ ∈ ℜ^D×m2 i.e., $h_1^{\left(j\right)}\in {\Re }^{D\times 1},h_2^{\left(j\right)}\in {\Re }^{D\times 1}$ are the hidden representation of the input x_j, with j ≤n. Input layer and output layer. Input layer is used to receive the input x_j. The output layer sends out the corresponding reconstructed input ${\hat{x}}_j$.

Through using the backpropagation manner to update the objective function J(w₁, w₂, b₁, b₂), the error between x_j and the corresponding reconstructed ${\hat{x}}_j$ can be minimized. Therefore, J(w₁, w₂, b₁, b₂) in the proposed sparse autoencoder is given in Eq. (2). (2)${\displaystyle \left\{\begin{array}{c}J\left({\mathbf{w}}_1,{\mathbf{w}}_2,b_1,b_2\right)=\frac{1}{n}\sum _{i=1}^n\left|\right|x_i-{\hat{x}}_i|{|}^2+{\Omega }_{weight}|\left|{\mathbf{w}}_1\right|{|}^2+{\Omega }_{weight}\left|\right|{\mathbf{w}}_2|{|}^2+\sum _{j=1}^D\left[p\ast log\frac{p}{{\hat{p}}_j}+\left(1-p\right)\ast log\frac{\left(1-p\right)}{\left(1-{\hat{p}}_j\right)}\right]\left(2\mathrm{a}\right)\hfill \\ {\hat{p}}_j=\frac{1}{n}\sum _{i=1,k=1,2}^n{h}_k^{\left(i\right)}\left(j\right)\left(2\mathrm{b}\right)\hfill \end{array}\right..}$ where ${\hat{p}}_j$, P are the actual activation and average activation for the j-th neuron in the hidden layer consisting of D neurons, respectively. w₁ and b ₁ are the weight and bias in the encoder. Similarly, w₂ and b ₂ are the weight and bias in the decoder. Ω_weight is the sparse item for constraining the weights. For selection of Ω_weight, we consider the empirical Ω_weight proposed in Olshausen & Field (1997).

Separation of abnormal features

In the reconstructed feature space, there is a certain difference between normal features and abnormal features. So we use a support vector machine (SVM) to separate the abnormal features from normal ones. Given a dataset in feature space $\left({\hat{x}}_1,y_1\right),\left({\hat{x}}_2,y_2\right),\dots ,\left({\hat{x}}_N,y_N\right)$, where, ${\hat{x}}_i\in {\Re }^n$, i = 1, 2,.., N, $y_i\in +1,-1.{\hat{x}}_i$ is the i-th eigenvector. y_i is the learned class labels. The proposed SVM can be described as following, (3)${\displaystyle \begin{array}{c}\underset{w,b,{\xi }_i}{min}\left(\frac{1}{2}\left|\right|\mathbf{w}|{|}^2+C\sum _{i=1}^N{\xi }_i\right)\hfill \\ s.t.y_i\left(w^{\mathbf{T}}{\hat{x}}_i+b\right)\ge 1-{\xi }_i\hfill \\ {\xi }_i\ge 0\left(i=1,2,\dots ,N\right)\hfill \end{array}.}$

where ξ_i isa slack variable. C is a penalty item, and C>0.

This is a convex quadratic programming problem with inequality constraints, so that the dual problem can be obtained by using the Lagrange multiplier. We then transform the constrained objective function in Eq. (3) into an unconstrained newly constructed Lagrange function (Peng & Xu, 2013), as following (4)${\displaystyle L\left(w,b,{\alpha }_i\right)=\frac{1}{2}\left|\right|w|{|}^2-\sum _{i=1}^N{\alpha }_i\left(y_i\left(w\ast {\hat{x}}_i+b\right)-1\right)}$ where α_i is the Lagrange multiplier, and α_i > 0. In order to minimize L(w, b, α_i), let the partial derivatives of L(w, b, α_i) be zero with respect to w, b, respectively, as following (5)${\displaystyle \left.\begin{array}{c}\frac{\partial L\left(w,b,\alpha \right)}{\partial w}=0\to w=\sum _{i=1}^N{\alpha }_i{y}_i{\hat{x}}_i\hfill \\ \frac{\partial L\left(w,b,\alpha \right)}{\partial b}=0\to \sum _{i=1}^N{\alpha }_i{y}_i=0\hfill \end{array}\right\}.}$

Using Eq. (5) to solve Eq. (4), as follows (6)${\displaystyle \begin{array}{c}\underset{\alpha }{max}\sum _{i=1}^N{\alpha }_i-\frac{1}{2}\sum _{i=1}^N\sum _{j=1}^N{\alpha }_i{\alpha }_j{y}_i{y}_j\left({\hat{x}}_i\ast {\hat{x}}_j\right)\hfill \\ s.t.\sum _{i=1}^N{\alpha }_i{y}_i=0,\hfill \\ 0\le {\alpha }_i\le C,i=1,2,\dots ,N\hfill \end{array}.}$

Equation (6) can be represented by the kernel function $\kappa \left({\hat{x}}_i,{\hat{x}}_j\right)$, as follows (7)${\displaystyle \begin{array}{c}\underset{\alpha }{max}\sum _{i=1}^N{\alpha }_i-\frac{1}{2}\sum _{i=1}^N\sum _{j=1}^N{\alpha }_i{\alpha }_j{y}_i{y}_j\kappa \left({\hat{x}}_i,{\hat{x}}_j\right)\hfill \\ s.t.\sum _{i=1}^N{\alpha }_i{y}_i=0,\hfill \\ 0\le {\alpha }_i\le C,i=1,2,\dots ,N\hfill \end{array}.}$

κ(,) is a kernel function satisfying the Mercer theorem. Clearly, Eq. (7) has inequality constraints so the solutions must satisfy the Karush–Kuhn–Tucker (KKT) (Peng & Xu, 2013) conditions. (8)${\displaystyle \left.\begin{array}{c}{\alpha }_i\ge 0\hfill \\ y_i\left(w_i\cdot {\hat{x}}_i+b\right)-1\ge 0\hfill \\ {\alpha }_i\left(y_i\left(w_i\cdot {\hat{x}}_i+b\right)-1\right)=0\hfill \end{array}\right\}.}$

According to Eq. (5), the form of the solutions is given in Eq. (9). (9)${\displaystyle \left.\begin{array}{c}w=\sum _{i=1}^N{\alpha }_i{y}_i{\hat{x}}_i\hfill \\ \sum _{i=1}^N{\alpha }_i{y}_i=0\hfill \end{array}\right\}.}$

It can be seen that in α^∗, there is at least one ${\alpha }_j^{\ast }>0$ for which j has (10)${\displaystyle y_j\left(w^{\ast }\cdot {\hat{x}}_j+b^{\ast }\right)-1=0.}$

Hence, the optimal weight vector w^∗ andthe optimal bias b^∗ can be obtained as following, (11)${\displaystyle \left.\begin{array}{c}{w}^{\ast }=\sum _{i=1}^N{\alpha }_{{}_i}^{\ast }{y}_i{\hat{x}}_i\hfill \\ b^{\ast }=y_j-\sum _{i=1}^N{\alpha }_{{}_i}^{\ast }{y}_i\left({\hat{x}}_i\cdot {\hat{x}}_j\right)\hfill \end{array}\right\}.}$

Through learning the decision function, the separation between abnormal and normal features can be achieved in the feature space, therefore, the decision function f(x) is given in Eq. (12): (12)${\displaystyle f\left(x\right)=sign\left(\sum _{i=1}^N{y}_i{\alpha }_i^{\ast }\kappa \left(x_1,x_2\right)+b^{\ast }\right).}$

The kernel

The kernel κ(x₁, x₂) in Eq. (12) is used for separating abnormal features from normal features. From the Mercer theorem in the low-dimensional feature extraction section we know that κ(x₁, x₂) is the kernel function satisfying Mercer theorem. Therefore, we use the cumulative distribution function (Jayasumana et al., 2013) to get a positive definite kernel (13)${\displaystyle {\kappa }_f=\left(1-{\left(1-x_1^{\alpha }\right)}^{\beta },1-{\left(1-x_2^{\alpha }\right)}^{\beta }|\alpha ,\beta \right).}$ where α, β are the non-negative kernel parameters.

In addition, to improve the precision of separation between abnormal and normal features, it requires a kernel to be able to perceive the location of data points, that is, the kernel needs to satisfy two properties (Snoek et al., 2014), i.e., non-stationarity and to be flexible in controlling searches in the normal data region. Indeed, the Matern52 kernel in Snoek et al. (2014) is a continuous kernel (see Lemma 2) satisfying the two properties. Because it can make the radius to be warping concave and non-decreasing (Jayasumana et al., 2013; Snoek et al., 2014), more areas with small radii can be observed. The Matern52 kernel is given in Eq. (14). (14)${\displaystyle K_{M52}=\theta \left(1+\sqrt{C{r}^2\left(x_1,x_2\right)}+A{r}^2\left(x_1,x_2\right)\right)exp-\sqrt{B{r}^2\left(x_1,x_2\right)}.}$

where θ, r are the kernel parameter and kernel radius, respectively. A, B, C are constants, respectively.

Based on the above derivation, the kernel κ(x₁, x₂) in Eq. (12) can be derived (15)${\displaystyle \kappa \left(x_1,x_2\right)={\kappa }_{M52}\circ {\kappa }_f.}$

Since kernel κ_M52 and kernel κ_f are positive definite kernels, according to Lemma 3, the kernel κ(x₁, x₂) is also a positive definite kernel. Consequently, the derivation of the kernel κ(x₁, x₂) is completed.

The upper boundary of the number of anomalies

When separating anomalies from the normal points, we want to count the number of anomalies in order to detect them accurately, which indeed, is difficult. However, the Chebyshev theorem (Lemma 4) can estimate the upper limit for the number of anomalous features in the feature space because it is capable of determining the upper boundary of the percentage of data that exists within k number of standard deviations from the mean, meanwhile, without assuming the data distribution. Using Eq. (16), it can be estimated that the percentage of the number of anomalies is lower than 1/k².

Item C in Eq. (3) is a predefined penalty item, which determines the tolerated ratio that normal points are mistaken as anomalies. If C is set too large, it may increase the penalty, and vice versa. For the setting of C, we apply Eq. (1) to estimate it, i.e., let C = P, which can reduce the probability those normal points being mistaken for anomalies.

Learning rate.

The item is responsible for whether the objective function can converge to a local minimum. A suitable learning rate can make the objective function converge to a local minimum within a certain time. For selection of learning rate, we considered the reference value in Qu et al. (2021), i.e., 1e-7.

Activation function.

We selected the sigmoid function as the activation function. Because the output of the sigmoid function is either zero or 1, it is suitable for representing anomalies and normal instances.

The number of neurons.

According to the number of input samples, we dynamically adjust the number of neurons δ within a certain range, i.e., let δ₁ = 20, _δ2 = 100, and Δδ = 20, then, δ is determined using cross-validation.

Training for SA-SVM.

The overall algorithm of SA-SVM training is given in Algorithm 1. The number of neurons δ is determined by cross-validation in step 1 and step 13. The data set X^Cro_train is used to train SA-SVM, then data set X^Cro_val is used as cross-validation of parameter δ in order to obtain the optimal configuration of δ. Once the optimal configuration for δ is obtained, which is denoted as Opt(δ), the training for SA-SVM is started again, as shown from step 14 to step 19. During training, the error between the input and the reconstructed input is minimized by iteratively learning the objective function O(J, f), meanwhile, the back propagation technique is used to update the hyper parameters. The training is stopped when SA-SVM converges. The process in step 20 and step 23 shows that SA-SVM is well trained then we save the trained SA-SVM, and the final training accuracy is outputted.

Algorithm 1. Training for SA-SVM.

Input: parameters, iteration epoch T, δ₁, δ₂, Δδ, training set T_set.

Output: training accuracy Max_TAcc.

Begin:

1. T_set is divided into X^Cro_train, X^Cro_val;

2. for t = 1 to T do:

3. for δ = δ₁ to _δ2 with step Δδ do:

4. Use data set X^Cro_train to train SA-SVM( X^Cro_train; δ);

5. Learn objective function O (J, f);

6. Use backpropagation technique to update hyper parameters until they converge;

7. Calculate training accuracy T_acc = SA-SVM (X^Cro_train; δ);

8. Use data set X^Cro_val to verify SA-SVM( X^Cro_val; δ);

9. Calculate validation accuracy Val_acc (X^Cro_val; δ) = SA-SVM (X^Cro_val; δ);

10. end for

11. Select δ so that δ(max) = arg max(T_acc (δ(max)));

12. Get the optimal for the number of neurons Opt(δ) = δ(max);

13. end for

14. for t = 1 to T do:

15. Use training set T_set to train SA-SVM(T_set; Opt(δ));

16. Learn objective function O(J, f);

17. Use backpropagation manner to update hyper parameters until they converge;

18. Calculate training accuracy Train_Acc (t) = SA-SVM(T_set; Opt(δ); t);

19. end for

20. Select the t so that t_max = arg max (Train_Acc (t));

21. Get the maximum training accuracy in t_max-th iteration Max_TAcc = SA-SVM(T_set; Opt(δ); t_max);

22. Save the trained SA-SVM (T_set; Opt(δ); t_max);

23. Output the maximum training accuracy Max_TAcc

End