Unsupervised Anomaly Detection

This topic introduces the unsupervised anomaly detection features for multivariate sample data available in Statistics and Machine Learning Toolbox™, and describes the workflows of the features for outlier detection (detecting anomalies in training data) and novelty detection (detecting anomalies in new data with uncontaminated training data).

For unlabeled multivariate sample data, you can detect anomalies by using isolation forest, one-class support vector machine (OCSVM), and Mahalanobis distance. These methods detect outliers either by training a model or by learning parameters. For novelty detection, you train a model or learn parameters with uncontaminated training data (data with no outliers) and detect anomalies in new data by using the trained model or learned parameters.

• Isolation forest — The Isolation Forest algorithm detects anomalies by isolating them from normal points using an ensemble of isolation trees. Detect outliers by using the iforest function, and detect novelties by using the object function isanomaly.

• One-class support vector machine (OCSVM) — One-Class Learning, or unsupervised SVM, tries to separate data from the origin in the transformed high-dimensional predictor space. Train an OCSVM model by using the fitcsvm function, and then detect outliers and novelties by using the resubPredict and predict object functions, respectively.

• Mahalanobis Distance — If sample data follows a multivariate normal distribution, then the squared Mahalanobis distances from samples to the distribution follow a chi-square distribution. Therefore, you can use the distances to detect anomalies based on the critical values of the chi-square distribution. For outlier detection, use the robustcov function to compute robust Mahalanobis distances. For novelty detection, you can compute Mahalanobis distances by using the robustcov and pdist2 functions.

Outlier Detection

This example illustrates the workflows of the three unsupervised anomaly detection methods (isolation forest, OCSVM, and Mahalanobis distance) for outlier detection.

Load Data

Load the humanactivity data set, which contains the variables feat and actid. The variable feat contains the predictor data matrix of 60 features for 24,075 observations, and the response variable actid contains the activity IDs for the observations as integers. This example uses the feat variable for anomaly detection.

load humanactivity

Find the size of the variable feat.

[N,D] = size(feat)

N = 24075

D = 60

Assume that the fraction of outliers in the data is 0.05.

contaminationFraction = 0.05;

Isolation Forest

Detect outliers by using the iforest function.

Train an isolation forest model by using the iforest function. Specify the fraction of outliers (ContaminationFraction) as 0.05.

rng("default") % For reproducibility
[forest,tf_forest,s_forest] = iforest(feat, ...
    ContaminationFraction=contaminationFraction);

forest is an IsolationForest object. iforest also returns the anomaly indicators (tf_forest) and anomaly scores (s_forest) for the data (feat). iforest determines the score threshold value (forest.ScoreThreshold) so that the function detects the specified fraction of observations as outliers.

Plot a histogram of the score values. Create a vertical line at the score threshold corresponding to the specified fraction.

figure
histogram(s_forest,Normalization="probability")
xline(forest.ScoreThreshold,"k-", ...
    join(["Threshold =" forest.ScoreThreshold]))
title("Histogram of Anomaly Scores for Isolation Forest")

Check the fraction of detected anomalies in the data.

OF_forest = sum(tf_forest)/N

OF_forest = 0.0496

The outlier fraction can be smaller than the specified fraction (0.05) when the scores can have tied values at the threshold.

One-Class Support Vector Machine (OCSVM)

Train an OCSVM model by using the fitcsvm function, and then detect outliers by using the resubPredict function.

Train a support vector machine model for one-class learning by using the fitcsvm function. The function trains a model for one-class learning if the class label variable is a vector of ones. Specify the fraction of outliers (OutlierFraction) as 0.05.

Mdl = fitcsvm(feat,ones(size(feat,1),1), ...
    OutlierFraction=contaminationFraction, ...
    KernelScale="auto",Standardize=true);

Mdl is a ClassificationSVM object.

Compute the outlier scores for feat by using the resubPredict function.

[~,s_OCSVM] = resubPredict(Mdl);

Negative score values indicate that the corresponding observations are outliers. Obtain the anomaly indicators.

tf_OCSVM = s_OCSVM < 0;

Plot a histogram of the score values.

figure
histogram(s_OCSVM,Normalization="probability");
xline(0,"k-","Threshold = 0")
title("Histogram of Anomaly Scores for OCSVM")

Check the fraction of detected anomalies in the data.

OF_OCSVM = sum(tf_OCSVM)/N

OF_OCSVM = 0.0500

fitcsvm trains the bias term of the SVM model so that the specified fraction of the training observations has negative scores. Therefore, OF_OCSVM is close to the specified fraction value.

Mahalanobis Distance

Use the robustcov function to compute robust Mahalanobis distances and robust estimates for the mean and covariance of the data.

Compute the Mahalanobis distance from feat to the distribution of feat by using the robustcov function. Specify the fraction of outliers (OutlierFraction) as 0.05. robustcov minimizes the covariance determinant over 95% of the observations.

[sigma,mu,s_robustcov,tf_robustcov_default] = robustcov(feat, ...
    OutlierFraction=contaminationFraction);

robustcov finds the robust covariance matrix estimate (sigma) and robust mean estimate (mu), which are less sensitive to outliers than the estimates from the cov and mean functions. The robustcov function also computes the Mahalanobis distances (s_robustcov) and the outlier indicators (tf_robustcov_default). By default, the function assumes that the data set follows a multivariate normal distribution, and identifies 2.5% of input observations as outliers based on the critical values of the chi-square distribution.

If the data set satisfies the normality assumption, then the squared Mahalanobis distance follows a chi-square distribution with D degrees of freedom, where D is the dimension of the data. In that case, you can find a new threshold by using the chi2inv function to detect the specified fraction of observations as outliers.

s_robustcov_threshold = sqrt(chi2inv(1-contaminationFraction,D));
tf_robustcov = s_robustcov > s_robustcov_threshold;

Create a distance-distance plot (DD plot) to check the multivariate normality of the data.

figure
d_classical = pdist2(feat,mean(feat),"mahalanobis");
gscatter(d_classical,s_robustcov,tf_robustcov,"kr",".x")
xline(s_robustcov_threshold,"k-")
yline(s_robustcov_threshold,"k-", ...
    join(["Threshold = " s_robustcov_threshold]));
l = refline([1 0]);
l.Color = "k";
xlabel("Mahalanobis Distance")
ylabel("Robust Distance")
legend("Normal Points","Outliers",Location="northwest")
title("Distance-Distance Plot")

Zoom in the axes to see the normal points.

xlim([0 10])
ylim([0 10])

If a data set follows a multivariate normal distribution, then data points cluster tightly around the 45 degree reference line. The DD plot indicates that the data set does not follow a multivariate normal distribution.

Because the data set does not satisfy the normality assumption, use the quantile of the distance values for the cumulative probability (1 — contaminationFraction) to find a threshold.

s_robustcov_threshold = quantile(s_robustcov,1-contaminationFraction);

Obtain the anomaly indicators for feat using the new threshold s_robustcov_threshold.

tf_robustcov = s_robustcov > s_robustcov_threshold;

Check the fraction of detected anomalies in the data.

OF_robustcov = sum(tf_robustcov)/N

OF_robustcov = 0.0500

Compare Detected Outliers

To visualize the detected outliers, reduce the data dimension by using the tsne function.

rng("default") % For reproducibility
T = tsne(feat,Standardize=true,Perplexity=100,Exaggeration=20);

Plot the normal points and outliers in the reduced dimension. Compare the results of the three methods: the isolation forest algorithm, OCSVM model, and robust Mahalanobis distance from robustcov.

figure
tiledlayout(2,2)
nexttile
gscatter(T(:,1),T(:,2),tf_forest,"kr",".x",[],"off")
title("Isolation Forest")
nexttile(3)
gscatter(T(:,1),T(:,2),tf_OCSVM,"kr",".x",[],"off")
title("OCSVM")
nexttile(4)
gscatter(T(:,1),T(:,2),tf_robustcov,"kr",".x",[],"off")
title("Robust Mahalanobis Distance")
l = legend("Normal Points","Novelties");
l.Layout.Tile = 2;

The novelties identified by the three methods are located near each other in the reduced dimension. Compute the fraction of outliers that the three methods have in common.

sum(tf_forest.*tf_OCSVM.*tf_robustcov)/N

ans = 0.0298

The three methods identify about 3% of the data (feat) as outliers.

You can also visualize observation values using the two most important features selected by the fsulaplacian function.

idx = fsulaplacian(feat);
figure
t = tiledlayout(2,2);
nexttile
gscatter(feat(:,idx(1)),feat(:,idx(2)),tf_forest,"kr",".x",[],"off")
title("Isolation Forest")
nexttile(3)
gscatter(feat(:,idx(1)),feat(:,idx(2)),tf_OCSVM,"kr",".x",[],"off")
title("OCSVM")
nexttile(4)
gscatter(feat(:,idx(1)),feat(:,idx(2)),tf_robustcov,"kr",".x",[],"off")
title("Mahalanobis Distance")
l = legend("Normal Points","Novelties");
l.Layout.Tile = 2;
xlabel(t,join(["Column" idx(1)]))
ylabel(t,join(["Column" idx(2)]))

Novelty Detection

This example illustrates the workflows of the three unsupervised anomaly detection methods (isolation forest, OCSVM, and Mahalanobis distance) for novelty detection.

Load Data

Load the humanactivity data set, which contains the variables feat and actid. The variable feat contains the predictor data matrix of 60 features for 24,075 observations, and the response variable actid contains the activity IDs for the observations as integers. This example uses the feat variable for anomaly detection.

load humanactivity

Partition the data into training and test sets by using the cvpartition function. Use 50% of the observations as training data and 50% of the observations as test data for novelty detection.

rng("default") % For reproducibility 
c = cvpartition(actid,Holdout=0.50);
trainingIndices = training(c); % Indices for the training set
testIndices = test(c); % Indices for the test set
XTrain = feat(trainingIndices,:);
XTest = feat(testIndices,:);

Assume that the training data is not contaminated (no outliers).

Find the size of the training and test sets.

[N,D] = size(XTrain)

N = 12038

D = 60

NTest = size(XTest,1)

NTest = 12037

Isolation Forest

Detect novelties using the object function isanomaly after training an isolation forest model by using the iforest function.

Train an isolation forest model.

[forest,tf_forest,s_forest] = iforest(XTrain);

forest is an IsolationForest object. iforest also returns the anomaly indicators (tf_forest) and anomaly scores (s_forest) for the training data (XTrain). By default, iforest treats all training observations as normal observations, and sets the score threshold (forest.ScoreThreshold) to the maximum score value.

Use the trained isolation forest model and the object function isanomaly to find novelties in XTest. The isanomaly function identifies observations with scores above the threshold (forest.ScoreThreshold) as novelties.

[tfTest_forest,sTest_forest] = isanomaly(forest,XTest);

The isanomaly function returns the anomaly indicators (tfTest_forest) and anomaly scores (sTest_forest) for the new data.

Plot histograms of the score values. Create a vertical line at the score threshold.

figure
histogram(s_forest,Normalization="probability")
hold on
histogram(sTest_forest,Normalization="probability")
xline(forest.ScoreThreshold,"k-", ...
    join(["Threshold =" forest.ScoreThreshold]))
legend("Training data","New data",Location="southeast")
title("Histograms of Anomaly Scores for Isolation Forest")
hold off

The anomaly score distribution of the test data is similar to that of the training data, so isanomaly detects a small number of anomalies in the test data.

Check the fraction of detected anomalies in the test data.

NF_forest = sum(tfTest_forest)/NTest

NF_forest = 8.3077e-05

Display the observation index of the anomalies in the test data.

idx_forest = find(tfTest_forest)

idx_forest = 3422

One-Class Support Vector Machine (OCSVM)

Train an OCSVM model by using the fitcsvm function, and then detect anomalies in new data by using the predict object function.

Train a support vector machine model for one-class learning by using the fitcsvm function. The function trains a model for one-class learning if the class label variable is a vector of ones. Specify the fraction of outliers (OutlierFraction) as 0.

Mdl = fitcsvm(XTrain,ones(size(XTrain,1),1),OutlierFraction=0, ...
    KernelScale="auto",Standardize=true);

Mdl is a ClassificationSVM object.

Compute the outlier scores for XTrain by using the resubPredict function.

[~,s_OCSVM] = resubPredict(Mdl);

Compute the novelty scores for XTest by using the predict function.

[~,sTest_OCSVM] = predict(Mdl,XTest);

Negative score values indicate that the corresponding observations are anomalies. Obtain the anomaly indicators.

tfTest_OCSVM = sTest_OCSVM < 0;

Plot histograms of the score values.

figure
histogram(s_OCSVM,Normalization="probability");
hold on
histogram(sTest_OCSVM,Normalization="probability");
xline(0,"k-","Threshold = 0")
legend("Training data","New Data",Location="best")
title("Histograms of Anomaly Scores for OCSVM")
hold off

Check the fraction of detected anomalies in the test data.

NF_OCSVM = sum(tfTest_OCSVM)/NTest

NF_OCSVM = 0.0707

Display the observation index of the anomalies in the test data.

idx_OCSVM = find(tfTest_OCSVM)

idx_OCSVM = 851×1

        1061
        1443
        1444
        1445
        1447
        1448
        1449
        1450
        1451
        1453
      ⋮

Mahalanobis Distance

Use the robustcov function to compute Mahalanobis distances of training data, and use the pdist2 function to compute Mahalanobis distances of new data.

Compute the Mahalanobis distance from XTrain to the distribution of XTrain by using the robustcov function. Specify the fraction of outliers (OutlierFraction) as 0.

[sigma,mu,s_mahal] = robustcov(XTrain,OutlierFraction=0);

robustcov also returns the estimates of covariance matrix (sigma) and mean (mu), which you can use to compute distances of new data.

Use the maximum value of s_mahal as the score threshold for novelty detection.

s_mahal_threshold = max(s_mahal);

Compute the Mahalanobis distance from XTest to the distribution of XTrain by using the pdist2 function.

sTest_mahal = pdist2(XTest,mu,"mahalanobis",sigma);

Obtain the anomaly indicators for XTest.

tfTest_mahal = sTest_mahal > s_mahal_threshold;

Plot histograms of the score values.

figure
histogram(s_mahal,Normalization="probability");
hold on
histogram(sTest_mahal,Normalization="probability");
xline(s_mahal_threshold,"k-", ...
    join(["Threshold =" s_mahal_threshold]))
legend("Training data","New Data",Location="southeast")
title("Histograms of Mahalanobis Distances")
hold off

Check the fraction of detected anomalies in the test data.

NF_mahal = sum(tfTest_mahal)/NTest

NF_mahal = 8.3077e-05

Display the observation index of the anomalies in the test data.

idx_mahal = find(tfTest_mahal)

idx_mahal = 3654

Compare Detected Anomalies

To visualize the detected anomalies, reduce the data dimension by using the tsne function.

rng("default") % For reproducibility 
T = tsne(feat,Standardize=true,Perplexity=100,Exaggeration=20);
XTest2D = T(testIndices,:);

Plot the normal points and novelties in the reduced dimension. Compare the results of the three methods: the isolation forest algorithm, OCSVM model, and Mahalanobis distance from mahal.

figure
tiledlayout(2,2)
nexttile
gscatter(XTest2D(:,1),XTest2D(:,2),tfTest_forest,"kr",".x",[],"off")
title("Isolation Forest")
nexttile(3)
gscatter(XTest2D(:,1),XTest2D(:,2),tfTest_OCSVM,"kr",".x",[],"off")
title("OCSVM")
nexttile(4)
gscatter(XTest2D(:,1),XTest2D(:,2),tfTest_mahal,"kr",".x",[],"off")
title("Mahalanobis Distance")
l = legend("Normal Points","Novelties");
l.Layout.Tile = 2;

The novelty identified by the isolation forest algorithm (idx_forest) and the novelty identified by the Mahalanobis distance (idx_mahal) are different, but they are located near each other in the reduced dimension.

Check if the novelties identified by OCSVM contain the novelties identified by the isolation forest algorithm and Mahalanobis distance.

ismember(idx_forest,idx_OCSVM)

ans = logical
   1

ismember(idx_mahal,idx_OCSVM)

ans = logical
   1

You can also visualize observation values using the two most important features selected by the fsulaplacian function.

idx = fsulaplacian([XTrain; XTest]);
figure
t = tiledlayout(2,2);
nexttile
gscatter(XTest(:,idx(1)),XTest(:,idx(2)),tfTest_forest,"kr",".x",[],"off")
title("Isolation Forest")
nexttile(3)
gscatter(XTest(:,idx(1)),XTest(:,idx(2)),tfTest_OCSVM,"kr",".x",[],"off")
title("OCSVM")
nexttile(4)
gscatter(XTest(:,idx(1)),XTest(:,idx(2)),tfTest_mahal,"kr",".x",[],"off")
title("Mahalanobis Distance")
l = legend("Normal Points","Novelties");
l.Layout.Tile = 2;
xlabel(t,join(["Column" idx(1)]))
ylabel(t,join(["Column" idx(2)]))