Learning low-dimensional manifolds under the L0-norm constraint for unsupervised outlier detection

Unsupervised outlier detection without the need for clean data has attracted great attention because it is suitable for real-world problems as a result of its low data collection costs. Reconstruction-based methods are popular approaches for unsupervised outlier detection. These methods decompose a data matrix into low-dimensional manifolds and an error matrix. Then, samples with a large error are detected as outliers. To achieve high outlier detection accuracy, when data are corrupted by large noise, the detection method should have the following two properties: (1) it should be able to decompose the data under the L0-norm constraint on the error matrix and (2) it should be able to reflect the nonlinear features of the data in the manifolds. Despite significant efforts, no method with both of these properties exists. To address this issue, we propose a novel reconstruction-based method: “L0-norm constrained autoencoders (L0-AE).” L0-AE uses autoencoders to learn low-dimensional manifolds that capture the nonlinear features of the data and uses a novel optimization algorithm that can decompose the data under the L0-norm constraints on the error matrix. This novel L0-AE algorithm provably guarantees the convergence of the optimization if the autoencoder is trained appropriately. The experimental results show that L0-AE is more robust, accurate and stable than other unsupervised outlier detection methods not only for artificial datasets with corrupted samples but also artificial datasets with well-known outlier distributions and real datasets. Additionally, the results show that the accuracy of L0-AE is moderately stable to changes in the parameter of the constrained term, and for real datasets, L0-AE achieves higher accuracy than the baseline non-robustified method for most parameter values.


Introduction
Outlier detection, which is the identification of samples that are different from the majority (inliers) in the data, is a fundamental problem in data mining and has been studied for a long time [1,16,26]. The demand for outlier detection has further as "corrupted data." Robust PCA (RPCA) [6], which is an improved method of PCA, linearly decomposes X into L and a sparse error matrix S. As the decomposition of RPCA is performed under an l 0 -norm constraint on S, RPCA is robust against corrupted data. However, because both methods decompose X linearly, there is the drawback of poor accuracy for data that have nonlinear features .
In recent years, the robust deep autoencoder (RDA) [34] has been proposed as an unsupervised outlier detection method that can accurately detect outliers, even for corrupted and nonlinear data. In the RDA, a low-dimensional nonlinear manifold that flexibly captures nonlinear features of input data is trained to reconstruct the data using autoencoders (AEs) [15] instead of the linear manifold used in RPCA. The RDA aims to decompose data as X = L D + S using an alternating optimization algorithm, where L D is a matrix whose elements are close on the low-dimensional nonlinear manifold and S is a sparse error matrix. To learn a low-dimensional nonlinear manifold that is not affected by corrupted elements, the RDA aims to decompose data under the l 0 -norm regularization on S; however, the RDA actually sparsifies S by solving the l 1 -norm relaxed regularization problem because of the computational intractability of l 0 -norm regularization. For corrupted data, l 0 -norm regularization, which can avoid the effect of the values of the corrupted elements, is more robust for outliers than l 1 -norm regularization, for which it is difficult to avoid the effects of the values of the corrupted elements (see Sect 5.3.1). Additionally, in the alternating optimization method of the RDA, no theoretical analysis of convergence is made. In practice, it has been confirmed that the progress of training the RDA model may be unstable (see Sect. 5.5).
In this paper, we propose L0-norm constrained AE (L0-AE), which is a novel reconstruction-based unsupervised outlier detection method that can learn low-dimensional manifolds under an l 0 -norm constraint for the error term using AE. In Table 1, we briefly summarize the comparison of the features of each reconstruction-based method. Compared with other reconstruction-based methods, L0-AE can provably guarantee the convergence of optimization under the l 0 -norm constraint in addition to considering nonlinear features. To the best of our knowledge, there is no reconstruction-based method that enforces neither relaxation of the l 0 -norm nor linearity. The key contributions of this work are as follows: 1. We propose a new alternating optimization algorithm that decomposes data nonlinearly under an l 0 -norm constraint for the error term (Sect. 3.1). 2. We prove that our alternating optimization algorithm always converges under the mild condition that AE is trained appropriately using gradient-based optimization instead of the RDA (Sect. 3.3).  Is trained appropriately 3. Through extensive experiments, we show that L0-AE is more robust and accurate than other reconstruction-based methods not only for corrupted data but also datasets that contain well-known outlier distributions. Additionally, for real datasets, we show that L0-AE achieves higher accuracy than the baseline non-robustified method for most parameter values and can learn more stably than the RDA (Sect. 5).

Preliminaries
In this section, we describe related reconstruction-based methods. Throughout the paper, we denote a given data matrix by X ∈ R N ×D , where N and D denote the number of samples and feature dimensions of X , respectively.

Robust principal component analysis
RPCA [6] is a modification of PCA [21]. First, we describe PCA. PCA decomposes a data matrix X into a low-rank matrix L and a small error matrix E such that This decomposition can be formulated as the following optimization problem where || · || 2 is the l 2 -norm and k is a parameter that determines the maximum value allowed for rank(L); that is, PCA seeks the best low-rank representation L of the given data matrix X with respect to (w.r.t.) the l 2 -norm. PCA is used for outlier detection because E • E indicates element-wise outlierness, where • is the element-wise product. Generally, the outlierness of each sample is obtained by adding E • E along the feature dimension. PCA is sensitive to outliers because of the use of the l 2 -norm, which makes the use of PCA for outlier detection difficult.
To avoid this drawback of PCA, RPCA decomposes X into a low-rank matrix L and a sparse error matrix S such that based on the assumption that there are not many corrupted elements in X . This decomposition problem can be represented as follows: minimize the rank of the matrix L satisfying the sparsity constraint ||S|| 0 ≤ k, where || · || 0 is the l 0 -norm that represents the number of nonzero elements and k is a parameter that determines the maximum value of ||S|| 0 . The Lagrangian reformulation of this optimization problem is where λ is a parameter that adjusts the sparsity of S instead of k. However, this non-convex optimization problem (4) is NP-hard [2]. A convex relaxation whose solution matches the solution of Eq. (4) under broad conditions has been proposed as follows: where || · || * is the nuclear norm and || · || 1 is the l 1 -norm (i.e., the sum of the absolute values of the elements). RPCA is also used for outlier detection because S • S indicates the outlierness of each element. RPCA is more robust against outliers than PCA because of the use of the l 0 -norm. However, RPCA cannot represent data distributed in a nonlinear manifold because the mapping from X to L is restricted to be linear.

Robust deep autoencoders
The RDA [34] is a method that relaxes the linearity assumption of RPCA. The RDA uses an AE for nonlinear mapping from X to L instead of the linear mapping of RPCA.
First, we describe the AE. The AE [15] is a special type of multilayer neural network in which the number of nodes in the output layer is the same as that in the input layer. Generally, the model parameters θ of AEs are trained to minimize the reconstruction error, which is defined by the following optimization problem: where f AE (X ; θ) is an output of the AE with an input X and parameters θ . If an AE with a bottleneck layer can be trained so that its reconstruction error becomes small, f AE (X ; θ) forms a low-dimensional nonlinear manifold; that is, the AE is a generalization of PCA and decomposes data such that whereL = f AE (X ; θ) denotes a matrix that forms a lowdimensional nonlinear manifold. The AE can be used as an outlier detection method that is a nonlinear version of PCA. Similar to PCA, the AE is sensitive to outliers because of the l 2 -error minimization.
Next, we describe the RDA. Similar to RPCA, the RDA aims to decompose X into X = L D + S, where S is a sparse error matrix whose nonzero elements indicate the reconstruction difficulty and L D is easily reconstructable data for the AE. This decomposition of the RDA is defined as the following optimization problem: This decomposition is more robust than Eq. (6) and can capture nonlinear features of X , unlike Eq. (4).
As Eq. (8) is difficult to optimize, an l 1 -relaxation of the following form has been also proposed [34]:

RDA for structured outliers
The outlier detection methods mentioned above assume that outliers are unstructured. However, in real applications, outliers are often structured [34]; that is, outliers are concentrated on a specific sample. To improve the detection accuracy for data that have structured outliers along the sample dimension, a method with grouped norm regularization has been proposed: where the l 2,1 -norm is defined as The RDA essentially uses Eq. (10) as an objective function for outlier detection.

Alternating optimization algorithm of the RDA
To optimize Eq. (10), the RDA applies an alternating optimization method to θ and S. In particular, the RDA uses a back-propagation method to train AE to minimize ||L D − f AE (L D ; θ)|| 2 with S fixed and a proximal-gradient-based method with a fixed step size to minimize the penalty term ||S T || 2,1 with θ fixed.

L0-norm constrained autoencoders
Although the RDA can detect outliers even for nonlinear data, there are several concerns with the RDA. First, because of the NP-hardness of the problem, the RDA uses the l 1 or l 2,1 -norm instead of the l 0 -norm, which makes the result highly sensitive to outliers. Second, the alternating optimization method of the RDA does not guarantee convergence. In practice, it has been experimentally confirmed that the progress of training the RDA model may be unstable (Sect. 5.5).
To address these issues, we propose an unsupervised outlier detection method that can decompose data nonlinearly using AE under an l 0 -norm constraint for the sparse matrix S. We prove that our algorithm always converges under a certain condition. For clarity, we first describe L0-AE for unstructured outliers and then extend it for structured outliers.

Formulation of the objective function
We consider decomposing a data matrix X into three terms: a low-dimensional nonlinear manifoldL = f AE (X ; θ), sparse error matrix S and small error matrix E similar to stable principal component pursuit [35], which overcomes the drawbacks of RPCA for noisy data, such that Note that we can calculate element-wise outlierness using (X − f AE (X ; θ)) • (X − f AE (X ; θ)) because X − AE θ (X ) equals S + E, which denotes the total error matrix. Equation (12) can be satisfied freely depending on E if E is not constrained. To train an AE that captures the features of X successfully, ||E|| 2 2 needs to be as small as possible. Therefore, we aim to obtain an AE model that successfully reconstructs the data by excluding the influence of outliers and a sparse matrix S by solving the following optimization problem that minimizes with l 0 -norm regularization of S: Specifically, this optimization problem (13) means "minimizing the reconstruction error of the AE without considering the sparse noisy elements S." Because of the l 0 -norm regularization of S, it is robust against outliers contained in X . However, Eq. (13) is difficult to optimize, as is Eq. (8) in the RDA. To facilitate the optimization described later, we consider the following l 0 -norm constrained optimization problem instead of Eq. (13): Equation (13) can be considered as the Lagrangian reformulation of this optimization problem (14). In unsupervised settings, it is not possible to know the optimal sparsity of S for a problem in advance; therefore, λ and k that adjust the sparsity of S are both hyperparameters. In practical situations, domain experts can often estimate the rate of outliers in the data. In Section 5.3.1, we experimentally show how to set the appropriate C p = k/N , which is the normalized value of k, for the data when the estimated rate of outliers is available. If we can solve Eq. (14), we can obtainL, which captures nonlinear features of X and can avoid the influence of outliers completely.

Alternating optimization algorithm
In the following, we propose an alternating optimization algorithm for θ and S for the l 0 -norm constrained optimization problem (14). We denote the reconstruction error terms X − f AE (X ; θ) by Z (θ ). Then the objective function can be expressed as ||Z (θ )−S|| 2 2 . In the optimization phase of θ with S fixed, we use a gradient-based method. In the optimization phase of S with θ fixed, the optimum S of ||Z (θ ) − S|| 2 2 can be obtained using the following method.
When ||S|| 0 ≤ 1 (S has at most one nonzero element), S that minimizes ||Z (θ ) − S|| 2 2 is the matrix that zeroes out the element with the largest absolute value of Z (θ ). Similarly, when ||S|| 0 ≤ k, S that minimizes ||Z (θ ) − S|| 2 2 is the matrix that zeroes out the elements with the top-k largest absolute values in Z (θ ); that is, S that minimizes Eq. (14) with θ fixed can be determined in a closed form as follows: where The optimality of Eq. (15) is proved by showing the following two facts.
(i) With H fixed, the optimal S is described by The optimal H consists of subscripts that have the top-k largest absolute values in (i) With H fixed, the objective function Eq. (14) is Note that, as |H | ≤ k, we do not need to consider the constraint. Hence, the optimal S is clearly given by H achieves an equal to or better objective value than S * H . This is because H can make more terms zero than H . This implies that if |H | is less than k, we can determine an equivalent or better H with k elements. 1 Hence, without loss of generality, we can assume |H | = k.
Consider an optimal H * . If we assume that there exists a subscript in H * that does not correspond to any of the top-k absolute values, then by replacing that subscript with a subscript corresponding to any of the top-k absolute values that is not already included in H * , we obtain an H that necessarily yields a lower value for the objective function. Therefore, H * is not optimal, which is a contradiction. Hence, the optimal H consists of subscripts that have the top-k largest absolute values in From the facts (i) and (ii), S that zeroes out the elements with the top-k largest absolute values in Z (θ ) with θ fixed is the optimal matrix that minimizes the objective function value; that is, S obtained by Eq. (15) always minimizes Eq. (14) with θ fixed.
Using the fact that the optimal S is described by s i j = z i j if (i, j) ∈ H * and s i j = 0 otherwise, we rewrite our proposed formulation (14) and alternating optimization method to make it algorithmically concise as follows: where A ∈ {0, 1} N ×D is a binary-valued matrix. In the alternating optimization of Eq. (17), θ is optimized using gradient-based optimization and A is optimized using If multiple elements in Z (θ ) match c, ties are broken arbitrarily so that ||A|| 0 ≥ ND − k.
The procedure of our proposed optimization algorithm is as follows: and E poch max ∈ N Initialize A ∈ R N ×D as a zero matrix, epoch counter E poch = 0 and an AE f AE ( · ; θ) with randomly initialized parameters. Repeat the following E poch max times: 2 2 using gradient-based optimization Return the element-wise outlierness R ∈ R N ×D , which is computed as follows: In step 3, the number of iterations in each gradient-based optimization process affects the performance of L0-AE. In practice, the values of a performance indicator and convergence of L0-AE are sufficient if the update occurs once based on a gradient-based optimization method (e.g., Adam [22]) (see Sect. 5). In this case, the total computational cost of L0-AE is the sum of the cost of a normal AE and sorting to obtain the top-k error value.

Convergence property
In this section, we prove that our alternating optimization algorithm always converges under the assumption that AE is trained appropriately using gradient-based optimization. Here, we denote the objective function ||A • Z (θ )|| 2 2 by K (A, θ) and the variables A and θ at the tth step of each alternating optimization phase by A t and θ t , respectively. Under this assumption, the convergence of the proposed alternating optimization method can be shown as follows: Proof By updating A using Eq. (18), the obtained A * minimizes Eq. (17) for any Z (θ ). Hence, for any θ t , we . This implies that a sequence {K (A t , θ t )} is a monotonically non-increasing and non-negative sequence. Therefore, by applying the monotone convergence theorem [4], there exists a value a ∞ = lim t→∞ K (A t , θ t ) ≥ 0.

Remark
The assumption K (A t+1 , θ t ) ≥ K (A t+1 , θ t+1 ) holds when the learning rate of the AE model is sufficiently small. Although this assumption might not hold for a fixed learning rate in practice, L0-AE results in better convergence than the RDA (see Sect. 5.5). Note that, by Theorem 2, our training algorithm is stable, although it does not guarantee that the parameters converge to a stationary point.

Algorithm for structured outliers
In the following, we describe an alternating optimization algorithm for data with structured outliers along the sample dimension. To detect structured outliers, Eqs. (17) and (18) are, respectively, reformulated as follows: where the l 2,0 -norm is defined as and c is the kth largest value of the vector D j=1 (z · j ) 2 . The sample-wise outlierness r is calculated using the R = [r i j ] defined by Eq. (19) as follows: L0-AE uses this version of the formulation and the alternating optimization method for outlier detection.
As with the update of A using Eq. (18), the update of A using Eq. (21) always minimizes the objective function (20) with θ fixed. The convergence of this algorithm using Eq. (21) is easily proved in a similar manner to Theorem 2.

Remark
The concept of our optimization methodology for structured outliers can be regarded as least trimmed squares (LTS) [27], in which the sum of all squared residuals except the largest k squared residuals is minimized.
have been proposed. However, they assume a different problem setting from ours; that is, the training data do not include anomalous data, and finding anomalies in test datasets is the target task. Therefore, these methods do not have a mechanism that excludes outliers during training. In [10], the equivalence of the global optimum of the VAE and RPCA is shown under the condition that a decoder has some type of affinity; however, connections between VAE and RPCA are not shown for general nonlinear activation functions.
The discriminative reconstruction AE (DRAE) [30] has been proposed for unsupervised outlier removal. The DRAE labels samples for which the reconstruction errors exceed a threshold as "outliers" and omits such samples for learning. To appropriately determine the threshold, the loss function of the DRAE has an additional term to separate the reconstruction errors of inliers and outliers. Because of this additional term, the DRAE does not solve an l 0 -norm constrained optimization problem; that is, the learned manifold is affected by outlier values, which degrades the detection performance (see Section 5).
RandNet [9] has been proposed as a method to increase robustness through an ensemble scheme. Although this ensemble may improve robustness, it does not completely avoid the adverse effects of outliers because each AE uses a non-regularized objective function. The deep structured energy-based model [31], which is a robust AE-based method that combines an energy-based model and non-regularized AE, has the same drawback. In [19], a method that simply combines AE and LTS was proposed; however, no theoretical analysis for the combined effects of AE and LTS was presented.

Datasets
We used both artificial and real datasets. Four types of artificial datasets were used to evaluate the detection performance against three types of well-known outliers [13,17]: global outliers, clustered outliers and local outliers, in addition to outliers contained in image data. Figure 1 illustrates the artificial datasets. In particular, the data in Fig. 1a have large noise, which can be regarded as an example of corrupted data. Figure 1a shows the dataset containing global outliers: We sampled 9,000 inlier samples (x, 2x 2 − 1) ∈ R 2 where x ∈ [−1, 1] was sampled uniformly. Furthermore, we sampled 1000 outliers uniformly from [−1, 1] × [−1, 1] that were regarded as corrupted samples. Figure 1b shows the two-dimensional dataset containing clustered outliers: We sampled 9000 inlier samples in the same manner as in Fig. 1a. Furthermore, we generated 1000 outliers with a normal distribution at the mean of (-0.1087, 0.6826) (randomly generated central coordinates) and the standard deviation of 0.04. Figure 1c shows the two-dimensional dataset containing the local outliers that were generated by referring to [17]. We sampled a total of 9750 inlier samples. Half of them were uniformly sampled from the center (0.5, 0.5), distance r ∈ [0, 0.5] and angle rad ∈ [0, 2π ] from the center. The other half were uniformly sampled from the center (-0.5, -0.5), r ∈ [0.4, 0.5] and rad ∈ [0, 2π ] from the center. Furthermore, we generated 250 outliers uniformly within a radius of 0.35 from the center (-0.5, -0.5). Figure 1d shows a part of the image dataset that was generated by referring to [34]. This dataset was generated using "MNIST": The inliers were 4750 randomly selected "4" images, and the outliers were 250 randomly selected images that were not "4." As real datasets, we used 11 datasets from Outlier Detection Datasets (ODDS) [25], which are commonly used as benchmarks for outlier detection methods. Additionally, we also used the "kdd99_rev" dataset introduced in [36]. Table 2 summarizes the 12 datasets. Before the experiments, we normalized the values of the datasets by dimension into the range of -1 to 1.

Evaluation method
Following the evaluation methods in [9,23,24], we compared the area under the receiver operating characteristic curves (AUCs) of the outlier detection accuracy. The evaluation was performed as follows: (1) all samples (whether inlier or outlier) were used for training; (2) the outlierness of each sample was calculated after training; and (3) the AUCs were calculated using outlierness and inlier/outlier labels. Note that we did not need to specify the detection threshold in this evaluation scheme.

Robust PCA (RPCA)
We used the RPCA implemented in [11] as a baseline linear method.

Normal autoencoders (N-AE)
We implemented N-AE with the loss function ||X − f AE (X ; θ)|| 2 2 as the baseline non-regularized detection method. For every AE-based method below, we used common network settings. We used three fully connected hidden layers (with a total of five layers), in which the number of neurons was ceil([D, D , where "linear" means that the layer is a fully connected layer with no activation function and "linear + relu" means that the layer is a fully connected layer with the rectified linear unit function. We set the mini-batch size to N /50 and applied Adam [22] (α = 0.001) for optimiza-tion with E poch max = 400. We used Eq. (23) to calculate the sample-wise outlierness. To prevent undue advantages to our method (L0-AE) and the other AE-based methods, we searched this architecture by maximizing the average AUC of N-AE.

RDA-Stbl
This baseline was used to confirm the effect of the l 0 -norm constraint of L0-AE. There are three differences between the objective functions of L0-AE and the RDA: (1) the l 0 -norm constrained problem versus the l 1 -norm regularized problem; (2) for the input to AEs, X versus X − S; and (3) in the first term of the objective function, || · || 2 2 versus || · || 2 . To bridge gaps (2) and (3), we introduce RDA-Stbl, which minimizes the objective function ||L D − f AE (X ; θ)|| 2 2 + λ||S T || 2,1 such that X − L D − S = 0, w.r.t. S and θ . By comparing the results of this model with those of L0-AE, we can confirm the effects that result from the difference between the l 0 -norm constraint and l 1 -norm regularization. We varied the λ values in 10 ways using the same values as those in the RDA and adopted the result with the highest average AUC.

L0-norm constrained autoencoders (L0-AE)
We used L0-AE for structured outliers (described in Sect. 3.4). The sample-wise outlierness of L0-AE was calculated using Eq. (23). We did not iterate the process to update the parameters of the AE at each gradient-based optimization step. Instead of k, we used C p = k/N (0 ≤ C p ≤ 1), which was normalized by the number of samples. This type of normalized parameter is often used in other methods such as the one-class support vector machine (OC-SVM) [29] and isolation forest (IForest) [24]. Note that L0-AE is equivalent to N-AE when C p = 0. We varied C p from 0.05 to 0.5 in 0.05 increments for each experiment and adopted the result with the highest average AUC.

Variational autoencoder (VAE)
We adopted VAE for our problem setting. Outlierness was computed using the reconstruction probability [3]. Note that the number of output dimensions of hidden layer1 and layer3 is twice that of the other AE-based methods.
We used Chainer (version 1.21.0) [8] for the implementation of the above AE-based methods. Additionally, we applied the following three conventional methods for a comparison of the AUCs for the real datasets.

One-class support vector machine (OC-SVM)
We used the OC-SVM implemented in scikit-learn and set kernel = "rbf."

Local outlier factor (LOF)
We used the LOF [5] implemented in scikit-learn and set "k" for the k-nearest neighbors to 100.

Isolation forest (IForest)
We used the IForest from a Python library pyod [33] with " n-estimators" (the number of trees) set to 100.
We tuned the above-mentioned parameters to achieve a high AUC on average over all the real datasets; for the other parameters, we used the recommended (default) values, unless otherwise noted.

Results for the global outlier dataset
First, we evaluated the AUCs and robustness against the outliers of L0-AE and the baseline reconstruction-based methods using the global outlier dataset that has highly corrupted samples (Fig. 1a) Table 3 presents the average of these measurements over 20 trials with different initial network weights, and Fig. 2 shows the distributions of outlierness for each method, except VAE, in the first run. Note that " params" in Table 3 indicates the parameters with which each method achieved the highest average AUC. As VAE uses the probability as outlierness, only the AUC was included in Table 3. These results show that L0-AE outperformed the other methods in terms of AUC and the distribution of outlierness between inliers and outliers (i.e., sparsity of the error matrix).
RPCA performed significantly poorer than the other methods because of its linearity. Among the AE-based methods, L0-AE performed best. In L0-AE, we observed that the learned manifold was almost entirely composed of inliers. Therefore, it can be confirmed that the l 0 -norm constraint of L0-AE functions as intended, and L0-AE can learn by almost completely eliminating the influence of corrupted samples while capturing nonlinear features. By contrast, the performances of the other AE-based methods were inferior to that of L0-AE because the other methods cannot completely exclude the influence of corrupted outliers. VAE was less accurate than the other AE-based methods; we consider that VAE was unable to demonstrate robustness because of the non-affinity of the decoder. For the DRAE, the reconstruction errors of inliers and outliers were relatively well separated, but the DRAE was more strongly affected by outliers than L0-AE because the DRAE objective function depends on how large outliers are, whereas the L0-AE objective function does not.
Next, we evaluated the parameter sensitivity of L0-AE using the artificial dataset with global outliers (Fig. 1a). We used 10 different datasets in which the outlier rate in the dataset was changed to {0.05, 0.1, 0.15, . . . , 0.45, 0.5}. The parameter sensitivity was evaluated based on the change in the average AUC over 10 trials with different random seeds when C p was changed from 0.05 to 0.55 in 0.05 increments for each dataset. Figure 3 shows the AUCs with different C p values for each dataset. From the results, we confirmed that, even if C p was changed significantly, the AUCs did not change significantly and remained above 80%. Therefore, we can say that the quality of the distribution of outlierness obtained by L0-AE is robust to changes in C p . For all results, the maximum AUC values occurred at C p values moderately greater than the true outlier rates. When C p was moderately greater than the true outlier rate, we considered that outliers were safely detected as outliers; therefore, we believe that the highest AUCs could be achieved. When C p was less than the true outlier rate, the AUCs were better than in the case of N-  Parameter sensitivity of L0-AE for the global outlier datasets with 10 outlier rates; each error bar represents a 95% confidence interval, gray vertical lines indicate the true outlier rates of each dataset, and gray horizontal lines indicate the AUCs of a normal AE (C p = 0) for each dataset AE (C p = 0) because some outliers were not reconstructed. Therefore, in practical use, when the true outlier rate can be approximately estimated by a domain expert, it is appropriate to set C p to a value moderately greater than the estimated outlier rate. On the basis of this experiment and experiments on various datasets presented below, we recommend setting C p to a value approximately 0.1 greater than the estimated outlier rate. By contrast, for the parameters in other robustified reconstruction-based methods, it is difficult to determine which value should be set even if the true outlier rate is accurately estimated by a domain expert because they are abstract in that they determine the ratio between the reconstruction error and the regularization term.

Results for the clustered, local and image outlier datasets
In the following, we evaluate L0-AE for the artificial datasets other than the global outlier dataset. First, we compare the AUCs and the distributions of outlierness. Table 4 presents the average of the measurements over 20 trials with different random seeds and the parameters for achieving the maximum average AUC for the three datasets (Figs. 1b-d). For the clustered and local outlier datasets, we used the same number of neurons used for the global outlier dataset. Figures 4, 5 and 6 show the distributions of outlierness for each method in the first run under the best parameters for the clustered, local and image outlier datasets, respectively.  For the clustered outlier dataset, we can see that L0-AE completely separated the distributions of outlierness of inliers and outliers. This is likely to be because of the appropriate C p settings and the l 0 -norm constraint, which allowed the manifold to be learned by excluding the effects of all clustered outliers. By contrast, the other methods had a lower average AUC than L0-AE because the learned manifolds were affected by the outliers as a result of the nonuse of the l 0 -norm constraint.
For the local and image outlier datasets, our L0-AE again had the best performance measures, as shown in Table 4, although the degree of improvement was smaller than that of the clustered outlierness experiments. A possible reason for this is that there was little room for improvement for the other robustified AE-based methods because the network structure, which was set up to provide the highest accuracy of unconstrained N-AE for fairness, was not suitable for the other methods. Thus, we verified that L0-AE can learn nonlinear manifolds better than other reconstruction-based methods for all cases, avoiding the influence of outliers.
Next, we evaluated the parameter sensitivity of L0-AE for the three datasets (Figs. 1b-d). Figure 7c shows the AUCs with different C p values for each dataset. As with the global outlier dataset, the AUCs were moderately robust to changes in C p , and the maximum AUC values were achieved at C p values moderately greater than the true outlier rates. A similar trend was observed for all well-known outliers, thereby confirming the utility of L0-AE. The development of an automatic optimal C p search method under the l 0 -norm constraint without ground truth outlier rates or labels is important future work. For the image outlier dataset, the AUC did not decrease as C p increased. Given the high AUC of RPCA, which is a linear method, for the image dataset, it is likely that the "4" images in the MNIST dataset were confined to a linear subspace and well separated from the images of the other classes. Therefore, even though C p was high and the number of inlier samples to be reconstructed reduced, we believe that L0-AE can learn manifolds such that the entire inlier sample can be successfully reconstructed.

Evaluation results for real datasets
First, we compared the AUCs for the real datasets. The AUC values were averaged over 50 trials with different random seeds. For the RDA, we used λ = 5.0 × 10 −5 ; for RDA-stbl, λ = 5.0 × 10 −4 ; for the DRAE, λ = 1.0 × 10 −1 ; and for L0-AE, C p = 0.3. Table 5 presents the average AUCs for each dataset; Avg. AUC, Avg. rank and Avg. time refer to the average AUC, average rank over the datasets and average run-time, respectively. L0-AE achieved the highest average AUC and average rank. Among the reconstruction-based methods, L0-AE achieved the highest AUCs for eight out of 12 datasets. Particularly on kdd99_rev, the AUC of L0-AE was considerably higher than those of the other AE-based methods. Because kdd99_rev had a high rate of outliers and they were distributed close to each other, the methods with l 1 -norm regularization and no regularization could not avoid reconstructing the outliers, whereas L0-AE almost completely avoided reconstruction because of its l 0 -norm constraint. Furthermore, we observed that the AUCs of the RDA and RDA-Stbl were nearly equal. This shows the importance of the l 0 -norm constraint. L0-AE outperformed the DRAE on average; we consider that L0-AE selectively reconstructs only the inliers, whereas the DRAE reconstructs inliers and reduces the variance of each label, thereby allowing outliers to affect manifolds. Additionally, the computational cost of the DRAE was higher than that of L0-AE, because of the calculation of the threshold. For VAE, training was unstable for some datasets. One possible reason is that VAE involves random sampling in the reparametrization trick, which increased the randomness of the results under these experimental settings. By contrast, among the AE-based methods, L0-AE achieved stable results. The RPCA results were relatively good for some datasets, which suggests that these datasets had linear features and l 0 -norm regularization worked; L0-AE performed well by capturing nonlinear features, even for the other datasets. The reason that RPCA outperformed some AE-based methods on average is that RPCA automatically detects the rank of the inlier, whereas the AE-based methods have a fixed latent dimension (there is no known method for obtaining an appropriate latent dimension in an unsupervised setting).
Next, we evaluated the parameter sensitivity of L0-AE using real datasets. Figure 8 shows the AUCs with different C p values for L0-AE (averaged over 50 trials). For most C p , the AUCs of L0-AE were higher than those of the normal AE (C p = 0), that is, the baseline of the non-robustified method. Therefore, even when the appropriate C p is unknown, L0-AE is still practical enough to be used. Additionally, as in the case of the artificial datasets, the maximum AUC values essentially occurred at C p values moderately greater than the true outlier rates. Therefore, in practical terms, it is easier to select an appropriate C p if the true outlier rate can be estimated.

Evaluation of convergence
We compared the convergence of L0-AE with that of the RDA. Here, we did not use mini-batch training to remove the effect of randomness. Table 6 presents the sum of the number of epochs in which the value of the objective function increased over the previous epoch during 20 trials (96,000 epochs in total). The results of N-AE are also included for reference. Additionally, Fig. 9 shows two transition examples of the values of the objective functions of the RDA and L0-AE. Among them, Fig. 9d shows the results of the only trial in which the objective function increased in L0-AE; the epochs in which the objective function increased were 294-302 when C p = 0.3, with an average increase of 0.23, which is considerably less than the value of the objective function. In Table 6 and Fig. 9, we observe that L0-AE converged well, regardless of the parameter C p , unlike the RDA. This empirically demonstrates the validity of Theorem 2, which states that our alternating optimization algorithm converges when the gradient-based optimization behaves ideally. For the RDA, when λ was small, the value of the objective function was unstable, but when λ was large, it was stable. We believe that this is because the sparsity of S improves as λ increases, thus reducing the gap in the objective function between phases of alternating optimization. We observe that, with N-AE, the values of the objective function did not increase, which implies that our gradient-based optimization generally satisfies the assumption in Theorem 2.

Conclusion
In this paper, we proposed L0-AE for unsupervised outlier detection. L0-AE decomposes data nonlinearly into lowdimensional manifolds that capture the nonlinear features of the data and a sparse error matrix under the l 0 -norm constraint. We proposed an efficient alternating optimization algorithm for training L0-AE and proved that this algorithm converges under a mild condition. We conducted extensive experiments with real and artificial datasets and confirmed that L0-AE is highly robust to outliers. We also confirmed that this high robustness leads to higher outlier detection accuracy, defined by the AUC, than those of existing unsupervised outlier detection methods.
Data Availability All data used in the paper are clearly marked with the sources or methods of preparation.

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://creativecomm ons.org/licenses/by/4.0/.